MakeNCCosRej_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: MakeNCCosRej
3 // Module Type: producer
4 // \file MakeNCCosRej_module.cc
5 // \brief A module to write out the NCCosRej object
6 //
7 // \author Shaokai Yang - yangs9uc@gmail.com
8 ////////////////////////////////////////////////////////////////////////
9 // system includes
10 #include <memory>
11 #include <string>
12 #include <vector>
13 
14 // art includes
23 #include "cetlib_except/exception.h"
24 #include "cetlib/filesystem.h"
25 #include "cetlib/search_path.h"
26 #include "fhiclcpp/ParameterSet.h"
28 
29 // ROOT includes
30 #include "TMath.h"
31 #include "TMVA/Reader.h"
32 #include "TVector3.h"
33 
34 // novasoft includes
35 #include "CosRej/NueCosRej.h"
36 #include "CVN/func/Result.h"
37 #include "Geometry/Geometry.h"
38 #include "MCCheater/BackTracker.h"
39 #include "NCID/NCCosRej.h"
41 #include "ShowerLID/ShowerLID.h"
42 #include "RecoBase/CellHit.h"
43 #include "RecoBase/Cluster.h"
44 #include "RecoBase/Prong.h"
45 #include "RecoBase/RecoHit.h"
46 #include "RecoBase/Shower.h"
48 #include "Utilities/AssociationUtil.h"
49 
50 
51 namespace ncid
52 {
54  {
55  public:
56 
57  explicit MakeNCCosRej(fhicl::ParameterSet const & p);
58  ~MakeNCCosRej();
59 
60  void beginRun(art::Run &evt) override;
61  void reconfigure(const fhicl::ParameterSet& p);
62  void produce(art::Event & e) override;
63 
64  private:
65 
78 
79  TMVA::Reader* fReaderDT;
80  TMVA::Reader* fReaderGBDT;
81 
82  // Variables need to send to TMVA::Reader
83  float TMVAvars[13];
84 
85  /// Particle ID alorithm (loglikelihoods and dE/dx)
86  //slid::ParticleIDAlg* fParticleIDAlg;
87  /// FCL parameters for particle ID alorithm (loglikelihoods and dE/dx)
89 
90  /// Return transverse momentum fraction of the event.
91  /// Calculation based on reconstructed prongs.
92  double TransMomFraction(const std::vector
93  < art::Ptr<rb::Prong> >& prongs) const;
94  };
95 
97  const art::Ptr<rb::Prong> b);
98 }
99 
100 //////////////////////////////////////////////////////////////////////////
101 
102 namespace ncid
103 {
105  fReaderDT(nullptr),
106  fReaderGBDT(nullptr),
107  //fParticleIDAlg(0),
108  fParticleIDAlgPSet(p.get< fhicl::ParameterSet >("ParticleIDAlgPSet"))
109  {
110  produces< std::vector<ncid::NCCosRej> > ();
111  produces< art::Assns<ncid::NCCosRej, rb::Cluster> > ();
112  this->reconfigure(p);
113  }
114 
116  {
117  if(fReaderDT) delete fReaderDT;
118  if(fReaderGBDT) delete fReaderGBDT;
119  }
120 
122  {
124 
125  // Make sure we can find the Real Adaptive Boosted Decision Trees weight file before we set up MVA
126  if(!cet::file_exists(pidlibpath+fNCCosRejPIDFile)){
127  throw cet::exception("MakeNCCosRej")
128  << "Couldn't find file " << pidlibpath+fNCCosRejPIDFile <<std::endl;
129  }
130  // Make sure we can find the Gradient Boosted Decision Trees weight file before we set up MVA
131  if(!cet::file_exists(pidlibpath+fNCCosRejPIDGFile)){
132  throw cet::exception("MakeNCCosRej")
133  << "Couldn't find file " << pidlibpath+fNCCosRejPIDGFile <<std::endl;
134  }
135 
137 
138  fReaderDT = new TMVA::Reader;
139  fReaderDT->AddVariable("cosmicid", &TMVAvars[0]); //cosmicid
140  fReaderDT->AddVariable("partptp", &TMVAvars[1]); //partptp
141  fReaderDT->AddVariable("shwnhit", &TMVAvars[2]); //shwnhit
142  fReaderDT->AddVariable("shwxminusy", &TMVAvars[3]); //shwminusy
143  fReaderDT->AddVariable("shwxplusy", &TMVAvars[4]); //shwxplusy
144  fReaderDT->AddVariable("shwxovery", &TMVAvars[5]); //shwxpvery
145  fReaderDT->AddVariable("shwcalE", &TMVAvars[6]); //shwcalE
146  fReaderDT->AddVariable("shwdirY", &TMVAvars[7]); //shwdirY
147  fReaderDT->AddVariable("shwlen", &TMVAvars[8]); //shwlen
148  fReaderDT->AddVariable("shwwwidth", &TMVAvars[9]); //shwwwidth
149  fReaderDT->AddVariable("shwGap", &TMVAvars[10]); //shwGap
150  fReaderDT->AddVariable("nshwlid", &TMVAvars[11]); //nshwlid
151  fReaderDT->AddVariable("nmiphit", &TMVAvars[12]); //nmiphit
152  fReaderDT->BookMVA("BDTA",pidlibpath+fNCCosRejPIDFile);
153 
154  fReaderGBDT = new TMVA::Reader;
155  fReaderGBDT->AddVariable("cosmicid", &TMVAvars[0]); //cosmicid
156  fReaderGBDT->AddVariable("partptp", &TMVAvars[1]); //partptp
157  fReaderGBDT->AddVariable("shwnhit", &TMVAvars[2]); //shwnhit
158  fReaderGBDT->AddVariable("shwxminusy", &TMVAvars[3]); //shwminusy
159  fReaderGBDT->AddVariable("shwxplusy", &TMVAvars[4]); //shwxplusy
160  fReaderGBDT->AddVariable("shwxovery", &TMVAvars[5]); //shwxpvery
161  fReaderGBDT->AddVariable("shwcalE", &TMVAvars[6]); //shwcalE
162  fReaderGBDT->AddVariable("shwdirY", &TMVAvars[7]); //shwdirY
163  fReaderGBDT->AddVariable("shwlen", &TMVAvars[8]); //shwlen
164  fReaderGBDT->AddVariable("shwwwidth", &TMVAvars[9]); //shwwwidth
165  fReaderGBDT->AddVariable("shwGap", &TMVAvars[10]); //shwGap
166  fReaderGBDT->AddVariable("nshwlid", &TMVAvars[11]); //nshwlid
167  fReaderGBDT->AddVariable("nmiphit", &TMVAvars[12]); //nmiphit
168  fReaderGBDT->BookMVA("BDTG", pidlibpath+fNCCosRejPIDGFile);
169 }
170 
172  {
173  //fMCTruthModuleLabel = p.get< std::string >("MCTruthModuleLabel");
174  fSliceLabel = p.get< std::string >("SliceLabel");
175  fCellHitLabel = p.get< std::string >("CellHitLabel");
176  fProngLabel = p.get< std::string >("ProngLabel");
177  fProngInstance = p.get< std::string >("ProngInstance");
178  fCVNLabel = p.get< std::string >("CVNLabel");
179  fShowerLabel = p.get< std::string >("ShowerLabel");
180  fShowerLIDLabel = p.get< std::string >("ShowerLIDLabel");
181  fPIDLibPath = p.get< std::string >("PIDLibPath");
182  fCosRejLabel = p.get< std::string >("CosRejLabel");
183  fNCCosRejPIDFile = p.get< std::string >("NCCosRejPIDFile");
184  fNCCosRejPIDGFile = p.get< std::string >("NCCosRejPIDGFile");
185  fParticleIDAlgPSet = p.get< fhicl::ParameterSet >("ParticleIDAlgPSet");
186  }
187 
188  //////////////////////////////////////////////////////////////////////////
190  {
193 
194  std::unique_ptr< std::vector<ncid::NCCosRej> >
195  ncCosRejjies(new std::vector< ncid::NCCosRej >);
196  std::unique_ptr< art::Assns<ncid::NCCosRej, rb::Cluster> >
198  // Get the slices in the event
200  evt.getByLabel( fSliceLabel, sliceHandle);
201 
202  std::vector< art::Ptr<rb::Cluster> > slices;
203  art::fill_ptr_vector(slices, sliceHandle);
204 
206  evt.getByLabel(fCellHitLabel, htHandle);
207 
209 
210  for(unsigned int m = 0; m < htHandle->size(); ++m){
211  allHits.push_back(art::Ptr<rb::CellHit>(htHandle,m));
212  }
213  // Get things associated with slices
214  art::FindManyP<rb::Prong> prongAssn(sliceHandle, evt, art::InputTag(fProngLabel, fProngInstance));
215  art::FindMany<slid::ShowerLID> showerLidAssn(sliceHandle, evt, fShowerLIDLabel);
217 
218  const int nslices = slices.size();
219  for(int isli = 0; isli < nslices; ++isli){
220  // Loop over slices
221  art::Ptr<rb::Cluster> thisSlice = slices.at(isli);
222 
223  if(thisSlice->IsNoise()) continue;
224  if(!prongAssn.isValid()) continue;
225 
226  ncid::NCCosRej nccosrej;
227  std::vector< art::Ptr<rb::Prong> > prongs = prongAssn.at(isli);
228 
229  if(prongs.empty()) continue;
230 
231  nccosrej.SetProngTransMom(TransMomFraction(prongs));
232  std::vector< art::Ptr<cosrej::NueCosRej> > cosrejs = fmcr.at(isli);
233 
234  if(cosrejs.empty()) continue;
235 
236  // Need CVN, input to Real Adaptive BDT
237  art::FindManyP<cvn::Result> fmCVN(slices, evt, fCVNLabel);
238 
239  if(!fmCVN.isValid()) continue;
240 
241  std::vector< art::Ptr<cvn::Result> > cvns = fmCVN.at(isli);
242  // CVN values
243  float cosmicid = 0.;
244  if(!cvns.empty()) cosmicid = cvns[0]->fOutput[14];
245 
246  // LIDShower values
247  float shwGap = -5.;
248  float shwnhit = -5.;
249  float shwnhitx = -5.;
250  float shwnhity = -5.;
251  float shwwwidth = -5.;
252  float shwcalE = -5.;
253  float shwdirY = -5.;
254  float shwlen = -5.;
255  float partptp = -5.;
256  /////////////////////// Fill the PIDs ////////////////////////////
257  // Only execute when there are showers associated with the shwlids
258  if(!showerLidAssn.isValid()) continue;
259 
260  std::vector<const slid::ShowerLID*> shwlids = showerLidAssn.at(isli);
261  if(shwlids.empty()) continue;
262 
263  sort(shwlids.begin(), shwlids.end(), slid::CompareByEnergy);
264  art::FindOneP<rb::Shower> foSh(shwlids, evt, fShowerLIDLabel);
265  float nshwlid = shwlids.size();
266 
267  if(!foSh.isValid()) continue;
268 
269  cet::maybe_ref< art::Ptr<rb::Shower> const > roSh(foSh.at(0));
270  art::Ptr<rb::Shower> shower = roSh.ref();
271  std::vector< art::Ptr<rb::Shower> > showers;
272  showers.push_back(shower);
273  art::FindOneP<slid::ShowerLID> foShLID(showers, evt,fShowerLIDLabel);
274 
275  if(!foShLID.isValid()) continue;
276 
277  cet::maybe_ref<art::Ptr<slid::ShowerLID> const> roShLID(foShLID.at(0));
278  art::Ptr<slid::ShowerLID> showerLID = roShLID.ref();
279 
280  int nMipHit = 0;
281  for(size_t hitIdx = 0; hitIdx < slices[isli]->NCell(); ++hitIdx){
282  const art::Ptr<rb::CellHit>& chptr = slices[isli]->Cell(hitIdx);
283  const rb::RecoHit rhit = slices[isli]->RecoHit(chptr);
284  if(!rhit.IsCalibrated()) continue;
285  if((rhit.PECorr() > 100.0) && (rhit.PECorr() < 245.0)) nMipHit++;
286  } // end of hits loop
287 
288  shwGap = showerLID->Gap();
289  shwnhit = shower->NCell();
290  shwwwidth = showerLID->Radius();
291  shwcalE = showerLID->ShowerEnergy();
292  shwnhitx = shower->NXCell();
293  shwdirY = (float)shower->Dir().Y();
294  shwnhity = shower->NYCell();
295  shwlen = shower->TotalLength();
296  partptp = cosrejs[0]->ParticleShowerTransMom();
297  float shwxminusy = shwnhitx - shwnhity;
298  float shwxplusy = shwnhitx + shwnhity;
299  float shwxovery = shwxminusy/shwxplusy;
300  float nmiphit = nMipHit;
301 
302 
303  TMVAvars[0] = cosmicid; //ncid
304  TMVAvars[1] = partptp; //partptp
305  //TMVAvars[1] = nccosrej.ProngTransMom(); //partptp
306  TMVAvars[2] = shwnhit; //shwnhit;
307  TMVAvars[3] = shwxminusy;
308  TMVAvars[4] = shwxplusy;
309  TMVAvars[5] = shwxovery;
310  TMVAvars[6] = shwcalE; //shwcalE
311  TMVAvars[7] = shwdirY; //shwdirY
312  TMVAvars[8] = shwlen; //shwlen
313  TMVAvars[9] = shwwwidth; //shwwidth
314  TMVAvars[10] = shwGap; //shwgap
315  TMVAvars[11] = nshwlid;
316  TMVAvars[12] = nmiphit; //nmiphit;
317 
318  nccosrej.SetCosPIDDT(fReaderDT->EvaluateMVA("BDTA"));
319  nccosrej.SetCosPIDDTG(fReaderGBDT->EvaluateMVA("BDTG"));
320  ncCosRejjies->push_back(nccosrej);
321  util::CreateAssn(*this, evt, *ncCosRejjies, thisSlice ,*sliceAssn);
322  } // end loop over slices
323 
324  evt.put(std::move(sliceAssn));
325  evt.put(std::move(ncCosRejjies));
326  } // end of producer
327 
328  ///////////////////////////////////////////////////////////////////////
329  double MakeNCCosRej::TransMomFraction(const std::vector< art::Ptr<rb::Prong> >& prongs) const
330  {
331  if(prongs.empty()) return 0;
332 
333  TVector3 ret;
334  for(const art::Ptr<rb::Prong>& prong: prongs){
335  if(!prong->Is3D()) continue;
336  const double w = prong->CalorimetricEnergy();
337  ret += w*prong->Dir();
338  }
339 
340  ret = ret.Unit();
341 
343  const TVector3 beamDir = geom->NuMIBeamDirection();
344  const double beamProj = ret.Dot(beamDir);
345 
346  return (ret-beamProj*beamDir).Mag();
347  }
348 
349  ///////////////////////////////////////////////////////////////////////
351  const art::Ptr<rb::Prong> b)
352  {
353  return a->CalorimetricEnergy() > b->CalorimetricEnergy();
354  }
355 
356 } // end namespace ncid
357 
359 ////////////////////////////////////////////////////////////////////////
float nshwlid
Definition: FillPIDs.h:18
back track the reconstruction to the simulation
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.
float partptp
Definition: NusVarsTemp.cxx:50
MakeNCCosRej(fhicl::ParameterSet const &p)
float shwnhity
float shwGap
void beginRun(art::Run &evt) override
const char * p
Definition: xmltok.h:285
std::string EnvExpansion(const std::string &inString)
Function to expand environment variables.
Definition: EnvExpand.cxx:8
float shwxplusy
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
void reconfigure(const fhicl::ParameterSet &p)
void SetProngTransMom(double input)
Definition: NCCosRej.h:28
DEFINE_ART_MODULE(TestTMapFile)
float shwnhitx
bool CompareByEnergy(const slid::ShowerLID *a, const slid::ShowerLID *b)
Definition: ShowerLID.cxx:51
float shwcalE
Definition: Run.h:31
void produce(art::Event &e) override
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
float cosmicid
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
unsigned short Cell() const
Definition: CellHit.h:40
double TransMomFraction(const std::vector< art::Ptr< rb::Prong > > &prongs) const
Result for CVN.
double CalorimetricEnergy(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple estimate of neutrino energy.
Definition: Cluster.cxx:439
float shwlen
float nmiphit
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
TMVA::Reader * fReaderGBDT
bool CompareByCalEnergy(const art::Ptr< rb::Prong > a, const art::Ptr< rb::Prong > b)
const double a
T get(std::string const &key) const
Definition: ParameterSet.h:231
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
void SetCosPIDDTG(double input)
Definition: NCCosRej.h:34
float shwdirY
Definition: run.py:1
TVector3 NuMIBeamDirection() const
Direction of neutrinos from the NuMI beam (unit vector)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
Calculate dE/dx and log likelihoods for different particle hypotheses. This information will be of us...
float shwxovery
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
float shwxminusy
void geom(int which=0)
Definition: geom.C:163
void SetCosPIDDT(double input)
Definition: NCCosRej.h:31
const hit & b
Definition: hits.cxx:21
bool file_exists(std::string const &qualified_filename)
float shwnhit
float PECorr() const
Definition: RecoHit.cxx:47
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464
Float_t e
Definition: plot.C:35
TMVA::Reader * fReaderDT
fhicl::ParameterSet fParticleIDAlgPSet
Particle ID alorithm (loglikelihoods and dE/dx)
Float_t w
Definition: plot.C:20
Definition: fwd.h:28
Encapsulate the geometry of one entire detector (near, far, ndos)