Public Member Functions | Static Public Member Functions | List of all members
LBNEImpWeight Class Reference

#include <LBNEImpWeight.hh>

Public Member Functions

 LBNEImpWeight ()
 
virtual ~LBNEImpWeight ()
 

Static Public Member Functions

static G4double CalculateImpWeight (const G4Track *aTrack)
 

Detailed Description

Definition at line 13 of file LBNEImpWeight.hh.

Constructor & Destructor Documentation

LBNEImpWeight::LBNEImpWeight ( )

Definition at line 18 of file LBNEImpWeight.cc.

19 {
20 }
LBNEImpWeight::~LBNEImpWeight ( )
virtual

Definition at line 22 of file LBNEImpWeight.cc.

23 {
24 }

Member Function Documentation

G4double LBNEImpWeight::CalculateImpWeight ( const G4Track *  aTrack)
static

Definition at line 25 of file LBNEImpWeight.cc.

26 {
27 
28  LBNERunManager* pRunManager =
29  dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
30 
31  LBNEPrimaryGeneratorAction* NPGA=(LBNEPrimaryGeneratorAction*)pRunManager->GetUserPrimaryGeneratorAction();
32  G4double impwt=NPGA->GetWeight();
33 
34  //Check if its a primary particle; if yes then weight either comes from external ntuple
35  // or its equal to 1.
36  if (aTrack->GetTrackID()==1) {
37  return impwt;
38  } else {
39  //not a primary then calculate weight:
40  G4int IMPWEIGHT_NBINS=40;
41  G4double psave;
42  G4ParticleDefinition * particleType = aTrack->GetDefinition();
43  G4double totE=aTrack->GetTotalEnergy()/CLHEP::GeV;
44 
45  //first find the importance weight of the parent
46  G4TrajectoryContainer* container =
47  pRunManager->GetCurrentEvent()->GetTrajectoryContainer();
48  if(container!=0) {
49  TrajectoryVector* vect = container->GetVector();
50  G4VTrajectory* tr;
51  G4int ii=0;
52  while (ii<G4int(vect->size())){
53  tr=(*vect)[ii];
54  LBNETrajectory* tr1=dynamic_cast<LBNETrajectory*>(tr);
55  if(tr1->GetTrackID()==aTrack->GetParentID()) impwt=tr1->GetNImpWt();
56  ii++;
57  }
58 // std::cerr << " LBNEImpWeight::CalculateImpWeight, ancestor weight for track ID "
59 // << aTrack->GetTrackID() << " is " << impwt << std::endl;
60  } else {
61  G4cout <<"LBNEImpWeight::Not a primary and no particle info"<<G4endl;
62  }
63  //Check if particle is neutrino
64  if (particleType==G4NeutrinoE::NeutrinoEDefinition()||
65  particleType==G4AntiNeutrinoE::AntiNeutrinoEDefinition()||
66  particleType==G4NeutrinoMu::NeutrinoMuDefinition()||
67  particleType==G4AntiNeutrinoMu::AntiNeutrinoMuDefinition()||
68  particleType==G4NeutrinoTau::NeutrinoTauDefinition()||
69  particleType==G4AntiNeutrinoTau::AntiNeutrinoTauDefinition())
70  {
71  return impwt;
72  }
73  //From gnumi/Impwght.F and beam/ffkeyImpWeight.F
74  G4int indx = (G4int)(totE+0.5); // Bins are 1 GeV
75  if (indx<1) indx = 1;
76  else if (indx>IMPWEIGHT_NBINS) indx = IMPWEIGHT_NBINS;
77 
78  G4double wEM = 0.01;
79  G4double wHad = (G4double)(indx)/30.0;
80  if (wHad>1.) wHad = 1.000;
81  if (wHad<0.001) wHad = 0.001;
82  //
83  // Get the survival probability from look up tables
84  //
85  if (particleType==G4Gamma::GammaDefinition()) psave = wEM;
86  else if (particleType==G4Electron::ElectronDefinition()) psave = wEM;
87  else if (particleType==G4Positron::PositronDefinition()) psave = wEM;
88  else if (particleType==G4MuonPlus::MuonPlusDefinition()) psave = 1.;
89  else if (particleType==G4MuonMinus::MuonMinusDefinition()) psave = 1.;
90  else if (particleType==G4PionZero::PionZeroDefinition()) psave = wHad;
91  else if (particleType==G4PionPlus::PionPlusDefinition()) psave = wHad;
92  else if (particleType==G4PionMinus::PionMinusDefinition()) psave = wHad;
93  else if (particleType==G4KaonZero::KaonZeroDefinition()) psave = wHad;
94  else if (particleType==G4KaonPlus::KaonPlusDefinition()) psave = wHad;
95  else if (particleType==G4KaonMinus::KaonMinusDefinition()) psave = wHad;
96  else if (particleType==G4Neutron::NeutronDefinition()) psave = wHad;
97  else if (particleType==G4Proton::ProtonDefinition()) psave = wHad;
98  else if (particleType==G4AntiProton::AntiProtonDefinition()) psave = wHad;
99  else if (particleType==G4KaonZeroShort::KaonZeroShortDefinition()) psave = wHad;
100  else psave = wHad;
101 
102  if (psave > 1. || impwt/psave>100.) psave = 1.0;
103  //don't let the weight be bigger than 100
104 
105 
106  if (G4UniformRand()<psave) { //G4UniformRand() returns random number between 0 and 1
107 // std::cerr << " LBNEImpWeight::CalculateImpWeight, final weight for track ID "
108 // << aTrack->GetTrackID() << " is " << impwt*1.0/psave <<
109 // " Track P " << aTrack->GetMomentum() << " E " << totE << " psave " << psave << " parent weight " << impwt << std::endl;
110  return impwt*1.0/psave;
111  } else {
112  return 0.;
113  }
114  }
115 }
G4int GetTrackID() const
virtual G4double GetNImpWt() const

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