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 9.3.p1)


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