source: trunk/source/processes/hadronic/models/im_r_matrix/src/G4Scatterer.cc @ 1201

Last change on this file since 1201 was 819, checked in by garnier, 16 years ago

import all except CVS
File size: 14.8 KB
Line 
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// $Id: G4Scatterer.cc,v 1.13.2.4 2006/06/29 20:40:31 gunter Exp $ //
27//
28
29#include "globals.hh"
30#include <vector>
31#include "G4ios.hh"
32#include "G4Scatterer.hh"
33#include "G4KineticTrack.hh"
34#include "G4ThreeVector.hh"
35#include "G4LorentzRotation.hh"
36#include "G4LorentzVector.hh"
37
38#include "G4CollisionNN.hh"
39#include "G4CollisionPN.hh"
40#include "G4CollisionMesonBaryon.hh"
41
42#include "G4CollisionInitialState.hh"
43#include "G4HadTmpUtil.hh"
44#include "G4Pair.hh"
45
46
47// Declare the categories of collisions the Scatterer can handle
48typedef GROUP2(G4CollisionNN, G4CollisionMesonBaryon) theChannels;
49
50
51G4Scatterer::G4Scatterer()
52{
53  Register aR;
54  G4ForEach<theChannels>::Apply(&aR, &collisions);
55}
56
57
58G4Scatterer::~G4Scatterer()
59{
60  std::for_each(collisions.begin(), collisions.end(), G4Delete());
61  collisions.clear();
62}
63
64G4double G4Scatterer::GetTimeToInteraction(const G4KineticTrack& trk1, 
65                                           const G4KineticTrack& trk2)
66{
67  G4double time = DBL_MAX;
68    G4double distance_fast;
69  G4LorentzVector mom1 = trk1.GetTrackingMomentum();
70//  G4cout << "zcomp=" << std::abs(mom1.vect().unit().z() -1 ) << G4endl;
71  G4double collisionTime;
72 
73  if ( std::abs(mom1.vect().unit().z() -1 ) < 1e-6 ) 
74  {
75     G4ThreeVector position = trk2.GetPosition() - trk1.GetPosition();
76     G4double deltaz=position.z();
77     G4double velocity = mom1.z()/mom1.e() * c_light;
78     
79     collisionTime=deltaz/velocity;
80     distance_fast=position.x()*position.x() + position.y()*position.y();
81  } else {
82 
83    //  The nucleons of the nucleus are FROZEN, ie. do not move..
84
85    G4ThreeVector position = trk2.GetPosition() - trk1.GetPosition();   
86
87    G4ThreeVector velocity = mom1.vect()/mom1.e() * c_light;  // mom1.boostVector() will exit on slightly negative mass
88    collisionTime = (position * velocity) / velocity.mag2();    // can't divide by /c_light;
89    position -= velocity * collisionTime;
90    distance_fast=position.mag2();
91   
92//    if ( collisionTime>0 ) G4cout << " dis1/2 square" << dis1 <<" "<< dis2 << G4endl;
93//     collisionTime = GetTimeToClosestApproach(trk1,trk2);
94  }
95     if (collisionTime > 0)
96        { 
97           static const G4double maxCrossSection = 500*millibarn;
98           if(0.7*pi*distance_fast>maxCrossSection) return time;
99
100       
101           G4LorentzVector mom2(0,0,0,trk2.Get4Momentum().mag());
102
103//         G4ThreeVector momLab = mom1.vect();// frozen Nucleus - mom2.vect();
104//         G4ThreeVector posLab = trk1.GetPosition() - trk2.GetPosition();
105//         G4double disLab=posLab * posLab - (posLab*momLab) * (posLab*momLab) /(momLab.mag2());
106
107           G4LorentzRotation toCMSFrame((-1)*(mom1 + mom2).boostVector());
108           mom1 = toCMSFrame * mom1;
109           mom2 = toCMSFrame * mom2;
110
111           G4LorentzVector coordinate1(trk1.GetPosition(), 100.);
112           G4LorentzVector coordinate2(trk2.GetPosition(), 100.);
113           G4ThreeVector pos = ((toCMSFrame * coordinate1).vect() - 
114                                (toCMSFrame * coordinate2).vect());
115
116           G4ThreeVector mom = mom1.vect() - mom2.vect();
117
118          // Calculate the impact parameter
119
120           G4double distance = pos * pos - (pos*mom) * (pos*mom) / (mom.mag2());
121
122//     G4cout << " disDiff " << distance-disLab << " " << disLab
123//            << " " << std::abs(distance-disLab)/distance << G4endl
124//          << " mom/Lab " << mom << " " << momLab << G4endl
125//          << " pos/Lab " << pos << " " << posLab
126//          << G4endl;
127
128           // global optimization
129//         static const G4double maxCrossSection = 500*millibarn;
130           if(pi*distance>maxCrossSection) return time;
131           
132           // charged particles special
133           static const G4double maxChargedCrossSection = 200*millibarn;
134           if(std::abs(trk1.GetDefinition()->GetPDGCharge())>0.1 && 
135              std::abs(trk2.GetDefinition()->GetPDGCharge())>0.1 &&
136              pi*distance>maxChargedCrossSection) return time;
137             
138           G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
139           // neutrons special   
140           if(( trk1.GetDefinition() == G4Neutron::Neutron() ||
141                trk1.GetDefinition() == G4Neutron::Neutron() ) &&
142                sqrtS>1.91*GeV && pi*distance>maxChargedCrossSection) return time;
143
144/*
145 *        if(distance <= sqr(1.14*fermi))
146 *        {
147 *          time = collisionTime;
148 *       
149 * *
150 *  *        G4cout << "Scatter distance/time: " << std::sqrt(distance)/fermi <<
151 *  *            " / "<< time/ns << G4endl;
152 *  *         G4ThreeVector pos1=trk1.GetPosition();
153 *  *         G4ThreeVector pos2=trk2.GetPosition();
154 *  *         G4LorentzVector xmom1 = trk1.Get4Momentum();
155 *  *         G4LorentzVector xmom2 = trk2.Get4Momentum(); 
156 *  *         G4cout << "position1: " <<  pos1.x() << " " << pos1.y() << " "
157 *  *                   << pos1.z();
158 *  *         pos1+=(collisionTime*c_light/xmom1.e())*xmom1.vect();
159 *  *         G4cout << " straight line trprt: "
160 *  *                   <<  pos1.x() << " " << pos1.y() << " "
161 *  *                   <<  pos1.z()  << G4endl;
162 *  *         G4cout << "position2: " <<  pos2.x() << " " << pos2.y() << " "
163 *  *                   << pos2.z()  << G4endl;
164 *  *         G4cout << "straight line distance 2 fixed:" << (pos1-pos2).mag()/fermi << G4endl;
165 *  *         pos2+= (collisionTime*c_light/xmom2.e())*xmom2.vect();
166 *  *         G4cout<< " straight line trprt: "
167 *  *                   <<  pos2.x() << " " << pos2.y() << " "
168 *  *                   <<  pos2.z() << G4endl;
169 *  *         G4cout << "straight line distance :" << (pos1-pos2).mag()/fermi << G4endl;
170 *  *
171 *        }
172 *       
173 *        if(1)
174 *          return time;
175 */
176           
177           if ((trk1.GetActualMass()+trk2.GetActualMass()) > sqrtS) return time;
178
179           
180         
181          G4VCollision* collision = FindCollision(trk1,trk2);
182          G4double totalCrossSection;
183          // The cross section is interpreted geometrically as an area
184          // Two particles are assumed to collide if their distance is < (totalCrossSection/pi)
185
186          if (collision != 0)
187            {
188              totalCrossSection = collision->CrossSection(trk1,trk2);
189              if ( totalCrossSection > 0 )
190                {
191/*                  G4cout << " totalCrossection = "<< totalCrossSection << ", trk1/2, s, e-m: "
192 *                         << trk1.GetDefinition()->GetParticleName()
193 *                         << " / "
194 *                         << trk2.GetDefinition()->GetParticleName()
195 *                         << ", "
196 *                         << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()
197 *                         << ", "
198 *                         << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()-
199 *                             trk1.Get4Momentum().mag() - trk2.Get4Momentum().mag()
200 *                         << G4endl;
201 */
202                 if (distance <= totalCrossSection / pi)
203                   {
204                     time = collisionTime;
205                   }
206                } else
207                {
208
209                 // For debugging...
210 //                 G4cout << " totalCrossection = 0, trk1/2, s, e-m: "
211 //                        << trk1.GetDefinition()->GetParticleName()
212 //                        << " / "
213 //                        << trk2.GetDefinition()->GetParticleName()
214 //                        << ", "
215 //                        << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()
216 //                        << ", "
217 //                        << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()-
218 //                            trk1.Get4Momentum().mag() - trk2.Get4Momentum().mag()
219 //                        << G4endl;
220
221                }
222/*
223 *            if(distance <= sqr(5.*fermi))
224 *             {
225 *                G4cout << " distance,xsect, std::sqrt(xsect/pi) : " << std::sqrt(distance)/fermi
226 *                       << " " << totalCrossSection/sqr(fermi)
227 *                       << " " << std::sqrt(totalCrossSection / pi)/fermi << G4endl;
228 *             }
229 */
230
231            }
232          else
233            {
234              time = DBL_MAX;
235//            /*
236              // For debugging
237//hpw                 G4cout << "G4Scatterer - collision not found: "
238//hpw                << trk1.GetDefinition()->GetParticleName()
239//hpw                << " - "
240//hpw                << trk2.GetDefinition()->GetParticleName()
241//hpw                << G4endl;
242              // End of debugging
243//            */
244            }
245        }
246
247      else
248        {
249              /*
250          // For debugging
251          G4cout << "G4Scatterer - negative collisionTime"
252                 << ": collisionTime = " << collisionTime
253                 << ", position = " << position
254                 << ", velocity = " << velocity
255                 << G4endl;
256          // End of debugging
257              */
258        }
259
260  return time;
261}
262
263G4KineticTrackVector* G4Scatterer::Scatter(const G4KineticTrack& trk1, 
264                                              const G4KineticTrack& trk2) 
265{
266//   G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
267   G4LorentzVector pInitial=trk1.Get4Momentum() + trk2.Get4Momentum();
268   G4double energyBalance = pInitial.t();
269   G4double pxBalance = pInitial.vect().x();
270   G4double pyBalance = pInitial.vect().y();
271   G4double pzBalance = pInitial.vect().z();
272   G4int chargeBalance = G4lrint(trk1.GetDefinition()->GetPDGCharge()
273                       + trk2.GetDefinition()->GetPDGCharge());
274   G4int baryonBalance = trk1.GetDefinition()->GetBaryonNumber()
275                       + trk2.GetDefinition()->GetBaryonNumber();
276   
277   G4VCollision* collision = FindCollision(trk1,trk2);
278   if (collision != 0)
279   {
280     G4double aCrossSection = collision->CrossSection(trk1,trk2);
281     if (aCrossSection > 0.0) 
282     {
283
284
285        #ifdef debug_G4Scatterer
286        G4cout << "be4 FinalState 1(p,e,m): "
287        << trk1.Get4Momentum() << " "
288        << trk1.Get4Momentum().mag()
289        << ", 2: "
290        << trk2.Get4Momentum()<< " "
291        << trk2.Get4Momentum().mag() << " "
292        << G4endl;
293        #endif
294
295
296       G4KineticTrackVector* products = collision->FinalState(trk1,trk2);
297       if(!products || products->size() == 0) return products;
298
299        #ifdef debug_G4Scatterer
300       G4cout << "size of FS: "<<products->size()<<G4endl; 
301        #endif
302
303       G4KineticTrack *final= products->operator[](0);
304
305
306        #ifdef debug_G4Scatterer
307        G4cout << "    FinalState 1: "
308                << final->Get4Momentum()<< " "
309                << final->Get4Momentum().mag() ;
310        #endif
311
312        if(products->size() == 1) return products;
313        final=products->operator[](1);
314        #ifdef debug_G4Scatterer
315        G4cout << ", 2: "
316                << final->Get4Momentum() << " "
317                << final->Get4Momentum().mag() << " " << G4endl;
318        #endif
319
320       final= products->operator[](0);
321       G4LorentzVector pFinal=final->Get4Momentum();
322       if(products->size()==2)
323       {
324         final=products->operator[](1);
325         pFinal +=final->Get4Momentum();
326       }
327
328       #ifdef debug_G4Scatterer
329       if ( (pInitial-pFinal).mag() > 0.1*MeV )
330       {
331          G4cout << "G4Scatterer: momentum imbalance, pInitial= " <<pInitial << " pFinal= " <<pFinal<< G4endl;
332       }
333       G4cout << "Scatterer costh= " << trk1.Get4Momentum().vect().unit() *(products->operator[](0))->Get4Momentum().vect().unit()<< G4endl;
334       #endif
335       
336       for(size_t hpw=0; hpw<products->size(); hpw++)
337       {
338         energyBalance-=products->operator[](hpw)->Get4Momentum().t();
339         pxBalance-=products->operator[](hpw)->Get4Momentum().vect().x();
340         pyBalance-=products->operator[](hpw)->Get4Momentum().vect().y();
341         pzBalance-=products->operator[](hpw)->Get4Momentum().vect().z();
342         chargeBalance-=G4lrint(products->operator[](hpw)->GetDefinition()->GetPDGCharge());
343         baryonBalance-=products->operator[](hpw)->GetDefinition()->GetBaryonNumber();
344       }
345       if(getenv("ScattererEnergyBalanceCheck"))
346         std::cout << "DEBUGGING energy balance A: "
347                   <<energyBalance<<" "
348                   <<pxBalance<<" "
349                   <<pyBalance<<" "
350                   <<pzBalance<<" "
351                   <<chargeBalance<<" "
352                   <<baryonBalance<<" "
353                   <<G4endl;
354       if(chargeBalance !=0 )
355       {
356         G4cout << "track 1"<<trk1.GetDefinition()->GetParticleName()<<G4endl;
357         G4cout << "track 2"<<trk2.GetDefinition()->GetParticleName()<<G4endl;
358         for(size_t hpw=0; hpw<products->size(); hpw++)
359         {
360           G4cout << products->operator[](hpw)->GetDefinition()->GetParticleName()<<G4endl;
361         }
362         G4Exception("We have the problem");
363       }
364       return products;
365     } 
366   } 
367   
368   return NULL;
369}
370
371
372G4VCollision* G4Scatterer::FindCollision(const G4KineticTrack& trk1, 
373                                         const G4KineticTrack& trk2)
374{
375  G4VCollision* collisionInCharge = 0;
376
377  size_t i;
378  for (i=0; i<collisions.size(); i++)
379    {
380      G4VCollision* component = collisions[i];
381      if (component->IsInCharge(trk1,trk2))
382        {
383          collisionInCharge = component;
384          break;
385        }
386    }
387//    if(collisionInCharge)
388//    {
389//      G4cout << "found collision : "
390//         << collisionInCharge->GetName()<< " "
391//      << "for "
392//      << trk1.GetDefinition()->GetParticleName()<<" + "
393//      << trk2.GetDefinition()->GetParticleName()<<" "
394//      << G4endl;;
395//    }
396  return collisionInCharge;
397}
398 
399G4double G4Scatterer::GetCrossSection(const G4KineticTrack& trk1, 
400                                      const G4KineticTrack& trk2) 
401{
402   G4VCollision* collision = FindCollision(trk1,trk2);
403   G4double aCrossSection = 0;
404   if (collision != 0)
405   {
406     aCrossSection = collision->CrossSection(trk1,trk2);
407   }
408   return aCrossSection;
409}
410
411
412const std::vector<G4CollisionInitialState *> & G4Scatterer::
413GetCollisions(G4KineticTrack * aProjectile, 
414              std::vector<G4KineticTrack *> & someCandidates,
415              G4double aCurrentTime)
416{
417  theCollisions.clear();
418  std::vector<G4KineticTrack *>::iterator j=someCandidates.begin();
419  for(; j != someCandidates.end(); ++j)
420  {
421    G4double collisionTime = GetTimeToInteraction(*aProjectile, **j);
422    if(collisionTime == DBL_MAX)  // no collision
423    {
424      continue;
425    }
426    G4KineticTrackVector aTarget;
427    aTarget.push_back(*j);
428    theCollisions.push_back(
429      new G4CollisionInitialState(collisionTime+aCurrentTime, aProjectile, aTarget, this) );
430//      G4cerr <<" !!!!!! debug collisions "<<collisionTime<<" "<<pkt->GetDefinition()->GetParticleName()<<G4endl;     
431   }
432   return theCollisions;
433}
434
435
436G4KineticTrackVector * G4Scatterer::
437GetFinalState(G4KineticTrack * aProjectile, 
438              std::vector<G4KineticTrack *> & theTargets)
439{
440    G4KineticTrack target_reloc(*(theTargets[0]));
441    return Scatter(*aProjectile, target_reloc);
442}
Note: See TracBrowser for help on using the repository browser.