Additional protection and change in the interface (Yu.Belikov)
[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,0,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.01&&abs(MC.fVDist[2])<0.01")
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 TCut citsin("citsin","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<5");
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.fParticle.Theta()/3.1415-0.5)<0.25");
50 TCut cteta05("cteta05","abs(MC.fParticle.Theta()/3.1415-0.5)<0.1");
51
52 TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1");
53 TCut csens("csens","abs(sqrt(fVDist[0]**2+fVDist[1]**2)-170)<50");
54 TCut cmuon("cmuon","abs(MC.fParticle.fPdgCode==-13)");
55 TCut cchi2("cchi2","fESDTrack.fITSchi2MIP[0]<7.&&fESDTrack.fITSchi2MIP[1]<5.&&fESDTrack.fITSchi2MIP[2]<7.&&fESDTrack.fITSchi2MIP[3]<7.5&&fESDTrack.fITSchi2MIP[4]<6.")
56
57 AliESDComparisonDraw comp;  
58 comp.SetIO(); 
59 TFile f("genHits.root");
60 TTree * treel = (TTree*)f.Get("HitLines");
61 if (treel) comp->fTree->AddFriend(treel,"L");
62
63 //
64 //example
65 comp.fTree->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)");
66 comp.fTree->SetAlias("direction","MC.fParticle.fVx*MC.fParticle.fPx+MC.fParticle.fVy*MC.fParticle.fPy");
67 comp.fTree->SetAlias("decaydir","MC.fTRdecay.fX*MC.fTRdecay.fPx+MC.fTRdecay.fY*MC.fTRdecay.fPy");
68 comp.fTree->SetAlias("theta","MC.fTrackRef.Theta()");
69 comp.fTree->SetAlias("primdca","sqrt(RC.fITStrack.fD[0]**2+RC.fITStrack.fD[1]**2)");
70 comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
71 comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
72
73
74 TH1F his("his","his",100,0,20);
75 TH1F hpools("hpools","hpools",100,-7,7);
76 TH1F hfake("hfake","hfake",1000,0,150);
77 TProfile profp0("profp0","profp0",20,-0.4,0.9)
78
79 comp.DrawXY("fTPCinP0[3]","fTPCDelta[4]/fTPCinP1[3]","fReconstructed==1"+cprim,"1",4,0.2,1.5,-0.06,0.06)
80 comp.fRes->Draw();
81 comp.fMean->Draw();  
82
83 comp.DrawXY("fITSinP0[3]","fITSDelta[4]/fITSinP1[3]","fReconstructed==1&&fITSOn"+cprim,"1",4,0.2,1.5,-0.06,0.06)
84 comp.fRes->Draw();
85
86 comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn",20,0.2,1.5)
87 comp.fRes->Draw();
88
89 comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn&&fITSOn&&fESDTrack.fITSFakeRatio<0.1",10,0.2,1.5)
90 comp.fRes->Draw();
91 comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1+cprim,"fTPCOn&&fITSOn&&fESDTrack.fITSFakeRatio>0.1",10,0.2,1.5)
92 comp.fRes->Draw();
93
94 comp.fTree->Draw("fESDTrack.fITSsignal/fESDTrack.fTPCsignal","fITSOn&&fTPCOn&&fESDTrack.fITSFakeRatio==0") 
95
96 TH1F his("his","his",100,0,20);
97 TH1F hpools("hpools","hpools",100,-7,7);
98
99 TH2F * hdedx0 = new TH2F("dEdx0","dEdx0",100, 0,2,200,0,550); hdedx0->SetMarkerColor(1);
100 TH2F * hdedx1 = new TH2F("dEdx1","dEdx1",100, 0,2,200,0,550); hdedx1->SetMarkerColor(4);
101 TH2F * hdedx2 = new TH2F("dEdx2","dEdx2",100, 0,2,200,0,550); hdedx2->SetMarkerColor(3);
102 TH2F * hdedx3 = new TH2F("dEdx3","dEdx3",100, 0,2,200,0,550); hdedx3->SetMarkerColor(2);
103
104 comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx0","fITSOn&&abs(fPdg)==211&&fITStrack.fN==6"+cprim) 
105 comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx1","fITSOn&&abs(fPdg)==2212&&fITStrack.fN==6"+cprim) 
106 comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx2","fITSOn&&abs(fPdg)==321&&fITStrack.fN==6"+cprim) 
107 comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx3","fITSOn&&abs(fPdg)==11&&fITStrack.fN==6"+cprim) 
108
109
110 comp.fTree->Draw("fESDTrack.fTRDsignal:MC.fParticle.P()>>dEdx0","fTRDOn&&abs(fPdg)==211&&fTRDtrack.fN>40&&fStatus[2]>1") 
111 comp.fTree->Draw("fESDTrack.fTRDsignal:MC.fParticle.P()>>dEdx1","fTRDOn&&abs(fPdg)==2212&&fTRDtrack.fN>40&&fStatus[2]>1") 
112 comp.fTree->Draw("fESDTrack.fTRDsignal:MC.fParticle.P()>>dEdx2","fTRDOn&&abs(fPdg)==321&&fTRDtrack.fN>40&&fStatus[2]>1") 
113 comp.fTree->Draw("fESDTrack.fTRDsignal:MC.fParticle.P()>>dEdx3","fTRDOn&&abs(fPdg)==11&&fTRDtrack.fN>40&&fStatus[2]>1") 
114
115 comp.fTree->Draw("fESDTrack.fTPCsignal:fTPCinP0[4]>>dEdx0","fTPCOn&&abs(fPdg)==211&&fESDTrack.fTPCncls>180&&fESDTrack.fTPCsignal>10"+cteta1); 
116 comp.fTree->Draw("fESDTrack.fTPCsignal:fTPCinP0[4]>>dEdx1","fTPCOn&&abs(fPdg)==2212&&fESDTrack.fTPCncls>180&&fESDTrack.fTPCsignal>10"+cteta1); 
117 comp.fTree->Draw("fESDTrack.fTPCsignal:fTPCinP0[4]>>dEdx2","fTPCOn&&abs(fPdg)==321&&fESDTrack.fTPCncls>180&&fESDTrack.fTPCsignal>10"+cteta1); 
118 comp.fTree->Draw("fESDTrack.fTPCsignal:fTPCinP0[4]>>dEdx3","fTPCOn&&abs(fPdg)==11&&fESDTrack.fTPCncls>180&&fESDTrack.fTPCsignal>10"+cteta1); 
119
120 hdedx3->SetXTitle("P(GeV/c)");
121 hdedx3->SetYTitle("dEdx(unit)");
122 hdedx3->Draw(); hdedx1->Draw("same"); hdedx2->Draw("same"); hdedx0->Draw("same");
123
124 comp.DrawXY("fITSinP0[3]","fITSPools[4]","fReconstructed==1&&fPdg==-211&&fITSOn"+cprim,"1",4,0.2,1.0,-8,8)
125
126 TProfile prof("prof","prof",10,0.5,5);
127
128
129
130
131 */
132
133
134 #if !defined(__CINT__) || defined(__MAKECINT__)
135 #include <stdio.h>
136 #include <string.h>
137 //ROOT includes
138 #include "Rtypes.h"
139 #include "TFile.h"
140 #include "TTree.h"
141 #include "TChain.h"
142 #include "TCut.h"
143 #include "TString.h"
144 #include "TBenchmark.h"
145 #include "TStopwatch.h"
146 #include "TParticle.h"
147 #include "TSystem.h"
148 #include "TTimer.h"
149 #include "TVector3.h"
150 #include "TPad.h"
151 #include "TCanvas.h"
152 #include "TH1F.h"
153 #include "TH2F.h"
154 #include "TF1.h"
155 #include "TText.h"
156 #include "Getline.h"
157 #include "TStyle.h"
158
159 //ALIROOT includes
160 #include "AliRun.h"
161 #include "AliStack.h"
162 #include "AliESDtrack.h"
163 #include "AliSimDigits.h"
164 #include "AliTPCParam.h"
165 #include "AliTPC.h"
166 #include "AliTPCLoader.h"
167 #include "AliDetector.h"
168 #include "AliTrackReference.h"
169 #include "AliRun.h"
170 #include "AliTPCParamSR.h"
171 #include "AliTracker.h"
172 #include "AliComplexCluster.h"
173 #include "AliMagF.h"
174 #include "AliESD.h"
175 #include "AliESDtrack.h"
176 #include "AliITStrackMI.h"
177 #include "AliTRDtrack.h"
178 #include "AliHelix.h"
179 #include "AliESDVertex.h"
180 #include "AliExternalTrackParam.h"
181 #include "AliESDkink.h"
182 #include "AliESDV0MI.h"
183
184 #endif
185 #include "AliGenInfo.h"
186 #include "AliESDComparisonMI.h"
187
188
189
190
191
192 void MakeAliases(AliESDComparisonDraw&comp)
193 {
194   //
195   // aliases definition
196   //
197   comp.fTree->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)");
198   comp.fTree->SetAlias("direction","MC.fParticle.fVx*MC.fParticle.fPx+MC.fParticle.fVy*MC.fParticle.fPy");
199   comp.fTree->SetAlias("decaydir","MC.fTRdecay.fX*MC.fTRdecay.fPx+MC.fTRdecay.fY*MC.fTRdecay.fPy");
200   comp.fTree->SetAlias("theta","MC.fTrackRef.Theta()");
201   comp.fTree->SetAlias("primdca","sqrt(RC.fITStrack.fD[0]**2+RC.fITStrack.fD[1]**2)");
202   comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
203   comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
204   
205   comp.fTree->SetAlias("trddedx","(RC.fESDTrack.fTRDsignals[0]+RC.fESDTrack.fTRDsignals[1]+RC.fESDTrack.fTRDsignals[2]+RC.fESDTrack.fTRDsignals[3]+RC.fESDTrack.fTRDsignals[4]+RC.fESDTrack.fTRDsignals[5])/6.");
206   
207   comp.fTree->SetAlias("dtofmc2","fESDTrack.fTrackTime[2]-(10^12*MC.fTOFReferences[0].fTime)");
208   comp.fTree->SetAlias("dtofrc2","(fESDTrack.fTrackTime[2]-fESDTrack.fTOFsignal)");
209
210   comp.fTree->SetAlias("psum","fESDTrack.fTOFr[4]+fESDTrack.fTOFr[3]+fESDTrack.fTOFr[2]+fESDTrack.fTOFr[1]+fESDTrack.fTOFr[0]");
211   comp.fTree->SetAlias("P0","fESDTrack.fTOFr[0]/psum");
212   comp.fTree->SetAlias("P1","fESDTrack.fTOFr[1]/psum");
213   comp.fTree->SetAlias("P2","fESDTrack.fTOFr[2]/psum");
214   comp.fTree->SetAlias("P3","fESDTrack.fTOFr[3]/psum");
215   comp.fTree->SetAlias("P4","fESDTrack.fTOFr[4]/psum");
216   comp.fTree->SetAlias("MaxP","max(max(max(P0,P1),max(P2,P3)),P4)");
217 }
218
219
220
221
222 void  AliESDRecInfo::UpdatePoints(AliESDtrack*track)
223 {
224   //
225   //
226   Int_t iclusters[200];
227   Float_t density[160];
228   for (Int_t i=0;i<160;i++) density[i]=-1.;
229   fTPCPoints[0]= 160;
230   fTPCPoints[1] = -1;
231   //
232   if (fTPCPoints[0]<fTPCPoints[1]) return;
233   //  Int_t nclusters=track->GetTPCclusters(iclusters);
234
235   Int_t ngood=0;
236   Int_t undeff=0;
237   Int_t nall =0;
238   Int_t range=20;
239   for (Int_t i=0;i<160;i++){
240     Int_t last = i-range;
241     if (nall<range) nall++;
242     if (last>=0){
243       if (iclusters[last]>0&& (iclusters[last]&0x8000)==0) ngood--;
244       if (iclusters[last]==-1) undeff--;
245     }
246     if (iclusters[i]>0&& (iclusters[i]&0x8000)==0)   ngood++;
247     if (iclusters[i]==-1) undeff++;
248     if (nall==range &&undeff<range/2) density[i-range/2] = Float_t(ngood)/Float_t(nall-undeff);
249   }
250   Float_t maxdens=0;
251   Int_t indexmax =0;
252   for (Int_t i=0;i<160;i++){
253     if (density[i]<0) continue;
254     if (density[i]>maxdens){
255       maxdens=density[i];
256       indexmax=i;
257     }
258   }
259   //
260   //max dens point
261   fTPCPoints[3] = maxdens;
262   fTPCPoints[1] = indexmax;
263   //
264   // last point
265   for (Int_t i=indexmax;i<160;i++){
266     if (density[i]<0) continue;
267     if (density[i]<maxdens/2.) {
268       break;
269     }
270     fTPCPoints[2]=i;
271   }
272   //
273   // first point
274   for (Int_t i=indexmax;i>0;i--){
275     if (density[i]<0) continue;
276     if (density[i]<maxdens/2.) {
277       break;
278     }
279     fTPCPoints[0]=i;
280   }
281   //
282   // Density at the last 30 padrows
283   //
284   // 
285   nall  = 0;
286   ngood = 0;
287   for (Int_t i=159;i>0;i--){
288     if (iclusters[i]==-1) continue; //dead zone
289     nall++;
290     if (iclusters[i]>0)   ngood++;
291     if (nall>20) break;
292   }
293   fTPCPoints[4] = Float_t(ngood)/Float_t(nall);
294   //
295   if ((track->GetStatus()&AliESDtrack::kITSrefit)>0) fTPCPoints[0]=-1;
296   //
297   //
298   // check TRDPoints
299   /*
300   nclusters=track->GetTRDclusters(iclusters);
301   for (Int_t i=nclusters;i>0;i--){
302     
303   }
304   */
305
306
307 }
308
309 //
310 //
311 void AliESDRecInfo::Update(AliMCInfo* info,AliTPCParam * /*par*/, Bool_t reconstructed)
312 {
313   //
314   //
315   //calculates derived variables
316   //  
317   //
318   UpdatePoints(&fESDTrack);
319   fBestTOFmatch=1000;
320   AliTrackReference * ref = &(info->fTrackRef);
321   fTPCinR0[0] = info->fTrackRef.X();    
322   fTPCinR0[1] = info->fTrackRef.Y();    
323   fTPCinR0[2] = info->fTrackRef.Z();
324   fTPCinR0[3] = TMath::Sqrt(fTPCinR0[0]*fTPCinR0[0]+fTPCinR0[1]*fTPCinR0[1]);
325   fTPCinR0[4] = TMath::ATan2(fTPCinR0[1],fTPCinR0[0]);
326   //
327   fTPCinP0[0] = ref->Px();
328   fTPCinP0[1] = ref->Py();
329   fTPCinP0[2] = ref->Pz();
330   fTPCinP0[3] = ref->Pt();
331   fTPCinP0[4] = ref->P();
332   fDeltaP     = (ref->P()-info->fParticle.P())/info->fParticle.P();
333   //
334   //
335   if (fTPCinP0[3]>0.0000001){
336     //
337     fTPCAngle0[0] = TMath::ATan2(fTPCinP0[1],fTPCinP0[0]);
338     fTPCAngle0[1] = TMath::ATan(fTPCinP0[2]/fTPCinP0[3]);
339   }
340   //
341   //
342   fITSinP0[0]=info->fParticle.Px();
343   fITSinP0[1]=info->fParticle.Py();
344   fITSinP0[2]=info->fParticle.Pz();
345   fITSinP0[3]=info->fParticle.Pt();    
346   //
347   fITSinR0[0]=info->fParticle.Vx();
348   fITSinR0[1]=info->fParticle.Vy();
349   fITSinR0[2]=info->fParticle.Vz();
350   fITSinR0[3] = TMath::Sqrt(fITSinR0[0]*fITSinR0[0]+fITSinR0[1]*fITSinR0[1]);
351   fITSinR0[4] = TMath::ATan2(fITSinR0[1],fITSinR0[0]);
352   //
353   //
354   if (fITSinP0[3]>0.0000001){
355     fITSAngle0[0] = TMath::ATan2(fITSinP0[1],fITSinP0[0]);
356     fITSAngle0[1] = TMath::ATan(fITSinP0[2]/fITSinP0[3]);
357   }
358   //
359   for (Int_t i=0;i<4;i++) fStatus[i] =0;
360   fReconstructed = kFALSE;
361   fTPCOn = kFALSE;
362   fITSOn = kFALSE;
363   fTRDOn = kFALSE;  
364   if (reconstructed==kFALSE) return;
365
366   fLabels[0] = info->fLabel;
367   fLabels[1] = info->fPrimPart;
368   fReconstructed = kTRUE;
369   fTPCOn = ((fESDTrack.GetStatus()&AliESDtrack::kTPCrefit)>0) ? kTRUE : kFALSE;
370   fITSOn = ((fESDTrack.GetStatus()&AliESDtrack::kITSrefit)>0) ? kTRUE : kFALSE;
371   fTRDOn = ((fESDTrack.GetStatus()&AliESDtrack::kTRDrefit)>0) ? kTRUE : kFALSE;
372   //
373   //  
374   if ((fESDTrack.GetStatus()&AliESDtrack::kTPCrefit)>0){
375     fStatus[1] =3;
376   }
377   else{
378     if ((fESDTrack.GetStatus()&AliESDtrack::kTPCout)>0){
379       fStatus[1] =2;
380     }
381     else{
382       fStatus[1]=1;
383     }      
384   }
385   //
386   if ((fESDTrack.GetStatus()&AliESDtrack::kITSout)>0){
387     fStatus[0] =2;
388   }
389   else{
390     if ((fESDTrack.GetStatus()&AliESDtrack::kITSrefit)>0){
391       fStatus[0] =1;
392     }
393     else{
394       fStatus[0]=0;
395     }      
396   }
397
398   //
399   //
400   if ((fESDTrack.GetStatus()&AliESDtrack::kTRDrefit)>0){
401     fStatus[2] =2;
402   }
403   else{
404     if ((fESDTrack.GetStatus()&AliESDtrack::kTRDout)>0){
405       fStatus[2] =1;
406     }
407   }
408   if ((fESDTrack.GetStatus()&AliESDtrack::kTRDStop)>0){
409     fStatus[2] =10;
410   }
411
412   //
413   //TOF 
414   // 
415   if (((fESDTrack.GetStatus()&AliESDtrack::kTOFout)>0)){
416     //
417     // best tof match
418     Double_t times[5];
419     fESDTrack.GetIntegratedTimes(times);    
420     for (Int_t i=0;i<5;i++){
421       if ( TMath::Abs(fESDTrack.GetTOFsignal()-times[i]) <TMath::Abs(fBestTOFmatch) ){
422         fBestTOFmatch = fESDTrack.GetTOFsignal()-times[i];
423       }
424     }
425     Int_t toflabel[3];
426     fESDTrack.GetTOFLabel(toflabel);
427     Bool_t toffake=kTRUE;
428     Bool_t tofdaughter=kFALSE;
429     for (Int_t i=0;i<3;i++){
430       if (toflabel[i]<0) continue;      
431       if (toflabel[i]== TMath::Abs(fESDTrack.GetLabel()))  toffake=kFALSE;      
432       if (toflabel[i]==info->fParticle.GetDaughter(0) || (toflabel[i]==info->fParticle.GetDaughter(1))) tofdaughter=kTRUE;  // decay product of original particle
433       fStatus[3]=1;
434     }
435     if (toffake) fStatus[3] =3;       //total fake
436     if (tofdaughter) fStatus[3]=2;    //fake because of decay
437   }else{
438     fStatus[3]=0;
439   }
440
441
442   if (fStatus[1]>0 &&info->fNTPCRef>0&&TMath::Abs(fTPCinP0[3])>0.0001){
443     //TPC
444     fESDTrack.GetInnerXYZ(fTPCinR1);
445     fTPCinR1[3] = TMath::Sqrt(fTPCinR1[0]*fTPCinR1[0]+fTPCinR1[1]*fTPCinR1[1]);
446     fTPCinR1[4] = TMath::ATan2(fTPCinR1[1],fTPCinR1[0]);        
447     fESDTrack.GetInnerPxPyPz(fTPCinP1);
448     fTPCinP1[3] = TMath::Sqrt(fTPCinP1[0]*fTPCinP1[0]+fTPCinP1[1]*fTPCinP1[1]);
449     fTPCinP1[4] = TMath::Sqrt(fTPCinP1[3]*fTPCinP1[3]+fTPCinP1[2]*fTPCinP1[2]);
450     //
451     //
452     if (fTPCinP1[3]>0.0000001){
453       fTPCAngle1[0] = TMath::ATan2(fTPCinP1[1],fTPCinP1[0]);
454       fTPCAngle1[1] = TMath::ATan(fTPCinP1[2]/fTPCinP1[3]);  
455     }    
456     Double_t cov[15], param[5],x;
457     fESDTrack.GetInnerExternalCovariance(cov);
458     fESDTrack.GetInnerExternalParameters(x,param);
459     if (x<50) return ;
460     //
461     fTPCDelta[0] = (fTPCinR0[4]-fTPCinR1[4])*fTPCinR1[3];  //delta rfi
462     fTPCPools[0] = fTPCDelta[0]/TMath::Sqrt(cov[0]);
463     fTPCDelta[1] = (fTPCinR0[2]-fTPCinR1[2]);              //delta z
464     fTPCPools[1] = fTPCDelta[1]/TMath::Sqrt(cov[2]);
465     fTPCDelta[2] = (fTPCAngle0[0]-fTPCAngle1[0]);
466     fTPCPools[2] = fTPCDelta[2]/TMath::Sqrt(cov[5]);
467     fTPCDelta[3] = (TMath::Tan(fTPCAngle0[1])-TMath::Tan(fTPCAngle1[1]));
468     fTPCPools[3] = fTPCDelta[3]/TMath::Sqrt(cov[9]);
469     fTPCDelta[4] = (fTPCinP0[3]-fTPCinP1[3]);
470     Double_t sign = (param[4]>0)? 1.:-1; 
471     fSign =sign;
472     fTPCPools[4] = sign*(1./fTPCinP0[3]-1./fTPCinP1[3])/TMath::Sqrt(TMath::Abs(cov[14]));    
473   }
474   if (fITSOn){
475     // ITS 
476     Double_t param[5],x;
477     fESDTrack.GetExternalParameters(x,param);   
478     //    fESDTrack.GetConstrainedExternalParameters(x,param);   
479     Double_t cov[15];
480     fESDTrack.GetExternalCovariance(cov);
481     //fESDTrack.GetConstrainedExternalCovariance(cov);
482     if (TMath::Abs(param[4])<0.0000000001) return;
483
484     fESDTrack.GetXYZ(fITSinR1);
485     fESDTrack.GetPxPyPz(fITSinP1);
486     fITSinP1[3] = TMath::Sqrt(fITSinP1[0]*fITSinP1[0]+fITSinP1[1]*fITSinP1[1]);
487     //
488     fITSinR1[3] = TMath::Sqrt(fITSinR1[0]*fITSinR1[0]+fITSinR1[1]*fITSinR1[1]);
489     fITSinR1[4] = TMath::ATan2(fITSinR1[1],fITSinR1[0]);
490     //
491     //
492     if (fITSinP1[3]>0.0000001){
493       fITSAngle1[0] = TMath::ATan2(fITSinP1[1],fITSinP1[0]);
494       fITSAngle1[1] = TMath::ATan(fITSinP1[2]/fITSinP1[3]);  
495     }
496     //
497     //
498     fITSDelta[0] = (fITSinR0[4]-fITSinR1[4])*fITSinR1[3];  //delta rfi
499     fITSPools[0] = fITSDelta[0]/TMath::Sqrt(cov[0]);
500     fITSDelta[1] = (fITSinR0[2]-fITSinR1[2]);              //delta z
501     fITSPools[1] = fITSDelta[1]/TMath::Sqrt(cov[2]);
502     fITSDelta[2] = (fITSAngle0[0]-fITSAngle1[0]);
503     fITSPools[2] = fITSDelta[2]/TMath::Sqrt(cov[5]);
504     fITSDelta[3] = (TMath::Tan(fITSAngle0[1])-TMath::Tan(fITSAngle1[1]));
505     fITSPools[3] = fITSDelta[3]/TMath::Sqrt(cov[9]);
506     fITSDelta[4] = (fITSinP0[3]-fITSinP1[3]);    
507     Double_t sign = (param[4]>0) ? 1:-1; 
508     fSign = sign;
509     fITSPools[4] = sign*(1./fITSinP0[3]-1./fITSinP1[3])/TMath::Sqrt(cov[14]);    
510   }
511   
512 }
513
514
515 void  AliESDRecV0Info::Update(Float_t vertex[3])
516
517
518   if ( (fT1.fStatus[1]>0)&& (fT2.fStatus[1]>0)){
519     Float_t distance1,distance2;
520     Double_t xx[3],pp[3];
521     //
522     Double_t xd[3],pd[3],signd;
523     Double_t xm[3],pm[3],signm;
524     //
525     //
526     if (fT1.fITSOn&&fT2.fITSOn){
527       for (Int_t i=0;i<3;i++){
528         xd[i] = fT2.fITSinR1[i];
529         pd[i] = fT2.fITSinP1[i];
530         xm[i] = fT1.fITSinR1[i];
531         pm[i] = fT1.fITSinP1[i];
532       }
533     }
534     else{
535       
536       for (Int_t i=0;i<3;i++){
537         xd[i] = fT2.fTPCinR1[i];
538         pd[i] = fT2.fTPCinP1[i];
539         xm[i] = fT1.fTPCinR1[i];
540         pm[i] = fT1.fTPCinP1[i];
541       }
542     }
543     //
544     //
545     signd =  fT2.fSign<0 ? -1:1;
546     signm =  fT1.fSign<0 ? -1:1;
547
548     AliHelix dhelix1(xd,pd,signd);
549     dhelix1.GetMomentum(0,pp,0);
550     dhelix1.Evaluate(0,xx);      
551     // 
552     //  Double_t x2[3],p2[3];
553     //            
554     AliHelix mhelix(xm,pm,signm);    
555     //
556     //find intersection linear
557     //
558     Double_t phase[2][2],radius[2];
559     Int_t  points = dhelix1.GetRPHIintersections(mhelix, phase, radius,200);
560     Double_t delta1=10000,delta2=10000;  
561
562     if (points==1){
563       fRs[0] = TMath::Sqrt(radius[0]);
564       fRs[1] = TMath::Sqrt(radius[0]);
565     }
566     if (points==2){
567       fRs[0] =TMath::Min(TMath::Sqrt(radius[0]),TMath::Sqrt(radius[1]));
568       fRs[1] =TMath::Max(TMath::Sqrt(radius[0]),TMath::Sqrt(radius[1]));
569     }
570     
571     if (points>0){
572       dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
573       dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
574       dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
575     }
576     if (points==2){    
577       dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
578       dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
579       dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
580     }
581     if (points==1){
582       fRs[0] = TMath::Sqrt(radius[0]);
583       fRs[1] = TMath::Sqrt(radius[0]);
584       fDistMinR = delta1;
585     }
586     if (points==2){
587       if (radius[0]<radius[1]){
588         fRs[0] = TMath::Sqrt(radius[0]);
589         fRs[1] = TMath::Sqrt(radius[1]);
590         fDistMinR = delta1;
591       }
592       else{
593         fRs[0] = TMath::Sqrt(radius[1]);
594         fRs[1] = TMath::Sqrt(radius[0]);
595         fDistMinR = delta2;
596       }
597     }
598     //
599     //
600     distance1 = TMath::Min(delta1,delta2);
601     //
602     //find intersection parabolic
603     //
604     points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
605     delta1=10000,delta2=10000;  
606     
607     if (points>0){
608       dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
609     }
610     if (points==2){    
611       dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
612     }
613     
614     distance2 = TMath::Min(delta1,delta2);
615     if (distance2>100) fDist2 =100;
616     return;
617     if (delta1<delta2){
618       //get V0 info
619       dhelix1.Evaluate(phase[0][0],fXr);
620       dhelix1.GetMomentum(phase[0][0],fPdr);
621       mhelix.GetMomentum(phase[0][1],fPm);
622       dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fAngle);
623       fRr = TMath::Sqrt(radius[0]);
624     }
625     else{
626       dhelix1.Evaluate(phase[1][0],fXr);
627       dhelix1.GetMomentum(phase[1][0], fPdr);
628       mhelix.GetMomentum(phase[1][1], fPm);
629       dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fAngle);
630       fRr = TMath::Sqrt(radius[1]);
631     }
632     fDist1 = TMath::Sqrt(distance1);
633     fDist2 = TMath::Sqrt(distance2);      
634     
635     if (fDist2<10.5){
636       Double_t x,alpha,param[5],cov[15];
637       //
638       fT1.fESDTrack.GetInnerExternalParameters(alpha,x,param);
639       fT1.fESDTrack.GetInnerExternalCovariance(cov);
640       AliExternalTrackParam paramm(x,alpha,param,cov);
641       //
642       fT2.fESDTrack.GetInnerExternalParameters(alpha,x,param);
643       fT2.fESDTrack.GetInnerExternalCovariance(cov);
644       AliExternalTrackParam paramd(x,alpha,param,cov);
645     }    
646     //            
647     //   
648     
649     Float_t v[3] = {fXr[0]-vertex[0],fXr[1]-vertex[1],fXr[2]-vertex[2]};
650     Float_t p[3] = {fPdr[0]+fPm[0], fPdr[1]+fPm[1],fPdr[2]+fPm[2]};
651     
652     Float_t vnorm2 = v[0]*v[0]+v[1]*v[1];
653     Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
654     vnorm2 = TMath::Sqrt(vnorm2);
655     Float_t pnorm2 = p[0]*p[0]+p[1]*p[1];
656     Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
657     pnorm2 = TMath::Sqrt(pnorm2);
658     
659     fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
660     fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);  
661     fPointAngle   = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
662   }
663 }
664
665 ////
666 void  AliESDRecKinkInfo::Update()
667 {
668
669   if ( (fT1.fTPCOn)&& (fT2.fTPCOn)){
670     //
671     // IF BOTH RECONSTRUCTED
672     Float_t distance1,distance2;
673     Double_t xx[3],pp[3];
674     //
675     Double_t xd[3],pd[3],signd;
676     Double_t xm[3],pm[3],signm;
677     for (Int_t i=0;i<3;i++){
678       xd[i] = fT2.fTPCinR1[i];
679       pd[i] = fT2.fTPCinP1[i];
680       xm[i] = fT1.fTPCinR1[i];
681       pm[i] = fT1.fTPCinP1[i];
682     }
683     signd =  fT2.fSign<0 ? -1:1;
684     signm =  fT1.fSign<0 ? -1:1;
685
686     AliHelix dhelix1(xd,pd,signd);
687     dhelix1.GetMomentum(0,pp,0);
688     dhelix1.Evaluate(0,xx);      
689     // 
690     //  Double_t x2[3],p2[3];
691     //            
692     AliHelix mhelix(xm,pm,signm);    
693     //
694     //find intersection linear
695     //
696     Double_t phase[2][2],radius[2];
697     Int_t  points = dhelix1.GetRPHIintersections(mhelix, phase, radius,200);
698     Double_t delta1=10000,delta2=10000;  
699
700     if (points==1){
701       fMinR = TMath::Sqrt(radius[0]);
702     }
703     if (points==2){
704       fMinR =TMath::Min(TMath::Sqrt(radius[0]),TMath::Sqrt(radius[1]));
705     }
706     
707     if (points>0){
708       dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
709       dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
710       dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
711     }
712     if (points==2){    
713       dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
714       dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
715       dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
716     }
717     if (points==1){
718       fMinR = TMath::Sqrt(radius[0]);
719       fDistMinR = delta1;
720     }
721     if (points==2){
722       if (radius[0]<radius[1]){
723         fMinR = TMath::Sqrt(radius[0]);
724         fDistMinR = delta1;
725       }
726       else{
727         fMinR = TMath::Sqrt(radius[1]);
728         fDistMinR = delta2;
729       }
730     }
731     //
732     //
733     distance1 = TMath::Min(delta1,delta2);
734     //
735     //find intersection parabolic
736     //
737     points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
738     delta1=10000,delta2=10000;  
739     
740     if (points>0){
741       dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
742     }
743     if (points==2){    
744       dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
745     }
746     
747     distance2 = TMath::Min(delta1,delta2);
748     if (delta1<delta2){
749       //get V0 info
750       dhelix1.Evaluate(phase[0][0],fXr);
751       dhelix1.GetMomentum(phase[0][0],fPdr);
752       mhelix.GetMomentum(phase[0][1],fPm);
753       dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fAngle);
754       fRr = TMath::Sqrt(radius[0]);
755     }
756     else{
757       dhelix1.Evaluate(phase[1][0],fXr);
758       dhelix1.GetMomentum(phase[1][0], fPdr);
759       mhelix.GetMomentum(phase[1][1], fPm);
760       dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fAngle);
761       fRr = TMath::Sqrt(radius[1]);
762     }
763     fDist1 = TMath::Sqrt(distance1);
764     fDist2 = TMath::Sqrt(distance2);      
765     
766     if (fDist2<10.5){
767       Double_t x,alpha,param[5],cov[15];
768       //
769       fT1.fESDTrack.GetInnerExternalParameters(alpha,x,param);
770       fT1.fESDTrack.GetInnerExternalCovariance(cov);
771       AliExternalTrackParam paramm(x,alpha,param,cov);
772       //
773       fT2.fESDTrack.GetInnerExternalParameters(alpha,x,param);
774       fT2.fESDTrack.GetInnerExternalCovariance(cov);
775       AliExternalTrackParam paramd(x,alpha,param,cov);
776       /*
777       AliESDkink kink;
778       kink.Update(&paramm,&paramd);
779       //      kink.Dump();
780       Double_t diff  = kink.fRr-fRr;
781       Double_t diff2 = kink.fDist2-fDist2;
782       printf("Diff\t%f\t%f\n",diff,diff2);
783       */
784     }
785     
786     //            
787     //
788   }
789
790 }
791
792
793 ////////////////////////////////////////////////////////////////////////
794 ESDCmpTr::ESDCmpTr()
795 {
796   Reset();
797 }
798
799 ////////////////////////////////////////////////////////////////////////
800 ESDCmpTr::ESDCmpTr(const char* fnGenTracks,
801                    const char* fnCmp,
802                    const char* fnGalice, Int_t direction,
803                    Int_t nEvents, Int_t firstEvent)
804 {
805   Reset();
806   //  fFnGenTracks = fnGenTracks;
807   //  fFnCmp = fnCmp;
808   sprintf(fFnGenTracks,"%s",fnGenTracks);
809   sprintf(fFnCmp,"%s",fnCmp);
810
811   fFirstEventNr = firstEvent;
812   fEventNr = firstEvent;
813   fNEvents = nEvents;
814   fDirection = direction;
815   //
816   fLoader = AliRunLoader::Open(fnGalice);
817   if (gAlice){
818     //delete gAlice->GetRunLoader();
819     delete gAlice;
820     gAlice = 0x0;
821   }
822   if (fLoader->LoadgAlice()){
823     cerr<<"Error occured while l"<<endl;
824   }
825   Int_t nall = fLoader->GetNumberOfEvents();
826   if (nEvents==0) {
827     nEvents =nall;
828     fNEvents=nall;
829     fFirstEventNr=0;
830   }    
831
832   if (nall<=0){
833     cerr<<"no events available"<<endl;
834     fEventNr = 0;
835     return;
836   }
837   if (firstEvent+nEvents>nall) {
838     fEventNr = nall-firstEvent;
839     cerr<<"restricted number of events availaible"<<endl;
840   }
841   AliMagF * magf = gAlice->Field();
842   AliTracker::SetFieldMap(magf,0);
843
844 }
845
846
847 ////////////////////////////////////////////////////////////////////////
848 ESDCmpTr::~ESDCmpTr()
849 {
850   if (fLoader) {
851     delete fLoader;
852   }
853 }
854
855 //////////////////////////////////////////////////////////////
856 Int_t ESDCmpTr::SetIO()
857 {
858   //
859   // 
860   CreateTreeCmp();
861   if (!fTreeCmp) return 1;
862   fParamTPC = GetTPCParam();
863   //
864   if (!ConnectGenTree()) {
865     cerr<<"Cannot connect tree with generated tracks"<<endl;
866     return 1;
867   }
868   return 0;
869 }
870
871 //////////////////////////////////////////////////////////////
872
873 Int_t ESDCmpTr::SetIO(Int_t eventNr)
874 {
875   //
876   // 
877   // SET INPUT
878   //
879   TFile f("AliESDs.root");
880   //
881  
882   TTree* tree = (TTree*) f.Get("esdTree");
883   if (!tree) { 
884     Char_t ename[100]; 
885     sprintf(ename,"%d",eventNr);
886     fEvent = (AliESD*)f.Get(ename);
887     if (!fEvent){
888       sprintf(ename,"ESD%d",eventNr);
889       fEvent = (AliESD*)f.Get(ename);
890     }
891   }
892   else{
893     tree->SetBranchAddress("ESD", &fEvent);
894     tree->GetEntry(eventNr);
895   }
896
897
898   /*
899   Char_t ename[100]; 
900   sprintf(ename,"%d",eventNr);
901   fEvent = (AliESD*)f.Get(ename);
902   if (!fEvent){
903     sprintf(ename,"ESD%d",eventNr);
904     fEvent = (AliESD*)f.Get(ename);
905   }
906   
907   TTree* tree = (TTree*) f.Get("esdTree");
908   if (!tree) {
909     Error("CheckESD", "no ESD tree found");
910     return kFALSE;
911   }
912   tree->SetBranchAddress("ESD", &fEvent);
913   tree->GetEntry(eventNr);
914   */
915
916   if (!fEvent) return 1;
917
918   return 0;
919 }
920
921
922
923 ////////////////////////////////////////////////////////////////////////
924 void ESDCmpTr::Reset()
925 {
926   fEventNr = 0;
927   fNEvents = 0;
928   fTreeCmp = 0;
929   fTreeCmpKinks =0;
930   fTreeCmpV0 =0;
931   //  fFnCmp = "cmpTracks.root";
932   fFileGenTracks = 0;
933   fDebug = 0;
934   //
935   fParamTPC = 0;
936   fEvent =0;
937 }
938
939 ////////////////////////////////////////////////////////////////////////
940 Int_t ESDCmpTr::Exec(Int_t nEvents, Int_t firstEventNr)
941 {
942   fNEvents = nEvents;
943   fFirstEventNr = firstEventNr;
944   return Exec();
945 }
946
947 ////////////////////////////////////////////////////////////////////////
948 Int_t ESDCmpTr::Exec()
949 {
950   TStopwatch timer;
951   timer.Start();
952
953   if (SetIO()==1) 
954     return 1;
955    
956   fNextTreeGenEntryToRead = 0;
957   fNextKinkToRead = 0;
958   fNextV0ToRead   =0;
959   cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl;
960   for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents;
961        eventNr++) {
962     fEventNr = eventNr;
963     SetIO(fEventNr);
964     fNParticles = gAlice->GetEvent(fEventNr);    
965
966     fIndexRecTracks = new Short_t[fNParticles*20];  //write at maximum 4 tracks corresponding to particle
967     fIndexRecKinks  = new Short_t[fNParticles*20];  //write at maximum 20 tracks corresponding to particle
968     fIndexRecV0  = new Short_t[fNParticles*20];  //write at maximum 20 tracks corresponding to particle
969
970     fFakeRecTracks = new Short_t[fNParticles];
971     fMultiRecTracks = new Short_t[fNParticles];
972     fMultiRecKinks = new Short_t[fNParticles];
973     fMultiRecV0 = new Short_t[fNParticles];
974
975     for (Int_t i = 0; i<fNParticles; i++) {
976       for (Int_t j=0;j<20;j++){
977         fIndexRecTracks[20*i+j] = -1;
978         fIndexRecKinks[20*i+j]  = -1;
979         fIndexRecV0[20*i+j]  = -1;
980       }
981       fFakeRecTracks[i] = 0;
982       fMultiRecTracks[i] = 0;
983       fMultiRecKinks[i] = 0;
984       fMultiRecV0[i] = 0;      
985     }
986   
987     cout<<"Start to process event "<<fEventNr<<endl;
988     cout<<"\tfNParticles = "<<fNParticles<<endl;
989     if (fDebug>2) cout<<"\tStart loop over TreeT"<<endl;
990     if (TreeTLoop()>0) return 1;
991
992     if (fDebug>2) cout<<"\tStart loop over tree genTracks"<<endl;
993     if (TreeGenLoop(eventNr)>0) return 1;
994     BuildKinkInfo0(eventNr);
995     BuildV0Info(eventNr);
996     fRecArray->Delete();
997
998     if (fDebug>2) cout<<"\tEnd loop over tree genTracks"<<endl;
999
1000     delete [] fIndexRecTracks;
1001     delete [] fIndexRecKinks;
1002     delete [] fIndexRecV0;
1003     delete [] fFakeRecTracks;
1004     delete [] fMultiRecTracks;
1005     delete [] fMultiRecKinks;
1006     delete [] fMultiRecV0;
1007   }
1008
1009   CloseOutputFile();
1010
1011   cerr<<"Exec finished"<<endl;
1012   timer.Stop();
1013   timer.Print();
1014   return 0;
1015
1016 }
1017 ////////////////////////////////////////////////////////////////////////
1018 Bool_t ESDCmpTr::ConnectGenTree()
1019 {
1020 //
1021 // connect all branches from the genTracksTree
1022 // use the same variables as for the new cmp tree, it may work
1023 //
1024   fFileGenTracks = TFile::Open(fFnGenTracks,"READ");
1025   if (!fFileGenTracks) {
1026     cerr<<"Error in ConnectGenTree: cannot open file "<<fFnGenTracks<<endl;
1027     return kFALSE;
1028   }
1029   fTreeGenTracks = (TTree*)fFileGenTracks->Get("genTracksTree");
1030   if (!fTreeGenTracks) {
1031     cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
1032         <<fFnGenTracks<<endl;
1033     return kFALSE;
1034   }
1035   //
1036   fMCInfo = new AliMCInfo;
1037   fTreeGenTracks->SetBranchAddress("MC",&fMCInfo);
1038   //
1039   //
1040   fTreeGenKinks = (TTree*)fFileGenTracks->Get("genKinksTree");
1041   if (!fTreeGenKinks) {
1042     cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
1043         <<fFnGenTracks<<endl;
1044     //return kFALSE;
1045   }
1046   else{
1047     fGenKinkInfo = new AliGenKinkInfo;
1048     fTreeGenKinks->SetBranchAddress("MC",&fGenKinkInfo);
1049   }
1050
1051   fTreeGenV0 = (TTree*)fFileGenTracks->Get("genV0Tree");
1052   if (!fTreeGenV0) {
1053     cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
1054         <<fFnGenTracks<<endl;
1055     //return kFALSE;
1056   }
1057   else{
1058     fGenV0Info = new AliGenV0Info;
1059     fTreeGenV0->SetBranchAddress("MC",&fGenV0Info);
1060   }
1061   //
1062   if (fDebug > 1) {
1063     cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl;
1064   }
1065   return kTRUE;
1066 }
1067
1068
1069 ////////////////////////////////////////////////////////////////////////
1070 void ESDCmpTr::CreateTreeCmp() 
1071 {
1072   fFileCmp = TFile::Open(fFnCmp,"RECREATE");
1073   if (!fFileCmp) {
1074     cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl;
1075     return;
1076   }
1077   //
1078   //
1079   fTreeCmp    = new TTree("ESDcmpTracks","ESDcmpTracks");
1080   fMCInfo = new AliMCInfo;
1081   fRecInfo = new AliESDRecInfo;
1082   AliESDtrack * esdTrack = new AliESDtrack;
1083   //  AliITStrackMI * itsTrack = new AliITStrackMI;  
1084   fTreeCmp->Branch("MC","AliMCInfo",&fMCInfo,256000);
1085   fTreeCmp->Branch("RC","AliESDRecInfo",&fRecInfo,256000);
1086   //  fTreeCmp->Branch("fESDTrack","AliESDtrack",&esdTrack);
1087   //  fTreeCmp->Branch("ITS","AliITStrackMI",&itsTrack);
1088   delete esdTrack;
1089   //
1090   //
1091   fTreeCmpKinks    = new TTree("ESDcmpKinks","ESDcmpKinks"); 
1092   fGenKinkInfo     = new AliGenKinkInfo;
1093   fRecKinkInfo     = new AliESDRecKinkInfo;
1094   fTreeCmpKinks->Branch("MC.","AliGenKinkInfo",&fGenKinkInfo,256000);
1095   fTreeCmpKinks->Branch("RC.","AliESDRecKinkInfo",&fRecKinkInfo,256000);
1096   //
1097   //
1098   fTreeCmpV0       = new TTree("ESDcmpV0","ESDcmpV0"); 
1099   fGenV0Info     = new AliGenV0Info;
1100   fRecV0Info     = new AliESDRecV0Info;
1101   fTreeCmpV0->Branch("MC.","AliGenV0Info",   &fGenV0Info,256000);
1102   fTreeCmpV0->Branch("RC.","AliESDRecV0Info",&fRecV0Info,256000);
1103   //
1104   fTreeCmp->AutoSave(); 
1105   fTreeCmpKinks->AutoSave(); 
1106   fTreeCmpV0->AutoSave(); 
1107 }
1108 ////////////////////////////////////////////////////////////////////////
1109 void ESDCmpTr::CloseOutputFile()  
1110 {
1111   if (!fFileCmp) {
1112     cerr<<"File "<<fFnCmp<<" not found as an open file."<<endl;
1113     return;
1114   }
1115   fFileCmp->cd();
1116   fTreeCmp->Write();    
1117   delete fTreeCmp;
1118   
1119   fFileCmp->Close();
1120   delete fFileCmp;
1121   return;
1122 }
1123 ////////////////////////////////////////////////////////////////////////
1124
1125 TVector3 ESDCmpTr::TR2Local(AliTrackReference *trackRef,
1126                             AliTPCParam *paramTPC) {
1127
1128   Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
1129   Int_t index[4];
1130   paramTPC->Transform0to1(x,index);
1131   paramTPC->Transform1to2(x,index);
1132   return TVector3(x);
1133 }
1134 ////////////////////////////////////////////////////////////////////////
1135
1136 Int_t ESDCmpTr::TreeTLoop()
1137 {
1138   //
1139   // loop over all ESD reconstructed tracks and store info in memory
1140   //
1141   // + loop over all reconstructed kinks
1142   TStopwatch  timer;
1143   timer.Start();
1144   //  
1145   Int_t nEntries = (Int_t)fEvent->GetNumberOfTracks();  
1146   Int_t nKinks = (Int_t) fEvent->GetNumberOfKinks();
1147   Int_t nV0MIs = (Int_t) fEvent->GetNumberOfV0MIs();
1148   fSignedKinks = new Short_t[nKinks];
1149   fSignedV0    = new Short_t[nV0MIs];
1150   //
1151   // load kinks to the memory
1152   for (Int_t i=0; i<nKinks;i++){
1153     AliESDkink * kink =fEvent->GetKink(i);
1154     fSignedKinks[i]=0;
1155     if (kink->fStatus<0) continue;
1156   }
1157   //
1158   for (Int_t i=0; i<nV0MIs;i++){
1159     AliESDV0MI * v0MI =fEvent->GetV0MI(i);
1160     fSignedV0[i]=0;
1161     if (v0MI->fStatus<0) continue;
1162   }
1163   
1164   //
1165   //
1166   AliESDtrack * track=0;
1167   for (Int_t iEntry=0; iEntry<nEntries;iEntry++){
1168     //track = (AliESDtrack*)fTracks->UncheckedAt(iEntry);
1169     track = (AliESDtrack*)fEvent->GetTrack(iEntry);
1170     //
1171     Int_t label = track->GetLabel();
1172     Int_t absLabel = abs(label);
1173     if (absLabel < fNParticles) {
1174       //      fIndexRecTracks[absLabel] =  iEntry;
1175       if (label < 0) fFakeRecTracks[absLabel]++;      
1176       if (fMultiRecTracks[absLabel]>0){
1177         if (fMultiRecTracks[absLabel]<20)
1178           fIndexRecTracks[absLabel*20+fMultiRecTracks[absLabel]] =  iEntry;     
1179       }
1180       else      
1181         fIndexRecTracks[absLabel*20] =  iEntry;
1182       fMultiRecTracks[absLabel]++;
1183     }
1184   }
1185   // sort reconstructed kinks  
1186   //
1187   AliESDkink * kink=0;
1188   for (Int_t iEntry=0; iEntry<nKinks;iEntry++){
1189     kink = (AliESDkink*)fEvent->GetKink(iEntry);
1190     if (!kink) continue;
1191     //
1192     Int_t label0 = TMath::Abs(kink->fLab[0]);
1193     Int_t label1 = TMath::Abs(kink->fLab[1]);
1194     Int_t absLabel = TMath::Min(label0,label1);
1195     if (absLabel < fNParticles) {
1196       if (fMultiRecKinks[absLabel]>0){
1197         if (fMultiRecKinks[absLabel]<20)
1198           fIndexRecKinks[absLabel*20+fMultiRecKinks[absLabel]] =  iEntry;       
1199       }
1200       else      
1201         fIndexRecKinks[absLabel*20] =  iEntry;
1202       fMultiRecKinks[absLabel]++;
1203     }
1204   }  
1205   // --sort reconstructed V0
1206   //
1207   AliESDV0MI * v0MI=0;
1208   for (Int_t iEntry=0; iEntry<nV0MIs;iEntry++){
1209     v0MI = fEvent->GetV0MI(iEntry);
1210     if (!v0MI) continue;
1211     //
1212     //    Int_t label0 = TMath::Abs(v0MI->fLab[0]);
1213     //Int_t label1 = TMath::Abs(v0MI->fLab[1]);
1214     //
1215     for (Int_t i=0;i<2;i++){
1216       // Int_t absLabel = TMath::Min(label0,label1);
1217       Int_t absLabel =  TMath::Abs(v0MI->fLab[i]);
1218       if (absLabel < fNParticles) {
1219         if (fMultiRecV0[absLabel]>0){
1220           if (fMultiRecV0[absLabel]<20)
1221             fIndexRecV0[absLabel*20+fMultiRecV0[absLabel]] =  iEntry;   
1222         }
1223         else      
1224           fIndexRecV0[absLabel*20] =  iEntry;
1225         fMultiRecV0[absLabel]++;
1226       }
1227     }
1228   }  
1229
1230
1231   printf("Time spended in TreeTLoop\n");
1232   timer.Print();
1233   
1234   if (fDebug > 2) cerr<<"end of TreeTLoop"<<endl;  
1235   return 0;
1236 }
1237
1238 ////////////////////////////////////////////////////////////////////////
1239 Int_t ESDCmpTr::TreeGenLoop(Int_t eventNr)
1240 {
1241 //
1242 // loop over all entries for a given event, find corresponding 
1243 // rec. track and store in the fTreeCmp
1244 //
1245   TStopwatch timer;
1246   timer.Start();
1247   Int_t entry = fNextTreeGenEntryToRead;
1248   Double_t nParticlesTR = fTreeGenTracks->GetEntriesFast();
1249   cerr<<"fNParticles, nParticlesTR, fNextTreeGenEntryToRead: "<<fNParticles<<" "
1250       <<nParticlesTR<<" "<<fNextTreeGenEntryToRead<<endl;
1251   TBranch * branch = fTreeCmp->GetBranch("RC");
1252   branch->SetAddress(&fRecInfo); // set all pointers
1253   fRecArray = new TObjArray(fNParticles);
1254
1255   while (entry < nParticlesTR) {
1256     fTreeGenTracks->GetEntry(entry);
1257     entry++;
1258     if (eventNr < fMCInfo->fEventNr) continue;
1259     if (eventNr > fMCInfo->fEventNr) continue;;
1260     //
1261     fNextTreeGenEntryToRead = entry-1;
1262     if (fDebug > 2 && fMCInfo->fLabel < 10) {
1263       cerr<<"Fill track with a label "<<fMCInfo->fLabel<<endl;
1264     }
1265     //    if (fMCInfo->fNTPCRef<1) continue; // not TPCref
1266     //
1267     fRecInfo->Reset();
1268     AliESDtrack * track=0;
1269     fRecInfo->fReconstructed =0;
1270     TVector3 local = TR2Local(&(fMCInfo->fTrackRef),fParamTPC);
1271     local.GetXYZ(fRecInfo->fTRLocalCoord);      
1272     //
1273     if (fIndexRecTracks[fMCInfo->fLabel*20] >= 0) {
1274       //track= (AliESDtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4]);
1275       track= (AliESDtrack*)fEvent->GetTrack(fIndexRecTracks[fMCInfo->fLabel*20]);
1276       //
1277       //
1278       // find nearest track if multifound
1279       //Int_t sign = Int_t(track->GetSign()*fMCInfo->fCharge);
1280       //
1281       Int_t status = 0;
1282       if  ((track->GetStatus()&AliESDtrack::kITSrefit)>0) status++;
1283       if  ((track->GetStatus()&AliESDtrack::kTPCrefit)>0) status++;
1284       if  ((track->GetStatus()&AliESDtrack::kTRDrefit)>0) status++;
1285
1286       //
1287       if (fIndexRecTracks[fMCInfo->fLabel*20+1]>0){
1288         //
1289         Double_t p[3];
1290         track->GetInnerPxPyPz(p);
1291         Float_t maxp = p[0]*p[0]+p[1]*p[1]+p[2]*p[2];
1292         //
1293         for (Int_t i=1;i<20;i++){
1294           if (fIndexRecTracks[fMCInfo->fLabel*20+i]>=0){
1295             AliESDtrack * track2 = (AliESDtrack*)fEvent->GetTrack(fIndexRecTracks[fMCInfo->fLabel*20+i]);
1296             if (!track2) continue;
1297             //Int_t sign2 = track->GetSign()*fMCInfo->fCharge; //           
1298             //if (sign2<0) continue;
1299             track2->GetInnerPxPyPz(p);
1300             Float_t mom = p[0]*p[0]+p[1]*p[1]+p[2]*p[2];
1301             /*
1302             if (sign<0){
1303               sign = sign2;
1304               track = track2;
1305               maxp = mom;
1306               continue;
1307             }
1308             */
1309             //
1310             Int_t status2 = 0;
1311             if  ((track2->GetStatus()&AliESDtrack::kITSrefit)>0) status2++;
1312             if  ((track2->GetStatus()&AliESDtrack::kTPCrefit)>0) status2++;
1313             if  ((track2->GetStatus()&AliESDtrack::kTRDrefit)>0) status2++;
1314             if (status2<status) continue;
1315             //
1316             if (mom<maxp) continue;
1317             maxp = mom;
1318             track = track2;
1319             //
1320           }
1321         }
1322       } 
1323       //
1324       fRecInfo->fESDTrack =*track; 
1325       if (track->GetITStrack())
1326         fRecInfo->fITStrack = *((AliITStrackMI*)track->GetITStrack());
1327       else{
1328         fRecInfo->fITStrack = *track;
1329       }
1330       if (track->GetTRDtrack()){
1331         fRecInfo->fTRDtrack = *((AliTRDtrack*)track->GetTRDtrack());
1332       }
1333       else{
1334         fRecInfo->fTRDtrack.SetdEdx(-1);
1335       }
1336       fRecInfo->fReconstructed = 1;
1337       fRecInfo->fFake     = fFakeRecTracks[fMCInfo->fLabel];
1338       fRecInfo->fMultiple = fMultiRecTracks[fMCInfo->fLabel];
1339       //
1340       fRecInfo->Update(fMCInfo,fParamTPC,kTRUE);          
1341     }
1342     else{
1343       fRecInfo->Update(fMCInfo,fParamTPC,kFALSE);
1344     }
1345     fRecArray->AddAt(new AliESDRecInfo(*fRecInfo),fMCInfo->fLabel);
1346     fTreeCmp->Fill();
1347   }
1348   fTreeCmp->AutoSave();
1349   //fTracks->Delete();
1350   printf("Time spended in TreeGenLoop\n");
1351   timer.Print();
1352   if (fDebug > 2) cerr<<"end of TreeGenLoop"<<endl;
1353
1354   return 0;
1355 }
1356
1357
1358
1359 ////////////////////////////////////////////////////////////////////////
1360 ////////////////////////////////////////////////////////////////////////
1361 ////////////////////////////////////////////////////////////////////////
1362 Int_t ESDCmpTr::BuildKinkInfo0(Int_t eventNr)
1363 {
1364 //
1365 // loop over all entries for a given event, find corresponding 
1366 // rec. track and store in the fTreeCmp
1367 //
1368   TStopwatch timer;
1369   timer.Start();
1370   Int_t entry = fNextKinkToRead;
1371   Double_t nParticlesTR = fTreeGenKinks->GetEntriesFast();
1372   cerr<<"fNParticles, nParticlesTR, fNextKinkToRead: "<<fNParticles<<" "
1373       <<nParticlesTR<<" "<<fNextKinkToRead<<endl;
1374   //
1375   TBranch * branch = fTreeCmpKinks->GetBranch("RC.");
1376   branch->SetAddress(&fRecKinkInfo); // set all pointers
1377   
1378   //
1379   while (entry < nParticlesTR) {
1380     fTreeGenKinks->GetEntry(entry);
1381     entry++;
1382     if (eventNr < fGenKinkInfo->fMCm.fEventNr) continue;
1383     if (eventNr > fGenKinkInfo->fMCm.fEventNr) continue;;
1384     //
1385     fNextKinkToRead = entry-1;
1386     //
1387     //
1388     AliESDRecInfo*  fRecInfo1 = (AliESDRecInfo*)fRecArray->At(fGenKinkInfo->fMCm.fLabel);
1389     AliESDRecInfo*  fRecInfo2 = (AliESDRecInfo*)fRecArray->At(fGenKinkInfo->fMCd.fLabel);
1390     fRecKinkInfo->fT1 = (*fRecInfo1);
1391     fRecKinkInfo->fT2 = (*fRecInfo2);
1392     fRecKinkInfo->fStatus =0;
1393     if (fRecInfo1 && fRecInfo1->fTPCOn) fRecKinkInfo->fStatus+=1;
1394     if (fRecInfo2 && fRecInfo2->fTPCOn) fRecKinkInfo->fStatus+=2;
1395     if (fRecKinkInfo->fStatus==3&&fRecInfo1->fSign!=fRecInfo2->fSign) fRecKinkInfo->fStatus*=-1;
1396     
1397     if (fRecKinkInfo->fStatus==3){
1398       fRecKinkInfo->Update();    
1399     }
1400     Int_t label =  TMath::Min(fGenKinkInfo->fMCm.fLabel,fGenKinkInfo->fMCd.fLabel);
1401     Int_t label2 = TMath::Max(fGenKinkInfo->fMCm.fLabel,fGenKinkInfo->fMCd.fLabel);
1402     
1403     AliESDkink *kink=0;
1404     fRecKinkInfo->fRecStatus   =0;
1405     fRecKinkInfo->fMultiple    = fMultiRecKinks[label];
1406     fRecKinkInfo->fKinkMultiple=0;
1407     //
1408     if (fMultiRecKinks[label]>0){
1409
1410       //      for (Int_t j=0;j<TMath::Min(fMultiRecKinks[label],100);j++){
1411       for (Int_t j=TMath::Min(fMultiRecKinks[label],Short_t(20))-1;j>=0;j--){
1412         Int_t index = fIndexRecKinks[label*20+j];
1413         //AliESDkink *kink2  = (AliESDkink*)fKinks->At(index);
1414         AliESDkink *kink2  = (AliESDkink*)fEvent->GetKink(index);
1415         if (TMath::Abs(kink2->fLab[0])==label &&TMath::Abs(kink2->fLab[1])==label2) {
1416           fRecKinkInfo->fKinkMultiple++;
1417           fSignedKinks[index]=1;
1418           Int_t c0=0;
1419           if (kink){
1420             //      if (kink->fTRDOn) c0++;
1421             //if (kink->fITSOn) c0++;
1422             if (kink->GetStatus(2)>0) c0++;
1423             if (kink->GetStatus(0)>0) c0++;
1424           }
1425           Int_t c2=0;
1426           //      if (kink2->fTRDOn) c2++;
1427           //if (kink2->fITSOn) c2++;
1428           if (kink2->GetStatus(2)>0) c2++;
1429           if (kink2->GetStatus(0)>0) c2++;
1430
1431           if (c2<c0) continue;
1432           kink =kink2;
1433         }
1434         if (TMath::Abs(kink2->fLab[1])==label &&TMath::Abs(kink2->fLab[0])==label2) {
1435           fRecKinkInfo->fKinkMultiple++;
1436           fSignedKinks[index]=1;
1437           Int_t c0=0;
1438           if (kink){
1439             //if (kink->fTRDOn) c0++;
1440             //if (kink->fITSOn) c0++;
1441             if (kink->GetStatus(2)>0) c0++;
1442             if (kink->GetStatus(0)>0) c0++;
1443
1444           }
1445           Int_t c2=0;
1446           //      if (kink2->fTRDOn) c2++;
1447           //if (kink2->fITSOn) c2++;
1448           if (kink2->GetStatus(2)>0) c2++;
1449           if (kink2->GetStatus(0)>0) c2++;
1450
1451           if (c2<c0) continue;
1452           kink =kink2;
1453         }
1454       }
1455     }
1456     if (kink){
1457       fRecKinkInfo->fKink = *kink;
1458       fRecKinkInfo->fRecStatus=1;
1459     }
1460     fTreeCmpKinks->Fill();
1461   }
1462   //  Int_t nkinks = fKinks->GetEntriesFast();
1463   Int_t nkinks = fEvent->GetNumberOfKinks();
1464   for (Int_t i=0;i<nkinks;i++){
1465     if (fSignedKinks[i]==0){
1466       //      AliESDkink *kink  = (AliESDkink*)fKinks->At(i);
1467       AliESDkink *kink  = (AliESDkink*)fEvent->GetKink(i);
1468       if (!kink) continue;
1469       //
1470       fRecKinkInfo->fKink = *kink;
1471       fRecKinkInfo->fRecStatus =-2;
1472       //
1473       AliESDRecInfo*  fRecInfo1 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(kink->fLab[0]));
1474       AliESDRecInfo*  fRecInfo2 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(kink->fLab[1]));
1475       if (fRecInfo1 && fRecInfo2){
1476         fRecKinkInfo->fT1 = (*fRecInfo1);
1477         fRecKinkInfo->fT2 = (*fRecInfo2);
1478         fRecKinkInfo->fRecStatus =-1;
1479       }
1480       fTreeCmpKinks->Fill();
1481     }
1482   }
1483
1484
1485   fTreeCmpKinks->AutoSave();
1486   printf("Time spended in BuilKinkInfo Loop\n");
1487   timer.Print();
1488   if (fDebug > 2) cerr<<"end of BuildKinkInfo Loop"<<endl;
1489   return 0;
1490 }
1491
1492
1493
1494 void   ESDCmpTr::MakePoints(AliESDtrack * track, AliPointsMI &points)
1495 {
1496   //
1497   // make points in global coordinate frame
1498   //
1499   return;
1500   /*
1501   UInt_t itscl[10];
1502   Int_t nits = track->GetITSclusters(itscl);
1503   Int_t tpccl[1000];
1504   Int_t ntpc = track->GetTPcclusters(tpccl);
1505   UInt_t trdcl[1000];
1506   Int_t ntrd = track->GetTRDclusters(trdcl);
1507   //
1508   AliLoader *itsloader = fLoader->GetLoader("ITSLoader");
1509   AliLoader *tpcloader = fLoader->GetLoader("TPCLoader");
1510   AliLoader *trdloader = fLoader->GetLoader("TRDLoader");
1511   //
1512   AliITStrackerMI itstracker(); 
1513   */
1514 }
1515
1516
1517 ////////////////////////////////////////////////////////////////////////
1518 ////////////////////////////////////////////////////////////////////////
1519 ////////////////////////////////////////////////////////////////////////
1520
1521
1522
1523 Int_t ESDCmpTr::BuildV0Info(Int_t eventNr)
1524 {
1525 //
1526 // loop over all entries for a given event, find corresponding 
1527 // rec. track and store in the fTreeCmp
1528 //
1529   TStopwatch timer;
1530   timer.Start();
1531   Int_t entry = fNextV0ToRead;
1532   Double_t nParticlesTR = fTreeGenV0->GetEntriesFast();
1533   cerr<<"fNParticles, nParticlesTR, fNextV0ToRead: "<<fNParticles<<" "
1534       <<nParticlesTR<<" "<<fNextV0ToRead<<endl;
1535   //
1536   TBranch * branch = fTreeCmpV0->GetBranch("RC.");
1537   branch->SetAddress(&fRecV0Info); // set all pointers
1538   const AliESDVertex * esdvertex = fEvent->GetVertex();
1539   Float_t vertex[3]= {esdvertex->GetXv(), esdvertex->GetYv(),esdvertex->GetZv()};
1540   
1541   //
1542   while (entry < nParticlesTR) {
1543     fTreeGenV0->GetEntry(entry);
1544     entry++;
1545     if (eventNr < fGenV0Info->fMCm.fEventNr) continue;
1546     if (eventNr > fGenV0Info->fMCm.fEventNr) continue;;
1547     //
1548     fNextV0ToRead = entry-1;
1549     //
1550     //
1551     AliESDRecInfo*  fRecInfo1 = (AliESDRecInfo*)fRecArray->At(fGenV0Info->fMCm.fLabel);
1552     AliESDRecInfo*  fRecInfo2 = (AliESDRecInfo*)fRecArray->At(fGenV0Info->fMCd.fLabel);
1553     if (fGenV0Info->fMCm.fCharge*fGenV0Info->fMCd.fCharge>0) continue;  // interactions
1554     if (!fRecInfo1 || !fRecInfo2) continue;
1555     fRecV0Info->fT1 = (*fRecInfo1);
1556     fRecV0Info->fT2 = (*fRecInfo2);
1557     fRecV0Info->fV0Status =0;
1558     if (fRecInfo1 && fRecInfo1->fStatus[1]>0) fRecV0Info->fV0Status+=1;
1559     if (fRecInfo2 && fRecInfo2->fStatus[1]>0) fRecV0Info->fV0Status+=2;
1560
1561     if (fRecV0Info->fV0Status==3&&fRecInfo1->fSign==fRecInfo2->fSign) fRecV0Info->fV0Status*=-1;
1562
1563
1564     if (abs(fRecV0Info->fV0Status)==3){
1565       fRecV0Info->Update(vertex);
1566       {
1567         //
1568         // TPC V0 Info
1569         Double_t x,alpha, param[5],cov[15];
1570         fRecV0Info->fT1.fESDTrack.GetInnerExternalParameters(alpha,x,param);
1571         fRecV0Info->fT1.fESDTrack.GetInnerExternalCovariance(cov);
1572         AliExternalTrackParam paramP(x,alpha,param,cov);
1573         //
1574         fRecV0Info->fT2.fESDTrack.GetInnerExternalParameters(alpha,x,param);
1575         fRecV0Info->fT2.fESDTrack.GetInnerExternalCovariance(cov);
1576         AliExternalTrackParam paramM(x,alpha,param,cov);
1577         //
1578         fRecV0Info->fV0tpc.SetM(paramM);
1579         fRecV0Info->fV0tpc.SetP(paramP);
1580         Double_t pid1[5],pid2[5];
1581         fRecV0Info->fT1.fESDTrack.GetESDpid(pid1);
1582         fRecV0Info->fT1.fESDTrack.GetESDpid(pid2);
1583         //
1584         fRecV0Info->fV0tpc.UpdatePID(pid1,pid2);
1585         fRecV0Info->fV0tpc.Update(vertex);
1586         //
1587         //
1588         fRecV0Info->fT1.fESDTrack.GetExternalParameters(x,param);
1589         fRecV0Info->fT1.fESDTrack.GetExternalCovariance(cov);
1590         alpha = fRecV0Info->fT1.fESDTrack.GetAlpha();
1591         new (&paramP) AliExternalTrackParam(x,alpha,param,cov);
1592         //
1593         fRecV0Info->fT2.fESDTrack.GetExternalParameters(x,param);
1594         fRecV0Info->fT2.fESDTrack.GetExternalCovariance(cov);
1595         alpha = fRecV0Info->fT2.fESDTrack.GetAlpha();
1596         new (&paramM) AliExternalTrackParam(x,alpha,param,cov);
1597         //
1598         fRecV0Info->fV0its.SetM(paramM);
1599         fRecV0Info->fV0its.SetP(paramP);
1600         fRecV0Info->fV0its.UpdatePID(pid1,pid2);
1601         fRecV0Info->fV0its.Update(vertex);
1602
1603       }
1604       if (TMath::Abs(fGenV0Info->fMCm.fPdg)==11 &&TMath::Abs(fGenV0Info->fMCd.fPdg)==11){
1605         if (fRecV0Info->fDist2>10){
1606           fRecV0Info->Update(vertex);
1607         }
1608         if (fRecV0Info->fDist2>10){
1609           fRecV0Info->Update(vertex);
1610         }
1611       }
1612     }   
1613     //
1614     // take the V0 from reconstruction
1615  
1616     Int_t label =  TMath::Min(fGenV0Info->fMCm.fLabel,fGenV0Info->fMCd.fLabel);
1617     Int_t label2 = TMath::Max(fGenV0Info->fMCm.fLabel,fGenV0Info->fMCd.fLabel);    
1618     AliESDV0MI *v0MI=0;
1619     fRecV0Info->fRecStatus   =0;
1620     fRecV0Info->fMultiple    = fMultiRecV0[label];
1621     fRecV0Info->fV0Multiple=0;
1622     //
1623     if (fMultiRecV0[label]>0 || fMultiRecV0[label2]>0){
1624
1625       //      for (Int_t j=0;j<TMath::Min(fMultiRecV0s[label],100);j++){
1626       for (Int_t j=TMath::Min(fMultiRecV0[label],Short_t(20))-1;j>=0;j--){
1627         Int_t index = fIndexRecV0[label*20+j];
1628         if (index<0) continue;
1629         AliESDV0MI *v0MI2  = fEvent->GetV0MI(index);
1630         if (TMath::Abs(v0MI2->fLab[0])==label &&TMath::Abs(v0MI2->fLab[1])==label2) {
1631           v0MI =v0MI2;
1632           fRecV0Info->fV0Multiple++;
1633           fSignedV0[index]=1;
1634         }
1635         if (TMath::Abs(v0MI2->fLab[1])==label &&TMath::Abs(v0MI2->fLab[0])==label2) {
1636           v0MI =v0MI2;
1637           fRecV0Info->fV0Multiple++;
1638           fSignedV0[index]=1;
1639         }
1640       }
1641     }
1642     if (v0MI){
1643       fRecV0Info->fV0rec = *v0MI;
1644       fRecV0Info->fRecStatus=1;
1645     }
1646
1647     fTreeCmpV0->Fill();
1648   }
1649   //
1650   // write fake v0s
1651
1652   Int_t nV0MIs = fEvent->GetNumberOfV0MIs();
1653   for (Int_t i=0;i<nV0MIs;i++){
1654     if (fSignedV0[i]==0){
1655       AliESDV0MI *v0MI  = fEvent->GetV0MI(i);
1656       if (!v0MI) continue;
1657       //
1658       fRecV0Info->fV0rec = *v0MI;
1659       fRecV0Info->fV0Status  =-10;
1660       fRecV0Info->fRecStatus =-2;
1661       //
1662       AliESDRecInfo*  fRecInfo1 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(v0MI->fLab[0]));
1663       AliESDRecInfo*  fRecInfo2 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(v0MI->fLab[1]));
1664       if (fRecInfo1 && fRecInfo2){
1665         fRecV0Info->fT1 = (*fRecInfo1);
1666         fRecV0Info->fT2 = (*fRecInfo2);
1667         fRecV0Info->fRecStatus =-1;
1668       }
1669       fRecV0Info->Update(vertex);
1670       fTreeCmpV0->Fill();
1671     }
1672   }
1673
1674
1675
1676   fTreeCmpV0->AutoSave();
1677   fRecArray->Delete();
1678   printf("Time spended in BuilV0Info Loop\n");
1679   timer.Print();
1680   if (fDebug > 2) cerr<<"end of BuildV0Info Loop"<<endl;
1681   return 0;
1682 }
1683 ////////////////////////////////////////////////////////////////////////
1684 ////////////////////////////////////////////////////////////////////////
1685
1686 void AliESDComparisonDraw::SetIO(const char *fname)
1687 {
1688   //
1689    TFile* file = TFile::Open(fname);
1690    //
1691    fTree = (TTree*) file->Get("ESDcmpTracks");
1692    if (!fTree) {
1693     printf("no track comparison tree found\n");
1694     file->Close();
1695     delete file;
1696   }
1697 }
1698
1699