All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FlashResAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: FlashResAna
3 // Plugin Type: analyzer (art v3_05_01)
4 // File: FlashResAna_module.cc
5 //
6 // Generated at Sun Jun 28 11:06:51 2020 by Christopher Hilgenberg using cetskelgen
7 // from cetlib version v3_10_00.
8 ////////////////////////////////////////////////////////////////////////
9 
10 //framework includes
11 #include "art/Framework/Core/EDAnalyzer.h"
12 #include "art/Framework/Core/ModuleMacros.h"
13 #include "art/Framework/Principal/Event.h"
14 #include "art/Framework/Principal/Handle.h"
15 #include "art/Framework/Principal/Run.h"
16 #include "art/Framework/Principal/SubRun.h"
17 #include "canvas/Utilities/InputTag.h"
18 #include "fhiclcpp/ParameterSet.h"
19 #include "messagefacility/MessageLogger/MessageLogger.h"
20 #include "art_root_io/TFileService.h"
21 #include "art/Framework/Services/Registry/ServiceHandle.h"
22 #include "canvas/Persistency/Common/FindManyP.h"
23 
24 //larsoft includes
25 //#include "nusimdata/SimulationBase/GTruth.h"
26 #include "nusimdata/SimulationBase/MCTruth.h"
27 #include "nusimdata/SimulationBase/MCParticle.h"
28 #include "larcore/Geometry/Geometry.h"
29 #include "larcore/CoreUtils/ServiceUtil.h" // lar::providerFrom()
30 #include "larcorealg/Geometry/GeometryCore.h"
31 #include "lardataobj/RecoBase/OpHit.h"
32 #include "lardataobj/RecoBase/OpFlash.h"
33 #include "lardataobj/Simulation/SimPhotons.h"
34 
35 //C++ includes
36 #include <vector>
37 #include <map>
38 
39 //ROOT includes
40 #include <TTree.h>
41 #include <TH1F.h>
42 
43 namespace icarus {
44  class FlashResAna;
45 }
46 
47 using std::vector;
48 using std::map;
49 
50 class icarus::FlashResAna : public art::EDAnalyzer {
51 public:
52  explicit FlashResAna(fhicl::ParameterSet const& p);
53  // The compiler-generated destructor is fine for non-base
54  // classes without bare pointers or other resource use.
55 
56  // Plugins should not be copied or assigned.
57  FlashResAna(FlashResAna const&) = delete;
58  FlashResAna(FlashResAna&&) = delete;
59  FlashResAna& operator=(FlashResAna const&) = delete;
60  FlashResAna& operator=(FlashResAna&&) = delete;
61 
62  // Required functions.
63  void analyze(art::Event const& e) override;
64 
65 private:
66 
67  art::InputTag fGenLabel;
68  art::InputTag fPhotLabel;
69  art::InputTag fFlashLabel0;
70  art::InputTag fFlashLabel1;
71  art::InputTag fFlashLabel2;
72  art::InputTag fFlashLabel3;
73  art::InputTag fHitLabel;
74 
75  map<int,art::InputTag> fFlashLabels;
76 
77  TTree* fTree;
78  int fNFlash;
79  bool fCC;
80  double fFlashPE;
81  int fFlashTPC;
82  double fFlashWidth[4];
83  int fMaxChan;
84  double fMaxPE;
85  double fNuE;
86  double fNuXYZT[4];
87  double fFlashXYZT[4];
88  double fFlashXYZTMin[4];
89  double fFlashXYZTMax[4];
90  double fDelta[4];
91  double fDeltaPmt;
92  bool fInFrame;
93  int fNPmtFlash;
94  double fHitXYZT[4];
95  double fHitDist;
96  double fHitTof;
97  int fNHit;
98  double fHitPE;
99  double fTPhot;
100 
101  TH1F* fDeltaTFlash;
102 
103  unsigned int GetMaxChan(vector<double> const& pes);
104 
105  geo::GeometryCore const* fGeometryService; ///< pointer to Geometry provider
106 
107 };
108 
109 
110 icarus::FlashResAna::FlashResAna(fhicl::ParameterSet const& p)
111  : EDAnalyzer{p},
112  fGenLabel(p.get<art::InputTag>("GenLabel","generator")),
113  fPhotLabel(p.get<art::InputTag>("PhotLabel","largeant")),
114  fFlashLabel0(p.get<art::InputTag>("FlashLabel0","opflashTPC0")),
115  fFlashLabel1(p.get<art::InputTag>("FlashLabel1","opflashTPC1")),
116  fFlashLabel2(p.get<art::InputTag>("FlashLabel2","opflashTPC2")),
117  fFlashLabel3(p.get<art::InputTag>("FlashLabel3","opflashTPC3")),
118  fHitLabel(p.get<art::InputTag>("HitLabel","ophit"))
119 {
120  fFlashLabels[0] = fFlashLabel0;
121  fFlashLabels[1] = fFlashLabel1;
122  fFlashLabels[2] = fFlashLabel2;
123  fFlashLabels[3] = fFlashLabel3;
124 
125  art::ServiceHandle<art::TFileService> tfs;
126 
127  fTree = tfs->make<TTree>("anatree","OpFlash resolution");
128  fTree->Branch("nFlash", &fNFlash, "nFlash/I");
129  fTree->Branch("nuE", &fNuE, "nuE/D");
130  fTree->Branch("cc", &fCC, "cc/O");
131  fTree->Branch("nuXYZT", fNuXYZT, "nuXYZT[4]/D");
132  fTree->Branch("flashXYZT", fFlashXYZT, "flashXYZT[4]/D");
133  fTree->Branch("flashXYZTMin", fFlashXYZTMin, "flashXYZTMin[4]/D");
134  fTree->Branch("flashXYZTMax", fFlashXYZTMax, "flashXYZTMax[4]/D");
135  fTree->Branch("flashWidth",fFlashWidth,"flashWidth[4]/D");
136  fTree->Branch("flashPE", &fFlashPE, "flashPE/D");
137  fTree->Branch("flashTPC", &fFlashTPC, "flashTPC/I");
138  fTree->Branch("maxChan", &fMaxChan, "maxChan/i");
139  fTree->Branch("maxPE", &fMaxPE, "maxPE/D");
140  fTree->Branch("delta", fDelta, "delta[4]/D");
141  fTree->Branch("inFrame", &fInFrame, "inFrame/O");
142  fTree->Branch("deltaPmt", &fDeltaPmt, "deltaPmt/D");
143  fTree->Branch("nPmtFlash", &fNPmtFlash,"nPmtFlash/I");
144  fTree->Branch("hitXYZT", fHitXYZT, "hittXYZT[4]/D");
145  fTree->Branch("hitDist", &fHitDist, "hitDist/D");
146  fTree->Branch("hitTof", &fHitTof, "hitTof/D");
147  fTree->Branch("nHit", &fNHit, "nHit/I");
148  fTree->Branch("hitPE", &fHitPE, "hitPE/D");
149  fTree->Branch("tPhot", &fTPhot, "tPhot/D");
150 
151  fDeltaTFlash = tfs->make<TH1F>("dtflash","delay between flashes;#Delta t [ns];",25000,0,500000);
152 
153  fGeometryService = lar::providerFrom<geo::Geometry>();
154 
155 }
156 
157 void icarus::FlashResAna::analyze(art::Event const& e)
158 {
159  //initialize for each event
160  fNFlash = 0;
161  fCC = false;
162  fFlashPE = 0.;
163  fFlashTPC = -1;
164  fMaxChan = -1;
165  fMaxPE = 0.;
166  fNuE = 0.;
167  fDeltaPmt = 0.;
168  fInFrame = false;;
169  fNPmtFlash= -1;
170  double dt = DBL_MAX;
171  geo::CryostatGeo const& cryo0 = fGeometryService->Cryostat(0);
172 
173  for(size_t i=0; i<4; i++) {
174  fNuXYZT[i] = dt;
175  fFlashXYZT[i] = dt;
176  fFlashXYZTMin[i] = dt;
177  fFlashXYZTMax[i] = dt;
178  fFlashWidth[i] = dt;
179  fDelta[i] = dt;
180  }
181 
182 
183  //neutrino truth
184  auto const& mctruths = //vector of MCTruths from GENIE
185  *e.getValidHandle<vector<simb::MCTruth>>(fGenLabel);
186 
187  auto const& nu = mctruths[0].GetNeutrino();
188  fNuE = nu.Nu().E();
189  if(nu.CCNC()==0)
190  fCC = true;
191 
192  const TLorentzVector nuxyzt = nu.Nu().Position(0);
193  fNuXYZT[0] = nuxyzt.X();
194  fNuXYZT[1] = nuxyzt.Y();
195  fNuXYZT[2] = nuxyzt.Z();
196  fNuXYZT[3] = nuxyzt.T();
197 
198  //SimPhot
199  art::Handle< std::vector<sim::SimPhotons>> photHandle;
200  std::vector< art::Ptr<sim::SimPhotons> > photList;
201  if( e.getByLabel(fPhotLabel, photHandle) )
202  art::fill_ptr_vector(photList, photHandle);
203 
204  double tphotmin = DBL_MAX;
205  for(auto const& phot : photList){
206  for(size_t i=0; i<phot->size(); i++) {
207  if ((phot->at(i)).Time < tphotmin)
208  tphotmin = (phot->at(i)).Time;
209  }
210  }
211  fTPhot = tphotmin;
212 
213  //OpFlash
214  map<int,art::Handle< std::vector<recob::OpFlash> > > flashHandles;
215  map<int,vector< art::Ptr<recob::OpFlash> >> flashLists;
216  if( e.getByLabel(fFlashLabel0,flashHandles[0]) )
217  art::fill_ptr_vector(flashLists[0], flashHandles[0]);
218  if( e.getByLabel(fFlashLabel1,flashHandles[1]) )
219  art::fill_ptr_vector(flashLists[1], flashHandles[1]);
220  if( e.getByLabel(fFlashLabel2,flashHandles[2]) )
221  art::fill_ptr_vector(flashLists[2], flashHandles[2]);
222  if( e.getByLabel(fFlashLabel3,flashHandles[3]) )
223  art::fill_ptr_vector(flashLists[3], flashHandles[3]);
224 
225  for(auto flashList : flashLists) {
226 
227  int tpc = flashList.first;
228  size_t nflashtpc = flashList.second.size();
229  if(nflashtpc==0) continue;
230  fNFlash += nflashtpc;
231 
232  std::sort(flashList.second.begin(),flashList.second.end(),
233  [](art::Ptr<recob::OpFlash> f1, art::Ptr<recob::OpFlash> f2)
234  {return f1->Time()<f2->Time();}
235  );
236 
237  art::FindManyP<recob::OpHit> findManyHits(
238  flashHandles[tpc], e, fFlashLabels[tpc]);
239 
240  //auto const& flash0 = flashList.second[0];
241  for(size_t iflash=0; iflash < nflashtpc; iflash++) {
242  auto const& flash0 = flashList.second[iflash]; //) {
243  if(abs(fNuXYZT[3]-flash0->Time()*1.0e3)<abs(dt)) {
244  dt = flash0->Time()*1.0e3 - fNuXYZT[3];
245  fFlashXYZT[0] = 0.;
246  fFlashXYZT[1] = flash0->YCenter();
247  fFlashXYZT[2] = flash0->ZCenter();
248  fFlashXYZT[3] = flash0->Time()*1.0e3;
249  fFlashWidth[0] = 0.;
250  fFlashWidth[1] = flash0->YWidth();
251  fFlashWidth[2] = flash0->ZWidth();
252  fFlashWidth[3] = flash0->TimeWidth();
253  fFlashPE = flash0->TotalPE();
254  fFlashTPC = tpc;
255  auto const& pes= flash0->PEs();
256  fNPmtFlash = pes.size();
257  fMaxChan = GetMaxChan(pes);
258  fMaxPE = flash0->PE(fMaxChan);
259  fInFrame = flash0->InBeamFrame();
260  for(size_t i=0; i<4; i++)
261  fDelta[i] = fFlashXYZT[i] - fNuXYZT[i];
262 
263  double tPmtMax=-DBL_MAX, tPmtMin=DBL_MAX;
264  vector<art::Ptr<recob::OpHit>> hits = findManyHits.at(iflash);
265  for(auto const& hit : hits) {
266  if(hit->OpChannel() == fMaxChan){
267  fDeltaPmt = hit->PeakTime()*1.0e3 - fNuXYZT[3];
268  }
269  if(hit->PeakTime()<tPmtMin){
270  tPmtMin = hit->PeakTime();
271  geo::OpDetGeo const& opDet = cryo0.OpDet(hit->OpChannel());
272  double pos[3];
273  opDet.GetCenter(pos);
274  for(int i=0; i<3; i++) fFlashXYZTMin[i] = pos[i];
275  fFlashXYZTMin[3] = tPmtMin;
276  }
277  if(hit->PeakTime()>tPmtMax){
278  tPmtMax = hit->PeakTime();
279  geo::OpDetGeo const& opDet = cryo0.OpDet(hit->OpChannel());
280  double pos[3];
281  opDet.GetCenter(pos);
282  for(int i=0; i<3; i++) fFlashXYZTMax[i] = pos[i];
283  fFlashXYZTMax[3] = tPmtMax;
284  }
285  }
286  }
287  }
288 
289  if(nflashtpc>1)
290  for(size_t i=0; i<nflashtpc-1; i++) {
291  auto const& flash = flashList.second.at(i);
292  auto const& flashnext = flashList.second.at(i+1);
293  fDeltaTFlash->Fill((flashnext->Time()-flash->Time())*1.0e3);
294  }
295  }//for flash lists
296 
297  //OpHits
298  art::Handle< std::vector<recob::OpHit> > opHitHandle;
299  std::vector< art::Ptr<recob::OpHit> > opHitList;
300 
301  if( e.getByLabel(fHitLabel,opHitHandle) ) {
302 
303  art::fill_ptr_vector(opHitList, opHitHandle);
304  fNHit = opHitList.size();
305  if(fNHit==0){
306  std::cout << "empty OpHit vector!" << std::endl;
307  for(int i=0; i<4; i++) fHitXYZT[i] = DBL_MAX;
308  fHitDist = DBL_MAX;
309  fHitTof = DBL_MAX;
310  fHitPE = DBL_MAX;
311  }
312 
313  else{
314  std::sort(opHitList.begin(),opHitList.end(),
315  [](art::Ptr<recob::OpHit> h1, art::Ptr<recob::OpHit> h2)
316  {return h1->PeakTime()<h2->PeakTime();}
317  );
318  if(opHitList.size()>1 && opHitList[0]->PeakTime()>opHitList.back()->PeakTime())
319  std::cout << "OpHits time sort failed!" << std::endl;
320 
321  geo::OpDetGeo const& opDet = cryo0.OpDet(opHitList[0]->OpChannel());
322  double pos[3];
323  opDet.GetCenter(pos);
324  for(int i=0; i<3; i++) fHitXYZT[i] = pos[i];
325  fHitXYZT[3] = opHitList[0]->PeakTime()*1.0e3; //time in ns
326  fHitDist = 0.;
327  for(int i=0; i<3; i++)
328  fHitDist+=pow(fHitXYZT[i]-fNuXYZT[i],2);
329  fHitDist = sqrt(fHitDist);
330  fHitTof = fHitXYZT[3] - fNuXYZT[3];
331  fHitPE = opHitList[0]->PE();
332  }
333  }
334  else{
335  std::cout << "no OpHit handle found!" << std::endl;
336  }
337 
338  fTree->Fill();
339 
340 }
341 
342 unsigned int icarus::FlashResAna::GetMaxChan(vector<double> const& pes){
343 
344  size_t imax = 0;
345  double max = 0.;
346  for(size_t i=0; i<pes.size(); i++){
347  if(pes[i]>max) {
348  max = pes[i];
349  imax = i;
350  }
351  }
352 
353  return imax;
354 }
355 
356 DEFINE_ART_MODULE(icarus::FlashResAna)
icarus::FlashResAna::analyze
void analyze(art::Event const &e) override
Definition: FlashResAna_module.cc:157
icarus::FlashResAna::fGeometryService
geo::GeometryCore const * fGeometryService
pointer to Geometry provider
Definition: FlashResAna_module.cc:105
ServiceUtil.h
Utilities related to art service access.
icarus::FlashResAna::fFlashLabel1
art::InputTag fFlashLabel1
Definition: FlashResAna_module.cc:70
p
pdgs p
Definition: selectors.fcl:22
icarus::FlashResAna::fFlashLabel0
art::InputTag fFlashLabel0
Definition: FlashResAna_module.cc:69
geo::OpDetGeo::GetCenter
void GetCenter(double *xyz, double localz=0.0) const
Definition: OpDetGeo.cxx:40
hit
process_name hit
Definition: cheaterreco.fcl:51
geo::CryostatGeo
Geometry information for a single cryostat.
Definition: CryostatGeo.h:43
lar::dump::vector
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
GeometryCore.h
Access the description of detector geometry.
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
SimPhotons.h
Simulation objects for optical detectors.
icarus::FlashResAna::fFlashLabels
map< int, art::InputTag > fFlashLabels
Definition: FlashResAna_module.cc:75
icarus::FlashResAna::fFlashLabel3
art::InputTag fFlashLabel3
Definition: FlashResAna_module.cc:72
geo::CryostatGeo::OpDet
const OpDetGeo & OpDet(unsigned int iopdet) const
Return the iopdet&#39;th optical detector in the cryostat.
Definition: CryostatGeo.cxx:107
geo::GeometryCore
Description of geometry of one entire detector.
Definition: GeometryCore.h:1476
icarus::FlashResAna::FlashResAna
FlashResAna(fhicl::ParameterSet const &p)
Definition: FlashResAna_module.cc:110
icarus::FlashResAna::operator=
FlashResAna & operator=(FlashResAna const &)=delete
e
do i e
Definition: updateVersionList.sh:68
icarus::FlashResAna::GetMaxChan
unsigned int GetMaxChan(vector< double > const &pes)
Definition: FlashResAna_module.cc:342
icarus::FlashResAna::fFlashLabel2
art::InputTag fFlashLabel2
Definition: FlashResAna_module.cc:71
tfs
art::ServiceHandle< art::TFileService > tfs
Definition: NuVertexChargeTree_module.cc:254
nu
BEGIN_PROLOG SN nu
Definition: clustermodules.fcl:121
Geometry.h
art framework interface to geometry description
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
Generated by   doxygen 1.8.5