Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MvtxQAHisto.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file MvtxQAHisto.C
1 
2 //
3 // This module is desgined to grab svtx tracks and put truth and cluster
4 // information into a TTree for GenFit testing
5 //
7 //
8 // Darren McGlinchey
9 // 1 Apr 2016
10 //
12 
13 
14 #include "MvtxQAHisto.h"
15 
16 #include <phool/phool.h>
17 #include <phool/getClass.h>
18 #include <g4main/PHG4HitContainer.h>
19 #include <g4main/PHG4TruthInfoContainer.h>
20 #include <g4main/PHG4Particle.h>
21 #include <g4main/PHG4Hit.h>
22 #include <g4main/PHG4VtxPoint.h>
23 #include <fun4all/PHTFileServer.h>
24 #include <fun4all/Fun4AllServer.h>
25 #include <phool/recoConsts.h>
26 
27 #include <trackbase/TrkrDefUtil.h>
28 #include <trackbase/TrkrClusterContainer.h>
29 #include <trackbase/TrkrCluster.h>
30 #include <trackbase/TrkrClusterv1.h>
31 #include <trackbase/TrkrHitSetContainer.h>
32 #include <trackbase/TrkrHitSet.h>
33 #include <trackbase/TrkrHitSetv1.h>
34 #include <mvtx/MvtxDefUtil.h>
35 #include <mvtx/MvtxHitSetv1.h>
36 
37 #include <TTree.h>
38 #include <TH1F.h>
39 #include <TH2F.h>
40 #include <TVector3.h>
41 
42 #include <iostream>
43 
44 using namespace std;
45 
46 //----------------------------------------------------------------------------//
47 //-- Constructor:
48 //-- simple initialization
49 //----------------------------------------------------------------------------//
50 MvtxQAHisto::MvtxQAHisto(const string &name):
51  SubsysReco( name )
52 {
53  //initialize
54  _outfile = "MvtxQAHisto.root";
55 
56  beam_x[0] = 826.693;
57  beam_x[1] = 822.267;
58  beam_x[2] = 818.413;
59  beam_x[3] = 830.190;
60 
61  beam_y[0] = 158.773;
62  beam_y[1] = 167.537;
63  beam_y[2] = 181.318;
64  beam_y[3] = 190.988;
65 
66 }
67 
68 //----------------------------------------------------------------------------//
69 //-- Init():
70 //-- Intialize all histograms, trees, and ntuples
71 //----------------------------------------------------------------------------//
72 int MvtxQAHisto::Init(PHCompositeNode *topNode)
73 {
74 
75  //recoConsts *rc = recoConsts::instance();
76  cout << PHWHERE << " Openning file " << _outfile << endl;
77  PHTFileServer::get().open( _outfile, "RECREATE");
78 
79  for (int ichip=0; ichip<4; ichip++){
80  h2d_hit[ichip] = new TH2F(Form("h2d_hit_chip%d",ichip),"",1024,-0.5,1023.5,512,-0.5,511.5);
81  h1d_hit_per_evt[ichip] = new TH1F(Form("h1d_hit_per_evt_chip%d",ichip),"",101,-0.5,100);
82  h1d_clus_per_evt[ichip] = new TH1F(Form("h1d_clus_per_evt_chip%d",ichip),"",51,-0.5,50.5);
83 
84  h2d_hit_beam[ichip] = new TH2F(Form("h2d_hit_beam_chip%d",ichip),"",1025,-512.5,512.5,513,-256.5,256.5);
85  h2d_hit_trk[ichip] = new TH2F(Form("h2d_hit_trk_chip%d",ichip),"",1025,-512.5,512.5,513,-256.5,256.5);
86  h2d_clus[ichip] = new TH2F(Form("h2d_clus_chip%d",ichip),"",1024,0,3.0,512,0,1.5);
87  h2d_clus_beam[ichip] = new TH2F(Form("h2d_clus_beam_chip%d",ichip),"",1024,-3.0/2,3.0/2,512,-1.5/2,1.5/2);
88 
89  h1d_clus_size_x[ichip] = new TH1F(Form("h1d_clus_size_x_chip%d",ichip),"",51,-0.5,50.5);
90  h1d_clus_size_z[ichip] = new TH1F(Form("h1d_clus_size_z_chip%d",ichip),"",51,-0.5,50.5);
91 
92  }
93  h1d_trk_finder_x = new TH1F("h1d_trk_finder_x","",513,-256.5,256.5);
94  h1d_trk_finder_z = new TH1F("h1d_trk_finder_z","",1025,-512.5,512.5);
95  h2d_trk_finder = new TH2F("h2d_trk_finder","",1025,-512.5,512.5,513,-256.5,256.5);
96 
97  h1d_clus_associated = new TH1F("h1d_clus_associated","",11,-0.5,10.5);
98  h1d_clus_eff = new TH1F("h1d_clus_eff","",5,-0.5,4.5);
99 
100  return 0;
101 }
102 
103 //----------------------------------------------------------------------------//
104 //-- process_event():
105 //-- Call user instructions for every event.
106 //-- This function contains the analysis structure.
107 //----------------------------------------------------------------------------//
108 int MvtxQAHisto::process_event(PHCompositeNode *topNode)
109 {
110  _event++;
111  //if (1)
112  if (_event % 1000 == 0)
113  cout << PHWHERE << "Events processed: " << _event << endl;
114 
115  GetNodes(topNode);
116 
117  TrkrDefUtil trkrdefutil;
118  MvtxDefUtil mvtxdefutil;
119 
120  int nhit_per_chip[4] = {0};
121  int nclus_per_chip[4] = {0};
122 
123  double avg_x = 0.0;
124  double avg_y = 0.0;
125  double norm = 0.0;
126 
127  h1d_trk_finder_x->Reset();
128  h1d_trk_finder_z->Reset();
129  h2d_trk_finder->Reset();
130 
131  TrkrHitSetContainer::ConstRange iter_range = hitsetcon->GetHitSets();
132  for ( TrkrHitSetContainer::ConstIterator iter = iter_range.first; iter!=iter_range.second; ++iter){
133 
134  int ichip = int(mvtxdefutil.GetChipId(iter->first));
135 
136  MvtxHitSetv1 *hitset = static_cast<MvtxHitSetv1 *>(iter->second);
137 
138  MvtxHitSetv1::ConstRange hit_iter_range = hitset->GetHits();
139  for ( MvtxHitSetv1::ConstIterator hit_iter = hit_iter_range.first; hit_iter!=hit_iter_range.second; ++hit_iter){
140 
141  int icol = int(hit_iter->first);
142  int irow = int(hit_iter->second);
143 
144  h2d_hit[ichip]->Fill(icol,irow);
145 
146  h1d_trk_finder_z->Fill(icol-beam_x[ichip]);
147  h1d_trk_finder_x->Fill(irow-beam_y[ichip]);
148  h2d_trk_finder->Fill(icol-beam_x[ichip],irow-beam_y[ichip]);
149 
150  nhit_per_chip[ichip]++;
151 
152  }//hit_iter
153  }//iter
154 
155  TrkrClusterContainer::ConstRange clus_iter_range = cluscon->GetClusters();
156  for (TrkrClusterContainer::ConstIterator clus_iter = clus_iter_range.first; clus_iter!=clus_iter_range.second; ++clus_iter){
157 
158  int ichip = trkrdefutil.GetLayer(clus_iter->first);
159 
160  TrkrCluster *clus = clus_iter->second;
161 
162  h2d_clus[ichip]->Fill(clus->GetZ(),clus->GetX());
163  h2d_clus_beam[ichip]->Fill(clus->GetZ() - beam_x[ichip]*(28e-4),clus->GetX() - beam_y[ichip]*(28e-4));
164 
165  h1d_clus_size_z[ichip]->Fill(clus->GetZSize());
166  h1d_clus_size_x[ichip]->Fill(clus->GetPhiSize());
167 
168  nclus_per_chip[ichip]++;
169  }//
170 
171  int nchip_w_clus = 0;
172 
173  for (int ichip=0; ichip<4; ichip++){
174  h1d_hit_per_evt[ichip]->Fill(nhit_per_chip[ichip]);
175  h1d_clus_per_evt[ichip]->Fill(nclus_per_chip[ichip]);
176 
177  if ( nclus_per_chip[ichip]>0 ) nchip_w_clus++;
178  }
179 
180  //if ( nchip_w_clus<3 ) return 0;
181  //if ( 1 ) return 0;
182 
183  //fast track reconstruction
184  /*
185  double trk_z = h1d_trk_finder_z->GetBinCenter(h1d_trk_finder_z->GetMaximumBin());
186  double trk_x = h1d_trk_finder_x->GetBinCenter(h1d_trk_finder_x->GetMaximumBin());
187  double trk_z_w = h1d_trk_finder_z->GetBinContent(h1d_trk_finder_z->GetMaximumBin());
188  double trk_x_w = h1d_trk_finder_x->GetBinContent(h1d_trk_finder_x->GetMaximumBin());
189  */
190 
191  /*
192  double trk_z = 0;
193  double trk_x = 0;
194  double trk_w = 0.1;
195  for (int ix=0; ix<h2d_trk_finder->GetNbinsX(); ix++){
196  for (int iy=0; iy<h2d_trk_finder->GetNbinsY(); iy++){
197  if ( h2d_trk_finder->GetBinContent(ix+1,iy+1)<trk_w ) continue;
198 
199  trk_z = h1d_trk_finder_z->GetBinCenter(ix+1);
200  trk_x = h1d_trk_finder_x->GetBinCenter(iy+1);
201  trk_w = h2d_trk_finder->GetBinContent(ix+1,iy+1);
202 
203  }
204  }
205 
206  if ( trk_w<3 ) return 0;
207 
208  int nclus_associated = 0;
209  int nclus_associated_per_chip[4] = {0};
210 
211  //TrkrClusterContainer::ConstRange clus_iter_range = cluscon->GetClusters();
212  for (TrkrClusterContainer::ConstIterator clus_iter = clus_iter_range.first; clus_iter!=clus_iter_range.second; ++clus_iter){
213 
214  int ichip = trkrdefutil.GetLayer(clus_iter->first);
215 
216  TrkrCluster *clus = clus_iter->second;
217 
218  double clusz = clus->GetZ()/(28e-4) - beam_x[ichip];
219  double clusx = clus->GetX()/(28e-4) - beam_y[ichip];
220 
221  if ( fabs(clusz-trk_z)<5 && fabs(clusx-trk_x)<5 ){
222  nclus_associated_per_chip[ichip]++;
223  nclus_associated++;
224  }
225 
226  }//
227 
228  h1d_clus_associated->Fill(nclus_associated);
229 
230  if ( nclus_associated<3 ) return 0;
231 
232  for (int ichip=0; ichip<4; ichip++){
233  if ( nclus_associated_per_chip[ichip]>0 ) h1d_clus_eff->Fill(ichip);
234  }
235 
236  //cout << "N Associated Clus: " << nclus_associated << endl;
237 
238 
239  //cout << "Z: " << trk_z << ", X: " << trk_x << ", S: " << trk_w << endl;
240 
241  //TrkrHitSetContainer::ConstRange iter_range = hitsetcon->GetHitSets();
242  for ( TrkrHitSetContainer::ConstIterator iter = iter_range.first; iter!=iter_range.second; ++iter){
243 
244  int ichip = int(mvtxdefutil.GetChipId(iter->first));
245 
246  MvtxHitSetv1 *hitset = static_cast<MvtxHitSetv1 *>(iter->second);
247 
248  MvtxHitSetv1::ConstRange hit_iter_range = hitset->GetHits();
249  for ( MvtxHitSetv1::ConstIterator hit_iter = hit_iter_range.first; hit_iter!=hit_iter_range.second; ++hit_iter){
250 
251  int icol = int(hit_iter->first);
252  int irow = int(hit_iter->second);
253 
254  h2d_hit_trk[ichip]->Fill(icol-beam_x[ichip]-trk_z,irow-beam_y[ichip]-trk_x);
255 
256  }//hit_iter
257  }//iter
258  */
259 
260 
261  return 0;
262 }
263 
264 //----------------------------------------------------------------------------//
265 //-- End():
266 //-- End method, wrap everything up
267 //----------------------------------------------------------------------------//
268 int MvtxQAHisto::End(PHCompositeNode *topNode)
269 {
270 
271  cout << "-----MvtxQAHisto::End------" << endl;
272 
273  PHTFileServer::get().cd( _outfile );
274  PHTFileServer::get().write( _outfile );
275  //PHTFileServer::get().close();
276 
277  return 0;
278 }
279 
280 
281 //----------------------------------------------------------------------------//
282 //-- GetNodes():
283 //-- Get all the all the required nodes off the node tree
284 //----------------------------------------------------------------------------//
285 void MvtxQAHisto::GetNodes(PHCompositeNode * topNode)
286 {
287  //DST objects
288  //
289  //
290  cluscon = findNode::getClass<TrkrClusterContainer>(topNode, "TrkrClusterContainer");
291  if (cluscon==0)
292  {
293  cout << "MvtxQAHisto::Process_Event - TrkrClusterContainer not found" << endl;
294  exit(-1);
295  }
296 
297  hitsetcon = findNode::getClass<TrkrHitSetContainer>(topNode, "TrkrHitSetContainer");
298  if (hitsetcon==0)
299  {
300  cout << "MvtxQAHisto::Process_Event - TrkrHitSetContainer not found" << endl;
301  exit(-1);
302  }
303 
304 }
305 
306 
307