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.00005&&abs(MC.fVDist[2])<0.0005")
  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
  45 TCut csec("csec","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)>0.5");
  46
  47
  48 TCut crec("crec","fReconstructed==1");
  49 TCut cteta1("cteta1","abs(MC.fTrackRef.Theta()/3.1415-0.5)<0.25");
  50 TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1");
  51 TCut csens("csens","abs(sqrt(fVDist[0]**2+fVDist[1]**2)-170)<50");
  52 TCut cmuon("cmuon","abs(MC.fParticle.fPdgCode==-13)");
  53 TCut cchi2("cchi2","fESDTrack.fITSchi2MIP[0]<7.&&fESDTrack.fITSchi2MIP[1]<5.&&fESDTrack.fITSchi2MIP[2]<7.&&fESDTrack.fITSchi2MIP[3]<7.5&&fESDTrack.fITSchi2MIP[4]<6.")
  54
  55 AliESDComparisonDraw comp;
  56 comp.SetIO();
  57
  58 //
  59 //example
  60 comp.fTree->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)");
  61
  62 TH1F his("his","his",100,0,20);
  63 TH1F hpools("hpools","hpools",100,-7,7);
  64 TH1F hfake("hfake","hfake",1000,0,150);
  65 TProfile profp0("profp0","profp0",20,-0.4,0.9)
  66
  67 comp.DrawLogXY("fTPCinP0[3]","fTPCDelta[4]/fTPCinP1[3]","fReconstructed==1&&abs(fPdg)==211"+cprim,"1",4,0.2,1.5,-0.06,0.06)
  68 comp.fRes->Draw();
  69 comp.fMean->Draw();
  70 comp.Eff("fTPCinP0[3]","fRowsWithDigits>120&&abs(fPdg)==211"+cteta1+cpos1+cprim,"fTPCOn",20,0.2,1.5)
  71 comp.fRes->Draw();
  72
  73 comp.fTree->Draw("fESDTrack.fITSsignal/fESDTrack.fTPCsignal","fITSOn&&fTPCOn&&fESDTrack.fITSFakeRatio==0")
  74
  75 TH1F his("his","his",100,0,20);
  76 TH1F hpools("hpools","hpools",100,-7,7);
  77
  78 TH2F * hdedx0 = new TH2F("dEdx0","dEdx0",100, 0,2,200,0,250); hdedx0->SetMarkerColor(2);
  79 TH2F * hdedx1 = new TH2F("dEdx1","dEdx1",100, 0,2,200,0,250); hdedx1->SetMarkerColor(3);
  80 TH2F * hdedx2 = new TH2F("dEdx2","dEdx2",100, 0,2,200,0,250); hdedx2->SetMarkerColor(4);
  81 TH2F * hdedx3 = new TH2F("dEdx3","dEdx3",100, 0,2,200,0,250); hdedx3->SetMarkerColor(5);
  82 comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx0","fITSOn&&MC.fParticle.P()<2&&abs(fPdg)==211&&fITStrack.fN==6"+cprim)
  83 comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx1","fITSOn&&MC.fParticle.P()<2&&abs(fPdg)==2212&&fITStrack.fN==6"+cprim)
  84 comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx2","fITSOn&&MC.fParticle.P()<2&&abs(fPdg)==321&&fITStrack.fN==6"+cprim)
  85 comp.fTree->Draw("fESDTrack.fITSsignal:MC.fParticle.P()>>dEdx3","fITSOn&&MC.fParticle.P()<2&&abs(fPdg)==11&&fITStrack.fN==6"+cprim)
  86 hdedx0->Draw(); hdedx1->Draw("same"); hdedx2->Draw("same"); hdedx3->Draw("same");
  87
  88 comp.DrawXY("fITSinP0[3]","fITSPools[4]","fReconstructed==1&&fPdg==-211&&fITSOn"+cprim,"1",4,0.2,1.0,-8,8)
  89
  90 */
  91
  92
  93 #if !defined(__CINT__) || defined(__MAKECINT__)
  94 #include <stdio.h>
  95 #include <string.h>
  96 //ROOT includes
  97 #include "Rtypes.h"
  98 #include "TFile.h"
  99 #include "TTree.h"
 100 #include "TChain.h"
 101 #include "TCut.h"
 102 #include "TString.h"
 103 #include "TBenchmark.h"
 104 #include "TStopwatch.h"
 105 #include "TParticle.h"
 106 #include "TSystem.h"
 107 #include "TTimer.h"
 108 #include "TVector3.h"
 109 #include "TPad.h"
 110 #include "TCanvas.h"
 111 #include "TH1F.h"
 112 #include "TH2F.h"
 113 #include "TF1.h"
 114 #include "TText.h"
 115 #include "Getline.h"
 116 #include "TStyle.h"
 117
 118 //ALIROOT includes
 119 #include "AliRun.h"
 120 #include "AliStack.h"
 121 #include "AliESDtrack.h"
 122 #include "AliSimDigits.h"
 123 #include "AliTPCParam.h"
 124 #include "AliTPC.h"
 125 #include "AliTPCLoader.h"
 126 #include "AliDetector.h"
 127 #include "AliTrackReference.h"
 128 #include "AliRun.h"
 129 #include "AliTPCParamSR.h"
 130 #include "AliTracker.h"
 131 #include "AliComplexCluster.h"
 132 #include "AliMagF.h"
 133 #include "AliESD.h"
 134 #include "AliESDtrack.h"
 135 #include "AliITStrackV2.h"
 136 #include "AliTRDtrack.h"
 137 #endif
 138 #include "AliGenInfo.h"
 139 #include "AliESDComparisonMI.h"
 140
 141
 142
 143 //
 144 //
 145 void AliESDRecInfo::Update(AliMCInfo* info,AliTPCParam * par, Bool_t reconstructed)
 146 {
 147   //
 148   //
 149   //calculates derived variables
 150   //
 151   //
 152   AliTrackReference * ref = &(info->fTrackRef);
 153   fTPCinR0[0] = info->fTrackRef.X();
 154   fTPCinR0[1] = info->fTrackRef.Y();
 155   fTPCinR0[2] = info->fTrackRef.Z();
 156   fTPCinR0[3] = TMath::Sqrt(fTPCinR0[0]*fTPCinR0[0]+fTPCinR0[1]*fTPCinR0[1]);
 157   fTPCinR0[4] = TMath::ATan2(fTPCinR0[1],fTPCinR0[0]);
 158   //
 159   fTPCinP0[0] = ref->Px();
 160   fTPCinP0[1] = ref->Py();
 161   fTPCinP0[2] = ref->Pz();
 162   fTPCinP0[3] = ref->Pt();
 163   fTPCinP0[4] = ref->P();
 164   fDeltaP     = (ref->P()-info->fParticle.P())/info->fParticle.P();
 165   //
 166   //
 167   if (fTPCinP0[3]>0.0000001){
 168     //
 169     fTPCAngle0[0] = TMath::ATan2(fTPCinP0[1],fTPCinP0[0]);
 170     fTPCAngle0[1] = TMath::ATan(fTPCinP0[2]/fTPCinP0[3]);
 171   }
 172   //
 173   //
 174   fITSinP0[0]=info->fParticle.Px();
 175   fITSinP0[1]=info->fParticle.Py();
 176   fITSinP0[2]=info->fParticle.Pz();
 177   fITSinP0[3]=info->fParticle.Pt();
 178   //
 179   fITSinR0[0]=info->fParticle.Vx();
 180   fITSinR0[1]=info->fParticle.Vy();
 181   fITSinR0[2]=info->fParticle.Vz();
 182   fITSinR0[3] = TMath::Sqrt(fITSinR0[0]*fITSinR0[0]+fITSinR0[1]*fITSinR0[1]);
 183   fITSinR0[4] = TMath::ATan2(fITSinR0[1],fITSinR0[0]);
 184   //
 185   //
 186   if (fITSinP0[3]>0.0000001){
 187     fITSAngle0[0] = TMath::ATan2(fITSinP0[1],fITSinP0[0]);
 188     fITSAngle0[1] = TMath::ATan(fITSinP0[2]/fITSinP0[3]);
 189   }
 190   //
 191   if (reconstructed==kFALSE){
 192     fReconstructed = kFALSE;
 193     fTPCOn = kFALSE;
 194     fITSOn = kFALSE;
 195     fTRDOn = kFALSE;
 196     return;
 197   }
 198   fReconstructed = kTRUE;
 199   fTPCOn = ((fESDTrack.GetStatus()&AliESDtrack::kTPCrefit)>0) ? kTRUE : kFALSE;
 200   fITSOn = ((fESDTrack.GetStatus()&AliESDtrack::kITSrefit)>0) ? kTRUE : kFALSE;
 201   fTRDOn = ((fESDTrack.GetStatus()&AliESDtrack::kTRDrefit)>0) ? kTRUE : kFALSE;
 202
 203   if (fTPCOn){
 204     //TPC
 205     fESDTrack.GetInnerXYZ(fTPCinR1);
 206     fTPCinR1[3] = TMath::Sqrt(fTPCinR1[0]*fTPCinR1[0]+fTPCinR1[1]*fTPCinR1[1]);
 207     fTPCinR1[4] = TMath::ATan2(fTPCinR1[1],fTPCinR1[0]);
 208     fESDTrack.GetInnerPxPyPz(fTPCinP1);
 209     fTPCinP1[3] = TMath::Sqrt(fTPCinP1[0]*fTPCinP1[0]+fTPCinP1[1]*fTPCinP1[1]);
 210     fTPCinP1[4] = TMath::Sqrt(fTPCinP1[3]*fTPCinP1[3]+fTPCinP1[2]*fTPCinP1[2]);
 211     //
 212     //
 213     if (fTPCinP1[3]>0.0000001){
 214       fTPCAngle1[0] = TMath::ATan2(fTPCinP1[1],fTPCinP1[0]);
 215       fTPCAngle1[1] = TMath::ATan(fTPCinP1[2]/fTPCinP1[3]);
 216     }
 217     Double_t cov[15], param[5],x;
 218     fESDTrack.GetInnerExternalCovariance(cov);
 219     fESDTrack.GetInnerExternalParameters(x,param);
 220     //
 221     fTPCDelta[0] = (fTPCinR0[4]-fTPCinR1[4])*fTPCinR1[3];  //delta rfi
 222     fTPCPools[0] = fTPCDelta[0]/TMath::Sqrt(cov[0]);
 223     fTPCDelta[1] = (fTPCinR0[2]-fTPCinR1[2]);              //delta z
 224     fTPCPools[1] = fTPCDelta[1]/TMath::Sqrt(cov[2]);
 225     fTPCDelta[2] = (fTPCAngle0[0]-fTPCAngle1[0]);
 226     fTPCPools[2] = fTPCDelta[2]/TMath::Sqrt(cov[5]);
 227     fTPCDelta[3] = (TMath::Tan(fTPCAngle0[1])-TMath::Tan(fTPCAngle1[1]));
 228     fTPCPools[3] = fTPCDelta[3]/TMath::Sqrt(cov[9]);
 229     fTPCDelta[4] = (fTPCinP0[3]-fTPCinP1[3]);
 230     Double_t sign = (param[4]>0)? 1.:-1;
 231     fTPCPools[4] = sign*(1./fTPCinP0[3]-1./fTPCinP1[3])/TMath::Sqrt(cov[14]);
 232   }
 233   if (fITSOn){
 234     // ITS
 235     Double_t param[5],x;
 236     //fESDTrack.GetExternalParameters(x,param);
 237     fESDTrack.GetConstrainedExternalParameters(x,param);
 238     Double_t cov[15];
 239     fESDTrack.GetExternalCovariance(cov);
 240     fESDTrack.GetConstrainedExternalCovariance(cov);
 241     if (TMath::Abs(param[4])<0.0000000001) return;
 242
 243     fESDTrack.GetXYZ(fITSinR1);
 244     fESDTrack.GetPxPyPz(fITSinP1);
 245     fITSinP1[3] = TMath::Sqrt(fITSinP1[0]*fITSinP1[0]+fITSinP1[1]*fITSinP1[1]);
 246     //
 247     fITSinR1[3] = TMath::Sqrt(fITSinR1[0]*fITSinR1[0]+fITSinR1[1]*fITSinR1[1]);
 248     fITSinR1[4] = TMath::ATan2(fITSinR1[1],fITSinR1[0]);
 249     //
 250     //
 251     if (fITSinP1[3]>0.0000001){
 252       fITSAngle1[0] = TMath::ATan2(fITSinP1[1],fITSinP1[0]);
 253       fITSAngle1[1] = TMath::ATan(fITSinP1[2]/fITSinP1[3]);
 254     }
 255     //
 256     //
 257     fITSDelta[0] = (fITSinR0[4]-fITSinR1[4])*fITSinR1[3];  //delta rfi
 258     fITSPools[0] = fITSDelta[0]/TMath::Sqrt(cov[0]);
 259     fITSDelta[1] = (fITSinR0[2]-fITSinR1[2]);              //delta z
 260     fITSPools[1] = fITSDelta[1]/TMath::Sqrt(cov[2]);
 261     fITSDelta[2] = (fITSAngle0[0]-fITSAngle1[0]);
 262     fITSPools[2] = fITSDelta[2]/TMath::Sqrt(cov[5]);
 263     fITSDelta[3] = (TMath::Tan(fITSAngle0[1])-TMath::Tan(fITSAngle1[1]));
 264     fITSPools[3] = fITSDelta[3]/TMath::Sqrt(cov[9]);
 265     fITSDelta[4] = (fITSinP0[3]-fITSinP1[3]);
 266     Double_t sign = (param[4]>0) ? 1:-1;
 267     fITSPools[4] = sign*(1./fITSinP0[3]-1./fITSinP1[3])/TMath::Sqrt(cov[14]);
 268   }
 269
 270 }
 271
 272
 273 void  AliESDRecV0Info::Update(Float_t vertex[3])
 274 {
 275   Float_t v[3] = {fXr[0]-vertex[0],fXr[1]-vertex[1],fXr[2]-vertex[2]};
 276   Float_t p[3] = {fPdr[0]+fPm[0], fPdr[1]+fPm[1],fPdr[2]+fPm[2]};
 277
 278   Float_t vnorm2 = v[0]*v[0]+v[1]*v[1];
 279   Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
 280   vnorm2 = TMath::Sqrt(vnorm2);
 281   Float_t pnorm2 = p[0]*p[0]+p[1]*p[1];
 282   Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
 283   pnorm2 = TMath::Sqrt(pnorm2);
 284
 285   fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
 286   fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);
 287   fPointAngle   = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
 288
 289 }
 290
 291
 292 ////////////////////////////////////////////////////////////////////////
 293 ESDCmpTr::ESDCmpTr()
 294 {
 295   Reset();
 296 }
 297
 298 ////////////////////////////////////////////////////////////////////////
 299 ESDCmpTr::ESDCmpTr(const char* fnGenTracks,
 300                    const char* fnCmp,
 301                    const char* fnGalice, Int_t direction,
 302                    Int_t nEvents, Int_t firstEvent)
 303 {
 304   Reset();
 305   //  fFnGenTracks = fnGenTracks;
 306   //  fFnCmp = fnCmp;
 307   sprintf(fFnGenTracks,"%s",fnGenTracks);
 308   sprintf(fFnCmp,"%s",fnCmp);
 309
 310   fFirstEventNr = firstEvent;
 311   fEventNr = firstEvent;
 312   fNEvents = nEvents;
 313   fDirection = direction;
 314   //
 315   fLoader = AliRunLoader::Open(fnGalice);
 316   if (gAlice){
 317     //delete gAlice->GetRunLoader();
 318     delete gAlice;
 319     gAlice = 0x0;
 320   }
 321   if (fLoader->LoadgAlice()){
 322     cerr<<"Error occured while l"<<endl;
 323   }
 324   Int_t nall = fLoader->GetNumberOfEvents();
 325   if (nEvents==0) {
 326     nEvents =nall;
 327     fNEvents=nall;
 328     fFirstEventNr=0;
 329   }
 330
 331   if (nall<=0){
 332     cerr<<"no events available"<<endl;
 333     fEventNr = 0;
 334     return;
 335   }
 336   if (firstEvent+nEvents>nall) {
 337     fEventNr = nall-firstEvent;
 338     cerr<<"restricted number of events availaible"<<endl;
 339   }
 340   AliMagF * magf = gAlice->Field();
 341   AliTracker::SetFieldMap(magf);
 342
 343 }
 344
 345
 346 ////////////////////////////////////////////////////////////////////////
 347 ESDCmpTr::~ESDCmpTr()
 348 {
 349   if (fLoader) {
 350     delete fLoader;
 351   }
 352 }
 353
 354 //////////////////////////////////////////////////////////////
 355 Int_t ESDCmpTr::SetIO()
 356 {
 357   //
 358   //
 359   CreateTreeCmp();
 360   if (!fTreeCmp) return 1;
 361   fParamTPC = GetTPCParam();
 362   //
 363   if (!ConnectGenTree()) {
 364     cerr<<"Cannot connect tree with generated tracks"<<endl;
 365     return 1;
 366   }
 367   return 0;
 368 }
 369
 370 //////////////////////////////////////////////////////////////
 371
 372 Int_t ESDCmpTr::SetIO(Int_t eventNr)
 373 {
 374   //
 375   //
 376   // SET INPUT
 377   //
 378   TFile f("AliESDs.root");
 379   //
 380
 381   TTree* tree = (TTree*) f.Get("esdTree");
 382   if (!tree) {
 383     Char_t ename[100];
 384     sprintf(ename,"%d",eventNr);
 385     fEvent = (AliESD*)f.Get(ename);
 386     if (!fEvent){
 387       sprintf(ename,"ESD%d",eventNr);
 388       fEvent = (AliESD*)f.Get(ename);
 389     }
 390   }
 391   else{
 392     tree->SetBranchAddress("ESD", &fEvent);
 393     tree->GetEntry(eventNr);
 394   }
 395
 396
 397   /*
 398   Char_t ename[100];
 399   sprintf(ename,"%d",eventNr);
 400   fEvent = (AliESD*)f.Get(ename);
 401   if (!fEvent){
 402     sprintf(ename,"ESD%d",eventNr);
 403     fEvent = (AliESD*)f.Get(ename);
 404   }
 405
 406   TTree* tree = (TTree*) f.Get("esdTree");
 407   if (!tree) {
 408     Error("CheckESD", "no ESD tree found");
 409     return kFALSE;
 410   }
 411   tree->SetBranchAddress("ESD", &fEvent);
 412   tree->GetEntry(eventNr);
 413   */
 414
 415   if (!fEvent) return 1;
 416
 417   return 0;
 418 }
 419
 420
 421
 422 ////////////////////////////////////////////////////////////////////////
 423 void ESDCmpTr::Reset()
 424 {
 425   fEventNr = 0;
 426   fNEvents = 0;
 427   fTreeCmp = 0;
 428   //  fFnCmp = "cmpTracks.root";
 429   fFileGenTracks = 0;
 430   fDebug = 0;
 431   //
 432   fParamTPC = 0;
 433   fEvent =0;
 434 }
 435
 436 ////////////////////////////////////////////////////////////////////////
 437 Int_t ESDCmpTr::Exec(Int_t nEvents, Int_t firstEventNr)
 438 {
 439   fNEvents = nEvents;
 440   fFirstEventNr = firstEventNr;
 441   return Exec();
 442 }
 443
 444 ////////////////////////////////////////////////////////////////////////
 445 Int_t ESDCmpTr::Exec()
 446 {
 447   TStopwatch timer;
 448   timer.Start();
 449
 450   if (SetIO()==1)
 451     return 1;
 452
 453   fNextTreeGenEntryToRead = 0;
 454   cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl;
 455   for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents;
 456        eventNr++) {
 457     fEventNr = eventNr;
 458     SetIO(fEventNr);
 459     fNParticles = gAlice->GetEvent(fEventNr);
 460
 461     fIndexRecTracks = new Int_t[fNParticles*4];  //write at maximum 4 tracks corresponding to particle
 462     fFakeRecTracks = new Int_t[fNParticles];
 463     fMultiRecTracks = new Int_t[fNParticles];
 464     for (Int_t i = 0; i<fNParticles; i++) {
 465       fIndexRecTracks[4*i] = -1;
 466       fIndexRecTracks[4*i+1] = -1;
 467       fIndexRecTracks[4*i+2] = -1;
 468       fIndexRecTracks[4*i+3] = -1;
 469
 470       fFakeRecTracks[i] = 0;
 471       fMultiRecTracks[i] = 0;
 472     }
 473
 474     cout<<"Start to process event "<<fEventNr<<endl;
 475     cout<<"\tfNParticles = "<<fNParticles<<endl;
 476     if (fDebug>2) cout<<"\tStart loop over TreeT"<<endl;
 477     if (TreeTLoop()>0) return 1;
 478
 479     if (fDebug>2) cout<<"\tStart loop over tree genTracks"<<endl;
 480     if (TreeGenLoop(eventNr)>0) return 1;
 481     if (fDebug>2) cout<<"\tEnd loop over tree genTracks"<<endl;
 482
 483     delete [] fIndexRecTracks;
 484     delete [] fFakeRecTracks;
 485     delete [] fMultiRecTracks;
 486   }
 487
 488   CloseOutputFile();
 489
 490   cerr<<"Exec finished"<<endl;
 491   timer.Stop();
 492   timer.Print();
 493   return 0;
 494
 495 }
 496 ////////////////////////////////////////////////////////////////////////
 497 Bool_t ESDCmpTr::ConnectGenTree()
 498 {
 499 //
 500 // connect all branches from the genTracksTree
 501 // use the same variables as for the new cmp tree, it may work
 502 //
 503   fFileGenTracks = TFile::Open(fFnGenTracks,"READ");
 504   if (!fFileGenTracks) {
 505     cerr<<"Error in ConnectGenTree: cannot open file "<<fFnGenTracks<<endl;
 506     return kFALSE;
 507   }
 508   fTreeGenTracks = (TTree*)fFileGenTracks->Get("genTracksTree");
 509   if (!fTreeGenTracks) {
 510     cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
 511         <<fFnGenTracks<<endl;
 512     return kFALSE;
 513   }
 514   //
 515   fMCInfo = new AliMCInfo;
 516   fTreeGenTracks->SetBranchAddress("MC",&fMCInfo);
 517   //
 518   if (fDebug > 1) {
 519     cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl;
 520   }
 521   return kTRUE;
 522 }
 523
 524
 525 ////////////////////////////////////////////////////////////////////////
 526 void ESDCmpTr::CreateTreeCmp()
 527 {
 528   fFileCmp = TFile::Open(fFnCmp,"RECREATE");
 529   if (!fFileCmp) {
 530     cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl;
 531     return;
 532   }
 533
 534
 535   fTreeCmp    = new TTree("ESDcmpTracks","ESDcmpTracks");
 536   //
 537   fMCInfo = new AliMCInfo;
 538   fRecInfo = new AliESDRecInfo;
 539   //
 540   AliESDtrack * esdTrack = new AliESDtrack;
 541   AliITStrackV2 * itsTrack = new AliITStrackV2;
 542   //
 543   fTreeCmp->Branch("MC","AliMCInfo",&fMCInfo);
 544   fTreeCmp->Branch("RC","AliESDRecInfo",&fRecInfo);
 545   fTreeCmp->Branch("fESDTrack","AliESDtrack",&esdTrack);
 546   //  fTreeCmp->Branch("ITS","AliITStrackV2",&itsTrack);
 547   delete esdTrack;
 548   //
 549   fTreeCmp->AutoSave();
 550
 551 }
 552 ////////////////////////////////////////////////////////////////////////
 553 void ESDCmpTr::CloseOutputFile()
 554 {
 555   if (!fFileCmp) {
 556     cerr<<"File "<<fFnCmp<<" not found as an open file."<<endl;
 557     return;
 558   }
 559   fFileCmp->cd();
 560   fTreeCmp->Write();
 561   delete fTreeCmp;
 562
 563   fFileCmp->Close();
 564   delete fFileCmp;
 565   return;
 566 }
 567 ////////////////////////////////////////////////////////////////////////
 568
 569 TVector3 ESDCmpTr::TR2Local(AliTrackReference *trackRef,
 570                             AliTPCParam *paramTPC) {
 571
 572   Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
 573   Int_t index[4];
 574   paramTPC->Transform0to1(x,index);
 575   paramTPC->Transform1to2(x,index);
 576   return TVector3(x);
 577 }
 578 ////////////////////////////////////////////////////////////////////////
 579
 580 Int_t ESDCmpTr::TreeTLoop()
 581 {
 582   //
 583   // loop over all ESD reconstructed tracks and store info in memory
 584   //
 585   TStopwatch  timer;
 586   timer.Start();
 587   //
 588   Int_t nEntries = (Int_t)fEvent->GetNumberOfTracks();
 589   //
 590   fTracks      = new TObjArray(nEntries);
 591   //
 592   //load tracks to the memory
 593   for (Int_t i=0; i<nEntries;i++){
 594     AliESDtrack * track =fEvent->GetTrack(i);
 595     fTracks->AddAt(track,i);
 596   }
 597   //
 598   //
 599   AliESDtrack * track=0;
 600   for (Int_t iEntry=0; iEntry<nEntries;iEntry++){
 601     track = (AliESDtrack*)fTracks->UncheckedAt(iEntry);
 602     //
 603     Int_t label = track->GetLabel();
 604     Int_t absLabel = abs(label);
 605     if (absLabel < fNParticles) {
 606       //      fIndexRecTracks[absLabel] =  iEntry;
 607       if (label < 0) fFakeRecTracks[absLabel]++;
 608       if (fMultiRecTracks[absLabel]>0){
 609         if (fMultiRecTracks[absLabel]<4)
 610           fIndexRecTracks[absLabel*4+fMultiRecTracks[absLabel]] =  iEntry;
 611       }
 612       else
 613         fIndexRecTracks[absLabel*4] =  iEntry;
 614       fMultiRecTracks[absLabel]++;
 615     }
 616   }
 617   printf("Time spended in TreeTLoop\n");
 618   timer.Print();
 619
 620   if (fDebug > 2) cerr<<"end of TreeTLoop"<<endl;
 621   return 0;
 622 }
 623
 624 ////////////////////////////////////////////////////////////////////////
 625 Int_t ESDCmpTr::TreeGenLoop(Int_t eventNr)
 626 {
 627 //
 628 // loop over all entries for a given event, find corresponding
 629 // rec. track and store in the fTreeCmp
 630 //
 631   TStopwatch timer;
 632   timer.Start();
 633   Int_t entry = fNextTreeGenEntryToRead;
 634   Double_t nParticlesTR = fTreeGenTracks->GetEntriesFast();
 635   cerr<<"fNParticles, nParticlesTR, fNextTreeGenEntryToRead: "<<fNParticles<<" "
 636       <<nParticlesTR<<" "<<fNextTreeGenEntryToRead<<endl;
 637   TBranch * branch = fTreeCmp->GetBranch("RC");
 638   branch->SetAddress(&fRecInfo); // set all pointers
 639
 640   while (entry < nParticlesTR) {
 641     fTreeGenTracks->GetEntry(entry);
 642     entry++;
 643     if (eventNr < fMCInfo->fEventNr) continue;
 644     if (eventNr > fMCInfo->fEventNr) continue;;
 645     //
 646     fNextTreeGenEntryToRead = entry-1;
 647     if (fDebug > 2 && fMCInfo->fLabel < 10) {
 648       cerr<<"Fill track with a label "<<fMCInfo->fLabel<<endl;
 649     }
 650     //
 651     fRecInfo->Reset();
 652     AliESDtrack * track=0;
 653     fRecInfo->fReconstructed =0;
 654     TVector3 local = TR2Local(&(fMCInfo->fTrackRef),fParamTPC);
 655     local.GetXYZ(fRecInfo->fTRLocalCoord);
 656     //
 657     if (fIndexRecTracks[fMCInfo->fLabel*4] >= 0) {
 658       track= (AliESDtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4]);
 659       //
 660       //
 661       // find nearest track if multifound
 662       if (fIndexRecTracks[fMCInfo->fLabel*4]+1){
 663         Double_t xyz[3];
 664         track->GetInnerXYZ(xyz);
 665         Float_t dz = TMath::Abs(local.Z()-xyz[2]);
 666         for (Int_t i=1;i<4;i++){
 667           if (fIndexRecTracks[fMCInfo->fLabel*4+i]>=0){
 668             AliESDtrack * track2 = (AliESDtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4+i]);
 669             track2->GetInnerXYZ(xyz);
 670             Float_t dz2 = TMath::Abs(local.Z()-xyz[2]);
 671             if  (TMath::Abs(dz2)<dz)
 672               track = track2;
 673           }
 674         }
 675       }
 676       //
 677       fRecInfo->fESDTrack =*track;
 678       if (track->GetITStrack())
 679         fRecInfo->fITStrack = *((AliITStrackV2*)track->GetITStrack());
 680       if (track->GetTRDtrack()){
 681         fRecInfo->fTRDtrack = *((AliTRDtrack*)track->GetTRDtrack());
 682       }
 683       else{
 684         fRecInfo->fTRDtrack.SetdEdx(-1);
 685       }
 686       fRecInfo->fReconstructed = 1;
 687       fRecInfo->fFake     = fFakeRecTracks[fMCInfo->fLabel];
 688       fRecInfo->fMultiple = fMultiRecTracks[fMCInfo->fLabel];
 689       //
 690       fRecInfo->Update(fMCInfo,fParamTPC,kTRUE);
 691     }
 692     else{
 693       fRecInfo->Update(fMCInfo,fParamTPC,kFALSE);
 694     }
 695
 696     fTreeCmp->Fill();
 697   }
 698   fTreeCmp->AutoSave();
 699   fTracks->Delete();
 700   printf("Time spended in TreeGenLoop\n");
 701   timer.Print();
 702   if (fDebug > 2) cerr<<"end of TreeGenLoop"<<endl;
 703
 704   return 0;
 705 }
 706 ////////////////////////////////////////////////////////////////////////
 707 ////////////////////////////////////////////////////////////////////////
 708
 709 void AliESDComparisonDraw::SetIO(const char *fname)
 710 {
 711   //
 712    TFile* file = TFile::Open(fname);
 713    //
 714    fTree = (TTree*) file->Get("ESDcmpTracks");
 715    if (!fTree) {
 716     printf("no track comparison tree found\n");
 717     file->Close();
 718     delete file;
 719   }
 720 }
 721
 722