larg4::MuNuclearSplittingProcessXSecBias Class Reference

#include <MuNuclearSplittingProcessXSecBias.h>

Inheritance diagram for larg4::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 24 of file MuNuclearSplittingProcessXSecBias.h.

Constructor & Destructor Documentation

larg4::MuNuclearSplittingProcessXSecBias::MuNuclearSplittingProcessXSecBias ( )

inline

Definition at line 27 of file MuNuclearSplittingProcessXSecBias.h.

{};

larg4::MuNuclearSplittingProcessXSecBias::~MuNuclearSplittingProcessXSecBias ( )

inline

Definition at line 28 of file MuNuclearSplittingProcessXSecBias.h.

{};

Member Function Documentation

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

Definition at line 168 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 larg4::MuNuclearSplittingProcessXSecBias::AlongStepGetPhysicalInteractionLength ( const G4Track &  track, G4double  previousStepSize, G4double  currentMinimumStep, G4double &  proposedSafety, G4GPILSelection *  selection  )

virtual

Definition at line 286 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 larg4::MuNuclearSplittingProcessXSecBias::GetTotalNumberOfInteractionLengthTraversed ( )

private

Definition at line 255 of file MuNuclearSplittingProcessXSecBias.cxx.

  {
    return theInitialNumberOfInteractionLength
      -G4VProcess::theNumberOfInteractionLengthLeft;
  }

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

Definition at line 45 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 larg4::MuNuclearSplittingProcessXSecBias::PostStepGetPhysicalInteractionLength ( const G4Track &  track, G4double  previousStepSize, G4ForceCondition *  condition  )

virtual

Definition at line 262 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 larg4::MuNuclearSplittingProcessXSecBias::ResetNumberOfInteractionLengthLeft ( )

inlineprotectedvirtual

Definition at line 52 of file MuNuclearSplittingProcessXSecBias.h.

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

void larg4::MuNuclearSplittingProcessXSecBias::SetIsActive ( G4bool  doIt )

inline

Definition at line 31 of file MuNuclearSplittingProcessXSecBias.h.

{fActive = doIt;};

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

inline

Definition at line 30 of file MuNuclearSplittingProcessXSecBias.h.

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

G4double larg4::MuNuclearSplittingProcessXSecBias::XBiasSecondaryWeight ( )

private

Definition at line 246 of file MuNuclearSplittingProcessXSecBias.cxx.

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

G4double larg4::MuNuclearSplittingProcessXSecBias::XBiasSurvivalProbability ( )

private

Definition at line 236 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 larg4::MuNuclearSplittingProcessXSecBias::eFactor

private

Definition at line 64 of file MuNuclearSplittingProcessXSecBias.h.

G4bool larg4::MuNuclearSplittingProcessXSecBias::fActive

private

Definition at line 62 of file MuNuclearSplittingProcessXSecBias.h.

G4int larg4::MuNuclearSplittingProcessXSecBias::fNSplit

private

Definition at line 61 of file MuNuclearSplittingProcessXSecBias.h.

G4VParticleChange larg4::MuNuclearSplittingProcessXSecBias::fParticleChange

private

Definition at line 66 of file MuNuclearSplittingProcessXSecBias.h.

G4double larg4::MuNuclearSplittingProcessXSecBias::theInitialNumberOfInteractionLength

private

Definition at line 69 of file MuNuclearSplittingProcessXSecBias.h.

G4double larg4::MuNuclearSplittingProcessXSecBias::wc

private

Definition at line 68 of file MuNuclearSplittingProcessXSecBias.h.

G4int larg4::MuNuclearSplittingProcessXSecBias::xBiasMode

private

Definition at line 63 of file MuNuclearSplittingProcessXSecBias.h.

---

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

• larsim/larsim/LegacyLArG4/MuNuclearSplittingProcessXSecBias.h
• larsim/larsim/LegacyLArG4/MuNuclearSplittingProcessXSecBias.cxx