Functions
protoana::ProtoDUNESelectionUtils Namespace Reference

Functions

TH1 * FillMCBackgroundHistogram_Pions (std::string filename, std::string treename, std::vector< double > recoBins, std::string channel, std::string topo, int toponum, double endZ_cut, double minval=0.0, double maxval=100000.0, bool doNegativeReco=false, int doSyst=0, std::string systName="", double weight=1.0, std::pair< double, double > PiMuScale={1., 1.})
 
TH1 * FillMCSignalHistogram_Pions (std::string filename, std::string treename, std::vector< double > recoBins, std::string channel, std::string topo, int toponum, double minval, double maxval, double endZ_cut, bool doNegativeReco=false, int doSyst=0, std::string systName="", double weight=1.0, std::pair< double, double > PiMuScale={1., 1.})
 
TH1 * FillMCIncidentHistogram_Pions (std::string filename, std::string treename, std::vector< double > recoBins, std::string topo, int toponum, double reco_beam_endZ_cut, double minval=0., double maxval=100000., bool doNegativeReco=false, int doSyst=0, std::string systName="", double weight=1.0, std::pair< double, double > PiMuScale={1., 1.})
 
TH1 * FillMCTruthSignalHistogram_Pions (std::string filename, std::string treename, std::vector< double > truthBins, std::string channel, double weight=1.0)
 
TH1 * FillDataHistogram_Pions (std::string filename, std::string treename, std::vector< double > recoBins, std::string channel, double reco_beam_endZ_cut, bool doNegativeReco=false, bool IsIncidentHisto=false)
 
TH1 * FillMCSidebandHistogram_Pions (std::string filename, std::string treename, std::string channel, std::string topo, int toponum, double endZ_cut, std::vector< double > binning, double weight=1.)
 
TH1 * FillDataSidebandHistogram_Pions (std::string filename, std::string treename, std::string channel, double endZ_cut, std::vector< double > binning, double weight=1.)
 
TH1 * FillMCFlux_Pions (std::string filename, std::string treename, std::vector< double > Bins, int mode=1, double weight=1.0)
 
int GetNTriggers_Pions (std::string filename, std::string treename, bool IsMC=true)
 
std::pair< TH1 *, TH1 * > GetMCIncidentEfficiency (std::string fileName, std::string treeName, std::vector< double > bins, double reco_beam_endZ_cut, int doSyst=0, double weight=1.)
 
std::pair< TH1 *, TH1 * > GetMCInteractingEfficiency (std::string fileName, std::string treeName, std::vector< double > bins, std::string channel, std::string topo, int toponum, double endZ_cut, int doSyst=0, double weight=1.)
 
int FindBin (double value, std::vector< double > binning)
 

Function Documentation

TH1 * protoana::ProtoDUNESelectionUtils::FillDataHistogram_Pions ( std::string  filename,
std::string  treename,
std::vector< double >  recoBins,
std::string  channel,
double  reco_beam_endZ_cut,
bool  doNegativeReco = false,
bool  IsIncidentHisto = false 
)

Definition at line 617 of file ProtoDUNESelectionUtils.cxx.

620  {
621  //********************************************************************
622 
623  TFile *file = new TFile(filename.c_str(), "READ");
624  TTree *defaultTree = (TTree*)file->Get(treename.c_str());
625 
626  Int_t reco_beam_type; // 13 -> track-like, 11 -> shower-like
627  Int_t reco_beam_nTrackDaughters, reco_beam_nShowerDaughters;
628  Double_t reco_beam_len, reco_beam_vtxX, reco_beam_vtxY, reco_beam_vtxZ, reco_beam_startX, reco_beam_startY, reco_beam_startZ, reco_beam_trackDirZ, reco_beam_interactingEnergy, reco_beam_Chi2_proton;
629  bool reco_beam_true_byHits_matched; // Does the true particle contributing most to the reconstructed beam track coincide with the actual beam particle that generated the event
630  Int_t reco_beam_true_byHits_origin, reco_beam_true_byHits_PDG; // Origin and PDG of the reconstructed beam track
631  Int_t true_beam_PDG;
632  Double_t true_beam_interactingEnergy;
633  std::string *true_beam_endProcess = 0;
634  std::string *reco_beam_true_byHits_endProcess = 0;
635  std::vector<double> *reco_beam_incidentEnergies = 0;
636 
637  //std::vector<int> * data_BI_PDG = 0x0;
638 
639  defaultTree->SetBranchAddress("reco_beam_type", &reco_beam_type);
640  defaultTree->SetBranchAddress("reco_beam_len", &reco_beam_len);
641  defaultTree->SetBranchAddress("reco_beam_vtxX", &reco_beam_vtxX);
642  defaultTree->SetBranchAddress("reco_beam_vtxY", &reco_beam_vtxY);
643  defaultTree->SetBranchAddress("reco_beam_vtxZ", &reco_beam_vtxZ);
644  defaultTree->SetBranchAddress("reco_beam_startX", &reco_beam_startX);
645  defaultTree->SetBranchAddress("reco_beam_startY", &reco_beam_startY);
646  defaultTree->SetBranchAddress("reco_beam_startZ", &reco_beam_startZ);
647  double reco_beam_endZ;
648  defaultTree->SetBranchAddress("reco_beam_endZ", &reco_beam_endZ);
649  defaultTree->SetBranchAddress("reco_beam_trackDirZ", &reco_beam_trackDirZ);
650  defaultTree->SetBranchAddress("reco_beam_nTrackDaughters", &reco_beam_nTrackDaughters);
651  defaultTree->SetBranchAddress("reco_beam_nShowerDaughters", &reco_beam_nShowerDaughters);
652  defaultTree->SetBranchAddress("reco_beam_interactingEnergy", &reco_beam_interactingEnergy);
653  defaultTree->SetBranchAddress("reco_beam_Chi2_proton", &reco_beam_Chi2_proton);
654  defaultTree->SetBranchAddress("reco_beam_incidentEnergies", &reco_beam_incidentEnergies);
655 
656  defaultTree->SetBranchAddress("reco_beam_true_byHits_matched", &reco_beam_true_byHits_matched);
657  defaultTree->SetBranchAddress("reco_beam_true_byHits_origin", &reco_beam_true_byHits_origin);
658  defaultTree->SetBranchAddress("reco_beam_true_byHits_PDG", &reco_beam_true_byHits_PDG);
659  defaultTree->SetBranchAddress("reco_beam_true_byHits_endProcess", &reco_beam_true_byHits_endProcess);
660 
661  defaultTree->SetBranchAddress("true_beam_interactingEnergy", &true_beam_interactingEnergy);
662  defaultTree->SetBranchAddress("true_beam_PDG", &true_beam_PDG);
663  defaultTree->SetBranchAddress("true_beam_endProcess", &true_beam_endProcess);
664  std::vector<double> *reco_beam_calo_wire = 0;
665  defaultTree->SetBranchAddress("reco_beam_calo_wire", &reco_beam_calo_wire);
666 
667  channel.erase(std::remove(channel.begin(), channel.end(), '.'), channel.end());
668  channel.erase(std::remove(channel.begin(), channel.end(), ' '), channel.end());
669 
670  const int nrecobins = recoBins.size();
671  TString hist_name = Form("Data_Channel%s_Histo",channel.c_str());
672  TString hist_title = Form("Data for channel %s",channel.c_str());
673  TH1D* datahisto = new TH1D(hist_name, hist_title,
674  (doNegativeReco ? nrecobins : nrecobins-1),
675  (doNegativeReco ? -1 : 0),
676  nrecobins-1);
677  if(IsIncidentHisto)
678  datahisto->SetNameTitle("Data_ChannelIncident_Histo", "Incident Data");
679  datahisto->SetDirectory(0);
680 
681  mf::LogInfo("FillDataHistogram_Pions") << "Filling data histogram " << datahisto->GetName() << " from file " << filename.c_str() << " for channel " << channel.c_str();
682 
683  double pitch = 0.4792;
684  double z0 = 0.56035;
685  int slice_cut = std::floor((reco_beam_endZ_cut - (z0 - pitch/2.)) / pitch);
686  for(Int_t k=0; k < defaultTree->GetEntries(); k++){
687  defaultTree->GetEntry(k);
688 
689  // Sometimes the reco energy at vertex is mis-reconstructed
690  if (!doNegativeReco && reco_beam_interactingEnergy < 0.0) continue;
691 
692  if (reco_beam_interactingEnergy == -999.) continue;
693 
694  if (IsIncidentHisto) {
695 
696  for (size_t l = 0; l < reco_beam_incidentEnergies->size(); ++l) {
697 
698  //Here -- need to add in a check for the reconstructed slice.
699  //Previously, was using an input file that already had these removed,
700  //but now I need to add this in.
701  if ((*reco_beam_calo_wire)[l] > slice_cut)
702  continue;
703 
704  double energy = (*reco_beam_incidentEnergies)[l];
705  if (doNegativeReco) {
706  if (energy < 0.) {
707  datahisto->AddBinContent(1, 1.);
708  }
709  else {
710  for (int m = 1; m < nrecobins; m++) {
711  if (energy > recoBins[m-1] &&
712  energy <= recoBins[m]) {
713  datahisto->AddBinContent(m+1, 1.);
714  break;
715  }
716  }
717  }
718  }
719  else {
720  for (int m = 1; m < nrecobins; m++) {
721  if (energy > recoBins[m-1] && energy <= recoBins[m]) {
722  datahisto->AddBinContent(m, 1);
723  break;
724  }
725  }
726  }
727  }
728  //continue;
729  }
730  else {
731  // Fill histogram in reco energy
732  if (doNegativeReco) {
733  if (reco_beam_interactingEnergy < 0.) {
734  datahisto->AddBinContent(1, 1);
735  }
736  else {
737  for (int l = 1; l < nrecobins; l++) {
738  if (reco_beam_interactingEnergy > recoBins[l-1] &&
739  reco_beam_interactingEnergy <= recoBins[l]) {
740  datahisto->AddBinContent(l+1, 1);
741  break;
742  }
743  }
744  }
745  }
746  else {
747  for (int l = 1; l < nrecobins; l++) {
748  if (reco_beam_interactingEnergy > recoBins[l-1] &&
749  reco_beam_interactingEnergy <= recoBins[l]) {
750  datahisto->AddBinContent(l, 1);
751  break;
752  }
753  }
754  }
755  }
756 
757  }
758 
759  file->Close();
760 
761  return datahisto;
762 
763 }
generate_CCQE_events.energy
string energy
Definition: generate_CCQE_events.py:15
python.larbatch_posix.remove
def remove(path)
Definition: larbatch_posix.py:715
string
std::string string
Definition: nybbler.cc:12
mf::LogInfo
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Definition: MessageLogger.h:209
channel
uint8_t channel
Definition: CRTFragment.hh:201
train.filename
string filename
Definition: train.py:213
l
static QStrList * l
Definition: config.cpp:1044
cache_state.file
file
Definition: cache_state.py:421
cache_state.m
m
Definition: cache_state.py:415
TH1 * protoana::ProtoDUNESelectionUtils::FillDataSidebandHistogram_Pions ( std::string  filename,
std::string  treename,
std::string  channel,
double  endZ_cut,
std::vector< double >  binning,
double  weight = 1. 
)

Definition at line 1779 of file ProtoDUNESelectionUtils.cxx.

