StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
qa.C
1 //usr/bin/env root4star -l -b -q $0'("'"$1"'")'; exit $?
2 
3 #include "TTree.h"
4 #include "TClonesArray.h"
5 #include <map>
6 // #include <__config>
7 
8 // TClonesArrays to read from the tree
9 TClonesArray *mcTracks = NULL;
10 TClonesArray *fttPoints = NULL;
11 TClonesArray *fttClusters = NULL;
12 TClonesArray *fstPoints = NULL;
13 TClonesArray *wcal = NULL;
14 TClonesArray *wcalHits = NULL;
15 TClonesArray *hcal = NULL;
16 TClonesArray *hcalHits = NULL;
17 TClonesArray *epdHits = NULL;
18 TClonesArray *fwdTracks = NULL;
19 TClonesArray *seeds = NULL;
20 TTree *t = NULL;
21 
22 std::map<std::string, TH1*> histograms;
23 TH1* addH1( string name, string title, int nx, float x1, float x2 ){
24  histograms[name] = new TH1F( name.c_str(), title.c_str(), nx, x1, x2 );
25  return histograms[name];
26 }
27 TH2* addH2( string name, string title, int nx, float x1, float x2, int ny, float y1, float y2 ){
28  histograms[name] = new TH2F( name.c_str(), title.c_str(), nx, x1, x2, ny, y1, y2 );
29  return (TH2*)histograms[name];
30 }
31 inline TH1 *getH1( string name ){
32  // printf( "Looking for histogram name=[%s]", name.c_str() );
33  assert( histograms.count( name ) && "Histogram cannot be found" && name.c_str() );
34  assert( histograms[name] && TString::Format( "Histogram %s is NULL", name.c_str() ) );
35  return histograms[name];
36 }
37 inline TH2 *getH2( string name ){
38  return (TH2*) getH1( name );
39 }
40 
41 void setupRead( TString filename = "fwdtree.root" ){
42  // setup and make sure libraries are loaded
43  gSystem->Load( "libStarRoot.so" );
44  gSystem->Load("libStarClassLibrary.so");
45  gROOT->SetMacroPath(".:/star-sw/StRoot/macros/:./StRoot/macros:./StRoot/macros/graphics:./StRoot/macros/analysis:./StRoot/macros/test:./StRoot/macros/examples:./StRoot/macros/html:./StRoot/macros/qa:./StRoot/macros/calib:./StRoot/macros/mudst:/afs/rhic.bnl.gov/star/packages/DEV/StRoot/macros:/afs/rhic.bnl.gov/star/packages/DEV/StRoot/macros/graphics:/afs/rhic.bnl.gov/star/packages/DEV/StRoot/macros/analysis:/afs/rhic.bnl.gov/star/packages/DEV/StRoot/macros/test:/afs/rhic.bnl.gov/star/packages/DEV/StRoot/macros/examples:/afs/rhic.bnl.gov/star/packages/DEV/StRoot/macros/html:/afs/rhic.bnl.gov/star/packages/DEV/StRoot/macros/qa:/afs/rhic.bnl.gov/star/packages/DEV/StRoot/macros/calib:/afs/rhic.bnl.gov/star/packages/DEV/StRoot/macros/mudst:/afs/rhic.bnl.gov/star/ROOT/36/5.34.38/.sl73_x8664_gcc485/rootdeb/macros:/afs/rhic.bnl.gov/star/ROOT/36/5.34.38/.sl73_x8664_gcc485/rootdeb/tutorials");
46  gROOT->LoadMacro("$STAR/StRoot/StMuDSTMaker/COMMON/macros/loadSharedLibraries.C");
47  loadSharedLibraries();
48 
49  gSystem->Load("libgenfit2.so");
50  gSystem->Load("libKiTrack.so");
51  gSystem->Load( "libStFwdTrackMaker.so" );
52 
53  // now open our data file
54  TFile *f = new TFile(filename);
55  t = (TTree*)f->Get("fwd");
56 
57  // create the readers for each branch
58  mcTracks = new TClonesArray("StMuMcTrack");
59  t->GetBranch("mcTracks")->SetAutoDelete(kFALSE);
60  t->SetBranchAddress("mcTracks",&mcTracks);
61 
62  fttPoints = new TClonesArray("StMuFttPoint");
63  t->GetBranch("fttPoints")->SetAutoDelete(kFALSE);
64  t->SetBranchAddress("fttPoints",&fttPoints);
65 
66  fttClusters = new TClonesArray("StMuFttCluster");
67  t->GetBranch("fttClusters")->SetAutoDelete(kFALSE);
68  t->SetBranchAddress("fttClusters",&fttClusters);
69 
70  fstPoints = new TClonesArray("StMuFstHit");
71  t->GetBranch("fstHits")->SetAutoDelete(kFALSE);
72  t->SetBranchAddress("fstHits",&fstPoints);
73 
74  wcal = new TClonesArray("FcsClusterWithStarXYZ");
75  t->GetBranch("wcalClusters")->SetAutoDelete(kFALSE);
76  t->SetBranchAddress("wcalClusters",&wcal);
77 
78  wcalHits = new TClonesArray("FcsHitWithStarXYZ");
79  t->GetBranch("wcalHits")->SetAutoDelete(kFALSE);
80  t->SetBranchAddress("wcalHits",&wcalHits);
81 
82  hcal = new TClonesArray("FcsClusterWithStarXYZ");
83  t->GetBranch("hcalClusters")->SetAutoDelete(kFALSE);
84  t->SetBranchAddress("hcalClusters",&hcal);
85 
86  hcalHits = new TClonesArray("FcsHitWithStarXYZ");
87  t->GetBranch("hcalHits")->SetAutoDelete(kFALSE);
88  t->SetBranchAddress("hcalHits",&hcalHits);
89 
90  epdHits = new TClonesArray("FcsHitWithStarXYZ");
91  t->GetBranch("epdHits")->SetAutoDelete(kFALSE);
92  t->SetBranchAddress("epdHits",&epdHits);
93 
94  fwdTracks = new TClonesArray("StMuFwdTrack");
95  t->GetBranch("reco")->SetAutoDelete(kFALSE);
96  t->SetBranchAddress("reco",&fwdTracks);
97 
98  seeds = new TClonesArray("StMuFwdTrackSeedPoint");
99  t->GetBranch("seeds")->SetAutoDelete(kFALSE);
100  t->SetBranchAddress("seeds",&seeds);
101 }
102 
103 void qaMomentumResolution(){
104 
105  if ( histograms.count( "curveRcVsMc" ) == 0 ) {
106  printf( "Creating Momentum Resolution Histograms\n" );
107  addH2( "curveRcVsMc", "Track curvature; MC; RC", 200, -10, 10, 200, -10, 10 );
108  addH2( "curveMcVsPtMc", ";MC Pt; MC Curve ", 200, 0, 10, 100, 0, 10 );
109  addH2( "ptRcVsMc", "Track Pt; MC; RC", 200, 0, 10, 200, 0, 10 );
110  addH1( "curveRes", "Curvature Resolution; (C^{MC}-C^{RC})/C^{MC}", 200, -2, 2 );
111  addH1( "transMomRes", "Pt Resolution; (Pt^{MC} - Pt^{RC}) / Pt^{MC}", 200, -2, 2 );
112  addH1( "deltaCharge", "deltaCharge; |q_{MC}-q_{RC}|;counts;", 5, 0, 5 );
113  }
114 
115  const TH2 * hCurveRcVsMc = getH2( "curveRcVsMc" );
116  const TH2 * hCurveMcVsPtMc = getH2( "curveMcVsPtMc" );
117  const TH2 * hPtRcVsMc = getH2( "ptRcVsMc" );
118 
119  const TH1 * hCurveRes = getH1( "curveRes" );
120  const TH1 * hTransMomRes = getH1( "transMomRes" );
121  const TH1 * hDeltaCharge = getH1( "deltaCharge" );
122 
123  for ( int j = 0; j < fwdTracks->GetEntries(); j++ ){
124 
125  StMuFwdTrack *fwt = fwdTracks->At(j);
126  UShort_t indexToMatchedMC = fwt->idTruth() - 1;
127  // // cout << "Processing track " << j << ", mcid = " << indexToMatchedMC << endl;
128  if (indexToMatchedMC >= mcTracks->GetEntries()) continue;
129  // // get corresponding MC track
130  StMuMcTrack *mct = mcTracks->At(indexToMatchedMC);
131 
132  float curveMc = fabs(mct->Charge() / mct->pT());
133  float curveRc = fabs(fwt->charge() / fwt->momentum().Pt());
134  // // cout << "mct->pT() = " << mct->pT() << endl;
135  if ( mct->pT() > 0.1 && fwt->pval() > 0.01){
136  hDeltaCharge->Fill( abs( mct->Charge() - fwt->charge() ) );
137  hCurveRes->Fill( (curveMc - curveRc) / curveMc );
138  hTransMomRes->Fill( (mct->pT() - fwt->momentum().Pt()) / mct->pT() );
139 
140  hCurveRcVsMc->Fill( curveMc, curveRc );
141  hCurveMcVsPtMc->Fill( curveMc, mct->pT() );
142  hPtRcVsMc->Fill( mct->pT(), fwt->momentum().Pt() );
143  }
144  }
145 }
146 
147 void qaFCSTrackMatch(){
148 
149  if (histograms.count("dxWCAL") == 0){
150  addH1( "dxWCAL", "WCAL; dx", 500, -100, 100 );
151  addH1( "dxWCAL2", "WCAL; dx", 500, -100, 100 );
152  addH1( "dyWCAL", "WCAL; dy", 500, -100, 100 );
153  addH1( "dxHCAL", "HCAL; dx", 500, -100, 100 );
154  addH1( "drWCAL", "WCAL; dr", 500, 0, 100 );
155  addH1( "drHCAL", "HCAL; dr", 500, 0, 100 );
156  addH2( "wcalTrackX", "; WCAL x; Track x", 500, -50, 50, 500, -50, 50 );
157  addH2( "wcalTrackY", "; WCAL y; Track y", 500, -50, 50, 500, -50, 50 );
158  }
159 
160  const TH1 * hdxWCAL = getH1("dxWCAL");
161  const TH1 * hdxWCAL2 = getH1("dxWCAL2");
162  const TH1 * hdyWCAL = getH1("dyWCAL");
163  const TH1 * hdxHCAL = getH1("dxHCAL");
164 
165  const TH1 * hdrWCAL = getH1("drWCAL");
166  const TH1 * hdrHCAL = getH1("drHCAL");
167 
168  const TH2 * hWCALX = getH1("wcalTrackX");
169  const TH2 * hWCALY = getH1("wcalTrackY");
170 
171  for ( int j = 0; j < fwdTracks->GetEntries(); j++ ){
172  StMuFwdTrack *track = (StMuFwdTrack *)fwdTracks->At(j);
173 
174  // printf("Track %d: pt=%f, eta=%f, phi=%f\n", j, track->momentum().Pt(), track->momentum().Eta(), track->momentum().Phi());
175  if ( track->pval() < 0.01 ) continue;
176  StMuFwdTrackProjection projWCAL;
177  track->getProjectionFor(kFcsWcalId, projWCAL);
178  // printf("Projection @ WCAL: det=%d, x=%f, y=%f, z=%f\n", projWCAL.mDetId, projWCAL.mXYZ.X(), projWCAL.mXYZ.Y(), projWCAL.mXYZ.Z());
179 
180  StMuFwdTrackProjection projHCAL;
181  track->getProjectionFor(kFcsHcalId, projHCAL);
182  // printf("Projection @ HCAL: det=%d, x=%f, y=%f, z=%f\n", projHCAL.mDetId, projHCAL.mXYZ.X(), projHCAL.mXYZ.Y(), projHCAL.mXYZ.Z());
183 
184  // loop over WCAL clusters
185  for ( int k = 0; k < wcal->GetEntries(); k++ ){
186  FcsClusterWithStarXYZ *cluster = (FcsClusterWithStarXYZ*)wcal->At(k);
187  double dx = projWCAL.mXYZ.X() - cluster->mXYZ.X();
188  double dy = projWCAL.mXYZ.Y() - cluster->mXYZ.Y();
189  double dr = sqrt( dx*dx + dy*dy );
190 
191  hdxWCAL2->Fill( dx );
192  if ( projWCAL.mXYZ.X() < 0 && cluster->mClu->detectorId() == 1 ) continue;
193  if ( projWCAL.mXYZ.X() > 0 && cluster->mClu->detectorId() == 0 ) continue;
194 
195  hdxWCAL->Fill( dx );
196  hdyWCAL->Fill( dy );
197  // printf("WCAL Cluster %d: x=%f, y=%f, z=%f\n", k, cluster->mXYZ.X(), cluster->mXYZ.Y(), cluster->mXYZ.Z());
198  hdrWCAL->Fill( dr );
199 
200  hWCALX->Fill( cluster->mXYZ.X(), projWCAL.mXYZ.X() );
201  hWCALY->Fill( cluster->mXYZ.Y(), projWCAL.mXYZ.Y() );
202  }
203 
204  // loop over WCAL clusters
205  for ( int k = 0; k < hcal->GetEntries(); k++ ){
206  FcsClusterWithStarXYZ *cluster = (FcsClusterWithStarXYZ*)hcal->At(k);
207  double dx = projHCAL.mXYZ.X() - cluster->mXYZ.X();
208  double dy = projHCAL.mXYZ.Y() - cluster->mXYZ.Y();
209  double dr = sqrt( dx*dx + dy*dy );
210 
211  hdxHCAL->Fill( dx );
212  hdrHCAL->Fill( dr );
213  }
214 
215  }
216 }
217 
218 
219 // Loop on events
220 void eventLoop( int numEventsLimit = -1, int reportEveryNthEvent = -1 ){
221  int lastEventIndex = (numEventsLimit > 0 ? numEventsLimit : t->GetEntries() );
222  for ( int i = 0; i < lastEventIndex; i++ ){
223  t->GetEntry(i);
224  if ( reportEveryNthEvent > 0 && i % reportEveryNthEvent == 0){
225  printf( "Processing Event %d...\n", i );
226  }
227  // run qa subroutines here
228  // qaMomentumResolution();
229  qaFCSTrackMatch();
230  }
231 }
232 
233 
234 
235 void qa( TString filename = "fwdtree.root" ){
236  setupRead( filename );
237  TFile * fOut = new TFile( "QuickQA.root", "RECREATE" );
238  fOut->cd();
239  eventLoop(5000, 100);
240  fOut->Write();
241  // writeHistograms();
242  return;
243  //loop over the events
244  int nEvents = t->GetEntries();
245  if (nEvents > 1000) nEvents = 1000;
246  for ( int i = 0; i < nEvents; i++ ){
247  t->GetEntry(i);
248  for ( int j = 0; j < fwdTracks->GetEntries(); j++ ){
249  StMuFwdTrack *track = fwdTracks->At(j);
250  printf("Track %d: pt=%f, eta=%f, phi=%f\n", j, track->momentum().Pt(), track->momentum().Eta(), track->momentum().Phi());
251 
252  StMuFwdTrackProjection projWCAL;
253  track->getProjectionFor(kFcsWcalId, projWCAL);
254  printf("Projection @ WCAL: det=%d, x=%f, y=%f, z=%f\n", projWCAL.mDetId, projWCAL.mXYZ.X(), projWCAL.mXYZ.Y(), projWCAL.mXYZ.Z());
255 
256 
257  StMuFwdTrackProjection projHCAL;
258  track->getProjectionFor(kFcsHcalId, projHCAL);
259  printf("Projection @ HCAL: det=%d, x=%f, y=%f, z=%f\n", projHCAL.mDetId, projHCAL.mXYZ.X(), projHCAL.mXYZ.Y(), projHCAL.mXYZ.Z());
260 
261  // loop over WCAL clusters
262  for ( int k = 0; k < wcal->GetEntries(); k++ ){
263  FcsClusterWithStarXYZ *cluster = (FcsClusterWithStarXYZ*)wcal->At(k);
264 
265  printf("WCAL Cluster %d: x=%f, y=%f, z=%f\n", k, cluster->mXYZ.X(), cluster->mXYZ.Y(), cluster->mXYZ.Z());
266 
267  }
268 
269 
270  // loop over WCAL clusters
271  for ( int k = 0; k < wcal->GetEntries(); k++ ){
272  FcsClusterWithStarXYZ *cluster = (FcsClusterWithStarXYZ*)wcal->At(k);
273 
274  printf("WCAL Cluster %d: x=%f, y=%f, z=%f\n", k, cluster->mXYZ.X(), cluster->mXYZ.Y(), cluster->mXYZ.Z());
275 
276  }
277 
278  }
279  }
280  cout << "Processed: " << t->GetEntries() << " entries" << endl;
281 }
FcsClusterWithStarXYZ
Store Cluster with STAR XYZ position.
Definition: StFwdQAMaker.h:119
Generated by   doxygen 1.8.5