Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/binary_cascade/src/G4Absorber.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/binary_cascade/src/G4Absorber.cc (Version 11.3.0) and /processes/hadronic/models/binary_cascade/src/G4Absorber.cc (Version 5.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 #include "G4Absorber.hh"                          
 27 #include "G4KineticTrack.hh"                      
 28 #include "G4PionPlus.hh"                          
 29 #include "G4PionMinus.hh"                         
 30 #include "G4PionZero.hh"                          
 31 #include "G4Proton.hh"                            
 32 #include "G4Neutron.hh"                           
 33                                                   
 34 #include "G4PhysicalConstants.hh"                 
 35 #include "G4SystemOfUnits.hh"                     
 36 #include "G4LorentzRotation.hh"                   
 37                                                   
 38 G4Absorber::G4Absorber(G4double cutOnP)           
 39 {                                                 
 40   theCutOnP = cutOnP;                             
 41   theAbsorbers = new G4KineticTrackVector;        
 42   theProducts = new G4KineticTrackVector;         
 43 }                                                 
 44                                                   
 45                                                   
 46 G4Absorber::~G4Absorber()                         
 47 {                                                 
 48   delete theAbsorbers;                            
 49   delete theProducts;                             
 50 }                                                 
 51                                                   
 52                                                   
 53 bool G4Absorber::WillBeAbsorbed(const G4Kineti    
 54 {                                                 
 55  // FixMe: actually only for pions                
 56 //  if(kt.Get4Momentum().vect().mag() < theCut    
 57 // Cut on kinetic Energy...                       
 58   if (kt.Get4Momentum().e() - kt.GetActualMass    
 59   {                                               
 60       if(kt.GetDefinition() == G4PionPlus::Pio    
 61    kt.GetDefinition() == G4PionZero::PionZero(    
 62    kt.GetDefinition() == G4PionMinus::PionMinu    
 63       {                                           
 64    return true;                                   
 65       }                                           
 66   }                                               
 67   return false;                                   
 68 }                                                 
 69                                                   
 70                                                   
 71                                                   
 72 G4bool G4Absorber::Absorb(G4KineticTrack & kt,    
 73 {                                                 
 74   if(!FindAbsorbers(kt, tgt))                     
 75     return false;                                 
 76   return FindProducts(kt);                        
 77 }                                                 
 78                                                   
 79                                                   
 80 G4bool G4Absorber::FindAbsorbers(G4KineticTrac    
 81          G4KineticTrackVector & tgt)              
 82 {                                                 
 83 //  Find a closest ( in space) pair of Nucleon    
 84 //    pi+ can be absorbed on np or nn resultin    
 85 //    pi- can be absorbed on np or pp resultin    
 86                                                   
 87 // @GF: FindAbsorbers is unused, logic is seri    
 88                                                   
 89   G4KineticTrack * kt1 = NULL;                    
 90   G4KineticTrack * kt2 = NULL;                    
 91   G4double dist1 = DBL_MAX;   // dist to close    
 92   G4double dist2 = DBL_MAX;   // dist to next     
 93   G4double charge1 = 0;                           
 94 //  G4double charge2 = 0;   // charge2 is only    
 95   G4double charge0 = kt.GetDefinition()->GetPD    
 96   G4ThreeVector pos = kt.GetPosition();           
 97                                                   
 98   std::vector<G4KineticTrack *>::iterator iter    
 99   for(iter = tgt.begin(); iter != tgt.end(); +    
100   {                                               
101     G4KineticTrack * curr = *iter;                
102     G4double dist = (pos-curr->GetPosition()).    
103     if(dist >= dist2)                             
104       continue;                                   
105     if(dist < dist1)                              
106     {                                             
107       if(dist1 == DBL_MAX) // accept 1st as a     
108       {                                           
109   kt1 = curr;                                     
110   charge1 = kt1->GetDefinition()->GetPDGCharge    
111   dist1 = dist;                                   
112   continue;                                       
113       }                                           
114       if(dist2 == DBL_MAX) // accept the candi    
115       {       // @GF: should'nt we check if co    
116   kt2 = kt1;                                      
117 //  charge2 = charge1;                            
118   dist2 = dist1;                                  
119   kt1 = curr;                                     
120   charge1 = kt1->GetDefinition()->GetPDGCharge    
121   dist1 = dist;                                   
122   continue;                                       
123       }                                           
124 // test the compatibility with charge conserva    
125       G4double charge = curr->GetDefinition()-    
126       if((charge0+charge1+charge < 0.) || //te    
127    (charge0+charge1+charge) > 2*eplus)            
128       {  // incompatible: change kt1 with curr    
129   kt1 = curr;                                     
130   charge1 = charge;                               
131   dist1 = dist;                                   
132       }                                           
133       else                                        
134       { // compatible: change kt1 with curr an    
135   kt2 = kt1;                                      
136 //  charge2 = charge1;                            
137   dist2 = dist1;                                  
138   kt1 = curr;                                     
139   charge1 = charge;                               
140   dist1 = dist;                                   
141       }                                           
142       continue;                                   
143     }                                             
144 // here if dist1 < dist < dist2                   
145     if(dist2 == DBL_MAX) // accept the candida    
146     {                                             
147       kt2 = curr;                                 
148 //      charge2 = kt2->GetDefinition()->GetPDG    
149       dist2 = dist;                               
150       continue;                                   
151     }                                             
152 // test the compatibility with charge conserva    
153     G4double charge = curr->GetDefinition()->G    
154     if((charge0+charge1+charge < 0.) ||           
155        (charge0+charge1+charge) > 2*eplus)        
156       continue;   // incomatible: do nothing      
157 // compatible: change kt2 with curr               
158     kt2 = curr;                                   
159 //    charge2 = charge;                           
160     dist2 = dist;                                 
161   }                                               
162                                                   
163   theAbsorbers->clear(); // do not delete trac    
164   if((kt1 == NULL) || (kt2 == NULL))              
165     return false;                                 
166                                                   
167   theAbsorbers->push_back(kt1);                   
168   theAbsorbers->push_back(kt2);                   
169   return true;                                    
170 }                                                 
171                                                   
172                                                   
173                                                   
174 G4bool G4Absorber::FindProducts(G4KineticTrack    
175 {                                                 
176 // Choose the products type                       
177   const G4ParticleDefinition * prod1;             
178   const G4ParticleDefinition * prod2;             
179   G4KineticTrack * abs1 = (*theAbsorbers)[0];     
180   G4KineticTrack * abs2 = (*theAbsorbers)[1];     
181                                                   
182   G4double charge = kt.GetDefinition()->GetPDG    
183   if(charge == eplus)                             
184   { // a neutron become proton                    
185     prod1 = G4Proton::Proton();                   
186     if(abs1->GetDefinition() == G4Neutron::Neu    
187       prod2 = abs2->GetDefinition();              
188     else                                          
189       prod2 = G4Proton::Proton();                 
190   }                                               
191   else if(charge == -eplus)                       
192   { // a proton become neutron                    
193     prod1 = G4Neutron::Neutron();                 
194     if(abs1->GetDefinition() == G4Proton::Prot    
195       prod2 = abs2->GetDefinition();              
196     else                                          
197       prod2 = G4Neutron::Neutron();               
198   }                                               
199   else  // charge = 0: leave particle types un    
200   {                                               
201     prod1 = abs1->GetDefinition();                
202     prod2 = abs2->GetDefinition();                
203   }                                               
204                                                   
205 // Translate to the CMS frame                     
206   G4LorentzVector momLab = kt.Get4Momentum()+a    
207     abs2->Get4Momentum();                         
208   G4LorentzRotation toCMSFrame((-1)*momLab.boo    
209   G4LorentzRotation toLabFrame(momLab.boostVec    
210   G4LorentzVector momCMS = toCMSFrame*momLab;     
211                                                   
212 // Evaluate the final momentum of products        
213   G4double ms1 = prod1->GetPDGMass();             
214   G4double ms2 = prod2->GetPDGMass();             
215   G4double e0 = momCMS.e();                       
216   G4double squareP = (e0*e0*e0*e0-2*e0*e0*(ms1    
217     (ms2*ms2-ms1*ms1)*(ms2*ms2-ms1*ms1))/(4*e0    
218 //  if(squareP < 0)  // should never happen       
219 //    squareP = 0;                                
220   G4ThreeVector mom1CMS = GetRandomDirection()    
221   mom1CMS = std::sqrt(squareP)*mom1CMS;           
222   G4LorentzVector final4Mom1CMS(mom1CMS, std::    
223   G4LorentzVector final4Mom2CMS((-1)*mom1CMS,     
224                                                   
225 // Go back to the lab frame                       
226   G4LorentzVector mom1 = toLabFrame*final4Mom1    
227   G4LorentzVector mom2 = toLabFrame*final4Mom2    
228                                                   
229 // ------ debug                                   
230 /*                                                
231   G4LorentzVector temp = mom1+mom2;               
232                                                   
233   cout << (1/MeV)*momLab.x() << " " << (1/MeV)    
234        << (1/MeV)*momLab.z() << " " << (1/MeV)    
235        << (1/MeV)*momLab.vect().mag() << " " <    
236        << (1/MeV)*temp.x() << " " << (1/MeV)*t    
237        << (1/MeV)*temp.z() << " " << (1/MeV)*t    
238        << (1/MeV)*temp.vect().mag() << " " <<     
239        << (1/MeV)*std::sqrt(squareP) << endl;     
240                                                   
241 */                                                
242 // ------ end debug                               
243                                                   
244 // Build two new kinetic tracks and add to pro    
245   G4KineticTrack * kt1 = new G4KineticTrack(pr    
246               mom1);                              
247   G4KineticTrack * kt2 = new G4KineticTrack(pr    
248               mom2);                              
249 // ------ debug                                   
250 /*                                                
251   G4LorentzVector initialMom1 = abs1->Get4Mome    
252   G4LorentzVector initialMom2 = abs2->Get4Mome    
253   G4LorentzVector pion4MomCMS = toCMSFrame*kt.    
254   cout << (1/MeV)*initialMom1.x() << " " << (1    
255        << (1/MeV)*initialMom1.z() << " " << (1    
256        << (1/MeV)*initialMom1.vect().mag() <<     
257        << (1/MeV)*initialMom2.x() << " " << (1    
258        << (1/MeV)*initialMom2.z() << " " << (1    
259        << (1/MeV)*initialMom2.vect().mag() <<     
260        << (1/MeV)*mom1.x() << " " << (1/MeV)*m    
261        << (1/MeV)*mom1.z() << " " << (1/MeV)*m    
262        << (1/MeV)*mom1.vect().mag() << " "        
263        << (1/MeV)*mom2.x() << " " << (1/MeV)*m    
264        << (1/MeV)*mom2.z() << " " << (1/MeV)*m    
265        << (1/MeV)*mom2.vect().mag() << " "        
266        << (1/MeV)*pion4MomCMS.x() << " " << (1    
267        << (1/MeV)*pion4MomCMS.z() << " " << (1    
268        << (1/MeV)*pion4MomCMS.vect().mag() <<     
269        << (1/MeV)*final4Mom1CMS.x() << " " <<     
270        << (1/MeV)*final4Mom1CMS.z() << " " <<     
271        << (1/MeV)*final4Mom1CMS.vect().mag() <    
272        << (1/MeV)*final4Mom2CMS.x() << " " <<     
273        << (1/MeV)*final4Mom2CMS.z() << " " <<     
274        << (1/MeV)*final4Mom2CMS.vect().mag() <    
275 */                                                
276 // ------ end debug                               
277                                                   
278   theProducts->clear();                           
279   theProducts->push_back(kt1);                    
280   theProducts->push_back(kt2);                    
281   return true;                                    
282 }                                                 
283                                                   
284                                                   
285                                                   
286 G4ThreeVector G4Absorber::GetRandomDirection()    
287 {                                                 
288   G4double theta = 2.0*G4UniformRand()-1.0;       
289   theta = std::acos(theta);                       
290   G4double phi = G4UniformRand()*2*pi;            
291   G4ThreeVector direction(std::sin(theta)*std:    
292   return direction;                               
293 }                                                 
294                                                   
295                                                   
296                                                   
297                                                   
298                                                   
299                                                   
300                                                   
301