ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLPionResonanceDecayChannel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLPionResonanceDecayChannel.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"
44 // #include <cassert>
45 
46 
47 namespace G4INCL {
48 
49  const G4double PionResonanceDecayChannel::angularSlope = 0.; // Slope to be defined, if needed.
50 
52  :theParticle(p), incidentDirection(dir)
53  { }
54 
56 
57 // Decay during the intranuclear cascade for Omega only
59  const G4double m = p->getMass();
60  const G4double geff = p->getEnergy()/m;
61 // const G4double geta = 1.31e-3;
62  const G4double gomega = 8.49;
63  G4double gg=0.;
64  switch (p->getType()) {
65 /* case Eta:
66  gg=geta;
67  break;*/
68  case Omega:
69  gg=gomega;
70  break;
71  default:
72  INCL_FATAL("Unrecognized pion resonance type; type=" << p->getType() << '\n');
73  break;
74  }
75  const G4double tpires = -G4INCL::PhysicalConstants::hc/gg*std::log(Random::shoot())*geff;
76  return tpires;
77  }
78 
79  void PionResonanceDecayChannel::sampleAngles(G4double *ctet_par, G4double *stet_par, G4double *phi_par) {
80 
81  (*ctet_par) = -1.0 + 2.0*Random::shoot();
82  if(std::abs(*ctet_par) > 1.0) (*ctet_par) = Math::sign(*ctet_par);
83  (*stet_par) = std::sqrt(1.-(*ctet_par)*(*ctet_par));
84  (*phi_par) = Math::twoPi * Random::shoot();
85  }
86 
88 
89  ParticleType createdType;
90  ParticleType pionType1=Neutron; // to avoid forgetting pionType definition when 3 particles are emitted
91  ParticleType pionType2=Neutron;
92 
93  const G4double sqrtS = theParticle->getMass();
94  G4int nbpart = 3; // number of emitted particles
95  G4double drnd=Random::shoot();
96  switch (theParticle->getType()) {
97  case Eta:
98  if (drnd < 0.3972) { // renormalized to the only four decays taken into account here
99 // 2 photons
100  nbpart=2;
102  createdType = Photon;
103  }
104  else if (drnd < 0.7265) {
105 // 3 pi0
107  pionType1 = PiZero;
108  pionType2 = PiZero;
109  }
110  else if (drnd < 0.9575) {
111 // pi+ pi- pi0
113  pionType1 = PiPlus;
114  pionType2 = PiMinus;
115  }
116  else {
117 // pi+ pi- photon
119  pionType1 = PiPlus;
120  pionType2 = PiMinus;
121  }
122  break;
123  case Omega:
124  if (drnd < 0.9009) { // renormalized to the only three decays taken into account here
125 // pi+ pi- pi0
127  pionType1 = PiPlus;
128  pionType2 = PiMinus;
129  }
130  else if (drnd < 0.9845) {
131 // pi0 photon
132  nbpart=2;
134  createdType = Photon;
135  }
136  else {
137 // pi+ pi-
138  nbpart=2;
140  createdType = PiMinus;
141  }
142  break;
143  default:
144  INCL_FATAL("Unrecognized pion resonance type; type=" << theParticle->getType() << '\n');
145  break;
146  }
147 
148  if (nbpart == 2) {
149  G4double fi, ctet, stet;
150  sampleAngles(&ctet, &stet, &fi);
151 
152  G4double cfi = std::cos(fi);
153  G4double sfi = std::sin(fi);
154  G4double beta = incidentDirection.mag();
155 
156  G4double q1, q2, q3;
157  G4double sal=0.0;
158  if (beta >= 1.0e-10)
159  sal = incidentDirection.perp()/beta;
160  if (sal >= 1.0e-6) {
164  G4double cal = b3/beta;
165  G4double t1 = ctet+cal*stet*sfi/sal;
166  G4double t2 = stet/sal;
167  q1=(b1*t1+b2*t2*cfi)/beta;
168  q2=(b2*t1-b1*t2*cfi)/beta;
169  q3=(b3*t1/beta-t2*sfi);
170  } else {
171  q1 = stet*cfi;
172  q2 = stet*sfi;
173  q3 = ctet;
174  }
175 
177  theParticle->getMass(),
178  ParticleTable::getINCLMass(createdType));
179  q1 *= xq;
180  q2 *= xq;
181  q3 *= xq;
182 
183  ThreeVector createdMomentum(q1, q2, q3);
184  ThreeVector createdPosition(theParticle->getPosition());
185  Particle *createdParticle = new Particle(createdType, createdMomentum, createdPosition);
186  theParticle->setMomentum(-createdMomentum);
188 
190  fs->addCreatedParticle(createdParticle);
191 
192  }
193  else if (nbpart == 3) {
194 // assert(pionType1!=Neutron && pionType2!=Neutron);
195  ParticleList list;
196  list.push_back(theParticle);
197  const ThreeVector &rposdecay = theParticle->getPosition();
198  const ThreeVector zero;
199  Particle *Pion1 = new Particle(pionType1,zero,rposdecay);
200  Particle *Pion2 = new Particle(pionType2,zero,rposdecay);
201  list.push_back(Pion1);
202  list.push_back(Pion2);
203 
205  fs->addCreatedParticle(Pion1);
206  fs->addCreatedParticle(Pion2);
207 
208  // PhaseSpaceGenerator::generateBiased(sqrtS, list, 0, angularSlope); Biasing?
209  PhaseSpaceGenerator::generate(sqrtS, list);
210  }
211 
212  }
213 }