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  EDProducer(p),
106  fReaderDT(nullptr),
107  fReaderGBDT(nullptr),
108  //fParticleIDAlg(0),
109  fParticleIDAlgPSet(p.get< fhicl::ParameterSet >("ParticleIDAlgPSet"))
110  {
111  produces< std::vector<ncid::NCCosRej> > ();
112  produces< art::Assns<ncid::NCCosRej, rb::Cluster> > ();
113  this->reconfigure(p);
114  }
115 
117  {
118  if(fReaderDT) delete fReaderDT;
119  if(fReaderGBDT) delete fReaderGBDT;
120  }
121 
123  {
125 
126  // Make sure we can find the Real Adaptive Boosted Decision Trees weight file before we set up MVA
127  if(!cet::file_exists(pidlibpath+fNCCosRejPIDFile)){
128  throw cet::exception("MakeNCCosRej")
129  << "Couldn't find file " << pidlibpath+fNCCosRejPIDFile <<std::endl;
130  }
131  // Make sure we can find the Gradient Boosted Decision Trees weight file before we set up MVA
132  if(!cet::file_exists(pidlibpath+fNCCosRejPIDGFile)){
133  throw cet::exception("MakeNCCosRej")
134  << "Couldn't find file " << pidlibpath+fNCCosRejPIDGFile <<std::endl;
135  }
136 
138 
139  fReaderDT = new TMVA::Reader;
140  fReaderDT->AddVariable("cosmicid", &TMVAvars[0]); //cosmicid
141  fReaderDT->AddVariable("partptp", &TMVAvars[1]); //partptp
142  fReaderDT->AddVariable("shwnhit", &TMVAvars[2]); //shwnhit
143  fReaderDT->AddVariable("shwxminusy", &TMVAvars[3]); //shwminusy
144  fReaderDT->AddVariable("shwxplusy", &TMVAvars[4]); //shwxplusy
145  fReaderDT->AddVariable("shwxovery", &TMVAvars[5]); //shwxpvery
146  fReaderDT->AddVariable("shwcalE", &TMVAvars[6]); //shwcalE
147  fReaderDT->AddVariable("shwdirY", &TMVAvars[7]); //shwdirY
148  fReaderDT->AddVariable("shwlen", &TMVAvars[8]); //shwlen
149  fReaderDT->AddVariable("shwwwidth", &TMVAvars[9]); //shwwwidth
150  fReaderDT->AddVariable("shwGap", &TMVAvars[10]); //shwGap
151  fReaderDT->AddVariable("nshwlid", &TMVAvars[11]); //nshwlid
152  fReaderDT->AddVariable("nmiphit", &TMVAvars[12]); //nmiphit
153  fReaderDT->BookMVA("BDTA",pidlibpath+fNCCosRejPIDFile);
154 
155  fReaderGBDT = new TMVA::Reader;
156  fReaderGBDT->AddVariable("cosmicid", &TMVAvars[0]); //cosmicid
157  fReaderGBDT->AddVariable("partptp", &TMVAvars[1]); //partptp
158  fReaderGBDT->AddVariable("shwnhit", &TMVAvars[2]); //shwnhit
159  fReaderGBDT->AddVariable("shwxminusy", &TMVAvars[3]); //shwminusy
160  fReaderGBDT->AddVariable("shwxplusy", &TMVAvars[4]); //shwxplusy
161  fReaderGBDT->AddVariable("shwxovery", &TMVAvars[5]); //shwxpvery
162  fReaderGBDT->AddVariable("shwcalE", &TMVAvars[6]); //shwcalE
163  fReaderGBDT->AddVariable("shwdirY", &TMVAvars[7]); //shwdirY
164  fReaderGBDT->AddVariable("shwlen", &TMVAvars[8]); //shwlen
165  fReaderGBDT->AddVariable("shwwwidth", &TMVAvars[9]); //shwwwidth
166  fReaderGBDT->AddVariable("shwGap", &TMVAvars[10]); //shwGap
167  fReaderGBDT->AddVariable("nshwlid", &TMVAvars[11]); //nshwlid
168  fReaderGBDT->AddVariable("nmiphit", &TMVAvars[12]); //nmiphit
169  fReaderGBDT->BookMVA("BDTG", pidlibpath+fNCCosRejPIDGFile);
170 }
171 
173  {
174  //fMCTruthModuleLabel = p.get< std::string >("MCTruthModuleLabel");
175  fSliceLabel = p.get< std::string >("SliceLabel");
176  fCellHitLabel = p.get< std::string >("CellHitLabel");
177  fProngLabel = p.get< std::string >("ProngLabel");
178  fProngInstance = p.get< std::string >("ProngInstance");
179  fCVNLabel = p.get< std::string >("CVNLabel");
180  fShowerLabel = p.get< std::string >("ShowerLabel");
181  fShowerLIDLabel = p.get< std::string >("ShowerLIDLabel");
182  fPIDLibPath = p.get< std::string >("PIDLibPath");
183  fCosRejLabel = p.get< std::string >("CosRejLabel");
184  fNCCosRejPIDFile = p.get< std::string >("NCCosRejPIDFile");
185  fNCCosRejPIDGFile = p.get< std::string >("NCCosRejPIDGFile");
186  fParticleIDAlgPSet = p.get< fhicl::ParameterSet >("ParticleIDAlgPSet");
187  }
188 
189  //////////////////////////////////////////////////////////////////////////
191  {
194 
195  std::unique_ptr< std::vector<ncid::NCCosRej> >
196  ncCosRejjies(new std::vector< ncid::NCCosRej >);
197  std::unique_ptr< art::Assns<ncid::NCCosRej, rb::Cluster> >
199  // Get the slices in the event
201  evt.getByLabel( fSliceLabel, sliceHandle);
202 
203  std::vector< art::Ptr<rb::Cluster> > slices;
204  art::fill_ptr_vector(slices, sliceHandle);
205 
207  evt.getByLabel(fCellHitLabel, htHandle);
208 
210 
211  for(unsigned int m = 0; m < htHandle->size(); ++m){
212  allHits.push_back(art::Ptr<rb::CellHit>(htHandle,m));
213  }
214  // Get things associated with slices
215  art::FindManyP<rb::Prong> prongAssn(sliceHandle, evt, art::InputTag(fProngLabel, fProngInstance));
216  art::FindMany<slid::ShowerLID> showerLidAssn(sliceHandle, evt, fShowerLIDLabel);
218 
219  const int nslices = slices.size();
220  for(int isli = 0; isli < nslices; ++isli){
221  // Loop over slices
222  art::Ptr<rb::Cluster> thisSlice = slices.at(isli);
223 
224  if(thisSlice->IsNoise()) continue;
225  if(!prongAssn.isValid()) continue;
226 
227  ncid::NCCosRej nccosrej;
228  std::vector< art::Ptr<rb::Prong> > prongs = prongAssn.at(isli);
229 
230  if(prongs.empty()) continue;
231 
232  nccosrej.SetProngTransMom(TransMomFraction(prongs));
233  std::vector< art::Ptr<cosrej::NueCosRej> > cosrejs = fmcr.at(isli);
234 
235  if(cosrejs.empty()) continue;
236 
237  // Need CVN, input to Real Adaptive BDT
238  art::FindManyP<cvn::Result> fmCVN(slices, evt, fCVNLabel);
239 
240  if(!fmCVN.isValid()) continue;
241 
242  std::vector< art::Ptr<cvn::Result> > cvns = fmCVN.at(isli);
243  // CVN values
244  float cosmicid = 0.;
245  if(!cvns.empty()) cosmicid = cvns[0]->fOutput[14];
246 
247  // LIDShower values
248  float shwGap = -5.;
249  float shwnhit = -5.;
250  float shwnhitx = -5.;
251  float shwnhity = -5.;
252  float shwwwidth = -5.;
253  float shwcalE = -5.;
254  float shwdirY = -5.;
255  float shwlen = -5.;
256  float partptp = -5.;
257  /////////////////////// Fill the PIDs ////////////////////////////
258  // Only execute when there are showers associated with the shwlids
259  if(!showerLidAssn.isValid()) continue;
260 
261  std::vector<const slid::ShowerLID*> shwlids = showerLidAssn.at(isli);
262  if(shwlids.empty()) continue;
263 
264  sort(shwlids.begin(), shwlids.end(), slid::CompareByEnergy);
265  art::FindOneP<rb::Shower> foSh(shwlids, evt, fShowerLIDLabel);
266  float nshwlid = shwlids.size();
267 
268  if(!foSh.isValid()) continue;
269 
270  cet::maybe_ref< art::Ptr<rb::Shower> const > roSh(foSh.at(0));
271  art::Ptr<rb::Shower> shower = roSh.ref();
272  std::vector< art::Ptr<rb::Shower> > showers;
273  showers.push_back(shower);
274  art::FindOneP<slid::ShowerLID> foShLID(showers, evt,fShowerLIDLabel);
275 
276  if(!foShLID.isValid()) continue;
277 
278  cet::maybe_ref<art::Ptr<slid::ShowerLID> const> roShLID(foShLID.at(0));
279  art::Ptr<slid::ShowerLID> showerLID = roShLID.ref();
280 
281  int nMipHit = 0;
282  for(size_t hitIdx = 0; hitIdx < slices[isli]->NCell(); ++hitIdx){
283  const art::Ptr<rb::CellHit>& chptr = slices[isli]->Cell(hitIdx);
284  const rb::RecoHit rhit = slices[isli]->RecoHit(chptr);
285  if(!rhit.IsCalibrated()) continue;
286  if((rhit.PECorr() > 100.0) && (rhit.PECorr() < 245.0)) nMipHit++;
287  } // end of hits loop
288 
289  shwGap = showerLID->Gap();
290  shwnhit = shower->NCell();
291  shwwwidth = showerLID->Radius();
292  shwcalE = showerLID->ShowerEnergy();
293  shwnhitx = shower->NXCell();
294  shwdirY = (float)shower->Dir().Y();
295  shwnhity = shower->NYCell();
296  shwlen = shower->TotalLength();
297  partptp = cosrejs[0]->ParticleShowerTransMom();
298  float shwxminusy = shwnhitx - shwnhity;
299  float shwxplusy = shwnhitx + shwnhity;
300  float shwxovery = shwxminusy/shwxplusy;
301  float nmiphit = nMipHit;
302 
303 
304  TMVAvars[0] = cosmicid; //ncid
305  TMVAvars[1] = partptp; //partptp
306  //TMVAvars[1] = nccosrej.ProngTransMom(); //partptp
307  TMVAvars[2] = shwnhit; //shwnhit;
308  TMVAvars[3] = shwxminusy;
309  TMVAvars[4] = shwxplusy;
310  TMVAvars[5] = shwxovery;
311  TMVAvars[6] = shwcalE; //shwcalE
312  TMVAvars[7] = shwdirY; //shwdirY
313  TMVAvars[8] = shwlen; //shwlen
314  TMVAvars[9] = shwwwidth; //shwwidth
315  TMVAvars[10] = shwGap; //shwgap
316  TMVAvars[11] = nshwlid;
317  TMVAvars[12] = nmiphit; //nmiphit;
318 
319  nccosrej.SetCosPIDDT(fReaderDT->EvaluateMVA("BDTA"));
320  nccosrej.SetCosPIDDTG(fReaderGBDT->EvaluateMVA("BDTG"));
321  ncCosRejjies->push_back(nccosrej);
322  util::CreateAssn(evt, *ncCosRejjies, thisSlice ,*sliceAssn);
323  } // end loop over slices
324 
325  evt.put(std::move(sliceAssn));
326  evt.put(std::move(ncCosRejjies));
327  } // end of producer
328 
329  ///////////////////////////////////////////////////////////////////////
330  double MakeNCCosRej::TransMomFraction(const std::vector< art::Ptr<rb::Prong> >& prongs) const
331  {
332  if(prongs.empty()) return 0;
333 
334  TVector3 ret;
335  for(const art::Ptr<rb::Prong>& prong: prongs){
336  if(!prong->Is3D()) continue;
337  const double w = prong->CalorimetricEnergy();
338  ret += w*prong->Dir();
339  }
340 
341  ret = ret.Unit();
342 
344  const TVector3 beamDir = geom->NuMIBeamDirection();
345  const double beamProj = ret.Dot(beamDir);
346 
347  return (ret-beamProj*beamDir).Mag();
348  }
349 
350  ///////////////////////////////////////////////////////////////////////
352  const art::Ptr<rb::Prong> b)
353  {
354  return a->CalorimetricEnergy() > b->CalorimetricEnergy();
355  }
356 
357 } // end namespace ncid
358 
360 ////////////////////////////////////////////////////////////////////////
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.
MakeNCCosRej(fhicl::ParameterSet const &p)
float shwnhity
float shwGap
void beginRun(art::Run &evt) override
const char * p
Definition: xmltok.h:285
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
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:21
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
unsigned short Cell() const
Definition: CellHit.h:40
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
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
float partptp
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:435
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
int evt
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)
Calculate dE/dx and log likelihoods for different particle hypotheses. This information will be of us...
float shwxovery
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:291
Float_t e
Definition: plot.C:35
TMVA::Reader * fReaderDT
ProductID put(std::unique_ptr< PROD > &&edp, FullSemantic< Level::Run > const semantic)
Definition: DataViewImpl.h:730
fhicl::ParameterSet fParticleIDAlgPSet
Particle ID alorithm (loglikelihoods and dE/dx)
Float_t w
Definition: plot.C:20
Definition: fwd.h:29
Encapsulate the geometry of one entire detector (near, far, ndos)
enum BeamMode string