LArSoft  v06_85_00
Liquid Argon Software toolkit - http://larsoft.org/
MuNuclearSplittingProcess.cxx
Go to the documentation of this file.
1 //
2 // Code adapted from EventBiasing from Jane Tinslay, SLAC. 2008. Now at Cisco,
3 // I think.
4 //
5 // This is a variance reduction technique.
6 //
7 // The point is to create Nsplit secondaries with each muNucl process, so
8 // that we may, e.g., create many, otherwise-rare, K0s which
9 // may charge exchange into pernicious K+s which then mimic proton dk.
10 //
11 // echurch@fnal.gov, May, 2011.
12 //
13 
14 
16 
17 #include <stdexcept> // std::runtime_error
18 
19 #include "Geant4/G4WrapperProcess.hh"
20 #include "Geant4/G4VParticleChange.hh"
21 
22 #include "Geant4/G4SDManager.hh"
23 #include "Geant4/G4RunManager.hh"
24 #include "Geant4/G4Event.hh"
25 #include "Geant4/G4HCofThisEvent.hh"
26 #include "Geant4/G4Track.hh"
27 #include "Geant4/G4Step.hh"
28 #include "Geant4/G4TrackStatus.hh"
29 #include "Geant4/G4ParticleDefinition.hh"
30 #include "Geant4/G4ParticleTypes.hh"
31 #include "Geant4/Randomize.hh"
32 #include "Geant4/G4KaonZeroLong.hh"
33 #include "Geant4/G4ios.hh"
34 
35 #include <TMath.h>
36 
37 namespace larg4 {
38 
39 G4VParticleChange* MuNuclearSplittingProcess::PostStepDoIt(const G4Track& track, const G4Step& step)
40 {
41 
42  G4VParticleChange* particleChange = new G4VParticleChange();
43  G4double weight = track.GetWeight()/fNSplit;
44  std::vector<G4Track*> secondaries; // Secondary store
45 
46  // Loop over PostStepDoIt method to generate multiple secondaries.
47  // The point is for each value of ii up to Nsplit we want to re-toss
48  // all secondaries for the muNucl process. Each toss gives back numSec
49  // secondaries, which will be a different sample of species/momenta each time.
50  for (unsigned int ii=0; ii<(unsigned int)fNSplit; ii++) {
51  particleChange = pRegProcess->PostStepDoIt(track, step);
52  if (!particleChange)
53  throw std::runtime_error("MuNuclearSplittingProcess::PostStepDoIt(): no particle change");
54  G4int j(0);
55  G4int numSec(particleChange->GetNumberOfSecondaries());
56  for (j=0; j<numSec; j++)
57  {
58  G4Track* newSec = new G4Track(*(particleChange->GetSecondary(j)));
59  G4String pdgstr = newSec->GetParticleDefinition()->GetParticleName();
60  G4double ke = newSec->GetKineticEnergy()/CLHEP::GeV;
61  G4int pdg = newSec->GetParticleDefinition()->GetPDGEncoding();
62  if (abs(pdg)==310 ||abs(pdg)==311 || abs(pdg)==3122 || abs(pdg)==2112)
63  {
64  // std::cout << "MuNuclSplProc: Creating " << pdgstr << " of Kinetic Energy " << ke << " GeV. numSec is " << j << std::endl;
65 
66  if (pdg==G4KaonZeroShort::KaonZeroShort()->GetPDGEncoding()&&G4UniformRand()<0.50)
67  {
68  pdg = G4KaonZeroLong::KaonZeroLong()->GetPDGEncoding();
69  pdgstr = G4KaonZeroLong::KaonZeroLong()->GetParticleName();
70  G4LorentzVector pK0L(newSec->GetMomentum(),TMath::Sqrt(TMath::Power(G4KaonZeroLong::KaonZeroLong()->GetPDGMass(),2)+TMath::Power(newSec->GetMomentum().mag(),2)));
71  G4DynamicParticle *newK0L = new G4DynamicParticle(G4KaonZeroLong::KaonZeroLong(),pK0L);
72 
73  G4Track* newSecK0L = new G4Track(newK0L,track.GetGlobalTime(),track.GetPosition());
74  secondaries.push_back(newSecK0L);
75  }
76  else
77  {
78  secondaries.push_back(newSec);
79  }
80 
81  if (abs(pdg)==130 ||abs(pdg)==310 ||abs(pdg)==311 || abs(pdg)==3122)
82  {
83  // WrappedmuNuclear always produces K0s if it produces a K0 at all.
84  // Let's make half of these K0Ls.
85  std::cout << "MuNuclSplProc: Creating " << pdgstr << " of Kinetic Energy " << ke << " GeV. numSec is " << j << std::endl;
86  }
87 
88  }
89  }
90  }
91 
92  particleChange->SetNumberOfSecondaries(secondaries.size());
93  particleChange->SetSecondaryWeightByProcess(true);
94  //int numSecAdd = secondaries.size();
95  std::vector<G4Track*>::iterator iter = secondaries.begin(); // Add all secondaries
96  while (iter != secondaries.end()) {
97  G4Track* myTrack = *iter;
98  G4String pdgstr = myTrack->GetParticleDefinition()->GetParticleName();
99  //G4double ke = myTrack->GetKineticEnergy()/GeV;
100  G4int pdg = myTrack->GetParticleDefinition()->GetPDGEncoding();
101  if (abs(pdg)==130 ||abs(pdg)==310 ||abs(pdg)==311 || abs(pdg)==3122 || abs(pdg)==2112)
102  // std::cout << "MuNuclSplProc: Adding " << pdgstr << " of Kinetic Energy " << ke << " GeV." << std::endl;
103  // Will ask for PdgCode and only Add K0s.
104  myTrack->SetWeight(weight);
105  particleChange->AddSecondary(myTrack);
106  iter++;
107  }
108  return particleChange;
109 }
110 
111 } // end namespace
intermediate_table::iterator iterator
Geant4 interface.
Check of Geant4 to run the LArSoft detector simulation.
double weight
Definition: plottest35.C:25
Float_t track
Definition: plot.C:34
G4VParticleChange * PostStepDoIt(const G4Track &track, const G4Step &step)