All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SinglePhoton_module.cc
Go to the documentation of this file.
1 #include "art/Framework/Core/ModuleMacros.h"
2 #include "art/Framework/Core/EDFilter.h"
3 #include "art/Framework/Core/EDAnalyzer.h"
4 #include "art/Framework/Principal/Event.h"
5 #include "art/Framework/Principal/SubRun.h"
6 #include "art/Framework/Principal/Handle.h"
7 #include "art/Framework/Services/Registry/ServiceHandle.h"
8 
9 #include "art_root_io/TFileService.h"
10 #include "art_root_io/TFileDirectory.h"
11 
30 
34 
37 
41 
42 #include "nusimdata/SimulationBase/MCParticle.h"
43 #include "nusimdata/SimulationBase/MCTruth.h"
44 #include "nusimdata/SimulationBase/simb.h"
45 #include "nusimdata/SimulationBase/MCFlux.h"
46 #include "nusimdata/SimulationBase/GTruth.h"
47 
48 #include "canvas/Utilities/ensurePointer.h"
49 #include "canvas/Persistency/Common/FindManyP.h"
50 #include "canvas/Persistency/Common/FindMany.h"
51 #include "canvas/Persistency/Common/FindOneP.h"
52 #include "canvas/Persistency/Common/FindOne.h"
53 
54 #include "fhiclcpp/ParameterSet.h"
55 #include "messagefacility/MessageLogger/MessageLogger.h"
56 #include "cetlib_except/exception.h"
57 
58 //#include "ubobj/Optical/UbooneOpticalFilter.h"
59 
60 // Helper function for PID stuff
61 //#include "ubana/ParticleID/Algorithms/uB_PlaneIDBitsetHelperFunctions.h"
62 
64 
65 #include "TCanvas.h"
66 #include "TTree.h"
67 #include "TFile.h"
68 #include "TGraph.h"
69 #include "TGraph2D.h"
70 #include "TGraphDelaunay.h"
71 #include "TRandom3.h"
72 #include "TGeoPolygon.h"
73 
74 //#include "Pandora/PdgTable.h"
75 #include <chrono>
76 
77 #include <utility>
78 #include <iostream>
79 #include <fstream>
80 #include <string>
81 #include <numeric>
82 #include <algorithm>
83 #include <map>
84 #include <vector>
85 #include <set>
86 #include <sys/stat.h>
87 
97 //#include "sbncode/SinglePhotonAnalysis/Libraries/DBSCAN.h"
99 
100 namespace single_photon
101 {
102 
103  /**
104  * @brief SinglePhoton class
105  */
106  class SinglePhoton : public art::EDFilter
107  {
108  public:
111  /**
112  * @brief Constructor
113  *
114  * @param pset the set of input fhicl parameters
115  */
116  SinglePhoton(fhicl::ParameterSet const &pset);
117 
118  /**
119  * @brief Configure memeber variables using FHiCL parameters
120  *
121  * @param pset the set of input fhicl parameters
122  */
123  void reconfigure(fhicl::ParameterSet const &pset);
124 
125  /**
126  * @brief Analyze an event!
127  *
128  * @param evt the art event to analyze
129  */
130  bool filter(art::Event &evt) override;
131 
132  /**
133  * @brief Begin the job, setting up !
134  *
135  */
136  void beginJob() override;
137 
138  /**
139  * @brief End the job, setting down !
140  *
141  */
142  void endJob() override;
143  /**
144  * @brief: grab run, subrun number, and subrun POT, fill the TTree */
145  bool beginSubRun(art::SubRun& sr) override;
146  bool endSubRun(art::SubRun& sr) override;
147 
148  };
149 
150  //Constructor from .fcl parameters
151  SinglePhoton::SinglePhoton(fhicl::ParameterSet const &pset) : art::EDFilter(pset)
152  {
153  std::cout<<"SinglePhoton::"<<__FUNCTION__<<" kicks off ------------------------------"<<std::endl;
154  this->reconfigure(pset);
155  //Set up some detector, timing, spacecharge and geometry services
156  //Keng, grab theDetector and detClocks in each event.
157  // theDetector = lar::providerFrom<detinfo::DetectorPropertiesService>();
158  // detClocks = lar::providerFrom<detinfo::DetectorClocksService>();
159  paras.s_SCE = lar::providerFrom<spacecharge::SpaceChargeService>();//Get space charge service
160  paras.s_geom = lar::providerFrom<geo::Geometry>();
161  std::cout<<"SinglePhoton::"<<__FUNCTION__<<" now can start jobs --------------------"<<std::endl;
162 
163  }
164 
165  //Reconfigure the internal class parameters from .fcl parameters
166  void SinglePhoton::reconfigure(fhicl::ParameterSet const &pset)
167  {
168  //Get Geometry
169  std::cout<<"SinglePhoton::"<<__FUNCTION__<<" starts ------------------------------"<<std::endl;
170 
171 
172  //input parameters for what file/mode were running in
173  paras.s_print_out_event = pset.get<bool>("PrintOut", false);
174  g_is_verbose = pset.get<bool>("Verbose",false);
175  paras.s_is_data = pset.get<bool>("isData",false);
176  paras.s_is_overlayed = pset.get<bool>("isOverlayed",false);
177  paras.s_is_textgen = pset.get<bool>("isTextGen",false);
178 
179  //some specific additonal info, default not include
180  paras.s_use_PID_algorithms = pset.get<bool>("usePID",false);
181  paras.s_use_delaunay = pset.get<bool>("useDelaunay",false);
182  paras.s_delaunay_max_hits = pset.get<int>("maxDelaunayHits",1000);
183 
184  //When we moved to filter instead of analyzer, kept ability to run 2g filter. A bit depreciated (not in code, but in our use case)
185  paras.s_fill_trees = pset.get<bool>("FillTrees",true);
186  paras.s_run_pi0_filter = pset.get<bool>("RunPi0Filter",false);
187  paras.s_run_pi0_filter_2g1p = pset.get<bool>("FilterMode2g1p",false);
188  paras.s_run_pi0_filter_2g0p = pset.get<bool>("FilterMode2g0p",false);
189 
190  if(paras.s_run_pi0_filter) paras.s_is_data = true;// If running in filter mode, treat all as data
191 
192  //Some output for logging
193  std::cout<<"SinglePhoton::reconfigure || whats configured? "<<std::endl;
194  std::cout<<"SinglePhoton::reconfigure || paras.s_is_data: "<<paras.s_is_data<<std::endl;
195  std::cout<<"SinglePhoton::reconfigure || paras.s_is_overlayed: "<<paras.s_is_overlayed<<std::endl;
196  std::cout<<"SinglePhoton::reconfigure || paras.s_is_textgen: "<<paras.s_is_textgen<<std::endl;
197 
198 
199  //Instead of running over ALL evnts in a file, can pass in a file to run over just the run:subrun:event listed
200  paras.s_runSelectedEvent = pset.get<bool>("SelectEvent", false);
201  paras.s_selected_event_list = pset.get<std::string>("SelectEventList", "");
202 
203  //Studies for photo nuclear EventWeights
204  paras.s_runPhotoNuTruth = pset.get<bool>("RunPhotoNu",false);
205 
206  //Ability to save some FULL eventweight components, rather than run later. Useful for systematic studies. Harcoded to two currently (TODO)
207  paras.s_runTrueEventweight = pset.get<bool>("RunTrueEventWeight",false);
208  paras.s_true_eventweight_label = pset.get<std::string>("true_eventweight_label","eventweight");
209  paras.s_Spline_CV_label = pset.get<std::string>("SplineCVLabel", "eventweight");//4to4aFix");
210 
211 
212  //Input ArtRoot data products
213  paras.s_pandoraLabel = pset.get<std::string>("PandoraLabel");
214  paras.s_trackLabel = pset.get<std::string>("TrackLabel");
215  //KENG paras.s_sliceLabel = pset.get<std::string>("SliceLabel","pandora");
216  paras.s_showerLabel = pset.get<std::string>("ShowerLabel");
217  paras.s_pidLabel = pset.get<std::string>("ParticleIDLabel","pandoracalipidSCE");
218  paras.s_caloLabel = pset.get<std::string>("CaloLabel");
219  paras.s_flashLabel = pset.get<std::string>("FlashLabel");
220  paras.s_potLabel = pset.get<std::string> ("POTLabel");
221  paras.s_hitfinderLabel = pset.get<std::string> ("HitFinderModule", "gaushit");
222  //KENG no such labels in sbnd?
223  //KENG paras.s_badChannelLabel = pset.get<std::string> ("BadChannelLabel","badmasks");
224  //KENG paras.s_showerKalmanLabel = pset.get<std::string> ("ShowerTrackFitter","pandoraKalmanShower");
225  //KENG paras.s_showerKalmanCaloLabel = pset.get<std::string> ("ShowerTrackFitterCalo","pandoraKalmanShowercali");
226  paras.s_generatorLabel = pset.get<std::string> ("GeneratorLabel","generator");
227  //KENG paras.s_mcTrackLabel = pset.get<std::string> ("MCTrackLabel","mcreco");
228  //KENG paras.s_mcShowerLabel = pset.get<std::string> ("MCShowerLabel","mcreco");
229  paras.s_geantModuleLabel = pset.get<std::string> ("GeantModule","largeant");
230  //KENG paras.s_backtrackerLabel = pset.get<std::string> ("BackTrackerModule","gaushitTruthMatch");
231  paras.s_hitMCParticleAssnsLabel = pset.get<std::string> ("HitMCParticleAssnLabel","gaushitTruthMatch");
232  paras.s_shower3dLabel = pset.get<std::string> ("Shower3DLabel","shrreco3d");
233 
234 
235 
236  //Flash related variables.
237  paras.s_beamgate_flash_start = pset.get<double>("beamgateStartTime",3.2); //Defaults to MC for now. Probably should change
238  paras.s_beamgate_flash_end = pset.get<double>("beamgateEndTime",4.8);
239 
240  // paras.s_badChannelProducer = pset.get<std::string>("BadChannelProducer","nfspl1");
241  // paras.s_badChannelProducer = pset.get<std::string>("BadChannelProducer","simnfspl1");
242 
243  //CRT related variables, should run only for RUN3+ enabled
244  paras.s_runCRT = pset.get<bool>("runCRT",false);
245  paras.s_CRTTzeroLabel = pset.get<std::string>("CRTTzeroLabel","crttzero");
246  paras.s_CRTVetoLabel = pset.get<std::string>("CRTVetoLabel","crtveto");
247  paras.s_CRTHitProducer = pset.get<std::string>("CRTHitProducer", "crthitcorr");
248  paras.s_DTOffset = pset.get<double>("DTOffset" , 68600.); //us, taken from ubcrt/UBCRTCosmicFilter/UBCRTCosmicFilter.fcl
249  paras.s_Resolution = pset.get<double>("Resolution" , 1.0); //us, taken from ubcrt/UBCRTCosmicFilter/UBCRTCosmicFilter.fcl
250  // paras.s_DAQHeaderProducer = pset.get<std::string>("DAQHeaderProducer" , "daq");
251 
252  //Some track calorimetry parameters
253  paras.s_track_calo_min_dEdx = pset.get<double>("Min_dEdx",0.005);
254  paras.s_track_calo_max_dEdx = pset.get<double>("Max_dEdx", 30);
255  paras.s_track_calo_min_dEdx_hits = pset.get<double>("Min_dEdx_hits",5); //might be good?
256  paras.s_track_calo_trunc_fraction = pset.get<double>("TruncMeanFraction",20.0);
257 
258  //Some shower calorimetry parameters
259  paras.s_work_function = pset.get<double>("work_function");
260  paras.s_recombination_factor = pset.get<double>("recombination_factor");
261  paras.s_gain_mc = pset.get<std::vector<double>>("gain_mc");
262  paras.s_gain_data = pset.get<std::vector<double>>("gain_data");
263  paras.s_wire_spacing = pset.get<double>("wire_spacing");
264  paras.s_width_dqdx_box = pset.get<double>("width_box");
265  paras.s_length_dqdx_box = pset.get<double>("length_box");
266  paras.s_truthmatching_signaldef = pset.get<std::string>("truthmatching_signaldef");
267 
268  //A seperate mode to run over AllPFPs and not just slice particles
269  paras.s_run_all_pfps = pset.get<bool>("runAllPFPs",false);
270 
271 
272  //Some paramaters for counting protons & photons
273  paras.s_exiting_photon_energy_threshold = pset.get<double>("exiting_photon_energy");
274  paras.s_exiting_proton_energy_threshold = pset.get<double>("exiting_proton_energy");
275  paras.s_max_conv_dist = pset.get<double>("convention_distance_cutoff", 999);
276  paras.s_mass_pi0_mev = 139.57;
277 
278  //SEAviwer Settings for shower clustering and proton stub finding
279  //Have two sets:
280  //Base SEAview is for Second Shower Veto
281  paras.s_runSEAview = pset.get<bool>("runSEAviewShower", false);
282  paras.s_SEAviewHitThreshold = pset.get<double>("SEAviewShowerHitThreshold",25);
283  paras.s_SEAviewPlotDistance = pset.get<double>("SEAviewShowerPlotDistance",80);
284  paras.s_SEAviewDbscanMinPts = pset.get<double>("SEAviewShowerDBSCANMinPts",8);
285  paras.s_SEAviewDbscanEps = pset.get<double>("SEAviewShowerDBSCANEps",4);
286  paras.s_SEAviewMaxPtsLinFit = pset.get<double>("SEAviewShowerMaxHitsLinFit",20.0);
287  paras.s_SEAviewMakePDF = pset.get<bool>("SEAviewShowerMakePDF",false);
288  paras.s_SEAviewNumRecoShower= pset.get<int>("SEAviewShowerNumRecoShower", -1);
289  paras.s_SEAviewNumRecoTrack = pset.get<int>("SEAviewShowerNumRecoTrack", -1);
290 
291  // Second set is for Proton Stub finding
292  paras.s_runSEAviewStub = pset.get<bool>("runSEAviewStub", false);
293  paras.s_SEAviewStubHitThreshold = pset.get<double>("SEAviewStubHitThreshold",25);
294  paras.s_SEAviewStubPlotDistance = pset.get<double>("SEAviewStubPlotDistance",80);
295  paras.s_SEAviewStubDbscanMinPts = pset.get<double>("SEAviewStubDBSCANMinPts",1);
296  paras.s_SEAviewStubDbscanEps = pset.get<double>("SEAviewStubDBSCANEps",1);
297  paras.s_SEAviewStubMakePDF = pset.get<bool>("SEAviewStubMakePDF",false);
298  paras.s_SEAviewStubNumRecoShower =pset.get<int>("SEAviewStubNumRecoShower", -1);
299  paras.s_SEAviewStubNumRecoTrack = pset.get<int>("SEAviewStubNumRecoTrack", -1);
300 
301  paras.s_make_sss_plots = true;
302 
303  //Misc setup CHECK
304  setTPCGeom(paras);//set activeTPC
305  paras.rangen = new TRandom3(22);
306 
307  std::cout<<"SinglePhoton::"<<__FUNCTION__<<" finishes ------------------------------"<<std::endl;
308 
309  }
310 
311  //------------------------------------------------------------------------------------------------------------------------------------------
312 
313 
314 
315  //--------------------------------------- Primary Filter------------------------------------------------------------------------------------
316  // Runs over every artroot event
317  bool SinglePhoton::filter(art::Event &evt)
318  {
319  // //Grab services
320  auto detClocks = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);//it is detinfo::DetectorClocksData
321  auto theDetector = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, detClocks);//it is detinfo::DetectorPropertiesData
322  std::cout<<"SinglePhoton::"<<__FUNCTION__<<" a new event ------------------------------"<<std::endl;
323 
324  //Clear all output branches
325  ClearMeta(vars);
329  ClearStubs(vars);
331  ClearTracks(vars);
336  ClearSlices(vars);
337 
338  Save_EventMeta( evt, vars);
339 
340  //if module is run in selected-event mode, and current event is not in the list, skip it
342  std::cout << "SinglePhoton::analyze()\t||\t event " << vars.m_run_number << "/" << vars.m_subrun_number << "/" << vars.m_event_number << " is not in the list, skip it" << std::endl;
343  return true;
344  }
345 
346  //Timing and TPC info
347  // auto const TPC = (*geom).begin_TPC();
348  // auto ID = TPC.ID();
349  // paras.s_Cryostat = ID.Cryostat;
350  // paras.s_TPC = ID.TPC;
351 
352  //use this in calculating dQdx
353  paras._time2cm = sampling_rate(detClocks) / 1000.0 * theDetector.DriftVelocity( theDetector.Efield(), theDetector.Temperature() );//found in ProtoShowerPandora_tool.cc
354 
355 
356  //******************************Setup*****************Setup**************************************/
357  //Use a new PandoraPFParticle class!
358  //Each class has a core PFParticle member and other associated recob objects;
359 
360 
361  // std::vector<art::Ptr<recob::Track>> trackVector;
362  // art::fill_ptr_vector(trackVector,trackHandle);
363 
364  //Collect the PFParticles from the event. This is the core!
365  art::ValidHandle<std::vector<recob::PFParticle>> const & pfParticleHandle = evt.getValidHandle<std::vector<recob::PFParticle>>(paras.s_pandoraLabel);
366  std::vector<art::Ptr<recob::PFParticle>> pfParticleVector;
367  art::fill_ptr_vector(pfParticleVector,pfParticleHandle);
368  //So a cross check
369  if (!pfParticleHandle.isValid())
370  { mf::LogDebug("SinglePhoton") << " Failed to find the PFParticles.\n";
371  return (paras.s_run_pi0_filter ? false : true) ;
372  }
373 
374  //get the cluster handle for the dQ/dx calc
375  art::ValidHandle<std::vector<recob::Cluster>> const & clusterHandle = evt.getValidHandle<std::vector<recob::Cluster>>(paras.s_pandoraLabel);
376  std::vector< art::Ptr<recob::Cluster> > clusterVector;
377  art::fill_ptr_vector(clusterVector,clusterHandle);
378 
379  // This is another pandora helper. I don't like PFParticle ID lookups but I guess lets keep for now;
380  // typedef std::map< size_t, art::Ptr<recob::PFParticle>>
381  // Produce a map of the PFParticle IDs for fast navigation through the hierarchy
382  // PFParticleIdMap pfParticleMap;
383  // this->GetPFParticleIdMap(pfParticleHandle, pfParticleMap);
384 
385 
386  //And some verticies.
387  art::ValidHandle<std::vector<recob::Vertex>> const & vertexHandle = evt.getValidHandle<std::vector<recob::Vertex>>(paras.s_pandoraLabel);
388  std::vector<art::Ptr<recob::Vertex>> vertexVector;
389  art::fill_ptr_vector(vertexVector,vertexHandle);
390  if(vertexVector.size()>0) vars.m_number_of_vertices++;
391 
392  //PFParticle to Vertices
393  art::FindManyP<recob::Vertex> vertices_per_pfparticle(pfParticleHandle, evt, paras.s_pandoraLabel);
394  std::map< art::Ptr<recob::PFParticle>, std::vector<art::Ptr<recob::Vertex>> > pfParticlesToVerticesMap;
395  for(size_t i=0; i< pfParticleVector.size(); ++i){
396  auto pfp = pfParticleVector[i];
397  pfParticlesToVerticesMap[pfp] =vertices_per_pfparticle.at(pfp.key());
398  }
399 
400  //------- 3D showers
401  art::FindOneP<recob::Shower> showerreco3D_per_pfparticle(pfParticleHandle, evt, paras.s_shower3dLabel);
402 
403  if(g_is_verbose) std::cout<<"SinglePhoton::analyze() \t||\t Get PandoraMetadata"<<std::endl;
404  //add the associaton between PFP and metadata, this is important to look at the slices and scores
405  art::FindManyP< larpandoraobj::PFParticleMetadata > pfPartToMetadataAssoc(pfParticleHandle, evt, paras.s_pandoraLabel);
406 
407  //PFPartciles -->Showers/Tracks
408  art::FindManyP< recob::Track > pfPartToTrackAssoc(pfParticleHandle, evt, paras.s_trackLabel);
409  art::FindManyP< recob::Shower > pfPartToShowerAssoc(pfParticleHandle, evt, paras.s_showerLabel);
410  //Assign Cluster here;
411  art::FindManyP<recob::Cluster> clusters_per_pfparticle(pfParticleHandle, evt, paras.s_pandoraLabel);
412  art::FindManyP<recob::Hit> hits_per_cluster(clusterHandle, evt, paras.s_pandoraLabel);
413 
414  std::map<art::Ptr<recob::PFParticle>, std::vector<art::Ptr<larpandoraobj::PFParticleMetadata>> > pfParticleToMetadataMap;
415  for(size_t i=0; i< pfParticleVector.size(); ++i){
416  const art::Ptr<recob::PFParticle> pfp = pfParticleVector[i];
417  pfParticleToMetadataMap[pfp] = pfPartToMetadataAssoc.at(pfp.key());
418  }
419 
420  //try the new class, Keng
421 
422  std::vector<PandoraPFParticle> allPPFParticles;
423  for(size_t index=0; index< pfParticleVector.size(); ++index){
424  auto pfp = pfParticleVector[index];
425  int temp_key = pfp.key();
426 
427  PandoraPFParticle temp_pf(
428  pfp,
429  pfPartToMetadataAssoc.at(temp_key),
430  vertices_per_pfparticle.at(temp_key),
431  clusters_per_pfparticle.at(temp_key),
432  pfPartToShowerAssoc.at(temp_key),
433  pfPartToTrackAssoc.at(temp_key),
434  hits_per_cluster
435  );
436 
437  allPPFParticles.push_back(temp_pf);
438 
439  }
440  PPFP_FindAncestor(allPPFParticles);
441 
442 
443  //Add slices & hits info.
444  //Slices
445  art::ValidHandle<std::vector<recob::Slice>> const & sliceHandle = evt.getValidHandle<std::vector<recob::Slice>>(paras.s_pandoraLabel);
446  std::vector<art::Ptr<recob::Slice>> sliceVector;
447  art::fill_ptr_vector(sliceVector,sliceHandle);
448 
449  //And some associations
450  art::FindManyP<recob::PFParticle> pfparticles_per_slice(sliceHandle, evt, paras.s_pandoraLabel);
451  art::FindManyP<recob::Hit> hits_per_slice(sliceHandle, evt, paras.s_pandoraLabel);
452 
453  //Slice to PFParticle
454  std::map< art::Ptr<recob::Slice>, std::vector<art::Ptr<recob::PFParticle>> > sliceToPFParticlesMap;
455  std::map<int, std::vector<art::Ptr<recob::PFParticle>> > sliceIDToPFParticlesMap;
456  for(size_t i=0; i< sliceVector.size(); ++i){
457  auto slice = sliceVector[i];
458  PPFP_FindSliceIDandHits(allPPFParticles, slice, pfparticles_per_slice.at(slice.key()), hits_per_slice.at(slice.key()) );
459  sliceToPFParticlesMap[slice] =pfparticles_per_slice.at(slice.key());
460  sliceIDToPFParticlesMap[slice->ID()] = pfparticles_per_slice.at(slice.key());
461  }
462 
463  //Mark all Pandora neutrino slices as neutrino slices and assign the highes neutrino score.
464  vars.pfp_w_bestnuID = DefineNuSlice(allPPFParticles);
465 
466  Save_PFParticleInfo( allPPFParticles, vars, paras);
468 
469  //Slice to Hits
470  std::map<int, std::vector<art::Ptr<recob::Hit>> > sliceIDToHitsMap;
471  for(size_t i=0; i< sliceVector.size(); ++i){
472 
473  auto slice = sliceVector[i];
474  sliceIDToHitsMap[slice->ID()] = hits_per_slice.at(slice.key());
475  }
476 
477 
478  // Once we have actual verticies, lets concentrate on JUST the neutrino PFParticles for now:
479  //--------------------------------
480  // Produce two PFParticle vectors containing final-state particles:
481  // 1. Particles identified as cosmic-rays - recontructed under cosmic-hypothesis
482  // 2. Daughters of the neutrino PFParticle - reconstructed under the neutrino hypothesis
483  std::vector< art::Ptr<recob::PFParticle> > nuParticles;
484  std::vector< art::Ptr<recob::PFParticle> > crParticles;
485 
486  for(size_t jndex=0; jndex< allPPFParticles.size(); ++jndex){
487  PandoraPFParticle* temp_pf = &allPPFParticles[jndex];
488  int temp_id = temp_pf->get_AncestorID();//indentify the jndex for the ancestor PPFParticle
489  PandoraPFParticle* temp_ancestor_ppfp = PPFP_GetPPFPFromPFID( allPPFParticles, temp_id);
490 
491  if ( temp_ancestor_ppfp->get_IsNeutrino() || paras.s_run_all_pfps){
492  nuParticles.push_back( temp_pf->pPFParticle);
493  // std::cout<<"Get nuParticles from pfpID "<<temp_pf->get_PFParticleID()<<std::endl;
494  }else{
495  crParticles.push_back( temp_pf->pPFParticle);
496  }
497  }
498 
499 
500  if(g_is_verbose) std::cout<<"SinglePhoton::analyze() \t||\t Get Spacepoints"<<std::endl;
501  //Spacepoint associaitions
502  art::FindManyP<recob::SpacePoint> spacePoints_per_pfparticle(pfParticleHandle, evt, paras.s_pandoraLabel);
503  std::map<art::Ptr<recob::PFParticle>, std::vector<art::Ptr<recob::SpacePoint>> > pfParticleToSpacePointsMap;
504  for(size_t i=0; i< nuParticles.size(); ++i){
505  const art::Ptr<recob::PFParticle> pfp = nuParticles[i];
506  pfParticleToSpacePointsMap[pfp] = spacePoints_per_pfparticle.at(pfp.key());
507  }
508 
509  if(g_is_verbose) std::cout<<"SinglePhoton::analyze() \t||\t Get Clusters"<<std::endl;
510  //Get a map between the PFP's and the clusters they're imporant for the shower dQ/dx
511  //Also need a map between clusters and hits
512  std::map<art::Ptr<recob::PFParticle>, std::vector<art::Ptr<recob::Cluster>> > pfParticleToClustersMap;
513  std::map<art::Ptr<recob::Cluster>, std::vector<art::Ptr<recob::Hit>> > clusterToHitsMap;
514  //fill map PFP to Clusters
515  for(size_t i=0; i< nuParticles.size(); ++i){
516  auto pfp = nuParticles[i];
517  pfParticleToClustersMap[pfp] = clusters_per_pfparticle.at(pfp.key());
518  }
519  //fill map Cluster to Hits
520  for(size_t i=0; i< clusterVector.size(); ++i){
521  auto cluster = clusterVector[i];
522  clusterToHitsMap[cluster] = hits_per_cluster.at(cluster.key());
523  }
524  if(g_is_verbose) std::cout<<"SinglePhoton::analyze() \t||\t Build hits to PFP Maps"<<std::endl;
525 
526 
527  //OK Here we build two IMPORTANT maps for the analysis, (a) given a PFParticle get a vector of hits..
528  //and (b) given a single hit, get the PFParticle it is in (MARK: is it only one? always? RE-MARK: Yes)
529  std::map<art::Ptr<recob::PFParticle>, std::vector<art::Ptr<recob::Hit>> > pfParticleToHitsMap;
530 
531 
532  //use pfp->cluster and cluster->hit to build pfp->hit map
533  //for each PFP
534  for(size_t i=0; i<nuParticles.size(); ++i){
535  auto pfp = nuParticles[i];
536 
537  //get the associated clusters
538  std::vector<art::Ptr<recob::Cluster>> clusters_vec = pfParticleToClustersMap[pfp] ;
539 
540  //make empty vector to store hits
541  std::vector<art::Ptr<recob::Hit>> hits_for_pfp = {};
542 
543 
544  //for each cluster, get the associated hits
545  for (art::Ptr<recob::Cluster> cluster: clusters_vec){
546  std::vector<art::Ptr<recob::Hit>> hits_vec = clusterToHitsMap[cluster];
547 
548  //insert hits into vector
549  hits_for_pfp.insert( hits_for_pfp.end(), hits_vec.begin(), hits_vec.end() );
550  }
551 
552  //fill the map
553  pfParticleToHitsMap[pfp] = hits_for_pfp;
554 
555  }//for each pfp
556 
557 
558 
559  /**************************************************************************
560  * For SEAview: grab cosmic-related PFPaticles and recob::Hits
561  *
562  **************************************************************************/
563  std::map<art::Ptr<recob::PFParticle>, std::vector<art::Ptr<recob::Cluster>> > cr_pfParticleToClustersMap;
564  std::map<art::Ptr<recob::PFParticle>, std::vector<art::Ptr<recob::Hit>> > cr_pfParticleToHitsMap;
565 
566  //first, collect all daughters of primary cosmic
567  // int nuvars.m_primary_cosmic_particle = crParticles.size();
568  // for(int i =0; i!=nuvars.m_primary_cosmic_particle; ++i){
569  // auto& pParticle = crParticles[i];
570  // for(const size_t daughterId : pParticle->Daughters())
571  // {
572  // if (pfParticleMap.find(daughterId) == pfParticleMap.end())
573  // throw cet::exception("SinglePhoton") << " Invalid PFParticle collection!";
574  //
575  // crParticles.push_back(pfParticleMap.at(daughterId));
576  // }
577  // }
578 
579  //second, build PFP to hits map for cosmic-related PFParticles
580  for(size_t i=0; i< crParticles.size(); ++i){
581  auto pfp = crParticles[i];
582  cr_pfParticleToClustersMap[pfp] = clusters_per_pfparticle.at(pfp.key());
583  }
584 
585  for(size_t i=0; i< crParticles.size(); ++i){
586  auto pfp = crParticles[i];
587 
588  // std::cout<<"starting to match to hits for pfp "<<pfp->Self()<<std::endl;
589  //get the associated clusters
590  std::vector<art::Ptr<recob::Cluster>> clusters_vec = cr_pfParticleToClustersMap[pfp] ;
591 
592  //make empty vector to store hits
593  std::vector<art::Ptr<recob::Hit>> hits_for_pfp = {};
594 
595  // std::cout<<"-- there are "<<clusters_vec.size()<<" associated clusters"<<std::endl;
596 
597  //for each cluster, get the associated hits
598  for (art::Ptr<recob::Cluster> cluster: clusters_vec){
599  std::vector<art::Ptr<recob::Hit>> hits_vec = clusterToHitsMap[cluster];
600 
601  // std::cout<<"looking at cluster in pfp "<<pfp->Self()<<" with "<<hits_vec.size() <<" hits"<<std::endl;
602  //insert hits into vector
603  hits_for_pfp.insert( hits_for_pfp.end(), hits_vec.begin(), hits_vec.end() );
604  }
605 
606  //fill the map
607  cr_pfParticleToHitsMap[pfp] = hits_for_pfp;
608  //std::cout<<"saving a total of "<<hits_for_pfp.size()<<" hits for pfp "<<pfp->Self()<<std::endl;
609 
610  }//for each pfp
611  // std::cout << "Guanqun SEAview test: initial crParticle size: " << nuvars.m_primary_cosmic_particle << ", current size: " << crParticles.size() <<" (including daughters)"<< std::endl;
612  /********************************************************************
613  * End of SEAview: grab cosmic-related PFPaticles and recob::Hits
614  *
615  **************************************************************************/
616 
617 
618  //****************************** Collecting Info and Analysis **************************************/
619 
620  // These are the vectors to hold the tracks and showers for the final-states of the reconstructed neutrino
621  //At this point, nuParticles is a std::vector< art::Ptr<recon::PFParticle>> of the PFParticles that we are interested in.
622  //tracks is a vector of recob::Tracks and same for showers.
623  //Implicitly, tracks.size() + showers.size() = nuParticles.size(); At this point I would like two things.
624  std::vector< art::Ptr<recob::Track> > tracks;
625  std::vector< art::Ptr<recob::Shower> > showers;
626  std::map< art::Ptr<recob::Track> , art::Ptr<recob::PFParticle >> trackToNuPFParticleMap;
627  std::map< art::Ptr<recob::Shower> , art::Ptr<recob::PFParticle>> showerToNuPFParticleMap;
628 
629  if(g_is_verbose) std::cout<<"SinglePhoton::analyze() \t||\t Get Tracks and Showers"<<std::endl;
630 
631  for(size_t jndex=0; jndex< allPPFParticles.size(); ++jndex){
632  PandoraPFParticle* temp_pf = &allPPFParticles[jndex];
633  if(!temp_pf->get_IsNuSlice() ) continue;//pass slice not defined as nu-slice in above lines
634  if(temp_pf->get_HasTrack()){
635  tracks.push_back(temp_pf->pTrack);
636  trackToNuPFParticleMap[temp_pf->pTrack] = temp_pf->pPFParticle;
637  }
638  if(temp_pf->get_HasShower()){
639  showers.push_back(temp_pf->pShower);
640  showerToNuPFParticleMap[temp_pf->pShower] = temp_pf->pPFParticle;
641  }
642  }
643 
644  //Helper function (can be found below) to collect tracks and showers in neutrino slice
645  //this->CollectTracksAndShowers(nuParticles, pfParticleMap, pfParticleHandle, evt, tracks, showers, trackToNuPFParticleMap, showerToNuPFParticleMap);
646 
647  //Track Calorimetry. Bit odd here but bear with me, good to match and fill here
648 
649  //tracks
650  art::ValidHandle<std::vector<recob::Track>> const & trackHandle = evt.getValidHandle<std::vector<recob::Track>>(paras.s_trackLabel);
651  //Keng, use pandoraCalo to get the correct Calorimetry;
652  art::FindManyP<anab::Calorimetry> calo_per_track(trackHandle, evt, paras.s_caloLabel);
653  // std::map<art::Ptr<recob::Track>, std::vector<art::Ptr<anab::Calorimetry>> > trackToCalorimetryMap;
654  //So a cross check
655  if (!calo_per_track.isValid())
656  {
657  mf::LogDebug("SinglePhoton") << " Failed to get Assns between recob::Track and anab::Calorimetry.\n";
658  return (paras.s_run_pi0_filter ? false : true);
659  }
660 
661  art::FindOneP<anab::ParticleID> pid_per_track(trackHandle, evt, paras.s_pidLabel);
662  std::map<art::Ptr<recob::Track>, art::Ptr<anab::ParticleID> > trackToPIDMap;
663  //Loop over all tracks we have to fill calorimetry map
664  for(size_t i=0; i< tracks.size(); ++i){
665  if(calo_per_track.at(tracks[i].key()).size() ==0){
666  std::cerr<<"Track Calorimetry Breaking! the vector of calo_per_track is of length 0 at this track."<<std::endl;
667  }
668  //size_t calo_size = calo_per_track.at(tracks[i].key()).size();
669  //std::cout<<"Track Calo from producer: "<<paras.s_caloLabel<<" has "<<calo_size<<" anab::Calorimetry objects associaed."<<std::endl;
670  //trackToCalorimetryMap[tracks[i]] = calo_per_track.at(tracks[i].key());
671  PandoraPFParticle * temp_ppfp = PPFP_GetPPFPFromTrack(allPPFParticles, tracks[i]);
672  temp_ppfp->set_Calorimetries(calo_per_track.at(tracks[i].key()));
673 
675  temp_ppfp->set_HasPID(true);
676  temp_ppfp->set_ParticleID(pid_per_track.at(tracks[i].key()));
677  }
678  }
679 
680 
681  // If we want PID algorithms to run. do so here
682  // Build a map to get PID from PFParticles, then call PID collection function
683  if(paras.s_use_PID_algorithms){//yes, as set in the FHiCL;
684  for(size_t i=0; i< tracks.size(); ++i){
685  trackToPIDMap[tracks[i]] = pid_per_track.at(tracks[i].key());
686  }
687  }
688 
689 
690 
691  //**********************************************************************************************/
692  //**********************************************************************************************/
693  //---------------------------------- MC TRUTH, MC Only---------------------------
694  //**********************************************************************************************/
695  //**********************************************************************************************/
696 
697  //Get the MCtruth handles and vectors
698  std::vector<art::Ptr<simb::MCTruth>> mcTruthVector;
699  std::vector<art::Ptr<simb::MCParticle>> mcParticleVector;
700 
701  //Then build a map from MCparticles to Hits and vice versa
702  std::map< art::Ptr<simb::MCParticle>, std::vector<art::Ptr<recob::Hit> > > mcParticleToHitsMap;
703  std::map< art::Ptr<recob::Hit>, art::Ptr<simb::MCParticle> > hitToMCParticleMap;
704 
705  //Apparrently a MCParticle doesn't know its origin (thanks Andy!)
706  //I would also like a map from MCparticle to MCtruth and then I will be done. and Vice Versa
707  //Note which map is which! //First is one-to-many. //Second is one-to-one
708  std::map< art::Ptr<simb::MCTruth>, std::vector<art::Ptr<simb::MCParticle>>> MCTruthToMCParticlesMap;
709  std::map< art::Ptr<simb::MCParticle>, art::Ptr<simb::MCTruth>> MCParticleToMCTruthMap;
710  std::map<int, art::Ptr<simb::MCParticle> > MCParticleToTrackIdMap;
711 
712  std::vector<art::Ptr<sim::MCTrack>> mcTrackVector;
713  std::vector<art::Ptr<sim::MCShower>> mcShowerVector;
714 
715  std::vector<art::Ptr<simb::MCParticle>> matchedMCParticleVector;
716  std::map<art::Ptr<recob::Track>, art::Ptr<simb::MCParticle> > trackToMCParticleMap;
717  std::map<art::Ptr<recob::Shower>, art::Ptr<simb::MCParticle> > showerToMCParticleMap;
718 
719  //Given a simb::MCParticle we would like a map to either a sim::MCTrack or sim::MCShower
720  std::map< art::Ptr<simb::MCParticle>, art::Ptr<sim::MCTrack> > MCParticleToMCTrackMap;
721  std::map< art::Ptr<simb::MCParticle>, art::Ptr<sim::MCShower> > MCParticleToMCShowerMap;
722 
723  if(g_is_verbose){
724  std::cout << "SinglePhoton::analyze()\t||\t Consolidated event summary:" << "\n";
725  std::cout << "SinglePhoton::analyze()\t||\t - Number of primary cosmic-ray PFParticles : " << crParticles.size() << "\n";
726  std::cout << "SinglePhoton::analyze()\t||\t - Number of neutrino final-state PFParticles : " << nuParticles.size() << "\n";
727  std::cout << "SinglePhoton::analyze()\t||\t ... of which are track-like : " << tracks.size() << "\n";
728  std::cout << "SinglePhoton::analyze()\t||\t ... of which are showers-like : " << showers.size() << "\n";
729  }
730 
731  //**********************************************************************************************/
732  //**********************************************************************************************/
733 
734  //and now get the simb::MCparticle to both MCtrack and MCshower maps (just for the MCparticles matched ok).
735  // if(!paras.s_run_pi0_filter) badChannelMatching<art::Ptr<recob::Track>>(badChannelVector, tracks, trackToNuPFParticleMap, pfParticleToHitsMap,geom,bad_channel_list_fixed_mcc9);
736  // badChannelMatching<art::Ptr<recob::Track>>(badChannelVector, tracks, trackToNuPFParticleMap, pfParticleToHitsMap,geom,bad_channel_list_fixed_mcc9);
737 
738 
739  //*******************************Slices***************************************************************/
740 
741  //these are all filled in analyze slice
742  // std::vector<std::pair<art::Ptr<recob::PFParticle>,int>> allPFPSliceIdVec; //stores a pair of all PFP's in the event and the slice ind
743  std::vector<std::pair<art::Ptr<recob::PFParticle>,int>> primaryPFPSliceIdVec; //stores a pair of only the primary PFP's in the event and the slice ind
744  std::map<int, double> sliceIdToNuScoreMap; //map between a slice Id and neutrino score
745  std::map<art::Ptr<recob::PFParticle>, bool> PFPToClearCosmicMap; //returns true for clear cosmic, false otherwise
746  // std::map<art::Ptr<recob::PFParticle>, int> PFPToSliceIdMap; //returns the slice id for all PFP's
747  std::map<art::Ptr<recob::PFParticle>,bool> PFPToNuSliceMap;
748  std::map<art::Ptr<recob::PFParticle>,double> PFPToTrackScoreMap;
749  std::map<int, int> sliceIdToNumPFPsMap;
750 
751  if(g_is_verbose) std::cout<<"\nSinglePhoton::AnalyzeSlice()\t||\t Starting"<<std::endl;
752 
753 
754 
755  //******************************* CRT CRT***************************************************************/
756  //Optical Flashes
757  art::ValidHandle<std::vector<recob::OpFlash>> const & flashHandle = evt.getValidHandle<std::vector<recob::OpFlash>>(paras.s_flashLabel);
758  std::vector<art::Ptr<recob::OpFlash>> flashVector;
759  art::fill_ptr_vector(flashVector,flashHandle);
760 
761  std::map<art::Ptr<recob::OpFlash>, std::vector< art::Ptr<sbn::crt::CRTHit>>> crtvetoToFlashMap;
762 
763  if(paras.s_runCRT){
764  art::FindManyP<sbn::crt::CRTHit> crtveto_per_flash(flashHandle, evt, paras.s_CRTVetoLabel);
765  for(size_t i=0; i< flashVector.size(); ++i){
766  crtvetoToFlashMap[flashVector[i]] = crtveto_per_flash.at(flashVector[i].key());
767  }
768  }
769 
770  art::Handle<std::vector<sbn::crt::CRTHit>> crthit_h; //only filled when there are hits, otherwise empty
771  //art::Handle<raw::DAQHeaderTimeUBooNE> rawHandle_DAQHeader;//Keng, this is to track CRT hit time;
772  double evt_timeGPS_nsec = -999 ;
773  if(paras.s_runCRT){
774  //evt.getByLabel(paras.s_DAQHeaderProducer, rawHandle_DAQHeader);
775 
776 
777  // raw::DAQHeaderTimeUBooNE const& my_DAQHeader(*rawHandle_DAQHeader);
778  //Keng, time() is set to the trigger time.
779  art::Timestamp evtTimeGPS = evt.time();
780  evt_timeGPS_nsec = evtTimeGPS.timeLow();
781 
782  evt.getByLabel(paras.s_CRTHitProducer, crthit_h);
783  std::cout<<"SinglePhoton::analyze \t||\t Got CRT hits"<<std::endl;
784  }
785 
786  //Analize the CRT flashes and store them in the vertex_tree.
787  //Found in analyze_OpFlashes.h
788  ResizeFlashes(flashVector.size(), vars);
789  AnalyzeFlashes(flashVector, crthit_h, evt_timeGPS_nsec, crtvetoToFlashMap, vars, paras);
790 
791 
792  //******************************* Common Optical Filter **************************************************************/
793  //Raw Optical fltr
794 
795  // art::Handle<uboone::UbooneOpticalFilter> uBooNE_common_optFltr;
796  // if(evt.getByLabel("opfiltercommon", uBooNE_common_optFltr)){
797  // vars.m_flash_optfltr_pe_beam = uBooNE_common_optFltr->PE_Beam();
798  // vars.m_flash_optfltr_pe_beavars.m_tot = uBooNE_common_optFltr->PE_Beavars.m_Total();
799  // vars.m_flash_optfltr_pe_veto = uBooNE_common_optFltr->PE_Veto();
800  // vars.m_flash_optfltr_pe_veto_tot = uBooNE_common_optFltr->PE_Veto_Total();
801  // }else{
802  // vars.m_flash_optfltr_pe_beam = -999;
803  // vars.m_flash_optfltr_pe_beavars.m_tot = -999;
804  // vars.m_flash_optfltr_pe_veto = -999;
805  // vars.m_flash_optfltr_pe_veto_tot = -999;
806  // std::cout<<"No opfiltercommon product:"<<std::endl;
807  // }
808 
809 
810  //******************************* Tracks **************************************************************/
811  std::cout<<"\nSinglePhoton::analyze \t||\t Start on Track Analysis "<<std::endl;
812 
813  //found in analyze_Tracks.h
815  allPPFParticles,
816  tracks,
817  pfParticleToSpacePointsMap, MCParticleToTrackIdMap, sliceIdToNuScoreMap,
818  vars, paras);
819  AnalyzeTrackCalo(tracks, allPPFParticles, vars, paras);
820 
821 
822  //Run over PID?
823  if(paras.s_use_PID_algorithms) CollectPID(tracks, allPPFParticles, vars);
824 
825 
826  //******************************* Showers **************************************************************/
827  std::cout<<"\nSinglePhoton::analyze \t||\t Start on Shower Analysis "<<std::endl;
828 
829 
830  //found in analyze_Showers.h
831  vars.m_reco_asso_showers=showers.size();
833  AnalyzeShowers(allPPFParticles, showers,
834  clusterToHitsMap,
835  trigger_offset(detClocks),
836  theDetector,
837  vars,
838  paras);
839 
840  std::cout<<"\nSinglePhoton::analyze \t||\t Finish Shower Analysis "<<std::endl;
841 
842  std::cout<<std::endl;
843  std::cout<<"SinglePhoton::Reconstruction Summary ---------------------------------------------- "<<std::endl;
844  std::vector<int> spacers = Printer_header({" pfpID", " Parent ", "Vertex(x, "," y, ",", z )", "slice "," nu_score"," trk_score "," prim? "," nu? "," nuSlice?"," cos? "," S# "," T# "});
845  for(size_t jndex=0; jndex< allPPFParticles.size(); ++jndex){
846  PandoraPFParticle temp_pf = allPPFParticles[jndex];
847  art::Ptr<recob::PFParticle> pfp = temp_pf.pPFParticle;
848 
850  {std::to_string(pfp->Self()),
851  std::to_string(pfp->Parent()) ,
852  std::to_string(temp_pf.get_Vertex_pos()[0]),
853  std::to_string(temp_pf.get_Vertex_pos()[1]),
854  std::to_string(temp_pf.get_Vertex_pos()[2]),
855  std::to_string(temp_pf.get_SliceID()),
856  std::to_string(temp_pf.get_NuScore()),
857  std::to_string(temp_pf.get_TrackScore()),
858  (pfp->IsPrimary())? "T":" ",
859  (temp_pf.get_IsNeutrino() )? "T":" ",
860  (temp_pf.get_IsNuSlice() )? "T":" ",
861  (temp_pf.get_IsClearCosmic() )? "T":" ",
862  (temp_pf.get_HasShower() >0)? "1":" ",
863  (temp_pf.get_HasTrack() >0)? "1":" "
864  }, spacers);
865 
866  }
867 
868  //KENG no KalmanShowers in SBND?
869  // this->AnalyzeKalmanShowers(showers,showerToNuPFParticleMap,pfParticlesToShowerKalmanMap, kalmanTrackToCaloMap, pfParticleToHitsMap, theDetector);
870 
871 
872  //Some misc things thrown in here rather than in a proper helper function. TODO. fix
873  //Calc a fake shower "end" distance. How to define an end distance? good question
874  for(size_t i_shr = 0; i_shr<showers.size();i_shr++){
875  const art::Ptr<recob::Shower> s = showers[i_shr];
876  const art::Ptr<recob::PFParticle> pfp = showerToNuPFParticleMap[s];
877  const std::vector< art::Ptr<recob::SpacePoint> > shr_spacepoints = pfParticleToSpacePointsMap[pfp];
878 
880  vars.m_reco_shower_end_dist_to_SCB[i_shr] = 99999;
881 
882  for(auto &sp: shr_spacepoints){
883  std::vector<double> tmp_spt = {sp->XYZ()[0],sp->XYZ()[1] , sp->XYZ()[2]};
885  double tmo;
886  distToSCB(tmo,tmp_spt, paras);
888 
889  //This section runs for only 1 shower events for purpose of testing delta specifics
890  if(showers.size()==1){
891  vars.m_reco_shower_spacepoint_x.push_back(sp->XYZ()[0]);
892  vars.m_reco_shower_spacepoint_y.push_back(sp->XYZ()[1]);
893  vars.m_reco_shower_spacepoint_z.push_back(sp->XYZ()[2]);
894  }
895 
896  }
897  }
898 
899  //******************************* MCTruth **************************************************************/
900  // Collect all the hits. We will need these. Lets grab both the handle as well as a vector of art::Ptr as I like both.
901  art::ValidHandle<std::vector<recob::Hit>> const & hitHandle = evt.getValidHandle<std::vector<recob::Hit>>(paras.s_hitfinderLabel);
902  std::vector<art::Ptr<recob::Hit>> hitVector;
903  art::fill_ptr_vector(hitVector,hitHandle);
904 
905 
906  //Grab the backtracker info for MCTruth Matching
907  art::FindManyP<simb::MCParticle,anab::BackTrackerHitMatchingData> mcparticles_per_hit(hitHandle, evt, paras.s_hitMCParticleAssnsLabel);
908 
909  // MCTruth, MCParticle, MCNeutrino information all comes directly from GENIE.
910  // MCShower and MCTrack come from energy depositions in GEANT4
911 
912  //Only run if its not data, i.e. MC events :)
913  if(!paras.s_is_data){
914 //GTruth information, don't need this at the moment;
915 // if(!paras.s_is_textgen){
916 // std::vector<art::Ptr<simb::GTruth>> gTruthVector;
917 // std::cout<<"\n Get some GTruth info."<<std::endl;
918 // // if Text is in the generator label, skip it. TODO this is a bit simple but works, maybe add a boolean
919 // art::ValidHandle<std::vector<simb::GTruth>> const & gTruthHandle= evt.getValidHandle<std::vector<simb::GTruth>>(paras.s_generatorLabel);
920 // art::fill_ptr_vector(gTruthVector,gTruthHandle);
921 // if(g_is_verbose){
922 // for(size_t p=0; p< gTruthVector.size();p++) std::cout<<gTruthVector[p]<<" "<<*gTruthVector[p]<<std::endl;
923 // }
924 // std::cout<<"done"<<std::endl;
925 // }
926 
927  //get MCTruth (GENIE)
928  art::ValidHandle<std::vector<simb::MCTruth>> const & mcTruthHandle= evt.getValidHandle<std::vector<simb::MCTruth>>(paras.s_generatorLabel);
929  art::fill_ptr_vector(mcTruthVector,mcTruthHandle);
930 
931  //get MCPartilces (GEANT4)
932  art::ValidHandle<std::vector<simb::MCParticle>> const & mcParticleHandle= evt.getValidHandle<std::vector<simb::MCParticle>>(paras.s_geantModuleLabel);
933  art::fill_ptr_vector(mcParticleVector,mcParticleHandle);
934 
935  //Found inanalyze_Geant4.h
936  //Currently just saves first 2 particles. TODO have a input list of things to save. Dont want to save everything!!
937 
938  if(g_is_verbose) std::cout<<"SinglePhoton::AnalyzeGeant4()\t||\t Begininning recob::Geant4 analysis suite"<<std::endl;;
939  AnalyzeGeant4(mcParticleVector, vars);
940 
941 
942  //Get the MCParticles (move to do this ourselves later)
943  CollectMCParticles(evt, paras.s_geantModuleLabel, MCTruthToMCParticlesMap, MCParticleToMCTruthMap, MCParticleToTrackIdMap,vars);
944 
945 
946  //mcc9 march miniretreat fix
947  std::vector<art::Ptr<simb::MCParticle>> particle_vec; //vector of all MCParticles associated with a given hit in the reco PFP
948  std::vector<anab::BackTrackerHitMatchingData const *> match_vec; //vector of some backtracker thing
950 
951  for(size_t j=0; j<hitVector.size();j++){
952  const art::Ptr<recob::Hit> hit = hitVector[j];
953 
954  particle_vec.clear(); match_vec.clear(); //only store per hit
955  mcparticles_per_hit.get(hit.key(), particle_vec, match_vec);
956  if(particle_vec.size() > 0){
958  }
959  }
960 
961 
962  //Important map, given a MCparticle, whats the "hits" associated
963  BuildMCParticleHitMaps(evt, paras.s_geantModuleLabel, hitVector, mcParticleToHitsMap, hitToMCParticleMap, lar_pandora::LArPandoraHelper::kAddDaughters, MCParticleToTrackIdMap, vars);
964 
965 
966  //The recoMC was originally templated for any track shower, but sufficient differences in showers emerged to have stand alone sadly
967  std::cout<<"\nSinglePhoton\t||\t Starting backtracker on recob::track"<<std::endl;
968  std::vector<double> trk_overlay_vec = trackRecoMCmatching( tracks,
969  trackToMCParticleMap,
970  trackToNuPFParticleMap,
971  pfParticleToHitsMap,
972  mcparticles_per_hit,
973  matchedMCParticleVector,
974  vars);
975 
976  std::cout<<"\nSinglePhoton\t||\t Starting backtracker on recob::shower"<<std::endl;
977  if(g_is_verbose) std::cout<<"Matching "<<showers.size()<<" showers"<<std::endl;
979  allPPFParticles,
980  showers,
981  showerToMCParticleMap,
982  mcparticles_per_hit,
983  matchedMCParticleVector,
984  MCParticleToTrackIdMap,
985  vars);
986 
987  //photoNuclearTesting(matchedMCParticleVector);
988 
989  std::cout<<"\nSinglePhoton\t||\tStarting RecoMCTracks "<<std::endl;
990  RecoMCTracks(
991  allPPFParticles,
992  tracks,
993  trackToMCParticleMap,
994  MCParticleToMCTruthMap,
995  mcParticleVector,
996  MCParticleToTrackIdMap,
997  trk_overlay_vec,
998  vars);
999 
1000  std::cout<<"\nSinglePhoton\t||\tStarting AnalyzeMCTruths"<<std::endl;
1001 
1002  ResizeMCTruths(mcTruthVector[0]->NParticles(), vars);//assume only one MCTruth
1003  AnalyzeMCTruths(mcTruthVector, mcParticleVector, vars, paras);
1004 
1005 
1006  if(!paras.s_is_textgen){// if Text is in the generator label, skip it. wont contain any
1007 
1008  //This is important for rebuilding MCFlux and GTruth to do eventweight
1009  std::cout<<"\nSinglePhoton\t||\tStarting AnalyzeEventWeight"<<std::endl;
1010  AnalyzeEventWeight(evt, vars);
1011 
1012 
1013  std::cout<<"\nSinglePhoton\t||\tStarting AnalyzeRecoMCSlices"<<std::endl;
1014  AnalyzeRecoMCSlices(paras.s_truthmatching_signaldef, allPPFParticles, MCParticleToTrackIdMap, showerToMCParticleMap, trackToMCParticleMap, vars, paras);
1015 
1016  std::cout<<"\nSinglePhoton\t||\tFinish AnalyzeRecoMCSlices"<<std::endl;
1017 
1018  //******************************* PhotoNuclear Absorption **************************************************************/
1019  vars.m_photonu_weight_low = -999;
1020  vars.m_photonu_weight_high = -999;
1022  art::Handle<std::vector<evwgh::MCEventWeight>> ev_evw_ph ;
1023  if( evt.getByLabel("eventweight",ev_evw_ph)){
1024  std::map<std::string, std::vector<double>> const & weight_map = ev_evw_ph->front().fWeight;
1025  for (auto const& x : weight_map){
1026  std::cout << x.first // string (key)
1027  << ':'
1028  << x.second.size() << std::endl ;
1029  if(x.first == "photonuclear_photon_PhotoNuclear"){
1030  auto vec = x.second;
1031  double ph_low = vec[1];
1032  double ph_high = vec[0];
1033  std::cout<<"PhotoNuBit: "<<ph_low<<" "<<ph_high<<std::endl;
1034  vars.m_photonu_weight_low = ph_low;
1035  vars.m_photonu_weight_high = ph_high;
1036  }
1037  }
1038 
1039  }
1040  }
1041 
1042 
1043  //******************************* True EventWeight **************************************************************/
1044  //Remove this
1045  // if(paras.s_runTrueEventweight){
1046  //
1047  // art::ValidHandle<std::vector<evwgh::MCEventWeight>> const & ev_evw_true = evt.getValidHandle<std::vector<evwgh::MCEventWeight>>(paras.s_true_eventweight_label);
1048  // std::map<std::string, std::vector<double>> const & weight_map = ev_evw_true->front().fWeight;
1049  // if(ev_evw_true->size() > 1) {
1050  // std::cout << __LINE__ << " " << __PRETTY_FUNCTION__ << "\n"
1051  // << "WARNING: eventweight has more than one entry\n";
1052  // }
1053  // fmcweight=weight_map;
1054  // }
1055  //
1056  }//end NOT textgen
1057 
1058 
1059 
1060  std::cout<<"SinglePhoton::analyze\t||\t finnished loop for this event"<<std::endl;
1061  }//end MC loop
1062 
1063 
1064 
1065  //******************************* Isolation (SSV precursor) **************************************************************/
1067  std::cout<<"\nSinglePhoton::analyze \t||\t Start IsolationStudy "<<std::endl;
1068  IsolationStudy(allPPFParticles, tracks, showers, theDetector, vars, paras);
1069  std::cout<<"\nSinglePhoton::analyze \t||\t Finish IsolationStudy "<<std::endl;
1070  }
1071  // if(!vars.m_run_all_pfps && ! paras.s_run_pi0_filter) this->IsolationStudy(allPPFParticles, tracks, trackToNuPFParticleMap, showers, showerToNuPFParticleMap, pfParticleToHitsMap, PFPToSliceIdMap, sliceIDToHitsMap, theDetector);
1072 
1073 
1074 
1075  // ################################################### SEAview SEAview #########################################################
1076  // ################################################### Proton Stub ###########################################
1077  // ------------- stub clustering ---------------------------
1078  std::cout << "----------------- Stub clustering --------------------------- " << std::endl;
1079  // std::cout << "SEAview Stub formation: " << (paras.s_runSEAviewStub ? "true" : "false" ) << " nshower requirement: " << paras.s_SEAviewStubNumRecoShower << ", actual num shower: " << showers.size() << " | ntrack requirement: " << paras.s_SEAviewStubNumRecoTrack << ", actual num track: " << tracks.size() << std::endl;
1080 
1081  if( !paras.s_run_pi0_filter &&
1083  (paras.s_SEAviewStubNumRecoShower== -1 || (int)showers.size()== paras.s_SEAviewStubNumRecoShower) &&
1084  (paras.s_SEAviewStubNumRecoTrack == -1 || (int)tracks.size() == paras.s_SEAviewStubNumRecoTrack)){
1085 
1086  // grab all hits in the slice of the reco shower
1087  art::Ptr<recob::Shower> p_shr = showers.front();
1088  PandoraPFParticle* ppfp = PPFP_GetPPFPFromShower(allPPFParticles, p_shr);
1089  art::Ptr<recob::PFParticle> p_pfp = ppfp->pPFParticle;//showerToNuPFParticleMap[p_shr];
1090  std::vector<art::Ptr<recob::Hit>> p_hits = pfParticleToHitsMap[p_pfp];
1091 
1092  int p_sliceid = ppfp->get_SliceID();//PFPToSliceIdMap[p_pfp];
1093  auto p_slice_hits = sliceIDToHitsMap[p_sliceid];
1094 
1095  std::string uniq_tag = "HitThres_"+ std::to_string(static_cast<int>(paras.s_SEAviewStubHitThreshold)) + "_" + std::to_string(vars.m_run_number)+"_"+std::to_string(vars.m_subrun_number)+"_"+std::to_string(vars.m_event_number);
1096 
1097  //Setup seaviewr object
1098  ////CHECK undefined reference? Found in header, but not defined;
1099  //Solution: add MODULE_LIBRARIES in cmake;
1100  seaview::SEAviewer sevd("Stub_"+uniq_tag,paras.s_geom, theDetector);
1101  //Pass in any bad channels you like
1102  // sevd.setBadChannelList(bad_channel_list_fixed_mcc9);
1103  //Give it a vertex to center around
1105 
1106  //Add the hits from just this slice, as well as hits within 150cm of the vertex
1107  sevd.addHitsToConsider(hitVector); // std::vector<art::Ptr<recob::Hit>>
1108  sevd.filterConsideredHits(150); //remve hits that're not within 150cm of the vertex on 2D view
1109  sevd.addHitsToConsider(p_slice_hits);
1110 
1111  sevd.addAllHits(hitVector); // std::vector<art::Ptr<recob::Hit>>
1113 
1114 
1115  //Add all the "nice " PFParticle Hits, as well as what to label
1116  //sevd.addPFParticleHits(p_hits, "Shower"); //std::vector<art::Ptr<recob::Hit>> and std::string
1117  sevd.addPFParticleHits(p_hits, "Shower", vars.m_reco_shower_energy_max[0], vars.m_reco_shower_conversion_distance[0], vars.m_reco_shower_impact_parameter[0]); //std::vector<art::Ptr<recob::Hit>> and std::string
1118 
1119  //and add the SingleShower we like
1120  sevd.addShower(p_shr); // art::Ptr<recob::Shower>
1121 
1122  //Add all track PFP
1123  int i_trk = 0;
1124  for(auto &trk: tracks){
1125  art::Ptr<recob::PFParticle> p_pfp_trk = trackToNuPFParticleMap[trk];
1126  std::vector<art::Ptr<recob::Hit>> p_hits_trk = pfParticleToHitsMap[p_pfp_trk];
1127  //sevd.addPFParticleHits(p_hits_trk,"track");
1128  sevd.addPFParticleHits(p_hits_trk,"track", vars.m_reco_track_length[i_trk], vars.m_reco_track_spacepoint_principal0[i_trk]);
1129  sevd.addTrack(trk);
1130  ++i_trk;
1131  }
1132 
1133  //Add all cosmic-relatd PFP
1134  for(auto &cr: crParticles){
1135  std::vector<art::Ptr<recob::Hit>> p_hits_cr = cr_pfParticleToHitsMap[cr];
1136  sevd.addPFParticleHits(p_hits_cr,"cosmic");
1137  }
1138 
1139 
1140  //We then calculate Unassociated hits, i.e the hits not associated to the "Shower" or tracksyou passed in.
1141  auto vnh= sevd.calcUnassociatedHits();
1142  vars.m_trackstub_num_unassociated_hits =vnh[1]+vnh[2];
1144  vars.m_trackstub_associated_hits = vnh[0]-vnh[1]-vnh[2];
1145 
1146  //Recluster, group unassociated hits into different clusters
1148 
1149  //And some plotting
1150  // If we want to plot pdfs again later, then we can't plot here
1151  //if(vars.m_SEAviewStubMakePDF) sevd.Print(vars.m_SEAviewStubPlotDistance);
1152 
1153  //Analyze formed clusters and save info
1154  std::vector<seaview::cluster> vec_SEAclusters ;
1155  sevd.analyzeTrackLikeClusters(paras.s_SEAviewStubDbscanEps, showerToNuPFParticleMap, pfParticleToHitsMap, vec_SEAclusters);
1156 
1157 
1158  //And save to file.
1159  std::cout<<"After SEAview we have "<<vec_SEAclusters.size()<<" Stub clusters to chat about"<<std::endl;
1160 
1162  for(size_t c=0; c< vec_SEAclusters.size(); c++){
1163  auto& clu = vec_SEAclusters.at(c); //type: seaview::cluster
1164  int pl = clu.getPlane();
1165  auto hitz = clu.getHits();
1166  double Ep = CalcEShowerPlane(hitz,pl, paras);
1167  int remerge = clu.getShowerRemerge();
1168  seaview::cluster_score * ssscorz = clu.getScore();
1169 
1170  std::cout<<c<<" "<<pl<<" "<<Ep<<" "<<clu.getMinHitImpactParam()<<" "<<clu.getMinHitConvDist()<<" "<<clu.getMinHitIOC()<<" "<<clu.getMeanADC()<<" "<<clu.getADCrms()<<" "<< clu.getLinearChi() << " " << remerge<<std::endl;
1171 
1172  //if the cluster is too close to the recob::shower, then do not include it
1173  if(remerge>=0 && remerge< (int)vars.m_reco_shower_reclustered_energy_plane2.size()){
1174  //decide not to add energy of the cluster to reco shower if it's matched
1175  //
1176  //if(pl==0)vars.m_reco_shower_reclustered_energy_plane0[remerge]+=Ep;
1177  //if(pl==1)vars.m_reco_shower_reclustered_energy_plane1[remerge]+=Ep;
1178  //if(pl==2)vars.m_reco_shower_reclustered_energy_plane2[remerge]+=Ep;
1179 
1180  continue;// Dont include this as a viable cluster!
1181  }
1182 
1184  //determine if this cluster is in neutrino slice
1185  vars.m_trackstub_candidate_in_nu_slice.push_back(clu.InNuSlice(sliceIDToHitsMap, p_sliceid));
1186 
1187  //Fill All the bits
1188  vars.m_trackstub_candidate_num_hits.push_back((int)hitz.size());
1189  vars.m_trackstub_candidate_num_wires.push_back((int)ssscorz->n_wires);
1190  vars.m_trackstub_candidate_num_ticks.push_back((int)ssscorz->n_ticks);
1191  vars.m_trackstub_candidate_plane.push_back(pl);
1192  vars.m_trackstub_candidate_PCA.push_back(ssscorz->pca_0);
1193  vars.m_trackstub_candidate_mean_tick.push_back(ssscorz->mean_tick);
1194  vars.m_trackstub_candidate_max_tick.push_back(ssscorz->max_tick);
1195  vars.m_trackstub_candidate_min_tick.push_back(ssscorz->min_tick);
1196  vars.m_trackstub_candidate_min_wire.push_back(ssscorz->min_wire);
1197  vars.m_trackstub_candidate_max_wire.push_back(ssscorz->max_wire);
1198  vars.m_trackstub_candidate_mean_wire.push_back(ssscorz->mean_wire);
1199  vars.m_trackstub_candidate_min_dist.push_back(ssscorz->min_dist);
1200  vars.m_trackstub_candidate_min_impact_parameter_to_shower.push_back(clu.getMinHitImpactParam());
1201  vars.m_trackstub_candidate_min_conversion_dist_to_shower_start.push_back(clu.getMinHitConvDist());
1202  vars.m_trackstub_candidate_min_ioc_to_shower_start.push_back(clu.getMinHitIOC());
1203  vars.m_trackstub_candidate_ioc_based_length.push_back(clu.getIOCbasedLength());
1204  vars.m_trackstub_candidate_wire_tick_based_length.push_back(clu.getWireTickBasedLength());
1205  vars.m_trackstub_candidate_mean_ADC_first_half.push_back(clu.getMeanADCFirstHalf());
1206  vars.m_trackstub_candidate_mean_ADC_second_half.push_back(clu.getMeanADCSecondHalf());
1207  vars.m_trackstub_candidate_mean_ADC_first_to_second_ratio.push_back(clu.getMeanADCRatio());
1208  vars.m_trackstub_candidate_track_angle_wrt_shower_direction.push_back(clu.getTrackAngleToShowerDirection());
1209  vars.m_trackstub_candidate_linear_fit_chi2.push_back(clu.getLinearChi());
1210  vars.m_trackstub_candidate_mean_ADC.push_back(clu.getMeanADC());
1211  vars.m_trackstub_candidate_ADC_RMS.push_back(clu.getADCrms());
1212  vars.m_trackstub_candidate_energy.push_back(Ep);
1213  vars.m_trackstub_candidate_remerge.push_back(remerge);
1214 
1215 
1216  //MCTruth matching for pi0's
1217  if(paras.s_is_data){
1218  vars.m_trackstub_candidate_matched.push_back(-1);
1219  vars.m_trackstub_candidate_pdg.push_back(-1);
1220  vars.m_trackstub_candidate_parent_pdg.push_back(-1);
1221  vars.m_trackstub_candidate_trackid.push_back(-1);
1225  }else{
1226 
1227  auto ssmatched = SecondShowerMatching(hitz, mcparticles_per_hit, mcParticleVector, MCParticleToTrackIdMap, vars);
1228  vars.m_trackstub_candidate_matched.push_back(ssmatched[0]);
1229  vars.m_trackstub_candidate_pdg.push_back(ssmatched[1]);
1230  vars.m_trackstub_candidate_parent_pdg.push_back(ssmatched[2]);
1231  vars.m_trackstub_candidate_trackid.push_back(ssmatched[3]);
1232  vars.m_trackstub_candidate_true_energy.push_back(ssmatched[4]);
1233  vars.m_trackstub_candidate_overlay_fraction.push_back(ssmatched[5]);
1235 
1236  //Guanqun: print out (best-matched) truth information of the cluster
1237  std::cout << "Cluster: " << vars.m_trackstub_num_candidates-1 << " plane: " << vars.m_trackstub_candidate_plane.back() << ", energy: " << vars.m_trackstub_candidate_energy.back() << ", min IOC of hit(wrt shower): " << vars.m_trackstub_candidate_min_ioc_to_shower_start.back() << "\n";
1238  std::cout << "Cluster is matched: " << vars.m_trackstub_candidate_matched.back() << ", matched PDG: " << vars.m_trackstub_candidate_pdg.back() << " track ID: " << vars.m_trackstub_candidate_trackid.back() << " overlay fraction: " << vars.m_trackstub_candidate_overlay_fraction.back() << std::endl;
1239  std::cout << "===============================================================" << std::endl;
1240  }
1242 
1243 
1244  } //end of cluster loop
1245 
1246  // Plot the event
1249  }
1250 
1251  //group clusters HERE
1252  std::pair<int, std::pair<std::vector<std::vector<double>>, std::vector<double>> > group_result = GroupClusterCandidate(vars.m_trackstub_num_candidates, vars.m_trackstub_candidate_plane, vars.m_trackstub_candidate_max_tick, vars.m_trackstub_candidate_min_tick);
1253  vars.m_trackstub_num_candidate_groups = group_result.first;
1254  vars.m_grouped_trackstub_candidate_indices = group_result.second.first;
1255  vars.m_trackstub_candidate_group_timeoverlap_fraction = group_result.second.second;
1256  }
1257 
1258  // --------------- shower clustering --------------------------
1259  std::cout << "------------- Shower clustering --------------------" << std::endl;
1260  std::cout << "SEAview Shower cluster formation: " << (paras.s_runSEAview ? "true" : "false" ) << " nshower requirement: " << paras.s_SEAviewNumRecoShower << ", actual num shower: " << showers.size() << " | ntrack requirement: " << paras.s_SEAviewNumRecoTrack << ", actual num track: " << tracks.size() << std::endl;
1261 
1262  if(!paras.s_run_pi0_filter &&
1263  paras.s_runSEAview &&
1264  (paras.s_SEAviewNumRecoShower == -1 || (int)showers.size()== paras.s_SEAviewNumRecoShower) &&
1265  (paras.s_SEAviewNumRecoTrack == -1 || (int)tracks.size() == paras.s_SEAviewNumRecoTrack) ){
1266 
1267  art::Ptr<recob::Shower> p_shr = showers.front();
1268 
1269  PandoraPFParticle* ppfp = PPFP_GetPPFPFromShower(allPPFParticles, p_shr);
1270  art::Ptr<recob::PFParticle> p_pfp = ppfp->pPFParticle;//showerToNuPFParticleMap[p_shr];
1271  std::vector<art::Ptr<recob::Hit>> p_hits = pfParticleToHitsMap[p_pfp];
1272 
1273 
1274  int p_sliceid = ppfp->get_SliceID();// PFPToSliceIdMap[p_pfp];
1275  auto p_slice_hits = sliceIDToHitsMap[p_sliceid];
1276 
1277  std::string uniq_tag = "HitThres_"+ std::to_string(static_cast<int>(paras.s_SEAviewHitThreshold)) + "_" + std::to_string(vars.m_run_number)+"_"+std::to_string(vars.m_subrun_number)+"_"+std::to_string(vars.m_event_number);
1278 
1279  //Setup seaviewr object
1280  seaview::SEAviewer sevd("Shower_"+uniq_tag, paras.s_geom, theDetector );
1281  //Pass in any bad channels you like
1282  // sevd.setBadChannelList(bad_channel_list_fixed_mcc9);
1283  //Give it a vertex to center around
1285 
1286  //Add hits to consider for clustering
1287  //sevd.addHitsToConsider(hitVector);// DONT do this yet, need something smarter for SSV
1288  //sevd.filterConsideredHits(150);
1289  sevd.addHitsToConsider(p_slice_hits);
1290 
1291  //Add all hits in the events
1292  sevd.addAllHits(hitVector); // std::vector<art::Ptr<recob::Hit>>
1294 
1295  //Add all the "nice " PFParticle Hits, as well as what to label
1296  //sevd.addPFParticleHits(p_hits, "Shower"); //std::vector<art::Ptr<recob::Hit>> and std::string
1297  sevd.addPFParticleHits(p_hits, "Shower", vars.m_reco_shower_energy_max[0], vars.m_reco_shower_conversion_distance[0], vars.m_reco_shower_impact_parameter[0]); //std::vector<art::Ptr<recob::Hit>> and std::string
1298 
1299  //and add the SingleShower we like
1300  sevd.addShower(p_shr); // art::Ptr<recob::Shower>
1301 
1302  //Add all track PFP
1303  int i_trk = 0;
1304  for(auto &trk: tracks){
1305  art::Ptr<recob::PFParticle> p_pfp_trk = trackToNuPFParticleMap[trk];
1306  std::vector<art::Ptr<recob::Hit>> p_hits_trk = pfParticleToHitsMap[p_pfp_trk];
1307  //sevd.addPFParticleHits(p_hits_trk,"track");
1308  sevd.addPFParticleHits(p_hits_trk,"track", vars.m_reco_track_length[i_trk], vars.m_reco_track_spacepoint_principal0[i_trk]);
1309  sevd.addTrack(trk);
1310  ++i_trk;
1311  }
1312 
1313  //Add all cosmic-relatd PFP // DONT do this yet, see line 1206
1314  /*for(auto &cr: crParticles){
1315  std::vector<art::Ptr<recob::Hit>> p_hits_cr = cr_pfParticleToHitsMap[cr];
1316  sevd.addPFParticleHits(p_hits_cr,"cosmic");
1317  }
1318  */
1319 
1320  //We then calculate Unassociated hits, i.e the hits not associated to the "Shower" or tracksyou passed in.
1321  auto vnh= sevd.calcUnassociatedHits();
1322  vars.m_sss_num_unassociated_hits =vnh[1]+vnh[2];
1324  vars.m_sss_num_associated_hits = vnh[0]-vnh[1]-vnh[2];
1325 
1326  //Recluster, group unassociated hits into different clusters
1328 
1329 
1330  //This is the place I will put the new Second Shower Search
1331  std::vector<seaview::cluster> vec_SEAclusters ;
1332  sevd.analyzeShowerLikeClusters(paras.s_SEAviewDbscanEps, showerToNuPFParticleMap, pfParticleToHitsMap, vec_SEAclusters);
1333 
1334  //And save to file.
1335  std::cout<<"After SEAview we have "<<vec_SEAclusters.size()<<" Shower clusters to chat about"<<std::endl;
1336 
1338  for(size_t c=0; c< vec_SEAclusters.size(); c++){
1339  auto clu = vec_SEAclusters.at(c); //type: seaview::cluster
1340  int pl = clu.getPlane();
1341  auto hitz = clu.getHits();
1342  double Ep = CalcEShowerPlane(hitz,pl, paras);
1343  int remerge = clu.getShowerRemerge();
1344  seaview::cluster_score * ssscorz = clu.getScore();
1345 
1346  std::cout<<c<<" "<<pl<<" "<<Ep<<" "<<clu.getImpactParam()<<" "<<clu.getFitSlope()<<" "<<clu.getFitCons()<<" "<<clu.getMeanADC() << " " << clu.getADCrms() << " "<<clu.getAngleWRTShower()<<" "<<remerge<<std::endl;
1347 
1348  //if the cluster is too close to the recob::shower, then do not include it
1349  if(remerge>=0 && remerge< (int)vars.m_reco_shower_reclustered_energy_plane2.size()){
1350  if(pl==0)vars.m_reco_shower_reclustered_energy_plane0[remerge]+=Ep;
1351  if(pl==1)vars.m_reco_shower_reclustered_energy_plane1[remerge]+=Ep;
1352  if(pl==2)vars.m_reco_shower_reclustered_energy_plane2[remerge]+=Ep;
1353 
1354  continue;// Dont include this as a viable cluster!
1355  }
1356 
1358 
1359  //determine if this cluster is in neutrino slice
1360  vars.m_sss_candidate_in_nu_slice.push_back(clu.InNuSlice(sliceIDToHitsMap, p_sliceid));
1361 
1362  //Fill All the bits
1363  vars.m_sss_candidate_num_hits.push_back((int)hitz.size());
1364  vars.m_sss_candidate_num_wires.push_back((int)ssscorz->n_wires);
1365  vars.m_sss_candidate_num_ticks.push_back((int)ssscorz->n_ticks);
1366  vars.m_sss_candidate_plane.push_back(pl);
1367  vars.m_sss_candidate_PCA.push_back(ssscorz->pca_0);
1368  vars.m_sss_candidate_impact_parameter.push_back(clu.getImpactParam());
1369  vars.m_sss_candidate_fit_slope.push_back(clu.getFitSlope());
1370  vars.m_sss_candidate_fit_constant.push_back(clu.getFitCons());
1371  vars.m_sss_candidate_mean_tick.push_back(ssscorz->mean_tick);
1372  vars.m_sss_candidate_max_tick.push_back(ssscorz->max_tick);
1373  vars.m_sss_candidate_min_tick.push_back(ssscorz->min_tick);
1374  vars.m_sss_candidate_min_wire.push_back(ssscorz->min_wire);
1375  vars.m_sss_candidate_max_wire.push_back(ssscorz->max_wire);
1376  vars.m_sss_candidate_mean_wire.push_back(ssscorz->mean_wire);
1377  vars.m_sss_candidate_min_dist.push_back(ssscorz->min_dist);
1378  vars.m_sss_candidate_wire_tick_based_length.push_back(clu.getWireTickBasedLength());
1379  vars.m_sss_candidate_mean_ADC.push_back(clu.getMeanADC());
1380  vars.m_sss_candidate_ADC_RMS.push_back(clu.getADCrms());
1381  vars.m_sss_candidate_energy.push_back(Ep);
1382  vars.m_sss_candidate_angle_to_shower.push_back(clu.getAngleWRTShower());
1383  vars.m_sss_candidate_remerge.push_back(remerge);
1384 
1385 
1386  //MCTruth matching for pi0's
1387  if(paras.s_is_data){
1388  vars.m_sss_candidate_matched.push_back(-1);
1389  vars.m_sss_candidate_pdg.push_back(-1);
1390  vars.m_sss_candidate_parent_pdg.push_back(-1);
1391  vars.m_sss_candidate_trackid.push_back(-1);
1392  vars.m_sss_candidate_true_energy.push_back(-1);
1393  vars.m_sss_candidate_overlay_fraction.push_back(-1);
1395  }else{
1396 
1397  auto ssmatched = SecondShowerMatching(hitz, mcparticles_per_hit, mcParticleVector, MCParticleToTrackIdMap, vars);
1398  vars.m_sss_candidate_matched.push_back(ssmatched[0]);
1399  vars.m_sss_candidate_pdg.push_back(ssmatched[1]);
1400  vars.m_sss_candidate_parent_pdg.push_back(ssmatched[2]);
1401  vars.m_sss_candidate_trackid.push_back(ssmatched[3]);
1402  vars.m_sss_candidate_true_energy.push_back(ssmatched[4]);
1403  vars.m_sss_candidate_overlay_fraction.push_back(ssmatched[5]);
1405 
1406  //Guanqun: print out (best-matched) truth information of the cluster
1407  std::cout << "Cluster: " << vars.m_sss_num_candidates-1 << " plane: " << vars.m_sss_candidate_plane.back() << ", energy: " << vars.m_sss_candidate_energy.back() << "\n";
1408  std::cout << "Cluster is matched: " << vars.m_sss_candidate_matched.back() << ", matched PDG: " << vars.m_sss_candidate_pdg.back() << " track ID: " << vars.m_sss_candidate_trackid.back() << " overlay fraction: " << vars.m_sss_candidate_overlay_fraction.back() << std::endl;
1409  std::cout << "===============================================================" << std::endl;
1410  }
1411 
1412 
1414 
1415  } //end of cluster loop
1416 
1417  // Plot the event
1418  if(paras.s_SEAviewMakePDF){
1420  }
1421 
1422  }
1423 
1424  for(int i =0; i<(int)showers.size(); i++){
1426  }
1427 
1428  // ################################################### END SEAview END SEAview #########################################################
1429  // #####################################################################################################################################
1430 
1431 
1432 
1433  // PandoraAllOutComes
1434  // I.e This runs over all 3D reco showers in the whole event and find second shower candidates
1435  if(!paras.s_run_pi0_filter){
1436  std::cout<<"------------ Shower3D --------------"<<std::endl;
1437  /*for(auto &s : showers){
1438  std::cout<<"shower pfp key : "<<showerToNuPFParticleMap[s].key()<<" self: "<<showerToNuPFParticleMap[s]->Self()<<std::endl;
1439  }
1440  for(auto &s : tracks){
1441  std::cout<<"track pfp key : "<<trackToNuPFParticleMap[s].key()<<" self: "<<trackToNuPFParticleMap[s]->Self()<<std::endl;
1442  }*/
1443 
1444  SecondShowerSearch3D(showers, showerToNuPFParticleMap, tracks,trackToNuPFParticleMap,evt, vars, paras);
1445 
1446 
1447  //And cluster the 2d and 3d second showers. Very simple TODO
1449 
1450  }
1451 
1452 
1453 
1454  //---------------------- END OF LOOP, fill vertex ---------------------
1455  bool filter_pass_2g1p = Pi0PreselectionFilter(vars,paras);
1456  bool filter_pass_2g0p = Pi0PreselectionFilter2g0p(vars,paras);
1457 
1458  if (paras.s_fill_trees && ( (filter_pass_2g1p && paras.s_run_pi0_filter_2g1p) || (filter_pass_2g0p && paras.s_run_pi0_filter_2g0p) || !paras.s_run_pi0_filter ) ) {
1459  vars.vertex_tree->Fill();
1460  vars.ncdelta_slice_tree->Fill();
1461  vars.eventweight_tree->Fill();
1462  vars.true_eventweight_tree->Fill();
1463  vars.geant4_tree->Fill();
1464  }
1465 
1466  //Rest the vertex after filling
1467  // this->ClearMeta();
1468  if(paras.s_run_pi0_filter_2g1p) return filter_pass_2g1p;
1469  else if(paras.s_run_pi0_filter_2g0p) return filter_pass_2g0p;
1470 
1471  //if not in filter mode pass all
1472  return true;
1473 
1474  }// end filter_module main class
1475 
1476 
1477 
1478  //-------------------------------------------------------------------------------------------
1479 
1480  //This runs ONCE at the start of the job and sets up all the necessary services and TTrees
1482  {
1483  mf::LogDebug("SinglePhoton") << " *** beginJob() *** " << "\n";
1484 
1485  art::ServiceHandle<art::TFileService> tfs;//output ROOT
1486 
1487  vars.run_subrun_tree = tfs->make<TTree>("run_subrun_tree","run_subrun_tree");
1488  vars.true_eventweight_tree = tfs->make<TTree>("true_eventweight_tree", "true_eventweight_tree");
1489  vars.pot_tree = tfs->make<TTree>("pot_tree", "pot_tree");
1490 
1491  vars.eventweight_tree = tfs->make<TTree>("eventweight_tree", "eventweight_tree");
1492  vars.ncdelta_slice_tree = tfs->make<TTree>("ncdelta_slice_tree", "ncdelta_slice_tree");
1493  vars.geant4_tree = tfs->make<TTree>("geant4_tree","geant4_tree");
1494  vars.vertex_tree = tfs->make<TTree>("vertex_tree", "vertex_tree");
1495 
1496  //run_subrun_tree, reset some POT
1497  vars.m_run = 0;
1498  vars.m_subrun = 0;
1499  vars.m_subrun_pot = 0;
1500 
1501  // --------------------- POT Releated variables -----------------
1504  vars.m_pot_count=0;
1505  vars.m_pot_per_event = 0;
1506  vars.m_pot_per_subrun = 0;
1508 
1509 
1519 
1523 
1524 
1525  //------------------- List of Selected Events to run --------
1527  std::cout << "SinglePhoton \t||\t Running in selected-event only mode " << std::endl;
1528 
1529  std::ifstream infile(paras.s_selected_event_list);
1530  if(!infile){
1531  std::cerr << "Fail to open file: " <<paras.s_selected_event_list << std::endl;
1532  return;
1533  }
1534 
1535  //read from file, run number, subrun number ,event number that should be run
1536  vars.m_selected_set.clear();
1537  std::string line;
1538  while(std::getline(infile, line)){
1539  std::istringstream ss(line);
1540 
1541  std::vector<int> event_info;
1542  for(int i; ss >> i; ) event_info.push_back(i);
1543 
1544  vars.m_selected_set.insert(event_info);
1545  }
1546 
1547  infile.close();
1548 
1549  if(g_is_verbose){
1550  std::cout << "Selected Events: " << std::endl;
1551  std::cout << "Run \t SubRun \t Event" << std::endl;
1552  for(auto & v: vars.m_selected_set){
1553  std::for_each(v.begin(), v.end(), [](int n){std::cout << n<<" \t "; });
1554  std::cout << std::endl;
1555  }
1556  }
1557  }
1558  }
1559 
1560 
1561 
1562  bool SinglePhoton::beginSubRun(art::SubRun& sr) {
1563 
1564  vars.m_run = sr.run();
1565  vars.m_subrun = sr.subRun();
1566 
1567  double this_pot = 0;
1568 
1569  //reset subrun count
1570  vars.m_subrun_counts = 0;
1571 
1572 
1573  if(paras.s_potLabel != ""){
1574  if(paras.s_potLabel == "generator"){
1575 
1576  art::Handle<sumdata::POTSummary> gen_pot_hand;
1577  if(sr.getByLabel(paras.s_potLabel,gen_pot_hand)){
1578  this_pot = gen_pot_hand->totgoodpot;
1579  vars.m_pot_count += this_pot;
1580  std::cout<<"SinglePhoton::beginSubRun()\t||\t SubRun POT: "<<this_pot<<" . Current total POT this file: "<<vars.m_pot_count<<" (label) "<<paras.s_potLabel<<std::endl;
1581  }
1582  }else{
1583 
1584  art::Handle<sumdata::POTSummary> potSummaryHandlebnbETOR875;
1585  if (sr.getByLabel("beamdata","bnbETOR875",potSummaryHandlebnbETOR875)){
1586  this_pot =potSummaryHandlebnbETOR875->totpot;
1587  vars.m_pot_count += this_pot;
1588  std::cout<<"SinglePhoton::beginSubRun()\t||\t SubRun POT: "<<potSummaryHandlebnbETOR875->totpot<<" . Current total POT this file: "<<vars.m_pot_count<<" (label) "<<paras.s_potLabel<<std::endl;
1589  }
1590  }
1591  }
1592  vars.m_subrun_pot = this_pot;
1593  return true;
1594  }
1595 
1596 
1597  bool SinglePhoton::endSubRun(art::SubRun&sr){
1598 
1599  vars.run_subrun_tree->Fill();
1600  return true;
1601  }
1602 
1604  {
1605  // if (vars.m_print_out_event){
1606  // out_stream.close();
1607  // }
1608  vars.pot_tree->Fill();
1609  }
1610 
1611  DEFINE_ART_MODULE(SinglePhoton)
1612 } //namespace single_photon
std::vector< double > m_sss_candidate_PCA
Definition: variables.h:207
std::vector< double > m_sss_candidate_overlay_fraction
Definition: variables.h:232
void AnalyzeTrackCalo(const std::vector< art::Ptr< recob::Track >> &tracks, std::vector< PandoraPFParticle > all_PPFPs, var_all &vars, para_all &paras)
TTree * ncdelta_slice_tree
Definition: variables.h:322
void CreateShowerBranches(var_all &vars)
std::vector< double > m_trackstub_candidate_PCA
Definition: variables.h:362
std::vector< double > SecondShowerMatching(std::vector< art::Ptr< recob::Hit >> &hitz, art::FindManyP< simb::MCParticle, anab::BackTrackerHitMatchingData > &mcparticles_per_hit, std::vector< art::Ptr< simb::MCParticle >> &mcParticleVector, std::map< int, art::Ptr< simb::MCParticle >> &MCParticleToTrackIdMap, var_all &vars)
std::string s_hitMCParticleAssnsLabel
Definition: variables.h:106
int loadVertex(double m_vertex_pos_x, double m_vertex_pos_y, double m_vertex_pos_z)
Definition: SEAviewer.cxx:350
std::vector< double > m_trackstub_candidate_mean_ADC_first_to_second_ratio
Definition: variables.h:380
std::vector< double > m_trackstub_candidate_max_wire
Definition: variables.h:370
std::vector< double > m_sss_candidate_wire_tick_based_length
Definition: variables.h:221
bool beginSubRun(art::SubRun &sr) override
: grab run, subrun number, and subrun POT, fill the TTree
bool endSubRun(art::SubRun &sr) override
double s_SEAviewPlotDistance
Definition: variables.h:76
std::vector< double > m_reco_track_spacepoint_principal0
Definition: variables.h:608
std::string s_showerLabel
The label for the shower producer from PFParticles.
Definition: variables.h:100
std::vector< double > m_trackstub_candidate_linear_fit_chi2
Definition: variables.h:382
std::vector< int > m_trackstub_candidate_in_nu_slice
Definition: variables.h:357
void CreateEventWeightBranches(var_all &vars)
PandoraPFParticle * PPFP_GetPPFPFromPFID(std::vector< PandoraPFParticle > &PPFPs, int id)
void CollectMCParticles(const art::Event &evt, const std::string &label, std::map< art::Ptr< simb::MCTruth >, std::vector< art::Ptr< simb::MCParticle >>> &truthToParticles, std::map< art::Ptr< simb::MCParticle >, art::Ptr< simb::MCTruth >> &particlesToTruth, std::map< int, art::Ptr< simb::MCParticle > > &MCParticleToTrackIdMap, var_all &vars)
Definition: Processors.cxx:96
void CreateFlashBranches(var_all &vars)
Utilities related to art service access.
void ClearFlashes(var_all &vars)
std::vector< double > m_sss_candidate_ADC_RMS
Definition: variables.h:209
std::set< std::vector< int > > m_selected_set
Definition: variables.h:314
void CreateSecondShowerBranches3D(var_all &vars)
void set_ParticleID(const art::Ptr< anab::ParticleID > input_ParticleID)
void ClearMCTruths(var_all &vars)
std::vector< double > m_trackstub_candidate_track_angle_wrt_shower_direction
Definition: variables.h:381
ClusterModuleLabel join with tracks
std::vector< double > m_sss_candidate_fit_constant
Definition: variables.h:213
int addPFParticleHits(std::vector< art::Ptr< recob::Hit >> &hits, std::string leg, double arg1=0.0, double arg2=0.0, double arg3=0.0)
Definition: SEAviewer.cxx:279
void AnalyzeEventWeight(art::Event const &e, var_all &vars)
Definition: analyze_MC.cxx:60
std::vector< int > m_sss_candidate_remerge
Definition: variables.h:225
std::vector< double > m_trackstub_candidate_mean_tick
Definition: variables.h:366
std::vector< int > m_sss_candidate_plane
Definition: variables.h:206
std::vector< double > m_trackstub_candidate_mean_ADC_first_half
Definition: variables.h:378
process_name opflash particleana ie x
BEGIN_PROLOG could also be cerr
process_name cluster
Definition: cheaterreco.fcl:51
std::vector< double > m_sss_candidate_impact_parameter
Definition: variables.h:210
std::vector< double > trackRecoMCmatching(std::vector< art::Ptr< recob::Track >> &objectVector, std::map< art::Ptr< recob::Track >, art::Ptr< simb::MCParticle >> &objectToMCParticleMap, std::map< art::Ptr< recob::Track >, art::Ptr< recob::PFParticle >> &objectToPFParticleMap, std::map< art::Ptr< recob::PFParticle >, std::vector< art::Ptr< recob::Hit >> > &pfParticleToHitsMap, art::FindManyP< simb::MCParticle, anab::BackTrackerHitMatchingData > &mcparticles_per_hit, std::vector< art::Ptr< simb::MCParticle >> &mcParticleVector, var_all &vars)
std::vector< int > m_sss_candidate_num_wires
Definition: variables.h:204
std::vector< double > m_reco_shower_impact_parameter
Definition: variables.h:801
std::vector< double > m_reco_shower_end_dist_to_SCB
Definition: variables.h:789
std::vector< double > m_trackstub_candidate_group_timeoverlap_fraction
Definition: variables.h:396
std::vector< int > m_trackstub_candidate_remerge
Definition: variables.h:384
void ClearGeant4Branches(var_all &vars)
: fill event weight related variables
std::vector< double > analyzeShowerLikeClusters(double eps, const std::map< art::Ptr< recob::Shower >, art::Ptr< recob::PFParticle >> &showerToPFParticleMap, const std::map< art::Ptr< recob::PFParticle >, std::vector< art::Ptr< recob::Hit >> > &pfParticleToHitsMap, std::vector< seaview::cluster > &vec_c)
Definition: SEAviewer.cxx:783
std::vector< int > m_trackstub_candidate_plane
Definition: variables.h:361
spacecharge::SpaceCharge const * s_SCE
Definition: variables.h:142
void CreateTrackBranches(var_all &vars)
std::string s_caloLabel
The label for calorimetry associations producer.
Definition: variables.h:101
int DefineNuSlice(std::vector< PandoraPFParticle > &PPFPs)
Declaration of signal hit object.
int m_trackstub_num_unassociated_hits
Definition: variables.h:351
std::vector< double > m_sss_candidate_angle_to_shower
Definition: variables.h:223
void BuildMCParticleHitMaps(const art::Event &evt, const std::string &label, const std::vector< art::Ptr< recob::Hit >> &hitVector, std::map< art::Ptr< simb::MCParticle >, std::vector< art::Ptr< recob::Hit > > > &particlesToHits, std::map< art::Ptr< recob::Hit >, art::Ptr< simb::MCParticle > > &hitsToParticles, const lar_pandora::LArPandoraHelper::DaughterMode daughterMode, std::map< int, art::Ptr< simb::MCParticle > > &MCParticleToTrackIdMap, var_all &vars)
Definition: Processors.cxx:160
void CreateSecondShowerBranches(var_all &vars)
void ClearSecondShowers(var_all &vars)
std::vector< int > m_trackstub_candidate_pdg
Definition: variables.h:387
int setHitThreshold(double)
Definition: SEAviewer.cxx:195
double s_SEAviewStubDbscanEps
Definition: variables.h:88
double s_track_calo_trunc_fraction
Definition: variables.h:157
std::vector< int > m_sss_candidate_trackid
Definition: variables.h:230
std::vector< int > m_sss_candidate_num_hits
Definition: variables.h:203
void CreateGeant4Branches(var_all &vars)
std::vector< int > m_trackstub_candidate_parent_pdg
Definition: variables.h:388
art::Ptr< recob::PFParticle > pPFParticle
std::string s_pandoraLabel
The label for the pandora producer.
Definition: variables.h:98
std::pair< int, std::pair< std::vector< std::vector< double > >, std::vector< double > > > GroupClusterCandidate(int num_clusters, const std::vector< int > &cluster_planes, const std::vector< double > &cluster_max_ticks, const std::vector< double > &cluster_min_ticks)
std::size_t size(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:561
int addTrack(art::Ptr< recob::Track > &trk)
Definition: SEAviewer.cxx:320
Contains data associated to particles from detector simulation.
int m_trackstub_num_candidate_groups
Definition: variables.h:394
process_name hit
Definition: cheaterreco.fcl:51
void AnalyzeGeant4(const std::vector< art::Ptr< simb::MCParticle >> &mcParticleVector, var_all &vars)
Definition: analyze_MC.cxx:11
void RecoMCTracks(std::vector< PandoraPFParticle > all_PPFPs, const std::vector< art::Ptr< recob::Track >> &tracks, std::map< art::Ptr< recob::Track >, art::Ptr< simb::MCParticle > > &trackToMCParticleMap, std::map< art::Ptr< simb::MCParticle >, art::Ptr< simb::MCTruth >> &MCParticleToMCTruthMap, std::vector< art::Ptr< simb::MCParticle >> &mcParticleVector, std::map< int, art::Ptr< simb::MCParticle > > &MCParticleToTrackIdMap, std::vector< double > &vfrac, var_all &vars)
std::string s_hitfinderLabel
Definition: variables.h:105
std::string s_potLabel
Definition: variables.h:107
int m_trackstub_unassociated_hits_below_threshold
Definition: variables.h:352
std::vector< int > m_sss_candidate_pdg
Definition: variables.h:228
void showerRecoMCmatching(std::vector< PandoraPFParticle > all_PPFPs, std::vector< art::Ptr< recob::Shower >> &showerVector, std::map< art::Ptr< recob::Shower >, art::Ptr< simb::MCParticle >> &showerToMCParticleMap, art::FindManyP< simb::MCParticle, anab::BackTrackerHitMatchingData > &mcparticles_per_hit, std::vector< art::Ptr< simb::MCParticle >> &mcParticleVector, std::map< int, art::Ptr< simb::MCParticle > > &MCParticleToTrackIdMap, var_all &vars)
std::string s_CRTHitProducer
Definition: variables.h:113
std::vector< double > m_trackstub_candidate_min_dist
Definition: variables.h:372
void beginJob() override
Begin the job, setting up !
std::vector< double > m_reco_shower_reclustered_energy_plane0
Definition: variables.h:993
void PPFP_FindAncestor(std::vector< PandoraPFParticle > &PPFPs)
void CreateStubBranches(var_all &vars)
void Save_EventMeta(art::Event &evt, var_all &vars)
std::vector< double > m_trackstub_candidate_ADC_RMS
Definition: variables.h:364
std::string s_flashLabel
Definition: variables.h:102
std::vector< double > m_trackstub_candidate_min_tick
Definition: variables.h:368
void ClearEventWeightBranches(var_all &vars)
std::vector< double > s_gain_data
Definition: variables.h:130
bool filter(art::Event &evt) override
Analyze an event!
std::vector< double > s_gain_mc
Definition: variables.h:129
void setTPCGeom(para_all &paras)
bool Pi0PreselectionFilter2g0p(var_all &vars, para_all &paras)
Definition: Processors.cxx:203
std::vector< int > m_trackstub_candidate_trackid
Definition: variables.h:389
double s_exiting_photon_energy_threshold
Definition: variables.h:160
int m_sss_num_unassociated_hits_below_threshold
Definition: variables.h:194
std::string s_trackLabel
The label for the track producer from PFParticles.
Definition: variables.h:99
std::vector< double > m_reco_shower_reclustered_energy_plane2
Definition: variables.h:996
std::vector< double > m_trackstub_candidate_mean_ADC_second_half
Definition: variables.h:379
std::vector< double > m_sss_candidate_true_energy
Definition: variables.h:231
double s_track_calo_min_dEdx_hits
Definition: variables.h:156
SinglePhoton(fhicl::ParameterSet const &pset)
Constructor.
bool Pi0PreselectionFilter(var_all &vars, para_all &paras)
Definition: Processors.cxx:185
std::vector< double > m_trackstub_candidate_matched_energy_fraction_best_plane
Definition: variables.h:386
std::vector< int > Printer_header(std::vector< std::string > headings)
std::vector< double > m_reco_shower_spacepoint_x
Definition: variables.h:1009
std::string s_generatorLabel
Definition: variables.h:108
std::vector< double > m_reco_shower_energy_max
Definition: variables.h:988
void IsolationStudy(std::vector< PandoraPFParticle > all_PPFPs, const std::vector< art::Ptr< recob::Track >> &tracks, const std::vector< art::Ptr< recob::Shower >> &showers, detinfo::DetectorPropertiesData const &theDetector, var_all &vars, para_all &paras)
std::vector< double > m_reco_shower_end_dist_to_active_TPC
Definition: variables.h:788
void CreateMCTruthBranches(var_all &vars)
std::vector< double > m_sss_candidate_energy
Definition: variables.h:222
std::vector< double > m_reco_shower_reclustered_energy_max
Definition: variables.h:992
std::string s_CRTVetoLabel
Definition: variables.h:111
std::string s_true_eventweight_label
Definition: variables.h:114
std::string s_selected_event_list
Definition: variables.h:94
std::vector< double > m_trackstub_candidate_overlay_fraction
Definition: variables.h:391
geo::GeometryCore const * s_geom
Definition: variables.h:143
std::string s_shower3dLabel
Definition: variables.h:95
void CreateSliceBranches(var_all &vars)
std::vector< double > m_sss_candidate_mean_tick
Definition: variables.h:214
std::vector< double > m_trackstub_candidate_max_tick
Definition: variables.h:367
std::vector< double > m_sss_candidate_max_wire
Definition: variables.h:218
std::vector< double > m_sss_candidate_min_dist
Definition: variables.h:220
std::vector< double > analyzeTrackLikeClusters(double eps, const std::map< art::Ptr< recob::Shower >, art::Ptr< recob::PFParticle >> &showerToPFParticleMap, const std::map< art::Ptr< recob::PFParticle >, std::vector< art::Ptr< recob::Hit >> > &pfParticleToHitsMap, std::vector< seaview::cluster > &vec_c)
Definition: SEAviewer.cxx:737
int distToSCB(double &dist, std::vector< double > &vec, para_all &paras)
std::vector< int > m_trackstub_candidate_num_ticks
Definition: variables.h:360
std::vector< std::vector< double > > m_grouped_trackstub_candidate_indices
Definition: variables.h:395
void ClearSecondShowers3D(var_all &vars)
Declaration of cluster object.
Class def header for mctrack data container.
int Print(double plot_distance)
Definition: SEAviewer.cxx:405
PandoraPFParticle * PPFP_GetPPFPFromTrack(std::vector< PandoraPFParticle > &PPFPs, art::Ptr< recob::Track > pTrack)
void AnalyzeTracks(std::vector< PandoraPFParticle > all_PPFPs, const std::vector< art::Ptr< recob::Track >> &tracks, std::map< art::Ptr< recob::PFParticle >, std::vector< art::Ptr< recob::SpacePoint >>> &pfParticleToSpacePointsMap, std::map< int, art::Ptr< simb::MCParticle > > &MCParticleToTrackIdMap, std::map< int, double > &sliceIdToNuScoreMap, var_all &vars, para_all &paras)
void set_Calorimetries(const std::vector< art::Ptr< anab::Calorimetry >> input_Calorimetries)
TTree * true_eventweight_tree
Definition: variables.h:326
std::vector< double > m_sss_candidate_matched_energy_fraction_best_plane
Definition: variables.h:227
std::vector< double > m_trackstub_candidate_wire_tick_based_length
Definition: variables.h:377
Provides recob::Track data product.
int addShower(art::Ptr< recob::Shower > &shr)
Definition: SEAviewer.cxx:312
double s_exiting_proton_energy_threshold
Definition: variables.h:161
int addAllHits(std::vector< art::Ptr< recob::Hit >> &hits)
Definition: SEAviewer.cxx:200
void PPFP_FindSliceIDandHits(std::vector< PandoraPFParticle > &PPFPs, art::Ptr< recob::Slice > slice, const std::vector< art::Ptr< recob::PFParticle > > PFP_in_slice, const std::vector< art::Ptr< recob::Hit > > Hit_inslice)
std::string s_Spline_CV_label
Definition: variables.h:116
void ResizeShowers(size_t size, var_all &vars)
std::vector< double > m_reco_shower_spacepoint_z
Definition: variables.h:1010
std::vector< double > m_sss_candidate_fit_slope
Definition: variables.h:211
void reconfigure(fhicl::ParameterSet const &pset)
Configure memeber variables using FHiCL parameters.
double CalcEShowerPlane(const std::vector< art::Ptr< recob::Hit >> &hits, int this_plane, para_all &paras)
Definition: Processors.cxx:293
void endJob() override
End the job, setting down !
void ClearMeta(var_all &vars)
: reset/clear data members
std::vector< int > m_sss_candidate_parent_pdg
Definition: variables.h:229
std::vector< double > m_trackstub_candidate_ioc_based_length
Definition: variables.h:376
void ClearStubs(var_all &vars)
double s_SEAviewStubDbscanMinPts
Definition: variables.h:87
std::vector< double > m_trackstub_candidate_min_impact_parameter_to_shower
Definition: variables.h:373
void ResizeMCTruths(size_t size, var_all &vars)
std::string s_CRTTzeroLabel
Definition: variables.h:112
void Printer_content(std::vector< std::string > nums, std::vector< int > spacers)
std::vector< double > m_trackstub_candidate_true_energy
Definition: variables.h:390
std::string to_string(WindowPattern const &pattern)
void ClearIsolation(var_all &vars)
then echo File list $list not found else cat $list while read file do echo $file sed s
Definition: file_to_url.sh:60
void Save_PFParticleInfo(std::vector< PandoraPFParticle > PPFPs, var_all &vars, para_all &paras)
double s_SEAviewStubPlotDistance
Definition: variables.h:86
PandoraPFParticle * PPFP_GetPPFPFromShower(std::vector< PandoraPFParticle > &PPFPs, art::Ptr< recob::Shower > pShower)
void SetClusterLegend(int cluster, double energy, int is_matched, int matched_pdg, double overlay_fraction)
Definition: SEAviewer.cxx:1262
double s_SEAviewMaxPtsLinFit
Definition: variables.h:80
double s_SEAviewStubHitThreshold
Definition: variables.h:85
object containing MC truth information necessary for making RawDigits and doing back tracking ...
int runseaDBSCAN(double min_pts, double eps)
Definition: SEAviewer.cxx:687
Class def header for MCShower data container.
int trigger_offset(DetectorClocksData const &data)
std::vector< double > m_trackstub_candidate_min_ioc_to_shower_start
Definition: variables.h:375
bool IsEventInList(int run, int subrun, int event, var_all &vars)
Definition: Processors.cxx:219
void SimpleSecondShowerCluster(var_all &vars, para_all &paras)
double s_SEAviewDbscanMinPts
Definition: variables.h:78
void set_HasPID(const bool input_bool)
double m_photonu_weight_high
Definition: variables.h:568
std::vector< double > m_reco_shower_spacepoint_y
Definition: variables.h:1011
std::vector< int > m_sss_candidate_matched
Definition: variables.h:226
art::ServiceHandle< art::TFileService > tfs
std::vector< int > m_trackstub_candidate_num_wires
Definition: variables.h:359
int addHitsToConsider(std::vector< art::Ptr< recob::Hit >> &hits)
Definition: SEAviewer.cxx:164
void CollectPID(std::vector< art::Ptr< recob::Track >> &tracks, std::vector< PandoraPFParticle > all_PPFPs, var_all &vars)
void ResizeSlices(size_t size, var_all &vars)
std::vector< double > m_trackstub_candidate_min_wire
Definition: variables.h:369
double distToTPCActive(std::vector< double > &vec, para_all &paras)
TCEvent evt
Definition: DataStructs.cxx:8
std::vector< int > m_sss_candidate_in_nu_slice
Definition: variables.h:202
std::vector< double > m_reco_shower_reclustered_energy_plane1
Definition: variables.h:995
std::string s_pidLabel
For PID stuff.
Definition: variables.h:110
void AnalyzeShowers(std::vector< PandoraPFParticle > all_PPFPs, const std::vector< art::Ptr< recob::Shower >> &showers, std::map< art::Ptr< recob::Cluster >, std::vector< art::Ptr< recob::Hit >> > &clusterToHitMap, double triggeroffset, detinfo::DetectorPropertiesData const &theDetector, var_all &vars, para_all &paras)
void ResizeFlashes(size_t size, var_all &vars)
std::string s_truthmatching_signaldef
Definition: variables.h:117
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
std::vector< double > m_trackstub_candidate_energy
Definition: variables.h:383
std::vector< int > calcUnassociatedHits()
Definition: SEAviewer.cxx:227
std::vector< double > m_trackstub_candidate_mean_wire
Definition: variables.h:371
std::vector< double > m_trackstub_candidate_min_conversion_dist_to_shower_start
Definition: variables.h:374
void SecondShowerSearch3D(std::vector< art::Ptr< recob::Shower >> &showers, std::map< art::Ptr< recob::Shower >, art::Ptr< recob::PFParticle >> &NormalShowerToPFParticleMap, std::vector< art::Ptr< recob::Track >> &tracks, std::map< art::Ptr< recob::Track >, art::Ptr< recob::PFParticle >> &NormalTrackToPFParticleMap, art::Event const &evt, var_all &vars, para_all &paras)
void CreateIsolationBranches(var_all &vars)
std::vector< int > m_trackstub_candidate_num_hits
Definition: variables.h:358
std::vector< double > m_reco_shower_conversion_distance
Definition: variables.h:799
void AnalyzeMCTruths(std::vector< art::Ptr< simb::MCTruth >> &mcTruthVector, std::vector< art::Ptr< simb::MCParticle >> &mcParticleVector, var_all &vars, para_all &paras)
Definition: analyze_MC.cxx:378
void AnalyzeRecoMCSlices(std::string signal_def, std::vector< PandoraPFParticle > all_PPFPs, std::map< int, art::Ptr< simb::MCParticle >> &MCParticleToTrackIDMap, std::map< art::Ptr< recob::Shower >, art::Ptr< simb::MCParticle > > &showerToMCParticleMap, std::map< art::Ptr< recob::Track >, art::Ptr< simb::MCParticle > > &trackToMCParticleMap, var_all &vars, para_all &paras)
Definition: analyze_MC.cxx:253
std::vector< double > m_sss_candidate_max_tick
Definition: variables.h:215
std::vector< double > m_trackstub_candidate_mean_ADC
Definition: variables.h:363
void AnalyzeFlashes(const std::vector< art::Ptr< recob::OpFlash >> &flashes, art::Handle< std::vector< sbn::crt::CRTHit >> crthit_h, double evt_timeGPS_nsec, std::map< art::Ptr< recob::OpFlash >, std::vector< art::Ptr< sbn::crt::CRTHit >>> crtvetoToFlashMap, var_all &vars, para_all &paras)
art framework interface to geometry description
BEGIN_PROLOG could also be cout
std::vector< int > m_sss_candidate_num_ticks
Definition: variables.h:205
std::vector< double > m_sss_candidate_mean_wire
Definition: variables.h:219
std::vector< int > m_trackstub_candidate_matched
Definition: variables.h:385
std::string s_geantModuleLabel
Definition: variables.h:103
void ClearSlices(var_all &vars)
std::vector< double > m_reco_track_length
Definition: variables.h:575
std::vector< double > m_sss_candidate_mean_ADC
Definition: variables.h:208
double s_SEAviewHitThreshold
Definition: variables.h:77
std::vector< double > m_sss_candidate_min_tick
Definition: variables.h:216
void ClearShowers(var_all &vars)
int filterConsideredHits(double dist_to_vertex)
Definition: SEAviewer.cxx:172
std::vector< double > m_sss_candidate_min_wire
Definition: variables.h:217
void CreateMetaBranches(var_all &vars)
void ClearTracks(var_all &vars)