Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/coherent_elastic/src/G4LEnp.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/hadronic/models/coherent_elastic/src/G4LEnp.cc (Version 11.3.0) and /processes/hadronic/models/coherent_elastic/src/G4LEnp.cc (Version 3.1)


  1 //                                                  1 
  2 // *******************************************    
  3 // * License and Disclaimer                       
  4 // *                                              
  5 // * The  Geant4 software  is  copyright of th    
  6 // * the Geant4 Collaboration.  It is provided    
  7 // * conditions of the Geant4 Software License    
  8 // * LICENSE and available at  http://cern.ch/    
  9 // * include a list of copyright holders.         
 10 // *                                              
 11 // * Neither the authors of this software syst    
 12 // * institutes,nor the agencies providing fin    
 13 // * work  make  any representation or  warran    
 14 // * regarding  this  software system or assum    
 15 // * use.  Please see the license in the file     
 16 // * for the full disclaimer and the limitatio    
 17 // *                                              
 18 // * This  code  implementation is the result     
 19 // * technical work of the GEANT4 collaboratio    
 20 // * By using,  copying,  modifying or  distri    
 21 // * any work based  on the software)  you  ag    
 22 // * use  in  resulting  scientific  publicati    
 23 // * acceptance of all terms of the Geant4 Sof    
 24 // *******************************************    
 25 //                                                
 26                                                   
 27 // G4 Low energy model: n-p scattering            
 28 // F.W. Jones, L.G. Greeniaus, H.P. Wellisch      
 29                                                   
 30 // 11-OCT-2007 F.W. Jones: removed erroneous c    
 31 //             exchange of particles.             
 32 // FWJ 27-AUG-2010: extended to 5 GeV by Tony     
 33                                                   
 34 #include "G4LEnp.hh"                              
 35 #include "G4PhysicalConstants.hh"                 
 36 #include "G4SystemOfUnits.hh"                     
 37 #include "Randomize.hh"                           
 38 #include "G4ios.hh"                               
 39                                                   
 40 // Initialization of static data arrays:          
 41 #include "G4LEnpData.hh"                          
 42 #include "Randomize.hh"                           
 43                                                   
 44 #include "G4PhysicsModelCatalog.hh"               
 45                                                   
 46                                                   
 47 G4LEnp::G4LEnp():                                 
 48  G4HadronElastic("G4LEnp")  // G4HadronicInter    
 49 {                                                 
 50   secID = G4PhysicsModelCatalog::GetModelID( "    
 51   //    theParticleChange.SetNumberOfSecondari    
 52                                                   
 53   //    SetMinEnergy(10.*MeV);                    
 54   //    SetMaxEnergy(1200.*MeV);                  
 55   SetMinEnergy(0.);                               
 56   SetMaxEnergy(5.*GeV);                           
 57 }                                                 
 58                                                   
 59 G4LEnp::~G4LEnp()                                 
 60 {                                                 
 61       theParticleChange.Clear();                  
 62 }                                                 
 63                                                   
 64 G4HadFinalState*                                  
 65 G4LEnp::ApplyYourself(const G4HadProjectile& a    
 66 {                                                 
 67     theParticleChange.Clear();                    
 68     const G4HadProjectile* aParticle = &aTrack    
 69                                                   
 70     G4double P = aParticle->GetTotalMomentum()    
 71     G4double Px = aParticle->Get4Momentum().x(    
 72     G4double Py = aParticle->Get4Momentum().y(    
 73     G4double Pz = aParticle->Get4Momentum().z(    
 74     G4double ek = aParticle->GetKineticEnergy(    
 75     G4ThreeVector theInitial = aParticle->Get4    
 76                                                   
 77     if (verboseLevel > 1) {                       
 78       G4double E = aParticle->GetTotalEnergy()    
 79       G4double E0 = aParticle->GetDefinition()    
 80       G4double Q = aParticle->GetDefinition()-    
 81       G4int A = targetNucleus.GetA_asInt();       
 82       G4int Z = targetNucleus.GetZ_asInt();       
 83       G4cout << "G4LEnp:ApplyYourself: inciden    
 84              << aParticle->GetDefinition()->Ge    
 85       G4cout << "P = " << P/GeV << " GeV/c"       
 86              << ", Px = " << Px/GeV << " GeV/c    
 87              << ", Py = " << Py/GeV << " GeV/c    
 88              << ", Pz = " << Pz/GeV << " GeV/c    
 89       G4cout << "E = " << E/GeV << " GeV"         
 90              << ", kinetic energy = " << ek/Ge    
 91              << ", mass = " << E0/GeV << " GeV    
 92              << ", charge = " << Q << G4endl;     
 93       G4cout << "G4LEnp:ApplyYourself: materia    
 94       G4cout << "A = " << A                       
 95              << ", Z = " << Z                     
 96              << ", atomic mass "                  
 97              <<  G4Proton::Proton()->GetPDGMas    
 98              << G4endl;                           
 99       //                                          
100       // GHEISHA ADD operation to get total en    
101       //                                          
102       E += proton_mass_c2;                        
103       G4double E02 = E*E - P*P;                   
104       E0 = std::sqrt(std::abs(E02));              
105       if (E02 < 0)E0 *= -1;                       
106       Q += Z;                                     
107       G4cout << "G4LEnp:ApplyYourself: total:"    
108       G4cout << "E = " << E/GeV << " GeV"         
109              << ", mass = " << E0/GeV << " GeV    
110              << ", charge = " << Q << G4endl;     
111     }                                             
112                                                   
113     // Find energy bin                            
114                                                   
115     G4int je1 = 0;                                
116     G4int je2 = NENERGY - 1;                      
117     ek = ek/GeV;                                  
118     do {                                          
119       G4int midBin = (je1 + je2)/2;               
120       if (ek < elab[midBin])                      
121         je2 = midBin;                             
122       else                                        
123         je1 = midBin;                             
124     } while (je2 - je1 > 1);  /* Loop checking    
125     G4double delab = elab[je2] - elab[je1];       
126                                                   
127     // Sample the angle                           
128                                                   
129     G4double sample = G4UniformRand();            
130     G4int ke1 = 0;                                
131     G4int ke2 = NANGLE - 1;                       
132     G4double dsig = sig[je2][0] - sig[je1][0];    
133     G4double rc = dsig/delab;                     
134     G4double b = sig[je1][0] - rc*elab[je1];      
135     G4double sigint1 = rc*ek + b;                 
136     G4double sigint2 = 0.;                        
137                                                   
138     if (verboseLevel > 1) {                       
139       G4cout << "sample=" << sample << G4endl     
140        << ke1 << " " << ke2 << " "                
141        << sigint1 << " " << sigint2 << G4endl;    
142     }                                             
143     do {                                          
144       G4int midBin = (ke1 + ke2)/2;               
145       dsig = sig[je2][midBin] - sig[je1][midBi    
146       rc = dsig/delab;                            
147       b = sig[je1][midBin] - rc*elab[je1];        
148       G4double sigint = rc*ek + b;                
149       if (sample < sigint) {                      
150         ke2 = midBin;                             
151         sigint2 = sigint;                         
152       }                                           
153       else {                                      
154         ke1 = midBin;                             
155         sigint1 = sigint;                         
156       }                                           
157       if (verboseLevel > 1) {                     
158   G4cout << ke1 << " " << ke2 << " "              
159          << sigint1 << " " << sigint2 << G4end    
160       }                                           
161     } while (ke2 - ke1 > 1);  /* Loop checking    
162                                                   
163     dsig = sigint2 - sigint1;                     
164     rc = 1./dsig;                                 
165     b = ke1 - rc*sigint1;                         
166     G4double kint = rc*sample + b;                
167     G4double theta = (0.5 + kint)*pi/180.;        
168                                                   
169     if (verboseLevel > 1) {                       
170       G4cout << "   energy bin " << je1 << " e    
171       G4cout << "   angle bin " << kint << " a    
172     }                                             
173                                                   
174     // Get the target particle                    
175                                                   
176     G4DynamicParticle* targetParticle = target    
177                                                   
178     G4double E1 = aParticle->GetTotalEnergy();    
179     G4double M1 = aParticle->GetDefinition()->    
180     G4double E2 = targetParticle->GetTotalEner    
181     G4double M2 = targetParticle->GetDefinitio    
182     G4double totalEnergy = E1 + E2;               
183     G4double pseudoMass = std::sqrt(totalEnerg    
184                                                   
185     // Transform into centre of mass system       
186                                                   
187     G4double px = (M2/pseudoMass)*Px;             
188     G4double py = (M2/pseudoMass)*Py;             
189     G4double pz = (M2/pseudoMass)*Pz;             
190     G4double p = std::sqrt(px*px + py*py + pz*    
191                                                   
192     if (verboseLevel > 1) {                       
193       G4cout << "  E1, M1 (GeV) " << E1/GeV <<    
194       G4cout << "  E2, M2 (GeV) " << E2/GeV <<    
195       G4cout << "  particle  1 momentum in CM     
196            << pz/GeV << " " << p/GeV << G4endl    
197     }                                             
198                                                   
199     // First scatter w.r.t. Z axis                
200     G4double phi = G4UniformRand()*twopi;         
201     G4double pxnew = p*std::sin(theta)*std::co    
202     G4double pynew = p*std::sin(theta)*std::si    
203     G4double pznew = p*std::cos(theta);           
204                                                   
205     // Rotate according to the direction of th    
206     if (px*px + py*py > 0) {                      
207       G4double cost, sint, ph, cosp, sinp;        
208       cost = pz/p;                                
209       sint = (std::sqrt(std::fabs((1-cost)*(1+    
210       py < 0 ? ph = 3*halfpi : ph = halfpi;       
211       if (std::abs(px) > 0.000001*GeV) ph = st    
212       cosp = std::cos(ph);                        
213       sinp = std::sin(ph);                        
214       px = (cost*cosp*pxnew - sinp*pynew + sin    
215       py = (cost*sinp*pxnew + cosp*pynew + sin    
216       pz = (-sint*pxnew                  + cos    
217     }                                             
218     else {                                        
219       px = pxnew;                                 
220       py = pynew;                                 
221       pz = pznew;                                 
222     }                                             
223                                                   
224     if (verboseLevel > 1) {                       
225       G4cout << "  AFTER SCATTER..." << G4endl    
226       G4cout << "  particle 1 momentum in CM "    
227            << pz/GeV << " " << p/GeV << G4endl    
228     }                                             
229                                                   
230     // Transform to lab system                    
231                                                   
232     G4double E1pM2 = E1 + M2;                     
233     G4double betaCM  = P/E1pM2;                   
234     G4double betaCMx = Px/E1pM2;                  
235     G4double betaCMy = Py/E1pM2;                  
236     G4double betaCMz = Pz/E1pM2;                  
237     G4double gammaCM = E1pM2/std::sqrt(E1pM2*E    
238                                                   
239     if (verboseLevel > 1) {                       
240       G4cout << "  betaCM " << betaCMx << " "     
241              << betaCMz << " " << betaCM << G4    
242       G4cout << "  gammaCM " << gammaCM << G4e    
243     }                                             
244                                                   
245     // Now following GLOREN...                    
246                                                   
247     G4double BETA[5], PA[5], PB[5];               
248     BETA[1] = -betaCMx;                           
249     BETA[2] = -betaCMy;                           
250     BETA[3] = -betaCMz;                           
251     BETA[4] = gammaCM;                            
252                                                   
253     //The incident particle...                    
254                                                   
255     PA[1] = px;                                   
256     PA[2] = py;                                   
257     PA[3] = pz;                                   
258     PA[4] = std::sqrt(M1*M1 + p*p);               
259                                                   
260     G4double BETPA  = BETA[1]*PA[1] + BETA[2]*    
261     G4double BPGAM  = (BETPA * BETA[4]/(BETA[4    
262                                                   
263     PB[1] = PA[1] + BPGAM  * BETA[1];             
264     PB[2] = PA[2] + BPGAM  * BETA[2];             
265     PB[3] = PA[3] + BPGAM  * BETA[3];             
266     PB[4] = (PA[4] - BETPA) * BETA[4];            
267                                                   
268     G4DynamicParticle* newP = new G4DynamicPar    
269     newP->SetDefinition(aParticle->GetDefiniti    
270     newP->SetMomentum(G4ThreeVector(PB[1], PB[    
271                                                   
272     //The target particle...                      
273                                                   
274     PA[1] = -px;                                  
275     PA[2] = -py;                                  
276     PA[3] = -pz;                                  
277     PA[4] = std::sqrt(M2*M2 + p*p);               
278                                                   
279     BETPA  = BETA[1]*PA[1] + BETA[2]*PA[2] + B    
280     BPGAM  = (BETPA * BETA[4]/(BETA[4] + 1.) -    
281                                                   
282     PB[1] = PA[1] + BPGAM  * BETA[1];             
283     PB[2] = PA[2] + BPGAM  * BETA[2];             
284     PB[3] = PA[3] + BPGAM  * BETA[3];             
285     PB[4] = (PA[4] - BETPA) * BETA[4];            
286                                                   
287     targetParticle->SetMomentum(G4ThreeVector(    
288                                                   
289     if (verboseLevel > 1) {                       
290       G4cout << "  particle 1 momentum in LAB     
291            << newP->GetMomentum()*(1./GeV)        
292            << " " << newP->GetTotalMomentum()/    
293       G4cout << "  particle 2 momentum in LAB     
294            << targetParticle->GetMomentum()*(1    
295            << " " << targetParticle->GetTotalM    
296       G4cout << "  TOTAL momentum in LAB "        
297            << (newP->GetMomentum()+targetParti    
298            << " "                                 
299            << (newP->GetMomentum()+targetParti    
300            << G4endl;                             
301     }                                             
302                                                   
303     theParticleChange.SetMomentumChange(newP->    
304     theParticleChange.SetEnergyChange(newP->Ge    
305     delete newP;                                  
306     theParticleChange.AddSecondary(targetParti    
307                                                   
308     return &theParticleChange;                    
309 }                                                 
310                                                   
311 //////////////////////////////////////////////    
312 //                                                
313 // sample momentum transfer using Lab. momentu    
314                                                   
315 G4double G4LEnp::SampleInvariantT(const G4Part    
316           G4double plab, G4int , G4int )          
317 {                                                 
318   G4double nMass = p->GetPDGMass(); // 939.565    
319   G4double ek = std::sqrt(plab*plab+nMass*nMas    
320                                                   
321     // Find energy bin                            
322                                                   
323   G4int je1 = 0;                                  
324   G4int je2 = NENERGY - 1;                        
325   ek = ek/GeV;                                    
326                                                   
327   do                                              
328   {                                               
329       G4int midBin = (je1 + je2)/2;               
330       if (ek < elab[midBin])                      
331         je2 = midBin;                             
332       else                                        
333         je1 = midBin;                             
334   } while (je2 - je1 > 1);  /* Loop checking,     
335                                                   
336   G4double delab = elab[je2] - elab[je1];         
337                                                   
338     // Sample the angle                           
339                                                   
340   G4double sample = G4UniformRand();              
341   G4int ke1 = 0;                                  
342   G4int ke2 = NANGLE - 1;                         
343   G4double dsig = sig[je2][0] - sig[je1][0];      
344   G4double rc = dsig/delab;                       
345   G4double b = sig[je1][0] - rc*elab[je1];        
346   G4double sigint1 = rc*ek + b;                   
347   G4double sigint2 = 0.;                          
348                                                   
349   do                                              
350   {                                               
351       G4int midBin = (ke1 + ke2)/2;               
352       dsig = sig[je2][midBin] - sig[je1][midBi    
353       rc = dsig/delab;                            
354       b = sig[je1][midBin] - rc*elab[je1];        
355       G4double sigint = rc*ek + b;                
356                                                   
357       if (sample < sigint)                        
358       {                                           
359         ke2 = midBin;                             
360         sigint2 = sigint;                         
361       }                                           
362       else                                        
363       {                                           
364         ke1 = midBin;                             
365         sigint1 = sigint;                         
366       }                                           
367   } while (ke2 - ke1 > 1);  /* Loop checking,     
368                                                   
369   dsig = sigint2 - sigint1;                       
370   rc = 1./dsig;                                   
371   b = ke1 - rc*sigint1;                           
372                                                   
373   G4double kint = rc*sample + b;                  
374   G4double theta = (0.5 + kint)*pi/180.;          
375   G4double t = 0.5*plab*plab*(1-std::cos(theta    
376                                                   
377   return t;                                       
378 }                                                 
379                                                   
380  // end of file                                   
381