FmmTriggerAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: FmmTriggerAna
3 // Module Type: analyzer
4 // File: FmmTriggerAna_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 ////////////////////////////////////////////////////////////////////////
14 
15 #include "MCCheater/BackTracker.h"
16 #include "Calibrator/Calibrator.h"
17 #include "Simulation/Particle.h"
21 #include "Utilities/func/MathUtil.h"
22 #include "fhiclcpp/ParameterSet.h"
24 
25 #include "Geometry/Geometry.h"
27 #include "GeometryObjects/Geo.h"
29 
30 #include "NovaDAQConventions/DAQConventions.h"
31 #include "RecoBase/CellHit.h"
32 #include "RecoBase/RecoHit.h"
33 #include "RecoBase/Cluster.h"
34 #include "RecoBase/Track.h"
35 
36 #include "TVector3.h"
37 #include "TTree.h"
38 #include "TProfile.h"
39 
40 #include <algorithm>
41 #include <utility>
42 #include <cmath>
43 #include <iostream>
44 
45 namespace zcl {
46  class FmmTriggerAna;
47 }
48 
50  public:
51  explicit FmmTriggerAna(fhicl::ParameterSet const & pset);
52  virtual ~FmmTriggerAna();
53 
54  void analyze(const art::Event & evt);
55  void beginJob();
56  void endJob(){};
57 
58  private:
59  std::string fClusterInput; ///< Input folder from cluster reco
60 
61  //Each entry is corresponding to 1 reconstructed slice
62  TTree * RC_Tree;
63  //TProfile * hitTdiff;
64  //TProfile * hitTdiff_log;
65  //TProfile * DCS_beta;
66  //TProfile * DCS_logbeta;
67 
69 
70  //Trigger decision branches
71  int _sliced;
73  int _trueHits;
74  int _nSatHits;
75  double _meanADC;
76  int _deltaP;
77  double _deltaTNS;
79 
80  double _tmin, _tmax; //The time window of the slice
81  int _pmin, _pmax;
82 
83  double beta_MC;
85 
86 };
87 
89  EDAnalyzer(p),
90  fClusterInput (p.get< std::string > ("ClusterInput") )
91 {
92 }
93 
95 {
96  // Clean up dynamic memory and other resources here.
97 }
98 
100 {
102 
103  //Combining cheater info and the RC info to make the data structure better
104  //One RC track info per entry, if no tracks reconstructed in this event, fill in 1 entry anyway tagging the
105  //reconstruction failed
106  RC_Tree = tfs->make<TTree>("RC_Tree","Information from RC track info");
107  RC_Tree -> Branch("event", &evtID, "event/I");
108  RC_Tree -> Branch("subRun", &subRunID, "subRun/I");
109  RC_Tree -> Branch("Run", &RunID, "Run/I");
110 
111  //Trigger threshold parameters
112  RC_Tree -> Branch("sliced", &_sliced, "sliced/I");
113  RC_Tree -> Branch("penetrated", &_penetrated, "penetrated/O");
114  RC_Tree -> Branch("nSatHits", &_nSatHits, "nSatHits/I");
115  RC_Tree -> Branch("meanADC", &_meanADC, "meanADC/D");
116  RC_Tree -> Branch("deltaP", &_deltaP, "deltaP/I");
117  RC_Tree -> Branch("deltaTNS", &_deltaTNS, "deltaTNS/D");
118  RC_Tree -> Branch("nxhS", &_nxhits_slice, "nxhS/I");
119  RC_Tree -> Branch("nyhS", &_nyhits_slice, "nyhS/I");
120  RC_Tree -> Branch("trueHits", &_trueHits, "trueHits/I");
121 
122  //Cheater info ntuple part
123  RC_Tree -> Branch("betaMC", &beta_MC, "betaMC/D");
124 }
125 
127 {
128  evtID = evt.id().event();
129  subRunID = evt.subRun();
130  RunID = evt.run();
131 
135 
136 
137  //Get MC info first:
138  const sim::ParticleNavigator& pnav = bt->ParticleNavigator();
139  for (sim::ParticleNavigator::const_iterator i = pnav.begin(); i != pnav.end(); ++i) {
140  nxhits_MC = 0;
141  nyhits_MC = 0;
142  const sim::Particle *p = (*i).second;
143  if (p->PdgCode()!=42 ) continue;
144  const std::vector<sim::FLSHit>& flshits = bt->ParticleToFLSHit(p->TrackId());
145  unsigned ntruehits_Tot = flshits.size();
146 
147  if ( ntruehits_Tot==0 ) continue;
148  double momentum2 = (p->P())*(p->P());
149  double mass2 = (p->Mass())*(p->Mass());
150  beta_MC = sqrt(momentum2 / (momentum2+mass2) );
151 
152 
153  for (unsigned h = 0; h!= flshits.size(); ++h) {
154  if (flshits[h].GetPlaneID() % 2) ++nxhits_MC;
155  else ++nyhits_MC;
156  }
157 
158  //1 monopole per event, so you can break out from the loop:
159  break;
160  }
161 
162  //=====================================================================================
163  //Now, start pseudo trigger:
164  const int _xMax = 383;
165  const int _yMax = 383;
166  const int _zMax = 895;
167  const int _xdelt = 35;
168  const int _ydelt = 35;
169  const int _zdelt = 15;
170 
171  int PX_max, PY_max, PX_min, PY_min;
172  double TX_max, TY_max, TX_min, TY_min;
173 
175  evt.getByLabel(fClusterInput, slices);
176 
177  _sliced = slices->size();
178  // std::cout<<"There are "<<slices->size()<<" slices in the event."<<std::endl;
179  //If no slice reconstructed, just return the MC info
180  if (_sliced == 0) {
181  _deltaP = 0;
182  _deltaTNS = 0;
183  _nxhits_slice = 0;
184  _nyhits_slice = 0;
185  _penetrated = false;
186  _trueHits = 0;
187  _nSatHits = 0;
188  _meanADC = 0;
189  RC_Tree->Fill();
190  }
191 
192  //Instead of splitting the hits into 2 views, just directly take the entire 3D slice:
193  for (unsigned sliceIdx = 0; sliceIdx != slices->size(); ++sliceIdx) {
194  const rb::Cluster& slice = (*slices)[sliceIdx];
195 
196  _nxhits_slice = 0;
197  _nyhits_slice = 0;
198 
199  PX_max = 0;
200  PY_max = 0;
201  PX_min = 999;
202  PY_min = 999;
203  TX_max = -9e7;
204  TY_max = -9e7;
205  TX_min = 9e9;
206  TY_min = 9e9;
207 
208  _penetrated = false;
209  _trueHits = 0;
210  _nSatHits = 0;
211  _meanADC = 0;
212 
213  bool Intrusion = false;
214  int entry_id = -1;
215  for (unsigned i = 0; i != slice.NCell(); ++i) {
216  const art::Ptr<rb::CellHit> chit = slice.Cell(i);
217  if (chit->ADC() < 501) continue;
218  if (chit->ADC() > 3100) ++_nSatHits;
219  _meanADC += chit->ADC();
220  float tns = chit->TNS();
221 
222  // int c = chit->Cell();
223  // int p = chit->Plane();
224 
225  if(chit->IsMC()==true) ++_trueHits;
226 
227  if (chit->View()==geo::kX) { //X view
228  ++_nxhits_slice;
229  if (PX_max < chit->Plane() ) PX_max = chit->Plane();
230  if (PX_min > chit->Plane() ) PX_min = chit->Plane();
231  if (TX_max < tns ) TX_max = tns;
232  if (TX_min > tns ) TX_min = tns;
233 
234  if (!Intrusion) {
235  if (chit->Cell() < _xdelt){
236  entry_id = 1;
237  Intrusion = true;
238  }
239  else if (chit->Cell() > _xMax - _xdelt){
240  entry_id = 2;
241  Intrusion = true;
242  }
243  }
244  else if (!_penetrated) {
245  if (chit->Cell() < _xdelt && entry_id!=1) _penetrated = true;
246  else if (chit->Cell() > _xMax-_xdelt && entry_id!=2) _penetrated = true;
247  }
248  }
249 
250  else { //Y view
251  ++_nyhits_slice;
252  if (PY_max < chit->Plane() ) PY_max = chit->Plane();
253  if (PY_min > chit->Plane() ) PY_min = chit->Plane();
254  if (TY_max < chit->TNS() ) TY_max = chit->TNS();
255  if (TY_min > chit->TNS() ) TY_min = chit->TNS();
256 
257  if (!Intrusion) {
258  if (chit->Cell() < _ydelt) {
259  entry_id = 3;
260  Intrusion = true;
261  }
262  else if (chit->Cell() > _yMax - _ydelt) {
263  entry_id = 4;
264  Intrusion = true;
265  }
266  }
267 
268  else if (!_penetrated) {
269  if (chit->Cell() < _ydelt && entry_id!=3) _penetrated = true;
270  else if (chit->Cell() > _yMax-_ydelt && entry_id!=4) _penetrated = true;
271  }
272  }
273 
274  if (!Intrusion) {
275  if (chit->Plane() < _zdelt) {
276  entry_id = 5;
277  Intrusion = true;
278  }
279  else if (chit->Plane() > _zMax - _zdelt) {
280  entry_id = 6;
281  Intrusion = true;
282  }
283  }
284  else if (!_penetrated) { //intrusion happened but not through yet
285  if (chit->Plane() < _zdelt && entry_id!=5) _penetrated=true;
286  else if (chit->Plane() > _zMax-_zdelt && entry_id!=6) _penetrated=true;
287  }
288  }
289  //Record the info
291  _deltaP = std::min(PX_max,PY_max) - std::max(PX_min,PY_min);
292  _deltaTNS = std::min(TX_max,TY_max) - std::max(TX_min,TY_min);
293  _tmin = std::min(TX_min,TY_min);
294  _tmax = std::max(TX_max,TY_max);
295  _pmin = std::min(PX_min,PY_min);
296  _pmax = std::max(PX_max,PY_max);
297 
298  // if (_nxhits_slice+_nyhits_slice < 10) _sliced = false;
299  // std::cout<<"xhits: "<<_nxhits_slice<<" yhits: "<<_nyhits_slice<<std::endl;
300 
301  // std::cout<<"There are "<<_trueHits<<" true hits in No."<<sliceIdx<<" slice in the event."<<std::endl;
302  //if(_trueHits>2) RC_Tree->Fill();
303  RC_Tree->Fill();
304  }
305 
306 } //end analyzing
307 
308 
float TNS() const
Definition: CellHit.h:46
T max(const caf::Proxy< T > &a, T b)
SubRunNumber_t subRun() const
Definition: Event.h:72
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 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
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
A collection of associated CellHits.
Definition: Cluster.h:47
DEFINE_ART_MODULE(TestTMapFile)
int TrackId() const
Definition: MCParticle.h:209
unsigned short Cell() const
Definition: CellHit.h:40
CDPStorage service.
double P(const int i=0) const
Definition: MCParticle.h:233
int evt
std::string fClusterInput
Input folder from cluster reco.
bool IsMC() const
Definition: RawDigit.h:108
Collect Geo headers and supply basic geometry functions.
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
FmmTriggerAna(fhicl::ParameterSet const &pset)
EventNumber_t event() const
Definition: EventID.h:116
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
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 analyze(const art::Event &evt)
void geom(int which=0)
Definition: geom.C:163
T min(const caf::Proxy< T > &a, T b)
RunNumber_t run() const
Definition: Event.h:77
EventID id() const
Definition: Event.h:56
Encapsulate the geometry of one entire detector (near, far, ndos)