MuNuclearSplittingProcessXSecBias.cxx

Go to the documentation of this file.

//
// This(these) is (are) a variance reduction technique(s).
//
//
// Secondary biasing code adapted from EventBiasing from Jane Tinslay, SLAC. 2008.
// (Who's now at Cisco, I think.) XSBiasing code from http://geant4.cern.ch/collaboration/workshops/workshop2002/slides/docs/flei/biasing.pdf. (And NOT from http://www.lcsim.org/software/geant4/doxygen/html/classG4HadronicProcess.html, where I *believe* the point is to bias primaries and secondaries from sub-classes of Cross-Sections in an already chosen
// broader class of XSs. This doesn't help us. We need more occurrences (albeit with small
// weights) of a specific process.
//
//
// The point is to create Nsplit secondaries with each muNucl process, so
// that we may, e.g., create many, otherwise-rare, K0s which
// may charge exchange into pernicious K+s which then mimic proton dk. And,
// additionally now to make the XSection for the muNucl process itself higher,
// keeping track of appropriate weights everywhere.
//
// echurch@fnal.gov, May, 2011.
//


#include "GArG4/MuNuclearSplittingProcessXSecBias.h"

#include "Geant4/G4WrapperProcess.hh"
#include "Geant4/G4VParticleChange.hh"
#include "Geant4/G4SDManager.hh"
#include "Geant4/G4RunManager.hh"
#include "Geant4/G4Event.hh"
#include "Geant4/G4HCofThisEvent.hh"
#include "Geant4/G4Track.hh"
#include "Geant4/G4Step.hh"
#include "Geant4/G4TrackStatus.hh"
#include "Geant4/G4ParticleDefinition.hh"
#include "Geant4/G4ParticleTypes.hh"
#include "Geant4/Randomize.hh"
#include "Geant4/G4KaonZeroLong.hh"
#include "Geant4/G4ios.hh"

#include "cetlib_except/exception.h"

#include <TMath.h>

namespace gar {
  namespace garg4 {
    
    
    
    G4VParticleChange* MuNuclearSplittingProcessXSecBias::PostStepDoIt(const G4Track& track, const G4Step& step)
    {
      
      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;
      
    }
    
    G4VParticleChange* MuNuclearSplittingProcessXSecBias::AlongStepDoIt(
                                                                        const G4Track& track,
                                                                        const G4Step& step
                                                                        )
    {
      
      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 MuNuclearSplittingProcessXSecBias::XBiasSurvivalProbability()
    {
      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;
    }
    
    G4double MuNuclearSplittingProcessXSecBias::XBiasSecondaryWeight()
    {
      G4double result = 0;
      G4double nLTraversed = GetTotalNumberOfInteractionLengthTraversed();
      result = 1./eFactor*std::exp(-nLTraversed/eFactor*(1-1./eFactor));
      return result;
    }
    
    
    G4double MuNuclearSplittingProcessXSecBias::GetTotalNumberOfInteractionLengthTraversed()
    {
      return theInitialNumberOfInteractionLength
      -G4VProcess::theNumberOfInteractionLengthLeft;
    }
    
    
    G4double MuNuclearSplittingProcessXSecBias::PostStepGetPhysicalInteractionLength( const G4Track& track,
                                                                                     G4double   previousStepSize,
                                                                                     G4ForceCondition* condition )
    {
        // 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 );
      
    }
    
    
    G4double MuNuclearSplittingProcessXSecBias::AlongStepGetPhysicalInteractionLength( const G4Track& track,
                                                                                      G4double previousStepSize,
                                                                                      G4double currentMinimumStep,
                                                                                      G4double& proposedSafety,
                                                                                      G4GPILSelection* selection )
    {
      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 ;
    }
    
    
  } // end namespace
} // gar