gar::garg4::MuNuclearSplittingProcessXSecBias Class Reference

#include <MuNuclearSplittingProcessXSecBias.h>

Inheritance diagram for gar::garg4::MuNuclearSplittingProcessXSecBias:

Public Member Functions
	MuNuclearSplittingProcessXSecBias ()
	~MuNuclearSplittingProcessXSecBias ()
void	SetNSplit (G4int nTrx, G4int xB=0, G4double xFac=1)
void	SetIsActive (G4bool doIt)
G4VParticleChange *	PostStepDoIt (const G4Track &track, const G4Step &step)
G4VParticleChange *	AlongStepDoIt (const G4Track &track, const G4Step &step)
virtual G4double	AlongStepGetPhysicalInteractionLength (const G4Track &track, G4double previousStepSize, G4double currentMinimumStep, G4double &proposedSafety, G4GPILSelection *selection)
virtual G4double	PostStepGetPhysicalInteractionLength (const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)

Protected Member Functions
virtual void	ResetNumberOfInteractionLengthLeft ()

Private Member Functions
G4double	XBiasSurvivalProbability ()
G4double	XBiasSecondaryWeight ()
G4double	GetTotalNumberOfInteractionLengthTraversed ()

Private Attributes
G4int	fNSplit
G4bool	fActive
G4int	xBiasMode
G4double	eFactor
G4VParticleChange	fParticleChange
G4double	wc
G4double	theInitialNumberOfInteractionLength

Detailed Description

Definition at line 31 of file MuNuclearSplittingProcessXSecBias.h.

Constructor & Destructor Documentation

gar::garg4::MuNuclearSplittingProcessXSecBias::MuNuclearSplittingProcessXSecBias ( )

inline

Definition at line 34 of file MuNuclearSplittingProcessXSecBias.h.

{};

gar::garg4::MuNuclearSplittingProcessXSecBias::~MuNuclearSplittingProcessXSecBias ( )

inline

Definition at line 35 of file MuNuclearSplittingProcessXSecBias.h.

{};

Member Function Documentation

G4VParticleChange * gar::garg4::MuNuclearSplittingProcessXSecBias::AlongStepDoIt	(	const G4Track &	track,
		const G4Step &	step
	)

Definition at line 170 of file MuNuclearSplittingProcessXSecBias.cxx.

    {
      
      wc = 1.0;
      
      if (xBiasMode==2)
      {
        /* I don't understand this. If I new it and Initialize(track) it works, but
         it goes off and grinds through the event.
         If I try to use pRegProcess->AlongStepDoIt() it crashes
         at the GetParentWeight() just below.
         */
        G4VParticleChange* pC  = new G4VParticleChange();
        pC->Initialize(track);
          //pC = pRegProcess->AlongStepDoIt( track, step);
        
        G4double nw = pC->GetParentWeight()*XBiasSecondaryWeight();
        pC->ProposeParentWeight(nw) ;
        for (G4int i = 0; i < pC->GetNumberOfSecondaries(); i++)
          pC->GetSecondary(i)->SetWeight(nw);
          //
        if(G4UniformRand()<XBiasSurvivalProbability())
        {
          pC->SetNumberOfSecondaries(pC->GetNumberOfSecondaries()+1);
          G4ThreeVector mDirection = track.GetDynamicParticle()->GetMomentumDirection();
          G4DynamicParticle* aD = new G4DynamicParticle((track.GetDynamicParticle())->GetDefinition(),
                                                        G4LorentzVector(mDirection,
                                                                        (track.GetDynamicParticle())->GetKineticEnergy()+step.GetTotalEnergyDeposit()));
          G4Track* secTrk = new G4Track(aD,track.GetGlobalTime(),track.GetPosition()); // +step.GetDeltaTime(), ,+step.GetDeltaPosition()
          secTrk->SetGoodForTrackingFlag(true);
          pC->AddSecondary(secTrk); // Add the primary right back.
          pC->GetSecondary(pC->GetNumberOfSecondaries()-1)->SetWeight(XBiasSurvivalProbability()*track.GetWeight());
            //pC->ProposeTrackStatus(fStopAndKill); // 8-June-2011-21:36 MDT, from Finland, I (EDC) added this. Seems right. Makes it go, anyway. Almost too fast.
        }
        
        return pC;
      }
      else if (xBiasMode==1)
      {
        fParticleChange.Initialize(track) ;
        G4double s = step.GetStepLength();
        G4double x = pRegProcess->GetCurrentInteractionLength();
          //fParticleChange = &(pRegProcess->AlongStepDoIt( track, step));
        wc =  exp(-(1. - eFactor)*s/x);
        G4double w = wc * fParticleChange.GetParentWeight();
        fParticleChange.SetParentWeightByProcess(false);
        fParticleChange.ProposeParentWeight(w) ;
        if (w>100)
        {
            //    G4cout << " MNSPXSB.AlongStepDoit(): GPW is " << w/wc << " wc = " << wc << ", CIL = "<< x << G4endl;
        }
        
        return &fParticleChange;
      }
      
      if (xBiasMode!=1 && xBiasMode!=2)
      {
        throw cet::exception("Incorrectly set Bias Mode. ")
        << "Set XBiasMode to 0 or 1. " << xBiasMode << " not allowed!\n";
      }
      fParticleChange.Initialize(track) ;
      return &fParticleChange;
    }

