MuonIDProd_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file MuonID_module.cc
3 // \brief BDT trained to identify muons and mubars
4 // \author Connor Johnson - Connor.Johnson@colostate.edu
5 ////////////////////////////////////////////////////////////////////////
6 #include <vector>
7 #include <string>
8 #include <limits>
9 
10 // // ROOT includes
11 
12 #include "TMVA/Reader.h"
13 
14 // // Framework includes
21 #include "Utilities/AssociationUtil.h"
24 #include "fhiclcpp/ParameterSet.h"
25 #include "cetlib/search_path.h"
26 #include "cetlib/filesystem.h"
27 
28 // // NOvA includes
29 #include "MuonID/MuonID.h"
30 #include "ReMId/ReMId.h"
31 #include "TrackInfo/TrackInfoObj.h"
32 #include "RecoBase/Track.h"
35 
36 
37 
38 /// A PID for muons
39 namespace muonid {
40  class MuonIDProd : public art::EDProducer {
41  public:
42 
43  explicit MuonIDProd(fhicl::ParameterSet const &pset);
44 
45  virtual ~MuonIDProd();
46 
47  void produce(art::Event &evt);
48 
49  void beginJob() override;
50 
51  private:
52  void Init();
53 
54  std::string fWeightFile; ///<BDT Weights from training
55  std::string fBDTLabel; ///<Title of BDT in weight training file
56  std::string fTrackLabel; ///<Where to find reconstructed tracks
57  std::string fTrackInfoLabel; ///<TrackInfoObj label
58  std::string fRemidLabel; ///<Where to find the remid
59 
61  TMVA::Reader *fReader;
62  float dedxllh, scatllh, avgdedxlast10cm, avgdedxlast40cm; ///<Objects for BDT calculation to use
63  };
64 
65  //....................................................................................
66 
68  trainfileloaded(false),
69  fReader(NULL),
70  dedxllh(-999),
71  scatllh(-999),
72  avgdedxlast10cm(-999),
73  avgdedxlast40cm(-999)
74  {
75  fWeightFile = pset.get<std::string>("WeightFile");
76  fBDTLabel = pset.get<std::string>("BDTLabel");
77  fTrackLabel = pset.get<std::string>("TrackLabel");
78  fTrackInfoLabel = pset.get<std::string>("TrackInfoLabel");
79  fRemidLabel = pset.get<std::string>("RemidLabel");
80 
81  produces < std::vector < muonid::MuonID > > ();
82  produces < art::Assns < muonid::MuonID, rb::Track > > ();
83  }
84 
85  //....................................................................................
86 
88  if (fReader) { delete fReader; }
89  }
90 
91  //....................................................................................
92 
94  {
95  // Book Method
96  if (!trainfileloaded) {
97  // Locate weight file
98  std::string xmlFileName = util::EnvExpansion("$SRT_PRIVATE_CONTEXT/") + fWeightFile;
99  mf::LogDebug("MuonID") << "Checking in SRT PRIVATE for " << xmlFileName << std::endl;
100  if (!cet::file_exists(xmlFileName)){
101  xmlFileName = util::EnvExpansion("$SRT_PUBLIC_CONTEXT/") + fWeightFile;
102  mf::LogDebug("MuonID") << "Checking in SRT PUBLIC for " << xmlFileName << std::endl;
103  if (!cet::file_exists(xmlFileName))
104  throw cet::exception("MuonIDProd") << "MuonID Train file " << fWeightFile << " not found" << std::endl;
105  }
106 
107  // Set up the TMVA::Reader to read in MVA weight file
108  fReader = new TMVA::Reader();
109 
110  // Add the pid variables
111  fReader->AddVariable("DedxLL", &dedxllh);
112  fReader->AddVariable("ScatLL", &scatllh);
113  fReader->AddVariable("Avededxlast10cm", &avgdedxlast10cm);
114  fReader->AddVariable("Avededxlast40cm", &avgdedxlast40cm);
115 
116  // Book the MVA
117  fReader->BookMVA("BDTG", xmlFileName);
118  trainfileloaded = true;
119  }
120  }
121  //....................................................................................
122 
123 
125  Init();
126 
127  std::unique_ptr<std::vector<muonid::MuonID>> pidcol(new std::vector <muonid::MuonID>);
128  std::unique_ptr<art::Assns<muonid::MuonID, rb::Track>> assncol(new art::Assns <muonid::MuonID, rb::Track>);
129 
130  // Load the tracks
132  evt.getByLabel(fTrackLabel, tracklist);
133 
134  // Load the remid and info
135  art::FindOneP<remid::ReMId> remidTracks(tracklist, evt, fRemidLabel);
136  art::FindOneP<trackinfo::TrackInfoObj> infoTracks(tracklist, evt, fTrackInfoLabel);
137 
138  for(unsigned int itrack = 0; itrack < tracklist->size(); ++itrack)
139  {
140  // MuonID Object to fill
141  MuonID muid(13, -999); // Default to PDG 13 (muon), -999 (value)
142 
143  // Track Ptr for Assns
144  const art::Ptr<rb::Track> track(tracklist, itrack);
145 
146  // Load remid and track info
147  const art::Ptr<remid::ReMId> trackremid = remidTracks.at(itrack);
148  const art::Ptr<trackinfo::TrackInfoObj> trackinfoobj = infoTracks.at(itrack);
149 
150  // Fill and evaluate the BDT
151  dedxllh = trackremid->DedxSeparation();
152  scatllh = trackremid->ScatSeparation();
153  avgdedxlast10cm = std::min(trackinfoobj->AvedEdxTrackLast10cm(), float(30.));
154  avgdedxlast40cm = std::min(trackinfoobj->AvedEdxTrackLast40cm(), float(30.));
155 
156  // Skip calculation if remid failed, or unphysical values
157  if(avgdedxlast10cm < 0 || avgdedxlast40cm < 0)
158  mf::LogDebug("Skipping MuonIDProd") << "Skipping muonid calculation due to invalid remid or unphysical values" << std::endl;
159  else
160  {
161  float muonidval = fReader->EvaluateMVA(fBDTLabel);
162  if(muonidval > -1.0 && muonidval < 1.0){
163  muid.SetVal(muonidval);
164  }
165  else
166  mf::LogWarning("Invalid MuonID") << "Invalid MuonID calculated (outside valid range)" <<
167  "\n\nMuonID = " << muonidval << std::endl;
168  }
169 
170  pidcol->push_back(muid);
171  util::CreateAssn(*this, evt, *(pidcol.get()), track, *(assncol.get()));
172 
173  } // end of loop over slice
174 
175  evt.put(std::move(pidcol));
176  evt.put(std::move(assncol));
177  }
178  //..........................................................................
179 
181  Init();
182  }
183  //..........................................................................
184 
186 
187 }
std::string fTrackInfoLabel
TrackInfoObj label.
std::string fTrackLabel
Where to find reconstructed 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.
void produce(art::Event &evt)
TMVA::Reader * fReader
double ScatSeparation() const
Return the scattering separation variable used as an input to the kNN that determines the pid value...
Definition: ReMId.cxx:206
std::string EnvExpansion(const std::string &inString)
Function to expand environment variables.
Definition: EnvExpand.cxx:8
std::string fWeightFile
BDT Weights from training.
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
DEFINE_ART_MODULE(TestTMapFile)
A PID for muons.
Definition: FillPIDs.h:11
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
T get(std::string const &key) const
Definition: ParameterSet.h:231
int evt
float AvedEdxTrackLast40cm() const
Definition: TrackInfoObj.h:32
MuonIDProd(fhicl::ParameterSet const &pset)
float avgdedxlast40cm
Objects for BDT calculation to use.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
std::string fRemidLabel
Where to find the remid.
void beginJob() override
bool file_exists(std::string const &qualified_filename)
T min(const caf::Proxy< T > &a, T b)
void SetVal(double val)
Definition: PID.h:24
Float_t track
Definition: plot.C:35
double DedxSeparation() const
Return the dE/dx separation variable used as an input to the kNN that determines the pid value...
Definition: ReMId.cxx:199
std::string fBDTLabel
Title of BDT in weight training file.
float AvedEdxTrackLast10cm() const
Definition: TrackInfoObj.h:29