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 
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  //----------------------------------------------------------------------
97  fSlicerToken(consumes<std::vector<rb::Cluster>>(pset.get<std::string>("SlicerLabel"))),
98  fGeneratorLabel (pset.get< std::string >("GeneratorLabel")),
99  fNuMIBeamLabel (pset.get< std::string >("NuMIBeamLabel")),
100  fTrackLabel (pset.get< std::string >("TrackLabel")),
101  fCosmicTrackLabel (pset.get< std::string >("CosmicTrackLabel")),
102  fPIDLabel (pset.get< std::string >("PIDLabel")),
103  fCVNLabel (pset.get< std::string >("CVNLabel")),
104  fTMVApath (pset.get< std::string >("tmvapath")),
105  fUseLongestTrack (pset.get< bool >("UseLongestTrack")),
106  fTrackContainmentOnly (pset.get< bool >("TrackContainmentOnly")),
107 
108  fReaderUCTuned (nullptr),
109  fReaderTAP1 (nullptr),
110  fReaderTAP2 (nullptr),
111  fReaderTAP3 (nullptr),
112 
113  fReaderProd5_Per1 (nullptr),
114  fReaderProd5_Per2 (nullptr),
115  fReaderProd5_FHC (nullptr),
116  fReaderProd5_RHC (nullptr)
117  {
118  produces< std::vector<cosrej::CosRejObj> >();
119  produces< art::Assns<cosrej::CosRejObj, rb::Cluster> >();
120  produces< std::vector<cosrej::TrkCntObj> >();
121  produces< art::Assns<cosrej::TrkCntObj, rb::Track> >();
122  }
123 
124 
125  //----------------------------------------------------------------------
127  {
128  if(fCosRejFxs) delete fCosRejFxs;
129  if(fReaderUCTuned) delete fReaderUCTuned;
130  if(fReaderTAP1) delete fReaderTAP1;
131  if(fReaderTAP2) delete fReaderTAP2;
132  if(fReaderTAP3) delete fReaderTAP3;
133 
138  }
139 
140  //-------------------------------------------------------------------
142  {
143  Init();
144  /*
145  fUncorrTA = tfs->make<TH1D>("UncorTA" , "Uncorrected TA; BDT Score, Number of slices" , 100, 0. , 1. );
146  fCorrTA = tfs->make<TH1D>("CorrTA" , "Corrected TA; BDT Score; Number of slices" , 100, 0. , 1. );
147  fDiffTA = tfs->make<TH1D>("DiffTA" , "Diff in Scores; Corr - Uncorr Score; Number of slices", 100, -0.5, 0.5 );
148  fCosRej20 = tfs->make<TH1D>("CosRej20", "CosRej 20; BDT Score; Number of slices" , 100, 0. , 1. );
149 
150  fCosRej20_CorrTA = tfs->make<TH2D>("CosRej20_CorrTA", "CosRej20 vs TA BDT Scores; Cos Rej 20 BDT Score; Corrected TA BDT Score",
151  100, 0., 1, 100, 0., 1. );
152  */
153  }
154 
155  //-------------------------------------------------------------------
157  {
158 
160 
161  // Easier to just declare regardless of detector. Memory footprint is small for BDTs
162  //if(geom->DetId() == novadaq::cnv::kFARDET) {
163 
164  fReaderUCTuned = new TMVA::Reader;
165  fReaderUCTuned->AddVariable("anglekal",&TMVAvarsUCTuned[0]);
166  fReaderUCTuned->AddVariable("kaldiry",&TMVAvarsUCTuned[1]);
167  fReaderUCTuned->AddVariable("tracklen",&TMVAvarsUCTuned[2]);
168  fReaderUCTuned->AddVariable("nhitcos/nhit",&TMVAvarsUCTuned[3]);
169  fReaderUCTuned->AddVariable("cdirscore",&TMVAvarsUCTuned[4]);
170  fReaderUCTuned->AddVariable("max(vty,endy)",&TMVAvarsUCTuned[5]);
171  fReaderUCTuned->AddVariable("cvncosmic",&TMVAvarsUCTuned[6]);
172  fReaderUCTuned->AddVariable("hadEPerNHit",&TMVAvarsUCTuned[7]);
173  fReaderUCTuned->AddVariable("slcCalEPerNHit",&TMVAvarsUCTuned[8]);
174  fReaderUCTuned->AddVariable("trkEPerNHit",&TMVAvarsUCTuned[9]);
175  fReaderUCTuned->AddVariable("vetoangl",&TMVAvarsUCTuned[10]);
176  fReaderUCTuned->AddVariable("scatt",&TMVAvarsUCTuned[11]);
177  fReaderUCTuned->AddSpectator("anglecos",&TMVAvarsUCTuned[12]);
178  fReaderUCTuned->AddSpectator("anglekal",&TMVAvarsUCTuned[13]);
179  fReaderUCTuned->AddSpectator("pngptp",&TMVAvarsUCTuned[14]);
180  fReaderUCTuned->AddSpectator("cosdiry",&TMVAvarsUCTuned[15]);
181  fReaderUCTuned->AddSpectator("scatt",&TMVAvarsUCTuned[16]);
182  fReaderUCTuned->AddSpectator("scatt/tracklen",&TMVAvarsUCTuned[17]);
183  fReaderUCTuned->AddSpectator("slcT",&TMVAvarsUCTuned[18]);
184  fReaderUCTuned->AddSpectator("remid",&TMVAvarsUCTuned[19]);
185  fReaderUCTuned->AddSpectator("oldbdt",&TMVAvarsUCTuned[20]);
186  fReaderUCTuned->AddSpectator("cvnnumu",&TMVAvarsUCTuned[21]);
187  fReaderUCTuned->AddSpectator("osc",&TMVAvarsUCTuned[22]);
188  fReaderUCTuned->AddSpectator("oscVal1",&TMVAvarsUCTuned[23]);
189  fReaderUCTuned->AddSpectator("oscVal2",&TMVAvarsUCTuned[24]);
190  fReaderUCTuned->AddSpectator("oscVal3",&TMVAvarsUCTuned[25]);
191  fReaderUCTuned->AddSpectator("oscVal4",&TMVAvarsUCTuned[26]);
192  fReaderUCTuned->AddSpectator("oscVal5",&TMVAvarsUCTuned[27]);
193  fReaderUCTuned->AddSpectator("nhit",&TMVAvarsUCTuned[28]);
194  fReaderUCTuned->AddSpectator("trueE",&TMVAvarsUCTuned[29]);
195  fReaderUCTuned->AddSpectator("visE",&TMVAvarsUCTuned[30]);
196  fReaderUCTuned->AddSpectator("ntracks3d",&TMVAvarsUCTuned[31]);
197  fReaderUCTuned->AddSpectator("npngs3d",&TMVAvarsUCTuned[32]);
198  fReaderUCTuned->AddSpectator("nprotons",&TMVAvarsUCTuned[33]);
199  fReaderUCTuned->AddSpectator("nneutrons",&TMVAvarsUCTuned[34]);
200  fReaderUCTuned->AddSpectator("nppions",&TMVAvarsUCTuned[35]);
201  fReaderUCTuned->AddSpectator("nnpions",&TMVAvarsUCTuned[36]);
202  fReaderUCTuned->AddSpectator("npions0",&TMVAvarsUCTuned[37]);
203  fReaderUCTuned->AddSpectator("nmuons",&TMVAvarsUCTuned[38]);
204  fReaderUCTuned->AddSpectator("namuons",&TMVAvarsUCTuned[39]);
205  fReaderUCTuned->AddSpectator("numuE",&TMVAvarsUCTuned[40]);
206  fReaderUCTuned->AddSpectator("calccE",&TMVAvarsUCTuned[41]);
207  fReaderUCTuned->AddSpectator("trkccE",&TMVAvarsUCTuned[42]);
208  fReaderUCTuned->AddSpectator("angleE",&TMVAvarsUCTuned[43]);
209  fReaderUCTuned->AddSpectator("recomuonE",&TMVAvarsUCTuned[44]);
210  fReaderUCTuned->AddSpectator("recohadtrkE",&TMVAvarsUCTuned[45]);
211  fReaderUCTuned->AddSpectator("hadcalE",&TMVAvarsUCTuned[46]);
212  fReaderUCTuned->AddSpectator("hadtrkE",&TMVAvarsUCTuned[47]);
213  fReaderUCTuned->AddSpectator("truemuonE",&TMVAvarsUCTuned[48]);
214  fReaderUCTuned->AddSpectator("truegoodmuon",&TMVAvarsUCTuned[49]);
215  fReaderUCTuned->AddSpectator("hadclustncalhits",&TMVAvarsUCTuned[50]);
216  fReaderUCTuned->AddSpectator("hadclustncontplanes",&TMVAvarsUCTuned[51]);
217  fReaderUCTuned->AddSpectator("hadclustfirstplane",&TMVAvarsUCTuned[52]);
218  fReaderUCTuned->AddSpectator("hadclustlastplane",&TMVAvarsUCTuned[53]);
219  fReaderUCTuned->AddSpectator("hadclustncellsedge",&TMVAvarsUCTuned[54]);
220  fReaderUCTuned->AddSpectator("hadclustcalE",&TMVAvarsUCTuned[55]);
221  fReaderUCTuned->AddSpectator("hadclustmeanposX",&TMVAvarsUCTuned[56]);
222  fReaderUCTuned->AddSpectator("hadclustmeanposY",&TMVAvarsUCTuned[57]);
223  fReaderUCTuned->AddSpectator("hadclustmeanposZ",&TMVAvarsUCTuned[58]);
224  fReaderUCTuned->AddSpectator("tracklen2nd",&TMVAvarsUCTuned[59]);
225  fReaderUCTuned->AddSpectator("trkE2nd",&TMVAvarsUCTuned[60]);
226  fReaderUCTuned->AddSpectator("angle2nd",&TMVAvarsUCTuned[61]);
227  fReaderUCTuned->AddSpectator("type",&TMVAvarsUCTuned[62]);
228  fReaderUCTuned->AddSpectator("vtx",&TMVAvarsUCTuned[63]);
229  fReaderUCTuned->AddSpectator("vty",&TMVAvarsUCTuned[64]);
230  fReaderUCTuned->AddSpectator("vtz",&TMVAvarsUCTuned[65]);
231  fReaderUCTuned->AddSpectator("endx",&TMVAvarsUCTuned[66]);
232  fReaderUCTuned->AddSpectator("endy",&TMVAvarsUCTuned[67]);
233  fReaderUCTuned->AddSpectator("endz",&TMVAvarsUCTuned[68]);
234  fReaderUCTuned->AddSpectator("vetoangl",&TMVAvarsUCTuned[69]);
235  fReaderUCTuned->AddSpectator("run",&TMVAvarsUCTuned[70]);
236  fReaderUCTuned->AddSpectator("subrun",&TMVAvarsUCTuned[71]);
237  fReaderUCTuned->AddSpectator("evt",&TMVAvarsUCTuned[72]);
238  fReaderUCTuned->AddSpectator("subevt",&TMVAvarsUCTuned[73]);
239  fReaderUCTuned->BookMVA("UCTuned_BDTA",pidpath+"UCTuned_prod3_BDTA_700.weights.xml");
240 
241  fReaderTAP1 = new TMVA::Reader;
242  fReaderTAP1->AddVariable("anglekal",&TMVAvarsTA[0]);
243  fReaderTAP1->AddVariable("kaldiry",&TMVAvarsTA[1]);
244  fReaderTAP1->AddVariable("tracklen",&TMVAvarsTA[2]);
245  fReaderTAP1->AddVariable("max(vty,endy)",&TMVAvarsTA[3]);
246  fReaderTAP1->AddVariable("cvncosmic",&TMVAvarsTA[4]);
247  fReaderTAP1->AddVariable("min(cosfwdcell+cosbakcell,kalfwdcell+kalbakcell)",&TMVAvarsTA[5]);
248  fReaderTAP1->AddVariable("nhitkal/nhit",&TMVAvarsTA[6]);
249  fReaderTAP1->BookMVA("TAP1_BDTA",pidpath+"TAP1_BDTA.weights.xml");
250 
251  fReaderTAP2 = new TMVA::Reader;
252  fReaderTAP2->AddVariable("anglekal",&TMVAvarsTA[0]);
253  fReaderTAP2->AddVariable("kaldiry",&TMVAvarsTA[1]);
254  fReaderTAP2->AddVariable("tracklen",&TMVAvarsTA[2]);
255  fReaderTAP2->AddVariable("max(vty,endy)",&TMVAvarsTA[3]);
256  fReaderTAP2->AddVariable("cvncosmic",&TMVAvarsTA[4]);
257  fReaderTAP2->AddVariable("min(cosfwdcell+cosbakcell,kalfwdcell+kalbakcell)",&TMVAvarsTA[5]);
258  fReaderTAP2->AddVariable("nhitkal/nhit",&TMVAvarsTA[6]);
259  fReaderTAP2->BookMVA("TAP2_BDTA",pidpath+"TAP2_BDTA.weights.xml");
260 
261  fReaderTAP3 = new TMVA::Reader;
262  fReaderTAP3->AddVariable("anglekal",&TMVAvarsTA[0]);
263  fReaderTAP3->AddVariable("kaldiry",&TMVAvarsTA[1]);
264  fReaderTAP3->AddVariable("tracklen",&TMVAvarsTA[2]);
265  fReaderTAP3->AddVariable("max(vty,endy)",&TMVAvarsTA[3]);
266  fReaderTAP3->AddVariable("cvncosmic",&TMVAvarsTA[4]);
267  fReaderTAP3->AddVariable("min(cosfwdcell+cosbakcell,kalfwdcell+kalbakcell)",&TMVAvarsTA[5]);
268  fReaderTAP3->AddVariable("nhitkal/nhit",&TMVAvarsTA[6]);
269  fReaderTAP3->BookMVA("TAP3_BDTA",pidpath+"TAP3_BDTA.weights.xml");
270 
271  // --- Prod5 FHC Period 1.
272  fReaderProd5_Per1 = new TMVA::Reader;
273  fReaderProd5_Per1 -> AddVariable( "min(CosFwdCell+CosBakCell,KalTrkFwdCell+KalTrkBakCell)", &TMVAvarsProd5[0 ]);
274  fReaderProd5_Per1 -> AddVariable( "KalTrkCosNumi" , &TMVAvarsProd5[1 ]);
275  fReaderProd5_Per1 -> AddVariable( "KalTrkLen" , &TMVAvarsProd5[2 ]);
276  fReaderProd5_Per1 -> AddVariable( "max(KalTrkStY,KalTrkEndY)" , &TMVAvarsProd5[3 ]);
277  fReaderProd5_Per1 -> AddVariable( "cos(KalTrkDirY)" , &TMVAvarsProd5[4 ]);
278  fReaderProd5_Per1 -> AddVariable( "KalTrkNHit/SliceNHit" , &TMVAvarsProd5[5 ]);
279  fReaderProd5_Per1 -> AddVariable( "KalTrkPtOverP" , &TMVAvarsProd5[6 ]);
280  // And my spectators
281  fReaderProd5_Per1 -> AddSpectator( "IsNu" , &TMVAvarsProd5[7 ]);
282  fReaderProd5_Per1 -> AddSpectator( "IsNuMu" , &TMVAvarsProd5[8 ]);
283  fReaderProd5_Per1 -> AddSpectator( "CVNMuon17" , &TMVAvarsProd5[9 ]);
284  fReaderProd5_Per1 -> AddSpectator( "CVNMuon18" , &TMVAvarsProd5[10]);
285  fReaderProd5_Per1 -> AddSpectator( "KalTrkRemID" , &TMVAvarsProd5[11]);
286  fReaderProd5_Per1 -> BookMVA("BDT", pidpath+"TMVA_MProd5_Kal_FHC_Per1_BDT.weights.xml");
287 
288  // --- Prod5 FHC Period 2.
289  fReaderProd5_Per2 = new TMVA::Reader;
290  fReaderProd5_Per2 -> AddVariable( "min(CosFwdCell+CosBakCell,KalTrkFwdCell+KalTrkBakCell)", &TMVAvarsProd5[0 ]);
291  fReaderProd5_Per2 -> AddVariable( "KalTrkCosNumi" , &TMVAvarsProd5[1 ]);
292  fReaderProd5_Per2 -> AddVariable( "KalTrkLen" , &TMVAvarsProd5[2 ]);
293  fReaderProd5_Per2 -> AddVariable( "max(KalTrkStY,KalTrkEndY)" , &TMVAvarsProd5[3 ]);
294  fReaderProd5_Per2 -> AddVariable( "cos(KalTrkDirY)" , &TMVAvarsProd5[4 ]);
295  fReaderProd5_Per2 -> AddVariable( "KalTrkNHit/SliceNHit" , &TMVAvarsProd5[5 ]);
296  fReaderProd5_Per2 -> AddVariable( "KalTrkPtOverP" , &TMVAvarsProd5[6 ]);
297  // And my spectators
298  fReaderProd5_Per2 -> AddSpectator( "IsNu" , &TMVAvarsProd5[7 ]);
299  fReaderProd5_Per2 -> AddSpectator( "IsNuMu" , &TMVAvarsProd5[8 ]);
300  fReaderProd5_Per2 -> AddSpectator( "CVNMuon17" , &TMVAvarsProd5[9 ]);
301  fReaderProd5_Per2 -> AddSpectator( "CVNMuon18" , &TMVAvarsProd5[10]);
302  fReaderProd5_Per2 -> AddSpectator( "KalTrkRemID" , &TMVAvarsProd5[11]);
303  fReaderProd5_Per2 -> BookMVA("BDT", pidpath+"TMVA_MProd5_Kal_FHC_Per2_BDT.weights.xml");
304 
305  // --- Prod5 FHC High Gain.
306  fReaderProd5_FHC = new TMVA::Reader;
307  fReaderProd5_FHC -> AddVariable( "min(CosFwdCell+CosBakCell,KalTrkFwdCell+KalTrkBakCell)", &TMVAvarsProd5[0 ]);
308  fReaderProd5_FHC -> AddVariable( "KalTrkCosNumi" , &TMVAvarsProd5[1 ]);
309  fReaderProd5_FHC -> AddVariable( "KalTrkLen" , &TMVAvarsProd5[2 ]);
310  fReaderProd5_FHC -> AddVariable( "max(KalTrkStY,KalTrkEndY)" , &TMVAvarsProd5[3 ]);
311  fReaderProd5_FHC -> AddVariable( "cos(KalTrkDirY)" , &TMVAvarsProd5[4 ]);
312  fReaderProd5_FHC -> AddVariable( "KalTrkNHit/SliceNHit" , &TMVAvarsProd5[5 ]);
313  fReaderProd5_FHC -> AddVariable( "KalTrkPtOverP" , &TMVAvarsProd5[6 ]);
314  // And my spectators
315  fReaderProd5_FHC -> AddSpectator( "IsNu" , &TMVAvarsProd5[7 ]);
316  fReaderProd5_FHC -> AddSpectator( "IsNuMu" , &TMVAvarsProd5[8 ]);
317  fReaderProd5_FHC -> AddSpectator( "CVNMuon17" , &TMVAvarsProd5[9 ]);
318  fReaderProd5_FHC -> AddSpectator( "CVNMuon18" , &TMVAvarsProd5[10]);
319  fReaderProd5_FHC -> AddSpectator( "KalTrkRemID" , &TMVAvarsProd5[11]);
320  fReaderProd5_FHC -> BookMVA("BDT", pidpath+"TMVA_MProd5_Kal_FHC_High_BDT.weights.xml");
321 
322  // --- Prod5 RHC High Gain.
323  fReaderProd5_RHC = new TMVA::Reader;
324  fReaderProd5_RHC -> AddVariable( "min(CosFwdCell+CosBakCell,KalTrkFwdCell+KalTrkBakCell)", &TMVAvarsProd5[0 ]);
325  fReaderProd5_RHC -> AddVariable( "KalTrkCosNumi" , &TMVAvarsProd5[1 ]);
326  fReaderProd5_RHC -> AddVariable( "KalTrkLen" , &TMVAvarsProd5[2 ]);
327  fReaderProd5_RHC -> AddVariable( "max(KalTrkStY,KalTrkEndY)" , &TMVAvarsProd5[3 ]);
328  fReaderProd5_RHC -> AddVariable( "cos(KalTrkDirY)" , &TMVAvarsProd5[4 ]);
329  fReaderProd5_RHC -> AddVariable( "KalTrkNHit/SliceNHit" , &TMVAvarsProd5[5 ]);
330  fReaderProd5_RHC -> AddVariable( "KalTrkPtOverP" , &TMVAvarsProd5[6 ]);
331  // And my spectators
332  fReaderProd5_RHC -> AddSpectator( "IsNu" , &TMVAvarsProd5[7 ]);
333  fReaderProd5_RHC -> AddSpectator( "IsNuMu" , &TMVAvarsProd5[8 ]);
334  fReaderProd5_RHC -> AddSpectator( "CVNMuon17" , &TMVAvarsProd5[9 ]);
335  fReaderProd5_RHC -> AddSpectator( "CVNMuon18" , &TMVAvarsProd5[10]);
336  fReaderProd5_RHC -> AddSpectator( "KalTrkRemID" , &TMVAvarsProd5[11]);
337  fReaderProd5_RHC -> BookMVA("BDT", pidpath+"TMVA_MProd5_Kal_RHC_High_BDT.weights.xml");
338 
339  }
340 
341  //-------------------------------------------------------------------
343  {
344  //Init();
345 
347  evt.getByToken(fSlicerToken,slicevec);
348 
349  if(slicevec->empty()) {
350  mf::LogWarning ("No Slices")<<"No Slices in the input file";
351  return;
352  }
353 
354  art::FindManyP<rb::Track> trackAssnList(slicevec, evt, fTrackLabel);
355  art::FindManyP<rb::Track> costrackAssnList(slicevec, evt, fCosmicTrackLabel);
356  std::unique_ptr< std::vector<cosrej::CosRejObj> > outputObjects(new std::vector<cosrej::CosRejObj>);
357  std::unique_ptr< art::Assns<cosrej::CosRejObj, rb::Cluster> > assoc(new art::Assns<cosrej::CosRejObj, rb::Cluster>);
358  std::unique_ptr< std::vector<cosrej::TrkCntObj> > outTrkObjects(new std::vector<cosrej::TrkCntObj>);
359  std::unique_ptr< art::Assns<cosrej::TrkCntObj, rb::Track> > trkAssoc(new art::Assns<cosrej::TrkCntObj, rb::Track>);
360 
361  // Karl::Get the horn current.
363  if ( !evt.isRealData() ) {
364  evt.getByLabel(fGeneratorLabel, spillPot);
365  } else {
366  evt.getByLabel(fNuMIBeamLabel , spillPot);
367  }
368  bool isRHC = false;
369  if ( spillPot.failedToGet() ) {
370  mf::LogError("CosRej") << "Spill Data not found, aborting without horn current information";
371  abort();
372  } else {
373  isRHC = spillPot->isRHC;
374  }
375  mf::LogDebug("CosRej") << "Horn Current set to isRHC: "<< isRHC;
376 
377  tf::TimingFitAlg tfalg;
378 
379  for(size_t i = 0; i < slicevec->size(); ++i){
380 
381  art::Ptr<rb::Cluster> slice(slicevec,i);
382 
383  if(slice->IsNoise()) {
384  continue;
385  }
386 
387  cosrej::CosRejObj cosrejobj; // initialization in constructor
388 
389  unsigned int bestIdx = 999; // for index in list of the highest ReMId PID track
390  TVector3 minp, maxp; // minimum and maximum extent vectors for slice
391  int nhadtracks = 0;
392  double fwddist = 0, bakdist = 0, scatt = 0;
393  int fwdcell = 0, bakcell = 0;
394  bool badtrack = false; // check for messed up tracks; a few kalman tracks have nans for directions
395 
396  const std::vector< art::Ptr<rb::Track> > cosTracks = costrackAssnList.at(i);
397 
398  minp = slice->MinXYZ();
399  maxp = slice->MaxXYZ();
400 
401  if(!costrackAssnList.isValid()) {
402  mf::LogError("CosRej") << "CosRej: No Cosmic Tracks! Do not go any further!";
403  outputObjects->push_back(cosrejobj);
404  util::CreateAssn(*this,evt,*(outputObjects.get()),slice,*(assoc.get()));
405  return;
406  }
407 
408  if(cosTracks.size()==1) { // for straight line fit version of Cosmic Track, should only be one track
409  const art::Ptr<rb::Track> track = cosTracks[0];
410  double angleOfBestTrack = fCosRejFxs->getAngle(track->Dir());
411  if(angleOfBestTrack>-6&&angleOfBestTrack<2) cosrejobj.SetAngleCos(angleOfBestTrack);
412 
413  TVector3 start = track->Start();
414  TVector3 end = track->Stop();
415  TVector3 dirbak = -track->Dir().Unit();
416  TVector3 dirfor = (track->TrajectoryPoint(track->NTrajectoryPoints()-1) - track->TrajectoryPoint(track->NTrajectoryPoints()-2)).Unit();
417 
418  cosrejobj.SetCosFwdDist(livegeom->ProjectedLiveDistanceToEdge(end,dirfor));
419  cosrejobj.SetCosBakDist(livegeom->ProjectedLiveDistanceToEdge(start,dirbak));
420  cosrejobj.SetCosFwdCell(livegeom->ProjectedLiveCellsToEdge(end,dirfor));
421  cosrejobj.SetCosBakCell(livegeom->ProjectedLiveCellsToEdge(start,dirbak));
422 
423  if(geom->DetId() == novadaq::cnv::kNEARDET) {
424  TVector3 dirfor = (track->TrajectoryPoint(track->NTrajectoryPoints()-1) - track->TrajectoryPoint(track->NTrajectoryPoints()-2)).Unit();
425 
426  cosrejobj.SetCosFwdSteel(livegeom->ProjectedSteelDist(track->Stop(), dirfor));
427  cosrejobj.SetCosBakSteel(livegeom->ProjectedSteelDist(track->Start(), -(track->Dir())));
428  cosrejobj.SetCosFwdAir(livegeom->ProjectedAirDist(track->Stop(), dirfor));
429  cosrejobj.SetCosBakAir(livegeom->ProjectedAirDist(track->Start(), -(track->Dir())));
430  cosrejobj.SetCosYPosAtTrans(livegeom->YPositionAtTransition(track->Start(), track->Dir()));
431  cosrejobj.SetCosFwdCellND(livegeom->FullNDProjectedCells(track->Stop(), dirfor));
432  cosrejobj.SetCosBakCellND(livegeom->FullNDProjectedCells(track->Start(), -(track->Dir())));
433  }
434 
435  TMVAvarsUCTuned[3] = (1.0*track->NCell() / (1.0*slice->NCell()));
436  double slope = 0, chisq = 0, chisqdiff = 0;
437  if(track->AllCells().size()>2)fCosRejFxs->getFits(track,slope,chisq,chisqdiff);
438  cosrejobj.SetCosSlope(slope);
439  cosrejobj.SetCosChisq(chisq);
440  cosrejobj.SetCosChiDiff(chisqdiff);
441  tf::TimingFitResult tfres = tfalg.HoughFit(track.get());
442  cosrejobj.SetCFitSpeed(tfres.bestInvSpeed);
443  cosrejobj.SetCDirScore(tfres.fwdScore-tfres.revScore);
444  cosrejobj.SetCScoreDiff(tfres.bestScore-std::max(tfres.fwdScore,tfres.revScore)); // Chris's timing fit (Hough)
445  cosrejobj.SetCosHitRatio(1.0*track->NCell()/(1.0*slice->NCell()));
446  TMVAvarsUCTuned[4] = cosrejobj.CDirScore();
447  TMVAvarsUCTuned[10] = fabs(angleOfBestTrack)*(track->Dir().Y() + 1.0);
448 
449  if(bt->HaveTruthInfo()) {
450  art::PtrVector<rb::CellHit> trkhits = track->AllCells();
451  const std::vector<const sim::Particle*> parts = bt->HitsToParticle(trkhits);
452 
453  if(!parts.empty()){
454  const TVector3 trueDir = parts[0]->Momentum().Vect().Unit();
455  cosrejobj.SetCosThetaTrue(trueDir.Dot(track->Dir()));
456  }
457  } // end cosmic truth
458  } // end cosmics
459 
460  if(!trackAssnList.isValid()) {
461  mf::LogWarning("CosRej") << "CosRej: No Kalman Tracks!";
462  outputObjects->push_back(cosrejobj);
463  util::CreateAssn(*this,evt,*(outputObjects.get()),slice,*(assoc.get()));
464  continue;
465  }
466 
467  const std::vector< art::Ptr<rb::Track> > sliceTracks = trackAssnList.at(i);
468 
469  //
470  // loop over tracks and fill per-track variables
471  //
472  for(size_t iTrack = 0; iTrack < sliceTracks.size(); ++iTrack)
473  {
474  cosrej::TrkCntObj trkcntobj;
475  art::Ptr<rb::Track> track = sliceTracks[iTrack];
476 
477  // tracks have to have more than 2 hits
478  if(track->NTrajectoryPoints() < 2){
479  outTrkObjects->push_back(trkcntobj);
480  util::CreateAssn(*this,evt,*(outTrkObjects.get()),track,*(trkAssoc.get()));
481  continue;
482  }
483 
484  // get track direction to check validity
485  double tdx = track->Dir().X();
486  double tdy = track->Dir().Y();
487  double tdz = track->Dir().Z();
488 
489  // ignore bad track
490  if(!(tdx>=-1&&tdx<=1&&tdy>=-1&&tdy<=1&&tdz>=-1&&tdz<=1)){
491  outTrkObjects->push_back(trkcntobj);
492  util::CreateAssn(*this,evt,*(outTrkObjects.get()),track,*(trkAssoc.get()));
493  continue;
494  }
495 
496  TVector3 start = track->Start();
497  TVector3 end = track->Stop();
498  TVector3 dirbak = -track->Dir().Unit();
499  TVector3 dirfor = (track->TrajectoryPoint(track->NTrajectoryPoints()-1) - track->TrajectoryPoint(track->NTrajectoryPoints()-2)).Unit();
500 
501  trkcntobj.SetTrkFwdCell(livegeom->ProjectedLiveCellsToEdge(end, dirfor));
502  trkcntobj.SetTrkBakCell(livegeom->ProjectedLiveCellsToEdge(start, dirbak));
503 
504  // if far detector, fill this way
506  trkcntobj.SetTrkLenInAct(track->TotalLength());
507 
508  // if near detector, fill this way
510  {
511  trkcntobj.SetTrkFwdCellND(livegeom->FullNDProjectedCells(end, dirfor));
512  trkcntobj.SetTrkBakCellND(livegeom->FullNDProjectedCells(start, dirbak));
513 
514  trkcntobj.SetTrkLenInAct(track->DistanceFromStart(livegeom->MCFrontZ()));
515  trkcntobj.SetTrkLenInCat(track->DistanceFromEnd(livegeom->MCFrontZ()));
516 
517  // get y opsition at transition plane
518  double xtran, ytran;
519  track->InterpolateXY(livegeom->MCFrontZ(), xtran, ytran);
520  trkcntobj.SetTrkYPosAtTrans(ytran);
521 
522  // Other containment variables
523  trkcntobj.SetTrkFwdSteel(livegeom->ProjectedSteelDist(end, dirfor));
524  trkcntobj.SetTrkBakSteel(livegeom->ProjectedSteelDist(start, dirbak));
525  trkcntobj.SetTrkFwdAir(livegeom->ProjectedAirDist(end, dirfor));
526  trkcntobj.SetTrkBakAir(livegeom->ProjectedAirDist(start, dirbak));
527  }
528 
529  // Other containment variables
530  trkcntobj.SetVtxDist(std::min(livegeom->DistToEdgeXY(start),livegeom->DistToEdgeZ(start)));
532  trkcntobj.SetTrkFwdDist(livegeom->ProjectedLiveDistanceToEdge(end,dirfor));
533  trkcntobj.SetTrkBakDist(livegeom->ProjectedLiveDistanceToEdge(start,dirbak));
534 
535  outTrkObjects->push_back(trkcntobj);
536  util::CreateAssn(*this,evt,*(outTrkObjects.get()),track,*(trkAssoc.get()));
537  } // end of loop over tracks on slice
538 
540  outputObjects->push_back(cosrejobj);
541  util::CreateAssn(*this,evt,*(outputObjects.get()),slice,*(assoc.get()));
542  continue;
543  }
544 
545  art::FindOneP<remid::ReMId> remidVec(sliceTracks, evt, fPIDLabel);
546 
547  if(fUseLongestTrack) {
548  bestIdx = 0;
549  if(sliceTracks.size() == 0) {
550  outputObjects->push_back(cosrejobj);
551  continue;
552  }
553  const art::Ptr<rb::Track> track = sliceTracks[0];
554  if(!track->Is3D()) {
555  outputObjects->push_back(cosrejobj);
556  continue;
557  }
558  cosrejobj.SetNTracks3D(sliceTracks.size());
559  }
560  else if(!remidVec.isValid()) {
561  mf::LogWarning("CosRej") << "CosRej: No ReMId information!";
562  outputObjects->push_back(cosrejobj);
563  util::CreateAssn(*this,evt,*(outputObjects.get()),slice,*(assoc.get()));
564  continue;
565  }
566  else {
567  bestIdx = remid::HighestPIDTrack(sliceTracks, fPIDLabel, evt);
568  for(size_t j = 0; j < sliceTracks.size(); ++j){ // loop over Kalman tracks
569  const art::Ptr<rb::Track> track = sliceTracks[j];
570  art::Ptr<remid::ReMId> trackRemid = remidVec.at(j);
571  if(!track->Is3D() || trackRemid->Value() < 0) continue;
572  if(cosrejobj.NTracks3D()==-5) cosrejobj.SetNTracks3D(1);
573  else cosrejobj.SetNTracks3D(cosrejobj.NTracks3D()+1);
574  if(trackRemid->Value() < 0.5) nhadtracks++;
575  }
576  }
577 
578  // Main condition
579  if(cosrejobj.NTracks3D()>0 && bestIdx!=999) {
580 
581  const art::Ptr<rb::Track> track = sliceTracks[bestIdx];
582 
583  art::PtrVector<rb::CellHit> trkhits = track->AllCells();
584 
585  if(cosrejobj.NTracks3D()>1) {
586  // Getting the 2nd highest ReMId PID track
587  unsigned int bestIdx2nd = 999;
588  double highestPID2nd = -999.0;
589  double bestTrkLen2nd = 0.0;
590  for(size_t j = 0; j < sliceTracks.size(); ++j){ // loop over Kalman tracks
591  if(j == bestIdx) continue;
592  art::Ptr<remid::ReMId> trackRemid2nd = remidVec.at(j);
593  if(trackRemid2nd->Value() < highestPID2nd) continue;
594  if(trackRemid2nd->Value() == highestPID2nd){
595  if(sliceTracks[j]->TotalLength() < bestTrkLen2nd) continue;
596  }
597  highestPID2nd = trackRemid2nd->Value();
598  bestTrkLen2nd = sliceTracks[j]->TotalLength();
599  bestIdx2nd = j;
600  }//end loop over kalman trks for 2nd best ReMId
601  const art::Ptr<rb::Track> track2nd = sliceTracks[bestIdx2nd];
602  cosrejobj.SetTrackLenKal2nd(track2nd->TotalLength());
603  cosrejobj.SetTrackE2nd(track2nd->CalorimetricEnergy());
604  double angleOf2ndBestTrack = fCosRejFxs->getAngle(track2nd->Dir());
605  if((angleOf2ndBestTrack>-6) && (angleOf2ndBestTrack<2)) cosrejobj.SetAngleKal2nd(angleOf2ndBestTrack);
606  TVector3 firstDir = track->Dir().Unit();
607  TVector3 secondDir = track2nd->Dir().Unit();
608  double angleBtwKalTrks = (firstDir*secondDir);
609  cosrejobj.SetAngleBtwKalTrks(angleBtwKalTrks);
610  }
611  else{
612  cosrejobj.SetTrackLenKal2nd(-5.0);
613  cosrejobj.SetTrackE2nd(-5.0);
614  cosrejobj.SetAngleKal2nd(-5.0);
615  cosrejobj.SetAngleBtwKalTrks(-5.0);
616  }
617 
618  if(bt->HaveTruthInfo()) {
619  const std::vector<const sim::Particle*> parts = bt->HitsToParticle(trkhits);
620  const std::vector< cheat::TrackIDE > trids = bt->HitsToTrackIDE(trkhits);
621 
622  if(!parts.empty()){
623  const TVector3 trueDir = parts[0]->Momentum().Vect().Unit();
624  cosrejobj.SetKalThetaTrue(trueDir.Dot(track->Dir()));
625  cosrejobj.SetPDGBest(parts[0]->PdgCode());
626  }
627  }
628 
629  double tdx = track->Dir().X();
630  double tdy = track->Dir().Y();
631  double tdz = track->Dir().Z();
632 
633  if(!(tdx>=-1&&tdx<=1&&tdy>=-1&&tdy<=1&&tdz>=-1&&tdz<=1)) badtrack = true;
634 
635  if(badtrack) continue;
636 
637  TVector3 start = track->Start();
638  TVector3 end = track->Stop();
639  TVector3 dirbak = -track->Dir().Unit();
640  TVector3 dirfor = (track->TrajectoryPoint(track->NTrajectoryPoints()-1) - track->TrajectoryPoint(track->NTrajectoryPoints()-2)).Unit();
641 
642  cosrejobj.SetKalFwdDist(livegeom->ProjectedLiveDistanceToEdge(end,dirfor));
643  cosrejobj.SetKalBakDist(livegeom->ProjectedLiveDistanceToEdge(start,dirbak));
644  cosrejobj.SetKalFwdCell(livegeom->ProjectedLiveCellsToEdge(end,dirfor));
645  cosrejobj.SetKalBakCell(livegeom->ProjectedLiveCellsToEdge(start,dirbak));
646 
647  cosrejobj.SetVtxDist(std::min(livegeom->DistToEdgeXY(start),livegeom->DistToEdgeZ(start)));
649 
650  if(geom->DetId() == novadaq::cnv::kNEARDET) {
651  cosrejobj.SetKalFwdSteel(livegeom->ProjectedSteelDist(track->Stop(), dirfor));
652  cosrejobj.SetKalBakSteel(livegeom->ProjectedSteelDist(track->Start(), -(track->Dir())));
653  cosrejobj.SetKalFwdAir(livegeom->ProjectedAirDist(track->Stop(), dirfor));
654  cosrejobj.SetKalBakAir(livegeom->ProjectedAirDist(track->Start(), -(track->Dir())));
655 
656  double xtran, ytran;
657  track->InterpolateXY(livegeom->MCFrontZ(), xtran, ytran);
658  cosrejobj.SetKalYPosAtTrans(ytran);
659  cosrejobj.SetKalFwdCellND(livegeom->FullNDProjectedCells(track->Stop(), dirfor));
660  cosrejobj.SetKalBakCellND(livegeom->FullNDProjectedCells(track->Start(), -(track->Dir())));
661  }
662 
663  if(std::min(fwddist,bakdist) < cosrejobj.MinDist()) cosrejobj.SetMinDist(std::min(fwddist,bakdist));
664  if(std::min(fwdcell,bakcell) < cosrejobj.MinCell()) cosrejobj.SetMinCell(std::min(fwdcell,bakcell));
665 
666  cosrejobj.SetERatio(track->TotalGeV()/slice->TotalGeV());
667  scatt = fCosRejFxs->getScatt(track);
668  if(track->TotalLength()>0) cosrejobj.SetScatt(scatt/track->TotalLength());
669 
670  double ferpars[4];
671  bool fdef = 0;
672  bool fFlag = fCosRejFxs->FScatterEstim(track,fdef,ferpars[0],ferpars[1],ferpars[2],ferpars[3]); // Fernandas scattering energy estimate
673  if(fFlag==1) {
674  cosrejobj.SetFScattExt(ferpars[0]);
675  cosrejobj.SetFScattSig(ferpars[1]);
676  cosrejobj.SetFScattMax(ferpars[2]);
677  cosrejobj.SetFScattSum(ferpars[3]);
678  }
679  double slope = 0, chisq = 0, chisqdiff = 0;
680  if(track->AllCells().size()>2)fCosRejFxs->getFits(track,slope,chisq,chisqdiff);
681  cosrejobj.SetKalSlope(slope);
682  cosrejobj.SetKalChisq(chisq);
683  cosrejobj.SetKalChiDiff(chisqdiff);
684 
685  tf::TimingFitResult tfres = tfalg.HoughFit(track.get());
686  cosrejobj.SetKFitSpeed(tfres.bestInvSpeed);
687  cosrejobj.SetKDirScore(tfres.fwdScore-tfres.revScore);
688  cosrejobj.SetKScoreDiff(tfres.bestScore-std::max(tfres.fwdScore,tfres.revScore));
689  cosrejobj.SetStartAct(fCosRejFxs->getActivity(track->AllCells(),slice->AllCells()));
690  cosrejobj.SetNHadTracks(nhadtracks);
691 
692  double angleOfBestTrack = fCosRejFxs->getAngle(sliceTracks[bestIdx]->Dir());
693  if(angleOfBestTrack>-6&&angleOfBestTrack<2) cosrejobj.SetAngleKal(angleOfBestTrack);
694 
695  cosrejobj.SetTrackLenKal(track->TotalLength());
696 
697  if(geom->DetId() == novadaq::cnv::kFARDET) {
698  const art::Ptr<rb::Track> track = sliceTracks[bestIdx];
699 
700  float fslcCalEPerNHit = 0.0;
701  float fhadEPerNHit = 0.0;
702  int frun = evt.run();
703 
704  art::FindManyP<cvn::Result> fmCVN(slicevec, evt, fCVNLabel); // copied from CAFMAker
705  if(fmCVN.isValid()) {
706  const std::vector<art::Ptr<cvn::Result>> results = fmCVN.at(i);
707  if(!results.empty()) {
708 
709  // --- Uncontained BDT
710  if(slice->NCell() > 0) fslcCalEPerNHit = ((slice->CalorimetricEnergy()) / (1.78*slice->NCell()));
711  if(track->NCell() > 0) cosrejobj.SetTrkEPerNHit((track->CalorimetricEnergy())/(track->NCell()));
712  if((slice->NCell()>0) && (track->NCell()>0)) fhadEPerNHit = (slice->CalorimetricEnergy() - track->CalorimetricEnergy())/(slice->NCell() - track->NCell());
713 
714  TMVAvarsUCTuned[0] = cosrejobj.AngleKal();
715  TMVAvarsUCTuned[1] = track->Dir().Y();
716  TMVAvarsUCTuned[2] = cosrejobj.TrackLenKal(); //track->TotalLength();
717  // [3] and [4] above.
718  TMVAvarsUCTuned[5] = std::max(track->Start().Y(),track->Stop().Y());
719  TMVAvarsUCTuned[6] = results[0]->fOutput[3]; // 2020 CVN has Cosmic at output 3.
720  TMVAvarsUCTuned[7] = fhadEPerNHit;
721  TMVAvarsUCTuned[8] = fslcCalEPerNHit;
722  TMVAvarsUCTuned[9] = cosrejobj.TrkEPerNHit(); //ftrkEPerNHit;
723  // [10] above.
724  TMVAvarsUCTuned[11] = cosrejobj.Scatt(); //track->TotalLength();
725  // --- Set the Uncontained BDT
726  cosrejobj.SetUnconTunedCosPID(fReaderUCTuned->EvaluateMVA("UCTuned_BDTA"));
727 
728  // --- Prod3 Trained BDT
729  TMVAvarsTA[0] = cosrejobj.AngleKal();
730  TMVAvarsTA[1] = track->Dir().Y();;
731  TMVAvarsTA[2] = track->TotalLength();
732  TMVAvarsTA[3] = std::max(track->Start().Y(),track->Stop().Y());
733  TMVAvarsTA[4] = results[0]->fOutput[3]; // 2020 CVN has Cosmic at output 3.
734  TMVAvarsTA[5] = std::min(cosrejobj.CosFwdCell()+cosrejobj.CosBakCell(),cosrejobj.KalFwdCell()+cosrejobj.KalBakCell());
735  TMVAvarsTA[6] = track->NCell() / (float)slice->NCell();
736  // --- Set the Prod3 Trained BDT
737  if (frun > 10000 && frun < 17300) {
738  cosrejobj.SetOldCosPID(fReaderTAP1->EvaluateMVA("TAP1_BDTA"));
739  } else if(frun > 17300 && frun < 20753) {
740  cosrejobj.SetOldCosPID(fReaderTAP2->EvaluateMVA("TAP2_BDTA"));
741  } else if (frun > 20752) {
742  cosrejobj.SetOldCosPID(fReaderTAP3->EvaluateMVA("TAP3_BDTA"));
743  }
744  } // end results not empty
745  } else {
746  mf::LogError("CosRej") << "CVN Data Product " << fCVNLabel << " not found, needed for CosRej2019, aborting.";
747  abort();
748  } // end CVN handle is valid
749 
750  // Karl::This is where we want to evaluate the Prod5 cosmic rejection...
751  // First, fill my BDTVars for this event...
752  TMVAvarsProd5[0] = std::min(cosrejobj.CosFwdCell()+cosrejobj.CosBakCell(),cosrejobj.KalFwdCell()+cosrejobj.KalBakCell());
753  TMVAvarsProd5[1] = track->Dir().Dot( geom->NuMIBeamDirection() ); // kCosNuMi
754  TMVAvarsProd5[2] = track->TotalLength()/100; // kTrkLength
755  TMVAvarsProd5[3] = std::max( track->Start().Y(), track->Stop().Y() )/100; // kMaxKalTrStEnY
756  TMVAvarsProd5[4] = cos( track->Dir().Y() ); // kcosDirY
757  TMVAvarsProd5[5] = track->NCell() / (float)slice->NCell(); // kKalHitRat
758  TMVAvarsProd5[6] = sqrt( 1 - (TMVAvarsProd5[1]*TMVAvarsProd5[1]) ); // kNumuMuonPtP
759  // Then, store the official number.
760  if (isRHC) { // RHC, so use the RHC BDT :O
761  cosrejobj.SetConCosPID ( fReaderProd5_RHC ->EvaluateMVA("BDT") );
762  } else { // FHC, so need to figure out if P1, P2, HighGain
763  if ( frun <= 17139 ) {
764  cosrejobj.SetConCosPID( fReaderProd5_Per1->EvaluateMVA("BDT") );
765  } else if ( frun <= 19746 ) {
766  cosrejobj.SetConCosPID( fReaderProd5_Per2->EvaluateMVA("BDT") );
767  } else {
768  cosrejobj.SetConCosPID( fReaderProd5_FHC ->EvaluateMVA("BDT" ) );
769  }
770  }
771  /*
772  std::cout << "\nLooking at Evt " << evt.event() << ", Slice " << i << "."
773  << I got: " << cosrejobj.ConCosPID() << ", Comp to " << cosrejobj.OldCosPID()
774  << "\n\t TMVAvarsProd5[0]: MinCellSum = " << TMVAvarsProd5[0]
775  << "\n\t TMVAvarsProd5[1]: kCosNuMi = " << TMVAvarsProd5[1]
776  << "\n\t TMVAvarsProd5[2]: kTrkLength = " << TMVAvarsProd5[2]
777  << "\n\t TMVAvarsProd5[3]: kMaxKalTrStEnY = " << TMVAvarsProd5[3]
778  << "\n\t TMVAvarsProd5[4]: kcosDirY = " << TMVAvarsProd5[4]
779  << "\n\t TMVAvarsProd5[5]: kKalHitRat = " << TMVAvarsProd5[5] << " .... " << track->NCell() / slice->NCell()
780  << " .... TMVAvarsUCTuned[3] = " << TMVAvarsUCTuned[3] << " .... " << TMVAvarsTA[6]
781  << "\n\t TMVAvarsProd5[6]: kNumuMuonPtP = " << TMVAvarsProd5[6]
782  << std::endl;
783  //*/
784  /*
785  fCosRej20 -> Fill( cosrejobj.ConCosPID() );
786  fCosRej20_CorrTA -> Fill( cosrejobj.ConCosPID(), cosrejobj.OldCosPID() );
787  //*/
788  } // end if FarDet
789  } // end if good track
790 
791  outputObjects->push_back(cosrejobj); // write out object with all information
792  util::CreateAssn(*this,evt,*(outputObjects.get()),slice,*(assoc.get()));
793 
794  } //end loop over slices
795 
796  evt.put(std::move(outputObjects));
797  evt.put(std::move(assoc));
798  evt.put(std::move(outTrkObjects));
799  evt.put(std::move(trkAssoc));
800 
801  } //end produce
802 
803 } //end namespace
804 
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
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
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
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
bool isRealData() const
Definition: Event.h:83
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:130
void produce(art::Event &evt)
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
float TrkEPerNHit() const
Definition: CosRejObj.cxx:827
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
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
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
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:308
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
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
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
T const * get() const
Definition: Ptr.h:321
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
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
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:133
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)
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
bool getByToken(ProductToken< PROD > const &token, Handle< PROD > &result) const
Definition: DataViewImpl.h:387
RunNumber_t run() const
Definition: Event.h:77
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
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
ProductToken< T > consumes(InputTag const &)
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:196
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]