STEER/AliESDComparisonMI.C

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