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 9.5.p2)


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