1782  {
1783 //********************************************************************
1784 
1785  TFile *file = new TFile(filename.c_str(), "READ");
1786  TTree *defaultTree = (TTree*)file->Get(treename.c_str());
1787 
1788  std::vector<double> *reco_beam_calo_wire = 0, *reco_beam_calo_wire_z = 0, *reco_beam_incidentEnergies = 0;
1789  double reco_beam_interactingEnergy;
1790  defaultTree->SetBranchAddress("reco_beam_calo_wire", &reco_beam_calo_wire);
1791  defaultTree->SetBranchAddress("reco_beam_calo_wire_z", &reco_beam_calo_wire_z);
1792  defaultTree->SetBranchAddress("reco_beam_interactingEnergy", &reco_beam_interactingEnergy);
1793  defaultTree->SetBranchAddress("reco_beam_incidentEnergies", &reco_beam_incidentEnergies);
1794 
1795  std::string hist_name = "Data_SidebandChannel" + channel + "_Histo";
1796  std::string hist_title = "Data Sideband topology " + channel;
1797 
1798  TH1D* hist = new TH1D(hist_name.c_str(), hist_title.c_str(), binning.size()-1, 0., binning.size()-1);
1799  hist->SetDirectory(0);
1800 
1801  double pitch = 0.4792;
1802  double z0 = 0.56035;
1803  int slice_cut = std::floor((endZ_cut - (z0 - pitch/2.)) / pitch);
1804 
1805  for (int i = 0; i < defaultTree->GetEntries(); ++i) {
1806 
1807  defaultTree->GetEntry(i);
1808  if (reco_beam_interactingEnergy == -999.) continue;
1809 
1810  for (size_t j = 0; j < reco_beam_incidentEnergies->size(); ++j) {
1811 
1812  if ((*reco_beam_calo_wire)[j] <= slice_cut)
1813  continue;
1814 
1815  //hist->Fill((*reco_beam_calo_wire_z)[j], weight);
1816  int bin = FindBin((*reco_beam_calo_wire_z)[j], binning);
1817  //std::cout << "Got bin: " << bin << std::endl;
1818  hist->AddBinContent(bin, weight);
1819  }
1820  }
1821 
1822  file->Close();
1823  return hist;
1824 }
string
std::string string
Definition: nybbler.cc:12
protoana::ProtoDUNESelectionUtils::FindBin
int FindBin(double value, std::vector< double > binning)
Definition: ProtoDUNESelectionUtils.cxx:1826
channel
uint8_t channel
Definition: CRTFragment.hh:201
train.filename
string filename
Definition: train.py:213
test.weight
weight
Definition: test.py:257
cache_state.file
file
Definition: cache_state.py:421
bin
QTextStream & bin(QTextStream &s)
Definition: qtextstream.cpp:2006
TH1 * protoana::ProtoDUNESelectionUtils::FillMCBackgroundHistogram_Pions ( std::string  filename,
std::string  treename,
std::vector< double >  recoBins,
std::string  channel,
std::string  topo,
int  toponum,
double  endZ_cut,
double  minval = 0.0,
double  maxval = 100000.0,
bool  doNegativeReco = false,
int  doSyst = 0,
std::string  systName = "",
double  weight = 1.0,
std::pair< double, double >  PiMuScale = {1., 1.} 
)

Definition at line 13 of file ProtoDUNESelectionUtils.cxx.

17  {
18 //********************************************************************
19 
20  TFile *file = new TFile(filename.c_str(), "READ");
21  TTree *defaultTree = (TTree*)file->Get(treename.c_str());
22 
23  int reco_beam_type; // 13 -> track-like, 11 -> shower-like
24  int reco_beam_nTrackDaughters, reco_beam_nShowerDaughters;
25  double reco_beam_len, reco_beam_vtxX, reco_beam_vtxY, reco_beam_vtxZ,
26  reco_beam_startX, reco_beam_startY, reco_beam_startZ, reco_beam_trackDirZ,
27  reco_beam_interactingEnergy, reco_beam_Chi2_proton;
28 
29  // Does the true particle contributing most to the reconstructed beam track
30  // coincide with the actual beam particle that generated the event
31  bool reco_beam_true_byHits_matched;
32 
33  // Origin and PDG of the reconstructed beam track
34  int reco_beam_true_byHits_origin, reco_beam_true_byHits_PDG;
35  int true_beam_PDG;
36  double true_beam_endZ;
37  double true_beam_interactingEnergy;
38  std::string *true_beam_endProcess = 0;
39  std::string *reco_beam_true_byHits_endProcess = 0;
40  std::vector<double> *reco_beam_incidentEnergies = 0;
41 
42  int true_chexSignal, true_absSignal, true_backGround, true_nPi0Signal;
43 
44  //new topology
45  int interaction_topology;
46  defaultTree->SetBranchAddress("interaction_topology", &interaction_topology);
47 
48  defaultTree->SetBranchAddress("reco_beam_type", &reco_beam_type);
49  defaultTree->SetBranchAddress("reco_beam_len", &reco_beam_len);
50  defaultTree->SetBranchAddress("reco_beam_vtxX", &reco_beam_vtxX);
51  defaultTree->SetBranchAddress("reco_beam_vtxY", &reco_beam_vtxY);
52  defaultTree->SetBranchAddress("reco_beam_vtxZ", &reco_beam_vtxZ);
53  defaultTree->SetBranchAddress("reco_beam_startX", &reco_beam_startX);
54  defaultTree->SetBranchAddress("reco_beam_startY", &reco_beam_startY);
55  defaultTree->SetBranchAddress("reco_beam_startZ", &reco_beam_startZ);
56  defaultTree->SetBranchAddress("reco_beam_trackDirZ", &reco_beam_trackDirZ);
57  defaultTree->SetBranchAddress("reco_beam_nTrackDaughters", &reco_beam_nTrackDaughters);
58  defaultTree->SetBranchAddress("reco_beam_nShowerDaughters", &reco_beam_nShowerDaughters);
59  defaultTree->SetBranchAddress("reco_beam_interactingEnergy", &reco_beam_interactingEnergy);
60  defaultTree->SetBranchAddress("reco_beam_Chi2_proton", &reco_beam_Chi2_proton);
61  defaultTree->SetBranchAddress("reco_beam_incidentEnergies", &reco_beam_incidentEnergies);
62 
63  defaultTree->SetBranchAddress("reco_beam_true_byHits_matched", &reco_beam_true_byHits_matched);
64  defaultTree->SetBranchAddress("reco_beam_true_byHits_origin", &reco_beam_true_byHits_origin);
65  defaultTree->SetBranchAddress("reco_beam_true_byHits_PDG", &reco_beam_true_byHits_PDG);
66  defaultTree->SetBranchAddress("reco_beam_true_byHits_endProcess", &reco_beam_true_byHits_endProcess);
67 
68  defaultTree->SetBranchAddress("true_beam_interactingEnergy", &true_beam_interactingEnergy);
69  defaultTree->SetBranchAddress("true_beam_PDG", &true_beam_PDG);
70  defaultTree->SetBranchAddress("true_beam_endZ", &true_beam_endZ);
71  defaultTree->SetBranchAddress("true_beam_endProcess", &true_beam_endProcess);
72 
73  defaultTree->SetBranchAddress("true_chexSignal", &true_chexSignal);
74  defaultTree->SetBranchAddress("true_nPi0Signal", &true_nPi0Signal);
75  defaultTree->SetBranchAddress("true_absSignal", &true_absSignal);
76  defaultTree->SetBranchAddress("true_backGround", &true_backGround);
77 
78  //Testing
79  std::vector<double> * reco_beam_calibrated_dEdX = 0x0;
80  defaultTree->SetBranchAddress("reco_beam_calibrated_dEdX", &reco_beam_calibrated_dEdX);
81  std::vector<double> * reco_beam_TrkPitch = 0x0;
82  defaultTree->SetBranchAddress("reco_beam_TrkPitch", &reco_beam_TrkPitch);
83 
84  double reco_beam_endX, reco_beam_endY, reco_beam_endZ;
85  defaultTree->SetBranchAddress("reco_beam_endX", &reco_beam_endX);
86  defaultTree->SetBranchAddress("reco_beam_endY", &reco_beam_endY);
87  defaultTree->SetBranchAddress("reco_beam_endZ", &reco_beam_endZ);
88 
89  int run, event;
90  defaultTree->SetBranchAddress("event", &event);
91  defaultTree->SetBranchAddress("run", &run);
92  /////////////////////////////////
93 
94  std::vector<int> * reco_beam_hit_true_origin = 0x0;
95  std::vector<int> * reco_beam_hit_true_ID = 0x0;
96  std::vector<int> * true_beam_daughter_ID = 0x0;
97  std::vector<int> * true_beam_grand_daughter_ID = 0x0;
98  int true_beam_ID;
99  int true_daughter_nPiPlus, true_daughter_nPiMinus, true_daughter_nPi0;
100  defaultTree->SetBranchAddress( "reco_beam_hit_true_origin", &reco_beam_hit_true_origin );
101  defaultTree->SetBranchAddress( "reco_beam_hit_true_ID", &reco_beam_hit_true_ID );
102  defaultTree->SetBranchAddress( "true_beam_ID", &true_beam_ID );
103  defaultTree->SetBranchAddress( "true_daughter_nPiPlus", &true_daughter_nPiPlus );
104  defaultTree->SetBranchAddress( "true_daughter_nPiMinus", &true_daughter_nPiMinus );
105  defaultTree->SetBranchAddress( "true_daughter_nPi0", &true_daughter_nPi0 );
106  defaultTree->SetBranchAddress( "true_beam_daughter_ID", &true_beam_daughter_ID );
107  defaultTree->SetBranchAddress( "true_beam_grand_daughter_ID", &true_beam_grand_daughter_ID );
108 
109  //For vertex type
110  std::vector< std::string > * true_beam_processes = 0x0;
111  defaultTree->SetBranchAddress( "true_beam_processes", &true_beam_processes );
112  //std::vector< std::vector < double > > * reco_beam_vertex_dRs = 0x0;
113  //std::vector< int > * reco_beam_vertex_hits_slices = 0x0;
114  //defaultTree->SetBranchAddress( "reco_beam_vertex_dRs", &reco_beam_vertex_dRs );
115  //defaultTree->SetBranchAddress( "reco_beam_vertex_hits_slices", &reco_beam_vertex_hits_slices );
116 
117  //For systs
118  std::vector<double> * g4rw_primary_plus_sigma_weight = 0x0;
119  std::vector<double> * g4rw_primary_minus_sigma_weight = 0x0;
120  std::vector<std::string> * g4rw_primary_var = 0x0;
121  if (doSyst != 0) {
122  defaultTree->SetBranchAddress("g4rw_primary_plus_sigma_weight", &g4rw_primary_plus_sigma_weight);
123  defaultTree->SetBranchAddress("g4rw_primary_minus_sigma_weight", &g4rw_primary_minus_sigma_weight);
124  defaultTree->SetBranchAddress("g4rw_primary_var", &g4rw_primary_var);
125  }
126 
127  channel.erase(std::remove(channel.begin(), channel.end(), '.'), channel.end());
128  channel.erase(std::remove(channel.begin(), channel.end(), ' '), channel.end());
129  topo.erase(std::remove(topo.begin(), topo.end(), '.'), topo.end());
130  topo.erase(std::remove(topo.begin(), topo.end(), ' '), topo.end());
131 
132  const int nrecobins = recoBins.size();
133  std::string hist_name = "MC_Channel" + channel + "_" + topo + "_Histo";
134  std::string hist_title = "MC Background for channel " + channel +
135  " and topology " + topo;
136 
137  if (doSyst == 1) {
138  hist_name += "_high_" + systName;
139  hist_title += " +1 sigma";
140  }
141  else if (doSyst == -1) {
142  hist_name += "_low_" + systName;
143  hist_title += " -1 sigma";
144  }
145 
146  TH1D* mchisto = new TH1D(hist_name.c_str(), hist_title.c_str(),
147  (doNegativeReco ? nrecobins : (nrecobins-1)),
148  (doNegativeReco ? -1 : 0), nrecobins-1);
149 
150  mchisto->SetDirectory(0);
151 
152  mf::LogInfo("FillMCBackgroundHistogram_Pions") <<
153  "Filling MC background histogram " << mchisto->GetName() <<
154  " from file " << filename.c_str() << " for channel " <<
155  channel.c_str() << " with topology " << topo.c_str();
156 
157  bool done_check = false;
158  for(Int_t k=0; k < defaultTree->GetEntries(); k++){
159 
160  defaultTree->GetEntry(k);
161 
162  if (!done_check && doSyst != 0) {
163  if (g4rw_primary_var->size() > 0) {
164  std::cout << "Checking" << std::endl;
165  auto syst_check = std::find(g4rw_primary_var->begin(),
166  g4rw_primary_var->end(), systName);
167  if (syst_check == g4rw_primary_var->end()) {
168  std::cout << "Error! Could not find syst named " << systName << std::endl;
169  std::exception e;
170  throw e;
171  }
172  done_check = true;
173  }
174  }
175 
176  double syst_weight = 1.;
177  std::map<std::string, double> weights;
178  if (doSyst == 1) {
179  for (size_t i = 0; i < g4rw_primary_var->size(); ++i) {
180  weights[(*g4rw_primary_var)[i]] = (*g4rw_primary_plus_sigma_weight)[i];
181  }
182  }
183  else if (doSyst == -1) {
184  for (size_t i = 0; i < g4rw_primary_var->size(); ++i) {
185  weights[(*g4rw_primary_var)[i]] = (*g4rw_primary_minus_sigma_weight)[i];
186  }
187  }
188 
189  // Different background topologies
190  //Int_t topology = -1;
191 
192  //Bin edges cases. If there are signal events that fall out the truth bins,
193  //but are present in the reco bins then count them as backgrounds
194  if (true_backGround == 0 && (true_chexSignal == 1 || true_absSignal == 1 ||
195  true_nPi0Signal == 1)){
196  if (true_beam_interactingEnergy < minval &&
197  reco_beam_interactingEnergy > recoBins[0]){
198  true_backGround = 1;
199  }
200 
201  if (true_beam_interactingEnergy > maxval &&
202  reco_beam_interactingEnergy < recoBins[nrecobins-1]){
203  true_backGround = 1;
204  }
205  }
206 
207  /////////////////////
208 
209  if (!reco_beam_hit_true_ID->size()) {
210  continue;
211  }
212 
213 
214  if (true_beam_PDG == 211) {
215  if (!reco_beam_hit_true_ID->size()) {
216  continue;
217  }
218  if (doSyst == 1 || doSyst == -1) {
219  syst_weight = weights[systName];
220  }
221  }
222 
223  /*
224  if( reco_beam_true_byHits_origin == 2 ){
225  topology = 7;
226  }
227  else if (true_beam_PDG == 211) {
228 
229  //make sure there is any hits at all
230  if( !reco_beam_hit_true_ID->size() ){
231  //std::cout << "no true_ID for hits " << reco_beam_interactingEnergy <<
232  // " " << reco_beam_calibrated_dEdX->size() << " " <<
233  // reco_beam_TrkPitch->size() << std::endl;
234  //std::cout << "Checking: " << event << " " << run << std::endl;
235  continue;
236  }
237 
238  //maybe no true_beam_slices
239  if (true_beam_endZ < 0.) {
240  topology = 4;
241  }
242  else if (true_beam_endZ > endZ_cut) {
243  topology = 7;
244  }
245  else if ( *true_beam_endProcess == "pi+Inelastic" ){
246  if ( true_daughter_nPiPlus == 0 && true_daughter_nPiMinus == 0 ) {
247  topology = 1;
248  }
249  else {
250 
251  // bool has_pion_above_threshold = false;
252  // for (size_t i = 0; i < true_daughter_startP->size(); ++i) {
253  // if (abs((*true_daughter_PDG)[i]) == 211 &&
254  // (*true_daughter_startP)[i] > .150) {
255  // has_pion_above_threshold = true;
256  // break;
257  // }
258  // }
259 
260  // if (has_pion_above_threshold) {
261  // topology = 3;
262  // }
263  // else {
264  // topology = 1;
265  // }
266  topology = 3;
267  }
268  }
269  else {
270  topology = 7;
271  }
272 
273  if (doSyst == 1 || doSyst == -1) { //Do +1 sigma
274  syst_weight = weights[systName];
275  }
276  }
277  else if ( true_beam_PDG == -13 ){
278  //First check that if the last hit is from a cosmic
279  if ( reco_beam_hit_true_origin->back() == 2 )
280  topology = 7; //Other for now
281  else if( reco_beam_hit_true_ID->back() == true_beam_ID )
282  topology = 5;
283  else
284  topology = 7;
285  }
286  else{ //Shouldn't be any but w/e
287  topology = 7;
288  }
289 
290 
291  // Select only the correct topology
292  if (topology != toponum) continue;
293 
294  if (topology != interaction_topology) {
295  std::cout << "ERROR topo mismatch " << topology << " " << interaction_topology << std::endl;
296  }
297 */
298  if (interaction_topology == -1)
299  std::cout << "Warning: Interaction topo -1" << std::endl;
300 
301  if (interaction_topology != toponum)
302  continue;
303 
304  double total_weight = weight*syst_weight;
305  if (true_beam_PDG == 211) {
306  total_weight *= PiMuScale.first;
307  }
308  else if (true_beam_PDG == -13) {
309  total_weight *= PiMuScale.second;
310  }
311 
312  if (doNegativeReco) {
313  if (reco_beam_interactingEnergy < 0.) {
314  mchisto->AddBinContent(1, total_weight);
315  }
316  else{
317  for (int l = 1; l < nrecobins; ++l) {
318  if (reco_beam_interactingEnergy > recoBins[l-1] &&
319  reco_beam_interactingEnergy <= recoBins[l]) {
320  //Fill +1 because the first bin is negative
321  mchisto->AddBinContent(l+1, total_weight);
322  break;
323  }
324  }
325  }
326  }
327  else {
328  // Sometimes the reco energy at vertex is mis-reconstructed
329  if (reco_beam_interactingEnergy < 0.0) continue;
330 
331  for (int l = 1; l < nrecobins; ++l) {
332  if (reco_beam_interactingEnergy > recoBins[l-1] &&
333  reco_beam_interactingEnergy <= recoBins[l]) {
334  mchisto->AddBinContent(l/*+1*/, total_weight);
335  break;
336  }
337  }
338  }
339  }
340 
341  file->Close();
342 
343  return mchisto;
344 
345 }
python.larbatch_posix.remove
def remove(path)
Definition: larbatch_posix.py:715
string
std::string string
Definition: nybbler.cc:12
mf::LogInfo
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Definition: MessageLogger.h:209
channel
uint8_t channel
Definition: CRTFragment.hh:201
train.filename
string filename
Definition: train.py:213
l
static QStrList * l
Definition: config.cpp:1044
test.weight
weight
Definition: test.py:257
e
const double e
Definition: gUpMuFluxGen.cxx:165
cache_state.file
file
Definition: cache_state.py:421
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
TH1 * protoana::ProtoDUNESelectionUtils::FillMCFlux_Pions ( std::string  filename,
std::string  treename,
std::vector< double >  Bins,
int  mode = 1,
double  weight = 1.0 
)

