ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLNDeltaToDeltaLKChannel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLNDeltaToDeltaLKChannel.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 
39 #include "G4INCLKinematicsUtils.hh"
41 #include "G4INCLRandom.hh"
42 #include "G4INCLGlobals.hh"
43 #include "G4INCLLogger.hh"
44 #include <algorithm>
46 
47 namespace G4INCL {
48 
50 
52  : particle1(p1), particle2(p2)
53  {}
54 
56 
59  const G4double maxDeltaMassRndm = std::atan((maxDeltaMass-ParticleTable::effectiveDeltaMass)*2./ParticleTable::effectiveDeltaWidth);
60  const G4double deltaMassRndmRange = maxDeltaMassRndm - ParticleTable::minDeltaMassRndm;
61 // assert(deltaMassRndmRange>0.);
62 
63  G4double y=ecm*ecm;
64  G4double q2=(y-1.157776E6)*(y-6.4E5)/y/4.0; // 1.157776E6 = 1076^2 = (mNucleon + mPion)^2, 6.4E5 = 800^2 = (mNucleon - mPion)^2
65  G4double q3=std::pow(std::sqrt(q2), 3.);
66  const G4double f3max=q3/(q3+5.832E6); // 5.832E6 = 180^3 = ???^3
67  G4double x;
68 
69  G4int nTries = 0;
70  G4bool success = false;
71  while(!success) { /* Loop checking, 10.07.2015, D.Mancusi */
72  if(++nTries >= 100000) {
73  INCL_WARN("NDeltaToDeltaLKChannel::sampleDeltaMass loop was stopped because maximum number of tries was reached. Minimum delta mass "
74  << ParticleTable::minDeltaMass << " MeV with CM energy " << ecm << " MeV may be unphysical." << '\n');
76  }
77 
78  G4double rndm = ParticleTable::minDeltaMassRndm + Random::shoot() * deltaMassRndmRange;
79  y = std::tan(rndm);
81 // assert(x>=ParticleTable::minDeltaMass && ecm >= x + ParticleTable::effectiveLambdaMass + ParticleTable::effectiveKaonMass + 1.0);
82 
83  // generation of the delta mass with the penetration factor
84  // (see prc56(1997)2431)
85  y=x*x;
86  q2=(y-1.157776E6)*(y-6.4E5)/y/4.0; // 1.157776E6 = 1076^2 = (mNucleon + mPion)^2, 6.4E5 = 800^2 = (mNucleon - mPion)^2
87  q3=std::pow(std::sqrt(q2), 3.);
88  const G4double f3=q3/(q3+5.832E6); // 5.832E6 = 180^3 = ???^3
89  rndm = Random::shoot();
90  if (rndm*f3max < f3)
91  success = true;
92  }
93  return x;
94  }
95 
97  // D++ p -> L K+ D++ (4)
98  //
99  // D++ n -> L K+ D+ (3)
100  // D++ n -> L K0 D++ (4)
101  //
102  // D+ p -> L K0 D++ (3)
103  // D+ p -> L K+ D+ (2)
104  //
105  // D+ n -> L K+ D0 (4)
106  // D+ n -> L K0 D+ (2)
107 
108  Particle *delta;
109  Particle *nucleon;
110 
111  if (particle1->isResonance()) {
112  delta = particle1;
113  nucleon = particle2;
114  }
115  else {
116  delta = particle2;
117  nucleon = particle1;
118  }
119 
120 
122 
124  const G4int iso_d = ParticleTable::getIsospin(delta->getType());
125  const G4double rdm = Random::shoot();
126 
127 /* const G4double m1 = particle1->getMass();
128  const G4double m2 = particle2->getMass();
129  const G4double pLab = KinematicsUtils::momentumInLab(particle1, particle2);*/
130 
131  ParticleType KaonType;
132  ParticleType DeltaType;
133  nucleon->setType(Lambda);
134 
135  if(std::abs(iso) == 4){// D++ p
136  KaonType = ParticleTable::getKaonType(iso/4);
137  DeltaType = ParticleTable::getDeltaType(3*iso/4);
138  }
139  else if(iso == 0){// D+ n
140  if(rdm*3 < 2){
141  KaonType = ParticleTable::getKaonType(iso_d);
142  DeltaType = ParticleTable::getDeltaType(-iso_d);
143  }
144  else{
145  KaonType = ParticleTable::getKaonType(-iso_d);
146  DeltaType = ParticleTable::getDeltaType(iso_d);
147  }
148  }
150  if(rdm*5 < 3){
151  KaonType = ParticleTable::getKaonType(-iso/2);
152  DeltaType = ParticleTable::getDeltaType(3*iso/2);
153  }
154  else{
155  KaonType = ParticleTable::getKaonType(iso/2);
156  DeltaType = ParticleTable::getDeltaType(iso/2);
157  }
158  }
159  else{// D++ n
160  if(rdm*7 < 3){
161  KaonType = ParticleTable::getKaonType(iso/2);
162  DeltaType = ParticleTable::getDeltaType(iso/2);
163  }
164  else{
165  KaonType = ParticleTable::getKaonType(-iso/2);
166  DeltaType = ParticleTable::getDeltaType(3*iso/2);
167  }
168  }
169 
170  delta->setType(DeltaType);
171  delta->setMass(sampleDeltaMass(sqrtS));
172 
173  ParticleList list;
174  list.push_back(delta);
175  list.push_back(nucleon);
176  const ThreeVector &rcol = nucleon->getPosition();
177  const ThreeVector zero;
178  Particle *kaon = new Particle(KaonType,zero,rcol);
179  list.push_back(kaon);
180 
183 
184 
185  fs->addModifiedParticle(delta);
186  fs->addModifiedParticle(nucleon);
187  fs->addCreatedParticle(kaon);
188 
189  }
190 }