G4BinaryCascade Class Reference

#include <G4BinaryCascade.hh>

Inheritance diagram for G4BinaryCascade:

Public Member Functions
	G4BinaryCascade (G4VPreCompoundModel *ptr=0)
virtual	~G4BinaryCascade ()
G4HadFinalState *	ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &theNucleus)
virtual G4ReactionProductVector *	Propagate (G4KineticTrackVector *, G4V3DNucleus *)
virtual void	ModelDescription (std::ostream &) const
virtual void	PropagateModelDescription (std::ostream &) const

---

Detailed Description

Definition at line 70 of file G4BinaryCascade.hh.

---

Constructor & Destructor Documentation

G4BinaryCascade::G4BinaryCascade

(

G4VPreCompoundModel *

ptr = 0

)

Definition at line 102 of file G4BinaryCascade.cc.

References G4ShortLivedConstructor::ConstructParticle(), G4HadronicInteractionRegistry::FindModel(), G4VIntraNuclearTransportModel::GetDeExcitation(), G4VPreCompoundModel::GetExcitationHandler(), G4HadronicInteractionRegistry::Instance(), G4VIntraNuclearTransportModel::SetDeExcitation(), G4HadronicInteraction::SetEnergyMomentumCheckLevels(), G4HadronicInteraction::SetMaxEnergy(), and G4HadronicInteraction::SetMinEnergy().

                                                         :
G4VIntraNuclearTransportModel("Binary Cascade", ptr)
{
    // initialise the resonance sector
    G4ShortLivedConstructor ShortLived;
    ShortLived.ConstructParticle();

    theCollisionMgr = new G4CollisionManager;
    theDecay=new G4BCDecay;
    theImR.push_back(theDecay);
    theLateParticle= new G4BCLateParticle;
    G4MesonAbsorption * aAb=new G4MesonAbsorption;
    theImR.push_back(aAb);
    G4Scatterer * aSc=new G4Scatterer;
    theH1Scatterer = new G4Scatterer;
    theImR.push_back(aSc);

    thePropagator = new G4RKPropagation;
    theCurrentTime = 0.;
    theBCminP = 45*MeV;
    theCutOnP = 90*MeV;
    theCutOnPAbsorb= 0*MeV;   // No Absorption of slow Mesons, other than above G4MesonAbsorption

    // reuse existing pre-compound model
    if(!ptr) {
      G4HadronicInteraction* p =
        G4HadronicInteractionRegistry::Instance()->FindModel("PRECO");
      G4VPreCompoundModel* pre = static_cast<G4VPreCompoundModel*>(p);
      if(!pre) { pre = new G4PreCompoundModel(); }
      SetDeExcitation(pre);
    }
    theExcitationHandler = GetDeExcitation()->GetExcitationHandler();
    SetMinEnergy(0.0*GeV);
    SetMaxEnergy(10.1*GeV);
    //PrintWelcomeMessage();
    thePrimaryEscape = true;
    thePrimaryType = 0;

    SetEnergyMomentumCheckLevels(1*perCent, 1*MeV);

    // init data members
    currentA=currentZ=0;
    lateA=lateZ=0;
    initialA=initialZ=0;
    projectileA=projectileZ=0;
    currentInitialEnergy=initial_nuclear_mass=0.;
    massInNucleus=0.;
    theOuterRadius=0.;
}

G4BinaryCascade::~G4BinaryCascade

(

)

[virtual]

Definition at line 159 of file G4BinaryCascade.cc.

{
    ClearAndDestroy(&theTargetList);
    ClearAndDestroy(&theSecondaryList);
    ClearAndDestroy(&theCapturedList);
    delete thePropagator;
    delete theCollisionMgr;
    std::for_each(theImR.begin(), theImR.end(), Delete<G4BCAction>());
    delete theLateParticle;
    //delete theExcitationHandler;
    delete theH1Scatterer;
}

---

Member Function Documentation

