]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG1/AliComparisonEff.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Histogram ranges changed to cut off saturation peak and noise
[u/mrichter/AliRoot.git] / PWG1 / AliComparisonEff.cxx
CommitLineData
09b20ad1 1//------------------------------------------------------------------------------
2// Implementation of AliComparisonEff class. It keeps information from
3// comparison of reconstructed and MC particle tracks. In addtion,
4// it keeps selection cuts used during comparison. The comparison
5// information is stored in the ROOT histograms. Analysis of these
6// histograms can be done by using Analyse() class function. The result of
7// the analysis (histograms) are stored in the output picture_eff.root file.
8//
9// Author: J.Otwinowski 04/02/2008
10//------------------------------------------------------------------------------
11
12/*
13 // after running analysis, read the file, and get component
14 gSystem->Load("libPWG1.so");
15 TFile f("Output.root");
16 AliComparisonEff * comp = (AliComparisonEff*)f.Get("AliComparisonEff");
17
18
19 // analyse comparison data (output stored in pictures_eff.root)
20 comp->Analyse();
21
22 // paramtetrisation of the TPC track length (for information only)
23 TF1 fl("fl","((min(250./(abs(x+0.000001)),250)-90))",0,2); // length function
24 TF1 fl2("fl2","[0]/((min(250./(abs(x+0.000001)),250)-90))^[1]",0,2);
25 fl2.SetParameter(1,1);
26 fl2.SetParameter(0,1);
27*/
28
29
30#include <iostream>
31
32#include "TFile.h"
33#include "TCint.h"
34#include "TH3F.h"
35#include "TH2F.h"
36#include "TF1.h"
37#include "TProfile.h"
38#include "TProfile2D.h"
39#include "TGraph2D.h"
40#include "TCanvas.h"
41#include "TGraph.h"
42//
43#include "AliESDEvent.h"
44#include "AliESD.h"
45#include "AliESDfriend.h"
46#include "AliESDfriendTrack.h"
47#include "AliRecInfoCuts.h"
48#include "AliMCInfoCuts.h"
49#include "AliLog.h"
50//
51#include "AliMathBase.h"
52#include "AliTreeDraw.h"
53#include "AliMagFMaps.h"
54#include "AliESDVertex.h"
55#include "AliExternalTrackParam.h"
56#include "AliTracker.h"
57
58#include "AliMCInfo.h"
59#include "AliESDRecInfo.h"
60#include "AliComparisonEff.h"
61
62using namespace std;
63
64
65ClassImp(AliComparisonEff)
66
67//_____________________________________________________________________________
68AliComparisonEff::AliComparisonEff():
69 TNamed("AliComparisonEff","AliComparisonEff"),
70
71 // histograms
72
73 fMCPt(0),
74 fMCRecPt(0),
75 fMCRecPrimPt(0),
76 fMCRecSecPt(0),
77
78 fEffTPCPt(0),
79 fEffTPCPtMC(0),
80 fEffTPCPtF(0),
81 //
82 fEffTPCPt_P(0),
83 fEffTPCPtMC_P(0),
84 fEffTPCPtF_P(0),
85 //
86 fEffTPCPt_Pi(0),
87 fEffTPCPtMC_Pi(0),
88 fEffTPCPtF_Pi(0),
89 //
90 fEffTPCPt_K(0),
91 fEffTPCPtMC_K(0),
92 fEffTPCPtF_K(0),
93
94 fEffTPCTan(0),
95 fEffTPCTanMC(0),
96 fEffTPCTanF(0),
97 //
98 fEffTPCPtTan(0),
99 fEffTPCPtTanMC(0),
100 fEffTPCPtTanF(0),
101
102 // Cuts
103 fCutsRC(0),
104 fCutsMC(0),
105
106 fVertex(0)
107{
108 // init vertex
109 fVertex = new AliESDVertex();
110 fVertex->SetXv(0.0); fVertex->SetYv(0.0); fVertex->SetZv(0.0);
111
112 for(Int_t i=0; i<4; ++i)
113 {
114 fTPCPtDCASigmaIdeal[i]=0;
115 fTPCPtDCASigmaFull[i]=0;
116 fTPCPtDCASigmaDay0[i]=0;
117
118 fTPCPtDCAXY[i]=0;
119 fTPCPtDCAZ[i]=0;
120
121 fTPCPtDCASigmaIdealPid[i]=0;
122 fTPCPtDCASigmaFullPid[i]=0;
123 fTPCPtDCASigmaDay0Pid[i]=0;
124
125 fTPCPtDCAXYPid[i]=0;
126 fTPCPtDCAZPid[i]=0;
127 }
128 InitHisto();
129 InitCuts();
130}
131
132//_____________________________________________________________________________
133AliComparisonEff::~AliComparisonEff(){
134
135 //
136 if(fMCPt) delete fEffTPCPt; fEffTPCPt=0;
137 if(fMCRecPt) delete fMCRecPt; fMCRecPt=0;
138 if(fMCRecPrimPt) delete fMCRecPrimPt; fMCRecPrimPt=0;
139 if(fMCRecSecPt) delete fMCRecSecPt; fMCRecSecPt=0;
140
141 //
142 if(fEffTPCPt) delete fEffTPCPt; fEffTPCPt=0;
143 if(fEffTPCPtMC) delete fEffTPCPtMC; fEffTPCPtMC=0;
144 if(fEffTPCPtF) delete fEffTPCPtF; fEffTPCPtF=0;
145
146 //
147 if(fEffTPCPt_P) delete fEffTPCPt_P; fEffTPCPt_P=0;
148 if(fEffTPCPtMC_P) delete fEffTPCPtMC_P; fEffTPCPtMC_P=0;
149 if(fEffTPCPtF_P) delete fEffTPCPtF_P; fEffTPCPtF_P=0;
150
151 //
152 if(fEffTPCPt_Pi) delete fEffTPCPt_Pi; fEffTPCPt_Pi=0;
153 if(fEffTPCPtMC_Pi) delete fEffTPCPtMC_Pi; fEffTPCPtMC_Pi=0;
154 if(fEffTPCPtF_Pi) delete fEffTPCPtF_Pi; fEffTPCPtF_Pi=0;
155
156 //
157 if(fEffTPCPt_K) delete fEffTPCPt_K; fEffTPCPt_K=0;
158 if(fEffTPCPtMC_K) delete fEffTPCPtMC_K; fEffTPCPtMC_K=0;
159 if(fEffTPCPtF_K) delete fEffTPCPtF_K; fEffTPCPtF_K=0;
160
161 //
162 if(fEffTPCTan) delete fEffTPCTan; fEffTPCTan=0;
163 if(fEffTPCTanMC) delete fEffTPCTanMC; fEffTPCTanMC=0;
164 if(fEffTPCTanF) delete fEffTPCTanF; fEffTPCTanF=0;
165
166 //
167 if(fEffTPCPtTan) delete fEffTPCPtTan; fEffTPCPtTan=0;
168 if(fEffTPCPtTanMC) delete fEffTPCPtTanMC; fEffTPCPtTanMC=0;
169 if(fEffTPCPtTanF) delete fEffTPCPtTanF; fEffTPCPtTanF=0;
170
171 for(Int_t i=0; i<4; ++i)
172 {
173 if(fTPCPtDCASigmaIdeal[i]) delete fTPCPtDCASigmaIdeal[i]; fTPCPtDCASigmaIdeal[i]=0;
174 if(fTPCPtDCASigmaFull[i]) delete fTPCPtDCASigmaFull[i]; fTPCPtDCASigmaFull[i]=0;
175 if(fTPCPtDCASigmaDay0[i]) delete fTPCPtDCASigmaDay0[i]; fTPCPtDCASigmaDay0[i]=0;
176
177 if(fTPCPtDCAXY[i]) delete fTPCPtDCAXY[i]; fTPCPtDCAXY[i]=0;
178 if(fTPCPtDCAZ[i]) delete fTPCPtDCAZ[i]; fTPCPtDCAZ[i]=0;
179
180 if(fTPCPtDCASigmaIdealPid[i]) delete fTPCPtDCASigmaIdealPid[i]; fTPCPtDCASigmaIdealPid[i]=0;
181 if(fTPCPtDCASigmaFullPid[i]) delete fTPCPtDCASigmaFullPid[i]; fTPCPtDCASigmaFullPid[i]=0;
182 if(fTPCPtDCASigmaDay0Pid[i]) delete fTPCPtDCASigmaDay0Pid[i]; fTPCPtDCASigmaDay0Pid[i]=0;
183
184 if(fTPCPtDCAXYPid[i]) delete fTPCPtDCAXYPid[i]; fTPCPtDCAXYPid[i]=0;
185 if(fTPCPtDCAZPid[i]) delete fTPCPtDCAZPid[i]; fTPCPtDCAZPid[i]=0;
186 }
187}
188
189//_____________________________________________________________________________
190void AliComparisonEff::InitHisto(){
191
192 // Init histograms
193 //
194 fMCPt = new TH1F("fMCPt","fMCPt",50,0.1,3);
195 fMCPt->SetXTitle("p_{t}");
196 fMCPt->SetYTitle("yield");
197
198 fMCRecPt = new TH1F("fMCRecPt","fMCRecPt",50,0.1,3);
199 fMCRecPt->SetXTitle("p_{t}");
200 fMCRecPt->SetYTitle("yield");
201
202 fMCRecPrimPt = new TH1F("fMCRecPrimPt","fMCRecPrimPt",50,0.1,3);
203 fMCRecPrimPt->SetXTitle("p_{t}");
204 fMCRecPrimPt->SetYTitle("yield");
205
206 fMCRecSecPt = new TH1F("fMCRecSecPt","fMCRecSecPt",50,0.1,3);
207 fMCRecSecPt->SetXTitle("p_{t}");
208 fMCRecSecPt->SetYTitle("yield");
209
210 // Efficiency as function of pt
211 fEffTPCPt = new TProfile("Eff_pt","Eff_Pt",50,0.1,3);
212 fEffTPCPt->SetXTitle("p_{t}");
213 fEffTPCPt->SetYTitle("TPC Efficiency");
214
215 fEffTPCPtMC = new TProfile("MC_Eff_pt","MC_Eff_Pt",50,0.1,3);
216 fEffTPCPtMC->SetXTitle("p_{t}");
217 fEffTPCPtMC->SetYTitle("MC TPC Efficiency");
218
219 fEffTPCPtF = new TProfile("F_Eff_pt","F_Eff_Pt",50,0.1,3);
220 fEffTPCPtF->SetXTitle("p_{t}");
221 fEffTPCPtF->SetYTitle("TPC Findable Efficiency");
222
223 // Efficiency as function of pt protons
224 fEffTPCPt_P = new TProfile("Eff_pt_P","Eff_Pt_P",50,0.1,3);
225 fEffTPCPt_P->SetXTitle("p_{t}");
226 fEffTPCPt_P->SetYTitle("TPC Efficiency");
227
228 fEffTPCPtMC_P = new TProfile("MC_Eff_pt_P","MC_Eff_Pt_P",50,0.1,3);
229 fEffTPCPtMC_P->SetXTitle("p_{t}");
230 fEffTPCPtMC_P->SetYTitle("MC TPC Efficiency");
231
232 fEffTPCPtF_P = new TProfile("F_Eff_pt_P","F_Eff_Pt_P",50,0.1,3);
233 fEffTPCPtF_P->SetXTitle("p_{t}");
234 fEffTPCPtF_P->SetYTitle("TPC Findable Efficiency");
235
236 // Efficiency as function of pt pions
237 fEffTPCPt_Pi = new TProfile("Eff_pt_Pi","Eff_Pit_Pi",50,0.1,3);
238 fEffTPCPt_Pi->SetXTitle("p_{t}");
239 fEffTPCPt_Pi->SetYTitle("TPC Efficiency");
240
241 fEffTPCPtMC_Pi = new TProfile("MC_Eff_pt_Pi","MC_Eff_Pit_Pi",50,0.1,3);
242 fEffTPCPtMC_Pi->SetXTitle("p_{t}");
243 fEffTPCPtMC_Pi->SetYTitle("MC TPC Efficiency");
244
245 fEffTPCPtF_Pi = new TProfile("F_Eff_pt_Pi","F_Eff_Pit_Pi",50,0.1,3);
246 fEffTPCPtF_Pi->SetXTitle("p_{t}");
247 fEffTPCPtF_Pi->SetYTitle("TPC Findable Efficiency");
248
249 // Efficiency as function of pt kaons
250 fEffTPCPt_K = new TProfile("Eff_pt_K","Eff_Kt_K",50,0.1,3);
251 fEffTPCPt_K->SetXTitle("p_{t}");
252 fEffTPCPt_K->SetYTitle("TPC Efficiency");
253
254 fEffTPCPtMC_K = new TProfile("MC_Eff_pt_K","MC_Eff_Kt_K",50,0.1,3);
255 fEffTPCPtMC_K->SetXTitle("p_{t}");
256 fEffTPCPtMC_K->SetYTitle("MC TPC Efficiency");
257
258 fEffTPCPtF_K = new TProfile("F_Eff_pt_K","F_Eff_Kt_K",50,0.1,3);
259 fEffTPCPtF_K->SetXTitle("p_{t}");
260 fEffTPCPtF_K->SetYTitle("TPC Findable Efficiency");
261
262 // Efficiency as function of tan(theta)
263 fEffTPCTan = new TProfile("Eff_tan","Eff_tan",50,-2.5,2.5);
264 fEffTPCTan->SetXTitle("tan(#theta)");
265 fEffTPCTan->SetYTitle("TPC Efficiency");
266
267 fEffTPCTanMC = new TProfile("MC_Eff_tan","MC_Eff_tan",50,-2.5,2.5);
268 fEffTPCTanMC->SetXTitle("tan(#theta)");
269 fEffTPCTanMC->SetYTitle("MC TPC Efficiency");
270
271 fEffTPCTanF = new TProfile("F_Eff_tan","F_Eff_tan",50,-2.5,2.5);
272 fEffTPCPtF->SetXTitle("tan(#theta)");
273 fEffTPCPtF->SetYTitle("TPC Findable Efficiency");
274
275 // Efficiency as function of pt and tan(theta)
276 fEffTPCPtTan = new TProfile2D("Eff_pt_tan","Eff_Pt_tan",10,0.1,3,20,-2.,2.);
277 fEffTPCPtTan->SetXTitle("tan(#theta)");
278 fEffTPCPtTan->SetYTitle("p_{t}");
279
280 fEffTPCPtTanMC = new TProfile2D("MC_Eff_pt_tan_MC","MC Eff Pt",10,0.1,3,20,-2.,2.);
281 fEffTPCPtTanMC->SetXTitle("tan(#theta)");
282 fEffTPCPtTanMC->SetYTitle("p_{t}");
283
284 fEffTPCPtTanF = new TProfile2D("MC_Eff_pt_tan_F","MC Eff Pt",10,0.1,3,20,-2.,2.);
285 fEffTPCPtTanF->SetXTitle("tan(#theta)");
286 fEffTPCPtTanF->SetYTitle("p_{t}");
287
288 char name[256];
289 for(Int_t i=0; i<4; ++i)
290 {
291 sprintf(name, "fTPCPtDCASigmaIdeal_%d",i);
292 fTPCPtDCASigmaIdeal[i] = new TH2F(name,name,50,0.1,3,100,0,100);
293
294 sprintf(name, "fTPCPtDCASigmaFull_%d",i);
295 fTPCPtDCASigmaFull[i] = new TH2F(name,name,50,0.1,3,100,0,100);
296
297 sprintf(name, "fTPCPtDCASigmaDay0_%d",i);
298 fTPCPtDCASigmaDay0[i] = new TH2F(name,name,50,0.1,3,100,0,100);
299
300 sprintf(name, "fTPCPtDCAXY_%d",i);
301 fTPCPtDCAXY[i]= new TH2F(name,name,50,0.1,3,100,0,100);
302 sprintf(name, "fTPCPtDCAZ_%d",i);
303 fTPCPtDCAZ[i]= new TH2F(name,name,50,0.1,3,100,0,100);
304
305 sprintf(name, "fTPCPtDCASigmaIdealPid_%d",i);
306 fTPCPtDCASigmaIdealPid[i] = new TH3F(name,name,50,0.1,3,100,0,100,5,0,5);
307
308 sprintf(name, "fTPCPtDCASigmaFullPid_%d",i);
309 fTPCPtDCASigmaFullPid[i]= new TH3F(name,name,50,0.1,3,100,0,100,5,0,5);
310
311 sprintf(name, "fTPCPtDCASigmaDay0Pid_%d",i);
312 fTPCPtDCASigmaDay0Pid[i]= new TH3F(name,name,50,0.1,3,100,0,100,5,0,5);
313
314 sprintf(name, "fTPCPtDCAXYPid_%d",i);
315 fTPCPtDCAXYPid[i]= new TH3F(name,name,50,0.1,3,100,0,100,5,0,5);
316
317 sprintf(name, "fTPCPtDCAZPid_%d",i);
318 fTPCPtDCAZPid[i]= new TH3F(name,name,50,0.1,3,100,0,100,5,0,5);
319 }
320}
321
322//_____________________________________________________________________________
323void AliComparisonEff::InitCuts()
324{
325
326 if(!fCutsMC)
327 AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
328 if(!fCutsRC)
329 AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
330}
331
332//_____________________________________________________________________________
333void AliComparisonEff::Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC)
334{
335 // Fill efficiency comparison information
336
337 AliExternalTrackParam *track = 0;
338 Double_t kRadius = 3.0; // beam pipe radius
339 Double_t kMaxStep = 5.0; // max step
340 Double_t field = 0.4; // mag. field
341 Double_t kMaxD = 123456.0; // max distance
342
343 // systematics
344 const Double_t kSigma2Full_xy = 0.25; // ExB full systematics [cm]
345 const Double_t kSigma2Full_z = 5.0; // [cm]
346
347 const Double_t kSigma2Day0_xy = 0.02; // ExB [cm]
348 const Double_t kSigma2Day0_z = 0.1; // [cm] goofie
349
350 //
351 Double_t DCASigmaIdeal=0;
352 Double_t DCASigmaFull=0;
353 Double_t DCASigmaDay0=0;
354
355 Double_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
356 const Double_t* dv;
357
358 // distance to Prim. vertex
359 dv = infoMC->GetVDist();
360
361 Float_t mcpt = infoMC->GetParticle().Pt();
362 Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);
363 //Bool_t isPrim = infoMC->GetParticle().R()<fCutsMC->GetMaxR() && TMath::Abs(infoMC->GetParticle().Vz())<fCutsMC->GetMaxVz();
364 Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz();
365
366 // calculate and set prim. vertex
367 fVertex->SetXv( infoMC->GetParticle().Vx() - dv[0] );
368 fVertex->SetYv( infoMC->GetParticle().Vy() - dv[1] );
369 fVertex->SetZv( infoMC->GetParticle().Vz() - dv[2] );
370
371 // Check selection cuts
372 if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return;
373
374 // transform Pdg to Pid
375 Double_t pid = -1;
376 if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetEP() ) pid = 0;
377 if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetMuP() ) pid = 1;
378 if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetPiP() ) pid = 2;
379 if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetKP() ) pid = 3;
380 if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetProt() ) pid = 4;
381
382 //cout << "dv[0] " << dv[0] << " dv[1] " << dv[1] << " dv[2] " << dv[2] << endl;
383 //cout << "v[0] " << fVertex->GetXv() << " v[1] " << fVertex->GetYv() << " v[2] " << fVertex->GetZv()<< endl;
384 if (TMath::Abs(tantheta)<fCutsRC->GetMaxAbsTanTheta())
385 {
386
387 // MC pt
388 fMCPt->Fill(mcpt);
389
390
391 if (infoRC->GetESDtrack()->GetTPCInnerParam())
392 {
393 if ((track = new AliExternalTrackParam(*infoRC->GetESDtrack()->GetTPCInnerParam())) != 0 )
394 {
395 AliTracker::PropagateTrackTo(track,kRadius,infoMC->GetMass(),kMaxStep,kTRUE);
396 track->PropagateToDCA(fVertex,field,kMaxD,dca,cov);
397
398 //
399 cov[2] = track->GetCovariance()[2];
400
401 // Eff = infoRC->GetStatus(1)==3 && isPrim / isPrim;
402 // Pt vs ( dca[0]^2/cov[0]^2 + dca[1]^2/cov[2]^2 )
403 // Pt vs ( dca[0]^2/(cov[0]^2 + kSigma2Full_xy) + dca[1]^2/( cov[2]^2 + kSigma2Full_z )
404 // Pt vs ( dca[0]^2/(cov[0]^2 + kSigma2_xy) + dca[1]^2/( cov[2]^2 + kSigma2_z )
405
406 if(cov[0]>0.0 && cov[2]>0.0)
407 {
408 DCASigmaIdeal = TMath::Power(dca[0],2)/cov[0]
409 + TMath::Power(dca[1],2)/cov[2];
410
411 DCASigmaFull = TMath::Power(dca[0],2)/(cov[0]+kSigma2Full_xy)
412 + TMath::Power(dca[1],2)/(cov[2]+kSigma2Full_z);
413
414 DCASigmaDay0 = TMath::Power(dca[0],2)/(cov[0]+kSigma2Day0_xy)
415 + TMath::Power(dca[1],2)/(cov[2]+kSigma2Day0_z);
416
417 //cout << "dca[0] " << dca[0] << " dca[1] " << dca[1] << endl;
418 //cout << "cov[0] " << cov[0] << " cov[2] " << cov[2] << endl;
419 //cout << "DCASigmaIdeal " << DCASigmaIdeal << " DCASigmaFull " << DCASigmaFull << " DCASigmaDay0 " <<DCASigmaDay0 << endl;
420 //cout << " -------------------------------------------------------- "<< endl;
421 }
422
423 if(infoRC->GetStatus(1)==3) fMCRecPt->Fill(mcpt);
424 if(infoRC->GetStatus(1)==3 && isPrim) fMCRecPrimPt->Fill(mcpt);
425 if(infoRC->GetStatus(1)==3 && !isPrim) fMCRecSecPt->Fill(mcpt);
426
427 if(isPrim)
428 {
429 fTPCPtDCASigmaIdeal[0]->Fill(mcpt,DCASigmaIdeal);
430 fTPCPtDCASigmaFull[0]->Fill(mcpt,DCASigmaFull);
431 fTPCPtDCASigmaDay0[0]->Fill(mcpt,DCASigmaDay0);
432
433 fTPCPtDCAXY[0]->Fill(mcpt,dca[0]);
434 fTPCPtDCAZ[0]->Fill(mcpt,dca[1]);
435
436 fTPCPtDCASigmaIdealPid[0]->Fill(mcpt,DCASigmaIdeal,pid);
437 fTPCPtDCASigmaFullPid[0]->Fill(mcpt,DCASigmaFull,pid);
438 fTPCPtDCASigmaDay0Pid[0]->Fill(mcpt,DCASigmaDay0,pid);
439
440 fTPCPtDCAXYPid[0]->Fill(mcpt,dca[0],pid);
441 fTPCPtDCAZPid[0]->Fill(mcpt,dca[1],pid);
442
443 if(infoRC->GetStatus(1)==3)
444 {
445 fTPCPtDCASigmaIdeal[1]->Fill(mcpt,DCASigmaIdeal);
446 fTPCPtDCASigmaFull[1]->Fill(mcpt,DCASigmaFull);
447 fTPCPtDCASigmaDay0[1]->Fill(mcpt,DCASigmaDay0);
448
449 fTPCPtDCAXY[1]->Fill(mcpt,dca[0]);
450 fTPCPtDCAZ[1]->Fill(mcpt,dca[1]);
451
452 fTPCPtDCASigmaIdealPid[1]->Fill(mcpt,DCASigmaIdeal,pid);
453 fTPCPtDCASigmaFullPid[1]->Fill(mcpt,DCASigmaFull,pid);
454 fTPCPtDCASigmaDay0Pid[1]->Fill(mcpt,DCASigmaDay0,pid);
455
456 fTPCPtDCAXYPid[1]->Fill(mcpt,dca[0],pid);
457 fTPCPtDCAZPid[1]->Fill(mcpt,dca[1],pid);
458 }
459 }
460
461 // Cont = infoRC->GetStatus(1)==3 && !isPrim / infoRC->GetStatus(1)==3
462 // Pt vs ( dz[0]^2/cov[0]^2 + dz[1]^2/cov[2]^2 )
463 // Pt vs ( dz[0]^2/(cov[0]^2 + kSigma2Full_xy) + dz[1]^2/( cov[2]^2 + kSigma2Full_z )
464 // Pt vs ( dz[0]^2/(cov[0]^2 + kSigma2_xy) + dz[1]^2/( cov[2]^2 + kSigma2_z )
465
466 if(infoRC->GetStatus(1)==3)
467 {
468 fTPCPtDCASigmaIdeal[2]->Fill(mcpt,DCASigmaIdeal);
469 fTPCPtDCASigmaFull[2]->Fill(mcpt,DCASigmaFull);
470 fTPCPtDCASigmaDay0[2]->Fill(mcpt,DCASigmaDay0);
471
472 fTPCPtDCAXY[2]->Fill(mcpt,dca[0]);
473 fTPCPtDCAZ[2]->Fill(mcpt,dca[1]);
474
475 fTPCPtDCASigmaIdealPid[2]->Fill(mcpt,DCASigmaIdeal,pid);
476 fTPCPtDCASigmaFullPid[2]->Fill(mcpt,DCASigmaFull,pid);
477 fTPCPtDCASigmaDay0Pid[2]->Fill(mcpt,DCASigmaDay0,pid);
478
479 fTPCPtDCAXYPid[2]->Fill(mcpt,dca[0],pid);
480 fTPCPtDCAZPid[2]->Fill(mcpt,dca[1],pid);
481
482 if(isPrim==0)
483 {
484 fTPCPtDCASigmaIdeal[3]->Fill(mcpt,DCASigmaIdeal);
485 fTPCPtDCASigmaFull[3]->Fill(mcpt,DCASigmaFull);
486 fTPCPtDCASigmaDay0[3]->Fill(mcpt,DCASigmaDay0);
487
488 fTPCPtDCAXY[3]->Fill(mcpt,dca[0]);
489 fTPCPtDCAZ[3]->Fill(mcpt,dca[1]);
490
491 fTPCPtDCASigmaIdealPid[3]->Fill(mcpt,DCASigmaIdeal,pid);
492 fTPCPtDCASigmaFullPid[3]->Fill(mcpt,DCASigmaFull,pid);
493 fTPCPtDCASigmaDay0Pid[3]->Fill(mcpt,DCASigmaDay0,pid);
494
495 fTPCPtDCAXYPid[3]->Fill(mcpt,dca[0],pid);
496 fTPCPtDCAZPid[3]->Fill(mcpt,dca[1],pid);
497 }
498 }
499 delete track;
500 }
501 }
502 else
503 {
504 if(isPrim)
505 {
506 fTPCPtDCASigmaIdeal[0]->Fill(mcpt,0.0);
507 fTPCPtDCASigmaFull[0]->Fill(mcpt,0.0);
508 fTPCPtDCASigmaDay0[0]->Fill(mcpt,0.0);
509
510 fTPCPtDCAXY[0]->Fill(mcpt,0.0);
511 fTPCPtDCAZ[0]->Fill(mcpt,0.0);
512
513 fTPCPtDCASigmaIdealPid[0]->Fill(mcpt,0.0,pid);
514 fTPCPtDCASigmaFullPid[0]->Fill(mcpt,0.0,pid);
515 fTPCPtDCASigmaDay0Pid[0]->Fill(mcpt,0.0,pid);
516
517 fTPCPtDCAXYPid[0]->Fill(mcpt,0.0,pid);
518 fTPCPtDCAZPid[0]->Fill(mcpt,0.0,pid);
519 }
520 }
521 }
522
523 // only primary particles
524 if (!isPrim) return;
525
526 // pt
527 if (TMath::Abs(tantheta)<fCutsRC->GetMaxAbsTanTheta()){
528
529 fEffTPCPt->Fill(mcpt, infoRC->GetStatus(1)==3);
530 fEffTPCPtMC->Fill(mcpt, infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits());
531 if (infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits()){
532 fEffTPCPtF->Fill(mcpt, infoRC->GetStatus(1)==3);
533 }
534
535 // protons
536 if(TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetProt()) {
537 fEffTPCPt_P->Fill(mcpt, infoRC->GetStatus(1)==3);
538 fEffTPCPtMC_P->Fill(mcpt, infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits());
539
540 if(infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits()) {
541 fEffTPCPtF_P->Fill(mcpt, infoRC->GetStatus(1)==3);
542 }
543 }
544
545 // pions
546 if(TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetPiP()) {
547 fEffTPCPt_Pi->Fill(mcpt, infoRC->GetStatus(1)==3);
548 fEffTPCPtMC_Pi->Fill(mcpt, infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits());
549
550 if(infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits()) {
551 fEffTPCPtF_Pi->Fill(mcpt, infoRC->GetStatus(1)==3);
552 }
553 }
554
555 // kaons
556 if(TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetKP()) {
557 fEffTPCPt_K->Fill(mcpt, infoRC->GetStatus(1)==3);
558 fEffTPCPtMC_K->Fill(mcpt, infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits());
559
560 if(infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits()) {
561 fEffTPCPtF_K->Fill(mcpt, infoRC->GetStatus(1)==3);
562 }
563 }
564 }
565
566 // theta
567 if (TMath::Abs(mcpt)>fCutsRC->GetPtMin()){
568 fEffTPCTan->Fill(tantheta, infoRC->GetStatus(1)==3);
569 fEffTPCTanMC->Fill(tantheta, infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits());
570 if (infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits()){
571 fEffTPCTanF->Fill(tantheta, infoRC->GetStatus(1)==3);
572 }
573 }
574
575 // pt-theta
576 fEffTPCPtTan->Fill(mcpt,tantheta,infoRC->GetStatus(1)==3);
577 fEffTPCPtTanMC->Fill(mcpt,tantheta,infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits());
578 if (infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits()){
579 fEffTPCPtTanF->Fill(mcpt,tantheta,infoRC->GetStatus(1)==3);
580 }
581}
582
583//_____________________________________________________________________________
584void AliComparisonEff::Exec(AliMCInfo* infoMC, AliESDRecInfo *infoRC)
585{
586 // Process comparison information
587 Process(infoMC,infoRC);
588}
589
590//_____________________________________________________________________________
591Long64_t AliComparisonEff::Merge(TCollection* list)
592{
593 // Merge list of objects (needed by PROOF)
594
595 if (!list)
596 return 0;
597
598 if (list->IsEmpty())
599 return 1;
600
601 TIterator* iter = list->MakeIterator();
602 TObject* obj = 0;
603
604 // collection of generated histograms
605
606 Int_t count=0;
607 while((obj = iter->Next()) != 0)
608 {
609 AliComparisonEff* entry = dynamic_cast<AliComparisonEff*>(obj);
610 if (entry == 0)
611 continue;
612
613 fMCPt->Add(entry->fMCPt);
614 fMCRecPt->Add(entry->fMCRecPt);
615 fMCRecPrimPt->Add(entry->fMCRecPrimPt);
616 fMCRecSecPt->Add(entry->fMCRecSecPt);
617
618 fEffTPCPt->Add(entry->fEffTPCPt);
619 fEffTPCPtMC->Add(entry->fEffTPCPtMC);
620 fEffTPCPtF->Add(entry->fEffTPCPtF);
621
622 fEffTPCPt_P->Add(entry->fEffTPCPt_P);
623 fEffTPCPtMC_P->Add(entry->fEffTPCPtMC_P);
624 fEffTPCPtF_P->Add(entry->fEffTPCPtF_P);
625
626 fEffTPCPt_Pi->Add(entry->fEffTPCPt_Pi);
627 fEffTPCPtMC_Pi->Add(entry->fEffTPCPtMC_Pi);
628 fEffTPCPtF_Pi->Add(entry->fEffTPCPtF_Pi);
629
630 fEffTPCPt_K->Add(entry->fEffTPCPt_K);
631 fEffTPCPtMC_K->Add(entry->fEffTPCPtMC_K);
632 fEffTPCPtF_K->Add(entry->fEffTPCPtF_K);
633
634 fEffTPCTan->Add(entry->fEffTPCTan);
635 fEffTPCTanMC->Add(entry->fEffTPCTanMC);
636 fEffTPCTanF->Add(entry->fEffTPCTanF);
637
638 fEffTPCPtTan->Add(entry->fEffTPCPtTan);
639 fEffTPCPtTanMC->Add(entry->fEffTPCPtTanMC);
640 fEffTPCPtTanF->Add(entry->fEffTPCPtTanF);
641
642 for(Int_t i=0; i<4; ++i)
643 {
644 fTPCPtDCASigmaIdeal[i]->Add(entry->fTPCPtDCASigmaIdeal[i]);
645 fTPCPtDCASigmaFull[i]->Add(entry->fTPCPtDCASigmaFull[i]);
646 fTPCPtDCASigmaDay0[i]->Add(entry->fTPCPtDCASigmaDay0[i]);
647
648 fTPCPtDCAXY[i]->Add(entry->fTPCPtDCAXY[i]);
649 fTPCPtDCAZ[i]->Add(entry->fTPCPtDCAXY[i]);
650
651 fTPCPtDCASigmaIdealPid[i]->Add(entry->fTPCPtDCASigmaIdealPid[i]);
652 fTPCPtDCASigmaFullPid[i]->Add(entry->fTPCPtDCASigmaFullPid[i]);
653 fTPCPtDCASigmaDay0Pid[i]->Add(entry->fTPCPtDCASigmaDay0Pid[i]);
654
655 fTPCPtDCAXYPid[i]->Add(entry->fTPCPtDCAXYPid[i]);
656 fTPCPtDCAZPid[i]->Add(entry->fTPCPtDCAXYPid[i]);
657 }
658
659 count++;
660 }
661
662return count;
663}
664
665//_____________________________________________________________________________
666void AliComparisonEff::Analyse()
667{
668 // Analyse output histograms
669
670 AliComparisonEff * comp=this;
671
672 // calculate efficiency and contamination (4 sigma)
673
674 TH1 *h_sigmaidealpid[20];
675 TH1 *h_sigmafullpid[20];
676 TH1 *h_sigmaday0pid[20];
677
678 //TH1 *h_sigmaday0pidclone[20];
679
680 TH1 *h_sigmaidealpidtot[4];
681 TH1 *h_sigmafullpidtot[4];
682 TH1 *h_sigmaday0pidtot[4];
683
684 //TH1 *h_sigmaday0pidtotclone[4];
685
686 char name[256];
687 char name1[256];
688 Int_t idx;
689
690 for(Int_t i=0; i<4; ++i)
691 {
692 //total
693 comp->fTPCPtDCASigmaIdealPid[i]->GetYaxis()->SetRange(1,4);
694 sprintf(name,"h_sigmaidealpidtot_%d",i);
695 h_sigmaidealpidtot[i] = comp->fTPCPtDCASigmaIdealPid[i]->Project3D();
696 h_sigmaidealpidtot[i]->SetName(name);
697
698 comp->fTPCPtDCASigmaFullPid[i]->GetYaxis()->SetRange(1,4);
699 sprintf(name,"h_sigmafullpidtot_%d",i);
700 h_sigmafullpidtot[i] = comp->fTPCPtDCASigmaFullPid[i]->Project3D();
701 h_sigmafullpidtot[i]->SetName(name);
702
703 comp->fTPCPtDCASigmaDay0Pid[i]->GetYaxis()->SetRange(1,4);
704 sprintf(name,"h_sigmaday0pidtot_%d",i);
705 h_sigmaday0pidtot[i] = comp->fTPCPtDCASigmaDay0Pid[i]->Project3D();
706 h_sigmaday0pidtot[i]->SetName(name);
707
708 // pid wise
709 for(Int_t j=0; j<5; ++j)
710 {
711 idx = i*5 + j;
712
713 comp->fTPCPtDCASigmaIdealPid[i]->GetYaxis()->SetRange(1,4);
714 comp->fTPCPtDCASigmaIdealPid[i]->GetZaxis()->SetRange(j+1,j+1);
715
716 sprintf(name,"h_sigmaidealpid_%d",idx);
717 h_sigmaidealpid[idx] = comp->fTPCPtDCASigmaIdealPid[i]->Project3D();
718 h_sigmaidealpid[idx]->SetName(name);
719
720
721 comp->fTPCPtDCASigmaFullPid[i]->GetYaxis()->SetRange(1,4);
722 comp->fTPCPtDCASigmaFullPid[i]->GetZaxis()->SetRange(j+1,j+1);
723
724 sprintf(name,"h_sigmafullpid_%d",idx);
725 h_sigmafullpid[idx] = comp->fTPCPtDCASigmaFullPid[i]->Project3D();
726 h_sigmafullpid[idx]->SetName(name);
727
728 comp->fTPCPtDCASigmaDay0Pid[i]->GetYaxis()->SetRange(1,4);
729 comp->fTPCPtDCASigmaDay0Pid[i]->GetZaxis()->SetRange(j+1,j+1);
730
731 sprintf(name,"h_sigmaday0pid_%d",idx);
732 h_sigmaday0pid[idx] = comp->fTPCPtDCASigmaDay0Pid[i]->Project3D();
733 h_sigmaday0pid[idx]->SetName(name);
734
735 }
736 }
737
738 // calculate efficiency and contamination (all pids)
739 h_sigmaidealpidtot[0]->Sumw2();
740 h_sigmaidealpidtot[1]->Divide(h_sigmaidealpidtot[0]);
741 h_sigmaidealpidtot[2]->Sumw2();
742 h_sigmaidealpidtot[3]->Divide(h_sigmaidealpidtot[2]);
743
744 h_sigmafullpidtot[0]->Sumw2();
745 h_sigmafullpidtot[1]->Divide(h_sigmafullpidtot[0]);
746 h_sigmafullpidtot[2]->Sumw2();
747 h_sigmafullpidtot[3]->Divide(h_sigmafullpidtot[2]);
748
749 h_sigmaday0pidtot[0]->Sumw2();
750 h_sigmaday0pidtot[1]->Divide(h_sigmaday0pidtot[0]);
751 h_sigmaday0pidtot[2]->Sumw2();
752 h_sigmaday0pidtot[3]->Divide(h_sigmaday0pidtot[2]);
753
754 // calculate efficiency pid wise
755 for(Int_t idx = 0; idx<5; idx++)
756 {
757 h_sigmaidealpid[idx]->Sumw2();
758 h_sigmaidealpid[idx+5]->Divide(h_sigmaidealpid[idx]);
759
760 h_sigmafullpid[idx]->Sumw2();
761 h_sigmafullpid[idx+5]->Divide(h_sigmafullpid[idx]);
762
763 h_sigmaday0pid[idx]->Sumw2();
764 h_sigmaday0pid[idx+5]->Divide(h_sigmaday0pid[idx]);
765 }
766
767 // calculate cont. pid wise
768 for(Int_t idx = 0; idx<5; idx++)
769 {
770 h_sigmaidealpid[idx+15]->Divide(h_sigmaidealpidtot[2]);
771 h_sigmafullpid[idx+15]->Divide(h_sigmafullpidtot[2]);
772 h_sigmaday0pid[idx+15]->Divide(h_sigmaday0pidtot[2]);
773 }
774
775 // write results
776 TFile *fp = new TFile("pictures_eff.root","recreate");
777 fp->cd();
778
779 TCanvas * c = new TCanvas("Efficiency","Track efficiency");
780 c->cd();
781
782 fMCPt->Write();
783 fMCRecPt->Write();
784 fMCRecPrimPt->Write();
785 fMCRecSecPt->Write();
786
787 for(int i = 0; i<4;i++)
788 {
789 comp->fTPCPtDCASigmaIdealPid[i]->GetYaxis()->SetRange();
790 comp->fTPCPtDCASigmaIdealPid[i]->GetZaxis()->SetRange();
791 comp->fTPCPtDCASigmaFullPid[i]->GetYaxis()->SetRange();
792 comp->fTPCPtDCASigmaFullPid[i]->GetZaxis()->SetRange();
793 comp->fTPCPtDCASigmaDay0Pid[i]->GetYaxis()->SetRange();
794 comp->fTPCPtDCASigmaDay0Pid[i]->GetZaxis()->SetRange();
795
796 comp->fTPCPtDCASigmaIdealPid[i]->Write();
797 comp->fTPCPtDCASigmaFullPid[i]->Write();
798 comp->fTPCPtDCASigmaDay0Pid[i]->Write();
799 }
800 //
801 comp->fEffTPCTanF->SetXTitle("Tan(#theta)");
802 comp->fEffTPCTanF->SetYTitle("eff_{findable}");
803 comp->fEffTPCTanF->Write("EffTanFindable");
804 //
805 comp->fEffTPCTan->SetXTitle("Tan(#theta)");
806 comp->fEffTPCTan->SetYTitle("eff_{all}");
807 comp->fEffTPCTan->Write("EffTanAll");
808
809 h_sigmaidealpidtot[1]->Write("Eff_SigmaIdeal");
810 h_sigmaidealpidtot[3]->Write("Cont_SigmaIdeal");
811
812 h_sigmafullpidtot[1]->Write("Eff_SigmaFull");
813 h_sigmafullpidtot[3]->Write("Cont_SigmaFull");
814
815 h_sigmaday0pidtot[1]->Write("Eff_SigmaDay0");
816 h_sigmaday0pidtot[3]->Write("Cont_SigmaDay0");
817
818 for(Int_t idx = 0; idx<5; idx++)
819 {
820 sprintf(name,"Eff_SigmaIdeal_%d",idx);
821 sprintf(name1,"Cont_SigmaIdeal_%d",idx);
822
823 h_sigmaidealpid[idx+5]->Write(name);
824 h_sigmaidealpid[idx+15]->Write(name1);
825
826 sprintf(name,"Eff_SigmaFull_%d",idx);
827 sprintf(name1,"Cont_SigmaFull_%d",idx);
828
829 h_sigmafullpid[idx+5]->Write(name);
830 h_sigmafullpid[idx+15]->Write(name1);
831
832 sprintf(name,"Eff_SigmaDay0_%d",idx);
833 sprintf(name1,"Cont_SigmaDay0_%d",idx);
834
835 h_sigmaday0pid[idx+5]->Write(name);
836 h_sigmaday0pid[idx+15]->Write(name1);
837 }
838
839 //
840 fp->Close();
841}