ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLNpiToMissingStrangenessChannel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLNpiToMissingStrangenessChannel.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 "G4INCLNpiToMissingStrangenessChannel.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 NpiToMissingStrangenessChannel::angularSlope = 1.;
50 
51  NpiToMissingStrangenessChannel::NpiToMissingStrangenessChannel(Particle *p1, Particle *p2)
52  : particle1(p1), particle2(p2)
53  {}
54 
55  NpiToMissingStrangenessChannel::~NpiToMissingStrangenessChannel(){}
56 
57  void NpiToMissingStrangenessChannel::fillFinalState(FinalState *fs) {
58 
59  const G4double sqrtS = KinematicsUtils::totalEnergyInCM(particle1, particle2); // MeV /!\/!\/!\.
60 // assert(sqrtS > 2.240); // ! > 2.109 Not supposed to be under 2.244 GeV.
61 
62  const G4int iso = ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
63 // assert(iso == -3 || iso == -1 || iso == 1 || iso == 3);
64  G4int iso_system = 0;
65  G4int available_iso = 0;
66  G4int nbr_pions = 0;
67  G4int min_pions = 0;
68  G4int max_pions = 0;
69 
70  Particle *pion_initial;
71  Particle *nucleon_initial;
72 
73  if(particle1->isPion()){
74  pion_initial = particle1;
75  nucleon_initial = particle2;
76  }
77  else{
78  pion_initial = particle2;
79  nucleon_initial = particle1;
80  }
81  const G4double pLab = 0.001*KinematicsUtils::momentumInLab(pion_initial, nucleon_initial); // GeV /!\/!\/!\.
82 
83  G4double rdm = Random::shoot();
84 
85  G4int nbr_particle = 2;
86 
87  if(rdm < 0.35){
88  // Lambda-K chosen
89  nucleon_initial->setType(Lambda);
90  available_iso = 1;
91  min_pions = 3;
92  max_pions = G4int((sqrtS-ParticleTable::getINCLMass(Lambda)-ParticleTable::getINCLMass(KZero)-10.)/ParticleTable::getINCLMass(PiPlus));
93  }
94  else if((iso == 0 && rdm < 0.55) || rdm < 0.5){
95  // N-K-Kb chosen
96  nbr_particle++;
97  available_iso = 3;
98  min_pions = 1;
99  max_pions = G4int((sqrtS-ParticleTable::getINCLMass(Proton)-2.*ParticleTable::getINCLMass(KZero)-10.)/ParticleTable::getINCLMass(PiPlus));
100  }
101  else{
102  // Sigma-K chosen
103  available_iso = 3;
104  min_pions = 3;
105  max_pions = G4int((sqrtS-ParticleTable::getINCLMass(SigmaMinus)-ParticleTable::getINCLMass(KZero)-10.)/ParticleTable::getINCLMass(PiPlus));
106  }
107  // Gaussian noise + mean value nbr pions fonction energy (choice)
108  G4double intermediaire = min_pions + Random::gauss(2) + std::sqrt(pLab-2.2);
109  nbr_pions = std::min(max_pions,std::max(min_pions,G4int(intermediaire )));
110 
111  available_iso += nbr_pions*2;
112  nbr_particle += nbr_pions;
113 
114  ParticleList list;
115  ParticleType PionType = PiZero;
116  const ThreeVector &rcol1 = pion_initial->getPosition();
117  const ThreeVector zero;
118 
119  // (pip piz pim) (sp sz sm) (L S Kb)
120  //pip_p 0.63 0.26 0.11 0.73 0.25 0.02 0.42 0.49 0.09 // inital
121  //pip_p 0.54 0.26 0.20 0.73 0.25 0.02 0.42 0.49 0.09 // choice
122  G4bool pip_p = (std::abs(iso) == 3);
123  //piz_p 0.32 0.45 0.23 0.52 0.40 0.08 0.40 0.41 0.19
124  G4bool piz_p = (ParticleTable::getIsospin(pion_initial->getType()) == 0);
125  //pim_p 0.18 0.37 0.45 0.20 0.63 0.17 0.39 0.33 0.28
126  G4bool pim_p = (!pip_p && !piz_p);
127 
128  for(Int_t i=0; i<nbr_pions; i++){
129  Particle *pion = new Particle(PionType,zero,rcol1);
130  if(available_iso-std::abs(iso-iso_system) >= 4){
131  rdm = Random::shoot();
132  if((pip_p && rdm < 0.54) || (piz_p && rdm < 0.32) || (pim_p && rdm < 0.45)){
133  pion->setType(ParticleTable::getPionType(G4int(Math::sign(iso))*2)); //pip/pip/pim
134  iso_system += 2*G4int(Math::sign(iso));
135  available_iso -= 2;
136  }
137  else if((pip_p && rdm < 0.80) || (piz_p && rdm < 0.77) || (pim_p && rdm < 0.82)){
138  pion->setType(PiZero);
139  available_iso -= 2;
140  }
141  else{
142  pion->setType(ParticleTable::getPionType(-G4int(Math::sign(iso))*2));
143  iso_system -= 2*G4int(Math::sign(iso));
144  available_iso -= 2;
145  }
146  }
147  else if(available_iso-std::abs(iso-iso_system) == 2){
148  rdm = Random::shoot();
149  if((pip_p && rdm < 0.26/0.37 && (Math::sign(iso)*Math::sign(iso-iso_system)+1)) || (pip_p && rdm < 0.26/0.89 && (Math::sign(iso)*Math::sign(iso-iso_system)-1)) ||
150  (piz_p && rdm < 0.45/0.68 && (Math::sign(iso)*Math::sign(iso-iso_system)+1)) || (piz_p && rdm < 0.45/0.77 && (Math::sign(iso)*Math::sign(iso-iso_system)-1)) ||
151  (pim_p && rdm < 0.37/0.82 && (Math::sign(iso)*Math::sign(iso-iso_system)+1)) || (piz_p && rdm < 0.37/0.55 && (Math::sign(iso)*Math::sign(iso-iso_system)-1))){
152  pion->setType(PiZero);
153  available_iso -= 2;
154  }
155  else{
156  pion->setType(ParticleTable::getPionType(Math::sign(iso-iso_system)*2));
157  iso_system += Math::sign(iso-iso_system)*2;
158  available_iso -= 2;
159  }
160  }
161  else if(available_iso-std::abs(iso-iso_system) == 0){
162  pion->setType(ParticleTable::getPionType(Math::sign(iso-iso_system)*2));
163  iso_system += Math::sign(iso-iso_system)*2;
164  available_iso -= 2;
165  }
166  else INCL_ERROR("Pion Generation Problem in NpiToMissingStrangenessChannel" << '\n');
167  list.push_back(pion);
168  }
169 
170  if(nucleon_initial->isLambda()){ // K-Lambda
171 // assert(available_iso == 1);
172  pion_initial->setType(ParticleTable::getKaonType(iso-iso_system));
173  }
174  else if(min_pions == 1){ // N-K-Kb chosen
175 // assert(available_iso == 3);
176  Particle *antikaon = new Particle(KMinus,zero,rcol1);
177  if(std::abs(iso-iso_system) == 3){
178  pion_initial->setType(ParticleTable::getKaonType((iso-iso_system)/3));
179  nucleon_initial->setType(ParticleTable::getNucleonType((iso-iso_system)/3));
180  antikaon->setType(ParticleTable::getAntiKaonType((iso-iso_system)/3));
181  }
182  else if(std::abs(iso-iso_system) == 1){ // equi-repartition
183  rdm = G4int(Random::shoot()*3.)-1;
184  nucleon_initial->setType(ParticleTable::getNucleonType((G4int(rdm+0.5)*2-1)*(iso_system-iso)));
185  pion_initial->setType(ParticleTable::getKaonType((std::abs(rdm*2)-1)*(iso-iso_system)));
186  antikaon->setType(ParticleTable::getAntiKaonType((G4int(rdm-0.5)*2+1)*(iso-iso_system)));
187  }
188  else INCL_ERROR("Isospin non-conservation in NNToMissingStrangenessChannel" << '\n');
189  list.push_back(antikaon);
190  nbr_pions += 1; // need for addCreatedParticle loop
191  }
192  else{// Sigma-K
193 // assert(available_iso == 3);
194  if(std::abs(iso-iso_system) == 3){
195  pion_initial->setType(ParticleTable::getKaonType((iso-iso_system)/3));
196  nucleon_initial->setType(ParticleTable::getSigmaType((iso-iso_system)*2/3));
197  }
198  else if(std::abs(iso-iso_system) == 1){
199  rdm = Random::shoot();
200  if((pip_p && rdm < 0.73) || (piz_p && rdm < 0.32) || (pim_p && rdm < 0.45)){
201  nucleon_initial->setType(SigmaZero);
202  pion_initial->setType(ParticleTable::getKaonType(iso-iso_system));
203  }
204  else{
205  nucleon_initial->setType(ParticleTable::getSigmaType((iso-iso_system)*2));
206  pion_initial->setType(ParticleTable::getKaonType(iso_system-iso));
207  }
208  }
209  else INCL_ERROR("Isospin non-conservation in NNToMissingStrangenessChannel" << '\n');
210  }
211 
212  list.push_back(pion_initial);
213  list.push_back(nucleon_initial);
214 
215  PhaseSpaceGenerator::generateBiased(sqrtS, list, list.size()-1, angularSlope);
216 
217  fs->addModifiedParticle(pion_initial);
218  fs->addModifiedParticle(nucleon_initial);
219  for(Int_t i=0; i<nbr_pions; i++) fs->addCreatedParticle(list[i]);
220 
221  }
222 }