ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLDeltaDecayChannel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLDeltaDecayChannel.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 "G4INCLDeltaDecayChannel.hh"
39 #include "G4INCLKinematicsUtils.hh"
40 #include "G4INCLBinaryCollisionAvatar.hh"
41 #include "G4INCLRandom.hh"
42 #include "G4INCLGlobals.hh"
43 
44 namespace G4INCL {
45 
46  DeltaDecayChannel::DeltaDecayChannel(Particle *p, ThreeVector const &dir)
47  :theParticle(p), incidentDirection(dir)
48  { }
49 
50  DeltaDecayChannel::~DeltaDecayChannel() {}
51 
52  G4double DeltaDecayChannel::computeDecayTime(Particle *p) {
53  const G4double m = p->getMass();
54  const G4double g0 = 115.0;
55  G4double gg = g0;
56  if(m > 1500.0) gg = 200.0;
57  const G4double geff = p->getEnergy()/m;
58  const G4double qqq = KinematicsUtils::momentumInCM(m, ParticleTable::effectiveNucleonMass, ParticleTable::effectivePionMass);
59  const G4double psf = std::pow(qqq, 3)/(std::pow(qqq, 3) + 5832000.0); // phase space factor 5.832E6 = 180^3
60  const G4double tdel = -G4INCL::PhysicalConstants::hc/(gg*psf)*std::log(Random::shoot())*geff; // fm
61  if( m > 1400) return tdel * 1./(1. + std::pow((m-1400)/g0,2)); // reduction of Delta life time for high masses.
62  return tdel; // fm
63  }
64 
65  void DeltaDecayChannel::sampleAngles(G4double *ctet_par, G4double *stet_par, G4double *phi_par) {
66  const G4double hel = theParticle->getHelicity();
67  unsigned long loopCounter = 0;
68  const unsigned long maxLoopCounter = 10000000;
69  do {
70  (*ctet_par) = -1.0 + 2.0*Random::shoot();
71  if(std::abs(*ctet_par) > 1.0) (*ctet_par) = Math::sign(*ctet_par);
72  ++loopCounter;
73  } while(loopCounter<maxLoopCounter && Random::shoot() > ((1.0 + 3.0 * hel * (*ctet_par) * (*ctet_par)) /(1.0 + 3.0 * hel))); /* Loop checking, 10.07.2015, D.Mancusi */
74  (*stet_par) = std::sqrt(1.-(*ctet_par)*(*ctet_par));
75  (*phi_par) = Math::twoPi * Random::shoot();
76  }
77 
78  void DeltaDecayChannel::fillFinalState(FinalState *fs) {
79  // SUBROUTINE DECAY2(P1,P2,P3,WP,ij,
80  // s X1,X2,hel,B1,B2,B3)
81 
82  // This routine describes the anisotropic decay of a particle of mass
83  // xi into 2 particles of masses x1,x2.
84  // The anisotropy is supposed to follow a 1+3*hel*(cos(theta))**2
85  // law with respect to the direction of the incoming particle.
86  // In the input, p1,p2,p3 is the momentum of particle xi.
87  // In the output, p1,p2,p3 is the momentum of particle x1 , while
88  // q1,q2,q3 is the momentum of particle x2.
89 
90  // COMMON/bl12/QQ1(200),QQ2(200),QQ3(200),QQ4(200),
91  // s YY1(200),YY2(200),YY3(200),YM(200),IPI(200)
92  // common/hazard/ial,IY1,IY2,IY3,IY4,IY5,IY6,IY7,IY8,IY9,IY10,
93  // s IY11,IY12,IY13,IY14,IY15,IY16,IY17,IY18,IY19
94 
95  // DATA IY8,IY9,IY10/82345,92345,45681/
96  // PCM(E,A,C)=0.5*SQRT((E**2-(A+C)**2)*(E**2-(A-C)**2))/E P-N20800
97  // XI=YM(ij)
98 
99  // XE=WP P-N20810
100  // B1=P1/XE P-N20820
101  // B2=P2/XE P-N20830
102  // B3=P3/XE
103  // XQ=PCM(XI,X1,X2)
104 
105  const G4double deltaMass = theParticle->getMass();
106 
107  G4double fi, ctet, stet;
108  sampleAngles(&ctet, &stet, &fi);
109 
110  G4double cfi = std::cos(fi);
111  G4double sfi = std::sin(fi);
112  G4double beta = incidentDirection.mag();
113 
114  G4double q1, q2, q3;
115  G4double sal=0.0;
116  if (beta >= 1.0e-10)
117  sal = incidentDirection.perp()/beta;
118  if (sal >= 1.0e-6) {
119  G4double b1 = incidentDirection.getX();
120  G4double b2 = incidentDirection.getY();
121  G4double b3 = incidentDirection.getZ();
122  G4double cal = b3/beta;
123  G4double t1 = ctet+cal*stet*sfi/sal;
124  G4double t2 = stet/sal;
125  q1=(b1*t1+b2*t2*cfi)/beta;
126  q2=(b2*t1-b1*t2*cfi)/beta;
127  q3=(b3*t1/beta-t2*sfi);
128  } else {
129  q1 = stet*cfi;
130  q2 = stet*sfi;
131  q3 = ctet;
132  }
133  theParticle->setHelicity(0.0);
134 
135  ParticleType pionType;
136  switch(theParticle->getType()) {
137  case DeltaPlusPlus:
138  theParticle->setType(Proton);
139  pionType = PiPlus;
140  break;
141  case DeltaPlus:
142  if(Random::shoot() < 1.0/3.0) {
143  theParticle->setType(Neutron);
144  pionType = PiPlus;
145  } else {
146  theParticle->setType(Proton);
147  pionType = PiZero;
148  }
149  break;
150  case DeltaZero:
151  if(Random::shoot() < 1.0/3.0) {
152  theParticle->setType(Proton);
153  pionType = PiMinus;
154  } else {
155  theParticle->setType(Neutron);
156  pionType = PiZero;
157  }
158  break;
159  case DeltaMinus:
160  theParticle->setType(Neutron);
161  pionType = PiMinus;
162  break;
163  default:
164  INCL_FATAL("Unrecognized delta type; type=" << theParticle->getType() << '\n');
165  pionType = UnknownParticle;
166  break;
167  }
168 
169  G4double xq = KinematicsUtils::momentumInCM(deltaMass,
170  theParticle->getMass(),
171  ParticleTable::getINCLMass(pionType));
172 
173  q1 *= xq;
174  q2 *= xq;
175  q3 *= xq;
176 
177  ThreeVector pionMomentum(q1, q2, q3);
178  ThreeVector pionPosition(theParticle->getPosition());
179  Particle *pion = new Particle(pionType, pionMomentum, pionPosition);
180  theParticle->setMomentum(-pionMomentum);
181  theParticle->adjustEnergyFromMomentum();
182 
183  fs->addModifiedParticle(theParticle);
184  fs->addCreatedParticle(pion);
185  // call loren(q1,q2,q3,b1,b2,b3,wq)
186  // call loren(p1,p2,p3,b1,b2,b3,wp)
187  // qq1(ij)=q1
188  // qq2(ij)=q2
189  // qq3(ij)=q3
190  // qq4(ij)=wq
191  // ym(ij)=xi
192  // RETURN P-N21120
193  // END P-N21130
194  }
195 }