All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PMTcoordinates_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: PMTcoordinates
3 // Plugin Type: analyzer (art v2_07_03)
4 // File: PMTcoordinates_module.cc
5 //
6 // Generated at Mon Sep 25 15:10:31 2017 by Andrea Falcone using cetskelgen
7 // from cetlib version v3_00_01.
8 ////////////////////////////////////////////////////////////////////////
9 
10 // ##########################
11 // ### Framework includes ###
12 // ##########################
13 #include "art/Framework/Core/EDAnalyzer.h"
14 #include "art/Framework/Core/ModuleMacros.h"
15 #include "art/Framework/Principal/Event.h"
16 #include "fhiclcpp/ParameterSet.h"
17 #include "art/Framework/Principal/Run.h"
18 #include "art/Framework/Principal/SubRun.h"
19 #include "art/Framework/Principal/Handle.h"
20 #include "canvas/Persistency/Common/Ptr.h"
21 #include "canvas/Persistency/Common/PtrVector.h"
22 #include "art/Framework/Services/Registry/ServiceHandle.h"
23 #include "art_root_io/TFileService.h"
24 #include "art_root_io/TFileDirectory.h"
25 #include "canvas/Persistency/Common/FindOneP.h"
26 #include "canvas/Persistency/Common/FindManyP.h"
27 #include "canvas/Utilities/InputTag.h"
28 #include "messagefacility/MessageLogger/MessageLogger.h"
29 //#include "cetlib/maybe_ref.h"
30 
31 // ########################
32 // ### LArSoft includes ###
33 // ########################
34 #include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
35 #include "larcoreobj/SimpleTypesAndConstants/RawTypes.h" // raw::ChannelID_t
36 #include "larcore/Geometry/Geometry.h"
37 #include "lardataobj/RecoBase/Wire.h"
38 #include "lardataobj/RecoBase/Hit.h"
39 #include "lardataobj/RecoBase/Cluster.h"
40 #include "lardataobj/RecoBase/Track.h"
41 #include "lardataobj/RecoBase/TrackHitMeta.h"
42 #include "lardataobj/RecoBase/Vertex.h"
43 #include "lardataobj/RecoBase/SpacePoint.h"
44 #include "lardata/ArtDataHelper/TrackUtils.h" // lar::util::TrackPitchInView()
45 #include "lardata/DetectorInfoServices/LArPropertiesService.h"
46 #include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
47 #include "lardata/Utilities/AssociationUtil.h"
48 
49 #include "larana/OpticalDetector/SimPhotonCounter.h"
50 #include "lardataobj/Simulation/SimPhotons.h"
51 
52 #include "lardataobj/Simulation/SimChannel.h"
53 
54 #include "nusimdata/SimulationBase/MCParticle.h"
55 #include "nusimdata/SimulationBase/MCTruth.h"
56 #include "nusimdata/SimulationBase/MCNeutrino.h"
57 
58 #include "nug4/ParticleNavigation/ParticleList.h"
59 
60 #include "lardataobj/Simulation/sim.h"
61 
62 // #####################
63 // ### ROOT includes ###
64 // #####################
65 #include "TComplex.h"
66 #include "TFile.h"
67 #include "TH2D.h"
68 #include "TF1.h"
69 #include "TTree.h"
70 #include "TTimeStamp.h"
71 
72 ////////////////////////////////// Define some constant variable //////////////////////////////
73 const int nPMTs = 360;
74 //const int PMTs_per_TPC = 90;
75 const int MaxPhotons = 10000;
76 
77 namespace icarus {
78  class PMTcoordinates;
79 }
80 
81 
82 class icarus::PMTcoordinates : public art::EDAnalyzer
83 {
84 public:
85 explicit PMTcoordinates(fhicl::ParameterSet const & p);
86 
87 // Plugins should not be copied or assigned.
88 PMTcoordinates(PMTcoordinates const &) = delete;
89 PMTcoordinates(PMTcoordinates &&) = delete;
90 PMTcoordinates & operator = (PMTcoordinates const &) = delete;
91 PMTcoordinates & operator = (PMTcoordinates &&) = delete;
92 
93 // Required functions.
94 void analyze(art::Event const & e) override;
95 
96 // Selected optional functions.
97 void beginJob() override;
98 //void reconfigure(fhicl::ParameterSet const & p);
99 
100 private:
101 
102 //TRandom* Ran;
103 
104 TTree* fTree;
105 
106 ////////////////////////////////// Variable in th tree//////////////////////////////
107 
108 int event;
109 
110 int event_type;
111 
112 int is_Neutrino;
113 int Neutrino_Interaction;
114 
115 int noPMT[nPMTs];
116 
117 int turned_PMT;
118 
119 double PMTx[nPMTs];
120 double PMTy[nPMTs];
121 double PMTz[nPMTs];
122 
123 int Cryostat[nPMTs];
124 int TPC[nPMTs];
125 
126 float photons_collected[nPMTs];
127 float QE_photons_collected[nPMTs];
128 float photon_time[nPMTs][MaxPhotons];
129 
130 float firstphoton_time[nPMTs];
131 
132 float true_barycentre_x;
133 float true_barycentre_y;
134 float true_barycentre_z;
135 
136 double vertex_x;
137 double vertex_y;
138 double vertex_z;
139 
140 float reco_barycentre_y;
141 float reco_barycentre_z;
142 
143 float total_quenched_energy;
144 float total_coll_photons;
145 
146 float PMT_error_y;
147 float PMT_error_z;
148 float PMT_total_error;
149 
150 art::InputTag photonLabel;
151 art::InputTag chargeLabel;
152 art::InputTag typoLabel;
153 
154 };
155 
156 
157 icarus::PMTcoordinates::PMTcoordinates(fhicl::ParameterSet const & p)
158  :
159  EDAnalyzer(p),
160  photonLabel(p.get<art::InputTag>("fottoni", "largeant")),
161  chargeLabel(p.get<art::InputTag>("carconi", "largeant")),
162  typoLabel (p.get<art::InputTag>("tiponi", "generator"))
163  // More initializers here.
164 {
165  std::cout << " PMT coordinates constructor " << std::endl;
166 }
167 
168 void icarus::PMTcoordinates::analyze(art::Event const & evt)
169 {
170 ////////////////////////////////// Create the LArsoft services and service handle//////////////////////////////
171 
172 art::ServiceHandle<geo::Geometry> geom;
173 
174 std::vector<sim::SimPhotons> const& optical = *(evt.getValidHandle<std::vector<sim::SimPhotons>>(photonLabel));
175 std::vector<sim::SimChannel> const& charge = *(evt.getValidHandle<std::vector<sim::SimChannel>>(chargeLabel));
176 //std::vector<simb::MCTruth> const& type = *(evt.getValidHandle<std::vector<simb::MCTruth>>(typoLabel));
177 
178 ////////////////////////////////// Event number//////////////////////////////
179 
180 event = evt.id().event();
181 
182 //std::vector< art::Handle< std::vector<simb::MCTruth> > > type;
183 //evt.getManyByType(type);
184 auto type = evt.getMany< std::vector<simb::MCTruth> >();
185 
186 for(size_t mcl = 0; mcl < type.size(); ++mcl)
187 {
188  art::Handle< std::vector<simb::MCTruth> > mclistHandle = type[mcl];
189 
190  for(size_t m = 0; m < mclistHandle->size(); ++m)
191  {
192  art::Ptr<simb::MCTruth> mct(mclistHandle, m);
193 // for(int ipart=0;ipart<mct->NParticles();ipart++)
194 // {
195 // int pdg=mct->GetParticle(ipart).PdgCode();
196 // double xx=mct->GetParticle(ipart).Vx();
197 // double yy=mct->GetParticle(ipart).Vy();
198 // double zz=mct->GetParticle(ipart).Vz();
199 
200  event_type=mct->GetParticle(0).PdgCode();
201  vertex_x=mct->GetParticle(0).Vx();
202  vertex_y=mct->GetParticle(0).Vy();
203  vertex_z=mct->GetParticle(0).Vz();
204 
205  if (event_type==12||event_type==-12||event_type==14||event_type==-14||event_type==16||event_type==-16)
206  {
207  is_Neutrino=1;
208  Neutrino_Interaction=mct->GetNeutrino().InteractionType();
209  }
210  else
211  {
212  is_Neutrino=0;
213  Neutrino_Interaction=-9999;
214  }
215 
216 // }
217  }
218 
219 }
220 
221 
222 ////////////////////////////////// Putting at 0 all the variables//////////////////////////////
223 
224 for (int g=0; g<MaxPhotons; g++)
225 {
226  for (int u=0; u<360; u++)
227  {
228  photon_time[u][g]=0;
229  }
230 }
231 
232 true_barycentre_x =0;
233 true_barycentre_y =0;
234 true_barycentre_z =0;
235 
236 total_quenched_energy =0;
237 
238 ////////////////////////////////// Charge part: identify the baricentre of the event //////////////////////////////
239 
240 for (std::size_t chargechannel = 0; chargechannel<charge.size(); ++chargechannel) //loop on SimChannel
241 {
242  auto const& channeltdcide = charge.at(chargechannel).TDCIDEMap();
243 
244  for (std::size_t TDCnu = 0; TDCnu<channeltdcide.size(); ++TDCnu) //loop on TDC
245  {
246 
247  sim::TDCIDE const& tdcide = channeltdcide.at(TDCnu);
248 
249  for (std::size_t IDEnu = 0; IDEnu<tdcide.second.size(); ++IDEnu) //loop on IDE
250  {
251  sim::IDE const& ida = tdcide.second.at(IDEnu);
252 
253 // std::cout << "IDA " << ida.x << '\t' << ida.y << '\t' << ida.z << std::endl;
254 
255  true_barycentre_x = true_barycentre_x + ida.x*ida.energy;
256  true_barycentre_y = true_barycentre_y + ida.y*ida.energy;
257  true_barycentre_z = true_barycentre_z + ida.z*ida.energy;
258  total_quenched_energy = total_quenched_energy + ida.energy;
259 
260  } //loop on IDE
261 
262  } //loop on TDC
263 
264 }//loop on SimChannel
265 
266 true_barycentre_x = true_barycentre_x/total_quenched_energy;
267 true_barycentre_y = true_barycentre_y/total_quenched_energy;
268 true_barycentre_z = true_barycentre_z/total_quenched_energy;
269 
270 total_quenched_energy = total_quenched_energy/3;
271 
272 ////////////////////////////////// Light part //////////////////////////////////////////////////
273 turned_PMT=0;
274 
275 reco_barycentre_y=0;
276 reco_barycentre_z=0;
277 
278 total_coll_photons=0;
279 
280 for (std::size_t channel = 0; channel < optical.size(); ++channel) {
281 
282  sim::SimPhotons const& photon_vec = optical[channel];
283 
284  noPMT[channel] = channel;
285 
286  photons_collected[channel]= photon_vec.size();
287  std::cout << " channel " << channel << " photons collected " << photons_collected[channel] << std::endl;
288 // double media = photons_collected[channel]*QE;
289 
290 // QE_photons_collected[channel]= Ran.Poisson(media);
291 
292  QE_photons_collected[channel]= 0.06*photons_collected[channel];
293 
294  if (photons_collected[channel]>0){
295 
296  turned_PMT++;
297  }
298 
299  double xyz[3];
300 
301  geom->OpDetGeoFromOpChannel(channel).GetCenter(xyz);
302 
303  PMTx[channel] = xyz[0];
304  PMTy[channel] = xyz[1];
305  PMTz[channel] = xyz[2];
306 
307  reco_barycentre_y = reco_barycentre_y + PMTy[channel]*photons_collected[channel];
308  reco_barycentre_z = reco_barycentre_z + PMTz[channel]*photons_collected[channel];
309  total_coll_photons= total_coll_photons + photons_collected[channel];
310  std::cout << " channel " << channel << " total photons " << total_coll_photons << std::endl;
311  firstphoton_time[channel] = 100000000;
312 
313  if (photons_collected[channel]>0)
314  {
315  for (size_t i = 0; i<photon_vec.size() && int(i)< MaxPhotons; ++i)
316  {
317  photon_time[channel][i]= photon_vec.at(i).Time;
318 
319  if (photon_time[channel][i]<firstphoton_time[channel])
320  {
321  firstphoton_time[channel]=photon_time[channel][i];
322  }
323  }
324 
325  }
326 
327 
328 // std::cout << PMTx[channel] << '\t' << PMTy[channel] << '\t' << PMTz[channel] << std::endl;
329 
330  if (PMTx[channel]<0){Cryostat[channel]=0;}
331  if (PMTx[channel]>0){Cryostat[channel]=1;}
332 
333  if (PMTx[channel]<-200){TPC[channel]=0;}
334  if (PMTx[channel]>-200 && PMTx[channel]<0){TPC[channel]=1;}
335  if (PMTx[channel]<200 && PMTx[channel]>0){TPC[channel]=2;}
336  if (PMTx[channel]>200){TPC[channel]=3;}
337 
338 }
339 
340 //total_coll_photons = total_coll_photons;
341 
342 std::cout << " fotoni finale = " <<total_coll_photons <<std::endl;
343 
344 reco_barycentre_y = reco_barycentre_y/total_coll_photons;
345 reco_barycentre_z = reco_barycentre_z/total_coll_photons;
346 
347 
348 PMT_error_y = reco_barycentre_y-true_barycentre_y;
349 PMT_error_z = reco_barycentre_z-true_barycentre_z;
350 PMT_total_error = sqrt((PMT_error_y*PMT_error_y)+(PMT_error_z*PMT_error_z));
351 
352 fTree->Fill();
353  std::cout << " after filling " << fTree << std::endl;
354 }
355 
356 void icarus::PMTcoordinates::beginJob()
357 {
358  std::cout << " PMTcoordinates beginjob " << std::endl;
359 
360 art::ServiceHandle<art::TFileService> tfs;
361 fTree = tfs->make<TTree>("lighttree","tree for the light response");
362 
363 fTree->Branch("event",&event,"event/I");
364 fTree->Branch("event_type",&event_type,"event_type/I");
365 fTree->Branch("is_Neutrino",&is_Neutrino,"is_Neutrino/I");
366 fTree->Branch("Neutrino_Interaction",&Neutrino_Interaction,"Neutrino_Interaction/I");
367 fTree->Branch("total_quenched_energy",&total_quenched_energy,"total_quenched_energy");
368 fTree->Branch("Cryostat",&Cryostat,("Cryostat[" + std::to_string(nPMTs) + "]/I").c_str());
369 fTree->Branch("TPC",&TPC,("TPC[" + std::to_string(nPMTs) + "]/I").c_str());
370 fTree->Branch("noPMT",&noPMT,("noPMT[" + std::to_string(nPMTs) + "]/I").c_str());
371 fTree->Branch("PMTx",&PMTx,("PMTx[" + std::to_string(nPMTs) + "]/D").c_str());
372 fTree->Branch("PMTy",&PMTy,("PMTy[" + std::to_string(nPMTs) + "]/D").c_str());
373 fTree->Branch("PMTz",&PMTz,("PMTz[" + std::to_string(nPMTs) + "]/D").c_str());
374 fTree->Branch("turned_PMT",&turned_PMT,"turned_PMT/I");
375 fTree->Branch("total_coll_photons",&total_coll_photons,"total_coll_photons/F");
376 fTree->Branch("photons_colleted",&photons_collected,("photons_collected[" + std::to_string(nPMTs) + "]/F").c_str());
377 fTree->Branch("QE_photons_colleted",&QE_photons_collected,("QE_photons_collected[" + std::to_string(nPMTs) + "]/F").c_str());
378 fTree->Branch("firstphoton_time",&firstphoton_time,("firstphoton_time[" + std::to_string(nPMTs) + "]/F").c_str());
379 fTree->Branch("photon_time",&photon_time,"photon_time[360][10000]/F");
380 fTree->Branch("vertex_x",&vertex_x,"vertex_x/D");
381 fTree->Branch("vertex_y",&vertex_y,"vertex_y/D");
382 fTree->Branch("vertex_z",&vertex_z,"vertex_z/D");
383 fTree->Branch("true_barycentre_x",&true_barycentre_x,"true_barycentre_x/F");
384 fTree->Branch("true_barycentre_y",&true_barycentre_y,"true_barycentre_y/F");
385 fTree->Branch("true_barycentre_z",&true_barycentre_z,"true_barycentre_z/F");
386 fTree->Branch("reco_barycentre_y",&reco_barycentre_y,"reco_barycentre_y/F");
387 fTree->Branch("reco_barycentre_z",&reco_barycentre_z,"reco_barycentre_z/F");
388 fTree->Branch("PMT_error_y",&PMT_error_y,"PMT_error_y/F");
389 fTree->Branch("PMT_error_z",&PMT_error_z,"PMT_error_z/F");
390 fTree->Branch("PMT_total_error",&PMT_total_error,"PMT_total_error/F");
391 }
392 
393 DEFINE_ART_MODULE(icarus::PMTcoordinates)
icarus::PMTcoordinates::photon_time
float photon_time[nPMTs][MaxPhotons]
Definition: PMTcoordinates_module.cc:128
sim::IDE::z
float z
z position of ionization [cm]
Definition: SimChannel.h:121
sim::TDCIDE
std::pair< unsigned short, std::vector< sim::IDE > > TDCIDE
List of energy deposits at the same time (on this channel)
Definition: SimChannel.h:127
Hit.h
Declaration of signal hit object.
get
get()
Definition: PMTconfigurationExtraction_module.cc:168
p
pdgs p
Definition: selectors.fcl:22
TrackHitMeta.h
Class to keep data related to recob::Hit associated with recob::Track.
icarus::PMTcoordinates::analyze
void analyze(art::Event const &e) override
Definition: PMTcoordinates_module.cc:168
g
BEGIN_PROLOG g
Definition: decoderTools_icarus.fcl:7
sim::IDE::x
float x
x position of ionization [cm]
Definition: SimChannel.h:119
gen_crt_frags.m
tuple m
now if test mode generate materials, CRT shell, world, gdml header else just generate CRT shell for u...
Definition: gen_crt_frags.py:1347
TPC
BEGIN_PROLOG TPC
Definition: sbnd_flashalgo.fcl:14
SimPhotons.h
Simulation objects for optical detectors.
sim::IDE
Ionization at a point of the TPC sensitive volume.
Definition: SimChannel.h:86
sim::IDE::energy
float energy
energy deposited by ionization by this track ID and time [MeV]
Definition: SimChannel.h:118
Cluster.h
Declaration of cluster object.
geo_types.h
Definition of data types for geometry description.
Track.h
Provides recob::Track data product.
sim::IDE::y
float y
y position of ionization [cm]
Definition: SimChannel.h:120
sim::SimPhotons
Collection of photons which recorded on one channel.
Definition: SimPhotons.h:136
icarus::PMTcoordinates::operator=
PMTcoordinates & operator=(PMTcoordinates const &)=delete
icarus::PMTcoordinates::PMTcoordinates
PMTcoordinates(fhicl::ParameterSet const &p)
Definition: PMTcoordinates_module.cc:157
icarus::trigger::to_string
std::string to_string(WindowPattern const &pattern)
Definition: WindowPattern.h:115
MaxPhotons
const int MaxPhotons
Definition: PMTcoordinates_module.cc:75
SimChannel.h
object containing MC truth information necessary for making RawDigits and doing back tracking ...
e
do i e
Definition: updateVersionList.sh:68
Wire.h
Declaration of basic channel signal object.
tfs
art::ServiceHandle< art::TFileService > tfs
Definition: NuVertexChargeTree_module.cc:254
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
TrackUtils.h
Utility functions to extract information from recob::Track
nPMTs
const int nPMTs
Definition: PMTcoordinates_module.cc:73
sim.h
Tools and modules for checking out the basics of the Monte Carlo.
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