| 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 "G4INCLNDeltaOmegaProductionChannel.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 NDeltaOmegaProductionChannel::angularSlope = 6.;
49 const G4int NDeltaOmegaProductionChannel::maxTries = 100000;
50
51 NDeltaOmegaProductionChannel::NDeltaOmegaProductionChannel(Particle *p1,Particle *p2)
52 : particle1(p1), particle2(p2)
53 {}
54
55 NDeltaOmegaProductionChannel::~NDeltaOmegaProductionChannel() {}
56
57 G4double NDeltaOmegaProductionChannel::sampleDeltaMass(G4double ecmorigin) {
58 const G4double ecm = ecmorigin - 783.437; // 783.437 MeV translation to open pion(delta) production in NNOmega
59 const G4double maxDeltaMass = ecm - ParticleTable::effectiveNucleonMass - 1.0;
60 const G4double maxDeltaMassRndm = std::atan((maxDeltaMass-ParticleTable::effectiveDeltaMass)*2./ParticleTable::effectiveDeltaWidth);
61 const G4double deltaMassRndmRange = maxDeltaMassRndm - ParticleTable::minDeltaMassRndm;
62 // assert(deltaMassRndmRange>0.);
63
64 G4double y=ecm*ecm;
65 G4double q2=(y-1.157776E6)*(y-6.4E5)/y/4.0; // 1.157776E6 = 1076^2, 6.4E5 = 800^2
66 G4double q3=std::pow(std::sqrt(q2), 3.);
67 const G4double f3max=q3/(q3+5.832E6); // 5.832E6 = 180^3
68 G4double x;
69
70 G4int nTries = 0;
71 G4bool success = false;
72 while(!success) { /* Loop checking, 10.07.2015, D.Mancusi */
73 if(++nTries >= maxTries) {
74 INCL_WARN("NDeltaOmegaProductionChannel::sampleDeltaMass loop was stopped because maximum number of tries was reached. Minimum delta mass "
75 << ParticleTable::minDeltaMass << " MeV with CM energy " << ecm << " MeV may be unphysical." << '\n');
76 return ParticleTable::minDeltaMass;
77 }
78
79 G4double rndm = ParticleTable::minDeltaMassRndm + Random::shoot() * deltaMassRndmRange;
80 y = std::tan(rndm);
81 x = ParticleTable::effectiveDeltaMass + 0.5*ParticleTable::effectiveDeltaWidth*y;
82 // assert(x>=ParticleTable::minDeltaMass && ecm >= x + ParticleTable::effectiveNucleonMass + 1.0);
83
84 // generation of the delta mass with the penetration factor
85 // (see prc56(1997)2431)
86 y=x*x;
87 q2=(y-1.157776E6)*(y-6.4E5)/y/4.0; // 1.157776E6 = 1076^2, 6.4E5 = 800^2
88 q3=std::pow(std::sqrt(q2), 3.);
89 const G4double f3=q3/(q3+5.832E6); // 5.832E6 = 180^3
90 rndm = Random::shoot();
91 if (rndm*f3max < f3)
92 success = true;
93 }
94 return x;
95 }
96
97 void NDeltaOmegaProductionChannel::fillFinalState(FinalState *fs) {
98
99 /**
100 *
101 * Unlike NN -> NDelta, NN -> NDeltaOmega is drawn from a phase-space generator
102 *
103 **/
104
105 G4int is1=ParticleTable::getIsospin(particle1->getType());
106 G4int is2=ParticleTable::getIsospin(particle2->getType());
107
108 ParticleList list;
109 list.push_back(particle1);
110 list.push_back(particle2);
111
112 // isospin Repartition of N and Delta;
113 G4double ecm = KinematicsUtils::totalEnergyInCM(particle1, particle2);
114 const G4int isospin = is1+is2;
115
116 G4double rndm = 0.0;
117 G4double xmdel = sampleDeltaMass(ecm);
118
119 G4int index2=0;
120 if (isospin == 0) { // pn case
121 rndm = Random::shoot();
122 if (rndm < 0.5) index2=1;
123 }
124
125 if (isospin == 0) {
126 if(index2 == 1) {
127 G4int isi=is1;
128 is1=is2;
129 is2=isi;
130 }
131 // particle1->setHelicity(0.0);
132 } else {
133 rndm = Random::shoot();
134 if (rndm >= 0.25) {
135 is1=3*is1;
136 is2=-is2;
137 }
138 // particle1->setHelicity(ctet*ctet);
139 }
140
141 if(is1 == ParticleTable::getIsospin(DeltaMinus)) {
142 particle1->setType(DeltaMinus);
143 } else if(is1 == ParticleTable::getIsospin(DeltaZero)) {
144 particle1->setType(DeltaZero);
145 } else if(is1 == ParticleTable::getIsospin(DeltaPlus)) {
146 particle1->setType(DeltaPlus);
147 } else if(is1 == ParticleTable::getIsospin(DeltaPlusPlus)) {
148 particle1->setType(DeltaPlusPlus);
149 }
150
151 if(is2 == ParticleTable::getIsospin(Proton)) {
152 particle2->setType(Proton);
153 } else if(is2 == ParticleTable::getIsospin(Neutron)) {
154 particle2->setType(Neutron);
155 }
156
157 if(particle1->isDelta()) particle1->setMass(xmdel);
158 if(particle2->isDelta()) particle2->setMass(xmdel);
159
160 const ThreeVector &rcolnucleon1 = particle1->getPosition();
161 const ThreeVector &rcolnucleon2 = particle2->getPosition();
162 const ThreeVector rcol = (rcolnucleon1+rcolnucleon2)*0.5;
163 const ThreeVector zero;
164 Particle *omega = new Particle(Omega,zero,rcol);
165 list.push_back(omega);
166 fs->addCreatedParticle(omega);
167
168 const G4double sqrtS = KinematicsUtils::totalEnergyInCM(particle1, particle2);
169 G4int biasIndex = ((Random::shoot()<0.5) ? 0 : 1);
170 PhaseSpaceGenerator::generateBiased(sqrtS, list, biasIndex, angularSlope);
171
172 const ThreeVector vz(0.0,0.0,1.0);
173 G4double ctet=(particle1->getMomentum().dot(vz))/particle1->getMomentum().mag();
174 if (isospin == 0)
175 particle1->setHelicity(0.0);
176 else
177 particle1->setHelicity(ctet*ctet);
178
179 fs->addModifiedParticle(particle1);
180 fs->addModifiedParticle(particle2);
181
182 }
183
184 }
185