PauliBlocker.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
8  University of Liverpool & STFC Rutherford Appleton Lab
9 
10  Joe Johnston (Univ of Pittsburgh) added code (Mar 18, 2016) to use
11  either local or relativistic Fermi Gas for Pauli blocking.
12 
13  Changes required to implement the GENIE Boosted Dark Matter module
14  were installed by Josh Berger (Univ. of Wisconsin)
15 */
16 //____________________________________________________________________________
17 
18 #include <TLorentzVector.h>
19 #include <TVector3.h>
20 
27 
41 
42 using namespace genie;
43 using namespace genie::constants;
44 
45 //___________________________________________________________________________
47 EventRecordVisitorI("genie::PauliBlocker")
48 {
49 
50 }
51 //___________________________________________________________________________
53 EventRecordVisitorI("genie::PauliBlocker", config)
54 {
55 
56 }
57 //___________________________________________________________________________
59 {
60 
61 }
62 //___________________________________________________________________________
64 {
65  // Return if the neutrino was not scatterred off a nuclear target
66  GHepParticle * nucltgt = evrec->TargetNucleus();
67  if (!nucltgt) {
68  LOG("PauliBlock", pINFO)
69  << "No nuclear target found - The Pauli Blocker exits";
70  return;
71  }
72 
73  // Handle only QEL for now...
74  // (can also be dark matter elastic)
75  Interaction * interaction = evrec->Summary();
76  const ProcessInfo & proc = interaction->ProcInfo();
77  if(!proc.IsQuasiElastic() && !proc.IsDarkMatterElastic()) {
78  LOG("PauliBlock", pINFO) << "Not a QEL event - The Pauli Blocker exits";
79  return;
80  }
81 
82  // Get the particle representing the initial hit nucleon
83  GHepParticle * hit = evrec->HitNucleon();
84  assert(hit);
85 
86  // Get the particle representing the recoiling final nucleon
87  GHepParticle * recoil = evrec->Particle(hit->FirstDaughter());
88  assert(recoil);
89  int nuc_pdgc = recoil->Pdg();
90 
91  const Target& tgt = interaction->InitState().Tgt();
92  double radius = hit->X4()->Vect().Mag();
93  double kf = this->GetFermiMomentum(tgt, nuc_pdgc, radius);
94 
95  LOG("PauliBlock", pINFO) << "KF = " << kf;
96 
97  // get the recoil momentum
98  double p = recoil->P4()->P(); // |p| for the recoil nucleon
99  LOG("PauliBlock", pINFO) << "Recoil nucleon |P| = " << p;
100 
101  // check for pauli blocking
102  bool is_blocked = (p < kf);
103 
104  // if it is blocked, set & thow an exception
105  if ( is_blocked ) {
106  LOG("PauliBlock", pNOTICE)
107  << " *** The generated event is Pauli-blocked ("
108  << "|p_{nucleon}| = " << p << " GeV < Fermi momentum = " << kf << " GeV) ***";
109 
110  evrec->EventFlags()->SetBitNumber(kPauliBlock, true);
112  exception.SetReason("Pauli-blocked event");
113 
114  // Include dark matter elastic
115  if(proc.IsQuasiElastic() || proc.IsDarkMatterElastic()) {
116  // nuclear suppression taken into account at the QEL cross
117  // section - should attempt to regenerate the event as QEL
118  exception.SwitchOnStepBack();
119  exception.SetReturnStep(0);
120  } else {
121  // end this event generation thread and start again at the
122  // interaction selection step
123  // - this is irrelevant for the time being as we only handle QEL-
124  exception.SwitchOnFastForward();
125  }
126  throw exception;
127  }
128 }
129 //___________________________________________________________________________
131 {
132  Algorithm::Configure(config);
133  this->LoadModelType();
134 }
135 //___________________________________________________________________________
136 void PauliBlocker::Configure(string param_set)
137 {
138  Algorithm::Configure(param_set);
139  this->LoadModelType();
140 }
141 //___________________________________________________________________________
144  const Registry * gc = confp->GlobalParameterList();
145 
146  // Create a nuclear model object to check the model type
147  RgKey nuclkey = "NuclearModel";
148  RgAlg nuclalg = gc->GetAlg(nuclkey);
149  AlgFactory * algf = AlgFactory::Instance();
150  const NuclearModelI* nuclModel =
151  dynamic_cast<const NuclearModelI*>(
152  algf->GetAlgorithm(nuclalg.name,nuclalg.config));
153  // Check if the model is a local Fermi gas
154  fLFG = (nuclModel && nuclModel->ModelType(Target()) == kNucmLocalFermiGas);
155 
156  if ( !fLFG ) {
157  // get the Fermi momentum table for relativistic Fermi gas
158  GetParam( "FermiMomentumTable", fKFTableName ) ;
159 
160  fKFTable = 0;
161 
163  fKFTable = kftp->GetTable(fKFTableName);
164  assert(fKFTable);
165  }
166 }
167 //___________________________________________________________________________
168 double PauliBlocker::GetFermiMomentum(const Target& tgt, int pdg_Nf,
169  double radius) const
170 {
171  // Pauli blocking should only be applied for nucleons
172  assert( pdg::IsProton(pdg_Nf) || pdg::IsNeutron(pdg_Nf) );
173  double kF = 0.;
174  if ( fLFG ) {
175  int A = tgt.A();
176  bool is_p = pdg::IsProton( pdg_Nf );
177  int numNuc = (is_p) ? tgt.Z() : tgt.N();
179  kF = TMath::Power(3 * kPi2 * numNuc *
180  genie::utils::nuclear::Density(radius, A), 1.0/3.0) * hbarc;
181  }
182  else {
183  kF = fKFTable->FindClosestKF(tgt.Pdg(), pdg_Nf);
184  }
185 
186  return kF;
187 }
Basic constants.
virtual GHepParticle * Particle(int position) const
Definition: GHepRecord.cxx:148
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
virtual Interaction * Summary(void) const
Definition: GHepRecord.cxx:135
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:79
Defines the EventRecordVisitorI interface. Concrete implementations of this interface use the &#39;Visito...
double Density(double r, int A, double ring=0.)
int FirstDaughter(void) const
Definition: GHepParticle.h:69
static const double fermi
Definition: Units.h:63
int A(void) const
Definition: Target.h:71
bool IsQuasiElastic(void) const
Definition: ProcessInfo.cxx:67
A singleton class holding all configuration registries built while parsing all loaded XML configurati...
Definition: AlgConfigPool.h:41
const char * p
Definition: xmltok.h:285
int Pdg(void) const
Definition: Target.h:72
static FermiMomentumTablePool * Instance(void)
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:42
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
Definition: config.py:1
void Configure(const Registry &config)
virtual NuclearModel_t ModelType(const Target &) const =0
int Pdg(void) const
Definition: GHepParticle.h:64
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:304
Summary information for an interaction.
Definition: Interaction.h:56
double GetFermiMomentum(const Target &tgt, int pdg_Nf, double radius=0.0) const
Get the Fermi momentum needed to check Pauli blocking.
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:299
Singleton class to load & serve tables of Fermi momentum constants.
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
const FermiMomentumTable * GetTable(string name)
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:44
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:41
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
const Algorithm * GetAlgorithm(const AlgId &algid)
Definition: AlgFactory.cxx:86
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:330
int Z(void) const
Definition: Target.h:69
bool IsDarkMatterElastic(void) const
Definition: ProcessInfo.cxx:72
Double_t radius
#define pINFO
Definition: Messenger.h:63
Registry * GlobalParameterList(void) const
static const double kLightSpeed
Definition: Constants.h:32
void LoadModelType(void)
static const double A
Definition: Units.h:82
int N(void) const
Definition: Target.h:70
virtual TBits * EventFlags(void) const
Definition: GHepRecord.h:118
static AlgFactory * Instance()
Definition: AlgFactory.cxx:75
string RgKey
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:80
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
Definition: structs.h:12
virtual GHepParticle * HitNucleon(void) const
Definition: GHepRecord.cxx:350
const FermiMomentumTable * fKFTable
Definition: PauliBlocker.h:59
assert(nhit_max >=nhit_nbins)
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
void ProcessEventRecord(GHepRecord *event_rec) const
double FindClosestKF(int target_pdgc, int nucleon_pdgc) const
#define pNOTICE
Definition: Messenger.h:62
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:67
The GENIE Algorithm Factory.
Definition: AlgFactory.h:40
static const double kPlankConstant
Definition: Constants.h:33
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:46
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:40
static const double kPi2
Definition: Constants.h:39
string config
RgAlg GetAlg(RgKey key) const
Definition: Registry.cxx:503
static AlgConfigPool * Instance()