Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/qmd/src/G4QMDReaction.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/G4QMDReaction.cc (Version 11.3.0) and /processes/hadronic/models/qmd/src/G4QMDReaction.cc (Version 10.0.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 // 080505 Fixed and changed sampling method of     26 // 080505 Fixed and changed sampling method of impact parameter by T. Koi 
 27 // 080602 Fix memory leaks by T. Koi               27 // 080602 Fix memory leaks by T. Koi 
 28 // 080612 Delete unnecessary dependency and un     28 // 080612 Delete unnecessary dependency and unused functions
 29 //        Change criterion of reaction by T. K     29 //        Change criterion of reaction by T. Koi
 30 // 081107 Add UnUseGEM (then use the default c     30 // 081107 Add UnUseGEM (then use the default channel of G4Evaporation)
 31 //            UseFrag (chage criterion of a in     31 //            UseFrag (chage criterion of a inelastic reaction)
 32 //        Fix bug in nucleon projectiles  by T     32 //        Fix bug in nucleon projectiles  by T. Koi    
 33 // 090122 Be8 -> Alpha + Alpha                     33 // 090122 Be8 -> Alpha + Alpha 
 34 // 090331 Change member shenXS and genspaXS ob     34 // 090331 Change member shenXS and genspaXS object to pointer 
 35 // 091119 Fix for incidence of neutral particl     35 // 091119 Fix for incidence of neutral particles 
 36 //                                                 36 //
 37 #include "G4QMDReaction.hh"                        37 #include "G4QMDReaction.hh"
 38 #include "G4QMDNucleus.hh"                         38 #include "G4QMDNucleus.hh"
 39 #include "G4QMDGroundStateNucleus.hh"              39 #include "G4QMDGroundStateNucleus.hh"
 40 #include "G4Pow.hh"                            <<  40 
 41 #include "G4PhysicalConstants.hh"                  41 #include "G4PhysicalConstants.hh"
 42 #include "G4SystemOfUnits.hh"                      42 #include "G4SystemOfUnits.hh"
 43 #include "G4NistManager.hh"                        43 #include "G4NistManager.hh"
 44                                                    44 
 45 #include "G4CrossSectionDataSetRegistry.hh"    << 
 46 #include "G4BGGPionElasticXS.hh"               << 
 47 #include "G4BGGPionInelasticXS.hh"             << 
 48 #include "G4VCrossSectionDataSet.hh"           << 
 49 #include "G4CrossSectionInelastic.hh"          << 
 50 #include "G4ComponentGGNuclNuclXsc.hh"         << 
 51 #include "G4PhysicsModelCatalog.hh"            << 
 52                                                << 
 53                                                << 
 54 G4QMDReaction::G4QMDReaction()                     45 G4QMDReaction::G4QMDReaction()
 55 : G4HadronicInteraction("QMDModel")                46 : G4HadronicInteraction("QMDModel")
 56 , system ( NULL )                                  47 , system ( NULL )
 57 , deltaT ( 1 ) // in fsec (c=1)                    48 , deltaT ( 1 ) // in fsec (c=1)
 58 , maxTime ( 100 ) // will have maxTime-th time     49 , maxTime ( 100 ) // will have maxTime-th time step
 59 , envelopF ( 1.05 ) // 10% for Peripheral reac     50 , envelopF ( 1.05 ) // 10% for Peripheral reactions
 60 , gem ( true )                                     51 , gem ( true )
 61 , frag ( false )                                   52 , frag ( false )
 62 , secID( -1 )                                  << 
 63 {                                                  53 {
 64    theXS = new G4CrossSectionInelastic( new G4 << 
 65    pipElNucXS = new G4BGGPionElasticXS(G4PionP << 
 66    pipElNucXS->BuildPhysicsTable(*(G4PionPlus: << 
 67                                                << 
 68    pimElNucXS = new G4BGGPionElasticXS(G4PionM << 
 69    pimElNucXS->BuildPhysicsTable(*(G4PionMinus << 
 70                                                << 
 71    pipInelNucXS = new G4BGGPionInelasticXS(G4P << 
 72    pipInelNucXS->BuildPhysicsTable(*(G4PionPlu << 
 73                                                << 
 74    pimInelNucXS = new G4BGGPionInelasticXS(G4P << 
 75    pimInelNucXS->BuildPhysicsTable(*(G4PionMin << 
 76                                                    54 
                                                   >>  55    //090331
                                                   >>  56    shenXS = new G4IonsShenCrossSection();
                                                   >>  57    //genspaXS = new G4GeneralSpaceNNCrossSection();
                                                   >>  58    piNucXS = new G4PiNuclearCrossSection();
 77    meanField = new G4QMDMeanField();               59    meanField = new G4QMDMeanField();
 78    collision = new G4QMDCollision();               60    collision = new G4QMDCollision();
 79                                                    61 
 80    excitationHandler = new G4ExcitationHandler <<  62    excitationHandler = new G4ExcitationHandler;
                                                   >>  63    evaporation = new G4Evaporation;
                                                   >>  64    excitationHandler->SetEvaporation( evaporation );
 81    setEvaporationCh();                             65    setEvaporationCh();
 82                                                << 
 83    coulomb_collision_gamma_proj = 0.0;         << 
 84    coulomb_collision_rx_proj = 0.0;            << 
 85    coulomb_collision_rz_proj = 0.0;            << 
 86    coulomb_collision_px_proj = 0.0;            << 
 87    coulomb_collision_pz_proj = 0.0;            << 
 88                                                << 
 89    coulomb_collision_gamma_targ = 0.0;         << 
 90    coulomb_collision_rx_targ = 0.0;            << 
 91    coulomb_collision_rz_targ = 0.0;            << 
 92    coulomb_collision_px_targ = 0.0;            << 
 93    coulomb_collision_pz_targ = 0.0;            << 
 94                                                << 
 95    secID = G4PhysicsModelCatalog::GetModelID(  << 
 96 }                                                  66 }
 97                                                    67 
 98                                                    68 
                                                   >>  69 
 99 G4QMDReaction::~G4QMDReaction()                    70 G4QMDReaction::~G4QMDReaction()
100 {                                                  71 {
                                                   >>  72    delete evaporation; 
101    delete excitationHandler;                       73    delete excitationHandler;
102    delete collision;                               74    delete collision;
103    delete meanField;                               75    delete meanField;
104 }                                                  76 }
105                                                    77 
106                                                    78 
                                                   >>  79 
107 G4HadFinalState* G4QMDReaction::ApplyYourself(     80 G4HadFinalState* G4QMDReaction::ApplyYourself( const G4HadProjectile & projectile , G4Nucleus & target )
108 {                                                  81 {
109    //G4cout << "G4QMDReaction::ApplyYourself"      82    //G4cout << "G4QMDReaction::ApplyYourself" << G4endl;
110                                                    83 
111    theParticleChange.Clear();                      84    theParticleChange.Clear();
112                                                    85 
113    system = new G4QMDSystem;                       86    system = new G4QMDSystem;
114                                                    87 
115    G4int proj_Z = 0;                               88    G4int proj_Z = 0;
116    G4int proj_A = 0;                               89    G4int proj_A = 0;
117    const G4ParticleDefinition* proj_pd = ( con <<  90    G4ParticleDefinition* proj_pd = ( G4ParticleDefinition* ) projectile.GetDefinition();
118    if ( proj_pd->GetParticleType() == "nucleus     91    if ( proj_pd->GetParticleType() == "nucleus" )
119    {                                               92    {
120       proj_Z = proj_pd->GetAtomicNumber();         93       proj_Z = proj_pd->GetAtomicNumber();
121       proj_A = proj_pd->GetAtomicMass();           94       proj_A = proj_pd->GetAtomicMass();
122    }                                               95    }
123    else                                            96    else
124    {                                               97    {
125       proj_Z = (int)( proj_pd->GetPDGCharge()/     98       proj_Z = (int)( proj_pd->GetPDGCharge()/eplus );
126       proj_A = 1;                                  99       proj_A = 1;
127    }                                              100    }
128    //G4int targ_Z = int ( target.GetZ() + 0.5     101    //G4int targ_Z = int ( target.GetZ() + 0.5 );
129    //G4int targ_A = int ( target.GetN() + 0.5     102    //G4int targ_A = int ( target.GetN() + 0.5 );
130    //migrate to integer A and Z (GetN_asInt re    103    //migrate to integer A and Z (GetN_asInt returns number of neutrons in the nucleus since this) 
131    G4int targ_Z = target.GetZ_asInt();            104    G4int targ_Z = target.GetZ_asInt();
132    G4int targ_A = target.GetA_asInt();            105    G4int targ_A = target.GetA_asInt();
133    const G4ParticleDefinition* targ_pd = G4Ion << 106    G4ParticleDefinition* targ_pd = G4IonTable::GetIonTable()->GetIon( targ_Z , targ_A , 0.0 );
134                                                   107 
135                                                   108 
136    //G4NistManager* nistMan = G4NistManager::I    109    //G4NistManager* nistMan = G4NistManager::Instance();
137 //   G4Element* G4NistManager::FindOrBuildElem    110 //   G4Element* G4NistManager::FindOrBuildElement( targ_Z );
138                                                   111 
139    const G4DynamicParticle* proj_dp = new G4Dy    112    const G4DynamicParticle* proj_dp = new G4DynamicParticle ( proj_pd , projectile.Get4Momentum() );
140    //const G4Element* targ_ele =  nistMan->Fin    113    //const G4Element* targ_ele =  nistMan->FindOrBuildElement( targ_Z ); 
141    //G4double aTemp = projectile.GetMaterial()    114    //G4double aTemp = projectile.GetMaterial()->GetTemperature();
142                                                   115 
143    // Glauber-Gribov nucleus-nucleus cross sec << 116      //090331
144    // therefore call GetElementCrossSection in << 117   
145    //G4double xs_0 = theXS->GetIsoCrossSection << 118    G4VCrossSectionDataSet* theXS = shenXS;
146    G4double xs_0 = theXS->GetElementCrossSecti << 119 
147                                                << 120    if ( proj_pd->GetParticleType() == "meson" ) theXS = piNucXS;
148    // When the projectile is a pion            << 
149    if (proj_pd == G4PionPlus::PionPlus() ) {   << 
150      xs_0 = pipElNucXS->GetElementCrossSection << 
151             pipInelNucXS->GetElementCrossSecti << 
152    } else if (proj_pd == G4PionMinus::PionMinu << 
153      xs_0 = pimElNucXS->GetElementCrossSection << 
154             pimInelNucXS->GetElementCrossSecti << 
155    }                                           << 
156                                                   121 
157    //G4double xs_0 = genspaXS->GetCrossSection    122    //G4double xs_0 = genspaXS->GetCrossSection ( proj_dp , targ_ele , aTemp );
158    //G4double xs_0 = theXS->GetCrossSection (     123    //G4double xs_0 = theXS->GetCrossSection ( proj_dp , targ_ele , aTemp );
159    //110822                                       124    //110822 
                                                   >> 125    G4double xs_0 = theXS->GetIsoCrossSection ( proj_dp , targ_Z , targ_A );
160                                                   126 
161      G4double bmax_0 = std::sqrt( xs_0 / pi );    127      G4double bmax_0 = std::sqrt( xs_0 / pi );
162      //std::cout << "bmax_0 in fm (fermi) " <<    128      //std::cout << "bmax_0 in fm (fermi) " <<  bmax_0/fermi << std::endl;
163                                                   129 
164      //delete proj_dp;                            130      //delete proj_dp; 
165                                                   131 
166    G4bool elastic = true;                         132    G4bool elastic = true;
167                                                   133    
168    std::vector< G4QMDNucleus* > nucleuses; //     134    std::vector< G4QMDNucleus* > nucleuses; // Secondary nuceluses 
169    G4ThreeVector boostToReac; // ReactionSyste    135    G4ThreeVector boostToReac; // ReactionSystem (CM or NN); 
170    G4ThreeVector boostBackToLAB; // Reaction S    136    G4ThreeVector boostBackToLAB; // Reaction System to LAB; 
171                                                   137 
172    G4LorentzVector targ4p( G4ThreeVector( 0.0     138    G4LorentzVector targ4p( G4ThreeVector( 0.0 ) , targ_pd->GetPDGMass()/GeV );
173    G4ThreeVector boostLABtoCM = targ4p.findBoo    139    G4ThreeVector boostLABtoCM = targ4p.findBoostToCM( proj_dp->Get4Momentum()/GeV ); // CM of target and proj; 
174                                                   140 
175    G4double p1 = proj_dp->GetMomentum().mag()/    141    G4double p1 = proj_dp->GetMomentum().mag()/GeV/proj_A; 
176    G4double m1 = proj_dp->GetDefinition()->Get    142    G4double m1 = proj_dp->GetDefinition()->GetPDGMass()/GeV/proj_A;
177    G4double e1 = std::sqrt( p1*p1 + m1*m1 );      143    G4double e1 = std::sqrt( p1*p1 + m1*m1 ); 
178    G4double e2 = targ_pd->GetPDGMass()/GeV/tar    144    G4double e2 = targ_pd->GetPDGMass()/GeV/targ_A;
179    G4double beta_nn = -p1 / ( e1+e2 );            145    G4double beta_nn = -p1 / ( e1+e2 );
180                                                   146 
181    G4ThreeVector boostLABtoNN ( 0. , 0. , beta    147    G4ThreeVector boostLABtoNN ( 0. , 0. , beta_nn ); // CM of NN; 
182                                                   148 
183    G4double beta_nncm = ( - boostLABtoCM.beta(    149    G4double beta_nncm = ( - boostLABtoCM.beta() + boostLABtoNN.beta() ) / ( 1 - boostLABtoCM.beta() * boostLABtoNN.beta() ) ;  
184                                                   150 
185    //std::cout << targ4p << std::endl;            151    //std::cout << targ4p << std::endl; 
186    //std::cout << proj_dp->Get4Momentum()<< st    152    //std::cout << proj_dp->Get4Momentum()<< std::endl; 
187    //std::cout << beta_nncm << std::endl;         153    //std::cout << beta_nncm << std::endl; 
188    G4ThreeVector boostNNtoCM( 0. , 0. , beta_n    154    G4ThreeVector boostNNtoCM( 0. , 0. , beta_nncm ); // 
189    G4ThreeVector boostCMtoNN( 0. , 0. , -beta_    155    G4ThreeVector boostCMtoNN( 0. , 0. , -beta_nncm ); // 
190                                                   156 
191    boostToReac = boostLABtoNN;                    157    boostToReac = boostLABtoNN; 
192    boostBackToLAB = -boostLABtoNN;                158    boostBackToLAB = -boostLABtoNN; 
193                                                   159 
194    delete proj_dp;                                160    delete proj_dp; 
195                                                   161 
196    G4int icounter = 0;                         << 162    while ( elastic ) 
197    G4int icounter_max = 1024;                  << 163    {
198    while ( elastic ) // Loop checking, 11.03.2 << 
199    {                                           << 
200       icounter++;                              << 
201       if ( icounter > icounter_max ) {         << 
202    G4cout << "Loop-counter exceeded the thresh << 
203          break;                                << 
204       }                                        << 
205                                                   164 
206 // impact parameter                               165 // impact parameter 
207       //G4double bmax = 1.05*(bmax_0/fermi);      166       //G4double bmax = 1.05*(bmax_0/fermi);  // 10% for Peripheral reactions
208       G4double bmax = envelopF*(bmax_0/fermi);    167       G4double bmax = envelopF*(bmax_0/fermi);
209       G4double b = bmax * std::sqrt ( G4Unifor    168       G4double b = bmax * std::sqrt ( G4UniformRand() );
210 //071112                                          169 //071112
211       //G4double b = 0;                           170       //G4double b = 0;
212       //G4double b = bmax;                        171       //G4double b = bmax;
213       //G4double b = bmax/1.05 * 0.7 * G4Unifo    172       //G4double b = bmax/1.05 * 0.7 * G4UniformRand();
214                                                   173 
215       //G4cout << "G4QMDRESULT bmax_0 = " << b    174       //G4cout << "G4QMDRESULT bmax_0 = " << bmax_0/fermi << " fm, bmax = " << bmax << " fm , b = " << b  << " fm " << G4endl; 
216                                                   175 
217       G4double plab = projectile.GetTotalMomen    176       G4double plab = projectile.GetTotalMomentum()/GeV;
218       G4double elab = ( projectile.GetKineticE    177       G4double elab = ( projectile.GetKineticEnergy() + proj_pd->GetPDGMass() + targ_pd->GetPDGMass() )/GeV;
219                                                   178 
220       calcOffSetOfCollision( b , proj_pd , tar    179       calcOffSetOfCollision( b , proj_pd , targ_pd , plab , elab , bmax , boostCMtoNN );
221                                                   180 
222 // Projectile                                     181 // Projectile
223       G4LorentzVector proj4pLAB = projectile.G    182       G4LorentzVector proj4pLAB = projectile.Get4Momentum()/GeV;
224                                                   183 
225       G4QMDGroundStateNucleus* proj(NULL);        184       G4QMDGroundStateNucleus* proj(NULL); 
226       if ( projectile.GetDefinition()->GetPart    185       if ( projectile.GetDefinition()->GetParticleType() == "nucleus" 
227         || projectile.GetDefinition()->GetPart    186         || projectile.GetDefinition()->GetParticleName() == "proton"
228         || projectile.GetDefinition()->GetPart    187         || projectile.GetDefinition()->GetParticleName() == "neutron" )
229       {                                           188       {
230                                                   189 
231          proj_Z = proj_pd->GetAtomicNumber();     190          proj_Z = proj_pd->GetAtomicNumber();
232          proj_A = proj_pd->GetAtomicMass();       191          proj_A = proj_pd->GetAtomicMass();
233                                                   192 
234          proj = new G4QMDGroundStateNucleus( p    193          proj = new G4QMDGroundStateNucleus( proj_Z , proj_A );
235          //proj->ShowParticipants();              194          //proj->ShowParticipants();
236                                                   195 
237                                                   196 
238          meanField->SetSystem ( proj );           197          meanField->SetSystem ( proj );
239          proj->SetTotalPotential( meanField->G    198          proj->SetTotalPotential( meanField->GetTotalPotential() );
240          proj->CalEnergyAndAngularMomentumInCM    199          proj->CalEnergyAndAngularMomentumInCM();
241                                                   200 
242       }                                           201       }
243                                                   202 
244 // Target                                         203 // Target
245       //G4int iz = int ( target.GetZ() );         204       //G4int iz = int ( target.GetZ() );
246       //G4int ia = int ( target.GetN() );         205       //G4int ia = int ( target.GetN() );
247       //migrate to integer A and Z (GetN_asInt    206       //migrate to integer A and Z (GetN_asInt returns number of neutrons in the nucleus since this) 
248       G4int iz = int ( target.GetZ_asInt() );     207       G4int iz = int ( target.GetZ_asInt() );
249       G4int ia = int ( target.GetA_asInt() );     208       G4int ia = int ( target.GetA_asInt() );
250                                                   209 
251       G4QMDGroundStateNucleus* targ = new G4QM    210       G4QMDGroundStateNucleus* targ = new G4QMDGroundStateNucleus( iz , ia );
252                                                   211 
253       meanField->SetSystem (targ );               212       meanField->SetSystem (targ );
254       targ->SetTotalPotential( meanField->GetT    213       targ->SetTotalPotential( meanField->GetTotalPotential() );
255       targ->CalEnergyAndAngularMomentumInCM();    214       targ->CalEnergyAndAngularMomentumInCM();
256                                                   215    
257       //G4LorentzVector targ4p( G4ThreeVector(    216       //G4LorentzVector targ4p( G4ThreeVector( 0.0 ) , targ->GetNuclearMass()/GeV );
258 // Boost Vector to CM                             217 // Boost Vector to CM
259       //boostToCM = targ4p.findBoostToCM( proj    218       //boostToCM = targ4p.findBoostToCM( proj4pLAB );
260                                                   219 
261 //    Target                                      220 //    Target 
262       for ( G4int i = 0 ; i < targ->GetTotalNu    221       for ( G4int i = 0 ; i < targ->GetTotalNumberOfParticipant() ; i++ )
263       {                                           222       {
264                                                   223 
265          G4ThreeVector p0 = targ->GetParticipa    224          G4ThreeVector p0 = targ->GetParticipant( i )->GetMomentum();
266          G4ThreeVector r0 = targ->GetParticipa    225          G4ThreeVector r0 = targ->GetParticipant( i )->GetPosition();
267                                                   226 
268          G4ThreeVector p ( p0.x() + coulomb_co    227          G4ThreeVector p ( p0.x() + coulomb_collision_px_targ 
269                          , p0.y()                 228                          , p0.y() 
270                          , p0.z() * coulomb_co    229                          , p0.z() * coulomb_collision_gamma_targ + coulomb_collision_pz_targ ); 
271                                                   230 
272          G4ThreeVector r ( r0.x() + coulomb_co    231          G4ThreeVector r ( r0.x() + coulomb_collision_rx_targ 
273                          , r0.y()                 232                          , r0.y() 
274                          , r0.z() / coulomb_co    233                          , r0.z() / coulomb_collision_gamma_targ + coulomb_collision_rz_targ ); 
275                                                   234      
276          system->SetParticipant( new G4QMDPart    235          system->SetParticipant( new G4QMDParticipant( targ->GetParticipant( i )->GetDefinition() , p , r ) );
277          system->GetParticipant( i )->SetTarge    236          system->GetParticipant( i )->SetTarget();
278                                                   237 
279       }                                           238       }
280                                                   239 
281       G4LorentzVector proj4pCM = CLHEP::boostO    240       G4LorentzVector proj4pCM = CLHEP::boostOf ( proj4pLAB , boostToReac );
282       G4LorentzVector targ4pCM = CLHEP::boostO    241       G4LorentzVector targ4pCM = CLHEP::boostOf ( targ4p , boostToReac );
283                                                   242 
284 //    Projectile                                  243 //    Projectile
285       //G4cout << "proj : " << proj << G4endl; << 244       if ( proj != NULL )
286       //if ( proj != NULL )                    << 
287       if ( proj_A != 1 )                       << 
288       {                                           245       {
289                                                   246 
290 //    projectile is nucleus                       247 //    projectile is nucleus
291                                                   248 
292          for ( G4int i = 0 ; i < proj->GetTota    249          for ( G4int i = 0 ; i < proj->GetTotalNumberOfParticipant() ; i++ )
293          {                                        250          {
294                                                   251 
295             G4ThreeVector p0 = proj->GetPartic    252             G4ThreeVector p0 = proj->GetParticipant( i )->GetMomentum();
296             G4ThreeVector r0 = proj->GetPartic    253             G4ThreeVector r0 = proj->GetParticipant( i )->GetPosition();
297                                                   254 
298             G4ThreeVector p ( p0.x() + coulomb    255             G4ThreeVector p ( p0.x() + coulomb_collision_px_proj 
299                             , p0.y()              256                             , p0.y() 
300                             , p0.z() * coulomb    257                             , p0.z() * coulomb_collision_gamma_proj + coulomb_collision_pz_proj ); 
301                                                   258 
302             G4ThreeVector r ( r0.x() + coulomb    259             G4ThreeVector r ( r0.x() + coulomb_collision_rx_proj 
303                             , r0.y()              260                             , r0.y() 
304                             , r0.z() / coulomb    261                             , r0.z() / coulomb_collision_gamma_proj + coulomb_collision_rz_proj ); 
305                                                   262      
306             system->SetParticipant( new G4QMDP    263             system->SetParticipant( new G4QMDParticipant( proj->GetParticipant( i )->GetDefinition() , p  , r ) );
307             system->GetParticipant ( i + targ-    264             system->GetParticipant ( i + targ->GetTotalNumberOfParticipant() )->SetProjectile();
308          }                                        265          }
309                                                   266 
310       }                                           267       }
311       else                                        268       else
312       {                                           269       {
313                                                   270 
314 //       projectile is particle                   271 //       projectile is particle
315                                                   272 
316          // avoid multiple set in "elastic" lo    273          // avoid multiple set in "elastic" loop
317          //G4cout << "system Total Participant << 
318          if ( system->GetTotalNumberOfParticip    274          if ( system->GetTotalNumberOfParticipant() == targ->GetTotalNumberOfParticipant() )
319          {                                        275          {
320                                                   276 
321             G4int i = targ->GetTotalNumberOfPa    277             G4int i = targ->GetTotalNumberOfParticipant(); 
322                                                   278       
323             G4ThreeVector p0( 0 );                279             G4ThreeVector p0( 0 ); 
324             G4ThreeVector r0( 0 );                280             G4ThreeVector r0( 0 );
325                                                   281 
326             G4ThreeVector p ( p0.x() + coulomb    282             G4ThreeVector p ( p0.x() + coulomb_collision_px_proj 
327                             , p0.y()              283                             , p0.y() 
328                             , p0.z() * coulomb    284                             , p0.z() * coulomb_collision_gamma_proj + coulomb_collision_pz_proj ); 
329                                                   285 
330             G4ThreeVector r ( r0.x() + coulomb    286             G4ThreeVector r ( r0.x() + coulomb_collision_rx_proj 
331                             , r0.y()              287                             , r0.y() 
332                             , r0.z() / coulomb    288                             , r0.z() / coulomb_collision_gamma_proj + coulomb_collision_rz_proj ); 
333                                                   289 
334             system->SetParticipant( new G4QMDP    290             system->SetParticipant( new G4QMDParticipant( (G4ParticleDefinition*)projectile.GetDefinition() , p , r ) );
335             // This is not important becase on    291             // This is not important becase only 1 projectile particle.
336             system->GetParticipant ( i )->SetP    292             system->GetParticipant ( i )->SetProjectile();
337          }                                        293          }
338                                                   294 
339       }                                           295       }
340       //system->ShowParticipants();               296       //system->ShowParticipants();
341                                                   297 
342       delete targ;                                298       delete targ;
343       delete proj;                                299       delete proj;
344                                                   300 
345    meanField->SetSystem ( system );               301    meanField->SetSystem ( system );
346    collision->SetMeanField ( meanField );         302    collision->SetMeanField ( meanField );
347                                                   303 
348 // Time Evolution                                 304 // Time Evolution 
349    //std::cout << "Start time evolution " << s    305    //std::cout << "Start time evolution " << std::endl;
350    //system->ShowParticipants();                  306    //system->ShowParticipants();
351    for ( G4int i = 0 ; i < maxTime ; i++ )        307    for ( G4int i = 0 ; i < maxTime ; i++ )
352    {                                              308    {
353       //G4cout << " do Paropagate " << i << "     309       //G4cout << " do Paropagate " << i << " th time step. " << G4endl;
354       meanField->DoPropagation( deltaT );         310       meanField->DoPropagation( deltaT );
355       //system->ShowParticipants();               311       //system->ShowParticipants();
356       collision->CalKinematicsOfBinaryCollisio    312       collision->CalKinematicsOfBinaryCollisions( deltaT );
357                                                   313 
358       if ( i / 10 * 10 == i )                     314       if ( i / 10 * 10 == i ) 
359       {                                           315       {
360          //G4cout << i << " th time step. " <<    316          //G4cout << i << " th time step. " << G4endl;
361          //system->ShowParticipants();            317          //system->ShowParticipants();
362       }                                           318       } 
363       //system->ShowParticipants();               319       //system->ShowParticipants();
364    }                                              320    }
365    //system->ShowParticipants();                  321    //system->ShowParticipants();
366                                                   322 
367                                                   323 
368    //std::cout << "Doing Cluster Judgment " <<    324    //std::cout << "Doing Cluster Judgment " << std::endl;
369                                                   325 
370    nucleuses = meanField->DoClusterJudgment();    326    nucleuses = meanField->DoClusterJudgment();
371                                                   327 
372 // Elastic Judgment                               328 // Elastic Judgment  
373                                                   329 
374    G4int numberOfSecondary = int ( nucleuses.s    330    G4int numberOfSecondary = int ( nucleuses.size() ) + system->GetTotalNumberOfParticipant(); 
375                                                   331 
376    G4int sec_a_Z = 0;                             332    G4int sec_a_Z = 0;
377    G4int sec_a_A = 0;                             333    G4int sec_a_A = 0;
378    const G4ParticleDefinition* sec_a_pd = NULL << 334    G4ParticleDefinition* sec_a_pd = NULL;
379    G4int sec_b_Z = 0;                             335    G4int sec_b_Z = 0;
380    G4int sec_b_A = 0;                             336    G4int sec_b_A = 0;
381    const G4ParticleDefinition* sec_b_pd = NULL << 337    G4ParticleDefinition* sec_b_pd = NULL;
382                                                   338 
383    if ( numberOfSecondary == 2 )                  339    if ( numberOfSecondary == 2 )
384    {                                              340    {
385                                                   341 
386       G4bool elasticLike_system = false;          342       G4bool elasticLike_system = false;
387       if ( nucleuses.size() == 2 )                343       if ( nucleuses.size() == 2 ) 
388       {                                           344       {
389                                                   345 
390          sec_a_Z = nucleuses[0]->GetAtomicNumb    346          sec_a_Z = nucleuses[0]->GetAtomicNumber();
391          sec_a_A = nucleuses[0]->GetMassNumber    347          sec_a_A = nucleuses[0]->GetMassNumber();
392          sec_b_Z = nucleuses[1]->GetAtomicNumb    348          sec_b_Z = nucleuses[1]->GetAtomicNumber();
393          sec_b_A = nucleuses[1]->GetMassNumber    349          sec_b_A = nucleuses[1]->GetMassNumber();
394                                                   350 
395          if ( ( sec_a_Z == proj_Z && sec_a_A =    351          if ( ( sec_a_Z == proj_Z && sec_a_A == proj_A && sec_b_Z == targ_Z && sec_b_A == targ_A )
396            || ( sec_a_Z == targ_Z && sec_a_A =    352            || ( sec_a_Z == targ_Z && sec_a_A == targ_A && sec_b_Z == proj_Z && sec_b_A == proj_A ) )
397          {                                        353          {
398             elasticLike_system = true;            354             elasticLike_system = true;
399          }                                        355          } 
400                                                   356 
401       }                                           357       }
402       else if ( nucleuses.size() == 1 )           358       else if ( nucleuses.size() == 1 ) 
403       {                                           359       {
404                                                   360 
405          sec_a_Z = nucleuses[0]->GetAtomicNumb    361          sec_a_Z = nucleuses[0]->GetAtomicNumber();
406          sec_a_A = nucleuses[0]->GetMassNumber    362          sec_a_A = nucleuses[0]->GetMassNumber();
407          sec_b_pd = system->GetParticipant( 0     363          sec_b_pd = system->GetParticipant( 0 )->GetDefinition();
408                                                   364 
409          if ( ( sec_a_Z == proj_Z && sec_a_A =    365          if ( ( sec_a_Z == proj_Z && sec_a_A == proj_A && sec_b_pd == targ_pd )
410            || ( sec_a_Z == targ_Z && sec_a_A =    366            || ( sec_a_Z == targ_Z && sec_a_A == targ_A && sec_b_pd == proj_pd ) )
411          {                                        367          {
412             elasticLike_system = true;            368             elasticLike_system = true;
413          }                                        369          } 
414                                                   370 
415       }                                           371       }  
416       else                                        372       else
417       {                                           373       {
418                                                   374 
419          sec_a_pd = system->GetParticipant( 0     375          sec_a_pd = system->GetParticipant( 0 )->GetDefinition();
420          sec_b_pd = system->GetParticipant( 1     376          sec_b_pd = system->GetParticipant( 1 )->GetDefinition();
421                                                   377  
422          if ( ( sec_a_pd == proj_pd && sec_b_p    378          if ( ( sec_a_pd == proj_pd && sec_b_pd == targ_pd ) 
423            || ( sec_a_pd == targ_pd && sec_b_p    379            || ( sec_a_pd == targ_pd && sec_b_pd == proj_pd ) ) 
424          {                                        380          {
425             elasticLike_system = true;            381             elasticLike_system = true;
426          }                                        382          } 
427                                                   383 
428       }                                           384       } 
429                                                   385 
430       if ( elasticLike_system == true )           386       if ( elasticLike_system == true )
431       {                                           387       {
432                                                   388 
433          G4bool elasticLike_energy = true;        389          G4bool elasticLike_energy = true;
434 //    Cal ExcitationEnergy                        390 //    Cal ExcitationEnergy 
435          for ( G4int i = 0 ; i < int ( nucleus    391          for ( G4int i = 0 ; i < int ( nucleuses.size() ) ; i++ )
436          {                                        392          { 
437                                                   393 
438             //meanField->SetSystem( nucleuses[    394             //meanField->SetSystem( nucleuses[i] );
439             meanField->SetNucleus( nucleuses[i    395             meanField->SetNucleus( nucleuses[i] );
440             //nucleuses[i]->SetTotalPotential(    396             //nucleuses[i]->SetTotalPotential( meanField->GetTotalPotential() );
441             //nucleuses[i]->CalEnergyAndAngula    397             //nucleuses[i]->CalEnergyAndAngularMomentumInCM();
442                                                   398 
443             if ( nucleuses[i]->GetExcitationEn    399             if ( nucleuses[i]->GetExcitationEnergy()*GeV > 1.0*MeV ) elasticLike_energy = false;  
444                                                   400 
445          }                                        401          } 
446                                                   402 
447 //    Check Collision                             403 //    Check Collision 
448          G4bool withCollision = true;             404          G4bool withCollision = true;
449          if ( system->GetNOCollision() == 0 )     405          if ( system->GetNOCollision() == 0 ) withCollision = false;
450                                                   406 
451 //    Final judegement for Inelasitc or Elasti    407 //    Final judegement for Inelasitc or Elastic;
452 //                                                408 //
453 //       ElasticLike without Collision            409 //       ElasticLike without Collision 
454          //if ( elasticLike_energy == true &&     410          //if ( elasticLike_energy == true && withCollision == false ) elastic = true;  // ielst = 0
455 //       ElasticLike with Collision               411 //       ElasticLike with Collision 
456          //if ( elasticLike_energy == true &&     412          //if ( elasticLike_energy == true && withCollision == true ) elastic = true;   // ielst = 1 
457 //       InelasticLike without Collision          413 //       InelasticLike without Collision 
458          //if ( elasticLike_energy == false )     414          //if ( elasticLike_energy == false ) elastic = false;                          // ielst = 2                
459          if ( frag == true )                      415          if ( frag == true )
460             if ( elasticLike_energy == false )    416             if ( elasticLike_energy == false ) elastic = false;
461 //       InelasticLike with Collision             417 //       InelasticLike with Collision 
462          if ( elasticLike_energy == false && w    418          if ( elasticLike_energy == false && withCollision == true ) elastic = false; // ielst = 3
463                                                   419 
464       }                                           420       }
465                                                   421 
466       }                                           422       }
467       else                                        423       else
468       {                                           424       {
469                                                   425 
470 //       numberOfSecondary != 2                   426 //       numberOfSecondary != 2 
471          elastic = false;                         427          elastic = false;
472                                                   428 
473       }                                           429       }
474                                                   430 
475 //071115                                          431 //071115
476       //G4cout << elastic << G4endl;              432       //G4cout << elastic << G4endl;
477       // if elastic is true try again from sam    433       // if elastic is true try again from sampling of impact parameter 
478                                                   434 
479       if ( elastic == true )                      435       if ( elastic == true )
480       {                                           436       {
481          // delete this nucleues                  437          // delete this nucleues
482          for ( std::vector< G4QMDNucleus* >::i    438          for ( std::vector< G4QMDNucleus* >::iterator
483                it = nucleuses.begin() ; it !=     439                it = nucleuses.begin() ; it != nucleuses.end() ; it++ )
484          {                                        440          {
485             delete *it;                           441             delete *it;
486          }                                        442          }
487          nucleuses.clear();                       443          nucleuses.clear();
488       }                                           444       }
489                                                << 
490    }                                              445    } 
491                                                   446 
492                                                   447 
493 // Statical Decay Phase                           448 // Statical Decay Phase
494                                                   449 
495    for ( std::vector< G4QMDNucleus* >::iterato    450    for ( std::vector< G4QMDNucleus* >::iterator it
496        = nucleuses.begin() ; it != nucleuses.e    451        = nucleuses.begin() ; it != nucleuses.end() ; it++ )
497    {                                              452    {
498                                                   453 
499 /*                                                454 /*
500       G4cout << "G4QMDRESULT "                    455       G4cout << "G4QMDRESULT "
501              << (*it)->GetAtomicNumber()          456              << (*it)->GetAtomicNumber() 
502              << " "                               457              << " " 
503              << (*it)->GetMassNumber()            458              << (*it)->GetMassNumber() 
504              << " "                               459              << " " 
505              << (*it)->Get4Momentum()             460              << (*it)->Get4Momentum() 
506              << " "                               461              << " " 
507              << (*it)->Get4Momentum().vect()      462              << (*it)->Get4Momentum().vect() 
508              << " "                               463              << " " 
509              << (*it)->Get4Momentum().restMass    464              << (*it)->Get4Momentum().restMass() 
510              << " "                               465              << " " 
511              << (*it)->GetNuclearMass()/GeV       466              << (*it)->GetNuclearMass()/GeV 
512              << G4endl;                           467              << G4endl;
513 */                                                468 */
514                                                   469 
515       meanField->SetNucleus ( *it );              470       meanField->SetNucleus ( *it );
516                                                   471 
517       if ( (*it)->GetAtomicNumber() == 0  // n    472       if ( (*it)->GetAtomicNumber() == 0  // neutron cluster
518         || (*it)->GetAtomicNumber() == (*it)->    473         || (*it)->GetAtomicNumber() == (*it)->GetMassNumber() ) // proton cluster
519       {                                           474       {
520          // push back system                      475          // push back system 
521          for ( G4int i = 0 ; i < (*it)->GetTot    476          for ( G4int i = 0 ; i < (*it)->GetTotalNumberOfParticipant() ; i++ )
522          {                                        477          {
523             G4QMDParticipant* aP = new G4QMDPa    478             G4QMDParticipant* aP = new G4QMDParticipant( ( (*it)->GetParticipant( i ) )->GetDefinition() , ( (*it)->GetParticipant( i ) )->GetMomentum() , ( (*it)->GetParticipant( i ) )->GetPosition() );  
524             system->SetParticipant ( aP );        479             system->SetParticipant ( aP );  
525          }                                        480          } 
526          continue;                                481          continue;  
527       }                                           482       }
528                                                   483 
529       G4double nucleus_e = std::sqrt ( G4Pow:: << 484       G4double nucleus_e = std::sqrt ( std::pow ( (*it)->GetNuclearMass()/GeV , 2 ) + std::pow ( (*it)->Get4Momentum().vect().mag() , 2 ) );
530       G4LorentzVector nucleus_p4CM ( (*it)->Ge    485       G4LorentzVector nucleus_p4CM ( (*it)->Get4Momentum().vect() , nucleus_e ); 
531                                                   486 
532 //      std::cout << "G4QMDRESULT nucleus delt    487 //      std::cout << "G4QMDRESULT nucleus deltaQ " << deltaQ << std::endl;
533                                                   488 
534       G4int ia = (*it)->GetMassNumber();          489       G4int ia = (*it)->GetMassNumber();
535       G4int iz = (*it)->GetAtomicNumber();        490       G4int iz = (*it)->GetAtomicNumber();
536                                                   491 
537       G4LorentzVector lv ( G4ThreeVector( 0.0     492       G4LorentzVector lv ( G4ThreeVector( 0.0 ) , (*it)->GetExcitationEnergy()*GeV + G4IonTable::GetIonTable()->GetIonMass( iz , ia ) );
538                                                   493 
539       G4Fragment* aFragment = new G4Fragment(     494       G4Fragment* aFragment = new G4Fragment( ia , iz , lv );
540                                                   495 
541       G4ReactionProductVector* rv;                496       G4ReactionProductVector* rv;
542       rv = excitationHandler->BreakItUp( *aFra    497       rv = excitationHandler->BreakItUp( *aFragment );
543       G4bool notBreak = true;                     498       G4bool notBreak = true;
544       for ( G4ReactionProductVector::iterator     499       for ( G4ReactionProductVector::iterator itt
545           = rv->begin() ; itt != rv->end() ; i    500           = rv->begin() ; itt != rv->end() ; itt++ )
546       {                                           501       {
547                                                   502 
548           notBreak = false;                       503           notBreak = false;
549           // Secondary from this nucleus (*it)    504           // Secondary from this nucleus (*it) 
550           const G4ParticleDefinition* pd = (*i << 505           G4ParticleDefinition* pd = (*itt)->GetDefinition();
551                                                   506 
552           G4LorentzVector p4 ( (*itt)->GetMome    507           G4LorentzVector p4 ( (*itt)->GetMomentum()/GeV , (*itt)->GetTotalEnergy()/GeV );  //in nucleus(*it) rest system
553           G4LorentzVector p4_CM = CLHEP::boost    508           G4LorentzVector p4_CM = CLHEP::boostOf( p4 , -nucleus_p4CM.findBoostToCM() );  // Back to CM
554           G4LorentzVector p4_LAB = CLHEP::boos    509           G4LorentzVector p4_LAB = CLHEP::boostOf( p4_CM , boostBackToLAB ); // Back to LAB  
555                                                   510 
556                                                   511 
557 //090122                                          512 //090122
558           //theParticleChange.AddSecondary( dp    513           //theParticleChange.AddSecondary( dp ); 
559           if ( !( pd->GetAtomicNumber() == 4 &    514           if ( !( pd->GetAtomicNumber() == 4 && pd->GetAtomicMass() == 8 ) )
560           {                                       515           {
561              //G4cout << "pd out of notBreak l << 
562              G4DynamicParticle* dp = new G4Dyn    516              G4DynamicParticle* dp = new G4DynamicParticle( pd , p4_LAB*GeV );  
563              theParticleChange.AddSecondary( d    517              theParticleChange.AddSecondary( dp ); 
564           }                                       518           }
565           else                                    519           else
566           {                                       520           {
567              //Be8 -> Alpha + Alpha + Q           521              //Be8 -> Alpha + Alpha + Q
568              G4ThreeVector randomized_directio    522              G4ThreeVector randomized_direction( G4UniformRand() , G4UniformRand() , G4UniformRand() );
569              randomized_direction = randomized    523              randomized_direction = randomized_direction.unit();
570              G4double q_decay = (*itt)->GetMas    524              G4double q_decay = (*itt)->GetMass() - 2*G4Alpha::Alpha()->GetPDGMass();
571              G4double p_decay = std::sqrt ( G4 << 525              G4double p_decay = std::sqrt ( std::pow(G4Alpha::Alpha()->GetPDGMass()+q_decay/2,2) - std::pow(G4Alpha::Alpha()->GetPDGMass() , 2 ) ); 
572              G4LorentzVector p4_a1 ( p_decay*r    526              G4LorentzVector p4_a1 ( p_decay*randomized_direction , G4Alpha::Alpha()->GetPDGMass()+q_decay/2 );  //in Be8 rest system
573                                                   527              
574              G4LorentzVector p4_a1_Be8 = CLHEP    528              G4LorentzVector p4_a1_Be8 = CLHEP::boostOf ( p4_a1/GeV , -p4.findBoostToCM() );
575              G4LorentzVector p4_a1_CM = CLHEP:    529              G4LorentzVector p4_a1_CM = CLHEP::boostOf ( p4_a1_Be8 , -nucleus_p4CM.findBoostToCM() );
576              G4LorentzVector p4_a1_LAB = CLHEP    530              G4LorentzVector p4_a1_LAB = CLHEP::boostOf ( p4_a1_CM , boostBackToLAB );
577                                                   531 
578              G4LorentzVector p4_a2 ( -p_decay*    532              G4LorentzVector p4_a2 ( -p_decay*randomized_direction , G4Alpha::Alpha()->GetPDGMass()+q_decay/2 );  //in Be8 rest system
579                                                   533              
580              G4LorentzVector p4_a2_Be8 = CLHEP    534              G4LorentzVector p4_a2_Be8 = CLHEP::boostOf ( p4_a2/GeV , -p4.findBoostToCM() );
581              G4LorentzVector p4_a2_CM = CLHEP:    535              G4LorentzVector p4_a2_CM = CLHEP::boostOf ( p4_a2_Be8 , -nucleus_p4CM.findBoostToCM() );
582              G4LorentzVector p4_a2_LAB = CLHEP    536              G4LorentzVector p4_a2_LAB = CLHEP::boostOf ( p4_a2_CM , boostBackToLAB );
583                                                   537              
584              G4DynamicParticle* dp1 = new G4Dy    538              G4DynamicParticle* dp1 = new G4DynamicParticle( G4Alpha::Alpha() , p4_a1_LAB*GeV );  
585              G4DynamicParticle* dp2 = new G4Dy    539              G4DynamicParticle* dp2 = new G4DynamicParticle( G4Alpha::Alpha() , p4_a2_LAB*GeV );  
586              theParticleChange.AddSecondary( d    540              theParticleChange.AddSecondary( dp1 ); 
587              theParticleChange.AddSecondary( d    541              theParticleChange.AddSecondary( dp2 ); 
588           }                                       542           }
589 //090122                                          543 //090122
590                                                   544 
591 /*                                                545 /*
592           G4cout                                  546           G4cout
593                 << "Regist Secondary "            547                 << "Regist Secondary "
594                 << (*itt)->GetDefinition()->Ge    548                 << (*itt)->GetDefinition()->GetParticleName()
595                 << " "                            549                 << " "
596                 << (*itt)->GetMomentum()/GeV      550                 << (*itt)->GetMomentum()/GeV
597                 << " "                            551                 << " "
598                 << (*itt)->GetKineticEnergy()/    552                 << (*itt)->GetKineticEnergy()/GeV
599                 << " "                            553                 << " "
600                 << (*itt)->GetMass()/GeV          554                 << (*itt)->GetMass()/GeV
601                 << " "                            555                 << " "
602                 << (*itt)->GetTotalEnergy()/Ge    556                 << (*itt)->GetTotalEnergy()/GeV
603                 << " "                            557                 << " "
604                 << (*itt)->GetTotalEnergy()/Ge    558                 << (*itt)->GetTotalEnergy()/GeV * (*itt)->GetTotalEnergy()/GeV
605                  - (*itt)->GetMomentum()/GeV *    559                  - (*itt)->GetMomentum()/GeV * (*itt)->GetMomentum()/GeV
606                 << " "                            560                 << " "
607                 << nucleus_p4CM.findBoostToCM(    561                 << nucleus_p4CM.findBoostToCM() 
608                 << " "                            562                 << " "
609                 << p4                             563                 << p4
610                 << " "                            564                 << " "
611                 << p4_CM                          565                 << p4_CM
612                 << " "                            566                 << " "
613                 << p4_LAB                         567                 << p4_LAB
614                 << G4endl;                        568                 << G4endl;
615 */                                                569 */
616                                                   570 
617       }                                           571       }
618       if ( notBreak == true )                     572       if ( notBreak == true )
619       {                                           573       {
620                                                   574 
621          const G4ParticleDefinition* pd = G4Io << 575          G4ParticleDefinition* pd = G4IonTable::GetIonTable()->GetIon( (*it)->GetAtomicNumber() , (*it)->GetMassNumber(), (*it)->GetExcitationEnergy()*GeV );
622              //G4cout << "pd in notBreak loop  << 
623          G4LorentzVector p4_CM = nucleus_p4CM;    576          G4LorentzVector p4_CM = nucleus_p4CM;
624          G4LorentzVector p4_LAB = CLHEP::boost    577          G4LorentzVector p4_LAB = CLHEP::boostOf( p4_CM , boostBackToLAB ); // Back to LAB  
625          G4DynamicParticle* dp = new G4Dynamic    578          G4DynamicParticle* dp = new G4DynamicParticle( pd , p4_LAB*GeV );  
626          theParticleChange.AddSecondary( dp );    579          theParticleChange.AddSecondary( dp ); 
627                                                   580 
628       }                                           581       }
629                                                   582 
630       for ( G4ReactionProductVector::iterator     583       for ( G4ReactionProductVector::iterator itt
631           = rv->begin() ; itt != rv->end() ; i    584           = rv->begin() ; itt != rv->end() ; itt++ )
632       {                                           585       {
633           delete *itt;                            586           delete *itt;
634       }                                           587       }
635       delete rv;                                  588       delete rv;
636                                                   589 
637       delete aFragment;                           590       delete aFragment;
638                                                   591 
639    }                                              592    }
640                                                   593 
641                                                   594 
642                                                   595 
643    for ( G4int i = 0 ; i < system->GetTotalNum    596    for ( G4int i = 0 ; i < system->GetTotalNumberOfParticipant() ; i++ )
644    {                                              597    {
                                                   >> 598 
645       // Secondary particles                      599       // Secondary particles 
646                                                   600 
647       const G4ParticleDefinition* pd = system- << 601       G4ParticleDefinition* pd = system->GetParticipant( i )->GetDefinition();
648       G4LorentzVector p4_CM = system->GetParti    602       G4LorentzVector p4_CM = system->GetParticipant( i )->Get4Momentum();
649       G4LorentzVector p4_LAB = CLHEP::boostOf(    603       G4LorentzVector p4_LAB = CLHEP::boostOf( p4_CM , boostBackToLAB );
650       G4DynamicParticle* dp = new G4DynamicPar    604       G4DynamicParticle* dp = new G4DynamicParticle( pd , p4_LAB*GeV );  
651       theParticleChange.AddSecondary( dp );       605       theParticleChange.AddSecondary( dp ); 
652       //G4cout << "In the last theParticleChan << 
653                                                   606 
654 /*                                                607 /*
655       G4cout << "G4QMDRESULT "                    608       G4cout << "G4QMDRESULT "
656       << "r" << i << " " << system->GetPartici    609       << "r" << i << " " << system->GetParticipant ( i ) -> GetPosition() << " "
657       << "p" << i << " " << system->GetPartici    610       << "p" << i << " " << system->GetParticipant ( i ) -> Get4Momentum()
658       << G4endl;                                  611       << G4endl;
659 */                                                612 */
660                                                   613 
661    }                                              614    }
662                                                   615 
663    for ( std::vector< G4QMDNucleus* >::iterato    616    for ( std::vector< G4QMDNucleus* >::iterator it
664        = nucleuses.begin() ; it != nucleuses.e    617        = nucleuses.begin() ; it != nucleuses.end() ; it++ )
665    {                                              618    {
666       delete *it;  // delete nulceuse             619       delete *it;  // delete nulceuse 
667    }                                              620    }
668    nucleuses.clear();                             621    nucleuses.clear();
669                                                   622 
670    system->Clear();                               623    system->Clear();
671    delete system;                                 624    delete system; 
672                                                   625 
673    theParticleChange.SetStatusChange( stopAndK    626    theParticleChange.SetStatusChange( stopAndKill );
674                                                   627 
675    for (G4int i = 0; i < G4int(theParticleChan << 
676    {                                           << 
677      //G4cout << "Particle : " << theParticleC << 
678      //G4cout << "KEnergy : " << theParticleCh << 
679      //G4cout << "modelID : " << theParticleCh << 
680      theParticleChange.GetSecondary(i)->SetCre << 
681    }                                           << 
682                                                << 
683    return &theParticleChange;                     628    return &theParticleChange;
684                                                   629 
685 }                                                 630 }
686                                                   631 
687                                                   632 
688                                                   633 
689 void G4QMDReaction::calcOffSetOfCollision( G4d    634 void G4QMDReaction::calcOffSetOfCollision( G4double b , 
690 const G4ParticleDefinition* pd_proj ,          << 635 G4ParticleDefinition* pd_proj , 
691 const G4ParticleDefinition* pd_targ ,          << 636 G4ParticleDefinition* pd_targ , 
692 G4double ptot , G4double etot , G4double bmax     637 G4double ptot , G4double etot , G4double bmax , G4ThreeVector boostToCM )
693 {                                                 638 {
694                                                   639 
695    G4double mass_proj = pd_proj->GetPDGMass()/    640    G4double mass_proj = pd_proj->GetPDGMass()/GeV;
696    G4double mass_targ = pd_targ->GetPDGMass()/    641    G4double mass_targ = pd_targ->GetPDGMass()/GeV;
697                                                   642   
698    G4double stot = std::sqrt ( etot*etot - pto    643    G4double stot = std::sqrt ( etot*etot - ptot*ptot );
699                                                   644 
700    G4double pstt = std::sqrt ( ( stot*stot - (    645    G4double pstt = std::sqrt ( ( stot*stot - ( mass_proj + mass_targ ) * ( mass_proj + mass_targ ) 
701                   ) * ( stot*stot - ( mass_pro    646                   ) * ( stot*stot - ( mass_proj - mass_targ ) * ( mass_proj - mass_targ ) ) ) 
702                  / ( 2.0 * stot );                647                  / ( 2.0 * stot );
703                                                   648 
704    G4double pzcc = pstt;                          649    G4double pzcc = pstt;
705    G4double eccm = stot - ( mass_proj + mass_t    650    G4double eccm = stot - ( mass_proj + mass_targ );
706                                                   651    
707    G4int zp = 1;                                  652    G4int zp = 1;
708    G4int ap = 1;                                  653    G4int ap = 1;
709    if ( pd_proj->GetParticleType() == "nucleus    654    if ( pd_proj->GetParticleType() == "nucleus" )
710    {                                              655    {
711       zp = pd_proj->GetAtomicNumber();            656       zp = pd_proj->GetAtomicNumber();
712       ap = pd_proj->GetAtomicMass();              657       ap = pd_proj->GetAtomicMass();
713    }                                              658    }
714    else                                           659    else 
715    {                                              660    {
716       // proton, neutron, mesons                  661       // proton, neutron, mesons
717       zp = int ( pd_proj->GetPDGCharge()/eplus    662       zp = int ( pd_proj->GetPDGCharge()/eplus + 0.5 );  
718       // ap = 1;                                  663       // ap = 1;
719    }                                              664    }
720                                                   665    
721                                                   666 
722    G4int zt = pd_targ->GetAtomicNumber();         667    G4int zt = pd_targ->GetAtomicNumber();
723    G4int at = pd_targ->GetAtomicMass();           668    G4int at = pd_targ->GetAtomicMass();
724                                                   669 
725                                                   670 
726    // Check the ramx0 value                    << 
727    //G4double rmax0 = 8.0;  // T.K dicide para    671    //G4double rmax0 = 8.0;  // T.K dicide parameter value  // for low energy
728    G4double rmax0 = bmax + 4.0;                   672    G4double rmax0 = bmax + 4.0;
729    G4double rmax = std::sqrt( rmax0*rmax0 + b*    673    G4double rmax = std::sqrt( rmax0*rmax0 + b*b );
730                                                   674 
731    G4double ccoul = 0.001439767;                  675    G4double ccoul = 0.001439767;
732    G4double pcca = 1.0 - double ( zp * zt ) *     676    G4double pcca = 1.0 - double ( zp * zt ) * ccoul / eccm / rmax - ( b / rmax )*( b / rmax );
733                                                   677 
734    G4double pccf = std::sqrt( pcca );             678    G4double pccf = std::sqrt( pcca );
735                                                   679 
736    //Fix for neutral particles                    680    //Fix for neutral particles
737    G4double aas1 = 0.0;                           681    G4double aas1 = 0.0;
738    G4double bbs = 0.0;                            682    G4double bbs = 0.0;
739                                                   683 
740    if ( zp != 0 )                                 684    if ( zp != 0 )
741    {                                              685    {
742       G4double aas = 2.0 * eccm * b / double (    686       G4double aas = 2.0 * eccm * b / double ( zp * zt ) / ccoul;
743       bbs = 1.0 / std::sqrt ( 1.0 + aas*aas );    687       bbs = 1.0 / std::sqrt ( 1.0 + aas*aas );
744       aas1 = ( 1.0 + aas * b / rmax ) * bbs;      688       aas1 = ( 1.0 + aas * b / rmax ) * bbs;
745    }                                              689    }
746                                                   690 
747    G4double cost = 0.0;                           691    G4double cost = 0.0;
748    G4double sint = 0.0;                           692    G4double sint = 0.0;
749    G4double thet1 = 0.0;                          693    G4double thet1 = 0.0;
750    G4double thet2 = 0.0;                          694    G4double thet2 = 0.0;
751    if ( 1.0 - aas1*aas1 <= 0 || 1.0 - bbs*bbs     695    if ( 1.0 - aas1*aas1 <= 0 || 1.0 - bbs*bbs <= 0.0 )   
752    {                                              696    {
753       cost = 1.0;                                 697       cost = 1.0;
754       sint = 0.0;                                 698       sint = 0.0;
755    }                                              699    } 
756    else                                           700    else 
757    {                                              701    {
758       G4double aat1 = aas1 / std::sqrt ( 1.0 -    702       G4double aat1 = aas1 / std::sqrt ( 1.0 - aas1*aas1 );
759       G4double aat2 = bbs / std::sqrt ( 1.0 -     703       G4double aat2 = bbs / std::sqrt ( 1.0 - bbs*bbs );
760                                                   704 
761       thet1 = std::atan ( aat1 );                 705       thet1 = std::atan ( aat1 );
762       thet2 = std::atan ( aat2 );                 706       thet2 = std::atan ( aat2 );
763                                                   707 
764 //    TK enter to else block                      708 //    TK enter to else block  
765       G4double theta = thet1 - thet2;             709       G4double theta = thet1 - thet2;
766       cost = std::cos( theta );                   710       cost = std::cos( theta );
767       sint = std::sin( theta );                   711       sint = std::sin( theta );
768    }                                              712    }
769                                                   713 
770    G4double rzpr = -rmax * cost * ( mass_targ     714    G4double rzpr = -rmax * cost * ( mass_targ ) / ( mass_proj + mass_targ );
771    G4double rzta =  rmax * cost * ( mass_proj     715    G4double rzta =  rmax * cost * ( mass_proj ) / ( mass_proj + mass_targ );
772                                                   716 
773    G4double rxpr = rmax / 2.0 * sint;             717    G4double rxpr = rmax / 2.0 * sint;
774                                                   718 
775    G4double rxta = -rxpr;                         719    G4double rxta = -rxpr;
776                                                   720 
777                                                   721 
778    G4double pzpc = pzcc * (  cost * pccf + sin    722    G4double pzpc = pzcc * (  cost * pccf + sint * b / rmax ); 
779    G4double pxpr = pzcc * ( -sint * pccf + cos    723    G4double pxpr = pzcc * ( -sint * pccf + cost * b / rmax ); 
780                                                   724 
781    G4double pztc = - pzpc;                        725    G4double pztc = - pzpc;
782    G4double pxta = - pxpr;                        726    G4double pxta = - pxpr;
783                                                   727 
784    G4double epc = std::sqrt ( pzpc*pzpc + pxpr    728    G4double epc = std::sqrt ( pzpc*pzpc + pxpr*pxpr + mass_proj*mass_proj );
785    G4double etc = std::sqrt ( pztc*pztc + pxta    729    G4double etc = std::sqrt ( pztc*pztc + pxta*pxta + mass_targ*mass_targ );
786                                                   730 
787    G4double pzpr = pzpc;                          731    G4double pzpr = pzpc;
788    G4double pzta = pztc;                          732    G4double pzta = pztc;
789    G4double epr = epc;                            733    G4double epr = epc;
790    G4double eta = etc;                            734    G4double eta = etc;
791                                                   735 
792 // CM -> NN                                       736 // CM -> NN
793    G4double gammacm = boostToCM.gamma();          737    G4double gammacm = boostToCM.gamma();
794    //G4double betacm = -boostToCM.beta();         738    //G4double betacm = -boostToCM.beta();
795    G4double betacm = boostToCM.z();               739    G4double betacm = boostToCM.z();
796    pzpr = pzpc + betacm * gammacm * ( gammacm     740    pzpr = pzpc + betacm * gammacm * ( gammacm / ( 1. + gammacm ) * pzpc * betacm + epc );
797    pzta = pztc + betacm * gammacm * ( gammacm     741    pzta = pztc + betacm * gammacm * ( gammacm / ( 1. + gammacm ) * pztc * betacm + etc );
798    epr = gammacm * ( epc + betacm * pzpc );       742    epr = gammacm * ( epc + betacm * pzpc );
799    eta = gammacm * ( etc + betacm * pztc );       743    eta = gammacm * ( etc + betacm * pztc );
800                                                   744 
801    //G4double betpr = pzpr / epr;                 745    //G4double betpr = pzpr / epr;
802    //G4double betta = pzta / eta;                 746    //G4double betta = pzta / eta;
803                                                   747 
804    G4double gammpr = epr / ( mass_proj );         748    G4double gammpr = epr / ( mass_proj );
805    G4double gammta = eta / ( mass_targ );         749    G4double gammta = eta / ( mass_targ );
806                                                   750       
807    pzta = pzta / double ( at );                   751    pzta = pzta / double ( at );
808    pxta = pxta / double ( at );                   752    pxta = pxta / double ( at );
809                                                   753       
810    pzpr = pzpr / double ( ap );                   754    pzpr = pzpr / double ( ap );
811    pxpr = pxpr / double ( ap );                   755    pxpr = pxpr / double ( ap );
812                                                   756 
813    G4double zeroz = 0.0;                          757    G4double zeroz = 0.0; 
814                                                   758 
815    rzpr = rzpr -zeroz;                            759    rzpr = rzpr -zeroz;
816    rzta = rzta -zeroz;                            760    rzta = rzta -zeroz;
817                                                   761 
818    // Set results                                 762    // Set results 
819    coulomb_collision_gamma_proj = gammpr;         763    coulomb_collision_gamma_proj = gammpr;
820    coulomb_collision_rx_proj = rxpr;              764    coulomb_collision_rx_proj = rxpr;
821    coulomb_collision_rz_proj = rzpr;              765    coulomb_collision_rz_proj = rzpr;
822    coulomb_collision_px_proj = pxpr;              766    coulomb_collision_px_proj = pxpr;
823    coulomb_collision_pz_proj = pzpr;              767    coulomb_collision_pz_proj = pzpr;
824                                                   768 
825    coulomb_collision_gamma_targ = gammta;         769    coulomb_collision_gamma_targ = gammta;
826    coulomb_collision_rx_targ = rxta;              770    coulomb_collision_rx_targ = rxta;
827    coulomb_collision_rz_targ = rzta;              771    coulomb_collision_rz_targ = rzta;
828    coulomb_collision_px_targ = pxta;              772    coulomb_collision_px_targ = pxta;
829    coulomb_collision_pz_targ = pzta;              773    coulomb_collision_pz_targ = pzta;
830                                                   774 
831 }                                                 775 }
832                                                   776 
                                                   >> 777 
                                                   >> 778 
833 void G4QMDReaction::setEvaporationCh()            779 void G4QMDReaction::setEvaporationCh()
834 {                                                 780 {
835   //fEvaporation - 8 default channels          << 
836   //fCombined    - 8 default + 60 GEM          << 
837   //fGEM         - 2 default + 66 GEM          << 
838   G4DeexChannelType ctype = gem ? fGEM : fComb << 
839   excitationHandler->SetDeexChannelsType(ctype << 
840 }                                              << 
841                                                   781 
842 void G4QMDReaction::ModelDescription(std::ostr << 782    if ( gem == true ) 
843 {                                              << 783       evaporation->SetGEMChannel();
844    outFile << "Lorentz covarianted Quantum Mol << 784    else
                                                   >> 785       evaporation->SetDefaultChannel();
                                                   >> 786 
845 }                                                 787 }
846                                                   788