TrackEva_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file TrackEva.cxx
3 //
4 // \brief Generating a root file containing information you need for getting
5 // track reconstruction efficiency.
6 //
7 // \version $Id: TrackEva.cxx,v 1.3 2012-09-25 06:30:12 gsdavies Exp $
8 // \author ricken@fnal.gov
9 ////////////////////////////////////////////////////////////////////////
10 #include "RecoBase/Track.h"
11 #include "Simulation/Particle.h"
14 #include "MCCheater/BackTracker.h"
15 
16 
17 // Framework includes
24 
25 #include "TVector3.h"
26 #include "TLorentzVector.h"
27 #include "TNtuple.h"
28 #include "TH1I.h"
29 
30 #include <string>
31 #include <iostream>
32 #include <cmath>
33 #include <algorithm>
34 #include <vector>
35 #include <fstream>
36 
37 
38 class TH1I;
39 class TNtuple;
40 
41 namespace simb { class MCTrajectory; }
42 
43 namespace vf {
44  class TrackEva : public art::EDAnalyzer
45  {
46  public:
47  explicit TrackEva(fhicl::ParameterSet const &pset);
48  virtual ~TrackEva();
49 
50  void analyze(const art::Event &evt);
51  void beginJob();
52  void reconfigure(fhicl::ParameterSet const &pset);
53 
54  static const int GAMMAID;
55  static const int PIZEROID;
56  static const int PIPLUSID;
57  static const int PIMINUSID;
58  static const int NUEID;
59  static const int NUMUID;
60  static const int NUTAUID;
61  static const int MUONID;
62  static const int PROTONID;
63 
64 
65  private:
66 
67  double Tracklength(const simb::MCTrajectory & trajectory);
68 
69  TNtuple* PROTON_RC;
70  TNtuple* MUON_RC;
71  TNtuple* PROTON_MC;
72  TNtuple* MUON_MC;
73 
74  TH1I* track_Nb;
75  TH1I* proton_Nb;
76  TH1I* proton_Tr;
77  TH1I* muon_Nb;
78  TH1I* muon_Tr;
79 
81  double fMatching;
82  double fEdgeX;
83  double fEdgeY;
84  double fEdgeZ;
85  double fVTX;
86  double fVTY;
87  double fVTZ;
88  };
89 
90 
91  const int TrackEva::GAMMAID = 22;
92  const int TrackEva::PIZEROID = 111;
93  const int TrackEva::PIPLUSID = 211;
94  const int TrackEva::PIMINUSID = -211;
95  const int TrackEva::NUEID = 12;
96  const int TrackEva::MUONID = 13;
97  const int TrackEva::NUMUID = 14;
98  const int TrackEva::NUTAUID = 16;
99  const int TrackEva::PROTONID = 2212;
100 
101 
102  TrackEva::TrackEva(fhicl::ParameterSet const& pset)
103  : EDAnalyzer(pset)
104  {
105  reconfigure(pset);
106  }
107 
108 
109 //......................................................................
111  {
112  Track_RC =pset.get< std::string >("Track_RC" );
113 
114  fMatching =pset.get< double >("MatchingPurity" );
115 
116  fEdgeX =pset.get< double >("EdgeX");
117  fEdgeY =pset.get< double >("EdgeY");
118  fEdgeZ =pset.get< double >("EdgeZ");
119 
120  }
121 
122 //......................................................................
123 
125 
126 
127 //......................................................................
128 
130  {
132 
133 
134  //The following two ntuples can give info:
135  // efficiency VS track length
136 
137  MUON_MC = tfs->make<TNtuple>("MUON_MC","muon_mc_info",
138  "energy:trackLength:Dir_Cos:ID");
139 
140  PROTON_MC = tfs->make<TNtuple>("PROTON_MC","proton_mc_info",
141  "energy:trackLength:Dir_Cos:ID");
142 
143  MUON_RC = tfs->make<TNtuple>("MUON_RC","muon_rc_info",
144  "tracklength_diff:energy_diff:Angle:TL_MC:E_MC:Dir_Cos:ID");
145 
146  PROTON_RC = tfs->make<TNtuple>("PROTON_RC","proton_rc_info",
147  "tracklength_diff:energy_diff:Angle:TL_MC:E_MC:Dir_Cos:ID");
148 
149 
150  track_Nb = tfs->make<TH1I>("track_Nb","Total Reconstructed Track-Number",
151  1001,-0.5,1000.5);
152 
153  proton_Nb = tfs->make<TH1I>("proton_Nb","Total Proton Number",
154  501,-0.5,500.5);
155  proton_Tr = tfs->make<TH1I>("proton_Tr","Reco Proton Track Number",
156  500,-0.5,500.5);
157 
158  muon_Nb = tfs->make<TH1I>("muon_Nb","Total Muon Number",
159  1001,-0.5,1000.5);
160  muon_Tr = tfs->make<TH1I>("muon_Tr","Reco Muon Track Number",
161  1001,-0.5,1000.5);
162 
163  }
164 
166  {
167  std::vector<double> Mu_E, P_E, LMC_mu, LMC_p;
168  std::vector<TVector3> P_Momenta, Mu_Momenta;
169 
170  TVector3 momentum_P, momentum_Mu;
171 
173  try{ evt.getByLabel(Track_RC,trackhandle); }
174  catch(...){ std::cout<<"couldn't get the trackhandle "<<std::endl; return;}
175 
176  //make a backtracker
178  const sim::ParticleNavigator& pnav = bt->ParticleNavigator();
179 
180 
181  //==============================================================
182 
183  int NbTrack = trackhandle->size();
184  track_Nb->Fill(NbTrack);
185 
186  // Get MC infor of the proton and muon:
187 
188  std::vector<int> MuonID, ProtonID;
189 
190  for (sim::ParticleNavigator::const_iterator i = pnav.begin(); i != pnav.end(); ++i){
191  const sim::Particle *p = (*i).second;
192 
193  if (p->PdgCode()==PROTONID ) {
194  LMC_p.push_back( Tracklength(p->Trajectory()));
195  momentum_P.SetXYZ(p->Px(),
196  p->Py(),
197  p->Pz());
198  P_Momenta.push_back(momentum_P);
199  P_E.push_back(p->E());
200  ProtonID.push_back(p->TrackId());
201  }
202  else if (p->PdgCode()==MUONID ){
203  LMC_mu.push_back( Tracklength(p->Trajectory()));
204  momentum_Mu.SetXYZ(p->Px(),
205  p->Py(),
206  p->Pz());
207  Mu_Momenta.push_back(momentum_Mu);
208  Mu_E.push_back( p->E());
209  MuonID.push_back(p->TrackId());
210  }
211  }
212 
213  //Filling proton mc ntuple
214  for (unsigned int n=0; n!=P_E.size();n++){
215  PROTON_MC->Fill(P_E[n],
216  LMC_p[n],
217  P_Momenta[n].Z()/P_Momenta[n].Mag(),
218  ProtonID[n]);
219  }
220 
221  //Filling muon mc ntuple
222  for (unsigned int n=0;n!=Mu_E.size();++n){
223  MUON_MC->Fill(Mu_E[n],
224  LMC_mu[n],
225  Mu_Momenta[n].Z()/Mu_Momenta[n].Mag(),
226  MuonID[n]);
227  }
228 
229  //Filling MC truth particle number info:
230 
231  proton_Nb->Fill(P_E.size());
232  muon_Nb ->Fill(Mu_E.size());
233 
234 
235  // Let's count qualified Tracks:
236 
237  int mu_Track=0;
238  int p_Track=0;
239 
240  for (int i=0; i!=NbTrack; i++){
241  art::Ptr<rb::Track> track(trackhandle, i);
242  int cell_Nb = track->NCell();
243 
244  //Muon Track?
245  for (unsigned int j=0; j!=MuonID.size(); j++){
246  int cell_count=0;
247 
248  for (int idx = 0; idx != cell_Nb; idx++){
249 
250  bool real_hit= false;
251 
252  std::vector<cheat::TrackIDE> trackides= bt->HitToTrackIDE(track->Cell(idx));
253 
254  for (unsigned int pts=0; pts!=trackides.size(); pts++){
255  if (trackides[pts].trackID==MuonID[j]) real_hit=true;
256  }
257 
258  if (real_hit) cell_count++;
259  }
260  if ((cell_count+0.0)/(cell_Nb+0.0)>=fMatching){
261  mu_Track++;
262  MUON_RC->Fill(
263  track->TotalLength()-LMC_mu[j],
264  track->TotalGeV()-Mu_E[j],
265  Mu_Momenta[j].Angle(track->Dir()),
266  LMC_mu[j],
267  Mu_E[j],
268  Mu_Momenta[j].Z()/Mu_Momenta[j].Mag(),
269  MuonID[j]
270  );
271  }
272  }
273 
274 
275  //Proton Track?
276  for (unsigned int j=0; j<ProtonID.size(); j++){
277  int cell_count=0;
278 
279  for (int idx = 0; idx < cell_Nb; idx++){
280 
281  bool real_hit= false;
282  std::vector<cheat::TrackIDE> trackides= bt->HitToTrackIDE(track->Cell(idx));
283  for (unsigned int pts=0; pts!=trackides.size(); pts++){
284  if (trackides[pts].trackID==ProtonID[j]) real_hit=true;
285  }
286 
287  if (real_hit) cell_count++;
288  }
289  if ((cell_count+0.0)/(cell_Nb+0.0)>fMatching){
290  p_Track++;
291  PROTON_RC->Fill(
292  track->TotalLength()-LMC_p[j],
293  track->TotalGeV()-P_E[j],
294  P_Momenta[j].Angle(track->Dir()),
295  LMC_p[j],
296  P_E[j],
297  P_Momenta[j].Z()/P_Momenta[j].Mag(),
298  ProtonID[j]
299  );
300  }
301  }
302  }
303 
304 
305  proton_Tr -> Fill(p_Track);
306  muon_Tr -> Fill(mu_Track);
307 
308  return;
309  }
310 
311  //Calculate the track length inside the detector
312  double TrackEva::Tracklength(const simb::MCTrajectory& trajectory){
313 
314  double length=0;
315 
316  for(simb::MCTrajectory::const_iterator it = trajectory.begin(); it!=trajectory.end()-1; ++it)
317  {
318  TLorentzVector pos1 = it->first;
319  TLorentzVector pos2 = (it+1)->first;
320 
321  if (fabs(pos2.X())<fEdgeX && fabs(pos2.Y())<fEdgeY){
322  if(pos2.Z()>9.96604 && pos2.Z()<fEdgeZ){
323  length+=(pos1.Vect()-pos2.Vect()).Mag();
324  }
325  }
326 
327  }
328  return (length);
329 
330  }
331 
332 
334 
335 }
double E(const int i=0) const
Definition: MCParticle.h:232
This is the version 2, works for multi-track event.
back track the reconstruction to the simulation
int PdgCode() const
Definition: MCParticle.h:211
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
static const int PIZEROID
double Py(const int i=0) const
Definition: MCParticle.h:230
Trajectory class.
set< int >::iterator it
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
std::vector< TrackIDE > HitToTrackIDE(const rb::CellHit &hit, bool useBirksE=false) const
Convenience function. HitsToTrackIDE but for a single hit.
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
Definition: BackTracker.h:744
const simb::MCTrajectory & Trajectory() const
Definition: MCParticle.h:252
const char * p
Definition: xmltok.h:285
TNtuple * PROTON_MC
list_type::const_iterator const_iterator
double Px(const int i=0) const
Definition: MCParticle.h:229
virtual ~TrackEva()
static const int PROTONID
double Tracklength(const simb::MCTrajectory &trajectory)
void reconfigure(fhicl::ParameterSet const &pset)
DEFINE_ART_MODULE(TestTMapFile)
int TrackId() const
Definition: MCParticle.h:209
static const int NUTAUID
Float_t Z
Definition: plot.C:38
virtual double TotalLength() const
Length (cm) of all the track segments.
Definition: Track.cxx:213
static const int NUMUID
length
Definition: demo0.py:21
void beginJob()
static const int MUONID
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
list_type::const_iterator const_iterator
Definition: MCTrajectory.h:66
T get(std::string const &key) const
Definition: ParameterSet.h:231
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
static const int GAMMAID
int evt
void analyze(const art::Event &evt)
static const int NUEID
const double j
Definition: BetheBloch.cxx:29
This class describes a particle created in the detector Monte Carlo simulation.
OStream cout
Definition: OStream.cxx:6
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
static const int PIPLUSID
std::string Track_RC
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
TNtuple * MUON_RC
float MuonID(const caf::SRProxy *sr)
Definition: MuonID.cxx:51
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
double Pz(const int i=0) const
Definition: MCParticle.h:231
float Mag() const
TNtuple * PROTON_RC
double TotalGeV(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple sum of the estimated GeV of all the hits.
Definition: Cluster.cxx:378
static const int PIMINUSID
Int_t trackID
Definition: plot.C:84
Float_t track
Definition: plot.C:35
TNtuple * MUON_MC