ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLPiNToMultiPionsChannel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLPiNToMultiPionsChannel.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France
33 //
34 #define INCLXX_IN_GEANT4_MODE 1
35 
36 #include "globals.hh"
37 
38 #include "G4INCLPiNToMultiPionsChannel.hh"
39 #include "G4INCLKinematicsUtils.hh"
40 #include "G4INCLBinaryCollisionAvatar.hh"
41 #include "G4INCLRandom.hh"
42 #include "G4INCLGlobals.hh"
43 #include "G4INCLLogger.hh"
44 #include <algorithm>
45 #include "G4INCLPhaseSpaceGenerator.hh"
46 
47 namespace G4INCL {
48 
49  const G4double PiNToMultiPionsChannel::angularSlope = 15.;
50 
51  PiNToMultiPionsChannel::PiNToMultiPionsChannel(const G4int npi, Particle *p1, Particle *p2)
52  : npion(npi),
53  ind2(0),
54  particle1(p1),
55  particle2(p2)
56  {
57  std::fill(isosp, isosp+4, 0);
58  }
59 
60  PiNToMultiPionsChannel::~PiNToMultiPionsChannel(){
61 
62  }
63 
64  void PiNToMultiPionsChannel::fillFinalState(FinalState *fs) {
65 
66 // assert(npion > 1 && npion < 5);
67 
68  Particle * nucleon;
69  Particle * pion;
70  if(particle1->isNucleon()) {
71  nucleon = particle1;
72  pion = particle2;
73  } else {
74  nucleon = particle2;
75  pion = particle1;
76  }
77  G4int ipi=ParticleTable::getIsospin(pion->getType());
78  ind2=ParticleTable::getIsospin(nucleon->getType());
79 
80  ParticleList list;
81  list.push_back(nucleon);
82  list.push_back(pion);
83  fs->addModifiedParticle(nucleon);
84  fs->addModifiedParticle(pion);
85 
86  isospinRepartition(ipi);
87 
88  const ParticleType tn=ParticleTable::getNucleonType(ind2);
89  nucleon->setType(tn);
90  ParticleType pionType=ParticleTable::getPionType(isosp[0]);
91  pion->setType(pionType);
92  const ThreeVector &rcolpion = pion->getPosition();
93  const ThreeVector zero;
94  for(G4int i=1; i<npion; ++i) {
95  pionType=ParticleTable::getPionType(isosp[i]);
96  Particle *newPion = new Particle(pionType,zero,rcolpion);
97  newPion->setType(pionType);
98  list.push_back(newPion);
99  fs->addCreatedParticle(newPion);
100  }
101 
102  const G4double sqrtS = KinematicsUtils::totalEnergyInCM(nucleon, pion);
103  PhaseSpaceGenerator::generateBiased(sqrtS, list, 0, angularSlope);
104 
105  }
106 
107  void PiNToMultiPionsChannel::isospinRepartition(G4int ipi) {
108  G4double rjcd=Random::shoot();
109  const G4int itot=ipi*ind2;
110 
111  isosp[1]=ipi;
112  if (npion != 3) {
113  if (npion == 4){
114  const G4double r2 = Random::shoot();
115  if (r2*3. > 2.) {
116  isosp[2]= 0;
117  isosp[3]= 0;
118  }
119  else {
120  isosp[2]= 2;
121  isosp[3]= -2;
122  }
123  }
124 
125  if (itot == 2) {
126  // CAS PI+ P ET PI- n
127  rjcd *= 5.;
128  if (rjcd > 3.) {
129  // PI+ PI+ N ET PI- PI- P
130  isosp[0]=2*ind2;
131  isosp[1]=ipi;
132  ind2=-ind2;
133  }
134  else {
135  // PI+ PI0 P ET PI- PI0 N
136  isosp[0]=0;
137  isosp[1]=ipi;
138  }
139  }
140  else if (itot == 0) {
141  // CAS PI0 P ET PI0 N
142  rjcd *= 90.;
143  if (rjcd > 13.) {
144  if (rjcd > 52.) {
145  // PI+ PI0 N ET PI- PI0 P
146  isosp[0]=2*ind2;
147  isosp[1]=0;
148  ind2=-ind2;
149  }
150  else {
151  // PI+ PI- P ET PI+ PI- N
152  isosp[1]=-2;
153  isosp[0]=2;
154  }
155  }
156  else {
157  // PI0 PI0 P ET PI0 PI0 N
158  isosp[0]=0;
159  isosp[1]=0;
160  }
161  }
162  else if (itot == -2) {
163  // CAS PI- P ET PI+ N
164  rjcd *= 45.;
165  if (rjcd > 17.) {
166  if (rjcd > 24.) {
167  // PI+ PI- N ET PI+ PI- P
168  isosp[0]=2*ind2;
169  ind2=-ind2;
170  }
171  else {
172  // PI0 PI0 N ET PI0 PI0 P
173  isosp[0]=0;
174  isosp[1]=0;
175  ind2=-ind2;
176  }
177  }
178  else
179  // PI- PI0 P ET PI+ PI0 N
180  isosp[0]=0;
181  }
182  } // if (npion != 3)
183  else {
184  // PI N -> PI PI PI N
185  // IF (IPI*IND2) 20,21,22
186  if (itot == -2) {
187  // CAS PI- P ET PI+ N
188  rjcd *= 135.;
189  if (rjcd <= 28.) {
190  // PI0 PI0 PI0 N ET PI0 PI0 PI0 P
191  isosp[0]=0;
192  isosp[1]=0;
193  isosp[2]=0;
194  ind2=-ind2;
195  }
196  else {
197  if (rjcd <= 84.) {
198  // PI+ PI- PI0 N ET PI- PI+ PI0 P
199  isosp[0]=2*ind2;
200  isosp[2]=0;
201  ind2=-ind2;
202  }
203  else {
204  if (rjcd <= 118.) {
205  // PI- PI- PI+ P ET PI- PI+ PI+ N
206  isosp[0]=ipi;
207  isosp[2]=-ipi;
208  }
209  else {
210  // PI- PI0 PI0 P ET PI0 PI0 PI+ N
211  isosp[0]=0;
212  isosp[2]=0;
213  }
214  }
215  }
216  }
217  else if (itot == 0) {
218  // CAS PI0 P ET PI0 N
219  rjcd *= 270.;
220  if (rjcd <= 39.) {
221  // PI0 PI0 PI0 P ET PI0 PI0 PI0 N
222  isosp[0]=0;
223  isosp[2]=0;
224  }
225  else {
226  if (rjcd <= 156.) {
227  // PI+ PI- PI0 P ET PI- PI+ PI0 N
228  isosp[0]=2;
229  isosp[2]=-2;
230  }
231  else {
232  if (rjcd <= 194.) {
233  // PI+ PI0 PI0 N ET PI0 PI0 PI- P
234  isosp[0]=0;
235  isosp[2]=2*ind2;
236  ind2=-ind2;
237  }
238  else {
239  // PI- PI+ PI+ N ET PI- PI- PI+ P
240  isosp[0]=2*ind2;
241  isosp[1]=2*ind2;
242  isosp[2]=-2*ind2;
243  ind2=-ind2;
244  }
245  }
246  }
247  }
248  else if (itot == 2) {
249  // CAS PI+ P ET PI- N
250  rjcd *= 5.;
251  if (rjcd <= 2.) {
252  // PI+ P PI0 PI0 ET PI- N PI0 PI0
253  isosp[0]=0;
254  isosp[2]=0;
255  }
256  else {
257  if (rjcd <= 3.) {
258  // PI+ P PI+ PI- ET PI- N PI+ PI-
259  isosp[0]=-2;
260  isosp[2]=2;
261  }
262  else {
263  // PI+ N PI+ PI0 ET PI- P PI0 PI-
264  isosp[0]=2*ind2;
265  isosp[2]=0;
266  ind2=-ind2;
267  }
268  }
269  }
270  }
271 
272  std::shuffle(isosp,isosp+npion,Random::getAdapter()); // isospin randomly distributed
273  }
274 
275 }