MakeXSecCCPi0Inc_module.cc
Go to the documentation of this file.
10 #include "cetlib/search_path.h"
11 
12 #include "Utilities/AssociationUtil.h"
13 #include "Utilities/func/MathUtil.h"
15 #include "ReMId/ReMId.h"
16 #include "RecoBase/Cluster.h"
17 #include "RecoBase/Prong.h"
18 #include "RecoBase/RecoHit.h"
19 #include "RecoBase/Vertex.h"
20 #include "RecoBase/Track.h"
21 
22 #include "TVector3.h"
23 #include "TH2.h"
24 #include "TFile.h"
25 
26 namespace xsrec
27 {
29  {
30  template<class T> using Atom = fhicl::Atom<T>;
31  template<class T> using Table = fhicl::Table<T>;
33  using Name = fhicl::Name;
34 
38  Atom<std::string> prong3DAssn {Name("Prong3DAssnInput")};
41 
43  };
44 
45  // Stores template histograms
46  struct CCPi0LLInput
47  {
48  TH2F *fBPIEPH_phot;
49  TH2F *fBPIEPH_oth;
50  TH2F *fGapNMP_phot;
51  TH2F *fGapNMP_oth;
52  };
53 
54  /// Calculates a prong level PID level to select photons
56  {
57  public:
59 
60  explicit MakeXSecCCPi0Inc(const Parameters& params);
61  virtual ~MakeXSecCCPi0Inc();
62 
63  void produce (art::Event& evt);
64  void beginJob();
65 
66  protected:
67  /// Calc \c \#Delta LL score input
68  double GetBPI(art::Ptr<rb::Prong> prng);
69  /// Calc \c \#Delta LL
70  double GetPhLL(art::Ptr<rb::Prong> prng, double bpi, double gap);
71  /// Find prong corresponding to muon
72  int GetMuIdx(std::vector< art::Ptr<rb::Prong> > ps, TVector3 trkdir);
73 
75 
77 
78  };
79 
81  : fParams(params())
82  {
83  produces< std::vector< xsrec::XSecCCPi0Inc> >();
84  produces< art::Assns< xsrec::XSecCCPi0Inc, rb::Prong> >();
85  }
86 
88 
90  {
91  // Load in template hists for #Delta LL
92  cet::search_path sp("XSECCCPI0INC_LIB_PATH");
93  std::string CCPi0IDTemplateName;
94  if (!sp.find_file(fParams.llInput(),CCPi0IDTemplateName))
95  throw cet::exception("MakeXSecCCPi0Inc")
96  << fParams.llInput()
97  << "could not be located";
98 
99  TFile *tempfile = TFile::Open(CCPi0IDTemplateName.c_str());
100  tempfile->cd();
101  fCCPi0LLInput.fBPIEPH_phot = (TH2F*)tempfile->Get("bpiephsig");
102  fCCPi0LLInput.fBPIEPH_oth = (TH2F*)tempfile->Get("bpiephbac");
103  fCCPi0LLInput.fGapNMP_phot = (TH2F*)tempfile->Get("gapnmpsig");
104  fCCPi0LLInput.fGapNMP_oth = (TH2F*)tempfile->Get("gapnmpbac");
105  tempfile->Close();
106  }
107 
108  ////////////////////////////////////////////////////////////////////////////
109  /// Determine which prong is the remid track, the most colinear
111  TVector3 trkdir)
112  {
113  double maxcos = -2;
114  unsigned int muIdx = -1;
115  for (unsigned int pIdx = 0; pIdx < prngs.size(); pIdx++){
116  TVector3 dir = prngs[pIdx]->Dir().Unit();
117  // If current prong is more parallel to trk dir, update cos and muIdx
118  if (dir.Dot(trkdir) > maxcos){
119  maxcos = dir.Dot(trkdir);
120  muIdx = pIdx;
121  }
122  }
123  return muIdx;
124  }
125 
126  ////////////////////////////////////////////////////////////////////////////
127  /// Get a BPI value for each prong - input to \c \#Delta LL
129  {
130  TVector3 start = prng->Start();
131  // Stores (dist to prong start, GeV) for each hit
132  std::vector< std::pair<double,double> > hitinfo;
133  // Loop over RecoHits
134  for (int rhIdx = 0; rhIdx < int(prng->NCell()); rhIdx++){
135  // Need to grap RecoHits for GeV / (x,y,z) estimate
136  rb::RecoHit rhit = prng->RecoHit(rhIdx);
137  if (!rhit.IsCalibrated()) continue;
138  double gev = rhit.GeV();
139  double l = util::pythag(start[0]-rhit.X(),
140  start[1]-rhit.Y(),
141  start[2]-rhit.Z()); // Cur hit to prong start
142  hitinfo.push_back(std::make_pair(l,gev));
143  }
144  // Default comparison is on first element (l), with second tie-breaking
145  std::sort(hitinfo.begin(),hitinfo.end());
146 
147  if (hitinfo.size() < 4) return -1; // Not defined for fewer than 4 hits
148  // Choose how many hits to consider "the end" and "the bulk"
149  int NHit = hitinfo.size();
150  int NEnd = std::min(6,NHit/2);
151  int NBulk = hitinfo.size()-NEnd;
152 
153  double endEPH = 0; double bulkEPH = 0;
154  for (int hIdx = 0; hIdx < NHit; hIdx++){
155  if (hIdx >= NHit-NEnd) // we're within NEnd cells of furthest
156  endEPH += hitinfo[hIdx].second / NEnd;
157  else // we're in the first NBulk cells
158  bulkEPH += hitinfo[hIdx].second / NBulk;
159  }
160  return bulkEPH > 0 ? endEPH / bulkEPH : -1; // guard against / by 0
161  }
162 
163  ////////////////////////////////////////////////////////////////////////////
164  /// PhLL - the \c \#Delta LL prong-level score calculator
166  double bpi, double gap)
167  {
168  if (prng->NCell() < 4) return -20;
169  double eph = prng->CalorimetricEnergy()/prng->NCell();
170  int nmp = prng->MostMissingPlanes(geo::kXorY);
171 
172  // Get bin #'s for each var
173  int bpibin = fCCPi0LLInput.fBPIEPH_phot->GetXaxis()->FindBin(bpi);
174  int ephbin = fCCPi0LLInput.fBPIEPH_phot->GetYaxis()->FindBin(eph);
175  int gapbin = fCCPi0LLInput.fGapNMP_phot->GetXaxis()->FindBin(gap);
176  int nmpbin = fCCPi0LLInput.fGapNMP_phot->GetYaxis()->FindBin(nmp);
177  if (nmpbin > 10) nmpbin = 10; // Cap at max
178 
179  // Use input hists
180  double pphot = fCCPi0LLInput.fBPIEPH_phot->GetBinContent(bpibin,ephbin) *
181  fCCPi0LLInput.fGapNMP_phot->GetBinContent(gapbin,nmpbin);
182  double poth = fCCPi0LLInput.fBPIEPH_oth ->GetBinContent(bpibin,ephbin) *
183  fCCPi0LLInput.fGapNMP_oth ->GetBinContent(gapbin,nmpbin);
184 
185  // If we don't have any information, set to low number
186  if (pphot <= 0 || poth <= 0) return -20;
187 
188  // return #Delta LL
189  return std::log(pphot)-std::log(poth);
190  }
191 
192 
193 
194  ////////////////////////////////////////////////////////////////////////////
196  {
197  // We make a vector of XSecCCPi0Inc, associating each with a prong
198  std::unique_ptr< std::vector<XSecCCPi0Inc> >
199  xsecccpi0col(new std::vector<XSecCCPi0Inc>);
200  std::unique_ptr< art::Assns<XSecCCPi0Inc, rb::Prong> >
201  xsecccpi0assncol(new art::Assns<XSecCCPi0Inc, rb::Prong>);
202 
203  // Pull out slice info, and associated tracks / vertices
205  evt.getByLabel(fParams.slicerInput(), slcs);
206  art::FindManyP<rb::Vertex> fmVertex(slcs, evt, fParams.vertexInput());
207  art::FindManyP<rb::Track> fmTrack(slcs, evt, fParams.trackInput());
208 
209  // Loop over slices
210  for (unsigned int sIdx = 0; sIdx < slcs->size(); sIdx++){
211  art::Ptr<rb::Cluster> slPtr(slcs,sIdx);
212  // Find vec of vertices, 3D prongs, trks
213  std::vector< art::Ptr<rb::Vertex> > vertexs = fmVertex.at(sIdx);
214  if (vertexs.size() == 0) continue; // We only get prongs if we see vertex
216  std::vector< art::Ptr<rb::Prong> > prngs = fmProng3d.at(0); // 1st entry
217 
218  std::vector< art::Ptr<rb::Track> > trks = fmTrack.at(sIdx);
219 
220  // Works on a prong-by-prong level, need to have a prong
221  if (prngs.size() == 0)
222  continue; // Don't bother calculating without a prong
223 
224  // Now, find the prong most parallel to muon
225  art::FindManyP<remid::ReMId> fmReMId(trks, evt, fParams.remidInput());
226  unsigned int bestReMIdIdx =
228 
229  // Start setting up BPI calculator //
230  ////////////////////////////////////////////////////////////////////
231 
232 
233  // Set muIdx = -1 if there is no track / suitable remid
234  int muIdx = bestReMIdIdx >= trks.size() ?
235  -1 : GetMuIdx(prngs,trks[bestReMIdIdx]->Dir().Unit());
236 
237  // Need to calculate vtx-prong gap
238  TVector3 vtx = vertexs[0]->GetXYZ();
239  for (int i = 0; i < int(prngs.size()); i++){
240 
241  TVector3 start = prngs[i]->Start();
242  double gap = (vtx-start).Mag(); // #Delta LL input
243  double bpi = GetBPI(prngs[i]); // #Delta LL input
244  double phll = GetPhLL(prngs[i],bpi,gap); // Calc #Delta LL
245  XSecCCPi0Inc curPID;
246 
247  curPID.SetBPI(bpi); // Can't calc this input var in CAFs
248  curPID.SetPhLL(phll);
249  curPID.SetIsMuProng(i==muIdx); // For convenience
250 
251  // Add to list of prong PID's / associations
252  xsecccpi0col->emplace_back(curPID);
253  util::CreateAssn(*this,evt,*(xsecccpi0col.get()),
254  prngs[i],*(xsecccpi0assncol.get()));
255 
256  } // End loop over prongs
257 
258  } // End loop over slices
259  evt.put(std::move(xsecccpi0col));
260  evt.put(std::move(xsecccpi0assncol));
261  } // End produce
262 
263 } // End namespace
264 
265 namespace xsrec
266 {
268 }
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
Calculates a prong level PID level to select photons.
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
X or Y views.
Definition: PlaneGeo.h:30
float Z() const
Definition: RecoHit.h:38
void SetBPI(double bpi)
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
std::string find_file(std::string const &filename) const
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
MakeXSecCCPi0Inc(const Parameters &params)
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
int MostMissingPlanes(geo::View_t view) const
Longest run of adjacent planes with no hits.
Definition: Cluster.cxx:668
double CalorimetricEnergy(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple estimate of neutrino energy.
Definition: Cluster.cxx:439
void beginJob()
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Definition: Cluster.cxx:259
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
double GetPhLL(art::Ptr< rb::Prong > prng, double bpi, double gap)
Calc #Delta LL.
TVector3 Unit() const
Vertex location in position and time.
double GetBPI(art::Ptr< rb::Prong > prng)
Calc #Delta LL score input.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
float GeV() const
Definition: RecoHit.cxx:69
TDirectory * dir
Definition: macro.C:5
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
void SetIsMuProng(bool isMuProng)
float X() const
Definition: RecoHit.h:36
float Mag() const
int GetMuIdx(std::vector< art::Ptr< rb::Prong > > ps, TVector3 trkdir)
Find prong corresponding to muon.
void SetPhLL(double phll)
float Y() const
Definition: RecoHit.h:37
T min(const caf::Proxy< T > &a, T b)
Definition: fwd.h:28
unsigned int HighestPIDTrack(const std::vector< art::Ptr< rb::Track > > &sliceTracks, const std::string &remidModuleLabel, const art::Event &e)
Definition: ReMId.cxx:249