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 "AliMagF.h"
 173 #include "AliESD.h"
 174 #include "AliESDtrack.h"
 175 #include "AliITStrackMI.h"
 176 #include "AliTRDtrack.h"
 177 #include "AliHelix.h"
 178 #include "AliESDVertex.h"
 179 #include "AliExternalTrackParam.h"
 180 #include "AliESDkink.h"
 181 #include "AliESDV0MI.h"
 182
 183 #endif
 184 #include "AliGenInfo.h"
 185 #include "AliESDComparisonMI.h"
 186
 187
 188
 189
 190
 191 void MakeAliases(AliESDComparisonDraw&comp)
 192 {
 193   //
 194   // aliases definition
 195   //
 196   comp.fTree->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)");
 197   comp.fTree->SetAlias("direction","MC.fParticle.fVx*MC.fParticle.fPx+MC.fParticle.fVy*MC.fParticle.fPy");
 198   comp.fTree->SetAlias("decaydir","MC.fTRdecay.fX*MC.fTRdecay.fPx+MC.fTRdecay.fY*MC.fTRdecay.fPy");
 199   comp.fTree->SetAlias("theta","MC.fTrackRef.Theta()");
 200   comp.fTree->SetAlias("primdca","sqrt(RC.fITStrack.fD[0]**2+RC.fITStrack.fD[1]**2)");
 201   comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
 202   comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
 203
 204   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.");
 205
 206   comp.fTree->SetAlias("dtofmc2","fESDTrack.fTrackTime[2]-(10^12*MC.fTOFReferences[0].fTime)");
 207   comp.fTree->SetAlias("dtofrc2","(fESDTrack.fTrackTime[2]-fESDTrack.fTOFsignal)");
 208
 209   comp.fTree->SetAlias("psum","fESDTrack.fTOFr[4]+fESDTrack.fTOFr[3]+fESDTrack.fTOFr[2]+fESDTrack.fTOFr[1]+fESDTrack.fTOFr[0]");
 210   comp.fTree->SetAlias("P0","fESDTrack.fTOFr[0]/psum");
 211   comp.fTree->SetAlias("P1","fESDTrack.fTOFr[1]/psum");
 212   comp.fTree->SetAlias("P2","fESDTrack.fTOFr[2]/psum");
 213   comp.fTree->SetAlias("P3","fESDTrack.fTOFr[3]/psum");
 214   comp.fTree->SetAlias("P4","fESDTrack.fTOFr[4]/psum");
 215   comp.fTree->SetAlias("MaxP","max(max(max(P0,P1),max(P2,P3)),P4)");
 216 }
 217
 218
 219
 220
 221 void  AliESDRecInfo::UpdatePoints(AliESDtrack*track)
 222 {
 223   //
 224   //
 225   Int_t iclusters[200];
 226   Float_t density[160];
 227   for (Int_t i=0;i<160;i++) density[i]=-1.;
 228   fTPCPoints[0]= 160;
 229   fTPCPoints[1] = -1;
 230   //
 231   if (fTPCPoints[0]<fTPCPoints[1]) return;
 232   //  Int_t nclusters=track->GetTPCclusters(iclusters);
 233
 234   Int_t ngood=0;
 235   Int_t undeff=0;
 236   Int_t nall =0;
 237   Int_t range=20;
 238   for (Int_t i=0;i<160;i++){
 239     Int_t last = i-range;
 240     if (nall<range) nall++;
 241     if (last>=0){
 242       if (iclusters[last]>0&& (iclusters[last]&0x8000)==0) ngood--;
 243       if (iclusters[last]==-1) undeff--;
 244     }
 245     if (iclusters[i]>0&& (iclusters[i]&0x8000)==0)   ngood++;
 246     if (iclusters[i]==-1) undeff++;
 247     if (nall==range &&undeff<range/2) density[i-range/2] = Float_t(ngood)/Float_t(nall-undeff);
 248   }
 249   Float_t maxdens=0;
 250   Int_t indexmax =0;
 251   for (Int_t i=0;i<160;i++){
 252     if (density[i]<0) continue;
 253     if (density[i]>maxdens){
 254       maxdens=density[i];
 255       indexmax=i;
 256     }
 257   }
 258   //
 259   //max dens point
 260   fTPCPoints[3] = maxdens;
 261   fTPCPoints[1] = indexmax;
 262   //
 263   // last point
 264   for (Int_t i=indexmax;i<160;i++){
 265     if (density[i]<0) continue;
 266     if (density[i]<maxdens/2.) {
 267       break;
 268     }
 269     fTPCPoints[2]=i;
 270   }
 271   //
 272   // first point
 273   for (Int_t i=indexmax;i>0;i--){
 274     if (density[i]<0) continue;
 275     if (density[i]<maxdens/2.) {
 276       break;
 277     }
 278     fTPCPoints[0]=i;
 279   }
 280   //
 281   // Density at the last 30 padrows
 282   //
 283   //
 284   nall  = 0;
 285   ngood = 0;
 286   for (Int_t i=159;i>0;i--){
 287     if (iclusters[i]==-1) continue; //dead zone
 288     nall++;
 289     if (iclusters[i]>0)   ngood++;
 290     if (nall>20) break;
 291   }
 292   fTPCPoints[4] = Float_t(ngood)/Float_t(nall);
 293   //
 294   if ((track->GetStatus()&AliESDtrack::kITSrefit)>0) fTPCPoints[0]=-1;
 295   //
 296   //
 297   // check TRDPoints
 298   /*
 299   nclusters=track->GetTRDclusters(iclusters);
 300   for (Int_t i=nclusters;i>0;i--){
 301
 302   }
 303   */
 304
 305
 306 }
 307
 308 //
 309 //
 310 void AliESDRecInfo::Update(AliMCInfo* info,AliTPCParam * /*par*/, Bool_t reconstructed, AliESD *event)
 311 {
 312   //
 313   //
 314   //calculates derived variables
 315   //
 316   //
 317   UpdatePoints(&fESDTrack);
 318   fBestTOFmatch=1000;
 319   AliTrackReference * ref = &(info->fTrackRef);
 320   fTPCinR0[0] = info->fTrackRef.X();
 321   fTPCinR0[1] = info->fTrackRef.Y();
 322   fTPCinR0[2] = info->fTrackRef.Z();
 323   fTPCinR0[3] = TMath::Sqrt(fTPCinR0[0]*fTPCinR0[0]+fTPCinR0[1]*fTPCinR0[1]);
 324   fTPCinR0[4] = TMath::ATan2(fTPCinR0[1],fTPCinR0[0]);
 325   //
 326   fTPCinP0[0] = ref->Px();
 327   fTPCinP0[1] = ref->Py();
 328   fTPCinP0[2] = ref->Pz();
 329   fTPCinP0[3] = ref->Pt();
 330   fTPCinP0[4] = ref->P();
 331   fDeltaP     = (ref->P()-info->fParticle.P())/info->fParticle.P();
 332   //
 333   //
 334   if (fTPCinP0[3]>0.0000001){
 335     //
 336     fTPCAngle0[0] = TMath::ATan2(fTPCinP0[1],fTPCinP0[0]);
 337     fTPCAngle0[1] = TMath::ATan(fTPCinP0[2]/fTPCinP0[3]);
 338   }
 339   //
 340   //
 341   fITSinP0[0]=info->fParticle.Px();
 342   fITSinP0[1]=info->fParticle.Py();
 343   fITSinP0[2]=info->fParticle.Pz();
 344   fITSinP0[3]=info->fParticle.Pt();
 345   //
 346   fITSinR0[0]=info->fParticle.Vx();
 347   fITSinR0[1]=info->fParticle.Vy();
 348   fITSinR0[2]=info->fParticle.Vz();
 349   fITSinR0[3] = TMath::Sqrt(fITSinR0[0]*fITSinR0[0]+fITSinR0[1]*fITSinR0[1]);
 350   fITSinR0[4] = TMath::ATan2(fITSinR0[1],fITSinR0[0]);
 351   //
 352   //
 353   if (fITSinP0[3]>0.0000001){
 354     fITSAngle0[0] = TMath::ATan2(fITSinP0[1],fITSinP0[0]);
 355     fITSAngle0[1] = TMath::ATan(fITSinP0[2]/fITSinP0[3]);
 356   }
 357   //
 358   for (Int_t i=0;i<4;i++) fStatus[i] =0;
 359   fReconstructed = kFALSE;
 360   fTPCOn = kFALSE;
 361   fITSOn = kFALSE;
 362   fTRDOn = kFALSE;
 363   if (reconstructed==kFALSE) return;
 364
 365   fLabels[0] = info->fLabel;
 366   fLabels[1] = info->fPrimPart;
 367   fReconstructed = kTRUE;
 368   fTPCOn = ((fESDTrack.GetStatus()&AliESDtrack::kTPCrefit)>0) ? kTRUE : kFALSE;
 369   fITSOn = ((fESDTrack.GetStatus()&AliESDtrack::kITSrefit)>0) ? kTRUE : kFALSE;
 370   fTRDOn = ((fESDTrack.GetStatus()&AliESDtrack::kTRDrefit)>0) ? kTRUE : kFALSE;
 371   //
 372   //
 373   if ((fESDTrack.GetStatus()&AliESDtrack::kTPCrefit)>0){
 374     fStatus[1] =3;
 375   }
 376   else{
 377     if ((fESDTrack.GetStatus()&AliESDtrack::kTPCout)>0){
 378       fStatus[1] =2;
 379     }
 380     else{
 381       if ((fESDTrack.GetStatus()&AliESDtrack::kTPCin)>0)
 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.000000000000001){
 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, alpha;
 457     fESDTrack.GetInnerExternalCovariance(cov);
 458     fESDTrack.GetInnerExternalParameters(alpha, 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->Transform1to2Ideal(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   AliESDtrack dummytrack;  //
1255
1256   while (entry < nParticlesTR) {
1257     fTreeGenTracks->GetEntry(entry);
1258     entry++;
1259     if (eventNr < fMCInfo->fEventNr) continue;
1260     if (eventNr > fMCInfo->fEventNr) continue;;
1261     //
1262     fNextTreeGenEntryToRead = entry-1;
1263     if (fDebug > 2 && fMCInfo->fLabel < 10) {
1264       cerr<<"Fill track with a label "<<fMCInfo->fLabel<<endl;
1265     }
1266     //    if (fMCInfo->fNTPCRef<1) continue; // not TPCref
1267     //
1268     fRecInfo->Reset();
1269     AliESDtrack * track=0;
1270     fRecInfo->fReconstructed =0;
1271     TVector3 local = TR2Local(&(fMCInfo->fTrackRef),fParamTPC);
1272     local.GetXYZ(fRecInfo->fTRLocalCoord);
1273     //
1274     if (fIndexRecTracks[fMCInfo->fLabel*20] >= 0) {
1275       //track= (AliESDtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4]);
1276       track= (AliESDtrack*)fEvent->GetTrack(fIndexRecTracks[fMCInfo->fLabel*20]);
1277       //
1278       //
1279       // find nearest track if multifound
1280       //Int_t sign = Int_t(track->GetSign()*fMCInfo->fCharge);
1281       //
1282       Int_t status = 0;
1283       if  ((track->GetStatus()&AliESDtrack::kITSrefit)>0) status++;
1284       if  ((track->GetStatus()&AliESDtrack::kTPCrefit)>0) status++;
1285       if  ((track->GetStatus()&AliESDtrack::kTRDrefit)>0) status++;
1286
1287       //
1288       if (fIndexRecTracks[fMCInfo->fLabel*20+1]>0){
1289         //
1290         Double_t p[3];
1291         track->GetInnerPxPyPz(p);
1292         Float_t maxp = p[0]*p[0]+p[1]*p[1]+p[2]*p[2];
1293         //
1294         for (Int_t i=1;i<20;i++){
1295           if (fIndexRecTracks[fMCInfo->fLabel*20+i]>=0){
1296             AliESDtrack * track2 = (AliESDtrack*)fEvent->GetTrack(fIndexRecTracks[fMCInfo->fLabel*20+i]);
1297             if (!track2) continue;
1298             //Int_t sign2 = track->GetSign()*fMCInfo->fCharge; //
1299             //if (sign2<0) continue;
1300             track2->GetInnerPxPyPz(p);
1301             Float_t mom = p[0]*p[0]+p[1]*p[1]+p[2]*p[2];
1302             /*
1303             if (sign<0){
1304               sign = sign2;
1305               track = track2;
1306               maxp = mom;
1307               continue;
1308             }
1309             */
1310             //
1311             Int_t status2 = 0;
1312             if  ((track2->GetStatus()&AliESDtrack::kITSrefit)>0) status2++;
1313             if  ((track2->GetStatus()&AliESDtrack::kTPCrefit)>0) status2++;
1314             if  ((track2->GetStatus()&AliESDtrack::kTRDrefit)>0) status2++;
1315             if (status2<status) continue;
1316             //
1317             if (mom<maxp) continue;
1318             maxp = mom;
1319             track = track2;
1320             //
1321           }
1322         }
1323       }
1324       //
1325       if (track) {
1326         new (&(fRecInfo->fESDTrack)) AliESDtrack(*track);
1327       }else{
1328         fRecInfo->fESDTrack = dummytrack;
1329       }
1330
1331       if (track->GetITStrack())
1332         fRecInfo->fITStrack = *((AliITStrackMI*)track->GetITStrack());
1333       else{
1334         fRecInfo->fITStrack = *track;
1335       }
1336       if (track->GetTRDtrack()){
1337         fRecInfo->fTRDtrack = *((AliTRDtrack*)track->GetTRDtrack());
1338       }
1339       else{
1340         fRecInfo->fTRDtrack.SetdEdx(-1);
1341       }
1342       fRecInfo->fReconstructed = 1;
1343       fRecInfo->fFake     = fFakeRecTracks[fMCInfo->fLabel];
1344       fRecInfo->fMultiple = fMultiRecTracks[fMCInfo->fLabel];
1345       //
1346       fRecInfo->Update(fMCInfo,fParamTPC,kTRUE, fEvent);
1347     }
1348     else{
1349       fRecInfo->fESDTrack = dummytrack;
1350       fRecInfo->Update(fMCInfo,fParamTPC,kFALSE, fEvent);
1351     }
1352     fRecArray->AddAt(new AliESDRecInfo(*fRecInfo),fMCInfo->fLabel);
1353     fTreeCmp->Fill();
1354   }
1355   fTreeCmp->AutoSave();
1356   //fTracks->Delete();
1357   printf("Time spended in TreeGenLoop\n");
1358   timer.Print();
1359   if (fDebug > 2) cerr<<"end of TreeGenLoop"<<endl;
1360
1361   return 0;
1362 }
1363
1364
1365
1366 ////////////////////////////////////////////////////////////////////////
1367 ////////////////////////////////////////////////////////////////////////
1368 ////////////////////////////////////////////////////////////////////////
1369 Int_t ESDCmpTr::BuildKinkInfo0(Int_t eventNr)
1370 {
1371 //
1372 // loop over all entries for a given event, find corresponding
1373 // rec. track and store in the fTreeCmp
1374 //
1375   TStopwatch timer;
1376   timer.Start();
1377   Int_t entry = fNextKinkToRead;
1378   Double_t nParticlesTR = fTreeGenKinks->GetEntriesFast();
1379   cerr<<"fNParticles, nParticlesTR, fNextKinkToRead: "<<fNParticles<<" "
1380       <<nParticlesTR<<" "<<fNextKinkToRead<<endl;
1381   //
1382   TBranch * branch = fTreeCmpKinks->GetBranch("RC.");
1383   branch->SetAddress(&fRecKinkInfo); // set all pointers
1384
1385   //
1386   while (entry < nParticlesTR) {
1387     fTreeGenKinks->GetEntry(entry);
1388     entry++;
1389     if (eventNr < fGenKinkInfo->fMCm.fEventNr) continue;
1390     if (eventNr > fGenKinkInfo->fMCm.fEventNr) continue;;
1391     //
1392     fNextKinkToRead = entry-1;
1393     //
1394     //
1395     AliESDRecInfo*  fRecInfo1 = (AliESDRecInfo*)fRecArray->At(fGenKinkInfo->fMCm.fLabel);
1396     AliESDRecInfo*  fRecInfo2 = (AliESDRecInfo*)fRecArray->At(fGenKinkInfo->fMCd.fLabel);
1397     fRecKinkInfo->fT1 = (*fRecInfo1);
1398     fRecKinkInfo->fT2 = (*fRecInfo2);
1399     fRecKinkInfo->fStatus =0;
1400     if (fRecInfo1 && fRecInfo1->fTPCOn) fRecKinkInfo->fStatus+=1;
1401     if (fRecInfo2 && fRecInfo2->fTPCOn) fRecKinkInfo->fStatus+=2;
1402     if (fRecKinkInfo->fStatus==3&&fRecInfo1->fSign!=fRecInfo2->fSign) fRecKinkInfo->fStatus*=-1;
1403
1404     if (fRecKinkInfo->fStatus==3){
1405       fRecKinkInfo->Update();
1406     }
1407     Int_t label =  TMath::Min(fGenKinkInfo->fMCm.fLabel,fGenKinkInfo->fMCd.fLabel);
1408     Int_t label2 = TMath::Max(fGenKinkInfo->fMCm.fLabel,fGenKinkInfo->fMCd.fLabel);
1409
1410     AliESDkink *kink=0;
1411     fRecKinkInfo->fRecStatus   =0;
1412     fRecKinkInfo->fMultiple    = fMultiRecKinks[label];
1413     fRecKinkInfo->fKinkMultiple=0;
1414     //
1415     if (fMultiRecKinks[label]>0){
1416
1417       //      for (Int_t j=0;j<TMath::Min(fMultiRecKinks[label],100);j++){
1418       for (Int_t j=TMath::Min(fMultiRecKinks[label],Short_t(20))-1;j>=0;j--){
1419         Int_t index = fIndexRecKinks[label*20+j];
1420         //AliESDkink *kink2  = (AliESDkink*)fKinks->At(index);
1421         AliESDkink *kink2  = (AliESDkink*)fEvent->GetKink(index);
1422         if (TMath::Abs(kink2->fLab[0])==label &&TMath::Abs(kink2->fLab[1])==label2) {
1423           fRecKinkInfo->fKinkMultiple++;
1424           fSignedKinks[index]=1;
1425           Int_t c0=0;
1426           if (kink){
1427             //      if (kink->fTRDOn) c0++;
1428             //if (kink->fITSOn) c0++;
1429             if (kink->GetStatus(2)>0) c0++;
1430             if (kink->GetStatus(0)>0) c0++;
1431           }
1432           Int_t c2=0;
1433           //      if (kink2->fTRDOn) c2++;
1434           //if (kink2->fITSOn) c2++;
1435           if (kink2->GetStatus(2)>0) c2++;
1436           if (kink2->GetStatus(0)>0) c2++;
1437
1438           if (c2<c0) continue;
1439           kink =kink2;
1440         }
1441         if (TMath::Abs(kink2->fLab[1])==label &&TMath::Abs(kink2->fLab[0])==label2) {
1442           fRecKinkInfo->fKinkMultiple++;
1443           fSignedKinks[index]=1;
1444           Int_t c0=0;
1445           if (kink){
1446             //if (kink->fTRDOn) c0++;
1447             //if (kink->fITSOn) c0++;
1448             if (kink->GetStatus(2)>0) c0++;
1449             if (kink->GetStatus(0)>0) c0++;
1450
1451           }
1452           Int_t c2=0;
1453           //      if (kink2->fTRDOn) c2++;
1454           //if (kink2->fITSOn) c2++;
1455           if (kink2->GetStatus(2)>0) c2++;
1456           if (kink2->GetStatus(0)>0) c2++;
1457
1458           if (c2<c0) continue;
1459           kink =kink2;
1460         }
1461       }
1462     }
1463     if (kink){
1464       fRecKinkInfo->fKink = *kink;
1465       fRecKinkInfo->fRecStatus=1;
1466     }
1467     fTreeCmpKinks->Fill();
1468   }
1469   //  Int_t nkinks = fKinks->GetEntriesFast();
1470   Int_t nkinks = fEvent->GetNumberOfKinks();
1471   for (Int_t i=0;i<nkinks;i++){
1472     if (fSignedKinks[i]==0){
1473       //      AliESDkink *kink  = (AliESDkink*)fKinks->At(i);
1474       AliESDkink *kink  = (AliESDkink*)fEvent->GetKink(i);
1475       if (!kink) continue;
1476       //
1477       fRecKinkInfo->fKink = *kink;
1478       fRecKinkInfo->fRecStatus =-2;
1479       //
1480       AliESDRecInfo*  fRecInfo1 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(kink->fLab[0]));
1481       AliESDRecInfo*  fRecInfo2 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(kink->fLab[1]));
1482       if (fRecInfo1 && fRecInfo2){
1483         fRecKinkInfo->fT1 = (*fRecInfo1);
1484         fRecKinkInfo->fT2 = (*fRecInfo2);
1485         fRecKinkInfo->fRecStatus =-1;
1486       }
1487       fTreeCmpKinks->Fill();
1488     }
1489   }
1490
1491
1492   fTreeCmpKinks->AutoSave();
1493   printf("Time spended in BuilKinkInfo Loop\n");
1494   timer.Print();
1495   if (fDebug > 2) cerr<<"end of BuildKinkInfo Loop"<<endl;
1496   return 0;
1497 }
1498
1499
1500
1501 void   ESDCmpTr::MakePoints(AliESDtrack * track, AliPointsMI &points)
1502 {
1503   //
1504   // make points in global coordinate frame
1505   //
1506   return;
1507   /*
1508   UInt_t itscl[10];
1509   Int_t nits = track->GetITSclusters(itscl);
1510   Int_t tpccl[1000];
1511   Int_t ntpc = track->GetTPcclusters(tpccl);
1512   UInt_t trdcl[1000];
1513   Int_t ntrd = track->GetTRDclusters(trdcl);
1514   //
1515   AliLoader *itsloader = fLoader->GetLoader("ITSLoader");
1516   AliLoader *tpcloader = fLoader->GetLoader("TPCLoader");
1517   AliLoader *trdloader = fLoader->GetLoader("TRDLoader");
1518   //
1519   AliITStrackerMI itstracker();
1520   */
1521 }
1522
1523
1524 ////////////////////////////////////////////////////////////////////////
1525 ////////////////////////////////////////////////////////////////////////
1526 ////////////////////////////////////////////////////////////////////////
1527
1528
1529
1530 Int_t ESDCmpTr::BuildV0Info(Int_t eventNr)
1531 {
1532 //
1533 // loop over all entries for a given event, find corresponding
1534 // rec. track and store in the fTreeCmp
1535 //
1536   TStopwatch timer;
1537   timer.Start();
1538   Int_t entry = fNextV0ToRead;
1539   Double_t nParticlesTR = fTreeGenV0->GetEntriesFast();
1540   cerr<<"fNParticles, nParticlesTR, fNextV0ToRead: "<<fNParticles<<" "
1541       <<nParticlesTR<<" "<<fNextV0ToRead<<endl;
1542   //
1543   TBranch * branch = fTreeCmpV0->GetBranch("RC.");
1544   branch->SetAddress(&fRecV0Info); // set all pointers
1545   const AliESDVertex * esdvertex = fEvent->GetVertex();
1546   Float_t vertex[3]= {esdvertex->GetXv(), esdvertex->GetYv(),esdvertex->GetZv()};
1547
1548   //
1549   while (entry < nParticlesTR) {
1550     fTreeGenV0->GetEntry(entry);
1551     entry++;
1552     if (eventNr < fGenV0Info->fMCm.fEventNr) continue;
1553     if (eventNr > fGenV0Info->fMCm.fEventNr) continue;;
1554     //
1555     fNextV0ToRead = entry-1;
1556     //
1557     //
1558     AliESDRecInfo*  fRecInfo1 = (AliESDRecInfo*)fRecArray->At(fGenV0Info->fMCm.fLabel);
1559     AliESDRecInfo*  fRecInfo2 = (AliESDRecInfo*)fRecArray->At(fGenV0Info->fMCd.fLabel);
1560     if (fGenV0Info->fMCm.fCharge*fGenV0Info->fMCd.fCharge>0) continue;  // interactions
1561     if (!fRecInfo1 || !fRecInfo2) continue;
1562     fRecV0Info->fT1 = (*fRecInfo1);
1563     fRecV0Info->fT2 = (*fRecInfo2);
1564     fRecV0Info->fV0Status =0;
1565     if (fRecInfo1 && fRecInfo1->fStatus[1]>0) fRecV0Info->fV0Status+=1;
1566     if (fRecInfo2 && fRecInfo2->fStatus[1]>0) fRecV0Info->fV0Status+=2;
1567
1568     if (fRecV0Info->fV0Status==3&&fRecInfo1->fSign==fRecInfo2->fSign) fRecV0Info->fV0Status*=-1;
1569
1570
1571     if (abs(fRecV0Info->fV0Status)==3){
1572       fRecV0Info->Update(vertex);
1573       {
1574         //
1575         // TPC V0 Info
1576         Double_t x,alpha, param[5],cov[15];
1577         if ( fRecV0Info->fT1.fESDTrack.GetInnerParam() && fRecV0Info->fT2.fESDTrack.GetInnerParam()){
1578           fRecV0Info->fT1.fESDTrack.GetInnerExternalParameters(alpha,x,param);
1579           fRecV0Info->fT1.fESDTrack.GetInnerExternalCovariance(cov);
1580           AliExternalTrackParam paramP(x,alpha,param,cov);
1581           //
1582           fRecV0Info->fT2.fESDTrack.GetInnerExternalParameters(alpha,x,param);
1583           fRecV0Info->fT2.fESDTrack.GetInnerExternalCovariance(cov);
1584           AliExternalTrackParam paramM(x,alpha,param,cov);
1585           //
1586           fRecV0Info->fV0tpc.SetM(paramM);
1587           fRecV0Info->fV0tpc.SetP(paramP);
1588           Double_t pid1[5],pid2[5];
1589           fRecV0Info->fT1.fESDTrack.GetESDpid(pid1);
1590           fRecV0Info->fT1.fESDTrack.GetESDpid(pid2);
1591           //
1592           fRecV0Info->fV0tpc.UpdatePID(pid1,pid2);
1593           fRecV0Info->fV0tpc.Update(vertex);
1594
1595           //
1596           //
1597           fRecV0Info->fT1.fESDTrack.GetExternalParameters(x,param);
1598           fRecV0Info->fT1.fESDTrack.GetExternalCovariance(cov);
1599           alpha = fRecV0Info->fT1.fESDTrack.GetAlpha();
1600           new (&paramP) AliExternalTrackParam(x,alpha,param,cov);
1601           //
1602           fRecV0Info->fT2.fESDTrack.GetExternalParameters(x,param);
1603           fRecV0Info->fT2.fESDTrack.GetExternalCovariance(cov);
1604           alpha = fRecV0Info->fT2.fESDTrack.GetAlpha();
1605           new (&paramM) AliExternalTrackParam(x,alpha,param,cov);
1606           //
1607           fRecV0Info->fV0its.SetM(paramM);
1608           fRecV0Info->fV0its.SetP(paramP);
1609           fRecV0Info->fV0its.UpdatePID(pid1,pid2);
1610           fRecV0Info->fV0its.Update(vertex);
1611         }
1612       }
1613       if (TMath::Abs(fGenV0Info->fMCm.fPdg)==11 &&TMath::Abs(fGenV0Info->fMCd.fPdg)==11){
1614         if (fRecV0Info->fDist2>10){
1615           fRecV0Info->Update(vertex);
1616         }
1617         if (fRecV0Info->fDist2>10){
1618           fRecV0Info->Update(vertex);
1619         }
1620       }
1621     }
1622     //
1623     // take the V0 from reconstruction
1624
1625     Int_t label =  TMath::Min(fGenV0Info->fMCm.fLabel,fGenV0Info->fMCd.fLabel);
1626     Int_t label2 = TMath::Max(fGenV0Info->fMCm.fLabel,fGenV0Info->fMCd.fLabel);
1627     AliESDV0MI *v0MI=0;
1628     fRecV0Info->fRecStatus   =0;
1629     fRecV0Info->fMultiple    = fMultiRecV0[label];
1630     fRecV0Info->fV0Multiple=0;
1631     //
1632     if (fMultiRecV0[label]>0 || fMultiRecV0[label2]>0){
1633
1634       //      for (Int_t j=0;j<TMath::Min(fMultiRecV0s[label],100);j++){
1635       for (Int_t j=TMath::Min(fMultiRecV0[label],Short_t(20))-1;j>=0;j--){
1636         Int_t index = fIndexRecV0[label*20+j];
1637         if (index<0) continue;
1638         AliESDV0MI *v0MI2  = fEvent->GetV0MI(index);
1639         if (TMath::Abs(v0MI2->fLab[0])==label &&TMath::Abs(v0MI2->fLab[1])==label2) {
1640           v0MI =v0MI2;
1641           fRecV0Info->fV0Multiple++;
1642           fSignedV0[index]=1;
1643         }
1644         if (TMath::Abs(v0MI2->fLab[1])==label &&TMath::Abs(v0MI2->fLab[0])==label2) {
1645           v0MI =v0MI2;
1646           fRecV0Info->fV0Multiple++;
1647           fSignedV0[index]=1;
1648         }
1649       }
1650     }
1651     if (v0MI){
1652       fRecV0Info->fV0rec = *v0MI;
1653       fRecV0Info->fRecStatus=1;
1654     }
1655
1656     fTreeCmpV0->Fill();
1657   }
1658   //
1659   // write fake v0s
1660
1661   Int_t nV0MIs = fEvent->GetNumberOfV0MIs();
1662   for (Int_t i=0;i<nV0MIs;i++){
1663     if (fSignedV0[i]==0){
1664       AliESDV0MI *v0MI  = fEvent->GetV0MI(i);
1665       if (!v0MI) continue;
1666       //
1667       fRecV0Info->fV0rec = *v0MI;
1668       fRecV0Info->fV0Status  =-10;
1669       fRecV0Info->fRecStatus =-2;
1670       //
1671       AliESDRecInfo*  fRecInfo1 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(v0MI->fLab[0]));
1672       AliESDRecInfo*  fRecInfo2 = (AliESDRecInfo*)fRecArray->At(TMath::Abs(v0MI->fLab[1]));
1673       if (fRecInfo1 && fRecInfo2){
1674         fRecV0Info->fT1 = (*fRecInfo1);
1675         fRecV0Info->fT2 = (*fRecInfo2);
1676         fRecV0Info->fRecStatus =-1;
1677       }
1678       fRecV0Info->Update(vertex);
1679       fTreeCmpV0->Fill();
1680     }
1681   }
1682
1683
1684
1685   fTreeCmpV0->AutoSave();
1686   printf("Time spended in BuilV0Info Loop\n");
1687   timer.Print();
1688   if (fDebug > 2) cerr<<"end of BuildV0Info Loop"<<endl;
1689   return 0;
1690 }
1691 ////////////////////////////////////////////////////////////////////////
1692 ////////////////////////////////////////////////////////////////////////
1693
1694 void AliESDComparisonDraw::SetIO(const char *fname)
1695 {
1696   //
1697    TFile* file = TFile::Open(fname);
1698    //
1699    fTree = (TTree*) file->Get("ESDcmpTracks");
1700    if (!fTree) {
1701     printf("no track comparison tree found\n");
1702     file->Close();
1703     delete file;
1704   }
1705 }
1706
1707