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