]>
Commit | Line | Data |
---|---|---|
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/graphs) are stored in the folder which is | |
8 | // a data member of AliComparisonRes. | |
9 | // | |
10 | // Author: J.Otwinowski 04/02/2008 | |
11 | //------------------------------------------------------------------------------ | |
12 | ||
13 | /* | |
14 | ||
15 | // after running comparison task, read the file, and get component | |
16 | gROOT->LoadMacro("$ALICE_ROOT/PWG1/Macros/LoadMyLibs.C"); | |
17 | LoadMyLibs(); | |
18 | ||
19 | TFile f("Output.root"); | |
20 | //AliComparisonRes * compObj = (AliComparisonRes*)f.Get("AliComparisonRes"); | |
21 | AliComparisonRes * compObj = (AliComparisonRes*)cOutput->FindObject("AliComparisonRes"); | |
22 | ||
23 | // analyse comparison data | |
24 | compObj->Analyse(); | |
25 | ||
26 | // the output histograms/graphs will be stored in the folder "folderRes" | |
27 | compObj->GetAnalysisFolder()->ls("*"); | |
28 | ||
29 | // user can save whole comparison object (or only folder with anlysed histograms) | |
30 | // in the seperate output file (e.g.) | |
31 | TFile fout("Analysed_Res.root","recreate"); | |
32 | compObj->Write(); // compObj->GetAnalysisFolder()->Write(); | |
33 | fout.Close(); | |
34 | ||
35 | */ | |
36 | ||
37 | #include "TCanvas.h" | |
38 | #include "TH1.h" | |
39 | #include "TAxis.h" | |
40 | ||
41 | #include "AliComparisonRes.h" | |
42 | #include "AliESDRecInfo.h" | |
43 | #include "AliESDVertex.h" | |
44 | #include "AliESDtrack.h" | |
45 | #include "AliLog.h" | |
46 | #include "AliMCInfo.h" | |
47 | #include "AliMCInfoCuts.h" | |
48 | #include "AliRecInfoCuts.h" | |
49 | #include "AliTracker.h" | |
50 | #include "AliTreeDraw.h" | |
51 | ||
52 | using namespace std; | |
53 | ||
54 | ClassImp(AliComparisonRes) | |
55 | ||
56 | //_____________________________________________________________________________ | |
57 | AliComparisonRes::AliComparisonRes(): | |
58 | AliComparisonObject("AliComparisonRes"), | |
59 | fResolHisto(0), | |
60 | fPullHisto(0), | |
61 | ||
62 | // Cuts | |
63 | fCutsRC(0), | |
64 | fCutsMC(0), | |
65 | ||
66 | // histogram folder | |
67 | fAnalysisFolder(0) | |
68 | { | |
69 | //Init(); | |
70 | } | |
71 | ||
72 | //_____________________________________________________________________________ | |
73 | AliComparisonRes::AliComparisonRes(Char_t* name="AliComparisonRes", Char_t* title="AliComparisonRes",Int_t analysisMode=0,Bool_t hptGenerator=kFALSE): | |
74 | //AliComparisonRes::AliComparisonRes(Char_t* name, Char_t* title,Int_t analysisMode,Bool_t hptGenerator): | |
75 | AliComparisonObject(name,title), | |
76 | fResolHisto(0), | |
77 | fPullHisto(0), | |
78 | ||
79 | // Cuts | |
80 | fCutsRC(0), | |
81 | fCutsMC(0), | |
82 | ||
83 | // histogram folder | |
84 | fAnalysisFolder(0) | |
85 | { | |
86 | // named constructor | |
87 | // | |
88 | SetAnalysisMode(analysisMode); | |
89 | SetHptGenerator(hptGenerator); | |
90 | ||
91 | Init(); | |
92 | } | |
93 | ||
94 | //_____________________________________________________________________________ | |
95 | AliComparisonRes::~AliComparisonRes() | |
96 | { | |
97 | // destructor | |
98 | ||
99 | if(fResolHisto) delete fResolHisto; fResolHisto=0; | |
100 | if(fPullHisto) delete fPullHisto; fPullHisto=0; | |
101 | ||
102 | if(fCutsRC) delete fCutsRC; fCutsRC=0; | |
103 | if(fCutsMC) delete fCutsMC; fCutsMC=0; | |
104 | ||
105 | if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0; | |
106 | } | |
107 | ||
108 | //_____________________________________________________________________________ | |
109 | void AliComparisonRes::Init(){ | |
110 | ||
111 | // | |
112 | // histogram bining | |
113 | // | |
114 | ||
115 | // set pt bins | |
116 | Int_t nPtBins = 31; | |
117 | Double_t binsPt[32] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.25,2.5,2.75,3.,3.5,4.,5.,6.,8.,10.}; | |
118 | Double_t ptMin = 0., ptMax = 10.; | |
119 | ||
120 | ||
121 | if(IsHptGenerator() == kTRUE) { | |
122 | nPtBins = 100; | |
123 | ptMin = 0.; ptMax = 100.; | |
124 | } | |
125 | ||
126 | // res_y:res_z:res_phi,res_lambda:res_1pt:y:z:eta:phi:pt | |
127 | Int_t binsResolHisto[10]={100,100,100,100,100,50,100,30,144,nPtBins}; | |
128 | Double_t minResolHisto[10]={-1.,-1.,-0.03,-0.03,-0.2,-0.03,-20.,-1.5,0.,ptMin}; | |
129 | Double_t maxResolHisto[10]={ 1., 1., 0.03, 0.03, 0.2, 0.03, 20., 1.5,2.*TMath::Pi(), ptMax}; | |
130 | ||
131 | ||
132 | fResolHisto = new THnSparseF("fResolHisto","res_y:res_z:res_phi:res_lambda:res_1pt:y:z:eta:phi:pt",10,binsResolHisto,minResolHisto,maxResolHisto); | |
133 | if(!IsHptGenerator()) fResolHisto->SetBinEdges(9,binsPt); | |
134 | ||
135 | fResolHisto->GetAxis(0)->SetTitle("res_y (cm)"); | |
136 | fResolHisto->GetAxis(1)->SetTitle("res_z (cm)"); | |
137 | fResolHisto->GetAxis(2)->SetTitle("res_phi (rad)"); | |
138 | fResolHisto->GetAxis(3)->SetTitle("res_lambda (rad)"); | |
139 | fResolHisto->GetAxis(4)->SetTitle("res_1pt (GeV/c)^{-1}"); | |
140 | fResolHisto->GetAxis(5)->SetTitle("y (cm)"); | |
141 | fResolHisto->GetAxis(6)->SetTitle("z (cm)"); | |
142 | fResolHisto->GetAxis(7)->SetTitle("eta"); | |
143 | fResolHisto->GetAxis(8)->SetTitle("phi (rad)"); | |
144 | fResolHisto->GetAxis(9)->SetTitle("pt (GeV/c)"); | |
145 | fResolHisto->Sumw2(); | |
146 | ||
147 | //pull_y:pull_z:pull_lambda:pull_phi:pull_1pt:y:z:eta:phi:pt | |
148 | Int_t binsPullHisto[10]={100,100,100,100,100,50,100,30,144,nPtBins}; | |
149 | Double_t minPullHisto[10]={-5.,-5.,-5.,-5.,-5.,-0.03, -20.,-1.5, 0., ptMin}; | |
150 | Double_t maxPullHisto[10]={ 5., 5., 5., 5., 5., 0.03, 20., 1.5, 2.*TMath::Pi(),ptMax}; | |
151 | ||
152 | fPullHisto = new THnSparseF("fPullHisto","pull_y:pull_z:pull_phi:pull_lambda:pull_1pt:y:z:eta:phi:pt",10,binsPullHisto,minPullHisto,maxPullHisto); | |
153 | if(!IsHptGenerator()) fPullHisto->SetBinEdges(9,binsPt); | |
154 | ||
155 | fPullHisto->GetAxis(0)->SetTitle("res_y (cm)"); | |
156 | fPullHisto->GetAxis(1)->SetTitle("res_z (cm)"); | |
157 | fPullHisto->GetAxis(2)->SetTitle("res_phi (rad)"); | |
158 | fPullHisto->GetAxis(3)->SetTitle("res_lambda (rad)"); | |
159 | fPullHisto->GetAxis(4)->SetTitle("res_1pt (GeV/c)^{-1}"); | |
160 | fPullHisto->GetAxis(5)->SetTitle("y (rad)"); | |
161 | fPullHisto->GetAxis(6)->SetTitle("z (rad)"); | |
162 | fPullHisto->GetAxis(7)->SetTitle("eta"); | |
163 | fPullHisto->GetAxis(8)->SetTitle("phi (rad)"); | |
164 | fPullHisto->GetAxis(9)->SetTitle("pt (GeV/c)"); | |
165 | fPullHisto->Sumw2(); | |
166 | ||
167 | // Init cuts | |
168 | if(!fCutsMC) | |
169 | AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object"); | |
170 | if(!fCutsRC) | |
171 | AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object"); | |
172 | ||
173 | // init folder | |
174 | fAnalysisFolder = CreateFolder("folderRes","Analysis Resolution Folder"); | |
175 | } | |
176 | ||
177 | //_____________________________________________________________________________ | |
178 | void AliComparisonRes::ProcessTPC(AliMCInfo* const infoMC, AliESDRecInfo *const infoRC) | |
179 | { | |
180 | // Fill resolution comparison information | |
181 | AliExternalTrackParam *track = 0; | |
182 | //Double_t field = AliTracker::GetBz(); // nominal Bz field [kG] | |
183 | Double_t kMaxD = 123456.0; // max distance | |
184 | AliESDVertex vertexMC; // MC primary vertex | |
185 | ||
186 | Double_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z | |
187 | ||
188 | Float_t delta1PtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC; | |
189 | Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC; | |
190 | //Float_t delta1Pt2TPC, deltaY1PtTPC, deltaZ1PtTPC, deltaPhi1PtTPC, deltaTheta1PtTPC; | |
191 | ||
192 | Float_t mceta = infoMC->GetParticle().Eta(); | |
193 | Float_t mcphi = infoMC->GetParticle().Phi(); | |
194 | if(mcphi<0) mcphi += 2.*TMath::Pi(); | |
195 | Float_t mcpt = infoMC->GetParticle().Pt(); | |
196 | ||
197 | // distance to Prim. vertex | |
198 | const Double_t* dv = infoMC->GetVDist(); | |
199 | Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz(); | |
200 | ||
201 | // Only 5 charged particle species (e,mu,pi,K,p) | |
202 | if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return; | |
203 | if (!isPrim) return; | |
204 | //if (infoRC->GetStatus(1)!=3) return; // TPC refit | |
205 | if (!infoRC->GetESDtrack()) return; | |
206 | if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; | |
207 | ||
208 | // calculate and set prim. vertex | |
209 | vertexMC.SetXv( infoMC->GetParticle().Vx() - dv[0] ); | |
210 | vertexMC.SetYv( infoMC->GetParticle().Vy() - dv[1] ); | |
211 | vertexMC.SetZv( infoMC->GetParticle().Vz() - dv[2] ); | |
212 | ||
213 | // calculate track parameters at vertex | |
214 | const AliExternalTrackParam *innerTPC = 0; | |
215 | if ((innerTPC = infoRC->GetESDtrack()->GetTPCInnerParam()) != 0) | |
216 | { | |
217 | if ((track = new AliExternalTrackParam(*infoRC->GetESDtrack()->GetTPCInnerParam())) != 0 ) | |
218 | { | |
219 | //Bool_t bDCAStatus = track->PropagateToDCA(&vertexMC,field,kMaxD,dca,cov); | |
220 | Double_t field[3]; track->GetBxByBz(field); | |
221 | Bool_t bDCAStatus = track->PropagateToDCABxByBz(&vertexMC,field,kMaxD,dca,cov); | |
222 | ||
223 | // Fill parametrisation histograms (only TPC track) | |
224 | if(bDCAStatus) | |
225 | { | |
226 | // select primaries | |
227 | if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ()) | |
228 | { | |
229 | ||
230 | deltaYTPC= track->GetY()-infoMC->GetParticle().Vy(); | |
231 | deltaZTPC = track->GetZ()-infoMC->GetParticle().Vz(); | |
232 | deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt()); | |
233 | deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px()); | |
234 | delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt; | |
235 | ||
236 | pullYTPC= (track->GetY()-infoMC->GetParticle().Vy()) / TMath::Sqrt(track->GetSigmaY2()); | |
237 | pullZTPC = (track->GetZ()-infoMC->GetParticle().Vz()) / TMath::Sqrt(track->GetSigmaZ2()); | |
238 | pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2()); | |
239 | pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2()); | |
240 | pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2()); | |
241 | ||
242 | Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt}; | |
243 | fResolHisto->Fill(vResolHisto); | |
244 | ||
245 | Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt}; | |
246 | fPullHisto->Fill(vPullHisto); | |
247 | } | |
248 | ||
249 | /* | |
250 | delta1Pt2TPC = (1/mcpt-innerTPC->OneOverPt())/TMath::Power(1+1/mcpt,2); | |
251 | deltaY1PtTPC= (infoMC->GetParticle().Vy()-innerTPC->GetY()) / (0.2+1/mcpt); | |
252 | deltaZ1PtTPC = (infoMC->GetParticle().Vz()-innerTPC->GetZ()) / (0.2+1/mcpt); | |
253 | deltaPhi1PtTPC = deltaPhiTPC / (0.1+1/mcpt); | |
254 | deltaTheta1PtTPC = deltaLambdaTPC / (0.1+1/mcpt); | |
255 | ||
256 | fPhiEtaPtTPC->Fill(TMath::ATan2(innerTPC->Py(),innerTPC->Px()), innerTPC->Eta(), innerTPC->Pt()); | |
257 | ||
258 | f1Pt2ResolS1PtTPC->Fill(s1mcpt,delta1Pt2TPC); | |
259 | fYResolS1PtTPC->Fill(s1mcpt,deltaY1PtTPC); | |
260 | fZResolS1PtTPC->Fill(s1mcpt,deltaZ1PtTPC); | |
261 | fPhiResolS1PtTPC->Fill(s1mcpt,deltaPhi1PtTPC); | |
262 | fThetaResolS1PtTPC->Fill(s1mcpt,deltaTheta1PtTPC); | |
263 | ||
264 | fPtResolPtTPC->Fill(mcpt,deltaPtTPC); | |
265 | fPtPullPtTPC->Fill(mcpt,pullPtTPC); | |
266 | fPhiResolEtaTPC->Fill(eta,deltaPhiTPC); | |
267 | fPhiPullEtaTPC->Fill(eta,pullPhiTPC); | |
268 | fLambdaResolEtaTPC->Fill(eta,deltaLambdaTPC); | |
269 | fLambdaPullEtaTPC->Fill(eta,pullLambdaTPC); | |
270 | */ | |
271 | } | |
272 | delete track; | |
273 | } | |
274 | } | |
275 | } | |
276 | ||
277 | //_____________________________________________________________________________ | |
278 | void AliComparisonRes::ProcessTPCITS(AliMCInfo* const infoMC, AliESDRecInfo *const infoRC) | |
279 | { | |
280 | Int_t clusterITS[200]; | |
281 | Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z | |
282 | ||
283 | Float_t delta1PtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC; | |
284 | Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC; | |
285 | ||
286 | Float_t mceta = infoMC->GetParticle().Eta(); | |
287 | Float_t mcphi = infoMC->GetParticle().Phi(); | |
288 | if(mcphi<0) mcphi += 2.*TMath::Pi(); | |
289 | Float_t mcpt = infoMC->GetParticle().Pt(); | |
290 | ||
291 | // distance to Prim. vertex | |
292 | const Double_t* dv = infoMC->GetVDist(); | |
293 | Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz(); | |
294 | ||
295 | // Only 5 charged particle species (e,mu,pi,K,p) | |
296 | if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return; | |
297 | if (!isPrim) return; | |
298 | if (infoRC->GetStatus(1)!=3) return; // TPC refit | |
299 | if (!infoRC->GetESDtrack()) return; | |
300 | if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; | |
301 | ||
302 | infoRC->GetESDtrack()->GetImpactParameters(dca,cov); | |
303 | ||
304 | // select primaries | |
305 | if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ()) | |
306 | { | |
307 | deltaYTPC= infoRC->GetESDtrack()->GetY()-infoMC->GetParticle().Vy(); | |
308 | deltaZTPC = infoRC->GetESDtrack()->GetZ()-infoMC->GetParticle().Vz(); | |
309 | deltaLambdaTPC = TMath::ATan2(infoRC->GetESDtrack()->Pz(),infoRC->GetESDtrack()->Pt())-TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt()); | |
310 | deltaPhiTPC = TMath::ATan2(infoRC->GetESDtrack()->Py(),infoRC->GetESDtrack()->Px())-TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px()); | |
311 | delta1PtTPC = (infoRC->GetESDtrack()->OneOverPt()-1./mcpt)*mcpt; | |
312 | ||
313 | pullYTPC= (infoRC->GetESDtrack()->GetY()-infoMC->GetParticle().Vy()) / TMath::Sqrt(infoRC->GetESDtrack()->GetSigmaY2()); | |
314 | pullZTPC = (infoRC->GetESDtrack()->GetZ()-infoMC->GetParticle().Vz()) / TMath::Sqrt(infoRC->GetESDtrack()->GetSigmaZ2()); | |
315 | pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(infoRC->GetESDtrack()->GetSigmaTgl2()); | |
316 | pullPhiTPC = deltaPhiTPC / TMath::Sqrt(infoRC->GetESDtrack()->GetSigmaSnp2()); | |
317 | pull1PtTPC = (infoRC->GetESDtrack()->OneOverPt()-1./mcpt) / TMath::Sqrt(infoRC->GetESDtrack()->GetSigma1Pt2()); | |
318 | ||
319 | Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt}; | |
320 | Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt}; | |
321 | ||
322 | // TPC and ITS clusters in the system | |
323 | if(infoRC->GetESDtrack()->GetITSclusters(clusterITS)>fCutsRC->GetMinNClustersITS()) | |
324 | { | |
325 | fResolHisto->Fill(vResolHisto); | |
326 | fPullHisto->Fill(vPullHisto); | |
327 | } | |
328 | } | |
329 | } | |
330 | ||
331 | //_____________________________________________________________________________ | |
332 | void AliComparisonRes::ProcessConstrained(AliMCInfo* const infoMC, AliESDRecInfo *const infoRC) | |
333 | { | |
334 | // Fill resolution comparison information (constarained parameters) | |
335 | // | |
336 | AliExternalTrackParam *track = 0; | |
337 | //Double_t field = AliTracker::GetBz(); // nominal Bz field [kG] | |
338 | Double_t kMaxD = 123456.0; // max distance | |
339 | AliESDVertex vertexMC; // MC primary vertex | |
340 | ||
341 | Double_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z | |
342 | ||
343 | Float_t delta1PtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC; | |
344 | Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC; | |
345 | ||
346 | Float_t mceta = infoMC->GetParticle().Eta(); | |
347 | Float_t mcphi = infoMC->GetParticle().Phi(); | |
348 | if(mcphi<0) mcphi += 2.*TMath::Pi(); | |
349 | Float_t mcpt = infoMC->GetParticle().Pt(); | |
350 | ||
351 | // distance to Prim. vertex | |
352 | const Double_t* dv = infoMC->GetVDist(); | |
353 | Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz(); | |
354 | ||
355 | // calculate and set prim. vertex | |
356 | vertexMC.SetXv( infoMC->GetParticle().Vx() - dv[0] ); | |
357 | vertexMC.SetYv( infoMC->GetParticle().Vy() - dv[1] ); | |
358 | vertexMC.SetZv( infoMC->GetParticle().Vz() - dv[2] ); | |
359 | ||
360 | // Check selection cuts | |
361 | if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return; | |
362 | if (!isPrim) return; | |
363 | if (infoRC->GetStatus(1)!=3) return; // TPC refit | |
364 | if (!infoRC->GetESDtrack()) return; | |
365 | if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; | |
366 | ||
367 | // constrained parameters resolution | |
368 | const AliExternalTrackParam * cparam = infoRC->GetESDtrack()->GetConstrainedParam(); | |
369 | if(!cparam) return; | |
370 | ||
371 | if ((track = new AliExternalTrackParam(*cparam)) != 0) | |
372 | { | |
373 | //Bool_t bDCAStatus = track->PropagateToDCA(&vertexMC,field,kMaxD,dca,cov); | |
374 | Double_t field[3]; track->GetBxByBz(field); | |
375 | Bool_t bDCAStatus = track->PropagateToDCABxByBz(&vertexMC,field,kMaxD,dca,cov); | |
376 | if(bDCAStatus) { | |
377 | if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ()) | |
378 | { | |
379 | deltaYTPC= track->GetY()-infoMC->GetParticle().Vy(); | |
380 | deltaZTPC = track->GetZ()-infoMC->GetParticle().Vz(); | |
381 | deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt()); | |
382 | deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px()); | |
383 | delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt; | |
384 | ||
385 | pullYTPC= (track->GetY()-infoMC->GetParticle().Vy()) / TMath::Sqrt(track->GetSigmaY2()); | |
386 | pullZTPC = (track->GetZ()-infoMC->GetParticle().Vz()) / TMath::Sqrt(track->GetSigmaZ2()); | |
387 | pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2()); | |
388 | pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2()); | |
389 | pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2()); | |
390 | ||
391 | Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt}; | |
392 | fResolHisto->Fill(vResolHisto); | |
393 | ||
394 | Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt}; | |
395 | fPullHisto->Fill(vPullHisto); | |
396 | } | |
397 | } | |
398 | delete track; | |
399 | } | |
400 | } | |
401 | ||
402 | //_____________________________________________________________________________ | |
403 | void AliComparisonRes::Exec(AliMCInfo* const infoMC, AliESDRecInfo* const infoRC){ | |
404 | ||
405 | // Process comparison information | |
406 | if(GetAnalysisMode() == 0) ProcessTPC(infoMC,infoRC); | |
407 | else if(GetAnalysisMode() == 1) ProcessTPCITS(infoMC,infoRC); | |
408 | else if(GetAnalysisMode() == 2) ProcessConstrained(infoMC,infoRC); | |
409 | else { | |
410 | printf("ERROR: AnalysisMode %d \n",fAnalysisMode); | |
411 | return; | |
412 | } | |
413 | } | |
414 | ||
415 | //_____________________________________________________________________________ | |
416 | TH1F* AliComparisonRes::MakeResol(TH2F * his, Int_t integ, Bool_t type){ | |
417 | // Create resolution histograms | |
418 | ||
419 | TH1F *hisr, *hism; | |
420 | if (!gPad) new TCanvas; | |
421 | hisr = AliTreeDraw::CreateResHistoI(his,&hism,integ); | |
422 | if (type) return hism; | |
423 | else | |
424 | return hisr; | |
425 | } | |
426 | ||
427 | //_____________________________________________________________________________ | |
428 | void AliComparisonRes::Analyse(){ | |
429 | // Analyse comparison information and store output histograms | |
430 | // in the folder "folderRes" | |
431 | // | |
432 | TH1::AddDirectory(kFALSE); | |
433 | TH1F *h=0; | |
434 | TH2F *h2D=0; | |
435 | TObjArray *aFolderObj = new TObjArray; | |
436 | ||
437 | // write results in the folder | |
438 | TCanvas * c = new TCanvas("Phi resol Tan","Phi resol Tan"); | |
439 | c->cd(); | |
440 | ||
441 | char name[256]; | |
442 | char res_title[256] = {"res_y:res_z:res_lambda:res_phi:res_1pt:y:z:eta:phi:pt"} ; | |
443 | char pull_title[256] = {"pull_y:pull_z:pull_lambda:pull_phi:pull_1pt:y:z:eta:phi:pt"}; | |
444 | ||
445 | for(Int_t i=0; i<5; i++) | |
446 | { | |
447 | for(Int_t j=5; j<10; j++) | |
448 | { | |
449 | if(j==7) fResolHisto->GetAxis(9)->SetRangeUser(0.2,10.); // pt threshold | |
450 | if(j==9) fResolHisto->GetAxis(7)->SetRangeUser(-0.9,0.9); // eta window | |
451 | ||
452 | h2D = (TH2F*)fResolHisto->Projection(i,j); | |
453 | h = AliComparisonRes::MakeResol(h2D,1,0); | |
454 | sprintf(name,"h_res_%d_vs_%d",i,j); | |
455 | h->SetName(name); | |
456 | h->SetTitle(res_title); | |
457 | ||
458 | aFolderObj->Add(h); | |
459 | ||
460 | h = AliComparisonRes::MakeResol(h2D,1,1); | |
461 | sprintf(name,"h_mean_res_%d_vs_%d",i,j); | |
462 | h->SetName(name); | |
463 | h->SetTitle(res_title); | |
464 | ||
465 | aFolderObj->Add(h); | |
466 | ||
467 | if(j==7) fResolHisto->GetAxis(9)->SetRangeUser(0.0,10.); | |
468 | if(j==9) fResolHisto->GetAxis(7)->SetRangeUser(-1.5,1.5); | |
469 | ||
470 | // | |
471 | if(j==7) fPullHisto->GetAxis(9)->SetRangeUser(0.2,10.); | |
472 | if(j==9) fPullHisto->GetAxis(7)->SetRangeUser(-0.9,0.9); | |
473 | ||
474 | h2D = (TH2F*)fPullHisto->Projection(i,j); | |
475 | h = AliComparisonRes::MakeResol(h2D,1,0); | |
476 | sprintf(name,"h_pull_%d_vs_%d",i,j); | |
477 | h->SetName(name); | |
478 | h->SetTitle(pull_title); | |
479 | ||
480 | aFolderObj->Add(h); | |
481 | ||
482 | h = AliComparisonRes::MakeResol(h2D,1,1); | |
483 | sprintf(name,"h_mean_pull_%d_vs_%d",i,j); | |
484 | h->SetName(name); | |
485 | h->SetTitle(pull_title); | |
486 | ||
487 | aFolderObj->Add(h); | |
488 | ||
489 | if(j==7) fPullHisto->GetAxis(9)->SetRangeUser(0.0,10.); | |
490 | if(j==9) fPullHisto->GetAxis(7)->SetRangeUser(-1.5,1.5); | |
491 | } | |
492 | } | |
493 | ||
494 | // export objects to analysis folder | |
495 | fAnalysisFolder = ExportToFolder(aFolderObj); | |
496 | ||
497 | // delete only TObjArray | |
498 | if(aFolderObj) delete aFolderObj; | |
499 | } | |
500 | ||
501 | //_____________________________________________________________________________ | |
502 | TFolder* AliComparisonRes::ExportToFolder(TObjArray * array) | |
503 | { | |
504 | // recreate folder avery time and export objects to new one | |
505 | // | |
506 | AliComparisonRes * comp=this; | |
507 | TFolder *folder = comp->GetAnalysisFolder(); | |
508 | ||
509 | TString name, title; | |
510 | TFolder *newFolder = 0; | |
511 | Int_t i = 0; | |
512 | Int_t size = array->GetSize(); | |
513 | ||
514 | if(folder) { | |
515 | // get name and title from old folder | |
516 | name = folder->GetName(); | |
517 | title = folder->GetTitle(); | |
518 | ||
519 | // delete old one | |
520 | delete folder; | |
521 | ||
522 | // create new one | |
523 | newFolder = CreateFolder(name.Data(),title.Data()); | |
524 | newFolder->SetOwner(); | |
525 | ||
526 | // add objects to folder | |
527 | while(i < size) { | |
528 | newFolder->Add(array->At(i)); | |
529 | i++; | |
530 | } | |
531 | } | |
532 | ||
533 | return newFolder; | |
534 | } | |
535 | ||
536 | //_____________________________________________________________________________ | |
537 | Long64_t AliComparisonRes::Merge(TCollection* const list) | |
538 | { | |
539 | // Merge list of objects (needed by PROOF) | |
540 | ||
541 | if (!list) | |
542 | return 0; | |
543 | ||
544 | if (list->IsEmpty()) | |
545 | return 1; | |
546 | ||
547 | TIterator* iter = list->MakeIterator(); | |
548 | TObject* obj = 0; | |
549 | ||
550 | // collection of generated histograms | |
551 | Int_t count=0; | |
552 | while((obj = iter->Next()) != 0) | |
553 | { | |
554 | AliComparisonRes* entry = dynamic_cast<AliComparisonRes*>(obj); | |
555 | if (entry == 0) continue; | |
556 | ||
557 | fResolHisto->Add(entry->fResolHisto); | |
558 | fPullHisto->Add(entry->fPullHisto); | |
559 | ||
560 | count++; | |
561 | } | |
562 | ||
563 | return count; | |
564 | } | |
565 | ||
566 | //_____________________________________________________________________________ | |
567 | TFolder* AliComparisonRes::CreateFolder(TString name,TString title) { | |
568 | // create folder for analysed histograms | |
569 | // | |
570 | TFolder *folder = 0; | |
571 | folder = new TFolder(name.Data(),title.Data()); | |
572 | ||
573 | return folder; | |
574 | } |