Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/binary_cascade/src/G4BinaryLightIonReaction.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/binary_cascade/src/G4BinaryLightIonReaction.cc (Version 11.3.0) and /processes/hadronic/models/binary_cascade/src/G4BinaryLightIonReaction.cc (Version 10.3)


  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 #include <algorithm>                               26 #include <algorithm>
 27 #include <vector>                                  27 #include <vector>
 28 #include <cmath>                                   28 #include <cmath>
 29 #include <numeric>                                 29 #include <numeric>
 30                                                    30 
 31 #include "G4BinaryLightIonReaction.hh"             31 #include "G4BinaryLightIonReaction.hh"
 32 #include "G4PhysicalConstants.hh"                  32 #include "G4PhysicalConstants.hh"
 33 #include "G4SystemOfUnits.hh"                      33 #include "G4SystemOfUnits.hh"
 34 #include "G4LorentzVector.hh"                      34 #include "G4LorentzVector.hh"
 35 #include "G4LorentzRotation.hh"                    35 #include "G4LorentzRotation.hh"
 36 #include "G4ReactionProductVector.hh"              36 #include "G4ReactionProductVector.hh"
 37 #include "G4ping.hh"                               37 #include "G4ping.hh"
 38 #include "G4Delete.hh"                             38 #include "G4Delete.hh"
 39 #include "G4Neutron.hh"                            39 #include "G4Neutron.hh"
 40 #include "G4VNuclearDensity.hh"                    40 #include "G4VNuclearDensity.hh"
 41 #include "G4FermiMomentum.hh"                      41 #include "G4FermiMomentum.hh"
 42 #include "G4HadTmpUtil.hh"                         42 #include "G4HadTmpUtil.hh"
 43 #include "G4PreCompoundModel.hh"                   43 #include "G4PreCompoundModel.hh"
 44 #include "G4HadronicInteractionRegistry.hh"        44 #include "G4HadronicInteractionRegistry.hh"
 45 #include "G4Log.hh"                                45 #include "G4Log.hh"
 46 #include "G4PhysicsModelCatalog.hh"            << 
 47 #include "G4HadronicParameters.hh"             << 
 48                                                << 
 49 G4int G4BinaryLightIonReaction::theBLIR_ID = - << 
 50                                                    46 
 51 //#define debug_G4BinaryLightIonReaction           47 //#define debug_G4BinaryLightIonReaction
 52 //#define debug_BLIR_finalstate                    48 //#define debug_BLIR_finalstate
 53 //#define debug_BLIR_result                        49 //#define debug_BLIR_result
 54                                                    50 
 55 G4BinaryLightIonReaction::G4BinaryLightIonReac     51 G4BinaryLightIonReaction::G4BinaryLightIonReaction(G4VPreCompoundModel* ptr)
 56 : G4HadronicInteraction("Binary Light Ion Casc     52 : G4HadronicInteraction("Binary Light Ion Cascade"),
 57   theProjectileFragmentation(ptr),                 53   theProjectileFragmentation(ptr),
 58   pA(0),pZ(0), tA(0),tZ(0),spectatorA(0),spect     54   pA(0),pZ(0), tA(0),tZ(0),spectatorA(0),spectatorZ(0),
 59   projectile3dNucleus(0),target3dNucleus(0)        55   projectile3dNucleus(0),target3dNucleus(0)
 60 {                                                  56 {
 61   if(!ptr) {                                       57   if(!ptr) {
 62     G4HadronicInteraction* p =                     58     G4HadronicInteraction* p =
 63       G4HadronicInteractionRegistry::Instance(     59       G4HadronicInteractionRegistry::Instance()->FindModel("PRECO");
 64     G4VPreCompoundModel* pre = static_cast<G4V     60     G4VPreCompoundModel* pre = static_cast<G4VPreCompoundModel*>(p);
 65     if(!pre) { pre = new G4PreCompoundModel();     61     if(!pre) { pre = new G4PreCompoundModel(); }
 66     theProjectileFragmentation = pre;              62     theProjectileFragmentation = pre;
 67   }                                                63   }
 68   theModel = new G4BinaryCascade(theProjectile     64   theModel = new G4BinaryCascade(theProjectileFragmentation);
 69   theHandler = theProjectileFragmentation->Get     65   theHandler = theProjectileFragmentation->GetExcitationHandler();
 70       theBLIR_ID = G4PhysicsModelCatalog::GetM <<  66 
 71   debug_G4BinaryLightIonReactionResults = G4Ha <<  67   debug_G4BinaryLightIonReactionResults=getenv("debug_G4BinaryLightIonReactionResults")!=0;
 72 }                                                  68 }
 73                                                    69 
 74 G4BinaryLightIonReaction::~G4BinaryLightIonRea     70 G4BinaryLightIonReaction::~G4BinaryLightIonReaction()
 75 {}                                                 71 {}
 76                                                    72 
 77 void G4BinaryLightIonReaction::ModelDescriptio     73 void G4BinaryLightIonReaction::ModelDescription(std::ostream& outFile) const
 78 {                                                  74 {
 79   outFile << "G4Binary Light Ion Cascade is an     75   outFile << "G4Binary Light Ion Cascade is an intra-nuclear cascade model\n"
 80       << "using G4BinaryCasacde to model the i     76       << "using G4BinaryCasacde to model the interaction of a light\n"
 81       << "nucleus with a nucleus.\n"               77       << "nucleus with a nucleus.\n"
 82       << "The lighter of the two nuclei is tre     78       << "The lighter of the two nuclei is treated like a set of projectiles\n"
 83       << "which are transported simultaneously <<  79       << "which are transported simultanously through the heavier nucleus.\n";
 84 }                                                  80 }
 85                                                    81 
 86 //--------------------------------------------     82 //--------------------------------------------------------------------------------
 87 struct ReactionProduct4Mom                         83 struct ReactionProduct4Mom
 88 {                                                  84 {
 89    G4LorentzVector operator()(G4LorentzVector      85    G4LorentzVector operator()(G4LorentzVector a,G4ReactionProduct* b) {return a + G4LorentzVector(b->GetMomentum(), b->GetTotalEnergy() );}
 90 };                                                 86 };
 91                                                    87 
 92 G4HadFinalState *G4BinaryLightIonReaction::        88 G4HadFinalState *G4BinaryLightIonReaction::
 93 ApplyYourself(const G4HadProjectile &aTrack, G     89 ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus & targetNucleus )
 94 {                                                  90 {
 95   if(debug_G4BinaryLightIonReactionResults) G4 <<  91   if(getenv("BLICDEBUG") ) G4cerr << " ######### Binary Light Ion Reaction starts ######### " << G4endl;
 96   G4ping debug("debug_G4BinaryLightIonReaction     92   G4ping debug("debug_G4BinaryLightIonReaction");
 97   pA=aTrack.GetDefinition()->GetBaryonNumber()     93   pA=aTrack.GetDefinition()->GetBaryonNumber();
 98   pZ=G4lrint(aTrack.GetDefinition()->GetPDGCha     94   pZ=G4lrint(aTrack.GetDefinition()->GetPDGCharge()/eplus);
 99   tA=targetNucleus.GetA_asInt();                   95   tA=targetNucleus.GetA_asInt();
100   tZ=targetNucleus.GetZ_asInt();                   96   tZ=targetNucleus.GetZ_asInt();
101   G4double timePrimary = aTrack.GetGlobalTime( <<  97 
102   G4LorentzVector mom(aTrack.Get4Momentum());      98   G4LorentzVector mom(aTrack.Get4Momentum());
103    //G4cout << "proj mom : " << mom << G4endl;     99    //G4cout << "proj mom : " << mom << G4endl;
104   G4LorentzRotation toBreit(mom.boostVector())    100   G4LorentzRotation toBreit(mom.boostVector());
105                                                   101 
106   G4bool swapped=SetLighterAsProjectile(mom, t    102   G4bool swapped=SetLighterAsProjectile(mom, toBreit);
107    //G4cout << "after swap, swapped? / mom " <    103    //G4cout << "after swap, swapped? / mom " << swapped << " / " << mom <<G4endl;
108   G4ReactionProductVector * result = 0;           104   G4ReactionProductVector * result = 0;
109   G4ReactionProductVector * cascaders=0; //new    105   G4ReactionProductVector * cascaders=0; //new G4ReactionProductVector;
110 //  G4double m_nucl(0);      // to check energ    106 //  G4double m_nucl(0);      // to check energy balance
111                                                   107 
112                                                   108 
113   //    G4double m1=G4ParticleTable::GetPartic    109   //    G4double m1=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(pZ,pA);
114   //    G4cout << "Entering the decision point    110   //    G4cout << "Entering the decision point "
115   //           << (mom.t()-mom.mag())/pA << "     111   //           << (mom.t()-mom.mag())/pA << " "
116   //     << pA<<" "<< pZ<<" "                     112   //     << pA<<" "<< pZ<<" "
117   //     << tA<<" "<< tZ<<G4endl                  113   //     << tA<<" "<< tZ<<G4endl
118   //     << " "<<mom.t()-mom.mag()<<" "           114   //     << " "<<mom.t()-mom.mag()<<" "
119   //     << mom.t()- m1<<G4endl;                  115   //     << mom.t()- m1<<G4endl;
120   if( (mom.t()-mom.mag())/pA < 50*MeV )           116   if( (mom.t()-mom.mag())/pA < 50*MeV )
121   {                                               117   {
122     //      G4cout << "Using pre-compound only    118     //      G4cout << "Using pre-compound only, E= "<<mom.t()-mom.mag()<<G4endl;
123     //      m_nucl = mom.mag();                   119     //      m_nucl = mom.mag();
124      cascaders=FuseNucleiAndPrompound(mom);    << 120       cascaders=FuseNucleiAndPrompound(mom);
125      if( !cascaders )                          << 121       if( !cascaders )
126      {                                         << 122       {
127                                                << 123 
128               // abort!! happens for too low e << 124          // abort!! happens for too low energy for nuclei to fuse
129                                                << 125 
130               theResult.Clear();               << 126          theResult.Clear();
131               theResult.SetStatusChange(isAliv << 127          theResult.SetStatusChange(isAlive);
132               theResult.SetEnergyChange(aTrack << 128          theResult.SetEnergyChange(aTrack.GetKineticEnergy());
133               theResult.SetMomentumChange(aTra << 129          theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
134               return &theResult;               << 130          return &theResult;
135      }                                         << 131       }
136   }                                               132   }
137   else                                            133   else
138   {                                               134   {
139      result=Interact(mom,toBreit);                135      result=Interact(mom,toBreit);
140                                                   136 
141      if(! result )                             << 137       if(! result )
142      {                                         << 138       {
143            // abort!!                          << 139             // abort!!
144                                                << 140 
145            G4cerr << "G4BinaryLightIonReaction << 141             G4cerr << "G4BinaryLightIonReaction no final state for: " << G4endl;
146            G4cerr << " Primary " << aTrack.Get << 142             G4cerr << " Primary " << aTrack.GetDefinition()
147               << ", (A,Z)=(" << aTrack.GetDefi << 143                << ", (A,Z)=(" << aTrack.GetDefinition()->GetBaryonNumber()
148               << "," << aTrack.GetDefinition() << 144                << "," << aTrack.GetDefinition()->GetPDGCharge()/eplus << ") "
149               << ", kinetic energy " << aTrack << 145                << ", kinetic energy " << aTrack.GetKineticEnergy()
150               << G4endl;                       << 146                << G4endl;
151            G4cerr << " Target nucleus (A,Z)=(" << 147             G4cerr << " Target nucleus (A,Z)=("
152                   <<  (swapped?pA:tA)  << ","  << 148                    <<  (swapped?pA:tA)  << ","
153                   << (swapped?pZ:tZ) << ")" << << 149                    << (swapped?pZ:tZ) << ")" << G4endl;
154            G4cerr << " if frequent, please sub << 150             G4cerr << " if frequent, please submit above information as bug report"
155                        << G4endl << G4endl;    << 151                   << G4endl << G4endl;
156                                                << 152 
157            theResult.Clear();                  << 153             theResult.Clear();
158            theResult.SetStatusChange(isAlive); << 154             theResult.SetStatusChange(isAlive);
159            theResult.SetEnergyChange(aTrack.Ge << 155             theResult.SetEnergyChange(aTrack.GetKineticEnergy());
160            theResult.SetMomentumChange(aTrack. << 156             theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
161            return &theResult;                  << 157             return &theResult;
162      }                                         << 158       }
163                                                << 159 
164          // Calculate excitation energy,       << 160     // Calculate excitation energy,
165      G4double theStatisticalExEnergy = GetProj << 161       G4double theStatisticalExEnergy = GetProjectileExcitation();
166                                                << 
167                                                << 
168      pInitialState = mom;                      << 
169         //G4cout << "BLIC: pInitialState from  << 
170      pInitialState.setT(pInitialState.getT() + << 
171     G4ParticleTable::GetParticleTable()->GetIo << 
172       //G4cout << "BLIC: target nucleus added  << 
173                                                << 
174      delete target3dNucleus;target3dNucleus=0; << 
175      delete projectile3dNucleus;projectile3dNu << 
176                                                << 
177      G4ReactionProductVector * spectators= new << 
178                                                << 
179      cascaders = new G4ReactionProductVector;  << 
180                                                   162 
181      G4LorentzVector pspectators=SortResult(re << 
182              // this also sets spectatorA and  << 
183                                                   163 
184      //      pFinalState=std::accumulate(casca << 164     pInitialState = mom;
                                                   >> 165        //G4cout << "BLIC: pInitialState from aTrack : " << pInitialState;
                                                   >> 166     pInitialState.setT(pInitialState.getT() +
                                                   >> 167          G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(tZ,tA));
                                                   >> 168            //G4cout << "BLIC: target nucleus added : " << pInitialState << G4endl;
185                                                   169 
186      std::vector<G4ReactionProduct *>::iterato << 170     delete target3dNucleus;target3dNucleus=0;
                                                   >> 171     delete projectile3dNucleus;projectile3dNucleus=0;
187                                                   172 
188              // G4cout << "pInitialState, pFin << 173       G4ReactionProductVector * spectators= new G4ReactionProductVector;
                                                   >> 174 
                                                   >> 175       cascaders = new G4ReactionProductVector;
                                                   >> 176 
                                                   >> 177       G4LorentzVector pspectators=SortResult(result,spectators,cascaders);
                                                   >> 178 
                                                   >> 179       //      pFinalState=std::accumulate(cascaders->begin(),cascaders->end(),pFinalState,ReactionProduct4Mom);
                                                   >> 180 
                                                   >> 181       std::vector<G4ReactionProduct *>::iterator iter;
                                                   >> 182 
                                                   >> 183         // G4cout << "pInitialState, pFinalState / pspectators"<< pInitialState << " / " << pFinalState << " / " << pspectators << G4endl;
189     //      if ( spectA-spectatorA !=0 || spec    184     //      if ( spectA-spectatorA !=0 || spectZ-spectatorZ !=0)
190     //      {                                     185     //      {
191     //          G4cout << "spect Nucl != spect    186     //          G4cout << "spect Nucl != spectators: nucl a,z; spect a,z" <<
192     //        spectatorA <<" "<< spectatorZ <<    187     //        spectatorA <<" "<< spectatorZ <<" ; " << spectA <<" "<< spectZ << G4endl;
193     //      }                                     188     //      }
194      delete result;                            << 189     delete result;
195      result=0;                                 << 190     result=0;
196      G4LorentzVector momentum(pInitialState-pF << 191     G4LorentzVector momentum(pInitialState-pFinalState);
197      G4int loopcount(0);                       << 192     G4int loopcount(0);
198         //G4cout << "BLIC: momentum, pspectato << 193        //G4cout << "BLIC: momentum, pspectators : " << momentum << " / " << pspectators << G4endl;
199      while (std::abs(momentum.e()-pspectators. << 194     while (std::abs(momentum.e()-pspectators.e()) > 10*MeV)                /* Loop checking, 31.08.2015, G.Folger */
200                                                << 195                                                                            // see if on loopcount 
201      {                                         << 
202        G4LorentzVector pCorrect(pInitialState- << 
203         //G4cout << "BLIC:: BIC nonconservatio << 
204        // Correct outgoing casacde particles.. << 
205        G4bool EnergyIsCorrect=EnergyAndMomentu << 
206        if ( ! EnergyIsCorrect && debug_G4Binar << 
207        {                                       << 
208          G4cout << "Warning - G4BinaryLightIon << 
209        }                                       << 
210        pFinalState=G4LorentzVector(0,0,0,0);   << 
211        for(iter=cascaders->begin(); iter!=casc << 
212        {                                       << 
213          pFinalState += G4LorentzVector( (*ite << 
214        }                                       << 
215        momentum=pInitialState-pFinalState;     << 
216        if (++loopcount > 10 )                  << 
217        {                                       << 
218            break;                              << 
219        }                                       << 
220      }                                         << 
221                                                << 
222 //      Check if Energy/Momemtum is now ok, if << 
223      if ( std::abs(momentum.e()-pspectators.e( << 
224      {                                         << 
225     for (iter=spectators->begin();iter!=specta << 
226     {                                             196     {
227         delete *iter;                          << 197       G4LorentzVector pCorrect(pInitialState-pspectators);
                                                   >> 198              //G4cout << "BLIC:: BIC nonconservation? (pInitialState-pFinalState) / spectators :" << momentum << " / " << pspectators << "pCorrect "<< pCorrect<< G4endl;
                                                   >> 199       // Correct outgoing casacde particles.... to have momentum of (initial state - spectators)
                                                   >> 200       G4bool EnergyIsCorrect=EnergyAndMomentumCorrector(cascaders, pCorrect);
                                                   >> 201       if ( ! EnergyIsCorrect && debug_G4BinaryLightIonReactionResults)
                                                   >> 202       {
                                                   >> 203         G4cout << "Warning - G4BinaryLightIonReaction E/P correction for cascaders failed" << G4endl;
                                                   >> 204       }
                                                   >> 205       pFinalState=G4LorentzVector(0,0,0,0);
                                                   >> 206         for(iter=cascaders->begin(); iter!=cascaders->end(); iter++)
                                                   >> 207       {
                                                   >> 208         pFinalState += G4LorentzVector( (*iter)->GetMomentum(), (*iter)->GetTotalEnergy() );
                                                   >> 209       }
                                                   >> 210       momentum=pInitialState-pFinalState;
                                                   >> 211       if (++loopcount > 10 )
                                                   >> 212       {
                                                   >> 213         if ( momentum.vect().mag() - momentum.e()> 10*keV  )
                                                   >> 214         {
                                                   >> 215           G4cerr << "G4BinaryLightIonReaction.cc: Cannot correct 4-momentum of cascade particles" << G4endl;
                                                   >> 216           throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCasacde::ApplyCollision()");
                                                   >> 217         } else {
                                                   >> 218           break;
                                                   >> 219         }
                                                   >> 220 
                                                   >> 221       }
228     }                                             222     }
229     delete spectators;                         << 223 
230      for(iter=cascaders->begin(); iter!=cascad << 224     if (spectatorA > 0 )
231      {                                         << 225     {
232          delete *iter;                         << 226        // check spectator momentum
233      }                                         << 227        if ( momentum.vect().mag() - momentum.e()> 10*keV )
234      delete cascaders;                         << 228        {
235                                                << 229 
236      G4cout << "G4BinaryLightIonReaction.cc: m << 230           for (iter=spectators->begin();iter!=spectators->end();iter++)
237            << " initial - final " << momentum  << 231           {
238            << " .. pInitialState/pFinalState/s << 232            delete *iter;
239            << pInitialState << G4endl          << 233           }
240            << pFinalState << G4endl            << 234           delete spectators;
241            << pspectators << G4endl            << 235           for(iter=cascaders->begin(); iter!=cascaders->end(); iter++)
242            << " .. A,Z " << spectatorA <<" "<< << 236           {
243      G4cout << "G4BinaryLightIonReaction inval << 237            delete *iter;
244      G4cout << " Primary " << aTrack.GetDefini << 238           }
245              << ", (A,Z)=(" << aTrack.GetDefin << 239           delete cascaders;
246              << "," << aTrack.GetDefinition()- << 240 
247              << ", kinetic energy " << aTrack. << 241           G4cout << "G4BinaryLightIonReaction.cc: mom check: " <<  momentum
248              << G4endl;                        << 242                 << " 3.mag "<< momentum.vect().mag() << G4endl
249      G4cout << " Target nucleus (A,Z)=(" <<  t << 243                 << " .. pInitialState/pFinalState/spectators " << pInitialState <<" "
250                  << "," << targetNucleus.GetZ_ << 244                 << pFinalState << " " << pspectators << G4endl
251      G4cout << " if frequent, please submit ab << 245                 << " .. A,Z " << spectatorA <<" "<< spectatorZ << G4endl;
252            << G4endl << G4endl;                << 246           G4cout << "G4BinaryLightIonReaction invalid final state for: " << G4endl;
                                                   >> 247           G4cout << " Primary " << aTrack.GetDefinition()
                                                   >> 248                   << ", (A,Z)=(" << aTrack.GetDefinition()->GetBaryonNumber()
                                                   >> 249                   << "," << aTrack.GetDefinition()->GetPDGCharge()/eplus << ") "
                                                   >> 250                   << ", kinetic energy " << aTrack.GetKineticEnergy()
                                                   >> 251                   << G4endl;
                                                   >> 252           G4cout << " Target nucleus (A,Z)=(" <<  targetNucleus.GetA_asInt()
                                                   >> 253                       << "," << targetNucleus.GetZ_asInt() << ")" << G4endl;
                                                   >> 254           G4cout << " if frequent, please submit above information as bug report"
                                                   >> 255                 << G4endl << G4endl;
253 #ifdef debug_G4BinaryLightIonReaction             256 #ifdef debug_G4BinaryLightIonReaction
254           G4ExceptionDescription ed;              257           G4ExceptionDescription ed;
255           ed << "G4BinaryLightIonreaction: Ter    258           ed << "G4BinaryLightIonreaction: Terminate for above error"  << G4endl;
256           G4Exception("G4BinaryLightIonreactio    259           G4Exception("G4BinaryLightIonreaction::ApplyYourSelf()", "BLIC001", FatalException,
257           ed);                                    260           ed);
258                                                   261 
259 #endif                                            262 #endif
260      theResult.Clear();                        << 263           theResult.Clear();
261      theResult.SetStatusChange(isAlive);       << 264           theResult.SetStatusChange(isAlive);
262      theResult.SetEnergyChange(aTrack.GetKinet << 265           theResult.SetEnergyChange(aTrack.GetKineticEnergy());
263      theResult.SetMomentumChange(aTrack.Get4Mo << 266           theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
264      return &theResult;                        << 267           return &theResult;
265                                                << 268        }
266      }                                         << 269 
267        if (spectatorA > 0 )                    << 270                    // DeExciteSpectatorNucleus() also handles also case of A=1, Z=0,1
268      {                                         << 271        DeExciteSpectatorNucleus(spectators, cascaders, theStatisticalExEnergy, momentum);
269            // DeExciteSpectatorNucleus() also  << 272     } else {
270                DeExciteSpectatorNucleus(specta << 273        delete spectators;
271      } else {    // no spectators              << 274     }
272          delete spectators;                    << 
273      }                                         << 
274   }                                               275   }
275   // Rotate to lab                                276   // Rotate to lab
276   G4LorentzRotation toZ;                          277   G4LorentzRotation toZ;
277   toZ.rotateZ(-1*mom.phi());                      278   toZ.rotateZ(-1*mom.phi());
278   toZ.rotateY(-1*mom.theta());                    279   toZ.rotateY(-1*mom.theta());
279   G4LorentzRotation toLab(toZ.inverse());         280   G4LorentzRotation toLab(toZ.inverse());
280                                                   281 
281   // Fill the particle change, while rotating.    282   // Fill the particle change, while rotating. Boost from projectile breit-frame in case we swapped.
282   // theResult.Clear();                           283   // theResult.Clear();
283   theResult.Clear();                              284   theResult.Clear();
284   theResult.SetStatusChange(stopAndKill);         285   theResult.SetStatusChange(stopAndKill);
285   G4LorentzVector ptot(0);                        286   G4LorentzVector ptot(0);
286   #ifdef debug_BLIR_result                     << 287     G4ReactionProductVector::iterator iter;
287      G4LorentzVector p_raw;                    << 288     #ifdef debug_BLIR_result
288   #endif                                       << 289        G4LorentzVector p_raw;
289   //G4int i=0;                                 << 290     #endif
                                                   >> 291     //G4int i=0;
290                                                   292 
291         G4ReactionProductVector::iterator iter << 
292   for(iter=cascaders->begin(); iter!=cascaders    293   for(iter=cascaders->begin(); iter!=cascaders->end(); iter++)
293   {                                               294   {
294     if((*iter)->GetNewlyAdded())                  295     if((*iter)->GetNewlyAdded())
295     {                                             296     {
296       G4DynamicParticle * aNewDP =             << 297       G4DynamicParticle * aNew =
297           new G4DynamicParticle((*iter)->GetDe    298           new G4DynamicParticle((*iter)->GetDefinition(),
298               (*iter)->GetTotalEnergy(),          299               (*iter)->GetTotalEnergy(),
299               (*iter)->GetMomentum() );           300               (*iter)->GetMomentum() );
300       G4LorentzVector tmp = aNewDP->Get4Moment << 301       G4LorentzVector tmp = aNew->Get4Momentum();
301              #ifdef debug_BLIR_result             302              #ifdef debug_BLIR_result
302            p_raw+= tmp;                           303            p_raw+= tmp;
303              #endif                               304              #endif
304       if(swapped)                                 305       if(swapped)
305       {                                           306       {
306         tmp*=toBreit.inverse();                   307         tmp*=toBreit.inverse();
307         tmp.setVect(-tmp.vect());                 308         tmp.setVect(-tmp.vect());
308       }                                           309       }
309       tmp *= toLab;                               310       tmp *= toLab;
310       aNewDP->Set4Momentum(tmp);               << 311       aNew->Set4Momentum(tmp);
311       G4HadSecondary aNew = G4HadSecondary(aNe << 
312             G4double time = 0;                 << 
313             //if(time < 0.0) { time = 0.0; }   << 
314             aNew.SetTime(timePrimary + time);  << 
315             //aNew.SetCreatorModelID((*iter)-> << 
316             aNew.SetCreatorModelID(theBLIR_ID) << 
317                                                << 
318       theResult.AddSecondary(aNew);               312       theResult.AddSecondary(aNew);
319       ptot += tmp;                                313       ptot += tmp;
320               //G4cout << "BLIC: Secondary " <    314               //G4cout << "BLIC: Secondary " << aNew->GetDefinition()->GetParticleName()
321               //       <<" "<<  aNew->GetMomen    315               //       <<" "<<  aNew->GetMomentum()<<" "<<  aNew->GetTotalEnergy() << G4endl;
322     }                                             316     }
323     delete *iter;                                 317     delete *iter;
324   }                                               318   }
325   delete cascaders;                               319   delete cascaders;
326                                                   320 
327 #ifdef debug_BLIR_result                          321 #ifdef debug_BLIR_result
328   //G4cout << "Result analysis, secondaries "     322   //G4cout << "Result analysis, secondaries " << theResult.GetNumberOfSecondaries() << G4endl;
329   //G4cout << "p_tot_raw " << p_raw << " sum p    323   //G4cout << "p_tot_raw " << p_raw << " sum p final " << ptot << G4endl;
330   G4double m_nucl=  G4ParticleTable::GetPartic    324   G4double m_nucl=  G4ParticleTable::GetParticleTable()->GetIonTable()->
331           GetIonMass(targetNucleus.GetZ_asInt(    325           GetIonMass(targetNucleus.GetZ_asInt(),targetNucleus.GetA_asInt());
332   // delete? tZ=targetNucleus.GetZ_asInt();       326   // delete? tZ=targetNucleus.GetZ_asInt();
333                                                   327 
334   //G4cout << "BLIC Energy conservation initia    328   //G4cout << "BLIC Energy conservation initial/primary/nucleus/final/delta(init-final) "
335    //     << aTrack.GetTotalEnergy()   + m_nuc    329    //     << aTrack.GetTotalEnergy()   + m_nucl <<" "<< aTrack.GetTotalEnergy() <<" "<< m_nucl <<" "<<ptot.e()
336    //     <<" "<< aTrack.GetTotalEnergy() + m_    330    //     <<" "<< aTrack.GetTotalEnergy() + m_nucl - ptot.e() << G4endl;
337   G4cout << "BLIC momentum conservation " << a    331   G4cout << "BLIC momentum conservation " << aTrack.Get4Momentum()+ G4LorentzVector(m_nucl)
338       << " ptot " << ptot << " delta " << aTra    332       << " ptot " << ptot << " delta " << aTrack.Get4Momentum()+ G4LorentzVector(m_nucl) - ptot
339       << "        3mom.mag() " << (aTrack.Get4    333       << "        3mom.mag() " << (aTrack.Get4Momentum()+ G4LorentzVector(m_nucl) - ptot).vect().mag() << G4endl;
340 #endif                                            334 #endif
341                                                   335 
342   if(debug_G4BinaryLightIonReactionResults) G4 << 336   if(getenv("BLICDEBUG") ) G4cerr << " ######### Binary Light Ion Reaction number ends ######### " << G4endl;
343                                                   337 
344   return &theResult;                              338   return &theResult;
345 }                                                 339 }
346                                                   340 
347 //--------------------------------------------    341 //--------------------------------------------------------------------------------
348                                                   342 
349 //********************************************    343 //****************************************************************************
350 G4bool G4BinaryLightIonReaction::EnergyAndMome    344 G4bool G4BinaryLightIonReaction::EnergyAndMomentumCorrector(
351     G4ReactionProductVector* Output, G4Lorentz    345     G4ReactionProductVector* Output, G4LorentzVector& TotalCollisionMom)
352 //********************************************    346 //****************************************************************************
353 {                                                 347 {
354   const int    nAttemptScale = 2500;              348   const int    nAttemptScale = 2500;
355   const double ErrLimit = 1.E-6;                  349   const double ErrLimit = 1.E-6;
356   if (Output->empty())                            350   if (Output->empty())
357     return TRUE;                                  351     return TRUE;
358   G4LorentzVector SumMom(0,0,0,0);                352   G4LorentzVector SumMom(0,0,0,0);
359   G4double        SumMass = 0;                    353   G4double        SumMass = 0;
360   G4double        TotalCollisionMass = TotalCo    354   G4double        TotalCollisionMass = TotalCollisionMom.m();
361   size_t i = 0;                                   355   size_t i = 0;
362   // Calculate sum hadron 4-momenta and summin    356   // Calculate sum hadron 4-momenta and summing hadron mass
363   for(i = 0; i < Output->size(); i++)             357   for(i = 0; i < Output->size(); i++)
364   {                                               358   {
365     SumMom  += G4LorentzVector((*Output)[i]->G    359     SumMom  += G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy());
366     SumMass += (*Output)[i]->GetDefinition()->    360     SumMass += (*Output)[i]->GetDefinition()->GetPDGMass();
367   }                                               361   }
368     // G4cout << " E/P corrector, SumMass, Sum    362     // G4cout << " E/P corrector, SumMass, SumMom.m2, TotalMass "
369     //       << SumMass <<" "<< SumMom.m2() <<    363     //       << SumMass <<" "<< SumMom.m2() <<" "<<TotalCollisionMass<< G4endl;
370   if (SumMass > TotalCollisionMass) return FAL    364   if (SumMass > TotalCollisionMass) return FALSE;
371   SumMass = SumMom.m2();                          365   SumMass = SumMom.m2();
372   if (SumMass < 0) return FALSE;                  366   if (SumMass < 0) return FALSE;
373   SumMass = std::sqrt(SumMass);                   367   SumMass = std::sqrt(SumMass);
374                                                   368 
375   // Compute c.m.s. hadron velocity and boost     369   // Compute c.m.s. hadron velocity and boost KTV to hadron c.m.s.
376   G4ThreeVector Beta = -SumMom.boostVector();     370   G4ThreeVector Beta = -SumMom.boostVector();
377         //G4cout << " == pre boost 2 "<< SumMo    371         //G4cout << " == pre boost 2 "<< SumMom.e()<< " "<< SumMom.mag()<<" "<< Beta <<G4endl;
378   //--old    Output->Boost(Beta);                 372   //--old    Output->Boost(Beta);
379   for(i = 0; i < Output->size(); i++)             373   for(i = 0; i < Output->size(); i++)
380   {                                               374   {
381     G4LorentzVector mom = G4LorentzVector((*Ou    375     G4LorentzVector mom = G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy());
382     mom *= Beta;                                  376     mom *= Beta;
383     (*Output)[i]->SetMomentum(mom.vect());        377     (*Output)[i]->SetMomentum(mom.vect());
384     (*Output)[i]->SetTotalEnergy(mom.e());        378     (*Output)[i]->SetTotalEnergy(mom.e());
385   }                                               379   }
386                                                   380 
387   // Scale total c.m.s. hadron energy (hadron     381   // Scale total c.m.s. hadron energy (hadron system mass).
388   // It should be equal interaction mass          382   // It should be equal interaction mass
389   G4double Scale = 0,OldScale=0;                  383   G4double Scale = 0,OldScale=0;
390   G4double factor = 1.;                           384   G4double factor = 1.;
391   G4int cAttempt = 0;                             385   G4int cAttempt = 0;
392   G4double Sum = 0;                               386   G4double Sum = 0;
393   G4bool success = false;                         387   G4bool success = false;
394   for(cAttempt = 0; cAttempt < nAttemptScale;     388   for(cAttempt = 0; cAttempt < nAttemptScale; cAttempt++)
395   {                                               389   {
396     Sum = 0;                                      390     Sum = 0;
397     for(i = 0; i < Output->size(); i++)           391     for(i = 0; i < Output->size(); i++)
398     {                                             392     {
399       G4LorentzVector HadronMom = G4LorentzVec    393       G4LorentzVector HadronMom = G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy());
400       HadronMom.setVect(HadronMom.vect()+ fact    394       HadronMom.setVect(HadronMom.vect()+ factor*Scale*HadronMom.vect());
401       G4double E = std::sqrt(HadronMom.vect().    395       G4double E = std::sqrt(HadronMom.vect().mag2() + sqr((*Output)[i]->GetDefinition()->GetPDGMass()));
402       HadronMom.setE(E);                          396       HadronMom.setE(E);
403       (*Output)[i]->SetMomentum(HadronMom.vect    397       (*Output)[i]->SetMomentum(HadronMom.vect());
404       (*Output)[i]->SetTotalEnergy(HadronMom.e    398       (*Output)[i]->SetTotalEnergy(HadronMom.e());
405       Sum += E;                                   399       Sum += E;
406     }                                             400     }
407     OldScale=Scale;                               401     OldScale=Scale;
408     Scale = TotalCollisionMass/Sum - 1;           402     Scale = TotalCollisionMass/Sum - 1;
409     //  G4cout << "E/P corr - " << cAttempt <<    403     //  G4cout << "E/P corr - " << cAttempt << " " << Scale << G4endl;
410     if (std::abs(Scale) <= ErrLimit               404     if (std::abs(Scale) <= ErrLimit
411         || OldScale == Scale)     // protect '    405         || OldScale == Scale)     // protect 'frozen' situation and divide by 0 in calculating new factor below
412     {                                             406     {
413       if (debug_G4BinaryLightIonReactionResult    407       if (debug_G4BinaryLightIonReactionResults) G4cout << "E/p corrector: " << cAttempt << G4endl;
414       success = true;                             408       success = true;
415       break;                                      409       break;
416     }                                             410     }
417     if ( cAttempt > 10 )                          411     if ( cAttempt > 10 )
418     {                                             412     {
419       //         G4cout << " speed it up? " <<    413       //         G4cout << " speed it up? " << std::abs(OldScale/(OldScale-Scale)) << G4endl;
420       factor=std::max(1.,G4Log(std::abs(OldSca    414       factor=std::max(1.,G4Log(std::abs(OldScale/(OldScale-Scale))));
421       //   G4cout << " ? factor ? " << factor     415       //   G4cout << " ? factor ? " << factor << G4endl;
422     }                                             416     }
423   }                                               417   }
424                                                   418 
425   if( (!success)  && debug_G4BinaryLightIonRea    419   if( (!success)  && debug_G4BinaryLightIonReactionResults)
426   {                                               420   {
427     G4cout << "G4G4BinaryLightIonReaction::Ene    421     G4cout << "G4G4BinaryLightIonReaction::EnergyAndMomentumCorrector - Warning"<<G4endl;
428     G4cout << "   Scale not unity at end of it    422     G4cout << "   Scale not unity at end of iteration loop: "<<TotalCollisionMass<<" "<<Sum<<" "<<Scale<<G4endl;
429     G4cout << "   Increase number of attempts     423     G4cout << "   Increase number of attempts or increase ERRLIMIT"<<G4endl;
430   }                                               424   }
431                                                   425 
432   // Compute c.m.s. interaction velocity and K    426   // Compute c.m.s. interaction velocity and KTV back boost
433   Beta = TotalCollisionMom.boostVector();         427   Beta = TotalCollisionMom.boostVector();
434   //--old    Output->Boost(Beta);                 428   //--old    Output->Boost(Beta);
435   for(i = 0; i < Output->size(); i++)             429   for(i = 0; i < Output->size(); i++)
436   {                                               430   {
437     G4LorentzVector mom = G4LorentzVector((*Ou    431     G4LorentzVector mom = G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy());
438     mom *= Beta;                                  432     mom *= Beta;
439     (*Output)[i]->SetMomentum(mom.vect());        433     (*Output)[i]->SetMomentum(mom.vect());
440     (*Output)[i]->SetTotalEnergy(mom.e());        434     (*Output)[i]->SetTotalEnergy(mom.e());
441   }                                               435   }
442   return TRUE;                                    436   return TRUE;
443 }                                                 437 }
444 G4bool G4BinaryLightIonReaction::SetLighterAsP    438 G4bool G4BinaryLightIonReaction::SetLighterAsProjectile(G4LorentzVector & mom,const G4LorentzRotation & toBreit)
445 {                                                 439 {
446    G4bool swapped = false;                        440    G4bool swapped = false;
447    if(tA<pA)                                      441    if(tA<pA)
448    {                                              442    {
449       swapped = true;                             443       swapped = true;
450       G4int tmp(0);                               444       G4int tmp(0);
451       tmp = tA; tA=pA; pA=tmp;                    445       tmp = tA; tA=pA; pA=tmp;
452       tmp = tZ; tZ=pZ; pZ=tmp;                    446       tmp = tZ; tZ=pZ; pZ=tmp;
453       G4double m1=G4ParticleTable::GetParticle    447       G4double m1=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(pZ,pA);
454       G4LorentzVector it(m1, G4ThreeVector(0,0    448       G4LorentzVector it(m1, G4ThreeVector(0,0,0));
455       mom = toBreit*it;                           449       mom = toBreit*it;
456    }                                              450    }
457    return swapped;                                451    return swapped;
458 }                                                 452 }
459 G4ReactionProductVector * G4BinaryLightIonReac    453 G4ReactionProductVector * G4BinaryLightIonReaction::FuseNucleiAndPrompound(const G4LorentzVector & mom)
460 {                                                 454 {
461    // Check if kinematically nuclei can fuse.     455    // Check if kinematically nuclei can fuse.
462    G4double mFused=G4ParticleTable::GetParticl    456    G4double mFused=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(pZ+tZ,pA+tA);
463    G4double mTarget=G4ParticleTable::GetPartic    457    G4double mTarget=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(tZ,tA);
464    G4LorentzVector pCompound(mom.e()+mTarget,m    458    G4LorentzVector pCompound(mom.e()+mTarget,mom.vect());
465    G4double m2Compound=pCompound.m2();            459    G4double m2Compound=pCompound.m2();
466    if (m2Compound < sqr(mFused) ) {               460    if (m2Compound < sqr(mFused) ) {
467      //G4cout << "G4BLIC: projectile p, mTarge    461      //G4cout << "G4BLIC: projectile p, mTarget, mFused, mCompound, delta: " <<mom <<  " " << mTarget <<  " " << mFused
468      //    <<  " " << sqrt(m2Compound)<<  " "     462      //    <<  " " << sqrt(m2Compound)<<  " " << sqrt(m2Compound) - mFused << G4endl;
469      return 0;                                    463      return 0;
470    }                                              464    }
471                                                   465 
472    G4Fragment aPreFrag;                           466    G4Fragment aPreFrag;
473    aPreFrag.SetZandA_asInt(pZ+tZ, pA+tA);         467    aPreFrag.SetZandA_asInt(pZ+tZ, pA+tA);
474    aPreFrag.SetNumberOfParticles(pA);             468    aPreFrag.SetNumberOfParticles(pA);
475    aPreFrag.SetNumberOfCharged(pZ);               469    aPreFrag.SetNumberOfCharged(pZ);
476    aPreFrag.SetNumberOfHoles(0);                  470    aPreFrag.SetNumberOfHoles(0);
477    //GF FIXME: whyusing plop in z direction? t    471    //GF FIXME: whyusing plop in z direction? this will not conserve momentum?
478    //G4ThreeVector plop(0.,0., mom.vect().mag(    472    //G4ThreeVector plop(0.,0., mom.vect().mag());
479    //G4LorentzVector aL(mom.t()+mTarget, plop)    473    //G4LorentzVector aL(mom.t()+mTarget, plop);
480    G4LorentzVector aL(mom.t()+mTarget,mom.vect    474    G4LorentzVector aL(mom.t()+mTarget,mom.vect());
481    aPreFrag.SetMomentum(aL);                      475    aPreFrag.SetMomentum(aL);
482                                                   476 
483                                                   477 
484          //G4cout << "Fragment INFO "<< pA+tA     478          //G4cout << "Fragment INFO "<< pA+tA <<" "<<pZ+tZ<<" "
485          //       << aL <<" "<<G4endl << aPreF    479          //       << aL <<" "<<G4endl << aPreFrag << G4endl;
486    G4ReactionProductVector * cascaders = thePr    480    G4ReactionProductVector * cascaders = theProjectileFragmentation->DeExcite(aPreFrag);
487    //G4double tSum = 0;                           481    //G4double tSum = 0;
488    for(size_t count = 0; count<cascaders->size    482    for(size_t count = 0; count<cascaders->size(); count++)
489    {                                              483    {
490       cascaders->operator[](count)->SetNewlyAd    484       cascaders->operator[](count)->SetNewlyAdded(true);
491       //tSum += cascaders->operator[](count)->    485       //tSum += cascaders->operator[](count)->GetKineticEnergy();
492    }                                              486    }
493    //       G4cout << "Exiting pre-compound on    487    //       G4cout << "Exiting pre-compound only, E= "<<tSum<<G4endl;
494    return cascaders;                              488    return cascaders;
495 }                                                 489 }
496 G4ReactionProductVector * G4BinaryLightIonReac    490 G4ReactionProductVector * G4BinaryLightIonReaction::Interact(G4LorentzVector & mom, const G4LorentzRotation & toBreit)
497 {                                                 491 {
498       G4ReactionProductVector * result = 0;       492       G4ReactionProductVector * result = 0;
499       G4double projectileMass(0);                 493       G4double projectileMass(0);
500       G4LorentzVector it;                         494       G4LorentzVector it;
501                                                   495 
502       G4int tryCount(0);                          496       G4int tryCount(0);
503       do                                          497       do
504       {                                           498       {
505          ++tryCount;                              499          ++tryCount;
506          projectile3dNucleus = new G4Fancy3DNu    500          projectile3dNucleus = new G4Fancy3DNucleus;
507          projectile3dNucleus->Init(pA, pZ);       501          projectile3dNucleus->Init(pA, pZ);
508          projectile3dNucleus->CenterNucleons()    502          projectile3dNucleus->CenterNucleons();
509          projectileMass=G4ParticleTable::GetPa    503          projectileMass=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(
510                projectile3dNucleus->GetCharge(    504                projectile3dNucleus->GetCharge(),projectile3dNucleus->GetMassNumber());
511          it=toBreit * G4LorentzVector(projecti    505          it=toBreit * G4LorentzVector(projectileMass,G4ThreeVector(0,0,0));
512                                                   506 
513          target3dNucleus = new G4Fancy3DNucleu    507          target3dNucleus = new G4Fancy3DNucleus;
514          target3dNucleus->Init(tA, tZ);           508          target3dNucleus->Init(tA, tZ);
515          G4double impactMax = target3dNucleus-    509          G4double impactMax = target3dNucleus->GetOuterRadius()+projectile3dNucleus->GetOuterRadius();
516          //        G4cout << "out radius - nuc    510          //        G4cout << "out radius - nucleus - projectile " << target3dNucleus->GetOuterRadius()/fermi << " - " << projectile3dNucleus->GetOuterRadius()/fermi << G4endl;
517          G4double aX=(2.*G4UniformRand()-1.)*i    511          G4double aX=(2.*G4UniformRand()-1.)*impactMax;
518          G4double aY=(2.*G4UniformRand()-1.)*i    512          G4double aY=(2.*G4UniformRand()-1.)*impactMax;
519          G4ThreeVector pos(aX, aY, -2.*impactM    513          G4ThreeVector pos(aX, aY, -2.*impactMax-5.*fermi);
520                                                   514 
521          G4KineticTrackVector * initalState =     515          G4KineticTrackVector * initalState = new G4KineticTrackVector;
522          projectile3dNucleus->StartLoop();        516          projectile3dNucleus->StartLoop();
523          G4Nucleon * aNuc;                        517          G4Nucleon * aNuc;
524          G4LorentzVector tmpV(0,0,0,0);           518          G4LorentzVector tmpV(0,0,0,0);
525          #ifdef debug_BLIR_finalstate             519          #ifdef debug_BLIR_finalstate
526              G4LorentzVector pinitial;            520              G4LorentzVector pinitial;
527          #endif                                   521          #endif
528          G4LorentzVector nucleonMom(1./pA*mom)    522          G4LorentzVector nucleonMom(1./pA*mom);
529          nucleonMom.setZ(nucleonMom.vect().mag    523          nucleonMom.setZ(nucleonMom.vect().mag());
530          nucleonMom.setX(0);                      524          nucleonMom.setX(0);
531          nucleonMom.setY(0);                      525          nucleonMom.setY(0);
532          theFermi.Init(pA,pZ);                    526          theFermi.Init(pA,pZ);
533          while( (aNuc=projectile3dNucleus->Get    527          while( (aNuc=projectile3dNucleus->GetNextNucleon()) )   /* Loop checking, 31.08.2015, G.Folger */
534          {                                        528          {
535             G4LorentzVector p4 = aNuc->GetMome    529             G4LorentzVector p4 = aNuc->GetMomentum();
536             tmpV+=p4;                             530             tmpV+=p4;
537             G4ThreeVector nucleonPosition(aNuc    531             G4ThreeVector nucleonPosition(aNuc->GetPosition());
538             G4double density=(projectile3dNucl    532             G4double density=(projectile3dNucleus->GetNuclearDensity())->GetDensity(nucleonPosition);
539             nucleonPosition += pos;               533             nucleonPosition += pos;
540             G4KineticTrack * it1 = new G4Kinet    534             G4KineticTrack * it1 = new G4KineticTrack(aNuc, nucleonPosition, nucleonMom );
541             it1->SetState(G4KineticTrack::outs    535             it1->SetState(G4KineticTrack::outside);
542             G4double pfermi= theFermi.GetFermi    536             G4double pfermi= theFermi.GetFermiMomentum(density);
543             G4double mass = aNuc->GetDefinitio    537             G4double mass = aNuc->GetDefinition()->GetPDGMass();
544             G4double Efermi= std::sqrt( sqr(ma    538             G4double Efermi= std::sqrt( sqr(mass) + sqr(pfermi)) - mass;
545             it1->SetProjectilePotential(-Eferm    539             it1->SetProjectilePotential(-Efermi);
546             initalState->push_back(it1);          540             initalState->push_back(it1);
547             #ifdef debug_BLIR_finalstate          541             #ifdef debug_BLIR_finalstate
548                pinitial += it1->Get4Momentum()    542                pinitial += it1->Get4Momentum();
549             #endif                                543             #endif
550          }                                        544          }
551                                                   545 
552          result=theModel->Propagate(initalStat    546          result=theModel->Propagate(initalState, target3dNucleus);
553          #ifdef debug_BLIR_finalstate             547          #ifdef debug_BLIR_finalstate
554            if( result && result->size()>0)     << 548          if( result && result->size()>0)
555            {                                   << 549          {
556      G4cout << "  Cascade result " << G4endl;  << 550            G4LorentzVector presult;
557              G4LorentzVector presult;          << 551            G4ReactionProductVector::iterator iter;
558              G4ReactionProductVector::iterator << 552            G4ReactionProduct xp;
559              G4ReactionProduct xp;             << 553            for (iter=result->begin(); iter !=result->end(); ++iter)
560              for (iter=result->begin(); iter ! << 554            {
561              {                                 << 555             presult += G4LorentzVector((*iter)->GetMomentum(),(*iter)->GetTotalEnergy());
562               presult += G4LorentzVector((*ite << 556            }
563         G4cout << (*iter)->GetDefinition()->Ge << 557 
564         << "("<< (*iter)->GetMomentum().x()<<" << 558                 G4cout << "BLIC check result :  initial " << pinitial << " mass tgt " << target3dNucleus->GetMass()
565         <<    (*iter)->GetMomentum().y()<<","  << 559                << " final " << presult
566         <<    (*iter)->GetMomentum().z()<<";"  << 560                << " IF - FF " << pinitial +G4LorentzVector(target3dNucleus->GetMass()) - presult << G4endl;
567         <<    (*iter)->GetTotalEnergy() <<")"< << 
568              }                                 << 
569                                                << 
570             G4cout << "BLIC check result :  in << 
571                  << " final " << presult       << 
572                  << " IF - FF " << pinitial +G << 
573                                                   561 
574            }                                   << 562          }
575          #endif                                   563          #endif
576          if( result && result->size()==0)         564          if( result && result->size()==0)
577          {                                        565          {
578             delete result;                        566             delete result;
579             result=0;                             567             result=0;
580          }                                        568          }
581          if ( ! result )                          569          if ( ! result )
582          {                                        570          {
583             delete target3dNucleus;               571             delete target3dNucleus;
584             delete projectile3dNucleus;           572             delete projectile3dNucleus;
585          }                                        573          }
586                                                   574 
587          // std::for_each(initalState->begin()    575          // std::for_each(initalState->begin(), initalState->end(), Delete<G4KineticTrack>());
588          // delete initalState;                   576          // delete initalState;
589                                                   577 
590       } while (! result && tryCount< 150);   /    578       } while (! result && tryCount< 150);   /* Loop checking, 31.08.2015, G.Folger */
591       return result;                              579       return result;
592 }                                                 580 }
593 G4double G4BinaryLightIonReaction::GetProjecti    581 G4double G4BinaryLightIonReaction::GetProjectileExcitation()
594 {                                                 582 {
595                                                   583 
596       G4Nucleon * aNuc;                           584       G4Nucleon * aNuc;
597       // the projectileNucleus excitation ener    585       // the projectileNucleus excitation energy estimate...
598       G4double theStatisticalExEnergy = 0;        586       G4double theStatisticalExEnergy = 0;
599       projectile3dNucleus->StartLoop();           587       projectile3dNucleus->StartLoop();
600       while( (aNuc=projectile3dNucleus->GetNex    588       while( (aNuc=projectile3dNucleus->GetNextNucleon()) )   /* Loop checking, 31.08.2015, G.Folger */
601       {                                           589       {
602                 //G4cout << " Nucleon : " << a    590                 //G4cout << " Nucleon : " << aNuc->GetDefinition()->GetParticleName() <<" "<< aNuc->AreYouHit() <<" "<<aNuc->GetMomentum()<<G4endl;
603          if(aNuc->AreYouHit()) {                  591          if(aNuc->AreYouHit()) {
604             G4ThreeVector aPosition(aNuc->GetP    592             G4ThreeVector aPosition(aNuc->GetPosition());
605             G4double localDensity = projectile    593             G4double localDensity = projectile3dNucleus->GetNuclearDensity()->GetDensity(aPosition);
606             G4double localPfermi = theFermi.Ge    594             G4double localPfermi = theFermi.GetFermiMomentum(localDensity);
607             G4double nucMass = aNuc->GetDefini    595             G4double nucMass = aNuc->GetDefinition()->GetPDGMass();
608             G4double localFermiEnergy = std::s    596             G4double localFermiEnergy = std::sqrt(nucMass*nucMass + localPfermi*localPfermi) - nucMass;
609             G4double deltaE = localFermiEnergy    597             G4double deltaE = localFermiEnergy - (aNuc->GetMomentum().t()-aNuc->GetMomentum().mag());
610             theStatisticalExEnergy += deltaE;     598             theStatisticalExEnergy += deltaE;
611          }                                        599          }
612       }                                           600       }
613       return theStatisticalExEnergy;              601       return theStatisticalExEnergy;
614 }                                                 602 }
615                                                   603 
616 G4LorentzVector G4BinaryLightIonReaction::Sort    604 G4LorentzVector G4BinaryLightIonReaction::SortResult(G4ReactionProductVector * result, G4ReactionProductVector * spectators,G4ReactionProductVector * cascaders)
617 {                                                 605 {
618    unsigned int i(0);                             606    unsigned int i(0);
619    spectatorA=spectatorZ=0;                       607    spectatorA=spectatorZ=0;
620    G4LorentzVector pspectators(0,0,0,0);          608    G4LorentzVector pspectators(0,0,0,0);
621    pFinalState=G4LorentzVector(0,0,0,0);          609    pFinalState=G4LorentzVector(0,0,0,0);
622    for(i=0; i<result->size(); i++)                610    for(i=0; i<result->size(); i++)
623    {                                              611    {
624       if( (*result)[i]->GetNewlyAdded() )         612       if( (*result)[i]->GetNewlyAdded() )
625       {                                           613       {
626          pFinalState += G4LorentzVector( (*res    614          pFinalState += G4LorentzVector( (*result)[i]->GetMomentum(), (*result)[i]->GetTotalEnergy() );
627          cascaders->push_back((*result)[i]);      615          cascaders->push_back((*result)[i]);
628       }                                           616       }
629       else {                                      617       else {
630          //          G4cout <<" spectator ...     618          //          G4cout <<" spectator ... ";
631          pspectators += G4LorentzVector( (*res    619          pspectators += G4LorentzVector( (*result)[i]->GetMomentum(), (*result)[i]->GetTotalEnergy() );
632          spectators->push_back((*result)[i]);     620          spectators->push_back((*result)[i]);
633          spectatorA++;                            621          spectatorA++;
634          spectatorZ+= G4lrint((*result)[i]->Ge    622          spectatorZ+= G4lrint((*result)[i]->GetDefinition()->GetPDGCharge()/eplus);
635       }                                           623       }
636                                                   624 
637       //       G4cout << (*result)[i]<< " "       625       //       G4cout << (*result)[i]<< " "
638       //        << (*result)[i]->GetDefinition    626       //        << (*result)[i]->GetDefinition()->GetParticleName() << " "
639       //        << (*result)[i]->GetMomentum()    627       //        << (*result)[i]->GetMomentum()<< " "
640       //        << (*result)[i]->GetTotalEnerg    628       //        << (*result)[i]->GetTotalEnergy() << G4endl;
641    }                                              629    }
642       //G4cout << "pFinalState / pspectators,     630       //G4cout << "pFinalState / pspectators, (A,Z), p " << pFinalState << " / " << spectators->size()
643       //    << " (" << spectatorA << ", "<< sp    631       //    << " (" << spectatorA << ", "<< spectatorZ  << "), 4-mom: " << pspectators << G4endl;
644                                                   632 
645    return pspectators;                            633    return pspectators;
646 }                                                 634 }
647                                                   635 
648 void G4BinaryLightIonReaction::DeExciteSpectat    636 void G4BinaryLightIonReaction::DeExciteSpectatorNucleus(G4ReactionProductVector * spectators, G4ReactionProductVector * cascaders,
649                                                   637                                                  G4double theStatisticalExEnergy, G4LorentzVector & pSpectators)
650 {                                                 638 {
651    // call precompound model                      639    // call precompound model
652    G4ReactionProductVector * proFrag = 0;         640    G4ReactionProductVector * proFrag = 0;
653    G4LorentzVector pFragment(0.,0.,0.,0.);        641    G4LorentzVector pFragment(0.,0.,0.,0.);
654    //      G4cout << " == pre boost 1 "<< mome    642    //      G4cout << " == pre boost 1 "<< momentum.e()<< " "<< momentum.mag()<<G4endl;
655    G4LorentzRotation boost_fragments;             643    G4LorentzRotation boost_fragments;
656    //      G4cout << " == post boost 1 "<< mom    644    //      G4cout << " == post boost 1 "<< momentum.e()<< " "<< momentum.mag()<<G4endl;
657    //    G4LorentzRotation boost_spectator_mom    645    //    G4LorentzRotation boost_spectator_mom(-momentum.boostVector());
658    //     G4cout << "- momentum " << boost_spe    646    //     G4cout << "- momentum " << boost_spectator_mom * momentum << G4endl;
659    G4LorentzVector pFragments(0,0,0,0);           647    G4LorentzVector pFragments(0,0,0,0);
660                                                   648 
661    if(spectatorZ>0 && spectatorA>1)               649    if(spectatorZ>0 && spectatorA>1)
662    {                                              650    {
663       //  Make the fragment                       651       //  Make the fragment
664       G4Fragment aProRes;                         652       G4Fragment aProRes;
665       aProRes.SetZandA_asInt(spectatorZ, spect    653       aProRes.SetZandA_asInt(spectatorZ, spectatorA);
666       aProRes.SetNumberOfParticles(0);            654       aProRes.SetNumberOfParticles(0);
667       aProRes.SetNumberOfCharged(0);              655       aProRes.SetNumberOfCharged(0);
668       aProRes.SetNumberOfHoles(pA-spectatorA);    656       aProRes.SetNumberOfHoles(pA-spectatorA);
669       G4double mFragment=G4ParticleTable::GetP    657       G4double mFragment=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(spectatorZ,spectatorA);
670       pFragment=G4LorentzVector(0,0,0,mFragmen    658       pFragment=G4LorentzVector(0,0,0,mFragment+std::max(0.,theStatisticalExEnergy) );
671       aProRes.SetMomentum(pFragment);             659       aProRes.SetMomentum(pFragment);
672                                                   660 
673       proFrag = theHandler->BreakItUp(aProRes)    661       proFrag = theHandler->BreakItUp(aProRes);
674                                                   662 
675       boost_fragments = G4LorentzRotation(pSpe    663       boost_fragments = G4LorentzRotation(pSpectators.boostVector());
676                                                   664 
677       //     G4cout << " Fragment a,z, Mass Fr    665       //     G4cout << " Fragment a,z, Mass Fragment, mass spect-mom, exitationE "
678       //       << spectatorA <<" "<< spectator    666       //       << spectatorA <<" "<< spectatorZ <<" "<< mFragment <<" "
679       //       << momentum.mag() <<" "<< momen    667       //       << momentum.mag() <<" "<< momentum.mag() - mFragment
680       //       << " "<<theStatisticalExEnergy     668       //       << " "<<theStatisticalExEnergy
681       //       << " "<< boost_fragments*pFragm    669       //       << " "<< boost_fragments*pFragment<< G4endl;
682       G4ReactionProductVector::iterator ispect    670       G4ReactionProductVector::iterator ispectator;
683       for (ispectator=spectators->begin();ispe    671       for (ispectator=spectators->begin();ispectator!=spectators->end();ispectator++)
684       {                                           672       {
685          delete *ispectator;                      673          delete *ispectator;
686       }                                           674       }
687    }                                              675    }
688    else if(spectatorA!=0)                         676    else if(spectatorA!=0)
689    {                                              677    {
690      G4ReactionProductVector::iterator ispecta    678      G4ReactionProductVector::iterator ispectator;
691      for (ispectator=spectators->begin();ispec    679      for (ispectator=spectators->begin();ispectator!=spectators->end();ispectator++)
692       {                                           680       {
693          (*ispectator)->SetNewlyAdded(true);      681          (*ispectator)->SetNewlyAdded(true);
694          cascaders->push_back(*ispectator);       682          cascaders->push_back(*ispectator);
695          pFinalState+=G4LorentzVector((*ispect    683          pFinalState+=G4LorentzVector((*ispectator)->GetMomentum(),(*ispectator)->GetTotalEnergy());
696                   //G4cout << "BLIC: spectator    684                   //G4cout << "BLIC: spectatorA>0, Z=0 from spectator "
697                   // << (*ispectator)->GetDefi    685                   // << (*ispectator)->GetDefinition()->GetParticleName() << " "
698                   // << (*ispectator)->GetMome    686                   // << (*ispectator)->GetMomentum()<< " "
699                   // << (*ispectator)->GetTota    687                   // << (*ispectator)->GetTotalEnergy() << G4endl;
700       }                                           688       }
701                                                   689 
702    }                                              690    }
703    // / if (spectators)                           691    // / if (spectators)
704    delete spectators;                             692    delete spectators;
705                                                   693 
706    // collect the evaporation part and boost t    694    // collect the evaporation part and boost to spectator frame
707    G4ReactionProductVector::iterator ii;          695    G4ReactionProductVector::iterator ii;
708    if(proFrag)                                    696    if(proFrag)
709    {                                              697    {
710       for(ii=proFrag->begin(); ii!=proFrag->en    698       for(ii=proFrag->begin(); ii!=proFrag->end(); ii++)
711       {                                           699       {
712          (*ii)->SetNewlyAdded(true);              700          (*ii)->SetNewlyAdded(true);
713          G4LorentzVector tmp((*ii)->GetMomentu    701          G4LorentzVector tmp((*ii)->GetMomentum(),(*ii)->GetTotalEnergy());
714          tmp *= boost_fragments;                  702          tmp *= boost_fragments;
715          (*ii)->SetMomentum(tmp.vect());          703          (*ii)->SetMomentum(tmp.vect());
716          (*ii)->SetTotalEnergy(tmp.e());          704          (*ii)->SetTotalEnergy(tmp.e());
717          //      result->push_back(*ii);          705          //      result->push_back(*ii);
718          pFragments += tmp;                       706          pFragments += tmp;
719       }                                           707       }
720    }                                              708    }
721                                                   709 
722    //    G4cout << "Fragmented p, momentum, de    710    //    G4cout << "Fragmented p, momentum, delta " << pFragments <<" "<<momentum
723    //            <<" "<< pFragments-momentum <    711    //            <<" "<< pFragments-momentum << G4endl;
724                                                   712 
725    //  correct p/E of Cascade secondaries         713    //  correct p/E of Cascade secondaries
726    G4LorentzVector pCas=pInitialState - pFragm    714    G4LorentzVector pCas=pInitialState - pFragments;
727                                                   715 
728        //G4cout <<"BLIC: Going to correct from    716        //G4cout <<"BLIC: Going to correct from " << pFinalState << " to " << pCas << G4endl;
729    //  the creation of excited fragment did vi    717    //  the creation of excited fragment did violate E/p, so correct cascaders to get overall conservation.
730    G4bool EnergyIsCorrect=EnergyAndMomentumCor    718    G4bool EnergyIsCorrect=EnergyAndMomentumCorrector(cascaders, pCas);
731    if ( ! EnergyIsCorrect && debug_G4BinaryLig    719    if ( ! EnergyIsCorrect && debug_G4BinaryLightIonReactionResults)
732    {                                              720    {
733       G4cout << "G4BinaryLightIonReaction E/P     721       G4cout << "G4BinaryLightIonReaction E/P correction for nucleus failed, will try to correct overall" << G4endl;
734    }                                              722    }
735                                                   723 
736    //  Add deexcitation secondaries               724    //  Add deexcitation secondaries
737    if(proFrag)                                    725    if(proFrag)
738    {                                              726    {
739       for(ii=proFrag->begin(); ii!=proFrag->en    727       for(ii=proFrag->begin(); ii!=proFrag->end(); ii++)
740       {                                           728       {
741          cascaders->push_back(*ii);               729          cascaders->push_back(*ii);
742       }                                           730       }
743       delete proFrag;                             731       delete proFrag;
744    }                                              732    }
745       //G4cout << "EnergyIsCorrect? " << Energ    733       //G4cout << "EnergyIsCorrect? " << EnergyIsCorrect << G4endl;
746    if ( ! EnergyIsCorrect )                       734    if ( ! EnergyIsCorrect )
747    {                                              735    {
748          // G4cout <<" ! EnergyIsCorrect " <<     736          // G4cout <<" ! EnergyIsCorrect " << pFinalState << " to " << pInitialState << G4endl;
749       if (! EnergyAndMomentumCorrector(cascade    737       if (! EnergyAndMomentumCorrector(cascaders,pInitialState))
750       {                                           738       {
751          if(debug_G4BinaryLightIonReactionResu    739          if(debug_G4BinaryLightIonReactionResults)
752             G4cout << "G4BinaryLightIonReactio    740             G4cout << "G4BinaryLightIonReaction E/P corrections failed" << G4endl;
753       }                                           741       }
754    }                                              742    }
755                                                   743 
756 }                                                 744 }
757                                                   745 
758                                                   746