G4HadFinalState * G4BinaryCascade::ApplyYourself	(	const G4HadProjectile &	aTrack,
		G4Nucleus &	theNucleus
	)			[virtual]

Implements G4HadronicInteraction.

Definition at line 208 of file G4BinaryCascade.cc.

References G4HadFinalState::AddSecondary(), G4VPreCompoundModel::ApplyYourself(), G4HadFinalState::Clear(), G4CollisionManager::ClearAndDestroy(), G4cerr, G4endl, G4HadProjectile::Get4Momentum(), G4Nucleus::GetA_asInt(), G4HadProjectile::GetDefinition(), G4HadProjectile::GetKineticEnergy(), G4V3DNucleus::GetOuterRadius(), G4Nucleus::GetZ_asInt(), G4VFieldPropagation::Init(), G4V3DNucleus::Init(), isAlive, G4Neutron::NeutronDefinition(), G4KineticTrack::outside, G4PionMinus::PionMinusDefinition(), G4PionPlus::PionPlusDefinition(), Propagate(), G4Proton::ProtonDefinition(), G4HadFinalState::SetEnergyChange(), G4HadFinalState::SetMomentumChange(), G4KineticTrack::SetState(), G4HadFinalState::SetStatusChange(), stopAndKill, G4VIntraNuclearTransportModel::the3DNucleus, G4VIntraNuclearTransportModel::theDeExcitation, and G4HadronicInteraction::theParticleChange.

