LArSoft  v08_35_01
Liquid Argon Software toolkit - http://larsoft.org/
ShowerReco_module.cc
Go to the documentation of this file.
1 //
3 // \file ShowerReco_module.cc
4 //
5 // biagio.rossi@lhep.unibe.ch (FWMK : argoslope.resize(fNPlanes);neut specific)
6 // thomas.strauss@lhep.unibe.ch (ART : general detector)
7 //
8 // andrzej.szelc@yale.edu (port to detector agnostic version)
9 // jasaadi@fnal.gov (Try to rewrite the code to make it more readable and to be able to handle
10 // multiple TPC's and cryostats and more than one cluster per plane.)
11 //
12 // This algorithm is designed to reconstruct showers
13 //
15 
16 // ### Generic C++ includes ###
17 #include <vector>
18 #include <string>
19 #include <cmath> // std::tan() ...
20 
21 // ### Framework includes ###
26 #include "fhiclcpp/ParameterSet.h"
31 #include "art_root_io/TFileService.h"
33 
34 // ### ROOT includes ###
35 #include "TMath.h"
36 #include "TTree.h"
37 
38 // ### LArSoft includes ###
39 
51 
52 
53 namespace shwf {
54 
55  class ShowerReco : public art::EDProducer {
56  public:
57  explicit ShowerReco(fhicl::ParameterSet const& pset);
58 
59  private:
60 
61  void beginJob();
62  void beginRun(art::Run& run);
63  void produce(art::Event& evt);
65  void GetVertexAndAnglesFromCluster(art::Ptr< recob::Cluster > clust,unsigned int plane); // Get shower vertex and slopes.
66  // void GetVertexN(art::Event& evt);
68 
69  void LongTransEnergy(unsigned int set, std::vector< art::Ptr<recob::Hit> > hitlist, bool isData=false); //Longtudinal
70 
71 
72 
73  void ClearandResizeVectors(unsigned int nPlanes);
74 
75 
77  // calo::CalorimetryAlg calorim;
78 //input labels:
79  float slope[3]; // in cm, cm
80  float angle[3]; // in radians
81 
82 
83  std::string fClusterModuleLabel;
84 
85  float ftimetick; // time sample in us
86  double xyz_vertex[3];
87 
88 
90  double fMean_wire_pitch ; // wire pitch in cm
92 
93  std::vector<double> fRMS_2cm;
94  std::vector<int> fNpoints_2cm;
95  std::vector<double> fCorr_MeV_2cm;
96  std::vector<double> fCorr_Charge_2cm;
97 
98  std::vector<int> fNpoints_corr_ADC_2cm;
99  std::vector<int> fNpoints_corr_MeV_2cm;
100 
101  std::vector<double> fTotChargeADC; //Total charge in ADC/cm for each plane
102  std::vector<double> fTotChargeMeV; //Total charge in MeV/cm for each plane
103  std::vector<double> fTotChargeMeV_MIPs; //Total charge in MeV/cm for each plane
104 
105  std::vector<double> fChargeADC_2cm; //Initial charge in ADC/cm for each plane first 4cm
106  std::vector<double> fChargeMeV_2cm; //initial charge in MeV/cm for each plane first 4cm
107 
108  std::vector<double> fChargeMeV_2cm_refined;
109  std::vector<double> fChargeMeV_2cm_axsum;
110 
111  std::vector<std::vector<double> > fDistribChargeADC; //vector with the first De/Dx points ADC
112  std::vector<std::vector<double> > fDistribChargeMeV; //vector with the first De/Dx points converted energy
113  std::vector<std::vector<double> > fDistribHalfChargeMeV;
114  std::vector<std::vector<double> > fDistribChargeposition; //vector with the first De/Dx points' positions
115 
116  std::vector<std::vector<double> > fSingleEvtAngle; //vector with the first De/Dx points
117  std::vector<std::vector<double> > fSingleEvtAngleVal; //vector with the first De/Dx points
118 
119 
120  std::vector<unsigned int> fWire_vertex; // wire coordinate of vertex for each plane
121  std::vector<double> fTime_vertex; // time coordinate of vertex for each plane
122 
123  std::vector<double> fWire_vertexError; // wire coordinate of vertex for each plane
124  std::vector<double> fTime_vertexError; // time coordinate of vertex for each plane
125 
126  std::vector<unsigned int> fWire_last; // wire coordinate of last point for each plane
127  std::vector<double> fTime_last; // time coordinate of last point for each plane
128 
129 
130 
131 
132  //reconstructed 3D position of shower start
133  std::vector<double> xyz_vertex_fit;
134 
135  std::vector< std::vector<double> > fNPitch; // double array, to use each plane for each set of angles
136 
137  //calorimetry variables
138  float Kin_En;
139  std::vector<float> vdEdx;
140  std::vector<float> vresRange;
141  std::vector<float> vdQdx;
142  std::vector<float> deadwire; //residual range for dead wires
143  float Trk_Length;
144  float fTrkPitchC;
145  float fdEdxlength; //distance that gets used to determine e/gamma separation
146  float fcalodEdxlength; // cutoff distance for hits saved to the calo object.
147  bool fUseArea;
148 
149  double xphi,xtheta; // new calculated angles.
150  unsigned int fTPC; //tpc type
151  unsigned int fNPlanes; // number of planes
152  unsigned int fNAngles;
153  TTree* ftree_shwf;
154 
155  //conversion and useful constants
156  double fWirePitch; //wire pitch in cm
157  double fTimeTick;
160 
161  std::vector< int > fNhitsperplane;
162  std::vector< double > fTotADCperplane;
163  //services
164 
166  const detinfo::DetectorProperties* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
167 
168 // //temporary
169 // int mcpdg;
170 // double mcenergy;
171 // double mcphi;
172 // double mctheta;
173 // std::vector< unsigned int> mcwirevertex; // wire coordinate of vertex for each plane
174 // std::vector< double> mctimevertex; // time coordinate of vertex for each plane
175 
176 
177 
178 
179  }; // class ShowerReco
180 
181 
182 
183 //------------------------------------------------------------------------------
185  EDProducer{pset}
186 {
187  fClusterModuleLabel = pset.get< std::string >("ClusterModuleLabel");
188 // fVertexCLusterModuleLabel=pset.get<std::string > ("VertexClusterModuleLabel");
189  fCaloPSet=pset.get< fhicl::ParameterSet >("CaloAlg");
190 
191  fdEdxlength= pset.get< double >("dEdxlength"); //distance that gets used to determine e/gamma separation
192  fcalodEdxlength= pset.get< double >("calodEdxlength"); // cutoff distance for hits saved to the calo object.
193  fUseArea= pset.get< bool >("UseArea");
194 
195  produces< std::vector<recob::Shower> >();
196  produces< art::Assns<recob::Shower, recob::Cluster> >();
197  produces< art::Assns<recob::Shower, recob::Hit> >();
198  produces< std::vector<anab::Calorimetry> >();
199  produces< art::Assns<recob::Shower, anab::Calorimetry> >();
200 }
201 
202 // struct SortByWire
203 // {
204 // bool operator() (art::Ptr<recob::Hit> const& h1, art::Ptr<recob::Hit> const& h2) const
205 // {
206 // return (h1->Wire()->RawDigit()->Channel() < h2->Wire()->RawDigit()->Channel());
207 // }
208 // };
209 
210 // ***************** //
212 {
213 
214 
215 
218 
219  detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
220 
223  fNPlanes = geo->Nplanes();
224  fMean_wire_pitch = geo->WirePitch(); //wire pitch in cm
225 
228 
229  // auto const* detp = lar::providerFrom<detinfo::DetectorPropertiesService>();
230  ftimetick=detprop->SamplingRate()/1000.;
231 
232 
233 
234 
235 
236 
237  ftree_shwf =tfs->make<TTree>("ShowerReco","Results");
239  // ftree_shwf->Branch("ftheta_Mean","std::vector<double>",&fTheta_Mean);
240  // ftree_shwf->Branch("ftheta_RMS","std::vector<double>",&fTheta_RMS );
241 
242  ftree_shwf->Branch("run",&fRun,"run/I");
243  ftree_shwf->Branch("subrun",&fSubRun,"subrun/I");
244  ftree_shwf->Branch("event",&fEvent,"event/I");
245  ftree_shwf->Branch("nplanes",&fNPlanes,"nplanes/I");
246  ftree_shwf->Branch("nangles",&fNAngles,"nangles/I");
247 
248 // ftree_shwf->Branch("fthetaN","std::vector<double>",&fThetaN_ang);
249 // ftree_shwf->Branch("fphiN","std::vector<double>",&fPhiN_ang);
250 
251  ftree_shwf->Branch("xtheta",&xtheta,"xtheta/D");
252  ftree_shwf->Branch("xphi",&xphi,"xphi/D");
253 
254  ftree_shwf->Branch("ftotChargeADC","std::vector<double>",&fTotChargeADC);
255  ftree_shwf->Branch("ftotChargeMeV","std::vector<double>",&fTotChargeMeV);
256  ftree_shwf->Branch("fTotChargeMeV_MIPs","std::vector<double>",&fTotChargeMeV_MIPs);
257 
258  ftree_shwf->Branch("NPitch","std::vector< std::vector<double> >", &fNPitch);
259  // ftree_shwf->Branch("Pitch","std::vector<double>", &fPitch);
260 
261  // this should be temporary - until the omega is sorted out.
262  ftree_shwf->Branch("RMS_2cm","std::vector<double>",&fRMS_2cm);
263  ftree_shwf->Branch("Npoints_2cm","std::vector<int>",&fNpoints_2cm);
264 // ftree_shwf->Branch("RMS_4cm","std::vector<double>",&fRMS_4cm);
265 // ftree_shwf->Branch("Npoints_4cm","std::vector<int>",&fNpoints_4cm);
266 
267  ftree_shwf->Branch("ChargeADC_2cm","std::vector<double>",&fChargeADC_2cm);
268  ftree_shwf->Branch("ChargeMeV_2cm","std::vector<double>",&fChargeMeV_2cm);
269 // ftree_shwf->Branch("ChargeADC_4cm","std::vector<double>",&fChargeADC_4cm);
270 // ftree_shwf->Branch("ChargeMeV_4cm","std::vector<double>",&fChargeMeV_4cm);
271 
272  ftree_shwf->Branch("ChargeMeV_2cm_refined","std::vector<double>",&fChargeMeV_2cm_refined);
273 // ftree_shwf->Branch("ChargeMeV_4cm_refined","std::vector<double>",&fChargeMeV_4cm_refined);
274 
275  ftree_shwf->Branch("ChargeMeV_2cm_axsum","std::vector<double>",&fChargeMeV_2cm_axsum);
276 // ftree_shwf->Branch("ChargeMeV_4cm_axsum","std::vector<double>",&fChargeMeV_4cm_axsum);
277 
278 
279  ftree_shwf->Branch("fNhitsperplane","std::vector<int>",&fNhitsperplane);
280  ftree_shwf->Branch("fTotADCperplane","std::vector<double>",&fTotADCperplane);
281 
282 
283 
284  ftree_shwf->Branch("ChargedistributionADC","std::vector<std::vector<double>>",&fDistribChargeADC);
285 
286  ftree_shwf->Branch("ChargedistributionMeV","std::vector<std::vector<double>>",&fDistribChargeMeV);
287  ftree_shwf->Branch("DistribHalfChargeMeV","std::vector<std::vector<double>>",&fDistribHalfChargeMeV);
288  ftree_shwf->Branch("ChargedistributionPosition","std::vector<std::vector<double>>",&fDistribChargeposition);
289  ftree_shwf->Branch("xyz_vertex_fit","std::vector<double>", &xyz_vertex_fit);
290 
291 }
292 
293 
295 {
296  detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
297 
298  fWirePitch = geom->WirePitch(); //wire pitch in cm
299  fTimeTick=detprop->SamplingRate()/1000.;
302 
303 }
304 
305 
306 
307 
308  void ShowerReco::ShowerReco::ClearandResizeVectors(unsigned int /*nPlanes*/) {
309  //calculate factorial for number of angles
310  int fact=1;
311  for (unsigned int i = 1; i <= fNPlanes; ++i) fact *= i;
312 
313  fNAngles=fact/2;
314 
315  fDistribChargeADC.clear();
316  fDistribChargeMeV.clear();
317  fDistribChargeposition.clear();
318 
319 
320  fDistribChargeADC.resize(fNPlanes);
321  fDistribChargeMeV.resize(fNPlanes);
323 
324 
325  fNPitch.clear();
326  fDistribChargeADC.clear();
327  fDistribChargeMeV.clear();
328  fDistribHalfChargeMeV.clear();
329  fDistribChargeposition.clear();
330 
331  fNPitch.resize(fNAngles);
332  fDistribChargeADC.resize(fNPlanes);
333  fDistribChargeMeV.resize(fNPlanes);
336 
337  for(unsigned int ii=0;ii<fNAngles;ii++) {
338  fNPitch[ii].resize(fNPlanes,-1);
339  }
340 
341  fNpoints_corr_ADC_2cm.clear();
342  fNpoints_corr_MeV_2cm.clear();
343 
344  fNpoints_corr_ADC_2cm.resize(fNAngles,-1);
345  fNpoints_corr_MeV_2cm.resize(fNAngles,-1);
346 // fNpoints_corr_ADC_4cm.resize(fNAngles,-1);
347 // fNpoints_corr_MeV_4cm.resize(fNAngles,-1);
348 
349 
350 
351  for(unsigned int ii=0;ii<fNPlanes;ii++)
352  { fDistribChargeADC[ii].resize(0); //vector with the first De/Dx points
353  fDistribChargeMeV[ii].resize(0); //vector with the first De/Dx points
354  fDistribHalfChargeMeV[ii].resize(0);
355  fDistribChargeposition[ii].resize(0); //vector with the first De/Dx points' positions
356  }
357 
358 
359  fWire_vertex.clear();
360  fTime_vertex.clear();
361  fWire_vertexError.clear();
362  fTime_vertexError.clear();
363  fWire_last.clear();
364  fTime_last.clear();
365  fTotChargeADC.clear();
366  fTotChargeMeV.clear();
367  fTotChargeMeV_MIPs.clear();
368  fRMS_2cm.clear();
369  fNpoints_2cm.clear();
370 
371  fNhitsperplane.clear();
372  fTotADCperplane.clear();
373 
374 
375 
376  fWire_vertex.resize(fNAngles,-1);
377  fTime_vertex.resize(fNAngles,-1);
378  fWire_vertexError.resize(fNPlanes,-1);
379  fTime_vertexError.resize(fNPlanes,-1);
380  fWire_last.resize(fNAngles,-1);
381  fTime_last.resize(fNAngles,-1);
382  fTotChargeADC.resize(fNAngles,0);
383  fTotChargeMeV.resize(fNAngles,0);
384  fTotChargeMeV_MIPs.resize(fNAngles,0);
385  fRMS_2cm.resize(fNAngles,0);
386  fNpoints_2cm.resize(fNAngles,0);
387 
388  fCorr_MeV_2cm.clear();
389  fCorr_Charge_2cm.clear();
390  xyz_vertex_fit.clear();
391 
392  fChargeADC_2cm.clear();
393  fChargeMeV_2cm.clear();
394  fChargeMeV_2cm_refined.clear();
395  fChargeMeV_2cm_refined.clear();
396  fChargeMeV_2cm_axsum.clear();
397 
398 
399  fCorr_MeV_2cm.resize(fNAngles,0);
400  fCorr_Charge_2cm.resize(fNAngles,0);
401  xyz_vertex_fit.resize(3);
402 
403  fChargeADC_2cm.resize(fNAngles,0); //Initial charge in ADC/cm for each plane angle calculation 4cm
404  fChargeMeV_2cm.resize(fNAngles,0); //initial charge in MeV/cm for each angle calculation first 4cm
405  fChargeMeV_2cm_refined.resize(0);
406  fChargeMeV_2cm_refined.resize(fNAngles,0);;
407  fChargeMeV_2cm_axsum.resize(fNAngles,0);;
408 
409  fNhitsperplane.resize(fNPlanes,-1);
410  fTotADCperplane.resize(fNPlanes,-1);
411 
412 
413  vdEdx.clear();
414  vresRange.clear();
415  vdQdx.clear();
416 
417 
418  }
419 
420 
421 
422 
423 
424 // ***************** //
426 {
427 
428  detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
429 
431  fNPlanes = geom->Nplanes();
432  //fdriftvelocity = detprop->DriftVelocity(Efield_SI,Temperature);
433  std::unique_ptr<std::vector<recob::Shower> > Shower3DVector(new std::vector<recob::Shower>);
434  std::unique_ptr< art::Assns<recob::Shower, recob::Cluster> > cassn(new art::Assns<recob::Shower, recob::Cluster>);
435  std::unique_ptr< art::Assns<recob::Shower, recob::Hit> > hassn(new art::Assns<recob::Shower, recob::Hit>);
436  std::unique_ptr< std::vector<anab::Calorimetry> > calorimetrycol(new std::vector<anab::Calorimetry>);
437  std::unique_ptr< art::Assns< anab::Calorimetry,recob::Shower> > calassn(new art::Assns<anab::Calorimetry,recob::Shower>);
438 
439 
445  art::Handle< std::vector<recob::Cluster> > clusterListHandle;
446  evt.getByLabel(fClusterModuleLabel,clusterListHandle);
447 
449  evt.getByLabel(fClusterModuleLabel,clusterAssociationHandle);
450 
451 
452  art::FindManyP<recob::Hit> fmh(clusterListHandle, evt,fClusterModuleLabel);
453 
454 // std::cout << " ---------- ShowerReco!!! -------------- " << std::endl;
455  fRun = evt.id().run();
456  fSubRun = evt.id().subRun();
457  fEvent = evt.id().event();
458 // unsigned int nCollections= clusterAssociationHandle->size();
459 // std::vector < art::PtrVector<recob::Cluster> >::const_iterator clusterSet = clusterAssociationHandle->begin();
460 // for(unsigned int iCol=0;iCol<nCollections;iCol++)
461 // {
462 // const art::PtrVector<recob::Cluster> pvcluster(*(clusterSet++));
463 // auto it = pvcluster.begin();
464 // int nClusts = pvcluster.size();
465 // ClearandResizeVectors(nClusts); // need to do this smarter (coutn at start and have an overall index?)
466 // for(int iClust = 0; iClust < nClusts; ++iClust) {
467 // const art::Ptr<recob::Cluster> pclust(*(it++));
468 // auto pcoll { pclust };
469 // art::FindManyP<recob::Hit> fs( pcoll, evt, fClusterModuleLabel);
470 // std::vector< art::Ptr<recob::Hit> > hitlist = fs.at(0);
471 // //std::cout << " hitlist size for coll " << iCol << " clust " << iClust << " " << hitlist.size() << std::endl;
472 // }
473 // }
474 
475 
476  // find all the hits associated to all the clusters (once and for all);
477  // the index of the query matches the index of the cluster in the collection
478  // (conveniently carried around in its art pointer)
479  art::FindManyP<recob::Hit> ClusterHits(clusterListHandle, evt, fClusterModuleLabel);
480 
481  std::vector < art::PtrVector<recob::Cluster> >::const_iterator clusterSet = clusterAssociationHandle->begin();
482  // loop over vector of vectors (each size of NPlanes) and reconstruct showers from each of those
483  for(size_t iClustSet = 0;iClustSet < clusterAssociationHandle->size(); iClustSet++){
484 
485  const art::PtrVector<recob::Cluster> CurrentClusters=(*(clusterSet++));
486 
487  // do some error checking - i.e. are the clusters themselves present.
488  if(clusterListHandle->size() < 2 || CurrentClusters.size() < 2) {
489  //std::cout << "not enough clusters to reconstruct" << std::endl;
490  //std::cout << "emergency filling tree @ run, evt," << fRun << " " << fEvent << std::endl;
491  ftree_shwf->Fill();
492  return;
493  }
494  //std::cout << " Cluster Set: " << iClustSet << " " << std::endl;
495 
497 
498  std::vector< std::vector< art::Ptr<recob::Hit> > > hitlist_all;
499  hitlist_all.resize(fNPlanes);
500 
501  // int nClusts = CurrentClusters.size();
502  for(size_t iClust = 0; iClust < CurrentClusters.size(); iClust++){
503  art::Ptr<recob::Cluster> const& pclust = CurrentClusters[iClust];
504  //size_t ii=0;
505  //std::cout << " clusterListHandle " << clusterListHandle->size() << " fNPlanes " << fNPlanes << " "<< CurrentClusters.size() << std::endl;
506 
507  // get all the hits for this cluster;
508  // pclust is a art::Ptr to the original cluster collection stored in the event;
509  // its key corresponds to its index in the collection
510  // (and therefore to the index in the query)
511  std::vector< art::Ptr<recob::Hit> > const& hitlist = ClusterHits.at(pclust.key());
512  //std::cout << " hitlist size for coll " << iClustSet << " clust " << iClust << " " << hitlist.size() << std::endl;
513 
514  unsigned int p(0); //c=channel, p=plane, w=wire
515 
516  //std::cout << " hitlist size " << hitlist.size() << std::endl;
517  if(hitlist.size() == 0) continue;
518 
519  p= (*hitlist.begin())->WireID().Plane;
520  //art::Ptr<recob::Cluster> cl(clusterListHandle, iClust);
521  //get vertex position and slope information to start with - ii is the posistion of the correct cluster:
523 
524 
525  double ADCcharge=0;
526  //loop over cluster hits
527  for(art::Ptr<recob::Hit> const& hit: hitlist){
528  p= hit->WireID().Plane;
529  hitlist_all[p].push_back(hit);
530  ADCcharge+= hit->PeakAmplitude();
531  }
532  fNhitsperplane[p]=hitlist_all[p].size();
533  fTotADCperplane[p]=ADCcharge;
534  // unsigned int nCollections= clusterAssociationHandle->size();
535 // std::vector < art::PtrVector<recob::Cluster> >::const_iterator clusterSet = clusterAssociationHandle->begin();
536 // for(unsigned int iCol=0;iCol<nCollections;iCol++)
537 // {
538 // const art::PtrVector<recob::Cluster> pvcluster(*(clusterSet++));
539 // auto it = pvcluster.begin();
540 // int nClusts = pvcluster.size();
541 // ClearandResizeVectors(nClusts); // need to do this smarter (coutn at start and have an overall index?)
542 // for(int iClust = 0; iClust < nClusts; ++iClust) {
543 // const art::Ptr<recob::Cluster> pclust(*(it++));
544 // auto pcoll { pclust };
545 // art::FindManyP<recob::Hit> fs( pcoll, evt, fClusterModuleLabel);
546 // std::vector< art::Ptr<recob::Hit> > hitlist = fs.at(0);
547 // //std::cout << " hitlist size for coll " << iCol << " clust " << iClust << " " << hitlist.size() << std::endl;
548 // }
549 // }
550 
551 
552 
553 
554 
555  // Find the right cluster from the standard cluster list -> to get the hitlist associations.
556 // for( ii = 0; ii < clusterListHandle->size(); ii++){
557 // art::Ptr<recob::Cluster> cl(clusterListHandle, ii);
558 // if((*cl).ID() == (*CurrentClusters[iClust]).ID() ) //find the right cluster out of the list of associated clusters
559 // break;
560 // }
561 //
562  //get vertex position and slope information to start with - ii is the posistion of the correct cluster:
563  // std::vector< art::Ptr<recob::Hit> > hitlist = fmh.at(ii);
564  // std::sort(hitlist.begin(), hitlist.end(), SortByWire());
565 
566 
567  } // End loop on clusters.
568  // Now I have the Hitlists and the relevent clusters parameters saved.
569 
570 
571 // for(unsigned int i = 0; i < fNPlanes; ++i){
572 // std::sort(hitlist_all[i].begin(), hitlist_all[i].end(),SortByWire());
573 //
574 // }
575 
576  //find best set:
577  unsigned int bp1 = 0,bp2 = 0;
578  double minerror1=99999999,minerror2=9999999;
579  for(unsigned int ii = 0; ii < fNPlanes; ++ii)
580  {
581  double locerror=fWire_vertexError[ii]*fWire_vertexError[ii]+ fTime_vertexError[ii]*fTime_vertexError[ii]; // time coordinate of vertex for each plane
582 
583  if(minerror1 >= locerror ) // the >= sign is to favorize collection
584  {
585  minerror1=locerror;
586  bp1=ii;
587  }
588  }
589  for(unsigned int ij = 0; ij < fNPlanes; ++ij)
590  {
591  double locerror=fWire_vertexError[ij]*fWire_vertexError[ij]+ fTime_vertexError[ij]*fTime_vertexError[ij]; // time coordinate of vertex for each plane
592 
593  if(minerror2 >= locerror && ij!=bp1 )
594  {
595  minerror2=locerror;
596  bp2=ij;
597  }
598  }
599  //bp1=0;
600  //bp1=2;
601  //std::cout << " best planes " << bp1 << " " << bp2 << std::endl;
602 for(unsigned int ij = 0; ij < fNPlanes; ++ij)
603  {
604  //std::cout << " wire distances: " << ij << " " << fabs((double)fWire_vertex[ij]-(double)fWire_last[ij]);
605 
606  }
607 
608 
609  //std::cout << " angles in: " << bp1 << " " << bp2 << " " << slope[bp1]*TMath::Pi()/180. << " " << slope[bp2]*TMath::Pi()/180. << " " << slope[bp1] << " " << slope[bp2] << std::endl;
610 // gser.Get3DaxisN(bp1,bp2,slope[bp1]*TMath::Pi()/180.,slope[bp2]*TMath::Pi()/180.,xphi,xtheta);
611  gser.Get3DaxisN(bp1,bp2,angle[bp1],angle[bp2],xphi,xtheta);
612 
614  const double origin[3] = {0.};
615  std::vector <std::vector < double > > position;
616  // auto const* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
617  double fTimeTick=detprop->SamplingRate()/1000.;
619  // get starting positions for all planes
620  for(unsigned int xx=0;xx<fNPlanes;xx++){
621  double pos1[3];
622  geom->Plane(xx).LocalToWorld(origin, pos1);
623  std::vector <double > pos2;
624  pos2.push_back(pos1[0]);
625  pos2.push_back(pos1[1]);
626  pos2.push_back(pos1[2]);
627  position.push_back(pos2);
628  }
629  // Assuming there is no problem ( and we found the best pair that comes close in time )
630  // we try to get the Y and Z coordinates for the start of the shower.
631  try{
632  int chan1=geom->PlaneWireToChannel(bp1,fWire_vertex[bp1], 0);
633  int chan2=geom->PlaneWireToChannel(bp2,fWire_vertex[bp2], 0);
634 
635  double y,z;
636 // bool wires_cross = geom->ChannelsIntersect(chan1,chan2,y,z);
637  geom->ChannelsIntersect(chan1,chan2,y,z);
638  // geom->ChannelsIntersect(chan1,chan2,y,z);
639 
640  xyz_vertex_fit[1]=y;
641  xyz_vertex_fit[2]=z;
642  xyz_vertex_fit[0]=(fTime_vertex[bp1]-detprop->TriggerOffset()) *fDriftVelocity*fTimeTick+position[0][0];
643 
644 
645  //std::cout << ":::::: found x,y,z vertex " << wires_cross << " " << xyz_vertex_fit[0] << " " << y << " " << z << " " << wires_cross << std::endl;
646  }
647  catch(cet::exception e) {
648  mf::LogWarning("ShowerReco") << "caught exception \n" << e;
649  xyz_vertex_fit[1]=0;
650  xyz_vertex_fit[2]=0;
651  xyz_vertex_fit[0]=0;
652  }
653 
654 
655  // if collection is not best plane, project starting point from that
656  if(bp1!=fNPlanes-1 && bp2!=fNPlanes-1)
657  {
658  double pos[3];
659  unsigned int wirevertex;
660 
661  geom->Plane(fNPlanes-1).LocalToWorld(origin, pos);
662  //planex[p] = pos[0];
663  //std::cout << "plane X positionp " << 2 << " " << pos[0] << std::endl;
664 
665  pos[1]=xyz_vertex_fit[1];
666  pos[2]=xyz_vertex_fit[2];
667  wirevertex = geom->NearestWire(pos,fNPlanes-1);
668  //geo->ChannelToWire(channel2,cs,t,p,wirevertex); //\fixme!
669  //wirevertex= (*a)->WireID().Wire;
670 
671 
672 
673 
674  double drifttick=(xyz_vertex_fit[0]/detprop->DriftVelocity(detprop->Efield(),detprop->Temperature()))*(1./fTimeTick);
675  fWire_vertex[fNPlanes-1]= wirevertex; // wire coordinate of vertex for each plane
676  fTime_vertex[fNPlanes-1] = drifttick-(pos[0]/detprop->DriftVelocity(detprop->Efield(),detprop->Temperature()))*(1./fTimeTick)+detprop->TriggerOffset();
677 
678 
679  }
680 
681 
682  // std::cout << "^^^^^^cross-check xphi and xtheta: " << xphi << " " << xtheta << std::endl;
683 
684 
685  if(fabs(xphi) < 5. )
686  { xtheta= gser.Get3DSpecialCaseTheta(bp1,bp2,fWire_last[bp1]-fWire_vertex[bp1], fWire_last[bp2]-fWire_vertex[bp2]);
687 
688  //std::cout << "xphi, xtheta1:,alt" << xphi << " " << xtheta <<std::endl;
689  }
690  //}
691 
692 // zero the arrays just to make sure
693  for(unsigned int i = 0; i < fNAngles; ++i){
694  fTotChargeADC[i]=0;
695  fTotChargeMeV[i]=0;
696  fTotChargeMeV_MIPs[i]=0;
699 
700  fRMS_2cm[i]=0;
701  fNpoints_2cm[i]=0;
702 
703  fCorr_MeV_2cm[i]=0;
704  fCorr_Charge_2cm[i]=0;
705 
706  fChargeADC_2cm[i]=0; //Initial charge in ADC/cm for each plane angle calculation 4cm
707  fChargeMeV_2cm[i]=0; //initial charge in MeV/cm for each angle calculation first 4cm
708 
709  }
710 
711 
712  // do loop and choose Collection. With ful calorimetry can do all.
713  //for(unsigned int set=0;set<fNAngles;set++)
714  if(!(fabs(xphi) >89 && fabs(xphi)<91)) // do not calculate pitch for extreme angles
715  LongTransEnergy(0,hitlist_all[fNPlanes-1]); //temporary only plane 2.
716 
717 
719 
720  //std::vector< recob::SpacePoint > spacepoints = std::vector<recob::SpacePoint>()
721 
722 
723 
724  // make an art::PtrVector of the clusters
726  for(unsigned int i = 0; i < clusterListHandle->size(); ++i){
727  art::Ptr<recob::Cluster> prod(clusterListHandle, i);
728  prodvec.push_back(prod);
729  }
730 
731  //create a singleSpacePoint at vertex.
732  std::vector< recob::SpacePoint > spcpts;
733 
734  //get direction cosines and set them for the shower
735  // TBD determine which angle to use for the actual shower
736  double fPhi=xphi;
737  double fTheta=xtheta;
738 
739  TVector3 dcosVtx(TMath::Cos(fPhi*TMath::Pi()/180)*TMath::Sin(fTheta*TMath::Pi()/180),
740  TMath::Cos(fTheta*TMath::Pi()/180),
741  TMath::Sin(fPhi*TMath::Pi()/180)*TMath::Sin(fTheta*TMath::Pi()/180));
743  // fill with bogus values for now
744  TVector3 dcosVtxErr(util::kBogusD, util::kBogusD, util::kBogusD);
745  //double maxTransWidth[2] = { util::kBogusD };
746  //double distMaxWidth = util::kBogusD;
747  //recob::Shower singShower(dcosVtx, dcosVtxErr, maxTransWidth, distMaxWidth, 1);
748  recob::Shower singShower;
749  singShower.set_direction(dcosVtx);
750  singShower.set_direction_err(dcosVtxErr);
751 
752  Shower3DVector->push_back(singShower);
753  // associate the shower with its clusters
754  util::CreateAssn(*this, evt, *Shower3DVector, prodvec, *cassn);
755 
756  // get the hits associated with each cluster and associate those with the shower
757  for(size_t p = 0; p < prodvec.size(); ++p){
758  std::vector< art::Ptr<recob::Hit> > hits = fmh.at(p);
759  util::CreateAssn(*this, evt, *Shower3DVector, hits, *hassn);
760  }
761 
762 
763  geo::PlaneID planeID(0,0,fNPlanes-1);
764  calorimetrycol->push_back(anab::Calorimetry(Kin_En,
765  vdEdx,
766  vdQdx,
767  vresRange,
768  deadwire,
769  Trk_Length,
770  fTrkPitchC,
771  planeID));
772 
774 
775  //for(unsigned int ip=0;ip<1;ip++) {
776  art::ProductID aid = evt.getProductID< std::vector < recob::Shower > >();
777  art::Ptr< recob::Shower > aptr(aid, 0, evt.productGetter(aid));
778  ssvec.push_back(aptr);
779  //}
780 
781 
782  //util::CreateAssn(*this, evt, *Shower3DVector, calorimetrycol, *calassn);
783  util::CreateAssn(*this, evt, *calorimetrycol,ssvec,*calassn);
785  //std::cout << " filling tree @ run, evt," << fRun << " " << fEvent << std::endl;
786  ftree_shwf->Fill();
787 
788  //for(unsigned int iplane = 0; iplane < fNPlanes; ++iplane)
789  //fh_theta[iplane]->Write(Form("fh_theta_%d_%d",iplane,evt.id().event()));
790  // This needs work, clearly.
791  //for(int p=0;p<2;p++)Shower3DVector->push_back(shower);
792 
793  } // end loop on Vectors of "Associated clusters"
794 
795  evt.put(std::move(Shower3DVector));
796  evt.put(std::move(cassn));
797  evt.put(std::move(hassn));
798  evt.put(std::move(calorimetrycol));
799  evt.put(std::move(calassn));
800 
801 
802 
803 }
804 
805 
806 
807 //------------------------------------------------------------------------------
808 void ShowerReco::LongTransEnergy(unsigned int set, std::vector < art::Ptr<recob::Hit> > hitlist, bool /*isData*/)
809 {
810  // alogorithm for energy vs dx of the shower (roto-translation) COLLECTION VIEW
811  // double wire_cm, time_cm;
812  // int loop_nrg = 0;
813 
816 
817 
818 
819  double totCnrg = 0,totCnrg_corr =0;//, totNewCnrg=0 ; // tot enegry of the shower in collection
820 // art::ServiceHandle<geo::Geometry const> geom;
821 // auto const* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
822 
823  double time;
824  unsigned int wire=0,plane=fNPlanes-1;
825 
826 
827  double mevav2cm=0.;
828  double sum=0.;
829  double npoints_calo=0;
830 
831  int direction=-1;
832 
833  //override direction if phi (XZ angle) is less than 90 degrees
834  if(fabs(xphi)<90)
835  direction=1;
836 
837  //variables to check whether a hit is close to the shower axis.
838  double ortdist,linedist;
839  double wire_on_line,time_on_line;
840 
841  //get effective pitch using 3D angles
842  double newpitch=gser.PitchInView(plane,xphi,xtheta);
843 
844 
845 
846  for(art::PtrVector<recob::Hit>::const_iterator hitIter = hitlist.begin(); hitIter != hitlist.end(); hitIter++){
847  art::Ptr<recob::Hit> theHit = (*hitIter);
848  time = theHit->PeakTime() ;
849 
850  wire= theHit->WireID().Wire;
851  plane= theHit->WireID().Plane;
852 
853  double dEdx_new;
854  // double dEdx_MIP;
855 
856  if(fUseArea)
857  { dEdx_new = calalg.dEdx_AREA((*hitIter), newpitch );
858  // dEdx_MIP = calalg.dEdx_AREA_forceMIP((*hitIter), newpitch );
859  }
860  else //this will hopefully go away, once all of the calibration factors are calculated.
861  {
862  dEdx_new = calalg.dEdx_AMP((*hitIter), newpitch );
863  // dEdx_MIP = calalg.dEdx_AMP_forceMIP((*hitIter), newpitch );
864  }
865 
866  //calculate total energy.
867  totCnrg_corr += dEdx_new;
868  //totNewCnrg+=dEdx_MIP;
869 
870  // calculate the wire,time coordinates of the hit projection on to the 2D shower axis
871  gser.GetPointOnLine(slope[plane]/fWireTimetoCmCm,fWire_vertex[plane],fTime_vertex[plane],wire,time,wire_on_line,time_on_line);
872  linedist=gser.Get2DDistance(wire_on_line,time_on_line,fWire_vertex[plane],fTime_vertex[plane]);
873  ortdist=gser.Get2DDistance(wire_on_line,time_on_line,wire,time);
874 
875 
876 
877  //calculate the distance from the vertex using the effective pitch metric
878  double wdist=(((double)wire-(double)fWire_vertex[plane])*newpitch)*direction; //wdist is always positive
879 
880 
881 // if( (fabs(wdist)<fcalodEdxlength)&&(fabs(wdist)>0.2)){
882  if( (wdist<fcalodEdxlength)&&(wdist>0.2)){
883 
884  vdEdx.push_back(dEdx_new);
885  vresRange.push_back(fabs(wdist));
886  vdQdx.push_back((*hitIter)->PeakAmplitude()/newpitch);
887  Trk_Length=wdist;
888  fTrkPitchC=fNPitch[set][plane];
889  Kin_En+=dEdx_new*newpitch;
890  npoints_calo++;
891  sum+=dEdx_new;
892 
893 
894 
895 //std::cout << " CALORIMETRY:" << " Pitch " <<newpitch << " dist: " << wdist << " dE/dx: " << dEdx_new << "MeV/cm " << " average: " << sum/npoints_calo << "hit: wire, time " << wire << " " << time << " line,ort " << linedist << " " << ortdist<< " direction " << direction << std::endl;
896 
897 
898  if(wdist<fdEdxlength
899  && ((direction==1 && wire>fWire_vertex[plane]) //take no hits before vertex (depending on direction)
900  || (direction==-1 && wire<fWire_vertex[plane]) )
901  && ortdist<4.5 && linedist < fdEdxlength ){
902  fChargeMeV_2cm[set]+= dEdx_new ;
903  fNpoints_2cm[set]++;
904  // std::cout << " CALORIMETRY:" << " Pitch " <<newpitch << " dist: " << wdist << " dE/dx: " << dEdx_new << "MeV/cm " << " average: " << sum/npoints_calo << "hit: wire, time " << wire << " " << time << " line,ort " << linedist << " " << ortdist<< " direction " << direction << std::endl;
905  }
906 
907  // fill out for 4cm preshower
908 
909  //fDistribChargeADC[set].push_back(ch_adc); //vector with the first De/Dx points
910  fDistribChargeMeV[set].push_back(dEdx_new); //vector with the first De/Dx points
911  //fDistribHalfChargeMeV[set].push_back(Bcorr_half);
912  fDistribChargeposition[set].push_back(wdist); //vector with the first De/Dx points' positions
913 
914  }//end inside range if statement
915 
916  }// end first loop on hits.
917 
918  auto const signalType
919  = hitlist.empty()? geo::kMysteryType: geom->SignalType(hitlist.front()->WireID());
920 
921  if(signalType == geo::kCollection)
922  {
923  fTotChargeADC[set]=totCnrg*newpitch;
924  fTotChargeMeV[set]=totCnrg_corr*newpitch;
925  //fTotChargeMeV_MIPs[set]=totNewCnrg*newpitch;
926  }
927 
928 
929  //calculate average dE/dx
930  if(fNpoints_2cm[set]>0) {
931  mevav2cm=fChargeMeV_2cm[set]/fNpoints_2cm[set];
932  }
933  //double RMS_ADC_2cm=0.;
934 
935 
936 
937  //second loop to calculate RMS
938  for(art::PtrVector<recob::Hit>::const_iterator hitIter = hitlist.begin(); hitIter != hitlist.end(); hitIter++){
939  art::Ptr<recob::Hit> theHit = (*hitIter);
940  time = theHit->PeakTime() ;
941  wire= theHit->WireID().Wire;
942  plane= theHit->WireID().Plane;
943  double dEdx=0;
944 
945  if(fUseArea)
946  { dEdx = calalg.dEdx_AREA((*hitIter), newpitch );
947  }
948  else //this will hopefully go away, once all of the calibration factors are calculated.
949  {
950  dEdx = calalg.dEdx_AMP((*hitIter), newpitch );
951  }
952 
953 
954 
955  gser.GetPointOnLine(slope[plane]/fWireTimetoCmCm,fWire_vertex[plane],fTime_vertex[plane],wire,time,wire_on_line,time_on_line);
956  linedist=gser.Get2DDistance(wire_on_line,time_on_line,fWire_vertex[plane],fTime_vertex[plane]);
957  ortdist=gser.Get2DDistance(wire_on_line,time_on_line,wire,time);
958 
959 
960 
961  double wdist=(((double)wire-(double)fWire_vertex[plane])*newpitch)*direction;
962 
963  // //std::cout << dEdx << " MeV, outside of if;; wd " << wdist << " ld " << linedist << " od " << ortdist << std::endl;
964 
965  if( (wdist<fcalodEdxlength)&&(wdist>0.2)){
966  if(wdist<fdEdxlength
967  && ((direction==1 && wire>fWire_vertex[plane]) ||
968  (direction==-1 && wire<fWire_vertex[plane]) )
969  && ortdist<4.5 && linedist < fdEdxlength)
970  {
971 // //std::cout << dEdx << " MeV " << std::endl;
972  fRMS_2cm[set]+= (dEdx-mevav2cm)*(dEdx-mevav2cm);
973  }
974 
975 
976  } //end if on correct hits.
977  } //end RMS_calculating loop.
978 
979  if(fNpoints_2cm[set]>0)
980  {
981  fRMS_2cm[set]=TMath::Sqrt(fRMS_2cm[set]/fNpoints_2cm[set]);
982  }
983 
984  //std::cout << " average dE/dx: " << mevav2cm << " RMS:: " << fRMS_2cm[set] << " " << fNpoints_2cm[set] << std::endl;
985 
987 
988  for(art::PtrVector<recob::Hit>::const_iterator hitIter = hitlist.begin(); hitIter != hitlist.end(); hitIter++){
989  art::Ptr<recob::Hit> theHit = (*hitIter);
990  time = theHit->PeakTime() ;
991  wire= theHit->WireID().Wire;
992  plane= theHit->WireID().Plane;
993 
994  double dEdx=0;
995  if(fUseArea)
996  { dEdx = calalg.dEdx_AREA((*hitIter), newpitch );
997  }
998  else //this will hopefully go away, once all of the calibration factors are calculated.
999  {
1000  dEdx = calalg.dEdx_AMP((*hitIter), newpitch );
1001  }
1002 
1003  gser.GetPointOnLine(slope[plane]/fWireTimetoCmCm,fWire_vertex[plane],fTime_vertex[plane],wire,time,wire_on_line,time_on_line);
1004  linedist=gser.Get2DDistance(wire_on_line,time_on_line,fWire_vertex[plane],fTime_vertex[plane]);
1005  ortdist=gser.Get2DDistance(wire_on_line,time_on_line,wire,time);
1006 
1007  double wdist=(((double)wire-(double)fWire_vertex[plane])*newpitch)*direction;
1008 
1009 
1010  if( (wdist < fcalodEdxlength) && (wdist > 0.2
1011  && ((direction==1 && wire>fWire_vertex[plane]) ||
1012  (direction==-1 && wire<fWire_vertex[plane]) )
1013  && ortdist<4.5 && linedist < fdEdxlength ))
1014  {
1015  if(wdist < fdEdxlength)
1016  {
1017  if( ((dEdx > (mevav2cm-fRMS_2cm[set]) )
1018  && (dEdx < (mevav2cm+fRMS_2cm[set]) ))
1019  || (newpitch > 0.3*fdEdxlength ) ) {
1020  fCorr_MeV_2cm[set]+= dEdx;
1021  fNpoints_corr_MeV_2cm[set]++;
1022  }
1023 
1024  } // end if on good hits
1025 
1026 
1027  }
1028  } //end of third loop on hits
1029 
1030  if(fNpoints_corr_MeV_2cm[set]>0){
1031  //std::cout << " ++ NPoints 2cm, ADC and MeV "
1032  // << fNpoints_corr_MeV_2cm[set] << " "
1033  // << fNpoints_corr_ADC_2cm[set]
1034  // << " corrected average De/Dx, charge, MeV: "
1035  // << fCorr_Charge_2cm[set]/fNpoints_corr_ADC_2cm[set]
1036  // << " " << fCorr_MeV_2cm[set]/fNpoints_corr_MeV_2cm[set] << std::endl;
1037  //fCorr_Charge_2cm[set]/=fNpoints_corr_ADC_2cm[set];
1040  }
1041 
1042 
1044 // std::cout << " total ENERGY, birks: " << fTotChargeMeV[set] << " MeV " << " |average: " << fChargeMeV_2cm_refined[set] << std::endl;
1045 }
1046 
1047 
1048 
1049 //------------------------------------------------------------------------------------//
1051 // Get shower vertex and slopes.
1052 {
1053  //convert to cm/cm units needed in the calculation
1054 // slope[plane]=clust->dTdW();//*(ftimetick*fdriftvelocity)/fMean_wire_pitch;
1055  angle[plane]=clust->StartAngle();
1056  slope[plane]=std::tan(clust->StartAngle());
1057  fWire_vertex[plane]=clust->StartWire();
1058  fTime_vertex[plane]=clust->StartTick();
1059 
1060  fWire_vertexError[plane]=clust->SigmaStartWire(); // wire coordinate of vertex for each plane
1061  fTime_vertexError[plane]=clust->SigmaStartTick(); // time coordinate of vertex for each plane
1062 
1063  fWire_last[plane]=clust->EndWire(); // wire coordinate of last point for each plane
1064  fTime_last[plane]=clust->EndTick();
1065 
1067 
1068  // std::cout << "======= setting slope for view: " << plane
1069 // << " " << slope[plane] << " " << fWire_vertex[plane]
1070 // << " " << fTime_vertex[plane] << " " << std::endl;
1071  // << fWire_vertex[plane]+50<< " "
1072  // << fTime_vertex[plane] + slope[plane]*(fWire_vertex[plane]+50)<< std::endl;
1073 
1074 }
1075 
1076 
1077 // //------------------------------------------------------------------------------
1078 // //to be moved to an _ana module
1079 // void ShowerReco::GetVertexN(art::Event& evt){
1080 //
1081 //
1082 //
1083 // double fTimeTick=detprop->SamplingRate()/1000.;
1084 //
1085 // art::Handle< std::vector<simb::MCTruth> > mctruthListHandle;
1086 // evt.getByLabel("generator",mctruthListHandle);
1087 //
1088 // art::PtrVector<simb::MCTruth> mclist;
1089 // for (unsigned int ii = 0; ii < mctruthListHandle->size(); ++ii)
1090 // {
1091 // art::Ptr<simb::MCTruth> mctparticle(mctruthListHandle,ii);
1092 // mclist.push_back(mctparticle);
1093 // }
1094 //
1095 //
1096 // mf::LogVerbatim("ShowerAngleClusterAna") << "%%%%%%% mc size size, "<<mclist.size() << std::endl;
1097 //
1098 //
1099 // art::Ptr<simb::MCTruth> mc(mclist[0]);
1100 // simb::MCParticle neut(mc->GetParticle(0));
1101 //
1102 // mcpdg=neut.PdgCode();
1103 // mcenergy=neut.E();
1104 // //std::cout << " ----------- mcenergy:: " << mcenergy << std::endl;
1105 // if (neut.P()){
1106 // double lep_dcosx_truth = neut.Px()/neut.P();
1107 // double lep_dcosy_truth = neut.Py()/neut.P();
1108 // double lep_dcosz_truth = neut.Pz()/neut.P();
1109 //
1110 // mf::LogVerbatim("ShowerAngleClusterAna") << "----- cx,cy,cz " << lep_dcosx_truth << " " << lep_dcosy_truth << " " << lep_dcosz_truth << std::endl;
1111 //
1112 //
1113 // mcphi= (lep_dcosx_truth == 0.0 && lep_dcosz_truth == 0.0) ? 0.0 : TMath::ATan2(lep_dcosx_truth,lep_dcosz_truth);
1114 // mctheta= (lep_dcosx_truth == 0.0 && lep_dcosy_truth == 0.0 && lep_dcosz_truth == 0.0) ? 0.0 : TMath::Pi()*0.5-TMath::ATan2(std::sqrt(lep_dcosx_truth*lep_dcosx_truth + lep_dcosz_truth*lep_dcosz_truth),lep_dcosy_truth);
1115 //
1116 //
1117 // mcphi=180*mcphi/TMath::Pi();
1118 // mctheta= 180*mctheta/TMath::Pi();
1119 // mf::LogVerbatim("ShowerAngleClusterAna") << "----- phi, theta " << mcphi << " " << mctheta << std::endl;
1120 //
1121 // }
1122 //
1123 // int npart=0;
1124 // // while(&& npart < mc->NParticles() )
1125 // // {
1126 // mf::LogVerbatim("ShowerAngleClusterAna") << "%%%%%%%####### is PDG: "<< npart <<" " << neut.PdgCode() << std::endl;
1127 // // neut=mc->GetParticle(npart++);
1128 //
1129 // // }
1130 //
1131 // mf::LogVerbatim("ShowerAngleClusterAna") << "%%%%%%%####### after loop is PDG: "<< npart <<" " << neut.PdgCode() << std::endl;
1132 // //if((neut.PdgCode()==11 || neut.PdgCode()==-11 )&& neut.StatusCode()==1){
1133 //
1134 //
1135 // xyz_vertex[0] =neut.Vx();
1136 // xyz_vertex[1] =neut.Vy();
1137 // xyz_vertex[2] =neut.Vz();
1138 //
1139 // mf::LogVerbatim("ShowerAngleClusterAna") <<"neut.Vx()= "<<neut.Vx()<<" ,y= "<<neut.Vy()<<" ,z= "<<neut.Vz()<<std::endl;
1140 // //if(((neut.PdgCode()==11 || neut.PdgCode()==-11 )&& neut.StatusCode()==1))
1141 // // break;
1142 //
1143 //
1144 // double drifttick=(xyz_vertex[0]/detprop->DriftVelocity(detprop->Efield(),detprop->Temperature()))*(1./fTimeTick);
1145 //
1146 // const double origin[3] = {0.};
1147 // for(unsigned int iplane=0;iplane<fNPlanes;iplane++)
1148 // {
1149 // double pos[3];
1150 // geom->Plane(iplane).LocalToWorld(origin, pos);
1151 // //planex[p] = pos[0];
1152 // mf::LogVerbatim("ShowerAngleClusterAna") << "plane X positionp " << iplane << " " << pos[0] << std::endl;
1153 //
1154 // pos[1]=xyz_vertex[1];
1155 // pos[2]=xyz_vertex[2];
1156 // ///\todo: have to change to use cryostat and TPC in NearestWire too
1157 // unsigned int wirevertex = geom->NearestWire(pos,iplane);
1158 //
1159 //
1160 // mcwirevertex[iplane]=wirevertex; // wire coordinate of vertex for each plane
1161 // mctimevertex[iplane]=drifttick-(pos[0]/detprop->DriftVelocity(detprop->Efield(),detprop->Temperature()))*(1./fTimeTick)+detprop->TriggerOffset(); // time coordinate of vertex for each plane
1162 //
1163 // //fWireVertex[p]=wirevertex;
1164 // //fTimeVertex[p]=drifttick-(pos[0]/detprop->DriftVelocity(detprop->Efield(),detprop->Temperature()))*(1./fTimeTick)+60;
1165 // mf::LogVerbatim("ShowerAngleClusterAna") << "wirevertex= "<< wirevertex
1166 // << " timevertex " << mctimevertex[iplane]
1167 // << " correction "
1168 // << (pos[0]/detprop->DriftVelocity(detprop->Efield(),
1169 // detprop->Temperature()))*(1./fTimeTick)
1170 // << " " << pos[0];
1171 //
1172 //
1173 // }
1174 //
1175 //
1176 //
1177 // return (void)0;
1178 // }
1179 
1180 
1181 
1182 
1183 
1184 
1185 
1187 
1188 
1189 }
std::vector< double > fTime_last
fhicl::ParameterSet fCaloPSet
std::vector< double > fWire_vertexError
ShowerReco(fhicl::ParameterSet const &pset)
std::vector< double > fTotChargeADC
std::vector< double > fCorr_Charge_2cm
void set_direction_err(const TVector3 &dir_e)
Definition: Shower.h:136
Double_t xx
Definition: macro.C:12
PlaneGeo const & Plane(unsigned int const p, unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified wire.
Who knows?
Definition: geo_types.h:143
ProductID getProductID(std::string const &instance_name="") const
Definition: DataViewImpl.h:349
virtual int TriggerOffset() const =0
Double_t Get2DDistance(Double_t wire1, Double_t time1, Double_t wire2, Double_t time2) const
void Get2DVariables(art::PtrVector< recob::Hit > hitlist)
double dEdx_AREA(art::Ptr< recob::Hit > hit, double pitch, double T0=0) const
geo::WireID WireID() const
Initial tdc tick for hit.
Definition: Hit.h:233
std::vector< int > fNhitsperplane
Declaration of signal hit object.
EDProducer()=default
Float_t y
Definition: compare.C:6
Double_t z
Definition: plot.C:279
The data type to uniquely identify a Plane.
Definition: geo_types.h:468
std::vector< double > fTotADCperplane
std::vector< std::vector< double > > fDistribChargeMeV
float StartWire() const
Returns the wire coordinate of the start of the cluster.
Definition: Cluster.h:286
bool ChannelsIntersect(raw::ChannelID_t c1, raw::ChannelID_t c2, double &y, double &z) const
Returns an intersection point of two channels.
virtual double SamplingRate() const =0
Returns the period of the TPC readout electronics clock.
std::vector< int > fNpoints_corr_MeV_2cm
std::vector< std::vector< double > > fSingleEvtAngle
shower finding
std::vector< double > fTime_vertex
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:576
float EndTick() const
Returns the tick coordinate of the end of the cluster.
Definition: Cluster.h:342
std::vector< unsigned int > fWire_last
SigType_t SignalType(geo::PlaneID const &pid) const
Returns the type of signal on the channels of specified TPC plane.
std::vector< std::vector< double > > fDistribHalfChargeMeV
float StartAngle() const
Returns the starting angle of the cluster.
Definition: Cluster.h:475
std::vector< double > fChargeMeV_2cm_refined
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
std::vector< float > vdEdx
RunNumber_t run() const
Definition: EventID.h:98
art::ServiceHandle< geo::Geometry const > geom
void GetVertexAndAnglesFromCluster(art::Ptr< recob::Cluster > clust, unsigned int plane)
Definition: Run.h:21
Double_t PitchInView(UInt_t plane, Double_t phi, Double_t theta) const
std::vector< float > vdQdx
std::vector< int > fNpoints_2cm
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
std::vector< double > fTime_vertexError
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:437
unsigned int Nplanes(unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wire planes in the specified TPC.
std::vector< std::vector< double > > fSingleEvtAngleVal
void LongTransEnergy(unsigned int set, std::vector< art::Ptr< recob::Hit > > hitlist, bool isData=false)
void hits()
Definition: readHits.C:15
void set_direction(const TVector3 &dir)
Definition: Shower.h:135
typename data_t::const_iterator const_iterator
Definition: PtrVector.h:55
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:56
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:435
std::vector< float > vresRange
EDProductGetter const * productGetter(ProductID const pid) const
virtual double Temperature() const =0
T get(std::string const &key) const
Definition: ParameterSet.h:231
geo::WireID::WireID_t NearestWire(geo::Point_t const &point, geo::PlaneID const &planeid) const
Returns the index of wire closest to position in the specified TPC.
key_type key() const
Definition: Ptr.h:238
std::vector< double > fRMS_2cm
void produce(art::Event &evt)
float SigmaStartWire() const
Returns the uncertainty on wire coordinate of the start of the cluster.
Definition: Cluster.h:306
float dEdx(TCSlice &slc, TP3D &tp3d)
Definition: PFPUtils.cxx:2191
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:489
Declaration of cluster object.
std::vector< std::vector< double > > fNPitch
size_type size() const
Definition: PtrVector.h:302
Detector simulation of raw signals on wires.
std::vector< double > fCorr_MeV_2cm
Int_t GetPointOnLine(Double_t slope, Double_t intercept, Double_t wire1, Double_t time1, Double_t &wireout, Double_t &timeout) const
std::string fClusterModuleLabel
raw::ChannelID_t PlaneWireToChannel(WireID const &wireid) const
Returns the ID of the TPC channel connected to the specified wire.
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:218
std::vector< double > xyz_vertex_fit
void beginRun(art::Run &run)
Utility object to perform functions of association.
bool CreateAssn(PRODUCER const &prod, art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t index=UINT_MAX)
Creates a single one-to-one association.
Encapsulate the construction of a single detector plane.
std::vector< double > fChargeMeV_2cm
std::vector< double > fTotChargeMeV_MIPs
virtual double DriftVelocity(double efield=0., double temperature=0.) const =0
const detinfo::DetectorProperties * detprop
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
float SigmaStartTick() const
Returns the uncertainty on tick coordinate of the start of the cluster.
Definition: Cluster.h:315
std::vector< double > fChargeADC_2cm
std::vector< float > deadwire
std::vector< std::vector< double > > fDistribChargeADC
EventNumber_t event() const
Definition: EventID.h:116
constexpr double kBogusD
obviously bogus double value
std::vector< double > fChargeMeV_2cm_axsum
void ClearandResizeVectors(unsigned int nPlanes)
virtual double Efield(unsigned int planegap=0) const =0
Returns the nominal electric field in the specified volume.
TCEvent evt
Definition: DataStructs.cxx:7
Float_t e
Definition: plot.C:34
std::vector< int > fNpoints_corr_ADC_2cm
ProductID put(std::unique_ptr< PROD > &&edp, FullSemantic< Level::Run > const semantic)
Definition: DataViewImpl.h:692
Namespace collecting geometry-related classes utilities.
std::vector< std::vector< double > > fDistribChargeposition
float StartTick() const
Returns the tick coordinate of the start of the cluster.
Definition: Cluster.h:297
std::vector< unsigned int > fWire_vertex
SubRunNumber_t subRun() const
Definition: EventID.h:110
void LocalToWorld(const double *plane, double *world) const
Transform point from local plane frame to world frame.
Definition: PlaneGeo.h:1221
double dEdx_AMP(art::Ptr< recob::Hit > hit, double pitch, double T0=0) const
EventID id() const
Definition: Event.cc:37
art framework interface to geometry description
constexpr Point origin()
Returns a origin position with a point of the specified type.
Definition: geo_vectors.h:227
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
std::vector< double > fTotChargeMeV
Signal from collection planes.
Definition: geo_types.h:142
float EndWire() const
Returns the wire coordinate of the end of the cluster.
Definition: Cluster.h:329