CosmicTrends_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 
19 #include "fhiclcpp/ParameterSet.h"
21 
22 #include "Geometry/Geometry.h"
23 #include "RecoBase/Track.h"
24 #include "MCCheater/BackTracker.h"
26 
27 #include "TStopwatch.h"
28 #include "TH1.h"
29 #include "TH2.h"
30 #include "TMath.h"
31 
32 namespace calib {
33  class CosmicTrends;
34 
35  enum _quality{
36  kAll = 0,
39  };
40 
51  };
52 
66  };
67 
68 }
69 
71 public:
72  explicit CosmicTrends(fhicl::ParameterSet const & p);
73  virtual ~CosmicTrends();
74 
75  virtual void analyze(art::Event const & e);
76 
77  virtual void beginRun(art::Run const&);
78  virtual void reconfigure(fhicl::ParameterSet const & p);
79 
80 private:
81 
82  void FillHistograms(rb::Track const& track,
83  unsigned int const& qualityLevel,
84  std::map<unsigned short, std::set<unsigned short> > const& planesAndCells,
85  bool mc);
86 
87  // useable tracks are defined as those that pass through a cell and its
88  // neighbors to either side in the same view
89  // useable tracks should also span a total number of planes in the detector
90 
91  double fTriCellFrac; ///< minimum fraction of tri-cells in a track
92  unsigned int fMinPlanes; ///< minimum number of planes a track must cross
93  ///< before being deemed useable
94  double fContainmentCut; ///< #cm in from the detector boundary to be called a stopper
95  std::string fTrackLabel; ///< label of module creating the tracks we use
96 
97  std::map<unsigned int, std::vector<TH1D*> > fOneDHists; ///< all 1D hists to make
98  std::map<unsigned int, std::vector<TH2D*> > fTwoDHists; ///< all 2D hists to make
99 
101 
102 };
103 
104 //---------------------------------------------------------------------
106  : EDAnalyzer(p)
107 {
108  fOneDHists.clear();
109  fTwoDHists.clear();
110 
111  this->reconfigure(p);
112 }
113 
114 //---------------------------------------------------------------------
116 {
117 }
118 
119 //---------------------------------------------------------------------
121 {
122 
123  // get the tracks for this event
125  e.getByLabel(fTrackLabel, th);
126 
127  int numTriCellTracks = 0;
128 
129  // loop over all the tracks, use only those that cross the minimum
130  // number of planes
131  for( size_t t = 0; t < th->size(); ++t){
132 
133  rb::Track const& track = th->at(t);
134 
135  // now loop over the hits in the track and make a
136  // map of plane and cells that are hit
137  std::map<unsigned short, std::set<unsigned short> > planesAndCells;
138 
139  // loop over the planes in the map and check out the cells
140  // that were hit
141  double useableCells[2] = {0.};
142  double totalCells[2] = {0.};
143 
144  for(size_t h = 0; h < track.NCell(); ++h)
145  planesAndCells[track.Cell(h)->Plane()].insert(track.Cell(h)->Cell());
146 
147  useableCells[0] = 0.;
148  useableCells[1] = 0.;
149  totalCells[0] = 0.;
150  totalCells[1] = 0.;
151 
152  for(auto const& pacitr : planesAndCells){
153 
154  unsigned short plane = pacitr.first;
155  // get the set for the current plane
156  // sets are presorted, go c++
157  std::set<unsigned short> cellSet = pacitr.second;
158 
159  // put the iterator over the set to the second entry
160  //auto citr = cellSet.cbegin();
161  //totalCells[fGeo->Plane(pacitr.first)->View()] += 1.;
162  //++citr;
163  //unsigned short cell = citr;
164  // loop over the cells and look to see if there is
165  // a cell in this plane that has both its neighbors also
166  // hit on this track
167  //while( citr != cellSet.end() ){
168  for( auto citr : cellSet ){
169  if( cellSet.find(citr - 1) != cellSet.end() &&
170  cellSet.find(citr + 1) != cellSet.end() ){
171  useableCells[fGeo->Plane(plane)->View()] += 1.;
172 
173  if(fGeo->Plane(plane)->View() == geo::kX){
174  fTwoDHists[kAll][kTotalCountX]->Fill(plane, citr);
175  }
176  else if(fGeo->Plane(plane)->View() == geo::kY){
177  fTwoDHists[kAll][kTotalCountY]->Fill(plane, citr);
178  }
179  }
180 
181  if(fGeo->Plane(plane)->View() == geo::kX){
182  fTwoDHists[kAll][kTotalHitCountX]->Fill(plane, citr);
183  }
184  else if(fGeo->Plane(plane)->View() == geo::kY){
185  fTwoDHists[kAll][kTotalHitCountY]->Fill(plane, citr);
186  }
187 
188  totalCells[fGeo->Plane(plane)->View()] += 1.;
189  ++citr;
190  }// end iterator over cells for this plane
191 
192  } // end iterator over plane/cell map
193 
194  double usefulFrac = 0.;
195  if(totalCells[0] + totalCells[1] > 0)
196  usefulFrac = (useableCells[0] + useableCells[1])/(totalCells[0] + totalCells[1]);
197 
198  // loop over the planesAndCells map again and fill histograms
199  // for cells that are passed tracks
200  unsigned int quality = (usefulFrac >= fTriCellFrac) ? kTriCell : kAll;
201 
202  if(totalCells[0] > 0){
203  fOneDHists[kAll][kFractionUsableX]->Fill(useableCells[0]/totalCells[0]);
204  if(quality == kTriCell)
205  fOneDHists[quality][kFractionUsableX]->Fill(useableCells[0]/totalCells[0]);
206  }
207  if(totalCells[1] > 0){
208  fOneDHists[kAll][kFractionUsableY]->Fill(useableCells[1]/totalCells[1]);
209  if(quality == kTriCell)
210  fOneDHists[quality][kFractionUsableY]->Fill(useableCells[1]/totalCells[1]);
211  }
212 
213  this->FillHistograms(track, kAll, planesAndCells, !e.isRealData());
214  if(quality == kTriCell){
215  ++numTriCellTracks;
216  this->FillHistograms(track, quality, planesAndCells, !e.isRealData());
217  }
218 
219  }// end loop over tracks in this event
220 
221  fOneDHists[kTriCell][kNumTracks]->Fill(numTriCellTracks);
222  fOneDHists[kAll][kNumTracks]->Fill(th->size());
223 
224  return;
225 }
226 
227 //---------------------------------------------------------------------
229 {
230  if(fOneDHists.size() > 0) return;
231 
233 
234  // use the geometry to get the number of cells for planes in each view
235  // also get the total number of planes
236 
237  unsigned int planeInView = *(fGeo->GetPlanesByView(geo::kX).begin());
238  double xCells = fGeo->Plane(planeInView)->Ncells();
239 
240  planeInView = *(fGeo->GetPlanesByView(geo::kY).begin());
241  double yCells = fGeo->Plane(planeInView)->Ncells();
242 
243  double totalPlanes = fGeo->NPlanes();
244  double maxCells = (yCells > xCells) ? yCells : xCells;
245 
246  for(int q = kAll; q < kMaxQuality; ++q){
247 
248  std::string qualString("AllTracks");
249  if(q == kTriCell) qualString = "TriCell";
250 
251  // # of useable stopping muons passing through each cell
252  fTwoDHists[q].push_back(tfs->make<TH2D>(("stopperCountX"+qualString).c_str(),
253  ";Plane;Cell",
254  (int)totalPlanes, 0., totalPlanes,
255  (int)maxCells, 0., maxCells));
256 
257  // # of useable muons passing through each cell
258  fTwoDHists[q].push_back(tfs->make<TH2D>(("totalCountX"+qualString).c_str(),
259  ";Plane;Cell",
260  (int)totalPlanes, 0., totalPlanes,
261  (int)maxCells, 0., maxCells));
262 
263  // # of useable stopping muons passing through each cell
264  fTwoDHists[q].push_back(tfs->make<TH2D>(("stopperCountY"+qualString).c_str(),
265  ";Plane;Cell",
266  (int)totalPlanes, 0., totalPlanes,
267  (int)maxCells, 0., maxCells));
268  // # of useable muons passing through each cell
269  fTwoDHists[q].push_back(tfs->make<TH2D>(("totalCountY"+qualString).c_str(),
270  ";Plane;Cell",
271  (int)totalPlanes, 0., totalPlanes,
272  (int)maxCells, 0., maxCells));
273 
274  // reconstructed cos(theta) vs reconstructed azimuth
275  fTwoDHists[q].push_back(tfs->make<TH2D>(("cosThetaVsAz"+qualString).c_str(),
276  ";cos(#theta);#phi (degrees)",
277  100, 0., 1.,
278  36, 0., 360.));
279 
280  // reconstructed cos(theta) vs reconstructed azimuth
281  fTwoDHists[q].push_back(tfs->make<TH2D>(("cosThetaVsAzTru"+qualString).c_str(),
282  ";True cos(#theta);True #phi (degrees)",
283  100, 0., 1.,
284  36, 0., 360.));
285 
286  // reconstructed cos(theta) vs planes crossed
287  fTwoDHists[q].push_back(tfs->make<TH2D>(("cosThetaVsPlanes"+qualString).c_str(),
288  ";cos(#theta);Planes Crossed",
289  100, 0., 1.,
290  (int)totalPlanes, 0., totalPlanes));
291 
292  // reconstructed azimuth vs planes crossed
293  fTwoDHists[q].push_back(tfs->make<TH2D>(("azimuthVsPlanes"+qualString).c_str(),
294  ";#phi (degrees);Planes Crossed",
295  36, 0., 360.,
296  (int)totalPlanes, 0., totalPlanes));
297 
298  // reco - true vs reconstructed azimuth
299  fTwoDHists[q].push_back(tfs->make<TH2D>(("deltaAzVsAz"+qualString).c_str(),
300  ";#phi (degrees);#Delta#phi (degrees)",
301  36, 0., 360.,
302  360, -180., 180.));
303 
304  // reco - true vs reconstructed cos(theta)
305  fTwoDHists[q].push_back(tfs->make<TH2D>(("deltaZenVsZen"+qualString).c_str(),
306  ";cos(#theta);#Deltacos(#theta)",
307  100, 0., 1.,
308  200, -1., 1.));
309 
310  // fraction of tri-cells in a muon track
311  fOneDHists[q].push_back(tfs->make<TH1D>(("fractionUsableX"+qualString).c_str(),
312  ";Usable Fraction;Number of Muons",
313  100, 0., 1.));
314 
315  // fraction of tri-cells in a muon track
316  fOneDHists[q].push_back(tfs->make<TH1D>(("fractionUsableY"+qualString).c_str(),
317  ";Usable Fraction;Number of Muons",
318  100, 0., 1.));
319 
320  // reconstructed cos(theta)
321  fOneDHists[q].push_back(tfs->make<TH1D>(("recoCosTheta"+qualString).c_str(),
322  ";cos(#theta);Number of Muons",
323  100, 0., 1. ));
324 
325  // reconstructed azimuth
326  fOneDHists[q].push_back(tfs->make<TH1D>(("recoAzimuth"+qualString).c_str(),
327  ";#phi (degrees);Number of Muons",
328  36, 0., 360. ));
329 
330  // reconstructed number of planes crossed by track
331  fOneDHists[q].push_back(tfs->make<TH1D>(("planesCrossed"+qualString).c_str(),
332  ";Planes Crossed;Number of Muons",
333  (int)totalPlanes*0.2, 0., (int)totalPlanes));
334 
335  // direction cosine in z
336  fOneDHists[q].push_back(tfs->make<TH1D>(("dcosZ"+qualString).c_str(),
337  ";dz/ds;Number of Tracks",
338  40, -1., 1.));
339 
340  // direction cosine in x
341  fOneDHists[q].push_back(tfs->make<TH1D>(("dcosX"+qualString).c_str(),
342  ";dx/ds;Number of Tracks",
343  40, -1., 1.));
344 
345  // direction cosine in y
346  fOneDHists[q].push_back(tfs->make<TH1D>(("dcosY"+qualString).c_str(),
347  ";dy/ds;Number of TriCell Tracks",
348  40, -1., 1.));
349 
350  // tracks in event
351  fOneDHists[q].push_back(tfs->make<TH1D>(("numTracks"+qualString).c_str(),
352  ";Number of TriCell Tracks;Events",
353  10, 0., 10.));
354 
355  // # of muons passing through each cell
356  fTwoDHists[q].push_back(tfs->make<TH2D>(("totalHitCountX"+qualString).c_str(),
357  ";Plane;Cell",
358  (int)totalPlanes, 0., totalPlanes,
359  (int)maxCells, 0., maxCells));
360  // # of muons passing through each cell
361  fTwoDHists[q].push_back(tfs->make<TH2D>(("totalHitCountY"+qualString).c_str(),
362  ";Plane;Cell",
363  (int)totalPlanes, 0., totalPlanes,
364  (int)maxCells, 0., maxCells));
365 
366 
367 
368  }
369 
370  return;
371 }
372 
373 //---------------------------------------------------------------------
375 {
376  fTriCellFrac = p.get< double >("TriCellFraction", 0.6);
377  fMinPlanes = p.get< unsigned int >("MinimumPlanesPerTrack", 8);
378  fTrackLabel = p.get< std::string >("TrackLabel" );
379  fContainmentCut = p.get< double >("ContainmentCut", 20.);
380 
381  return;
382 }
383 
384 //---------------------------------------------------------------------
386  unsigned int const& qualityLevel,
387  std::map<unsigned short, std::set<unsigned short> > const& planesAndCells,
388  bool mc)
389 {
390 
391  fOneDHists[qualityLevel][kPlanesCrossed]->Fill(track.ExtentPlane());
392 
393  if(track.ExtentPlane() < fMinPlanes) return;
394 
395  TVector3 startDir = track.Dir();
396  double cosTheta = -startDir.y()/startDir.Mag();
397  double azimuth = std::atan2(startDir.x()/startDir.Mag(), -startDir.z()/startDir.Mag()) * 180./TMath::Pi();
398  if(azimuth < 0) azimuth += 360.;
399 
400  fTwoDHists[qualityLevel][kCosThetaVsAz] ->Fill(cosTheta, azimuth);
401  fTwoDHists[qualityLevel][kCosThetaVsPlanes]->Fill(cosTheta, track.ExtentPlane());
402  fTwoDHists[qualityLevel][kAzVsPlanes] ->Fill(azimuth, track.ExtentPlane());
403 
404  fOneDHists[qualityLevel][kRecoCosTheta]->Fill(cosTheta);
405  fOneDHists[qualityLevel][kRecoAz] ->Fill(azimuth);
406 
407  // add code to determine reconstruction qualityLevel of the
408  // tracking in terms of azimuth and cos(theta) when
409  // back tracker is able to do that.
410  if(mc){
412 
413  // get the MCParticle that contributes the most information to the
414  // current track
415  const sim::Particle* truth = bt->HitsToParticle(track.AllCells()).at(0);
416  if(!truth)
417  LOG_WARNING("CosmicTrends") << "no simulated particle matches current track, skip it";
418  else{
419  float trueCosTheta = -truth->Py()/truth->P();
420  float trueAzimuth = std::atan2(truth->Px()/truth->P(), -truth->Pz()/truth->P()) * 180./TMath::Pi();
421  if(trueAzimuth < 0) trueAzimuth += 360.;
422  fTwoDHists[qualityLevel][kCosThetaVsAzTru]->Fill(trueCosTheta, trueAzimuth);
423  fTwoDHists[qualityLevel][kDeltaAzVsAz] ->Fill(azimuth, azimuth-trueAzimuth);
424  fTwoDHists[qualityLevel][kDeltaZenVsZen] ->Fill(cosTheta, cosTheta-trueCosTheta);
425  }
426  }
427 
428  // check if this track stops in the detector
429  TVector3 startPoint(track.Start());
430  TVector3 endPoint(track.Stop());
431  bool stopper = false;
432  if(std::abs(endPoint.X()) < fGeo->DetHalfWidth() - fContainmentCut &&
433  std::abs(endPoint.Y()) < fGeo->DetHalfHeight() - fContainmentCut &&
434  endPoint.Z() < fGeo->DetLength() - fContainmentCut &&
435  endPoint.Z() > fContainmentCut )
436  stopper = true;
437 
438  for(auto const& pacitr : planesAndCells){
439 
440  unsigned short plane = pacitr.first;
441  // get the set for the current plane
442  // sets are presorted, go c++
443  std::set<unsigned short> cellSet = pacitr.second;
444 
445  // loop over the cells and look to see if there is
446  // a cell in this plane that has both its neighbors also
447  // hit on this track
448  for( auto citr : cellSet ){
449 
450 // if(fGeo->Plane(plane)->View() == geo::kX)
451 // fTwoDHists[qualityLevel][kTotalCountX]->Fill(plane, citr);
452 // else if(fGeo->Plane(plane)->View() == geo::kY)
453 // fTwoDHists[qualityLevel][kTotalCountY]->Fill(plane, citr);
454 
455  if(stopper){
456  if(fGeo->Plane(plane)->View() == geo::kX)
457  fTwoDHists[qualityLevel][kStopperCountX]->Fill(plane, citr);
458  else if(fGeo->Plane(plane)->View() == geo::kY)
459  fTwoDHists[qualityLevel][kStopperCountY]->Fill(plane, citr);
460  }
461 
462  ++citr;
463  }// end iterator over cells for this plane
464  } // end iterator over plane/cell map
465 
466  fOneDHists[qualityLevel][kDCosZ]->Fill(startDir.Z()/startDir.Mag());
467  fOneDHists[qualityLevel][kDCosX]->Fill(startDir.X()/startDir.Mag());
468  fOneDHists[qualityLevel][kDCosY]->Fill(startDir.Y()/startDir.Mag());
469 
470  return;
471 }
472 
474 
475 #endif /* CosmicTrends_h */
back track the reconstruction to the simulation
double fContainmentCut
#cm in from the detector boundary to be called a stopper
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
virtual void analyze(art::Event const &e)
double Py(const int i=0) const
Definition: MCParticle.h:230
unsigned short Plane() const
Definition: CellHit.h:39
art::ServiceHandle< geo::Geometry > fGeo
std::map< unsigned int, std::vector< TH2D * > > fTwoDHists
all 2D hists to make
const char * p
Definition: xmltok.h:285
Vertical planes which measure X.
Definition: PlaneGeo.h:28
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
double Px(const int i=0) const
Definition: MCParticle.h:229
Definition: event.h:19
double DetLength() const
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
const PlaneGeo * Plane(unsigned int i) const
bool isRealData() const
Definition: Event.h:83
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
float abs(float number)
Definition: d0nt_math.hpp:39
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
Definition: Run.h:31
std::map< unsigned int, std::vector< TH1D * > > fOneDHists
all 1D hists to make
unsigned short Cell() const
Definition: CellHit.h:40
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
CDPStorage service.
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
double P(const int i=0) const
Definition: MCParticle.h:233
T get(std::string const &key) const
Definition: ParameterSet.h:231
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
double fTriCellFrac
minimum fraction of tri-cells in a track
double DetHalfHeight() const
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
#define LOG_WARNING(category)
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
double DetHalfWidth() const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
std::string fTrackLabel
label of module creating the tracks we use
const std::set< unsigned int > & GetPlanesByView(View_t v=kXorY) const
void FillHistograms(rb::Track const &track, unsigned int const &qualityLevel, std::map< unsigned short, std::set< unsigned short > > const &planesAndCells, bool mc)
double Pz(const int i=0) const
Definition: MCParticle.h:231
unsigned int ExtentPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.h:250
virtual void reconfigure(fhicl::ParameterSet const &p)
unsigned int NPlanes() const
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
Float_t e
Definition: plot.C:35
virtual void beginRun(art::Run const &)
Encapsulate the geometry of one entire detector (near, far, ndos)
CosmicTrends(fhicl::ParameterSet const &p)
T atan2(T number)
Definition: d0nt_math.hpp:72
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.