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 167 of file MuNuclearSplittingProcessXSecBias.cxx.

References eFactor, fParticleChange, w, x, xBiasMode, XBiasSecondaryWeight(), and XBiasSurvivalProbability().

Referenced by SetIsActive().

  {

    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 273 of file MuNuclearSplittingProcessXSecBias.cxx.

References eFactor, and xBiasMode.

Referenced by SetIsActive().

  {
    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 251 of file MuNuclearSplittingProcessXSecBias.cxx.

References theInitialNumberOfInteractionLength.

Referenced by XBiasSecondaryWeight(), and XBiasSurvivalProbability().

  {
    return theInitialNumberOfInteractionLength - G4VProcess::theNumberOfInteractionLengthLeft;
  }

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

Definition at line 42 of file MuNuclearSplittingProcessXSecBias.cxx.

References util::abs(), e, eFactor, fNSplit, x, xBiasMode, XBiasSecondaryWeight(), and XBiasSurvivalProbability().

Referenced by SetIsActive().

  {

    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 256 of file MuNuclearSplittingProcessXSecBias.cxx.

References eFactor, and xBiasMode.

Referenced by SetIsActive().

  {
    // 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 51 of file MuNuclearSplittingProcessXSecBias.h.

References theInitialNumberOfInteractionLength.

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

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

inline

Definition at line 34 of file MuNuclearSplittingProcessXSecBias.h.

References AlongStepDoIt(), AlongStepGetPhysicalInteractionLength(), fActive, PostStepDoIt(), PostStepGetPhysicalInteractionLength(), and track.

Referenced by larg4::ModularPhysicsList::ConstructProcess().

{ fActive = doIt; };

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

inline

Definition at line 30 of file MuNuclearSplittingProcessXSecBias.h.

References eFactor, fNSplit, and xBiasMode.

Referenced by larg4::ModularPhysicsList::ConstructProcess().

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

G4double larg4::MuNuclearSplittingProcessXSecBias::XBiasSecondaryWeight ( )

private

Definition at line 243 of file MuNuclearSplittingProcessXSecBias.cxx.

References eFactor, and GetTotalNumberOfInteractionLengthTraversed().

Referenced by AlongStepDoIt(), and PostStepDoIt().

  {
    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 233 of file MuNuclearSplittingProcessXSecBias.cxx.

References eFactor, and GetTotalNumberOfInteractionLengthTraversed().

Referenced by AlongStepDoIt(), and PostStepDoIt().

  {
    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 62 of file MuNuclearSplittingProcessXSecBias.h.

Referenced by AlongStepDoIt(), AlongStepGetPhysicalInteractionLength(), PostStepDoIt(), PostStepGetPhysicalInteractionLength(), SetNSplit(), XBiasSecondaryWeight(), and XBiasSurvivalProbability().

G4bool larg4::MuNuclearSplittingProcessXSecBias::fActive

private

Definition at line 60 of file MuNuclearSplittingProcessXSecBias.h.

Referenced by SetIsActive().

G4int larg4::MuNuclearSplittingProcessXSecBias::fNSplit

private

Definition at line 59 of file MuNuclearSplittingProcessXSecBias.h.

Referenced by PostStepDoIt(), and SetNSplit().

G4VParticleChange larg4::MuNuclearSplittingProcessXSecBias::fParticleChange

private

Definition at line 64 of file MuNuclearSplittingProcessXSecBias.h.

Referenced by AlongStepDoIt().

G4double larg4::MuNuclearSplittingProcessXSecBias::theInitialNumberOfInteractionLength

private

Definition at line 67 of file MuNuclearSplittingProcessXSecBias.h.

Referenced by GetTotalNumberOfInteractionLengthTraversed(), and ResetNumberOfInteractionLengthLeft().

G4double larg4::MuNuclearSplittingProcessXSecBias::wc

private

Definition at line 66 of file MuNuclearSplittingProcessXSecBias.h.

G4int larg4::MuNuclearSplittingProcessXSecBias::xBiasMode

private

Definition at line 61 of file MuNuclearSplittingProcessXSecBias.h.

Referenced by AlongStepDoIt(), AlongStepGetPhysicalInteractionLength(), PostStepDoIt(), PostStepGetPhysicalInteractionLength(), and SetNSplit().

---

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

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