| Geant4 Cross Reference |
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 "G4INCLNDeltaEtaProductionChannel.hh"
39 #include "G4INCLKinematicsUtils.hh"
40 #include "G4INCLBinaryCollisionAvatar.hh"
41 #include "G4INCLRandom.hh"
42 #include "G4INCLGlobals.hh"
43 #include "G4INCLLogger.hh"
44 #include "G4INCLPhaseSpaceGenerator.hh"
45
46 namespace G4INCL {
47
48 const G4double NDeltaEtaProductionChannel::angularSlope = 6.;
49 const G4int NDeltaEtaProductionChannel::maxTries = 100000;
50
51 NDeltaEtaProductionChannel::NDeltaEtaProductionChannel(Particle *p1,Particle *p2)
52 : particle1(p1), particle2(p2)
53 {}
54
55 NDeltaEtaProductionChannel::~NDeltaEtaProductionChannel() {}
56
57 G4double NDeltaEtaProductionChannel::sampleDeltaMass(G4double ecmorigin) {
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');
77 return ParticleTable::minDeltaMass;
78 }
79
80 G4double rndm = ParticleTable::minDeltaMassRndm + Random::shoot() * deltaMassRndmRange;
81 y = std::tan(rndm);
82 x = ParticleTable::effectiveDeltaMass + 0.5*ParticleTable::effectiveDeltaWidth*y;
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
98 void NDeltaEtaProductionChannel::fillFinalState(FinalState *fs) {
99
100 /**
101 *
102 * Unlike NN -> NDelta, NN -> NDeltaEta is drawn from a phase-space generator
103 *
104 **/
105
106 G4int is1=ParticleTable::getIsospin(particle1->getType());
107 G4int is2=ParticleTable::getIsospin(particle2->getType());
108
109 ParticleList list;
110 list.push_back(particle1);
111 list.push_back(particle2);
112
113 // isospin Repartition of N and Delta;
114 G4double ecm = KinematicsUtils::totalEnergyInCM(particle1, particle2);
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
142 if(is1 == ParticleTable::getIsospin(DeltaMinus)) {
143 particle1->setType(DeltaMinus);
144 } else if(is1 == ParticleTable::getIsospin(DeltaZero)) {
145 particle1->setType(DeltaZero);
146 } else if(is1 == ParticleTable::getIsospin(DeltaPlus)) {
147 particle1->setType(DeltaPlus);
148 } else if(is1 == ParticleTable::getIsospin(DeltaPlusPlus)) {
149 particle1->setType(DeltaPlusPlus);
150 }
151
152 if(is2 == ParticleTable::getIsospin(Proton)) {
153 particle2->setType(Proton);
154 } else if(is2 == ParticleTable::getIsospin(Neutron)) {
155 particle2->setType(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
169 const G4double sqrtS = KinematicsUtils::totalEnergyInCM(particle1, particle2);
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);
174 G4double ctet=(particle1->getMomentum().dot(vz))/particle1->getMomentum().mag();
175 if (isospin == 0)
176 particle1->setHelicity(0.0);
177 else
178 particle1->setHelicity(ctet*ctet);
179
180 fs->addModifiedParticle(particle1);
181 fs->addModifiedParticle(particle2);
182
183 }
184
185 }
186