ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLNNToMissingStrangenessChannel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLNNToMissingStrangenessChannel.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 "G4INCLNNToMissingStrangenessChannel.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 NNToMissingStrangenessChannel::angularSlope = 1.;
50 
51  NNToMissingStrangenessChannel::NNToMissingStrangenessChannel(Particle *p1, Particle *p2)
52  : particle1(p1), particle2(p2)
53  {}
54 
55  NNToMissingStrangenessChannel::~NNToMissingStrangenessChannel(){}
56 
57  void NNToMissingStrangenessChannel::fillFinalState(FinalState *fs) {
58 
59  const G4double sqrtS = KinematicsUtils::totalEnergyInCM(particle1, particle2); // MeV
60  const G4double pLab = 0.001*KinematicsUtils::momentumInLab(particle1, particle2); // GeV
61 // assert(sqrtS > 3100); // ! > 3047.34. Not supposed to be under 3,626 GeV.
62 
63  const G4int iso = ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
64 // assert(iso == -2 || iso == 0 || iso == 2);
65  G4int iso_system = 0;
66  G4int available_iso = 0;
67  G4int nbr_pions = 0;
68  G4int min_pions = 0;
69  G4int max_pions = 0;
70 
71  G4double rdm = Random::shoot();
72 
73  G4int nbr_particle = 3;
74 
75  if(rdm < 0.35){
76  // N-K-Lambda chosen
77  particle2->setType(Lambda);
78  available_iso = 2;
79  min_pions = 3;
80  max_pions = G4int((sqrtS-ParticleTable::getINCLMass(Lambda)-ParticleTable::getINCLMass(Proton)-ParticleTable::getINCLMass(KZero))/ParticleTable::getINCLMass(PiPlus));
81  }
82  else if((iso == 0 && rdm < 0.55) || rdm < 0.5){
83  // N-N-K-Kb chosen
84  nbr_particle++;
85  available_iso = 4;
86  min_pions = 1;
87  max_pions = G4int((sqrtS-2.*ParticleTable::getINCLMass(Proton)-2.*ParticleTable::getINCLMass(KZero))/ParticleTable::getINCLMass(PiPlus));
88  }
89  else{
90  // N-K-Sigma chosen
91  available_iso = 4;
92  min_pions = 3;
93  max_pions = G4int((sqrtS-ParticleTable::getINCLMass(SigmaMinus)-ParticleTable::getINCLMass(Proton)-ParticleTable::getINCLMass(KZero))/ParticleTable::getINCLMass(PiPlus));
94  }
95  /* Gaussian noise + mean value nbr pions fonction energy (choice)*/
96  G4double intermediaire = min_pions + Random::gauss(2) + std::sqrt(pLab-5);
97  nbr_pions = std::min(max_pions,std::max(min_pions,G4int(intermediaire )));
98 
99  available_iso += nbr_pions*2;
100  nbr_particle += nbr_pions;
101 
102  ParticleList list;
103  ParticleType PionType = PiZero;
104  const ThreeVector &rcol1 = particle1->getPosition();
105  const ThreeVector zero;
106  Particle *kaon = new Particle(KPlus,zero,rcol1);
107 
108  for(Int_t i=0; i<nbr_pions; i++){
109  Particle *pion = new Particle(PionType,zero,rcol1);
110  if(available_iso-std::abs(iso-iso_system) >= 4){ // pp(pn) pip:0.40(0.3) piz:0.35(0.4) pim:0.25(0.3)
111  rdm = Random::shoot();
112  if(((iso == 0) && rdm < 0.3) || ((iso == -2) && rdm < 0.40) || rdm < 0.25){
113  pion->setType(PiMinus);
114  iso_system -= 2;
115  available_iso -= 2;
116  }
117  else if(((iso == 0) && rdm < 0.7) || ((iso == -2) && rdm < 0.75) || rdm < 0.60){
118  pion->setType(PiZero);
119  available_iso -= 2;
120  }
121  else{
122  pion->setType(PiPlus);
123  iso_system += 2;
124  available_iso -= 2;
125  }
126  }
127  else if(available_iso-std::abs(iso-iso_system) == 2){
128  rdm = Random::shoot();
129  // pn pp too high (nn too low) -> PiMinus(PiPlus) pp too low (nn too high) -> PiPlus(PiMinus)
130  if(((iso == 0) && (rdm*0.7 < 0.3)) || ((rdm*0.60 < 0.25) && (Math::sign(iso-iso_system)*2-iso != 0)) || ((rdm*0.75 < 0.40) && (Math::sign(iso-iso_system)*2+iso != 0) && (iso != 0))){
131  pion->setType(ParticleTable::getPionType(Math::sign(iso-iso_system)*2));
132  iso_system += Math::sign(iso-iso_system)*2;
133  available_iso -= 2;
134  }
135  else{
136  PionType = PiZero;
137  available_iso -= 2;
138  }
139  }
140  else if(available_iso-std::abs(iso-iso_system) == 0){
141  pion->setType(ParticleTable::getPionType(Math::sign(iso-iso_system)*2));
142  iso_system += Math::sign(iso-iso_system)*2;
143  available_iso -= 2;
144  }
145  else INCL_ERROR("Pion Generation Problem in NNToMissingStrangenessChannel" << '\n');
146  list.push_back(pion);
147  }
148 
149  if(particle2->isLambda()){ // N-K-Lambda
150 // assert(available_iso == 2);
151  if(std::abs(iso-iso_system) == 2){
152  particle1->setType(ParticleTable::getNucleonType((iso-iso_system)/2));
153  kaon->setType(ParticleTable::getKaonType((iso-iso_system)/2));
154  }
155  else if(std::abs(iso-iso_system) == 0){
156  rdm = G4int(Random::shoot()*2.)*2-1;
157  particle1->setType(ParticleTable::getNucleonType(G4int(rdm)));
158  kaon->setType(ParticleTable::getKaonType(G4int(-rdm)));
159  }
160  else INCL_ERROR("Isospin non-conservation in NNToMissingStrangenessChannel" << '\n');
161  }
162  else if(min_pions == 1){ // N-N-K-Kb chosen
163 // assert(available_iso == 4);
164  Particle *antikaon = new Particle(KMinus,zero,rcol1);
165  if(std::abs(iso-iso_system) == 4){
166  particle1->setType(ParticleTable::getNucleonType((iso-iso_system)/4));
167  particle2->setType(ParticleTable::getNucleonType((iso-iso_system)/4));
168  kaon->setType(ParticleTable::getKaonType((iso-iso_system)/4));
169  antikaon->setType(ParticleTable::getAntiKaonType((iso-iso_system)/4));
170  }
171  else if(std::abs(iso-iso_system) == 2){ // choice: kaon type free, nucleon type fixed
172  rdm = G4int(Random::shoot()*2.)*2-1;
173  particle1->setType(ParticleTable::getNucleonType((iso-iso_system)/2));
174  particle2->setType(ParticleTable::getNucleonType((iso-iso_system)/2));
175  kaon->setType(ParticleTable::getKaonType(G4int(rdm)));
176  antikaon->setType(ParticleTable::getAntiKaonType(G4int(-rdm)));
177  }
178  else if(std::abs(iso-iso_system) == 0){ // particle1 3/4 proton, 1/4 neutron; particle2 1/4 proton, 3/4 neutron
179  rdm = G4int(Random::shoot()*2.)*2-1;
180  G4double rdm2 = G4int(Random::shoot()*2.)*2-1;
181  particle1->setType(ParticleTable::getNucleonType(std::max(rdm,rdm2)));
182  particle2->setType(ParticleTable::getNucleonType(std::min(rdm,rdm2)));
183  kaon->setType(ParticleTable::getKaonType(-G4int(rdm)));
184  antikaon->setType(ParticleTable::getAntiKaonType(-G4int(rdm2)));
185  }
186  else INCL_ERROR("Isospin non-conservation in NNToMissingStrangenessChannel" << '\n');
187  list.push_back(antikaon);
188  nbr_pions += 1; // need for addCreatedParticle loop
189  }
190  else{// N-K-Sigma
191 // assert(available_iso == 4);
192  if(std::abs(iso-iso_system) == 4){
193  particle1->setType(ParticleTable::getNucleonType((iso-iso_system)/4));
194  particle2->setType(ParticleTable::getSigmaType((iso-iso_system)/2));
195  kaon->setType(ParticleTable::getKaonType((iso-iso_system)/4));
196  }
197  else if(std::abs(iso-iso_system) == 2){ // choice: sigma quasi-free, kaon type semi-free, nucleon type fixed
198  rdm = Random::shoot();
199  // pp(pn) sp:0.42(0.23) sz:0.51(0.54) sm:0.07(0.23)
200  if(((iso == 0) && (rdm*0.77 < 0.23)) || ((rdm*0.58 < 0.07) && (Math::sign(iso-iso_system)*2-iso != 0)) || ((rdm*0.93 < 0.42) && (Math::sign(iso-iso_system)*2+iso != 0) && (iso != 0))){
201  particle2->setType(ParticleTable::getSigmaType(Math::sign(iso-iso_system)*2));
202  rdm = G4int(Random::shoot()*2.)*2-1;
203  particle1->setType(ParticleTable::getNucleonType(G4int(rdm)));
204  kaon->setType(ParticleTable::getKaonType(-G4int(rdm)));
205  }
206  else{
207  particle2->setType(SigmaZero);
208  particle1->setType(ParticleTable::getNucleonType((iso-iso_system)/2));
209  kaon->setType(ParticleTable::getKaonType((iso-iso_system)/2));
210  }
211  }
212  else if(std::abs(iso-iso_system) == 0){ // choice: sigma free, kaontype semi-free, nucleon type fixed
213  if(((iso == 0) && rdm < 0.23) || ((iso == 2) && rdm < 0.42) || rdm < 0.07){
214  particle2->setType(SigmaPlus);
215  particle1->setType(Neutron);
216  kaon->setType(KZero);
217  }
218  else if(((iso == 0) && rdm < 0.77) || ((iso == 2) && rdm < 0.93) || rdm < 0.58){
219  particle2->setType(SigmaZero);
220  rdm = G4int(Random::shoot()*2.)*2-1;
221  particle1->setType(ParticleTable::getNucleonType(G4int(rdm)));
222  kaon->setType(ParticleTable::getKaonType(-G4int(rdm)));
223  }
224  else{
225  particle2->setType(SigmaMinus);
226  particle1->setType(Proton);
227  kaon->setType(KPlus);
228  }
229  }
230  else INCL_ERROR("Isospin non-conservation in NNToMissingStrangenessChannel" << '\n');
231  }
232 
233  list.push_back(particle1);
234  list.push_back(particle2);
235  list.push_back(kaon);
236 
237  PhaseSpaceGenerator::generateBiased(sqrtS, list, list.size()-3, angularSlope);
238 
239  fs->addModifiedParticle(particle1);
240  fs->addModifiedParticle(particle2);
241  fs->addCreatedParticle(kaon);
242  for(Int_t i=0; i<nbr_pions; i++) fs->addCreatedParticle(list[i]);
243 
244  }
245 }