- changes by Markus Koehler
[u/mrichter/AliRoot.git] / PWG1 / TPC / AliPerformanceMC.cxx
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
60 using namespace std;
61
62 ClassImp(AliPerformanceMC)
63
64 //_____________________________________________________________________________
65 AliPerformanceMC::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 //_____________________________________________________________________________
81 AliPerformanceMC::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 //_____________________________________________________________________________
102 AliPerformanceMC::~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 //_____________________________________________________________________________
113 void 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 = 20.;
122
123   Double_t *binsPt = 0;
124
125   if (IsHptGenerator())  { 
126         ptMax = 100.;
127   } 
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,90,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 //_____________________________________________________________________________
187 void AliPerformanceMC::Exec(AliMCEvent* const mcEvent, AliESDEvent *const /*esdEvent*/, AliESDfriend *const /*esdFriend*/, const Bool_t bUseMC, const Bool_t /*bUseESDfriend*/)
188 {
189   // Process pure MC information 
190   //
191   AliHeader* header = 0;
192   AliGenEventHeader* genHeader = 0;
193   AliStack* stack = 0;
194   TArrayF vtxMC(3);
195   
196   if(bUseMC)
197   {
198     if(!mcEvent) {
199       Error("Exec","mcEvent not available");
200       return;
201     }
202     // get MC event header
203     header = mcEvent->Header();
204     if (!header) {
205       Error("Exec","Header not available");
206       return;
207     }
208     // MC particle stack
209     stack = mcEvent->Stack();
210     if (!stack) {
211       Error("Exec","Stack not available");
212       return;
213     }
214     // get MC vertex
215     genHeader = header->GenEventHeader();
216     if (!genHeader) {
217       Error("Exec","Could not retrieve genHeader from Header");
218       return;
219     }
220     genHeader->PrimaryVertex(vtxMC);
221   } 
222   else {
223     Error("Exec","MC information required!");
224     return;
225   }
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
243     AliMCParticle *mcParticle = (AliMCParticle*) mcEvent->GetTrack(iPart);
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::PropagateTrackToBxByBz()
326       //
327       if(refTPCIn && refTPCOut) 
328         ProcessInnerTPC(refTPCIn,refTPCOut,particle); 
329     }
330
331     if(GetAnalysisMode() == 4) { 
332       // propagate refTPCIn to refTPCOut using AliExternalTrackParam::PropagateToBxByBz()
333       ProcessOuterTPCExt(part,trefs); 
334     }
335   }
336
337   //if(part) delete part;
338   //if(trefs) delete trefs;
339 }
340
341 //_____________________________________________________________________________
342 void 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   AliTracker::PropagateTrackToBxByBz(track, radius+step, mass, step, kTRUE,0.99);
360   AliTracker::PropagateTrackToBxByBz(track, radius+0.5, mass, step*0.1, kTRUE,0.99);
361   Double_t xyz[3] = {particle->Vx(),particle->Vy(),particle->Vz()};
362   Double_t sxyz[3]={0.0,0.0,0.0};
363   AliESDVertex vertex(xyz,sxyz);
364   Double_t dca[2], cov[3];
365   //Bool_t isOK = track->PropagateToDCA(&vertex,AliTracker::GetBz(),10,dca,cov);
366   Double_t field[3];  track->GetBxByBz(field); 
367   Bool_t isOK = track->PropagateToDCABxByBz(&vertex,field,10,dca,cov);
368   if(!isOK) return;
369
370   // Fill histogram
371   Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC; 
372   Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC; 
373
374   Float_t mceta =  particle->Eta();
375   Float_t mcphi =  particle->Phi();
376   if(mcphi<0) mcphi += 2.*TMath::Pi();
377   Float_t mcpt = particle->Pt();
378   Float_t mcsnp = TMath::Sin(TMath::ATan2(particle->Py(),particle->Px()));
379   Float_t mctgl = TMath::Tan(TMath::ATan2(particle->Pz(),particle->Pt()));
380
381   //if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ()) 
382   //{ 
383     if(mcpt == 0) return;
384     
385     //deltaYTPC= track->GetY()-particle->Vy();
386     deltaYTPC= track->GetY(); // it is already rotated
387     deltaZTPC = track->GetZ()-particle->Vz();
388     deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
389     deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(particle->Py(),particle->Px());
390     deltaPtTPC = (track->Pt()-mcpt) / mcpt;
391
392     /*
393     printf("------------------------------ \n");
394     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() ); 
395     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 ); 
396     */
397
398
399     pullYTPC= (track->GetY()-particle->Vy()) / TMath::Sqrt(track->GetSigmaY2());
400     pullZTPC = (track->GetZ()-particle->Vz()) / TMath::Sqrt(track->GetSigmaZ2());
401  
402     pullPhiTPC = (track->GetSnp() - mcsnp) / TMath::Sqrt(track->GetSigmaSnp2());
403     pullLambdaTPC = (track->GetTgl() - mctgl) / TMath::Sqrt(track->GetSigmaTgl2());
404
405     if (mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
406     else pull1PtTPC = 0.;
407
408     Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,particle->Vy(),particle->Vz(),mcphi,mceta,mcpt};
409     fResolHisto->Fill(vResolHisto);
410
411     Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mcsnp,mctgl,1./mcpt};
412     fPullHisto->Fill(vPullHisto);
413   //}
414   if(track) delete track;
415 }
416
417 //_____________________________________________________________________________
418 void AliPerformanceMC::ProcessInnerTPC(AliTrackReference* const refIn,  AliTrackReference *const  refOut, TParticle *const particle) {
419   //
420   // Test propagation from Out to In
421   //
422   if(!refIn) return;
423   if(!refOut) return;
424   if(!particle) return;
425
426   AliExternalTrackParam *track=MakeTrack(refOut,particle);
427   if (!track) return;
428
429   Double_t mass = particle->GetMass();
430   Double_t step=1;
431
432   Float_t alphaIn= TMath::ATan2(refIn->Y(),refIn->X());
433   //Float_t alphaOut= TMath::ATan2(refOut->Y(),refOut->X());
434   //track->Rotate(alphaIn-alphaOut);
435   track->Rotate(alphaIn);
436   //printf("alphaIn %f, alphaOut %f \n",alphaIn,alphaOut);
437
438   //
439   //Bool_t isOK = AliTracker::PropagateTrackTo(track, refIn->R(), mass, step, kTRUE,0.99);
440   //Bool_t isOK = AliTracker::PropagateTrackTo(track, refIn->R(), mass, step, kFALSE,0.99);
441   Bool_t isOK = AliTracker::PropagateTrackToBxByBz(track, refIn->R(), mass, step, kFALSE,0.99);
442   if(!isOK) return;
443
444   // calculate alpha angle
445   //Double_t xyz[3] = {refIn->X(),refIn->Y(),refIn0->Z()};
446   //Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
447   //if (alpha<0) alpha += TMath::TwoPi();
448
449   // rotate outer track to inner 
450   // and propagate to the radius of the first track referenco of TPC
451   //Double_t trRadius = TMath::Sqrt(xyz[1] * xyz[1] + xyz[0] * xyz[0]);
452   //Bool_t isOK = track->Propagate(alpha,trRadius,AliTracker::GetBz());
453   //if(!isOK) return;
454   
455   Float_t mceta =  -TMath::Log(TMath::Tan(0.5 * refIn->Theta()));
456   Float_t mcphi =  refIn->Phi();
457   if(mcphi<0) mcphi += 2.*TMath::Pi();
458   Float_t mcpt = refIn->Pt();
459   Float_t mcsnp = TMath::Sin(TMath::ATan2(refIn->Py(),refIn->Px()));
460   Float_t mctgl = TMath::Tan(TMath::ATan2(refIn->Pz(),refIn->Pt()));
461
462   Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC; 
463   Float_t pull1PtTPC=0., pullYTPC=0., pullZTPC=0., pullPhiTPC=0., pullLambdaTPC=0.; 
464
465   if(mcpt == 0) return;
466
467   //deltaYTPC = track->GetY()-refIn->Y();
468   deltaYTPC = track->GetY();
469   deltaZTPC = track->GetZ()-refIn->Z();
470   deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(refIn->Pz(),refIn->Pt());
471   deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(refIn->Py(),refIn->Px());
472   deltaPtTPC = (track->Pt()-mcpt) / mcpt;
473
474   /*
475   printf("------------------------------ \n");
476   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() ); 
477   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 ); 
478   */
479
480   if(TMath::Sqrt(track->GetSigmaY2())) pullYTPC = track->GetY() / TMath::Sqrt(track->GetSigmaY2());
481   if(TMath::Sqrt(track->GetSigmaZ2())) pullZTPC = (track->GetZ()-refIn->Z()) / TMath::Sqrt(track->GetSigmaZ2());
482   if(TMath::Sqrt(track->GetSigmaSnp2())) pullPhiTPC = (track->GetSnp() - mcsnp) / TMath::Sqrt(track->GetSigmaSnp2());
483   if(TMath::Sqrt(track->GetSigmaTgl2())) pullLambdaTPC = (track->GetTgl() - mctgl) / TMath::Sqrt(track->GetSigmaTgl2());
484   if(TMath::Sqrt(track->GetSigma1Pt2())) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
485
486   //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);
487
488   Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,refIn->Y(),refIn->Z(),mcphi,mceta,mcpt};
489   fResolHisto->Fill(vResolHisto);
490
491   Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,refIn->Y(),refIn->Z(),mcsnp,mctgl,1./mcpt};
492   fPullHisto->Fill(vPullHisto);
493
494   if(track) delete track;
495 }
496
497 //_____________________________________________________________________________
498 void AliPerformanceMC::ProcessOuterTPCExt(TParticle *const part, TClonesArray * const trefs)
499 {
500   const Int_t kMinRefs=5;
501   Int_t nrefs = trefs->GetEntries();
502   if (nrefs<kMinRefs) return; // we should have enough references
503   Int_t iref0 =-1;
504   Int_t iref1 =-1;
505   
506   for (Int_t iref=0; iref<nrefs; iref++){
507     AliTrackReference * ref = (AliTrackReference*)trefs->At(iref);
508     if (!ref) continue;    
509     Float_t dir = ref->X()*ref->Px()+ref->Y()*ref->Py();
510     if (dir<0) break;
511     if (ref->DetectorId()!=AliTrackReference::kTPC) continue;
512     if (iref0<0) iref0 = iref;
513     iref1 = iref;    
514   }
515   if (iref1-iref0<kMinRefs) return;
516   Double_t covar[15];
517   for (Int_t icov=0; icov<15; icov++) covar[icov]=0;
518   covar[0]=1; 
519   covar[2]=1; 
520   covar[5]=1;
521   covar[9]=1;
522   covar[14]=1;
523
524   AliTrackReference * refIn = (AliTrackReference*)trefs->At(iref0);
525   AliTrackReference * refOut = (AliTrackReference*)trefs->At(iref1);
526   AliExternalTrackParam *paramPropagate= MakeTrack(refIn,part);
527   AliExternalTrackParam *paramUpdate   = MakeTrack(refIn,part);
528   paramUpdate->AddCovariance(covar);
529
530   Bool_t isOKP=kTRUE;
531   Bool_t isOKU=kTRUE;
532   AliMagF * field = (AliMagF*) TGeoGlobalMagField::Instance()->GetField();
533   if(!field) return;
534   for (Int_t iref = iref0; iref<=iref1; iref++){
535     AliTrackReference * ref = (AliTrackReference*)trefs->At(iref);
536     if(!ref) continue;
537     Float_t alphaC= TMath::ATan2(ref->Y(),ref->X());
538     Double_t pos[3] = {ref->X(), ref->Y(), ref->Z()};
539     Double_t mag[3];
540     field->Field(pos,mag);
541     isOKP&=paramPropagate->Rotate(alphaC);
542     isOKU&=paramUpdate->Rotate(alphaC);
543     for (Float_t xref= paramPropagate->GetX(); xref<ref->R(); xref++){
544       //isOKP&=paramPropagate->PropagateTo(xref, mag[2]);
545       //isOKU&=paramUpdate->PropagateTo(xref, mag[2]);
546       isOKP&=paramPropagate->PropagateToBxByBz(xref, mag);
547       isOKU&=paramUpdate->PropagateToBxByBz(xref, mag);
548     }
549     //isOKP&=paramPropagate->PropagateTo(ref->R(), mag[2]);
550     //isOKU&=paramUpdate->PropagateTo(ref->R(), mag[2]);
551     isOKP&=paramPropagate->PropagateToBxByBz(ref->R(), mag);
552     isOKU&=paramUpdate->PropagateToBxByBz(ref->R(), mag);
553     Double_t clpos[2] = {0, ref->Z()};
554     Double_t clcov[3] = { 0.005,0,0.005};
555     isOKU&= paramUpdate->Update(clpos, clcov);  
556   }
557
558   Float_t mceta =  -TMath::Log(TMath::Tan(0.5 * refOut->Theta()));
559   Float_t mcphi =  refOut->Phi();
560   if(mcphi<0) mcphi += 2.*TMath::Pi();
561   Float_t mcpt = refOut->Pt();
562   Float_t mcsnp = TMath::Sin(TMath::ATan2(refOut->Py(),refOut->Px()));
563   Float_t mctgl = TMath::Tan(TMath::ATan2(refOut->Pz(),refOut->Pt()));
564
565   Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC; 
566   Float_t pull1PtTPC=0., pullYTPC=0., pullZTPC=0., pullPhiTPC=0., pullLambdaTPC=0.; 
567
568   if(mcpt == 0) return;
569
570   //deltaYTPC = track->GetY()-refIn->Y();
571   deltaYTPC = paramPropagate->GetY();
572   deltaZTPC = paramPropagate->GetZ()-refOut->Z();
573   deltaLambdaTPC = TMath::ATan2(paramPropagate->Pz(),paramPropagate->Pt())-TMath::ATan2(refOut->Pz(),refOut->Pt());
574   deltaPhiTPC = TMath::ATan2(paramPropagate->Py(),paramPropagate->Px())-TMath::ATan2(refOut->Py(),refOut->Px());
575   deltaPtTPC = (paramPropagate->Pt()-mcpt) / mcpt;
576
577   if(TMath::Sqrt(paramPropagate->GetSigmaY2())) pullYTPC = paramPropagate->GetY() / TMath::Sqrt(paramPropagate->GetSigmaY2());
578   if(TMath::Sqrt(paramPropagate->GetSigmaZ2())) pullZTPC = (paramPropagate->GetZ()-refOut->Z()) / TMath::Sqrt(paramPropagate->GetSigmaZ2());
579   if(TMath::Sqrt(paramPropagate->GetSigmaSnp2())) pullPhiTPC = (paramPropagate->GetSnp() - mcsnp) / TMath::Sqrt(paramPropagate->GetSigmaSnp2());
580   if(TMath::Sqrt(paramPropagate->GetSigmaTgl2())) pullLambdaTPC = (paramPropagate->GetTgl() - mctgl) / TMath::Sqrt(paramPropagate->GetSigmaTgl2());
581   if(TMath::Sqrt(paramPropagate->GetSigma1Pt2())) pull1PtTPC = (paramPropagate->OneOverPt()-1./mcpt) / TMath::Sqrt(paramPropagate->GetSigma1Pt2());
582
583   //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);
584
585   Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,refIn->Y(),refIn->Z(),mcphi,mceta,mcpt};
586   fResolHisto->Fill(vResolHisto);
587
588   Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,refIn->Y(),refIn->Z(),mcsnp,mctgl,1./mcpt};
589   fPullHisto->Fill(vPullHisto);
590 }
591  
592 //_____________________________________________________________________________
593 AliExternalTrackParam * AliPerformanceMC::MakeTrack(const AliTrackReference* const ref, TParticle* const part)
594 {
595   //
596   // Make track out of the track ref
597   // part - TParticle used to determine chargr
598   // the covariance matrix - equal 0 - starting from ideal MC position
599   Double_t xyz[3]={ref->X(),ref->Y(),ref->Z()};
600   Double_t pxyz[3]={ref->Px(),ref->Py(),ref->Pz()};
601   Double_t cv[21];
602   for (Int_t i=0; i<21;i++) cv[i]=0;
603   if (!part->GetPDG()) return 0;
604   AliExternalTrackParam * param = new AliExternalTrackParam(xyz,pxyz,cv,TMath::Nint(part->GetPDG()->Charge()/3.));
605   return param;
606 }
607
608 //_____________________________________________________________________________
609 TH1F* AliPerformanceMC::MakeResol(TH2F * his, Int_t integ, Bool_t type, Int_t cut){
610   // Create resolution histograms
611  
612   TH1F *hisr, *hism;
613   if (!gPad) new TCanvas;
614   hisr = AliTreeDraw::CreateResHistoII(his,&hism,integ,kTRUE,cut);
615   if (type) return hism;
616   else 
617     return hisr;
618 }
619
620 //_____________________________________________________________________________
621 void AliPerformanceMC::Analyse() {
622   // Analyse comparison information and store output histograms
623   // in the folder "folderRes"
624   //
625   TH1::AddDirectory(kFALSE);
626   TH1F *h=0;
627   TH2F *h2D=0;
628   TObjArray *aFolderObj = new TObjArray;
629   if(!aFolderObj) return;
630
631   // write results in the folder 
632   TCanvas * c = new TCanvas("Phi resol Tan","Phi resol Tan");
633   c->cd();
634
635   char name[256];
636   char title[256];
637
638   for(Int_t i=0; i<5; i++) 
639   {
640     for(Int_t j=5; j<10; j++) 
641     {
642       //if(j!=8) fResolHisto->GetAxis(8)->SetRangeUser(-0.9,0.89); // eta window
643       fResolHisto->GetAxis(9)->SetRangeUser(0.16,10.); // pt threshold
644       if(j!=8) fResolHisto->GetAxis(8)->SetRangeUser(0.,0.9); // eta window
645       //fResolHisto->GetAxis(9)->SetRangeUser(0.16,3.); // pt window
646       if(GetAnalysisMode() == 3) fResolHisto->GetAxis(5)->SetRangeUser(-80.,80.); // y range
647
648       h2D = (TH2F*)fResolHisto->Projection(i,j);
649
650       h = AliPerformanceMC::MakeResol(h2D,1,0,20);
651       snprintf(name,256,"h_res_%d_vs_%d",i,j);
652       h->SetName(name);
653
654       h->GetXaxis()->SetTitle(fResolHisto->GetAxis(j)->GetTitle());
655       snprintf(title,256,"%s %s",fResolHisto->GetAxis(i)->GetTitle(),"(resolution)");
656       h->GetYaxis()->SetTitle(title);
657       snprintf(title,256,"%s vs %s",title,fResolHisto->GetAxis(j)->GetTitle());
658       h->SetTitle(title);
659
660       if(j==9) h->SetBit(TH1::kLogX);    
661       aFolderObj->Add(h);
662
663       h = AliPerformanceMC::MakeResol(h2D,1,1,20);
664       //h = (TH1F*)arr->At(1);
665       snprintf(name,256,"h_mean_res_%d_vs_%d",i,j);
666       h->SetName(name);
667
668       h->GetXaxis()->SetTitle(fResolHisto->GetAxis(j)->GetTitle());
669       snprintf(title,256,"%s %s",fResolHisto->GetAxis(i)->GetTitle(),"(mean)");
670       h->GetYaxis()->SetTitle(title);
671
672       snprintf(title,256,"%s vs %s",title,fResolHisto->GetAxis(j)->GetTitle());
673       h->SetTitle(title);
674
675       if(j==9) h->SetBit(TH1::kLogX);    
676       aFolderObj->Add(h);
677
678       fResolHisto->GetAxis(8)->SetRangeUser(-1.5,1.5);
679       fResolHisto->GetAxis(9)->SetRangeUser(0.0,10.);
680
681       //
682       if(j!=8) fPullHisto->GetAxis(8)->SetRangeUser(-1.0,0.99); // theta window
683       else  fPullHisto->GetAxis(8)->SetRangeUser(-1.5,1.5); // theta window
684       fPullHisto->GetAxis(9)->SetRangeUser(0.0,10.);  // pt threshold
685       if(GetAnalysisMode() == 3) fPullHisto->GetAxis(5)->SetRangeUser(-80.,80.); // y range
686
687       h2D = (TH2F*)fPullHisto->Projection(i,j);
688
689       h = AliPerformanceMC::MakeResol(h2D,1,0,20);
690       snprintf(name,256,"h_pull_%d_vs_%d",i,j);
691       h->SetName(name);
692
693       h->GetXaxis()->SetTitle(fPullHisto->GetAxis(j)->GetTitle());
694       snprintf(title,256,"%s %s",fPullHisto->GetAxis(i)->GetTitle(),"(resolution)");
695       h->GetYaxis()->SetTitle(title);
696       snprintf(title,256,"%s vs %s",title,fPullHisto->GetAxis(j)->GetTitle());
697       h->SetTitle(title);
698
699       //if(j==9) h->SetBit(TH1::kLogX);    
700       aFolderObj->Add(h);
701
702       h = AliPerformanceMC::MakeResol(h2D,1,1,20);
703       snprintf(name,256,"h_mean_pull_%d_vs_%d",i,j);
704       h->SetName(name);
705
706       h->GetXaxis()->SetTitle(fPullHisto->GetAxis(j)->GetTitle());
707       snprintf(title,256,"%s %s",fPullHisto->GetAxis(i)->GetTitle(),"(mean)");
708       h->GetYaxis()->SetTitle(title);
709       snprintf(title,256,"%s vs %s",title,fPullHisto->GetAxis(j)->GetTitle());
710       h->SetTitle(title);
711
712       //if(j==9) h->SetBit(TH1::kLogX);    
713       aFolderObj->Add(h);
714     }
715   }
716
717   // export objects to analysis folder
718   fAnalysisFolder = ExportToFolder(aFolderObj);
719
720   // delete only TObjArray
721   if(aFolderObj) delete aFolderObj;
722 }
723
724 //_____________________________________________________________________________
725 TFolder* AliPerformanceMC::ExportToFolder(TObjArray * array) 
726 {
727   // recreate folder avery time and export objects to new one
728   //
729   AliPerformanceMC * comp=this;
730   TFolder *folder = comp->GetAnalysisFolder();
731
732   TString name, title;
733   TFolder *newFolder = 0;
734   Int_t i = 0;
735   Int_t size = array->GetSize();
736
737   if(folder) { 
738      // get name and title from old folder
739      name = folder->GetName();  
740      title = folder->GetTitle();  
741
742          // delete old one
743      delete folder;
744
745          // create new one
746      newFolder = CreateFolder(name.Data(),title.Data());
747      newFolder->SetOwner();
748
749          // add objects to folder
750      while(i < size) {
751            newFolder->Add(array->At(i));
752            i++;
753          }
754   }
755
756 return newFolder;
757 }
758
759 //_____________________________________________________________________________
760 Long64_t AliPerformanceMC::Merge(TCollection* const list) 
761 {
762   // Merge list of objects (needed by PROOF)
763
764   if (!list)
765   return 0;
766
767   if (list->IsEmpty())
768   return 1;
769
770   TIterator* iter = list->MakeIterator();
771   TObject* obj = 0;
772
773   // collection of generated histograms
774   Int_t count=0;
775   while((obj = iter->Next()) != 0) 
776   {
777   AliPerformanceMC* entry = dynamic_cast<AliPerformanceMC*>(obj);
778   if (entry == 0) continue; 
779
780   fResolHisto->Add(entry->fResolHisto);
781   fPullHisto->Add(entry->fPullHisto);
782
783   count++;
784   }
785
786 return count;
787 }
788
789 //_____________________________________________________________________________
790 TFolder* AliPerformanceMC::CreateFolder(TString name,TString title) { 
791 // create folder for analysed histograms
792 //
793 TFolder *folder = 0;
794   folder = new TFolder(name.Data(),title.Data());
795
796   return folder;
797 }