G4double gar::garg4::MuNuclearSplittingProcessXSecBias::AlongStepGetPhysicalInteractionLength ( const G4Track &  track, G4double  previousStepSize, G4double  currentMinimumStep, G4double &  proposedSafety, G4GPILSelection *  selection  )

virtual

Definition at line 288 of file MuNuclearSplittingProcessXSecBias.cxx.

    {
      if (xBiasMode==2)
      {
          // I am pretty sure this mode does not work because MuNuclear doesn't interact
          // AlongStep. Below often (always?) returns a negative number, which I have to
          // swizzle.
        G4double intLen;
        intLen = (1./eFactor) *
        pRegProcess->AlongStepGetPhysicalInteractionLength( track,previousStepSize,
                                                           currentMinimumStep,proposedSafety,selection);
        if (intLen<0) return currentMinimumStep;
        else return intLen;
      }
      else
        return DBL_MAX ;
    }

G4double gar::garg4::MuNuclearSplittingProcessXSecBias::GetTotalNumberOfInteractionLengthTraversed ( )

private

Definition at line 257 of file MuNuclearSplittingProcessXSecBias.cxx.

    {
      return theInitialNumberOfInteractionLength
      -G4VProcess::theNumberOfInteractionLengthLeft;
    }

G4VParticleChange * gar::garg4::MuNuclearSplittingProcessXSecBias::PostStepDoIt	(	const G4Track &	track,
		const G4Step &	step
	)