Definition at line 1173 of file ProtoDUNESelectionUtils.cxx.

1173  {
1174  //********************************************************************
1175 
1176  TFile *file = new TFile(filename.c_str(), "READ");
1177  TTree *truthTree = (TTree*)file->Get(treename.c_str());
1178 
1179  Int_t reco_beam_type; // 13 -> track-like, 11 -> shower-like
1180  Double_t reco_beam_len, reco_beam_startX, reco_beam_startY, reco_beam_startZ, reco_beam_trackDirZ, reco_beam_interactingEnergy, reco_beam_Chi2_proton;
1181  Int_t true_beam_PDG, reco_beam_true_byHits_origin, reco_beam_true_byHits_PDG, reco_beam_nTrackDaughters, reco_beam_nShowerDaughters;
1182  Double_t true_beam_interactingEnergy, true_beam_endZ;
1183  std::vector<double> *true_beam_incidentEnergies = 0; std::vector<double> *reco_beam_incidentEnergies = 0;
1184 
1185  truthTree->SetBranchAddress("reco_beam_type", &reco_beam_type);
1186  truthTree->SetBranchAddress("reco_beam_len", &reco_beam_len);
1187  truthTree->SetBranchAddress("reco_beam_startX", &reco_beam_startX);
1188  truthTree->SetBranchAddress("reco_beam_startY", &reco_beam_startY);
1189  truthTree->SetBranchAddress("reco_beam_startZ", &reco_beam_startZ);
1190  truthTree->SetBranchAddress("reco_beam_trackDirZ", &reco_beam_trackDirZ);
1191  truthTree->SetBranchAddress("reco_beam_interactingEnergy", &reco_beam_interactingEnergy);
1192  truthTree->SetBranchAddress("reco_beam_Chi2_proton", &reco_beam_Chi2_proton);
1193  truthTree->SetBranchAddress("reco_beam_true_byHits_origin", &reco_beam_true_byHits_origin);
1194  truthTree->SetBranchAddress("reco_beam_true_byHits_PDG", &reco_beam_true_byHits_PDG);
1195  truthTree->SetBranchAddress("reco_beam_incidentEnergies", &reco_beam_incidentEnergies);
1196  truthTree->SetBranchAddress("reco_beam_nTrackDaughters", &reco_beam_nTrackDaughters);
1197  truthTree->SetBranchAddress("reco_beam_nShowerDaughters", &reco_beam_nShowerDaughters);
1198 
1199  truthTree->SetBranchAddress("true_beam_interactingEnergy", &true_beam_interactingEnergy);
1200  truthTree->SetBranchAddress("true_beam_PDG", &true_beam_PDG);
1201  truthTree->SetBranchAddress("true_beam_endZ", &true_beam_endZ);
1202  truthTree->SetBranchAddress("true_beam_incidentEnergies", &true_beam_incidentEnergies);
1203 
1204  const int nbins = Bins.size();
1205  TH1D* mchisto = new TH1D(Form("MC_PionFlux%i_Histo",mode), Form("MC Pion Flux - %i",mode), nbins-1, 0, nbins-1);
1206  mchisto->SetDirectory(0);
1207 
1208  mf::LogInfo("FillMCFlux_Pions") << "Filling MC pion flux histogram " << mchisto->GetName() << " from file " << filename.c_str() << " for mode " << mode;
1209 
1210  for(Int_t k=0; k < truthTree->GetEntries(); k++){
1211  truthTree->GetEntry(k);
1212 
1213  // Pion beam
1214  if(true_beam_PDG != 211) continue;
1215 
1216  // If the pion does not make at the front face of the TPC, then don't count. Probably need a better estimate.
1217  if(true_beam_endZ < 0.0) continue;
1218 
1219  if(true_beam_incidentEnergies->size() == 0 && reco_beam_incidentEnergies->size() == 0) continue;
1220 
1221  if(mode == 1){
1222  for(UInt_t l = 0; l < true_beam_incidentEnergies->size(); l++){
1223  for(Int_t m = 1; m <= nbins; m++){
1224  if(true_beam_incidentEnergies->at(l) > Bins[m-1] && true_beam_incidentEnergies->at(l) <= Bins[m]){
1225  mchisto->SetBinContent(m, mchisto->GetBinContent(m) + weight);
1226  break;
1227  }
1228  }
1229  }
1230  }
1231  else{
1232  if(reco_beam_true_byHits_origin != 4) continue;
1233  if(reco_beam_true_byHits_PDG != 211) continue;
1234  if(reco_beam_incidentEnergies->size() == 0) continue;
1235 
1236  // Track-like beam particles
1237  if(reco_beam_type != 13) continue;
1238 
1239  // Position cuts
1240  if(reco_beam_startZ < 28.0 || reco_beam_startZ > 32.0) continue;
1241  if(reco_beam_startY < 410.0 || reco_beam_startY > 445.0) continue;
1242  if(reco_beam_startX < -45.0 || reco_beam_startX > 0.0) continue;
1243 
1244  // Direction cuts
1245  if(reco_beam_trackDirZ < 0.9) continue;
1246 
1247  // Remove secondary protons from pion interacting outside the TPC
1248  if(reco_beam_Chi2_proton < 450.0) continue;
1249 
1250  // Minimum and maximum track length
1251  if(reco_beam_len < 5.0 || reco_beam_len > 270.0) continue;
1252 
1253  // Sometimes the reco energy at vertex is mis-reconstructed
1254  if(reco_beam_interactingEnergy < 0.0) continue;
1255 
1256  // For inelastic scattering at least one daughter is expected
1257  // Jake: This is deprecated in later versions
1258  Int_t allDaughters = reco_beam_nTrackDaughters + reco_beam_nShowerDaughters;
1259  if(allDaughters <= 0) continue;
1260 
1261  if(mode == 2){
1262  for(UInt_t l = 0; l < true_beam_incidentEnergies->size(); l++){
1263  for(Int_t m = 1; m < nbins; m++){
1264  if(true_beam_incidentEnergies->at(l) > Bins[m-1] && true_beam_incidentEnergies->at(l) <= Bins[m]){
1265  mchisto->SetBinContent(m, mchisto->GetBinContent(m) + weight);
1266  break;
1267  }
1268  }
1269  }
1270  }
1271  else if(mode == 3){
1272  for(UInt_t l = 0; l < reco_beam_incidentEnergies->size(); l++){
1273  for(Int_t m = 1; m < nbins; m++){
1274  if(reco_beam_incidentEnergies->at(l) > Bins[m-1] && reco_beam_incidentEnergies->at(l) <= Bins[m]){
1275  mchisto->SetBinContent(m, mchisto->GetBinContent(m) + weight);
1276  break;
1277  }
1278  }
1279  }
1280  }
1281  }
1282  }
1283 
1284  file->Close();
1285 
1286  return mchisto;
1287 
1288 }
mf::LogInfo
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Definition: MessageLogger.h:209
train.filename
string filename
Definition: train.py:213
l
static QStrList * l
Definition: config.cpp:1044
test.weight
weight
Definition: test.py:257
cache_state.file
file
Definition: cache_state.py:421
cache_state.m
m
Definition: cache_state.py:415
TH1 * protoana::ProtoDUNESelectionUtils::FillMCIncidentHistogram_Pions ( std::string  filename,
std::string  treename,
std::vector< double >  recoBins,
std::string  topo,
int  toponum,
double  reco_beam_endZ_cut,
double  minval = 0.,
double  maxval = 100000.,
bool  doNegativeReco = false,
int  doSyst = 0,
std::string  systName = "",
double  weight = 1.0,
std::pair< double, double >  PiMuScale = {1., 1.} 
)

