ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLNDeltaToDeltaSKChannel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLNDeltaToDeltaSKChannel.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 "G4INCLNDeltaToDeltaSKChannel.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 NDeltaToDeltaSKChannel::angularSlope = 2.;
50 
51  NDeltaToDeltaSKChannel::NDeltaToDeltaSKChannel(Particle *p1, Particle *p2)
52  : particle1(p1), particle2(p2)
53  {}
54 
55  NDeltaToDeltaSKChannel::~NDeltaToDeltaSKChannel(){}
56 
57  G4double NDeltaToDeltaSKChannel::sampleDeltaMass(G4double ecm) {
58  const G4double maxDeltaMass = ecm - ParticleTable::effectiveSigmaMass - ParticleTable::effectiveKaonMass - 1.0;
59  const G4double maxDeltaMassRndm = std::atan((maxDeltaMass-ParticleTable::effectiveDeltaMass)*2./ParticleTable::effectiveDeltaWidth); // atan((mass-1232)*2/130)
60  const G4double deltaMassRndmRange = maxDeltaMassRndm - ParticleTable::minDeltaMassRndm; // atan
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 // (prc56(1997)2431) (eq 3.7)^2
65  G4double q3=std::pow(std::sqrt(q2), 3.);
66  const G4double f3max=q3/(q3+5.832E6); // 5.832E6 = 180^3 = cut_parameter^3 // (prc56(1997)2431) (cf eq 3.6)
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("NDeltaToDeltaSKChannel::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');
75  return ParticleTable::minDeltaMass;
76  }
77 
78  G4double rndm = ParticleTable::minDeltaMassRndm + Random::shoot() * deltaMassRndmRange; // atan in order to avec a distribution in 1/(1+x^2)
79  y = std::tan(rndm); // (mass-1232)*2/130
80  x = ParticleTable::effectiveDeltaMass + 0.5*ParticleTable::effectiveDeltaWidth*y; // probability to have the mass M = 1/(1+(M-1232)^2)/Pi cut with min and max mass
81 // assert(x>=ParticleTable::minDeltaMass && ecm >= x + ParticleTable::effectiveSigmaMass + 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 // (prc56(1997)2431) (eq 3.7)^2
87  q3=std::pow(std::sqrt(q2), 3.);
88  const G4double f3=q3/(q3+5.832E6); // 5.832E6 = 180^3 = cut_parameter^3 // (prc56(1997)2431) (eq 3.6)
89  rndm = Random::shoot();
90  if (rndm*f3max < f3) success = true; // promoting high masses
91  }
92  return x;
93  }
94 
95  void NDeltaToDeltaSKChannel::fillFinalState(FinalState *fs) {
96  //
97  // D++ p -> S+ K+ D+ (2)
98  // D++ p -> S0 K+ D++ (1)
99  // D++ p -> S+ K0 D++ (6)
100  //
101  // D++ n -> S+ K+ D0 (2)
102  // D++ n -> S0 K+ D+ (4)
103  // D++ n -> S- K+ D++ (6)
104  // D++ n -> S+ K0 D+ (2)
105  // D++ n -> S0 K0 D++ (1)
106  //
107  // D+ p -> S+ K+ D0 (2)
108  // D+ p -> S0 K+ D+ (1)
109  // D+ p -> S- K+ D++ (2)
110  // D+ p -> S+ K0 D+ (2)
111  // D+ p -> S0 K0 D++ (4)
112  //
113  // D+ n -> S+ K+ D- (2)
114  // D+ n -> S0 K+ D0 (4)
115  // D+ n -> S- K+ D+ (2)
116  // D+ n -> S+ K0 D0 (2)
117  // D+ n -> S0 K0 D+ (1)
118  // D+ n -> S- K0 D++ (2)
119 
120 
121  Particle *delta;
122 
123  if (particle1->isResonance()) {
124  delta = particle1;
125  }
126  else {
127  delta = particle2;
128  }
129 
130  const G4double sqrtS = KinematicsUtils::totalEnergyInCM(particle1, particle2);
131 
132  const G4int iso = ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
133  const G4int iso_d = ParticleTable::getIsospin(delta->getType());
134 
135  ParticleType KaonType;
136  ParticleType DeltaType;
137  ParticleType SigmaType;
138 
139  const G4double rdm = Random::shoot();
140 
141  if(std::abs(iso) == 4){// D++ p
142  if(rdm*9 < 2){
143  KaonType = ParticleTable::getKaonType(iso/4);
144  DeltaType = ParticleTable::getDeltaType(iso/4);
145  SigmaType = ParticleTable::getSigmaType(iso/2);
146  }
147  else if(rdm*9 < 3){
148  KaonType = ParticleTable::getKaonType(iso/4);
149  DeltaType = ParticleTable::getDeltaType(3*iso/4);
150  SigmaType = SigmaZero;
151  }
152  else{
153  KaonType = ParticleTable::getKaonType(-iso/4);
154  DeltaType = ParticleTable::getDeltaType(3*iso/4);
155  SigmaType = ParticleTable::getSigmaType(iso/2);
156  }
157  }
158  else if(iso == 0){// D+ n
159  if(rdm*13 < 2){
160  KaonType = ParticleTable::getKaonType(iso_d);
161  DeltaType = ParticleTable::getDeltaType(-3*iso_d);
162  SigmaType = ParticleTable::getSigmaType(2*iso_d);
163  }
164  else if(rdm*13 < 6){
165  KaonType = ParticleTable::getKaonType(iso_d);
166  DeltaType = ParticleTable::getDeltaType(-iso_d);
167  SigmaType = SigmaZero;
168  }
169  else if(rdm*13 < 8){
170  KaonType = ParticleTable::getKaonType(iso_d);
171  DeltaType = ParticleTable::getDeltaType(iso_d);
172  SigmaType = ParticleTable::getSigmaType(-2*iso_d);
173  }
174  else if(rdm*13 < 10){
175  KaonType = ParticleTable::getKaonType(-iso_d);
176  DeltaType = ParticleTable::getDeltaType(-iso_d);
177  SigmaType = ParticleTable::getSigmaType(2*iso_d);
178  }
179  else if(rdm*13 < 11){
180  KaonType = ParticleTable::getKaonType(-iso_d);
181  DeltaType = ParticleTable::getDeltaType(iso_d);
182  SigmaType = SigmaZero;
183  }
184  else{
185  KaonType = ParticleTable::getKaonType(-iso_d);
186  DeltaType = ParticleTable::getDeltaType(3*iso_d);
187  SigmaType = ParticleTable::getSigmaType(-2*iso_d);
188  }
189  }
190  else if(ParticleTable::getIsospin(particle1->getType()) == ParticleTable::getIsospin(particle2->getType())){// D+ p
191  if(rdm*11 < 2){
192  KaonType = ParticleTable::getKaonType(iso/2);
193  DeltaType = ParticleTable::getDeltaType(-iso/2);
194  SigmaType = ParticleTable::getSigmaType(iso);
195  }
196  else if(rdm*11 < 3){
197  KaonType = ParticleTable::getKaonType(iso/2);
198  DeltaType = ParticleTable::getDeltaType(iso/2);
199  SigmaType = SigmaZero;
200  }
201  else if(rdm*11 < 5){
202  KaonType = ParticleTable::getKaonType(iso/2);
203  DeltaType = ParticleTable::getDeltaType(3*iso/2);
204  SigmaType = ParticleTable::getSigmaType(-iso);
205  }
206  else if(rdm*11 < 7){
207  KaonType = ParticleTable::getKaonType(-iso/2);
208  DeltaType = ParticleTable::getDeltaType(iso/2);
209  SigmaType = ParticleTable::getSigmaType(iso);
210  }
211  else{
212  KaonType = ParticleTable::getKaonType(-iso/2);
213  DeltaType = ParticleTable::getDeltaType(3*iso/2);
214  SigmaType = SigmaZero;
215  }
216  }
217  else{// D++ n
218  if(rdm*15 < 2){
219  KaonType = ParticleTable::getKaonType(iso/2);
220  DeltaType = ParticleTable::getDeltaType(-iso/2);
221  SigmaType = ParticleTable::getSigmaType(iso);
222  }
223  else if(rdm*15 < 6){
224  KaonType = ParticleTable::getKaonType(iso/2);
225  DeltaType = ParticleTable::getDeltaType(iso/2);
226  SigmaType = SigmaZero;
227  }
228  else if(rdm*15 < 12){
229  KaonType = ParticleTable::getKaonType(iso/2);
230  DeltaType = ParticleTable::getDeltaType(3*iso/2);
231  SigmaType = ParticleTable::getSigmaType(-iso);
232  }
233  else if(rdm*15 < 14){
234  KaonType = ParticleTable::getKaonType(-iso/2);
235  DeltaType = ParticleTable::getDeltaType(iso/2);
236  SigmaType = ParticleTable::getSigmaType(iso);
237  }
238  else{
239  KaonType = ParticleTable::getKaonType(-iso/2);
240  DeltaType = ParticleTable::getDeltaType(3*iso/2);
241  SigmaType = SigmaZero;
242  }
243  }
244 
245 
246  particle1->setType(DeltaType);
247  particle1->setMass(sampleDeltaMass(sqrtS));
248  particle2->setType(SigmaType);
249 
250  ParticleList list;
251  list.push_back(particle1);
252  list.push_back(particle2);
253  const ThreeVector &rcol = particle2->getPosition();
254  const ThreeVector zero;
255  Particle *kaon = new Particle(KaonType,zero,rcol);
256  list.push_back(kaon);
257 
258  if(Random::shoot()<0.5) PhaseSpaceGenerator::generateBiased(sqrtS, list, 0, angularSlope);
259  else PhaseSpaceGenerator::generateBiased(sqrtS, list, 1, angularSlope);
260 
261  fs->addModifiedParticle(particle1);
262  fs->addModifiedParticle(particle2);
263  fs->addCreatedParticle(kaon);
264 
265  }
266 }