ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLNDeltaEtaProductionChannel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLNDeltaEtaProductionChannel.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"
45 
46 namespace G4INCL {
47 
50 
52  : particle1(p1), particle2(p2)
53  {}
54 
56 
58 // const G4double ecm = ecmorigin - 686.987; // 686.987 MeV translation to open pion(delta) production in NNEta
59  const G4double ecm = ecmorigin - 581.437; // 581.437 MeV translation to open pion(delta) production in NNEta
60  const G4double maxDeltaMass = ecm - ParticleTable::effectiveNucleonMass - 1.0;
61  const G4double maxDeltaMassRndm = std::atan((maxDeltaMass-ParticleTable::effectiveDeltaMass)*2./ParticleTable::effectiveDeltaWidth);
62  const G4double deltaMassRndmRange = maxDeltaMassRndm - ParticleTable::minDeltaMassRndm;
63 // assert(deltaMassRndmRange>0.);
64 
65  G4double y=ecm*ecm;
66  G4double q2=(y-1.157776E6)*(y-6.4E5)/y/4.0; // 1.157776E6 = 1076^2, 6.4E5 = 800^2
67  G4double q3=std::pow(std::sqrt(q2), 3.);
68  const G4double f3max=q3/(q3+5.832E6); // 5.832E6 = 180^3
69  G4double x;
70 
71  G4int nTries = 0;
72  G4bool success = false;
73  while(!success) { /* Loop checking, 10.07.2015, D.Mancusi */
74  if(++nTries >= maxTries) {
75  INCL_WARN("NDeltaEtaProductionChannel::sampleDeltaMass loop was stopped because maximum number of tries was reached. Minimum delta mass "
76  << ParticleTable::minDeltaMass << " MeV with CM energy " << ecm << " MeV may be unphysical." << '\n');
78  }
79 
80  G4double rndm = ParticleTable::minDeltaMassRndm + Random::shoot() * deltaMassRndmRange;
81  y = std::tan(rndm);
83 // assert(x>=ParticleTable::minDeltaMass && ecm >= x + ParticleTable::effectiveNucleonMass + 1.0);
84 
85  // generation of the delta mass with the penetration factor
86  // (see prc56(1997)2431)
87  y=x*x;
88  q2=(y-1.157776E6)*(y-6.4E5)/y/4.0; // 1.157776E6 = 1076^2, 6.4E5 = 800^2
89  q3=std::pow(std::sqrt(q2), 3.);
90  const G4double f3=q3/(q3+5.832E6); // 5.832E6 = 180^3
91  rndm = Random::shoot();
92  if (rndm*f3max < f3)
93  success = true;
94  }
95  return x;
96  }
97 
99 
108 
109  ParticleList list;
110  list.push_back(particle1);
111  list.push_back(particle2);
112 
113 // isospin Repartition of N and Delta;
115  const G4int isospin = is1+is2;
116 
117  G4double rndm = 0.0;
118  G4double xmdel = sampleDeltaMass(ecm);
119 
120  G4int index2=0;
121  if (isospin == 0) { // pn case
122  rndm = Random::shoot();
123  if (rndm < 0.5) index2=1;
124  }
125 
126  if (isospin == 0) {
127  if(index2 == 1) {
128  G4int isi=is1;
129  is1=is2;
130  is2=isi;
131  }
132 // particle1->setHelicity(0.0);
133  } else {
134  rndm = Random::shoot();
135  if (rndm >= 0.25) {
136  is1=3*is1;
137  is2=-is2;
138  }
139 // particle1->setHelicity(ctet*ctet);
140  }
141 
144  } else if(is1 == ParticleTable::getIsospin(DeltaZero)) {
146  } else if(is1 == ParticleTable::getIsospin(DeltaPlus)) {
148  } else if(is1 == ParticleTable::getIsospin(DeltaPlusPlus)) {
150  }
151 
152  if(is2 == ParticleTable::getIsospin(Proton)) {
154  } else if(is2 == ParticleTable::getIsospin(Neutron)) {
156  }
157 
158  if(particle1->isDelta()) particle1->setMass(xmdel);
159  if(particle2->isDelta()) particle2->setMass(xmdel);
160 
161  const ThreeVector &rcolnucleon1 = particle1->getPosition();
162  const ThreeVector &rcolnucleon2 = particle2->getPosition();
163  const ThreeVector rcol = (rcolnucleon1+rcolnucleon2)*0.5;
164  const ThreeVector zero;
165  Particle *eta = new Particle(Eta,zero,rcol);
166  list.push_back(eta);
167  fs->addCreatedParticle(eta);
168 
170  G4int biasIndex = ((Random::shoot()<0.5) ? 0 : 1);
171  PhaseSpaceGenerator::generateBiased(sqrtS, list, biasIndex, angularSlope);
172 
173  const ThreeVector vz(0.0,0.0,1.0);
175  if (isospin == 0)
176  particle1->setHelicity(0.0);
177  else
178  particle1->setHelicity(ctet*ctet);
179 
182 
183  }
184 
185 }