RockAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file RockAna.cxx
3 // \brief Analyse monopole monte carlo events
4 // \version $Id: RockAna.cxx,v 1.2 2012-10-09 03:56:03 ricken Exp $
5 // \author $Author: ricken $
6 ////////////////////////////////////////////////////////////////////////
7 #include "Geometry/Geometry.h"
10 #include "RecoBase/CellHit.h"
11 #include "RecoBase/Cluster.h"
14 #include "Utilities/func/MathUtil.h"
15 
16 
23 #include "fhiclcpp/ParameterSet.h"
29 
30 #include "TLorentzVector.h"
31 #include "TTree.h"
32 #include "TH1F.h"
33 
34 #include <iostream>
35 #include <string>
36 
37 class TNtuple;
38 class TTree;
39 
40 namespace mcchk {
41  class RockAna : public art::EDAnalyzer
42  {
43  public:
44  explicit RockAna(fhicl::ParameterSet const &pset);
45  virtual ~RockAna();
46 
47  void analyze(const art::Event &evt);
48  void beginJob();
49  void reconfigure(fhicl::ParameterSet const &pset);
50 
51  private:
52 
54 
55  TTree* Rock_MC;
56  TTree* Rock_Cheater;
57 
58  int nu_pdg;
59  float vx,vy,vz,px,py,pz,Enu;
60 
61  float mcPE,cx,cy,cz,ct;
62  int cview;
63 
64 
65  };
66 
67 
68 }
69 
70 
71 namespace mcchk{
72 
74  : EDAnalyzer(pset)
75  {
76  reconfigure(pset);
77  }
78 
80 
81  {
82  fClusterInput = pset.get< std::string >("ClusterInput") ;
83  }
84 
85 
86 
87 //......................................................................
88 
90 
91 
92 //......................................................................
93 
95  {
96 
98 
99  Rock_MC = tfs->make<TTree>("Rock_MC","Neutrino Info");
100  Rock_MC -> Branch("pdg",&nu_pdg,"pdg/I");
101  Rock_MC -> Branch("vx",&vx,"vx/F");
102  Rock_MC -> Branch("vy",&vy,"vy/F");
103  Rock_MC -> Branch("vz",&vz,"vz/F");
104  Rock_MC -> Branch("px",&px,"px/F");
105  Rock_MC -> Branch("py",&py,"py/F");
106  Rock_MC -> Branch("pz",&pz,"pz/F");
107  Rock_MC -> Branch("E", &Enu,"E/F");
108 
109  Rock_Cheater = tfs->make<TTree>("Rock_Cheater","all hits");
110  Rock_Cheater -> Branch("PE",&mcPE,"PE/F");
111  Rock_Cheater -> Branch("x",&cx,"x/F");
112  Rock_Cheater -> Branch("y",&cy,"y/F");
113  Rock_Cheater -> Branch("z",&cz,"z/F");
114  Rock_Cheater -> Branch("t",&ct,"t/F");
115  Rock_Cheater -> Branch("view",&cview,"view/I");
116 
117  }
118 
120  {
121 
123 
124 
125  //Now, get the vertex position and energy of neutrino
127  try{ evt.getByLabel("generator",mcTruths); }
128  catch(...){ std::cout<<"couldn't get the mcTruths"<<std::endl; return;}
129 
130  for (unsigned idx =0; idx!=mcTruths->size();++idx) {
131  art::Ptr<simb::MCTruth> mc (mcTruths,idx);
132 
133  const simb::MCNeutrino& mc_neutrino = mc->GetNeutrino();
134  const simb::MCParticle& nu = mc_neutrino.Nu();
135  nu_pdg = nu.PdgCode();
136  vx = nu.Vx();
137  vy = nu.Vy();
138  vz = nu.Vz();
139  px = nu.Px();
140  py = nu.Py();
141  pz = nu.Pz();
142  Enu = nu.P();
143 
144  Rock_MC->Fill();
145  }
146 
147  //make a backtracker
149  // const sim::ParticleNavigator& pnav = bt->ParticleNavigator();
150 
151 
152  // now we collect all hits from the slice
154  try{ evt.getByLabel(fClusterInput, slices);}
155  catch(...){ std::cout<<"couldn't get slices"<<std::endl; return;}
156  const int sliceMax = slices->size();
157 
158  for(int sliceIdx = 0; sliceIdx < sliceMax; ++sliceIdx){
159  const rb::Cluster& slice = (*slices)[sliceIdx];
160  const int hitMax = slice.NCell();
161 
162  for (int hitIdx=0; hitIdx<hitMax; ++hitIdx)
163  {
164  const art::Ptr<rb::CellHit> chit = slice.Cell(hitIdx);
165 
166  unsigned short plane = chit->Plane();
167  unsigned short cell = chit->Cell();
168  if (chit->View() == geo::kX) cview = 0;
169  else cview =1;
170  try {
171  if (!(bt->IsNoise(chit))){
172  //Get the center position of each hit
173  double xyz[3];
174  geom->Plane(plane)->Cell(cell)->GetCenter(xyz);
175 
176  cx = xyz[0];
177  cy = xyz[1];
178  cz = xyz[2];
179 
180  mcPE=0;
181  ct =0;
182  for (unsigned int i=0; i!=(bt->HitToPhotonSignal(chit)).size();++i){
183  mcPE += (bt->HitToPhotonSignal(chit))[i].NPhoton();
184  ct += (bt->HitToPhotonSignal(chit))[i].TimeMean();
185  }
186  ct = ct/((bt->HitToPhotonSignal(chit)).size()+0.0);
187 
188 
189  Rock_Cheater->Fill();
190 
191  }
192 
193 
194  }
195 
196 
197  catch (...) {
198  continue;
199  }
200 
201 
202  }
203 
204 
205  }
206 
207 
208 
209 
210  return;
211  }
212 
214 
215 }
void analyze(const art::Event &evt)
std::string fClusterInput
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
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:76
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
double Py(const int i=0) const
Definition: MCParticle.h:230
Simple module to analyze MC cosmics distributions.
void reconfigure(fhicl::ParameterSet const &pset)
unsigned short Plane() const
Definition: CellHit.h:39
geo::View_t View() const
Definition: CellHit.h:41
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
Vertical planes which measure X.
Definition: PlaneGeo.h:28
double Px(const int i=0) const
Definition: MCParticle.h:229
A collection of associated CellHits.
Definition: Cluster.h:47
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
bool IsNoise(const art::Ptr< rb::CellHit > &hit) const
Is this hit not associated with any particles?
unsigned short Cell() const
Definition: CellHit.h:40
double P(const int i=0) const
Definition: MCParticle.h:233
T get(std::string const &key) const
Definition: ParameterSet.h:231
int evt
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
virtual ~RockAna()
RockAna(fhicl::ParameterSet const &pset)
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
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
void geom(int which=0)
Definition: geom.C:163
double Pz(const int i=0) const
Definition: MCParticle.h:231
TTree * Rock_Cheater
double Vz(const int i=0) const
Definition: MCParticle.h:222
Event generator information.
Definition: MCNeutrino.h:18
double Vy(const int i=0) const
Definition: MCParticle.h:221
Encapsulate the geometry of one entire detector (near, far, ndos)