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