Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/inclxx/incl_physics/src/G4INCLNNToMissingStrangenessChannel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 }
246