NDRockMuon_module.cc
Go to the documentation of this file.
1 ///
2 /// \brief
3 /// \author Xuebing Bu - xbbu@fnal.gov
4 ///
5 /// select rock muons in ND to study DQ and calibration
6 ///
7 
9 #include "Geometry/Geometry.h"
11 
12 // ROOT includes
13 #include "TVector3.h"
14 #include "TH1F.h"
15 #include "TMath.h"
16 #include "TTree.h"
17 #include "TProfile2D.h"
18 #include "time.h"
19 #include "TString.h"
20 
21 #include "RecoBase/CellHit.h"
22 #include "RecoBase/Cluster.h"
23 #include "RecoBase/RecoHit.h"
24 #include "RecoBase/Track.h"
25 #include "RecoBase/Vertex.h"
26 #include "RecoBase/Prong.h"
27 
28 #include "Utilities/AssociationUtil.h"
29 #include "Utilities/func/MathUtil.h"
30 
31 // Framework includes
40 #include "fhiclcpp/ParameterSet.h"
46 
47 #include <cmath>
48 #include <vector>
49 #include <string>
50 #include <vector>
51 #include <sstream>
52 
53 class TVector3;
54 class TH1F;
55 class TTree;
56 
57 namespace rb{class Cluster; class Track;}
58 
59 namespace rockmuon {
60  class NDRockMuon : public art::EDAnalyzer {
61  public:
62  explicit NDRockMuon(fhicl::ParameterSet const& pset);
63  virtual ~NDRockMuon();
64 
65  void beginJob();
66  void analyze(const art::Event& evt);
67  void reconfigure(const fhicl::ParameterSet& p);
68  void endJob();
69 
70  int Nevts_percell[214][96];
71 
72  private:
73 
74  //folders to get objects from
75  std::string fSliceLabel; ///label for slices
76  std::string fKTrackLabel; ///label for kalman tracks
77  std::string fCellHitLabel; ///label for cell hits
78  std::string fVertexLabel; ///label for vertex
79 
80  //tunable parameters
81  int fMinHit; //slices must have at least this many hits
82  float fMinTrackLength; //tracks must have at least this length in cm
83  float fMinTrackCosTheta; //tracks must have at least this cosTheta
84  float fMaxTrackStartZ; //tracks must start within this Z in cm
85  float fMaxVertexX; //verties must be within this X
86  float fMaxVertexY; //verties must be within this Y
87 
88  //define histograms
89  TH1F *vertexX;
90  TH1F *vertexY;
91  TH1F *vertexZ;
92 
93  TH1F *trackStartX;
94  TH1F *trackStartY;
95  TH1F *trackStartZ;
96  TH1F *trackLength;
97 
98  TH1F *Ncells_track[6]; //# of track associated cells in all(0), 3DBs(1), DB1(2), DB2(3),DB3(4), and muon-catcher(5);
99  TH1F *Ecells_track[6]; //sum energy of track associated cells
100 
101  TH1F *PECorr_track[6];//PECorr distributions
102  TH1F *PECorr_track_MCH;//PECorr distributions for MC H modules
103  TH1F *PECorr_track_MCV;//PECorr distributions for MC V modules
104 
105  TH1F *Ecell_track[6];//Ecell distributions
106  TH1F *Ecell_track_MCH;//Ecell distributions for MC H modules
107  TH1F *Ecell_track_MCV;//Ecell distributions for MC V modules
108 
109  TH1F *planeEnergy[214];//energy per plane
110  TH1F *moduleADC[214][3];//ADC for each module
111  TH1F *modulePECorr[214][3];//PECorr for each module
112  TH1F *modulePE[214][3];//PE for each module
113  TH1F *moduleEnergy[214][3];//energy for each module
114  };
115 }
116 
117 namespace rockmuon{
118 
119  //......................................................................
120  NDRockMuon::NDRockMuon(fhicl::ParameterSet const& pset)
121  : EDAnalyzer(pset)
122  {
123  mf::LogInfo("NDRockMuon")<<__PRETTY_FUNCTION__<<"\n";
124  reconfigure(pset);
125  }
126 
127  //......................................................................
129  {
130  }
131 
132  //......................................................................
134  {
135  #define SET(T,N) { f##N = pset.get<T>(#N); }
136  SET(std::string, SliceLabel);
137  SET(std::string, CellHitLabel);
138  SET(std::string, KTrackLabel);
139  SET(std::string, VertexLabel);
140 
141  SET(int, MinHit);
142  SET(float, MinTrackLength);
143  SET(float, MinTrackCosTheta);
144  SET(float, MaxTrackStartZ);
145  SET(float, MaxVertexX);
146  SET(float, MaxVertexY);
147  #undef SET
148  }
149 
150  //......................................................................
152  {
154 
155  for( int ip=0; ip<214; ++ip ){
156  for( int ic=0; ic<96; ++ic ){
157  Nevts_percell[ip][ic]=0;
158  }
159  }
160 
161  vertexX=tfs->make<TH1F>("vertexX","",40,-200,200);
162  vertexY=tfs->make<TH1F>("vertexY","",40,-200,200);
163  vertexZ=tfs->make<TH1F>("vertexZ","",100,0,100);
164  vertexX->Sumw2();
165  vertexY->Sumw2();
166  vertexZ->Sumw2();
167 
168  trackStartX=tfs->make<TH1F>("trackStartX","",40,-200,200);
169  trackStartY=tfs->make<TH1F>("trackStartY","",40,-200,200);
170  trackStartZ=tfs->make<TH1F>("trackStartZ","",100,0,100);
171  trackStartX->Sumw2();
172  trackStartY->Sumw2();
173  trackStartZ->Sumw2();
174 
175  trackLength=tfs->make<TH1F>("trackLength","",160,0,1600);
176  trackLength->Sumw2();
177 
178  TString detname[6]={"all","DB123","DB1","DB2","DB3","MC"};
179  for( int i=0; i<6; ++i ){
180  TString name1="Ncells_track_"+detname[i];
181  Ncells_track[i]=tfs->make<TH1F>(name1,"",300,0,300);
182  Ncells_track[i]->Sumw2();
183 
184  TString name2="Ecells_track_"+detname[i];
185  Ecells_track[i]=tfs->make<TH1F>(name2,"",300,0,3);
186  Ecells_track[i]->Sumw2();
187 
188  TString name3="PECorr_track_"+detname[i];
189  PECorr_track[i]=tfs->make<TH1F>(name3,"",5000,0,5000);
190  PECorr_track[i]->Sumw2();
191 
192  TString name4="Ecell_track_"+detname[i];
193  Ecell_track[i]=tfs->make<TH1F>(name4,"",1000,0,0.1);
194  Ecell_track[i]->Sumw2();
195  }//loop det parts
196 
197  PECorr_track_MCH=tfs->make<TH1F>("PECorr_track_MCH","",5000,0,5000);
198  PECorr_track_MCH->Sumw2();
199  PECorr_track_MCV=tfs->make<TH1F>("PECorr_track_MCV","",5000,0,5000);
200  PECorr_track_MCV->Sumw2();
201 
202  Ecell_track_MCH=tfs->make<TH1F>("Ecell_track_MCH","",1000,0,0.1);
203  Ecell_track_MCH->Sumw2();
204  Ecell_track_MCV=tfs->make<TH1F>("Ecell_track_MCV","",1000,0,0.1);
205  Ecell_track_MCV->Sumw2();
206 
207  for( int i=0; i<214; ++i){
208  TString perplaneEnergy=Form("Ecells_plane%i",i);
209  planeEnergy[i]=tfs->make<TH1F>(perplaneEnergy,"",1000,0,0.1);
210  planeEnergy[i]->Sumw2();
211  }
212 
213  for( int ip=0; ip<214; ++ip ){
214  std::stringstream planenumber;
215  planenumber << ip;
216  std::string str_planenumber=planenumber.str();
217 
218  for( int ic=0; ic<3; ++ic ){
219  std::stringstream modulenumber;
220  modulenumber << ic;
221  std::string str_modulenumber=modulenumber.str();
222 
223  std::string namePECorr="modulePECorr_plane"+str_planenumber+"_module"+str_modulenumber;
224  TString str_namePECorr=namePECorr;
225  modulePECorr[ip][ic]=tfs->make<TH1F>(str_namePECorr,"",5000,0,5000);
226  modulePECorr[ip][ic]->Sumw2();
227 
228  std::string namePE="modulePE_plane"+str_planenumber+"_module"+str_modulenumber;
229  TString str_namePE=namePE;
230  modulePE[ip][ic]=tfs->make<TH1F>(str_namePE,"",5000,0,5000);
231  modulePE[ip][ic]->Sumw2();
232 
233  std::string nameADC="moduleADC_plane"+str_planenumber+"_module"+str_modulenumber;
234  TString str_nameADC=nameADC;
235  moduleADC[ip][ic]=tfs->make<TH1F>(str_nameADC,"",5000,0,5000);
236  moduleADC[ip][ic]->Sumw2();
237 
238  std::string nameEnergy="moduleEnergy_plane"+str_planenumber+"_module"+str_modulenumber;
239  TString str_nameEnergy=nameEnergy;
240  moduleEnergy[ip][ic]=tfs->make<TH1F>(str_nameEnergy,"",1000,0,0.1);
241  moduleEnergy[ip][ic]->Sumw2();
242  }//loop module index
243  }//loop plane index
244 
245  }
246 
247  //......................................................................
249  {
252 
253  //get all hits info
255  evt.getByLabel(fCellHitLabel, chits);
256 
257  //get reco slicer info
259  evt.getByLabel(fSliceLabel, slices);
260 
261  //get kalman track
262  art::FindManyP<rb::Track> fmKTrack(slices, evt, fKTrackLabel);
263 
264  //get reco vertex
265  art::FindManyP<rb::Vertex> fmVertex(slices, evt, fVertexLabel);
266 
267  for(unsigned int sliceIdx = 0; sliceIdx < slices->size(); ++sliceIdx){
268 
269  const rb::Cluster& slice = (*slices)[sliceIdx];
270  if(slice.IsNoise()) continue; //veto noise slice
271  double nCells = slice.NCell();
272  if( nCells<fMinHit ) continue; //minimum cells associated to slice
273 
274  //select reco vertex
275  if (!(fmVertex.isValid())) continue;
276  std::vector<art::Ptr<rb::Vertex>> vert = fmVertex.at(sliceIdx);
277  if (vert.size() != 1) continue;
278  if( fabs(vert[0]->GetX())>fMaxVertexX ) continue;
279  if( fabs(vert[0]->GetY())>fMaxVertexY ) continue;
280 
281  //select the kalman track
282  std::vector<art::Ptr<rb::Track>> Ktracks = fmKTrack.at(sliceIdx);
283  if( !(fmKTrack.isValid()) ) continue;
284 
285  //require there is only one track in the slice
286  if( Ktracks.size() != 1 ) continue;
287 
288  double Ltrack = Ktracks[0]->TotalLength();
289  if( Ltrack<fMinTrackLength ) continue; //track go through the full active detector
290 
291  TVector3 dir = Ktracks[0]->Dir();
292  double track_theta=dir.Theta();
293  if( cos(track_theta)<fMinTrackCosTheta ) continue;//track point forward
294 
295  TVector3 start = Ktracks[0]->Start();
296  double track_startz = start.Z();
297  if( track_startz>fMaxTrackStartZ ) continue;//track starts within first 4 planes
298 
299  vertexX->Fill(vert[0]->GetX());
300  vertexY->Fill(vert[0]->GetY());
301  vertexZ->Fill(vert[0]->GetZ());
302  trackLength->Fill(Ltrack);
303  trackStartX->Fill(start.X());
304  trackStartY->Fill(start.Y());
305  trackStartZ->Fill(start.Z());
306 
307  //save track associated cell information to vectors
308  std::vector<float> trkcellE;
309  std::vector<int> trkcellPlane;
310  std::vector<int> trkcellCell;
311  std::vector<int> trkcellADC;
312  std::vector<float> trkcellPECorr;
313  std::vector<float> trkcellPE;
314 
315  for( unsigned int icell=0; icell<Ktracks[0]->NCell(); ++ icell){
316  const art::Ptr<rb::CellHit>& chit = Ktracks[0]->Cell(icell);
317  const rb::RecoHit rhit(cal->MakeRecoHit(*chit, Ktracks[0]->W(chit.get())));
318 
319  if( !rhit.IsCalibrated() ) continue;
320 
321  trkcellE.push_back(rhit.GeV());
322  trkcellADC.push_back(chit->ADC());
323  trkcellPlane.push_back(chit->Plane());
324  trkcellCell.push_back(chit->Cell());
325  trkcellPECorr.push_back(rhit.PECorr());
326  trkcellPE.push_back(chit->PE());
327  }
328 
329  int ncells_track[6]={0};
330  double ecells_track[6]={0};
331  for( unsigned int ic=0; ic<trkcellE.size(); ++ic ){
332  ncells_track[0]++;
333  ecells_track[0] += trkcellE[ic];
334 
335  planeEnergy[trkcellPlane[ic]]->Fill(trkcellE[ic]);
336  Nevts_percell[trkcellPlane[ic]][trkcellCell[ic]]++;
337 
338  if(trkcellCell[ic]<32 ){
339  modulePE[trkcellPlane[ic]][0]->Fill(trkcellPE[ic]);
340  modulePECorr[trkcellPlane[ic]][0]->Fill(trkcellPECorr[ic]);
341  moduleADC[trkcellPlane[ic]][0]->Fill(trkcellADC[ic]);
342  moduleEnergy[trkcellPlane[ic]][0]->Fill(trkcellE[ic]);
343  }//H bottom or V east
344  else if( trkcellCell[ic]<64 ){
345  modulePE[trkcellPlane[ic]][1]->Fill(trkcellPE[ic]);
346  modulePECorr[trkcellPlane[ic]][1]->Fill(trkcellPECorr[ic]);
347  moduleADC[trkcellPlane[ic]][1]->Fill(trkcellADC[ic]);
348  moduleEnergy[trkcellPlane[ic]][1]->Fill(trkcellE[ic]);
349  }//middle
350  else{
351  modulePE[trkcellPlane[ic]][2]->Fill(trkcellPE[ic]);
352  modulePECorr[trkcellPlane[ic]][2]->Fill(trkcellPECorr[ic]);
353  moduleADC[trkcellPlane[ic]][2]->Fill(trkcellADC[ic]);
354  moduleEnergy[trkcellPlane[ic]][2]->Fill(trkcellE[ic]);
355  }//H top or V west
356 
357  PECorr_track[0]->Fill(trkcellPECorr[ic]);
358  Ecell_track[0]->Fill(trkcellE[ic]);
359 
360  if( trkcellPlane[ic]<192 ){
361  ncells_track[1]++;
362  ecells_track[1] += trkcellE[ic];
363  PECorr_track[1]->Fill(trkcellPECorr[ic]);
364  Ecell_track[1]->Fill(trkcellE[ic]);
365 
366  if( trkcellPlane[ic]<64 ){
367  ncells_track[2]++;
368  ecells_track[2] += trkcellE[ic];
369  PECorr_track[2]->Fill(trkcellPECorr[ic]);
370  Ecell_track[2]->Fill(trkcellE[ic]);
371  }//DB1
372  else if( trkcellPlane[ic]<128 ){
373  ncells_track[3]++;
374  ecells_track[3] += trkcellE[ic];
375  PECorr_track[3]->Fill(trkcellPECorr[ic]);
376  Ecell_track[3]->Fill(trkcellE[ic]);
377  }//DB2
378  else {
379  ncells_track[4]++;
380  ecells_track[4] += trkcellE[ic];
381  PECorr_track[4]->Fill(trkcellPECorr[ic]);
382  Ecell_track[4]->Fill(trkcellE[ic]);
383  }//DB3
384  }//DB1,2,3
385  else{
386  ncells_track[5]++;
387  ecells_track[5] += trkcellE[ic];
388  PECorr_track[5]->Fill(trkcellPECorr[ic]);
389  Ecell_track[5]->Fill(trkcellE[ic]);
390 
391  if( trkcellPlane[ic]%2==0 ){
392  PECorr_track_MCH->Fill(trkcellPECorr[ic]);
393  Ecell_track_MCH->Fill(trkcellE[ic]);
394  }
395  else{
396  PECorr_track_MCV->Fill(trkcellPECorr[ic]);
397  Ecell_track_MCV->Fill(trkcellE[ic]);
398  }
399 
400  }//muon catcher
401  }//end loop of track associated cells
402 
403  for( int id=0; id<6; ++id ){
404  Ncells_track[id]->Fill(ncells_track[id]);
405  Ecells_track[id]->Fill(ecells_track[id]);
406  }
407 
408  }//end loop of slice
409 
410  return;
411  }//end analyze
412 
414  {
416  TH1F *Nevts_incell_perplane[214];//# of events in each cell per each plane
417  for( int i=0; i<214; ++i ){
418  TString perplaneEvts=Form("Nevts_plane%i",i);
419  Nevts_incell_perplane[i]=tfs->make<TH1F>(perplaneEvts,"",96,-0.5,95.5);
420  }
421 
422  for( int ip=0; ip<214; ++ip ){
423  for( int ic=0; ic<96; ++ic ){
424  Nevts_incell_perplane[ip]->SetBinContent(ic+1,Nevts_percell[ip][ic]);
425  }//loop cell number
426  }//loop plane number
427 
428  }
429 
430 }//end namespace
431 
432 namespace rockmuon{
433 
435 
436 }
::xsd::cxx::tree::id< char, ncname > id
Definition: Database.h:165
std::string fKTrackLabel
label for slices
TH1F * moduleEnergy[214][3]
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
#define SET(T, N)
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
std::string fVertexLabel
label for cell hits
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
unsigned short Plane() const
Definition: CellHit.h:39
const char * p
Definition: xmltok.h:285
void reconfigure(const fhicl::ParameterSet &p)
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
A collection of associated CellHits.
Definition: Cluster.h:47
TString ip
Definition: loadincs.C:5
DEFINE_ART_MODULE(TestTMapFile)
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
TH1F * modulePECorr[214][3]
void beginJob()
int evt
float PE() const
Definition: CellHit.h:42
double GetY(int dcm=1, int feb=0, int pix=0)
Definition: PlotOnMon.C:119
std::string fCellHitLabel
label for kalman tracks
Vertex location in position and time.
Perform a "2 point" Hough transform on a collection of hits.
::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
T const * get() const
Definition: Ptr.h:321
TDirectory * dir
Definition: macro.C:5
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
T cos(T number)
Definition: d0nt_math.hpp:78
int fMinHit
label for vertex
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
void analyze(const art::Event &evt)
Encapsulate the geometry of one entire detector (near, far, ndos)
double GetX(int ndb=14, int db=1, int feb=0, int pix=0)
Definition: PlotOnMon.C:111