Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/im_r_matrix/src/G4Scatterer.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/im_r_matrix/src/G4Scatterer.cc (Version 11.3.0) and /processes/hadronic/models/im_r_matrix/src/G4Scatterer.cc (Version 11.2)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 //                                                 26 //
 27                                                    27 
 28 #include <vector>                                  28 #include <vector>
 29                                                    29 
 30 #include "globals.hh"                              30 #include "globals.hh"
 31 #include "G4PhysicalConstants.hh"                  31 #include "G4PhysicalConstants.hh"
 32 #include "G4SystemOfUnits.hh"                      32 #include "G4SystemOfUnits.hh"
 33 #include "G4ios.hh"                                33 #include "G4ios.hh"
 34 #include "G4Scatterer.hh"                          34 #include "G4Scatterer.hh"
 35 #include "G4KineticTrack.hh"                       35 #include "G4KineticTrack.hh"
 36 #include "G4ThreeVector.hh"                        36 #include "G4ThreeVector.hh"
 37 #include "G4LorentzRotation.hh"                    37 #include "G4LorentzRotation.hh"
 38 #include "G4LorentzVector.hh"                      38 #include "G4LorentzVector.hh"
 39                                                    39 
 40 #include "G4CollisionNN.hh"                        40 #include "G4CollisionNN.hh"
 41 #include "G4CollisionPN.hh"                        41 #include "G4CollisionPN.hh"
 42 #include "G4CollisionMesonBaryon.hh"               42 #include "G4CollisionMesonBaryon.hh"
 43                                                    43 
 44 #include "G4CollisionInitialState.hh"              44 #include "G4CollisionInitialState.hh"
 45 #include "G4HadTmpUtil.hh"                         45 #include "G4HadTmpUtil.hh"
 46 #include "G4Pair.hh"                               46 #include "G4Pair.hh"
 47 #include "G4AutoLock.hh"                           47 #include "G4AutoLock.hh"
 48                                                    48 
 49 //Mutex for control of shared resource             49 //Mutex for control of shared resource
 50 namespace  {                                       50 namespace  {
 51     G4Mutex collisions_mutex = G4MUTEX_INITIAL     51     G4Mutex collisions_mutex = G4MUTEX_INITIALIZER;
 52     G4bool setupDone = false;                      52     G4bool setupDone = false;
 53 }                                                  53 }
 54                                                    54 
 55 // Declare the categories of collisions the Sc     55 // Declare the categories of collisions the Scatterer can handle
 56 typedef GROUP2(G4CollisionNN, G4CollisionMeson     56 typedef GROUP2(G4CollisionNN, G4CollisionMesonBaryon) theChannels;
 57                                                    57 
 58 G4CollisionVector G4Scatterer::collisions;         58 G4CollisionVector G4Scatterer::collisions;
 59                                                    59 
 60 //--------------------------------------------     60 //----------------------------------------------------------------------------
 61                                                    61 
 62 G4Scatterer::G4Scatterer()                         62 G4Scatterer::G4Scatterer()
 63 {                                                  63 {
 64   G4AutoLock l(&collisions_mutex);                 64   G4AutoLock l(&collisions_mutex);
 65   if ( ! setupDone )                               65   if ( ! setupDone )
 66   {                                                66   {
 67       Register aR;                                 67       Register aR;
 68       G4ForEach<theChannels>::Apply(&aR, &coll     68       G4ForEach<theChannels>::Apply(&aR, &collisions);
 69       setupDone = true;                            69       setupDone = true;
 70   }                                                70   }
 71 }                                                  71 }
 72                                                    72 
 73 //--------------------------------------------     73 //----------------------------------------------------------------------------
 74                                                    74 
 75 G4Scatterer::~G4Scatterer()                        75 G4Scatterer::~G4Scatterer()
 76 {                                                  76 {
 77   G4AutoLock l(&collisions_mutex);                 77   G4AutoLock l(&collisions_mutex);
 78   std::for_each(collisions.begin(), collisions     78   std::for_each(collisions.begin(), collisions.end(), G4Delete());
 79   collisions.clear();                              79   collisions.clear();
 80 }                                                  80 }
 81                                                    81 
 82 //--------------------------------------------     82 //----------------------------------------------------------------------------
 83                                                    83 
 84 G4double G4Scatterer::GetTimeToInteraction(con     84 G4double G4Scatterer::GetTimeToInteraction(const G4KineticTrack& trk1,
 85              const G4KineticTrack& trk2) const     85              const G4KineticTrack& trk2) const
 86 {                                                  86 {
 87   G4double time = DBL_MAX;                         87   G4double time = DBL_MAX;
 88     G4double distance_fast;                        88     G4double distance_fast;
 89   G4LorentzVector mom1 = trk1.GetTrackingMomen     89   G4LorentzVector mom1 = trk1.GetTrackingMomentum();
 90 //  G4cout << "zcomp=" << std::abs(mom1.vect()     90 //  G4cout << "zcomp=" << std::abs(mom1.vect().unit().z() -1 ) << G4endl;
 91   G4double collisionTime;                          91   G4double collisionTime;
 92                                                    92 
 93   if ( std::abs(mom1.vect().unit().z() -1 ) <      93   if ( std::abs(mom1.vect().unit().z() -1 ) < 1e-6 )
 94   {                                                94   {
 95      G4ThreeVector position = trk2.GetPosition     95      G4ThreeVector position = trk2.GetPosition() - trk1.GetPosition();
 96      G4double deltaz=position.z();                 96      G4double deltaz=position.z();
 97      G4double velocity = mom1.z()/mom1.e() * c     97      G4double velocity = mom1.z()/mom1.e() * c_light;
 98                                                    98 
 99      collisionTime=deltaz/velocity;                99      collisionTime=deltaz/velocity;
100      distance_fast=position.x()*position.x() +    100      distance_fast=position.x()*position.x() + position.y()*position.y();
101   } else {                                        101   } else {
102                                                   102 
103     //  The nucleons of the nucleus are FROZEN    103     //  The nucleons of the nucleus are FROZEN, ie. do not move..
104                                                   104 
105     G4ThreeVector position = trk2.GetPosition(    105     G4ThreeVector position = trk2.GetPosition() - trk1.GetPosition();
106                                                   106 
107     G4ThreeVector velocity = mom1.vect()/mom1.    107     G4ThreeVector velocity = mom1.vect()/mom1.e() * c_light;  // mom1.boostVector() will exit on slightly negative mass
108     collisionTime = (position * velocity) / ve    108     collisionTime = (position * velocity) / velocity.mag2();    // can't divide by /c_light;
109     position -= velocity * collisionTime;         109     position -= velocity * collisionTime;
110     distance_fast=position.mag2();                110     distance_fast=position.mag2();
111                                                   111 
112 //    if ( collisionTime>0 ) G4cout << " dis1/    112 //    if ( collisionTime>0 ) G4cout << " dis1/2 square" << dis1 <<" "<< dis2 << G4endl;
113 //     collisionTime = GetTimeToClosestApproac    113 //     collisionTime = GetTimeToClosestApproach(trk1,trk2);
114   }                                               114   }
115      if (collisionTime > 0)                       115      if (collisionTime > 0)
116   {                                               116   {
117      static const G4double maxCrossSection = 5    117      static const G4double maxCrossSection = 500*millibarn;
118      if(0.7*pi*distance_fast>maxCrossSection)     118      if(0.7*pi*distance_fast>maxCrossSection) return time;
119                                                   119 
120                                                   120 
121            G4LorentzVector mom2(0,0,0,trk2.Get    121            G4LorentzVector mom2(0,0,0,trk2.Get4Momentum().mag());
122                                                   122 
123 //     G4ThreeVector momLab = mom1.vect();// f    123 //     G4ThreeVector momLab = mom1.vect();// frozen Nucleus - mom2.vect();
124 //     G4ThreeVector posLab = trk1.GetPosition    124 //     G4ThreeVector posLab = trk1.GetPosition() - trk2.GetPosition();
125 //     G4double disLab=posLab * posLab - (posL    125 //     G4double disLab=posLab * posLab - (posLab*momLab) * (posLab*momLab) /(momLab.mag2());
126                                                   126 
127      G4LorentzRotation toCMSFrame((-1)*(mom1 +    127      G4LorentzRotation toCMSFrame((-1)*(mom1 + mom2).boostVector());
128      mom1 = toCMSFrame * mom1;                    128      mom1 = toCMSFrame * mom1;
129      mom2 = toCMSFrame * mom2;                    129      mom2 = toCMSFrame * mom2;
130                                                   130 
131      G4LorentzVector coordinate1(trk1.GetPosit    131      G4LorentzVector coordinate1(trk1.GetPosition(), 100.);
132      G4LorentzVector coordinate2(trk2.GetPosit    132      G4LorentzVector coordinate2(trk2.GetPosition(), 100.);
133      G4ThreeVector pos = ((toCMSFrame * coordi    133      G4ThreeVector pos = ((toCMSFrame * coordinate1).vect() -
134         (toCMSFrame * coordinate2).vect());       134         (toCMSFrame * coordinate2).vect());
135                                                   135 
136      G4ThreeVector mom = mom1.vect() - mom2.ve    136      G4ThreeVector mom = mom1.vect() - mom2.vect();
137                                                   137 
138     // Calculate the impact parameter             138     // Calculate the impact parameter
139                                                   139 
140      G4double distance = pos * pos - (pos*mom)    140      G4double distance = pos * pos - (pos*mom) * (pos*mom) / (mom.mag2());
141                                                   141 
142 //     G4cout << " disDiff " << distance-disLa    142 //     G4cout << " disDiff " << distance-disLab << " " << disLab
143 //            << " " << std::abs(distance-disL    143 //            << " " << std::abs(distance-disLab)/distance << G4endl
144 //      << " mom/Lab " << mom << " " << momLab    144 //      << " mom/Lab " << mom << " " << momLab << G4endl
145 //      << " pos/Lab " << pos << " " << posLab    145 //      << " pos/Lab " << pos << " " << posLab
146 //      << G4endl;                                146 //      << G4endl;
147                                                   147 
148      if(pi*distance>maxCrossSection) return ti    148      if(pi*distance>maxCrossSection) return time;
149                                                   149 
150      // charged particles special                 150      // charged particles special
151      static const G4double maxChargedCrossSect    151      static const G4double maxChargedCrossSection = 200*millibarn;
152      if(std::abs(trk1.GetDefinition()->GetPDGC    152      if(std::abs(trk1.GetDefinition()->GetPDGCharge())>0.1 &&
153         std::abs(trk2.GetDefinition()->GetPDGC    153         std::abs(trk2.GetDefinition()->GetPDGCharge())>0.1 &&
154         pi*distance>maxChargedCrossSection) re    154         pi*distance>maxChargedCrossSection) return time;
155                                                   155 
156            G4double sqrtS = (trk1.Get4Momentum    156            G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
157      // neutrons special  pn is largest cross-    157      // neutrons special  pn is largest cross-section, but above 1.91 GeV is less than 200 mb
158      if(( trk1.GetDefinition() == G4Neutron::N    158      if(( trk1.GetDefinition() == G4Neutron::Neutron() ||
159           trk2.GetDefinition() == G4Neutron::N    159           trk2.GetDefinition() == G4Neutron::Neutron() ) &&
160          sqrtS>1.91*GeV && pi*distance>maxChar    160          sqrtS>1.91*GeV && pi*distance>maxChargedCrossSection) return time;
161                                                   161 
162 /*                                                162 /*
163  *    if(distance <= sqr(1.14*fermi))             163  *    if(distance <= sqr(1.14*fermi))
164  *    {                                           164  *    {
165  *      time = collisionTime;                     165  *      time = collisionTime;
166  *                                                166  *
167  * *                                              167  * *
168  *  *        G4cout << "Scatter distance/time:    168  *  *        G4cout << "Scatter distance/time: " << std::sqrt(distance)/fermi <<
169  *  *            " / "<< time/ns << G4endl;       169  *  *            " / "<< time/ns << G4endl;
170  *  *         G4ThreeVector pos1=trk1.GetPosit    170  *  *         G4ThreeVector pos1=trk1.GetPosition();
171  *  *         G4ThreeVector pos2=trk2.GetPosit    171  *  *         G4ThreeVector pos2=trk2.GetPosition();
172  *  *         G4LorentzVector xmom1 = trk1.Get    172  *  *         G4LorentzVector xmom1 = trk1.Get4Momentum();
173  *  *         G4LorentzVector xmom2 = trk2.Get    173  *  *         G4LorentzVector xmom2 = trk2.Get4Momentum();
174  *  *         G4cout << "position1: " <<  pos1    174  *  *         G4cout << "position1: " <<  pos1.x() << " " << pos1.y() << " "
175  *  *             << pos1.z();                    175  *  *             << pos1.z();
176  *  *         pos1+=(collisionTime*c_light/xmo    176  *  *         pos1+=(collisionTime*c_light/xmom1.e())*xmom1.vect();
177  *  *         G4cout << " straight line trprt:    177  *  *         G4cout << " straight line trprt: "
178  *  *             <<  pos1.x() << " " << pos1.    178  *  *             <<  pos1.x() << " " << pos1.y() << " "
179  *  *       <<  pos1.z()  << G4endl;              179  *  *       <<  pos1.z()  << G4endl;
180  *  *         G4cout << "position2: " <<  pos2    180  *  *         G4cout << "position2: " <<  pos2.x() << " " << pos2.y() << " "
181  *  *             << pos2.z()  << G4endl;         181  *  *             << pos2.z()  << G4endl;
182  *  *         G4cout << "straight line distanc    182  *  *         G4cout << "straight line distance 2 fixed:" << (pos1-pos2).mag()/fermi << G4endl;
183  *  *         pos2+= (collisionTime*c_light/xm    183  *  *         pos2+= (collisionTime*c_light/xmom2.e())*xmom2.vect();
184  *  *         G4cout<< " straight line trprt:     184  *  *         G4cout<< " straight line trprt: "
185  *  *             <<  pos2.x() << " " << pos2.    185  *  *             <<  pos2.x() << " " << pos2.y() << " "
186  *  *       <<  pos2.z() << G4endl;               186  *  *       <<  pos2.z() << G4endl;
187  *  *         G4cout << "straight line distanc    187  *  *         G4cout << "straight line distance :" << (pos1-pos2).mag()/fermi << G4endl;
188  *  *                                             188  *  *
189  *    }                                           189  *    }
190  *                                                190  *
191  *    if(1)                                       191  *    if(1)
192  *      return time;                              192  *      return time;
193  */                                               193  */
194                                                   194 
195      if ((trk1.GetActualMass()+trk2.GetActualM    195      if ((trk1.GetActualMass()+trk2.GetActualMass()) > sqrtS) return time;
196                                                   196 
197                                                   197 
198                                                   198 
199     const G4VCollision* collision = FindCollis    199     const G4VCollision* collision = FindCollision(trk1,trk2);
200     G4double totalCrossSection;                   200     G4double totalCrossSection;
201     // The cross section is interpreted geomet    201     // The cross section is interpreted geometrically as an area
202     // Two particles are assumed to collide if    202     // Two particles are assumed to collide if their distance is < (totalCrossSection/pi)
203                                                   203 
204     if (collision != 0)                           204     if (collision != 0)
205       {                                           205       {
206         totalCrossSection = collision->CrossSe    206         totalCrossSection = collision->CrossSection(trk1,trk2);
207         if ( totalCrossSection > 0 )              207         if ( totalCrossSection > 0 )
208           {                                       208           {
209 /*        G4cout << " totalCrossection = "<< t    209 /*        G4cout << " totalCrossection = "<< totalCrossSection << ", trk1/2, s, e-m: "
210  *               << trk1.GetDefinition()->GetP    210  *               << trk1.GetDefinition()->GetParticleName()
211  *         << " / "                               211  *         << " / "
212  *               << trk2.GetDefinition()->GetP    212  *               << trk2.GetDefinition()->GetParticleName()
213  *         << ", "                                213  *         << ", "
214  *         << (trk1.Get4Momentum()+trk2.Get4Mo    214  *         << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()
215  *         << ", "                                215  *         << ", "
216  *         << (trk1.Get4Momentum()+trk2.Get4Mo    216  *         << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()-
217  *             trk1.Get4Momentum().mag() - trk    217  *             trk1.Get4Momentum().mag() - trk2.Get4Momentum().mag()
218  *         << G4endl;                             218  *         << G4endl;
219  */                                               219  */
220      if (distance <= totalCrossSection / pi)      220      if (distance <= totalCrossSection / pi)
221        {                                          221        {
222          time = collisionTime;                    222          time = collisionTime;
223        }                                          223        }
224           } else                                  224           } else
225     {                                             225     {
226                                                   226 
227      // For debugging...                          227      // For debugging...
228  //       G4cout << " totalCrossection = 0, tr    228  //       G4cout << " totalCrossection = 0, trk1/2, s, e-m: "
229  //              << trk1.GetDefinition()->GetP    229  //              << trk1.GetDefinition()->GetParticleName()
230  //        << " / "                               230  //        << " / "
231  //              << trk2.GetDefinition()->GetP    231  //              << trk2.GetDefinition()->GetParticleName()
232  //        << ", "                                232  //        << ", "
233  //        << (trk1.Get4Momentum()+trk2.Get4Mo    233  //        << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()
234  //        << ", "                                234  //        << ", "
235  //        << (trk1.Get4Momentum()+trk2.Get4Mo    235  //        << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()-
236  //            trk1.Get4Momentum().mag() - trk    236  //            trk1.Get4Momentum().mag() - trk2.Get4Momentum().mag()
237  //        << G4endl;                             237  //        << G4endl;
238                                                   238 
239     }                                             239     }
240 /*                                                240 /*
241  *        if(distance <= sqr(5.*fermi))           241  *        if(distance <= sqr(5.*fermi))
242  *         {                                      242  *         {
243  *      G4cout << " distance,xsect, std::sqrt(    243  *      G4cout << " distance,xsect, std::sqrt(xsect/pi) : " << std::sqrt(distance)/fermi
244  *       << " " << totalCrossSection/sqr(fermi    244  *       << " " << totalCrossSection/sqr(fermi)
245  *       << " " << std::sqrt(totalCrossSection    245  *       << " " << std::sqrt(totalCrossSection / pi)/fermi << G4endl;
246  *         }                                      246  *         }
247  */                                               247  */
248                                                   248 
249       }                                           249       }
250     else                                          250     else
251       {                                           251       {
252         time = DBL_MAX;                           252         time = DBL_MAX;
253 //        /*                                      253 //        /*
254         // For debugging                          254         // For debugging
255 //hpw               G4cout << "G4Scatterer - c    255 //hpw               G4cout << "G4Scatterer - collision not found: "
256 //hpw        << trk1.GetDefinition()->GetParti    256 //hpw        << trk1.GetDefinition()->GetParticleName()
257 //hpw        << " - "                             257 //hpw        << " - "
258 //hpw        << trk2.GetDefinition()->GetParti    258 //hpw        << trk2.GetDefinition()->GetParticleName()
259 //hpw        << G4endl;                           259 //hpw        << G4endl;
260         // End of debugging                       260         // End of debugging
261 //        */                                      261 //        */
262       }                                           262       }
263   }                                               263   }
264                                                   264 
265       else                                        265       else
266   {                                               266   {
267         /*                                        267         /*
268     // For debugging                              268     // For debugging
269     G4cout << "G4Scatterer - negative collisio    269     G4cout << "G4Scatterer - negative collisionTime"
270      << ": collisionTime = " << collisionTime     270      << ": collisionTime = " << collisionTime
271      << ", position = " << position               271      << ", position = " << position
272      << ", velocity = " << velocity               272      << ", velocity = " << velocity
273      << G4endl;                                   273      << G4endl;
274     // End of debugging                           274     // End of debugging
275         */                                        275         */
276   }                                               276   }
277                                                   277 
278   return time;                                    278   return time;
279 }                                                 279 }
280                                                   280 
281 //--------------------------------------------    281 //----------------------------------------------------------------------------
282                                                   282 
283 G4KineticTrackVector* G4Scatterer::Scatter(con    283 G4KineticTrackVector* G4Scatterer::Scatter(const G4KineticTrack& trk1,
284                 const G4KineticTrack& trk2) co    284                 const G4KineticTrack& trk2) const
285 {                                                 285 {
286 //   G4double sqrtS = (trk1.Get4Momentum() + t    286 //   G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
287    G4LorentzVector pInitial=trk1.Get4Momentum(    287    G4LorentzVector pInitial=trk1.Get4Momentum() + trk2.Get4Momentum();
288    G4double energyBalance = pInitial.t();         288    G4double energyBalance = pInitial.t();
289    G4double pxBalance = pInitial.vect().x();      289    G4double pxBalance = pInitial.vect().x();
290    G4double pyBalance = pInitial.vect().y();      290    G4double pyBalance = pInitial.vect().y();
291    G4double pzBalance = pInitial.vect().z();      291    G4double pzBalance = pInitial.vect().z();
292    G4int chargeBalance = G4lrint(trk1.GetDefin    292    G4int chargeBalance = G4lrint(trk1.GetDefinition()->GetPDGCharge()
293                        + trk2.GetDefinition()-    293                        + trk2.GetDefinition()->GetPDGCharge());
294    G4int baryonBalance = trk1.GetDefinition()-    294    G4int baryonBalance = trk1.GetDefinition()->GetBaryonNumber()
295                        + trk2.GetDefinition()-    295                        + trk2.GetDefinition()->GetBaryonNumber();
296                                                   296 
297    const G4VCollision* collision = FindCollisi    297    const G4VCollision* collision = FindCollision(trk1,trk2);
298    if (collision != 0)                            298    if (collision != 0)
299    {                                              299    {
300      G4double aCrossSection = collision->Cross    300      G4double aCrossSection = collision->CrossSection(trk1,trk2);
301      if (aCrossSection > 0.0)                     301      if (aCrossSection > 0.0)
302      {                                            302      {
303                                                   303 
304                                                   304 
305     #ifdef debug_G4Scatterer                      305     #ifdef debug_G4Scatterer
306   G4cout << "be4 FinalState 1(p,e,m): "           306   G4cout << "be4 FinalState 1(p,e,m): "
307         << trk1.Get4Momentum() << " "             307         << trk1.Get4Momentum() << " "
308         << trk1.Get4Momentum().mag()              308         << trk1.Get4Momentum().mag()
309   << ", 2: "                                      309   << ", 2: "
310         << trk2.Get4Momentum()<< " "              310         << trk2.Get4Momentum()<< " "
311         << trk2.Get4Momentum().mag() << " "       311         << trk2.Get4Momentum().mag() << " "
312         << G4endl;                                312         << G4endl;
313   #endif                                          313   #endif
314                                                   314 
315                                                   315 
316        G4KineticTrackVector* products = collis    316        G4KineticTrackVector* products = collision->FinalState(trk1,trk2);
317        if(!products || products->size() == 0)     317        if(!products || products->size() == 0) return products;
318                                                   318 
319     #ifdef debug_G4Scatterer                      319     #ifdef debug_G4Scatterer
320        G4cout << "size of FS: "<<products->siz    320        G4cout << "size of FS: "<<products->size()<<G4endl;
321   #endif                                          321   #endif
322                                                   322 
323        G4KineticTrack *final= products->operat    323        G4KineticTrack *final= products->operator[](0);
324                                                   324 
325                                                   325 
326     #ifdef debug_G4Scatterer                      326     #ifdef debug_G4Scatterer
327         G4cout << "    FinalState 1: "            327         G4cout << "    FinalState 1: "
328     << final->Get4Momentum()<< " "                328     << final->Get4Momentum()<< " "
329     << final->Get4Momentum().mag() ;              329     << final->Get4Momentum().mag() ;
330   #endif                                          330   #endif
331                                                   331 
332         if(products->size() == 1) return produ    332         if(products->size() == 1) return products;
333   final=products->operator[](1);                  333   final=products->operator[](1);
334     #ifdef debug_G4Scatterer                      334     #ifdef debug_G4Scatterer
335   G4cout << ", 2: "                               335   G4cout << ", 2: "
336     << final->Get4Momentum() << " "               336     << final->Get4Momentum() << " "
337           << final->Get4Momentum().mag() << "     337           << final->Get4Momentum().mag() << " " << G4endl;
338   #endif                                          338   #endif
339                                                   339 
340        final= products->operator[](0);            340        final= products->operator[](0);
341        G4LorentzVector pFinal=final->Get4Momen    341        G4LorentzVector pFinal=final->Get4Momentum();
342        if(products->size()==2)                    342        if(products->size()==2)
343        {                                          343        {
344          final=products->operator[](1);           344          final=products->operator[](1);
345          pFinal +=final->Get4Momentum();          345          pFinal +=final->Get4Momentum();
346        }                                          346        }
347                                                   347 
348        #ifdef debug_G4Scatterer                   348        #ifdef debug_G4Scatterer
349        if ( (pInitial-pFinal).mag() > 0.1*MeV     349        if ( (pInitial-pFinal).mag() > 0.1*MeV )
350        {                                          350        {
351           G4cout << "G4Scatterer: momentum imb    351           G4cout << "G4Scatterer: momentum imbalance, pInitial= " <<pInitial << " pFinal= " <<pFinal<< G4endl;
352        }                                          352        }
353        G4cout << "Scatterer costh= " << trk1.G    353        G4cout << "Scatterer costh= " << trk1.Get4Momentum().vect().unit() *(products->operator[](0))->Get4Momentum().vect().unit()<< G4endl;
354        #endif                                     354        #endif
355                                                   355 
356        for(size_t hpw=0; hpw<products->size();    356        for(size_t hpw=0; hpw<products->size(); hpw++)
357        {                                          357        {
358          energyBalance-=products->operator[](h    358          energyBalance-=products->operator[](hpw)->Get4Momentum().t();
359          pxBalance-=products->operator[](hpw)-    359          pxBalance-=products->operator[](hpw)->Get4Momentum().vect().x();
360          pyBalance-=products->operator[](hpw)-    360          pyBalance-=products->operator[](hpw)->Get4Momentum().vect().y();
361          pzBalance-=products->operator[](hpw)-    361          pzBalance-=products->operator[](hpw)->Get4Momentum().vect().z();
362    chargeBalance-=G4lrint(products->operator[]    362    chargeBalance-=G4lrint(products->operator[](hpw)->GetDefinition()->GetPDGCharge());
363          baryonBalance-=products->operator[](h    363          baryonBalance-=products->operator[](hpw)->GetDefinition()->GetBaryonNumber();
364        }                                          364        }
365        if(std::getenv("ScattererEnergyBalanceC    365        if(std::getenv("ScattererEnergyBalanceCheck"))
366          std::cout << "DEBUGGING energy balanc    366          std::cout << "DEBUGGING energy balance A: "
367              <<energyBalance<<" "                 367              <<energyBalance<<" "
368              <<pxBalance<<" "                     368              <<pxBalance<<" "
369              <<pyBalance<<" "                     369              <<pyBalance<<" "
370              <<pzBalance<<" "                     370              <<pzBalance<<" "
371        <<chargeBalance<<" "                       371        <<chargeBalance<<" "
372        <<baryonBalance<<" "                       372        <<baryonBalance<<" "
373        <<G4endl;                                  373        <<G4endl;
374        if(chargeBalance !=0 )                     374        if(chargeBalance !=0 )
375        {                                          375        {
376          G4cout << "track 1"<<trk1.GetDefiniti    376          G4cout << "track 1"<<trk1.GetDefinition()->GetParticleName()<<G4endl;
377          G4cout << "track 2"<<trk2.GetDefiniti    377          G4cout << "track 2"<<trk2.GetDefinition()->GetParticleName()<<G4endl;
378          for(size_t hpw=0; hpw<products->size(    378          for(size_t hpw=0; hpw<products->size(); hpw++)
379          {                                        379          {
380             G4cout << products->operator[](hpw    380             G4cout << products->operator[](hpw)->GetDefinition()->GetParticleName()<<G4endl;
381          }                                        381          }
382          G4Exception("G4Scatterer", "im_r_matr    382          G4Exception("G4Scatterer", "im_r_matrix001", FatalException,
383              "Problem in ChargeBalance");         383              "Problem in ChargeBalance");
384        }                                          384        }
385        return products;                           385        return products;
386      }                                            386      }
387    }                                              387    }
388                                                   388 
389    return NULL;                                   389    return NULL;
390 }                                                 390 }
391                                                   391 
392 //--------------------------------------------    392 //----------------------------------------------------------------------------
393                                                   393 
394 const G4VCollision* G4Scatterer::FindCollision    394 const G4VCollision* G4Scatterer::FindCollision(const G4KineticTrack& trk1,
395            const G4KineticTrack& trk2) const      395            const G4KineticTrack& trk2) const
396 {                                                 396 {
397   G4VCollision* collisionInCharge = 0;            397   G4VCollision* collisionInCharge = 0;
398                                                   398 
399   size_t i;                                       399   size_t i;
400   for (i=0; i<collisions.size(); i++)             400   for (i=0; i<collisions.size(); i++)
401     {                                             401     {
402       G4VCollision* component = collisions[i];    402       G4VCollision* component = collisions[i];
403       if (component->IsInCharge(trk1,trk2))       403       if (component->IsInCharge(trk1,trk2))
404   {                                               404   {
405     collisionInCharge = component;                405     collisionInCharge = component;
406     break;                                        406     break;
407   }                                               407   }
408     }                                             408     }
409 //    if(collisionInCharge)                       409 //    if(collisionInCharge)
410 //    {                                           410 //    {
411 //      G4cout << "found collision : "            411 //      G4cout << "found collision : "
412 //         << collisionInCharge->GetName()<< "    412 //         << collisionInCharge->GetName()<< " "
413 //  << "for "                                     413 //  << "for "
414 //  << trk1.GetDefinition()->GetParticleName()    414 //  << trk1.GetDefinition()->GetParticleName()<<" + "
415 //  << trk2.GetDefinition()->GetParticleName()    415 //  << trk2.GetDefinition()->GetParticleName()<<" "
416 //  << G4endl;;                                   416 //  << G4endl;;
417 //    }                                           417 //    }
418   return collisionInCharge;                       418   return collisionInCharge;
419 }                                                 419 }
420                                                   420 
421 //--------------------------------------------    421 //----------------------------------------------------------------------------
422                                                   422 
423 G4double G4Scatterer::GetCrossSection(const G4    423 G4double G4Scatterer::GetCrossSection(const G4KineticTrack& trk1,
424               const G4KineticTrack& trk2) cons    424               const G4KineticTrack& trk2) const
425 {                                                 425 {
426    const G4VCollision* collision = FindCollisi    426    const G4VCollision* collision = FindCollision(trk1,trk2);
427    G4double aCrossSection = 0;                    427    G4double aCrossSection = 0;
428    if (collision != 0)                            428    if (collision != 0)
429    {                                              429    {
430      aCrossSection = collision->CrossSection(t    430      aCrossSection = collision->CrossSection(trk1,trk2);
431    }                                              431    }
432    return aCrossSection;                          432    return aCrossSection;
433 }                                                 433 }
434                                                   434 
435 //--------------------------------------------    435 //----------------------------------------------------------------------------
436                                                   436 
437 const std::vector<G4CollisionInitialState *> &    437 const std::vector<G4CollisionInitialState *> & G4Scatterer::
438 GetCollisions(G4KineticTrack * aProjectile,       438 GetCollisions(G4KineticTrack * aProjectile,
439               std::vector<G4KineticTrack *> &     439               std::vector<G4KineticTrack *> & someCandidates,
440         G4double aCurrentTime)                    440         G4double aCurrentTime)
441 {                                                 441 {
442   theCollisions.clear();                          442   theCollisions.clear();
443   std::vector<G4KineticTrack *>::iterator j=so    443   std::vector<G4KineticTrack *>::iterator j=someCandidates.begin();
444   for(; j != someCandidates.end(); ++j)           444   for(; j != someCandidates.end(); ++j)
445   {                                               445   {
446     G4double collisionTime = GetTimeToInteract    446     G4double collisionTime = GetTimeToInteraction(*aProjectile, **j);
447     if(collisionTime == DBL_MAX)  // no collis    447     if(collisionTime == DBL_MAX)  // no collision
448     {                                             448     {
449       continue;                                   449       continue;
450     }                                             450     }
451     G4KineticTrackVector aTarget;                 451     G4KineticTrackVector aTarget;
452     aTarget.push_back(*j);                        452     aTarget.push_back(*j);
453     theCollisions.push_back(                      453     theCollisions.push_back(
454       new G4CollisionInitialState(collisionTim    454       new G4CollisionInitialState(collisionTime+aCurrentTime, aProjectile, aTarget, this) );
455 //      G4cerr <<" !!!!!! debug collisions "<<    455 //      G4cerr <<" !!!!!! debug collisions "<<collisionTime<<" "<<pkt->GetDefinition()->GetParticleName()<<G4endl;
456    }                                              456    }
457    return theCollisions;                          457    return theCollisions;
458 }                                                 458 }
459                                                   459 
460                                                   460 
461 G4KineticTrackVector * G4Scatterer::              461 G4KineticTrackVector * G4Scatterer::
462 GetFinalState(G4KineticTrack * aProjectile,       462 GetFinalState(G4KineticTrack * aProjectile,
463         std::vector<G4KineticTrack *> & theTar    463         std::vector<G4KineticTrack *> & theTargets)
464 {                                                 464 {
465     G4KineticTrack target_reloc(*(theTargets[0    465     G4KineticTrack target_reloc(*(theTargets[0]));
466     return Scatter(*aProjectile, target_reloc)    466     return Scatter(*aProjectile, target_reloc);
467 }                                                 467 }
468 //--------------------------------------------    468 //----------------------------------------------------------------------------
469                                                   469