Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/qmd/src/G4QMDCollision.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/qmd/src/G4QMDCollision.cc (Version 11.3.0) and /processes/hadronic/models/qmd/src/G4QMDCollision.cc (Version 10.1.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 // 080602 Fix memory leaks by T. Koi               26 // 080602 Fix memory leaks by T. Koi 
 27 // 081120 Add deltaT in signature of CalKinema     27 // 081120 Add deltaT in signature of CalKinematicsOfBinaryCollisions
 28 //        Add several required updating of Mea     28 //        Add several required updating of Mean Filed 
 29 //        Modified handling of absorption case     29 //        Modified handling of absorption case by T. Koi
 30 // 090126 Fix in absorption case by T. Koi         30 // 090126 Fix in absorption case by T. Koi  
 31 // 090331 Fix for gamma participant by T. Koi      31 // 090331 Fix for gamma participant by T. Koi 
 32 //                                                 32 //
 33 #include "G4QMDCollision.hh"                       33 #include "G4QMDCollision.hh"
 34 #include "G4Scatterer.hh"                          34 #include "G4Scatterer.hh"
 35 #include "G4Pow.hh"                            << 
 36 #include "G4Exp.hh"                            << 
 37 #include "G4Log.hh"                            << 
 38 #include "G4PhysicalConstants.hh"                  35 #include "G4PhysicalConstants.hh"
 39 #include "G4SystemOfUnits.hh"                      36 #include "G4SystemOfUnits.hh"
 40 #include "Randomize.hh"                            37 #include "Randomize.hh"
 41                                                    38 
 42 G4QMDCollision::G4QMDCollision()                   39 G4QMDCollision::G4QMDCollision()
 43 : fdeltar ( 4.0 )                                  40 : fdeltar ( 4.0 )
 44 , fbcmax0 ( 1.323142 ) // NN maximum impact pa     41 , fbcmax0 ( 1.323142 ) // NN maximum impact parameter
 45 , fbcmax1 ( 2.523 )    // others maximum impac     42 , fbcmax1 ( 2.523 )    // others maximum impact parameter
 46 // , sig0 ( 55 )   // NN cross section             43 // , sig0 ( 55 )   // NN cross section
 47 //110617 fix for gcc 4.6 compilation warnings      44 //110617 fix for gcc 4.6 compilation warnings 
 48 //, sig1 ( 200 )  // others cross section          45 //, sig1 ( 200 )  // others cross section
 49 , fepse ( 0.0001 )                                 46 , fepse ( 0.0001 )
 50 {                                                  47 {
 51    //These two pointers will be set through Se     48    //These two pointers will be set through SetMeanField method
 52    theSystem=NULL;                                 49    theSystem=NULL;
 53    theMeanField=NULL;                              50    theMeanField=NULL;
 54    theScatterer = new G4Scatterer();               51    theScatterer = new G4Scatterer();
 55 }                                                  52 }
 56                                                    53 
 57 /*                                                 54 /*
 58 G4QMDCollision::G4QMDCollision( const G4QMDCol     55 G4QMDCollision::G4QMDCollision( const G4QMDCollision& obj )
 59 : fdeltar ( obj.fdeltar )                          56 : fdeltar ( obj.fdeltar )
 60 , fbcmax0 ( obj.fbcmax0 ) // NN maximum impact     57 , fbcmax0 ( obj.fbcmax0 ) // NN maximum impact parameter
 61 , fbcmax1 ( obj.fbcmax1 )    // others maximum     58 , fbcmax1 ( obj.fbcmax1 )    // others maximum impact parameter
 62 , fepse ( obj.fepse )                              59 , fepse ( obj.fepse )
 63 {                                                  60 {
 64                                                    61    
 65    if ( obj.theSystem != NULL ) {                  62    if ( obj.theSystem != NULL ) {
 66       theSystem = new G4QMDSystem;                 63       theSystem = new G4QMDSystem;
 67       *theSystem = *obj.theSystem;                 64       *theSystem = *obj.theSystem;
 68    } else {                                        65    } else {
 69       theSystem = NULL;                            66       theSystem = NULL;
 70    }                                               67    }
 71    if ( obj.theMeanField != NULL ) {               68    if ( obj.theMeanField != NULL ) {
 72       theMeanField = new G4QMDMeanField;           69       theMeanField = new G4QMDMeanField;
 73       *theMeanField = *obj.theMeanField;           70       *theMeanField = *obj.theMeanField;
 74    } else {                                        71    } else {
 75       theMeanField = NULL;                         72       theMeanField = NULL;
 76    }                                               73    }
 77    theScatterer = new G4Scatterer();               74    theScatterer = new G4Scatterer();
 78    *theScatterer = *obj.theScatterer;              75    *theScatterer = *obj.theScatterer;
 79 }                                                  76 }
 80                                                    77 
 81 G4QMDCollision & G4QMDCollision::operator= ( c     78 G4QMDCollision & G4QMDCollision::operator= ( const G4QMDCollision& obj)
 82 {                                                  79 {
 83    fdeltar = obj.fdeltar;                          80    fdeltar = obj.fdeltar;
 84    fbcmax0 = obj.fbcmax1;                          81    fbcmax0 = obj.fbcmax1;
 85    fepse = obj.fepse;                              82    fepse = obj.fepse;
 86                                                    83 
 87    if ( obj.theSystem != NULL ) {                  84    if ( obj.theSystem != NULL ) {
 88       delete theSystem;                            85       delete theSystem;
 89       theSystem = new G4QMDSystem;                 86       theSystem = new G4QMDSystem;
 90       *theSystem = *obj.theSystem;                 87       *theSystem = *obj.theSystem;
 91    } else {                                        88    } else {
 92       theSystem = NULL;                            89       theSystem = NULL;
 93    }                                               90    }
 94    if ( obj.theMeanField != NULL ) {               91    if ( obj.theMeanField != NULL ) {
 95       delete theMeanField;                         92       delete theMeanField;
 96       theMeanField = new G4QMDMeanField;           93       theMeanField = new G4QMDMeanField;
 97       *theMeanField = *obj.theMeanField;           94       *theMeanField = *obj.theMeanField;
 98    } else {                                        95    } else {
 99       theMeanField = NULL;                         96       theMeanField = NULL;
100    }                                               97    }
101    delete theScatterer;                            98    delete theScatterer;
102    theScatterer = new G4Scatterer();               99    theScatterer = new G4Scatterer();
103    *theScatterer = *obj.theScatterer;             100    *theScatterer = *obj.theScatterer;
104                                                   101 
105    return *this;                                  102    return *this;
106 }                                                 103 }
107 */                                                104 */
108                                                   105 
109                                                   106 
110 G4QMDCollision::~G4QMDCollision()                 107 G4QMDCollision::~G4QMDCollision()
111 {                                                 108 {
112    //if ( theSystem != NULL ) delete theSystem    109    //if ( theSystem != NULL ) delete theSystem;
113    //if ( theMeanField != NULL ) delete theMea    110    //if ( theMeanField != NULL ) delete theMeanField;
114    delete theScatterer;                           111    delete theScatterer;
115 }                                                 112 }
116                                                   113 
117                                                   114 
118 void G4QMDCollision::CalKinematicsOfBinaryColl    115 void G4QMDCollision::CalKinematicsOfBinaryCollisions( G4double dt )
119 {                                                 116 {
120    G4double deltaT = dt;                          117    G4double deltaT = dt; 
121                                                   118 
122    G4int n = theSystem->GetTotalNumberOfPartic    119    G4int n = theSystem->GetTotalNumberOfParticipant(); 
123 //081118                                          120 //081118
124    //G4int nb = 0;                                121    //G4int nb = 0;
125    for ( G4int i = 0 ; i < n ; i++ )              122    for ( G4int i = 0 ; i < n ; i++ )
126    {                                              123    {
127       theSystem->GetParticipant( i )->UnsetHit    124       theSystem->GetParticipant( i )->UnsetHitMark();
128       theSystem->GetParticipant( i )->UnsetHit    125       theSystem->GetParticipant( i )->UnsetHitMark();
129       //nb += theSystem->GetParticipant( i )->    126       //nb += theSystem->GetParticipant( i )->GetBaryonNumber();
130    }                                              127    }
131    //G4cout << "nb = " << nb << " n = " << n <    128    //G4cout << "nb = " << nb << " n = " << n << G4endl;
132                                                   129 
133                                                   130 
134 //071101                                          131 //071101
135    for ( G4int i = 0 ; i < n ; i++ )              132    for ( G4int i = 0 ; i < n ; i++ )
136    {                                              133    {
137                                                   134 
138       //std::cout << i << " " << theSystem->Ge    135       //std::cout << i << " " << theSystem->GetParticipant( i )->GetDefinition()->GetParticleName() << " " << theSystem->GetParticipant( i )->GetPosition() << std::endl;
139                                                   136 
140       if ( theSystem->GetParticipant( i )->Get    137       if ( theSystem->GetParticipant( i )->GetDefinition()->IsShortLived() )
141       {                                           138       {
142                                                   139 
143          G4bool decayed = false;                  140          G4bool decayed = false; 
144                                                   141 
145          const G4ParticleDefinition* pd0 = the    142          const G4ParticleDefinition* pd0 = theSystem->GetParticipant( i )->GetDefinition();
146          G4ThreeVector p0 = theSystem->GetPart    143          G4ThreeVector p0 = theSystem->GetParticipant( i )->GetMomentum();
147          G4ThreeVector r0 = theSystem->GetPart    144          G4ThreeVector r0 = theSystem->GetParticipant( i )->GetPosition();
148                                                   145 
149          G4LorentzVector p40 = theSystem->GetP    146          G4LorentzVector p40 = theSystem->GetParticipant( i )->Get4Momentum();
150                                                   147 
151          G4double eini = theMeanField->GetTota    148          G4double eini = theMeanField->GetTotalPotential() + p40.e();
152                                                   149 
153          G4int n0 = theSystem->GetTotalNumberO    150          G4int n0 = theSystem->GetTotalNumberOfParticipant(); 
154          G4int i0 = 0;                            151          G4int i0 = 0; 
155                                                   152 
156          G4bool isThisEnergyOK = false;           153          G4bool isThisEnergyOK = false;
157                                                   154 
158          G4int maximumNumberOfTrial=4;         << 155          for ( G4int ii = 0 ; ii < 4 ; ii++ )
159          for ( G4int ii = 0 ; ii < maximumNumb << 
160          {                                        156          { 
161                                                   157 
162             //G4LorentzVector p4 = theSystem->    158             //G4LorentzVector p4 = theSystem->GetParticipant( i )->Get4Momentum();
163             G4LorentzVector p400 = p40;           159             G4LorentzVector p400 = p40;
164                                                   160 
165             p400 *= GeV;                          161             p400 *= GeV;
166             //G4KineticTrack kt( theSystem->Ge    162             //G4KineticTrack kt( theSystem->GetParticipant( i )->GetDefinition() , 0.0 , (theSystem->GetParticipant( i )->GetPosition())*fermi , p4 );
167             G4KineticTrack kt( pd0 , 0.0 , r0*    163             G4KineticTrack kt( pd0 , 0.0 , r0*fermi , p400 );
168             //std::cout << "G4KineticTrack " <    164             //std::cout << "G4KineticTrack " << i << " " <<  kt.GetDefinition()->GetParticleName() <<  kt.GetPosition() << std::endl;
169             G4KineticTrackVector* secs = NULL;    165             G4KineticTrackVector* secs = NULL;
170             secs = kt.Decay();                    166             secs = kt.Decay();
171             G4int id = 0;                         167             G4int id = 0;
172             G4double et = 0;                      168             G4double et = 0;
173             if ( secs )                           169             if ( secs )
174             {                                     170             {
175                for ( G4KineticTrackVector::ite    171                for ( G4KineticTrackVector::iterator it 
176                      = secs->begin() ; it != s    172                      = secs->begin() ; it != secs->end() ; it++ )
177                {                                  173                {
178 /*                                                174 /*
179                   G4cout << "G4KineticTrack"      175                   G4cout << "G4KineticTrack" 
180                   << " " << (*it)->GetDefiniti    176                   << " " << (*it)->GetDefinition()->GetParticleName()
181                   << " " << (*it)->Get4Momentu    177                   << " " << (*it)->Get4Momentum()
182                   << " " << (*it)->GetPosition    178                   << " " << (*it)->GetPosition()/fermi
183                   << G4endl;                      179                   << G4endl;
184 */                                                180 */
185                    if ( id == 0 )                 181                    if ( id == 0 ) 
186                    {                              182                    {
187                       theSystem->GetParticipan    183                       theSystem->GetParticipant( i )->SetDefinition( (*it)->GetDefinition() );
188                       theSystem->GetParticipan    184                       theSystem->GetParticipant( i )->SetMomentum( (*it)->Get4Momentum().v()/GeV );
189                       theSystem->GetParticipan    185                       theSystem->GetParticipant( i )->SetPosition( (*it)->GetPosition()/fermi );
190                       //theMeanField->Cal2Body    186                       //theMeanField->Cal2BodyQuantities( i ); 
191                       et += (*it)->Get4Momentu    187                       et += (*it)->Get4Momentum().e()/GeV;
192                    }                              188                    }
193                    if ( id > 0 )                  189                    if ( id > 0 )
194                    {                              190                    {
195                       // Append end;              191                       // Append end;
196                       theSystem->SetParticipan    192                       theSystem->SetParticipant ( new G4QMDParticipant ( (*it)->GetDefinition() , (*it)->Get4Momentum().v()/GeV , (*it)->GetPosition()/fermi ) );
197                       et += (*it)->Get4Momentu    193                       et += (*it)->Get4Momentum().e()/GeV;
198                       if ( id > 1 )               194                       if ( id > 1 )
199                       {                           195                       {
200                          //081118                 196                          //081118
201                          //G4cout << "G4QMDCol    197                          //G4cout << "G4QMDCollision id >2; id= " << id  << G4endl; 
202                       }                           198                       }
203                    }                              199                    }
204                    id++; // number of daughter    200                    id++; // number of daughter particles
205                                                   201 
206                    delete *it;                    202                    delete *it;
207                }                                  203                }
208                                                   204 
209                theMeanField->Update();            205                theMeanField->Update();
210                i0 = id-1; // 0 enter to i         206                i0 = id-1; // 0 enter to i
211                                                   207 
212                delete secs;                       208                delete secs;
213             }                                     209             }
214                                                   210 
215 //          EnergyCheck                           211 //          EnergyCheck  
216                                                   212 
217             G4double efin = theMeanField->GetT    213             G4double efin = theMeanField->GetTotalPotential() + et; 
218             //std::cout <<  std::abs ( eini -     214             //std::cout <<  std::abs ( eini - efin ) - fepse << std::endl; 
219 //            std::cout <<  std::abs ( eini -     215 //            std::cout <<  std::abs ( eini - efin ) - fepse*10 << std::endl; 
220 //                                           *    216 //                                           *10 TK  
221             if ( std::abs ( eini - efin ) < fe    217             if ( std::abs ( eini - efin ) < fepse*10 ) 
222             {                                     218             {
223                // Energy OK                       219                // Energy OK 
224                isThisEnergyOK = true;             220                isThisEnergyOK = true;
225                break;                             221                break; 
226             }                                     222             }
227             else                                  223             else
228             {                                     224             {
229                                                   225 
230                theSystem->GetParticipant( i )-    226                theSystem->GetParticipant( i )->SetDefinition( pd0 );
231                theSystem->GetParticipant( i )-    227                theSystem->GetParticipant( i )->SetPosition( r0 );
232                theSystem->GetParticipant( i )-    228                theSystem->GetParticipant( i )->SetMomentum( p0 );
233                                                   229 
234                //for ( G4int i0i = 0 ; i0i < i << 230                for ( G4int i0i = 0 ; i0i < id-1 ; i0i++ )
235                //160210 deletion must be done  << 231                {
236                for ( G4int i0i = id-2 ; 0 <= i << 
237                   //081118                        232                   //081118
238                   //std::cout << "Decay Energi    233                   //std::cout << "Decay Energitically Blocked deleteing " << i0i+n0  << std::endl;
239                   theSystem->DeleteParticipant    234                   theSystem->DeleteParticipant( i0i+n0 );
240                }                                  235                }
241                //081103                           236                //081103
242                theMeanField->Update();            237                theMeanField->Update();
243             }                                     238             }
244                                                   239 
245          }                                        240          }
246                                                   241 
247                                                   242 
248 //       Pauli Check                              243 //       Pauli Check 
249          if ( isThisEnergyOK == true )            244          if ( isThisEnergyOK == true )
250          {                                        245          {
251             if ( theMeanField->IsPauliBlocked     246             if ( theMeanField->IsPauliBlocked ( i ) != true ) 
252             {                                     247             {
253                                                   248 
254                G4bool allOK = true;               249                G4bool allOK = true; 
255                for ( G4int i0i = 0 ; i0i < i0     250                for ( G4int i0i = 0 ; i0i < i0 ; i0i++ )
256                {                                  251                {
257                   if ( theMeanField->IsPauliBl    252                   if ( theMeanField->IsPauliBlocked ( i0i+n0 ) == true )
258                   {                               253                   {
259                      allOK = false;               254                      allOK = false;
260                      break;                       255                      break;
261                   }                               256                   } 
262                }                                  257                }
263                                                   258 
264                if ( allOK )                       259                if ( allOK ) 
265                {                                  260                {
266                   decayed = true; //Decay Succ    261                   decayed = true; //Decay Succeeded
267                }                                  262                }
268             }                                     263             }
269                                                   264 
270          }                                        265          }
271 //                                                266 //       
272                                                   267 
273          if ( decayed )                           268          if ( decayed )
274          {                                        269          {
275             //081119                              270             //081119
276             //G4cout << "Decay Suceeded! " <<     271             //G4cout << "Decay Suceeded! " << std::endl;
277             theSystem->GetParticipant( i )->Se    272             theSystem->GetParticipant( i )->SetHitMark();
278             for ( G4int i0i = 0 ; i0i < i0 ; i    273             for ( G4int i0i = 0 ; i0i < i0 ; i0i++ )
279             {                                     274             {
280                 theSystem->GetParticipant( i0i    275                 theSystem->GetParticipant( i0i+n0 )->SetHitMark();
281             }                                     276             }
282                                                   277 
283          }                                        278          }
284          else                                     279          else
285          {                                        280          {
286                                                   281 
287 //          Decay Blocked and re-enter orginal    282 //          Decay Blocked and re-enter orginal participant;
288                                                   283 
289             if ( isThisEnergyOK == true )  //     284             if ( isThisEnergyOK == true )  // for false case already done
290             {                                     285             {
291                                                   286 
292                theSystem->GetParticipant( i )-    287                theSystem->GetParticipant( i )->SetDefinition( pd0 );
293                theSystem->GetParticipant( i )-    288                theSystem->GetParticipant( i )->SetPosition( r0 );
294                theSystem->GetParticipant( i )-    289                theSystem->GetParticipant( i )->SetMomentum( p0 );
295                                                   290 
296                for ( G4int i0i = 0 ; i0i < i0     291                for ( G4int i0i = 0 ; i0i < i0 ; i0i++ )
297                {                                  292                {
298                   //081118                        293                   //081118
299                   //std::cout << "Decay Blocke    294                   //std::cout << "Decay Blocked deleteing " << i0i+n0  << std::endl;
300                   //160210 adding commnet: del << 295                   theSystem->DeleteParticipant( i0i+n0 );
301                   theSystem->DeleteParticipant << 
302                }                                  296                }
303                //081103                           297                //081103
304                theMeanField->Update();            298                theMeanField->Update();
305             }                                     299             }
306                                                   300 
307          }                                        301          }
308                                                   302 
309       }  //shortlive                              303       }  //shortlive
310    }  // go next participant                      304    }  // go next participant 
311 //071101                                          305 //071101
312                                                   306 
313                                                   307 
314    n = theSystem->GetTotalNumberOfParticipant(    308    n = theSystem->GetTotalNumberOfParticipant(); 
315                                                   309 
316 //081118                                          310 //081118
317    //for ( G4int i = 1 ; i < n ; i++ )            311    //for ( G4int i = 1 ; i < n ; i++ )
318    for ( G4int i = 1 ; i < theSystem->GetTotal    312    for ( G4int i = 1 ; i < theSystem->GetTotalNumberOfParticipant() ; i++ )
319    {                                              313    {
320                                                   314 
321       //std::cout << "Collision i " << i << st    315       //std::cout << "Collision i " << i << std::endl;
322       if ( theSystem->GetParticipant( i )->IsT    316       if ( theSystem->GetParticipant( i )->IsThisHit() ) continue; 
323                                                   317 
324       G4ThreeVector ri =  theSystem->GetPartic    318       G4ThreeVector ri =  theSystem->GetParticipant( i )->GetPosition();
325       G4LorentzVector p4i =  theSystem->GetPar    319       G4LorentzVector p4i =  theSystem->GetParticipant( i )->Get4Momentum();
326       G4double rmi =  theSystem->GetParticipan    320       G4double rmi =  theSystem->GetParticipant( i )->GetMass();
327       const G4ParticleDefinition* pdi =  theSy    321       const G4ParticleDefinition* pdi =  theSystem->GetParticipant( i )->GetDefinition();
328 //090331 gamma                                    322 //090331 gamma 
329       if ( pdi->GetPDGMass() == 0.0 ) continue    323       if ( pdi->GetPDGMass() == 0.0 ) continue;
330                                                   324 
331       //std::cout << " p4i00 " << p4i << std::    325       //std::cout << " p4i00 " << p4i << std::endl;
332       for ( G4int j = 0 ; j < i ; j++ )           326       for ( G4int j = 0 ; j < i ; j++ )
333       {                                           327       {
334                                                   328 
335                                                   329 
336 /*                                                330 /*
337          G4cout << "Collision " << i << " " <<    331          G4cout << "Collision " << i << " " << theSystem->GetParticipant( i )->IsThisProjectile() << G4endl;
338          G4cout << "Collision " << j << " " <<    332          G4cout << "Collision " << j << " " << theSystem->GetParticipant( j )->IsThisProjectile() << G4endl;
339          G4cout << "Collision " << i << " " <<    333          G4cout << "Collision " << i << " " << theSystem->GetParticipant( i )->IsThisTarget() << G4endl;
340          G4cout << "Collision " << j << " " <<    334          G4cout << "Collision " << j << " " << theSystem->GetParticipant( j )->IsThisTarget() << G4endl;
341 */                                                335 */
342                                                   336 
343          // Only 1 Collision allowed for each     337          // Only 1 Collision allowed for each particle in a time step. 
344          //081119                                 338          //081119
345          if ( theSystem->GetParticipant( i )->    339          if ( theSystem->GetParticipant( i )->IsThisHit() ) continue; 
346          if ( theSystem->GetParticipant( j )->    340          if ( theSystem->GetParticipant( j )->IsThisHit() ) continue; 
347                                                   341 
348          //std::cout << "Collision " << i << "    342          //std::cout << "Collision " << i << " " << j << std::endl;
349                                                   343 
350          // Do not allow collision between nuc    344          // Do not allow collision between nucleons in target/projectile til its first collision.
351          if ( theSystem->GetParticipant( i )->    345          if ( theSystem->GetParticipant( i )->IsThisProjectile() )
352          {                                        346          {
353             if ( theSystem->GetParticipant( j     347             if ( theSystem->GetParticipant( j )->IsThisProjectile() ) continue;
354          }                                        348          }
355          else if ( theSystem->GetParticipant(     349          else if ( theSystem->GetParticipant( i )->IsThisTarget() )
356          {                                        350          {
357             if ( theSystem->GetParticipant( j     351             if ( theSystem->GetParticipant( j )->IsThisTarget() ) continue;
358          }                                        352          }
359                                                   353 
360                                                   354 
361          G4ThreeVector rj =  theSystem->GetPar    355          G4ThreeVector rj =  theSystem->GetParticipant( j )->GetPosition();
362          G4LorentzVector p4j =  theSystem->Get    356          G4LorentzVector p4j =  theSystem->GetParticipant( j )->Get4Momentum();
363          G4double rmj =  theSystem->GetPartici    357          G4double rmj =  theSystem->GetParticipant( j )->GetMass();
364          const G4ParticleDefinition* pdj =  th    358          const G4ParticleDefinition* pdj =  theSystem->GetParticipant( j )->GetDefinition();
365 //090331 gamma                                    359 //090331 gamma 
366          if ( pdj->GetPDGMass() == 0.0 ) conti    360          if ( pdj->GetPDGMass() == 0.0 ) continue;
367                                                   361 
368          G4double rr2 = theMeanField->GetRR2(     362          G4double rr2 = theMeanField->GetRR2( i , j );
369                                                   363 
370 //       Here we assume elab (beam momentum le    364 //       Here we assume elab (beam momentum less than 5 GeV/n )
371          if ( rr2 > fdeltar*fdeltar ) continue    365          if ( rr2 > fdeltar*fdeltar ) continue;
372                                                   366 
373          //G4double s = (p4i+p4j)*(p4i+p4j);      367          //G4double s = (p4i+p4j)*(p4i+p4j);
374          //G4double srt = std::sqrt ( s );        368          //G4double srt = std::sqrt ( s );
375                                                   369 
376          G4double srt = std::sqrt( (p4i+p4j)*(    370          G4double srt = std::sqrt( (p4i+p4j)*(p4i+p4j) );
377                                                   371 
378          G4double cutoff = 0.0;                   372          G4double cutoff = 0.0;
379          G4double fbcmax = 0.0;                   373          G4double fbcmax = 0.0;
380          //110617 fix for gcc 4.6 compilation     374          //110617 fix for gcc 4.6 compilation warnings 
381          //G4double sig = 0.0;                    375          //G4double sig = 0.0;
382                                                   376 
383          if ( rmi < 0.94 && rmj < 0.94 )          377          if ( rmi < 0.94 && rmj < 0.94 ) 
384          {                                        378          {
385 //          nucleon or pion case                  379 //          nucleon or pion case
386             cutoff = rmi + rmj + 0.02;            380             cutoff = rmi + rmj + 0.02; 
387             fbcmax = fbcmax0;                     381             fbcmax = fbcmax0;
388             //110617 fix for gcc 4.6 compilati    382             //110617 fix for gcc 4.6 compilation warnings 
389             //sig = sig0;                         383             //sig = sig0;
390          }                                        384          }
391          else                                     385          else
392          {                                        386          {
393             cutoff = rmi + rmj;                   387             cutoff = rmi + rmj; 
394             fbcmax = fbcmax1;                     388             fbcmax = fbcmax1;
395             //110617 fix for gcc compilation w    389             //110617 fix for gcc compilation warnings 
396             //sig = sig1;                         390             //sig = sig1;
397          }                                        391          }
398                                                   392 
399          //std::cout << "Collision cutoff " <<    393          //std::cout << "Collision cutoff " << i << " " << j << " " << cutoff << std::endl;
400          if ( srt < cutoff ) continue;            394          if ( srt < cutoff ) continue; 
401                                                   395         
402          G4ThreeVector dr = ri - rj;              396          G4ThreeVector dr = ri - rj;
403          G4double rsq = dr*dr;                    397          G4double rsq = dr*dr;
404                                                   398 
405          G4double pij = p4i*p4j;                  399          G4double pij = p4i*p4j; 
406          G4double pidr = p4i.vect()*dr;           400          G4double pidr = p4i.vect()*dr;
407          G4double pjdr = p4j.vect()*dr;           401          G4double pjdr = p4j.vect()*dr;
408                                                   402 
409          G4double aij = 1.0 - ( rmi*rmj /pij )    403          G4double aij = 1.0 - ( rmi*rmj /pij ) * ( rmi*rmj /pij ); 
410          G4double bij = pidr / rmi - pjdr*rmi/    404          G4double bij = pidr / rmi - pjdr*rmi/pij;
411          G4double cij = rsq + ( pidr / rmi ) *    405          G4double cij = rsq + ( pidr / rmi ) * ( pidr / rmi );
412          G4double brel = std::sqrt ( std::abs     406          G4double brel = std::sqrt ( std::abs ( cij - bij*bij/aij ) );
413                                                   407  
414          if ( brel > fbcmax ) continue;           408          if ( brel > fbcmax ) continue;
415          //std::cout << "collisions3 " << std:    409          //std::cout << "collisions3 " << std::endl;
416                                                   410      
417          G4double bji = -pjdr/rmj + pidr * rmj    411          G4double bji = -pjdr/rmj + pidr * rmj /pij;
418                                                   412  
419          G4double ti = ( pidr/rmi - bij / aij     413          G4double ti = ( pidr/rmi - bij / aij ) * p4i.e() / rmi;
420          G4double tj = (-pjdr/rmj - bji / aij     414          G4double tj = (-pjdr/rmj - bji / aij ) * p4j.e() / rmj;
421                                                   415 
422                                                   416 
423 /*                                                417 /*
424          G4cout << "collisions4  p4i " << p4i     418          G4cout << "collisions4  p4i " << p4i << G4endl;
425          G4cout << "collisions4  ri " << ri <<    419          G4cout << "collisions4  ri " << ri << G4endl;
426          G4cout << "collisions4  p4j " << p4j     420          G4cout << "collisions4  p4j " << p4j << G4endl;
427          G4cout << "collisions4  rj " << rj <<    421          G4cout << "collisions4  rj " << rj << G4endl;
428          G4cout << "collisions4  dr " << dr <<    422          G4cout << "collisions4  dr " << dr << G4endl;
429          G4cout << "collisions4  pij " << pij     423          G4cout << "collisions4  pij " << pij << G4endl;
430          G4cout << "collisions4  aij " << aij     424          G4cout << "collisions4  aij " << aij << G4endl;
431          G4cout << "collisions4  bij bji " <<     425          G4cout << "collisions4  bij bji " << bij << " " << bji << G4endl;
432          G4cout << "collisions4  pidr pjdr " <    426          G4cout << "collisions4  pidr pjdr " << pidr << " " << pjdr << G4endl;
433          G4cout << "collisions4  p4i.e() p4j.e    427          G4cout << "collisions4  p4i.e() p4j.e() " << p4i.e() << " " << p4j.e() << G4endl;
434          G4cout << "collisions4  rmi rmj " <<     428          G4cout << "collisions4  rmi rmj " << rmi << " " << rmj << G4endl;
435          G4cout << "collisions4 " << ti << " "    429          G4cout << "collisions4 " << ti << " " << tj << G4endl;
436 */                                                430 */
437          if ( std::abs ( ti + tj ) > deltaT )     431          if ( std::abs ( ti + tj ) > deltaT ) continue;
438          //std::cout << "collisions4 " << std:    432          //std::cout << "collisions4 " << std::endl;
439                                                   433 
440          G4ThreeVector beta = ( p4i + p4j ).bo    434          G4ThreeVector beta = ( p4i + p4j ).boostVector();
441                                                   435 
442          G4LorentzVector p = p4i;                 436          G4LorentzVector p = p4i;
443          G4LorentzVector p4icm = p.boost( p.fi    437          G4LorentzVector p4icm = p.boost( p.findBoostToCM ( p4j ) );
444          G4ThreeVector pcm = p4icm.vect();        438          G4ThreeVector pcm = p4icm.vect();
445                                                   439          
446          G4double prcm = pcm.mag();               440          G4double prcm = pcm.mag();
447                                                   441 
448          if ( prcm <= 0.00001 ) continue;         442          if ( prcm <= 0.00001 ) continue; 
449          //std::cout << "collisions5 " << std:    443          //std::cout << "collisions5 " << std::endl;
450                                                   444 
451          G4bool energetically_forbidden = !( C    445          G4bool energetically_forbidden = !( CalFinalStateOfTheBinaryCollision ( i , j ) ); // Use Geant4 Collision Library
452          //G4bool energetically_forbidden = !(    446          //G4bool energetically_forbidden = !( CalFinalStateOfTheBinaryCollisionJQMD ( sig , cutoff , pcm , prcm , srt, beta , gamma , i , j ) ); // JQMD Elastic 
453                                                   447 
454 /*                                                448 /*
455          G4bool pauli_blocked = false;            449          G4bool pauli_blocked = false;
456          if ( energetically_forbidden == false    450          if ( energetically_forbidden == false ) // result true 
457          {                                        451          { 
458             if ( theMeanField->IsPauliBlocked     452             if ( theMeanField->IsPauliBlocked ( i ) == true || theMeanField->IsPauliBlocked ( j ) == true ) 
459             {                                     453             {
460                pauli_blocked = true;              454                pauli_blocked = true;
461                //std::cout << "G4QMDRESULT Col    455                //std::cout << "G4QMDRESULT Collsion Pauli Blocked " << std::endl;
462             }                                     456             }
463          }                                        457          }
464          else                                     458          else
465          {                                        459          {
466             if ( theMeanField->IsPauliBlocked     460             if ( theMeanField->IsPauliBlocked ( i ) == true || theMeanField->IsPauliBlocked ( j ) == true ) 
467                pauli_blocked = false;             461                pauli_blocked = false;
468             //std::cout << "G4QMDRESULT Collsi    462             //std::cout << "G4QMDRESULT Collsion Blocked " << std::endl;
469          }                                        463          } 
470 */                                                464 */
471                                                   465 
472 /*                                                466 /*
473             G4cout << "G4QMDRESULT Collsion in    467             G4cout << "G4QMDRESULT Collsion initial p4 i and j " 
474                       << p4i << " " << p4j        468                       << p4i << " " << p4j
475                       << G4endl;                  469                       << G4endl;
476 */                                                470 */
477 //       081118                                   471 //       081118
478          //if ( energetically_forbidden == tru    472          //if ( energetically_forbidden == true || pauli_blocked == true )
479          if ( energetically_forbidden == true     473          if ( energetically_forbidden == true )
480          {                                        474          {
481                                                   475 
482             //G4cout << " energetically_forbid    476             //G4cout << " energetically_forbidden  " << G4endl;
483 //          Collsion not allowed then re enter    477 //          Collsion not allowed then re enter orginal participants 
484 //          Now only momentum, becasuse we onl    478 //          Now only momentum, becasuse we only consider elastic scattering of nucleons
485                                                   479 
486             theSystem->GetParticipant( i )->Se    480             theSystem->GetParticipant( i )->SetMomentum( p4i.vect() );
487             theSystem->GetParticipant( i )->Se    481             theSystem->GetParticipant( i )->SetDefinition( pdi );
488             theSystem->GetParticipant( i )->Se    482             theSystem->GetParticipant( i )->SetPosition( ri );
489                                                   483 
490             theSystem->GetParticipant( j )->Se    484             theSystem->GetParticipant( j )->SetMomentum( p4j.vect() );
491             theSystem->GetParticipant( j )->Se    485             theSystem->GetParticipant( j )->SetDefinition( pdj );
492             theSystem->GetParticipant( j )->Se    486             theSystem->GetParticipant( j )->SetPosition( rj );
493                                                   487 
494             theMeanField->Cal2BodyQuantities(     488             theMeanField->Cal2BodyQuantities( i ); 
495             theMeanField->Cal2BodyQuantities(     489             theMeanField->Cal2BodyQuantities( j ); 
496                                                   490 
497          }                                        491          }
498          else                                     492          else 
499          {                                        493          {
500                                                   494 
501                                                   495             
502            G4bool absorption = false;             496            G4bool absorption = false; 
503            if ( n == theSystem->GetTotalNumber    497            if ( n == theSystem->GetTotalNumberOfParticipant()+1 ) absorption = true;
504            if ( absorption )                      498            if ( absorption ) 
505            {                                      499            {
506               //G4cout << "Absorption happend     500               //G4cout << "Absorption happend " << G4endl; 
507               i = i-1;                            501               i = i-1; 
508               n = n-1;                            502               n = n-1;
509            }                                      503            } 
510                                                   504               
511 //          Collsion allowed (really happened)    505 //          Collsion allowed (really happened) 
512                                                   506 
513             // Unset Projectile/Target flag       507             // Unset Projectile/Target flag
514             theSystem->GetParticipant( i )->Un    508             theSystem->GetParticipant( i )->UnsetInitialMark();
515             if ( !absorption ) theSystem->GetP    509             if ( !absorption ) theSystem->GetParticipant( j )->UnsetInitialMark();
516                                                   510 
517             theSystem->GetParticipant( i )->Se    511             theSystem->GetParticipant( i )->SetHitMark();
518             if ( !absorption ) theSystem->GetP    512             if ( !absorption ) theSystem->GetParticipant( j )->SetHitMark();
519                                                   513 
520             theSystem->IncrementCollisionCount    514             theSystem->IncrementCollisionCounter();
521                                                   515 
522 /*                                                516 /*
523             G4cout << "G4QMDRESULT Collsion Re    517             G4cout << "G4QMDRESULT Collsion Really Happened between " 
524                       << i << " and " << j        518                       << i << " and " << j 
525                       << G4endl;                  519                       << G4endl;
526             G4cout << "G4QMDRESULT Collsion in    520             G4cout << "G4QMDRESULT Collsion initial p4 i and j " 
527                       << p4i << " " << p4j        521                       << p4i << " " << p4j
528                       << G4endl;                  522                       << G4endl;
529             G4cout << "G4QMDRESULT Collsion af    523             G4cout << "G4QMDRESULT Collsion after p4 i and j " 
530                       << theSystem->GetPartici    524                       << theSystem->GetParticipant( i )->Get4Momentum()
531                       << " "                      525                       << " " 
532                       << theSystem->GetPartici    526                       << theSystem->GetParticipant( j )->Get4Momentum()
533                       << G4endl;                  527                       << G4endl;
534             G4cout << "G4QMDRESULT Collsion Di    528             G4cout << "G4QMDRESULT Collsion Diff " 
535                       << p4i + p4j - theSystem    529                       << p4i + p4j - theSystem->GetParticipant( i )->Get4Momentum() - theSystem->GetParticipant( j )->Get4Momentum()
536                       << G4endl;                  530                       << G4endl;
537             G4cout << "G4QMDRESULT Collsion in    531             G4cout << "G4QMDRESULT Collsion initial r i and j " 
538                       << ri << " " << rj          532                       << ri << " " << rj
539                       << G4endl;                  533                       << G4endl;
540             G4cout << "G4QMDRESULT Collsion af    534             G4cout << "G4QMDRESULT Collsion after r i and j " 
541                       << theSystem->GetPartici    535                       << theSystem->GetParticipant( i )->GetPosition()
542                       << " "                      536                       << " " 
543                       << theSystem->GetPartici    537                       << theSystem->GetParticipant( j )->GetPosition()
544                       << G4endl;                  538                       << G4endl;
545 */                                                539 */
546                                                   540              
547                                                   541 
548          }                                        542          }
549                                                   543 
550       }                                           544       }
551                                                   545 
552    }                                              546    }
553                                                   547 
554                                                   548 
555 }                                                 549 }
556                                                   550 
557                                                   551 
558                                                   552 
559 G4bool G4QMDCollision::CalFinalStateOfTheBinar    553 G4bool G4QMDCollision::CalFinalStateOfTheBinaryCollision( G4int i , G4int j )
560 {                                                 554 {
561                                                   555 
562 //081103                                          556 //081103
563    //G4cout << "CalFinalStateOfTheBinaryCollis    557    //G4cout << "CalFinalStateOfTheBinaryCollision " << i << " " << j << " " << theSystem->GetTotalNumberOfParticipant() << G4endl;
564                                                   558 
565    G4bool result = false;                         559    G4bool result = false;
566    G4bool energyOK = false;                       560    G4bool energyOK = false; 
567    G4bool pauliOK = false;                        561    G4bool pauliOK = false; 
568    G4bool abs = false;                            562    G4bool abs = false;
569    G4QMDParticipant* absorbed = NULL;             563    G4QMDParticipant* absorbed = NULL;  
570                                                   564 
571    G4LorentzVector p4i = theSystem->GetPartici    565    G4LorentzVector p4i = theSystem->GetParticipant( i )->Get4Momentum();
572    G4LorentzVector p4j = theSystem->GetPartici    566    G4LorentzVector p4j = theSystem->GetParticipant( j )->Get4Momentum();
573                                                   567 
574 //071031                                          568 //071031
575                                                   569 
576    G4double epot = theMeanField->GetTotalPoten    570    G4double epot = theMeanField->GetTotalPotential();
577                                                   571 
578    G4double eini = epot + p4i.e() + p4j.e();      572    G4double eini = epot + p4i.e() + p4j.e();
579                                                   573 
580 //071031                                          574 //071031
581    // will use KineticTrack                       575    // will use KineticTrack
582    const G4ParticleDefinition* pdi0 =theSystem    576    const G4ParticleDefinition* pdi0 =theSystem->GetParticipant( i )->GetDefinition();
583    const G4ParticleDefinition* pdj0 =theSystem    577    const G4ParticleDefinition* pdj0 =theSystem->GetParticipant( j )->GetDefinition();
584    G4LorentzVector p4i0 = p4i*GeV;                578    G4LorentzVector p4i0 = p4i*GeV;
585    G4LorentzVector p4j0 = p4j*GeV;                579    G4LorentzVector p4j0 = p4j*GeV;
586    G4ThreeVector ri0 = ( theSystem->GetPartici    580    G4ThreeVector ri0 = ( theSystem->GetParticipant( i )->GetPosition() )*fermi;
587    G4ThreeVector rj0 = ( theSystem->GetPartici    581    G4ThreeVector rj0 = ( theSystem->GetParticipant( j )->GetPosition() )*fermi;
588                                                   582 
589    for ( G4int iitry = 0 ; iitry < 4 ; iitry++    583    for ( G4int iitry = 0 ; iitry < 4 ; iitry++ )
590    {                                              584    {
591                                                   585 
592       abs = false;                                586       abs = false;
593                                                   587 
594       G4KineticTrack kt1( pdi0 , 0.0 , ri0 , p    588       G4KineticTrack kt1( pdi0 , 0.0 , ri0 , p4i0 );
595       G4KineticTrack kt2( pdj0 , 0.0 , rj0 , p    589       G4KineticTrack kt2( pdj0 , 0.0 , rj0 , p4j0 );
596                                                   590 
597       G4LorentzVector p4ix_new;                   591       G4LorentzVector p4ix_new; 
598       G4LorentzVector p4jx_new;                   592       G4LorentzVector p4jx_new; 
599       G4KineticTrackVector* secs = NULL;          593       G4KineticTrackVector* secs = NULL;
600       secs = theScatterer->Scatter( kt1 , kt2     594       secs = theScatterer->Scatter( kt1 , kt2 );
601                                                   595 
602       //std::cout << "G4QMDSCATTERER BEFORE "     596       //std::cout << "G4QMDSCATTERER BEFORE " << kt1.GetDefinition()->GetParticleName() << " " << kt1.Get4Momentum()/GeV << " " << kt1.GetPosition()/fermi << std::endl;
603       //std::cout << "G4QMDSCATTERER BEFORE "     597       //std::cout << "G4QMDSCATTERER BEFORE " << kt2.GetDefinition()->GetParticleName() << " " << kt2.Get4Momentum()/GeV << " " << kt2.GetPosition()/fermi << std::endl;
604       //std::cout << "THESCATTERER " << theSca    598       //std::cout << "THESCATTERER " << theScatterer->GetCrossSection ( kt1 , kt2 )/millibarn << " " << elastic << " " << sig << std::endl;
605                                                   599 
606                                                   600 
607       if ( secs )                                 601       if ( secs )
608       {                                           602       {
609          G4int iti = 0;                           603          G4int iti = 0;
610          if (  secs->size() == 2 )                604          if (  secs->size() == 2 )
611          {                                        605          {
612             for ( G4KineticTrackVector::iterat    606             for ( G4KineticTrackVector::iterator it 
613                 = secs->begin() ; it != secs->    607                 = secs->begin() ; it != secs->end() ; it++ )
614             {                                     608             {
615                if ( iti == 0 )                    609                if ( iti == 0 ) 
616                {                                  610                {
617                   theSystem->GetParticipant( i    611                   theSystem->GetParticipant( i )->SetDefinition( (*it)->GetDefinition() );
618                   p4ix_new = (*it)->Get4Moment    612                   p4ix_new = (*it)->Get4Momentum()/GeV;
619                   //std::cout << "THESCATTERER    613                   //std::cout << "THESCATTERER " << (*it)->GetDefinition()->GetParticleName() << std::endl;
620                   theSystem->GetParticipant( i    614                   theSystem->GetParticipant( i )->SetMomentum( p4ix_new.v() );
621                }                                  615                } 
622                if ( iti == 1 )                    616                if ( iti == 1 )
623                {                                  617                {
624                   theSystem->GetParticipant( j    618                   theSystem->GetParticipant( j )->SetDefinition( (*it)->GetDefinition() );
625                   p4jx_new = (*it)->Get4Moment    619                   p4jx_new = (*it)->Get4Momentum()/GeV;
626                   //std::cout << "THESCATTERER    620                   //std::cout << "THESCATTERER " << p4jx_new.e()-p4jx_new.m() << std::endl;
627                   theSystem->GetParticipant( j    621                   theSystem->GetParticipant( j )->SetMomentum( p4jx_new.v() );
628                }                                  622                }
629                //std::cout << "G4QMDSCATTERER     623                //std::cout << "G4QMDSCATTERER AFTER " << (*it)->GetDefinition()->GetParticleName() << " " << (*it)->Get4Momentum()/GeV << std::endl;
630                iti++;                             624                iti++;
631             }                                     625             }
632          }                                        626          }
633          else if ( secs->size() == 1 )            627          else if ( secs->size() == 1 )
634          {                                        628          {
635 //081118                                          629 //081118
636             abs = true;                           630             abs = true;
637             //G4cout << "G4QMDCollision pion a    631             //G4cout << "G4QMDCollision pion absrorption " << secs->front()->GetDefinition()->GetParticleName() << G4endl;
638             //secs->front()->Decay();             632             //secs->front()->Decay();
639             theSystem->GetParticipant( i )->Se    633             theSystem->GetParticipant( i )->SetDefinition( secs->front()->GetDefinition() );
640             p4ix_new = secs->front()->Get4Mome    634             p4ix_new = secs->front()->Get4Momentum()/GeV;
641             theSystem->GetParticipant( i )->Se    635             theSystem->GetParticipant( i )->SetMomentum( p4ix_new.v() );
642                                                   636 
643          }                                        637          } 
644                                                   638 
645 //081118                                          639 //081118
646          if ( secs->size() > 2 )                  640          if ( secs->size() > 2 ) 
647          {                                        641          {
648                                                   642 
649             G4cout << "G4QMDCollision secs siz    643             G4cout << "G4QMDCollision secs size > 2;  " << secs->size() << G4endl;
650                                                   644 
651             for ( G4KineticTrackVector::iterat    645             for ( G4KineticTrackVector::iterator it 
652                 = secs->begin() ; it != secs->    646                 = secs->begin() ; it != secs->end() ; it++ )
653             {                                     647             {
654                G4cout << "G4QMDSCATTERER AFTER    648                G4cout << "G4QMDSCATTERER AFTER " << (*it)->GetDefinition()->GetParticleName() << " " << (*it)->Get4Momentum()/GeV << G4endl;
655             }                                     649             }
656                                                   650 
657          }                                        651          }
658                                                   652 
659          // deleteing KineticTrack                653          // deleteing KineticTrack
660          for ( G4KineticTrackVector::iterator     654          for ( G4KineticTrackVector::iterator it 
661                = secs->begin() ; it != secs->e    655                = secs->begin() ; it != secs->end() ; it++ )
662          {                                        656          {  
663             delete *it;                           657             delete *it;
664          }                                        658          }
665                                                   659 
666          delete secs;                             660          delete secs;
667       }                                           661       }
668 //071031                                          662 //071031
669                                                   663 
670       if ( !abs )                                 664       if ( !abs )
671       {                                           665       { 
672          theMeanField->Cal2BodyQuantities( i )    666          theMeanField->Cal2BodyQuantities( i ); 
673          theMeanField->Cal2BodyQuantities( j )    667          theMeanField->Cal2BodyQuantities( j ); 
674       }                                           668       } 
675       else                                        669       else
676       {                                           670       {
677          absorbed = theSystem->EraseParticipan    671          absorbed = theSystem->EraseParticipant( j ); 
678          theMeanField->Update();                  672          theMeanField->Update();
679       }                                           673       }
680                                                   674 
681       epot = theMeanField->GetTotalPotential()    675       epot = theMeanField->GetTotalPotential();
682                                                   676 
683       G4double efin = epot + p4ix_new.e() + p4    677       G4double efin = epot + p4ix_new.e() + p4jx_new.e(); 
684                                                   678 
685       //std::cout << "Collision NEW epot " <<     679       //std::cout << "Collision NEW epot " << i << " " << j << " " << epot << " " << std::abs ( eini - efin ) - fepse << std::endl;
686                                                   680 
687 /*                                                681 /*
688       G4cout << "Collision efin " << i << " "     682       G4cout << "Collision efin " << i << " " << j << " " << efin << G4endl;
689       G4cout << "Collision " << i << " " << j     683       G4cout << "Collision " << i << " " << j << " " << std::abs ( eini - efin ) << " " << fepse << G4endl;
690       G4cout << "Collision " << std::abs ( ein    684       G4cout << "Collision " << std::abs ( eini - efin ) << " " << fepse << G4endl;
691 */                                                685 */
692                                                   686 
693 //071031                                          687 //071031
694       if ( std::abs ( eini - efin ) < fepse )     688       if ( std::abs ( eini - efin ) < fepse ) 
695       {                                           689       {
696          // Collison OK                           690          // Collison OK 
697          //std::cout << "collisions6" << std::    691          //std::cout << "collisions6" << std::endl;
698          //std::cout << "collisions before " <    692          //std::cout << "collisions before " << p4i << " " << p4j << std::endl;
699          //std::cout << "collisions after " <<    693          //std::cout << "collisions after " << theSystem->GetParticipant( i )->Get4Momentum() << " " << theSystem->GetParticipant( j )->Get4Momentum() << std::endl;
700          //std::cout << "collisions dif " << (    694          //std::cout << "collisions dif " << ( p4i + p4j ) - ( theSystem->GetParticipant( i )->Get4Momentum() + theSystem->GetParticipant( j )->Get4Momentum() ) << std::endl;
701          //std::cout << "collisions before " <    695          //std::cout << "collisions before " << ri0/fermi << " " << rj0/fermi << std::endl;
702          //std::cout << "collisions after " <<    696          //std::cout << "collisions after " << theSystem->GetParticipant( i )->GetPosition() << " " << theSystem->GetParticipant( j )->GetPosition() << std::endl;
703          energyOK = true;                         697          energyOK = true;
704          break;                                   698          break;
705       }                                           699       }
706       else                                        700       else
707       {                                           701       {
708          //G4cout << "Energy Not OK " << G4end    702          //G4cout << "Energy Not OK " << G4endl;
709          if ( abs )                               703          if ( abs )
710          {                                        704          {
711             //G4cout << "TKDB reinsert j " <<     705             //G4cout << "TKDB reinsert j " << G4endl;
712             theSystem->InsertParticipant( abso    706             theSystem->InsertParticipant( absorbed , j );   
713             theMeanField->Update();               707             theMeanField->Update();
714          }                                        708          }
715          // do not need reinsert in no absropt    709          // do not need reinsert in no absroption case 
716       }                                           710       }
717 //071031                                          711 //071031
718    }                                              712    }
719                                                   713 
720 // Energetically forbidden collision              714 // Energetically forbidden collision
721                                                   715 
722    if ( energyOK )                                716    if ( energyOK )
723    {                                              717    {
724       // Pauli Check                              718       // Pauli Check 
725       //G4cout << "Pauli Checking " << theSyst    719       //G4cout << "Pauli Checking " << theSystem->GetTotalNumberOfParticipant() << G4endl;
726       if ( !abs )                                 720       if ( !abs ) 
727       {                                           721       {
728          if ( !( theMeanField->IsPauliBlocked     722          if ( !( theMeanField->IsPauliBlocked ( i ) == true || theMeanField->IsPauliBlocked ( j ) == true ) ) 
729          {                                        723          {
730             //G4cout << "Binary Collision Happ    724             //G4cout << "Binary Collision Happen " << theSystem->GetTotalNumberOfParticipant() << G4endl;
731             pauliOK = true;                       725             pauliOK = true;
732          }                                        726          }
733       }                                           727       }
734       else                                        728       else 
735       {                                           729       {
736          //if ( theMeanField->IsPauliBlocked (    730          //if ( theMeanField->IsPauliBlocked ( i ) == false ) 
737          //090126                            i    731          //090126                            i-1 cause jth is erased
738          if ( theMeanField->IsPauliBlocked ( i    732          if ( theMeanField->IsPauliBlocked ( i-1 ) == false ) 
739          {                                        733          { 
740             //G4cout << "Absorption Happen " <    734             //G4cout << "Absorption Happen " << theSystem->GetTotalNumberOfParticipant() << G4endl;
741             delete absorbed;                      735             delete absorbed;
742             pauliOK = true;                       736             pauliOK = true;
743          }                                        737          }
744       }                                           738       }
745                                                   739       
746                                                   740 
747       if ( pauliOK )                              741       if ( pauliOK ) 
748       {                                           742       {
749          result = true;                           743          result = true;
750       }                                           744       }
751       else                                        745       else
752       {                                           746       {
753          //G4cout << "Pauli Blocked" << G4endl    747          //G4cout << "Pauli Blocked" << G4endl;
754          if ( abs )                               748          if ( abs )
755          {                                        749          {
756             //G4cout << "TKDB reinsert j pauli    750             //G4cout << "TKDB reinsert j pauli block" << G4endl;
757             theSystem->InsertParticipant( abso    751             theSystem->InsertParticipant( absorbed , j );   
758             theMeanField->Update();               752             theMeanField->Update(); 
759          }                                        753          }
760       }                                           754       }
761    }                                              755    }
762                                                   756 
763    return result;                                 757    return result;
764                                                   758 
765 }                                                 759 } 
766                                                   760 
767                                                   761 
768                                                   762 
769 G4bool G4QMDCollision::CalFinalStateOfTheBinar    763 G4bool G4QMDCollision::CalFinalStateOfTheBinaryCollisionJQMD( G4double sig , G4double cutoff , G4ThreeVector pcm , G4double prcm , G4double srt , G4ThreeVector beta , G4double gamma , G4int i , G4int j )
770 {                                                 764 {
771                                                   765 
772    //G4cout << "CalFinalStateOfTheBinaryCollis    766    //G4cout << "CalFinalStateOfTheBinaryCollisionJQMD" << G4endl;
773                                                   767 
774    G4bool result = true;                          768    G4bool result = true;
775                                                   769 
776    G4LorentzVector p4i =  theSystem->GetPartic    770    G4LorentzVector p4i =  theSystem->GetParticipant( i )->Get4Momentum();
777    G4double rmi =  theSystem->GetParticipant(     771    G4double rmi =  theSystem->GetParticipant( i )->GetMass();
778    G4int zi =  theSystem->GetParticipant( i )-    772    G4int zi =  theSystem->GetParticipant( i )->GetChargeInUnitOfEplus();
779                                                   773 
780    G4LorentzVector p4j =  theSystem->GetPartic    774    G4LorentzVector p4j =  theSystem->GetParticipant( j )->Get4Momentum();
781    G4double rmj =  theSystem->GetParticipant(     775    G4double rmj =  theSystem->GetParticipant( j )->GetMass();
782    G4int zj =  theSystem->GetParticipant( j )-    776    G4int zj =  theSystem->GetParticipant( j )->GetChargeInUnitOfEplus();
783                                                   777 
784    G4double pr = prcm;                            778    G4double pr = prcm;
785                                                   779 
786    G4double c2  = pcm.z()/pr;                     780    G4double c2  = pcm.z()/pr;
787                                                   781     
788    G4double csrt = srt - cutoff;                  782    G4double csrt = srt - cutoff;
789                                                   783 
790    //G4double pri = prcm;                         784    //G4double pri = prcm;
791    //G4double prf = sqrt ( 0.25 * srt*srt -rm2    785    //G4double prf = sqrt ( 0.25 * srt*srt -rm2 );
792                                                   786     
793    G4double asrt = srt - rmi - rmj;               787    G4double asrt = srt - rmi - rmj;
794    G4double pra = prcm;                           788    G4double pra = prcm;
795                                                   789    
796                                                   790 
797                                                   791 
798    G4double elastic = 0.0;                        792    G4double elastic = 0.0; 
799                                                   793 
800    if ( zi == zj )                                794    if ( zi == zj )
801    {                                              795    {
802       if ( csrt < 0.4286 )                        796       if ( csrt < 0.4286 )
803       {                                           797       {
804          elastic = 35.0 / ( 1. + csrt * 100.0     798          elastic = 35.0 / ( 1. + csrt * 100.0 )  +  20.0;
805       }                                           799       }
806       else                                        800       else
807       {                                           801       {
808          elastic = ( - std::atan( ( csrt - 0.4    802          elastic = ( - std::atan( ( csrt - 0.4286 ) * 1.5 - 0.8 )
809                  *   2. / pi + 1.0 ) * 9.65 +     803                  *   2. / pi + 1.0 ) * 9.65 + 7.0;
810       }                                           804       }         
811    }                                              805    }
812    else                                           806    else
813    {                                              807    {
814       if ( csrt < 0.4286 )                        808       if ( csrt < 0.4286 )
815       {                                           809       {
816          elastic = 28.0 / ( 1. + csrt * 100.0     810          elastic = 28.0 / ( 1. + csrt * 100.0 )  +  27.0;
817       }                                           811       }
818       else                                        812       else
819       {                                           813       {
820          elastic = ( - std::atan( ( csrt - 0.4    814          elastic = ( - std::atan( ( csrt - 0.4286 ) * 1.5 - 0.8 )
821                  *   2. / pi + 1.0 ) * 12.34 +    815                  *   2. / pi + 1.0 ) * 12.34 + 10.0;
822       }                                           816       }         
823    }                                              817    }
824                                                   818 
825 //   std::cout << "Collision csrt " << i << "     819 //   std::cout << "Collision csrt " << i << " " << j << " " << csrt << std::endl;
826 //   std::cout << "Collision elstic " << i <<     820 //   std::cout << "Collision elstic " << i << " " << j << " " << elastic << std::endl;
827                                                   821 
828                                                   822 
829 //   std::cout << "Collision sig " << i << " "    823 //   std::cout << "Collision sig " << i << " " << j  << " " << sig << std::endl;
830    if ( G4UniformRand() > elastic / sig )         824    if ( G4UniformRand() > elastic / sig ) 
831    {                                              825    { 
832       //std::cout << "Inelastic " << std::endl    826       //std::cout << "Inelastic " << std::endl; 
833       //std::cout << "elastic/sig " << elastic    827       //std::cout << "elastic/sig " << elastic/sig << std::endl; 
834       return result;                              828       return result; 
835    }                                              829    }
836    else                                           830    else
837    {                                              831    {
838       //std::cout << "Elastic " << std::endl;     832       //std::cout << "Elastic " << std::endl; 
839    }                                              833    } 
840 //   std::cout << "Collision ELSTIC " << i <<     834 //   std::cout << "Collision ELSTIC " << i << " " << j << std::endl;
841                                                   835 
842                                                   836    
843    G4double as = G4Pow::GetInstance()->powN (  << 837    G4double as = std::pow ( 3.65 * asrt , 6 );
844    G4double a = 6.0 * as / (1.0 + as);            838    G4double a = 6.0 * as / (1.0 + as);
845    G4double ta = -2.0 * pra*pra;                  839    G4double ta = -2.0 * pra*pra;
846    G4double x = G4UniformRand();                  840    G4double x = G4UniformRand(); 
847    G4double t1 = G4Log( (1-x) * G4Exp(2.*a*ta) << 841    G4double t1 = std::log( (1-x) * std::exp(2.*a*ta) + x )  /  a;
848    G4double c1 = 1.0 - t1/ta;                     842    G4double c1 = 1.0 - t1/ta;
849                                                   843  
850    if( std::abs(c1) > 1.0 ) c1 = 2.0 * x - 1.0    844    if( std::abs(c1) > 1.0 ) c1 = 2.0 * x - 1.0;
851                                                   845 
852 /*                                                846 /*
853    G4cout << "Collision as " << i << " " << j     847    G4cout << "Collision as " << i << " " << j << " " << as << G4endl;
854    G4cout << "Collision a " << i << " " << j <    848    G4cout << "Collision a " << i << " " << j << " " << a << G4endl;
855    G4cout << "Collision ta " << i << " " << j     849    G4cout << "Collision ta " << i << " " << j << " " << ta << G4endl;
856    G4cout << "Collision x " << i << " " << j <    850    G4cout << "Collision x " << i << " " << j << " " << x << G4endl;
857    G4cout << "Collision t1 " << i << " " << j     851    G4cout << "Collision t1 " << i << " " << j << " " << t1 << G4endl;
858    G4cout << "Collision c1 " << i << " " << j     852    G4cout << "Collision c1 " << i << " " << j << " " << c1 << G4endl;
859 */                                                853 */
860    t1 = 2.0*pi*G4UniformRand();                   854    t1 = 2.0*pi*G4UniformRand(); 
861 //   std::cout << "Collision t1 " << i << " "     855 //   std::cout << "Collision t1 " << i << " " << j << " " << t1 << std::endl;
862    G4double t2 = 0.0;                             856    G4double t2 = 0.0;
863    if ( pcm.x() == 0.0 && pcm.y() == 0 )          857    if ( pcm.x() == 0.0 && pcm.y() == 0 )
864    {                                              858    {
865       t2 = 0.0;                                   859       t2 = 0.0;
866    }                                              860    }
867    else                                           861    else 
868    {                                              862    {
869       t2 = std::atan2( pcm.y() , pcm.x() );       863       t2 = std::atan2( pcm.y() , pcm.x() );
870    }                                              864    }
871 //      std::cout << "Collision t2 " << i << "    865 //      std::cout << "Collision t2 " << i << " " << j << " " << t2 << std::endl;
872                                                   866 
873    G4double s1 = std::sqrt ( 1.0 - c1*c1 );       867    G4double s1 = std::sqrt ( 1.0 - c1*c1 );
874    G4double s2 = std::sqrt ( 1.0 - c2*c2 );       868    G4double s2 = std::sqrt ( 1.0 - c2*c2 );
875                                                   869  
876    G4double ct1 = std::cos(t1);                   870    G4double ct1 = std::cos(t1);      
877    G4double st1 = std::sin(t1);                   871    G4double st1 = std::sin(t1);      
878                                                   872 
879    G4double ct2 = std::cos(t2);                   873    G4double ct2 = std::cos(t2);      
880    G4double st2 = std::sin(t2);                   874    G4double st2 = std::sin(t2);      
881                                                   875 
882    G4double ss = c2*s1*ct1 + s2*c1;               876    G4double ss = c2*s1*ct1 + s2*c1;
883                                                   877 
884    pcm.setX( pr * ( ss*ct2 - s1*st1*st2) );       878    pcm.setX( pr * ( ss*ct2 - s1*st1*st2) );
885    pcm.setY( pr * ( ss*st2 + s1*st1*ct2) );       879    pcm.setY( pr * ( ss*st2 + s1*st1*ct2) );
886    pcm.setZ( pr * ( c1*c2 - s1*s2*ct1) );         880    pcm.setZ( pr * ( c1*c2 - s1*s2*ct1) );
887                                                   881 
888 // std::cout << "Collision pcm " << i << " " <    882 // std::cout << "Collision pcm " << i << " " << j << " " << pcm << std::endl;
889                                                   883 
890    G4double epot = theMeanField->GetTotalPoten    884    G4double epot = theMeanField->GetTotalPotential();
891                                                   885 
892    G4double eini = epot + p4i.e() + p4j.e();      886    G4double eini = epot + p4i.e() + p4j.e();
893    G4double etwo = p4i.e() + p4j.e();             887    G4double etwo = p4i.e() + p4j.e();
894                                                   888 
895 /*                                                889 /*
896    G4cout << "Collision epot " << i << " " <<     890    G4cout << "Collision epot " << i << " " << j << " " << epot << G4endl;
897    G4cout << "Collision eini " << i << " " <<     891    G4cout << "Collision eini " << i << " " << j << " " << eini << G4endl;
898    G4cout << "Collision etwo " << i << " " <<     892    G4cout << "Collision etwo " << i << " " << j << " " << etwo << G4endl;
899 */                                                893 */
900                                                   894 
901                                                   895 
902    for ( G4int itry = 0 ; itry < 4 ; itry++ )     896    for ( G4int itry = 0 ; itry < 4 ; itry++ )
903    {                                              897    {
904                                                   898 
905       G4double eicm = std::sqrt ( rmi*rmi + pc    899       G4double eicm = std::sqrt ( rmi*rmi + pcm*pcm );
906       G4double pibeta = pcm*beta;                 900       G4double pibeta = pcm*beta;
907                                                   901 
908       G4double trans = gamma * ( gamma * pibet    902       G4double trans = gamma * ( gamma * pibeta / ( gamma + 1 ) + eicm );
909                                                   903    
910       G4ThreeVector pi_new = beta*trans + pcm;    904       G4ThreeVector pi_new = beta*trans + pcm;
911                                                   905    
912       G4double ejcm = std::sqrt ( rmj*rmj + pc    906       G4double ejcm = std::sqrt ( rmj*rmj + pcm*pcm );
913       trans = gamma * ( gamma * pibeta / ( gam    907       trans = gamma * ( gamma * pibeta / ( gamma + 1 ) + ejcm );
914                                                   908 
915       G4ThreeVector pj_new = beta*trans - pcm;    909       G4ThreeVector pj_new = beta*trans - pcm;
916                                                   910 
917 //                                                911 //
918 // Delete old                                     912 // Delete old 
919 // Add new Particitipants                         913 // Add new Particitipants
920 //                                                914 //
921 // Now only change momentum ( Beacuse we only     915 // Now only change momentum ( Beacuse we only have elastic sctter of nucleon
922 // In future Definition also will be change       916 // In future Definition also will be change 
923 //                                                917 //
924                                                   918 
925       theSystem->GetParticipant( i )->SetMomen    919       theSystem->GetParticipant( i )->SetMomentum( pi_new );
926       theSystem->GetParticipant( j )->SetMomen    920       theSystem->GetParticipant( j )->SetMomentum( pj_new );
927                                                   921 
928       G4double pi_new_e = (theSystem->GetParti    922       G4double pi_new_e = (theSystem->GetParticipant( i )->Get4Momentum()).e();
929       G4double pj_new_e = (theSystem->GetParti    923       G4double pj_new_e = (theSystem->GetParticipant( j )->Get4Momentum()).e();
930                                                   924 
931       theMeanField->Cal2BodyQuantities( i );      925       theMeanField->Cal2BodyQuantities( i ); 
932       theMeanField->Cal2BodyQuantities( j );      926       theMeanField->Cal2BodyQuantities( j ); 
933                                                   927 
934       epot = theMeanField->GetTotalPotential()    928       epot = theMeanField->GetTotalPotential();
935                                                   929 
936       G4double efin = epot + pi_new_e + pj_new    930       G4double efin = epot + pi_new_e + pj_new_e ; 
937                                                   931 
938       //std::cout << "Collision NEW epot " <<     932       //std::cout << "Collision NEW epot " << i << " " << j << " " << epot << " " << std::abs ( eini - efin ) - fepse << std::endl;
939 /*                                                933 /*
940       G4cout << "Collision efin " << i << " "     934       G4cout << "Collision efin " << i << " " << j << " " << efin << G4endl;
941       G4cout << "Collision " << i << " " << j     935       G4cout << "Collision " << i << " " << j << " " << std::abs ( eini - efin ) << " " << fepse << G4endl;
942       G4cout << "Collision " << std::abs ( ein    936       G4cout << "Collision " << std::abs ( eini - efin ) << " " << fepse << G4endl;
943 */                                                937 */
944                                                   938 
945 //071031                                          939 //071031
946       if ( std::abs ( eini - efin ) < fepse )     940       if ( std::abs ( eini - efin ) < fepse ) 
947       {                                           941       {
948    // Collison OK                                 942    // Collison OK 
949          //std::cout << "collisions6" << std::    943          //std::cout << "collisions6" << std::endl;
950          //std::cout << "collisions before " <    944          //std::cout << "collisions before " << p4i << " " << p4j << std::endl;
951          //std::cout << "collisions after " <<    945          //std::cout << "collisions after " << theSystem->GetParticipant( i )->Get4Momentum() << " " << theSystem->GetParticipant( j )->Get4Momentum() << std::endl;
952          //std::cout << "collisions dif " << (    946          //std::cout << "collisions dif " << ( p4i + p4j ) - ( theSystem->GetParticipant( i )->Get4Momentum() + theSystem->GetParticipant( j )->Get4Momentum() ) << std::endl;
953          //std::cout << "collisions before " <    947          //std::cout << "collisions before " << rix/fermi << " " << rjx/fermi << std::endl;
954          //std::cout << "collisions after " <<    948          //std::cout << "collisions after " << theSystem->GetParticipant( i )->GetPosition() << " " << theSystem->GetParticipant( j )->GetPosition() << std::endl;
955       }                                           949       }
956 //071031                                          950 //071031
957                                                   951 
958          if ( std::abs ( eini - efin ) < fepse    952          if ( std::abs ( eini - efin ) < fepse ) return result;  // Collison OK 
959                                                   953       
960          G4double cona = ( eini - efin + etwo     954          G4double cona = ( eini - efin + etwo ) / gamma;
961          G4double fac2 = 1.0 / ( 4.0 * cona*co    955          G4double fac2 = 1.0 / ( 4.0 * cona*cona * pr*pr ) *
962                        ( ( cona*cona - ( rmi*r    956                        ( ( cona*cona - ( rmi*rmi + rmj*rmj ) )*( cona*cona - ( rmi*rmi + rmj*rmj ) )
963                        - 4.0 * rmi*rmi * rmj*r    957                        - 4.0 * rmi*rmi * rmj*rmj );
964                                                   958 
965          if ( fac2 > 0 )                          959          if ( fac2 > 0 ) 
966          {                                        960          {
967             G4double fact = std::sqrt ( fac2 )    961             G4double fact = std::sqrt ( fac2 ); 
968             pcm = fact*pcm;                       962             pcm = fact*pcm;
969          }                                        963          }
970                                                   964 
971                                                   965 
972    }                                              966    }
973                                                   967 
974 // Energetically forbidden collision              968 // Energetically forbidden collision
975    result = false;                                969    result = false;
976                                                   970 
977    return result;                                 971    return result;
978                                                   972 
979 }                                                 973 } 
980                                                   974