Definition at line 766 of file ProtoDUNESelectionUtils.cxx.

771  {
772  //********************************************************************
773 
774  TFile *file = new TFile(filename.c_str(), "READ");
775  TTree *defaultTree = (TTree*)file->Get(treename.c_str());
776 
777  Int_t reco_beam_type; // 13 -> track-like, 11 -> shower-like
778  Int_t reco_beam_nTrackDaughters, reco_beam_nShowerDaughters;
779  Double_t reco_beam_len, reco_beam_vtxX, reco_beam_vtxY, reco_beam_vtxZ, reco_beam_startX, reco_beam_startY, reco_beam_startZ, reco_beam_trackDirZ, reco_beam_interactingEnergy, reco_beam_Chi2_proton;
780  bool reco_beam_true_byHits_matched; // Does the true particle contributing most to the reconstructed beam track coincide with the actual beam particle that generated the event
781  Int_t reco_beam_true_byHits_origin, reco_beam_true_byHits_PDG; // Origin and PDG of the reconstructed beam track
782  Int_t true_beam_PDG, true_beam_ID;
783  double true_beam_endZ;
784  Double_t true_beam_interactingEnergy;
785  std::string *true_beam_endProcess = 0;
786  std::string *reco_beam_true_byHits_endProcess = 0;
787  std::vector<double> *reco_beam_incidentEnergies = 0;
788  std::vector<double> *reco_beam_calo_wire = 0;
789 
790  //Int_t true_chexSignal, true_absSignal, true_backGround;
791 
792  defaultTree->SetBranchAddress("reco_beam_type", &reco_beam_type);
793  defaultTree->SetBranchAddress("reco_beam_len", &reco_beam_len);
794  defaultTree->SetBranchAddress("reco_beam_vtxX", &reco_beam_vtxX);
795  defaultTree->SetBranchAddress("reco_beam_vtxY", &reco_beam_vtxY);
796  defaultTree->SetBranchAddress("reco_beam_vtxZ", &reco_beam_vtxZ);
797  defaultTree->SetBranchAddress("reco_beam_startX", &reco_beam_startX);
798  defaultTree->SetBranchAddress("reco_beam_startY", &reco_beam_startY);
799  defaultTree->SetBranchAddress("reco_beam_startZ", &reco_beam_startZ);
800  double reco_beam_endZ;
801  defaultTree->SetBranchAddress("reco_beam_endZ", &reco_beam_endZ);
802  defaultTree->SetBranchAddress("reco_beam_trackDirZ", &reco_beam_trackDirZ);
803  defaultTree->SetBranchAddress("reco_beam_nTrackDaughters", &reco_beam_nTrackDaughters);
804  defaultTree->SetBranchAddress("reco_beam_nShowerDaughters", &reco_beam_nShowerDaughters);
805  defaultTree->SetBranchAddress("reco_beam_interactingEnergy", &reco_beam_interactingEnergy);
806  defaultTree->SetBranchAddress("reco_beam_Chi2_proton", &reco_beam_Chi2_proton);
807  defaultTree->SetBranchAddress("reco_beam_incidentEnergies", &reco_beam_incidentEnergies);
808  defaultTree->SetBranchAddress("reco_beam_calo_wire", &reco_beam_calo_wire);
809 
810  defaultTree->SetBranchAddress("reco_beam_true_byHits_matched", &reco_beam_true_byHits_matched);
811  defaultTree->SetBranchAddress("reco_beam_true_byHits_origin", &reco_beam_true_byHits_origin);
812  defaultTree->SetBranchAddress("reco_beam_true_byHits_PDG", &reco_beam_true_byHits_PDG);
813  defaultTree->SetBranchAddress("reco_beam_true_byHits_endProcess", &reco_beam_true_byHits_endProcess);
814 
815  defaultTree->SetBranchAddress("true_beam_interactingEnergy", &true_beam_interactingEnergy);
816  defaultTree->SetBranchAddress("true_beam_PDG", &true_beam_PDG);
817  defaultTree->SetBranchAddress("true_beam_endZ", &true_beam_endZ);
818  defaultTree->SetBranchAddress("true_beam_ID", &true_beam_ID);
819  defaultTree->SetBranchAddress("true_beam_endProcess", &true_beam_endProcess);
820 
821  //defaultTree->SetBranchAddress("true_chexSignal", &true_chexSignal);
822  //defaultTree->SetBranchAddress("true_absSignal", &true_absSignal);
823  //defaultTree->SetBranchAddress("true_backGround", &true_backGround);
824 
825  std::vector< int > * reco_beam_hit_true_ID = 0x0;
826  std::vector< int > * reco_beam_hit_true_origin = 0x0;
827  std::vector< int > * reco_beam_hit_true_slice = 0x0;
828 
829  defaultTree->SetBranchAddress("reco_beam_hit_true_ID", &reco_beam_hit_true_ID);
830  defaultTree->SetBranchAddress("reco_beam_hit_true_slice", &reco_beam_hit_true_slice);
831  defaultTree->SetBranchAddress("reco_beam_hit_true_origin", &reco_beam_hit_true_origin);
832 
833  std::vector< std::string > * true_beam_processes = 0x0;
834  defaultTree->SetBranchAddress("true_beam_processes", &true_beam_processes);
835 
836  std::vector< int > * true_beam_daughter_ID = 0x0;
837  std::vector< int > * true_beam_grand_daughter_ID = 0x0;
838  defaultTree->SetBranchAddress( "true_beam_daughter_ID", &true_beam_daughter_ID );
839  defaultTree->SetBranchAddress( "true_beam_grand_daughter_ID", &true_beam_grand_daughter_ID );
840  //std::vector< int > * reco_beam_vertex_hits_slices = 0x0;
841  //std::vector< std::vector<double> > * reco_beam_vertex_dRs = 0x0;
842  //defaultTree->SetBranchAddress( "reco_beam_vertex_hits_slices", &reco_beam_vertex_hits_slices );
843  //defaultTree->SetBranchAddress( "reco_beam_vertex_dRs", &reco_beam_vertex_dRs );
844 
845  //For systs
846  std::vector<double> * g4rw_primary_plus_sigma_weight = 0x0;
847  std::vector<double> * g4rw_primary_minus_sigma_weight = 0x0;
848  std::vector<std::string> * g4rw_primary_var = 0x0;
849  if (doSyst != 0) {
850  defaultTree->SetBranchAddress("g4rw_primary_plus_sigma_weight", &g4rw_primary_plus_sigma_weight);
851  defaultTree->SetBranchAddress("g4rw_primary_minus_sigma_weight", &g4rw_primary_minus_sigma_weight);
852  defaultTree->SetBranchAddress("g4rw_primary_var", &g4rw_primary_var);
853  }
854 
855  //Splitting by true energy
856  std::vector<int> * true_beam_slices = 0x0;
857  std::vector<double> * true_beam_incidentEnergies = 0x0;
858  defaultTree->SetBranchAddress( "true_beam_slices", &true_beam_slices );
859  defaultTree->SetBranchAddress("true_beam_incidentEnergies",
860  &true_beam_incidentEnergies);
861 
862  //std::replace(channel.begin(), channel.end(), ' ', '-');
863  //channel.erase(std::remove(channel.begin(), channel.end(), '.'), channel.end());
864  //channel.erase(std::remove(channel.begin(), channel.end(), ' '), channel.end());
865  //std::replace(topo.begin(), topo.end(), ' ', '-');
866  topo.erase(std::remove(topo.begin(), topo.end(), '.'), topo.end());
867  topo.erase(std::remove(topo.begin(), topo.end(), ' '), topo.end());
868 
869  size_t nrecobins = recoBins.size();
870 
871  TString hist_name = Form("MC_ChannelIncident_%s_Histo", topo.c_str());
872  //Form("MC_ChannelIncident%s_%s_Histo",
873  // channel.c_str(),topo.c_str());
874  TString hist_title = Form("Incident MC for topology %s", topo.c_str());
875  //Form("Incident MC for channel %s and topology %s",
876  // channel.c_str(),topo.c_str());
877 
878  if (doSyst == 1) {
879  hist_name += "_high_" + systName;
880  hist_title += " +1 sigma";
881  }
882  else if (doSyst == -1) {
883  hist_name += "_low_" + systName;
884  hist_title += " -1 sigma";
885  }
886 
887  //TH1D* mchisto = new TH1D(hist_name, hist_title, nrecobins-1, 0, nrecobins-1);
888  TH1D* mchisto = new TH1D(hist_name, hist_title,
889  (doNegativeReco ? nrecobins : (nrecobins-1)),
890  (doNegativeReco ? -1 : 0), nrecobins-1);
891  mchisto->SetDirectory(0);
892 
893  mf::LogInfo("FillMCIncidentHistogram_Pions") <<
894  "Filling MC incident histogram " << mchisto->GetName() <<
895  " from file " << filename.c_str() << /*" for channel " <<
896  channel.c_str() <<*/ " with topology " << topo.c_str();
897 
898  double pitch = 0.4792;
899  double z0 = 0.56035;
900  int slice_cut = std::floor((reco_beam_endZ_cut - (z0 - pitch/2.)) / pitch);
901 
902  bool done_check = false;
903  for(Int_t k=0; k < defaultTree->GetEntries(); k++){
904 
905  defaultTree->GetEntry(k);
906 
907  ////For first attempt at g4rw systematics
908  if (!done_check && doSyst != 0) {
909  if (g4rw_primary_var->size() > 0) {
910  std::cout << "Checking" << std::endl;
911  auto syst_check = std::find(g4rw_primary_var->begin(), g4rw_primary_var->end(), systName);
912  if (syst_check == g4rw_primary_var->end()) {
913  std::cout << "Error! Could not find syst named " << systName << std::endl;
914  std::exception e;
915  throw e;
916  }
917  done_check = true;
918  }
919  }
920 
921  double syst_weight = 1.;
922  std::map<std::string, double> weights;
923  if (doSyst == 1) {
924  for (size_t i = 0; i < g4rw_primary_var->size(); ++i) {
925  weights[(*g4rw_primary_var)[i]] = (*g4rw_primary_plus_sigma_weight)[i];
926  }
927  }
928  else if (doSyst == -1) {
929  for (size_t i = 0; i < g4rw_primary_var->size(); ++i) {
930  weights[(*g4rw_primary_var)[i]] = (*g4rw_primary_minus_sigma_weight)[i];
931  }
932  }
933  //////////////////////////////////////
934 
935  // Different background topologies
936  Int_t topology = -1;
937 
938  // Sometimes the reco energy at vertex is mis-reconstructed
939  if (!doNegativeReco && reco_beam_interactingEnergy < 0.0) continue;
940 
941  if (reco_beam_interactingEnergy == -999.) continue;
942 
943  //bool check_cosmics = ( reco_beam_true_byHits_origin == 2 );
944 
945  // Go through the incident energies from the beam particle
946  for(size_t l = 0; l < reco_beam_incidentEnergies->size(); ++l){
947 
948  //Here -- need to add in a check for the reconstructed slice.
949  //Previously, was using an input file that already had these removed,
950  //but now I need to add this in.
951  if ((*reco_beam_calo_wire)[l] > slice_cut)
952  continue;
953 
954  // Check the ID, origin, true_slice
955  int true_id = (*reco_beam_hit_true_ID)[l];
956  int true_origin = (*reco_beam_hit_true_origin)[l];
957  int true_slice = (*reco_beam_hit_true_slice)[l];
958 
959  /*
960  if ( check_cosmics && true_origin != 2 ) {
961  std::cout << "Notice! Beam matched to cosmic, with non-cosmic hit"
962  << std::endl;
963  }
964  */
965  double true_energy = 0.;
966 
967  // Cosmic
968  if( true_origin == 2 ){
969  topology = 5;
970  }
971  else{
972  if( true_id != true_beam_ID ){
973  if ( true_id == -999 )
974  topology = 8;
975  else {
976  //maybe no true_beam_slices instead?
977  if ( true_beam_endZ < 0. ) //Beam particle ends before TPC
978  topology = 7; //Just consider it downstream
979  else {
980  // Consider it downstream. This is slightly different than just above
981  // i.e. if the beam particle interacts/decays before the TPC
982  auto daughter_ID_check = std::find(
983  true_beam_daughter_ID->begin(), true_beam_daughter_ID->end(),
984  true_id);
985 
986  auto g_daughter_ID_check = std::find(
987  true_beam_grand_daughter_ID->begin(),
988  true_beam_grand_daughter_ID->end(),
989  true_id);
990 
991  if (daughter_ID_check != true_beam_daughter_ID->end()) {
992  topology = 7;
993  }
994  else if (g_daughter_ID_check != true_beam_grand_daughter_ID->end()){
995  topology = 7;
996  }
997  else {
998  topology = 8;
999  }
1000 
1001  }
1002  }
1003  }
1004  else{
1005  // True id matched to hit, but no slice matched
1006  // i.e. This is a 'messy' event
1007  if (true_slice == -999) {
1008  topology = 6;
1009  }
1010  else{
1011  if (true_beam_PDG == 211) {
1012  if (true_slice > slice_cut) {
1013  topology = 8; // Past endZ cut
1014  }
1015  else {
1016  topology = 3; // Is Pion
1017 
1018  //Go through the true slices to find the
1019  //true incident energy
1020  bool found_true_slice = false;
1021  for (size_t i = 0; i < true_beam_slices->size(); ++i) {
1022  int check_slice = (*true_beam_slices)[i];
1023  true_energy = (*true_beam_incidentEnergies)[i];
1024  if (true_slice == check_slice) {
1025  //std::cout << "Reco inc energy: " <<
1026  // (*reco_beam_incidentEnergies)[l] << " " <<
1027  // "True inc energy: " << true_energy <<
1028  // std::endl;
1029  found_true_slice = true;
1030  break;
1031  }
1032  }
1033  if (!found_true_slice) {
1034  std::cout << "Could not find true slice" << std::endl;
1035  }
1036  }
1037  }
1038  else if (true_beam_PDG == -13) {
1039  topology = 4; // Is muon
1040  }
1041 
1042  if (doSyst == 1 || doSyst == -1) { //Do +1 sigma
1043  syst_weight = weights[systName];
1044  }
1045  }
1046  }
1047  }
1048 
1049 
1050  if (topology == -1)
1051  std::cout << "Warning: incident topology == -1" <<std::endl;
1052  if (topology != toponum) {
1053  continue;
1054  }
1055 
1056  double total_weight = weight*syst_weight;
1057  if (true_beam_PDG == 211) {
1058  total_weight *= PiMuScale.first;
1059  }
1060  else if (true_beam_PDG == -13) {
1061  total_weight *= PiMuScale.second;
1062  }
1063 
1064  //Patch
1065  //if (true_energy < 0.)
1066  // true_energy = 10.;
1067  if (true_energy < minval || true_energy >= maxval) {
1068  //if (true_energy > maxval) {
1069  // std::cout << "TrueEnergy G: " << true_energy << " " << minval <<
1070  // " " << maxval << " " << (*reco_beam_incidentEnergies)[l] << std::endl;
1071  //}
1072  //else if (true_energy == maxval) {
1073  // std::cout << "TrueEnergy E: " << true_energy << " " << minval <<
1074  // " " << maxval << " " << (*reco_beam_incidentEnergies)[l] << std::endl;
1075  //}
1076  //else if (true_energy < minval) {
1077  // std::cout << "TrueEnergy L: " << true_energy << " " << minval <<
1078  // " " << maxval << " " << (*reco_beam_incidentEnergies)[l] << std::endl;
1079  //}
1080  continue;
1081  }
1082 
1083  double energy = (*reco_beam_incidentEnergies)[l];
1084  if (doNegativeReco) {
1085  if (energy < 0.) {
1086  mchisto->AddBinContent(1, total_weight);
1087  }
1088  else{
1089  for (size_t m = 1; m < nrecobins; ++m) {
1090  if (energy > recoBins[m-1] &&
1091  energy <= recoBins[m]) {
1092  //Fill +1 because the first bin is negative
1093  mchisto->AddBinContent(m+1, total_weight);
1094  break;
1095  }
1096  }
1097  }
1098  }
1099  else {
1100  // Sometimes the reco energy at vertex is mis-reconstructed
1101  if (energy < 0.0) continue;
1102 
1103  for (size_t m = 1; m < nrecobins; ++m) {
1104  if (energy > recoBins[m-1] &&
1105  energy <= recoBins[m]) {
1106  mchisto->AddBinContent(m/*+1*/, total_weight);
1107  break;
1108  }
1109  }
1110  }
1111 
1112  }
1113  }
1114 
1115  file->Close();
1116 
1117  return mchisto;
1118 
1119 }
generate_CCQE_events.energy
string energy
Definition: generate_CCQE_events.py:15
python.larbatch_posix.remove
def remove(path)
Definition: larbatch_posix.py:715
string
std::string string
Definition: nybbler.cc:12
mf::LogInfo
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Definition: MessageLogger.h:209
train.filename
string filename
Definition: train.py:213
l
static QStrList * l
Definition: config.cpp:1044
test.weight
weight
Definition: test.py:257
e
const double e
Definition: gUpMuFluxGen.cxx:165
cache_state.file
file
Definition: cache_state.py:421
cache_state.m
m
Definition: cache_state.py:415
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
TH1 * protoana::ProtoDUNESelectionUtils::FillMCSidebandHistogram_Pions ( std::string  filename,
std::string  treename,
std::string  channel,
std::string  topo,
int  toponum,
double  endZ_cut,
std::vector< double >  binning,
double  weight = 1. 
)

