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