PWG1/AliESDComparisonMI.cxx

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