]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG1/AliPerformanceMC.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
adding cut on relative uncertainty on pt (Bastian)
[u/mrichter/AliRoot.git] / PWG1 / AliPerformanceMC.cxx
CommitLineData
29755947 1//------------------------------------------------------------------------------
2// Implementation of AliPerformanceMC 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 AliPerformanceMC.
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 AliPerformanceMC * compObj = (AliPerformanceMC*)coutput->FindObject("AliPerformanceMC");
21
22 // analyse comparison data
23 compObj->Analyse();
24
25 // the output histograms/graphs will be stored in the folder "folderRes"
26 compObj->GetAnalysisFolder()->ls("*");
27
28 // user can save whole comparison object (or only folder with anlysed histograms)
29 // in the seperate output file (e.g.)
30 TFile fout("Analysed_MC.root","recreate");
31 compObj->Write(); // compObj->GetAnalysisFolder()->Write();
32 fout.Close();
33
34*/
35
36#include "TCanvas.h"
37#include "TH1.h"
38#include "TH2.h"
39#include "TAxis.h"
40#include "TF1.h"
41
42#include "AliPerformanceMC.h"
43#include "AliESDEvent.h"
44#include "AliESDVertex.h"
45#include "AliESDtrack.h"
46#include "AliESDfriendTrack.h"
47#include "AliESDfriend.h"
48#include "AliLog.h"
49#include "AliMCEvent.h"
50#include "AliMCParticle.h"
51#include "AliHeader.h"
52#include "AliGenEventHeader.h"
53#include "AliStack.h"
54#include "AliMCInfoCuts.h"
55#include "AliRecInfoCuts.h"
56#include "AliTracker.h"
57#include "AliTreeDraw.h"
58#include "AliTPCseed.h"
59
60using namespace std;
61
62ClassImp(AliPerformanceMC)
63
64//_____________________________________________________________________________
65AliPerformanceMC::AliPerformanceMC():
66 AliPerformanceObject("AliPerformanceMC"),
67 fResolHisto(0),
68 fPullHisto(0),
69
70 // Cuts
71 fCutsRC(0),
72 fCutsMC(0),
73
74 // histogram folder
75 fAnalysisFolder(0)
76{
77 Init();
78}
79
80//_____________________________________________________________________________
81AliPerformanceMC::AliPerformanceMC(Char_t* name="AliPerformanceMC", Char_t* title="AliPerformanceMC",Int_t analysisMode=0,Bool_t hptGenerator=kFALSE):
82 AliPerformanceObject(name,title),
83 fResolHisto(0),
84 fPullHisto(0),
85
86 // Cuts
87 fCutsRC(0),
88 fCutsMC(0),
89
90 // histogram folder
91 fAnalysisFolder(0)
92{
93 // named constructor
94 //
95 SetAnalysisMode(analysisMode);
96 SetHptGenerator(hptGenerator);
97
98 Init();
99}
100
101//_____________________________________________________________________________
102AliPerformanceMC::~AliPerformanceMC()
103{
104 // destructor
105
106 if(fResolHisto) delete fResolHisto; fResolHisto=0;
107 if(fPullHisto) delete fPullHisto; fPullHisto=0;
108
109 if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
110}
111
112//_____________________________________________________________________________
113void AliPerformanceMC::Init(){
114
115 //
116 // histogram bining
117 //
118
119 // set pt bins
120 Int_t nPtBins = 50;
121 Double_t ptMin = 1.e-2, ptMax = 10.;
122
123 Double_t *binsPt = 0;
124 if (IsHptGenerator()) {
125 nPtBins = 100; ptMax = 100.;
126 binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
127 } else {
128 binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
129 }
130
131 Double_t yMin = -0.02, yMax = 0.02;
132 Double_t zMin = -12.0, zMax = 12.0;
133 if(GetAnalysisMode() == 3) { // TrackRef coordinate system
134 yMin = -100.; yMax = 100.;
135 zMin = -100.; zMax = 100.;
136 }
137
138 // res_y:res_z:res_phi,res_lambda:res_pt:y:z:eta:phi:pt
139 Int_t binsResolHisto[10]={100,100,100,100,100,25,50,144,30,nPtBins};
140 Double_t minResolHisto[10]={-0.2,-0.2,-0.002,-0.002,-0.002, yMin, zMin, 0., -1.5, ptMin};
141 Double_t maxResolHisto[10]={ 0.2, 0.2, 0.002, 0.002, 0.002, yMax, zMax, 2.*TMath::Pi(), 1.5, ptMax};
142
143 fResolHisto = new THnSparseF("fResolHisto","res_y:res_z:res_phi:res_lambda:res_pt:y:z:phi:eta:pt",10,binsResolHisto,minResolHisto,maxResolHisto);
144 fResolHisto->SetBinEdges(9,binsPt);
145
146 fResolHisto->GetAxis(0)->SetTitle("y-y_{mc} (cm)");
147 fResolHisto->GetAxis(1)->SetTitle("z-z_{mc} (cm)");
148 fResolHisto->GetAxis(2)->SetTitle("#phi-#phi_{mc} (rad)");
149 fResolHisto->GetAxis(3)->SetTitle("#lambda-#lambda_{mc} (rad)");
150 fResolHisto->GetAxis(4)->SetTitle("(p_{T}/p_{Tmc}-1)");
151 fResolHisto->GetAxis(5)->SetTitle("y_{mc} (cm)");
152 fResolHisto->GetAxis(6)->SetTitle("z_{mc} (cm)");
153 fResolHisto->GetAxis(7)->SetTitle("#phi_{mc} (rad)");
154 fResolHisto->GetAxis(8)->SetTitle("#eta_{mc}");
155 fResolHisto->GetAxis(9)->SetTitle("p_{Tmc} (GeV/c)");
156 fResolHisto->Sumw2();
157
158 //pull_y:pull_z:pull_snp:pull_tgl:pull_1pt:y:z:snp:tgl:1pt
159 Int_t binsPullHisto[10]={100,100,100,100,100,50,50,50,50,nPtBins};
160 Double_t minPullHisto[10]={-5.,-5.,-5.,-5.,-5., yMin, zMin,-1., -2.0, ptMin};
161 Double_t maxPullHisto[10]={ 5., 5., 5., 5., 5., yMax, zMax, 1., 2.0, ptMax};
162 fPullHisto = new THnSparseF("fPullHisto","pull_y:pull_z:pull_y:pull_z:pull_snp:pull_tgl:pull_1pt:y:z:snp:tgl:1pt",10,binsPullHisto,minPullHisto,maxPullHisto);
163
164 fPullHisto->GetAxis(0)->SetTitle("(y-y_{mc})/#sigma");
165 fPullHisto->GetAxis(1)->SetTitle("(z-z_{mc})/#sigma");
166 fPullHisto->GetAxis(2)->SetTitle("(sin#phi-sin#phi_{mc})/#sigma");
167 fPullHisto->GetAxis(3)->SetTitle("(tan#lambda-tan#lambda_{mc})/#sigma");
168 fPullHisto->GetAxis(4)->SetTitle("(p_{Tmc}/p_{T}-1)/#sigma");
169 fPullHisto->GetAxis(5)->SetTitle("y_{mc} (cm)");
170 fPullHisto->GetAxis(6)->SetTitle("z_{mc} (cm)");
171 fPullHisto->GetAxis(7)->SetTitle("sin#phi_{mc}");
172 fPullHisto->GetAxis(8)->SetTitle("tan#lambda_{mc}");
173 fPullHisto->GetAxis(9)->SetTitle("1/p_{Tmc} (GeV/c)^{-1}");
174 fPullHisto->Sumw2();
175
176 // Init cuts
177 if(!fCutsMC)
178 AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
179 if(!fCutsRC)
180 AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
181
182 // init folder
183 fAnalysisFolder = CreateFolder("folderMC","Analysis MC Folder");
184}
185
186//_____________________________________________________________________________
187void AliPerformanceMC::Exec(AliMCEvent* const mcEvent, AliESDEvent *const /*esdEvent*/, AliESDfriend *const /*esdFriend*/, const Bool_t bUseMC, const Bool_t /*bUseESDfriend*/)
188{
189 //
190 // Process pure MC information
191 //
192 AliHeader* header = 0;
193 AliGenEventHeader* genHeader = 0;
194 AliStack* stack = 0;
195 TArrayF vtxMC(3);
196
197 if(bUseMC)
198 {
199 if(!mcEvent) {
200 AliDebug(AliLog::kError, "mcEvent not available");
201 return;
202 }
203
204 // get MC event header
205 header = mcEvent->Header();
206 if (!header) {
207 AliDebug(AliLog::kError, "Header not available");
208 return;
209 }
210 // MC particle stack
211 stack = mcEvent->Stack();
212 if (!stack) {
213 AliDebug(AliLog::kError, "Stack not available");
214 return;
215 }
216
217 // get MC vertex
218 genHeader = header->GenEventHeader();
219 if (!genHeader) {
220 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
221 return;
222 }
223 genHeader->PrimaryVertex(vtxMC);
224
225 } // end bUseMC
226
227 Int_t nPart = mcEvent->GetNumberOfTracks();
228 if (nPart==0) return;
229
230 //TParticle * part = new TParticle;
231 //TClonesArray * trefs = new TClonesArray("AliTrackReference");
232 TParticle *part=0;
233 TClonesArray *trefs=0;
234
235 // Process events
236 for (Int_t iPart = 0; iPart < nPart; iPart++)
237 {
238 Int_t status = mcEvent->GetParticleAndTR(iPart, part, trefs);
239 if (status<0) continue;
240 if(!part) continue;
241 if(!trefs) continue;
242
7aad0c47 243 AliMCParticle *mcParticle = (AliMCParticle*) mcEvent->GetTrack(iPart);
29755947 244 if(!mcParticle) continue;
245
246 TParticle *particle = mcParticle->Particle();
247 if(!particle) continue;
248
249 Bool_t prim = stack->IsPhysicalPrimary(iPart);
250
251 // Only 5 charged particle species (e,mu,pi,K,p)
252 if (fCutsMC->IsPdgParticle(TMath::Abs(particle->GetPdgCode())) == kFALSE) break;
253
254 // exclude electrons
255 if (fCutsMC->GetEP()==TMath::Abs(particle->GetPdgCode())) return;
256
257 AliTrackReference *refTPCIn = 0;
258 AliTrackReference *refTPCOut = 0;
259 AliTrackReference *refITSIn = 0;
260 AliTrackReference *refITSOut = 0;
261 AliTrackReference *refTRDIn = 0;
262 AliTrackReference *refTRDOut = 0;
263
264 Float_t rmax = 0.;
265 Int_t nTrackRef = mcParticle->GetNumberOfTrackReferences();
266 if(nTrackRef < 5) continue;
267
268 for (Int_t iref = 0; iref < nTrackRef; iref++)
269 {
270 AliTrackReference *ref = mcParticle->GetTrackReference(iref);
271 if(!ref) continue;
272
273 Float_t dir = ref->X()*ref->Px()+ref->Y()*ref->Py();
274 if(dir < 0.) break; // looping tracks (return back)
275 if (ref->R()<rmax) break;
276
277 // pt threshold
278 if(ref->Pt() < fCutsRC->GetPtMin()) break;
279
280 // TPC
281 if(ref->DetectorId()==AliTrackReference::kTPC)
282 {
283 if(!refTPCIn) {
284 refTPCIn = ref;
285 }
286 else {
287 if(ref->R() > refTPCIn->R())
288 refTPCOut = ref;
289 }
290 }
291 // ITS
292 if(ref->DetectorId()==AliTrackReference::kITS)
293 {
294 if(!refITSIn) {
295 refITSIn = ref;
296 }
297 else {
298 if(ref->R() > refITSIn->R())
299 refITSOut = ref;
300 }
301 }
302 // TRD
303 if(ref->DetectorId()==AliTrackReference::kTRD)
304 {
305 if(!refTRDIn) {
306 refTRDIn = ref;
307 }
308 else {
309 if(ref->R() > refTRDIn->R())
310 refTRDOut = ref;
311 }
312 }
313 if (ref->R()>rmax) rmax=ref->R();
314 }
315
316 if(GetAnalysisMode() == 0) {
317 // Propagate refTPCIn to DCA to prim. vertex and make comparison
318 // only primary particles
319 if(refTPCIn) {
320 if(prim) ProcessTPC(refTPCIn,particle);
321 }
322 }
323
324 if(GetAnalysisMode() == 3) {
325 // Propagate refTPCOut to refTPCIn using AliTrack::PropagateTrackTo()
326 //
327 if(refTPCIn && refTPCOut)
328 ProcessInnerTPC(refTPCIn,refTPCOut,particle);
329 }
330
331 if(GetAnalysisMode() == 4) {
332 // propagate refTPCIn to refTPCOut using AliExternalTrackParam::PropagateTo()
333 ProcessOuterTPCExt(part,trefs);
334 }
335 }
336
337 //if(part) delete part;
338 //if(trefs) delete trefs;
339}
340
341//_____________________________________________________________________________
342void AliPerformanceMC::ProcessTPC(AliTrackReference* const refIn, TParticle *const particle) {
343 //
344 // Test propagation from refIn to DCA
345 //
346 if(!refIn) return;
347 if(!particle) return;
348
349 AliExternalTrackParam *track = 0;
350 Double_t radius = 2.8; // cm
351 Double_t mass = particle->GetMass();
352 Double_t step=1;
353 //
354 track=MakeTrack(refIn,particle);
355 if (!track) return;
356 //
357 AliTracker::PropagateTrackTo(track, radius+step, mass, step, kTRUE,0.99);
358 AliTracker::PropagateTrackTo(track, radius+0.5, mass, step*0.1, kTRUE,0.99);
359 Double_t xyz[3] = {particle->Vx(),particle->Vy(),particle->Vz()};
360 Double_t sxyz[3]={0.0,0.0,0.0};
361 AliESDVertex vertex(xyz,sxyz);
362 Double_t dca[2], cov[3];
363 Bool_t isOK = track->PropagateToDCA(&vertex,AliTracker::GetBz(),10,dca,cov);
364 if(!isOK) return;
365
366 // Fill histogram
367 Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
368 Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
369
370 Float_t mceta = particle->Eta();
371 Float_t mcphi = particle->Phi();
372 if(mcphi<0) mcphi += 2.*TMath::Pi();
373 Float_t mcpt = particle->Pt();
374 Float_t mcsnp = TMath::Sin(TMath::ATan2(particle->Py(),particle->Px()));
375 Float_t mctgl = TMath::Tan(TMath::ATan2(particle->Pz(),particle->Pt()));
376
377 //if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
378 //{
379 if(mcpt == 0) return;
380
381 //deltaYTPC= track->GetY()-particle->Vy();
382 deltaYTPC= track->GetY(); // it is already rotated
383 deltaZTPC = track->GetZ()-particle->Vz();
384 deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
385 deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(particle->Py(),particle->Px());
386 deltaPtTPC = (track->Pt()-mcpt) / mcpt;
387
388 /*
389 printf("------------------------------ \n");
390 printf("trY %f, trZ %f, trTheta %f ,trPhi %f, trPt %f \n", track->GetY(), track->GetZ(), TMath::ATan2(track->Pz(),track->Pt()), TMath::ATan2(track->Py(),track->Px()), track->Pt() );
391 printf("partY %f, partZ %f, partTheta %f ,partPhi %f, partPt %f \n", particle->Vy(), particle->Vz(), TMath::ATan2(particle->Pz(),particle->Pt()), TMath::ATan2(particle->Py(),particle->Px()), mcpt );
392 */
393
394
395 pullYTPC= (track->GetY()-particle->Vy()) / TMath::Sqrt(track->GetSigmaY2());
396 pullZTPC = (track->GetZ()-particle->Vz()) / TMath::Sqrt(track->GetSigmaZ2());
397
398 pullPhiTPC = (track->GetSnp() - mcsnp) / TMath::Sqrt(track->GetSigmaSnp2());
399 pullLambdaTPC = (track->GetTgl() - mctgl) / TMath::Sqrt(track->GetSigmaTgl2());
400
401 if (mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
402 else pull1PtTPC = 0.;
403
404 Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,particle->Vy(),particle->Vz(),mcphi,mceta,mcpt};
405 fResolHisto->Fill(vResolHisto);
406
407 Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mcsnp,mctgl,1./mcpt};
408 fPullHisto->Fill(vPullHisto);
409 //}
410 if(track) delete track;
411}
412
413//_____________________________________________________________________________
414void AliPerformanceMC::ProcessInnerTPC(AliTrackReference* const refIn, AliTrackReference *const refOut, TParticle *const particle) {
415 //
416 // Test propagation from Out to In
417 //
418 if(!refIn) return;
419 if(!refOut) return;
420 if(!particle) return;
421
422 AliExternalTrackParam *track=MakeTrack(refOut,particle);
423 if (!track) return;
424
425 Double_t mass = particle->GetMass();
426 Double_t step=1;
427
428 Float_t alphaIn= TMath::ATan2(refIn->Y(),refIn->X());
429 //Float_t alphaOut= TMath::ATan2(refOut->Y(),refOut->X());
430 //track->Rotate(alphaIn-alphaOut);
431 track->Rotate(alphaIn);
432 //printf("alphaIn %f, alphaOut %f \n",alphaIn,alphaOut);
433
434 //
435 //Bool_t isOK = AliTracker::PropagateTrackTo(track, refIn->R(), mass, step, kTRUE,0.99);
436 Bool_t isOK = AliTracker::PropagateTrackTo(track, refIn->R(), mass, step, kFALSE,0.99);
437 if(!isOK) return;
438
439 // calculate alpha angle
440 //Double_t xyz[3] = {refIn->X(),refIn->Y(),refIn0->Z()};
441 //Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
442 //if (alpha<0) alpha += TMath::TwoPi();
443
444 // rotate outer track to inner
445 // and propagate to the radius of the first track referenco of TPC
446 //Double_t trRadius = TMath::Sqrt(xyz[1] * xyz[1] + xyz[0] * xyz[0]);
447 //Bool_t isOK = track->Propagate(alpha,trRadius,AliTracker::GetBz());
448 //if(!isOK) return;
449
450 Float_t mceta = -TMath::Log(TMath::Tan(0.5 * refIn->Theta()));
451 Float_t mcphi = refIn->Phi();
452 if(mcphi<0) mcphi += 2.*TMath::Pi();
453 Float_t mcpt = refIn->Pt();
454 Float_t mcsnp = TMath::Sin(TMath::ATan2(refIn->Py(),refIn->Px()));
455 Float_t mctgl = TMath::Tan(TMath::ATan2(refIn->Pz(),refIn->Pt()));
456
457 Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
458 Float_t pull1PtTPC=0., pullYTPC=0., pullZTPC=0., pullPhiTPC=0., pullLambdaTPC=0.;
459
460 if(mcpt == 0) return;
461
462 //deltaYTPC = track->GetY()-refIn->Y();
463 deltaYTPC = track->GetY();
464 deltaZTPC = track->GetZ()-refIn->Z();
465 deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(refIn->Pz(),refIn->Pt());
466 deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(refIn->Py(),refIn->Px());
467 deltaPtTPC = (track->Pt()-mcpt) / mcpt;
468
469 /*
470 printf("------------------------------ \n");
471 printf("trX %f, trY %f, trZ %f, trTheta %f ,trPhi %f, trPt %f \n", track->GetX(), track->GetY(), track->GetZ(), TMath::ATan2(track->Pz(),track->Pt()), TMath::ATan2(track->Py(),track->Px()), track->Pt() );
472 printf("partX %f, partY %f, partZ %f, partTheta %f ,partPhi %f, partPt %f \n",refIn->X(), refIn->Y(), refIn->Z(), TMath::ATan2(refIn->Pz(),refIn->Pt()), TMath::ATan2(refIn->Py(),refIn->Px()), mcpt );
473 */
474
475 if(TMath::Sqrt(track->GetSigmaY2())) pullYTPC = track->GetY() / TMath::Sqrt(track->GetSigmaY2());
476 if(TMath::Sqrt(track->GetSigmaZ2())) pullZTPC = (track->GetZ()-refIn->Z()) / TMath::Sqrt(track->GetSigmaZ2());
477 if(TMath::Sqrt(track->GetSigmaSnp2())) pullPhiTPC = (track->GetSnp() - mcsnp) / TMath::Sqrt(track->GetSigmaSnp2());
478 if(TMath::Sqrt(track->GetSigmaTgl2())) pullLambdaTPC = (track->GetTgl() - mctgl) / TMath::Sqrt(track->GetSigmaTgl2());
479 if(TMath::Sqrt(track->GetSigma1Pt2())) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
480
481 //printf("pullYTPC %f,pullZTPC %f,pullPhiTPC %f,pullLambdaTPC %f,pull1PtTPC %f,refIn->Y() %f,refIn->Z() %f,mcsnp %f,mctgl %f,1./mcpt %f \n",pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,refIn->Y(),refIn->Z(),mcsnp,mctgl,1./mcpt);
482
483 Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,refIn->Y(),refIn->Z(),mcphi,mceta,mcpt};
484 fResolHisto->Fill(vResolHisto);
485
486 Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,refIn->Y(),refIn->Z(),mcsnp,mctgl,1./mcpt};
487 fPullHisto->Fill(vPullHisto);
488
489 if(track) delete track;
490}
491
492//_____________________________________________________________________________
493void AliPerformanceMC::ProcessOuterTPCExt(TParticle *const part, TClonesArray * const trefs)
494{
495 const Int_t kMinRefs=5;
496 Int_t nrefs = trefs->GetEntries();
497 if (nrefs<kMinRefs) return; // we should have enough references
498 Int_t iref0 =-1;
499 Int_t iref1 =-1;
500
501 for (Int_t iref=0; iref<nrefs; iref++){
502 AliTrackReference * ref = (AliTrackReference*)trefs->At(iref);
503 if (!ref) continue;
504 Float_t dir = ref->X()*ref->Px()+ref->Y()*ref->Py();
505 if (dir<0) break;
506 if (ref->DetectorId()!=AliTrackReference::kTPC) continue;
507 if (iref0<0) iref0 = iref;
508 iref1 = iref;
509 }
510 if (iref1-iref0<kMinRefs) return;
511 Double_t covar[14];
512 for (Int_t icov=0; icov<14; icov++) covar[icov]=0;
513 covar[0]=1;
514 covar[2]=1;
515 covar[5]=1;
516 covar[9]=1;
517 covar[14]=1;
518
519 AliTrackReference * refIn = (AliTrackReference*)trefs->At(iref0);
520 AliTrackReference * refOut = (AliTrackReference*)trefs->At(iref1);
521 AliExternalTrackParam *paramPropagate= MakeTrack(refIn,part);
522 AliExternalTrackParam *paramUpdate = MakeTrack(refIn,part);
523 paramUpdate->AddCovariance(covar);
524
525 Bool_t isOKP=kTRUE;
526 Bool_t isOKU=kTRUE;
527 AliMagF * field = (AliMagF*) TGeoGlobalMagField::Instance()->GetField();
528 if(!field) return;
529 for (Int_t iref = iref0; iref<=iref1; iref++){
530 AliTrackReference * ref = (AliTrackReference*)trefs->At(iref);
531 if(!ref) continue;
532 Float_t alphaC= TMath::ATan2(ref->Y(),ref->X());
533 Double_t pos[3] = {ref->X(), ref->Y(), ref->Z()};
534 Double_t mag[3];
535 field->Field(pos,mag);
536 isOKP&=paramPropagate->Rotate(alphaC);
537 isOKU&=paramUpdate->Rotate(alphaC);
538 for (Float_t xref= paramPropagate->GetX(); xref<ref->R(); xref++){
539 isOKP&=paramPropagate->PropagateTo(xref, mag[2]);
540 isOKU&=paramUpdate->PropagateTo(xref, mag[2]);
541 }
542 isOKP&=paramPropagate->PropagateTo(ref->R(), mag[2]);
543 isOKU&=paramUpdate->PropagateTo(ref->R(), mag[2]);
544 Double_t clpos[2] = {0, ref->Z()};
545 Double_t clcov[3] = { 0.005,0,0.005};
546 isOKU&= paramUpdate->Update(clpos, clcov);
547 }
548
549 Float_t mceta = -TMath::Log(TMath::Tan(0.5 * refOut->Theta()));
550 Float_t mcphi = refOut->Phi();
551 if(mcphi<0) mcphi += 2.*TMath::Pi();
552 Float_t mcpt = refOut->Pt();
553 Float_t mcsnp = TMath::Sin(TMath::ATan2(refOut->Py(),refOut->Px()));
554 Float_t mctgl = TMath::Tan(TMath::ATan2(refOut->Pz(),refOut->Pt()));
555
556 Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
557 Float_t pull1PtTPC=0., pullYTPC=0., pullZTPC=0., pullPhiTPC=0., pullLambdaTPC=0.;
558
559 if(mcpt == 0) return;
560
561 //deltaYTPC = track->GetY()-refIn->Y();
562 deltaYTPC = paramPropagate->GetY();
563 deltaZTPC = paramPropagate->GetZ()-refOut->Z();
564 deltaLambdaTPC = TMath::ATan2(paramPropagate->Pz(),paramPropagate->Pt())-TMath::ATan2(refOut->Pz(),refOut->Pt());
565 deltaPhiTPC = TMath::ATan2(paramPropagate->Py(),paramPropagate->Px())-TMath::ATan2(refOut->Py(),refOut->Px());
566 deltaPtTPC = (paramPropagate->Pt()-mcpt) / mcpt;
567
568 if(TMath::Sqrt(paramPropagate->GetSigmaY2())) pullYTPC = paramPropagate->GetY() / TMath::Sqrt(paramPropagate->GetSigmaY2());
569 if(TMath::Sqrt(paramPropagate->GetSigmaZ2())) pullZTPC = (paramPropagate->GetZ()-refOut->Z()) / TMath::Sqrt(paramPropagate->GetSigmaZ2());
570 if(TMath::Sqrt(paramPropagate->GetSigmaSnp2())) pullPhiTPC = (paramPropagate->GetSnp() - mcsnp) / TMath::Sqrt(paramPropagate->GetSigmaSnp2());
571 if(TMath::Sqrt(paramPropagate->GetSigmaTgl2())) pullLambdaTPC = (paramPropagate->GetTgl() - mctgl) / TMath::Sqrt(paramPropagate->GetSigmaTgl2());
572 if(TMath::Sqrt(paramPropagate->GetSigma1Pt2())) pull1PtTPC = (paramPropagate->OneOverPt()-1./mcpt) / TMath::Sqrt(paramPropagate->GetSigma1Pt2());
573
574 //printf("pullYTPC %f,pullZTPC %f,pullPhiTPC %f,pullLambdaTPC %f,pull1PtTPC %f,refIn->Y() %f,refIn->Z() %f,mcsnp %f,mctgl %f,1./mcpt %f \n",pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,refIn->Y(),refIn->Z(),mcsnp,mctgl,1./mcpt);
575
576 Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,refIn->Y(),refIn->Z(),mcphi,mceta,mcpt};
577 fResolHisto->Fill(vResolHisto);
578
579 Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,refIn->Y(),refIn->Z(),mcsnp,mctgl,1./mcpt};
580 fPullHisto->Fill(vPullHisto);
581}
582
583//_____________________________________________________________________________
584AliExternalTrackParam * AliPerformanceMC::MakeTrack(const AliTrackReference* const ref, TParticle* const part)
585{
586 //
587 // Make track out of the track ref
588 // part - TParticle used to determine chargr
589 // the covariance matrix - equal 0 - starting from ideal MC position
590 Double_t xyz[3]={ref->X(),ref->Y(),ref->Z()};
591 Double_t pxyz[3]={ref->Px(),ref->Py(),ref->Pz()};
592 Double_t cv[21];
593 for (Int_t i=0; i<21;i++) cv[i]=0;
594 if (!part->GetPDG()) return 0;
595 AliExternalTrackParam * param = new AliExternalTrackParam(xyz,pxyz,cv,TMath::Nint(part->GetPDG()->Charge()/3.));
596 return param;
597}
598
599//_____________________________________________________________________________
600TH1F* AliPerformanceMC::MakeResol(TH2F * his, Int_t integ, Bool_t type, Int_t cut){
601 // Create resolution histograms
602
603 TH1F *hisr, *hism;
604 if (!gPad) new TCanvas;
605 hisr = AliTreeDraw::CreateResHistoII(his,&hism,integ,kTRUE,cut);
606 if (type) return hism;
607 else
608 return hisr;
609}
610
611//_____________________________________________________________________________
612void AliPerformanceMC::Analyse() {
613 // Analyse comparison information and store output histograms
614 // in the folder "folderRes"
615 //
616 TH1::AddDirectory(kFALSE);
617 TH1F *h=0;
618 TH2F *h2D=0;
619 TObjArray *aFolderObj = new TObjArray;
620
621 // write results in the folder
622 TCanvas * c = new TCanvas("Phi resol Tan","Phi resol Tan");
623 c->cd();
624
625 char name[256];
626 char title[256];
627
628 for(Int_t i=0; i<5; i++)
629 {
630 for(Int_t j=5; j<10; j++)
631 {
632 //if(j!=8) fResolHisto->GetAxis(8)->SetRangeUser(-0.9,0.89); // eta window
633 fResolHisto->GetAxis(9)->SetRangeUser(0.16,10.); // pt threshold
634 if(j!=8) fResolHisto->GetAxis(8)->SetRangeUser(0.,0.9); // eta window
635 //fResolHisto->GetAxis(9)->SetRangeUser(0.16,3.); // pt window
636 if(GetAnalysisMode() == 3) fResolHisto->GetAxis(5)->SetRangeUser(-80.,80.); // y range
637
638 h2D = (TH2F*)fResolHisto->Projection(i,j);
639
640 h = AliPerformanceMC::MakeResol(h2D,1,0,20);
641 sprintf(name,"h_res_%d_vs_%d",i,j);
642 h->SetName(name);
643
644 h->GetXaxis()->SetTitle(fResolHisto->GetAxis(j)->GetTitle());
645 sprintf(title,"%s %s",fResolHisto->GetAxis(i)->GetTitle(),"(resolution)");
646 h->GetYaxis()->SetTitle(title);
647 sprintf(title,"%s vs %s",title,fResolHisto->GetAxis(j)->GetTitle());
648 h->SetTitle(title);
649
650 if(j==9) h->SetBit(TH1::kLogX);
651 aFolderObj->Add(h);
652
653 h = AliPerformanceMC::MakeResol(h2D,1,1,20);
654 //h = (TH1F*)arr->At(1);
655 sprintf(name,"h_mean_res_%d_vs_%d",i,j);
656 h->SetName(name);
657
658 h->GetXaxis()->SetTitle(fResolHisto->GetAxis(j)->GetTitle());
659 sprintf(title,"%s %s",fResolHisto->GetAxis(i)->GetTitle(),"(mean)");
660 h->GetYaxis()->SetTitle(title);
661
662 sprintf(title,"%s vs %s",title,fResolHisto->GetAxis(j)->GetTitle());
663 h->SetTitle(title);
664
665 if(j==9) h->SetBit(TH1::kLogX);
666 aFolderObj->Add(h);
667
668 fResolHisto->GetAxis(8)->SetRangeUser(-1.5,1.5);
669 fResolHisto->GetAxis(9)->SetRangeUser(0.0,10.);
670
671 //
672 if(j!=8) fPullHisto->GetAxis(8)->SetRangeUser(-1.0,0.99); // theta window
673 else fPullHisto->GetAxis(8)->SetRangeUser(-1.5,1.5); // theta window
674 fPullHisto->GetAxis(9)->SetRangeUser(0.0,10.); // pt threshold
675 if(GetAnalysisMode() == 3) fPullHisto->GetAxis(5)->SetRangeUser(-80.,80.); // y range
676
677 h2D = (TH2F*)fPullHisto->Projection(i,j);
678
679 h = AliPerformanceMC::MakeResol(h2D,1,0,20);
680 sprintf(name,"h_pull_%d_vs_%d",i,j);
681 h->SetName(name);
682
683 h->GetXaxis()->SetTitle(fPullHisto->GetAxis(j)->GetTitle());
684 sprintf(title,"%s %s",fPullHisto->GetAxis(i)->GetTitle(),"(resolution)");
685 h->GetYaxis()->SetTitle(title);
686 sprintf(title,"%s vs %s",title,fPullHisto->GetAxis(j)->GetTitle());
687 h->SetTitle(title);
688
689 //if(j==9) h->SetBit(TH1::kLogX);
690 aFolderObj->Add(h);
691
692 h = AliPerformanceMC::MakeResol(h2D,1,1,20);
693 sprintf(name,"h_mean_pull_%d_vs_%d",i,j);
694 h->SetName(name);
695
696 h->GetXaxis()->SetTitle(fPullHisto->GetAxis(j)->GetTitle());
697 sprintf(title,"%s %s",fPullHisto->GetAxis(i)->GetTitle(),"(mean)");
698 h->GetYaxis()->SetTitle(title);
699 sprintf(title,"%s vs %s",title,fPullHisto->GetAxis(j)->GetTitle());
700 h->SetTitle(title);
701
702 //if(j==9) h->SetBit(TH1::kLogX);
703 aFolderObj->Add(h);
704 }
705 }
706
707 // export objects to analysis folder
708 fAnalysisFolder = ExportToFolder(aFolderObj);
709
710 // delete only TObjArray
711 if(aFolderObj) delete aFolderObj;
712}
713
714//_____________________________________________________________________________
715TFolder* AliPerformanceMC::ExportToFolder(TObjArray * array)
716{
717 // recreate folder avery time and export objects to new one
718 //
719 AliPerformanceMC * comp=this;
720 TFolder *folder = comp->GetAnalysisFolder();
721
722 TString name, title;
723 TFolder *newFolder = 0;
724 Int_t i = 0;
725 Int_t size = array->GetSize();
726
727 if(folder) {
728 // get name and title from old folder
729 name = folder->GetName();
730 title = folder->GetTitle();
731
732 // delete old one
733 delete folder;
734
735 // create new one
736 newFolder = CreateFolder(name.Data(),title.Data());
737 newFolder->SetOwner();
738
739 // add objects to folder
740 while(i < size) {
741 newFolder->Add(array->At(i));
742 i++;
743 }
744 }
745
746return newFolder;
747}
748
749//_____________________________________________________________________________
750Long64_t AliPerformanceMC::Merge(TCollection* const list)
751{
752 // Merge list of objects (needed by PROOF)
753
754 if (!list)
755 return 0;
756
757 if (list->IsEmpty())
758 return 1;
759
760 TIterator* iter = list->MakeIterator();
761 TObject* obj = 0;
762
763 // collection of generated histograms
764 Int_t count=0;
765 while((obj = iter->Next()) != 0)
766 {
767 AliPerformanceMC* entry = dynamic_cast<AliPerformanceMC*>(obj);
768 if (entry == 0) continue;
769
770 fResolHisto->Add(entry->fResolHisto);
771 fPullHisto->Add(entry->fPullHisto);
772
773 count++;
774 }
775
776return count;
777}
778
779//_____________________________________________________________________________
780TFolder* AliPerformanceMC::CreateFolder(TString name,TString title) {
781// create folder for analysed histograms
782//
783TFolder *folder = 0;
784 folder = new TFolder(name.Data(),title.Data());
785
786 return folder;
787}