bc333e93fd065c679f18b1153f661ce321ba79bc
[u/mrichter/AliRoot.git] / STEER / AliESDComparisonMI.C
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16
17 ///////////////////////////////////////////////////////////////////////////////
18 //                                                                           //
19 //  Time Projection Chamber                                                  //
20 //  Comparison macro for ESD                                                 //
21 //  responsible: 
22 //  marian.ivanov@cern.ch                                                    //
23 //
24 //
25
26 /*
27 marian.ivanov@cern.ch
28 Usage:
29
30
31 .L $ALICE_ROOT/STEER/AliGenInfo.C+
32 //be sure you created genTracks file before
33 .L $ALICE_ROOT/STEER/AliESDComparisonMI.C+
34
35 //
36 ESDCmpTr *t2 = new ESDCmpTr("genTracks.root","cmpESDTracks.root","galice.root",-1,1,0);
37 t2->Exec();
38
39 //
40 //some cuts definition
41 TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.00005&&abs(MC.fVDist[2])<0.0005")
42 //TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.5&&abs(MC.fVDist[2])<0.5")
43 TCut citsin("citsin","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<3.9");
44
45 TCut csec("csec","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)>0.5");
46
47
48 TCut crec("crec","fReconstructed==1");
49 TCut cteta1("cteta1","abs(MC.fTrackRef.Theta()/3.1415-0.5)<0.25");
50 TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1");
51 TCut csens("csens","abs(sqrt(fVDist[0]**2+fVDist[1]**2)-170)<50");
52 TCut cmuon("cmuon","abs(MC.fParticle.fPdgCode==-13)");
53 TCut cchi2("cchi2","fESDTrack.fITSchi2MIP[0]<7.&&fESDTrack.fITSchi2MIP[1]<5.&&fESDTrack.fITSchi2MIP[2]<7.&&fESDTrack.fITSchi2MIP[3]<7.5&&fESDTrack.fITSchi2MIP[4]<6.")
54
55 AliESDComparisonDraw comp;  
56 comp.SetIO(); 
57
58 //
59 //example
60 comp.fTree->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)");
61
62 TH1F his("his","his",100,0,20);
63 TH1F hpools("hpools","hpools",100,-7,7);
64 TH1F hfake("hfake","hfake",1000,0,150);
65 TProfile profp0("profp0","profp0",20,-0.4,0.9)
66
67 comp.DrawLogXY("fTPCinP0[3]","fTPCDelta[4]/fTPCinP1[3]","fReconstructed==1&&abs(fPdg)==211"+cprim,"1",4,0.2,1.5,-0.06,0.06)
68 comp.fRes->Draw();
69 comp.fMean->Draw();  
70 comp.Eff("fTPCinP0[3]","fRowsWithDigits>120&&abs(fPdg)==211"+cteta1+cpos1+cprim,"fTPCOn",20,0.2,1.5)
71 comp.fRes->Draw();
72
73 comp.fTree->Draw("fESDTrack.fITSsignal/fESDTrack.fTPCsignal","fITSOn&&fTPCOn&&fESDTrack.fITSFakeRatio==0") 
74
75 TH1F his("his","his",100,0,20);
76 TH1F hpools("hpools","hpools",100,-7,7);
77
78 TH2F * hdedx0 = new TH2F("dEdx0","dEdx0",100, 0,2,200,0,250); hdedx0->SetMarkerColor(2);
79 TH2F * hdedx1 = new TH2F("dEdx1","dEdx1",100, 0,2,200,0,250); hdedx1->SetMarkerColor(3);
80 TH2F * hdedx2 = new TH2F("dEdx2","dEdx2",100, 0,2,200,0,250); hdedx2->SetMarkerColor(4);
81 TH2F * hdedx3 = new TH2F("dEdx3","dEdx3",100, 0,2,200,0,250); hdedx3->SetMarkerColor(5);
82 comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx0","fITSOn&&MC.fParticle.P()<2&&abs(fPdg)==211&&fITStrack.fN==6"+cprim) 
83 comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx1","fITSOn&&MC.fParticle.P()<2&&abs(fPdg)==2212&&fITStrack.fN==6"+cprim) 
84 comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx2","fITSOn&&MC.fParticle.P()<2&&abs(fPdg)==321&&fITStrack.fN==6"+cprim) 
85 comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx3","fITSOn&&MC.fParticle.P()<2&&abs(fPdg)==11&&fITStrack.fN==6"+cprim) 
86 hdedx0->Draw(); hdedx1->Draw("same"); hdedx2->Draw("same"); hdedx3->Draw("same");
87
88 comp.DrawXY("fITSinP0[3]","fITSPools[4]","fReconstructed==1&&fPdg==-211&&fITSOn"+cprim,"1",4,0.2,1.0,-8,8)
89
90 */
91
92
93 #if !defined(__CINT__) || defined(__MAKECINT__)
94 #include <stdio.h>
95 #include <string.h>
96 //ROOT includes
97 #include "Rtypes.h"
98 #include "TFile.h"
99 #include "TTree.h"
100 #include "TChain.h"
101 #include "TCut.h"
102 #include "TString.h"
103 #include "TBenchmark.h"
104 #include "TStopwatch.h"
105 #include "TParticle.h"
106 #include "TSystem.h"
107 #include "TTimer.h"
108 #include "TVector3.h"
109 #include "TPad.h"
110 #include "TCanvas.h"
111 #include "TH1F.h"
112 #include "TH2F.h"
113 #include "TF1.h"
114 #include "TText.h"
115 #include "Getline.h"
116 #include "TStyle.h"
117
118 //ALIROOT includes
119 #include "AliRun.h"
120 #include "AliStack.h"
121 #include "AliESDtrack.h"
122 #include "AliSimDigits.h"
123 #include "AliTPCParam.h"
124 #include "AliTPC.h"
125 #include "AliTPCLoader.h"
126 #include "AliDetector.h"
127 #include "AliTrackReference.h"
128 #include "AliRun.h"
129 #include "AliTPCParamSR.h"
130 #include "AliTracker.h"
131 #include "AliComplexCluster.h"
132 #include "AliMagF.h"
133 #include "AliESD.h"
134 #include "AliESDtrack.h"
135 #include "AliITStrackV2.h"
136 #include "AliTRDtrack.h"
137 #endif
138 #include "AliGenInfo.h"
139 #include "AliESDComparisonMI.h"
140
141
142
143 //
144 //
145 void AliESDRecInfo::Update(AliMCInfo* info,AliTPCParam * par, Bool_t reconstructed)
146 {
147   //
148   //
149   //calculates derived variables
150   //  
151   //
152   AliTrackReference * ref = &(info->fTrackRef);
153   fTPCinR0[0] = info->fTrackRef.X();    
154   fTPCinR0[1] = info->fTrackRef.Y();    
155   fTPCinR0[2] = info->fTrackRef.Z();
156   fTPCinR0[3] = TMath::Sqrt(fTPCinR0[0]*fTPCinR0[0]+fTPCinR0[1]*fTPCinR0[1]);
157   fTPCinR0[4] = TMath::ATan2(fTPCinR0[1],fTPCinR0[0]);
158   //
159   fTPCinP0[0] = ref->Px();
160   fTPCinP0[1] = ref->Py();
161   fTPCinP0[2] = ref->Pz();
162   fTPCinP0[3] = ref->Pt();
163   fTPCinP0[4] = ref->P();
164   fDeltaP     = (ref->P()-info->fParticle.P())/info->fParticle.P();
165   //
166   //
167   if (fTPCinP0[3]>0.0000001){
168     //
169     fTPCAngle0[0] = TMath::ATan2(fTPCinP0[1],fTPCinP0[0]);
170     fTPCAngle0[1] = TMath::ATan(fTPCinP0[2]/fTPCinP0[3]);
171   }
172   //
173   //
174   fITSinP0[0]=info->fParticle.Px();
175   fITSinP0[1]=info->fParticle.Py();
176   fITSinP0[2]=info->fParticle.Pz();
177   fITSinP0[3]=info->fParticle.Pt();    
178   //
179   fITSinR0[0]=info->fParticle.Vx();
180   fITSinR0[1]=info->fParticle.Vy();
181   fITSinR0[2]=info->fParticle.Vz();
182   fITSinR0[3] = TMath::Sqrt(fITSinR0[0]*fITSinR0[0]+fITSinR0[1]*fITSinR0[1]);
183   fITSinR0[4] = TMath::ATan2(fITSinR0[1],fITSinR0[0]);
184   //
185   //
186   if (fITSinP0[3]>0.0000001){
187     fITSAngle0[0] = TMath::ATan2(fITSinP0[1],fITSinP0[0]);
188     fITSAngle0[1] = TMath::ATan(fITSinP0[2]/fITSinP0[3]);
189   }
190   //
191   if (reconstructed==kFALSE){
192     fReconstructed = kFALSE;
193     fTPCOn = kFALSE;
194     fITSOn = kFALSE;
195     fTRDOn = kFALSE;
196     return;
197   }
198   fReconstructed = kTRUE;
199   fTPCOn = ((fESDTrack.GetStatus()&AliESDtrack::kTPCrefit)>0) ? kTRUE : kFALSE;
200   fITSOn = ((fESDTrack.GetStatus()&AliESDtrack::kITSrefit)>0) ? kTRUE : kFALSE;
201   fTRDOn = ((fESDTrack.GetStatus()&AliESDtrack::kTRDrefit)>0) ? kTRUE : kFALSE;
202   
203   if (fTPCOn){
204     //TPC
205     fESDTrack.GetInnerXYZ(fTPCinR1);
206     fTPCinR1[3] = TMath::Sqrt(fTPCinR1[0]*fTPCinR1[0]+fTPCinR1[1]*fTPCinR1[1]);
207     fTPCinR1[4] = TMath::ATan2(fTPCinR1[1],fTPCinR1[0]);        
208     fESDTrack.GetInnerPxPyPz(fTPCinP1);
209     fTPCinP1[3] = TMath::Sqrt(fTPCinP1[0]*fTPCinP1[0]+fTPCinP1[1]*fTPCinP1[1]);
210     fTPCinP1[4] = TMath::Sqrt(fTPCinP1[3]*fTPCinP1[3]+fTPCinP1[2]*fTPCinP1[2]);
211     //
212     //
213     if (fTPCinP1[3]>0.0000001){
214       fTPCAngle1[0] = TMath::ATan2(fTPCinP1[1],fTPCinP1[0]);
215       fTPCAngle1[1] = TMath::ATan(fTPCinP1[2]/fTPCinP1[3]);  
216     }    
217     Double_t cov[15], param[5],x;
218     fESDTrack.GetInnerExternalCovariance(cov);
219     fESDTrack.GetInnerExternalParameters(x,param);    
220     //
221     fTPCDelta[0] = (fTPCinR0[4]-fTPCinR1[4])*fTPCinR1[3];  //delta rfi
222     fTPCPools[0] = fTPCDelta[0]/TMath::Sqrt(cov[0]);
223     fTPCDelta[1] = (fTPCinR0[2]-fTPCinR1[2]);              //delta z
224     fTPCPools[1] = fTPCDelta[1]/TMath::Sqrt(cov[2]);
225     fTPCDelta[2] = (fTPCAngle0[0]-fTPCAngle1[0]);
226     fTPCPools[2] = fTPCDelta[2]/TMath::Sqrt(cov[5]);
227     fTPCDelta[3] = (TMath::Tan(fTPCAngle0[1])-TMath::Tan(fTPCAngle1[1]));
228     fTPCPools[3] = fTPCDelta[3]/TMath::Sqrt(cov[9]);
229     fTPCDelta[4] = (fTPCinP0[3]-fTPCinP1[3]);
230     Double_t sign = (param[4]>0)? 1.:-1; 
231     fTPCPools[4] = sign*(1./fTPCinP0[3]-1./fTPCinP1[3])/TMath::Sqrt(cov[14]);    
232   }
233   if (fITSOn){
234     // ITS 
235     Double_t param[5],x;
236     //fESDTrack.GetExternalParameters(x,param);   
237     fESDTrack.GetConstrainedExternalParameters(x,param);   
238     Double_t cov[15];
239     fESDTrack.GetExternalCovariance(cov);
240     fESDTrack.GetConstrainedExternalCovariance(cov);
241     if (TMath::Abs(param[4])<0.0000000001) return;
242
243     fESDTrack.GetXYZ(fITSinR1);
244     fESDTrack.GetPxPyPz(fITSinP1);
245     fITSinP1[3] = TMath::Sqrt(fITSinP1[0]*fITSinP1[0]+fITSinP1[1]*fITSinP1[1]);
246     //
247     fITSinR1[3] = TMath::Sqrt(fITSinR1[0]*fITSinR1[0]+fITSinR1[1]*fITSinR1[1]);
248     fITSinR1[4] = TMath::ATan2(fITSinR1[1],fITSinR1[0]);
249     //
250     //
251     if (fITSinP1[3]>0.0000001){
252       fITSAngle1[0] = TMath::ATan2(fITSinP1[1],fITSinP1[0]);
253       fITSAngle1[1] = TMath::ATan(fITSinP1[2]/fITSinP1[3]);  
254     }
255     //
256     //
257     fITSDelta[0] = (fITSinR0[4]-fITSinR1[4])*fITSinR1[3];  //delta rfi
258     fITSPools[0] = fITSDelta[0]/TMath::Sqrt(cov[0]);
259     fITSDelta[1] = (fITSinR0[2]-fITSinR1[2]);              //delta z
260     fITSPools[1] = fITSDelta[1]/TMath::Sqrt(cov[2]);
261     fITSDelta[2] = (fITSAngle0[0]-fITSAngle1[0]);
262     fITSPools[2] = fITSDelta[2]/TMath::Sqrt(cov[5]);
263     fITSDelta[3] = (TMath::Tan(fITSAngle0[1])-TMath::Tan(fITSAngle1[1]));
264     fITSPools[3] = fITSDelta[3]/TMath::Sqrt(cov[9]);
265     fITSDelta[4] = (fITSinP0[3]-fITSinP1[3]);    
266     Double_t sign = (param[4]>0) ? 1:-1; 
267     fITSPools[4] = sign*(1./fITSinP0[3]-1./fITSinP1[3])/TMath::Sqrt(cov[14]);    
268   }
269   
270 }
271
272
273 void  AliESDRecV0Info::Update(Float_t vertex[3])
274 {
275   Float_t v[3] = {fXr[0]-vertex[0],fXr[1]-vertex[1],fXr[2]-vertex[2]};
276   Float_t p[3] = {fPdr[0]+fPm[0], fPdr[1]+fPm[1],fPdr[2]+fPm[2]};
277
278   Float_t vnorm2 = v[0]*v[0]+v[1]*v[1];
279   Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
280   vnorm2 = TMath::Sqrt(vnorm2);
281   Float_t pnorm2 = p[0]*p[0]+p[1]*p[1];
282   Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
283   pnorm2 = TMath::Sqrt(pnorm2);
284   
285   fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
286   fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);  
287   fPointAngle   = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
288
289 }
290
291
292 ////////////////////////////////////////////////////////////////////////
293 ESDCmpTr::ESDCmpTr()
294 {
295   Reset();
296 }
297
298 ////////////////////////////////////////////////////////////////////////
299 ESDCmpTr::ESDCmpTr(const char* fnGenTracks,
300                    const char* fnCmp,
301                    const char* fnGalice, Int_t direction,
302                    Int_t nEvents, Int_t firstEvent)
303 {
304   Reset();
305   //  fFnGenTracks = fnGenTracks;
306   //  fFnCmp = fnCmp;
307   sprintf(fFnGenTracks,"%s",fnGenTracks);
308   sprintf(fFnCmp,"%s",fnCmp);
309
310   fFirstEventNr = firstEvent;
311   fEventNr = firstEvent;
312   fNEvents = nEvents;
313   fDirection = direction;
314   //
315   fLoader = AliRunLoader::Open(fnGalice);
316   if (gAlice){
317     //delete gAlice->GetRunLoader();
318     delete gAlice;
319     gAlice = 0x0;
320   }
321   if (fLoader->LoadgAlice()){
322     cerr<<"Error occured while l"<<endl;
323   }
324   Int_t nall = fLoader->GetNumberOfEvents();
325   if (nEvents==0) {
326     nEvents =nall;
327     fNEvents=nall;
328     fFirstEventNr=0;
329   }    
330
331   if (nall<=0){
332     cerr<<"no events available"<<endl;
333     fEventNr = 0;
334     return;
335   }
336   if (firstEvent+nEvents>nall) {
337     fEventNr = nall-firstEvent;
338     cerr<<"restricted number of events availaible"<<endl;
339   }
340   AliMagF * magf = gAlice->Field();
341   AliTracker::SetFieldMap(magf);
342
343 }
344
345
346 ////////////////////////////////////////////////////////////////////////
347 ESDCmpTr::~ESDCmpTr()
348 {
349   if (fLoader) {
350     delete fLoader;
351   }
352 }
353
354 //////////////////////////////////////////////////////////////
355 Int_t ESDCmpTr::SetIO()
356 {
357   //
358   // 
359   CreateTreeCmp();
360   if (!fTreeCmp) return 1;
361   fParamTPC = GetTPCParam();
362   //
363   if (!ConnectGenTree()) {
364     cerr<<"Cannot connect tree with generated tracks"<<endl;
365     return 1;
366   }
367   return 0;
368 }
369
370 //////////////////////////////////////////////////////////////
371
372 Int_t ESDCmpTr::SetIO(Int_t eventNr)
373 {
374   //
375   // 
376   // SET INPUT
377   //
378   TFile f("AliESDs.root");
379   //
380  
381   TTree* tree = (TTree*) f.Get("esdTree");
382   if (!tree) { 
383     Char_t ename[100]; 
384     sprintf(ename,"%d",eventNr);
385     fEvent = (AliESD*)f.Get(ename);
386     if (!fEvent){
387       sprintf(ename,"ESD%d",eventNr);
388       fEvent = (AliESD*)f.Get(ename);
389     }
390   }
391   else{
392     tree->SetBranchAddress("ESD", &fEvent);
393     tree->GetEntry(eventNr);
394   }
395
396
397   /*
398   Char_t ename[100]; 
399   sprintf(ename,"%d",eventNr);
400   fEvent = (AliESD*)f.Get(ename);
401   if (!fEvent){
402     sprintf(ename,"ESD%d",eventNr);
403     fEvent = (AliESD*)f.Get(ename);
404   }
405   
406   TTree* tree = (TTree*) f.Get("esdTree");
407   if (!tree) {
408     Error("CheckESD", "no ESD tree found");
409     return kFALSE;
410   }
411   tree->SetBranchAddress("ESD", &fEvent);
412   tree->GetEntry(eventNr);
413   */
414
415   if (!fEvent) return 1;
416
417   return 0;
418 }
419
420
421
422 ////////////////////////////////////////////////////////////////////////
423 void ESDCmpTr::Reset()
424 {
425   fEventNr = 0;
426   fNEvents = 0;
427   fTreeCmp = 0;
428   //  fFnCmp = "cmpTracks.root";
429   fFileGenTracks = 0;
430   fDebug = 0;
431   //
432   fParamTPC = 0;
433   fEvent =0;
434 }
435
436 ////////////////////////////////////////////////////////////////////////
437 Int_t ESDCmpTr::Exec(Int_t nEvents, Int_t firstEventNr)
438 {
439   fNEvents = nEvents;
440   fFirstEventNr = firstEventNr;
441   return Exec();
442 }
443
444 ////////////////////////////////////////////////////////////////////////
445 Int_t ESDCmpTr::Exec()
446 {
447   TStopwatch timer;
448   timer.Start();
449
450   if (SetIO()==1) 
451     return 1;
452    
453   fNextTreeGenEntryToRead = 0;
454   cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl;
455   for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents;
456        eventNr++) {
457     fEventNr = eventNr;
458     SetIO(fEventNr);
459     fNParticles = gAlice->GetEvent(fEventNr);    
460
461     fIndexRecTracks = new Int_t[fNParticles*4];  //write at maximum 4 tracks corresponding to particle
462     fFakeRecTracks = new Int_t[fNParticles];
463     fMultiRecTracks = new Int_t[fNParticles];
464     for (Int_t i = 0; i<fNParticles; i++) {
465       fIndexRecTracks[4*i] = -1;
466       fIndexRecTracks[4*i+1] = -1;
467       fIndexRecTracks[4*i+2] = -1;
468       fIndexRecTracks[4*i+3] = -1;
469
470       fFakeRecTracks[i] = 0;
471       fMultiRecTracks[i] = 0;
472     }
473   
474     cout<<"Start to process event "<<fEventNr<<endl;
475     cout<<"\tfNParticles = "<<fNParticles<<endl;
476     if (fDebug>2) cout<<"\tStart loop over TreeT"<<endl;
477     if (TreeTLoop()>0) return 1;
478
479     if (fDebug>2) cout<<"\tStart loop over tree genTracks"<<endl;
480     if (TreeGenLoop(eventNr)>0) return 1;
481     if (fDebug>2) cout<<"\tEnd loop over tree genTracks"<<endl;
482
483     delete [] fIndexRecTracks;
484     delete [] fFakeRecTracks;
485     delete [] fMultiRecTracks;
486   }
487
488   CloseOutputFile();
489
490   cerr<<"Exec finished"<<endl;
491   timer.Stop();
492   timer.Print();
493   return 0;
494
495 }
496 ////////////////////////////////////////////////////////////////////////
497 Bool_t ESDCmpTr::ConnectGenTree()
498 {
499 //
500 // connect all branches from the genTracksTree
501 // use the same variables as for the new cmp tree, it may work
502 //
503   fFileGenTracks = TFile::Open(fFnGenTracks,"READ");
504   if (!fFileGenTracks) {
505     cerr<<"Error in ConnectGenTree: cannot open file "<<fFnGenTracks<<endl;
506     return kFALSE;
507   }
508   fTreeGenTracks = (TTree*)fFileGenTracks->Get("genTracksTree");
509   if (!fTreeGenTracks) {
510     cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
511         <<fFnGenTracks<<endl;
512     return kFALSE;
513   }
514   //
515   fMCInfo = new AliMCInfo;
516   fTreeGenTracks->SetBranchAddress("MC",&fMCInfo);
517   //
518   if (fDebug > 1) {
519     cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl;
520   }
521   return kTRUE;
522 }
523
524
525 ////////////////////////////////////////////////////////////////////////
526 void ESDCmpTr::CreateTreeCmp() 
527 {
528   fFileCmp = TFile::Open(fFnCmp,"RECREATE");
529   if (!fFileCmp) {
530     cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl;
531     return;
532   }
533
534
535   fTreeCmp    = new TTree("ESDcmpTracks","ESDcmpTracks");
536   //
537   fMCInfo = new AliMCInfo;
538   fRecInfo = new AliESDRecInfo;
539   //
540   AliESDtrack * esdTrack = new AliESDtrack;
541   AliITStrackV2 * itsTrack = new AliITStrackV2;  
542   //
543   fTreeCmp->Branch("MC","AliMCInfo",&fMCInfo);
544   fTreeCmp->Branch("RC","AliESDRecInfo",&fRecInfo);
545   fTreeCmp->Branch("fESDTrack","AliESDtrack",&esdTrack);
546   //  fTreeCmp->Branch("ITS","AliITStrackV2",&itsTrack);
547   delete esdTrack;
548   //
549   fTreeCmp->AutoSave(); 
550
551 }
552 ////////////////////////////////////////////////////////////////////////
553 void ESDCmpTr::CloseOutputFile()  
554 {
555   if (!fFileCmp) {
556     cerr<<"File "<<fFnCmp<<" not found as an open file."<<endl;
557     return;
558   }
559   fFileCmp->cd();
560   fTreeCmp->Write();    
561   delete fTreeCmp;
562   
563   fFileCmp->Close();
564   delete fFileCmp;
565   return;
566 }
567 ////////////////////////////////////////////////////////////////////////
568
569 TVector3 ESDCmpTr::TR2Local(AliTrackReference *trackRef,
570                             AliTPCParam *paramTPC) {
571
572   Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
573   Int_t index[4];
574   paramTPC->Transform0to1(x,index);
575   paramTPC->Transform1to2(x,index);
576   return TVector3(x);
577 }
578 ////////////////////////////////////////////////////////////////////////
579
580 Int_t ESDCmpTr::TreeTLoop()
581 {
582   //
583   // loop over all ESD reconstructed tracks and store info in memory
584   //
585   TStopwatch  timer;
586   timer.Start();
587   //  
588   Int_t nEntries = (Int_t)fEvent->GetNumberOfTracks();  
589   //
590   fTracks      = new TObjArray(nEntries);
591   //
592   //load tracks to the memory
593   for (Int_t i=0; i<nEntries;i++){
594     AliESDtrack * track =fEvent->GetTrack(i);
595     fTracks->AddAt(track,i);
596   }
597   //
598   //
599   AliESDtrack * track=0;
600   for (Int_t iEntry=0; iEntry<nEntries;iEntry++){
601     track = (AliESDtrack*)fTracks->UncheckedAt(iEntry);
602     //
603     Int_t label = track->GetLabel();
604     Int_t absLabel = abs(label);
605     if (absLabel < fNParticles) {
606       //      fIndexRecTracks[absLabel] =  iEntry;
607       if (label < 0) fFakeRecTracks[absLabel]++;      
608       if (fMultiRecTracks[absLabel]>0){
609         if (fMultiRecTracks[absLabel]<4)
610           fIndexRecTracks[absLabel*4+fMultiRecTracks[absLabel]] =  iEntry;      
611       }
612       else      
613         fIndexRecTracks[absLabel*4] =  iEntry;
614       fMultiRecTracks[absLabel]++;
615     }
616   }  
617   printf("Time spended in TreeTLoop\n");
618   timer.Print();
619   
620   if (fDebug > 2) cerr<<"end of TreeTLoop"<<endl;  
621   return 0;
622 }
623
624 ////////////////////////////////////////////////////////////////////////
625 Int_t ESDCmpTr::TreeGenLoop(Int_t eventNr)
626 {
627 //
628 // loop over all entries for a given event, find corresponding 
629 // rec. track and store in the fTreeCmp
630 //
631   TStopwatch timer;
632   timer.Start();
633   Int_t entry = fNextTreeGenEntryToRead;
634   Double_t nParticlesTR = fTreeGenTracks->GetEntriesFast();
635   cerr<<"fNParticles, nParticlesTR, fNextTreeGenEntryToRead: "<<fNParticles<<" "
636       <<nParticlesTR<<" "<<fNextTreeGenEntryToRead<<endl;
637   TBranch * branch = fTreeCmp->GetBranch("RC");
638   branch->SetAddress(&fRecInfo); // set all pointers
639
640   while (entry < nParticlesTR) {
641     fTreeGenTracks->GetEntry(entry);
642     entry++;
643     if (eventNr < fMCInfo->fEventNr) continue;
644     if (eventNr > fMCInfo->fEventNr) continue;;
645     //
646     fNextTreeGenEntryToRead = entry-1;
647     if (fDebug > 2 && fMCInfo->fLabel < 10) {
648       cerr<<"Fill track with a label "<<fMCInfo->fLabel<<endl;
649     }
650     //
651     fRecInfo->Reset();
652     AliESDtrack * track=0;
653     fRecInfo->fReconstructed =0;
654     TVector3 local = TR2Local(&(fMCInfo->fTrackRef),fParamTPC);
655     local.GetXYZ(fRecInfo->fTRLocalCoord);      
656     //
657     if (fIndexRecTracks[fMCInfo->fLabel*4] >= 0) {
658       track= (AliESDtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4]);
659       //
660       //
661       // find nearest track if multifound
662       if (fIndexRecTracks[fMCInfo->fLabel*4]+1){
663         Double_t xyz[3];
664         track->GetInnerXYZ(xyz);
665         Float_t dz = TMath::Abs(local.Z()-xyz[2]);
666         for (Int_t i=1;i<4;i++){
667           if (fIndexRecTracks[fMCInfo->fLabel*4+i]>=0){
668             AliESDtrack * track2 = (AliESDtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4+i]);
669             track2->GetInnerXYZ(xyz);
670             Float_t dz2 = TMath::Abs(local.Z()-xyz[2]);
671             if  (TMath::Abs(dz2)<dz)
672               track = track2;              
673           }
674         }
675       } 
676       //
677       fRecInfo->fESDTrack =*track; 
678       if (track->GetITStrack())
679         fRecInfo->fITStrack = *((AliITStrackV2*)track->GetITStrack());
680       if (track->GetTRDtrack()){
681         fRecInfo->fTRDtrack = *((AliTRDtrack*)track->GetTRDtrack());
682       }
683       else{
684         fRecInfo->fTRDtrack.SetdEdx(-1);
685       }
686       fRecInfo->fReconstructed = 1;
687       fRecInfo->fFake     = fFakeRecTracks[fMCInfo->fLabel];
688       fRecInfo->fMultiple = fMultiRecTracks[fMCInfo->fLabel];
689       //
690       fRecInfo->Update(fMCInfo,fParamTPC,kTRUE);          
691     }
692     else{
693       fRecInfo->Update(fMCInfo,fParamTPC,kFALSE);
694     }
695
696     fTreeCmp->Fill();
697   }
698   fTreeCmp->AutoSave();
699   fTracks->Delete();
700   printf("Time spended in TreeGenLoop\n");
701   timer.Print();
702   if (fDebug > 2) cerr<<"end of TreeGenLoop"<<endl;
703
704   return 0;
705 }
706 ////////////////////////////////////////////////////////////////////////
707 ////////////////////////////////////////////////////////////////////////
708
709 void AliESDComparisonDraw::SetIO(const char *fname)
710 {
711   //
712    TFile* file = TFile::Open(fname);
713    //
714    fTree = (TTree*) file->Get("ESDcmpTracks");
715    if (!fTree) {
716     printf("no track comparison tree found\n");
717     file->Close();
718     delete file;
719   }
720 }
721
722