Definition at line 47 of file MuNuclearSplittingProcessXSecBias.cxx.

    {
      
      G4VParticleChange* pChange = pRegProcess->PostStepDoIt( track, step );
      pChange->SetVerboseLevel(0);
      pChange->SetSecondaryWeightByProcess(true);
      
      G4double secWeight = 1.0;
      if (fNSplit>=1) secWeight = 1.0/fNSplit;
      std::vector<G4Track*> secondaries; // Secondary store
      
      G4double nw=1.0;
        // Let's go bias the underlying xsection.
      if(xBiasMode==1)
      {
        G4double s = step.GetStepLength();
        G4double x = pRegProcess->GetCurrentInteractionLength();
        nw =  pChange->GetParentWeight();
          // need to undo the change made to the parent track in
          // most recent step (1/wc).
          // It seems to be possible to get here w.o. first having hit the
          // AlongStep Method that initialied wc. So protect against division y 0.
        if (wc < 1.0e-10) wc = 1.0;
        nw *= 1/wc * (1. - exp(-s/x))/(1 - exp(-eFactor*s/x)) ;
        G4cout << " MNSPXSB.PostStepDoit(): original wt = " <<  pChange->GetParentWeight() << " eF = " << eFactor << G4endl;
        G4cout << " MNSPXSB.PostStepDoit(): New weight = " << nw << " wc = " << wc << " s= " << s <<" x = " << x <<  G4endl;
        ((G4ParticleChange*)pChange)->ProposeParentWeight(nw) ;
      }
      else if (xBiasMode==2)
      {
        nw = pChange->GetParentWeight()*XBiasSecondaryWeight();
        ((G4ParticleChange*)pChange)->ProposeParentWeight(nw) ;
          //  for (G4int i = 0; i < pChange->GetNumberOfSecondaries(); i++)
          //pChange->GetSecondary(i)->SetWeight(nw);
        if(G4UniformRand()<XBiasSurvivalProbability())
        { // need to add the primary back in to balance the weight
          // the position of the extra track should be moved to the next volume, but I don't know how yet.
          // need to change number of secondary first before adding an additional track
          pChange->SetNumberOfSecondaries(pChange->GetNumberOfSecondaries()+1);
          G4ThreeVector mDirection = track.GetDynamicParticle()->GetMomentumDirection();
          G4DynamicParticle* aD = new G4DynamicParticle((track.GetDynamicParticle())->GetDefinition(),
                                                        G4LorentzVector(mDirection,
                                                                        (track.GetDynamicParticle())->GetKineticEnergy()+step.GetTotalEnergyDeposit()));
          G4Track* secTrk = new G4Track(aD,step.GetDeltaTime(),step.GetDeltaPosition());
          pChange->AddSecondary(secTrk);
          pChange->GetSecondary(pChange->GetNumberOfSecondaries()-1)->
          SetWeight(XBiasSurvivalProbability()*track.GetWeight());
        }
        else
        {
          nw = track.GetWeight();
          pChange->ProposeParentWeight(nw);
        }
      }
      
      secWeight = secWeight*nw;
      
        // Loop over PostStepDoIt method to generate multiple secondaries.
        // The point is for each value of ii up to Nsplit we want to re-toss
        // all secondaries for the muNucl process. Each toss gives back numSec
        // secondaries, which will be a different sample of species/momenta each time.
      G4VParticleChange* particleChange = new G4VParticleChange();
      
      for (unsigned int ii=0; ii<(unsigned int)fNSplit; ii++) {
        particleChange = pRegProcess->PostStepDoIt(track, step);
        if (!particleChange)
          throw std::runtime_error("MuNuclearSplittingProcessXSecBias::PostStepDoIt(): no particle change");
        G4int j(0);
        G4int numSec(particleChange->GetNumberOfSecondaries());
          // Don't change weight of last secondary. It's the just-added primary in this mode.
        for (j=0; j<numSec; j++)
        {
          G4Track* newSec = new G4Track(*(particleChange->GetSecondary(j)));
          G4String pdgstr = newSec->GetParticleDefinition()->GetParticleName();
            //G4double ke = newSec->GetKineticEnergy()/GeV;
          G4int pdg = newSec->GetParticleDefinition()->GetPDGEncoding();
          if (abs(pdg)==310 ||abs(pdg)==311 || abs(pdg)==3122 || abs(pdg)==2112)
          {
              //      std::cout << "MuNuclSplProc: Creating " << pdgstr << " of Kinetic Energy " << ke << " GeV. numSec is " << j << std::endl;
            
            if (pdg==G4KaonZeroShort::KaonZeroShort()->GetPDGEncoding()&&G4UniformRand()<0.50)
            {
              pdg = G4KaonZeroLong::KaonZeroLong()->GetPDGEncoding();
              pdgstr = G4KaonZeroLong::KaonZeroLong()->GetParticleName();
              G4LorentzVector pK0L(newSec->GetMomentum(),TMath::Sqrt(TMath::Power(G4KaonZeroLong::KaonZeroLong()->GetPDGMass(),2)+TMath::Power(newSec->GetMomentum().mag(),2)));
              G4DynamicParticle *newK0L = new G4DynamicParticle(G4KaonZeroLong::KaonZeroLong(),pK0L);
              
              G4Track* newSecK0L = new G4Track(newK0L,track.GetGlobalTime(),track.GetPosition());
              secondaries.push_back(newSecK0L);
            }
            else
            {
              secondaries.push_back(newSec);
            }
          }
        }
      }
      
      pChange->SetNumberOfSecondaries(secondaries.size());
      pChange->SetSecondaryWeightByProcess(true);
        //pChange->ProposeTrackStatus(fStopAndKill); // End of the story for this muon.
      
        //int numSecAdd = secondaries.size();
      std::vector<G4Track*>::iterator iter = secondaries.begin(); // Add all secondaries
      while (iter != secondaries.end()) {
        G4Track* myTrack = *iter;
        G4String pdgstr = myTrack->GetParticleDefinition()->GetParticleName();
        G4double ke = myTrack->GetKineticEnergy()/CLHEP::GeV;
        G4int pdg = myTrack->GetParticleDefinition()->GetPDGEncoding();
        if (abs(pdg)==130 ||abs(pdg)==310 ||abs(pdg)==311 || abs(pdg)==3122 || abs(pdg)==2112)
        {
          if (pdg!=2112)
            std::cout << "MuNuclSplXSB(): Adding " << pdgstr << " of Kinetic Energy and weight " << ke << " GeV and " << secWeight << "." << std::endl;
          
          myTrack->SetWeight(secWeight);
          pChange->AddSecondary(myTrack);
        }
        iter++;
      }
      return pChange;
      
    }

