MMCheater_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: MMCheater
3 // Module Type: analyzer
4 // File: MMCheater_module.cc
5 //
6 // Generated at Wed Dec 18 15:41:32 2013 by Zukai Wang using artmod
7 // from cetpkgsupport v1_04_02.
8 ////////////////////////////////////////////////////////////////////////
15 
16 #include "MCCheater/BackTracker.h"
17 #include "Simulation/Particle.h"
21 #include "fhiclcpp/ParameterSet.h"
22 
23 
24 #include "Geometry/Geometry.h"
27 
28 #include "NovaDAQConventions/DAQConventions.h"
29 #include "RecoBase/CellHit.h"
30 #include "RecoBase/Cluster.h"
31 
32 #include "TH3.h"
33 #include "TVector3.h"
34 #include "TTree.h"
35 
36 #include <algorithm>
37 #include <utility>
38 #include <cmath>
39 #include <iostream>
40 
41 namespace zcl {
42  class MMCheater;
43 }
44 
46 public:
47  explicit MMCheater(fhicl::ParameterSet const & pset);
48  virtual ~MMCheater();
49 
50  void beginJob();
51  double distance(float c1, float c2, float p1, float p2);
52  void analyze(const art::Event & evt);
53 
54 private:
56 
57  TTree* Cheater_Tree;
58  double beta, cosx, cosy, cosz;
59  double edep, pl, tof;
60  int NXhits, NYhits;
62  int ADC_min;
63  double ADC_mean;
67 
68 
69 };
70 
72  EDAnalyzer(p),
73  fHitInput (p.get< std::string > ("CellHitInput") )
74 {
75 }
76 
78 {
79  // Clean up dynamic memory and other resources here.
80 }
81 
83 {
85  Cheater_Tree = tfs->make<TTree>("Cheater_Tree","Information from MC Cheater on event basis");
86  Cheater_Tree -> Branch("beta", &beta, "beta/D");
87  Cheater_Tree -> Branch("cosx", &cosx, "cosx/D");
88  Cheater_Tree -> Branch("cosy", &cosy, "cosy/D");
89  Cheater_Tree -> Branch("cosz", &cosz, "cosz/D");
90  Cheater_Tree -> Branch("edep", &edep, "edep/D");
91  Cheater_Tree -> Branch("PathLength", &pl, "PathLength/D");
92  Cheater_Tree -> Branch("tof", &tof, "tof/D");
93  Cheater_Tree -> Branch("nxhits", &NXhits, "nxhits/I");
94  Cheater_Tree -> Branch("nyhits", &NYhits, "nyhits/I");
95  Cheater_Tree -> Branch("D_TDC_max", &D_TDC_max, "D_TDC_max/I");
96  Cheater_Tree -> Branch("D_Z_max", &D_Z_max, "D_Z_max/I");
97  Cheater_Tree -> Branch("ADC_min", &ADC_min, "ADC_min/I");
98  Cheater_Tree -> Branch("ADC_mean", &ADC_mean, "ADC_mean/D");
99  Cheater_Tree -> Branch("D_P_min", &D_P_min, "D_P_min/I");
100  Cheater_Tree -> Branch("D_C_min", &D_C_min, "D_C_min/I");
101  Cheater_Tree -> Branch("Lx", &Lx, "Lx/D"); //In DDT unit
102  Cheater_Tree -> Branch("D_lx_max", &D_lx_max, "D_lx_max/D");
103  Cheater_Tree -> Branch("vpro_x", &vpro_x, "vpro_x/D");
104  Cheater_Tree -> Branch("D_Tx_max_l", &D_Tx_max_l, "D_Tx_max_l/D");
105  Cheater_Tree -> Branch("D_Tx_max_s", &D_Tx_max_s, "D_Tx_max_s/D");
106  Cheater_Tree -> Branch("Ly", &Ly, "Ly/D");
107  Cheater_Tree -> Branch("D_ly_max", &D_ly_max, "D_ly_max/D");
108  Cheater_Tree -> Branch("vpro_y", &vpro_y, "vpro_y/D");
109  Cheater_Tree -> Branch("D_Ty_max_l", &D_Ty_max_l, "D_Ty_max_l/D");
110  Cheater_Tree -> Branch("D_Ty_max_s", &D_Ty_max_s, "D_Ty_max_s/D");
111 }
112 
114 {
117  //Get MC info first:
118  const sim::ParticleNavigator& pnav = bt->ParticleNavigator();
119  for (sim::ParticleNavigator::const_iterator i = pnav.begin(); i != pnav.end(); ++i) {
120  const sim::Particle *p = (*i).second;
121  if (p->PdgCode()!=42 ) continue;
122  const std::vector<sim::FLSHit>& flshits = bt->ParticleToFLSHit(p->TrackId());
123 
124  if (flshits.size()==0 ) continue;
125  cosx = (p->Px())/(p->P());
126  cosy = (p->Py())/(p->P());
127  cosz = (p->Pz())/(p->P());
128  double momentum2 = (p->P())*(p->P());
129  double mass2 = (p->Mass())*(p->Mass());
130  beta = sqrt(momentum2 / (momentum2+mass2) );
131 
132  double ti_MC = flshits[0].GetEntryT();
133  double tf_MC = flshits[0].GetExitT();
134 
135  edep = flshits[0].GetEdep();
136  pl = flshits[0].GetTotalPathLength();
137 
138  for (unsigned h =1; h!= flshits.size(); ++h) {
139  double ti = flshits[h].GetEntryT();
140  double tf = flshits[h].GetExitT();
141  if (ti < ti_MC)
142  ti_MC = ti;
143 
144  if (tf > tf_MC)
145  tf_MC = tf;
146 
147  pl += flshits[h].GetTotalPathLength();
148  edep += flshits[h].GetEdep();
149  }
150  tof = tf_MC-ti_MC;
151 
152  //1 monopole per event, so you can break out from the loop:
153  break;
154  }
155 
156  // std::cout<<"particle info collected: FLS finish!"<<std::endl;
157 
158  std::vector<art::Ptr<rb::CellHit> > TrueHitsList;
159  std::vector<art::Ptr<rb::CellHit> > Xhits;
160  std::vector<art::Ptr<rb::CellHit> > Yhits;
162  evt.getByLabel(fHitInput, hitcol);
163 
164  for (unsigned int i = 0; i != hitcol->size(); ++i) {
165  art::Ptr<rb::CellHit> hit_ptr(hitcol, i);
166  if (!bt->IsNoise(hit_ptr))
167  TrueHitsList.emplace_back(hit_ptr);
168  }
169 
170  if (TrueHitsList.size()<6) return;
171  rb::SortByTime(TrueHitsList);
172 
173  D_TDC_max = 0;
174  D_Z_max = 0;
175  ADC_mean = 0;
176  ADC_min = 4096;
177  D_P_min = 888;
178  D_C_min = 384;
179  int plane_tmp = TrueHitsList[0]->Plane();
180  int TDC_tmp = TrueHitsList[0]->TDC();
181 
182  int xti,xci,xpi,xtf,xcf,xpf;
183  int yti,yci,ypi,ytf,ycf,ypf;
184  xci = 384;
185  yci = 384;
186  xcf = -1;
187  ycf = -1;
188  xpi = 888;
189  xpf = -1;
190  ypi = 888;
191  ypf = -1;
192  xti = TDC_tmp;
193  yti = TDC_tmp;
194  xtf = TDC_tmp;
195  ytf = TDC_tmp;
196 
197  for (unsigned i =0; i!= TrueHitsList.size(); ++i) {
198  art::Ptr<rb::CellHit> hit_ptr = TrueHitsList[i];
199  int ADC = hit_ptr->ADC();
200  int plane = hit_ptr->Plane();
201  int cell = hit_ptr->Cell();
202  int TDC = hit_ptr->TDC();
203 
204  if (hit_ptr->View()==geo::kX) {
205  Xhits.emplace_back(hit_ptr);
206  if (cell<xci) {
207  xpi = plane;
208  xci = cell;
209  xti = TDC;
210  }
211  if (cell>xcf) {
212  xpf = plane;
213  xcf = cell;
214  xtf = TDC;
215  }
216  }
217  else {
218  Yhits.emplace_back(hit_ptr);
219  if (cell<yci) {
220  ypi = plane;
221  yci = cell;
222  yti = TDC;
223  }
224  if (cell>ycf) {
225  ypf = plane;
226  ycf = cell;
227  ytf = TDC;
228  }
229  }
230  int delta_tdc = abs(TDC - TDC_tmp);
231  TDC_tmp = TDC;
232  if (delta_tdc > D_TDC_max) D_TDC_max = delta_tdc;
233 
234  int delta_z = abs(plane - plane_tmp);
235  plane_tmp = plane;
236  if (delta_z > D_Z_max) D_Z_max = delta_z;
237 
238  ADC_mean += ADC;
239  if (ADC_min > ADC) ADC_min = ADC;
240 
241  int delta_p = std::min(887-plane, plane);
242  int delta_c = std::min(383-cell, cell);
243  if (D_P_min > delta_p) D_P_min = delta_p;
244  if (D_C_min > delta_c) D_C_min = delta_c;
245  }
246 
247  ADC_mean = ADC_mean/( TrueHitsList.size()+0.);
248 
249  NXhits = Xhits.size();
250  NYhits = Yhits.size();
251 
252  // std::cout<<"Cosmic Slice info collected!"<<std::endl;
253 
254  if (NXhits<3 || NYhits<3) return;
255 
256  //OK, now the slow MM rc and trigger tunes:
257  // rb::SortByPlaneAndCell(Xhits);
258  Lx = distance(xci, xcf, xpi, xpf);
259  // std::cout<<"xci: "<<xci<<" xcf: "<<xcf<<" xpi: "<<xpi<<" xpf: "<<xpf<<std::endl;
260  vpro_x = (xtf-xti+0.)/Lx;
261  D_lx_max = 0;
262  D_Tx_max_l = 0;
263  D_Tx_max_s = 0;
264  for (unsigned i = 1; i!= Xhits.size()-1; ++i) {
265  auto hit_ptr = Xhits[i];
266  int pi = hit_ptr->Plane();
267  int ci = hit_ptr->Cell();
268  int ti = hit_ptr->TDC();
269  double dl1 = distance(xci, ci, xpi, pi);
270  double dl2 = distance(xcf, ci, xpf, pi);
271  if (dl1+dl2-Lx > D_lx_max) D_lx_max=dl1+dl2-Lx;
272 
273  double dtl, dts;
274  if (dl1>dl2) {
275  dtl = fabs(ti-(xti+vpro_x*dl1));
276  dts = fabs(ti-(xtf-vpro_x*dl2));
277  }
278 
279  else {
280  dts = fabs(ti-(xti+vpro_x*dl1));
281  dtl = fabs(ti-(xtf-vpro_x*dl2));
282  }
283 
284  if (dts>D_Tx_max_s) D_Tx_max_s = dts;
285  if (dtl>D_Tx_max_l) D_Tx_max_l = dtl;
286  }
287  // std::cout<<"X SMM finished!"<<std::endl;
288  //Yview :
289  // rb::SortByPlaneAndCell(Yhits);
290  Ly = distance(yci, ycf, ypi, ypf);
291  vpro_y = (ytf-yti+0.)/Ly;
292  D_ly_max = 0;
293  D_Ty_max_l = 0;
294  D_Ty_max_s = 0;
295  for (unsigned i = 1; i!= Yhits.size()-1; ++i) {
296  auto hit_ptr = Yhits[i];
297  int pi = hit_ptr->Plane();
298  int ci = hit_ptr->Cell();
299  int ti = hit_ptr->TDC();
300  double dl1 = distance(yci, ci, ypi, pi);
301  double dl2 = distance(ycf, ci, ypf, pi);
302  if (dl1+dl2-Ly > D_ly_max) D_ly_max=dl1+dl2-Ly;
303 
304  double dtl, dts;
305  if (dl1>dl2) {
306  dtl = fabs(ti-(yti+vpro_y*dl1));
307  dts = fabs(ti-(ytf-vpro_y*dl2));
308  }
309 
310  else {
311  dts = fabs(ti-(yti+vpro_y*dl1));
312  dtl = fabs(ti-(ytf-vpro_y*dl2));
313  }
314 
315  if (dts>D_Ty_max_s) D_Ty_max_s = dts;
316  if (dtl>D_Ty_max_l) D_Ty_max_l = dtl;
317  }
318 
319  // std::cout<<"Delta TDC: "<<D_TDC_max<<std::endl;
320  Cheater_Tree->Fill();
321  return;
322 } //end production function
323 
324 double zcl::MMCheater::distance(float c1, float c2,
325  float p1, float p2)
326 {
327  return sqrt((p1-p2)*(p1-p2)*2.82072+(c1-c2)*(c1-c2));
328 }
329 
back track the reconstruction to the simulation
int PdgCode() const
Definition: MCParticle.h:211
double Py(const int i=0) const
Definition: MCParticle.h:230
int32_t TDC() const
The time of the last baseline sample.
Definition: RawDigit.h:94
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
Definition: BackTracker.h:744
unsigned short Plane() const
Definition: CellHit.h:39
MMCheater(fhicl::ParameterSet const &pset)
geo::View_t View() const
Definition: CellHit.h:41
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
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...
double Mass() const
Definition: MCParticle.h:238
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
void abs(TH1 *hist)
void SortByTime(std::vector< art::Ptr< rb::CellHit > > &c)
Sort c in time order (earliest to latest).
Definition: CellHit.cxx:134
DEFINE_ART_MODULE(TestTMapFile)
Timing fit.
bool IsNoise(const art::Ptr< rb::CellHit > &hit) const
Is this hit not associated with any particles?
int TrackId() const
Definition: MCParticle.h:209
double distance(float c1, float c2, float p1, float p2)
c2
Definition: demo5.py:33
unsigned short Cell() const
Definition: CellHit.h:40
double P(const int i=0) const
Definition: MCParticle.h:233
int evt
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
virtual ~MMCheater()
T * make(ARGS...args) const
void analyze(const art::Event &evt)
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
void geom(int which=0)
Definition: geom.C:163
double Pz(const int i=0) const
Definition: MCParticle.h:231
std::string fHitInput
c1
Definition: demo5.py:24
T min(const caf::Proxy< T > &a, T b)
Encapsulate the geometry of one entire detector (near, far, ndos)
enum BeamMode string