CosRej_module.cc
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////
2 /// \brief Cosmic Rejection PIDs for Numu analysis
3 /// \author Kirk Bays bays@caltech.edu
4 ////////////////////////////////////////////////////////
5 #include <string>
6 
11 #include "art_root_io/TFileService.h"
12 #include "fhiclcpp/ParameterSet.h"
14 
15 #include "CosRej/CosRejFxs.h"
16 #include "CosRej/CosRejObj.h"
17 #include "CosRej/TrkCntObj.h"
18 #include "Geometry/Geometry.h"
19 #include "Geometry/LiveGeometry.h"
20 #include "MCCheater/BackTracker.h"
21 #include "NovaDAQConventions/DAQConventions.h"
22 #include "RecoBase/CellHit.h"
23 #include "RecoBase/Cluster.h"
24 #include "RecoBase/Track.h"
25 #include "ReMId/ReMId.h"
26 #include "SummaryData/SpillData.h"
27 #include "TimingFit/TimingFitAlg.h"
28 #include "Utilities/AssociationUtil.h"
30 #include "CVN/func/Result.h"
31 
32 #include "TFile.h"
33 #include "TMVA/Reader.h"
34 
35 #include "TH1.h"
36 #include "TH2.h"
37 namespace cosrej {
38 
39  class CosRej : public art::EDProducer {
40  public:
41  explicit CosRej(fhicl::ParameterSet const & pset);
42  virtual ~CosRej();
43 
44  void produce (art::Event & evt);
45  void beginJob();
46 
47  protected:
48  void Init();
49 
61 
62  TMVA::Reader* fReaderUCTuned;
63  TMVA::Reader* fReaderTAP1;
64  TMVA::Reader* fReaderTAP2;
65  TMVA::Reader* fReaderTAP3;
66 
67  TMVA::Reader* fReaderProd5_Per1;
68  TMVA::Reader* fReaderProd5_Per2;
69  TMVA::Reader* fReaderProd5_FHC;
70  TMVA::Reader* fReaderProd5_RHC;
71 
72  float TMVAvarsTA [7];
73  float TMVAvarsUCTuned [74] = {0.};
74  float TMVAvarsProd5 [12] = {0.};
75 
79  /*
80  art::ServiceHandle<art::TFileService> tfs;
81  TH1D* fUncorrTA;
82  TH1D* fCorrTA;
83  TH1D* fDiffTA;
84  TH1D* fCosRej20;
85  TH2D* fCosRej20_CorrTA;
86  */
87  };
88 }
89 
90 
91 namespace cosrej
92 {
93 
94  //----------------------------------------------------------------------
96  EDProducer(pset),
98  fSlicerToken(consumes<std::vector<rb::Cluster>>(pset.get<std::string>("SlicerLabel"))),
99  fGeneratorLabel (pset.get< std::string >("GeneratorLabel")),
100  fNuMIBeamLabel (pset.get< std::string >("NuMIBeamLabel")),
101  fTrackLabel (pset.get< std::string >("TrackLabel")),
102  fCosmicTrackLabel (pset.get< std::string >("CosmicTrackLabel")),
103  fPIDLabel (pset.get< std::string >("PIDLabel")),
104  fCVNLabel (pset.get< std::string >("CVNLabel")),
105  fTMVApath (pset.get< std::string >("tmvapath")),
106  fUseLongestTrack (pset.get< bool >("UseLongestTrack")),
107  fTrackContainmentOnly (pset.get< bool >("TrackContainmentOnly")),
108 
109  fReaderUCTuned (nullptr),
110  fReaderTAP1 (nullptr),
111  fReaderTAP2 (nullptr),
112  fReaderTAP3 (nullptr),
113 
114  fReaderProd5_Per1 (nullptr),
115  fReaderProd5_Per2 (nullptr),
116  fReaderProd5_FHC (nullptr),
117  fReaderProd5_RHC (nullptr)
118  {
119  produces< std::vector<cosrej::CosRejObj> >();
120  produces< art::Assns<cosrej::CosRejObj, rb::Cluster> >();
121  produces< std::vector<cosrej::TrkCntObj> >();
122  produces< art::Assns<cosrej::TrkCntObj, rb::Track> >();
123  }
124 
125 
126  //----------------------------------------------------------------------
128  {
129  if(fCosRejFxs) delete fCosRejFxs;
130  if(fReaderUCTuned) delete fReaderUCTuned;
131  if(fReaderTAP1) delete fReaderTAP1;
132  if(fReaderTAP2) delete fReaderTAP2;
133  if(fReaderTAP3) delete fReaderTAP3;
134 
139  }
140 
141  //-------------------------------------------------------------------
143  {
144  Init();
145  /*
146  fUncorrTA = tfs->make<TH1D>("UncorTA" , "Uncorrected TA; BDT Score, Number of slices" , 100, 0. , 1. );
147  fCorrTA = tfs->make<TH1D>("CorrTA" , "Corrected TA; BDT Score; Number of slices" , 100, 0. , 1. );
148  fDiffTA = tfs->make<TH1D>("DiffTA" , "Diff in Scores; Corr - Uncorr Score; Number of slices", 100, -0.5, 0.5 );
149  fCosRej20 = tfs->make<TH1D>("CosRej20", "CosRej 20; BDT Score; Number of slices" , 100, 0. , 1. );
150 
151  fCosRej20_CorrTA = tfs->make<TH2D>("CosRej20_CorrTA", "CosRej20 vs TA BDT Scores; Cos Rej 20 BDT Score; Corrected TA BDT Score",
152  100, 0., 1, 100, 0., 1. );
153  */
154  }
155 
156  //-------------------------------------------------------------------
158  {
159 
161 
162  // Easier to just declare regardless of detector. Memory footprint is small for BDTs
163  //if(geom->DetId() == novadaq::cnv::kFARDET) {
164 
165  fReaderUCTuned = new TMVA::Reader;
166  fReaderUCTuned->AddVariable("anglekal",&TMVAvarsUCTuned[0]);
167  fReaderUCTuned->AddVariable("kaldiry",&TMVAvarsUCTuned[1]);
168  fReaderUCTuned->AddVariable("tracklen",&TMVAvarsUCTuned[2]);
169  fReaderUCTuned->AddVariable("nhitcos/nhit",&TMVAvarsUCTuned[3]);
170  fReaderUCTuned->AddVariable("cdirscore",&TMVAvarsUCTuned[4]);
171  fReaderUCTuned->AddVariable("max(vty,endy)",&TMVAvarsUCTuned[5]);
172  fReaderUCTuned->AddVariable("cvncosmic",&TMVAvarsUCTuned[6]);
173  fReaderUCTuned->AddVariable("hadEPerNHit",&TMVAvarsUCTuned[7]);
174  fReaderUCTuned->AddVariable("slcCalEPerNHit",&TMVAvarsUCTuned[8]);
175  fReaderUCTuned->AddVariable("trkEPerNHit",&TMVAvarsUCTuned[9]);
176  fReaderUCTuned->AddVariable("vetoangl",&TMVAvarsUCTuned[10]);
177  fReaderUCTuned->AddVariable("scatt",&TMVAvarsUCTuned[11]);
178  fReaderUCTuned->AddSpectator("anglecos",&TMVAvarsUCTuned[12]);
179  fReaderUCTuned->AddSpectator("anglekal",&TMVAvarsUCTuned[13]);
180  fReaderUCTuned->AddSpectator("pngptp",&TMVAvarsUCTuned[14]);
181  fReaderUCTuned->AddSpectator("cosdiry",&TMVAvarsUCTuned[15]);
182  fReaderUCTuned->AddSpectator("scatt",&TMVAvarsUCTuned[16]);
183  fReaderUCTuned->AddSpectator("scatt/tracklen",&TMVAvarsUCTuned[17]);
184  fReaderUCTuned->AddSpectator("slcT",&TMVAvarsUCTuned[18]);
185  fReaderUCTuned->AddSpectator("remid",&TMVAvarsUCTuned[19]);
186  fReaderUCTuned->AddSpectator("oldbdt",&TMVAvarsUCTuned[20]);
187  fReaderUCTuned->AddSpectator("cvnnumu",&TMVAvarsUCTuned[21]);
188  fReaderUCTuned->AddSpectator("osc",&TMVAvarsUCTuned[22]);
189  fReaderUCTuned->AddSpectator("oscVal1",&TMVAvarsUCTuned[23]);
190  fReaderUCTuned->AddSpectator("oscVal2",&TMVAvarsUCTuned[24]);
191  fReaderUCTuned->AddSpectator("oscVal3",&TMVAvarsUCTuned[25]);
192  fReaderUCTuned->AddSpectator("oscVal4",&TMVAvarsUCTuned[26]);
193  fReaderUCTuned->AddSpectator("oscVal5",&TMVAvarsUCTuned[27]);
194  fReaderUCTuned->AddSpectator("nhit",&TMVAvarsUCTuned[28]);
195  fReaderUCTuned->AddSpectator("trueE",&TMVAvarsUCTuned[29]);
196  fReaderUCTuned->AddSpectator("visE",&TMVAvarsUCTuned[30]);
197  fReaderUCTuned->AddSpectator("ntracks3d",&TMVAvarsUCTuned[31]);
198  fReaderUCTuned->AddSpectator("npngs3d",&TMVAvarsUCTuned[32]);
199  fReaderUCTuned->AddSpectator("nprotons",&TMVAvarsUCTuned[33]);
200  fReaderUCTuned->AddSpectator("nneutrons",&TMVAvarsUCTuned[34]);
201  fReaderUCTuned->AddSpectator("nppions",&TMVAvarsUCTuned[35]);
202  fReaderUCTuned->AddSpectator("nnpions",&TMVAvarsUCTuned[36]);
203  fReaderUCTuned->AddSpectator("npions0",&TMVAvarsUCTuned[37]);
204  fReaderUCTuned->AddSpectator("nmuons",&TMVAvarsUCTuned[38]);
205  fReaderUCTuned->AddSpectator("namuons",&TMVAvarsUCTuned[39]);
206  fReaderUCTuned->AddSpectator("numuE",&TMVAvarsUCTuned[40]);
207  fReaderUCTuned->AddSpectator("calccE",&TMVAvarsUCTuned[41]);
208  fReaderUCTuned->AddSpectator("trkccE",&TMVAvarsUCTuned[42]);
209  fReaderUCTuned->AddSpectator("angleE",&TMVAvarsUCTuned[43]);
210  fReaderUCTuned->AddSpectator("recomuonE",&TMVAvarsUCTuned[44]);
211  fReaderUCTuned->AddSpectator("recohadtrkE",&TMVAvarsUCTuned[45]);
212  fReaderUCTuned->AddSpectator("hadcalE",&TMVAvarsUCTuned[46]);
213  fReaderUCTuned->AddSpectator("hadtrkE",&TMVAvarsUCTuned[47]);
214  fReaderUCTuned->AddSpectator("truemuonE",&TMVAvarsUCTuned[48]);
215  fReaderUCTuned->AddSpectator("truegoodmuon",&TMVAvarsUCTuned[49]);
216  fReaderUCTuned->AddSpectator("hadclustncalhits",&TMVAvarsUCTuned[50]);
217  fReaderUCTuned->AddSpectator("hadclustncontplanes",&TMVAvarsUCTuned[51]);
218  fReaderUCTuned->AddSpectator("hadclustfirstplane",&TMVAvarsUCTuned[52]);
219  fReaderUCTuned->AddSpectator("hadclustlastplane",&TMVAvarsUCTuned[53]);
220  fReaderUCTuned->AddSpectator("hadclustncellsedge",&TMVAvarsUCTuned[54]);
221  fReaderUCTuned->AddSpectator("hadclustcalE",&TMVAvarsUCTuned[55]);
222  fReaderUCTuned->AddSpectator("hadclustmeanposX",&TMVAvarsUCTuned[56]);
223  fReaderUCTuned->AddSpectator("hadclustmeanposY",&TMVAvarsUCTuned[57]);
224  fReaderUCTuned->AddSpectator("hadclustmeanposZ",&TMVAvarsUCTuned[58]);
225  fReaderUCTuned->AddSpectator("tracklen2nd",&TMVAvarsUCTuned[59]);
226  fReaderUCTuned->AddSpectator("trkE2nd",&TMVAvarsUCTuned[60]);
227  fReaderUCTuned->AddSpectator("angle2nd",&TMVAvarsUCTuned[61]);
228  fReaderUCTuned->AddSpectator("type",&TMVAvarsUCTuned[62]);
229  fReaderUCTuned->AddSpectator("vtx",&TMVAvarsUCTuned[63]);
230  fReaderUCTuned->AddSpectator("vty",&TMVAvarsUCTuned[64]);
231  fReaderUCTuned->AddSpectator("vtz",&TMVAvarsUCTuned[65]);
232  fReaderUCTuned->AddSpectator("endx",&TMVAvarsUCTuned[66]);
233  fReaderUCTuned->AddSpectator("endy",&TMVAvarsUCTuned[67]);
234  fReaderUCTuned->AddSpectator("endz",&TMVAvarsUCTuned[68]);
235  fReaderUCTuned->AddSpectator("vetoangl",&TMVAvarsUCTuned[69]);
236  fReaderUCTuned->AddSpectator("run",&TMVAvarsUCTuned[70]);
237  fReaderUCTuned->AddSpectator("subrun",&TMVAvarsUCTuned[71]);
238  fReaderUCTuned->AddSpectator("evt",&TMVAvarsUCTuned[72]);
239  fReaderUCTuned->AddSpectator("subevt",&TMVAvarsUCTuned[73]);
240  fReaderUCTuned->BookMVA("UCTuned_BDTA",pidpath+"UCTuned_prod3_BDTA_700.weights.xml");
241 
242  fReaderTAP1 = new TMVA::Reader;
243  fReaderTAP1->AddVariable("anglekal",&TMVAvarsTA[0]);
244  fReaderTAP1->AddVariable("kaldiry",&TMVAvarsTA[1]);
245  fReaderTAP1->AddVariable("tracklen",&TMVAvarsTA[2]);
246  fReaderTAP1->AddVariable("max(vty,endy)",&TMVAvarsTA[3]);
247  fReaderTAP1->AddVariable("cvncosmic",&TMVAvarsTA[4]);
248  fReaderTAP1->AddVariable("min(cosfwdcell+cosbakcell,kalfwdcell+kalbakcell)",&TMVAvarsTA[5]);
249  fReaderTAP1->AddVariable("nhitkal/nhit",&TMVAvarsTA[6]);
250  fReaderTAP1->BookMVA("TAP1_BDTA",pidpath+"TAP1_BDTA.weights.xml");
251 
252  fReaderTAP2 = new TMVA::Reader;
253  fReaderTAP2->AddVariable("anglekal",&TMVAvarsTA[0]);
254  fReaderTAP2->AddVariable("kaldiry",&TMVAvarsTA[1]);
255  fReaderTAP2->AddVariable("tracklen",&TMVAvarsTA[2]);
256  fReaderTAP2->AddVariable("max(vty,endy)",&TMVAvarsTA[3]);
257  fReaderTAP2->AddVariable("cvncosmic",&TMVAvarsTA[4]);
258  fReaderTAP2->AddVariable("min(cosfwdcell+cosbakcell,kalfwdcell+kalbakcell)",&TMVAvarsTA[5]);
259  fReaderTAP2->AddVariable("nhitkal/nhit",&TMVAvarsTA[6]);
260  fReaderTAP2->BookMVA("TAP2_BDTA",pidpath+"TAP2_BDTA.weights.xml");
261 
262  fReaderTAP3 = new TMVA::Reader;
263  fReaderTAP3->AddVariable("anglekal",&TMVAvarsTA[0]);
264  fReaderTAP3->AddVariable("kaldiry",&TMVAvarsTA[1]);
265  fReaderTAP3->AddVariable("tracklen",&TMVAvarsTA[2]);
266  fReaderTAP3->AddVariable("max(vty,endy)",&TMVAvarsTA[3]);
267  fReaderTAP3->AddVariable("cvncosmic",&TMVAvarsTA[4]);
268  fReaderTAP3->AddVariable("min(cosfwdcell+cosbakcell,kalfwdcell+kalbakcell)",&TMVAvarsTA[5]);
269  fReaderTAP3->AddVariable("nhitkal/nhit",&TMVAvarsTA[6]);
270  fReaderTAP3->BookMVA("TAP3_BDTA",pidpath+"TAP3_BDTA.weights.xml");
271 
272  // --- Prod5 FHC Period 1.
273  fReaderProd5_Per1 = new TMVA::Reader;
274  fReaderProd5_Per1 -> AddVariable( "min(CosFwdCell+CosBakCell,KalTrkFwdCell+KalTrkBakCell)", &TMVAvarsProd5[0 ]);
275  fReaderProd5_Per1 -> AddVariable( "KalTrkCosNumi" , &TMVAvarsProd5[1 ]);
276  fReaderProd5_Per1 -> AddVariable( "KalTrkLen" , &TMVAvarsProd5[2 ]);
277  fReaderProd5_Per1 -> AddVariable( "max(KalTrkStY,KalTrkEndY)" , &TMVAvarsProd5[3 ]);
278  fReaderProd5_Per1 -> AddVariable( "cos(KalTrkDirY)" , &TMVAvarsProd5[4 ]);
279  fReaderProd5_Per1 -> AddVariable( "KalTrkNHit/SliceNHit" , &TMVAvarsProd5[5 ]);
280  fReaderProd5_Per1 -> AddVariable( "KalTrkPtOverP" , &TMVAvarsProd5[6 ]);
281  // And my spectators
282  fReaderProd5_Per1 -> AddSpectator( "IsNu" , &TMVAvarsProd5[7 ]);
283  fReaderProd5_Per1 -> AddSpectator( "IsNuMu" , &TMVAvarsProd5[8 ]);
284  fReaderProd5_Per1 -> AddSpectator( "CVNMuon17" , &TMVAvarsProd5[9 ]);
285  fReaderProd5_Per1 -> AddSpectator( "CVNMuon18" , &TMVAvarsProd5[10]);
286  fReaderProd5_Per1 -> AddSpectator( "KalTrkRemID" , &TMVAvarsProd5[11]);
287  fReaderProd5_Per1 -> BookMVA("BDT", pidpath+"TMVA_MProd5_Kal_FHC_Per1_BDT.weights.xml");
288 
289  // --- Prod5 FHC Period 2.
290  fReaderProd5_Per2 = new TMVA::Reader;
291  fReaderProd5_Per2 -> AddVariable( "min(CosFwdCell+CosBakCell,KalTrkFwdCell+KalTrkBakCell)", &TMVAvarsProd5[0 ]);
292  fReaderProd5_Per2 -> AddVariable( "KalTrkCosNumi" , &TMVAvarsProd5[1 ]);
293  fReaderProd5_Per2 -> AddVariable( "KalTrkLen" , &TMVAvarsProd5[2 ]);
294  fReaderProd5_Per2 -> AddVariable( "max(KalTrkStY,KalTrkEndY)" , &TMVAvarsProd5[3 ]);
295  fReaderProd5_Per2 -> AddVariable( "cos(KalTrkDirY)" , &TMVAvarsProd5[4 ]);
296  fReaderProd5_Per2 -> AddVariable( "KalTrkNHit/SliceNHit" , &TMVAvarsProd5[5 ]);
297  fReaderProd5_Per2 -> AddVariable( "KalTrkPtOverP" , &TMVAvarsProd5[6 ]);
298  // And my spectators
299  fReaderProd5_Per2 -> AddSpectator( "IsNu" , &TMVAvarsProd5[7 ]);
300  fReaderProd5_Per2 -> AddSpectator( "IsNuMu" , &TMVAvarsProd5[8 ]);
301  fReaderProd5_Per2 -> AddSpectator( "CVNMuon17" , &TMVAvarsProd5[9 ]);
302  fReaderProd5_Per2 -> AddSpectator( "CVNMuon18" , &TMVAvarsProd5[10]);
303  fReaderProd5_Per2 -> AddSpectator( "KalTrkRemID" , &TMVAvarsProd5[11]);
304  fReaderProd5_Per2 -> BookMVA("BDT", pidpath+"TMVA_MProd5_Kal_FHC_Per2_BDT.weights.xml");
305 
306  // --- Prod5 FHC High Gain.
307  fReaderProd5_FHC = new TMVA::Reader;
308  fReaderProd5_FHC -> AddVariable( "min(CosFwdCell+CosBakCell,KalTrkFwdCell+KalTrkBakCell)", &TMVAvarsProd5[0 ]);
309  fReaderProd5_FHC -> AddVariable( "KalTrkCosNumi" , &TMVAvarsProd5[1 ]);
310  fReaderProd5_FHC -> AddVariable( "KalTrkLen" , &TMVAvarsProd5[2 ]);
311  fReaderProd5_FHC -> AddVariable( "max(KalTrkStY,KalTrkEndY)" , &TMVAvarsProd5[3 ]);
312  fReaderProd5_FHC -> AddVariable( "cos(KalTrkDirY)" , &TMVAvarsProd5[4 ]);
313  fReaderProd5_FHC -> AddVariable( "KalTrkNHit/SliceNHit" , &TMVAvarsProd5[5 ]);
314  fReaderProd5_FHC -> AddVariable( "KalTrkPtOverP" , &TMVAvarsProd5[6 ]);
315  // And my spectators
316  fReaderProd5_FHC -> AddSpectator( "IsNu" , &TMVAvarsProd5[7 ]);
317  fReaderProd5_FHC -> AddSpectator( "IsNuMu" , &TMVAvarsProd5[8 ]);
318  fReaderProd5_FHC -> AddSpectator( "CVNMuon17" , &TMVAvarsProd5[9 ]);
319  fReaderProd5_FHC -> AddSpectator( "CVNMuon18" , &TMVAvarsProd5[10]);
320  fReaderProd5_FHC -> AddSpectator( "KalTrkRemID" , &TMVAvarsProd5[11]);
321  fReaderProd5_FHC -> BookMVA("BDT", pidpath+"TMVA_MProd5_Kal_FHC_High_BDT.weights.xml");
322 
323  // --- Prod5 RHC High Gain.
324  fReaderProd5_RHC = new TMVA::Reader;
325  fReaderProd5_RHC -> AddVariable( "min(CosFwdCell+CosBakCell,KalTrkFwdCell+KalTrkBakCell)", &TMVAvarsProd5[0 ]);
326  fReaderProd5_RHC -> AddVariable( "KalTrkCosNumi" , &TMVAvarsProd5[1 ]);
327  fReaderProd5_RHC -> AddVariable( "KalTrkLen" , &TMVAvarsProd5[2 ]);
328  fReaderProd5_RHC -> AddVariable( "max(KalTrkStY,KalTrkEndY)" , &TMVAvarsProd5[3 ]);
329  fReaderProd5_RHC -> AddVariable( "cos(KalTrkDirY)" , &TMVAvarsProd5[4 ]);
330  fReaderProd5_RHC -> AddVariable( "KalTrkNHit/SliceNHit" , &TMVAvarsProd5[5 ]);
331  fReaderProd5_RHC -> AddVariable( "KalTrkPtOverP" , &TMVAvarsProd5[6 ]);
332  // And my spectators
333  fReaderProd5_RHC -> AddSpectator( "IsNu" , &TMVAvarsProd5[7 ]);
334  fReaderProd5_RHC -> AddSpectator( "IsNuMu" , &TMVAvarsProd5[8 ]);
335  fReaderProd5_RHC -> AddSpectator( "CVNMuon17" , &TMVAvarsProd5[9 ]);
336  fReaderProd5_RHC -> AddSpectator( "CVNMuon18" , &TMVAvarsProd5[10]);
337  fReaderProd5_RHC -> AddSpectator( "KalTrkRemID" , &TMVAvarsProd5[11]);
338  fReaderProd5_RHC -> BookMVA("BDT", pidpath+"TMVA_MProd5_Kal_RHC_High_BDT.weights.xml");
339 
340  }
341 
342  //-------------------------------------------------------------------
344  {
345  //Init();
346 
348  evt.getByToken(fSlicerToken,slicevec);
349 
350  if(slicevec->empty()) {
351  mf::LogWarning ("No Slices")<<"No Slices in the input file";
352  return;
353  }
354 
355  art::FindManyP<rb::Track> trackAssnList(slicevec, evt, fTrackLabel);
356  art::FindManyP<rb::Track> costrackAssnList(slicevec, evt, fCosmicTrackLabel);
357  std::unique_ptr< std::vector<cosrej::CosRejObj> > outputObjects(new std::vector<cosrej::CosRejObj>);
358  std::unique_ptr< art::Assns<cosrej::CosRejObj, rb::Cluster> > assoc(new art::Assns<cosrej::CosRejObj, rb::Cluster>);
359  std::unique_ptr< std::vector<cosrej::TrkCntObj> > outTrkObjects(new std::vector<cosrej::TrkCntObj>);
360  std::unique_ptr< art::Assns<cosrej::TrkCntObj, rb::Track> > trkAssoc(new art::Assns<cosrej::TrkCntObj, rb::Track>);
361 
362  // Karl::Get the horn current.
364  if ( !evt.isRealData() ) {
365  evt.getByLabel(fGeneratorLabel, spillPot);
366  } else {
367  evt.getByLabel(fNuMIBeamLabel , spillPot);
368  }
369  bool isRHC = false;
370  if ( spillPot.failedToGet() ) {
371  mf::LogError("CosRej") << "Spill Data not found, aborting without horn current information";
372  abort();
373  } else {
374  isRHC = spillPot->isRHC;
375  }
376  mf::LogDebug("CosRej") << "Horn Current set to isRHC: "<< isRHC;
377 
378  tf::TimingFitAlg tfalg;
379 
380  for(size_t i = 0; i < slicevec->size(); ++i){
381 
382  art::Ptr<rb::Cluster> slice(slicevec,i);
383 
384  if(slice->IsNoise()) {
385  continue;
386  }
387 
388  cosrej::CosRejObj cosrejobj; // initialization in constructor
389 
390  unsigned int bestIdx = 999; // for index in list of the highest ReMId PID track
391  TVector3 minp, maxp; // minimum and maximum extent vectors for slice
392  int nhadtracks = 0;
393  double fwddist = 0, bakdist = 0, scatt = 0;
394  int fwdcell = 0, bakcell = 0;
395  bool badtrack = false; // check for messed up tracks; a few kalman tracks have nans for directions
396 
397  const std::vector< art::Ptr<rb::Track> > cosTracks = costrackAssnList.at(i);
398 
399  minp = slice->MinXYZ();
400  maxp = slice->MaxXYZ();
401 
402  if(!costrackAssnList.isValid()) {
403  mf::LogError("CosRej") << "CosRej: No Cosmic Tracks! Do not go any further!";
404  outputObjects->push_back(cosrejobj);
405  util::CreateAssn(evt,*(outputObjects.get()),slice,*(assoc.get()));
406  return;
407  }
408 
409  if(cosTracks.size()==1) { // for straight line fit version of Cosmic Track, should only be one track
410  const art::Ptr<rb::Track> track = cosTracks[0];
411  double angleOfBestTrack = fCosRejFxs->getAngle(track->Dir());
412  if(angleOfBestTrack>-6&&angleOfBestTrack<2) cosrejobj.SetAngleCos(angleOfBestTrack);
413 
414  TVector3 start = track->Start();
415  TVector3 end = track->Stop();
416  TVector3 dirbak = -track->Dir().Unit();
417  TVector3 dirfor = (track->TrajectoryPoint(track->NTrajectoryPoints()-1) - track->TrajectoryPoint(track->NTrajectoryPoints()-2)).Unit();
418 
419  cosrejobj.SetCosFwdDist(livegeom->ProjectedLiveDistanceToEdge(end,dirfor));
420  cosrejobj.SetCosBakDist(livegeom->ProjectedLiveDistanceToEdge(start,dirbak));
421  cosrejobj.SetCosFwdCell(livegeom->ProjectedLiveCellsToEdge(end,dirfor));
422  cosrejobj.SetCosBakCell(livegeom->ProjectedLiveCellsToEdge(start,dirbak));
423 
424  if(geom->DetId() == novadaq::cnv::kNEARDET) {
425  TVector3 dirfor = (track->TrajectoryPoint(track->NTrajectoryPoints()-1) - track->TrajectoryPoint(track->NTrajectoryPoints()-2)).Unit();
426 
427  cosrejobj.SetCosFwdSteel(livegeom->ProjectedSteelDist(track->Stop(), dirfor));
428  cosrejobj.SetCosBakSteel(livegeom->ProjectedSteelDist(track->Start(), -(track->Dir())));
429  cosrejobj.SetCosFwdAir(livegeom->ProjectedAirDist(track->Stop(), dirfor));
430  cosrejobj.SetCosBakAir(livegeom->ProjectedAirDist(track->Start(), -(track->Dir())));
431  cosrejobj.SetCosYPosAtTrans(livegeom->YPositionAtTransition(track->Start(), track->Dir()));
432  cosrejobj.SetCosFwdCellND(livegeom->FullNDProjectedCells(track->Stop(), dirfor));
433  cosrejobj.SetCosBakCellND(livegeom->FullNDProjectedCells(track->Start(), -(track->Dir())));
434  }
435 
436  TMVAvarsUCTuned[3] = (1.0*track->NCell() / (1.0*slice->NCell()));
437  double slope = 0, chisq = 0, chisqdiff = 0;
438  if(track->AllCells().size()>2)fCosRejFxs->getFits(track,slope,chisq,chisqdiff);
439  cosrejobj.SetCosSlope(slope);
440  cosrejobj.SetCosChisq(chisq);
441  cosrejobj.SetCosChiDiff(chisqdiff);
442  tf::TimingFitResult tfres = tfalg.HoughFit(track.get());
443  cosrejobj.SetCFitSpeed(tfres.bestInvSpeed);
444  cosrejobj.SetCDirScore(tfres.fwdScore-tfres.revScore);
445  cosrejobj.SetCScoreDiff(tfres.bestScore-std::max(tfres.fwdScore,tfres.revScore)); // Chris's timing fit (Hough)
446  cosrejobj.SetCosHitRatio(1.0*track->NCell()/(1.0*slice->NCell()));
447  TMVAvarsUCTuned[4] = cosrejobj.CDirScore();
448  TMVAvarsUCTuned[10] = fabs(angleOfBestTrack)*(track->Dir().Y() + 1.0);
449 
450  if(bt->HaveTruthInfo()) {
451  art::PtrVector<rb::CellHit> trkhits = track->AllCells();
452  const std::vector<const sim::Particle*> parts = bt->HitsToParticle(trkhits);
453 
454  if(!parts.empty()){
455  const TVector3 trueDir = parts[0]->Momentum().Vect().Unit();
456  cosrejobj.SetCosThetaTrue(trueDir.Dot(track->Dir()));
457  }
458  } // end cosmic truth
459  } // end cosmics
460 
461  if(!trackAssnList.isValid()) {
462  mf::LogWarning("CosRej") << "CosRej: No Kalman Tracks!";
463  outputObjects->push_back(cosrejobj);
464  util::CreateAssn(evt,*(outputObjects.get()),slice,*(assoc.get()));
465  continue;
466  }
467 
468  const std::vector< art::Ptr<rb::Track> > sliceTracks = trackAssnList.at(i);
469 
470  //
471  // loop over tracks and fill per-track variables
472  //
473  for(size_t iTrack = 0; iTrack < sliceTracks.size(); ++iTrack)
474  {
475  cosrej::TrkCntObj trkcntobj;
476  art::Ptr<rb::Track> track = sliceTracks[iTrack];
477 
478  // tracks have to have more than 2 hits
479  if(track->NTrajectoryPoints() < 2){
480  outTrkObjects->push_back(trkcntobj);
481  util::CreateAssn(evt,*(outTrkObjects.get()),track,*(trkAssoc.get()));
482  continue;
483  }
484 
485  // get track direction to check validity
486  double tdx = track->Dir().X();
487  double tdy = track->Dir().Y();
488  double tdz = track->Dir().Z();
489 
490  // ignore bad track
491  if(!(tdx>=-1&&tdx<=1&&tdy>=-1&&tdy<=1&&tdz>=-1&&tdz<=1)){
492  outTrkObjects->push_back(trkcntobj);
493  util::CreateAssn(evt,*(outTrkObjects.get()),track,*(trkAssoc.get()));
494  continue;
495  }
496 
497  TVector3 start = track->Start();
498  TVector3 end = track->Stop();
499  TVector3 dirbak = -track->Dir().Unit();
500  TVector3 dirfor = (track->TrajectoryPoint(track->NTrajectoryPoints()-1) - track->TrajectoryPoint(track->NTrajectoryPoints()-2)).Unit();
501 
502  trkcntobj.SetTrkFwdCell(livegeom->ProjectedLiveCellsToEdge(end, dirfor));
503  trkcntobj.SetTrkBakCell(livegeom->ProjectedLiveCellsToEdge(start, dirbak));
504 
505  // if far detector, fill this way
507  trkcntobj.SetTrkLenInAct(track->TotalLength());
508 
509  // if near detector, fill this way
511  {
512  trkcntobj.SetTrkFwdCellND(livegeom->FullNDProjectedCells(end, dirfor));
513  trkcntobj.SetTrkBakCellND(livegeom->FullNDProjectedCells(start, dirbak));
514 
515  trkcntobj.SetTrkLenInAct(track->DistanceFromStart(livegeom->MCFrontZ()));
516  trkcntobj.SetTrkLenInCat(track->DistanceFromEnd(livegeom->MCFrontZ()));
517 
518  // get y opsition at transition plane
519  double xtran, ytran;
520  track->InterpolateXY(livegeom->MCFrontZ(), xtran, ytran);
521  trkcntobj.SetTrkYPosAtTrans(ytran);
522 
523  // Other containment variables
524  trkcntobj.SetTrkFwdSteel(livegeom->ProjectedSteelDist(end, dirfor));
525  trkcntobj.SetTrkBakSteel(livegeom->ProjectedSteelDist(start, dirbak));
526  trkcntobj.SetTrkFwdAir(livegeom->ProjectedAirDist(end, dirfor));
527  trkcntobj.SetTrkBakAir(livegeom->ProjectedAirDist(start, dirbak));
528  }
529 
530  // Other containment variables
531  trkcntobj.SetVtxDist(std::min(livegeom->DistToEdgeXY(start),livegeom->DistToEdgeZ(start)));
533  trkcntobj.SetTrkFwdDist(livegeom->ProjectedLiveDistanceToEdge(end,dirfor));
534  trkcntobj.SetTrkBakDist(livegeom->ProjectedLiveDistanceToEdge(start,dirbak));
535 
536  outTrkObjects->push_back(trkcntobj);
537  util::CreateAssn(evt,*(outTrkObjects.get()),track,*(trkAssoc.get()));
538  } // end of loop over tracks on slice
539 
541  outputObjects->push_back(cosrejobj);
542  util::CreateAssn(evt,*(outputObjects.get()),slice,*(assoc.get()));
543  continue;
544  }
545 
546  art::FindOneP<remid::ReMId> remidVec(sliceTracks, evt, fPIDLabel);
547 
548  if(fUseLongestTrack) {
549  bestIdx = 0;
550  if(sliceTracks.size() == 0) {
551  outputObjects->push_back(cosrejobj);
552  continue;
553  }
554  const art::Ptr<rb::Track> track = sliceTracks[0];
555  if(!track->Is3D()) {
556  outputObjects->push_back(cosrejobj);
557  continue;
558  }
559  cosrejobj.SetNTracks3D(sliceTracks.size());
560  }
561  else if(!remidVec.isValid()) {
562  mf::LogWarning("CosRej") << "CosRej: No ReMId information!";
563  outputObjects->push_back(cosrejobj);
564  util::CreateAssn(evt,*(outputObjects.get()),slice,*(assoc.get()));
565  continue;
566  }
567  else {
568  bestIdx = remid::HighestPIDTrack(sliceTracks, fPIDLabel, evt);
569  for(size_t j = 0; j < sliceTracks.size(); ++j){ // loop over Kalman tracks
570  const art::Ptr<rb::Track> track = sliceTracks[j];
571  art::Ptr<remid::ReMId> trackRemid = remidVec.at(j);
572  if(!track->Is3D() || trackRemid->Value() < 0) continue;
573  if(cosrejobj.NTracks3D()==-5) cosrejobj.SetNTracks3D(1);
574  else cosrejobj.SetNTracks3D(cosrejobj.NTracks3D()+1);
575  if(trackRemid->Value() < 0.5) nhadtracks++;
576  }
577  }
578 
579  // Main condition
580  if(cosrejobj.NTracks3D()>0 && bestIdx!=999) {
581 
582  const art::Ptr<rb::Track> track = sliceTracks[bestIdx];
583 
584  art::PtrVector<rb::CellHit> trkhits = track->AllCells();
585 
586  if(cosrejobj.NTracks3D()>1) {
587  // Getting the 2nd highest ReMId PID track
588  unsigned int bestIdx2nd = 999;
589  double highestPID2nd = -999.0;
590  double bestTrkLen2nd = 0.0;
591  for(size_t j = 0; j < sliceTracks.size(); ++j){ // loop over Kalman tracks
592  if(j == bestIdx) continue;
593  art::Ptr<remid::ReMId> trackRemid2nd = remidVec.at(j);
594  if(trackRemid2nd->Value() < highestPID2nd) continue;
595  if(trackRemid2nd->Value() == highestPID2nd){
596  if(sliceTracks[j]->TotalLength() < bestTrkLen2nd) continue;
597  }
598  highestPID2nd = trackRemid2nd->Value();
599  bestTrkLen2nd = sliceTracks[j]->TotalLength();
600  bestIdx2nd = j;
601  }//end loop over kalman trks for 2nd best ReMId
602  const art::Ptr<rb::Track> track2nd = sliceTracks[bestIdx2nd];
603  cosrejobj.SetTrackLenKal2nd(track2nd->TotalLength());
604  cosrejobj.SetTrackE2nd(track2nd->CalorimetricEnergy());
605  double angleOf2ndBestTrack = fCosRejFxs->getAngle(track2nd->Dir());
606  if((angleOf2ndBestTrack>-6) && (angleOf2ndBestTrack<2)) cosrejobj.SetAngleKal2nd(angleOf2ndBestTrack);
607  TVector3 firstDir = track->Dir().Unit();
608  TVector3 secondDir = track2nd->Dir().Unit();
609  double angleBtwKalTrks = (firstDir*secondDir);
610  cosrejobj.SetAngleBtwKalTrks(angleBtwKalTrks);
611  }
612  else{
613  cosrejobj.SetTrackLenKal2nd(-5.0);
614  cosrejobj.SetTrackE2nd(-5.0);
615  cosrejobj.SetAngleKal2nd(-5.0);
616  cosrejobj.SetAngleBtwKalTrks(-5.0);
617  }
618 
619  if(bt->HaveTruthInfo()) {
620  const std::vector<const sim::Particle*> parts = bt->HitsToParticle(trkhits);
621  const std::vector< cheat::TrackIDE > trids = bt->HitsToTrackIDE(trkhits);
622 
623  if(!parts.empty()){
624  const TVector3 trueDir = parts[0]->Momentum().Vect().Unit();
625  cosrejobj.SetKalThetaTrue(trueDir.Dot(track->Dir()));
626  cosrejobj.SetPDGBest(parts[0]->PdgCode());
627  }
628  }
629 
630  double tdx = track->Dir().X();
631  double tdy = track->Dir().Y();
632  double tdz = track->Dir().Z();
633 
634  if(!(tdx>=-1&&tdx<=1&&tdy>=-1&&tdy<=1&&tdz>=-1&&tdz<=1)) badtrack = true;
635 
636  if(badtrack) continue;
637 
638  TVector3 start = track->Start();
639  TVector3 end = track->Stop();
640  TVector3 dirbak = -track->Dir().Unit();
641  TVector3 dirfor = (track->TrajectoryPoint(track->NTrajectoryPoints()-1) - track->TrajectoryPoint(track->NTrajectoryPoints()-2)).Unit();
642 
643  cosrejobj.SetKalFwdDist(livegeom->ProjectedLiveDistanceToEdge(end,dirfor));
644  cosrejobj.SetKalBakDist(livegeom->ProjectedLiveDistanceToEdge(start,dirbak));
645  cosrejobj.SetKalFwdCell(livegeom->ProjectedLiveCellsToEdge(end,dirfor));
646  cosrejobj.SetKalBakCell(livegeom->ProjectedLiveCellsToEdge(start,dirbak));
647 
648  cosrejobj.SetVtxDist(std::min(livegeom->DistToEdgeXY(start),livegeom->DistToEdgeZ(start)));
650 
651  if(geom->DetId() == novadaq::cnv::kNEARDET) {
652  cosrejobj.SetKalFwdSteel(livegeom->ProjectedSteelDist(track->Stop(), dirfor));
653  cosrejobj.SetKalBakSteel(livegeom->ProjectedSteelDist(track->Start(), -(track->Dir())));
654  cosrejobj.SetKalFwdAir(livegeom->ProjectedAirDist(track->Stop(), dirfor));
655  cosrejobj.SetKalBakAir(livegeom->ProjectedAirDist(track->Start(), -(track->Dir())));
656 
657  double xtran, ytran;
658  track->InterpolateXY(livegeom->MCFrontZ(), xtran, ytran);
659  cosrejobj.SetKalYPosAtTrans(ytran);
660  cosrejobj.SetKalFwdCellND(livegeom->FullNDProjectedCells(track->Stop(), dirfor));
661  cosrejobj.SetKalBakCellND(livegeom->FullNDProjectedCells(track->Start(), -(track->Dir())));
662  }
663 
664  if(std::min(fwddist,bakdist) < cosrejobj.MinDist()) cosrejobj.SetMinDist(std::min(fwddist,bakdist));
665  if(std::min(fwdcell,bakcell) < cosrejobj.MinCell()) cosrejobj.SetMinCell(std::min(fwdcell,bakcell));
666 
667  cosrejobj.SetERatio(track->TotalGeV()/slice->TotalGeV());
668  scatt = fCosRejFxs->getScatt(track);
669  if(track->TotalLength()>0) cosrejobj.SetScatt(scatt/track->TotalLength());
670 
671  double ferpars[4];
672  bool fdef = 0;
673  bool fFlag = fCosRejFxs->FScatterEstim(track,fdef,ferpars[0],ferpars[1],ferpars[2],ferpars[3]); // Fernandas scattering energy estimate
674  if(fFlag==1) {
675  cosrejobj.SetFScattExt(ferpars[0]);
676  cosrejobj.SetFScattSig(ferpars[1]);
677  cosrejobj.SetFScattMax(ferpars[2]);
678  cosrejobj.SetFScattSum(ferpars[3]);
679  }
680  double slope = 0, chisq = 0, chisqdiff = 0;
681  if(track->AllCells().size()>2)fCosRejFxs->getFits(track,slope,chisq,chisqdiff);
682  cosrejobj.SetKalSlope(slope);
683  cosrejobj.SetKalChisq(chisq);
684  cosrejobj.SetKalChiDiff(chisqdiff);
685 
686  tf::TimingFitResult tfres = tfalg.HoughFit(track.get());
687  cosrejobj.SetKFitSpeed(tfres.bestInvSpeed);
688  cosrejobj.SetKDirScore(tfres.fwdScore-tfres.revScore);
689  cosrejobj.SetKScoreDiff(tfres.bestScore-std::max(tfres.fwdScore,tfres.revScore));
690  cosrejobj.SetStartAct(fCosRejFxs->getActivity(track->AllCells(),slice->AllCells()));
691  cosrejobj.SetNHadTracks(nhadtracks);
692 
693  double angleOfBestTrack = fCosRejFxs->getAngle(sliceTracks[bestIdx]->Dir());
694  if(angleOfBestTrack>-6&&angleOfBestTrack<2) cosrejobj.SetAngleKal(angleOfBestTrack);
695 
696  cosrejobj.SetTrackLenKal(track->TotalLength());
697 
698  if(geom->DetId() == novadaq::cnv::kFARDET) {
699  const art::Ptr<rb::Track> track = sliceTracks[bestIdx];
700 
701  float fslcCalEPerNHit = 0.0;
702  float fhadEPerNHit = 0.0;
703  int frun = evt.run();
704 
705  art::FindManyP<cvn::Result> fmCVN(slicevec, evt, fCVNLabel); // copied from CAFMAker
706  if(fmCVN.isValid()) {
707  const std::vector<art::Ptr<cvn::Result>> results = fmCVN.at(i);
708  if(!results.empty()) {
709 
710  // --- Uncontained BDT
711  if(slice->NCell() > 0) fslcCalEPerNHit = ((slice->CalorimetricEnergy()) / (1.78*slice->NCell()));
712  if(track->NCell() > 0) cosrejobj.SetTrkEPerNHit((track->CalorimetricEnergy())/(track->NCell()));
713  if((slice->NCell()>0) && (track->NCell()>0)) fhadEPerNHit = (slice->CalorimetricEnergy() - track->CalorimetricEnergy())/(slice->NCell() - track->NCell());
714 
715  TMVAvarsUCTuned[0] = cosrejobj.AngleKal();
716  TMVAvarsUCTuned[1] = track->Dir().Y();
717  TMVAvarsUCTuned[2] = cosrejobj.TrackLenKal(); //track->TotalLength();
718  // [3] and [4] above.
719  TMVAvarsUCTuned[5] = std::max(track->Start().Y(),track->Stop().Y());
720  TMVAvarsUCTuned[6] = results[0]->fOutput[3]; // 2020 CVN has Cosmic at output 3.
721  TMVAvarsUCTuned[7] = fhadEPerNHit;
722  TMVAvarsUCTuned[8] = fslcCalEPerNHit;
723  TMVAvarsUCTuned[9] = cosrejobj.TrkEPerNHit(); //ftrkEPerNHit;
724  // [10] above.
725  TMVAvarsUCTuned[11] = cosrejobj.Scatt(); //track->TotalLength();
726  // --- Set the Uncontained BDT
727  cosrejobj.SetUnconTunedCosPID(fReaderUCTuned->EvaluateMVA("UCTuned_BDTA"));
728 
729  // --- Prod3 Trained BDT
730  TMVAvarsTA[0] = cosrejobj.AngleKal();
731  TMVAvarsTA[1] = track->Dir().Y();;
732  TMVAvarsTA[2] = track->TotalLength();
733  TMVAvarsTA[3] = std::max(track->Start().Y(),track->Stop().Y());
734  TMVAvarsTA[4] = results[0]->fOutput[3]; // 2020 CVN has Cosmic at output 3.
735  TMVAvarsTA[5] = std::min(cosrejobj.CosFwdCell()+cosrejobj.CosBakCell(),cosrejobj.KalFwdCell()+cosrejobj.KalBakCell());
736  TMVAvarsTA[6] = track->NCell() / (float)slice->NCell();
737  // --- Set the Prod3 Trained BDT
738  if (frun > 10000 && frun < 17300) {
739  cosrejobj.SetOldCosPID(fReaderTAP1->EvaluateMVA("TAP1_BDTA"));
740  } else if(frun > 17300 && frun < 20753) {
741  cosrejobj.SetOldCosPID(fReaderTAP2->EvaluateMVA("TAP2_BDTA"));
742  } else if (frun > 20752) {
743  cosrejobj.SetOldCosPID(fReaderTAP3->EvaluateMVA("TAP3_BDTA"));
744  }
745  } // end results not empty
746  } else {
747  mf::LogError("CosRej") << "CVN Data Product " << fCVNLabel << " not found, needed for CosRej2019, aborting.";
748  abort();
749  } // end CVN handle is valid
750 
751  // Karl::This is where we want to evaluate the Prod5 cosmic rejection...
752  // First, fill my BDTVars for this event...
753  TMVAvarsProd5[0] = std::min(cosrejobj.CosFwdCell()+cosrejobj.CosBakCell(),cosrejobj.KalFwdCell()+cosrejobj.KalBakCell());
754  TMVAvarsProd5[1] = track->Dir().Dot( geom->NuMIBeamDirection() ); // kCosNuMi
755  TMVAvarsProd5[2] = track->TotalLength()/100; // kTrkLength
756  TMVAvarsProd5[3] = std::max( track->Start().Y(), track->Stop().Y() )/100; // kMaxKalTrStEnY
757  TMVAvarsProd5[4] = cos( track->Dir().Y() ); // kcosDirY
758  TMVAvarsProd5[5] = track->NCell() / (float)slice->NCell(); // kKalHitRat
759  TMVAvarsProd5[6] = sqrt( 1 - (TMVAvarsProd5[1]*TMVAvarsProd5[1]) ); // kNumuMuonPtP
760  // Then, store the official number.
761  if (isRHC) { // RHC, so use the RHC BDT :O
762  cosrejobj.SetConCosPID ( fReaderProd5_RHC ->EvaluateMVA("BDT") );
763  } else { // FHC, so need to figure out if P1, P2, HighGain
764  if ( frun <= 17139 ) {
765  cosrejobj.SetConCosPID( fReaderProd5_Per1->EvaluateMVA("BDT") );
766  } else if ( frun <= 19746 ) {
767  cosrejobj.SetConCosPID( fReaderProd5_Per2->EvaluateMVA("BDT") );
768  } else {
769  cosrejobj.SetConCosPID( fReaderProd5_FHC ->EvaluateMVA("BDT" ) );
770  }
771  }
772  /*
773  std::cout << "\nLooking at Evt " << evt.event() << ", Slice " << i << "."
774  << I got: " << cosrejobj.ConCosPID() << ", Comp to " << cosrejobj.OldCosPID()
775  << "\n\t TMVAvarsProd5[0]: MinCellSum = " << TMVAvarsProd5[0]
776  << "\n\t TMVAvarsProd5[1]: kCosNuMi = " << TMVAvarsProd5[1]
777  << "\n\t TMVAvarsProd5[2]: kTrkLength = " << TMVAvarsProd5[2]
778  << "\n\t TMVAvarsProd5[3]: kMaxKalTrStEnY = " << TMVAvarsProd5[3]
779  << "\n\t TMVAvarsProd5[4]: kcosDirY = " << TMVAvarsProd5[4]
780  << "\n\t TMVAvarsProd5[5]: kKalHitRat = " << TMVAvarsProd5[5] << " .... " << track->NCell() / slice->NCell()
781  << " .... TMVAvarsUCTuned[3] = " << TMVAvarsUCTuned[3] << " .... " << TMVAvarsTA[6]
782  << "\n\t TMVAvarsProd5[6]: kNumuMuonPtP = " << TMVAvarsProd5[6]
783  << std::endl;
784  //*/
785  /*
786  fCosRej20 -> Fill( cosrejobj.ConCosPID() );
787  fCosRej20_CorrTA -> Fill( cosrejobj.ConCosPID(), cosrejobj.OldCosPID() );
788  //*/
789  } // end if FarDet
790  } // end if good track
791 
792  outputObjects->push_back(cosrejobj); // write out object with all information
793  util::CreateAssn(evt,*(outputObjects.get()),slice,*(assoc.get()));
794 
795  } //end loop over slices
796 
797  evt.put(std::move(outputObjects));
798  evt.put(std::move(assoc));
799  evt.put(std::move(outTrkObjects));
800  evt.put(std::move(trkAssoc));
801 
802  } //end produce
803 
804 } //end namespace
805 
void SetKalChisq(float kalchisq)
chisq result associated with slope fit. Mostly not useful
Definition: CosRejObj.cxx:165
T max(const caf::Proxy< T > &a, T b)
void SetCDirScore(float cdirscore)
chisq diff between assuming forwards going and backwards going for Cosmic Track
Definition: CosRejObj.cxx:141
size_t NTrajectoryPoints() const
Definition: Track.h:83
bool isRHC
is the beam in antineutrino mode, aka RHC
Definition: SpillData.h:28
void SetKalFwdDist(float kalfwddist)
distance (dist,cm) of best track from end to det wall projected forwards (Fwd), based on Kalman Track...
Definition: CosRejObj.cxx:219
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.
void SetKFitSpeed(float kfitspeed)
inverse speed (ns/cm) of track fit by TimingFit module (negative indicates backwards track) for best ...
Definition: CosRejObj.cxx:123
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
void SetTrkLenInCat(double trklenincat)
Definition: TrkCntObj.h:28
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
Definition: Track.cxx:158
bool getByToken(ProductToken< PROD > const &, Handle< PROD > &result) const
Definition: DataViewImpl.h:462
virtual ~CosRej()
void SetKalSlope(float kalslope)
slope of timing fit; a positive slope indicates track is properly directed
Definition: CosRejObj.cxx:159
double ProjectedSteelDist(TVector3 vertex, TVector3 dir)
TMVA::Reader * fReaderProd5_FHC
float TMVAvarsTA[7]
void SetTrackE2nd(float tracke2nd)
Calorimetric energy of 2nd highest ReMId kalman track.
Definition: CosRejObj.cxx:429
void SetKalBakSteel(float kalbaksteel)
distance (dist,cm) through muon catcher traveled, when projected forwards to detector edge...
Definition: CosRejObj.cxx:249
double DistToEdgeXY(TVector3 vertex)
void SetCosBakAir(float cosbakair)
distance (dist,cm) through air traveled, when projected backwards to detector edge. For ND only, muon catcher air gap. NYI - just use ytrans
Definition: CosRejObj.cxx:207
void SetKalChiDiff(float kalchidiff)
chisq difference of assuming + or - c as slope. A positive chisqdiff means the track is backwards...
Definition: CosRejObj.cxx:363
void SetOldCosPID(float oldcospid)
PID for contained cosmic rejection (Old Alg, Prod 3 training)
Definition: CosRejObj.cxx:98
void SetKScoreDiff(float kscorediff)
chisq diff between best directed fit and free fit (may indicate if things are wonky) for best ReMId K...
Definition: CosRejObj.cxx:147
T sqrt(T number)
Definition: d0nt_math.hpp:156
void SetCosFwdCell(int cosfwdcell)
cells in that line segment instead of distance of line segment
Definition: CosRejObj.cxx:255
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
CosRej(fhicl::ParameterSet const &pset)
Store results of TimingFitAlg.
void SetAngleKal(float angleb)
cos of angle of best Kalman track wrt beam
Definition: CosRejObj.cxx:80
void SetTrkBakCellND(int trkbakcellnd)
Definition: TrkCntObj.h:26
void SetNHadTracks(int nhadtracks)
number of 3D Kalman tracks with ReMId value of < 0.4 (not mu-like, or hadronic-like) ...
Definition: CosRejObj.cxx:117
void SetCosFwdDist(float cosfwddist)
distance (dist,cm) of track from end to det wall projected forwards (Fwd), based on Cosmic Tracker (C...
Definition: CosRejObj.cxx:183
void SetTrackLenKal2nd(float tracklenkal2nd)
Total length of 2nd highest ReMId kalman track.
Definition: CosRejObj.cxx:423
void SetScatt(float scatt)
basic scattering variable; sum of all angular deviations between best kal track trajectory points / t...
Definition: CosRejObj.cxx:333
Definition: event.h:19
float Scatt() const
Definition: CosRejObj.cxx:706
void SetCosFwdAir(float cosfwdair)
distance (dist,cm) through air traveled, when projected forwards to detector edge. For ND only, muon catcher air gap. NYI
Definition: CosRejObj.cxx:189
void SetVtxDist(float vtxdist)
Definition: TrkCntObj.h:30
Algorithm to determine the direction of tracks using timing.
Definition: TimingFitAlg.h:26
void SetCosYPosAtTrans(float cosyposattrans)
Y position at transition to muon catcher, for determining if track went through air gap (ND only) ...
Definition: CosRejObj.cxx:303
void SetTrkBakAir(float trkbakair)
distance (dist,cm) through muon catcher traveled, when projected forwards to detector edge...
Definition: TrkCntObj.h:36
std::string fNuMIBeamLabel
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
void SetTrkFwdAir(float trkfwdair)
cosmic track distance projected backwards from track start
Definition: TrkCntObj.h:34
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
void SetMinDist(float mindist)
minimum projected distance, forwards or backwards, of ANY Kalman track with more than 15 hits ...
Definition: CosRejObj.cxx:315
void SetTrkLenInAct(double trkleninact)
Definition: TrkCntObj.h:27
void SetKalBakCellND(int kalbakcellnd)
note that fwd<->bak values are interchanged if track direction is backwards, including MC...
Definition: CosRejObj.cxx:297
void SetCosBakCell(int cosbakcell)
note that fwd<->bak values are interchanged if track direction is backwards
Definition: CosRejObj.cxx:267
void SetUnconTunedCosPID(float uncontunedcospid)
Tuned PID for uncontained cosmic rejection (TMVA) - Jose&#39;s.
Definition: CosRejObj.cxx:105
int KalFwdCell() const
Definition: CosRejObj.cxx:652
TimingFitResult HoughFit(const rb::Track *track) const
double Value() const
Definition: PID.h:22
int ProjectedLiveCellsToEdge(TVector3 vertex, TVector3 dir)
recommended projection to use
void SetNTracks3D(int ntracks3d)
Definition: CosRejObj.cxx:111
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
TVector3 MaxXYZ() const
Definition: Cluster.cxx:492
void SetTrkYPosAtTrans(float trkyposattrans)
Definition: TrkCntObj.h:29
void SetStartAct(float startact)
measure of activity near best Kalman track start position; start pos is a good measure of vertex ...
Definition: CosRejObj.cxx:339
double MCFrontZ() const
Z position at which the muon catcher starts.
Definition: LiveGeometry.h:132
void produce(art::Event &evt)
virtual double TotalLength() const
Length (cm) of all the track segments.
Definition: Track.cxx:213
float TrkEPerNHit() const
Definition: CosRejObj.cxx:827
bool isRealData() const
bool FScatterEstim(art::Ptr< rb::Track > track, bool &fromTrkEnd, double &angleExt, double &angleSigma, double &angleMax, double &angleSum)
Definition: CosRejFxs.cxx:228
double ProjectedLiveDistanceToEdge(TVector3 vertex, TVector3 dir)
void SetKalFwdCell(int kalfwdcell)
cells in that line segment instead of distance of line segment
Definition: CosRejObj.cxx:279
Far Detector at Ash River, MN.
void SetAngleBtwKalTrks(float anglebtw)
cos of angle of 2nd best Kalman track wrt beam
Definition: CosRejObj.cxx:441
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
void SetTrkFwdCellND(int trkfwdcellnd)
Definition: TrkCntObj.h:24
void SetCosChiDiff(float coschidiff)
same for Cosmic track
Definition: CosRejObj.cxx:369
int PdgCode(Resonance_t res, int Q)
(resonance id, charge) -> PDG code
Result for CVN.
double CalorimetricEnergy(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple estimate of neutrino energy.
Definition: Cluster.cxx:439
void SetTrkFwdSteel(float trkfwdsteel)
distance (dist,cm) through air traveled, when projected forwards to detector edge. For ND only, muon catcher air gap. NYI
Definition: TrkCntObj.h:35
double DistToEdgeZ(TVector3 vertex)
void SetCosThetaTrue(float cosqtrue)
cosine of angle between true particle (bt most contributing) and reco track dir; use to tell if track...
Definition: CosRejObj.cxx:345
void getFits(const art::Ptr< rb::Track > track, double &slope, double &chisq, double &chisqdiff)
Definition: CosRejFxs.cxx:44
void SetEndDist(float enddist)
Definition: TrkCntObj.h:31
double YPositionAtTransition(TVector3 vertex, TVector3 dir)
float getScatt(const art::Ptr< rb::Track > track)
Definition: CosRejFxs.cxx:189
void SetCosBakCellND(int cosbakcellnd)
note that fwd<->bak values are interchanged if track direction is backwards, including MC...
Definition: CosRejObj.cxx:273
void SetTrkFwdCell(int trkfwdcell)
Definition: TrkCntObj.h:23
float getAngle(TVector3 trackdir)
Definition: CosRejFxs.cxx:33
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
int evt
int KalBakCell() const
Definition: CosRejObj.cxx:664
int FullNDProjectedCells(TVector3 vertex, TVector3 dir)
Near Detector in the NuMI cavern.
std::string fTrackLabel
TMVA::Reader * fReaderTAP3
void SetCFitSpeed(float cfitspeed)
inverse speed (ns/cm) of track fit by TimingFit module (negative indicates backwards track) for Cosmi...
Definition: CosRejObj.cxx:129
art::ServiceHandle< geo::LiveGeometry > livegeom
int CosBakCell() const
Definition: CosRejObj.cxx:640
int CosFwdCell() const
Definition: CosRejObj.cxx:628
const double j
Definition: BetheBloch.cxx:29
void SetTrackLenKal(float tracklenkal)
Total length of kalman track.
Definition: CosRejObj.cxx:417
void SetKalYPosAtTrans(float kalyposattrans)
Y position at transition to muon catcher, for determining if track went through air gap (ND only) ...
Definition: CosRejObj.cxx:309
RunNumber_t run() const
std::string fTMVApath
Perform a "2 point" Hough transform on a collection of hits.
float CDirScore() const
Definition: CosRejObj.cxx:508
void SetCosFwdSteel(float cosfwdsteel)
distance (dist,cm) through muon catcher traveled, when projected forwards to detector edge...
Definition: CosRejObj.cxx:195
float chisq
Definition: plotTbData.C:11
void SetMinCell(int mincell)
minimum projected number of cell, forwards or backwards, of ANY Kalman track with more than 15 hits ...
Definition: CosRejObj.cxx:321
size_type size() const
Definition: PtrVector.h:302
int MinCell() const
Definition: CosRejObj.cxx:688
TMVA::Reader * fReaderProd5_RHC
std::string fCosmicTrackLabel
void SetCScoreDiff(float cscorediff)
chisq diff between best directed fit and free fit (may indicate if things are wonky) for Cosmic Track...
Definition: CosRejObj.cxx:153
virtual double DistanceFromStart(double z) const
Definition: Track.cxx:229
void SetKDirScore(float kdirscore)
chisq diff between assuming forwards going and backwards going for best ReMId Kalman track ...
Definition: CosRejObj.cxx:135
void SetKalThetaTrue(float kalqtrue)
cosine of angle between true particle (bt most contributing) and reco track dir; use to tell if track...
Definition: CosRejObj.cxx:351
void SetPDGBest(int pdgbest)
pdg code of true most contributing particle to reco track hits. Calculate for best Kalman track...
Definition: CosRejObj.cxx:357
TMVA::Reader * fReaderProd5_Per2
TVector3 NuMIBeamDirection() const
Direction of neutrinos from the NuMI beam (unit vector)
void SetFScattSum(float fscattsum)
sum of scattering angles (Fernanda&#39;s)
Definition: CosRejObj.cxx:393
std::vector< TrackIDE > HitsToTrackIDE(const std::vector< const rb::CellHit * > &hits, bool useBirksE=false) const
Returns vector of TrackIDE structs contributing to the given collection of hits.
Cosmic Rejection PIDs for Numu analysis.
Definition: FillParentInfo.h:9
float TrackLenKal() const
Definition: CosRejObj.cxx:797
float TMVAvarsUCTuned[74]
void SetKalBakAir(float kalbakair)
distance (dist,cm) through air traveled, when projected forwards to detector edge. For ND only, muon catcher air gap. NYI
Definition: CosRejObj.cxx:243
void SetKalFwdAir(float kalfwdair)
distance (dist,cm) through air traveled, when projected forwards to detector edge. For ND only, muon catcher air gap. NYI
Definition: CosRejObj.cxx:225
TVector3 MinXYZ() const
Definition: Cluster.cxx:446
void SetTrkBakCell(int trkbakcell)
Definition: TrkCntObj.h:25
void SetVtxDist(float vtxdist)
minimum distance of start position of best Kalman track to detector wall. Negative indicates outside ...
Definition: CosRejObj.cxx:375
int NTracks3D() const
Definition: CosRejObj.cxx:484
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
void SetKalBakCell(int kalbakcell)
note that fwd<->bak values are interchanged if track direction is backwards
Definition: CosRejObj.cxx:291
void SetCosFwdCellND(int cosfwdcellnd)
cells in that line segment instead of distance of line segment, including MC, ND only
Definition: CosRejObj.cxx:261
std::string fCVNLabel
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
Definition: BackTracker.h:137
float getActivity(art::PtrVector< rb::CellHit > const &trackHits, art::PtrVector< rb::CellHit > const &sliceHits)
Definition: CosRejFxs.cxx:210
float MinDist() const
Definition: CosRejObj.cxx:694
void SetKalFwdCellND(int kalfwdcellnd)
cells in that line segment instead of distance of line segment, including MC, ND only
Definition: CosRejObj.cxx:285
T cos(T number)
Definition: d0nt_math.hpp:78
double ProjectedAirDist(TVector3 vertex, TVector3 dir)
ProductToken< T > consumes(InputTag const &)
Definition: ModuleBase.h:55
double TotalGeV(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple sum of the estimated GeV of all the hits.
Definition: Cluster.cxx:378
void SetFScattMax(float fscattmax)
maximum scattering angle (Fernanda&#39;s)
Definition: CosRejObj.cxx:387
void SetEndDist(float enddist)
minimum distance of end position of best Kalman track to detector wall. Negative indicates outside of...
Definition: CosRejObj.cxx:381
void SetCosHitRatio(float coshitratio)
ratio of hits in cosmic track to hits in slice
Definition: CosRejObj.cxx:411
art::ServiceHandle< cheat::BackTracker > bt
void SetAngleCos(float anglel)
cos of angle of Cosmic track wrt beam
Definition: CosRejObj.cxx:86
void SetTrkBakSteel(float trkbaksteel)
distance (dist,cm) through air traveled, when projected forwards to detector edge. For ND only, muon catcher air gap. NYI
Definition: TrkCntObj.h:37
const art::ProductToken< std::vector< rb::Cluster > > fSlicerToken
virtual void InterpolateXY(double z, double &x, double &y) const
Definition: Track.cxx:325
void SetTrkBakDist(float trkbakdist)
distance (dist,cm) of best track from end to det wall projected forwards (Fwd), based on Kalman Track...
Definition: TrkCntObj.h:33
TMVA::Reader * fReaderTAP1
cosrej::CosRejFxs * fCosRejFxs
TMVA::Reader * fReaderUCTuned
TMVA::Reader * fReaderTAP2
T min(const caf::Proxy< T > &a, T b)
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
T const * get() const
Definition: Ptr.h:149
void SetCosChisq(float coschisq)
same for Cosmic track
Definition: CosRejObj.cxx:177
void SetTrkEPerNHit(float trkepernhit)
Track energy per kalman track hit.
Definition: CosRejObj.cxx:447
ProductID put(std::unique_ptr< PROD > &&edp, FullSemantic< Level::Run > const semantic)
Definition: DataViewImpl.h:730
double bestInvSpeed
ns/cm
void SetCosSlope(float cosslope)
same for Cosmic track
Definition: CosRejObj.cxx:171
virtual double DistanceFromEnd(double z) const
Definition: Track.cxx:255
void SetKalBakDist(float kalbakdist)
cosmic track distance projected backwards from track start
Definition: CosRejObj.cxx:237
void SetCosBakDist(float cosbakdist)
cosmic track distance projected backwards from track start
Definition: CosRejObj.cxx:201
void SetFScattSig(float fscattsig)
sigma of scattering angles (Fernanda&#39;s);
Definition: CosRejObj.cxx:405
void SetAngleKal2nd(float angle2nd)
cos of angle of 2nd best Kalman track wrt beam
Definition: CosRejObj.cxx:435
Encapsulate the geometry of one entire detector (near, far, ndos)
virtual bool Is3D() const
Definition: Prong.h:71
std::string fGeneratorLabel
unsigned int HighestPIDTrack(const std::vector< art::Ptr< rb::Track > > &sliceTracks, const std::string &remidModuleLabel, const art::Event &e)
Definition: ReMId.cxx:249
bool failedToGet() const
Definition: Handle.h:190
art::ServiceHandle< geo::Geometry > geom
std::string fPIDLabel
void SetConCosPID(float concospid)
PID for contained cosmic rejection (The one to use - P5 train)
Definition: CosRejObj.cxx:92
void SetFScattExt(float fscattext)
scattering variable (Fernanda&#39;s)
Definition: CosRejObj.cxx:399
bool fTrackContainmentOnly
void SetTrkFwdDist(float trkfwddist)
Definition: TrkCntObj.h:32
std::vector< const sim::Particle * > HitsToParticle(const std::vector< const rb::CellHit * > &hits) const
Returns vector of sim::Particle objects contributing to the given collection of hits.
void SetERatio(float eratio)
ratio of visible E in best track/ visible E in whole slice
Definition: CosRejObj.cxx:327
void SetKalFwdSteel(float kalfwdsteel)
distance (dist,cm) through muon catcher traveled, when projected forwards to detector edge...
Definition: CosRejObj.cxx:231
void SetCosBakSteel(float cosbaksteel)
distance (dist,cm) through muon catcher traveled, when projected backwards to detector edge...
Definition: CosRejObj.cxx:213
TMVA::Reader * fReaderProd5_Per1
float AngleKal() const
Definition: CosRejObj.cxx:453
float TMVAvarsProd5[12]
enum BeamMode string