StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StTruthTestMaker.cxx
1 #include "StTruthTestMaker.h"
3 
4 #include "TH1Helper.h"
5 #include "StMuDSTMaker/COMMON/StMuDst.h"
6 #include "StMuDSTMaker/COMMON/StMuTrack.h"
7 #include "StMuDSTMaker/COMMON/StMuDebug.h"
8 
9 #include "TDataSetIter.h"
10 #include "tables/St_g2t_track_Table.h"
11 
12 #include "TString.h"
13 
14 #include "StarGenerator/EVENT/StarGenEvent.h"
15 #include "StarGenerator/EVENT/StarGenParticle.h"
16 
17 #include <map>
18 using namespace std;
19 
20 // Few histograms in global scope
21 #include "TH1.h"
22 #include "TH1F.h"
23 #include "TH2F.h"
24 
25 #include "TVector2.h"
26 
27 // Few cuts also
28 const Float_t invCut = 0.12000 * 3.0;
29 const Float_t etaCut = 0.01200 * 3.0;
30 const Float_t phiCut = 0.01900 * 3.0;
31 const Float_t absEtaCut = 1.0;
32 
33 #define Res(x) TMath::Abs( x##R - x##T )
34 
35 #define MUDST_CUT qaTruth < 95 || track -> nHits() < 10
36 
37 // ------------------------------------------------------------------------------------------ Init //
38 Int_t StTruthTestMaker::Init()
39 {
40  hMatchedEta = new TH2F("hMatchedEta",";#eta thrown;#eta reco",201,-1.005,1.005,201,-1.005,1.005);
41  hMatchedPhi = new TH2F("hMatchedPhi",";#phi thrown;#phi reco",201,-1.005*TMath::Pi(),1.005*TMath::Pi(),201,-1.005*TMath::Pi(),1.005*TMath::Pi());
42  hMatchedPt = new TH2F("hMatchedPt", ";p_{T} thrown;p_{T} reco",200,0.,10.,200,0.,10.);
43  hMatchedInv = new TH2F("hMatchedInv",";1/p_{T} thrown;1/p_{T} reco",200,0.,10.,200,0.,10.);
44  hMatchedPID = new TH1F("hMatchedPID","PID of matched particles",2112*2+1.0,-2112.5,+2112.5);
45  TH1Helper::SetCanRebin(hMatchedPID);
46  hMatchedQA = new TH1F("hMatchedQA", "QA truth of matched tracks",101,-0.5,100.5);
47 
48  hMatchedEtaRes = new TH1F("hMatchedEtaRes",";#eta resolution",201,-0.2*1.005,0.2*1.005);
49  hMatchedPhiRes = new TH1F("hMatchedPhiRes",";#phi resolution",201,-0.2*1.005,0.2*1.005);
50  hMatchedPtRes = new TH1F("hMatchedPtRes", ";p_{T} resolution",201,-1.005,1.005);
51  hMatchedInvRes = new TH1F("hMatchedInvRes",";1/p_{T} resolution",201,-1.005,1.005);
52 
53  hNumMismatched = new TH1F("hNumMismatched","Number of mismatched tracks / event;N mismatched",11,-0.5,10.5);
54  hPerMismatched = new TH1F("hPerMismatched","Percentage of mismatched globals / event;N mismatched/total",21,-0.025,1.025);
55  hPidMismatched = new TH1F("hPidMismatched","PID of mismatched globals",2112*2+1.0,-2112.5,+2112.5);
56  TH1Helper::SetCanRebin(hPidMismatched);
57 
58  TH1 *hList[]={hMatchedEta, hMatchedPhi, hMatchedPt, hMatchedPID,hMatchedEtaRes, hMatchedPhiRes, hMatchedPtRes, hNumMismatched,hPerMismatched, hPidMismatched, hMatchedInv, hMatchedInvRes, hMatchedQA};
59  // const Char_t *xtitle[]={"#eta thrown","#phi thrown","pt thrown","pid", "#eta resolution","#phi resolution","pt resolution","N mismatched","fraction mismatched","pid mismatched","1/pt thrown","1/pt resolution","QA truth"};
60  // const Char_t *ytitle[]={"#eta reco", "#phi reco", "pt reco", "", "", "", "", "", "", "", "1/pt reco", "", "" };
61 
62  for ( UInt_t i=0;i<sizeof(hList)/sizeof(TH1*);i++ ) {
63  // hList[i]->GetXaxis()->SetTitle(xtitle[i]);
64  // hList[i]->GetYaxis()->SetTitle(ytitle[i]);
65  AddHist( hList[i] );
66  // cout << "Added histogram " << hList[i]->GetName() << endl;
67  }
68 
69 
70 
71  TString name=GetName();
72  if ( name=="testGeant" ) mDoGeant = true;
73  else mDoGeant = false;
74 
75  return kStOK;
76 }
77 // ------------------------------------------------------------------------------------------ Init //
78 
79 Int_t StTruthTestMaker::Make()
80 {
81 
82  if ( mDoGeant ) MakeGeant();
83  else MakeRecord();
84  return kStOK;
85 
86 }
87 
88 Int_t StTruthTestMaker::MakeGeant()
89 {
90 
91  static Bool_t first = true;
92 
93  /*
94  *****************************************************************
95 
96  Obtain the MuDst and map all tracks to their ID truth value
97 
98  *****************************************************************
99  */
100 
101  map< Int_t, StMuTrack * > muTrackMap;
102 
103  StMuDst *muDst = (StMuDst *)GetInputDS("MuDst"); assert(muDst);
104  StMuEvent *muEvent = muDst -> event(); assert(muEvent);
105 
106  { Int_t ntracks = muDst -> globalTracks() -> GetEntries();
107  for ( Int_t itrack = 0; itrack < ntracks; itrack++ )
108  {
109 
110  StMuTrack *track = muDst -> globalTracks(itrack);
111 
112  Int_t idTruth = track->idTruth();
113  Int_t qaTruth = track->qaTruth();
114 // Int_t id = track->id();
115 
116  if ( MUDST_CUT ) continue;
117 
118  if ( idTruth ) muTrackMap[idTruth] = track;
119  // cout << Form( "Mu Track %03i has idtruth=%03i w/ qa=%03i/100", id, idTruth, qaTruth ) << endl;
120 
121  }
122  }
123 
124  /*
125  *************************************************************************
126 
127  Obtain the geant dataset and loop over all tracks. Map g2t_track to the
128  primary key
129 
130  *************************************************************************
131  */
132 
133  TDataSet *geant = GetDataSet("geant"); assert(geant);
134  TDataSetIter geantIter(geant);
135  St_g2t_track *g2t_trackTable = (St_g2t_track *)geantIter("g2t_track");
136  g2t_track_st *trackTable = (g2t_track_st *)g2t_trackTable->GetTable();
137 
138  if ( first ) { // old fortran habits die hard
139  g2t_trackTable->Print();
140  first = false;
141  }
142 
143  cout << "-- G2T Table Event Record ----------------------------------------" << endl;
144 
145  Int_t count = 0;
146  Int_t missed = 0;
147  Int_t matched = 0;
148 
149  { Int_t ntracks = g2t_trackTable->GetNRows();
150  for ( Int_t itrack = 0; itrack<ntracks; itrack++ )
151  {
152 
153  Int_t key = trackTable[itrack].id;
154  Int_t pid = trackTable[itrack].eg_pid;
155  Int_t lab = trackTable[itrack].eg_label;
156 
157  // Skip MC track if it's a geant track
158  if ( !lab ) continue;
159 
160  // Retrieve corresponding MuDst track, skip if there is none
161  StMuTrack *muTrack = muTrackMap[ key ];
162  if ( !muTrack ) continue;
163  count++;
164 
165 // Int_t id = muTrack->id();
166 // Int_t type = muTrack->type();
167 // Int_t flag = muTrack->flag();
168  Int_t qaTruth = muTrack->qaTruth();
169  Int_t idTruth = muTrack->idTruth(); assert(idTruth==key);
170  // cout << Form(" + matched to Mu Track with id=%03i type=%03i flag=%03i qa=%03i",id,type,flag,qaTruth) << endl;
171 
172  // Get MuTrack kinematics
173  Float_t etaR = muTrack->eta();
174  Float_t phiR = muTrack->phi();
175  Float_t ptR = muTrack->pt();
176  Float_t invR = (ptR>0)? 1.0/ptR : 999;
177 
178  // Get g2tTrack kinematics
179  Float_t etaT = trackTable[itrack].eta;
180  Float_t ptT = trackTable[itrack].pt;
181  Float_t px = trackTable[itrack].p[0];
182  Float_t py = trackTable[itrack].p[1];
183  Float_t pz = trackTable[itrack].p[2];
184  Float_t invT = (ptT>0)? 1.0/ptT : 999;
185  TVector2 pT(px,py);
186  Float_t phiT = TVector2::Phi_mpi_pi( pT.Phi() );
187 
188  cout << Form("g2tTrack idtruth=%03i pid=%i px=%6.3f py=%6.3f pz=%7.3f",key,pid,px,py,pz) << endl;
189 
190  hMatchedEta -> Fill( etaT, etaR );
191  hMatchedPhi -> Fill( phiT, phiR );
192  hMatchedPt -> Fill( ptT, ptR );
193  hMatchedInv -> Fill( invT, invR );
194  hMatchedQA -> Fill( qaTruth );
195 
196  Float_t etaD = (etaR - etaT);
197  Float_t phiD = TVector2::Phi_mpi_pi(phiR - phiT);
198  Float_t ptD = ptR - ptT;
199  Float_t invD = invR - invT;
200 
201  hMatchedEtaRes -> Fill( etaD );
202  hMatchedPhiRes -> Fill( phiD );
203  hMatchedPtRes -> Fill( ptD );
204  hMatchedInvRes -> Fill( invD );
205 
206  if ( TMath::Abs( etaD ) > etaCut &&
207  TMath::Abs( phiD ) > phiCut &&
208  TMath::Abs( invD ) > invCut ) goto MISMATCH;
209 
210  matched++;
211 
212  continue;
213 
214  MISMATCH:
215 
216  missed++;
217 
218  continue;
219 
220  }
221  }
222 
223 
224  hNumMismatched -> Fill( missed );
225  if ( count > 0 ) {
226  hPerMismatched -> Fill( float(missed / count) );
227  };
228 
229  return kStOK;
230 }
231 
232 Int_t StTruthTestMaker::MakeRecord()
233 {
234 
235  static Bool_t first = true;
236 
237  /*
238  *****************************************************************
239 
240  Obtain the MuDst and map all tracks to their ID truth value
241 
242  *****************************************************************
243  */
244 
245  map< Int_t, StMuTrack * > muTrackMap;
246 
247  StMuDst *muDst = (StMuDst *)GetInputDS("MuDst"); assert(muDst);
248  StMuEvent *muEvent = muDst -> event(); assert(muEvent);
249 
250  { Int_t ntracks = muDst -> globalTracks() -> GetEntries();
251  for ( Int_t itrack = 0; itrack < ntracks; itrack++ )
252  {
253 
254  StMuTrack *track = muDst -> globalTracks(itrack);
255 
256  Int_t idTruth = track->idTruth();
257  Int_t qaTruth = track->qaTruth();
258 // Int_t id = track->id();
259 
260  if ( MUDST_CUT ) continue;
261 
262  if ( idTruth ) muTrackMap[idTruth] = track;
263  // cout << Form( "Mu Track %03i has idtruth=%03i w/ qa=%03i/100", id, idTruth, qaTruth ) << endl;
264 
265  }
266  }
267 
268 
269  /*
270  *****************************************************************
271 
272  Obtain the generator event record
273 
274  *****************************************************************
275  */
276  TObjectSet *oset = (TObjectSet *)GetData("GenEvent",".event");
277  StarGenEvent *event = (StarGenEvent *)oset->GetObject();
278  if ( first ) { // old fortran habits die har
279  event->Print();
280  first = false;
281  }
282 
283  Int_t count = 0;
284  Int_t missed = 0;
285  Int_t matched = 0;
286 
287 
288  // Loop over all event generator tracks and print those with ID truth info
289  // (e.g. the tracks which were fed out to geant for simulation)
290  TIter Next( event -> IterAll() );
291  // TIter Next( event -> IterStacked() );
292  StarGenParticle *particle = 0;
293  cout << "-- Event Generator Record ----------------------------------------" << endl;
294  while ( (particle=(StarGenParticle *)Next()) )
295  {
296  Int_t idtruth = particle->GetStack();
297 
298  if ( particle->GetPrimaryKey() != 0 ) particle->Print();
299 
300 
301  StMuTrack *muTrack = muTrackMap[idtruth];
302  if ( !muTrack ) continue;
303  count++;
304 
305  // particle->Print();
306 
307 // Int_t id = muTrack->id();
308 // Int_t type = muTrack->type();
309 // Int_t flag = muTrack->flag();
310  Int_t qaTruth = muTrack->qaTruth();
311 // Int_t idTruth = muTrack->idTruth(); // assert(idTruth==key);
312 
313  // Get MuTrack kinematics
314  Float_t etaR = muTrack->eta();
315  Float_t phiR = muTrack->phi();
316  Float_t ptR = muTrack->pt();
317  Float_t invR = (ptR>0)? 1.0/ptR : 999;
318 
319  // Get generator kinematics
320  Float_t etaT = particle->momentum().Eta();
321  Float_t ptT = particle->pt();
322  Float_t px = particle->GetPx();
323  Float_t py = particle->GetPy();
324  Float_t pz = particle->GetPz();
325  Float_t invT = (ptT>0)? 1.0/ptT : 999;
326  TVector2 pT(px,py);
327  Float_t phiT = TVector2::Phi_mpi_pi( pT.Phi() );
328  Int_t pid = particle->GetId();
329 
330  cout << Form("evtTrack idtruth=%03i pid=%i px=%6.3f py=%6.3f pz=%7.3f",idtruth ,pid,px,py,pz) << endl;
331 
332  hMatchedEta -> Fill( etaT, etaR );
333  hMatchedPhi -> Fill( phiT, phiR );
334  hMatchedPt -> Fill( ptT, ptR );
335  hMatchedInv -> Fill( invT, invR );
336  hMatchedQA -> Fill( qaTruth );
337 
338  Float_t etaD = (etaR - etaT);
339  Float_t phiD = TVector2::Phi_mpi_pi(phiR - phiT);
340  Float_t ptD = ptR - ptT;
341  Float_t invD = invR - invT;
342 
343  hMatchedEtaRes -> Fill( etaD );
344  hMatchedPhiRes -> Fill( phiD );
345  hMatchedPtRes -> Fill( ptD );
346  hMatchedInvRes -> Fill( invD );
347 
348  if ( TMath::Abs( etaD ) > etaCut &&
349  TMath::Abs( phiD ) > phiCut &&
350  TMath::Abs( invD ) > invCut ) goto MISMATCH;
351 
352  matched++;
353 
354  continue;
355 
356  MISMATCH:
357 
358  missed++;
359 
360  continue;
361 
362  }
363 
364  hNumMismatched -> Fill( missed );
365  if ( count > 0 ) {
366  hPerMismatched -> Fill( float(missed / count) );
367  }
368 
369 
370  return kStOK;
371 }
372 
StMuTrack::pt
Double_t pt() const
Returns pT at point of dca to primary vertex.
Definition: StMuTrack.h:256
StarGenParticle::pt
Float_t pt()
Returns the transverse momentum of the particle.
Definition: StarGenParticle.h:131
StarGenParticle
Yet another particle class.
Definition: StarGenParticle.h:69
StarGenParticle::GetPz
Float_t GetPz()
Get the z-component of the momentum.
Definition: StarGenParticle.h:128
StarGenParticle::GetId
Int_t GetId()
Get the id code of the particle according to the PDG standard.
Definition: StarGenParticle.h:95
StarGenParticle::GetStack
Int_t GetStack()
Get the line number in the particle stack.
Definition: StarGenParticle.h:178
StMuTrack::eta
Double_t eta() const
Returns pseudo rapidity at point of dca to primary vertex.
Definition: StMuTrack.h:257
StarGenParticle::GetPx
Float_t GetPx()
Get the x-component of the momentum.
Definition: StarGenParticle.h:120
StarGenParticle::momentum
TLorentzVector momentum()
Return the 4-momentum of the particle.
Definition: StarGenParticle.h:140
StarGenParticle::GetPy
Float_t GetPy()
Get the y-component of the momentum.
Definition: StarGenParticle.h:124
StarGenEvent
Base class for event records.
Definition: StarGenEvent.h:81
kStOK
Definition: Stypes.h:40
StMuTrack::phi
Double_t phi() const
Returns phi at point of dca to primary vertex.
Definition: StMuTrack.h:258
StarGenParticle::Print
void Print(const Option_t *opts="") const
Print the particle.
Definition: StarGenParticle.cxx:39
TObjectSet::GetObject
virtual TObject * GetObject() const
The depricated method (left here for the sake of the backward compatibility)
Definition: TObjectSet.h:56
Generated by   doxygen 1.8.5