G4double gar::garg4::MuNuclearSplittingProcessXSecBias::PostStepGetPhysicalInteractionLength ( const G4Track &  track, G4double  previousStepSize, G4ForceCondition *  condition  )

virtual

Definition at line 264 of file MuNuclearSplittingProcessXSecBias.cxx.

    {
        // I believe this is the money line for the whole thing. By shrinking the
        // interaction length in MuNuclear's process, it's more likely that
        // this process will be chosen in the competition among processes.
      if(xBiasMode)
      {
        G4double psgPIL = (1./eFactor) * pRegProcess->PostStepGetPhysicalInteractionLength(
                                                                                           track,
                                                                                           previousStepSize,
                                                                                           condition );
        return psgPIL;
      }
      
      return pRegProcess->PostStepGetPhysicalInteractionLength(
                                                               track,
                                                               previousStepSize,
                                                               condition );
      
    }

virtual void gar::garg4::MuNuclearSplittingProcessXSecBias::ResetNumberOfInteractionLengthLeft ( )

inlineprotectedvirtual

Definition at line 59 of file MuNuclearSplittingProcessXSecBias.h.

      {
        G4VProcess::theNumberOfInteractionLengthLeft =  -std::log( G4UniformRand() );
        theInitialNumberOfInteractionLength = G4VProcess::theNumberOfInteractionLengthLeft;
      }

void gar::garg4::MuNuclearSplittingProcessXSecBias::SetIsActive ( G4bool  doIt )

inline

Definition at line 38 of file MuNuclearSplittingProcessXSecBias.h.

{fActive = doIt;};

void gar::garg4::MuNuclearSplittingProcessXSecBias::SetNSplit ( G4int  nTrx, G4int  xB = 0, G4double  xFac = 1  )

inline

Definition at line 37 of file MuNuclearSplittingProcessXSecBias.h.

{fNSplit = nTrx, eFactor = (G4double) xFac, xBiasMode = xB;};

G4double gar::garg4::MuNuclearSplittingProcessXSecBias::XBiasSecondaryWeight ( )

private

Definition at line 248 of file MuNuclearSplittingProcessXSecBias.cxx.

    {
      G4double result = 0;
      G4double nLTraversed = GetTotalNumberOfInteractionLengthTraversed();
      result = 1./eFactor*std::exp(-nLTraversed/eFactor*(1-1./eFactor));
      return result;
    }

G4double gar::garg4::MuNuclearSplittingProcessXSecBias::XBiasSurvivalProbability ( )

private

Definition at line 238 of file MuNuclearSplittingProcessXSecBias.cxx.

    {
      G4double result = 0;
      G4double nLTraversed = GetTotalNumberOfInteractionLengthTraversed();
      G4double biasedProbability = 1.-std::exp(-nLTraversed);
      G4double realProbability = 1-std::exp(-nLTraversed/eFactor);
      result = (biasedProbability-realProbability)/biasedProbability;
      return result;
    }

Member Data Documentation

G4double gar::garg4::MuNuclearSplittingProcessXSecBias::eFactor

private

Definition at line 71 of file MuNuclearSplittingProcessXSecBias.h.

G4bool gar::garg4::MuNuclearSplittingProcessXSecBias::fActive

private

Definition at line 69 of file MuNuclearSplittingProcessXSecBias.h.

G4int gar::garg4::MuNuclearSplittingProcessXSecBias::fNSplit

private

Definition at line 68 of file MuNuclearSplittingProcessXSecBias.h.

G4VParticleChange gar::garg4::MuNuclearSplittingProcessXSecBias::fParticleChange

private

Definition at line 73 of file MuNuclearSplittingProcessXSecBias.h.

G4double gar::garg4::MuNuclearSplittingProcessXSecBias::theInitialNumberOfInteractionLength

private

Definition at line 76 of file MuNuclearSplittingProcessXSecBias.h.

G4double gar::garg4::MuNuclearSplittingProcessXSecBias::wc

private

Definition at line 75 of file MuNuclearSplittingProcessXSecBias.h.

G4int gar::garg4::MuNuclearSplittingProcessXSecBias::xBiasMode

private

Definition at line 70 of file MuNuclearSplittingProcessXSecBias.h.

---

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

• garsoft/GArG4/MuNuclearSplittingProcessXSecBias.h
• garsoft/GArG4/MuNuclearSplittingProcessXSecBias.cxx