MonopoleAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file MonopoleAna.cxx
3 // \brief Analyse monopole monte carlo events
4 // \version $Id: MonopoleAna.cxx,v 1.6 2012-08-23 01:52:17 bckhouse Exp $
5 // \author $Author: bckhouse $
6 ////////////////////////////////////////////////////////////////////////
7 #include "Geometry/Geometry.h"
9 #include "Simulation/Particle.h"
11 #include "Utilities/func/MathUtil.h"
13 #include "Calibrator/Calibrator.h"
14 #include "MCCheater/BackTracker.h"
15 #include "RecoBase/RecoHit.h"
16 
26 #include "fhiclcpp/ParameterSet.h"
29 
30 #include "TLorentzVector.h"
31 #include "TTree.h"
32 #include "TH1I.h"
33 #include "TH1F.h"
34 
35 #include <iostream>
36 #include <string>
37 class TTree;
38 
39 namespace mcchk {
41  {
42  public:
43  explicit MonopoleAna(fhicl::ParameterSet const &pset);
44  virtual ~MonopoleAna();
45 
46  void analyze(const art::Event &evt);
47  void beginJob();
48  void reconfigure(fhicl::ParameterSet const &pset);
49 
50  private:
51 
52  TTree* MM_MC;
53  TTree* MMH_cheater;
54 
58 
59  double x0,y0,z0;
61  double pathLength; //pathlength inside FD
62  double mcTi, mcTf; //MC entry and exit time of monopole to FD
63  double rcTi, rcTf; //RC entry and exit time of monopole to FD
65  double dE, dPE;
66 
68  double entry_t;
69  double exit_t;
71  double dEnergy, chpl;
72  };
73 }
74 
75 
76 namespace mcchk{
77 
79  : EDAnalyzer(pset)
80  {
81  reconfigure(pset);
82  }
83 
84 
85  //......................................................................
87  {
88  fClusterInput = pset.get< std::string >("ClusterInput") ;
89  fFLSHitListLabel = pset.get< std::string >("FLSHitListLabel") ;
90  fRawDataLabel = pset.get< std::string >("RawDataLabel");
91  }
92 
93  //......................................................................
94 
96 
97 
98  //......................................................................
99 
101  {
103  std::cout<<"begin job!"<<std::endl;
104  //The following two ntuples can give info:
105  // efficiency VS track length
106 
107  MM_MC = tfs->make<TTree>("MM_MC","monopole_mc_info");
108  MM_MC -> Branch("x0",&x0,"x0/D");
109  MM_MC -> Branch("y0",&y0,"y0/D");
110  MM_MC -> Branch("z0",&z0,"z0/D");
111  MM_MC -> Branch("pathLength",&pathLength,"pathLength/D");
112  MM_MC -> Branch("cosx",&cosx,"cosx/D");
113  MM_MC -> Branch("cosy",&cosy,"cosy/D");
114  MM_MC -> Branch("cosz",&cosz,"cosz/D");
115  MM_MC -> Branch("momentum",&momentum,"momentum/D");
116  MM_MC -> Branch("mcTi",&mcTi,"mcTi/D");
117  MM_MC -> Branch("mcTf",&mcTf,"mcTf/D");
118  MM_MC -> Branch("rcTi",&rcTi,"rcTi/D");
119  MM_MC -> Branch("dE",&dE,"dE/D");
120  MM_MC -> Branch("dPE",&dPE,"dPE/D");
121  MM_MC -> Branch("nModule",&nModule,"nModule/I");
122  MM_MC -> Branch("nCell",&nCell,"nCell/I");
123  MM_MC -> Branch("nSCell",&nSCell,"nSCell/I");
124 
125 
126  MMH_cheater = tfs->make<TTree>("MMH_cheater","monopole_hits_mc_info");
127  MMH_cheater -> Branch("plane",&plane_id,"plane/I");
128  MMH_cheater -> Branch("cell",&cell_id,"cell/I");
129  MMH_cheater -> Branch("x0",&entry_x,"x0/D");
130  MMH_cheater -> Branch("y0",&entry_y,"y0/D");
131  MMH_cheater -> Branch("z0",&entry_z,"z0/D");
132  MMH_cheater -> Branch("t0",&entry_t,"t0/D");
133  MMH_cheater -> Branch("x1",&exit_x,"x1/D");
134  MMH_cheater -> Branch("y1",&exit_y,"y1/D");
135  MMH_cheater -> Branch("z1",&exit_z,"z1/D");
136  MMH_cheater -> Branch("t1",&exit_t,"t1/D");
137  MMH_cheater -> Branch("energydump",&dEnergy,"dEnergy/D");
138  MMH_cheater -> Branch("pathlength",&chpl,"pathlength/D");
139 
140  }
141 
143  {
144  //Initialization of MM_MC (each entry corresponds to a single MM):
145  mcTi = 9999999;
146  mcTf = -9999999;
147  pathLength = 0;
148  nCell = 0;
149  nSCell = 0;
150  nModule = 0;
151  dE =0;
152  dPE =0;
153  x0 = 0;
154  y0 = 0;
155  z0 =0;
156  momentum = 0;
157  cosx = 0;
158  cosy = 0;
159  cosz = 0;
160 
161 
162  //make a backtracker
164  const sim::ParticleNavigator& pnav = bt->ParticleNavigator();
165 
166  for (sim::ParticleNavigator::const_iterator i = pnav.begin(); i != pnav.end(); ++i){
167  const sim::Particle *p = (*i).second;
168  if (p->PdgCode()!=42 ) continue;
169 
170  momentum = p->P();
171  x0 = p->Vx();
172  y0 = p->Vy();
173  z0 = p->Vz();
174  cosx = (p->Px())/momentum;
175  cosy = (p->Py())/momentum;
176  cosz = (p->Pz())/momentum;
177 
178  const std::vector<sim::FLSHit>& flshits = bt->ParticleToFLSHit(p->TrackId());
179  for (unsigned h =0; h!= flshits.size(); ++h) {
180  plane_id = flshits[h].GetPlaneID();
181  cell_id = flshits[h].GetCellID();
182  entry_x = flshits[h].GetEntryX();
183  exit_x = flshits[h].GetExitX();
184  entry_y = flshits[h].GetEntryY();
185  exit_y = flshits[h].GetExitY();
186  entry_z = flshits[h].GetEntryZ();
187  exit_z = flshits[h].GetExitZ();
188  entry_t = flshits[h].GetEntryT();
189  exit_t = flshits[h].GetExitT();
190  dEnergy = flshits[h].GetEdep();
191  chpl = flshits[h].GetTotalPathLength();
192  MMH_cheater->Fill();
193 
194  pathLength += chpl;
195  if (mcTi>entry_t) mcTi = entry_t;
196  if (mcTf<exit_t) mcTf = exit_t;
197  }
198  }
199 
200  // define vector holding the clusters
202  evt.getByLabel(fClusterInput, slices);
204 
205  const int sliceMax = slices->size();
206  for(int sliceIdx = 0; sliceIdx < sliceMax; ++sliceIdx){
207 
208  const rb::Cluster& slice = (*slices)[sliceIdx];
209 
210  if (slice.IsNoise()) continue;
211 
212  nCell = slice.NCell();
213  nModule =0;
214  unsigned short plane,cell;
215  int pStore = -99;
216  int mStore = -99;
217 
218  for(int hitIdx = 0; hitIdx < nCell; ++hitIdx){
219  const art::Ptr<rb::CellHit> chit = slice.Cell(hitIdx);
220  plane = chit->Plane();
221  cell = chit->Cell();
222  int module = cell/32;
223  if (plane!=pStore or module!=mStore) nModule++;
224 
225  try {
226  if(!bt->IsNoise(chit) || (chit->ADC()>3100 && chit->TNS() > mcTi-1000 && chit->TNS() < mcTf+1000) ) {
227 
228  std::vector<cheat::TrackIDE> trackides= bt->HitToTrackIDE(chit);
229  TVector3 position = bt->HitToXYZ(chit);
230 
231  int mcPE=0;
232 
233  for (unsigned int i=0; i!=(bt->HitToPhotonSignal(chit)).size();++i)
234  mcPE += (bt->HitToPhotonSignal(chit))[i].NPhoton();
235 
236  double w = position[0];
237  if (chit->View() == geo::kX) w= position[1];
238 
239  rb::RecoHit recohit(calib->MakeRecoHit(*chit,w));
240 
241  const double rbTime = recohit.T();
242 
243  if (!hitIdx) rcTi = rbTime;
244  rcTf = rbTime;
245 
246  dE += (trackides[0].energy) / (trackides[0].energyFrac);
247  dPE += mcPE;
248 
249  if (chit->ADC(3)==4095) nSCell++;
250 
251  }
252 
253  pStore = plane;
254  mStore = module;
255 
256 
257  }
258 
259  catch (...) {continue;}
260 
261  }
262  }
263 
264  MM_MC->Fill();
265 
266  return;
267  }
268 
269 
270 
271 
272 
274 
275 }
float T() const
Definition: RecoHit.h:39
void analyze(const art::Event &evt)
float TNS() const
Definition: CellHit.h:46
back track the reconstruction to the simulation
std::string fFLSHitListLabel
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
double Py(const int i=0) const
Definition: MCParticle.h:230
Trajectory class.
Simple module to analyze MC cosmics distributions.
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
unsigned short Plane() const
Definition: CellHit.h:39
geo::View_t View() const
Definition: CellHit.h:41
const char * p
Definition: xmltok.h:285
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...
void reconfigure(fhicl::ParameterSet const &pset)
list_type::const_iterator const_iterator
Vertical planes which measure X.
Definition: PlaneGeo.h:28
double Px(const int i=0) const
Definition: MCParticle.h:229
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
A collection of associated CellHits.
Definition: Cluster.h:47
DEFINE_ART_MODULE(TestTMapFile)
bool IsNoise(const art::Ptr< rb::CellHit > &hit) const
Is this hit not associated with any particles?
int TrackId() const
Definition: MCParticle.h:209
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
unsigned short Cell() const
Definition: CellHit.h:40
CDPStorage service.
double P(const int i=0) const
Definition: MCParticle.h:233
T get(std::string const &key) const
Definition: ParameterSet.h:231
MonopoleAna(fhicl::ParameterSet const &pset)
int evt
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
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
double Vx(const int i=0) const
Definition: MCParticle.h:220
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const std::vector< sim::PhotonSignal > HitToPhotonSignal(const art::Ptr< rb::CellHit > &hit) const
Returns the PhotonSignals contributing the signal in the specified hit.
T * make(ARGS...args) const
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
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
double Vz(const int i=0) const
Definition: MCParticle.h:222
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
TVector3 HitToXYZ(art::Ptr< rb::CellHit > const &hit, bool useBirksE=false) const
Returns the XYZ position of the energy deposition for a given hit.
Float_t w
Definition: plot.C:20
double Vy(const int i=0) const
Definition: MCParticle.h:221
Encapsulate the geometry of one entire detector (near, far, ndos)