LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
ReweightAna_module.cc
Go to the documentation of this file.
1 //
8 #include <vector>
9 #include <cmath>
10 #include "TH1.h"
11 #include "TH2.h"
12 #include "TH3.h"
13 #include "TDatabasePDG.h"
14 #include "TSystem.h"
19 
22 #include "fhiclcpp/ParameterSet.h"
27 #include "art_root_io/TFileService.h"
28 #include "art_root_io/TFileDirectory.h"
30 
32 
34 
36 
37 
38 namespace simb { class MCTruth; }
39 namespace simb { class GTruth; }
40 namespace rwgt {class NuReweight; }
41 
42 class TH1F;
43 class TH2F;
44 
45 namespace rwgt {
46 
47  //class GENIEReweight;
49 
50  class ReweightAna : public art::EDAnalyzer {
51 
52  public:
53  explicit ReweightAna(fhicl::ParameterSet const& pset);
54  virtual ~ReweightAna();
55 
56  void analyze(art::Event const& evt);
57  void beginSubRun(art::SubRun const& sr);
58  void beginJob();
59  void endJob();
60  void endSubRun(art::SubRun const& sr);
61  void reconfigure(const fhicl::ParameterSet& p);
62 
63  private:
64  void LoadMCInfo(art::Event const& evt);
65 
66  private:
69  TH1F* fWgtQE[3];
70  TH1F* fWgtRES[3];
71  TH1F* fWgtDIS[3];
72  rwgt::NuReweight* fGrwgt[3];
73 
74  std::string fMCTruthModuleLabel;
75  std::string fPotLabel;
76 
77  };
78 }
79 
80 namespace rwgt{
81 
82  static int cntEvent = 0; //event number
83 
84  //......................................................................
85  ReweightAna::ReweightAna(fhicl::ParameterSet const& p)
86  : EDAnalyzer(p)
87  {
88  this->reconfigure(p);
89  }
90 
91 
92  //......................................................................
94 
95  //......................................................................
97  {
99 
100  mf::LogVerbatim("ReweightAna") << "make histograms" ;
101 
102  fEnergyNeutrino = tfs->make<TH1F>("fEnergyneutrino", "Total number of events", 50, 0., 25);
103  fNeventsSubrun = tfs->make<TH1F>("fNeventsSubrun", "Total number of events", 1, 0., 1.);
104 
105  char name[300];
106  for(int i = 0; i < 3; i++) {
107  sprintf(name, "fWgtQE_%dsigma", i+1);
108  fWgtQE[i] = tfs->make<TH1F>(name, "Evt Wgts", 100, 0., 2.0);
109  sprintf(name, "fWgtRES_%dsigma", i+1);
110  fWgtRES[i] = tfs->make<TH1F>(name, "Evt Wgts", 100, 0., 2.0);
111  sprintf(name, "fWgtDIS_%dsigma", i+1);
112  fWgtDIS[i] = tfs->make<TH1F>(name, "Evt Wgts", 100, 0., 2.0);
113 
114  double sigma = (double)(i+1);
115  fGrwgt[i] = new rwgt::NuReweight();
127 
128  fGrwgt[i]->Configure();
129  }
130  }
131 
132  //......................................................................
134 
135  }
136 
137  //......................................................................
139  {
140 
141  // // Pull the MC generator information out of the event
142  mf::LogVerbatim("ReweightAna") << "Start analyze" ;
144  evt.getByLabel(fMCTruthModuleLabel, mclist);
145  if (mclist->empty()) {
146  mf::LogWarning("ReweightAna") << "Error retrieving MCTruth list" ;
147  return;
148  }
149 
151  evt.getByLabel(fMCTruthModuleLabel, gtlist);
152  if (gtlist->empty()) {
153  mf::LogWarning("ReweightAna") << "Error retrieving GTruth list" ;
154  return;
155  }
156 
157  MF_LOG_DEBUG("ReweightAna")<<"MC List sizes:" << mclist->size() << " " << gtlist->size() << "\n";
158 
159  // // Loop over neutrino interactions
160  for(size_t i_intx = 0; i_intx < mclist->size(); ++i_intx){
161  MF_LOG_DEBUG("ReweightAna") << "start loop";
162 
163  // // Link to the MCNeutrino class.
164  // // The class contains information not only about
165  // // the incoming neutrino, but about the products of the decay
166  simb::MCTruth const& truth = mclist->at(i_intx);
167  simb::GTruth const& gtruth = gtlist->at(i_intx);
168  simb::MCNeutrino const& mc_neutrino = truth.GetNeutrino();
169 
170  fEnergyNeutrino->Fill(mc_neutrino.Nu().E());
171  for(int i = 0; i < 3; i++) {
172  double wgt = fGrwgt[i]->CalcWeight(truth, gtruth);
173  //double wgt = 1.;
174  if(mc_neutrino.Mode()==0 && mc_neutrino.CCNC()==0) {
175  fWgtQE[i]->Fill(wgt);
176  }
177  else if(mc_neutrino.Mode()==1 && mc_neutrino.CCNC()==0) {
178  fWgtRES[i]->Fill(wgt);
179  }
180  else if(mc_neutrino.Mode()==2 && mc_neutrino.CCNC()==0) {
181  fWgtDIS[i]->Fill(wgt);
182  }
183  }
184 
185  MF_LOG_DEBUG("ReweightAna") << "end loop" ;
186  }//end loop over interactions
187 
188 
189  return;
190  }
191 
192 
193  //......................................................................
195  {
196 
197 
198  }
199 
200 
201  //......................................................................
203  {
204  fMCTruthModuleLabel = p.get< std::string>("MCTruthModuleLabel");
205  }
206 
207  //......................................................................
209  {
210  fNeventsSubrun->Fill(sr.subRun(), cntEvent);
211  cntEvent = 0;
212 
213  }
214 
215  //......................................................................
217  {
218 
219  }
220 
221 
223 }
224 //Module Definition
225 
226 namespace rwgt{
227 
229 
230 }
231 
232 
double E(const int i=0) const
Definition: MCParticle.h:234
void AddReweightValue(ReweightLabel_t rLabel, double value)
Change a reweight parameter. If it hasn&#39;t been added yet add it.
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
int CCNC() const
Definition: MCNeutrino.h:148
tweak the 2pi non-RES bkg in the RES region, for v+n CC
TH1F * fNeventsSubrun
Total number of events per subrun.
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:77
void Configure()
Reconfigure the weight calculators.
void endSubRun(art::SubRun const &sr)
tweak the 2pi non-RES bkg in the RES region, for v+p CC
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
static int cntEvent
tweak the 2pi non-RES bkg in the RES region, for v+p NC
object containing MC flux information
void beginSubRun(art::SubRun const &sr)
void analyze(art::Event const &evt)
std::string fMCTruthModuleLabel
label for module producing mc truth information
tweak the 2pi non-RES bkg in the RES region, for v+n NC
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
tweak Ma NCRES, affects d2sigma(NCRES)/dWdQ2 both in shape and normalization
void beginJob()
Definition: Breakpoints.cc:14
tweak the 1pi non-RES bkg in the RES region, for v+n CC
T get(std::string const &key) const
Definition: ParameterSet.h:314
std::string fPotLabel
Module that produced the POTSum object.
tweak the 1pi non-RES bkg in the RES region, for v+p NC
ART objects.
TH1F * fEnergyNeutrino
Total number of events.
void reconfigure(const fhicl::ParameterSet &p)
SubRunNumber_t subRun() const
Definition: SubRun.cc:37
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
#define MF_LOG_DEBUG(id)
void LoadMCInfo(art::Event const &evt)
tweak Ma CCRES, affects d2sigma(CCRES)/dWdQ2 both in shape and normalization
tweak Ma CCQE, affects dsigma(CCQE)/dQ2 both in shape and normalization
rwgt::NuReweight * fGrwgt[3]
X-sec weight calculator.
TCEvent evt
Definition: DataStructs.cxx:8
Event generator information.
Definition: MCTruth.h:32
Event generator information.
Definition: MCNeutrino.h:18
int Mode() const
Definition: MCNeutrino.h:149
double CalcWeight(const simb::MCTruth &truth, const simb::GTruth &gtruth) const
Definition: NuReweight.cxx:112
A module to check the results from the Monte Carlo generator.
tweak the 1pi non-RES bkg in the RES region, for v+p CC
tweak the 1pi non-RES bkg in the RES region, for v+n NC
Wrapper for reweightings neutrino interactions within the ART framework.