{
    static G4int eventcounter=0;

    //   if ( eventcounter == 0 ) {
    //      SetEpReportLevel(3);   // report non conservation with model etc.
    //      G4double relativeLevel = 1*perCent;
    //      G4double absoluteLevel = 2*MeV;
    //      SetEnergyMomentumCheckLevels(relativeLevel,absoluteLevel);
    //   }

    //if(eventcounter == 100*(eventcounter/100) )
    eventcounter++;
    if(getenv("BCDEBUG") ) G4cerr << " ######### Binary Cascade Reaction number starts ######### "<<eventcounter<<G4endl;

    G4LorentzVector initial4Momentum = aTrack.Get4Momentum();
    G4ParticleDefinition * definition = const_cast<G4ParticleDefinition *>(aTrack.GetDefinition());

    if(initial4Momentum.e()-initial4Momentum.m()<theBCminP &&
            ( definition==G4Neutron::NeutronDefinition() || definition==G4Proton::ProtonDefinition() ) )
    {
        return theDeExcitation->ApplyYourself(aTrack, aNucleus);
    }

    theParticleChange.Clear();
    // initialize the G4V3DNucleus from G4Nucleus
    the3DNucleus = new G4Fancy3DNucleus;

    // Build a KineticTrackVector with the G4Track
    G4KineticTrackVector * secondaries;// = new G4KineticTrackVector;
    G4ThreeVector initialPosition(0., 0., 0.); // will be set later

    if(!getenv("I_Am_G4BinaryCascade_Developer") )
    {
        if(definition!=G4Neutron::NeutronDefinition() &&
                definition!=G4Proton::ProtonDefinition() &&
                definition!=G4PionPlus::PionPlusDefinition() &&
                definition!=G4PionMinus::PionMinusDefinition() )
        {
            G4cerr << "You are using G4BinaryCascade for projectiles other than nucleons or pions."<<G4endl;
            G4cerr << "If you want to continue, please switch on the developer environment: "<<G4endl;
            G4cerr << "setenv I_Am_G4BinaryCascade_Developer 1 "<<G4endl<<G4endl;
            throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade - used for unvalid particle type - Fatal");
        }
    }

    // keep primary
    thePrimaryType = definition;
    thePrimaryEscape = false;

    // try until an interaction will happen
    G4ReactionProductVector * products = 0;
    G4int interactionCounter = 0;
    do
    {
        // reset status that could be changed in previous loop event
        theCollisionMgr->ClearAndDestroy();

        if(products != 0)
        {  // free memory from previous loop event
            ClearAndDestroy(products);
            delete products;
            products=0;
        }

        G4int massNumber=aNucleus.GetA_asInt();
        the3DNucleus->Init(massNumber, aNucleus.GetZ_asInt());
        thePropagator->Init(the3DNucleus);
        //      GF Leak on kt??? but where to delete?
        G4KineticTrack * kt;// = new G4KineticTrack(definition, 0., initialPosition, initial4Momentum);
        do                  // sample impact parameter until collisions are found
        {
            theCurrentTime=0;
            G4double radius = the3DNucleus->GetOuterRadius()+3*fermi;
            initialPosition=GetSpherePoint(1.1*radius, initial4Momentum);  // get random position
            kt = new G4KineticTrack(definition, 0., initialPosition, initial4Momentum);
            kt->SetState(G4KineticTrack::outside);
            // secondaries has been cleared by Propagate() in the previous loop event
            secondaries= new G4KineticTrackVector;
            secondaries->push_back(kt);
            if(massNumber > 1) // 1H1 is special case
            {
                products = Propagate(secondaries, the3DNucleus);
            } else {
                products = Propagate1H1(secondaries,the3DNucleus);
            }
        } while(! products );  // until we FIND a collision...

        if(++interactionCounter>99) break;
    } while(products->size() == 0);  // ...until we find an ALLOWED collision

    if(products->size()>0)
    {
        //  G4cout << "BIC Applyyourself: number of products " << products->size() << G4endl;

        // Fill the G4ParticleChange * with products
        theParticleChange.SetStatusChange(stopAndKill);
        G4ReactionProductVector::iterator iter;

        for(iter = products->begin(); iter != products->end(); ++iter)
        {
            G4DynamicParticle * aNew =
                    new G4DynamicParticle((*iter)->GetDefinition(),
                            (*iter)->GetTotalEnergy(),
                            (*iter)->GetMomentum());
            theParticleChange.AddSecondary(aNew);
        }

        // DebugEpConservation(track, products);


    } else {  // no interaction, return primary
        if(getenv("BCDEBUG") ) G4cerr << " ######### Binary Cascade Reaction number void ######### "<<eventcounter<<G4endl;
        theParticleChange.SetStatusChange(isAlive);
        theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
        theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
    }

    ClearAndDestroy(products);
    delete products;

    delete the3DNucleus;
    the3DNucleus = NULL;

    if(getenv("BCDEBUG") ) G4cerr << " ######### Binary Cascade Reaction number ends ######### "<<eventcounter<<G4endl;

    return &theParticleChange;
}

void G4BinaryCascade::ModelDescription

(

std::ostream &

)

const [virtual]

Reimplemented from G4VIntraNuclearTransportModel.

Definition at line 172 of file G4BinaryCascade.cc.

{
    outFile << "G4BinaryCascade is an intra-nuclear cascade model in which\n"
            << "an incident hadron collides with a nucleon, forming two\n"
            << "final-state particles, one or both of which may be resonances.\n"
            << "The resonances then decay hadronically and the decay products\n"
            << "are then propagated through the nuclear potential along curved\n"
            << "trajectories until they re-interact or leave the nucleus.\n"
            << "This model is valid for incident pions up to 1.5 GeV and\n"
            << "nucleons up to 10 GeV.\n";
}

G4ReactionProductVector * G4BinaryCascade::Propagate	(	G4KineticTrackVector *	,
		G4V3DNucleus *
	)			[virtual]

Implements G4VIntraNuclearTransportModel.

Definition at line 339 of file G4BinaryCascade.cc.

References _CheckChargeAndBaryonNumber_, G4CollisionManager::ClearAndDestroy(), G4CollisionManager::Entries(), G4cerr, G4cout, G4endl, G4V3DNucleus::GetMass(), G4CollisionManager::GetNextCollision(), G4V3DNucleus::GetOuterRadius(), G4VFieldPropagation::Init(), G4CollisionManager::RemoveCollision(), and G4VIntraNuclearTransportModel::the3DNucleus.

