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