StopperThreshold_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: CosmicTrends
3 // Module Type: analyzer
4 // File: CosmicTrends_module.cc
5 //
6 // Generated at Tue Aug 21 13:38:36 2012 by Brian_Rebel using artmod
7 // from art v1_00_11.
8 ////////////////////////////////////////////////////////////////////////
9 #ifndef CosmicTrends_h
10 #define CosmicTrends_h
11 
12 #include <map>
13 
14 #include "art/Framework/Core/ModuleMacros.h"
15 #include "canvas/Persistency/Common/FindOne.h"
16 #include "art/Framework/Core/EDAnalyzer.h"
17 #include "art/Framework/Services/Optional/TFileService.h"
18 #include "art/Framework/Services/Optional/TFileDirectory.h"
19 #include "art/Framework/Services/Registry/ServiceHandle.h"
20 #include "art/Framework/Principal/Event.h"
21 #include "art/Framework/Principal/Handle.h"
22 #include "canvas/Persistency/Common/Ptr.h"
23 #include "canvas/Persistency/Common/PtrVector.h"
24 #include "fhiclcpp/ParameterSet.h"
25 #include "messagefacility/MessageLogger/MessageLogger.h"
26 #include "MCCheater/BackTracker.h"
27 
28 #include "Geometry/Geometry.h"
29 #include "GeometryObjects/PlaneGeo.h"
30 #include "GeometryObjects/CellGeo.h"
31 #include "RecoBase/Track.h"
32 #include "RecoBase/Cluster.h"
33 #include "RecoBase/RecoHit.h"
34 #include "RecoBase/CellHit.h"
35 #include "RunHistory/service/RunHistoryService.h"
36 #include "Geometry/LiveGeometry.h"
37 #include "CalibrationUtils/CalibUtil.h"
38 
39 #include "TStopwatch.h"
40 #include "TH1.h"
41 #include "TH2.h"
42 #include "TVector3.h"
43 #include "TTree.h"
44 #include <string>
45 #include <sstream>
46 
47 namespace calib {
48  class StopperThreshold;
49 
50  enum _oneDHists{
51  };
52 
53  enum _twoDHists{
54  };
55 
56 }
57 
58 class calib::StopperThreshold : public art::EDAnalyzer {
59 public:
60  explicit StopperThreshold(fhicl::ParameterSet const & p);
61  virtual ~StopperThreshold();
62 
63  virtual void analyze(art::Event const & e);
64 
65  virtual void beginRun(art::Run const&);
66  virtual void endRun(art::Run const&);
67  virtual void reconfigure(fhicl::ParameterSet const & p);
68 
69 private:
70 
71  void FindTrajectoryPoints(rb::Track const& track, std::set<TVector3, bool(*)(TVector3 const&, TVector3 const&)> & trajPoints);
72 
73  void FillTree(rb::Track const& track, std::set<TVector3, bool(*)(TVector3 const&, TVector3 const&) > const& trajPoints);
74  void FillHist(rb::Track const& track, std::set<TVector3, bool(*)(TVector3 const&, TVector3 const&) > const& trajPoints);
75  void testPath(rb::Track const& track, std::set<TVector3, bool(*)(TVector3 const&, TVector3 const&) > const& trajPoints);
76 
77  int CalculateWall(float x0, float z0, float dx, float dz, float xlow, float xhigh, float zlow, float zhigh);
78 
79  art::ServiceHandle<geo::Geometry> fGeo; ///< handle to the geometry
80 
81  std::string fTrackLabel; ///< label of module creating the tracks used
82  double fMinDistToEdge;
83  double fMinDistFromEnd; ///< minimum distance from end of track to use a point
84  double fMaxDistFromEnd; ///< maximum distance from end of track to use a point
85  //int histfill;
86 
87  // histograms for shadow calculation
88  TH2F* fADCvsYTot[4];
89  TH2F* fADCratevsYTot[4];
90  TH2F* fADCratevsYlongTot[4];
91  //TH2F* fADCvsYmod[4][12];
92  //TH2F* fADCratevsYmod[4][12];
93  TH2* ftrueShadowVsYTot[4];
94 
95  // histograms for threshold calculation
96  TH2F* fADCvsWTot[4];
97  TH2F* fADCratevsWTot[4];
98  TH2F* fADCratevsWlongTot[4];
99  TH2F* fLvsWTot[4];
100  TH2* fLvsWnothreshTot[4];
101  //TH2F* fADCvsWmod[4][12];
102  //TH2F* fADCratevsWmod[4][12];
103  //TH2F* fLvsWmod[4][12];
104  TH2* ftrueThredVsWTot[4];
105  TH2* flambdaVsWTot[4];
106  TH2* fPEVslambdaTot[4];
107  TH2* fADCratevsLTot[4];
108  TH2* fPEvsLTot[4];
109  TH2* fPECorrvsLTot[4];
110  TH1* fADCrateideal[4];
111  TH1* hrecopath[4];
112  TH1* hrecopathTest[4]; //tmp
113  TH1* htruepath[4];
114  TH1* htruepathTest[4]; //tmp
115  TH1* hrecow[4];
116  TH1* htruew[4];
117  TH1* hrecodx[4];
118  TH1* hrecody[4];
119  TH1* hrecodz[4];
120  //TH2* htruepathvstruemev[4];
121  //TH2* htruewvslambda[4];
122  TH2F* fLvsdxTot[4];
123  TH2F* fLvsdyTot[4];
124  TH2F* fLvsdzTot[4];
125  TH2F* fLvswallsTot[4];
126  TH1F* fpathtest[4];
127 
128  //NEW
129  TTree* tInfoTree;
130  int tView;
131  int tCell;
132  int tPlane;
133  float tADC;
134  float tPE;
135  float tPECorr;
136  float tW;
137  float tTruePath;
138  float tPath;
139  float tRecodx;
140  float tRecody;
141  float tRecodz;
142  float tWall;
143  float tYTot;
144  float tPoissonLambda;
145  float tTrueMeV;
146  float tTrueW;
147  int tIsTruePath;
148  int tAngleRestr;
149  int tLongSample;
150  int tLongNear;
151  int tIsStopper;
152 
153  art::ServiceHandle<geo::LiveGeometry> fLivegeo;
154  art::ServiceHandle<cheat::BackTracker> fBT;
155  bool isstopper;
156 
157 };
158 
159 bool cellSortZ(art::Ptr<rb::CellHit> const& a, art::Ptr<rb::CellHit> const& b)
160 {
161  return a->Plane() < b->Plane();
162 }
163 
164 bool ltTVec(TVector3 const& a, TVector3 const& b)
165 {
166  return a.Z() < b.Z();
167 }
168 
169 //---------------------------------------------------------------------
170 calib::StopperThreshold::StopperThreshold(fhicl::ParameterSet const & p)
171  : EDAnalyzer(p)
172 {
173  this->reconfigure(p);
174 }
175 
176 //---------------------------------------------------------------------
177 calib::StopperThreshold::~StopperThreshold()
178 {
179 }
180 
181 //---------------------------------------------------------------------
182 void calib::StopperThreshold::beginRun(art::Run const&)
183 {
184  art::ServiceHandle<art::TFileService> tfs;
185 
186  // Info Tree
187  tInfoTree = tfs->make<TTree>("InfoTree","Shadow and Threshold variables");
188  tInfoTree->Branch("tView",&tView,"tView/I");
189  tInfoTree->Branch("tCell",&tCell,"tCell/I");
190  tInfoTree->Branch("tPlane",&tPlane,"tPlane/I");
191  tInfoTree->Branch("tADC",&tADC,"tADC/F");
192  tInfoTree->Branch("tPE",&tPE,"tPE/F");
193  tInfoTree->Branch("tPECorr",&tPECorr,"tPECorr/F");
194  tInfoTree->Branch("tW",&tW,"tW/F");
195  tInfoTree->Branch("tTruePath",&tTruePath,"tTruePath/F");
196  tInfoTree->Branch("tPath",&tPath,"tPath/F");
197  tInfoTree->Branch("tRecodx",&tRecodx,"tRecodx/F");
198  tInfoTree->Branch("tRecody",&tRecody,"tRecody/F");
199  tInfoTree->Branch("tRecodz",&tRecodz,"tRecodz/F");
200  tInfoTree->Branch("tWall",&tWall,"tWall/F");
201  tInfoTree->Branch("tYTot",&tYTot,"tYTot/F");
202  tInfoTree->Branch("tPoissonLambda",&tPoissonLambda,"tPoissonLambda/F");
203  tInfoTree->Branch("tTrueMeV",&tTrueMeV,"tTrueMeV/F");
204  tInfoTree->Branch("tTrueW",&tTrueW,"tTrueW/F");
205  tInfoTree->Branch("tIsTruePath",&tIsTruePath,"tIsTruePath/I");
206  tInfoTree->Branch("tAngleRestr",&tAngleRestr,"tAngleRestr/I");
207  tInfoTree->Branch("tLongSample",&tLongSample,"tLongSample/I");
208  tInfoTree->Branch("tLongNear",&tLongNear,"tLongNear/I");
209  tInfoTree->Branch("tIsStopper",&tIsStopper,"tIsStopper/I");
210 
211  //unsigned int planeInView = *(fGeo->GetPlanesByView(geo::kX).begin());
212  //double xCells = fGeo->Plane(planeInView)->Ncells();
213  //planeInView = *(fGeo->GetPlanesByView(geo::kY).begin());
214  //double yCells = fGeo->Plane(planeInView)->Ncells();
215  //double maxCells = std::max(yCells, xCells);
216  //int modules = maxCells/32.;
217  int nADCbin = 1000;
218  int maxADC = 1000;
219  float maxADCrate = 200.0;
220 
221  double detHalfWidth = std::max(fGeo->DetHalfWidth(), fGeo->DetHalfHeight());
222  int wBins = detHalfWidth/5; // 10 cm/bin
223 
224  TString sview[4] = {"Xall","Yall","Xmip","Ymip"};
225  TString name = "";
226  for(int i = 0; i < 4; ++i){
227 
228  fADCvsYTot[i] = tfs->make<TH2F>(sview[i]+"ADCvsYTot",";Y ;ADC",
229  wBins,-detHalfWidth,detHalfWidth,nADCbin,0,maxADC);
230 
231  fADCratevsYTot[i] = tfs->make<TH2F>(sview[i]+"ADCratevsYTot",";Y ;ADC/Path",
232  wBins,-detHalfWidth,detHalfWidth,nADCbin,0,maxADCrate);
233 
234  fADCratevsYlongTot[i] = tfs->make<TH2F>(sview[i]+"ADCratevsYlongTot",";Y ;ADC/Path",
235  wBins,-detHalfWidth,detHalfWidth,nADCbin,0,maxADCrate);
236 
237  fADCvsWTot[i] = tfs->make<TH2F>(sview[i]+"ADCvsWTot",";W ;ADC",
238  wBins,-detHalfWidth,detHalfWidth,nADCbin,0,maxADC);
239 
240  fADCratevsWTot[i] = tfs->make<TH2F>(sview[i]+"ADCratevsWTot",";W ;ADC/Path",
241  wBins,-detHalfWidth,detHalfWidth,nADCbin,0,maxADCrate);
242 
243  fADCratevsWlongTot[i] = tfs->make<TH2F>(sview[i]+"ADCratevsWlongTot",";W ;ADC/Path",
244  wBins,-detHalfWidth,detHalfWidth,nADCbin,0,maxADCrate);
245 
246  fLvsWTot[i] = tfs->make<TH2F>(sview[i]+"LvsWTot",";W ;Path",
247  wBins,-detHalfWidth,detHalfWidth,nADCbin,0,20);
248 
249  fLvsdxTot[i] = tfs->make<TH2F>(sview[i]+"LvsdxTot",";DirX ;Path",
250  wBins,-1,1,nADCbin,0,20);
251 
252  fLvsdyTot[i] = tfs->make<TH2F>(sview[i]+"LvsdyTot",";DirY ;Path",
253  wBins,-1,1,nADCbin,0,20);
254 
255  fLvsdzTot[i] = tfs->make<TH2F>(sview[i]+"LvsdzTot",";DirZ ;Path",
256  wBins,-1,1,nADCbin,0,20);
257 
258  fLvswallsTot[i] = tfs->make<TH2F>(sview[i]+"LvswallsTot",";walls ;Path",
259  5,-0.5,4.5,nADCbin,0,20);
260 
261  fLvsWnothreshTot[i] = tfs->make<TH2F>(sview[i]+"LvsWnothreshTot",";W ;Path",
262  wBins,-detHalfWidth,detHalfWidth,nADCbin,0,20);
263 
264  ftrueThredVsWTot[i] = tfs->make<TH2F>(sview[i]+"trueThredVsWTot", ";W ;Threshold Correction",
265  wBins, -detHalfWidth, detHalfWidth, nADCbin, 0, 5);
266 
267  ftrueShadowVsYTot[i] = tfs->make<TH2F>(sview[i]+"trueShadowVsYTot", ";Y ;Shadow Correction",
268  wBins, -detHalfWidth, detHalfWidth, nADCbin, 0, 5);
269 
270  flambdaVsWTot[i] = tfs->make<TH2F>(sview[i]+"lambdavsWTot",";W;lambda",
271  wBins,-detHalfWidth,detHalfWidth,nADCbin,0,maxADC);
272 
273  fPEVslambdaTot[i] = tfs->make<TH2F>(sview[i]+"PEvslambdaTot",";Lambda;ADC",
274  nADCbin,0,maxADC,nADCbin,0,maxADC);
275 
276  fADCratevsLTot[i] = tfs->make<TH2F>(sview[i]+"ADCratevsLTot",";Path ;ADC/Path",
277  wBins,0,20.0,nADCbin,0,maxADCrate);
278  //new
279  fPEvsLTot[i] = tfs->make<TH2F>(sview[i]+"PEvsLTot",";Path ;PE/Path",
280  wBins,0,20.0,nADCbin,0,maxADC);
281  //new
282  fPECorrvsLTot[i] = tfs->make<TH2F>(sview[i]+"PECorrvsLTot",";Path ;PECorr/Path",
283  wBins,0,20.0,nADCbin,0,maxADC);
284 
285  fADCrateideal[i] = tfs->make<TH1F>(sview[i]+"ADCrateideal",";ADC/Path",
286  nADCbin,0,maxADCrate);
287 
288  htruepath[i] = tfs->make<TH1F>(sview[i]+"truepath",";Path",
289  nADCbin,0,20);
290  //tmp
291  htruepathTest[i] = tfs->make<TH1F>(sview[i]+"truepathTest",";Path",
292  nADCbin,0,20);
293 
294  hrecopath[i] = tfs->make<TH1F>(sview[i]+"recopath",";Path",
295  nADCbin,0,20);
296  //tmp
297  hrecopathTest[i] = tfs->make<TH1F>(sview[i]+"recopathTest",";Path",
298  nADCbin,0,20);
299 
300  fpathtest[i] = tfs->make<TH1F>(sview[i]+"testpath",";Path",
301  nADCbin,0,20);
302 
303  htruew[i] = tfs->make<TH1F>(sview[i]+"truew",";W",
304  wBins,-detHalfWidth,detHalfWidth);
305 
306  hrecow[i] = tfs->make<TH1F>(sview[i]+"recow",";W",
307  wBins,-detHalfWidth,detHalfWidth);
308 
309  hrecodx[i] = tfs->make<TH1F>(sview[i]+"recodx",";DirX",
310  210,-1.05,1.05);
311 
312  hrecody[i] = tfs->make<TH1F>(sview[i]+"recody",";DirY",
313  210,-1.05,1.05);
314  hrecodz[i] = tfs->make<TH1F>(sview[i]+"recodz",";DirZ",
315  210,-1.05,1.05);
316  /*
317  for(int ii = 0; ii<modules;ii++){
318  name = std::to_string(ii);
319 
320  fADCvsYmod[i][ii] = tfs->make<TH2F>(sview[i]+"ADCvsYmod"+name,";Y;ADC",
321  wBins,-detHalfWidth,detHalfWidth,nADCbin,0,maxADC);
322 
323  fADCratevsYmod[i][ii] = tfs->make<TH2F>(sview[i]+"ADCratevsYmod"+name,";Y;ADC/Path",
324  wBins,-detHalfWidth,detHalfWidth,nADCbin,0,maxADCrate);
325 
326  fADCvsWmod[i][ii] = tfs->make<TH2F>(sview[i]+"ADCvsWmod"+name,";W;ADC",
327  wBins,-detHalfWidth,detHalfWidth,nADCbin,0,maxADC);
328 
329  fADCratevsWmod[i][ii] = tfs->make<TH2F>(sview[i]+"ADCratevsWmod"+name,";W;ADC/Path",
330  wBins,-detHalfWidth,detHalfWidth,nADCbin,0,maxADCrate);
331 
332  fLvsWmod[i][ii] = tfs->make<TH2F>(sview[i]+"LvsWmod"+name,";W;Path",
333  wBins,-detHalfWidth,detHalfWidth,nADCbin,0,20);
334  }
335  */
336  }// end loop to make histograms
337 
338  return;
339 }
340 
341 void calib::StopperThreshold::endRun(art::Run const&)
342 {
343 }
344 //---------------------------------------------------------------------
345 void calib::StopperThreshold::reconfigure(fhicl::ParameterSet const & p)
346 {
347  fTrackLabel = p.get< std::string >("TrackLabel" );
348  fMinDistToEdge = p.get< double >("MinDistToEdge" ); //50.0
349  fMinDistFromEnd = p.get< double >("MinDistFromEnd" ); //100.0
350  fMaxDistFromEnd = p.get< double >("MaxDistFromEnd" ); //200.0
351 
352  return;
353 }
354 
355 //---------------------------------------------------------------------
356 void calib::StopperThreshold::analyze(art::Event const & e)
357 {
358 
359  art::Handle< std::vector<rb::Track> > th;
360  e.getByLabel(fTrackLabel, th);
361 
362  for(size_t t = 0; t < th->size(); ++t){
363  rb::Track const& track = th->at(t);
364  if(!track.Is3D()) continue;
365  std::set<TVector3, bool(*)(TVector3 const&, TVector3 const&)> trajPoints(ltTVec);
366  this->FindTrajectoryPoints(track, trajPoints);
367  if(trajPoints.size() == 0) continue;
368  // look at only stopping muon
369  //if(trajPoints.size() == 0 || !isstopper) continue;
370 
371  this->FillTree(track, trajPoints);
372  this->testPath(track, trajPoints);
373  this->FillHist(track, trajPoints);
374  }
375 
376  return;
377 }
378 
379 //-------------------------------------------------------------
380 void calib::StopperThreshold::FindTrajectoryPoints(rb::Track const& track,
381  std::set<TVector3, bool(*)(TVector3 const&, TVector3 const&) > & trajPoints)
382 {
383  // For stopping muons
384  if(fLivegeo->DistToEdgeXY(track.Stop()) > fMinDistToEdge &&
385  fLivegeo->DistToEdgeZ (track.Stop()) > fMinDistToEdge ){ //checks for 50 cm
386  isstopper = true;
387  double totalDist = 0.;
388  for(size_t tp = track.NTrajectoryPoints()-1; tp > 1; --tp){
389  totalDist += (track.TrajectoryPoint(tp) - track.TrajectoryPoint(tp - 1)).Mag();
390  if(totalDist > fMinDistFromEnd && totalDist < fMaxDistFromEnd)
391  trajPoints.insert(track.TrajectoryPoint(tp));
392  if(totalDist > fMaxDistFromEnd) break;
393  }
394  }
395  else{
396  isstopper = false;
397  for(size_t itp = track.NTrajectoryPoints() - 1; itp > 1; --itp){
398  trajPoints.insert(track.TrajectoryPoint(itp));
399  }
400  }// end loop over track trajectorypoints
401 
402  return;
403 }
404 
405 //-------------------------------------------------------------
406 void calib::StopperThreshold::FillHist( rb::Track const& track,
407  std::set<TVector3, bool(*)(TVector3 const&, TVector3 const&)> const& trajPoints)
408 {
409  float minZ = trajPoints.begin()->Z();
410  float maxZ = trajPoints.rbegin()->Z();
411  float cellHalfDepth = fGeo->Plane(0)->Cell(0)->HalfD();
412  if(maxZ < minZ) std::swap(minZ, maxZ);
413 
414  //std::vector<EPathQuality> quals;
415  //const std::vector<double> paths = BestPathEstimates(track.OfflineChans(),track.PlaneDirMap(), quals);
416  //if(track.NCell() != paths.size() || track.NCell() != quals.size()) return;
417  //if(track.NCell() != paths.size()) continue;
418 
419  minZ -= cellHalfDepth;
420  maxZ += cellHalfDepth;
421  float xi = trajPoints.begin()->X();
422  float yi = trajPoints.begin()->Y();
423  float zi = trajPoints.begin()->Z();
424  float xf = trajPoints.rbegin()->X();
425  float yf = trajPoints.rbegin()->Y();
426  float zf = trajPoints.rbegin()->Z();
427  float dr = sqrt((xi-xf)*(xi-xf)+(yi-yf)*(yi-yf)+(zi-zf)*(zi-zf));
428  if(dr == 0) return;
429  float dx = (xf-xi)/dr;
430  float dy = (yf-yi)/dr;
431  float dz = (zf-zi)/dr;
432  double xyz[3]= {0.};
433  unsigned int ntp = track.NTrajectoryPoints();
434  std::pair<uint32_t, uint32_t> pc;
435  TVector3 point;
436  std::set<double> planeZBounds;
437  calib::FindZBoundaries(planeZBounds);
438  auto boundaryMap = calib::MakeZBoundaryMap(planeZBounds, track.Trajectory());
439  if(boundaryMap.size() < 1) return;
440  for(size_t c = 0; c < track.NCell(); ++c){
441  const art::Ptr<rb::CellHit> ch = track.Cell(c);
442 
443  float fpoissonlambda = 0.0;
444  float ftruemev = 0.0;
445  float ftruepath = -10.0;
446  float ftruew = -9999.0;
447  // Check PCHit code later (needs update)
448  if (fBT->HaveTruthInfo()){
449  // count poissonlambda using photonsignal
450  const std::vector<sim::PhotonSignal> pes = fBT->HitToPhotonSignal(track.Cell(c));
451  for(sim::PhotonSignal pe: pes){
452  fpoissonlambda += pe.PoissonLambda();
453  }
454 
455  // count truemev using fls??
456  const std::vector<cheat::TrackIDE> ides = fBT->HitToTrackIDE(track.Cell(c));
457  for(size_t n = 0; n < ides.size(); ++n)
458  ftruemev += ides[n].energy*1000;
459 
460  // calculate true path length and W using fls
461  const std::vector<const sim::Particle*> parts = fBT->HitsToParticle(track.AllCells());
462  if(!parts.empty()){
463  // front() is to select the particle which contributes the most to the reco track.
464  auto particle = parts.front();
465  std::vector<sim::FLSHit> flsHits = fBT->ParticleToFLSHit(particle->TrackId(), particle->PdgCode());
466  for(auto fls : flsHits){
467  if(fls.GetPlaneID() == track.Cell(c)->Plane() &&
468  fls.GetCellID() == track.Cell(c)->Cell() ){
469  ftruepath = fls.GetTotalPathLength();
470  ftruew = fls.GetZAverage();
471  break;
472  }
473  }
474  }
475  } //end HaveTruthInfo
476  if((ftruemev == 0.0) || (fpoissonlambda == 0.0)) continue;
477  //if(ftruepath == 0 || ftruew == 0 || ftruemev == 0.0 || fpoissonlambda == 0.0) continue;
478 
479  // tricell condition
480  /*
481  if(c == 0 || c == track.NCell()-1) continue;
482  const art::Ptr<rb::CellHit> chm1 = track.Cell(c-1);
483  const art::Ptr<rb::CellHit> chp1 = track.Cell(c+1);
484  if(chm1->Plane() != ch->Plane() ||
485  chp1->Plane() != ch->Plane() ||
486  chm1->Cell() != ch->Cell()-1 ||
487  chp1->Cell() != ch->Cell()+1) continue;
488  */
489  fGeo->Plane(ch->Plane())->Cell(ch->Cell())->GetCenter(xyz);
490  double detHalfWidth = std::max(fGeo->DetHalfWidth(), fGeo->DetHalfHeight());
491  // to find the direction of segment_i
492  TVector3 vtp0;
493  TVector3 vtp1;
494  // if cell is near the edge of a traj point
495  if(track.TrajectoryPoint(0).Z() < track.TrajectoryPoint(ntp-1).Z()){
496  if(track.TrajectoryPoint(0).Z() >= xyz[2]){
497  vtp0 = track.TrajectoryPoint(0);
498  vtp1 = track.TrajectoryPoint(1);
499  }
500  else if(track.TrajectoryPoint(ntp-1).Z() <= xyz[2]){
501  vtp0 = track.TrajectoryPoint(ntp-2);
502  vtp1 = track.TrajectoryPoint(ntp-1);
503  }
504  else{
505  for(unsigned int tp=0;tp<ntp;tp++){
506  if(track.TrajectoryPoint(tp).Z() > xyz[2]){
507  vtp0 = track.TrajectoryPoint(tp-1);
508  vtp1 = track.TrajectoryPoint(tp);
509  break;
510  }
511  }
512  }
513  }
514  else{// inverse of the previous case
515  if(track.TrajectoryPoint(0).Z() <= xyz[2]){
516  vtp0 = track.TrajectoryPoint(0);
517  vtp1 = track.TrajectoryPoint(1);
518  }
519  else if(track.TrajectoryPoint(ntp-1).Z() >= xyz[2]){
520  vtp0 = track.TrajectoryPoint(ntp-2);
521  vtp1 = track.TrajectoryPoint(ntp-1);
522  }
523  else{
524  for(unsigned int tp=0;tp<ntp;tp++){
525  if(track.TrajectoryPoint(tp).Z() < xyz[2]){
526  vtp0 = track.TrajectoryPoint(tp-1);
527  vtp1 = track.TrajectoryPoint(tp);
528  break;
529  }
530  }
531  }
532  }
533  TVector3 vdir;
534  vdir = (vtp1 - vtp0).Unit();
535  //auto const vdir = DirectionAtPlane(track.PlaneDirMap(),ch->Plane());
536  float cellwidth = fGeo->Plane(ch->Plane())->Cell(ch->Cell())->HalfW();
537  pc.first = ch->Plane();
538  pc.second = ch->Cell();
539  const rb::RecoHit rh = track.RecoHit(track.Cell(c));
540  point = TVector3(rh.X(), rh.Y(), rh.Z());
541 
542  // whether the track pass tricells.
543  /*
544  if(ch->View() == geo::kX && vdir.Z() != 0.0){
545  double xlow = xyz[0] - cellwidth - vtp0.X();
546  double xhigh = xyz[0] + cellwidth - vtp0.X();
547  double zlow = (xyz[2] - cellHalfDepth - vtp0.Z())*vdir.X()/vdir.Z();
548  double zhigh = (xyz[2] + cellHalfDepth - vtp0.Z())*vdir.X()/vdir.Z();
549  if(zlow > xlow || zlow > xhigh || zhigh < xlow || zhigh < xhigh) continue;
550  }
551  else if(ch->View() == geo::kY && vdir.Z() != 0.0){
552  double ylow = xyz[1] - cellwidth - vtp0.Y();
553  double yhigh = xyz[1] + cellwidth - vtp0.Y();
554  double zlow = (xyz[2] - cellHalfDepth - vtp0.Z())*vdir.Y()/vdir.Z();
555  double zhigh = (xyz[2] + cellHalfDepth - vtp0.Z())*vdir.Y()/vdir.Z();
556  if(zlow > ylow || zlow > yhigh || zhigh < ylow || zhigh < yhigh) continue;
557  }
558  */
559  double fadc = ch->ADC();
560  double fpe = ch->PE();
561 
562  double fpecorr = -5.0;
563  if(rh.IsCalibrated()){
564  fpecorr = rh.PECorr();
565  }
566 
567  double fw = 0.0;
568  int fwall = 0;
569  double fpath = calib::PathLengthInCell(boundaryMap,point,pc);
570  if((ch->View() == geo::kX) && (vdir.X() != 0 || vdir.Z() != 0)){
571 
572  fwall = (int)CalculateWall(rh.X(), rh.Z(), vdir.X(), vdir.Z(),
573  (xyz[0]-cellwidth), (xyz[0]+cellwidth), (xyz[2]-cellHalfDepth), (xyz[2]+cellHalfDepth));
574 
575  if(vdir.Z() != 0.0){
576  fw = (xyz[2]-vtp0.Z())*vdir.Y()/vdir.Z() + vtp0.Y();
577  }
578  else{
579  fw = (xyz[2]-zi)*dy/dz + yi;
580  }
581  }//end if ch->View
582 
583  if((ch->View() == geo::kY) && (vdir.Y() != 0 || vdir.Z() != 0)){
584 
585  fwall = (int)CalculateWall(rh.Y(), rh.Z(), vdir.Y(), vdir.Z(),
586  (xyz[1]-cellwidth), (xyz[1]+cellwidth), (xyz[2]-cellHalfDepth), (xyz[2]+cellHalfDepth));
587 
588  if(vdir.Z() != 0.0){
589  fw = (xyz[2]-vtp0.Z())*vdir.X()/vdir.Z() + vtp0.X();
590  }
591  else{
592  fw = (xyz[2]-zi)*dx/dz + xi;
593  }
594  }//end if ch->View
595 
596  if((fpath < 1.0) || (fpath > 25.0) || (abs(fw) > detHalfWidth) ) continue;
597 
598  double ftempy = 0.0;
599  //int imodulew = 99;
600  //int imoduley = 99;
601  if(ch->View() == geo::kX){
602  ftempy = fw;
603  //imodulew = (int)ch->Cell()/32.0;
604  //imoduley = imodulew;
605  //imoduley = (int)(fw + detHalfWidth)*(fXADC_W.size())/(2.0*detHalfWidth);
606  }
607  else if(ch->View() == geo::kY){
608  ftempy = xyz[1];
609  //imodulew = (int)ch->Cell()/32.0;
610  //imoduley = (int)(fw + detHalfWidth)*12/(2.0*detHalfWidth);
611  }
612 
613  //if(imodulew > 11 || imoduley > 11) continue;
614 
615  float fpathverylong[2] = {6.0,4.5};
616  float fpathlong[2] = {6.0,4.5};
617  int iview = (int)ch->View();
618  fADCvsWTot[iview]->Fill(fw,fadc);
619  fADCratevsWTot[iview]->Fill(fw,fadc/fpath);
620  fADCratevsLTot[iview]->Fill(fpath,fadc/fpath);
621  fPEvsLTot[iview]->Fill(fpath,fpe/fpath);
622  fPECorrvsLTot[iview]->Fill(fpath,fpecorr/fpath);
623  fLvsWTot[iview]->Fill(fw,fpath);
624  fLvsdxTot[iview]->Fill(vdir.X(),fpath);
625  fLvsdyTot[iview]->Fill(vdir.Y(),fpath);
626  fLvsdzTot[iview]->Fill(vdir.Z(),fpath);
627  fLvswallsTot[iview]->Fill(fwall,fpath);
628  hrecodx[iview]->Fill(vdir.X());
629  hrecody[iview]->Fill(vdir.Y());
630  hrecodz[iview]->Fill(vdir.Z());
631  fADCvsYTot[iview]->Fill(ftempy,fadc);
632  fADCratevsYTot[iview]->Fill(ftempy,fadc/fpath);
633  ftrueThredVsWTot[iview]->Fill(fw,fadc/fpoissonlambda);
634  ftrueShadowVsYTot[iview]->Fill(ftempy,ftruemev/(1.78*fpath));
635  fPEVslambdaTot[iview]-> Fill(fpoissonlambda,fadc);
636  flambdaVsWTot[iview]-> Fill(fw,fpoissonlambda);
637 
638  if(ftruepath > 0){
639  htruepath[iview]->Fill(ftruepath);
640  hrecopath[iview]->Fill(fpath);
641  //std::cout<<"ftruepath: "<<ftruepath<<", fpath: "<<fpath<<std::endl;
642 
643  if(((ftruepath > 3.7) && (ftruepath <= 3.9)) //||
644  //((fpath > 3.7) && (fpath <= 3.9))
645  ){
646  htruepathTest[iview]->Fill(ftruepath);
647  hrecopathTest[iview]->Fill(fpath);
648  //std::cout<<"FILLED !!"<<std::endl;
649  //std::cout<<"xi,yi,zi: "<<xi<<", "<<yi<<", "<<zi<<" xf,yf,zf: "<<xf<<", "<<yf<<", "<<zf<<std::endl;
650  }
651 
652  }
653 
654  //angle restricted sample
655  if( (abs(vdir.Y()) < 0.5) && (ftruew > -1000)){
656  htruew[iview]->Fill(ftruew);
657  hrecow[iview]->Fill(fw);
658  }
659 
660  // long sample
661  if( fpath > fpathverylong[iview]){
662  fADCratevsYlongTot[iview]->Fill(ftempy,fadc/fpath); //THIS is used for shadow corr
663  fADCratevsWlongTot[iview]->Fill(fw,fadc/fpath); //this is for threshold corr
664  }
665  // long and near sample
666  if( (fw > 550.0) && (fpath > fpathlong[iview])){
667  fADCrateideal[iview]-> Fill(fadc/fpath); // DEBUGIN
668  }
669 
670  /*
671  fADCvsWmod[iview][imodulew] ->Fill(fw,fadc);
672  fADCratevsWmod[iview][imodulew] ->Fill(fw,fadc/fpath);
673  fLvsWmod[iview][imodulew] ->Fill(fw,fpath);
674  fADCvsYmod[iview][imoduley] ->Fill(ftempy,fadc);
675  fADCratevsYmod[iview][imoduley] ->Fill(ftempy,fadc/fpath);
676  */
677 
678  // MIP region histograms only
679  if(isstopper && xyz[2] > minZ && xyz[2] < maxZ){
680  fADCvsWTot[iview+2]->Fill(fw,fadc);
681  fADCratevsWTot[iview+2]->Fill(fw,fadc/fpath);
682  fADCratevsLTot[iview+2]->Fill(fpath,fadc/fpath);
683  fLvsWTot[iview+2]->Fill(fw,fpath);
684  fLvsdxTot[iview+2]->Fill(vdir.X(),fpath);
685  fLvsdyTot[iview+2]->Fill(vdir.Y(),fpath);
686  fLvsdzTot[iview+2]->Fill(vdir.Z(),fpath);
687  hrecodx[iview+2]->Fill(vdir.X());
688  hrecody[iview+2]->Fill(vdir.Y());
689  hrecodz[iview+2]->Fill(vdir.Z());
690  fLvswallsTot[iview+2]->Fill(fwall,fpath);
691  fADCvsYTot[iview+2]->Fill(ftempy,fadc);
692  fADCratevsYTot[iview+2]->Fill(ftempy,fadc/fpath);
693  ftrueThredVsWTot[iview+2]->Fill(fw,fadc/fpoissonlambda);
694  ftrueShadowVsYTot[iview+2]-> Fill(ftempy,ftruemev/(1.78*fpath));
695  fPEVslambdaTot[iview+2]-> Fill(fpoissonlambda,fadc);
696  flambdaVsWTot[iview+2]-> Fill(fw,fpoissonlambda);
697 
698  if(ftruepath > 0){
699  htruepath[iview+2]->Fill(ftruepath);
700  hrecopath[iview+2]->Fill(fpath);
701  //std::cout<<"ftruepath: "<<ftruepath<<", fpath: "<<fpath<<std::endl;
702 
703  if(((ftruepath > 3.7) && (ftruepath <= 3.9)) //||
704  //((fpath > 3.7) && (fpath <= 3.9))
705  ){
706  htruepathTest[iview]->Fill(ftruepath);
707  hrecopathTest[iview]->Fill(fpath);
708  //std::cout<<"FILLED stoppers!!"<<std::endl;
709  //std::cout<<"xi,yi,zi: "<<xi<<", "<<yi<<", "<<zi<<" xf,yf,zf: "<<xf<<", "<<yf<<", "<<zf<<std::endl;
710  }
711 
712  }
713 
714  if(abs(vdir.Y()) < 0.5 && (ftruew > -1000.0)){
715  htruew[iview+2]->Fill(ftruew);
716  hrecow[iview+2]->Fill(fw);
717  }
718 
719  if(fpath > fpathverylong[iview]){
720  fADCratevsYlongTot[iview+2]->Fill(ftempy,fadc/fpath);
721  fADCratevsWlongTot[iview+2]->Fill(fw,fadc/fpath);
722  }
723 
724  if((fw > 550.0) && (fpath > fpathlong[iview])){
725  fADCrateideal[iview+2] -> Fill(fadc/fpath);
726  }
727  /*
728  fADCvsWmod[iview+2][imodulew] ->Fill(fw,fadc);
729  fADCratevsWmod[iview+2][imodulew]->Fill(fw,fadc/fpath);
730  fLvsWmod[iview+2][imodulew] ->Fill(fw,fpath);
731  fADCvsYmod[iview+2][imoduley] ->Fill(ftempy,fadc);
732  fADCratevsYmod[iview+2][imoduley]->Fill(ftempy,fadc/fpath);
733  */
734  }//end isstopper
735 
736  } // loop over on-track cells
737 
738  return;
739 }
740 
741 //-------------------------------------------------------------
742 void calib::StopperThreshold::FillTree( rb::Track const& track,
743  std::set<TVector3, bool(*)(TVector3 const&, TVector3 const&)> const& trajPoints)
744 {
745  float minZ = trajPoints.begin()->Z();
746  float maxZ = trajPoints.rbegin()->Z();
747  float cellHalfDepth = fGeo->Plane(0)->Cell(0)->HalfD();
748  if(maxZ < minZ) std::swap(minZ, maxZ);
749 
750  minZ -= cellHalfDepth;
751  maxZ += cellHalfDepth;
752  float xi = trajPoints.begin()->X();
753  float yi = trajPoints.begin()->Y();
754  float zi = trajPoints.begin()->Z();
755  float xf = trajPoints.rbegin()->X();
756  float yf = trajPoints.rbegin()->Y();
757  float zf = trajPoints.rbegin()->Z();
758  float dr = sqrt((xi-xf)*(xi-xf)+(yi-yf)*(yi-yf)+(zi-zf)*(zi-zf));
759  if(dr == 0) return;
760  float dx = (xf-xi)/dr;
761  float dy = (yf-yi)/dr;
762  float dz = (zf-zi)/dr;
763  double xyz[3]= {0.};
764  unsigned int ntp = track.NTrajectoryPoints();
765  std::pair<uint32_t, uint32_t> pc;
766  TVector3 point;
767  std::set<double> planeZBounds;
768  calib::FindZBoundaries(planeZBounds);
769  auto boundaryMap = calib::MakeZBoundaryMap(planeZBounds, track.Trajectory());
770  if(boundaryMap.size() < 1) return;
771 
772  // Loop over all selected cells in track
773  for(size_t c = 0; c < track.NCell(); ++c){
774 
775  //Define pointer to a given cell
776  const art::Ptr<rb::CellHit> ch = track.Cell(c);
777 
778  float fpoissonlambda = 0.0;
779  float ftruemev = 0.0;
780  float ftruepath = -10.0;
781  float ftruew = -9999.0;
782  // Check PCHit code later (needs update)
783  if (fBT->HaveTruthInfo()){
784  // count poissonlambda using photonsignal
785  const std::vector<sim::PhotonSignal> pes = fBT->HitToPhotonSignal(track.Cell(c));
786  for(sim::PhotonSignal pe: pes){
787  fpoissonlambda += pe.PoissonLambda();
788  }
789 
790  // count truemev using fls??
791  const std::vector<cheat::TrackIDE> ides = fBT->HitToTrackIDE(track.Cell(c));
792  for(size_t n = 0; n < ides.size(); ++n)
793  ftruemev += ides[n].energy*1000;
794 
795  // calculate true path length and W using fls
796  const std::vector<const sim::Particle*> parts = fBT->HitsToParticle(track.AllCells());
797  if(!parts.empty()){
798  // front() is to select the particle which contributes the most to the reco track.
799  auto particle = parts.front();
800  std::vector<sim::FLSHit> flsHits = fBT->ParticleToFLSHit(particle->TrackId(), particle->PdgCode());
801  for(auto fls : flsHits){
802  if(fls.GetPlaneID() == track.Cell(c)->Plane() &&
803  fls.GetCellID() == track.Cell(c)->Cell() ){
804  ftruepath = fls.GetTotalPathLength();
805  ftruew = fls.GetZAverage();
806  break;
807  }
808  }
809  }
810  }
811  if((ftruemev == 0.0) || (fpoissonlambda == 0.0)) continue;
812 
813  fGeo->Plane(ch->Plane())->Cell(ch->Cell())->GetCenter(xyz);
814  double detHalfWidth = std::max(fGeo->DetHalfWidth(), fGeo->DetHalfHeight());
815  // to find the direction of segment_i
816  TVector3 vtp0;
817  TVector3 vtp1;
818  // if cell is near the edge of a traj point
819  if(track.TrajectoryPoint(0).Z() < track.TrajectoryPoint(ntp-1).Z()){
820  if(track.TrajectoryPoint(0).Z() >= xyz[2]){
821  vtp0 = track.TrajectoryPoint(0);
822  vtp1 = track.TrajectoryPoint(1);
823  }
824  else if(track.TrajectoryPoint(ntp-1).Z() <= xyz[2]){
825  vtp0 = track.TrajectoryPoint(ntp-2);
826  vtp1 = track.TrajectoryPoint(ntp-1);
827  }
828  else{
829  for(unsigned int tp=0;tp<ntp;tp++){
830  if(track.TrajectoryPoint(tp).Z() > xyz[2]){
831  vtp0 = track.TrajectoryPoint(tp-1);
832  vtp1 = track.TrajectoryPoint(tp);
833  break;
834  }
835  }
836  }
837  }
838  else{// inverse of the previous case
839  if(track.TrajectoryPoint(0).Z() <= xyz[2]){
840  vtp0 = track.TrajectoryPoint(0);
841  vtp1 = track.TrajectoryPoint(1);
842  }
843  else if(track.TrajectoryPoint(ntp-1).Z() >= xyz[2]){
844  vtp0 = track.TrajectoryPoint(ntp-2);
845  vtp1 = track.TrajectoryPoint(ntp-1);
846  }
847  else{
848  for(unsigned int tp=0;tp<ntp;tp++){
849  if(track.TrajectoryPoint(tp).Z() < xyz[2]){
850  vtp0 = track.TrajectoryPoint(tp-1);
851  vtp1 = track.TrajectoryPoint(tp);
852  break;
853  }
854  }
855  }
856  }
857  TVector3 vdir;
858  vdir = (vtp1 - vtp0).Unit();
859 
860  float cellwidth = fGeo->Plane(ch->Plane())->Cell(ch->Cell())->HalfW();
861  pc.first = ch->Plane();
862  pc.second = ch->Cell();
863  const rb::RecoHit rh = track.RecoHit(track.Cell(c));
864  point = TVector3(rh.X(), rh.Y(), rh.Z());
865 
866  double fadc = ch->ADC();
867  double fpe = ch->PE();
868 
869  double fpecorr = -5.0;
870  if(rh.IsCalibrated()){
871  fpecorr = rh.PECorr();
872  }
873 
874  double fw = 0.0;
875  int fwall = 0;
876  double fpath = calib::PathLengthInCell(boundaryMap,point,pc);
877  if((ch->View() == geo::kX) && (vdir.X() != 0 || vdir.Z() != 0)){
878 
879  fwall = (int)CalculateWall(rh.X(), rh.Z(), vdir.X(), vdir.Z(),
880  (xyz[0]-cellwidth), (xyz[0]+cellwidth), (xyz[2]-cellHalfDepth), (xyz[2]+cellHalfDepth));
881 
882  if(vdir.Z() != 0.0){
883  fw = (xyz[2]-vtp0.Z())*vdir.Y()/vdir.Z() + vtp0.Y();
884  }
885  else{
886  fw = (xyz[2]-zi)*dy/dz + yi;
887  }
888  }//end if ch->View
889 
890  if((ch->View() == geo::kY) && (vdir.Y() != 0 || vdir.Z() != 0)){
891 
892  fwall = (int)CalculateWall(rh.Y(), rh.Z(), vdir.Y(), vdir.Z(),
893  (xyz[1]-cellwidth), (xyz[1]+cellwidth), (xyz[2]-cellHalfDepth), (xyz[2]+cellHalfDepth));
894 
895  if(vdir.Z() != 0.0){
896  fw = (xyz[2]-vtp0.Z())*vdir.X()/vdir.Z() + vtp0.X();
897  }
898  else{
899  fw = (xyz[2]-zi)*dx/dz + xi;
900  }
901  }//end if ch->View
902 
903  if((fpath < 1.0) || (fpath > 25.0) || (abs(fw) > detHalfWidth) ) continue;
904 
905  double ftempy = 0.0;
906 
907  if(ch->View() == geo::kX){
908  ftempy = fw;
909  }
910  else if(ch->View() == geo::kY){
911  ftempy = xyz[1];
912  }
913 
914  float fpathverylong[2] = {6.0,4.5};
915  float fpathlong[2] = {6.0,4.5};
916  int iview = (int)ch->View();
917  int icell = (int)ch->Cell();
918  int iplane = (int)ch->Plane();
919 
920  tView = iview;
921  tCell = icell;
922  tPlane = iplane;
923  tADC = fadc;
924  tPE = fpe;
925  tPECorr = fpecorr;
926  tW = fw;
927  tTruePath = ftruepath;
928  tPath = fpath;
929  tRecodx = vdir.X();
930  tRecody = vdir.Y();
931  tRecodz = vdir.Z();
932  tWall = fwall;
933  tYTot = ftempy;
934  tPoissonLambda = fpoissonlambda;
935  tTrueMeV = ftruemev;
936  tTrueW = ftruew;
937 
938  bool istruepath = false;
939  if(ftruepath > 0){
940  istruepath = true;
941  }
942  tIsTruePath = istruepath;
943 
944  bool isanglerestr = false;
945  //angle restricted sample
946  if( (abs(vdir.Y()) < 0.5) && (ftruew > -1000)){
947  isanglerestr = true;
948  }
949  tAngleRestr = isanglerestr;
950 
951  // long sample
952  bool islong = false;
953  if( fpath > fpathverylong[iview]){
954  islong = true;
955  }
956  tLongSample = islong;
957 
958  bool islongnear = false;
959  // long and near sample
960  if( (fw > 550.0) && (fpath > fpathlong[iview])){
961  islongnear = true;
962  }
963  tLongNear = islongnear;
964 
965  bool isMIPstopper = false;
966  // MIP region histograms only
967  if(isstopper && xyz[2] > minZ && xyz[2] < maxZ){
968  isMIPstopper = true;
969  }//end isstopper
970  tIsStopper = isMIPstopper;
971 
972  // Filling info tree
973  tInfoTree->Fill();
974  //std::cout<<"Tree Info FILLED!!"<<std::endl;
975 
976  } // loop over on-track cells
977 
978  return;
979 }
980 
981 //-------------------------------------------------------------
982 void calib::StopperThreshold::testPath( rb::Track const& track,
983  std::set<TVector3, bool(*)(TVector3 const&, TVector3 const&)> const& trajPoints)
984 {
985  float minZ = trajPoints.begin()->Z();
986  float maxZ = trajPoints.rbegin()->Z();
987  float cellHalfDepth = fGeo->Plane(0)->Cell(0)->HalfD();
988  if(maxZ < minZ) std::swap(minZ, maxZ);
989  minZ -= cellHalfDepth;
990  maxZ += cellHalfDepth;
991 
992  std::pair<uint32_t, uint32_t> pc;
993  TVector3 point;
994  std::set<double> planeZBounds;
995  calib::FindZBoundaries(planeZBounds);
996  auto boundaryMap = calib::MakeZBoundaryMap(planeZBounds, track.Trajectory());
997  if(boundaryMap.size() < 1) return;
998  for(size_t c = 0; c < track.NCell(); ++c){
999  const art::Ptr<rb::CellHit> ch = track.Cell(c);
1000  pc.first = ch->Plane();
1001  pc.second = ch->Cell();
1002  const rb::RecoHit rh = track.RecoHit(track.Cell(c));
1003  point = TVector3(rh.X(), rh.Y(), rh.Z());
1004  double fpath = calib::PathLengthInCell(boundaryMap,point,pc);
1005  int iview = (int)ch->View();
1006  fpathtest[iview]->Fill(fpath);
1007  double xyz[3]= {0.};
1008  fGeo->Plane(ch->Plane())->Cell(ch->Cell())->GetCenter(xyz);
1009  if(isstopper && xyz[2] > minZ && xyz[2] < maxZ)
1010  fpathtest[iview+2]->Fill(fpath);
1011 
1012  } // loop over on-track cells
1013 
1014  return;
1015 }
1016 
1017 //-------------------------------------------------------------
1018 int calib::StopperThreshold::CalculateWall(float x0, float z0, float dx, float dz,
1019  float xlow, float xhigh, float zlow, float zhigh)
1020 {
1021  int fwalls = -5;
1022  float lambda[4] = {0.0};
1023  float otherview[4] = {0.0};
1024  float sign[4] = {0.0};
1025  if(dx == 0){
1026  if((x0 - xlow)*(x0 - xhigh) < 0){
1027  fwalls = 1;
1028  }
1029  else if((x0 - xlow)*(x0 - xhigh) == 0){
1030  fwalls = 4;
1031  }
1032  else{
1033  fwalls = 0;
1034  }
1035  }//end if dx==0
1036  else if(dz == 0){
1037  if((z0 - zlow)*(z0 - zhigh) < 0){
1038  fwalls = 2;
1039  }
1040  else if((z0 - zlow)*(z0 - zhigh) == 0){
1041  fwalls = 4;
1042  }
1043  else {
1044  fwalls = 0;
1045  }
1046 
1047  }//end else if
1048  else{
1049  lambda[0] = (xlow - x0)/dx;
1050  lambda[1] = (xhigh - x0)/dx;
1051  lambda[2] = (zlow - z0)/dz;
1052  lambda[3] = (zhigh - z0)/dz;
1053  otherview[0] = lambda[0]*dz + z0;
1054  otherview[1] = lambda[1]*dz + z0;
1055  otherview[2] = lambda[2]*dx + x0;
1056  otherview[3] = lambda[3]*dx + x0;
1057  sign[0] = (otherview[0] - zlow)*(otherview[0] - zhigh);
1058  sign[1] = (otherview[1] - zlow)*(otherview[1] - zhigh);
1059  sign[2] = (otherview[2] - xlow)*(otherview[2] - xhigh);
1060  sign[3] = (otherview[3] - xlow)*(otherview[3] - xhigh);
1061 
1062  if((sign[0] > 0) && (sign[1] > 0) &&
1063  (sign[2] > 0) && (sign[3] > 0)){
1064  fwalls = 0;
1065  }
1066  else if((sign[0] == 0) || (sign[1] == 0) || (sign[2] == 0) || (sign[3] == 0)){
1067  fwalls = 4;
1068  }
1069  else if((sign[0] < 0) && (sign[1] < 0)){
1070  fwalls = 1;
1071  }
1072  else if((sign[2] < 0) && (sign[3] < 0)){
1073  fwalls = 2;
1074  }
1075  else {
1076  fwalls = 3;
1077  }
1078 
1079  }//end else
1080 
1081  return fwalls;
1082 }
1083 
1084 //-------------------------------------------------------------
1085 
1086 DEFINE_ART_MODULE(calib::StopperThreshold)
1087 
1088 #endif /* StopperThreshold_h */
std::max
T max(const caf::Proxy< T > &a, T b)
Definition: BasicTypesProxy.h:416
rb::Track::NTrajectoryPoints
size_t NTrajectoryPoints() const
Definition: Track.h:83
calib::StopperThreshold::FillTree
void FillTree(rb::Track const &track, std::set< TVector3, bool(*)(TVector3 const &, TVector3 const &) > const &trajPoints)
Definition: StopperThreshold_module.cc:742
name
const XML_Char * name
Definition: expat.h:151
calib::StopperThreshold::FillHist
void FillHist(rb::Track const &track, std::set< TVector3, bool(*)(TVector3 const &, TVector3 const &) > const &trajPoints)
Definition: StopperThreshold_module.cc:406
BackTracker.h
back track the reconstruction to the simulation
calib::StopperThreshold::analyze
virtual void analyze(art::Event const &e)
Definition: StopperThreshold_module.cc:356
geo::CellGeo::HalfD
double HalfD() const
Definition: CellGeo.cxx:205
rb::Cluster::NCell
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
calib::StopperThreshold::StopperThreshold
StopperThreshold(fhicl::ParameterSet const &p)
Definition: StopperThreshold_module.cc:170
rb::Track::TrajectoryPoint
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
Definition: Track.cxx:158
rh
TH2 * rh
Definition: drawXsec.C:5
calib::MakeZBoundaryMap
zBoundMap MakeZBoundaryMap(std::set< double > const &planeZBounds, std::vector< TVector3 > const &trajectory)
Return a map of the z position of each trajectory point on a track to the bounding positions of where...
Definition: CalibUtil.cxx:354
cheat::BackTracker::HitToTrackIDE
std::vector< TrackIDE > HitToTrackIDE(const rb::CellHit &hit, bool useBirksE=false) const
Convenience function. HitsToTrackIDE but for a single hit.
Definition: BackTracker_service.cc:619
rb::CellHit::Plane
unsigned short Plane() const
Definition: CellHit.h:39
geo::LiveGeometry::DistToEdgeXY
double DistToEdgeXY(TVector3 vertex)
Definition: LiveGeometry_service.cc:568
calib::StopperThreshold::fTrackLabel
std::string fTrackLabel
label of module creating the tracks used
Definition: StopperThreshold_module.cc:81
calib::StopperThreshold::testPath
void testPath(rb::Track const &track, std::set< TVector3, bool(*)(TVector3 const &, TVector3 const &) > const &trajPoints)
Definition: StopperThreshold_module.cc:982
rb::CellHit::View
geo::View_t View() const
Definition: CellHit.h:41
geo::PlaneGeo::Cell
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
p
const char * p
Definition: xmltok.h:285
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
cheat::BackTracker::ParticleToFLSHit
std::vector< sim::FLSHit > ParticleToFLSHit(const int &trackID) const
All the FLSHits that were created by the track id trackID, sorted from most to least light...
Definition: BackTracker_service.cc:1514
geo::kX
Vertical planes which measure X.
Definition: PlaneGeo.h:28
rb::RecoHit::Z
float Z() const
Definition: RecoHit.h:38
track
Definition: event.h:19
abs
void abs(TH1 *hist)
Definition: wrong_sign_uncertainty.C:716
rb::Track
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
geo::GeometryBase::Plane
const PlaneGeo * Plane(unsigned int i) const
Definition: GeometryBase.cxx:749
DEFINE_ART_MODULE
DEFINE_ART_MODULE(TestTMapFile)
geo::kY
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
rb::Cluster::AllCells
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
art::Run
Definition: Run.h:31
rb::RecoHit
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
std::swap
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
Definition: HLTGlobalStatus.h:199
calib::StopperThreshold::fLivegeo
art::ServiceHandle< geo::LiveGeometry > fLivegeo
Definition: StopperThreshold_module.cc:153
rb::CellHit::Cell
unsigned short Cell() const
Definition: CellHit.h:40
calib::StopperThreshold::endRun
virtual void endRun(art::Run const &)
Definition: StopperThreshold_module.cc:341
calib
CDPStorage service.
Definition: AbsCalib_module.cc:25
dy
double dy[NP][NC]
Definition: checkMeasured_Target.C:19
dx
double dx[NP][NC]
Definition: checkMeasured_Target.C:18
geo::LiveGeometry::DistToEdgeZ
double DistToEdgeZ(TVector3 vertex)
Definition: LiveGeometry_service.cc:480
dz
double dz[NP][NC]
Definition: checkMeasured_Target.C:20
Fill
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
a
const double a
Definition: gUpMuFluxGen.cxx:164
fhicl::ParameterSet::get
T get(std::string const &key) const
Definition: ParameterSet.h:231
plot_validation_datamc.c
c
Definition: plot_validation_datamc.py:198
energy
double energy
Definition: plottest35.C:25
rb::CellHit::PE
float PE() const
Definition: CellHit.h:42
rb::Cluster::RecoHit
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Definition: Cluster.cxx:259
rb::RecoHit::IsCalibrated
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
geo::GeometryBase::DetHalfHeight
double DetHalfHeight() const
Definition: GeometryBase.cxx:474
art::EDAnalyzer::EDAnalyzer
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
calib::StopperThreshold::fMaxDistFromEnd
double fMaxDistFromEnd
maximum distance from end of track to use a point
Definition: StopperThreshold_module.cc:84
submit_concat_project.parts
parts
Definition: submit_concat_project.py:367
Unit
TVector3 Unit() const
Definition: SRVector3DAddons.h:14
Event.h
rb::Cluster::Cell
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
calib::StopperThreshold::FindTrajectoryPoints
void FindTrajectoryPoints(rb::Track const &track, std::set< TVector3, bool(*)(TVector3 const &, TVector3 const &)> &trajPoints)
Definition: StopperThreshold_module.cc:380
cheat::BackTracker::HitToPhotonSignal
const std::vector< sim::PhotonSignal > HitToPhotonSignal(const art::Ptr< rb::CellHit > &hit) const
Returns the PhotonSignals contributing the signal in the specified hit.
Definition: BackTracker_service.cc:294
art::TFileDirectory::make
T * make(ARGS...args) const
Definition: TFileDirectory.h:82
geo::GeometryBase::DetHalfWidth
double DetHalfWidth() const
Definition: GeometryBase.cxx:467
rawdata::RawDigit::ADC
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
calib::StopperThreshold::beginRun
virtual void beginRun(art::Run const &)
Definition: StopperThreshold_module.cc:182
art::DataViewImpl::getByLabel
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
ltTVec
bool ltTVec(TVector3 const &a, TVector3 const &b)
Definition: StopperThreshold_module.cc:164
rb::RecoHit::X
float X() const
Definition: RecoHit.h:36
cheat::BackTracker::HaveTruthInfo
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
Definition: BackTracker.h:133
rb::Track::Trajectory
std::vector< TVector3 > const & Trajectory() const
return a constant reference to the track trajectory points
Definition: Track.h:86
Mag
float Mag() const
Definition: SRLorentzVectorAddons.h:10
b
const hit & b
Definition: hits.cxx:21
cellSortZ
bool cellSortZ(art::Ptr< rb::CellHit > const &a, art::Ptr< rb::CellHit > const &b)
Definition: StopperThreshold_module.cc:159
calib::StopperThreshold::fBT
art::ServiceHandle< cheat::BackTracker > fBT
Definition: StopperThreshold_module.cc:154
rb::RecoHit::Y
float Y() const
Definition: RecoHit.h:37
rb::RecoHit::PECorr
float PECorr() const
Definition: RecoHit.cxx:47
rb::Track::Stop
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
calib::FindZBoundaries
void FindZBoundaries(std::set< double > &planeZBounds)
Find the boundaries in the z direction of planes in the detector.
Definition: CalibUtil.cxx:332
e
Float_t e
Definition: plot.C:35
calib::StopperThreshold::fGeo
art::ServiceHandle< geo::Geometry > fGeo
handle to the geometry
Definition: StopperThreshold_module.cc:79
calib::StopperThreshold::fMinDistFromEnd
double fMinDistFromEnd
minimum distance from end of track to use a point
Definition: StopperThreshold_module.cc:83
calib::StopperThreshold::reconfigure
virtual void reconfigure(fhicl::ParameterSet const &p)
Definition: StopperThreshold_module.cc:345
Geometry.h
Encapsulate the geometry of one entire detector (near, far, ndos)
rb::Prong::Is3D
virtual bool Is3D() const
Definition: Prong.h:71
canMan.sign
def sign(x)
Definition: canMan.py:197
cheat::BackTracker::HitsToParticle
std::vector< const sim::Particle * > HitsToParticle(const std::vector< const rb::CellHit * > &hits) const
Returns vector of sim::Particle objects contributing to the given collection of hits.
Definition: BackTracker_service.cc:542
calib::StopperThreshold::CalculateWall
int CalculateWall(float x0, float z0, float dx, float dz, float xlow, float xhigh, float zlow, float zhigh)
Definition: StopperThreshold_module.cc:1018
calib::PathLengthInCell
float PathLengthInCell(zBoundMap const &zBounds, TVector3 const &recoHitLoc, std::pair< uint32_t, uint32_t > const &pc)
Return the path length of a track in the cell in question.
Definition: CalibUtil.cxx:498
string
enum BeamMode string
Generated on Wed Feb 24 2021 14:12:19 for NOvA by   doxygen 1.8.11