ProtoDUNE_opdet_eventdisplay_module.cc
Go to the documentation of this file.
1 // -*- mode: c++; c-basic-offset: 2; -*-
2 // This analyzer writes out a folder containing the
3 // PDS event display. For every OpFLsh, the OpHits are
4 // plotted.
5 //
6 //
7 
8 
9 #ifndef ProtoDUNE_opdet_eventdisplay_H
10 #define ProtoDUNE_opdet_eventdisplay_H
11 
12 // ROOT includes
13 #include "TH1.h"
14 #include "TH2.h"
15 #include "TLorentzVector.h"
16 #include "TVector3.h"
17 #include "TTree.h"
18 #include "TLine.h"
19 #include "TLegend.h"
20 
21 // C++ includes
22 #include <map>
23 #include <vector>
24 #include <iostream>
25 #include <cstring>
26 #include <sstream>
27 #include "math.h"
28 #include <climits>
29 #include <functional>
30 #include <algorithm>
31 #include <map>
32 
33 
34 // LArSoft includes
42 
44 //#include "OpticalDetector/OpDigiProperties.h"
45 
47 #include <numeric>
48 
49 // ART includes.
52 #include "fhiclcpp/ParameterSet.h"
57 #include "art_root_io/TFileService.h"
58 #include "art_root_io/TFileDirectory.h"
61 #include "canvas/Persistency/Common/FindManyP.h"
62 
63 namespace opdet {
64 
66  public:
67 
68  // Standard constructor and destructor for an ART module.
71 
72  // This method is called once, at the start of the job. In this
73  // example, it will define the histogram we'll write.
74  void beginJob();
75 
76  // The analyzer routine, called once per event.
77  void analyze (const art::Event&);
78 
79  private:
80 
81  // The stuff below is the part you'll most likely have to change to
82  // go from this custom example to your own task.
83 
84  // The parameters we'll read from the .fcl file.
85  std::string fOpFlashModuleLabel; // Input tag for OpFlash collection
86  std::string fOpHitModuleLabel; // Input tag for OpHit collection
87  float fSampleFreq; // in MHz
88  float fTimeBegin; // in us
89  float fTimeEnd; // in us
90 
91  float fYMin, fYMax, fZMin, fZMax;
92 
94 
95 // bool fMakeFlashTimeHist;
96 // bool fMakeFlashPosHist;
97 // bool fMakePerFlashHists;
98 // bool fMakePerOpHitTree;
99 
100 // TTree * fPerOpHitTree;
101 
102  Int_t fEventID;
103  Int_t fFlashID;
104  Int_t fHitID;
105  Float_t fFlashTime;
106 // Float_t fAbsTime;
107 // bool fInBeamFrame;
108 // int fOnBeamTime;
109 // Float_t fTotalPE;
110 // Int_t fFlashFrame;
111 
112 // Float_t fNPe;
113 // Float_t fYCenter;
114 // Float_t fYWidth;
115 // Float_t fZCenter;
116 // Float_t fZWidth;
117 
118  Int_t fOpChannel;
119  Float_t fPeakTimeAbs;
120  Float_t fPeakTime;
121 // Int_t fFrame;
122 // Float_t fWidth;
123 // Float_t fArea;
124 // Float_t fAmplitude;
125  Float_t fPE;
126 // Float_t fFastToTotal;
127  Float_t fYOpCenter;
128  Float_t fZOpCenter;
129  Float_t fXOpCenter;
130 
131 // int fNFlashes;
132  std::vector< int > fFlashIDVector;
133  std::vector< float > fYCenterVector;
134  std::vector< float > fZCenterVector;
135  std::vector< float > fYWidthVector;
136  std::vector< float > fZWidthVector;
137  std::vector< float > fFlashTimeVector;
138  std::vector< float > fAbsTimeVector;
139  std::vector< int > fFlashFrameVector;
140  std::vector< bool > fInBeamFrameVector;
141  std::vector< int > fOnBeamTimeVector;
142  std::vector< float > fTotalPEVector;
143 // int fNChannels;
144  std::vector< float > fPEsPerFlashPerChannelVector;
145  };
146 
147 }
148 
149 #endif // ProtoDUNE_opdet_eventdisplay_H
150 
151 namespace opdet {
152 
153  //-----------------------------------------------------------------------
154  // Constructor
156  : EDAnalyzer(pset)
157  {
158 
159  // Indicate that the Input Module comes from .fcl
160  fOpFlashModuleLabel = pset.get<std::string>("OpFlashModuleLabel");
161  fOpHitModuleLabel = pset.get<std::string>("OpHitModuleLabel");
162 
163 // art::ServiceHandle<OpDigiProperties> odp;
164 // fTimeBegin = odp->TimeBegin();
165 // fTimeEnd = odp->TimeEnd();
166 // fSampleFreq = odp->SampleFreq();
167 
168  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataForJob();
169  fTimeBegin = clockData.OpticalClock().Time();
170  fTimeEnd = clockData.OpticalClock().FramePeriod();
171  fSampleFreq = clockData.OpticalClock().Frequency();
172 
173  fYMin = pset.get<float>("YMin");
174  fYMax = pset.get<float>("YMax");
175  fZMin = pset.get<float>("ZMin");
176  fZMax = pset.get<float>("ZMax");
177 
178  PosHistYRes = 20;
179  PosHistZRes = 42;
180 
182 
183  fFlashID = 0;
184  }
185 
186  //-----------------------------------------------------------------------
187  // Destructor
189  {}
190 
191  //-----------------------------------------------------------------------
193  {}
194 
195  //-----------------------------------------------------------------------
197  {
198 
199  float displace[288]={0};
200 
201  for(int j=0;j<287;j++) displace[j]=j;
202  for(int j=0;j<284;j+=4){
203  displace[j]=-100;
204  displace[j+1]=-40;
205  displace[j+2]=5;
206  displace[j+3]=60;
207  if((j>=132 && j<=143) || (j>=264 && j<=275)){
208  displace[j]=-100;
209  displace[j+1]=-75;
210  displace[j+2]=-60;
211  displace[j+3]=-40;
212  displace[j+4]=-20;
213  displace[j+5]=-5;
214  displace[j+6]=5;
215  displace[j+7]=20;
216  displace[j+8]=40;
217  displace[j+9]=60;
218  displace[j+10]=75;
219  displace[j+11]=100;
220  j+=8;
221  }
222 
223  }
224 
225  Float_t binsz[] = {0,10, 29, 48, 66, 84, 102, 120, 138, 156, 174, 192, 211, 230, 240, 257, 278, 296, 314, 332, 350, 368, 386, 404, 422, 441, 460, 470, 489, 508, 526, 544, 562, 580, 598, 616, 634, 650, 671, 690, 700};
226  Int_t binnumz = sizeof(binsz)/sizeof(Float_t) - 1;
227 
228  Float_t binsy[]={-600,-540, -480, -420, -360, -300, -240, -180, -120, -60, 0, 60, 120, 180, 240, 300, 360, 420, 480, 540, 600};
229  Int_t binnumy = sizeof(binsy)/sizeof(Float_t) - 1;
230 
231 
232 
233  // Get flashes from event
234  auto FlashHandle = evt.getHandle< std::vector< recob::OpFlash > >(fOpFlashModuleLabel);
235 
236  // Get assosciations between flashes and hits
237  art::FindManyP< recob::OpHit > Assns(FlashHandle, evt, fOpFlashModuleLabel);
238 
239  // Create string for histogram name
240  char HistName[50];
241 
242  fFlashID = 0;
243 
244  // Access ART's TFileService, which will handle creating and writing
245  // histograms for us.
247 
248  //std::vector<TH1D*> FlashHist;
249  std::vector<TH2D*> DispHist;
250 
251  fEventID = evt.id().event();
252 
253  sprintf(HistName, "Event %d Disps", evt.id().event());
254 
256  //unsigned int NOpChannels = geom->NOpChannels();
257 
258  // For every OpFlash in the vector
259  for(unsigned int i = 0; i < FlashHandle->size(); ++i)
260  {
261 
262  // Get OpFlash
263  art::Ptr< recob::OpFlash > TheFlashPtr(FlashHandle, i);
264  recob::OpFlash TheFlash = *TheFlashPtr;
265 
266  fFlashTime = TheFlash.Time();
267  fFlashID = i; //++;
268 
269 
270  TH2D * ThisEventDisp = nullptr;
271  ThisEventDisp = tfs->make<TH2D>(HistName,"PE from OpHits in OpFlash;z ;DaS RaS (y)",
272  binnumz, binsz,
273  binnumy, binsy);
274 
275 
276  double xyz[3];
277  auto OpHitHandle = evt.getHandle< std::vector< recob::OpHit > >(fOpHitModuleLabel);
278 
279  for(size_t i=0; i!=OpHitHandle->size(); ++i)
280  {
281  fOpChannel = OpHitHandle->at(i).OpChannel();
282  geom->OpDetGeoFromOpChannel(fOpChannel).GetCenter(xyz);
283  fXOpCenter = xyz[0];
284  fYOpCenter = xyz[1];
285  fZOpCenter = xyz[2];
286  fPeakTimeAbs = OpHitHandle->at(i).PeakTimeAbs();
287  fPeakTime = OpHitHandle->at(i).PeakTime();
288  fPE = OpHitHandle->at(i).PE();
289  fHitID = i;
290  //fPerOpHitTree->Fill();
291  //ThisEventDisp->Fill(fZOpCenter, fYOpCenter,fPE);
292  if(fXOpCenter >0){
293  if(fOpChannel>=132 && fOpChannel<=143)
294  ThisEventDisp->Fill(fZOpCenter+displace[fOpChannel], fYOpCenter,fPE);
295  else{
296  ThisEventDisp->Fill(fZOpCenter+displace[fOpChannel], fYOpCenter,fPE);
297  ThisEventDisp->Fill(fZOpCenter+displace[fOpChannel]+18, fYOpCenter,fPE);
298  ThisEventDisp->Fill(fZOpCenter+displace[fOpChannel]+36, fYOpCenter,fPE);
299  }
300  }
301  else if(fXOpCenter<0){
302  if(fOpChannel>=264 && fOpChannel<=275)
303  ThisEventDisp->Fill(fZOpCenter+displace[fOpChannel], fYOpCenter-600,fPE);
304  else{
305  ThisEventDisp->Fill(fZOpCenter+displace[fOpChannel], fYOpCenter-600,fPE);
306  ThisEventDisp->Fill(fZOpCenter+displace[fOpChannel]+18, fYOpCenter-600,fPE);
307  ThisEventDisp->Fill(fZOpCenter+displace[fOpChannel]+38, fYOpCenter-600,fPE);
308  }
309  }
310  TLine *line = new TLine(0,0,600,0);
311  line->SetLineColor(kBlack);
312  line->Draw();
313  }
314 
315 
316 
317 
318 
319  }
320 
321 
322  }
323 
324 } // namespace opdet
325 
326 namespace opdet {
328 }
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
std::string string
Definition: nybbler.cc:12
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:25
art framework interface to geometry description
double Time() const
Definition: OpFlash.h:106
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
T get(std::string const &key) const
Definition: ParameterSet.h:271
Encapsulate the geometry of an optical detector.
void line(double t, double *p, double &x, double &y, double &z)
ProtoDUNE_opdet_eventdisplay(const fhicl::ParameterSet &)
EventNumber_t event() const
Definition: EventID.h:116
Access the description of detector geometry.
TCEvent evt
Definition: DataStructs.cxx:7
EventID id() const
Definition: Event.cc:34
pure virtual base interface for detector clocks