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