Definition at line 1643 of file ProtoDUNESelectionUtils.cxx.

1646  {
1647 //********************************************************************
1648 
1649  TFile *file = new TFile(filename.c_str(), "READ");
1650  TTree *defaultTree = (TTree*)file->Get(treename.c_str());
1651 
1652  int reco_beam_type; // 13 -> track-like, 11 -> shower-like
1653  int reco_beam_nTrackDaughters, reco_beam_nShowerDaughters;
1654  double reco_beam_len, reco_beam_vtxX, reco_beam_vtxY, reco_beam_vtxZ,
1655  reco_beam_startX, reco_beam_startY, reco_beam_startZ, reco_beam_trackDirZ,
1656  reco_beam_interactingEnergy, reco_beam_Chi2_proton;
1657 
1658  // Origin and PDG of the reconstructed beam track
1659  int true_beam_PDG;
1660  double true_beam_endZ;
1661  double true_beam_interactingEnergy;
1662  std::string *true_beam_endProcess = 0;
1663  std::vector<double> *reco_beam_incidentEnergies = 0;
1664 
1665  defaultTree->SetBranchAddress("reco_beam_type", &reco_beam_type);
1666  defaultTree->SetBranchAddress("reco_beam_len", &reco_beam_len);
1667  defaultTree->SetBranchAddress("reco_beam_vtxX", &reco_beam_vtxX);
1668  defaultTree->SetBranchAddress("reco_beam_vtxY", &reco_beam_vtxY);
1669  defaultTree->SetBranchAddress("reco_beam_vtxZ", &reco_beam_vtxZ);
1670  defaultTree->SetBranchAddress("reco_beam_startX", &reco_beam_startX);
1671  defaultTree->SetBranchAddress("reco_beam_startY", &reco_beam_startY);
1672  defaultTree->SetBranchAddress("reco_beam_startZ", &reco_beam_startZ);
1673  defaultTree->SetBranchAddress("reco_beam_trackDirZ", &reco_beam_trackDirZ);
1674  defaultTree->SetBranchAddress("reco_beam_nTrackDaughters", &reco_beam_nTrackDaughters);
1675  defaultTree->SetBranchAddress("reco_beam_nShowerDaughters", &reco_beam_nShowerDaughters);
1676  defaultTree->SetBranchAddress("reco_beam_interactingEnergy", &reco_beam_interactingEnergy);
1677  defaultTree->SetBranchAddress("reco_beam_Chi2_proton", &reco_beam_Chi2_proton);
1678  defaultTree->SetBranchAddress("reco_beam_incidentEnergies", &reco_beam_incidentEnergies);
1679 
1680  defaultTree->SetBranchAddress("true_beam_interactingEnergy", &true_beam_interactingEnergy);
1681  defaultTree->SetBranchAddress("true_beam_PDG", &true_beam_PDG);
1682  defaultTree->SetBranchAddress("true_beam_endZ", &true_beam_endZ);
1683  defaultTree->SetBranchAddress("true_beam_endProcess", &true_beam_endProcess);
1684 
1685  std::vector<int> * reco_beam_hit_true_origin = 0x0;
1686  std::vector<int> * reco_beam_hit_true_ID = 0x0;
1687  int true_beam_ID;
1688  defaultTree->SetBranchAddress( "true_beam_ID", &true_beam_ID );
1689  defaultTree->SetBranchAddress( "reco_beam_hit_true_origin", &reco_beam_hit_true_origin );
1690  defaultTree->SetBranchAddress( "reco_beam_hit_true_ID", &reco_beam_hit_true_ID );
1691 
1692  std::vector<int> * true_beam_daughter_ID = 0, * true_beam_daughter_PDG = 0;
1693  std::vector<int> * true_beam_grand_daughter_ID = 0x0;
1694  std::vector<double> *reco_beam_calo_wire = 0, *reco_beam_calo_wire_z = 0;
1695  defaultTree->SetBranchAddress("reco_beam_calo_wire", &reco_beam_calo_wire);
1696  defaultTree->SetBranchAddress("true_beam_daughter_ID",
1697  &true_beam_daughter_ID);
1698  defaultTree->SetBranchAddress("true_beam_daughter_PDG",
1699  &true_beam_daughter_PDG);
1700  defaultTree->SetBranchAddress("true_beam_grand_daughter_ID",
1701  &true_beam_grand_daughter_ID);
1702  defaultTree->SetBranchAddress("reco_beam_calo_wire_z",
1703  &reco_beam_calo_wire_z);
1704 
1705 
1706  std::string hist_name = "MC_SidebandChannel" + channel + "_" + topo + "_Histo";
1707  std::string hist_title = "MC Sideband Channel " + channel +
1708  " topology " + topo;
1709 
1710  TH1D* hist = new TH1D(hist_name.c_str(), hist_title.c_str(), binning.size()-1, 0., binning.size()-1);
1711  hist->SetDirectory(0);
1712 
1713  //TH1 * hist = new TH1D("MC_SidebandChannelMuons_Hist", "", 1, 0, 1);
1714 
1715  //hist->SetDirectory(0);
1716 
1717 
1718  double pitch = 0.4792;
1719  double z0 = 0.56035;
1720  int slice_cut = std::floor((endZ_cut - (z0 - pitch/2.)) / pitch);
1721 
1722  for (int i = 0; i < defaultTree->GetEntries(); ++i) {
1723 
1724  defaultTree->GetEntry(i);
1725  if (reco_beam_interactingEnergy == -999.) continue;
1726 
1727  for (size_t j = 0; j < reco_beam_incidentEnergies->size(); ++j) {
1728 
1729  if ((*reco_beam_calo_wire)[j] <= slice_cut)
1730  continue;
1731 
1732  int topology = -1;
1733 
1734  // Check the ID, origin
1735  int true_id = (*reco_beam_hit_true_ID)[j];
1736  //int true_origin = (*reco_beam_hit_true_origin)[j];
1737 
1738  auto daughter_find = std::find(true_beam_daughter_ID->begin(),
1739  true_beam_daughter_ID->end(), true_id);
1740  int daughter_index = daughter_find - true_beam_daughter_ID->begin();
1741  if (true_id == true_beam_ID && abs(true_beam_PDG) == 13) {
1742  topology = 1;
1743  }
1744  else if (true_id == true_beam_ID && abs(true_beam_PDG) == 211) {
1745  topology = 2;
1746  }
1747  else if (daughter_find != true_beam_daughter_ID->end()) {
1748  int daughter_PDG = (*true_beam_daughter_PDG)[daughter_index];
1749  if (abs(daughter_PDG) == 13) {
1750  topology = 3;
1751  }
1752  else {
1753  topology = 4;
1754  }
1755  }
1756  else if (std::find(true_beam_grand_daughter_ID->begin(),
1757  true_beam_grand_daughter_ID->end(), true_id) !=
1758  true_beam_grand_daughter_ID->end()) {
1759  topology = 4;
1760  }
1761  else {
1762  topology = 5;
1763  }
1764 
1765  if (topology == toponum) {
1766  //hist->Fill((*reco_beam_calo_wire_z)[j], weight);
1767  int bin = FindBin((*reco_beam_calo_wire_z)[j], binning);
1768  //std::cout << "Got bin: " << bin << std::endl;
1769  hist->AddBinContent(bin, weight);
1770  }
1771  }
1772  }
1773 
1774  file->Close();
1775  return hist;
1776 }
string
std::string string
Definition: nybbler.cc:12
protoana::ProtoDUNESelectionUtils::FindBin
int FindBin(double value, std::vector< double > binning)
Definition: ProtoDUNESelectionUtils.cxx:1826
channel
uint8_t channel
Definition: CRTFragment.hh:201
train.filename
string filename
Definition: train.py:213
test.weight
weight
Definition: test.py:257
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
cache_state.file
file
Definition: cache_state.py:421
bin
QTextStream & bin(QTextStream &s)
Definition: qtextstream.cpp:2006
TH1 * protoana::ProtoDUNESelectionUtils::FillMCSignalHistogram_Pions ( std::string  filename,
std::string  treename,
std::vector< double >  recoBins,
std::string  channel,
std::string  topo,
int  toponum,
double  minval,
double  maxval,
double  endZ_cut,
bool  doNegativeReco = false,
int  doSyst = 0,
std::string  systName = "",
double  weight = 1.0,
std::pair< double, double >  PiMuScale = {1., 1.} 
)