Referenced by ApplyYourself().

{
    G4ping debug("debug_G4BinaryCascade");
#ifdef debug_BIC_Propagate
    G4cout << "G4BinaryCascade Propagate starting -------------------------------------------------------" <<G4endl;
#endif

    the3DNucleus=aNucleus;
    G4ReactionProductVector * products = new G4ReactionProductVector;
    theOuterRadius = the3DNucleus->GetOuterRadius();
    theCurrentTime=0;
    theProjectile4Momentum=G4LorentzVector(0,0,0,0);
    // build theSecondaryList, theProjectileList and theCapturedList
    ClearAndDestroy(&theCapturedList);
    ClearAndDestroy(&theSecondaryList);
    theSecondaryList.clear();
    ClearAndDestroy(&theFinalState);
    std::vector<G4KineticTrack *>::iterator iter;
    theCollisionMgr->ClearAndDestroy();

    theCutOnP=90*MeV;
    if(the3DNucleus->GetMass()>30) theCutOnP = 70*MeV;
    if(the3DNucleus->GetMass()>60) theCutOnP = 50*MeV;
    if(the3DNucleus->GetMass()>120) theCutOnP = 45*MeV;


    BuildTargetList();

#ifdef debug_BIC_GetExcitationEnergy
    G4cout << "ExcitationEnergy0 " << GetExcitationEnergy() << G4endl;
#endif

    thePropagator->Init(the3DNucleus);

    G4bool success = BuildLateParticleCollisions(secondaries);
    if (! success )   // fails if no excitation energy left....
    {
       products=HighEnergyModelFSProducts(products, secondaries);
       ClearAndDestroy(secondaries);
       delete secondaries;

#ifdef debug_G4BinaryCascade
       G4cout << "G4BinaryCascade::Propagate: warning - high energy model failed energy conservation, returning unchanged high energy final state" << G4endl;
#endif

       return products;
    }
    // check baryon and charge ...

    _CheckChargeAndBaryonNumber_("lateparticles");

    // if called stand alone find first collisions
    FindCollisions(&theSecondaryList);


    if(theCollisionMgr->Entries() == 0 )      //late particles ALWAYS create Entries
    {
        //G4cout << " no collsions -> return 0" << G4endl;
        delete products;
#ifdef debug_BIC_return
        G4cout << "return @ begin2,  no collisions "<< G4endl;
#endif
        return 0;
    }

    // end of initialization: do the job now
    // loop until there are no more collisions


    G4bool haveProducts = false;
    G4int collisionCount=0;
    theMomentumTransfer=G4ThreeVector(0,0,0);
    while(theCollisionMgr->Entries() > 0 && currentZ)
    {
        if(Absorb()) {  // absorb secondaries, pions only
            haveProducts = true;
        }
        _CheckChargeAndBaryonNumber_("absorbed");
        if(Capture()) { // capture secondaries, nucleons only
            haveProducts = true;
        }
        _CheckChargeAndBaryonNumber_("captured");

        // propagate to the next collision if any (collisions could have been deleted
        // by previous absorption or capture)
        if(theCollisionMgr->Entries() > 0)
        {
            G4CollisionInitialState *
            nextCollision = theCollisionMgr->GetNextCollision();
#ifdef debug_BIC_Propagate_Collisions
            G4cout << " NextCollision  * , Time, curtime = " << nextCollision << " "
                    <<nextCollision->GetCollisionTime()<< " " <<
                    theCurrentTime<< G4endl;
#endif
            if (!DoTimeStep(nextCollision->GetCollisionTime()-theCurrentTime) )
            {
                // Check if nextCollision is still valid, ie. particle did not leave nucleus
                if (theCollisionMgr->GetNextCollision() != nextCollision )
                {
                    nextCollision = 0;
                }
            }

            if( nextCollision )
            {
                if (ApplyCollision(nextCollision))
                {
                    //G4cerr << "ApplyCollision success " << G4endl;
                    haveProducts = true;
                    collisionCount++;
                    _CheckChargeAndBaryonNumber_("ApplyCollision");

                } else {
                    //G4cerr << "ApplyCollision failure " << G4endl;
                    theCollisionMgr->RemoveCollision(nextCollision);
                }
            }
        }
    }

    //--------- end of while on Collisions
    //G4cout << "currentZ @ end loop " << currentZ << G4endl;
    if ( ! currentZ ){
        // nucleus completely destroyed, fill in ReactionProductVector
        products = FillVoidNucleusProducts(products);
#ifdef debug_BIC_return
        G4cout << "return @ Z=0 after collision loop "<< G4endl;
        PrintKTVector(&theSecondaryList,std::string(" theSecondaryList"));
        G4cout << "theTargetList size: " << theTargetList.size() << G4endl;
        PrintKTVector(&theTargetList,std::string(" theTargetList"));
        PrintKTVector(&theCapturedList,std::string(" theCapturedList"));

        G4cout << " ExcitE be4 Correct : " <<GetExcitationEnergy() << G4endl;
        G4cout << " Mom Transfered to nucleus : " << theMomentumTransfer << " " << theMomentumTransfer.mag() << G4endl;
        PrintKTVector(&theFinalState,std::string(" FinalState uncorrected"));
        G4cout << "returned products: " << products->size() << G4endl;
#endif
        // @@GF FixMe cannot return here.
        return products;
    }

    // No more collisions: absorb, capture and propagate the secondaries out of the nucleus
    if(Absorb()) {
        haveProducts = true;
        // G4cout << "Absorb sucess " << G4endl;
    }

    if(Capture()) {
        haveProducts = true;
        // G4cout << "Capture sucess " << G4endl;
    }

    if(!haveProducts)  // no collisions happened. Return an empty vector.
    {
#ifdef debug_BIC_return
        G4cout << "return 3, no products "<< G4endl;
#endif
        return products;
    }


#ifdef debug_BIC_Propagate
    G4cout << " Momentum transfer to Nucleus " << theMomentumTransfer << " " << theMomentumTransfer.mag() << G4endl;
    G4cout << "  Stepping particles out...... " << G4endl;
#endif

    StepParticlesOut();


    if ( theSecondaryList.size() > 0 )
    {
#ifdef debug_G4BinaryCascade
        G4cerr << "G4BinaryCascade: Warning, have active particles at end" << G4endl;
        PrintKTVector(&theSecondaryList, "active particles @ end  added to theFinalState");
#endif
        //  add left secondaries to FinalSate
        for ( iter =theSecondaryList.begin(); iter != theSecondaryList.end(); ++iter)
        {
            theFinalState.push_back(*iter);
        }
        theSecondaryList.clear();

    }
    while ( theCollisionMgr->Entries() > 0 )
    {
#ifdef debug_G4BinaryCascade
        G4cerr << " Warning: remove left over collision(s) " << G4endl;
#endif
        theCollisionMgr->RemoveCollision(theCollisionMgr->GetNextCollision());
    }

#ifdef debug_BIC_Propagate_Excitation

    PrintKTVector(&theSecondaryList,std::string(" theSecondaryList"));
    G4cout << "theTargetList size: " << theTargetList.size() << G4endl;
    //  PrintKTVector(&theTargetList,std::string(" theTargetList"));
    PrintKTVector(&theCapturedList,std::string(" theCapturedList"));

    G4cout << " ExcitE be4 Correct : " <<GetExcitationEnergy() << G4endl;
    G4cout << " Mom Transfered to nucleus : " << theMomentumTransfer << " " << theMomentumTransfer.mag() << G4endl;
    PrintKTVector(&theFinalState,std::string(" FinalState uncorrected"));
#endif

    //


    G4double ExcitationEnergy=GetExcitationEnergy();

#ifdef debug_BIC_Propagate_finals
    PrintKTVector(&theFinalState,std::string(" FinalState be4 corr"));
    G4cout << " Excitation Energy prefinal,  #collisions:, out, captured  "
            << ExcitationEnergy << " "
            << collisionCount << " "
            << theFinalState.size() << " "
            << theCapturedList.size()<<G4endl;
#endif

    if (ExcitationEnergy < 0 )
    {
        G4int maxtry=5, ntry=0;
        do {
            CorrectFinalPandE();
            ExcitationEnergy=GetExcitationEnergy();
        } while ( ++ntry < maxtry && ExcitationEnergy < 0 );
    }

#ifdef debug_BIC_Propagate_finals
    PrintKTVector(&theFinalState,std::string(" FinalState corrected"));
    G4cout << " Excitation Energy final,  #collisions:, out, captured  "
            << ExcitationEnergy << " "
            << collisionCount << " "
            << theFinalState.size() << " "
            << theCapturedList.size()<<G4endl;
#endif


    if ( ExcitationEnergy < 0. )
    {
        //      if ( ExcitationEnergy < 0. )
        {
            //#ifdef debug_G4BinaryCascade
            //            G4cerr << "G4BinaryCascade-Warning: negative excitation energy ";
            //            G4cerr <<ExcitationEnergy<<G4endl;
            //     PrintKTVector(&theFinalState,std::string("FinalState"));
            //    PrintKTVector(&theCapturedList,std::string("captured"));
            //    G4cout << "negative ExE:Final 4Momentum .mag: " << GetFinal4Momentum()
            //            << " "<< GetFinal4Momentum().mag()<< G4endl
            //            << "negative ExE:FinalNucleusMom  .mag: " << GetFinalNucleusMomentum()
            //            << " "<< GetFinalNucleusMomentum().mag()<< G4endl;
            //#endif
        }
        ClearAndDestroy(products);
        //G4cout << "  negative Excitation E return empty products " << products << G4endl;
#ifdef debug_BIC_return
        G4cout << "return 4, excit < 0 "<< G4endl;
#endif
        return products;   // return empty products- FIXME
    }

    G4ReactionProductVector * precompoundProducts=DeExcite();


    G4DecayKineticTracks decay(&theFinalState);

    products= ProductsAddFinalState(products, theFinalState);

    products= ProductsAddPrecompound(products, precompoundProducts);

//    products=ProductsAddFakeGamma(products);


    thePrimaryEscape = true;

    #ifdef debug_BIC_return
    G4cout << "return @end, all ok "<< G4endl;
    //G4cout << "  return products " << products << G4endl;
    #endif

    return products;
}

void G4BinaryCascade::PropagateModelDescription

(

std::ostream &

)

const [virtual]

Reimplemented from G4VIntraNuclearTransportModel.

Definition at line 183 of file G4BinaryCascade.cc.

{
    outFile << "G4BinaryCascade propagtes secondaries produced by a high\n"
            << "energy model through the wounded nucleus.\n"
            << "Secondaries are followed after the formation time and if\n"
            << "within the nucleus are propagated through the nuclear\n"
            << "potential along curved trajectories until they interact\n"
            << "with a nucleon, decay, or leave the nucleus.\n"
            << "An interaction of a secondary with a nucleon produces two\n"
            << "final-state particles, one or both of which may be resonances.\n"
            << "Resonances decay hadronically and the decay products\n"
            << "are in turn propagated through the nuclear potential along curved\n"
            << "trajectories until they re-interact or leave the nucleus.\n"
            << "This model is valid for pions up to 1.5 GeV and\n"
            << "nucleons up to about 3.5 GeV.\n";
}

---

The documentation for this class was generated from the following files:

• G4BinaryCascade.hh
• G4BinaryCascade.cc

---