CVNCosmicEventDump_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file CVNCosmicEventDump_module.cc
3 // \brief Analyzer module for creating CVN Cosmic Rejection PixelMap objects
4 // \original author Dominick Rocco - rocco@physics.umn.edu
5 // \modified author Adam Moren - amoren@fnal.gov
6 ////////////////////////////////////////////////////////////////////////
7 
8 // C/C++ includes
9 #include <iostream>
10 #include <sstream>
11 #include <vector>
12 
13 // ROOT includes
14 #include "TTree.h"
15 #include "TH2F.h"
16 
17 // Framework includes
24 #include "fhiclcpp/ParameterSet.h"
28 
29 
30 // NOvASoft includes
31 #include "MCCheater/BackTracker.h"
32 #include "RecoBase/Cluster.h"
33 #include "RecoBase/FilterList.h"
37 #include "Utilities/AssociationUtil.h"
38 
39 #include "CVN/func/AssignLabels.h"
40 #include "CVN/func/TrainingData.h"
45 
47 #include "CVN/func/ParticlesType.h"
48 
49 // Required for CAFCutter
52 
53 
54 namespace cvn {
56  public:
57  explicit CVNCosmicEventDump(fhicl::ParameterSet const& pset);
59 
60  void analyze(const art::Event& evt);
61  void reconfigure(const fhicl::ParameterSet& pset);
62  void beginJob();
63  void endJob();
64 
65 
66 
67  private:
68 
73 
74  bool interactionTraining = true;
75  bool particlesTraining = false;
76 
79  TTree* fTrainTree;
80 
82 
83  /// Check rb::IsFiltered?
85  std::vector<std::string> fPreselectionLabels;
86 
87  /// Function to extract TH2 from PixelMap and write to TFile
88  void WriteMapTH2(const art::Event& evt, int slice, const PixelMap& pm);
89 
90  std::string fCAFLabel; ///< label for the process that made the standard records
91  bool fApplyCAFCuts; ///< should we apply the CAF-level cuts?
92  int fCutType; ///< what cuts to apply (see CAFCutter.h)
93 
94  recovalid::CAFCutter fCAFCutter; ///< the CAFCutter helper class
95  };
96 
97  //.......................................................................
99  : EDAnalyzer(pset)
100  {
101  this->reconfigure(pset);
102  }
103 
104  //......................................................................
106  {
107  }
108 
109  //......................................................................
111  {
112  fClusterLabel = pset.get<std::string> ("ClusterLabel");
113  fPixelMapInput = pset.get<std::string> ("PixelMapInput");
114  fChooseLabels = pset.get<std::string> ("ChooseLabels");
115  fWriteMapTH2 = pset.get<bool> ("WriteMapTH2");
116  fObeyPreselection = pset.get<bool> ("ObeyPreselection");
117  fPreselectionLabels = pset.get< std::vector<std::string> >("PreselectionLabels");
118  fCAFLabel = pset.get< std::string >("CAFLabel");
119  fApplyCAFCuts = pset.get< bool > ("ApplyCAFCuts");
120  fCutType = pset.get< int > ("CutType");
121  }
122 
123  //......................................................................
125  {
126 
128 
129  if(fChooseLabels == "particleType"){
130  particlesTraining = true;
131  interactionTraining = false;
132  }
133 
134  fTrainTree = tfs->make<TTree>("CVNTrainTree", "Training records");
136  fTrainTree->Branch("train", "cvn::TrainingData", &fTrain);
138  fTrainTree->Branch("train", "cvn::ParticlesTrainingData", &fParticlesTrain);
139 
140  }
141 
142  //......................................................................
144  {
145  }
146 
147  //......................................................................
149  {
150 
151  // get the slices
153  evt.getByLabel(art::InputTag(fClusterLabel, "cvnmap"), clusterHandle);
154  const std::vector<rb::Cluster> & clusters = *clusterHandle;
156  for(unsigned int i = 0; i < clusterHandle->size(); ++i){
157  slices.push_back(art::Ptr<rb::Cluster>(clusterHandle, i));
158  }
159 
160 
161 // art::FindManyP<PixelMap> fmPixelMap(clusterHandle, evt, art::InputTag("cvnmap", fPixelMapInput));
162  art::FindManyP<PixelMap> fmPixelMap(clusterHandle, evt, "cvnmap");
163  // Get the Standard Records
164  art::FindManyP<caf::StandardRecord> recordFMP(slices, evt, fCAFLabel);
165 
166  std::vector<int> bestNuId;
167  std::vector<cheat::NeutrinoWithIndex> nus;
168  std::vector< std::vector<cheat::NeutrinoEffPur> > sEffPur;
169  if(fBT->HaveTruthInfo()){
170  nus = fBT->allMCTruth();
171  sEffPur = fBT->SlicesToMCTruthsTable(slices);
172  bestNuId = fBT->SliceToOrderedNuIdsByEnergy(nus, sEffPur);
173  }
174 
175  for(size_t iClust = 0; iClust < clusters.size(); ++iClust){
176  if(fObeyPreselection && rb::IsFiltered(evt, clusterHandle, iClust,
177  fPreselectionLabels)) continue;
178  // get record associated with the slice
179  std::vector< art::Ptr<caf::StandardRecord> > records;
180  bool pass = true;
181  // Apply the CAF-level cuts
182  if(fApplyCAFCuts && recordFMP.isValid()){
183  records = recordFMP.at(iClust);
184  if(records.size() > 0) pass = fCAFCutter.passCuts(fCutType, records[0].get());
185  }
186 
187 // if(sEffPur[iClust][bestNuId[iClust]].nSliceHits <= 70) continue;
188  if(!pass) continue;
189  const rb::Cluster& cluster = clusters[iClust];
190  if(!fmPixelMap.isValid()) continue;
191 
192  const std::vector< art::Ptr<PixelMap> > pixelMaps = fmPixelMap.at(iClust);
193  if(pixelMaps.empty()) continue;
194  InteractionType interaction = kCosmic;
195  FinalStateType finalstate = kOtherFS;
196  FinalStateProngsType finalstateprongs = kOtherFSp;
197  ParentParticleType parentparticle = kUnknownPP;
198 
199  ParticlesType particles_type = kCosmic_PT;
200 
201  float nuEnergy = 0;
202  float lepEnergy = 0;
203  double vtxx = std::numeric_limits<double>::lowest();
204  double vtxy = std::numeric_limits<double>::lowest();
205  double vtxz = std::numeric_limits<double>::lowest();
207  if(fBT->HaveTruthInfo()){
208  if(bestNuId[iClust] != -1) {
209 
210  truth = sEffPur[iClust][bestNuId[iClust]].neutrinoInt;
211 
212  interaction = SliceClassify(truth, cluster, nuEnergy, lepEnergy, vtxx, vtxy, vtxz);
213  particles_type = ParticlesSliceClassify(truth, cluster, nuEnergy, lepEnergy, vtxx, vtxy, vtxz);
214  finalstate = GetFinalStateType(truth);
215  finalstateprongs = GetFinalStateProngsType(truth);
216  parentparticle = GetParentParticleType(truth, evt);
217 
218  if (interaction != kCosmic) {
219  printf("\ncluster: %li Particles_type: %i Interaction: %i, nSliceHits %i",
220  iClust, particles_type, interaction, sEffPur[iClust][bestNuId[iClust]].nSliceHits);
221  if(sEffPur[iClust][bestNuId[iClust]].nSliceHits <= 10) {printf(" **** Skipping ***\n\n");}
222  else printf("\n\n");
223 // if(sEffPur[iClust][bestNuId[iClust]].nSliceHits >= 70) continue;
224 // printf("\n it continued with hits = %i \n", sEffPur[iClust][bestNuId[iClust]].nSliceHits);
225  }
226 
227  //If the cluster contains the Neutrino event, but has less than 10 hits from it, cut it out
228  if(sEffPur[iClust][bestNuId[iClust]].nSliceHits <= 10 && interaction != kCosmic) continue;
229  }
230  }
231  else particles_type = kCosmic_PT; // ! Careful This will only work for cosmic data !
232 
233 //if (interaction != 14) {
234 // printf("\ncluster: %li Particles_type: %i Interaction: %i \n", iClust, particles_type, interaction);
235 // if(sEffPur[iClust][bestNuId[iClust]].nSliceHits >= 70) continue;
236 // printf("\n** Saving **");
237 // }
238 
239 
240 
242  TrainingData train(interaction, parentparticle, finalstate, finalstateprongs,
243  particles_type, nuEnergy, lepEnergy,
244  vtxx, vtxy, vtxz, *pixelMaps[0]);
245  if(fWriteMapTH2) WriteMapTH2(evt, (int)iClust, train.fPMap);
246  fTrain = &train;
247  fTrainTree->Fill();
248  }
249 
250  if(particlesTraining){
251  ParticlesTrainingData train(particles_type, nuEnergy, lepEnergy,
252  vtxx, vtxy, vtxz, *pixelMaps[0]);
253  if(fWriteMapTH2) WriteMapTH2(evt, (int)iClust, train.fPMap);
255  fTrainTree->Fill();
256  }
257  } // end loop over slices
258 
259  } // analyze
260 
261  //----------------------------------------------------------------------
262  void CVNCosmicEventDump::WriteMapTH2(const art::Event& evt, int slice, const PixelMap& pm)
263  {
264  TString mapid = "_r"+std::to_string(evt.run())+"_s"+std::to_string(evt.subRun())+"_e"
265  +std::to_string(evt.event())+"_sl"+std::to_string(slice);
266  TString data = "PixelMap";
267  TString label = "TruthPixelMap";
268  TString object= "ObjectPixelMap";
269  TString time = "TimePixelMap";
270  TH2F* hist = pm.ToTH2( data+mapid );
271  TH2F* histL = pm.ToLabTH2( label+mapid);
272  TH2F* histO = pm.ToObjTH2( object+mapid);
273  TH2F* histX = pm.SingleViewToTH2(0, data+"ViewX"+mapid);
274  TH2F* histY = pm.SingleViewToTH2(1, data+"ViewY"+mapid);
275  TH2F* histXL = pm.SingleViewToLabTH2(0, label+"ViewX"+mapid, false);
276  TH2F* histYL = pm.SingleViewToLabTH2(1, label+"ViewY"+mapid, false);
277  TH2F* histXT = pm.SingleViewToLabTH2(0, time+"ViewX"+mapid, true);
278  TH2F* histYT = pm.SingleViewToLabTH2(1, time+"ViewY"+mapid, true);
279  TH2F* histXO = pm.SingleViewToObjTH2(0, object+"ViewX"+mapid, false);
280  TH2F* histYO = pm.SingleViewToObjTH2(1, object+"ViewY"+mapid, false);
281 
283  TH2F* histWrite = tfs->make<TH2F>(*hist);
284  histWrite->Write();
285  TH2F* histWriteL = tfs->make<TH2F>(*histL);
286 // histWriteL->GetZaxis()->SetRangeUser(0,10);
287  histWriteL->Write();
288  TH2F* histWriteO = tfs->make<TH2F>(*histO);
289  histWriteL->GetZaxis()->SetRangeUser(0,10);
290  histWriteO->Write();
291  TH2F* histWriteX = tfs->make<TH2F>(*histX);
292  histWriteX->Write();
293  TH2F* histWriteY = tfs->make<TH2F>(*histY);
294  histWriteY->Write();
295  TH2F* histWriteXL = tfs->make<TH2F>(*histXL);
296  histWriteXL->Write();
297  TH2F* histWriteYL = tfs->make<TH2F>(*histYL);
298  histWriteYL->Write();
299  TH2F* histWriteXO = tfs->make<TH2F>(*histXO);
300  histWriteXO->Write();
301  TH2F* histWriteYO = tfs->make<TH2F>(*histYO);
302  histWriteYO->Write();
303  TH2F* histWriteXT = tfs->make<TH2F>(*histXT);
304  histWriteXT->Write();
305  TH2F* histWriteYT = tfs->make<TH2F>(*histYT);
306  histWriteYT->Write();
307 
308  delete hist;
309  delete histWrite;
310  delete histL;
311  delete histWriteL;
312  delete histO;
313  delete histXO;
314  delete histYO;
315  delete histWriteO;
316  delete histWriteXO;
317  delete histWriteYO;
318  delete histX;
319  delete histWriteX;
320  delete histY;
321  delete histWriteY;
322  delete histXL;
323  delete histWriteXL;
324  delete histWriteXT;
325  delete histYL;
326  delete histWriteYL;
327  delete histWriteYT;
328 
329  } // WriteMapTH2
330 
332 } // end namespace cvn
333 ////////////////////////////////////////////////////////////////////////
art::ServiceHandle< cheat::BackTracker > fBT
SubRunNumber_t subRun() const
Definition: Event.h:72
back track the reconstruction to the simulation
InteractionType SliceClassify(const art::Ptr< simb::MCTruth > truth, const rb::Cluster &slice, float &nuEnergy, float &lepEnergy, double &vtxx, double &vtxy, double &vtxz)
The TrainingData objects contains a PixelMap and the output class type, and any other bit that goes i...
FinalStateType GetFinalStateType(const art::Ptr< simb::MCTruth > truth)
void analyze(const art::Event &evt)
TH2F * SingleViewToTH2(const unsigned int &view, TString name) const
Definition: PixelMap.cxx:312
std::vector< NeutrinoWithIndex > allMCTruth() const
PixelMap fPMap
PixelMap for the event.
float vtxx
A collection of associated CellHits.
Definition: Cluster.h:47
enum cvn::Interaction InteractionType
std::vector< std::string > fPreselectionLabels
ParticlesType ParticlesSliceClassify(const art::Ptr< simb::MCTruth > truth, const rb::Cluster &slice, float &nuEnergy, float &lepEnergy, double &vtxx, double &vtxy, double &vtxz)
DEFINE_ART_MODULE(TestTMapFile)
Defines an enumeration for prong classification.
Particle class.
std::vector< std::vector< cheat::NeutrinoEffPur > > SlicesToMCTruthsTable(const std::vector< const rb::Cluster * > &sliceList) const
Given ALL the slices in an event, including the noise slice, returns a vector of vector of structures...
void WriteMapTH2(const art::Event &evt, int slice, const PixelMap &pm)
Function to extract TH2 from PixelMap and write to TFile.
Cosmic ray background.
Something else. Tau? Hopefully we don&#39;t use this.
const char * label
PixelMap fPMap
PixelMap for the event.
Definition: TrainingData.h:69
TH2F * ToLabTH2(TString name) const
Definition: PixelMap.cxx:258
const XML_Char const XML_Char * data
Definition: expat.h:268
std::vector< int > SliceToOrderedNuIdsByEnergy(const std::vector< cheat::NeutrinoWithIndex > &nusWithIdx, const std::vector< std::vector< cheat::NeutrinoEffPur >> &slTruthTable) const
TH2F * ToTH2(TString name) const
Definition: PixelMap.cxx:236
bool passCuts(int cut, const caf::StandardRecord *sr)
Definition: CAFCutter.cxx:37
ParticlesTrainingData * fParticlesTrain
ParentParticleType GetParentParticleType(const art::Ptr< simb::MCTruth > truth, const art::Event &evt)
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
Something else. Tau? Hopefully we don&#39;t use this.
T get(std::string const &key) const
Definition: ParameterSet.h:231
int evt
printf("%d Experimental points found\n", nlines)
recovalid::CAFCutter fCAFCutter
the CAFCutter helper class
float vtxy
FinalStateProngsType GetFinalStateProngsType(const art::Ptr< simb::MCTruth > truth)
bool IsFiltered(const art::Event &evt, art::Ptr< T > x, const std::vector< std::string > &labels)
Is this Ptr marked "filtered out"?
Definition: FilterList.h:96
EventNumber_t event() const
Definition: Event.h:67
TH2F * ToObjTH2(TString name) const
Definition: PixelMap.cxx:285
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
enum cvn::FinalState FinalStateType
TH2F * SingleViewToObjTH2(const unsigned int &view, TString name, bool nuID) const
Definition: PixelMap.cxx:359
enum cvn::Particles ParticlesType
std::string fCAFLabel
label for the process that made the standard records
CVNCosmicEventDump(fhicl::ParameterSet const &pset)
Catch all term, shouldn&#39;t be used.
float vtxz
T * make(ARGS...args) const
void reconfigure(const fhicl::ParameterSet &pset)
enum cvn::ParentParticle ParentParticleType
void train()
Definition: train.C:24
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
The TrainingData objects contains a PixelMap and the output class type, and any other bit that goes i...
Definition: TrainingData.h:23
PixelMap, basic input to CVN neural net.
Definition: PixelMap.h:23
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
Definition: BackTracker.h:133
enum cvn::FinalStateProngs FinalStateProngsType
TH2F * SingleViewToLabTH2(const unsigned int &view, TString name, bool time) const
Definition: PixelMap.cxx:331
Helper class for Reco Validation modules.
Definition: CAFCutter.h:58
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
RunNumber_t run() const
Definition: Event.h:77
bool fObeyPreselection
Check rb::IsFiltered?
int fCutType
what cuts to apply (see CAFCutter.h)
bool fApplyCAFCuts
should we apply the CAF-level cuts?
enum BeamMode string