]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG1/AliComparisonRes.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Update of Macros
[u/mrichter/AliRoot.git] / PWG1 / AliComparisonRes.cxx
CommitLineData
96b3109a 1//------------------------------------------------------------------------------
2// Implementation of AliComparisonRes 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_res.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 AliComparisonRes * comp = (AliComparisonRes*)f.Get("AliComparisonRes");
17
18 // analyse comparison data (output stored in pictures_res.root)
19 comp->Analyse();
20
21 // paramtetrisation of the TPC track length (for information only)
22 TF1 fl("fl","((min(250./(abs(x+0.000001)),250)-90))",0,2); // TPC track length function
23 TF1 fl2("fl2","[0]/((min(250./(abs(x+0.000001)),250)-90))^[1]",0,2);
24 fl2.SetParameter(1,1);
25 fl2.SetParameter(0,1);
26*/
27
28#include <iostream>
29
30#include "TFile.h"
31#include "TCint.h"
32#include "TH3F.h"
33#include "TH2F.h"
34#include "TF1.h"
35#include "TProfile.h"
36#include "TProfile2D.h"
37#include "TGraph2D.h"
38#include "TCanvas.h"
39#include "TGraph.h"
40
41#include "AliESDEvent.h"
42#include "AliESD.h"
43#include "AliESDfriend.h"
44#include "AliESDfriendTrack.h"
45#include "AliESDVertex.h"
46#include "AliRecInfoCuts.h"
47#include "AliMCInfoCuts.h"
48#include "AliLog.h"
49#include "AliTracker.h"
50
51#include "AliMathBase.h"
52#include "AliTreeDraw.h"
53
54#include "AliMCInfo.h"
55#include "AliESDRecInfo.h"
56#include "AliComparisonRes.h"
57
58using namespace std;
59
60ClassImp(AliComparisonRes)
61
62//_____________________________________________________________________________
63AliComparisonRes::AliComparisonRes():
64 TNamed("AliComparisonRes","AliComparisonRes"),
65
66 // Resolution
67 fPtResolLPT(0), // pt resolution - low pt
68 fPtResolHPT(0), // pt resolution - high pt
69 fPtPullLPT(0), // pt resolution - low pt
70 fPtPullHPT(0), // pt resolution - high pt
71 //
72 // Resolution constrained param
73 //
74 fCPhiResolTan(0), // angular resolution - constrained
75 fCTanResolTan(0), // angular resolution - constrained
76 fCPtResolTan(0), // pt resolution - constrained
77 fCPhiPullTan(0), // angular resolution - constrained
78 fCTanPullTan(0), // angular resolution - constrained
79 fCPtPullTan(0), // pt resolution - constrained
80
81 //
82 // Parametrisation histograms
83 //
84
85 f1Pt2Resol1PtTPC(0),
86 f1Pt2Resol1PtTPCITS(0),
87 fYResol1PtTPC(0),
88 fYResol1PtTPCITS(0),
89 fZResol1PtTPC(0),
90 fZResol1PtTPCITS(0),
91 fPhiResol1PtTPC(0),
92 fPhiResol1PtTPCITS(0),
93 fThetaResol1PtTPC(0),
94 fThetaResol1PtTPCITS(0),
95
96 // constrained
97 fC1Pt2Resol1PtTPC(0),
98 fC1Pt2Resol1PtTPCITS(0),
99 fCYResol1PtTPC(0),
100 fCYResol1PtTPCITS(0),
101 fCZResol1PtTPC(0),
102 fCZResol1PtTPCITS(0),
103 fCPhiResol1PtTPC(0),
104 fCPhiResol1PtTPCITS(0),
105 fCThetaResol1PtTPC(0),
106 fCThetaResol1PtTPCITS(0),
107
108 // vertex
109 fVertex(0),
110
111 // Cuts
112 fCutsRC(0),
113 fCutsMC(0)
114{
115 InitHisto();
116 InitCuts();
117
118 // vertex (0,0,0)
119 fVertex = new AliESDVertex();
120 fVertex->SetXv(0.0);
121 fVertex->SetYv(0.0);
122 fVertex->SetZv(0.0);
123}
124
125//_____________________________________________________________________________
126AliComparisonRes::~AliComparisonRes(){
127
128 // Resolution histograms
129 if(fPtResolLPT) delete fPtResolLPT; fPtResolLPT=0;
130 if(fPtResolHPT) delete fPtResolHPT; fPtResolHPT=0;
131 if(fPtPullLPT) delete fPtPullLPT; fPtPullLPT=0;
132 if(fPtPullHPT) delete fPtPullHPT; fPtPullHPT=0;
133
134 // Resolution histograms (constrained param)
135 if(fCPhiResolTan) delete fCPhiResolTan; fCPhiResolTan=0;
136 if(fCTanResolTan) delete fCTanResolTan; fCTanResolTan=0;
137 if(fCPtResolTan) delete fCPtResolTan; fCPtResolTan=0;
138 if(fCPhiPullTan) delete fCPhiPullTan; fCPhiPullTan=0;
139 if(fCTanPullTan) delete fCTanPullTan; fCTanPullTan=0;
140 if(fCPtPullTan) delete fCPtPullTan; fCPtPullTan=0;
141
142 // Parametrisation histograms
143 //
144 if(f1Pt2Resol1PtTPC) delete f1Pt2Resol1PtTPC; f1Pt2Resol1PtTPC=0;
145 if(f1Pt2Resol1PtTPCITS) delete f1Pt2Resol1PtTPCITS; f1Pt2Resol1PtTPCITS=0;
146 if(fYResol1PtTPC) delete fYResol1PtTPC; fYResol1PtTPC=0;
147 if(fYResol1PtTPCITS) delete fYResol1PtTPCITS; fYResol1PtTPCITS=0;
148 if(fZResol1PtTPC) delete fZResol1PtTPC; fZResol1PtTPC=0;
149 if(fZResol1PtTPCITS) delete fZResol1PtTPCITS; fZResol1PtTPCITS=0;
150 if(fPhiResol1PtTPC) delete fPhiResol1PtTPC; fPhiResol1PtTPC=0;
151 if(fPhiResol1PtTPCITS) delete fPhiResol1PtTPCITS; fPhiResol1PtTPCITS=0;
152 if(fThetaResol1PtTPC) delete fThetaResol1PtTPC; fThetaResol1PtTPC=0;
153 if(fThetaResol1PtTPCITS) delete fThetaResol1PtTPCITS; fThetaResol1PtTPCITS=0;
154
155 // constrained
156 if(fC1Pt2Resol1PtTPC) delete fC1Pt2Resol1PtTPC; fC1Pt2Resol1PtTPC=0;
157 if(fC1Pt2Resol1PtTPCITS) delete fC1Pt2Resol1PtTPCITS; fC1Pt2Resol1PtTPCITS=0;
158 if(fCYResol1PtTPC) delete fCYResol1PtTPC; fCYResol1PtTPC=0;
159 if(fCYResol1PtTPCITS) delete fCYResol1PtTPCITS; fCYResol1PtTPCITS=0;
160 if(fCZResol1PtTPC) delete fCZResol1PtTPC; fCZResol1PtTPC=0;
161 if(fCZResol1PtTPCITS) delete fCZResol1PtTPCITS; fCZResol1PtTPCITS=0;
162 if(fCPhiResol1PtTPC) delete fCPhiResol1PtTPC; fCPhiResol1PtTPC=0;
163 if(fCPhiResol1PtTPCITS) delete fCPhiResol1PtTPCITS; fCPhiResol1PtTPCITS=0;
164 if(fCThetaResol1PtTPC) delete fCThetaResol1PtTPC; fCThetaResol1PtTPC=0;
165 if(fCThetaResol1PtTPCITS) delete fCThetaResol1PtTPCITS; fCThetaResol1PtTPCITS=0;
166
167 if(fVertex) delete fVertex; fVertex=0;
168
169}
170
171//_____________________________________________________________________________
172void AliComparisonRes::InitHisto(){
173
174 // Init histograms
175 fCPhiResolTan = new TH2F("CPhiResolTan","CPhiResolTan",50, -2,2,200,-0.025,0.025);
176 fCPhiResolTan->SetXTitle("tan(#theta)");
177 fCPhiResolTan->SetYTitle("#Delta#phi");
178
179 fCTanResolTan = new TH2F("CTanResolTan","CTanResolTan",50, -2,2,200,-0.025,0.025);
180 fCTanResolTan->SetXTitle("tan(#theta)");
181 fCTanResolTan->SetYTitle("#Delta#theta");
182
183 fCPtResolTan=new TH2F("CPtResol","CPtResol",50, -2,2,200,-0.2,0.2);
184 fCPtResolTan->SetXTitle("Tan(#theta)");
185 fCPtResolTan->SetYTitle("#Deltap_{t}/p_{t}");
186
187 fCPhiPullTan = new TH2F("CPhiPullTan","CPhiPullTan",50, -2,2,200,-5,5);
188 fCPhiPullTan->SetXTitle("Tan(#theta)");
189 fCPhiPullTan->SetYTitle("#Delta#phi/#Sigma");
190
191 fCTanPullTan = new TH2F("CTanPullTan","CTanPullTan",50, -2,2,200,-5,5);
192 fCTanPullTan->SetXTitle("Tan(#theta)");
193 fCTanPullTan->SetYTitle("#Delta#theta/#Sigma");
194
195 fCPtPullTan=new TH2F("CPtPull","CPtPull",50, -2,2,200,-5,5);
196 fCPtPullTan->SetXTitle("Tan(#theta)");
197 fCPtPullTan->SetYTitle("(1/mcp_{t}-1/p_{t})/#Sigma");
198
199 fPtResolLPT = new TH2F("Pt_resol_lpt","pt resol",10, 0.1,3,200,-0.2,0.2);
200 fPtResolLPT->SetXTitle("p_{t}");
201 fPtResolLPT->SetYTitle("#Deltap_{t}/p_{t}");
202
203 fPtResolHPT = new TH2F("Pt_resol_hpt","pt resol",10, 2,100,200,-0.3,0.3);
204 fPtResolHPT->SetXTitle("p_{t}");
205 fPtResolHPT->SetYTitle("#Deltap_{t}/p_{t}");
206
207 fPtPullLPT = new TH2F("Pt_pull_lpt","pt pull",10, 0.1,3,200,-6,6);
208 fPtPullLPT->SetXTitle("p_{t}");
209 fPtPullLPT->SetYTitle("#Deltap_{t}/#Sigma");
210
211 fPtPullHPT = new TH2F("Pt_pull_hpt","pt pull",10,2,100,200,-6,6);
212 fPtPullHPT->SetXTitle("p_{t}");
213 fPtPullHPT->SetYTitle("#Deltap_{t}/#Sigma");
214
215 //
216 // Parametrisation histograms
217 //
218
219 f1Pt2Resol1PtTPC = new TH2F("f1Pt2Resol1PtTPC","(1/mcpt-1/pt)/(1+1/mcpt)^2 vs 1/pt)",100,0,10,200,-0.010,0.010);
220 f1Pt2Resol1PtTPC->SetXTitle("1/mcp_{t}");
221 f1Pt2Resol1PtTPC->SetYTitle("(1/mcp_{t}-1/p_{t})/(1+1/mcp_{t})^2)");
222
223 f1Pt2Resol1PtTPCITS = new TH2F("f1Pt2Resol1PtTPCITS","(1/mcpt-1/pt)/(1+1/mcpt)^2 vs 1/pt)",100,0,10,200,-0.010,0.010);
224 f1Pt2Resol1PtTPCITS->SetXTitle("1/mcp_{t}");
225 f1Pt2Resol1PtTPCITS->SetYTitle("(1/mcp_{t}-1/p_{t})/(1+1/mcp_{t})^2)");
226
227 fYResol1PtTPC = new TH2F("fYResol1PtTPC","fYResol1PtTPC",100, 0,10,200,-1.0,1.0);
228 fYResol1PtTPC->SetXTitle("1/mcpt");
229 fYResol1PtTPC->SetYTitle("#DeltaY");
230
231 fYResol1PtTPCITS = new TH2F("fYResol1PtTPCITS","fYResol1PtTPCITS",100, 0,10,200,-0.05,0.05);
232 fYResol1PtTPCITS->SetXTitle("1/mcpt");
233 fYResol1PtTPCITS->SetYTitle("#DeltaY");
234
235 fZResol1PtTPC = new TH2F("fZResol1PtTPC","fZResol1PtTPC",100, 0,10,200,-1.0,1.0);
236 fZResol1PtTPC->SetXTitle("1/mcpt");
237 fZResol1PtTPC->SetYTitle("#DeltaZ");
238
239 fZResol1PtTPCITS = new TH2F("fZResol1PtTPCITS","fZResol1PtTPCITS",100, 0,10,200,-0.05,0.05);
240 fZResol1PtTPCITS->SetXTitle("1/mcpt");
241 fZResol1PtTPCITS->SetYTitle("#DeltaZ");
242
243 fPhiResol1PtTPC = new TH2F("fPhiResol1PtTPC","fPhiResol1PtTPC",100, 0,10,200,-0.025,0.025);
244 fPhiResol1PtTPC->SetXTitle("1/mcpt");
245 fPhiResol1PtTPC->SetYTitle("#Delta#phi");
246
247 fPhiResol1PtTPCITS = new TH2F("fPhiResol1PtTPCITS","fPhiResol1PtTPCITS",100, 0,10,200,-0.01,0.01);
248 fPhiResol1PtTPCITS->SetXTitle("1/mcpt");
249 fPhiResol1PtTPCITS->SetYTitle("#Delta#phi");
250
251 fThetaResol1PtTPC = new TH2F("fThetaResol1PtTPC","fThetaResol1PtTPC",100, 0,10,200,-0.025,0.025);
252 fThetaResol1PtTPC->SetXTitle("1/mcpt");
253 fThetaResol1PtTPC->SetYTitle("#Delta#theta");
254
255 fThetaResol1PtTPCITS = new TH2F("fThetaResol1PtTPCITS","fThetaResol1PtTPCITS",100, 0,10,200,-0.01,0.01);
256 fThetaResol1PtTPCITS->SetXTitle("1/mcpt");
257 fThetaResol1PtTPCITS->SetYTitle("#Delta#theta");
258
259 // constrained
260 fC1Pt2Resol1PtTPC = new TH2F("fC1Pt2Resol1PtTPC","(1/mcpt-1/pt)/(1+1/mcpt)^2 vs 1/pt)",100,0,10,200,-0.010,0.010);
261 fC1Pt2Resol1PtTPC->SetXTitle("1/mcp_{t}");
262 fC1Pt2Resol1PtTPC->SetYTitle("(1/mcp_{t}-1/p_{t})/(1+1/mcp_{t})^2)");
263
264 fC1Pt2Resol1PtTPCITS = new TH2F("fC1Pt2Resol1PtTPCITS","(1/mcpt-1/pt)/(1+1/mcpt)^2 vs 1/pt)",100,0,10,200,-0.010,0.010);
265 fC1Pt2Resol1PtTPCITS->SetXTitle("1/mcp_{t}");
266 fC1Pt2Resol1PtTPCITS->SetYTitle("(1/mcp_{t}-1/p_{t})/(1+1/mcp_{t})^2)");
267
268 fCYResol1PtTPC = new TH2F("fCYResol1PtTPC","fCYResol1PtTPC",100, 0,10,200,-1.0,1.0);
269 fCYResol1PtTPC->SetXTitle("1/mcpt");
270 fCYResol1PtTPC->SetYTitle("#DeltaY");
271
272 fCYResol1PtTPCITS = new TH2F("fCYResol1PtTPCITS","fCYResol1PtTPCITS",100, 0,10,200,-0.05,0.05);
273 fCYResol1PtTPCITS->SetXTitle("1/mcpt");
274 fCYResol1PtTPCITS->SetYTitle("#DeltaY");
275
276 fCZResol1PtTPC = new TH2F("fCZResol1PtTPC","fCZResol1PtTPC",100, 0,10,200,-1.0,1.0);
277 fCZResol1PtTPC->SetXTitle("1/mcpt");
278 fCZResol1PtTPC->SetYTitle("#DeltaZ");
279
280 fCZResol1PtTPCITS = new TH2F("fCZResol1PtTPCITS","fCZResol1PtTPCITS",100, 0,10,200,-0.05,0.05);
281 fCZResol1PtTPCITS->SetXTitle("1/mcpt");
282 fCZResol1PtTPCITS->SetYTitle("#DeltaZ");
283
284 fCPhiResol1PtTPC = new TH2F("fCPhiResol1PtTPC","fCPhiResol1PtTPC",100, 0,10,200,-0.025,0.025);
285 fCPhiResol1PtTPC->SetXTitle("1/mcpt");
286 fCPhiResol1PtTPC->SetYTitle("#Delta#phi");
287
288 fCPhiResol1PtTPCITS = new TH2F("fCPhiResol1PtTPCITS","fCPhiResol1PtTPCITS",100, 0,10,200,-0.01,0.01);
289 fCPhiResol1PtTPCITS->SetXTitle("1/mcpt");
290 fCPhiResol1PtTPCITS->SetYTitle("#Delta#phi");
291
292 fCThetaResol1PtTPC = new TH2F("fCThetaResol1PtTPC","fCThetaResol1PtTPC",100, 0,10,200,-0.025,0.025);
293 fCThetaResol1PtTPC->SetXTitle("1/mcpt");
294 fCThetaResol1PtTPC->SetYTitle("#Delta#theta");
295
296 fCThetaResol1PtTPCITS = new TH2F("fCThetaResol1PtTPCITS","fCThetaResol1PtTPCITS",100, 0,10,200,-0.01,0.01);
297 fCThetaResol1PtTPCITS->SetXTitle("1/mcpt");
298 fCThetaResol1PtTPCITS->SetYTitle("#Delta#theta");
299}
300
301//_____________________________________________________________________________
302void AliComparisonRes::InitCuts()
303{
304 // Init cuts
305 if(!fCutsMC)
306 AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
307 if(!fCutsRC)
308 AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
309}
310
311//_____________________________________________________________________________
312void AliComparisonRes::Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC)
313{
314 // Fill resolution comparison information
315 AliExternalTrackParam *track = 0;
316 Double_t kRadius = 3.0; // beam pipe radius
317 Double_t kMaxStep = 5.0; // max step
318 Double_t field = AliTracker::GetBz(); // nominal Bz field [kG]
319 Double_t kMaxD = 123456.0; // max distance
320
321 Int_t clusterITS[200];
322 Double_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
323
324 Float_t deltaPt, pullPt, deltaPhi, deltaTan, delta1Pt2, deltaY1Pt, deltaZ1Pt, deltaPhi1Pt, deltaTheta1Pt;
325 Float_t deltaPtTPC, pullPtTPC, deltaPhiTPC, deltaTanTPC, delta1Pt2TPC, deltaY1PtTPC, deltaZ1PtTPC, deltaPhi1PtTPC, deltaTheta1PtTPC;
326
327 Float_t mcpt = infoMC->GetParticle().Pt();
328
329 // distance to Prim. vertex
330 const Double_t* dv = infoMC->GetVDist();
331 Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz();
332
333 // Check selection cuts
334 if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return;
335 if (!isPrim) return;
336 //if (infoRC->GetStatus(1)==0) return;
337 if (infoRC->GetStatus(1)!=3) return; // TPC refit
338 if (!infoRC->GetESDtrack()) return;
339 if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
340
341 // calculate and set prim. vertex
342 fVertex->SetXv( infoMC->GetParticle().Vx() - dv[0] );
343 fVertex->SetYv( infoMC->GetParticle().Vy() - dv[1] );
344 fVertex->SetZv( infoMC->GetParticle().Vz() - dv[2] );
345
346 deltaPt= (mcpt-infoRC->GetESDtrack()->Pt())/mcpt;
347 pullPt= (1/mcpt-infoRC->GetESDtrack()->OneOverPt())/TMath::Sqrt(infoRC->GetESDtrack()->GetSigma1Pt2());
348 deltaPhi = TMath::ATan2(infoRC->GetESDtrack()->Py(),infoRC->GetESDtrack()->Px())-
349 TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px());
350
351 deltaTan = TMath::ATan2(infoRC->GetESDtrack()->Pz(),infoRC->GetESDtrack()->Pt())-
352 TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt());
353
354 delta1Pt2 = (1/mcpt-infoRC->GetESDtrack()->OneOverPt())/TMath::Power(1+1/mcpt,2);
355 deltaY1Pt = (infoMC->GetParticle().Vy()-infoRC->GetESDtrack()->GetY()) / (0.2+1/mcpt);
356 deltaZ1Pt = (infoMC->GetParticle().Vz()-infoRC->GetESDtrack()->GetZ()) / (0.2+1/mcpt);
357 deltaPhi1Pt = deltaPhi / (0.1+1/mcpt);
358 deltaTheta1Pt = deltaTan / (0.1+1/mcpt);
359
360 //track parameters at the first measured point (TPC)
361 //Double_t param[5],x,alpha; // [0]-Y [cm],[1]-Z [cm],[2]-sin(phi),[3]-tan(theta),[4]-1/pt [1/GeV]
362 //infoRC->GetESDtrack()->GetInnerExternalParameters(alpha,x,param);
363 //const AliExternalTrackParam *innerTPC = infoRC->GetESDtrack()->GetInnerParam();
364 //const AliExternalTrackParam *innerTPC = infoRC->GetESDtrack()->GetTPCInnerParam();
365
366
367 // calculate track parameters at vertex
368 const AliExternalTrackParam *innerTPC = 0;
369 if ((innerTPC = infoRC->GetESDtrack()->GetTPCInnerParam()) != 0)
370 {
371 if ((track = new AliExternalTrackParam(*infoRC->GetESDtrack()->GetTPCInnerParam())) != 0 )
372 {
373 Bool_t bStatus = AliTracker::PropagateTrackTo(track,kRadius,infoMC->GetMass(),kMaxStep,kTRUE);
374 Bool_t bDCAStatus = track->PropagateToDCA(fVertex,field,kMaxD,dca,cov);
375
376 // Fill parametrisation histograms (only TPC track)
377 if(bStatus && bDCAStatus)
378 {
379 deltaPtTPC= (mcpt-innerTPC->Pt())/mcpt;
380 pullPtTPC= (1/mcpt-innerTPC->OneOverPt())/TMath::Sqrt(innerTPC->GetSigma1Pt2());
381 deltaPhiTPC = TMath::ATan2(innerTPC->Py(),innerTPC->Px())-
382 TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px());
383
384 deltaTanTPC = TMath::ATan2(innerTPC->Pz(),innerTPC->Pt())-
385 TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt());
386
387 delta1Pt2TPC = (1/mcpt-innerTPC->OneOverPt())/TMath::Power(1+1/mcpt,2);
388 deltaY1PtTPC= (infoMC->GetParticle().Vy()-innerTPC->GetY()) / (0.2+1/mcpt);
389 deltaZ1PtTPC = (infoMC->GetParticle().Vz()-innerTPC->GetZ()) / (0.2+1/mcpt);
390 deltaPhi1PtTPC = deltaPhiTPC / (0.1+1/mcpt);
391 deltaTheta1PtTPC = deltaTanTPC / (0.1+1/mcpt);
392
393 f1Pt2Resol1PtTPC->Fill(1/mcpt,delta1Pt2TPC);
394 fYResol1PtTPC->Fill(1/mcpt,deltaY1PtTPC);
395 fZResol1PtTPC->Fill(1/mcpt,deltaZ1PtTPC);
396 fPhiResol1PtTPC->Fill(1/mcpt,deltaPhi1PtTPC);
397 fThetaResol1PtTPC->Fill(1/mcpt,deltaTheta1PtTPC);
398 }
399 delete track;
400 }
401 }
402
403 // TPC and ITS (nb. of clusters >2) in the system
404 if(infoRC->GetESDtrack()->GetITSclusters(clusterITS)>2)
405 {
406 f1Pt2Resol1PtTPCITS->Fill(1/mcpt,delta1Pt2);
407 fYResol1PtTPCITS->Fill(1/mcpt,deltaY1Pt);
408 fZResol1PtTPCITS->Fill(1/mcpt,deltaZ1Pt);
409 fPhiResol1PtTPCITS->Fill(1/mcpt,deltaPhi1Pt);
410 fThetaResol1PtTPCITS->Fill(1/mcpt,deltaTheta1Pt);
411 }
412
413 // Fill histograms
414 fPtResolLPT->Fill(mcpt,deltaPt);
415 fPtResolHPT->Fill(mcpt,deltaPt);
416 fPtPullLPT->Fill(mcpt,pullPt);
417 fPtPullHPT->Fill(mcpt,pullPt);
418}
419
420//_____________________________________________________________________________
421void AliComparisonRes::ProcessConstrained(AliMCInfo* infoMC, AliESDRecInfo *infoRC)
422{
423 // Fill resolution comparison information (constarained parameters)
424 //
425 AliExternalTrackParam *track = 0;
426 Double_t kRadius = 3.0; // beam pipe radius
427 Double_t kMaxStep = 5.0; // max step
428 Double_t field = AliTracker::GetBz(); // nominal Bz field [kG]
429 Double_t kMaxD = 123456.0; // max distance
430
431 Int_t clusterITS[200];
432 Double_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
433
434 Float_t deltaPt, pullPt, deltaPhi, pullPhi, deltaTan, pullTan, delta1Pt2, deltaY1Pt, deltaZ1Pt, deltaPhi1Pt, deltaTheta1Pt;
435 Float_t deltaPtTPC, pullPtTPC, deltaPhiTPC, deltaTanTPC, delta1Pt2TPC, deltaY1PtTPC, deltaZ1PtTPC, deltaPhi1PtTPC, deltaTheta1PtTPC;
436
437 Float_t mcpt = infoMC->GetParticle().Pt();
438 Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);
439
440 // distance to Prim. vertex
441 const Double_t* dv = infoMC->GetVDist();
442 Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz();
443
444 // Check selection cuts
445 if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return;
446 if (!isPrim) return;
447 if (infoRC->GetStatus(1)!=3) return;
448 if (!infoRC->GetESDtrack()) return;
449 if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
450 if (!infoRC->GetESDtrack()->GetConstrainedParam()) return;
451
452// calculate and set prim. vertex
453 fVertex->SetXv( infoMC->GetParticle().Vx() - dv[0] );
454 fVertex->SetYv( infoMC->GetParticle().Vy() - dv[1] );
455 fVertex->SetZv( infoMC->GetParticle().Vz() - dv[2] );
456
457 // constrained parameters resolution
458 const AliExternalTrackParam * cparam = infoRC->GetESDtrack()->GetConstrainedParam();
459 deltaPt= (mcpt-cparam->Pt())/mcpt;
460 pullPt= (1/mcpt-cparam->OneOverPt())/TMath::Sqrt(cparam->GetSigma1Pt2());
461 deltaPhi = TMath::ATan2(cparam->Py(),cparam->Px())-
462 TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px());
463 pullPhi = deltaPhi/TMath::Sqrt(cparam->GetSigmaSnp2());
464 deltaTan = TMath::ATan2(cparam->Pz(),cparam->Pt())-TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt());
465 pullTan = deltaPhi/TMath::Sqrt(cparam->GetSigmaSnp2());
466
467
468 delta1Pt2 = (1/mcpt-cparam->OneOverPt())/TMath::Power(1+1/mcpt,2);
469
470 deltaY1Pt = (infoMC->GetParticle().Vy()-cparam->GetY()) / (0.2+1/mcpt);
471 deltaZ1Pt = (infoMC->GetParticle().Vz()-cparam->GetZ()) / (0.2+1/mcpt);
472 deltaPhi1Pt = deltaPhi / (0.1+1/mcpt);
473 deltaTheta1Pt = deltaTan / (0.1+1/mcpt);
474
475 // calculate track parameters at vertex
476 const AliExternalTrackParam *innerTPC = 0;
477 if ((innerTPC = infoRC->GetESDtrack()->GetTPCInnerParam()) != 0)
478 {
479 if ((track = new AliExternalTrackParam(*infoRC->GetESDtrack()->GetTPCInnerParam())) != 0 )
480 {
481 Bool_t bStatus = AliTracker::PropagateTrackTo(track,kRadius,infoMC->GetMass(),kMaxStep,kTRUE);
482 Bool_t bDCAStatus = track->PropagateToDCA(fVertex,field,kMaxD,dca,cov);
483
484 // Fill parametrisation histograms (only TPC track)
485 if(bStatus && bDCAStatus)
486 {
487 deltaPtTPC= (mcpt-innerTPC->Pt())/mcpt;
488 pullPtTPC= (1/mcpt-innerTPC->OneOverPt())/TMath::Sqrt(innerTPC->GetSigma1Pt2());
489 deltaPhiTPC = TMath::ATan2(innerTPC->Py(),innerTPC->Px())-
490 TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px());
491
492 deltaTanTPC = TMath::ATan2(innerTPC->Pz(),innerTPC->Pt())-
493 TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt());
494
495 delta1Pt2TPC = (1/mcpt-innerTPC->OneOverPt())/TMath::Power(1+1/mcpt,2);
496 deltaY1PtTPC= (infoMC->GetParticle().Vy()-innerTPC->GetY()) / (0.2+1/mcpt);
497 deltaZ1PtTPC = (infoMC->GetParticle().Vz()-innerTPC->GetZ()) / (0.2+1/mcpt);
498 deltaPhi1PtTPC = deltaPhiTPC / (0.1+1/mcpt);
499 deltaTheta1PtTPC = deltaTanTPC / (0.1+1/mcpt);
500
501 fC1Pt2Resol1PtTPC->Fill(1/mcpt,delta1Pt2TPC);
502 fCYResol1PtTPC->Fill(1/mcpt,deltaY1PtTPC);
503 fCZResol1PtTPC->Fill(1/mcpt,deltaZ1PtTPC);
504 fCPhiResol1PtTPC->Fill(1/mcpt,deltaPhi1PtTPC);
505 fCThetaResol1PtTPC->Fill(1/mcpt,deltaTheta1PtTPC);
506 }
507 delete track;
508 }
509 }
510
511 // TPC and ITS (nb. of clusters >2) in the system
512 if(infoRC->GetESDtrack()->GetITSclusters(clusterITS)>2)
513 {
514 fC1Pt2Resol1PtTPCITS->Fill(1/mcpt,delta1Pt2);
515 fCYResol1PtTPCITS->Fill(1/mcpt,deltaY1Pt);
516 fCZResol1PtTPCITS->Fill(1/mcpt,deltaZ1Pt);
517 fCPhiResol1PtTPCITS->Fill(1/mcpt,deltaPhi1Pt);
518 fCThetaResol1PtTPCITS->Fill(1/mcpt,deltaTheta1Pt);
519 }
520
521 // Fill histograms
522 fCPtResolTan->Fill(tantheta,deltaPt);
523 fCPtPullTan->Fill(tantheta,pullPt);
524 fCPhiResolTan->Fill(tantheta,deltaPhi);
525 fCPhiPullTan->Fill(tantheta,pullPhi);
526 fCTanResolTan->Fill(tantheta,deltaTan);
527 fCTanPullTan->Fill(tantheta,pullTan);
528}
529
530//_____________________________________________________________________________
531void AliComparisonRes::Exec(AliMCInfo* infoMC, AliESDRecInfo *infoRC){
532
533 // Process comparison information
534 Process(infoMC,infoRC);
535 ProcessConstrained(infoMC,infoRC);
536}
537
538//_____________________________________________________________________________
539TH1F* AliComparisonRes::MakeResol(TH2F * his, Int_t integ, Bool_t type){
540 // Create resolution histograms
541
542 TH1F *hisr, *hism;
543 if (!gPad) new TCanvas;
544 hisr = AliTreeDraw::CreateResHistoI(his,&hism,integ);
545 if (type) return hism;
546 else
547 return hisr;
548}
549
550//_____________________________________________________________________________
551void AliComparisonRes::Analyse(){
552 // Analyse comparison information and store output histograms
553 // in the "pictures_res.root" file
554
555 AliComparisonRes * comp=this;
556 TH1F *hiss=0;
557
558 TFile *fp = new TFile("pictures_res.root","recreate");
559 fp->cd();
560
561 TCanvas * c = new TCanvas("Phi resol Tan","Phi resol Tan");
562 c->cd();
563 //
564 hiss = comp->MakeResol(comp->fCPtResolTan,1,0);
565 hiss->SetXTitle("Tan(#theta)");
566 hiss->SetYTitle("#sigmap_{t}/p_{t}");
567 hiss->Draw();
568 hiss->Write("CptResolTan");
569 //
570 hiss = comp->MakeResol(comp->fCPhiResolTan,1,0);
571 hiss->SetXTitle("Tan(#theta)");
572 hiss->SetYTitle("#sigma#phi (rad)");
573 hiss->Draw();
574 hiss->Write("PhiResolTan");
575 //
576 hiss = comp->MakeResol(comp->fCTanResolTan,1,0);
577 hiss->SetXTitle("Tan(#theta)");
578 hiss->SetYTitle("#sigma#theta (rad)");
579 hiss->Draw();
580 hiss->Write("ThetaResolTan");
581 //
582 hiss = comp->MakeResol(comp->fCPtPullTan,1,0);
583 hiss->SetXTitle("Tan(#theta)");
584 hiss->SetYTitle("1/mcp_{t}-1/p_{t}/#Sigma(1/p_{t})");
585 hiss->Draw();
586 hiss->Write("CptPullTan");
587 //
588 hiss = comp->MakeResol(comp->fC1Pt2Resol1PtTPC,1,0);
589 hiss->SetXTitle("1/mcp_{t}");
590 hiss->SetYTitle("1/mcp_{t}-1/p_{t}/(1+1/p_{t})^2");
591 hiss->Draw();
592 hiss->Write("C1Pt2Resol1PtTPC");
593 fC1Pt2Resol1PtTPC->Write();
594
595 hiss = comp->MakeResol(comp->fC1Pt2Resol1PtTPCITS,1,0);
596 hiss->SetXTitle("1/mcp_{t}");
597 hiss->SetYTitle("1/mcp_{t}-1/p_{t}/(1+1/p_{t})^2");
598 hiss->Draw();
599 hiss->Write("C1Pt2Resol1PtTPCITS");
600 fC1Pt2Resol1PtTPCITS->Write();
601 //
602 hiss = comp->MakeResol(comp->fCYResol1PtTPC,1,0);
603 hiss->SetXTitle("1/mcp_{t}");
604 hiss->SetYTitle("(mcy-y)/(0.2+1/mcp_{t})");
605 hiss->Draw();
606 hiss->Write("CYResol1PtTPC");
607 fCYResol1PtTPC->Write();
608
609 hiss = comp->MakeResol(comp->fCYResol1PtTPCITS,1,0);
610 hiss->SetXTitle("1/mcp_{t}");
611 hiss->SetYTitle("(mcy-y)/(0.2+1/mcp_{t})");
612 hiss->Draw();
613 hiss->Write("CYResol1PtTPCITS");
614 fCYResol1PtTPCITS->Write();
615 //
616 hiss = comp->MakeResol(comp->fCZResol1PtTPC,1,0);
617 hiss->SetXTitle("1/mcp_{t}");
618 hiss->SetYTitle("(mcz-z)/(0.2+1/mcp_{t})");
619 hiss->Draw();
620 hiss->Write("CZResol1PtTPC");
621 fCZResol1PtTPC->Write();
622
623 hiss = comp->MakeResol(comp->fCZResol1PtTPCITS,1,0);
624 hiss->SetXTitle("1/mcp_{t}");
625 hiss->SetYTitle("(mcz-z)/(0.2+1/mcp_{t})");
626 hiss->Draw();
627 hiss->Write("CZResol1PtTPCITS");
628 fCZResol1PtTPCITS->Write();
629 //
630 hiss = comp->MakeResol(comp->fCPhiResol1PtTPC,1,0);
631 hiss->SetXTitle("1/mcp_{t}");
632 hiss->SetYTitle("(mc#phi-#phi)/(0.1+1/mcp_{t})");
633 hiss->Draw();
634 hiss->Write("CPhiResol1PtTPC");
635 fCPhiResol1PtTPC->Write();
636
637 hiss = comp->MakeResol(comp->fCPhiResol1PtTPCITS,1,0);
638 hiss->SetXTitle("1/mcp_{t}");
639 hiss->SetYTitle("(mc#phi-#phi)/(0.1+1/mcp_{t})");
640 hiss->Draw();
641 hiss->Write("CPhiResol1PtTPCITS");
642 fCPhiResol1PtTPCITS->Write();
643 //
644 hiss = comp->MakeResol(comp->fCThetaResol1PtTPC,1,0);
645 hiss->SetXTitle("1/mcp_{t}");
646 hiss->SetYTitle("(mc#theta-#theta)/(0.1+1/mcp_{t})");
647 hiss->Draw();
648 hiss->Write("CThetaResol1PtTPC");
649 fCThetaResol1PtTPC->Write();
650
651 hiss = comp->MakeResol(comp->fCThetaResol1PtTPCITS,1,0);
652 hiss->SetXTitle("1/mcp_{t}");
653 hiss->SetYTitle("(mc#theta-#theta)/(0.1+1/mcp_{t})");
654 hiss->Draw();
655 hiss->Write("CThetaResol1PtTPCITS");
656 fCThetaResol1PtTPCITS->Write();
657
658 //
659 hiss = comp->MakeResol(comp->f1Pt2Resol1PtTPC,1,0);
660 hiss->SetXTitle("1/mcp_{t}");
661 hiss->SetYTitle("1/mcp_{t}-1/p_{t}/(1+1/p_{t})^2");
662 hiss->Draw();
663 hiss->Write("OnePt2Resol1PtTPC");
664 f1Pt2Resol1PtTPC->Write();
665
666 hiss = comp->MakeResol(comp->f1Pt2Resol1PtTPCITS,1,0);
667 hiss->SetXTitle("1/mcp_{t}");
668 hiss->SetYTitle("1/mcp_{t}-1/p_{t}/(1+1/p_{t})^2");
669 hiss->Draw();
670 hiss->Write("OnePt2Resol1PtTPCITS");
671 f1Pt2Resol1PtTPCITS->Write();
672 //
673 hiss = comp->MakeResol(comp->fYResol1PtTPC,1,0);
674 hiss->SetXTitle("1/mcp_{t}");
675 hiss->SetYTitle("(mcy-y)/(0.2+1/mcp_{t})");
676 hiss->Draw();
677 hiss->Write("YResol1PtTPC");
678 fYResol1PtTPC->Write();
679
680 hiss = comp->MakeResol(comp->fYResol1PtTPCITS,1,0);
681 hiss->SetXTitle("1/mcp_{t}");
682 hiss->SetYTitle("(mcy-y)/(0.2+1/mcp_{t})");
683 hiss->Draw();
684 hiss->Write("YResol1PtTPCITS");
685 fYResol1PtTPCITS->Write();
686 //
687 hiss = comp->MakeResol(comp->fZResol1PtTPC,1,0);
688 hiss->SetXTitle("1/mcp_{t}");
689 hiss->SetYTitle("(mcz-z)/(0.2+1/mcp_{t})");
690 hiss->Draw();
691 hiss->Write("ZResol1PtTPC");
692 fZResol1PtTPC->Write();
693
694 hiss = comp->MakeResol(comp->fZResol1PtTPCITS,1,0);
695 hiss->SetXTitle("1/mcp_{t}");
696 hiss->SetYTitle("(mcz-z)/(0.2+1/mcp_{t})");
697 hiss->Draw();
698 hiss->Write("ZResol1PtTPCITS");
699 fZResol1PtTPCITS->Write();
700 //
701 hiss = comp->MakeResol(comp->fPhiResol1PtTPC,1,0);
702 hiss->SetXTitle("1/mcp_{t}");
703 hiss->SetYTitle("(mc#phi-#phi)/(0.1+1/mcp_{t})");
704 hiss->Draw();
705 hiss->Write("PhiResol1PtTPC");
706 fPhiResol1PtTPC->Write();
707
708 hiss = comp->MakeResol(comp->fPhiResol1PtTPCITS,1,0);
709 hiss->SetXTitle("1/mcp_{t}");
710 hiss->SetYTitle("(mc#phi-#phi)/(0.1+1/mcp_{t})");
711 hiss->Draw();
712 hiss->Write("PhiResol1PtTPCITS");
713 fPhiResol1PtTPCITS->Write();
714 //
715 hiss = comp->MakeResol(comp->fThetaResol1PtTPC,1,0);
716 hiss->SetXTitle("1/mcp_{t}");
717 hiss->SetYTitle("(mc#theta-#theta)/(0.1+1/mcp_{t})");
718 hiss->Draw();
719 hiss->Write("ThetaResol1PtTPC");
720 fThetaResol1PtTPC->Write();
721
722 hiss = comp->MakeResol(comp->fThetaResol1PtTPCITS,1,0);
723 hiss->SetXTitle("1/mcp_{t}");
724 hiss->SetYTitle("(mc#theta-#theta)/(0.1+1/mcp_{t})");
725 hiss->Draw();
726 hiss->Write("ThetaResol1PtTPCITS");
727 fThetaResol1PtTPCITS->Write();
728
729 fp->Close();
730}
731
732//_____________________________________________________________________________
733Long64_t AliComparisonRes::Merge(TCollection* list)
734{
735 // Merge list of objects (needed by PROOF)
736
737 if (!list)
738 return 0;
739
740 if (list->IsEmpty())
741 return 1;
742
743 TIterator* iter = list->MakeIterator();
744 TObject* obj = 0;
745
746 // collection of generated histograms
747 Int_t count=0;
748 while((obj = iter->Next()) != 0)
749 {
750 AliComparisonRes* entry = dynamic_cast<AliComparisonRes*>(obj);
751 if (entry == 0) continue;
752
753 fPtResolLPT->Add(entry->fPtResolLPT);
754 fPtResolHPT->Add(entry->fPtResolHPT);
755 fPtPullLPT->Add(entry->fPtPullLPT);
756 fPtPullHPT->Add(entry->fPtPullHPT);
757
758 // Histograms for 1/pt parameterisation
759 f1Pt2Resol1PtTPC->Add(entry->f1Pt2Resol1PtTPC);
760 fYResol1PtTPC->Add(entry->fYResol1PtTPC);
761 fZResol1PtTPC->Add(entry->fZResol1PtTPC);
762 fPhiResol1PtTPC->Add(entry->fPhiResol1PtTPC);
763 fThetaResol1PtTPC->Add(entry->fThetaResol1PtTPC);
764
765 f1Pt2Resol1PtTPCITS->Add(entry->f1Pt2Resol1PtTPCITS);
766 fYResol1PtTPCITS->Add(entry->fYResol1PtTPCITS);
767 fZResol1PtTPCITS->Add(entry->fZResol1PtTPCITS);
768 fPhiResol1PtTPCITS->Add(entry->fPhiResol1PtTPCITS);
769 fThetaResol1PtTPCITS->Add(entry->fThetaResol1PtTPCITS);
770
771 // Resolution histograms (constrained param)
772 fCPhiResolTan->Add(entry->fCPhiResolTan);
773 fCTanResolTan->Add(entry->fCTanResolTan);
774 fCPtResolTan->Add(entry->fCPtResolTan);
775 fCPhiPullTan->Add(entry->fCPhiPullTan);
776 fCTanPullTan->Add(entry->fCTanPullTan);
777 fCPtPullTan->Add(entry->fCPtPullTan);
778
779 // Histograms for 1/pt parameterisation (constrained)
780 fC1Pt2Resol1PtTPC->Add(entry->fC1Pt2Resol1PtTPC);
781 fCYResol1PtTPC->Add(entry->fCYResol1PtTPC);
782 fCZResol1PtTPC->Add(entry->fCZResol1PtTPC);
783 fCPhiResol1PtTPC->Add(entry->fCPhiResol1PtTPC);
784 fCThetaResol1PtTPC->Add(entry->fCThetaResol1PtTPC);
785
786 fC1Pt2Resol1PtTPCITS->Add(entry->fC1Pt2Resol1PtTPCITS);
787 fCYResol1PtTPCITS->Add(entry->fCYResol1PtTPCITS);
788 fCZResol1PtTPCITS->Add(entry->fCZResol1PtTPCITS);
789 fCPhiResol1PtTPCITS->Add(entry->fCPhiResol1PtTPCITS);
790 fCThetaResol1PtTPCITS->Add(entry->fCThetaResol1PtTPCITS);
791
792 count++;
793 }
794
795return count;
796}