RecVarPID_module.cc
Go to the documentation of this file.
1 // NOvA includes
2 #include "RecVarPID/RVP.h"
4 #include "RecoBase/Track.h"
5 #include "RecoBase/Prong.h"
6 #include "RecoBase/Cluster.h"
7 #include "RecoBase/CellHit.h"
8 #include "RecoBase/FilterList.h"
9 #include "RecoBase/RecoHit.h"
10 #include "RecoBase/Vertex.h"
12 #include "Utilities/AssociationUtil.h"
13 
14 // Framework includes
18 #include "fhiclcpp/ParameterSet.h"
25 
26 #include <vector>
27 #include <string>
28 
29 // ROOT includes
30 #include "TObject.h"
31 #include "TNtuple.h"
32 #include "TMVA/Factory.h"
33 #include "TMVA/Tools.h"
34 #include "TMVA/Reader.h"
35 #include "TString.h"
36 #include "TMath.h"
37 
38 
39 #include <boost/algorithm/string.hpp>
40 #include <boost/foreach.hpp>
41 #include <boost/tokenizer.hpp>
42 
43 /// for std::out
44 #define MAX_BUFFER 100000
45 
46 /// In case this will be changed to whatever in future
47 //#define FUZZYKVERTEXOUTPUT rb::Track
48 //#define FUZZYKVERTEXOUTPUT rb::Shower
49 #define FUZZYKVERTEXOUTPUT rb::Prong
50 
51 
52 namespace rvp {
53  class RecVarPID : public art::EDProducer {
54  public:
55 
56  explicit RecVarPID(fhicl::ParameterSet const& pset);
57  virtual ~RecVarPID();
58  void beginJob();
59  void produce(art::Event& evt);
60 
61  private:
62 
64  double const& w);
65 
66  /// Book rvp reader
67  bool bookVariables(const std::string& weight_file, TMVA::Reader& rvp, std::string& method_name);
68 
69 
70  /// Method to get Variable names from XML
71  std::vector<std::string> getVarNamesFromXML(const std::string weights_file_name) const;
72 
73 
74  /// Get TMVA method name from XML
75  std::string getMethodNameFromXML(const std::string weights_file_name, const bool get_full_name=true) const;
76 
77  protected:
78  std::string fRVPWeightFile; ///< RVP weight file name
79  std::string fRVP12WeightFile; ///< RVP12 weight file name
80  // std::vector<std::string> fExtendedRVPWeightFile; ///< Extended RVP weight file name
81  std::string fSliceLabel; ///< Label for slices
82  std::string fTrackLabel; ///< Label for tracks
83  std::string fProngLabel; ///< Label for prong
84  std::string fInstLabel3D; ///< Instance label for 3d prongs
85  std::string fInstLabel2D; ///< Instance label for 2d prongs
86  std::string fVertexLabel; ///< Label for vertex
87 
88  std::vector<Float_t> fInputVarsRVP; ///< Vector of the variables that are input for TMVA
89 
90  TMVA::Reader* fbdtRVP; ///< Reader for RVP15
91  TMVA::Reader* fbdtRVP12; ///< Reader for RVP15
94 
95  bool fRunTMVA; ///< Run the actual (memory-hungry) TMVA step?
96  };
97 
98  //----------------------------------------------------------------------------
100  fbdtRVP (new TMVA::Reader()),
101  fbdtRVP12 (new TMVA::Reader())
102  {
103  /// Reserve the number of variables to be something big
104  fInputVarsRVP.resize(10000);
105 
106  const std::string bdt = pset.get< std::string >("BDTweight");
107  fRVPWeightFile = pset.get< std::string >("RVPWeightFile");
108  fRVP12WeightFile = pset.get< std::string >("RVP12WeightFile");
109  fSliceLabel = pset.get< std::string >("SliceLabel");
110  fTrackLabel = pset.get< std::string >("TrackLabel");
111  fProngLabel = pset.get< std::string >("ProngLabel");
112  fInstLabel3D = pset.get< std::string >("InstLabel3D");
113  fInstLabel2D = pset.get< std::string >("InstLabel2D");
114  fVertexLabel = pset.get< std::string >("VertexLabel");
115 
116  fRunTMVA = pset.get< bool >("RunTMVA");
117 
118  /// Expand fcl input path to BDTs
119  std::string bdtpath = util::EnvExpansion(bdt);
120 
121  /// Add BDT path to the files
122  fRVPWeightFile = bdtpath + fRVPWeightFile;
123  fRVP12WeightFile = bdtpath + fRVP12WeightFile;
124 
125 
126  produces< std::vector<rvp::RVP> >();
127  produces< art::Assns<rvp::RVP,rb::Cluster> >();
128 
129  }
130 
131  //...........................................................................
133  {
134  if(fbdtRVP) delete fbdtRVP;
135  if(fbdtRVP12) delete fbdtRVP12;
136 
137  }
138 
139  //...........................................................................
141  {
142  if(fRunTMVA){
143  /// Need to set Weight file
144  fbdtRVPMethodName = "RVP ";
146 
147  /// Need to set Weight file
148  fbdtRVP12MethodName = "RVP12 ";
150  }
151  }
152 
153  //......................................................................
155  double const& w)
156  {
158  rb::RecoHit rhit(cal->MakeRecoHit(*chit, w));
159 
160  return rhit;
161  }
162 
163  //...........................................................................
165  {
167 
168  typedef std::unique_ptr<std::vector<rvp::RVP> > RVPs_t;
169  typedef std::unique_ptr<art::Assns<rvp::RVP,rb::Cluster> > RVPAssn_t;
170  RVPs_t pidcol (new std::vector<rvp::RVP>);
171  RVPAssn_t assncol(new art::Assns<rvp::RVP,rb::Cluster>);
172 
173  double BDT12 = -1.0;//set the default 12 variable BDT output to -1
174  double BDT15 = -1.0;//set the default 15 variable BDT output to -1
175 
177  evt.getByLabel(fSliceLabel, slices);
178 
179 
180  const art::FindManyP<rb::Track> fmTrack (slices, evt, fTrackLabel);
181  const art::FindManyP<rb::Vertex> fmVertex (slices, evt, fVertexLabel);
182 
183  //loop over all slices
184  for(unsigned int sliceIdx = 0; sliceIdx < slices->size(); ++sliceIdx){
185 
186  if(rb::IsFiltered(evt, slices, sliceIdx)) continue;
187  const rb::Cluster& slice = (*slices)[sliceIdx];
188 
189  if(slice.IsNoise()) continue;
190 
191  const unsigned int nCells = slice.NCell();//# of cells
192  const int lengthPlanes = slice.ExtentPlane();//# of extended planes
193 
194  //select the longest track
195  double longestTrack = 0.0;//length of longest track
196  int longestTrack_ncells = 0; //# of cells associated to the longest track
197 
198  if (fmTrack.isValid()){
199  std::vector<art::Ptr<rb::Track>> tracks = fmTrack.at(sliceIdx);
200 
201  for(unsigned int n = 0; n < tracks.size(); ++n){
202  const art::Ptr<rb::Track> current_track = tracks[n];
203  const double L = current_track->TotalLength();
204  if(L > longestTrack){
205  longestTrack = L;
206  longestTrack_ncells = current_track->NCell();
207  }
208  }
209  }
210 
211  //mip cells and reco energy for slicer
212  double recoEnergy = 0.0;//reco energy of slicer
213  double nmip = 0.0;//# of mip cells
214  double Ntotal = 0.0;//# of total cells
215 
216  const unsigned int ncell_hits = slice.NCell();
217  double wx = (slice.NXCell() > 0) ? slice.W(slice.XCell(0).get()) : 0;
218  double wy = (slice.NYCell() > 0) ? slice.W(slice.YCell(0).get()) : 0;
219  std::vector<rb::RecoHit> all_reco_hits;
220 
221  for( unsigned int icell=0; icell<ncell_hits; ++icell){
222  const art::Ptr<rb::CellHit>& chit = slice.Cell(icell);
223  all_reco_hits.push_back(this->MakeRecoHit(chit, chit->View() == geo::kX ? wx : wy));
224  const rb::RecoHit& rhit = all_reco_hits[icell];
225 
226  if( !rhit.IsCalibrated() ) continue;
227  recoEnergy += rhit.GeV();
228 
229  if( rhit.PECorr()>100. && rhit.PECorr()<245. ) ++nmip;
230  ++Ntotal;
231  }
232 
233  //reco vertex
234  //only continue if this get by label succeeded
235  if (fmVertex.isValid()){
236  std::vector<art::Ptr<rb::Vertex>> vert = fmVertex.at(sliceIdx);
237 
238  //pass over slice if there was no vertex
239  //NOTE: right now it is impossible for a slice to have more then one vertex, so the size
240  //should be exactly one
241  if (vert.size() != 1) continue;
242  //this is extra preselection that is being removed so nothing is hidden from caf
243  //if( recoEnergy>0.05 && nCells>5 && lengthPlanes<140 && vert.size()>0 ){
244  //if( recoEnergy>0.05 && nCells>5 && lengthPlanes<140 ){
245  //calculate per plane cell energy
246  std::vector<double> PlaneSumEnergy;
247  double totalE_road2sigma = 0.0; //total energy in 2 sigma road
248  double Esig_3 = 0.0; //energy not in a 3 sigma window of shower core
249 
250  for( int iplane=0; iplane<lengthPlanes; ++iplane ){
251  double cellX = 0.0;//energy-weighted cell X in each plane
252  double cellY = 0.0;//energy-weighted cell Y in each plane
253  double Ecell = 0.0;//sum cell energy in each plane
254  double EXcell = 0.0;
255  double EYcell = 0.0;
256 
257  std::vector<double> Xxview;
258  std::vector<double> Exview;
259  std::vector<double> Yyview;
260  std::vector<double> Eyview;
261  Xxview.clear();
262  Exview.clear();
263  Yyview.clear();
264  Eyview.clear();
265 
266  for( unsigned int icell=0; icell<ncell_hits; ++icell){
267  const art::Ptr<rb::CellHit>& chit = slice.Cell(icell);
268  const rb::RecoHit& rhit = all_reco_hits[icell];
269  if( !rhit.IsCalibrated() ) continue;
270 
271  const int Dplane=chit->Plane()-slice.MinPlane();
272  if( Dplane<0 ) std::cout<<"plane # of this cell is beyond the minimal plane."<<std::endl;
273 
274  if( Dplane==iplane ){
275  Ecell += rhit.GeV();
276  EXcell += rhit.GeV()*rhit.X();
277  EYcell += rhit.GeV()*rhit.Y();
278 
279  if( chit->View() == geo::kX ){
280  Xxview.push_back(rhit.X());
281  Exview.push_back(rhit.GeV());
282  }//in X view
283  else {
284  Yyview.push_back(rhit.Y());
285  Eyview.push_back(rhit.GeV());
286  }//in Y view
287 
288  }
289  }//loop all cells to calculate the centroid cell position in each plane
290 
291  if( Ecell>0. ){
292  cellX = EXcell/Ecell;
293  cellY = EYcell/Ecell;
294  }
295 
296  PlaneSumEnergy.push_back(Ecell);
297 
298  const int Nxviews=Xxview.size();
299  const int Nyviews=Yyview.size();
300  if( Nxviews>0 ){
301  for( int ixv=0; ixv<Nxviews; ++ixv ){
302  if( fabs(Xxview[ixv]-cellX) < 2.*2. ) totalE_road2sigma+=Exview[ixv];
303  if( fabs(Xxview[ixv]-cellX) < 2.*3. ) Esig_3+=Exview[ixv];
304  }
305  }//in X view
306  if( Nyviews>0 ){
307  for( int iyv=0; iyv<Nyviews; ++iyv ){
308  if( fabs(Yyview[iyv]-cellY) < 2.*2. ) totalE_road2sigma+=Eyview[iyv];
309  if( fabs(Yyview[iyv]-cellY) < 2.*3. ) Esig_3+=Eyview[iyv];
310  }
311  }//in Y view
312 
313  }//loop all planes
314 
315  double efrac_2slides = -1.0;//energy fraction in 2 sliding window
316  double efrac_6slides = -1.0;//energy fraction in 6 sliding window
317  double Eplane20 = 0.0;
318  int nplanes=PlaneSumEnergy.size();
319  if (recoEnergy > 0.0){
320  for( int i=0; i<nplanes; ++i ){
321  if( i<20 ) Eplane20 += PlaneSumEnergy[i];
322 
323  //maximal energy fraction in 2 sliding window
324  if( i<nplanes-2 ){
325  const double e_frac2=(PlaneSumEnergy[i]+PlaneSumEnergy[i+1])/recoEnergy;
326  if( efrac_2slides<e_frac2 ) efrac_2slides=e_frac2; }
327  else if( nplanes<=2 ) efrac_2slides=1.;
328  //maximal energy fraction in 6 sliding window
329  if( i<nplanes-6 ){
330  const double e_frac6=(PlaneSumEnergy[i]+PlaneSumEnergy[i+1]+PlaneSumEnergy[i+2]+PlaneSumEnergy[i+3]+PlaneSumEnergy[i+4]+PlaneSumEnergy[i+5])/recoEnergy;
331  if( efrac_6slides<e_frac6 ) efrac_6slides=e_frac6; }
332  else if( nplanes<=6 ) efrac_6slides=1.;
333  }//loop all planes
334  }//end if
335  //reco prongs
336  //Pulling Tracks from event since this is what FuzzyK makes for others
337  //Tracks are then converted to prongs for the rest of this code
338  //art::Handle< std::vector<rb::Track> > pngcol;
339  //evt.getByLabel(fProngLabel, fInstLabel, pngcol);
340  //std::vector<const rb::Prong*> SelectedProng;
343  std::vector<art::Ptr<FUZZYKVERTEXOUTPUT>> prongs = fmProng3D.at(0);
344  std::vector<art::Ptr<FUZZYKVERTEXOUTPUT>> SelectedProng;
345  std::vector<art::Ptr<FUZZYKVERTEXOUTPUT>> SelectedProng2D;
346  //check if instance label is defined to handle sorting of prongs
347  if ((fInstLabel3D != "")&&(fInstLabel2D != "")){
348  SelectedProng = fmProng3D.at(0);
349  SelectedProng2D = fmProng2D.at(0);
350  }
351  else{
352  for( unsigned int ip=0; ip<prongs.size(); ++ip ){
353  art::Ptr<FUZZYKVERTEXOUTPUT> current_prong = prongs[ip];
354  if(current_prong->IsNoise() ) continue;
355  if(current_prong->View() == geo::kXorY)//3D prongs
356  SelectedProng.push_back(current_prong);
357  else
358  SelectedProng2D.push_back(current_prong);
359  }
360  }
361 
362  int nprongs = SelectedProng.size();
363  double Ebalance = 1.0;//energy balance between 2 most energetic prongs
364  std::vector<double> ProngEnergy;
365  std::vector<double> ProngEnergyXView;
366  std::vector<double> ProngEnergyYView;
367  for(int ip=0; ip<nprongs; ++ip ){
368  double Eprong = 0.0;
369  double EprongX = 0.0;
370  double EprongY = 0.0;
371  const art::Ptr<FUZZYKVERTEXOUTPUT> current_prong = SelectedProng[ip];
372 
373  const unsigned int nselectedprong_cells = current_prong->NCell();
374 
375  for( unsigned int icell=0; icell<nselectedprong_cells; ++icell){
376  const art::Ptr<rb::CellHit>& chit = current_prong->Cell(icell);
377  const rb::RecoHit rhit = current_prong->RecoHit(icell);
378 
379  if( !rhit.IsCalibrated() ) continue;
380  Eprong += rhit.GeV();
381  if( chit->View() == geo::kX ) EprongX += rhit.GeV();
382  if( chit->View() == geo::kY ) EprongY += rhit.GeV();
383  }
384  ProngEnergy.push_back(Eprong);
385  ProngEnergyXView.push_back(EprongX);
386  ProngEnergyYView.push_back(EprongY);
387  }//loop all prongs
388 
389  if (nprongs > 1){
390  int Index1 = -1;
391  int Index2 = -1;
392  double MaxEnergy1 = -999.0;
393  double MaxEnergy2 = -999.0;
394  for( int i=0;i<nprongs; ++i ){
395  if( MaxEnergy1<ProngEnergy[i] ){
396  MaxEnergy1 = ProngEnergy[i];
397  Index1 = i;
398  }
399  }//leading
400  for( int i=0;i<nprongs; ++i ){
401  if( i==Index1 ) continue;
402  if( MaxEnergy2<ProngEnergy[i] ){
403  MaxEnergy2=ProngEnergy[i];
404  Index2=i;
405  }
406  }//sub-leading
407  double Eden = ProngEnergy[Index1]+ProngEnergy[Index2];
408  double Enum = ProngEnergy[Index1]-ProngEnergy[Index2];
409  if( Eden>0.0 ) Ebalance=Enum/Eden;
410  }//has at least 2 prongs
411  //study 2D prongs
412 
413  std::vector<int> Prong2DInXView;
414  std::vector<double> ProngEnergy2D;
415  int nprongs2D=SelectedProng2D.size();
416  for(int ip=0; ip<nprongs2D; ++ip ){
417  const art::Ptr<FUZZYKVERTEXOUTPUT> current_prong = SelectedProng2D[ip];
418 
419  int prongview=SelectedProng2D[ip]->View();
420  Prong2DInXView.push_back(prongview);
421 
422  //reco momentum
423  double Eprong = 0.0;
424 
425  const unsigned int nselectedprong2d_cells = current_prong->NCell();
426  for( unsigned int icell=0; icell<nselectedprong2d_cells; ++icell){
427  const rb::RecoHit rhit = current_prong->RecoHit(icell);
428 
429  if( !rhit.IsCalibrated() ) continue;
430  Eprong += rhit.GeV();
431  }
432  ProngEnergy2D.push_back(Eprong);
433  }
434 
435  int index3D = -1;
436  double max3D = -99.0;
437  for( int i=0; i<nprongs; ++i ){
438  if( max3D<ProngEnergy[i] ){
439  max3D = ProngEnergy[i];
440  index3D = i;
441  }
442  }//loop 3D prongs
443  int index2D = -1;
444  double max2D = -99.;
445  for( int i=0;i<nprongs2D; ++i ){
446  if( max2D<ProngEnergy2D[i] ){
447  max2D = ProngEnergy2D[i];
448  index2D = i;
449  }
450  }//loop 2D prongs
451  double ebalance2D = 1.0;
452  if( nprongs>0 && nprongs2D>0 ){
453  double Eden = 0.0;
454  double Enum = 0.0;
455  if( Prong2DInXView[index2D]==0 ){//XZ view
456  Eden = ProngEnergy2D[index2D]+ProngEnergyXView[index3D];
457  Enum = ProngEnergy2D[index2D]-ProngEnergyXView[index3D];
458  }
459  else {//YZ view
460  Eden = ProngEnergy2D[index2D]+ProngEnergyYView[index3D];
461  Enum = ProngEnergy2D[index2D]-ProngEnergyYView[index3D];
462  }
463  if( Eden > 0.0 ) ebalance2D=fabs(Enum)/Eden;
464  }
465 
466  fInputVarsRVP[0] = (Float_t)nCells;
467  fInputVarsRVP[1] = (Float_t)recoEnergy;
468  fInputVarsRVP[2] = (Float_t)longestTrack;
469  if (nCells > 0) fInputVarsRVP[3] = (Float_t)longestTrack_ncells/nCells;
470  else fInputVarsRVP[3] = -1.0;
471  if (Ntotal > 0) fInputVarsRVP[4] = (Float_t)nmip/Ntotal;
472  else fInputVarsRVP[4] = -1.0;
473  fInputVarsRVP[5] = (Float_t)nmip;
474  //make sure total rec energy is nonzero
475  if (recoEnergy > 0.0) fInputVarsRVP[6] = (Float_t)Eplane20/recoEnergy;
476  else fInputVarsRVP[6] = -1.0;
477  fInputVarsRVP[7] = (Float_t)efrac_2slides;
478  fInputVarsRVP[8] = (Float_t)efrac_6slides;
479  if (recoEnergy > 0.0) fInputVarsRVP[9] = (Float_t)totalE_road2sigma/recoEnergy;
480  else fInputVarsRVP[9] = -1.0;
481  fInputVarsRVP[10] = (Float_t)nprongs;
482  fInputVarsRVP[11] = (Float_t)Ebalance;
483  fInputVarsRVP[12] = (Float_t)nprongs2D;
484  fInputVarsRVP[13] = (Float_t)ebalance2D;
485  if (Esig_3 > 0.0) fInputVarsRVP[14] = (Float_t)(recoEnergy - Esig_3)/Esig_3;
486  else fInputVarsRVP[14] = -1.0;
487 
488  if(fRunTMVA){
489  BDT15 = fbdtRVP->EvaluateMVA(fbdtRVPMethodName);
490  BDT12 = fbdtRVP12->EvaluateMVA(fbdtRVP12MethodName);
491  }
492 
493  /// Slice for the PID
494  const art::Ptr<rb::Cluster> bestSlice(slices,sliceIdx);
495 
496  RVP pid;
497  pid.SetVal (BDT15);
498  pid.SetPdg (12);
499  pid.SetNCells (fInputVarsRVP[0]);
500  pid.SetRecoE (fInputVarsRVP[1]);
501  pid.SetLongTr (fInputVarsRVP[2]);
502  pid.SetLongTrFrac (fInputVarsRVP[3]);
503  pid.SetMipFrac (fInputVarsRVP[4]);
504  pid.SetMipHits (fInputVarsRVP[5]);
505  pid.SetEPl20Frac (fInputVarsRVP[6]);
506  pid.SetEFrac2PlWin (fInputVarsRVP[7]);
507  pid.SetEFrac6PlWin (fInputVarsRVP[8]);
508  pid.SetEFrac2SigRd (fInputVarsRVP[9]);
509  pid.SetProngs3D (fInputVarsRVP[10]);
510  pid.SetProngEBal3D (fInputVarsRVP[11]);
511  pid.SetProngs2D (fInputVarsRVP[12]);
512  pid.SetProngEBal2D (fInputVarsRVP[13]);
513  pid.SetEIso3Sig (fInputVarsRVP[14]);
514  pid.SetRVP12 (BDT12);
515 
516  pidcol->push_back(pid);
517  util::CreateAssn(*this,evt,*(pidcol.get()),bestSlice,*(assncol.get()));
518 
519 
520 
521  }//check Assn between Vertex and Slice exists
522  }//has reco slicer
523 
524  /// If producing RVP, move products
525  evt.put(std::move(pidcol));
526  evt.put(std::move(assncol));
527 
528  }// end of produce
529 
530 
531  //...........................................................................
532  bool RecVarPID::bookVariables(const std::string& weight_file, TMVA::Reader& rvp, std::string& method_name)
533  {
534  /// Get all variable names from the weight XML
535  const std::vector<std::string> varnames_from_xml = getVarNamesFromXML(weight_file);
536 
537  /// Number of variables
538  const unsigned int nvars = varnames_from_xml.size();
539 
540  /// Add variables to TMVA
541  for(unsigned int i=0; i<nvars; ++i){
542  rvp.AddVariable(varnames_from_xml[i], &fInputVarsRVP[i]);
543  }// end of adding variables to TMVA
544 
545  /// Extract method name and add it to the current name
546  method_name += getMethodNameFromXML(weight_file);
547 
548  /// Book TMVA
549  rvp.BookMVA(method_name, weight_file);
550 
551  return true;
552  }// end of RecVarPID::bookVariables(TMVA::Reader& rvp)
553 
554 
555  //---------------------------------------------------------------------------
556  // Parse XML file to extract TMVA method name
558  const bool get_full_name
559  ) const
560  {
561 
562  /// System command to read out variables
563  std::string cmd = "cat " + weights_file_name;
564  cmd += "| grep \"<MethodSetup\"";
565 
566  /// d
567  std::string data; ///< stream from system
568 
569  /// Perform a system call of cmd with the output stream going to data
570  {
571  FILE *stream;
572  char buffer[MAX_BUFFER];
573 
574  stream = popen(cmd.c_str(), "r");
575 
576  while ( fgets(buffer, MAX_BUFFER, stream) != NULL ){
577  data.append(buffer);
578  }
579  pclose(stream);
580  }// end of the system call
581 
582  /// Replace stuff
583  boost::replace_all(data, "<MethodSetup Method=", "");
584  boost::replace_all(data, ">" , "");
585  boost::replace_all(data, "\"" , "");
586  boost::replace_all(data, " " , "");
587  boost::replace_all(data, "\n" , "");
588 
589  /// If getting a full name, return now
590  /// Full name is, for instance, "BDT::BDTG"
591  /// Short name would be "BDTG"
592  if(get_full_name) return data;
593 
594 
595  /// Getting the short name from the full name
596 
597  std::vector<std::string> names;
598 
599  /// Tokenize the string, with a separator ";"
600  const boost::char_separator<char> sep(":");
601  boost::tokenizer< boost::char_separator<char> > tokens(data, sep);
602  BOOST_FOREACH (const std::string& t, tokens) {
603 
604  /// If the size is non-zero, add the string to the vector of arguments
605  if(t.size() > 0)
606  names.push_back(t);
607  }// end of loop over tokens
608 
609  std::string out_string;
610 
611  if(names.size() > 0)
612  out_string = names[names.size() - 1];
613 
614  return out_string;
615  }
616 
617  //---------------------------------------------------------------------------
618  // Parse XML file to extract variable names
619  std::vector<std::string> RecVarPID::getVarNamesFromXML(const std::string file_name) const
620  {
621  /// System command to read out variables
622  std::string cmd = "cat " + file_name;
623  cmd += "| grep \"Variable VarIndex=\" | awk \'{print $3}\'";
624 
625  /// d
626  std::string data; ///< stream from system
627 
628  /// Perform a system call of cmd with the output stream going to data
629  {
630  FILE *stream;
631  char buffer[MAX_BUFFER];
632 
633  stream = popen(cmd.c_str(), "r");
634 
635  while ( fgets(buffer, MAX_BUFFER, stream) != NULL ){
636  data.append(buffer);
637  }
638  pclose(stream);
639  }// end of the system call
640 
641  /// Replace stuff
642  boost::replace_all(data, "Expression=", ";");
643  boost::replace_all(data, "\"" , "");
644  boost::replace_all(data, "\n" , "");
645  boost::replace_all(data, "," , "");
646  boost::replace_all(data, " " , "");
647  /// After the replacement the string of variables have variables separated by ";"
648 
649  /// Output vector of variables in the XML file
650  std::vector<std::string> vars;
651 
652  /// Tokenize the string, with a separator ";"
653  const boost::char_separator<char> sep(";");
654  boost::tokenizer< boost::char_separator<char> > tokens(data, sep);
655  BOOST_FOREACH (const std::string& t, tokens) {
656 
657  /// If the size is non-zero, add the string to the vector of arguments
658  if(t.size() > 0)
659  vars.push_back(t);
660  }// end of loop over tokens
661 
662  return vars;
663  }
664 
666 
667 }//end namespace
668 
std::string fTrackLabel
Label for tracks.
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.
virtual double W(const rb::CellHit *chit) const
Estimate the unmeasured coordinate of chit.
Definition: Cluster.cxx:121
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
void SetRecoE(float recE)
Definition: RVP.h:43
std::string getMethodNameFromXML(const std::string weights_file_name, const bool get_full_name=true) const
Get TMVA method name from XML.
art::Ptr< rb::CellHit > XCell(unsigned int xIdx) const
Get the ith cell in the x-view.
Definition: Cluster.cxx:157
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
RecVarPID(fhicl::ParameterSet const &pset)
X or Y views.
Definition: PlaneGeo.h:30
unsigned short Plane() const
Definition: CellHit.h:39
geo::View_t View() const
Definition: CellHit.h:41
void SetEFrac2SigRd(float eFrac2Sig)
Definition: RVP.h:51
std::string EnvExpansion(const std::string &inString)
Function to expand environment variables.
Definition: EnvExpand.cxx:8
Vertical planes which measure X.
Definition: PlaneGeo.h:28
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
A collection of associated CellHits.
Definition: Cluster.h:47
TString ip
Definition: loadincs.C:5
std::string fInstLabel3D
Instance label for 3d prongs.
DEFINE_ART_MODULE(TestTMapFile)
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
void SetProngEBal3D(float pngEBal)
Definition: RVP.h:53
std::string fbdtRVP12MethodName
::xsd::cxx::tree::buffer< char > buffer
Definition: Database.h:179
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
std::string fInstLabel2D
Instance label for 2d prongs.
const XML_Char const XML_Char * data
Definition: expat.h:268
void SetNCells(unsigned int ncell)
Definition: RVP.h:42
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
virtual double TotalLength() const
Length (cm) of all the track segments.
Definition: Track.cxx:213
bool replace_all(std::string &in, std::string const &from, std::string const &to)
Replace all occurrences of from in string with to.
TMVA::Reader * fbdtRVP12
Reader for RVP15.
string cmd
Definition: run_hadd.py:52
void SetProngEBal2D(float pngEBal2D)
Definition: RVP.h:55
void SetProngs3D(int png3D)
Definition: RVP.h:52
static constexpr double L
void SetEFrac6PlWin(float eFrac6Plane)
Definition: RVP.h:50
bool fRunTMVA
Run the actual (memory-hungry) TMVA step?
art::Ptr< rb::CellHit > YCell(unsigned int yIdx) const
Get the ith cell in the y-view.
Definition: Cluster.cxx:165
void SetProngs2D(int png2D)
Definition: RVP.h:54
T get(std::string const &key) const
Definition: ParameterSet.h:231
Definition: RVP.h:16
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
bool IsFiltered(const art::Event &evt, art::Ptr< T > x, const std::vector< std::string > &labels)
Is this Ptr marked "filtered out"?
Definition: FilterList.h:96
const std::map< std::pair< std::string, std::string >, Variable > vars
std::string fVertexLabel
Label for vertex.
void SetEIso3Sig(float eIso3Sig)
Definition: RVP.h:56
Vertex location in position and time.
void SetLongTr(float lTrack)
Definition: RVP.h:44
std::string fbdtRVPMethodName
TMVA::Reader * fbdtRVP
Reader for RVP15.
void produce(art::Event &evt)
const char sep
void SetMipHits(float mHits)
Definition: RVP.h:47
OStream cout
Definition: OStream.cxx:6
int nplanes
Definition: geom.C:145
void SetMipFrac(float mFrac)
Definition: RVP.h:46
void SetPdg(int pdg)
Definition: PID.h:23
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::string fRVP12WeightFile
RVP12 weight file name.
std::vector< std::string > getVarNamesFromXML(const std::string weights_file_name) const
Method to get Variable names from XML.
float GeV() const
Definition: RecoHit.cxx:69
bool bookVariables(const std::string &weight_file, TMVA::Reader &rvp, std::string &method_name)
Book rvp reader.
#define MAX_BUFFER
for std::out
const int nvars
T const * get() const
Definition: Ptr.h:321
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
void SetEFrac2PlWin(float eFrac2Plane)
Definition: RVP.h:49
float X() const
Definition: RecoHit.h:36
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
unsigned int ExtentPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.h:250
std::string fSliceLabel
Label for slices.
std::string fProngLabel
Label for prong.
Definition: tmvaglob.h:28
void SetEPl20Frac(float ePlane20)
Definition: RVP.h:48
float Y() const
Definition: RecoHit.h:37
void SetLongTrFrac(float ltFrac)
Definition: RVP.h:45
float PECorr() const
Definition: RecoHit.cxx:47
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
void SetVal(double val)
Definition: PID.h:24
Float_t w
Definition: plot.C:20
std::string fRVPWeightFile
RVP weight file name.
Definition: FillPIDs.h:15
rb::RecoHit MakeRecoHit(art::Ptr< rb::CellHit > const &chit, double const &w)
void SetRVP12(float rvp12)
Definition: RVP.h:57
std::vector< Float_t > fInputVarsRVP
Vector of the variables that are input for TMVA.