Definition at line 348 of file ProtoDUNESelectionUtils.cxx.

352  {
353 //********************************************************************
354 
355  TFile *file = new TFile(filename.c_str(), "READ");
356  TTree *defaultTree = (TTree*)file->Get(treename.c_str());
357 
358  int reco_beam_type; // 13 -> track-like, 11 -> shower-like
359  int reco_beam_nTrackDaughters, reco_beam_nShowerDaughters;
360  double reco_beam_len, reco_beam_vtxX, reco_beam_vtxY, reco_beam_vtxZ, reco_beam_startX, reco_beam_startY, reco_beam_startZ, reco_beam_trackDirZ, reco_beam_interactingEnergy, reco_beam_Chi2_proton;
361  bool reco_beam_true_byHits_matched; // Does the true particle contributing most to the reconstructed beam track coincide with the actual beam particle that generated the event
362  int reco_beam_true_byHits_origin, reco_beam_true_byHits_PDG; // Origin and PDG of the reconstructed beam track
363  int true_beam_PDG;
364  double true_beam_endZ;
365  double true_beam_interactingEnergy;
366  std::string *true_beam_endProcess = 0;
367  std::string *reco_beam_true_byHits_endProcess = 0;
368  std::vector<double> *reco_beam_incidentEnergies = 0;
369  int true_chexSignal, true_absSignal, true_backGround, true_nPi0Signal;
370 
371  int event, run;
372 
373  defaultTree->SetBranchAddress("reco_beam_type", &reco_beam_type);
374  defaultTree->SetBranchAddress("reco_beam_len", &reco_beam_len);
375  defaultTree->SetBranchAddress("reco_beam_vtxX", &reco_beam_vtxX);
376  defaultTree->SetBranchAddress("reco_beam_vtxY", &reco_beam_vtxY);
377  defaultTree->SetBranchAddress("reco_beam_vtxZ", &reco_beam_vtxZ);
378  defaultTree->SetBranchAddress("reco_beam_startX", &reco_beam_startX);
379  defaultTree->SetBranchAddress("reco_beam_startY", &reco_beam_startY);
380  defaultTree->SetBranchAddress("reco_beam_startZ", &reco_beam_startZ);
381  defaultTree->SetBranchAddress("reco_beam_trackDirZ", &reco_beam_trackDirZ);
382  defaultTree->SetBranchAddress("reco_beam_nTrackDaughters", &reco_beam_nTrackDaughters);
383  defaultTree->SetBranchAddress("reco_beam_nShowerDaughters", &reco_beam_nShowerDaughters);
384  defaultTree->SetBranchAddress("reco_beam_interactingEnergy", &reco_beam_interactingEnergy);
385  defaultTree->SetBranchAddress("reco_beam_Chi2_proton", &reco_beam_Chi2_proton);
386  defaultTree->SetBranchAddress("reco_beam_incidentEnergies", &reco_beam_incidentEnergies);
387 
388  defaultTree->SetBranchAddress("reco_beam_true_byHits_matched", &reco_beam_true_byHits_matched);
389  defaultTree->SetBranchAddress("reco_beam_true_byHits_origin", &reco_beam_true_byHits_origin);
390  defaultTree->SetBranchAddress("reco_beam_true_byHits_PDG", &reco_beam_true_byHits_PDG);
391  defaultTree->SetBranchAddress("reco_beam_true_byHits_endProcess", &reco_beam_true_byHits_endProcess);
392 
393  defaultTree->SetBranchAddress("true_beam_interactingEnergy", &true_beam_interactingEnergy);
394  defaultTree->SetBranchAddress("true_beam_PDG", &true_beam_PDG);
395  defaultTree->SetBranchAddress("true_beam_endZ", &true_beam_endZ);
396  defaultTree->SetBranchAddress("true_beam_endProcess", &true_beam_endProcess);
397 
398  defaultTree->SetBranchAddress("true_chexSignal", &true_chexSignal);
399  defaultTree->SetBranchAddress("true_nPi0Signal", &true_nPi0Signal);
400  defaultTree->SetBranchAddress("true_absSignal", &true_absSignal);
401  defaultTree->SetBranchAddress("true_backGround", &true_backGround);
402 
403  defaultTree->SetBranchAddress("event", &event);
404  defaultTree->SetBranchAddress("run", &run);
405 
406  int interaction_topology;
407  defaultTree->SetBranchAddress("interaction_topology", &interaction_topology);
408 
409  //Testing
410  std::vector<double> * reco_beam_calibrated_dEdX = 0x0;
411  defaultTree->SetBranchAddress("reco_beam_calibrated_dEdX", &reco_beam_calibrated_dEdX);
412  std::vector<double> * reco_beam_TrkPitch = 0x0;
413  defaultTree->SetBranchAddress("reco_beam_TrkPitch", &reco_beam_TrkPitch);
414 
415  /////////////////////////////////
416 
417 
418  int true_daughter_nPiPlus, true_daughter_nPiMinus, true_daughter_nPi0, true_beam_ID;
419  std::vector< std::string > * true_beam_processes = 0x0;
420  std::vector< int > * reco_beam_hit_true_ID = 0x0;
421  //double true_beam_interactingEnergy;
422  defaultTree->SetBranchAddress( "true_daughter_nPiPlus", &true_daughter_nPiPlus );
423  defaultTree->SetBranchAddress( "true_daughter_nPiMinus", &true_daughter_nPiMinus );
424  defaultTree->SetBranchAddress( "true_daughter_nPi0", &true_daughter_nPi0 );
425  defaultTree->SetBranchAddress( "true_beam_ID", &true_beam_ID );
426  defaultTree->SetBranchAddress( "true_beam_processes", &true_beam_processes );
427  defaultTree->SetBranchAddress( "reco_beam_hit_true_ID", &reco_beam_hit_true_ID );
428  //defaultTree->SetBranchAddress("true_beam_interactingEnergy", &true_beam_interactingEnergy);
429 
430 
431  //std::vector< int > * reco_beam_vertex_hits_slices = 0x0;
432  //std::vector< std::vector< double > > * reco_beam_vertex_dRs = 0x0;
433  //defaultTree->SetBranchAddress( "reco_beam_vertex_hits_slices", &reco_beam_vertex_hits_slices );
434  //defaultTree->SetBranchAddress( "reco_beam_vertex_dRs", &reco_beam_vertex_dRs );
435 
436  //For systs
437  std::vector<double> * g4rw_primary_plus_sigma_weight = 0x0;
438  std::vector<double> * g4rw_primary_minus_sigma_weight = 0x0;
439  std::vector<std::string> * g4rw_primary_var = 0x0;
440  if (doSyst != 0) {
441  defaultTree->SetBranchAddress("g4rw_primary_plus_sigma_weight", &g4rw_primary_plus_sigma_weight);
442  defaultTree->SetBranchAddress("g4rw_primary_minus_sigma_weight", &g4rw_primary_minus_sigma_weight);
443  defaultTree->SetBranchAddress("g4rw_primary_var", &g4rw_primary_var);
444  }
445 
446  channel.erase(std::remove(channel.begin(), channel.end(), '.'), channel.end());
447  channel.erase(std::remove(channel.begin(), channel.end(), ' '), channel.end());
448  topo.erase(std::remove(topo.begin(), topo.end(), '.'), topo.end());
449  topo.erase(std::remove(topo.begin(), topo.end(), ' '), topo.end());
450 
451  const int nrecobins = recoBins.size();
452  TString hist_name = Form("MC_Channel%s_%s_%.1f-%.1f_Histo",
453  channel.c_str(), topo.c_str(), minval, maxval);
454  TString hist_title = Form("MC Signal for channel %s and topology %s and true region %.1f-%.1f",
455  channel.c_str(), topo.c_str(), minval, maxval);
456 
457  if (doSyst == 1) {
458  hist_name += "_high_" + systName;
459  hist_title += " +1 sigma";
460  }
461  else if (doSyst == -1) {
462  hist_name += "_low_" + systName;
463  hist_title += " -1 sigma";
464  }
465 
466  TH1D* mchisto = new TH1D(hist_name, hist_title,
467  (doNegativeReco ? nrecobins : (nrecobins-1)),
468  (doNegativeReco ? -1 : 0), nrecobins-1);
469  mchisto->SetDirectory(0);
470 
471  mf::LogInfo("FillMCSignalHistogram_Pions") << "Filling MC signal histogram " <<
472  mchisto->GetName() << " from file " << filename.c_str() <<
473  " for channel " << channel.c_str() << " with topology " <<
474  topo.c_str() << " in the true region " << minval << "-" << maxval;
475 
476  bool done_check = false;
477  for (int k=0; k < defaultTree->GetEntries(); k++) {
478 
479  defaultTree->GetEntry(k);
480 
481  if (!done_check && doSyst != 0) {
482  if (g4rw_primary_var->size() > 0) {
483  auto syst_check = std::find(g4rw_primary_var->begin(), g4rw_primary_var->end(), systName);
484  if (syst_check == g4rw_primary_var->end()) {
485  std::cout << "Error! Could not find syst named " << systName << std::endl;
486  std::exception e;
487  throw e;
488  }
489  done_check = true;
490  }
491  }
492 
493  double syst_weight = 1.;
494  std::map<std::string, double> weights;
495  if (doSyst == 1) {
496  for (size_t i = 0; i < g4rw_primary_var->size(); ++i) {
497  weights[(*g4rw_primary_var)[i]] = (*g4rw_primary_plus_sigma_weight)[i];
498  }
499  }
500  else if (doSyst == -1) {
501  for (size_t i = 0; i < g4rw_primary_var->size(); ++i) {
502  weights[(*g4rw_primary_var)[i]] = (*g4rw_primary_minus_sigma_weight)[i];
503  }
504  }
505 
506  if ( reco_beam_true_byHits_origin == 2 )
507  continue;
508 
509 /*
510  int topology = -1;
511  //First determine if this is a pion ending in abs/cex/nPi0
512  if ( true_beam_PDG == 211 && *true_beam_endProcess == "pi+Inelastic" &&
513  true_daughter_nPiPlus == 0 && true_daughter_nPiMinus == 0 ){
514 
515  //make sure there is any hits at all
516  if ( !reco_beam_hit_true_ID->size() ) {
517  //std::cout << "no true_ID for hits " << reco_beam_interactingEnergy <<
518  // " " << reco_beam_calibrated_dEdX->size() << " " <<
519  // reco_beam_TrkPitch->size() << std::endl;
520  continue;
521  }
522 
523  //Instead of true_beam_endZ < 0. --> no true_beam_slices?
524  if (true_beam_endZ < 0. || true_beam_endZ > endZ_cut) continue;
525 
526  if (doSyst == 1 || doSyst == -1 ) {
527  syst_weight = weights[systName];
528  }
529 
530  // Absorption
531  if ( true_daughter_nPi0 == 0 )
532  topology = 1;
533  // Charge Exchange + nPi0
534  else
535  topology = 2;
536  }
537  else continue;
538  */
539 
540 
541  if (interaction_topology == -1)
542  std::cout << "Warning: Interaction topo -1" << std::endl;
543 
544  if (interaction_topology == 1 || interaction_topology ==2) {
545  //make sure there is any hits at all
546  if ( !reco_beam_hit_true_ID->size() ) {
547  //std::cout << "no true_ID for hits " << reco_beam_interactingEnergy <<
548  // " " << reco_beam_calibrated_dEdX->size() << " " <<
549  // reco_beam_TrkPitch->size() << std::endl;
550  continue;
551  }
552  if (doSyst == 1 || doSyst == -1 ) {
553  syst_weight = weights[systName];
554  }
555  }
556  else {
557  continue;
558  }
559 
560  // Remove events with the wrong topology
561  //if (topology != interaction_topology) {
562  // std::cout << "ERROR topo mismatch " << topology << " " << interaction_topology << std::endl;
563  //}
564  //if(topology != toponum) continue;
565  if(interaction_topology != toponum) continue;
566 
567  double total_weight = weight*syst_weight;
568  if (true_beam_PDG == 211) {
569  total_weight *= PiMuScale.first;
570  }
571  else if (true_beam_PDG == -13) {
572  total_weight *= PiMuScale.second;
573  }
574 
575  // True energy bin
576  //if(true_beam_interactingEnergy < minval) continue;
577  //if(true_beam_interactingEnergy >= maxval) continue;
578  if (true_beam_interactingEnergy < minval ||
579  true_beam_interactingEnergy >= maxval) continue;
580 
581  if (doNegativeReco) {
582  if (reco_beam_interactingEnergy < 0.0){
583  mchisto->AddBinContent(1, total_weight);
584  }
585  else{
586  for (int l = 1; l < nrecobins; l++) {
587  if (reco_beam_interactingEnergy > recoBins[l-1] &&
588  reco_beam_interactingEnergy <= recoBins[l]) {
589  mchisto->AddBinContent(l+1, total_weight);
590  break;
591  }
592  }
593  }
594  }
595  else {
596  // Sometimes the reco energy at vertex is mis-reconstructed
597  if (reco_beam_interactingEnergy < 0.0) continue;
598 
599  // Fill histogram in reco energy
600  for (int l = 1; l < nrecobins; l++) {
601  if (reco_beam_interactingEnergy > recoBins[l-1] &&
602  reco_beam_interactingEnergy <= recoBins[l]) {
603  mchisto->AddBinContent(l, total_weight);
604  break;
605  }
606  }
607  }
608  }
609 
610  file->Close();
611 
612  return mchisto;
613 
614 }
python.larbatch_posix.remove
def remove(path)
Definition: larbatch_posix.py:715
string
std::string string
Definition: nybbler.cc:12
mf::LogInfo
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Definition: MessageLogger.h:209
channel
uint8_t channel
Definition: CRTFragment.hh:201
train.filename
string filename
Definition: train.py:213
l
static QStrList * l
Definition: config.cpp:1044
test.weight
weight
Definition: test.py:257
e
const double e
Definition: gUpMuFluxGen.cxx:165
cache_state.file
file
Definition: cache_state.py:421
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
TH1 * protoana::ProtoDUNESelectionUtils::FillMCTruthSignalHistogram_Pions ( std::string  filename,
std::string  treename,
std::vector< double >  truthBins,
std::string  channel,
double  weight = 1.0 
)

