G4BinaryCascade.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 
00027 #include "globals.hh"
00028 #include "G4PhysicalConstants.hh"
00029 #include "G4SystemOfUnits.hh"
00030 #include "G4Proton.hh"
00031 #include "G4Neutron.hh"
00032 #include "G4LorentzRotation.hh"
00033 #include "G4BinaryCascade.hh"
00034 #include "G4KineticTrackVector.hh"
00035 #include "G4DecayKineticTracks.hh"
00036 #include "G4ReactionProductVector.hh"
00037 #include "G4Track.hh"
00038 #include "G4V3DNucleus.hh"
00039 #include "G4Fancy3DNucleus.hh"
00040 #include "G4Scatterer.hh"
00041 #include "G4MesonAbsorption.hh"
00042 #include "G4ping.hh"
00043 #include "G4Delete.hh"
00044 
00045 #include "G4CollisionManager.hh"
00046 #include "G4Absorber.hh"
00047 
00048 #include "G4CollisionInitialState.hh"
00049 #include "G4ListOfCollisions.hh"
00050 #include "G4Fragment.hh"
00051 #include "G4RKPropagation.hh"
00052 
00053 #include "G4NuclearShellModelDensity.hh"
00054 #include "G4NuclearFermiDensity.hh"
00055 #include "G4FermiMomentum.hh"
00056 
00057 #include "G4PreCompoundModel.hh"
00058 #include "G4ExcitationHandler.hh"
00059 #include "G4HadronicInteractionRegistry.hh"
00060 
00061 #include "G4FermiPhaseSpaceDecay.hh"
00062 
00063 #include "G4PreCompoundModel.hh"
00064 
00065 #include <algorithm>
00066 #include "G4ShortLivedConstructor.hh"
00067 #include <typeinfo>
00068 
00069 //   turn on general debugging info, and consistency checks
00070 //#define debug_G4BinaryCascade 1
00071 
00072 //  more detailed debugging -- deprecated
00073 //#define debug_H1_BinaryCascade 1
00074 
00075 //  specific debugging info per method or functionality
00076 //#define debug_BIC_ApplyCollision 1
00077 //#define debug_BIC_CheckPauli 1
00078 //#define debug_BIC_CorrectFinalPandE 1
00079 //#define debug_BIC_Propagate 1
00080 //#define debug_BIC_Propagate_Excitation 1
00081 //#define debug_BIC_Propagate_Collisions 1
00082 //#define debug_BIC_Propagate_finals 1
00083 //#define debug_BIC_DoTimeStep 1
00084 //#define debug_BIC_CorrectBarionsOnBoundary 1
00085 //#define debug_BIC_GetExcitationEnergy 1
00086 //#define debug_BIC_FinalNucleusMomentum 1
00087 //#define debug_BIC_Final4Momentum 1
00088 //#define debug_BIC_FindFragments 1
00089 //#define debug_BIC_BuildTargetList 1
00090 //#define debug_BIC_FindCollision 1
00091 //#define debug_BIC_return 1
00092 
00093 #if defined(debug_G4BinaryCascade)
00094   #define _CheckChargeAndBaryonNumber_(val) CheckChargeAndBaryonNumber(val)
00095 #else
00096   #define _CheckChargeAndBaryonNumber_(val)
00097 #endif
00098 //
00099 //  C O N S T R U C T O R S   A N D   D E S T R U C T O R S
00100 //
00101 
00102 G4BinaryCascade::G4BinaryCascade(G4VPreCompoundModel* ptr) :
00103 G4VIntraNuclearTransportModel("Binary Cascade", ptr)
00104 {
00105     // initialise the resonance sector
00106     G4ShortLivedConstructor ShortLived;
00107     ShortLived.ConstructParticle();
00108 
00109     theCollisionMgr = new G4CollisionManager;
00110     theDecay=new G4BCDecay;
00111     theImR.push_back(theDecay);
00112     theLateParticle= new G4BCLateParticle;
00113     G4MesonAbsorption * aAb=new G4MesonAbsorption;
00114     theImR.push_back(aAb);
00115     G4Scatterer * aSc=new G4Scatterer;
00116     theH1Scatterer = new G4Scatterer;
00117     theImR.push_back(aSc);
00118 
00119     thePropagator = new G4RKPropagation;
00120     theCurrentTime = 0.;
00121     theBCminP = 45*MeV;
00122     theCutOnP = 90*MeV;
00123     theCutOnPAbsorb= 0*MeV;   // No Absorption of slow Mesons, other than above G4MesonAbsorption
00124 
00125     // reuse existing pre-compound model
00126     if(!ptr) {
00127       G4HadronicInteraction* p =
00128         G4HadronicInteractionRegistry::Instance()->FindModel("PRECO");
00129       G4VPreCompoundModel* pre = static_cast<G4VPreCompoundModel*>(p);
00130       if(!pre) { pre = new G4PreCompoundModel(); }
00131       SetDeExcitation(pre);
00132     }
00133     theExcitationHandler = GetDeExcitation()->GetExcitationHandler();
00134     SetMinEnergy(0.0*GeV);
00135     SetMaxEnergy(10.1*GeV);
00136     //PrintWelcomeMessage();
00137     thePrimaryEscape = true;
00138     thePrimaryType = 0;
00139 
00140     SetEnergyMomentumCheckLevels(1*perCent, 1*MeV);
00141 
00142     // init data members
00143     currentA=currentZ=0;
00144     lateA=lateZ=0;
00145     initialA=initialZ=0;
00146     projectileA=projectileZ=0;
00147     currentInitialEnergy=initial_nuclear_mass=0.;
00148     massInNucleus=0.;
00149     theOuterRadius=0.;
00150 }
00151 
00152 /*
00153 G4BinaryCascade::G4BinaryCascade(const G4BinaryCascade& )
00154 : G4VIntraNuclearTransportModel("Binary Cascade")
00155 {
00156 }
00157  */
00158 
00159 G4BinaryCascade::~G4BinaryCascade()
00160 {
00161     ClearAndDestroy(&theTargetList);
00162     ClearAndDestroy(&theSecondaryList);
00163     ClearAndDestroy(&theCapturedList);
00164     delete thePropagator;
00165     delete theCollisionMgr;
00166     std::for_each(theImR.begin(), theImR.end(), Delete<G4BCAction>());
00167     delete theLateParticle;
00168     //delete theExcitationHandler;
00169     delete theH1Scatterer;
00170 }
00171 
00172 void G4BinaryCascade::ModelDescription(std::ostream& outFile) const
00173 {
00174     outFile << "G4BinaryCascade is an intra-nuclear cascade model in which\n"
00175             << "an incident hadron collides with a nucleon, forming two\n"
00176             << "final-state particles, one or both of which may be resonances.\n"
00177             << "The resonances then decay hadronically and the decay products\n"
00178             << "are then propagated through the nuclear potential along curved\n"
00179             << "trajectories until they re-interact or leave the nucleus.\n"
00180             << "This model is valid for incident pions up to 1.5 GeV and\n"
00181             << "nucleons up to 10 GeV.\n";
00182 }
00183 void G4BinaryCascade::PropagateModelDescription(std::ostream& outFile) const
00184 {
00185     outFile << "G4BinaryCascade propagtes secondaries produced by a high\n"
00186             << "energy model through the wounded nucleus.\n"
00187             << "Secondaries are followed after the formation time and if\n"
00188             << "within the nucleus are propagated through the nuclear\n"
00189             << "potential along curved trajectories until they interact\n"
00190             << "with a nucleon, decay, or leave the nucleus.\n"
00191             << "An interaction of a secondary with a nucleon produces two\n"
00192             << "final-state particles, one or both of which may be resonances.\n"
00193             << "Resonances decay hadronically and the decay products\n"
00194             << "are in turn propagated through the nuclear potential along curved\n"
00195             << "trajectories until they re-interact or leave the nucleus.\n"
00196             << "This model is valid for pions up to 1.5 GeV and\n"
00197             << "nucleons up to about 3.5 GeV.\n";
00198 }
00199 
00200 //----------------------------------------------------------------------------
00201 
00202 //
00203 //      I M P L E M E N T A T I O N
00204 //
00205 
00206 
00207 //----------------------------------------------------------------------------
00208 G4HadFinalState * G4BinaryCascade::ApplyYourself(const G4HadProjectile & aTrack,
00209         G4Nucleus & aNucleus)
00210 //----------------------------------------------------------------------------
00211 {
00212     static G4int eventcounter=0;
00213 
00214     //   if ( eventcounter == 0 ) {
00215     //      SetEpReportLevel(3);   // report non conservation with model etc.
00216     //      G4double relativeLevel = 1*perCent;
00217     //      G4double absoluteLevel = 2*MeV;
00218     //      SetEnergyMomentumCheckLevels(relativeLevel,absoluteLevel);
00219     //   }
00220 
00221     //if(eventcounter == 100*(eventcounter/100) )
00222     eventcounter++;
00223     if(getenv("BCDEBUG") ) G4cerr << " ######### Binary Cascade Reaction number starts ######### "<<eventcounter<<G4endl;
00224 
00225     G4LorentzVector initial4Momentum = aTrack.Get4Momentum();
00226     G4ParticleDefinition * definition = const_cast<G4ParticleDefinition *>(aTrack.GetDefinition());
00227 
00228     if(initial4Momentum.e()-initial4Momentum.m()<theBCminP &&
00229             ( definition==G4Neutron::NeutronDefinition() || definition==G4Proton::ProtonDefinition() ) )
00230     {
00231         return theDeExcitation->ApplyYourself(aTrack, aNucleus);
00232     }
00233 
00234     theParticleChange.Clear();
00235     // initialize the G4V3DNucleus from G4Nucleus
00236     the3DNucleus = new G4Fancy3DNucleus;
00237 
00238     // Build a KineticTrackVector with the G4Track
00239     G4KineticTrackVector * secondaries;// = new G4KineticTrackVector;
00240     G4ThreeVector initialPosition(0., 0., 0.); // will be set later
00241 
00242     if(!getenv("I_Am_G4BinaryCascade_Developer") )
00243     {
00244         if(definition!=G4Neutron::NeutronDefinition() &&
00245                 definition!=G4Proton::ProtonDefinition() &&
00246                 definition!=G4PionPlus::PionPlusDefinition() &&
00247                 definition!=G4PionMinus::PionMinusDefinition() )
00248         {
00249             G4cerr << "You are using G4BinaryCascade for projectiles other than nucleons or pions."<<G4endl;
00250             G4cerr << "If you want to continue, please switch on the developer environment: "<<G4endl;
00251             G4cerr << "setenv I_Am_G4BinaryCascade_Developer 1 "<<G4endl<<G4endl;
00252             throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade - used for unvalid particle type - Fatal");
00253         }
00254     }
00255 
00256     // keep primary
00257     thePrimaryType = definition;
00258     thePrimaryEscape = false;
00259 
00260     // try until an interaction will happen
00261     G4ReactionProductVector * products = 0;
00262     G4int interactionCounter = 0;
00263     do
00264     {
00265         // reset status that could be changed in previous loop event
00266         theCollisionMgr->ClearAndDestroy();
00267 
00268         if(products != 0)
00269         {  // free memory from previous loop event
00270             ClearAndDestroy(products);
00271             delete products;
00272             products=0;
00273         }
00274 
00275         G4int massNumber=aNucleus.GetA_asInt();
00276         the3DNucleus->Init(massNumber, aNucleus.GetZ_asInt());
00277         thePropagator->Init(the3DNucleus);
00278         //      GF Leak on kt??? but where to delete?
00279         G4KineticTrack * kt;// = new G4KineticTrack(definition, 0., initialPosition, initial4Momentum);
00280         do                  // sample impact parameter until collisions are found
00281         {
00282             theCurrentTime=0;
00283             G4double radius = the3DNucleus->GetOuterRadius()+3*fermi;
00284             initialPosition=GetSpherePoint(1.1*radius, initial4Momentum);  // get random position
00285             kt = new G4KineticTrack(definition, 0., initialPosition, initial4Momentum);
00286             kt->SetState(G4KineticTrack::outside);
00287             // secondaries has been cleared by Propagate() in the previous loop event
00288             secondaries= new G4KineticTrackVector;
00289             secondaries->push_back(kt);
00290             if(massNumber > 1) // 1H1 is special case
00291             {
00292                 products = Propagate(secondaries, the3DNucleus);
00293             } else {
00294                 products = Propagate1H1(secondaries,the3DNucleus);
00295             }
00296         } while(! products );  // until we FIND a collision...
00297 
00298         if(++interactionCounter>99) break;
00299     } while(products->size() == 0);  // ...until we find an ALLOWED collision
00300 
00301     if(products->size()>0)
00302     {
00303         //  G4cout << "BIC Applyyourself: number of products " << products->size() << G4endl;
00304 
00305         // Fill the G4ParticleChange * with products
00306         theParticleChange.SetStatusChange(stopAndKill);
00307         G4ReactionProductVector::iterator iter;
00308 
00309         for(iter = products->begin(); iter != products->end(); ++iter)
00310         {
00311             G4DynamicParticle * aNew =
00312                     new G4DynamicParticle((*iter)->GetDefinition(),
00313                             (*iter)->GetTotalEnergy(),
00314                             (*iter)->GetMomentum());
00315             theParticleChange.AddSecondary(aNew);
00316         }
00317 
00318         // DebugEpConservation(track, products);
00319 
00320 
00321     } else {  // no interaction, return primary
00322         if(getenv("BCDEBUG") ) G4cerr << " ######### Binary Cascade Reaction number void ######### "<<eventcounter<<G4endl;
00323         theParticleChange.SetStatusChange(isAlive);
00324         theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
00325         theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
00326     }
00327 
00328     ClearAndDestroy(products);
00329     delete products;
00330 
00331     delete the3DNucleus;
00332     the3DNucleus = NULL;
00333 
00334     if(getenv("BCDEBUG") ) G4cerr << " ######### Binary Cascade Reaction number ends ######### "<<eventcounter<<G4endl;
00335 
00336     return &theParticleChange;
00337 }
00338 //----------------------------------------------------------------------------
00339 G4ReactionProductVector * G4BinaryCascade::Propagate(
00340         G4KineticTrackVector * secondaries, G4V3DNucleus * aNucleus)
00341 //----------------------------------------------------------------------------
00342 {
00343     G4ping debug("debug_G4BinaryCascade");
00344 #ifdef debug_BIC_Propagate
00345     G4cout << "G4BinaryCascade Propagate starting -------------------------------------------------------" <<G4endl;
00346 #endif
00347 
00348     the3DNucleus=aNucleus;
00349     G4ReactionProductVector * products = new G4ReactionProductVector;
00350     theOuterRadius = the3DNucleus->GetOuterRadius();
00351     theCurrentTime=0;
00352     theProjectile4Momentum=G4LorentzVector(0,0,0,0);
00353     // build theSecondaryList, theProjectileList and theCapturedList
00354     ClearAndDestroy(&theCapturedList);
00355     ClearAndDestroy(&theSecondaryList);
00356     theSecondaryList.clear();
00357     ClearAndDestroy(&theFinalState);
00358     std::vector<G4KineticTrack *>::iterator iter;
00359     theCollisionMgr->ClearAndDestroy();
00360 
00361     theCutOnP=90*MeV;
00362     if(the3DNucleus->GetMass()>30) theCutOnP = 70*MeV;
00363     if(the3DNucleus->GetMass()>60) theCutOnP = 50*MeV;
00364     if(the3DNucleus->GetMass()>120) theCutOnP = 45*MeV;
00365 
00366 
00367     BuildTargetList();
00368 
00369 #ifdef debug_BIC_GetExcitationEnergy
00370     G4cout << "ExcitationEnergy0 " << GetExcitationEnergy() << G4endl;
00371 #endif
00372 
00373     thePropagator->Init(the3DNucleus);
00374 
00375     G4bool success = BuildLateParticleCollisions(secondaries);
00376     if (! success )   // fails if no excitation energy left....
00377     {
00378        products=HighEnergyModelFSProducts(products, secondaries);
00379        ClearAndDestroy(secondaries);
00380        delete secondaries;
00381 
00382 #ifdef debug_G4BinaryCascade
00383        G4cout << "G4BinaryCascade::Propagate: warning - high energy model failed energy conservation, returning unchanged high energy final state" << G4endl;
00384 #endif
00385 
00386        return products;
00387     }
00388     // check baryon and charge ...
00389 
00390     _CheckChargeAndBaryonNumber_("lateparticles");
00391 
00392     // if called stand alone find first collisions
00393     FindCollisions(&theSecondaryList);
00394 
00395 
00396     if(theCollisionMgr->Entries() == 0 )      //late particles ALWAYS create Entries
00397     {
00398         //G4cout << " no collsions -> return 0" << G4endl;
00399         delete products;
00400 #ifdef debug_BIC_return
00401         G4cout << "return @ begin2,  no collisions "<< G4endl;
00402 #endif
00403         return 0;
00404     }
00405 
00406     // end of initialization: do the job now
00407     // loop until there are no more collisions
00408 
00409 
00410     G4bool haveProducts = false;
00411     G4int collisionCount=0;
00412     theMomentumTransfer=G4ThreeVector(0,0,0);
00413     while(theCollisionMgr->Entries() > 0 && currentZ)
00414     {
00415         if(Absorb()) {  // absorb secondaries, pions only
00416             haveProducts = true;
00417         }
00418         _CheckChargeAndBaryonNumber_("absorbed");
00419         if(Capture()) { // capture secondaries, nucleons only
00420             haveProducts = true;
00421         }
00422         _CheckChargeAndBaryonNumber_("captured");
00423 
00424         // propagate to the next collision if any (collisions could have been deleted
00425         // by previous absorption or capture)
00426         if(theCollisionMgr->Entries() > 0)
00427         {
00428             G4CollisionInitialState *
00429             nextCollision = theCollisionMgr->GetNextCollision();
00430 #ifdef debug_BIC_Propagate_Collisions
00431             G4cout << " NextCollision  * , Time, curtime = " << nextCollision << " "
00432                     <<nextCollision->GetCollisionTime()<< " " <<
00433                     theCurrentTime<< G4endl;
00434 #endif
00435             if (!DoTimeStep(nextCollision->GetCollisionTime()-theCurrentTime) )
00436             {
00437                 // Check if nextCollision is still valid, ie. particle did not leave nucleus
00438                 if (theCollisionMgr->GetNextCollision() != nextCollision )
00439                 {
00440                     nextCollision = 0;
00441                 }
00442             }
00443 
00444             if( nextCollision )
00445             {
00446                 if (ApplyCollision(nextCollision))
00447                 {
00448                     //G4cerr << "ApplyCollision success " << G4endl;
00449                     haveProducts = true;
00450                     collisionCount++;
00451                     _CheckChargeAndBaryonNumber_("ApplyCollision");
00452 
00453                 } else {
00454                     //G4cerr << "ApplyCollision failure " << G4endl;
00455                     theCollisionMgr->RemoveCollision(nextCollision);
00456                 }
00457             }
00458         }
00459     }
00460 
00461     //--------- end of while on Collisions
00462     //G4cout << "currentZ @ end loop " << currentZ << G4endl;
00463     if ( ! currentZ ){
00464         // nucleus completely destroyed, fill in ReactionProductVector
00465         products = FillVoidNucleusProducts(products);
00466 #ifdef debug_BIC_return
00467         G4cout << "return @ Z=0 after collision loop "<< G4endl;
00468         PrintKTVector(&theSecondaryList,std::string(" theSecondaryList"));
00469         G4cout << "theTargetList size: " << theTargetList.size() << G4endl;
00470         PrintKTVector(&theTargetList,std::string(" theTargetList"));
00471         PrintKTVector(&theCapturedList,std::string(" theCapturedList"));
00472 
00473         G4cout << " ExcitE be4 Correct : " <<GetExcitationEnergy() << G4endl;
00474         G4cout << " Mom Transfered to nucleus : " << theMomentumTransfer << " " << theMomentumTransfer.mag() << G4endl;
00475         PrintKTVector(&theFinalState,std::string(" FinalState uncorrected"));
00476         G4cout << "returned products: " << products->size() << G4endl;
00477 #endif
00478         // @@GF FixMe cannot return here.
00479         return products;
00480     }
00481 
00482     // No more collisions: absorb, capture and propagate the secondaries out of the nucleus
00483     if(Absorb()) {
00484         haveProducts = true;
00485         // G4cout << "Absorb sucess " << G4endl;
00486     }
00487 
00488     if(Capture()) {
00489         haveProducts = true;
00490         // G4cout << "Capture sucess " << G4endl;
00491     }
00492 
00493     if(!haveProducts)  // no collisions happened. Return an empty vector.
00494     {
00495 #ifdef debug_BIC_return
00496         G4cout << "return 3, no products "<< G4endl;
00497 #endif
00498         return products;
00499     }
00500 
00501 
00502 #ifdef debug_BIC_Propagate
00503     G4cout << " Momentum transfer to Nucleus " << theMomentumTransfer << " " << theMomentumTransfer.mag() << G4endl;
00504     G4cout << "  Stepping particles out...... " << G4endl;
00505 #endif
00506 
00507     StepParticlesOut();
00508 
00509 
00510     if ( theSecondaryList.size() > 0 )
00511     {
00512 #ifdef debug_G4BinaryCascade
00513         G4cerr << "G4BinaryCascade: Warning, have active particles at end" << G4endl;
00514         PrintKTVector(&theSecondaryList, "active particles @ end  added to theFinalState");
00515 #endif
00516         //  add left secondaries to FinalSate
00517         for ( iter =theSecondaryList.begin(); iter != theSecondaryList.end(); ++iter)
00518         {
00519             theFinalState.push_back(*iter);
00520         }
00521         theSecondaryList.clear();
00522 
00523     }
00524     while ( theCollisionMgr->Entries() > 0 )
00525     {
00526 #ifdef debug_G4BinaryCascade
00527         G4cerr << " Warning: remove left over collision(s) " << G4endl;
00528 #endif
00529         theCollisionMgr->RemoveCollision(theCollisionMgr->GetNextCollision());
00530     }
00531 
00532 #ifdef debug_BIC_Propagate_Excitation
00533 
00534     PrintKTVector(&theSecondaryList,std::string(" theSecondaryList"));
00535     G4cout << "theTargetList size: " << theTargetList.size() << G4endl;
00536     //  PrintKTVector(&theTargetList,std::string(" theTargetList"));
00537     PrintKTVector(&theCapturedList,std::string(" theCapturedList"));
00538 
00539     G4cout << " ExcitE be4 Correct : " <<GetExcitationEnergy() << G4endl;
00540     G4cout << " Mom Transfered to nucleus : " << theMomentumTransfer << " " << theMomentumTransfer.mag() << G4endl;
00541     PrintKTVector(&theFinalState,std::string(" FinalState uncorrected"));
00542 #endif
00543 
00544     //
00545 
00546 
00547     G4double ExcitationEnergy=GetExcitationEnergy();
00548 
00549 #ifdef debug_BIC_Propagate_finals
00550     PrintKTVector(&theFinalState,std::string(" FinalState be4 corr"));
00551     G4cout << " Excitation Energy prefinal,  #collisions:, out, captured  "
00552             << ExcitationEnergy << " "
00553             << collisionCount << " "
00554             << theFinalState.size() << " "
00555             << theCapturedList.size()<<G4endl;
00556 #endif
00557 
00558     if (ExcitationEnergy < 0 )
00559     {
00560         G4int maxtry=5, ntry=0;
00561         do {
00562             CorrectFinalPandE();
00563             ExcitationEnergy=GetExcitationEnergy();
00564         } while ( ++ntry < maxtry && ExcitationEnergy < 0 );
00565     }
00566 
00567 #ifdef debug_BIC_Propagate_finals
00568     PrintKTVector(&theFinalState,std::string(" FinalState corrected"));
00569     G4cout << " Excitation Energy final,  #collisions:, out, captured  "
00570             << ExcitationEnergy << " "
00571             << collisionCount << " "
00572             << theFinalState.size() << " "
00573             << theCapturedList.size()<<G4endl;
00574 #endif
00575 
00576 
00577     if ( ExcitationEnergy < 0. )
00578     {
00579         //      if ( ExcitationEnergy < 0. )
00580         {
00581             //#ifdef debug_G4BinaryCascade
00582             //            G4cerr << "G4BinaryCascade-Warning: negative excitation energy ";
00583             //            G4cerr <<ExcitationEnergy<<G4endl;
00584             //     PrintKTVector(&theFinalState,std::string("FinalState"));
00585             //    PrintKTVector(&theCapturedList,std::string("captured"));
00586             //    G4cout << "negative ExE:Final 4Momentum .mag: " << GetFinal4Momentum()
00587             //            << " "<< GetFinal4Momentum().mag()<< G4endl
00588             //            << "negative ExE:FinalNucleusMom  .mag: " << GetFinalNucleusMomentum()
00589             //            << " "<< GetFinalNucleusMomentum().mag()<< G4endl;
00590             //#endif
00591         }
00592         ClearAndDestroy(products);
00593         //G4cout << "  negative Excitation E return empty products " << products << G4endl;
00594 #ifdef debug_BIC_return
00595         G4cout << "return 4, excit < 0 "<< G4endl;
00596 #endif
00597         return products;   // return empty products- FIXME
00598     }
00599 
00600     G4ReactionProductVector * precompoundProducts=DeExcite();
00601 
00602 
00603     G4DecayKineticTracks decay(&theFinalState);
00604 
00605     products= ProductsAddFinalState(products, theFinalState);
00606 
00607     products= ProductsAddPrecompound(products, precompoundProducts);
00608 
00609 //    products=ProductsAddFakeGamma(products);
00610 
00611 
00612     thePrimaryEscape = true;
00613 
00614     #ifdef debug_BIC_return
00615     G4cout << "return @end, all ok "<< G4endl;
00616     //G4cout << "  return products " << products << G4endl;
00617     #endif
00618 
00619     return products;
00620 }
00621 
00622 //----------------------------------------------------------------------------
00623 G4double G4BinaryCascade::GetExcitationEnergy()
00624 //----------------------------------------------------------------------------
00625 {
00626 
00627     // get A and Z for the residual nucleus
00628 #if defined(debug_G4BinaryCascade) || defined(debug_BIC_GetExcitationEnergy)
00629     G4int finalA = theTargetList.size()+theCapturedList.size();
00630     G4int finalZ = GetTotalCharge(theTargetList)+GetTotalCharge(theCapturedList);
00631     if ( (currentA - finalA) != 0 || (currentZ - finalZ) != 0 )
00632     {
00633         G4cerr << "G4BIC:GetExcitationEnergy(): Nucleon counting error current/final{A,Z} "
00634                 << currentA << " " << finalA << " "<< currentZ << " " << finalZ << G4endl;
00635     }
00636 
00637 #endif
00638 
00639     G4double excitationE(0);
00640     G4double nucleusMass(0);
00641     if(currentZ>.5)
00642     {
00643         nucleusMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(currentZ,currentA);
00644     }
00645     else if (currentZ==0 )     // Uzhi && currentA==1 )                     // Uzhi
00646     {                                                                       // Uzhi
00647         if(currentA == 1) {nucleusMass = G4Neutron::Neutron()->GetPDGMass();}// Uzhi
00648         else              {nucleusMass = GetFinalNucleusMomentum().mag()     // Uzhi
00649             - 3.*MeV*currentA;}                 // Uzhi
00650     }                                                                       // Uzhi
00651     else
00652     {
00653 #ifdef debug_G4BinaryCascade
00654         G4cout << "G4BinaryCascade::GetExcitationEnergy(): Warning - invalid nucleus (A,Z)=("
00655                 << currentA << "," << currentZ << ")" << G4endl;
00656 #endif
00657         return 0;
00658     }
00659 
00660 #ifdef debug_BIC_GetExcitationEnergy
00661     G4ping debug("debug_ExcitationEnergy");
00662     debug.push_back("====> current A, Z");
00663     debug.push_back(currentZ);
00664     debug.push_back(currentA);
00665     debug.push_back("====> final A, Z");
00666     debug.push_back(finalZ);
00667     debug.push_back(finalA);
00668     debug.push_back(nucleusMass);
00669     debug.push_back(GetFinalNucleusMomentum().mag());
00670     debug.dump();
00671     //  PrintKTVector(&theTargetList, std::string(" current target list info"));
00672     //PrintKTVector(&theCapturedList, std::string(" current captured list info"));
00673 #endif
00674 
00675     excitationE = GetFinalNucleusMomentum().mag() - nucleusMass;
00676 
00677     //G4double exE2 = GetFinal4Momentum().mag() - nucleusMass;
00678 
00679     //G4cout << "old/new excitE " << excitationE << " / "<< exE2 << G4endl;
00680 
00681 #ifdef debug_BIC_GetExcitationEnergy
00682     // ------ debug
00683     if ( excitationE < 0 )
00684     {
00685         G4cout << "negative ExE final Ion mass " <<nucleusMass<< G4endl;
00686         G4LorentzVector Nucl_mom=GetFinalNucleusMomentum();
00687         if(finalZ>.5) G4cout << " Final nuclmom/mass " << Nucl_mom << " " << Nucl_mom.mag()
00688                                            << " (A,Z)=("<< finalA <<","<<finalZ <<")"
00689                                            << " mass " << nucleusMass << " "
00690                                            << " excitE " << excitationE << G4endl;
00691 
00692 
00693         G4int A = the3DNucleus->GetMassNumber();
00694         G4int Z = the3DNucleus->GetCharge();
00695         G4double initialExc(0);
00696         if(Z>.5)
00697         {
00698             initialExc = theInitial4Mom.mag()-
00699                     G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z, A);
00700             G4cout << "GetExcitationEnergy: Initial nucleus A Z " << A << " " << Z << " " << initialExc << G4endl;
00701         }
00702     }
00703 
00704 #endif
00705 
00706     return excitationE;
00707 }
00708 
00709 
00710 //----------------------------------------------------------------------------
00711 //
00712 //       P R I V A T E   M E M B E R   F U N C T I O N S
00713 //
00714 //----------------------------------------------------------------------------
00715 
00716 //----------------------------------------------------------------------------
00717 void G4BinaryCascade::BuildTargetList()
00718 //----------------------------------------------------------------------------
00719 {
00720 
00721     if(!the3DNucleus->StartLoop())
00722     {
00723         //    G4cerr << "G4BinaryCascade::BuildTargetList(): StartLoop() error!"
00724         //         << G4endl;
00725         return;
00726     }
00727 
00728     ClearAndDestroy(&theTargetList);  // clear theTargetList before rebuilding
00729 
00730     G4Nucleon * nucleon;
00731     G4ParticleDefinition * definition;
00732     G4ThreeVector pos;
00733     G4LorentzVector mom;
00734     // if there are nucleon hit by higher energy models, then SUM(momenta) != 0
00735     initialZ=the3DNucleus->GetCharge();
00736     initialA=the3DNucleus->GetMassNumber();
00737     initial_nuclear_mass=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(initialZ,initialA);
00738     theInitial4Mom = G4LorentzVector(0,0,0,initial_nuclear_mass);
00739     currentA=0;
00740     currentZ=0;
00741     while((nucleon = the3DNucleus->GetNextNucleon()) != NULL)
00742     {
00743         // check if nucleon is hit by higher energy model.
00744         if ( ! nucleon->AreYouHit() )
00745         {
00746             definition = nucleon->GetDefinition();
00747             pos = nucleon->GetPosition();
00748             mom = nucleon->GetMomentum();
00749             //    G4cout << "Nucleus " << pos.mag()/fermi << " " << mom.e() << G4endl;
00750             //theInitial4Mom += mom;
00751             //        the potential inside the nucleus is taken into account, and nucleons are on mass shell.
00752             mom.setE( std::sqrt( mom.vect().mag2() + sqr(definition->GetPDGMass()) ) );
00753             G4KineticTrack * kt = new G4KineticTrack(definition, 0., pos, mom);
00754             kt->SetState(G4KineticTrack::inside);
00755             kt->SetNucleon(nucleon);
00756             theTargetList.push_back(kt);
00757             ++currentA;
00758             if (definition->GetPDGCharge() > .5 ) ++currentZ;
00759         }
00760 
00761 #ifdef debug_BIC_BuildTargetList
00762         else { G4cout << "nucleon is hit" << nucleon << G4endl;}
00763 #endif
00764     }
00765     massInNucleus = 0;
00766     if(currentZ>.5)
00767     {
00768         massInNucleus = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(currentZ,currentA);
00769     } else if (currentZ==0 && currentA>=1 )
00770     {
00771         massInNucleus = currentA * G4Neutron::Neutron()->GetPDGMass();
00772     } else
00773     {
00774         G4cerr << "G4BinaryCascade::BuildTargetList(): Fatal Error - invalid nucleus (A,Z)=("
00775                 << currentA << "," << currentZ << ")" << G4endl;
00776         throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCasacde::BuildTargetList()");
00777     }
00778     currentInitialEnergy=       theInitial4Mom.e() + theProjectile4Momentum.e();
00779 
00780 #ifdef debug_BIC_BuildTargetList
00781     G4cout << "G4BinaryCascade::BuildTargetList():  nucleus (A,Z)=("
00782             << currentA << "," << currentZ << ") mass: " << massInNucleus <<
00783             ", theInitial4Mom " << theInitial4Mom <<
00784             ", currentInitialEnergy " << currentInitialEnergy << G4endl;
00785 #endif
00786 
00787 }
00788 
00789 //----------------------------------------------------------------------------
00790 G4bool  G4BinaryCascade::BuildLateParticleCollisions(G4KineticTrackVector * secondaries)
00791 //----------------------------------------------------------------------------
00792 {
00793    G4bool success(false);
00794    std::vector<G4KineticTrack *>::iterator iter;
00795 
00796    lateA=lateZ=0;
00797    projectileA=projectileZ=0;
00798 
00799    G4double StartingTime=DBL_MAX;        // Search for minimal formation time
00800    for(iter = secondaries->begin(); iter != secondaries->end(); ++iter)
00801    {
00802       if((*iter)->GetFormationTime() < StartingTime)
00803          StartingTime = (*iter)->GetFormationTime();
00804    }
00805 
00806    //PrintKTVector(secondaries, "initial late particles ");
00807    G4LorentzVector lateParticles4Momentum(0,0,0,0);
00808    for(iter = secondaries->begin(); iter != secondaries->end(); ++iter)
00809    {
00810       //  G4cout << " Formation time : " << (*iter)->GetDefinition()->GetParticleName() << " "
00811       //   << (*iter)->GetFormationTime() << G4endl;
00812       G4double FormTime = (*iter)->GetFormationTime() - StartingTime;
00813       (*iter)->SetFormationTime(FormTime);
00814       if( (*iter)->GetState() == G4KineticTrack::undefined  )   // particles from high energy generator
00815       {
00816          FindLateParticleCollision(*iter);
00817          lateParticles4Momentum += (*iter)->GetTrackingMomentum();
00818          lateA += (*iter)->GetDefinition()->GetBaryonNumber();
00819          lateZ += G4lrint((*iter)->GetDefinition()->GetPDGCharge()/eplus);
00820          //PrintKTVector(*iter, "late particle ");
00821       } else
00822       {
00823          theSecondaryList.push_back(*iter);
00824          //PrintKTVector(*iter, "incoming particle ");
00825          theProjectile4Momentum += (*iter)->GetTrackingMomentum();
00826          projectileA += (*iter)->GetDefinition()->GetBaryonNumber();
00827          projectileZ += G4lrint((*iter)->GetDefinition()->GetPDGCharge()/eplus);
00828 #ifdef debug_BIC_Propagate
00829          G4cout << " Adding initial secondary " << *iter
00830                << " time" << (*iter)->GetFormationTime()
00831                << ", state " << (*iter)->GetState() << G4endl;
00832 #endif
00833       }
00834    }
00835    //theCollisionMgr->Print();
00836 
00837    const G4HadProjectile * primary = GetPrimaryProjectile();  // check for primary from TheoHE model
00838    if (primary){
00839       G4LorentzVector mom=primary->Get4Momentum();
00840       theProjectile4Momentum += mom;
00841       projectileA = primary->GetDefinition()->GetBaryonNumber();
00842       projectileZ = G4lrint(primary->GetDefinition()->GetPDGCharge()/eplus);
00843       // now check if "excitation" energy left by TheoHE model
00844       G4double excitation= theProjectile4Momentum.e() + initial_nuclear_mass - lateParticles4Momentum.e() - massInNucleus;
00845 #ifdef debug_BIC_GetExcitationEnergy
00846       G4cout << "BIC: Proj.e, nucl initial, nucl final, lateParticles"
00847             << theProjectile4Momentum << ",  "
00848             << initial_nuclear_mass<< ",  " << massInNucleus << ",  "
00849             << lateParticles4Momentum << G4endl;
00850       G4cout << "BIC: Proj.e / initial excitation: " << theProjectile4Momentum.e() << " / " << excitation << G4endl;
00851 #endif
00852       success = excitation > 0;
00853 #ifdef debug_G4BinaryCascade
00854       if ( ! success ) {
00855          G4cout << "G4BinaryCascade::BuildLateParticleCollisions(): Proj.e / initial excitation: " << theProjectile4Momentum.e() << " / " << excitation << G4endl;
00856          //PrintKTVector(secondaries);
00857       }
00858 #endif
00859    } else {
00860       // no primary from HE model -> cascade
00861       success=true;
00862    }
00863 
00864    if (success) {
00865       secondaries->clear(); // Don't leave "G4KineticTrack *"s in two vectors
00866       delete secondaries;
00867    }
00868    return success;
00869 }
00870 
00871 //----------------------------------------------------------------------------
00872 G4ReactionProductVector *  G4BinaryCascade::DeExcite()
00873 //----------------------------------------------------------------------------
00874 {
00875    // find a fragment and call the precompound model.
00876    G4Fragment * fragment = 0;
00877    G4ReactionProductVector * precompoundProducts = 0;
00878 
00879    G4LorentzVector pFragment(0);
00880    // G4cout << " final4mon " << GetFinal4Momentum() /MeV << G4endl;
00881 
00882    //   if ( ExcitationEnergy >= 0 )                                         // closed by Uzhi
00883    //   {                                                                    // closed by Uzhi
00884    fragment = FindFragments();
00885    if(fragment)                                                            // Uzhi
00886    {                                                                       // Uzhi
00887       if(fragment->GetA() >1.5)                                          // Uzhi
00888       {
00889          pFragment=fragment->GetMomentum();
00890          // G4cout << " going to preco with fragment 4 mom " << pFragment << G4endl;
00891          if (theDeExcitation)                // pre-compound
00892          {
00893             precompoundProducts= theDeExcitation->DeExcite(*fragment);
00894          }
00895          else if (theExcitationHandler)    // de-excitation
00896          {
00897             precompoundProducts=theExcitationHandler->BreakItUp(*fragment);
00898          }
00899 
00900       } else
00901       {                                   // fragment->GetA() < 1.5, so a single proton, as a fragment must have Z>0
00902          if (theTargetList.size() + theCapturedList.size() > 1 ) {
00903             throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCasacde:: Invalid Fragment");
00904          }
00905 
00906          std::vector<G4KineticTrack *>::iterator i;
00907          if ( theTargetList.size() == 1 )  {i=theTargetList.begin();}
00908          if ( theCapturedList.size() == 1 ) {i=theCapturedList.begin();}                             // Uzhi
00909          G4ReactionProduct * aNew = new G4ReactionProduct((*i)->GetDefinition());
00910          aNew->SetTotalEnergy((*i)->GetDefinition()->GetPDGMass());
00911          aNew->SetMomentum(G4ThreeVector(0));// see boost for preCompoundProducts below..
00912          precompoundProducts = new G4ReactionProductVector();
00913          precompoundProducts->push_back(aNew);
00914       }                            // End of fragment->GetA() < 1.5
00915       delete fragment;
00916       fragment=0;
00917 
00918    } else                            // End of if(fragment)
00919    {                                 // No fragment, can be neutrons only  // Uzhi
00920 
00921       precompoundProducts = DecayVoidNucleus();
00922    }
00923    return precompoundProducts;
00924 }
00925 
00926 //----------------------------------------------------------------------------
00927 G4ReactionProductVector *  G4BinaryCascade::DecayVoidNucleus()
00928 //----------------------------------------------------------------------------
00929 {
00930    G4ReactionProductVector * result=0;
00931    if ( (theTargetList.size()+theCapturedList.size()) > 0 )
00932    {
00933       result = new G4ReactionProductVector;
00934       std::vector<G4KineticTrack *>::iterator aNuc;
00935       G4LorentzVector aVec;
00936       std::vector<G4double> masses;
00937       G4double sumMass(0);
00938 
00939       if ( theTargetList.size() != 0)                                      // Uzhi
00940             {
00941          for ( aNuc=theTargetList.begin(); aNuc != theTargetList.end(); aNuc++)
00942          {
00943             G4double mass=(*aNuc)->GetDefinition()->GetPDGMass();
00944             masses.push_back(mass);
00945             sumMass += mass;
00946          }
00947             }                                                                    // Uzhi
00948 
00949       if ( theCapturedList.size() != 0)                                    // Uzhi
00950       {                                                                    // Uzhi
00951          for(aNuc = theCapturedList.begin();                               // Uzhi
00952                aNuc != theCapturedList.end(); aNuc++)                        // Uzhi
00953          {                                                                 // Uzhi
00954             G4double mass=(*aNuc)->GetDefinition()->GetPDGMass();          // Uzhi
00955             masses.push_back(mass);                                        // Uzhi
00956             sumMass += mass;                                               // Uzhi
00957          }
00958       }
00959 
00960       G4LorentzVector finalP=GetFinal4Momentum();
00961       G4FermiPhaseSpaceDecay decay;
00962       // G4cout << " some neutrons? " << masses.size() <<" " ;
00963       // G4cout<< theTargetList.size()<<" "<<finalP <<" " << finalP.mag()<<G4endl;
00964 
00965       G4double eCMS=finalP.mag();
00966       if ( eCMS < sumMass )                    // @@GF --- Cheat!!
00967       {
00968          eCMS=sumMass + 2*MeV*masses.size();
00969          finalP.setE(std::sqrt(finalP.vect().mag2() + sqr(eCMS)));
00970       }
00971 
00972       precompoundLorentzboost.set(finalP.boostVector());
00973       std::vector<G4LorentzVector*> * momenta=decay.Decay(eCMS,masses);
00974       std::vector<G4LorentzVector*>::iterator aMom=momenta->begin();
00975 
00976 
00977       if ( theTargetList.size() != 0)
00978       {
00979          for ( aNuc=theTargetList.begin();
00980                (aNuc != theTargetList.end()) && (aMom!=momenta->end());
00981                aNuc++, aMom++ )
00982          {
00983             G4ReactionProduct * aNew = new G4ReactionProduct((*aNuc)->GetDefinition());
00984             aNew->SetTotalEnergy((*aMom)->e());
00985             aNew->SetMomentum((*aMom)->vect());
00986             result->push_back(aNew);
00987 
00988             delete *aMom;
00989          }
00990       }
00991 
00992       if ( theCapturedList.size() != 0)                                    // Uzhi
00993             {                                                                    // Uzhi
00994          for ( aNuc=theCapturedList.begin();                                // Uzhi
00995                (aNuc != theCapturedList.end()) && (aMom!=momenta->end());    // Uzhi
00996                aNuc++, aMom++ )                                             // Uzhi
00997          {                                                                  // Uzhi
00998             G4ReactionProduct * aNew = new G4ReactionProduct(               // Uzhi
00999                   (*aNuc)->GetDefinition());            // Uzhi
01000             aNew->SetTotalEnergy((*aMom)->e());                             // Uzhi
01001             aNew->SetMomentum((*aMom)->vect());                             // Uzhi
01002             result->push_back(aNew);                           // Uzhi
01003             delete *aMom;                                                   // Uzhi
01004          }                                                                  // Uzhi
01005             }                                                                    // Uzhi
01006 
01007       delete momenta;
01008    }
01009    return result;
01010 }                   // End if(!fragment)
01011 
01012 //----------------------------------------------------------------------------
01013 G4ReactionProductVector * G4BinaryCascade::ProductsAddFinalState(G4ReactionProductVector * products, G4KineticTrackVector & fs)
01014 //----------------------------------------------------------------------------
01015 {
01016 // fill in products the outgoing particles
01017     size_t i(0);
01018     for(i = 0; i< fs.size(); i++)
01019     {
01020         G4KineticTrack * kt = fs[i];
01021         G4ReactionProduct * aNew = new G4ReactionProduct(kt->GetDefinition());
01022         aNew->SetMomentum(kt->Get4Momentum().vect());
01023         aNew->SetTotalEnergy(kt->Get4Momentum().e());
01024         aNew->SetNewlyAdded(kt->IsParticipant());
01025         products->push_back(aNew);
01026 
01027 #ifdef debug_BIC_Propagate_finals
01028         G4cout <<kt->GetDefinition()->GetParticleName();
01029         G4cout << " Particle Ekin " << aNew->GetKineticEnergy();
01030         G4cout << "final is " << kt->GetDefinition()->GetPDGStable() ? "stable" :
01031                 ( kt->GetDefinition()->IsShortLived() ? "short lived " : "non stable") << G4endl;;
01032 #endif
01033 
01034     }
01035     return products;
01036 }
01037 //----------------------------------------------------------------------------
01038 G4ReactionProductVector * G4BinaryCascade::ProductsAddPrecompound(G4ReactionProductVector * products, G4ReactionProductVector * precompoundProducts)
01039 //----------------------------------------------------------------------------
01040 {
01041    G4LorentzVector pSumPreco(0), pPreco(0);
01042 
01043    if ( precompoundProducts )
01044    {
01045       std::vector<G4ReactionProduct *>::iterator j;
01046       for(j = precompoundProducts->begin(); j != precompoundProducts->end(); ++j)
01047       {
01048          // boost back to system of moving nucleus
01049          G4LorentzVector pProduct((*j)->GetMomentum(),(*j)->GetTotalEnergy());
01050          pPreco+= pProduct;
01051 #ifdef debug_BIC_Propagate_finals
01052          G4cout << " pProduct be4 boost " <<pProduct << G4endl;
01053 #endif
01054          pProduct *= precompoundLorentzboost;
01055 #ifdef debug_BIC_Propagate_finals
01056          G4cout << " pProduct aft boost " <<pProduct << G4endl;
01057 #endif
01058          pSumPreco += pProduct;
01059          (*j)->SetTotalEnergy(pProduct.e());
01060          (*j)->SetMomentum(pProduct.vect());
01061          (*j)->SetNewlyAdded(true);
01062          products->push_back(*j);
01063       }
01064       // G4cout << " unboosted preco result mom " << pPreco / MeV << "  ..- fragmentMom " << (pPreco - pFragment)/MeV<< G4endl;
01065       // G4cout << " preco result mom " << pSumPreco / MeV << "  ..-file4Mom " << (pSumPreco - GetFinal4Momentum())/MeV<< G4endl;
01066       precompoundProducts->clear();
01067       delete precompoundProducts;
01068    }
01069    return products;
01070 }
01071 //----------------------------------------------------------------------------
01072 void  G4BinaryCascade::FindCollisions(G4KineticTrackVector * secondaries)
01073 //----------------------------------------------------------------------------
01074 {
01075     for(std::vector<G4KineticTrack *>::iterator i = secondaries->begin();
01076             i != secondaries->end(); ++i)
01077     {
01078 
01079         for(std::vector<G4BCAction *>::iterator j = theImR.begin();
01080                 j!=theImR.end(); j++)
01081         {
01082             //      G4cout << "G4BinaryCascade::FindCollisions: at action " << *j << G4endl;
01083             const std::vector<G4CollisionInitialState *> & aCandList
01084             = (*j)->GetCollisions(*i, theTargetList, theCurrentTime);
01085             for(size_t count=0; count<aCandList.size(); count++)
01086             {
01087                 theCollisionMgr->AddCollision(aCandList[count]);
01088                 //4cout << "====================== New Collision ================="<<G4endl;
01089                 //theCollisionMgr->Print();
01090             }
01091         }
01092     }
01093 }
01094 
01095 
01096 //----------------------------------------------------------------------------
01097 void  G4BinaryCascade::FindDecayCollision(G4KineticTrack * secondary)
01098 //----------------------------------------------------------------------------
01099 {
01100     const std::vector<G4CollisionInitialState *> & aCandList
01101     = theDecay->GetCollisions(secondary, theTargetList, theCurrentTime);
01102     for(size_t count=0; count<aCandList.size(); count++)
01103     {
01104         theCollisionMgr->AddCollision(aCandList[count]);
01105     }
01106 }
01107 
01108 //----------------------------------------------------------------------------
01109 void  G4BinaryCascade::FindLateParticleCollision(G4KineticTrack * secondary)
01110 //----------------------------------------------------------------------------
01111 {
01112 
01113     G4double tin=0., tout=0.;
01114     if (((G4RKPropagation*)thePropagator)->GetSphereIntersectionTimes(secondary,tin,tout))
01115     {
01116         if ( tin > 0 )
01117         {
01118             secondary->SetState(G4KineticTrack::outside);
01119         } else if ( tout > 0 )
01120         {
01121             secondary->SetState(G4KineticTrack::inside);
01122         } else {
01123             //G4cout << "G4BC set miss , tin, tout " << tin << " , " << tout <<G4endl;
01124             secondary->SetState(G4KineticTrack::miss_nucleus);
01125         }
01126     } else {
01127         secondary->SetState(G4KineticTrack::miss_nucleus);
01128         //G4cout << "G4BC set miss ,no intersect tin, tout " << tin << " , " << tout <<G4endl;
01129     }
01130 
01131 
01132 #ifdef debug_BIC_FindCollision
01133     G4cout << "FindLateP Particle, 4-mom, times newState "
01134             << secondary->GetDefinition()->GetParticleName() << " "
01135             << secondary->Get4Momentum()
01136             << " times " <<  tin << " " << tout << " "
01137             << secondary->GetState() << G4endl;
01138 #endif
01139 
01140     const std::vector<G4CollisionInitialState *> & aCandList
01141     = theLateParticle->GetCollisions(secondary, theTargetList, theCurrentTime);
01142     for(size_t count=0; count<aCandList.size(); count++)
01143     {
01144 #ifdef debug_BIC_FindCollision
01145         G4cout << " Adding a late Col : " << aCandList[count] << G4endl;
01146 #endif
01147         theCollisionMgr->AddCollision(aCandList[count]);
01148     }
01149 }
01150 
01151 
01152 //----------------------------------------------------------------------------
01153 G4bool G4BinaryCascade::ApplyCollision(G4CollisionInitialState * collision)
01154 //----------------------------------------------------------------------------
01155 {
01156     G4KineticTrack * primary = collision->GetPrimary();
01157 
01158 #ifdef debug_BIC_ApplyCollision
01159     G4cerr << "G4BinaryCascade::ApplyCollision start"<<G4endl;
01160     theCollisionMgr->Print();
01161     G4cout << "ApplyCollisions : projte 4mom " << primary->GetTrackingMomentum()<< G4endl;
01162 #endif
01163 
01164     G4KineticTrackVector target_collection=collision->GetTargetCollection();
01165     G4bool haveTarget=target_collection.size()>0;
01166     if( haveTarget && (primary->GetState() != G4KineticTrack::inside) )
01167     {
01168 #ifdef debug_G4BinaryCascade
01169         G4cout << "G4BinaryCasacde::ApplyCollision(): StateError " << primary << G4endl;
01170         PrintKTVector(primary,std::string("primay- ..."));
01171         PrintKTVector(&target_collection,std::string("... targets"));
01172         collision->Print();
01173         G4cout << G4endl;
01174         theCollisionMgr->Print();
01175         //*GF*     throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCasacde::ApplyCollision()");
01176 #endif
01177         return false;
01178 //    } else {
01179 //       G4cout << "G4BinaryCasacde::ApplyCollision(): decay " << G4endl;
01180 //       PrintKTVector(primary,std::string("primay- ..."));
01181 //       G4double tin=0., tout=0.;
01182 //       if (((G4RKPropagation*)thePropagator)->GetSphereIntersectionTimes(primary,tin,tout))
01183 //       {
01184 //           G4cout << "tin tout: " << tin << " " << tout << G4endl;
01185 //       }
01186 
01187     }
01188 
01189     G4LorentzVector mom4Primary=primary->Get4Momentum();
01190 
01191     G4int initialBaryon(0);
01192     G4int initialCharge(0);
01193     if (primary->GetState() == G4KineticTrack::inside)
01194     {
01195         initialBaryon = primary->GetDefinition()->GetBaryonNumber();
01196         initialCharge = G4lrint(primary->GetDefinition()->GetPDGCharge()/eplus);
01197     }
01198 
01199     // for primary resonances, subtract neutron ( = proton) field ( ie. add std::abs(field))
01200     G4double initial_Efermi=CorrectShortlivedPrimaryForFermi(primary,target_collection);
01201     //****************************************
01202 
01203 
01204     G4KineticTrackVector * products = collision->GetFinalState();
01205 
01206 #ifdef debug_BIC_ApplyCollision
01207     DebugApplyCollisionFail(collision, products);
01208 #endif
01209 
01210     // reset primary to initial state, in case there is a veto...
01211     primary->Set4Momentum(mom4Primary);
01212 
01213     G4bool lateParticleCollision= (!haveTarget) && products && products->size() == 1;
01214     G4bool decayCollision= (!haveTarget) && products && products->size() > 1;
01215     G4bool Success(true);
01216 
01217 
01218 #ifdef debug_G4BinaryCascade
01219     G4int lateBaryon(0), lateCharge(0);
01220 #endif
01221 
01222     if ( lateParticleCollision )
01223     {  // for late particles, reset charges
01224         //G4cout << "lateP, initial B C state " << initialBaryon << " "
01225         //        << initialCharge<< " " << primary->GetState() << " "<< primary->GetDefinition()->GetParticleName()<< G4endl;
01226 #ifdef debug_G4BinaryCascade
01227         lateBaryon = initialBaryon;
01228         lateCharge = initialCharge;
01229 #endif
01230         initialBaryon=initialCharge=0;
01231         lateA -= primary->GetDefinition()->GetBaryonNumber();
01232         lateZ -= G4lrint(primary->GetDefinition()->GetPDGCharge()/eplus);
01233     }
01234 
01235     initialBaryon += collision->GetTargetBaryonNumber();
01236     initialCharge += G4lrint(collision->GetTargetCharge());
01237     if (!lateParticleCollision)
01238     {
01239        if( !products || products->size()==0 || !CheckPauliPrinciple(products) )
01240        {
01241 #ifdef debug_BIC_ApplyCollision
01242           if (products) G4cout << " ======Failed Pauli =====" << G4endl;
01243           G4cerr << "G4BinaryCascade::ApplyCollision blocked"<<G4endl;
01244 #endif
01245           Success=false;
01246        }
01247 
01248 
01249 
01250        if (Success && primary->GetState() == G4KineticTrack::inside ) {   // if the primary was outside, nothing to correct
01251           if (! CorrectShortlivedFinalsForFermi(products, initial_Efermi)){
01252              Success=false;
01253           }
01254        }
01255     }
01256 
01257 #ifdef debug_BIC_ApplyCollision
01258     DebugApplyCollision(collision, products);
01259 #endif
01260 
01261     if ( ! Success ){
01262         if (products) ClearAndDestroy(products);
01263         if ( decayCollision ) FindDecayCollision(primary);  // for decay, sample new decay
01264         delete products;
01265         products=0;
01266         return false;
01267     }
01268 
01269     G4int finalBaryon(0);
01270     G4int finalCharge(0);
01271     G4KineticTrackVector toFinalState;
01272     for(std::vector<G4KineticTrack *>::iterator i =products->begin(); i != products->end(); i++)
01273     {
01274         if ( ! lateParticleCollision )
01275         {
01276             (*i)->SetState(primary->GetState());  // decay may be anywhere!
01277             if ( (*i)->GetState() == G4KineticTrack::inside ){
01278                 finalBaryon+=(*i)->GetDefinition()->GetBaryonNumber();
01279                 finalCharge+=G4lrint((*i)->GetDefinition()->GetPDGCharge()/eplus);
01280             } else {
01281                G4double tin=0., tout=0.;
01282                if (((G4RKPropagation*)thePropagator)->GetSphereIntersectionTimes((*i),tin,tout) &&
01283                      tin < 0 && tout > 0 )
01284                {
01285                   PrintKTVector((*i),"particle inside marked not-inside");
01286                    G4cout << "tin tout: " << tin << " " << tout << G4endl;
01287                }
01288             }
01289         } else {
01290             G4double tin=0., tout=0.;
01291             if (((G4RKPropagation*)thePropagator)->GetSphereIntersectionTimes((*i),tin,tout))
01292             {
01293                 //G4cout << "tin tout: " << tin << " " << tout << G4endl;
01294                 if ( tin > 0 )
01295                 {
01296                     (*i)->SetState(G4KineticTrack::outside);
01297                 }
01298                 else if ( tout > 0 )
01299                 {
01300                     (*i)->SetState(G4KineticTrack::inside);
01301                     finalBaryon+=(*i)->GetDefinition()->GetBaryonNumber();
01302                     finalCharge+=G4lrint((*i)->GetDefinition()->GetPDGCharge()/eplus);
01303                 }
01304                 else
01305                 {
01306                     (*i)->SetState(G4KineticTrack::gone_out);
01307                     toFinalState.push_back((*i));
01308                 }
01309             } else
01310             {
01311                 (*i)->SetState(G4KineticTrack::miss_nucleus);
01312                 //G4cout << " G4BC - miss -late Part- no intersection found " << G4endl;
01313                 toFinalState.push_back((*i));
01314             }
01315 
01316         }
01317     }
01318     if(!toFinalState.empty())
01319     {
01320         theFinalState.insert(theFinalState.end(),
01321                 toFinalState.begin(),toFinalState.end());
01322         std::vector<G4KineticTrack *>::iterator iter1, iter2;
01323         for(iter1 = toFinalState.begin(); iter1 != toFinalState.end();
01324                 ++iter1)
01325         {
01326             iter2 = std::find(products->begin(), products->end(),
01327                     *iter1);
01328             if ( iter2 != products->end() ) products->erase(iter2);
01329         }
01330         theCollisionMgr->RemoveTracksCollisions(&toFinalState);
01331     }
01332 
01333     //G4cout << " currentA, Z be4: " << currentA << " " << currentZ << G4endl;
01334     currentA += finalBaryon-initialBaryon;
01335     currentZ += finalCharge-initialCharge;
01336     //G4cout << " ApplyCollision currentA, Z aft: " << currentA << " " << currentZ << G4endl;
01337 
01338     G4KineticTrackVector oldSecondaries;
01339     oldSecondaries.push_back(primary);
01340     primary->Hit();
01341 
01342 #ifdef debug_G4BinaryCascade
01343     if ( (finalBaryon-initialBaryon-lateBaryon) != 0 || (finalCharge-initialCharge-lateCharge) != 0 )
01344     {
01345         G4cout << "G4BinaryCascade: Error in Balancing: " << G4endl;
01346         G4cout << "initial/final baryon number, initial/final Charge "
01347                 << initialBaryon <<" "<< finalBaryon <<" "
01348                 << initialCharge <<" "<< finalCharge <<" "
01349                 << " in Collision type: "<< typeid(*collision->GetGenerator()).name()
01350                 << ", with number of products: "<< products->size() <<G4endl;
01351         G4cout << G4endl<<"Initial condition are these:"<<G4endl;
01352         G4cout << "proj: "<<collision->GetPrimary()->GetDefinition()->GetParticleName()<<G4endl;
01353         for(size_t it=0; it<collision->GetTargetCollection().size(); it++)
01354         {
01355             G4cout << "targ: "
01356                     <<collision->GetTargetCollection()[it]->GetDefinition()->GetParticleName()<<G4endl;
01357         }
01358         PrintKTVector(&collision->GetTargetCollection(),std::string(" Target particles"));
01359         G4cout << G4endl<<G4endl;
01360     }
01361 #endif
01362 
01363     G4KineticTrackVector oldTarget = collision->GetTargetCollection();
01364     for(size_t ii=0; ii< oldTarget.size(); ii++)
01365     {
01366         oldTarget[ii]->Hit();
01367     }
01368 
01369     UpdateTracksAndCollisions(&oldSecondaries, &oldTarget, products);
01370     std::for_each(oldSecondaries.begin(), oldSecondaries.end(), Delete<G4KineticTrack>());
01371     std::for_each(oldTarget.begin(), oldTarget.end(), Delete<G4KineticTrack>());
01372 
01373     delete products;
01374     return true;
01375 }
01376 
01377 
01378 //----------------------------------------------------------------------------
01379 G4bool G4BinaryCascade::Absorb()
01380 //----------------------------------------------------------------------------
01381 {
01382     // Do it in two step: cannot change theSecondaryList inside the first loop.
01383     G4Absorber absorber(theCutOnPAbsorb);
01384 
01385     // Build the vector of G4KineticTracks that must be absorbed
01386     G4KineticTrackVector absorbList;
01387     std::vector<G4KineticTrack *>::iterator iter;
01388     //  PrintKTVector(&theSecondaryList, " testing for Absorb" );
01389     for(iter = theSecondaryList.begin();
01390             iter != theSecondaryList.end(); ++iter)
01391     {
01392         G4KineticTrack * kt = *iter;
01393         if(kt->GetState() == G4KineticTrack::inside)// absorption happens only inside the nucleus
01394         {
01395             if(absorber.WillBeAbsorbed(*kt))
01396             {
01397                 absorbList.push_back(kt);
01398             }
01399         }
01400     }
01401 
01402     if(absorbList.empty())
01403         return false;
01404 
01405     G4KineticTrackVector toDelete;
01406     for(iter = absorbList.begin(); iter != absorbList.end(); ++iter)
01407     {
01408         G4KineticTrack * kt = *iter;
01409         if(!absorber.FindAbsorbers(*kt, theTargetList))
01410             throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade::Absorb(): Cannot absorb a particle.");
01411 
01412         if(!absorber.FindProducts(*kt))
01413             throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade::Absorb(): Cannot absorb a particle.");
01414 
01415         G4KineticTrackVector * products = absorber.GetProducts();
01416         // ------ debug
01417         G4int count = 0;
01418         // ------ end debug
01419         while(!CheckPauliPrinciple(products))
01420         {
01421             // ------ debug
01422             count++;
01423             // ------ end debug
01424             ClearAndDestroy(products);
01425             if(!absorber.FindProducts(*kt))
01426                 throw G4HadronicException(__FILE__, __LINE__,
01427                         "G4BinaryCascade::Absorb(): Cannot absorb a particle.");
01428         }
01429         // ------ debug
01430         //    G4cerr << "Absorb CheckPauliPrinciple count= " <<  count << G4endl;
01431         // ------ end debug
01432         G4KineticTrackVector toRemove;  // build a vector for UpdateTrack...
01433         toRemove.push_back(kt);
01434         toDelete.push_back(kt);  // delete the track later
01435         G4KineticTrackVector * absorbers = absorber.GetAbsorbers();
01436         UpdateTracksAndCollisions(&toRemove, absorbers, products);
01437         ClearAndDestroy(absorbers);
01438     }
01439     ClearAndDestroy(&toDelete);
01440     return true;
01441 }
01442 
01443 
01444 
01445 // Capture all p and n with Energy < theCutOnP
01446 //----------------------------------------------------------------------------
01447 G4bool G4BinaryCascade::Capture(G4bool verbose)
01448 //----------------------------------------------------------------------------
01449 {
01450     G4KineticTrackVector captured;
01451     G4bool capture = false;
01452     std::vector<G4KineticTrack *>::iterator i;
01453 
01454     G4RKPropagation * RKprop=(G4RKPropagation *)thePropagator;
01455 
01456     G4double capturedEnergy = 0;
01457     G4int particlesAboveCut=0;
01458     G4int particlesBelowCut=0;
01459     if ( verbose ) G4cout << " Capture: secondaries " << theSecondaryList.size() << G4endl;
01460     for(i = theSecondaryList.begin(); i != theSecondaryList.end(); ++i)
01461     {
01462         G4KineticTrack * kt = *i;
01463         if (verbose) G4cout << "Capture position, radius, state " <<kt->GetPosition().mag()<<" "<<theOuterRadius<<" "<<kt->GetState()<<G4endl;
01464         if(kt->GetState() == G4KineticTrack::inside) // capture happens only inside the nucleus
01465         {
01466             if((kt->GetDefinition() == G4Proton::Proton()) ||
01467                     (kt->GetDefinition() == G4Neutron::Neutron()))
01468             {
01469                 //GF cut on kinetic energy    if(kt->Get4Momentum().vect().mag() >= theCutOnP)
01470                 G4double field=RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition())
01471                                            -RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding());
01472                 G4double energy= kt->Get4Momentum().e() - kt->GetActualMass() + field;
01473                 if (verbose ) G4cout << "Capture: .e(), mass, field, energy" << kt->Get4Momentum().e() <<" "<<kt->GetActualMass()<<" "<<field<<" "<<energy<< G4endl;
01474                 //       if( energy < theCutOnP )
01475                 //       {
01476                 capturedEnergy+=energy;
01477                 ++particlesBelowCut;
01478                 //       } else
01479                 //       {
01480                 //          ++particlesAboveCut;
01481                 //       }
01482             }
01483         }
01484     }
01485     if (verbose) G4cout << "Capture particlesAboveCut,particlesBelowCut, capturedEnergy,capturedEnergy/particlesBelowCut <? 0.2*theCutOnP "
01486             << particlesAboveCut << " " << particlesBelowCut << " " << capturedEnergy
01487             << " "  << G4endl;
01488 //    << " " << (particlesBelowCut>0) ? (capturedEnergy/particlesBelowCut) : (capturedEnergy) << " " << 0.2*theCutOnP << G4endl;
01489     //  if(particlesAboveCut==0 && particlesBelowCut>0 && capturedEnergy/particlesBelowCut<0.2*theCutOnP)
01490     if(particlesBelowCut>0 && capturedEnergy/particlesBelowCut<0.2*theCutOnP)
01491     {
01492         capture=true;
01493         for(i = theSecondaryList.begin(); i != theSecondaryList.end(); ++i)
01494         {
01495             G4KineticTrack * kt = *i;
01496             if(kt->GetState() == G4KineticTrack::inside) // capture happens only inside the nucleus
01497             {
01498                 if((kt->GetDefinition() == G4Proton::Proton()) ||
01499                         (kt->GetDefinition() == G4Neutron::Neutron()))
01500                 {
01501                     captured.push_back(kt);
01502                     kt->Hit();                          //
01503                     theCapturedList.push_back(kt);
01504                 }
01505             }
01506         }
01507         UpdateTracksAndCollisions(&captured, NULL, NULL);
01508     }
01509 
01510     return capture;
01511 }
01512 
01513 
01514 //----------------------------------------------------------------------------
01515 G4bool G4BinaryCascade::CheckPauliPrinciple(G4KineticTrackVector * products)
01516 //----------------------------------------------------------------------------
01517 {
01518     G4int A = the3DNucleus->GetMassNumber();
01519     G4int Z = the3DNucleus->GetCharge();
01520 
01521     G4FermiMomentum fermiMom;
01522     fermiMom.Init(A, Z);
01523 
01524     const G4VNuclearDensity * density=the3DNucleus->GetNuclearDensity();
01525 
01526     G4KineticTrackVector::iterator i;
01527     G4ParticleDefinition * definition;
01528 
01529     // ------ debug
01530     G4bool myflag = true;
01531     // ------ end debug
01532     //  G4ThreeVector xpos(0);
01533     for(i = products->begin(); i != products->end(); ++i)
01534     {
01535         definition = (*i)->GetDefinition();
01536         if((definition == G4Proton::Proton()) ||
01537                 (definition == G4Neutron::Neutron()))
01538         {
01539             G4ThreeVector pos = (*i)->GetPosition();
01540             G4double d = density->GetDensity(pos);
01541             // energy correspondiing to fermi momentum
01542             G4double eFermi = std::sqrt( sqr(fermiMom.GetFermiMomentum(d)) + (*i)->Get4Momentum().mag2() );
01543             if( definition == G4Proton::Proton() )
01544             {
01545                 eFermi -= the3DNucleus->CoulombBarrier();
01546             }
01547             G4LorentzVector mom = (*i)->Get4Momentum();
01548             // ------ debug
01549             /*
01550              *        G4cout << "p:[" << (1/MeV)*mom.x() << " " << (1/MeV)*mom.y() << " "
01551              *            << (1/MeV)*mom.z() << "] |p3|:"
01552              *            << (1/MeV)*mom.vect().mag() << " E:" << (1/MeV)*mom.t() << " ] m: "
01553              *            << (1/MeV)*mom.mag() << "  pos["
01554              *            << (1/fermi)*pos.x() << " "<< (1/fermi)*pos.y() << " "
01555              *            << (1/fermi)*pos.z() << "] |Dpos|: "
01556              *            << (1/fermi)*(pos-xpos).mag() << " Pfermi:"
01557              *            << (1/MeV)*p << G4endl;
01558              *         xpos=pos;
01559              */
01560             // ------ end debug
01561             if(mom.e() < eFermi )
01562             {
01563                 // ------ debug
01564                 myflag = false;
01565                 // ------ end debug
01566                 //      return false;
01567             }
01568         }
01569     }
01570 #ifdef debug_BIC_CheckPauli
01571     if ( myflag  )
01572     {
01573         for(i = products->begin(); i != products->end(); ++i)
01574         {
01575             definition = (*i)->GetDefinition();
01576             if((definition == G4Proton::Proton()) ||
01577                     (definition == G4Neutron::Neutron()))
01578             {
01579                 G4ThreeVector pos = (*i)->GetPosition();
01580                 G4double d = density->GetDensity(pos);
01581                 G4double pFermi = fermiMom.GetFermiMomentum(d);
01582                 G4LorentzVector mom = (*i)->Get4Momentum();
01583                 G4double field =((G4RKPropagation*)thePropagator)->GetField(definition->GetPDGEncoding(),pos);
01584                 if ( mom.e()-mom.mag()+field > 160*MeV )
01585                 {
01586                     G4cout << "momentum problem pFermi=" <<  pFermi
01587                             << " mom, mom.m " << mom << " " << mom.mag()
01588                             << " field " << field << G4endl;
01589                 }
01590             }
01591         }
01592     }
01593 #endif
01594 
01595     return myflag;
01596 }
01597 
01598 //----------------------------------------------------------------------------
01599 void G4BinaryCascade::StepParticlesOut()
01600 //----------------------------------------------------------------------------
01601 {
01602     G4int counter=0;
01603     G4int countreset=0;
01604     //G4cout << " nucl. Radius " << radius << G4endl;
01605     // G4cerr <<"pre-while- theSecondaryList "<<G4endl;
01606     while( theSecondaryList.size() > 0 )
01607     {
01608         G4int nsec=0;
01609         G4double minTimeStep = 1.e-12*ns;   // about 30*fermi/(0.1*c_light);1.e-12*ns
01610         // i.e. a big step
01611         std::vector<G4KineticTrack *>::iterator i;
01612         for(i = theSecondaryList.begin(); i != theSecondaryList.end(); ++i)
01613         {
01614             G4KineticTrack * kt = *i;
01615             if( kt->GetState() == G4KineticTrack::inside )
01616             {
01617                 nsec++;
01618                 G4double tStep(0), tdummy(0);
01619                 G4bool intersect =
01620                         ((G4RKPropagation*)thePropagator)->GetSphereIntersectionTimes(kt,tdummy,tStep);
01621 #ifdef debug_BIC_StepParticlesOut
01622                 G4cout << " minTimeStep, tStep Particle " <<minTimeStep << " " <<tStep
01623                         << " " <<kt->GetDefinition()->GetParticleName()
01624                         << " 4mom " << kt->GetTrackingMomentum()<<G4endl;
01625                 if ( ! intersect );
01626                 {
01627                     PrintKTVector(&theSecondaryList, std::string(" state ERROR....."));
01628                     throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade::StepParticlesOut() particle not in nucleus");
01629                 }
01630 #endif
01631                 if(intersect && tStep<minTimeStep && tStep> 0 )
01632                 {
01633                     minTimeStep = tStep;
01634                 }
01635             } else if ( kt->GetState() != G4KineticTrack::outside ){
01636                 PrintKTVector(&theSecondaryList, std::string(" state ERROR....."));
01637                 throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade::StepParticlesOut() particle not in nucleus");
01638             }
01639         }
01640         minTimeStep *= 1.2;
01641         //    G4cerr << "CaptureCount = "<<counter<<" "<<nsec<<" "<<minTimeStep<<" "<<1*ns<<G4endl;
01642         G4double timeToCollision=DBL_MAX;
01643         G4CollisionInitialState * nextCollision=0;
01644         if(theCollisionMgr->Entries() > 0)
01645         {
01646             nextCollision = theCollisionMgr->GetNextCollision();
01647             timeToCollision = nextCollision->GetCollisionTime()-theCurrentTime;
01648             G4cout << " NextCollision  * , Time= " << nextCollision << " "
01649                     <<timeToCollision<< G4endl;
01650         }
01651         if ( timeToCollision > minTimeStep )
01652         {
01653             DoTimeStep(minTimeStep);
01654             ++counter;
01655         } else
01656         {
01657             if (!DoTimeStep(timeToCollision) )
01658             {
01659                 // Check if nextCollision is still valid, ie. partcile did not leave nucleus
01660                 if (theCollisionMgr->GetNextCollision() != nextCollision )
01661                 {
01662                     nextCollision = 0;
01663                 }
01664             }
01665             // G4cerr <<"post- DoTimeStep 3"<<G4endl;
01666 
01667             if(nextCollision)
01668             {
01669                 if  ( ApplyCollision(nextCollision))
01670                 {
01671                     // G4cout << "ApplyCollision sucess " << G4endl;
01672                 } else
01673                 {
01674                     theCollisionMgr->RemoveCollision(nextCollision);
01675                 }
01676             }
01677         }
01678 
01679         if(countreset>100)
01680         {
01681 #ifdef debug_G4BinaryCascade
01682             G4cerr << "G4BinaryCascade.cc: Warning - aborting looping particle(s)" << G4endl;
01683             PrintKTVector(&theSecondaryList," looping particles added to theFinalState");
01684 #endif
01685 
01686             //  add left secondaries to FinalSate
01687             std::vector<G4KineticTrack *>::iterator iter;
01688             for ( iter =theSecondaryList.begin(); iter != theSecondaryList.end(); ++iter)
01689             {
01690                 theFinalState.push_back(*iter);
01691             }
01692             theSecondaryList.clear();
01693 
01694             break;
01695         }
01696 
01697         if(Absorb())
01698         {
01699             //       haveProducts = true;
01700             // G4cout << "Absorb sucess " << G4endl;
01701         }
01702 
01703         if(Capture(false))
01704         {
01705             //       haveProducts = true;
01706 #ifdef debug_BIC_StepParticlesOut
01707             G4cout << "Capture sucess " << G4endl;
01708 #endif
01709         }
01710         if ( counter > 100 && theCollisionMgr->Entries() == 0)   // no collision, and stepping a while....
01711         {
01712 #ifdef debug_BIC_StepParticlesOut
01713             PrintKTVector(&theSecondaryList,std::string("stepping 100 steps"));
01714 #endif
01715             FindCollisions(&theSecondaryList);
01716             counter=0;
01717             ++countreset;
01718         }
01719     }
01720     //  G4cerr <<"Finished capture loop "<<G4endl;
01721 
01722     //G4cerr <<"pre- DoTimeStep 4"<<G4endl;
01723     DoTimeStep(DBL_MAX);
01724     //G4cerr <<"post- DoTimeStep 4"<<G4endl;
01725 
01726 
01727 }
01728 
01729 //----------------------------------------------------------------------------
01730 G4double G4BinaryCascade::CorrectShortlivedPrimaryForFermi(
01731         G4KineticTrack* primary,G4KineticTrackVector target_collection)
01732 //----------------------------------------------------------------------------
01733 {
01734     G4double Efermi(0);
01735     if (primary->GetState() == G4KineticTrack::inside ) {
01736         G4int PDGcode=primary->GetDefinition()->GetPDGEncoding();
01737         Efermi=((G4RKPropagation *)thePropagator)->GetField(PDGcode,primary->GetPosition());
01738 
01739         if ( std::abs(PDGcode > 1000) && PDGcode != 2112 && PDGcode != 2212 )
01740         {
01741             Efermi = ((G4RKPropagation *)thePropagator)->GetField(G4Neutron::Neutron()->GetPDGEncoding(),primary->GetPosition());
01742             G4LorentzVector mom4Primary=primary->Get4Momentum();
01743             primary->Update4Momentum(mom4Primary.e() - Efermi);
01744         }
01745 
01746         std::vector<G4KineticTrack *>::iterator titer;
01747         for ( titer=target_collection.begin() ; titer!=target_collection.end(); ++titer)
01748         {
01749             G4ParticleDefinition * aDef=(*titer)->GetDefinition();
01750             G4int aCode=aDef->GetPDGEncoding();
01751             G4ThreeVector aPos=(*titer)->GetPosition();
01752             Efermi+= ((G4RKPropagation *)thePropagator)->GetField(aCode, aPos);
01753         }
01754     }
01755     return Efermi;
01756 }
01757 
01758 //----------------------------------------------------------------------------
01759 G4bool G4BinaryCascade::CorrectShortlivedFinalsForFermi(G4KineticTrackVector * products,
01760         G4double initial_Efermi)
01761 //----------------------------------------------------------------------------
01762 {
01763     G4double final_Efermi(0);
01764     G4KineticTrackVector resonances;
01765     for ( std::vector<G4KineticTrack *>::iterator i =products->begin(); i != products->end(); i++)
01766     {
01767         G4int PDGcode=(*i)->GetDefinition()->GetPDGEncoding();
01768         //       G4cout << " PDGcode, state " << PDGcode << " " << (*i)->GetState()<<G4endl;
01769         final_Efermi+=((G4RKPropagation *)thePropagator)->GetField(PDGcode,(*i)->GetPosition());
01770         if ( std::abs(PDGcode) > 1000 && PDGcode != 2112 && PDGcode != 2212 )
01771         {
01772             resonances.push_back(*i);
01773         }
01774     }
01775     if ( resonances.size() > 0 )
01776     {
01777         G4double delta_Fermi= (initial_Efermi-final_Efermi)/resonances.size();
01778         for (std::vector<G4KineticTrack *>::iterator res=resonances.begin(); res != resonances.end(); res++)
01779         {
01780             G4LorentzVector mom=(*res)->Get4Momentum();
01781             G4double mass2=mom.mag2();
01782             G4double newEnergy=mom.e() + delta_Fermi;
01783             G4double newEnergy2= newEnergy*newEnergy;
01784             //G4cout << "mom = " << mom <<" newE " << newEnergy<< G4endl;
01785             if ( newEnergy2 < mass2 )
01786             {
01787                 return false;
01788             }
01789             //    G4cout << " correct resonance from /to " << mom.e() << " / " << newEnergy<< G4endl;
01790             G4ThreeVector mom3=std::sqrt(newEnergy2 - mass2) * mom.vect().unit();
01791             (*res)->Set4Momentum(G4LorentzVector(mom3,newEnergy));
01792         }
01793     }
01794     return true;
01795 }
01796 
01797 //----------------------------------------------------------------------------
01798 void G4BinaryCascade::CorrectFinalPandE()
01799 //----------------------------------------------------------------------------
01800 //
01801 //  Modify momenta of outgoing particles.
01802 //   Assume two body decay, nucleus(@nominal mass) + sum of final state particles(SFSP).
01803 //   momentum of SFSP shall be less than momentum for two body decay.
01804 //
01805 {
01806 #ifdef debug_BIC_CorrectFinalPandE
01807     G4cerr << "BIC: -CorrectFinalPandE called" << G4endl;
01808 #endif
01809 
01810     if ( theFinalState.size() == 0 ) return;
01811 
01812     G4KineticTrackVector::iterator i;
01813     G4LorentzVector pNucleus=GetFinal4Momentum();
01814     if ( pNucleus.e() == 0 ) return;    // check against explicit 0 from GetNucleus4Momentum()
01815 #ifdef debug_BIC_CorrectFinalPandE
01816     G4cerr << " -CorrectFinalPandE 3" << G4endl;
01817 #endif
01818     G4LorentzVector pFinals(0);
01819     G4int nFinals(0);
01820     for(i = theFinalState.begin(); i != theFinalState.end(); ++i)
01821     {
01822         pFinals += (*i)->Get4Momentum();
01823         ++nFinals;
01824 #ifdef debug_BIC_CorrectFinalPandE
01825         G4cout <<"CorrectFinalPandE a final " << (*i)->GetDefinition()->GetParticleName()
01826                                    << " 4mom " << (*i)->Get4Momentum()<< G4endl;
01827 #endif
01828     }
01829 #ifdef debug_BIC_CorrectFinalPandE
01830     G4cout << "CorrectFinalPandE pN pF: " <<pNucleus << " " <<pFinals << G4endl;
01831 #endif
01832     G4LorentzVector pCM=pNucleus + pFinals;
01833 
01834     G4LorentzRotation toCMS(-pCM.boostVector());
01835     pFinals *=toCMS;
01836 #ifdef debug_BIC_CorrectFinalPandE
01837     G4cout << "CorrectFinalPandE pCM, CMS pCM " << pCM << " " <<toCMS*pCM<< G4endl;
01838     G4cout << "CorrectFinal CMS pN pF " <<toCMS*pNucleus << " "
01839             <<pFinals << G4endl
01840             << " nucleus initial mass : " <<GetFinal4Momentum().mag()
01841             <<" massInNucleus m(nucleus) m(finals) std::sqrt(s): " << massInNucleus << " " <<pNucleus.mag()<< " "
01842             << pFinals.mag() << " " << pCM.mag() << G4endl;
01843 #endif
01844 
01845     G4LorentzRotation toLab = toCMS.inverse();
01846 
01847     G4double s0 = pCM.mag2();
01848     G4double m10 = GetIonMass(currentZ,currentA);
01849     G4double m20 = pFinals.mag();
01850     if( s0-(m10+m20)*(m10+m20) < 0 )
01851     {
01852 #ifdef debug_BIC_CorrectFinalPandE
01853         G4cout << "G4BinaryCascade::CorrectFinalPandE() : error! " << G4endl;
01854 
01855         G4cout << "not enough mass to correct: mass, A,Z, mass(nucl), mass(finals) "
01856                 << std::sqrt(s0-(m10+m20)*(m10+m20)) << " "
01857                 << currentA << " " << currentZ << " "
01858                 << m10 << " " << m20
01859                 << G4endl;
01860         G4cerr << " -CorrectFinalPandE 4" << G4endl;
01861 
01862         PrintKTVector(&theFinalState," mass problem");
01863 #endif
01864         return;
01865     }
01866 
01867     // Three momentum in cm system
01868     G4double pInCM = std::sqrt((s0-(m10+m20)*(m10+m20))*(s0-(m10-m20)*(m10-m20))/(4.*s0));
01869 #ifdef debug_BIC_CorrectFinalPandE
01870     G4cout <<" CorrectFinalPandE pInCM  new, CURRENT, ratio : " << pInCM
01871             << " " << (pFinals).vect().mag()<< " " <<  pInCM/(pFinals).vect().mag() << G4endl;
01872 #endif
01873     if ( pFinals.vect().mag() > pInCM )
01874     {
01875         G4ThreeVector p3finals=pInCM*pFinals.vect().unit();
01876 
01877         //    G4ThreeVector deltap=(p3finals - pFinals.vect() ) / nFinals;
01878         G4double factor=std::max(0.98,pInCM/pFinals.vect().mag());   // small correction
01879         G4LorentzVector qFinals(0);
01880         for(i = theFinalState.begin(); i != theFinalState.end(); ++i)
01881         {
01882             //      G4ThreeVector p3((toCMS*(*i)->Get4Momentum()).vect() + deltap);
01883             G4ThreeVector p3(factor*(toCMS*(*i)->Get4Momentum()).vect());
01884             G4LorentzVector p(p3,std::sqrt((*i)->Get4Momentum().mag2() + p3.mag2()));
01885             qFinals += p;
01886             p *= toLab;
01887 #ifdef debug_BIC_CorrectFinalPandE
01888             G4cout << " final p corrected: " << p << G4endl;
01889 #endif
01890             (*i)->Set4Momentum(p);
01891         }
01892 #ifdef debug_BIC_CorrectFinalPandE
01893         G4cout << "CorrectFinalPandE nucleus corrected mass : " << GetFinal4Momentum() << " "
01894                 <<GetFinal4Momentum().mag() << G4endl
01895                 << " CMS pFinals , mag, 3.mag : " << qFinals << " " << qFinals.mag() << " " << qFinals.vect().mag()<< G4endl;
01896         G4cerr << " -CorrectFinalPandE 5 " << factor <<  G4endl;
01897 #endif
01898     }
01899 #ifdef debug_BIC_CorrectFinalPandE
01900     else { G4cerr << " -CorrectFinalPandE 6 - no correction done" << G4endl; }
01901 #endif
01902 
01903 }
01904 
01905 //----------------------------------------------------------------------------
01906 void G4BinaryCascade::UpdateTracksAndCollisions(
01907         //----------------------------------------------------------------------------
01908         G4KineticTrackVector * oldSecondaries,
01909         G4KineticTrackVector * oldTarget,
01910         G4KineticTrackVector * newSecondaries)
01911 {
01912     std::vector<G4KineticTrack *>::iterator iter1, iter2;
01913 
01914     // remove old secondaries from the secondary list
01915     if(oldSecondaries)
01916     {
01917         if(!oldSecondaries->empty())
01918         {
01919             for(iter1 = oldSecondaries->begin(); iter1 != oldSecondaries->end();
01920                     ++iter1)
01921             {
01922                 iter2 = std::find(theSecondaryList.begin(), theSecondaryList.end(),
01923                         *iter1);
01924                 if ( iter2 != theSecondaryList.end() ) theSecondaryList.erase(iter2);
01925             }
01926             theCollisionMgr->RemoveTracksCollisions(oldSecondaries);
01927         }
01928     }
01929 
01930     // remove old target from the target list
01931     if(oldTarget)
01932     {
01933         // G4cout << "################## Debugging 0 "<<G4endl;
01934         if(oldTarget->size()!=0)
01935         {
01936 
01937             // G4cout << "################## Debugging 1 "<<oldTarget->size()<<G4endl;
01938             for(iter1 = oldTarget->begin(); iter1 != oldTarget->end(); ++iter1)
01939             {
01940                 iter2 = std::find(theTargetList.begin(), theTargetList.end(),
01941                         *iter1);
01942                 theTargetList.erase(iter2);
01943             }
01944             theCollisionMgr->RemoveTracksCollisions(oldTarget);
01945         }
01946     }
01947 
01948     if(newSecondaries)
01949     {
01950         if(!newSecondaries->empty())
01951         {
01952             // insert new secondaries in the secondary list
01953             for(iter1 = newSecondaries->begin(); iter1 != newSecondaries->end();
01954                     ++iter1)
01955             {
01956                 theSecondaryList.push_back(*iter1);
01957                 if ((*iter1)->GetState() == G4KineticTrack::undefined)
01958                 {
01959                    PrintKTVector(*iter1, "undefined in FindCollisions");
01960                 }
01961 
01962 
01963             }
01964             // look for collisions of new secondaries
01965             FindCollisions(newSecondaries);
01966         }
01967     }
01968     // G4cout << "Exiting ... "<<oldTarget<<G4endl;
01969 }
01970 
01971 
01972 class SelectFromKTV
01973 {
01974 private:
01975     G4KineticTrackVector * ktv;
01976     G4KineticTrack::CascadeState wanted_state;
01977 public:
01978     SelectFromKTV(G4KineticTrackVector * out, G4KineticTrack::CascadeState astate)
01979     :
01980         ktv(out), wanted_state(astate)
01981     {};
01982     void operator() (G4KineticTrack *& kt) const
01983     {
01984         if ( (kt)->GetState() == wanted_state ) ktv->push_back(kt);
01985     };
01986 };
01987 
01988 
01989 
01990 //----------------------------------------------------------------------------
01991 G4bool G4BinaryCascade::DoTimeStep(G4double theTimeStep)
01992 //----------------------------------------------------------------------------
01993 {
01994 
01995 #ifdef debug_BIC_DoTimeStep
01996     G4ping debug("debug_G4BinaryCascade");
01997     debug.push_back("======> DoTimeStep 1"); debug.dump();
01998     G4cerr <<"G4BinaryCascade::DoTimeStep: enter step="<< theTimeStep
01999             << " , time="<<theCurrentTime << G4endl;
02000     PrintKTVector(&theSecondaryList, std::string("DoTimeStep - theSecondaryList"));
02001     //PrintKTVector(&theTargetList, std::string("DoTimeStep - theTargetList"));
02002 #endif
02003 
02004     G4bool success=true;
02005     std::vector<G4KineticTrack *>::iterator iter;
02006 
02007     G4KineticTrackVector * kt_outside = new G4KineticTrackVector;
02008     std::for_each( theSecondaryList.begin(),theSecondaryList.end(),
02009             SelectFromKTV(kt_outside,G4KineticTrack::outside));
02010     //PrintKTVector(kt_outside, std::string("DoTimeStep - found outside"));
02011 
02012     G4KineticTrackVector * kt_inside = new G4KineticTrackVector;
02013     std::for_each( theSecondaryList.begin(),theSecondaryList.end(),
02014             SelectFromKTV(kt_inside, G4KineticTrack::inside));
02015     //  PrintKTVector(kt_inside, std::string("DoTimeStep - found inside"));
02016     //-----
02017     G4KineticTrackVector dummy;   // needed for re-usability
02018 #ifdef debug_BIC_DoTimeStep
02019     G4cout << "NOW WE ARE ENTERING THE TRANSPORT"<<G4endl;
02020 #endif
02021 
02022     // =================== Here we move the particles  ===================
02023 
02024     thePropagator->Transport(theSecondaryList, dummy, theTimeStep);
02025 
02026     // =================== Here we move the particles  ===================
02027 
02028     //------
02029 
02030     theMomentumTransfer += thePropagator->GetMomentumTransfer();
02031 #ifdef debug_BIC_DoTimeStep
02032     G4cout << "DoTimeStep : theMomentumTransfer = " << theMomentumTransfer << G4endl;
02033     PrintKTVector(&theSecondaryList, std::string("DoTimeStep - secondaries aft trsprt"));
02034 #endif
02035 
02036     // Partclies which went INTO nucleus
02037 
02038     G4KineticTrackVector * kt_gone_in = new G4KineticTrackVector;
02039     std::for_each( kt_outside->begin(),kt_outside->end(),
02040             SelectFromKTV(kt_gone_in,G4KineticTrack::inside));
02041     //  PrintKTVector(kt_gone_in, std::string("DoTimeStep - gone in"));
02042 
02043 
02044     // Partclies which  went OUT OF nucleus
02045     G4KineticTrackVector * kt_gone_out = new G4KineticTrackVector;
02046     std::for_each( kt_inside->begin(),kt_inside->end(),
02047             SelectFromKTV(kt_gone_out, G4KineticTrack::gone_out));
02048 
02049     //  PrintKTVector(kt_gone_out, std::string("DoTimeStep - gone out"));
02050 
02051     G4KineticTrackVector *fail=CorrectBarionsOnBoundary(kt_gone_in,kt_gone_out);
02052 
02053     if ( fail )
02054     {
02055         // some particle(s) supposed to enter/leave were miss_nucleus/captured by the correction
02056         kt_gone_in->clear();
02057         std::for_each( kt_outside->begin(),kt_outside->end(),
02058                 SelectFromKTV(kt_gone_in,G4KineticTrack::inside));
02059 
02060         kt_gone_out->clear();
02061         std::for_each( kt_inside->begin(),kt_inside->end(),
02062                 SelectFromKTV(kt_gone_out, G4KineticTrack::gone_out));
02063 
02064 #ifdef debug_BIC_DoTimeStep
02065         PrintKTVector(fail,std::string(" Failed to go in/out -> miss_nucleus/captured"));
02066         PrintKTVector(kt_gone_in, std::string("recreated kt_gone_in"));
02067         PrintKTVector(kt_gone_out, std::string("recreated kt_gone_out"));
02068 #endif
02069         delete fail;
02070     }
02071 
02072     // Add tracks missing nucleus and tracks going straight though  to addFinals
02073     std::for_each( kt_outside->begin(),kt_outside->end(),
02074             SelectFromKTV(kt_gone_out,G4KineticTrack::miss_nucleus));
02075     //PrintKTVector(kt_gone_out, std::string("miss to append to final state.."));
02076     std::for_each( kt_outside->begin(),kt_outside->end(),
02077             SelectFromKTV(kt_gone_out,G4KineticTrack::gone_out));
02078 
02079 #ifdef debug_BIC_DoTimeStep
02080     PrintKTVector(kt_gone_out, std::string("append gone_outs to final state.. theFinalState"));
02081 #endif
02082 
02083     theFinalState.insert(theFinalState.end(),
02084             kt_gone_out->begin(),kt_gone_out->end());
02085 
02086     // Partclies which could not leave nucleus,  captured...
02087     G4KineticTrackVector * kt_captured = new G4KineticTrackVector;
02088     std::for_each( theSecondaryList.begin(),theSecondaryList.end(),
02089             SelectFromKTV(kt_captured, G4KineticTrack::captured));
02090 
02091     // Check no track is part in next collision, ie.
02092     //  this step was to far, and collisions should not occur any more
02093 
02094     if ( theCollisionMgr->Entries()> 0 )
02095     {
02096         if (kt_gone_out->size() )
02097         {
02098             G4KineticTrack * nextPrimary = theCollisionMgr->GetNextCollision()->GetPrimary();
02099             iter = std::find(kt_gone_out->begin(),kt_gone_out->end(),nextPrimary);
02100             if ( iter !=  kt_gone_out->end() )
02101             {
02102                 success=false;
02103 #ifdef debug_BIC_DoTimeStep
02104                 G4cout << " DoTimeStep - WARNING: deleting current collision!" << G4endl;
02105 #endif
02106             }
02107         }
02108         if ( kt_captured->size() )
02109         {
02110             G4KineticTrack * nextPrimary = theCollisionMgr->GetNextCollision()->GetPrimary();
02111             iter = std::find(kt_captured->begin(),kt_captured->end(),nextPrimary);
02112             if ( iter !=  kt_captured->end() )
02113             {
02114                 success=false;
02115 #ifdef debug_BIC_DoTimeStep
02116                 G4cout << " DoTimeStep - WARNING: deleting current collision!" << G4endl;
02117 #endif
02118             }
02119         }
02120 
02121     }
02122     // PrintKTVector(kt_gone_out," kt_gone_out be4 updatetrack...");
02123     UpdateTracksAndCollisions(kt_gone_out,0 ,0);
02124 
02125 
02126     if ( kt_captured->size() )
02127     {
02128         theCapturedList.insert(theCapturedList.end(),
02129                 kt_captured->begin(),kt_captured->end());
02130         //should be      std::for_each(kt_captured->begin(),kt_captured->end(),
02131         //              std::mem_fun(&G4KineticTrack::Hit));
02132         // but VC 6 requires:
02133         std::vector<G4KineticTrack *>::iterator i_captured;
02134         for(i_captured=kt_captured->begin();i_captured!=kt_captured->end();i_captured++)
02135         {
02136             (*i_captured)->Hit();
02137         }
02138         //     PrintKTVector(kt_captured," kt_captured be4 updatetrack...");
02139         UpdateTracksAndCollisions(kt_captured, NULL, NULL);
02140     }
02141 
02142 #ifdef debug_G4BinaryCascade
02143     delete kt_inside;
02144     kt_inside = new G4KineticTrackVector;
02145     std::for_each( theSecondaryList.begin(),theSecondaryList.end(),
02146             SelectFromKTV(kt_inside, G4KineticTrack::inside));
02147     if ( currentZ != (GetTotalCharge(theTargetList)
02148             + GetTotalCharge(theCapturedList)
02149             + GetTotalCharge(*kt_inside)) )
02150     {
02151         G4cout << " error-DoTimeStep aft, A, Z: " << currentA << " " << currentZ
02152                 << " sum(tgt,capt,active) "
02153                 << GetTotalCharge(theTargetList) + GetTotalCharge(theCapturedList) + GetTotalCharge(*kt_inside)
02154                 << " targets: "  << GetTotalCharge(theTargetList)
02155                 << " captured: " << GetTotalCharge(theCapturedList)
02156                 << " active: "   << GetTotalCharge(*kt_inside)
02157                 << G4endl;
02158     }
02159 #endif
02160 
02161     delete kt_inside;
02162     delete kt_outside;
02163     delete kt_captured;
02164     delete kt_gone_in;
02165     delete kt_gone_out;
02166 
02167     //  G4cerr <<"G4BinaryCascade::DoTimeStep: exit "<<G4endl;
02168     theCurrentTime += theTimeStep;
02169 
02170     //debug.push_back("======> DoTimeStep 2"); debug.dump();
02171     return success;
02172 
02173 }
02174 
02175 //----------------------------------------------------------------------------
02176 G4KineticTrackVector* G4BinaryCascade::CorrectBarionsOnBoundary(
02177         G4KineticTrackVector *in,
02178         G4KineticTrackVector *out)
02179 //----------------------------------------------------------------------------
02180 {
02181     G4KineticTrackVector * kt_fail(0);
02182     std::vector<G4KineticTrack *>::iterator iter;
02183     //  G4cout << "CorrectBarionsOnBoundary,currentZ,currentA,"
02184     //         << currentZ << " "<< currentA << G4endl;
02185     if (in->size())
02186     {
02187         G4int secondaries_in(0);
02188         G4int secondaryBarions_in(0);
02189         G4int secondaryCharge_in(0);
02190         G4double secondaryMass_in(0);
02191 
02192         for ( iter =in->begin(); iter != in->end(); ++iter)
02193         {
02194             ++secondaries_in;
02195             secondaryCharge_in += G4lrint((*iter)->GetDefinition()->GetPDGCharge()/eplus);
02196             if ((*iter)->GetDefinition()->GetBaryonNumber()!=0 )
02197             {
02198                 secondaryBarions_in += (*iter)->GetDefinition()->GetBaryonNumber();
02199                 if((*iter)->GetDefinition() == G4Neutron::Neutron() ||
02200                         (*iter)->GetDefinition() == G4Proton::Proton()  )
02201                 {
02202                     secondaryMass_in += (*iter)->GetDefinition()->GetPDGMass();
02203                 } else    {
02204                     secondaryMass_in += G4Proton::Proton()->GetPDGMass();
02205                 }
02206             }
02207         }
02208         G4double mass_initial= GetIonMass(currentZ,currentA);
02209 
02210         currentZ += secondaryCharge_in;
02211         currentA += secondaryBarions_in;
02212 
02213         //  G4cout << "CorrectBarionsOnBoundary,secondaryCharge_in, secondaryBarions_in "
02214         //         <<    secondaryCharge_in << " "<<  secondaryBarions_in << G4endl;
02215 
02216         G4double mass_final= GetIonMass(currentZ,currentA);
02217 
02218         G4double correction= secondaryMass_in + mass_initial - mass_final;
02219         if (secondaries_in>1)
02220         {correction /= secondaries_in;}
02221 
02222 #ifdef debug_BIC_CorrectBarionsOnBoundary
02223         G4cout << "CorrectBarionsOnBoundary,currentZ,currentA,"
02224                 << "secondaryCharge_in,secondaryBarions_in,"
02225                 << "energy correction,m_secondry,m_nucl_init,m_nucl_final "
02226                 << currentZ << " "<< currentA <<" "
02227                 << secondaryCharge_in<<" "<<secondaryBarions_in<<" "
02228                 << correction << " "
02229                 << secondaryMass_in << " "
02230                 << mass_initial << " "
02231                 << mass_final << " "
02232                 << G4endl;
02233         PrintKTVector(in,std::string("in be4 correction"));
02234 #endif
02235 
02236         for ( iter = in->begin(); iter != in->end(); ++iter)
02237         {
02238             if ((*iter)->GetTrackingMomentum().e()+correction > (*iter)->GetActualMass())
02239             {
02240                 (*iter)->UpdateTrackingMomentum((*iter)->GetTrackingMomentum().e() + correction);
02241             } else {
02242                 //particle cannot go in, put to miss_nucleus
02243                 G4RKPropagation * RKprop=(G4RKPropagation *)thePropagator;
02244                 (*iter)->SetState(G4KineticTrack::miss_nucleus);
02245                 // Undo correction for Colomb Barrier
02246                 G4double barrier=RKprop->GetBarrier((*iter)->GetDefinition()->GetPDGEncoding());
02247                 (*iter)->UpdateTrackingMomentum((*iter)->GetTrackingMomentum().e() + barrier);
02248                 if ( ! kt_fail ) kt_fail=new G4KineticTrackVector;
02249                 kt_fail->push_back(*iter);
02250                 currentZ -= G4lrint((*iter)->GetDefinition()->GetPDGCharge()/eplus);
02251                 currentA -= (*iter)->GetDefinition()->GetBaryonNumber();
02252 
02253             }
02254 
02255         }
02256 #ifdef debug_BIC_CorrectBarionsOnBoundary
02257         G4cout << " CorrectBarionsOnBoundary, aft, A, Z, sec-Z,A,m,m_in_nucleus "
02258                 << currentA << " " << currentZ << " "
02259                 << secondaryCharge_in << " " << secondaryBarions_in << " "
02260                 << secondaryMass_in  << " "
02261                 << G4endl;
02262         PrintKTVector(in,std::string("in AFT correction"));
02263 #endif
02264 
02265     }
02266     //----------------------------------------------
02267     if (out->size())
02268     {
02269         G4int secondaries_out(0);
02270         G4int secondaryBarions_out(0);
02271         G4int secondaryCharge_out(0);
02272         G4double secondaryMass_out(0);
02273 
02274         for ( iter =out->begin(); iter != out->end(); ++iter)
02275         {
02276             ++secondaries_out;
02277             secondaryCharge_out += G4lrint((*iter)->GetDefinition()->GetPDGCharge()/eplus);
02278             if ((*iter)->GetDefinition()->GetBaryonNumber() !=0 )
02279             {
02280                 secondaryBarions_out += (*iter)->GetDefinition()->GetBaryonNumber();
02281                 if((*iter)->GetDefinition() == G4Neutron::Neutron() ||
02282                         (*iter)->GetDefinition() == G4Proton::Proton()  )
02283                 {
02284                     secondaryMass_out += (*iter)->GetDefinition()->GetPDGMass();
02285                 } else {
02286                     secondaryMass_out += G4Neutron::Neutron()->GetPDGMass();
02287                 }
02288             }
02289         }
02290 
02291         G4double mass_initial=  GetIonMass(currentZ,currentA);
02292         currentA -=secondaryBarions_out;
02293         currentZ -=secondaryCharge_out;
02294 
02295         //  G4cout << "CorrectBarionsOnBoundary,secondaryCharge_out, secondaryBarions_out"
02296         //         <<    secondaryCharge_out << " "<<  secondaryBarions_out << G4endl;
02297 
02298         //                        a delta minus will do currentZ < 0 in light nuclei
02299         //     if (currentA < 0 || currentZ < 0 )
02300         if (currentA < 0 )
02301         {
02302             G4cerr << "G4BinaryCascade - secondaryBarions_out,secondaryCharge_out " <<
02303                     secondaryBarions_out << " " << secondaryCharge_out << G4endl;
02304             PrintKTVector(&theTargetList,"CorrectBarionsOnBoundary Target");
02305             PrintKTVector(&theCapturedList,"CorrectBarionsOnBoundary Captured");
02306             PrintKTVector(&theSecondaryList,"CorrectBarionsOnBoundary Secondaries");
02307             G4cerr << "G4BinaryCascade - currentA, currentZ " << currentA << " " << currentZ << G4endl;
02308             throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade::CorrectBarionsOnBoundary() - fatal error");
02309         }
02310         G4double mass_final=GetIonMass(currentZ,currentA);
02311         G4double correction= mass_initial - mass_final - secondaryMass_out;
02312 
02313         if (secondaries_out>1) correction /= secondaries_out;
02314 #ifdef debug_BIC_CorrectBarionsOnBoundary
02315         G4cout << "DoTimeStep,currentZ,currentA,"
02316                 << "secondaries_out,"
02317                 <<"secondaryCharge_out,secondaryBarions_out,"
02318                 <<"energy correction,m_secondry,m_nucl_init,m_nucl_final "
02319                 << " "<< currentZ << " "<< currentA <<" "
02320                 << secondaries_out << " "
02321                 << secondaryCharge_out<<" "<<secondaryBarions_out<<" "
02322                 << correction << " "
02323                 << secondaryMass_out << " "
02324                 << mass_initial << " "
02325                 << mass_final << " "
02326                 << G4endl;
02327         PrintKTVector(out,std::string("out be4 correction"));
02328 #endif
02329 
02330         for ( iter = out->begin(); iter != out->end(); ++iter)
02331         {
02332             if ((*iter)->GetTrackingMomentum().e()+correction > (*iter)->GetActualMass())
02333             {
02334                 (*iter)->UpdateTrackingMomentum((*iter)->GetTrackingMomentum().e() + correction);
02335             } else
02336             {
02337                 // particle cannot go out due to change of nuclear potential!
02338                 //  capture protons and neutrons;
02339                 if(((*iter)->GetDefinition() == G4Proton::Proton()) ||
02340                         ((*iter)->GetDefinition() == G4Neutron::Neutron()))
02341                 {
02342                     (*iter)->SetState(G4KineticTrack::captured);
02343                     // Undo correction for Colomb Barrier
02344                     G4double barrier=((G4RKPropagation *)thePropagator)->GetBarrier((*iter)->GetDefinition()->GetPDGEncoding());
02345                     (*iter)->UpdateTrackingMomentum((*iter)->GetTrackingMomentum().e() - barrier);
02346                     if ( kt_fail == 0 ) kt_fail=new G4KineticTrackVector;
02347                     kt_fail->push_back(*iter);
02348                     currentZ += G4lrint((*iter)->GetDefinition()->GetPDGCharge()/eplus);
02349                     currentA += (*iter)->GetDefinition()->GetBaryonNumber();
02350                 }
02351 #ifdef debug_BIC_CorrectBarionsOnBoundary
02352                 else
02353                 {
02354                     G4cout << "Not correcting outgoing " << *iter << " "
02355                             << (*iter)->GetDefinition()->GetPDGEncoding() << " "
02356                             << (*iter)->GetDefinition()->GetParticleName() << G4endl;
02357                     PrintKTVector(out,std::string("outgoing, one not corrected"));
02358                 }
02359 #endif
02360             }
02361         }
02362 
02363 #ifdef debug_BIC_CorrectBarionsOnBoundary
02364         PrintKTVector(out,std::string("out AFTER correction"));
02365         G4cout << " DoTimeStep, nucl-update, A, Z, sec-Z,A,m,m_in_nucleus, table-mass, delta "
02366                 << currentA << " "<< currentZ << " "
02367                 << secondaryCharge_out << " "<< secondaryBarions_out << " "<<
02368                 secondaryMass_out << " "
02369                 << massInNucleus << " "
02370                 << G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(currentZ,currentA)
02371                 << " " << massInNucleus -G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(currentZ,currentA)
02372                 << G4endl;
02373 #endif
02374     }
02375 
02376     return kt_fail;
02377 }
02378 
02379 
02380 //----------------------------------------------------------------------------
02381 
02382 G4Fragment * G4BinaryCascade::FindFragments()
02383 //----------------------------------------------------------------------------
02384 {
02385 
02386 #ifdef debug_BIC_FindFragments
02387     G4cout << "target, captured, secondary: "
02388             << theTargetList.size() << " "
02389             << theCapturedList.size()<< " "
02390             << theSecondaryList.size()
02391             << G4endl;
02392 #endif
02393 
02394     G4int a = theTargetList.size()+theCapturedList.size();
02395     G4int zTarget = 0;
02396     G4KineticTrackVector::iterator i;
02397     for(i = theTargetList.begin(); i != theTargetList.end(); ++i)
02398     {
02399         if(G4lrint((*i)->GetDefinition()->GetPDGCharge()/eplus) == 1 )
02400         {
02401             zTarget++;
02402         }
02403     }
02404 
02405     G4int zCaptured = 0;
02406     G4LorentzVector CapturedMomentum(0.,0.,0.,0.);
02407     for(i = theCapturedList.begin(); i != theCapturedList.end(); ++i)
02408     {
02409         CapturedMomentum += (*i)->Get4Momentum();
02410         if(G4lrint((*i)->GetDefinition()->GetPDGCharge()/eplus) == 1 )
02411         {
02412             zCaptured++;
02413         }
02414     }
02415 
02416     G4int z = zTarget+zCaptured;
02417 
02418 #ifdef debug_G4BinaryCascade
02419     if ( z != (GetTotalCharge(theTargetList) + GetTotalCharge(theCapturedList)) )
02420     {
02421         G4cout << " FindFragment Counting error z a " << z << " " <<a << " "
02422                 << GetTotalCharge(theTargetList) << " " <<  GetTotalCharge(theCapturedList)<<
02423                 G4endl;
02424         PrintKTVector(&theTargetList, std::string("theTargetList"));
02425         PrintKTVector(&theCapturedList, std::string("theCapturedList"));
02426     }
02427 #endif
02428     //debug
02429     /*
02430      *   G4cout << " Fragment mass table / real "
02431      *          << G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(z, a)
02432      *   << " / " << GetFinal4Momentum().mag()
02433      *   << " difference "
02434      *   <<  GetFinal4Momentum().mag() - G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(z, a)
02435      *   << G4endl;
02436      */
02437     //
02438     //  if(getenv("BCDEBUG") ) G4cerr << "Fragment A, Z "<< a <<" "<< z<<G4endl;
02439     if ( z < 1 ) return 0;
02440 
02441     G4int holes = the3DNucleus->GetMassNumber() - theTargetList.size();
02442     G4int excitons = theCapturedList.size();
02443 #ifdef debug_BIC_FindFragments
02444     G4cout << "Fragment: a= " << a << " z= " << z << " particles= " <<  excitons
02445             << " Charged= " << zCaptured << " holes= " << holes
02446             << " excitE= " <<GetExcitationEnergy()
02447             << " Final4Momentum= " << GetFinalNucleusMomentum() << " capturMomentum= " << CapturedMomentum
02448             << G4endl;
02449 #endif
02450 
02451     G4Fragment * fragment = new G4Fragment(a,z,GetFinalNucleusMomentum());
02452     fragment->SetNumberOfHoles(holes);
02453 
02454     //GF  fragment->SetNumberOfParticles(excitons-holes);
02455     fragment->SetNumberOfParticles(excitons);
02456     fragment->SetNumberOfCharged(zCaptured);
02457     G4ParticleDefinition * aIonDefinition =
02458             G4ParticleTable::GetParticleTable()->FindIon(a,z,0,z);
02459     fragment->SetParticleDefinition(aIonDefinition);
02460 
02461     return fragment;
02462 }
02463 
02464 //----------------------------------------------------------------------------
02465 
02466 G4LorentzVector G4BinaryCascade::GetFinal4Momentum()
02467 //----------------------------------------------------------------------------
02468 // Return momentum of reminder nulceus;
02469 //  ie. difference of (initial state(primaries+nucleus) - final state) particles, ignoring remnant nucleus
02470 {
02471     G4LorentzVector final4Momentum = theInitial4Mom + theProjectile4Momentum;
02472     G4LorentzVector finals(0,0,0,0);
02473     for(G4KineticTrackVector::iterator i = theFinalState.begin(); i != theFinalState.end(); ++i)
02474     {
02475         final4Momentum -= (*i)->Get4Momentum();
02476         finals             += (*i)->Get4Momentum();
02477     }
02478 
02479     if((final4Momentum.vect()/final4Momentum.e()).mag()>1.0 && currentA > 0)
02480     {
02481 #ifdef debug_BIC_Final4Momentum
02482         G4cerr << G4endl;
02483         G4cerr << "G4BinaryCascade::GetFinal4Momentum - Fatal"<<G4endl;
02484         G4KineticTrackVector::iterator i;
02485         G4cerr <<"Total initial 4-momentum " << theProjectile4Momentum << G4endl;
02486         G4cerr <<" GetFinal4Momentum: Initial nucleus "<<theInitial4Mom<<G4endl;
02487         for(i = theFinalState.begin(); i != theFinalState.end(); ++i)
02488         {
02489             G4cerr <<" Final state: "<<(*i)->Get4Momentum()<<(*i)->GetDefinition()->GetParticleName()<<G4endl;
02490         }
02491         G4cerr << "Sum( 4-mon ) finals " << finals << G4endl;
02492         G4cerr<< " Final4Momentum = "<<final4Momentum <<" "<<final4Momentum.m()<<G4endl;
02493         G4cerr <<" current A, Z = "<< currentA<<", "<<currentZ<<G4endl;
02494         G4cerr << G4endl;
02495 #endif
02496 
02497         final4Momentum=G4LorentzVector(0,0,0,0);
02498     }
02499     return final4Momentum;
02500 }
02501 
02502 //----------------------------------------------------------------------------
02503 G4LorentzVector G4BinaryCascade::GetFinalNucleusMomentum()
02504 //----------------------------------------------------------------------------
02505 {
02506     // return momentum of nucleus for use with precompound model; also keep transformation to
02507     //   apply to precompoud products.
02508 
02509     G4LorentzVector CapturedMomentum(0,0,0,0);
02510     G4KineticTrackVector::iterator i;
02511     //  G4cout << "GetFinalNucleusMomentum Captured size: " <<theCapturedList.size() << G4endl;
02512     for(i = theCapturedList.begin(); i != theCapturedList.end(); ++i)
02513     {
02514         CapturedMomentum += (*i)->Get4Momentum();
02515     }
02516     //G4cout << "GetFinalNucleusMomentum CapturedMomentum= " <<CapturedMomentum << G4endl;
02517     //  G4cerr << "it 9"<<G4endl;
02518 
02519     G4LorentzVector NucleusMomentum = GetFinal4Momentum();
02520     if ( NucleusMomentum.e() > 0 )
02521     {
02522         // G4cout << "GetFinalNucleusMomentum GetFinal4Momentum= " <<NucleusMomentum <<" "<<NucleusMomentum.mag()<<G4endl;
02523         // boost nucleus to a frame such that the momentum of nucleus == momentum of Captured
02524         G4ThreeVector boost= (NucleusMomentum.vect() -CapturedMomentum.vect())/NucleusMomentum.e();
02525         if(boost.mag2()>1.0)
02526         {
02527 #     ifdef debug_BIC_FinalNucleusMomentum
02528             G4cerr << "G4BinaryCascade::GetFinalNucleusMomentum - Fatal"<<G4endl;
02529             G4cerr << "it 0"<<boost <<G4endl;
02530             G4cerr << "it 01"<<NucleusMomentum<<" "<<CapturedMomentum<<" "<<G4endl;
02531             G4cout <<" testing boost "<<boost<<" "<<boost.mag()<<G4endl;
02532 #      endif
02533             boost=G4ThreeVector(0,0,0);
02534             NucleusMomentum=G4LorentzVector(0,0,0,0);
02535         }
02536         G4LorentzRotation  nucleusBoost( -boost );
02537         precompoundLorentzboost.set( boost );
02538 #ifdef debug_debug_BIC_FinalNucleusMomentum
02539         G4cout << "GetFinalNucleusMomentum be4 boostNucleusMomentum, CapturedMomentum"<<NucleusMomentum<<" "<<CapturedMomentum<<" "<<G4endl;
02540 #endif
02541         NucleusMomentum *= nucleusBoost;
02542 #ifdef debug_BIC_FinalNucleusMomentum
02543         G4cout << "GetFinalNucleusMomentum aft boost GetFinal4Momentum= " <<NucleusMomentum <<G4endl;
02544 #endif
02545     }
02546     return NucleusMomentum;
02547 }
02548 
02549 //----------------------------------------------------------------------------
02550 G4ReactionProductVector * G4BinaryCascade::Propagate1H1(
02551         //----------------------------------------------------------------------------
02552         G4KineticTrackVector * secondaries, G4V3DNucleus * nucleus)
02553 {
02554     G4ReactionProductVector * products = new G4ReactionProductVector;
02555     G4ParticleDefinition * aHTarg = G4Proton::ProtonDefinition();
02556     G4double mass = aHTarg->GetPDGMass();
02557     if (nucleus->GetCharge() == 0) aHTarg = G4Neutron::NeutronDefinition();
02558     mass = aHTarg->GetPDGMass();
02559     G4KineticTrackVector * secs = 0;
02560     G4ThreeVector pos(0,0,0);
02561     G4LorentzVector mom(mass);
02562     G4KineticTrack aTarget(aHTarg, 0., pos, mom);
02563     G4bool done(false);
02564     std::vector<G4KineticTrack *>::iterator iter, jter;
02565     // data member    static G4Scatterer theH1Scatterer;
02566     //G4cout << " start 1H1 for " << (*secondaries).front()->GetDefinition()->GetParticleName()
02567     //       << " on " << aHTarg->GetParticleName() << G4endl;
02568     G4int tryCount(0);
02569     while(!done && tryCount++ <200)
02570     {
02571         if(secs)
02572         {
02573             std::for_each(secs->begin(), secs->end(), DeleteKineticTrack());
02574             delete secs;
02575         }
02576         secs = theH1Scatterer->Scatter(*(*secondaries).front(), aTarget);
02577 #ifdef debug_H1_BinaryCascade
02578         PrintKTVector(secs," From Scatter");
02579 #endif
02580         for(size_t ss=0; secs && ss<secs->size(); ss++)
02581         {
02582             // must have one resonance in final state, or it was elastic, not allowed here.
02583             if((*secs)[ss]->GetDefinition()->IsShortLived()) done = true;
02584         }
02585     }
02586 
02587     size_t current(0);
02588     ClearAndDestroy(&theFinalState);
02589     for(current=0; secs && current<secs->size(); current++)
02590     {
02591         if((*secs)[current]->GetDefinition()->IsShortLived())
02592         {
02593             done = true;        // must have one resonance in final state, elastic not allowed here!
02594             G4KineticTrackVector * dec = (*secs)[current]->Decay();
02595             for(jter=dec->begin(); jter != dec->end(); jter++)
02596             {
02597                 //G4cout << "Decay"<<G4endl;
02598                 secs->push_back(*jter);
02599                 //G4cout << "decay "<<(*jter)->GetDefinition()->GetParticleName()<<G4endl;
02600             }
02601             delete (*secs)[current];
02602             delete dec;
02603         }
02604         else
02605         {
02606             //G4cout << "FS "<<G4endl;
02607             //G4cout << "FS "<<(*secs)[current]->GetDefinition()->GetParticleName()<<G4endl;
02608             theFinalState.push_back((*secs)[current]);
02609         }
02610     }
02611 
02612     delete secs;
02613 #ifdef debug_H1_BinaryCascade
02614     PrintKTVector(&theFinalState," FinalState");
02615 #endif
02616     for(iter = theFinalState.begin(); iter != theFinalState.end(); ++iter)
02617     {
02618         G4KineticTrack * kt = *iter;
02619         G4ReactionProduct * aNew = new G4ReactionProduct(kt->GetDefinition());
02620         aNew->SetMomentum(kt->Get4Momentum().vect());
02621         aNew->SetTotalEnergy(kt->Get4Momentum().e());
02622         products->push_back(aNew);
02623 #ifdef debug_H1_BinaryCascade
02624         if (! kt->GetDefinition()->GetPDGStable() )
02625         {
02626             if (kt->GetDefinition()->IsShortLived())
02627             {
02628                 G4cout << "final shortlived : ";
02629             } else
02630             {
02631                 G4cout << "final un stable : ";
02632             }
02633             G4cout <<kt->GetDefinition()->GetParticleName()<< G4endl;
02634         }
02635 #endif
02636         delete kt;
02637     }
02638     theFinalState.clear();
02639     return products;
02640 
02641 }
02642 
02643 //----------------------------------------------------------------------------
02644 G4ThreeVector G4BinaryCascade::GetSpherePoint(
02645         G4double r, const G4LorentzVector & mom4)
02646 //----------------------------------------------------------------------------
02647 {
02648     //  Get a point outside radius.
02649     //     point is random in plane (circle of radius r) orthogonal to mom,
02650     //      plus -1*r*mom->vect()->unit();
02651     G4ThreeVector o1, o2;
02652     G4ThreeVector mom = mom4.vect();
02653 
02654     o1= mom.orthogonal();  // we simply need any vector non parallel
02655     o2= mom.cross(o1);     //  o2 is now orthogonal to mom and o1, ie. o1 and o2  define plane.
02656 
02657     G4double x2, x1;
02658 
02659     do
02660     {
02661         x1=(G4UniformRand()-.5)*2;
02662         x2=(G4UniformRand()-.5)*2;
02663     } while (sqr(x1) +sqr(x2) > 1.);
02664 
02665     return G4ThreeVector(r*(x1*o1.unit() + x2*o2.unit() - 1.5* mom.unit()));
02666 
02667 
02668 
02669     /*
02670      * // Get a point uniformly distributed on the surface of a sphere,
02671      * // with z < 0.
02672      *   G4double b = r*G4UniformRand();  // impact parameter
02673      *   G4double phi = G4UniformRand()*2*pi;
02674      *   G4double x = b*std::cos(phi);
02675      *   G4double y = b*std::sin(phi);
02676      *   G4double z = -std::sqrt(r*r-b*b);
02677      *   z *= 1.001; // Get position a little bit out of the sphere...
02678      *   point.setX(x);
02679      *   point.setY(y);
02680      *   point.setZ(z);
02681      */
02682 }
02683 
02684 //----------------------------------------------------------------------------
02685 
02686 void G4BinaryCascade::ClearAndDestroy(G4KineticTrackVector * ktv)
02687 //----------------------------------------------------------------------------
02688 {
02689     std::vector<G4KineticTrack *>::iterator i;
02690     for(i = ktv->begin(); i != ktv->end(); ++i)
02691         delete (*i);
02692     ktv->clear();
02693 }
02694 
02695 //----------------------------------------------------------------------------
02696 void G4BinaryCascade::ClearAndDestroy(G4ReactionProductVector * rpv)
02697 //----------------------------------------------------------------------------
02698 {
02699     std::vector<G4ReactionProduct *>::iterator i;
02700     for(i = rpv->begin(); i != rpv->end(); ++i)
02701         delete (*i);
02702     rpv->clear();
02703 }
02704 
02705 //----------------------------------------------------------------------------
02706 void G4BinaryCascade::PrintKTVector(G4KineticTrackVector * ktv, std::string comment)
02707 //----------------------------------------------------------------------------
02708 {
02709     if (comment.size() > 0 ) G4cout << "G4BinaryCascade::PrintKTVector() " << comment << G4endl;
02710     if (ktv) {
02711         G4cout << "  vector: " << ktv << ", number of tracks: " << ktv->size()
02712                                  << G4endl;
02713         std::vector<G4KineticTrack *>::iterator i;
02714         G4int count;
02715 
02716         for(count = 0, i = ktv->begin(); i != ktv->end(); ++i, ++count)
02717         {
02718             G4KineticTrack * kt = *i;
02719             G4cout << "  track n. " << count;
02720             PrintKTVector(kt);
02721         }
02722     } else {
02723         G4cout << "G4BinaryCascade::PrintKTVector():No KineticTrackVector given " << G4endl;
02724     }
02725 }
02726 //----------------------------------------------------------------------------
02727 void G4BinaryCascade::PrintKTVector(G4KineticTrack * kt, std::string comment)
02728 //----------------------------------------------------------------------------
02729 {
02730     if (comment.size() > 0 ) G4cout << "G4BinaryCascade::PrintKTVector() "<< comment << G4endl;
02731     if ( kt ){
02732         G4cout << ", id: " << kt << G4endl;
02733         G4ThreeVector pos = kt->GetPosition();
02734         G4LorentzVector mom = kt->Get4Momentum();
02735         G4LorentzVector tmom = kt->GetTrackingMomentum();
02736         G4ParticleDefinition * definition = kt->GetDefinition();
02737         G4cout << "    definition: " << definition->GetPDGEncoding() << " pos: "
02738                 << 1/fermi*pos << " R: " << 1/fermi*pos.mag() << " 4mom: "
02739                 << 1/MeV*mom <<"Tr_mom" <<  1/MeV*tmom << " P: " << 1/MeV*mom.vect().mag()
02740                 << " M: " << 1/MeV*mom.mag() << G4endl;
02741         G4cout <<"    trackstatus: "<<kt->GetState()<<G4endl;
02742     } else {
02743         G4cout << "G4BinaryCascade::PrintKTVector(): No Kinetictrack given" << G4endl;
02744     }
02745 }
02746 
02747 
02748 //----------------------------------------------------------------------------
02749 G4double G4BinaryCascade::GetIonMass(G4int Z, G4int A)
02750 //----------------------------------------------------------------------------
02751 {
02752     G4double mass(0);
02753     if ( Z > 0 && A >= Z )
02754     {
02755         mass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z,A);
02756 
02757     } else if ( A > 0 && Z>0 )
02758     {
02759         // charge Z > A; will happen for light nuclei with pions involved.
02760         mass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(A,A);
02761 
02762     } else if ( A >= 0 && Z<=0 )
02763     {
02764         // all neutral, or empty nucleus
02765         mass = A * G4Neutron::Neutron()->GetPDGMass();
02766 
02767     } else if ( A == 0 && std::abs(Z)<2 )
02768     {
02769         // empty nucleus, except maybe pions
02770         mass = 0;
02771 
02772     } else
02773     {
02774         G4cerr << "G4BinaryCascade::GetIonMass() - invalid (A,Z) = ("
02775                 << A << "," << Z << ")" <<G4endl;
02776         throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade::GetIonMass() - giving up");
02777 
02778     }
02779     return mass;
02780 }
02781 G4ReactionProductVector * G4BinaryCascade::FillVoidNucleusProducts(G4ReactionProductVector * products)
02782 {
02783     // return product when nucleus is destroyed, i.e. charge=0
02784     G4double Esecondaries(0.);
02785     std::vector<G4KineticTrack *>::iterator iter;
02786     decayKTV.Decay(&theFinalState);
02787     //PrintKTVector(&theFinalState, " theFinalState @ void Nucl");
02788     for(iter = theFinalState.begin(); iter != theFinalState.end(); ++iter)
02789     {
02790         G4ReactionProduct * aNew = new G4ReactionProduct((*iter)->GetDefinition());
02791         aNew->SetMomentum((*iter)->Get4Momentum().vect());
02792         aNew->SetTotalEnergy((*iter)->Get4Momentum().e());
02793         Esecondaries +=(*iter)->Get4Momentum().e();
02794         aNew->SetNewlyAdded(true);
02795         //G4cout << " Particle Ekin " << aNew->GetKineticEnergy() << G4endl;
02796         products->push_back(aNew);
02797     }
02798 
02799     //PrintKTVector(&theCapturedList, " theCapturedList @ void Nucl");
02800     for(iter = theCapturedList.begin(); iter != theCapturedList.end(); ++iter)
02801     {
02802         G4ReactionProduct * aNew = new G4ReactionProduct((*iter)->GetDefinition());
02803         aNew->SetMomentum((*iter)->Get4Momentum().vect());
02804         aNew->SetTotalEnergy((*iter)->Get4Momentum().e());
02805         Esecondaries +=(*iter)->Get4Momentum().e();
02806         aNew->SetNewlyAdded(true);
02807         //G4cout << " Particle Ekin " << aNew->GetKineticEnergy() << G4endl;
02808         products->push_back(aNew);
02809     }
02810     // pull out late particles from collisions
02811     //theCollisionMgr->Print();
02812     while(theCollisionMgr->Entries() > 0)
02813     {
02814         G4CollisionInitialState *
02815         collision = theCollisionMgr->GetNextCollision();
02816 
02817         if ( ! collision->GetTargetCollection().size() ){
02818             G4KineticTrackVector * lates = collision->GetFinalState();
02819             if ( lates->size() == 1 ) {
02820                 G4KineticTrack * atrack=*(lates->begin());
02821                 G4ReactionProduct * aNew = new G4ReactionProduct(atrack->GetDefinition());
02822                 aNew->SetMomentum(atrack->Get4Momentum().vect());
02823                 aNew->SetTotalEnergy(atrack->Get4Momentum().e());
02824                 Esecondaries +=atrack->Get4Momentum().e();
02825                 aNew->SetNewlyAdded(true);
02826                 //G4cout << " Particle Ekin " << aNew->GetKineticEnergy() << G4endl;
02827                 products->push_back(aNew);
02828 
02829             }
02830         }
02831         theCollisionMgr->RemoveCollision(collision);
02832 
02833     }
02834 
02835     // decay must be after loop on Collisions, and Decay() will delete entries in theSecondaryList, refered
02836     //   to by Collisions.
02837     decayKTV.Decay(&theSecondaryList);
02838     //PrintKTVector(&theSecondaryList, " secondaires @ void Nucl");
02839     for(iter = theSecondaryList.begin(); iter != theSecondaryList.end(); ++iter)
02840     {
02841         G4ReactionProduct * aNew = new G4ReactionProduct((*iter)->GetDefinition());
02842         aNew->SetMomentum((*iter)->Get4Momentum().vect());
02843         aNew->SetTotalEnergy((*iter)->Get4Momentum().e());
02844         Esecondaries +=(*iter)->Get4Momentum().e();
02845         aNew->SetNewlyAdded(true);
02846         //G4cout << " Particle Ekin " << aNew->GetKineticEnergy() << G4endl;
02847         products->push_back(aNew);
02848     }
02849 
02850     G4double SumMassNucleons(0.);
02851     for(iter = theTargetList.begin(); iter != theTargetList.end(); ++iter)
02852     {
02853         SumMassNucleons += (*iter)->GetDefinition()->GetPDGMass();
02854     }
02855 
02856     G4double Ekinetic=theProjectile4Momentum.e() + initial_nuclear_mass - Esecondaries - SumMassNucleons;
02857     //G4cout << "Void nucleus nucleons : "<<theTargetList.size() << " ,  T: " << Ekinetic << G4endl;
02858     if (Ekinetic > 0.){
02859         if (theTargetList.size() ) Ekinetic /= theTargetList.size();
02860     } else {
02861         //G4cout << " zero or negative kinetic energy, use random: " << Ekinetic<< G4endl;
02862         Ekinetic= ( 0.1 + G4UniformRand()*5.) * MeV;    // violate Energy conservation by small amount here
02863     }
02864     //G4cout << " using T per nucleon: " << Ekinetic << G4endl;
02865 
02866     for(iter = theTargetList.begin(); iter != theTargetList.end(); ++iter)
02867     {
02868         G4ReactionProduct * aNew = new G4ReactionProduct((*iter)->GetDefinition());
02869         G4double mass=(*iter)->GetDefinition()->GetPDGMass();
02870         G4double p=std::sqrt(sqr(Ekinetic) + 2.*Ekinetic*mass);
02871         aNew->SetMomentum(p*(*iter)->Get4Momentum().vect().unit());
02872         aNew->SetTotalEnergy(mass+Ekinetic);
02873         aNew->SetNewlyAdded(true);
02874         //G4cout << " Particle Ekin " << aNew->GetKineticEnergy() << G4endl;
02875         products->push_back(aNew);
02876     }
02877 
02878     return products;
02879 }
02880 G4ReactionProductVector * G4BinaryCascade::HighEnergyModelFSProducts(G4ReactionProductVector * products,
02881         G4KineticTrackVector * secondaries)
02882 {
02883     std::vector<G4KineticTrack *>::iterator iter;
02884     for(iter = secondaries->begin(); iter != secondaries->end(); ++iter)
02885     {
02886         G4ReactionProduct * aNew = new G4ReactionProduct((*iter)->GetDefinition());
02887         aNew->SetMomentum((*iter)->Get4Momentum().vect());
02888         aNew->SetTotalEnergy((*iter)->Get4Momentum().e());
02889         aNew->SetNewlyAdded(true);
02890         //G4cout << " Particle Ekin " << aNew->GetKineticEnergy() << G4endl;
02891         products->push_back(aNew);
02892     }
02893     G4ParticleDefinition* fragment = 0;
02894     if (currentA == 1 && currentZ == 0) {
02895         fragment = G4Neutron::NeutronDefinition();
02896     } else if (currentA == 1 && currentZ == 1) {
02897         fragment = G4Proton::ProtonDefinition();
02898     } else if (currentA == 2 && currentZ == 1) {
02899         fragment = G4Deuteron::DeuteronDefinition();
02900     } else if (currentA == 3 && currentZ == 1) {
02901         fragment = G4Triton::TritonDefinition();
02902     } else if (currentA == 3 && currentZ == 2) {
02903         fragment = G4He3::He3Definition();
02904     } else if (currentA == 4 && currentZ == 2) {
02905         fragment = G4Alpha::AlphaDefinition();;
02906     } else {
02907         fragment =
02908                 G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(currentZ,currentA,0.0);
02909     }
02910     if (fragment != 0) {
02911         G4ReactionProduct * theNew = new G4ReactionProduct(fragment);
02912         theNew->SetMomentum(G4ThreeVector(0,0,0));
02913         theNew->SetTotalEnergy(massInNucleus);
02914         //theNew->SetFormationTime(??0.??);
02915         //G4cout << " Nucleus (" << currentZ << ","<< currentA << "), mass "<< massInNucleus << G4endl;
02916         products->push_back(theNew);
02917     }
02918     return products;
02919 }
02920 
02921 void G4BinaryCascade::PrintWelcomeMessage()
02922 {
02923     G4cout <<"Thank you for using G4BinaryCascade. "<<G4endl;
02924 }
02925 
02926 //----------------------------------------------------------------------------
02927 void G4BinaryCascade::DebugApplyCollisionFail(G4CollisionInitialState * collision,
02928       G4KineticTrackVector * products)
02929 {
02930    G4bool havePion=false;
02931    if (products)
02932    {
02933       for ( std::vector<G4KineticTrack *>::iterator i =products->begin(); i != products->end(); i++)
02934       {
02935          G4int PDGcode=std::abs((*i)->GetDefinition()->GetPDGEncoding());
02936          if (std::abs(PDGcode)==211 || PDGcode==111 ) havePion=true;
02937       }
02938    }
02939    if ( !products  || havePion)
02940    {
02941       G4cout << " Collision " << collision << ", type: "<< typeid(*collision->GetGenerator()).name()
02942                                                                             << ", with NO products! " <<G4endl;
02943       G4cout << G4endl<<"Initial condition are these:"<<G4endl;
02944       G4cout << "proj: "<<collision->GetPrimary()->GetDefinition()->GetParticleName()<<G4endl;
02945       PrintKTVector(collision->GetPrimary());
02946       for(size_t it=0; it<collision->GetTargetCollection().size(); it++)
02947       {
02948          G4cout << "targ: "
02949                <<collision->GetTargetCollection()[it]->GetDefinition()->GetParticleName()<<G4endl;
02950       }
02951       PrintKTVector(&collision->GetTargetCollection(),std::string(" Target particles"));
02952    }
02953    //  if ( lateParticleCollision ) G4cout << " Added late particle--------------------------"<<G4endl;
02954    //  if ( lateParticleCollision && products ) PrintKTVector(products, " reaction products");
02955 }
02956 
02957 //----------------------------------------------------------------------------
02958 
02959 G4bool G4BinaryCascade::CheckChargeAndBaryonNumber(G4String where)
02960 {
02961    static G4int lastdA(0), lastdZ(0);
02962    G4int iStateA = the3DNucleus->GetMassNumber() + projectileA;
02963    G4int iStateZ = the3DNucleus->GetCharge()  + projectileZ;
02964 
02965    G4int fStateA(0);
02966    G4int fStateZ(0);
02967 
02968    std::vector<G4KineticTrack *>::iterator i;
02969    G4int CapturedA(0), CapturedZ(0);
02970    G4int secsA(0), secsZ(0);
02971    for ( i=theCapturedList.begin(); i!=theCapturedList.end(); ++i) {
02972       CapturedA += (*i)->GetDefinition()->GetBaryonNumber();
02973       CapturedZ += G4lrint((*i)->GetDefinition()->GetPDGCharge()/eplus);
02974    }
02975 
02976    for ( i=theSecondaryList.begin(); i!=theSecondaryList.end(); ++i) {
02977       if ( (*i)->GetState() != G4KineticTrack::inside ) {
02978          secsA += (*i)->GetDefinition()->GetBaryonNumber();
02979          secsZ += G4lrint((*i)->GetDefinition()->GetPDGCharge()/eplus);
02980       }
02981    }
02982 
02983    for ( i=theFinalState.begin(); i!=theFinalState.end(); ++i) {
02984       fStateA += (*i)->GetDefinition()->GetBaryonNumber();
02985       fStateZ += G4lrint((*i)->GetDefinition()->GetPDGCharge()/eplus);
02986    }
02987 
02988    G4int deltaA= iStateA -  secsA - fStateA -currentA - lateA;
02989    G4int deltaZ= iStateZ -  secsZ - fStateZ -currentZ - lateZ;
02990 
02991    if (deltaA != 0  || deltaZ!=0 ) {
02992       if (deltaA != lastdA || deltaZ != lastdZ ) {
02993          G4cout << "baryon/charge imbalance - " << where << G4endl
02994                << "deltaA " <<deltaA<<", iStateA "<<iStateA<< ",  CapturedA "<<CapturedA <<",  secsA "<<secsA
02995                << ", fStateA "<<fStateA << ", currentA "<<currentA << ", lateA "<<lateA << G4endl
02996                << "deltaZ "<<deltaZ<<", iStateZ "<<iStateZ<< ",  CapturedZ "<<CapturedZ <<",  secsZ "<<secsZ
02997                << ", fStateZ "<<fStateZ << ", currentZ "<<currentZ << ", lateZ "<<lateZ << G4endl<< G4endl;
02998          lastdA=deltaA;
02999          lastdZ=deltaZ;
03000       }
03001    } else { lastdA=lastdZ=0;}
03002 
03003    return true;
03004 }
03005 //----------------------------------------------------------------------------
03006 void G4BinaryCascade::DebugApplyCollision(G4CollisionInitialState * collision,
03007         G4KineticTrackVector * products)
03008 {
03009 
03010     PrintKTVector(collision->GetPrimary(),std::string(" Primary particle"));
03011     PrintKTVector(&collision->GetTargetCollection(),std::string(" Target particles"));
03012     PrintKTVector(products,std::string(" Scatterer products"));
03013 
03014 #ifdef dontUse
03015     G4double thisExcitation(0);
03016     //  excitation energy from this collision
03017     //  initial state:
03018     G4double initial(0);
03019     G4KineticTrack * kt=collision->GetPrimary();
03020     initial +=  kt->Get4Momentum().e();
03021 
03022     G4RKPropagation * RKprop=(G4RKPropagation *)thePropagator;
03023 
03024     initial +=  RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition());
03025     initial -=  RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding());
03026     G4cout << "prim. E/field/Barr/Sum " << kt->Get4Momentum().e()
03027                                   << " " << RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition())
03028                                   << " " << RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding())
03029                                   << " " << initial << G4endl;;
03030 
03031     G4KineticTrackVector ktv=collision->GetTargetCollection();
03032     for ( unsigned int it=0; it < ktv.size(); it++)
03033     {
03034         kt=ktv[it];
03035         initial +=  kt->Get4Momentum().e();
03036         thisExcitation += kt->GetDefinition()->GetPDGMass()
03037                                          - kt->Get4Momentum().e()
03038                                          - RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition());
03039         //     initial +=  RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition());
03040         //     initial -=  RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding());
03041         G4cout << "Targ. def/E/field/Barr/Sum " <<  kt->GetDefinition()->GetPDGEncoding()
03042                                   << " " << kt->Get4Momentum().e()
03043                                   << " " << RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition())
03044                                   << " " << RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding())
03045                                   << " " << initial <<" Excit " << thisExcitation << G4endl;;
03046     }
03047 
03048     G4double final(0);
03049     G4double mass_out(0);
03050     G4int product_barions(0);
03051     if ( products )
03052     {
03053         for ( unsigned int it=0; it < products->size(); it++)
03054         {
03055             kt=(*products)[it];
03056             final +=  kt->Get4Momentum().e();
03057             final +=  RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition());
03058             final +=  RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding());
03059             if ( kt->GetDefinition()->GetBaryonNumber()==1 ) product_barions++;
03060             mass_out += kt->GetDefinition()->GetPDGMass();
03061             G4cout << "sec. def/E/field/Barr/Sum " << kt->GetDefinition()->GetPDGEncoding()
03062                                  << " " << kt->Get4Momentum().e()
03063                                  << " " << RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition())
03064                                  << " " << RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding())
03065                                  << " " << final << G4endl;;
03066         }
03067     }
03068 
03069 
03070     G4int finalA = currentA;
03071     G4int finalZ = currentZ;
03072     if ( products )
03073     {
03074         finalA -= product_barions;
03075         finalZ -= GetTotalCharge(*products);
03076     }
03077     G4double delta = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(currentZ,currentA)
03078                                            - (G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(finalZ,finalA)
03079                                                    + mass_out);
03080     G4cout << " current/final a,z " << currentA << " " << currentZ << " "<< finalA<< " "<< finalZ
03081             <<  " delta-mass " << delta<<G4endl;
03082     final+=delta;
03083     mass_out  = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(finalZ,finalA);
03084     G4cout << " initE/ E_out/ Mfinal/ Excit " << currentInitialEnergy
03085             << " " <<   final << " "
03086             <<  mass_out<<" "
03087             <<  currentInitialEnergy - final - mass_out
03088             << G4endl;
03089     currentInitialEnergy-=final;
03090 #endif
03091 }
03092 
03093 //----------------------------------------------------------------------------
03094 G4bool G4BinaryCascade::DebugEpConservation(const G4HadProjectile & aTrack,
03095         G4ReactionProductVector* products)
03096 //----------------------------------------------------------------------------
03097 {
03098     G4ReactionProductVector::iterator iter;
03099     G4double Efinal(0);
03100     G4ThreeVector pFinal(0);
03101     if (std::abs(theParticleChange.GetWeightChange() -1 ) > 1e-5 )
03102     {
03103         G4cout <<" BIC-weight change " << theParticleChange.GetWeightChange()<< G4endl;
03104     }
03105 
03106     for(iter = products->begin(); iter != products->end(); ++iter)
03107     {
03108 
03109         //   G4cout << " Secondary E - Ekin / p " <<
03110         //      (*iter)->GetDefinition()->GetParticleName() << " " <<
03111         //      (*iter)->GetTotalEnergy() << " - " <<
03112         //      (*iter)->GetKineticEnergy()<< " / " <<
03113         //      (*iter)->GetMomentum().x() << " " <<
03114         //      (*iter)->GetMomentum().y() << " " <<
03115         //      (*iter)->GetMomentum().z() << G4endl;
03116         Efinal += (*iter)->GetTotalEnergy();
03117         pFinal += (*iter)->GetMomentum();
03118     }
03119 
03120     //  G4cout << "e outgoing/ total : " << Efinal << " " << Efinal+GetFinal4Momentum().e()<< G4endl;
03121     G4cout << "BIC E/p delta " <<
03122             (aTrack.Get4Momentum().e()+ the3DNucleus->GetMass() - Efinal)/MeV <<
03123             " MeV / mom " << (aTrack.Get4Momentum()  - pFinal ) /MeV << G4endl;
03124 
03125     return (aTrack.Get4Momentum().e() + the3DNucleus->GetMass() - Efinal)/aTrack.Get4Momentum().e() < perCent;
03126 
03127 }
03128 
03129 //----------------------------------------------------------------------------
03130 G4ReactionProductVector * G4BinaryCascade::ProductsAddFakeGamma(G4ReactionProductVector *products )
03131 //----------------------------------------------------------------------------
03132 {
03133    //    else
03134 //    {
03135 //        G4ReactionProduct * aNew=0;
03136 //        // return nucleus e and p
03137 //        if  (fragment != 0 ) {
03138 //            aNew = new G4ReactionProduct(G4Gamma::GammaDefinition());   // we only want to pass e/p
03139 //            aNew->SetMomentum(fragment->GetMomentum().vect());
03140 //            aNew->SetTotalEnergy(fragment->GetMomentum().e());
03141 //            delete fragment;
03142 //            fragment=0;
03143 //        } else if (products->size() == 0) {
03144 //            // FixMe GF: for testing without precompound, return 1 gamma of 0.01 MeV in +x
03145 //#include "G4Gamma.hh"
03146 //            aNew = new G4ReactionProduct(G4Gamma::GammaDefinition());
03147 //            aNew->SetMomentum(G4ThreeVector(0.01*MeV,0,0));
03148 //            aNew->SetTotalEnergy(0.01*MeV);
03149 //        }
03150 //        if ( aNew != 0 ) products->push_back(aNew);
03151 //    }
03152    return products;
03153 }
03154 
03155 //----------------------------------------------------------------------------

Generated on Mon May 27 17:47:45 2013 for Geant4 by  doxygen 1.4.7