Definition at line 1122 of file ProtoDUNESelectionUtils.cxx.

1122  {
1123  //********************************************************************
1124 
1125  TFile *file = new TFile(filename.c_str(), "READ");
1126  TTree *truthTree = (TTree*)file->Get(treename.c_str());
1127 
1128  Int_t true_beam_PDG;
1129  Double_t true_beam_interactingEnergy, true_beam_endZ;
1130  std::string *true_beam_endProcess = 0;
1131 
1132  truthTree->SetBranchAddress("true_beam_interactingEnergy", &true_beam_interactingEnergy);
1133  truthTree->SetBranchAddress("true_beam_PDG", &true_beam_PDG);
1134  truthTree->SetBranchAddress("true_beam_endZ", &true_beam_endZ);
1135  truthTree->SetBranchAddress("true_beam_endProcess", &true_beam_endProcess);
1136 
1137  channel.erase(std::remove(channel.begin(), channel.end(), '.'), channel.end());
1138  channel.erase(std::remove(channel.begin(), channel.end(), ' '), channel.end());
1139 
1140  const int ntruthbins = truthBins.size();
1141  TH1D* mchisto = new TH1D(Form("MC_Channel%s_TruthSig_Histo",channel.c_str()), Form("MC Truth Signal for channel %s",channel.c_str()), ntruthbins-1, 0, ntruthbins-1);
1142  mchisto->SetDirectory(0);
1143 
1144  mf::LogInfo("FillMCTruthSignalHistogram_Pions") << "Filling MC truth signal histogram " << mchisto->GetName() << " from file " << filename.c_str() << " for channel " << channel.c_str();
1145 
1146  for(Int_t k=0; k < truthTree->GetEntries(); k++){
1147  truthTree->GetEntry(k);
1148 
1149  // Pion beam
1150  if(true_beam_PDG != 211) continue;
1151 
1152  // If the pion does not make at the front face of the TPC, then don't count. Probably need a better estimate.
1153  if(true_beam_endZ < 0.0) continue;
1154 
1155  // True energy bin
1156  //if(true_beam_interactingEnergy < minval) continue;
1157  //if(true_beam_interactingEnergy >= maxval) continue;
1158 
1159  for(Int_t l = 1; l < ntruthbins; l++){
1160  if(true_beam_interactingEnergy > truthBins[l-1] && true_beam_interactingEnergy <= truthBins[l]){
1161  mchisto->AddBinContent(l, weight);
1162  break;
1163  }
1164  }
1165  }
1166 
1167  file->Close();
1168 
1169  return mchisto;
1170 }
python.larbatch_posix.remove
def remove(path)
Definition: larbatch_posix.py:715
string
std::string string
Definition: nybbler.cc:12
mf::LogInfo
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Definition: MessageLogger.h:209
channel
uint8_t channel
Definition: CRTFragment.hh:201
train.filename
string filename
Definition: train.py:213
l
static QStrList * l
Definition: config.cpp:1044
test.weight
weight
Definition: test.py:257
cache_state.file
file
Definition: cache_state.py:421
int protoana::ProtoDUNESelectionUtils::FindBin ( double  value,
std::vector< double >  binning 
)

Definition at line 1826 of file ProtoDUNESelectionUtils.cxx.

1826  {
1827  if (binning.size() < 2)
1828  return -1;
1829 
1830  for (size_t i = 1; i < binning.size(); ++i) {
1831  if ((binning[i-1] < value) && (value < binning[i])) {
1832  return i;
1833  }
1834  }
1835 
1836  return -1;
1837 }
submit_mcc.value
value
Definition: submit_mcc.py:159
std::pair< TH1 *, TH1 * > protoana::ProtoDUNESelectionUtils::GetMCIncidentEfficiency ( std::string  fileName,
std::string  treeName,
std::vector< double >  bins,
double  reco_beam_endZ_cut,
int  doSyst = 0,
double  weight = 1. 
)

Definition at line 1319 of file ProtoDUNESelectionUtils.cxx.

1322  {
1323 //********************************************************************
1324 
1325  TFile * file = new TFile(fileName.c_str(), "READ");
1326  TTree * defaultTree = (TTree*)file->Get(treeName.c_str());
1327 
1328 
1329  const size_t nBins = bins.size();
1330  TString hist_name = "MC_Incident_Efficiency_Denominator";
1331  TString hist_title = "Incident MC Efficiency Denominator";
1332 
1333  if (doSyst == 1) {
1334  hist_name += "_high";
1335  hist_title += " +1 sigma";
1336  }
1337  else if (doSyst == -1) {
1338  hist_name += "_low";
1339  hist_title += " -1 sigma";
1340  }
1341 
1342  TH1D * denominator = new TH1D(hist_name, hist_title, nBins-1, 0, nBins-1 );
1343  denominator->SetDirectory(0);
1344 
1345  hist_name = "MC_Incident_Efficiency_Numerator";
1346  hist_title = "Incident MC Efficiency Numerator";
1347 
1348  if (doSyst == 1) {
1349  hist_name += "_high";
1350  hist_title += " +1 sigma";
1351  }
1352  else if (doSyst == -1) {
1353  hist_name += "_low";
1354  hist_title += " -1 sigma";
1355  }
1356 
1357  TH1D * numerator = new TH1D(hist_name, hist_title, nBins-1, 0, nBins-1 );
1358  numerator->SetDirectory(0);
1359 
1360 
1361  double reco_beam_interactingEnergy, true_beam_endZ;
1362  int true_beam_ID, true_beam_PDG;
1363  std::vector< int > * reco_beam_hit_true_ID = 0x0;
1364  std::vector< int > * reco_beam_hit_true_slice = 0x0;
1365 
1366  std::vector< std::string > * true_beam_processes = 0x0;
1367  std::vector< int > * true_beam_slices = 0x0;
1368  std::vector< double > * true_beam_incidentEnergies = 0x0;
1369 
1370  defaultTree->SetBranchAddress( "reco_beam_interactingEnergy", &reco_beam_interactingEnergy );
1371  defaultTree->SetBranchAddress( "reco_beam_hit_true_ID", &reco_beam_hit_true_ID );
1372  defaultTree->SetBranchAddress( "reco_beam_hit_true_slice", &reco_beam_hit_true_slice );
1373  defaultTree->SetBranchAddress( "true_beam_processes", &true_beam_processes );
1374  defaultTree->SetBranchAddress( "true_beam_slices", &true_beam_slices );
1375  defaultTree->SetBranchAddress( "true_beam_ID", &true_beam_ID );
1376  defaultTree->SetBranchAddress( "true_beam_PDG", &true_beam_PDG );
1377  defaultTree->SetBranchAddress( "true_beam_endZ", &true_beam_endZ );
1378  defaultTree->SetBranchAddress( "true_beam_incidentEnergies", &true_beam_incidentEnergies );
1379 
1380  //std::vector< int > * reco_beam_vertex_hits_slices = 0x0;
1381  //std::vector< std::vector< double > > * reco_beam_vertex_dRs = 0x0;
1382  //defaultTree->SetBranchAddress( "reco_beam_vertex_hits_slices", &reco_beam_vertex_hits_slices );
1383  //defaultTree->SetBranchAddress( "reco_beam_vertex_dRs", &reco_beam_vertex_dRs );
1384 
1385  double g4rw_primary_plus_sigma_weight, g4rw_primary_minus_sigma_weight;
1386  if (doSyst != 0) {
1387  defaultTree->SetBranchAddress("g4rw_primary_plus_sigma_weight", &g4rw_primary_plus_sigma_weight);
1388  defaultTree->SetBranchAddress("g4rw_primary_minus_sigma_weight", &g4rw_primary_minus_sigma_weight);
1389  }
1390 
1391  for( int k=0; k < defaultTree->GetEntries(); ++k ){
1392  double syst_weight = 1.;
1393  defaultTree->GetEntry(k);
1394 
1395  if ( true_beam_PDG != 211 )
1396  continue;
1397 
1398  if (doSyst == 1) {
1399  syst_weight = g4rw_primary_plus_sigma_weight;
1400  }
1401  else if (doSyst == -1) {
1402  syst_weight = g4rw_primary_minus_sigma_weight;
1403  }
1404 
1405  //if ( true_beam_endZ < 0. && true_beam_slices->size() )
1406  // std::cout << "NOTICE: endZ < 0. but has true slices "
1407  // << true_beam_slices->size() << std::endl;
1408 
1409  //if ( true_beam_endZ > 225. )
1410  // continue;
1411 
1412  if ((true_beam_incidentEnergies->size() && !true_beam_slices->size()) ||
1413  (!true_beam_incidentEnergies->size() && true_beam_slices->size())) {
1414  std::cout << "NOTICE! Energies: " <<
1415  true_beam_incidentEnergies->size() << " slices: " <<
1416  true_beam_slices->size() << std::endl;
1417  }
1418 
1419  double pitch = 0.4792;
1420  double z0 = 0.56035;
1421  int slice_cut = std::floor((reco_beam_endZ_cut - (z0 - pitch/2.)) / pitch);
1422 
1423  std::vector< int > denom_slices;
1424 
1425  //First, fill the denominator histogram
1426  //
1427  for (size_t i = 0; i < true_beam_incidentEnergies->size(); ++i) {
1428 
1429  int the_slice = (*true_beam_slices)[i];
1430  if (the_slice > slice_cut) continue;
1431 
1432  for ( size_t j = 1; j < nBins; ++j ) {
1433  if ( true_beam_incidentEnergies->at(i) > bins[j-1] &&
1434  true_beam_incidentEnergies->at(i) <= bins[j] ) {
1435  denominator->AddBinContent(j, weight*syst_weight);
1436  break;
1437  }
1438  }
1439  denom_slices.push_back(the_slice);
1440  }
1441 
1442 
1443  // Go through the incident energies from the beam particle
1444  for ( size_t i = 0; i < reco_beam_hit_true_ID->size(); ++i ) {
1445 
1446  // Check the ID, origin, true_slice
1447  int true_id = (*reco_beam_hit_true_ID)[i];
1448  int true_slice = (*reco_beam_hit_true_slice)[i];
1449 
1450  if ( true_id == true_beam_ID ) {
1451  if (true_slice > slice_cut) continue;
1452 
1453  //search through the true beam slices
1454  for ( size_t j = 0; j < true_beam_slices->size(); ++j ) {
1455  if ( true_slice == (*true_beam_slices)[j] ) {
1456 
1457  auto slice_check = std::find(denom_slices.begin(),
1458  denom_slices.end(),
1459  true_slice);
1460 
1461  if ( slice_check == denom_slices.end() ) {
1462  std::cout << "Warning found true slice from hit"
1463  << "that was not in denominator" << std::endl;
1464  }
1465  for ( size_t m = 1; m < nBins; ++m ) {
1466  if ( (*true_beam_incidentEnergies)[j] > bins[m-1]
1467  && (*true_beam_incidentEnergies)[j] <= bins[m] ) {
1468  numerator->AddBinContent(m, weight*syst_weight);
1469  break;
1470  }
1471  }
1472  }
1473  }
1474  }
1475  }
1476  }
1477 
1478  for (size_t i = 1; i < bins.size(); ++i) {
1479  TString label = Form("%.1f-%.1f", bins[i-1], bins[i]);
1480  numerator->GetXaxis()->SetBinLabel(i, label.Data());
1481  denominator->GetXaxis()->SetBinLabel(i, label.Data());
1482  }
1483  numerator->SetTitle(";E_{True} (MeV)");
1484  denominator->SetTitle(";E_{True} (MeV)");
1485 
1486  file->Close();
1487 
1488  return {numerator, denominator};
1489 }
test_nxdot.label
label
Definition: test_nxdot.py:6
test.weight
weight
Definition: test.py:257
dumpTree.fileName
fileName
Definition: dumpTree.py:9
cache_state.file
file
Definition: cache_state.py:421
make_THn_beam_input.nBins
nBins
Definition: make_THn_beam_input.py:20
cache_state.m
m
Definition: cache_state.py:415
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
std::pair< TH1 *, TH1 * > protoana::ProtoDUNESelectionUtils::GetMCInteractingEfficiency ( std::string  fileName,
std::string  treeName,
std::vector< double >  bins,
std::string  channel,
std::string  topo,
int  toponum,
double  endZ_cut,
int  doSyst = 0,
double  weight = 1. 
)

Definition at line 1493 of file ProtoDUNESelectionUtils.cxx.

1497  {
1498 //********************************************************************
1499 
1500  TFile * file = new TFile(fileName.c_str(), "READ");
1501  TTree * defaultTree = (TTree*)file->Get(treeName.c_str());
1502 
1503 
1504  const size_t nBins = bins.size();
1505  TString hist_name = Form("MC_Channel%s_%s_Interacting_Denominator",channel.c_str(),topo.c_str());
1506  TString hist_title = Form( "Interacting MC Efficiency Denominator for channel %s and topology %s", channel.c_str(), topo.c_str());
1507 
1508  if (doSyst == 1) {
1509  hist_name += "_high";
1510  hist_title += " +1 sigma";
1511  }
1512  else if (doSyst == -1) {
1513  hist_name += "_low";
1514  hist_title += " -1 sigma";
1515  }
1516 
1517  TH1D * denominator = new TH1D(hist_name, hist_title, nBins-1, 0, nBins-1);
1518  denominator->SetDirectory(0);
1519 
1520  hist_name = Form("MC_Channel%s_%s_Interacting_Numerator",channel.c_str(),topo.c_str());
1521  hist_title = Form("Interacting MC Efficiency Numerator for channel %s and topology %s", channel.c_str(), topo.c_str());
1522 
1523  if (doSyst == 1) {
1524  hist_name += "_high";
1525  hist_title += " +1 sigma";
1526  }
1527  else if (doSyst == -1) {
1528  hist_name += "_low";
1529  hist_title += " -1 sigma";
1530  }
1531 
1532 
1533  TH1D * numerator = new TH1D(hist_name, hist_title, nBins-1, 0, nBins-1 );
1534  numerator->SetDirectory(0);
1535 
1536  //Init stuff here
1537  int true_beam_PDG, true_daughter_nPiPlus, true_daughter_nPiMinus, true_daughter_nPi0;
1538  std::string * true_beam_endProcess = 0x0;
1539  std::vector< std::string > * true_beam_processes = 0x0;
1540  double true_beam_interactingEnergy, true_beam_endZ;
1541  bool has_noPion_daughter, has_shower_nHits_distance;
1542 
1543  std::vector< int > * true_beam_slices = 0x0;
1544  defaultTree->SetBranchAddress( "true_beam_slices", &true_beam_slices );
1545 
1546  defaultTree->SetBranchAddress( "true_beam_PDG", &true_beam_PDG );
1547  defaultTree->SetBranchAddress( "true_beam_endZ", &true_beam_endZ );
1548  defaultTree->SetBranchAddress( "true_beam_endProcess", &true_beam_endProcess );
1549  defaultTree->SetBranchAddress( "true_daughter_nPiPlus", &true_daughter_nPiPlus );
1550  defaultTree->SetBranchAddress( "true_daughter_nPiMinus", &true_daughter_nPiMinus );
1551  defaultTree->SetBranchAddress( "true_daughter_nPi0", &true_daughter_nPi0 );
1552  defaultTree->SetBranchAddress( "true_beam_processes", &true_beam_processes );
1553  defaultTree->SetBranchAddress( "true_beam_interactingEnergy", &true_beam_interactingEnergy );
1554  defaultTree->SetBranchAddress( "has_noPion_daughter", &has_noPion_daughter );
1555  defaultTree->SetBranchAddress( "has_shower_nHits_distance", &has_shower_nHits_distance );
1556 
1557  //std::vector< int > * reco_beam_vertex_hits_slices = 0x0;
1558  //std::vector< std::vector< double > > * reco_beam_vertex_dRs = 0x0;
1559  //defaultTree->SetBranchAddress( "reco_beam_vertex_hits_slices", &reco_beam_vertex_hits_slices );
1560  //defaultTree->SetBranchAddress( "reco_beam_vertex_dRs", &reco_beam_vertex_dRs );
1561 
1562  double g4rw_primary_plus_sigma_weight, g4rw_primary_minus_sigma_weight;
1563  if (doSyst != 0) {
1564  defaultTree->SetBranchAddress("g4rw_primary_plus_sigma_weight", &g4rw_primary_plus_sigma_weight);
1565  defaultTree->SetBranchAddress("g4rw_primary_minus_sigma_weight", &g4rw_primary_minus_sigma_weight);
1566  }
1567 
1568  int topology = -1;
1569 
1570  for( int k=0; k < defaultTree->GetEntries(); ++k ){
1571  double syst_weight = 1.;
1572  defaultTree->GetEntry(k);
1573 
1574  //First, check if it's a pion
1575  if ( true_beam_PDG != 211 )
1576  continue;
1577 
1578  if (!true_beam_slices->size())//( true_beam_endZ < 0. )
1579  continue;
1580 
1581  //Then if it ends in abs/cex/nPi0 in the FV
1582  if (!(*true_beam_endProcess == "pi+Inelastic" &&
1583  true_daughter_nPiPlus == 0 && true_daughter_nPiMinus == 0 &&
1584  true_beam_endZ < endZ_cut)){
1585  continue;
1586  }
1587 
1588  //Abs
1589  if (true_daughter_nPi0 == 0) {
1590  topology = 1;
1591  }
1592  //Cex/nPi0
1593  else {
1594  topology = 2;
1595  }
1596 
1597  if (topology != toponum){
1598  continue;
1599  }
1600 
1601  if (doSyst == 1) {
1602  syst_weight = g4rw_primary_plus_sigma_weight;
1603  }
1604  else if (doSyst == -1) {
1605  syst_weight = g4rw_primary_minus_sigma_weight;
1606  }
1607 
1608  //Fill the denominator
1609  for ( size_t i = 1; i < nBins; ++i ) {
1610  if (true_beam_interactingEnergy > bins[i-1]
1611  && true_beam_interactingEnergy <= bins[i]){
1612  denominator->AddBinContent(i, weight*syst_weight);
1613 
1614  if ( has_noPion_daughter ) { //Abs/Cex selection
1615  if (topology == 1 && !has_shower_nHits_distance) {
1616  numerator->AddBinContent(i, weight*syst_weight); //Abs
1617  }
1618  else if (topology == 2 && has_shower_nHits_distance) {
1619  numerator->AddBinContent(i, weight*syst_weight); //Cex
1620  }
1621  }
1622  break;
1623 
1624  }
1625  }
1626  }
1627 
1628  for (size_t i = 1; i < bins.size(); ++i) {
1629  TString label = Form("%.1f-%.1f", bins[i-1], bins[i]);
1630  numerator->GetXaxis()->SetBinLabel(i, label.Data());
1631  denominator->GetXaxis()->SetBinLabel(i, label.Data());
1632  }
1633 
1634  numerator->SetTitle(";E_{True} (MeV)");
1635  denominator->SetTitle(";E_{True} (MeV)");
1636 
1637  file->Close();
1638  return {numerator, denominator};
1639 }
string
std::string string
Definition: nybbler.cc:12
test_nxdot.label
label
Definition: test_nxdot.py:6
channel
uint8_t channel
Definition: CRTFragment.hh:201
test.weight
weight
Definition: test.py:257
dumpTree.fileName
fileName
Definition: dumpTree.py:9
cache_state.file
file
Definition: cache_state.py:421
make_THn_beam_input.nBins
nBins
Definition: make_THn_beam_input.py:20
int protoana::ProtoDUNESelectionUtils::GetNTriggers_Pions ( std::string  filename,
std::string  treename,
bool  IsMC = true 
)

Definition at line 1291 of file ProtoDUNESelectionUtils.cxx.

1291  {
1292 //********************************************************************
1293 
1294  TFile *file = new TFile(filename.c_str(), "READ");
1295  TTree *defaultTree = (TTree*)file->Get(treename.c_str());
1296 
1297  int ntriggers = 0;
1298 
1299  Int_t true_beam_PDG, reco_beam_true_byHits_origin;
1300 
1301  defaultTree->SetBranchAddress("true_beam_PDG", &true_beam_PDG);
1302  defaultTree->SetBranchAddress("reco_beam_true_byHits_origin", &reco_beam_true_byHits_origin);
1303 
1304  for(Int_t k=0; k < defaultTree->GetEntries(); k++){
1305  defaultTree->GetEntry(k);
1306 
1307  if(true_beam_PDG != 211) continue; // should be removed when proper data input is used
1308  if(reco_beam_true_byHits_origin != 4) continue;
1309  ntriggers++;
1310 
1311  }
1312 
1313  return ntriggers;
1314 
1315 }
train.filename
string filename
Definition: train.py:213
cache_state.file
file
Definition: cache_state.py:421
Generated by   doxygen 1.8.11