HLT/TPCLib/macros/AliHLTTPCTrackerEvaluation.C

   1 // $Id$
   2 /*
   3  * Modified version of $ALICE_ROOT/TPC/AliTPCComparison.C for evaluating the
   4  * performance of the CA tracker.
   5  *
   6  * Usage:
   7  * <pre>
   8  *   aliroot AliHLTTPCTrackerEvaluation.C
   9  * </pre>
  10  *
  11  * documentation to be filled in when the macro is finalized
  12  *
  13  *
  14  * @ingroup alihlt_tpc
  15  */
  16 #if defined(__CINT__) && !defined(__MAKECINT__)
  17
  18 //Int_t AliHLTTPCTrackerEvaluation(const Char_t *dir=".", const char* input="AliHLTESDs.root", Float_t ptLowerCut=0.1, Float_t ptUpperCut=10., Bool_t bDedxAndClus=kFALSE)
  19 {
  20   gSystem->AddIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/include -I$ALICE_ROOT -I$ALICE_ROOT/TPC -I$ALICE/geant3/TGeant3");
  21   gROOT->LoadMacro("AliHLTTPCTrackerEvaluation.C++");
  22
  23   AliHLTTPCTrackerEvaluation();
  24   // the attempt to pass the arguments with a self-called macro hasn't worked so far
  25   //AliHLTTPCTrackerEvaluation(dir, input, ptLowerCut, ptUpperCut, bDedxAndClus);
  26
  27 }
  28 #else
  29
  30 #include <AliLog.h>
  31 #include <TMath.h>
  32 #include <TError.h>
  33 #include <Riostream.h>
  34 #include <TH1F.h>
  35 #include <TH2F.h>
  36 #include <TF1.h>
  37 #include <TTree.h>
  38 #include <TParticle.h>
  39 #include <TCanvas.h>
  40 #include <TLegend.h>
  41 #include <TLine.h>
  42 #include <TText.h>
  43 #include <TBenchmark.h>
  44 #include <TStyle.h>
  45 #include <TFile.h>
  46 #include <TROOT.h>
  47
  48 #include "AliStack.h"
  49 #include "AliHeader.h"
  50 #include "AliTrackReference.h"
  51 #include "AliRunLoader.h"
  52 #include "AliRun.h"
  53 #include "AliESDEvent.h"
  54 #include "AliESDtrack.h"
  55 #include "AliTPCtrack.h"
  56
  57 #include "AliSimDigits.h"
  58 #include "AliTPC.h"
  59 #include "AliTPCParamSR.h"
  60 #include "AliTPCClustersArray.h"
  61 #include "AliTPCClustersRow.h"
  62 #include "AliTPCcluster.h"
  63 #include "AliTPCLoader.h"
  64 #include "TParticlePDG.h"
  65 #include "TDatabasePDG.h"
  66
  67 #include "AliCDBManager.h"
  68 #include "AliTPCcalibDB.h"
  69 #include "TProfile.h"
  70
  71 //_______________________________________________________________________________________
  72
  73
  74 Int_t GoodTPCTracks(const Char_t *dir = ".");
  75
  76 extern AliRun     *gAlice;
  77 extern TBenchmark *gBenchmark;
  78 extern TROOT      *gROOT;
  79
  80 static Int_t allgood     = 0;
  81 static Int_t allselected = 0;
  82 static Int_t allfound    = 0;
  83
  84
  85
  86 Int_t AliHLTTPCTrackerEvaluation(const Char_t *dir=".", const char* input="AliHLTESDs.root", Float_t ptLowerCut=0.1, Float_t ptUpperCut=10., Bool_t bDedxAndClus=kFALSE){
  87
  88    if(!input){
  89        cerr << "please specify an input file" << endl;
  90        return 1;
  91    }
  92
  93    gBenchmark->Start("AliHLTTPCTrackerEvaluation");
  94
  95    ::Info("AliHLTTPCTrackerEvaluation.C","Calculating reconstruction efficiency...");
  96
  97    gSystem->SetIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/include -I$ALICE_ROOT -I$ALICE/geant3/TGeant3");
  98
  99    gStyle->SetTitleOffset(1.2,"X");
 100    gStyle->SetTitleOffset(1.1,"Y");
 101
 102
 103    // --------------- Efficiency histograms ------------------------
 104
 105    TH1F *hClus      = new TH1F("hClus","# clusters on reconstructed tracks",300,0,300);
 106
 107    TH1F *hGood      = new TH1F("hGood",     "P_{t} distribution of MC tracks (selected for efficiency calculation)", 30, 0, 6);
 108    TH1F *hFake      = new TH1F("hFake",     "P_{t} distribution of fake tracks",                                     30, 0, 6);
 109    TH1F *hFound     = new TH1F("hFound",    "P_{t} distribution of reconstructed tracks",                            30, 0, 6);
 110    TH1F *hFoundMult = new TH1F("hFoundMult","P_{t} distribution of reconstructed tracks with clones",                30, 0, 6);
 111    TH1F *hClone     = new TH1F("hClone",    "P_{t} distribution of clone tracks",                                    30, 0, 6);
 112
 113    TH1F *hEff      = new TH1F("hEff",     "Reconstruction efficiency based on selected MC tracks", 30, 0, 6);
 114    TH1F *hFakeEff  = new TH1F("hFakeEff", "Efficiency for fake tracks",                            30, 0, 6);
 115    TH1F *hCloneEff = new TH1F("hCloneEff","Efficiency for clone tracks",                           30, 0, 6);
 116
 117    hEff->SetLineWidth(2);
 118    hEff->SetMaximum(1.4);
 119    hEff->SetYTitle("efficiency");
 120    hEff->GetYaxis()->CenterTitle();
 121    hEff->SetXTitle("P_{t} (GeV/c)");
 122
 123    hFakeEff->SetLineColor(kRed);
 124    hFakeEff->SetLineWidth(2);
 125
 126    hCloneEff->SetLineColor(kBlue);
 127    hCloneEff->SetLineWidth(2);
 128
 129    TH2F *hDedx = new TH2F("hDedx","dE/dx vs. momentum", 50,0.,2.,50,0.,400.);
 130    hDedx->SetMarkerStyle(8);
 131    hDedx->SetMarkerSize(0.4);
 132    hDedx->SetXTitle("P (Gev/c)");
 133    hDedx->SetYTitle("dE/dx (a.u.)");
 134
 135    TH1F *hnhit_ref = new TH1F("hnhit_ref","# clusters on reco tracks, used in fit performance",300,0,300);
 136    TH1F *he = new TH1F("he","dE/dx for pions with 0.4<p<0.5 GeV/c",50,0.,100.);
 137
 138    // --------------- Resolution histograms ---------------------------
 139
 140    TH2F *hPt     = new TH2F("hPt",     "", 30, 0, 6, 50,  -3,  3);
 141    TH2F *hPhi    = new TH2F("hPhi",    "", 30, 0, 6, 50, -16, 16);
 142    TH2F *hLambda = new TH2F("hLambda", "", 30, 0, 6, 50, -10, 10);
 143    TH2F *hY      = new TH2F("hY",      "", 30, 0, 6, 30,  -5,  5);
 144    TH2F *hZ      = new TH2F("hZ",      "", 30, 0, 250, 30,  -5,  5); // fix the x axis, if this is z dependent
 145
 146    TProfile *hResPt     = new TProfile("hResPt",     "P_{t} resolution (%)",  30, 0, 6);
 147    TProfile *hResPhi    = new TProfile("hResPhi",    "#phi resolution (mrad)",   30, 0, 6);
 148    TProfile *hResLambda = new TProfile("hResLambda", "#lambda resolution (mrad)",30, 0, 6);
 149    TProfile *hResY      = new TProfile("hResY",      "Y resolution (mm)",      30, 0, 6);
 150    TProfile *hResZ      = new TProfile("hResZ",      "Z resolution (mm)",      30, 0, 250); // select appropriate x axis, KKK
 151
 152    hResPt->SetXTitle("P_{t} (GeV/c)");
 153    hResPt->SetYTitle("#sigma_{(P_{t}-P_{t_{MC}})/P_{t_{MC}}} (%)");
 154    hResPt->GetYaxis()->CenterTitle(true);
 155    hResPt->GetYaxis()->CenterTitle(true);
 156
 157    hResPhi->SetXTitle("P_{t} (GeV/c)");
 158    hResPhi->SetYTitle("#sigma_{(#phi_{rec}-#phi_{MC})} (mrad)");
 159    hResPhi->GetYaxis()->CenterTitle(true);
 160
 161    hResLambda->SetXTitle("P_{t} (GeV/c)");
 162    hResLambda->SetYTitle("#sigma_{(#lambda_{rec}-#lambda_{MC})} (mrad)"); // KKK perhaps we should add the definition on the histo pad
 163    hResLambda->GetYaxis()->CenterTitle(true);
 164
 165    hResY->SetXTitle("P_{t} (GeV/c)");
 166    hResY->SetYTitle("#sigma_{(Y_{rec}-Y_{MC})} (mm)");
 167    hResY->GetYaxis()->CenterTitle(true);
 168
 169    hResZ->SetXTitle("Z (mm)");
 170    hResZ->SetYTitle("#sigma_{(Z_{rec}-Z_{MC})} (mm)");
 171    hResZ->GetYaxis()->CenterTitle(true);
 172
 173    // KKK I leave the Z coordinate up to you. You can make the hResZ Z dependent.
 174
 175    /*
 176
 177    TH1D *hProjPt[15], *hProjPhi[15], *hProjLambda[15], *hProjY[15], *hProjZ[15]; // I am not sure we need these histos to be pointers
 178    Char_t name[15];  Char_t title[100];
 179
 180    for(Int_t i=0; i<15; i++){
 181        sprintf(name,"hProjPt%i",i+1);
 182        sprintf(title,"(Pt_{MC}-Pt_{Rec})/Pt_{MC} @ Pt#in[%f, %f] GeV/c", i-0.5, i+0.5);
 183        hProjPt[i] = new TH1D(name, title, 50, -3, -3);
 184
 185        sprintf(name,"hProjPhi%i",i+1);
 186        sprintf(title,"(#phi_{MC}-#phi_{Rec}) @ #phi#in[%f, %f] GeV/c", i-0.5, i+0.5);
 187        hProjPhi[i] = new TH1D(name, title, 50, -16, -16);
 188
 189        sprintf(name,"hProjLambda%i",i+1);
 190        sprintf(title,"(#lambda_{MC}-#lambda_{Rec}) @ #lambda#in[%f, %f] GeV/c", i-0.5, i+0.5);
 191        hProjLambda[i] = new TH1D(name, title, 50, -16, -16);
 192
 193        sprintf(name,"hProjY%i",i+1);
 194        sprintf(title,"(Y_{MC}-Y_{Rec}) @ Y#in[%f, %f] GeV/c", i-0.5, i+0.5);
 195        hProjY[i] = new TH1D(name, title, 50, -16, -16);
 196
 197        // here you add the hProjZ[] histograms    KKK
 198    }
 199    */
 200
 201
 202    // --------------- Pull variable histograms ------------------------
 203
 204    TH1F *hpullPhi  = new TH1F("hpullPhi", "SinPhi pull", 30,-10.,10.);
 205    TH1F *hpullY    = new TH1F("hpullY",   "Y pull",      30,-10.,10.);
 206    TH1F *hpullZ    = new TH1F("hpullZ",   "Z pull",      30,-10.,10.);
 207    TH1F *hpullDzds = new TH1F("hpullDzds","Dzds pull",   30,-10.,10.);
 208    TH1F *hpullK    = new TH1F("hpullK",   "Kappa pull",  30,-10.,10.);
 209
 210 //---------------------------------------------------------------------------------------
 211
 212    Char_t fname[100];
 213    sprintf(fname,"%s/GoodTPCTracks.root",dir);
 214
 215    TFile *refFile = TFile::Open(fname,"old");
 216    if(!refFile || !refFile->IsOpen()){
 217       ::Info("AliHLTTPCTrackerEvaluation.C","Marking good tracks (will take a while)...");
 218       if(GoodTPCTracks(dir)){
 219          ::Error("AliHLTTPCTrackerEvaluation.C","Cannot generate the reference file!");
 220          return 1;
 221      }
 222    }
 223
 224    refFile = TFile::Open(fname,"old");
 225    if(!refFile || !refFile->IsOpen()){
 226       ::Error("AliHLTTPCTrackerEvaluation.C","Cannot open the reference file!");
 227       return 2;
 228    }
 229
 230    //------------ access the contents of the tree created by the function GoodTracks() -------------
 231
 232    TTree *tpcTree = (TTree*)refFile->Get("tpcTree");
 233    if(!tpcTree){
 234       ::Error("AliHLTTPCTrackerEvaluation.C","Cannot get the reference tree!");
 235       return 3;
 236    }
 237
 238    TBranch *branch = tpcTree->GetBranch("TPC");
 239    if(!branch){
 240       ::Error("AliHLTTPCTrackerEvaluation.C","Cannot get the TPC branch!");
 241       return 4;
 242    }
 243
 244    TClonesArray dummy("AliTrackReference",1000), *refs = &dummy;
 245    branch->SetAddress(&refs);
 246
 247    //------------ access the HLT output --------------------
 248
 249    sprintf(fname,"%s/%s", dir, input);
 250    TFile *ef = TFile::Open(fname);
 251
 252    if( !ef || !ef->IsOpen() ){
 253       ::Error("AliHLTTPCTrackerEvaluation.C","Cannot open [%s]!",fname);
 254       return 5;
 255    }
 256
 257    AliESDEvent *event = new AliESDEvent();
 258    TTree *esdTree = NULL;
 259    TString tsinput = input;
 260
 261    if(tsinput.CompareTo("AliESDs.root")==0){
 262       esdTree = (TTree*)ef->Get("HLTesdTree");
 263       if(!esdTree){
 264          ::Error("AliHLTTPCTrackerEvaluation.C", "no HLTESD tree found");
 265          return 6;
 266       }
 267    event->ReadFromTree(esdTree);
 268    } else if(tsinput.CompareTo("AliHLTESDs.root")==0){
 269       esdTree = (TTree*)ef->Get("esdTree");
 270       if(!esdTree){
 271          ::Error("AliHLTTPCTrackerEvaluation.C", "no ESD tree found");
 272          return 7;
 273       }
 274    event->ReadFromTree(esdTree);
 275    } else return 8;
 276
 277
 278    //---------- Loop over events -------------------
 279
 280    Int_t iEvent = 0;
 281    while(esdTree->GetEvent(iEvent)){
 282
 283       Int_t nEntries = event->GetNumberOfTracks();
 284       cout << "********* Processing event number: " << iEvent << ", " << nEntries << " reconstructed tracks *******\n" << endl;
 285
 286       allfound += nEntries;
 287
 288       if(tpcTree->GetEvent(iEvent++)==0){
 289          cerr << "No reconstructable tracks !\n";
 290          continue;
 291       }
 292
 293       Int_t ngood = refs->GetEntriesFast(); // access the track references
 294       cout << ngood << " good MC tracks" << endl;
 295       allgood += ngood;
 296
 297       const Int_t MAX = 15000;
 298
 299       Int_t track_notfound[MAX],                            itrack_notfound   = 0;
 300       Int_t track_fake[MAX],                                itrack_fake       = 0;
 301       Int_t track_multifound[MAX], track_multifound_n[MAX], itrack_multifound = 0;
 302
 303       for(Int_t i=0; i<nEntries; i++){
 304
 305           hClus->Fill(event->GetTrack(i)->GetTPCNcls()); // filled but not displayed
 306           cout << "TPC label for track " << i << " = " << event->GetTrack(i)->GetTPCLabel() << ", # clusters " << event->GetTrack(i)->GetTPCNcls() << endl;
 307       }
 308
 309       Int_t i;
 310       for(Int_t k=0; k<ngood; k++){ // read the MC information
 311
 312           AliTrackReference *ref = (AliTrackReference*)refs->UncheckedAt(k);
 313
 314           Int_t  label = ref->Label();
 315           Int_t tlabel = -1;
 316           Float_t ptMC = TMath::Sqrt(ref->Px()*ref->Px() + ref->Py()*ref->Py());
 317           Float_t pMC  = TMath::Sqrt(ref->Px()*ref->Px() + ref->Py()*ref->Py()+ref->Pz()*ref->Pz());
 318
 319           if (ptMC<1e-33 || ptMC<ptLowerCut || ptMC>ptUpperCut) continue; // first condition is for tracks not crossing padrow 0
 320
 321           allselected++;
 322
 323           hGood->Fill(ptMC);
 324
 325           for(i=0; i<nEntries; i++){
 326               Int_t ttlabel = event->GetTrack(i)->GetTPCLabel();
 327               if( label==TMath::Abs(ttlabel) && ttlabel<0 ){
 328                  track_fake[itrack_fake++] = label;
 329                  hFake->Fill(ptMC);
 330               }
 331           }
 332
 333           AliESDtrack *esdtrack = 0;
 334           for(i=0; i<nEntries; i++){
 335             esdtrack = event->GetTrack(i);
 336             tlabel   = esdtrack->GetTPCLabel();
 337             if(label==tlabel) break;
 338           }
 339
 340           if(i==nEntries){
 341             track_notfound[itrack_notfound++] = label;
 342             cout << "MC track " << label << " not reconstructed" << endl;
 343             continue;
 344           }
 345
 346           Int_t micount = 0;
 347           Int_t mi;
 348           AliESDtrack *mitrack;
 349
 350           for(mi=0; mi<nEntries; mi++){
 351               mitrack = event->GetTrack(mi);
 352               if(label==mitrack->GetTPCLabel()) micount++;
 353           }
 354
 355           if(micount>1){
 356              track_multifound[itrack_multifound]   = label;
 357              track_multifound_n[itrack_multifound] = micount;
 358              itrack_multifound++;
 359              hClone->Fill(ptMC,micount-1);
 360           }
 361
 362         //if((mitrack->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
 363         //if((mitrack->GetStatus()&AliESDtrack::kTPCin)==0)    continue;
 364           cout << "MC track " << label << " reconstructed "<<micount<<" times" << endl;
 365
 366         if(label==tlabel){
 367            hFound->Fill(ptMC);
 368            //if( micount<1 ) cout<<"ERROR!!!!"<<endl;
 369            hFoundMult->Fill(ptMC,micount);
 370         }
 371
 372         AliExternalTrackParam tpctrack =*(esdtrack->GetInnerParam());
 373
 374         if( TMath::Abs(tpctrack.GetSigned1Pt()) <1.e-8 ) continue;
 375
 376         Double_t mcxyz[3]   = { ref->X(), ref->Y(), ref->Z() };
 377         Double_t mclocal[3] = { ref->X(), ref->Y(), ref->Z() };
 378
 379         if(1){
 380           Double_t c = TMath::Cos(tpctrack.GetAlpha());
 381           Double_t s = TMath::Sin(tpctrack.GetAlpha());
 382           Double_t x = mclocal[0];
 383           mclocal[0] = x*c + mclocal[1]*s;
 384           mclocal[1] =-x*s + mclocal[1]*c;
 385         }
 386
 387         if(0){
 388           Double_t local[3] = { tpctrack.GetX(),tpctrack.GetY(),tpctrack.GetZ() };
 389 //        cout << "label: "      << label << endl;
 390 //        cout << "Z diff = "    << local[2] - mclocal[2] << endl;
 391 //        cout << "orig rec: "   << local[0]<<" "<<local[1]<<" "<<local[2] <<" "<<TMath::Sqrt(local[0]*local[0]+local[1]*local[1])<<endl;
 392 //        cout << "rotated mc: " << mclocal[0]<<" "<<mclocal[1]<<" "<<mclocal[2] <<" "<<TMath::Sqrt(mclocal[0]*mclocal[0]+mclocal[1]*mclocal[1])<<endl;
 393 //        cout << "orig mc: "    << mcxyz[0]<<" "<<mcxyz[1]<<" "<<mcxyz[2]<<endl;
 394 //        cout << "alpha = "     << tpctrack.GetAlpha()/TMath::Pi()*180.<<" "<<ref->Alpha()/TMath::Pi()*180.<<endl;
 395         }
 396
 397         //cout << "X = " << mclocal[0] << " / " << tpctrack.GetX() << endl;
 398         tpctrack.AliExternalTrackParam::PropagateTo(mclocal[0],5.);
 399
 400         Double_t xyz[3], pxpypz[3];
 401
 402         tpctrack.GetXYZ(xyz);       // GetInnerXYZ
 403         tpctrack.GetPxPyPz(pxpypz); // GetInnerPxPyPz
 404
 405         Double_t local[3] = { tpctrack.GetX(),tpctrack.GetY(),tpctrack.GetZ() };
 406         Float_t phiRec = TMath::ATan2(pxpypz[1],pxpypz[0]);
 407
 408         if(phiRec <= TMath::Pi()) phiRec += 2*TMath::Pi();
 409         if(phiRec >= TMath::Pi()) phiRec -= 2*TMath::Pi();
 410
 411         Double_t ptRec = TMath::Sqrt(pxpypz[0]*pxpypz[0]+pxpypz[1]*pxpypz[1]);
 412         Double_t pRec  = TMath::Sqrt(pxpypz[0]*pxpypz[0]+pxpypz[1]*pxpypz[1]+pxpypz[2]*pxpypz[2]);
 413
 414         Float_t  lambdaRec = TMath::ATan2(pxpypz[2],ptRec);
 415         Float_t pt_1 = 1./ptRec;
 416         //Float_t pts_1 = tpctrack->GetSigned1Pt();
 417         //Float_t ptMCs_1 = ref->/ptMC;
 418
 419
 420
 421         Int_t pdg = (Int_t)ref->GetLength();  // particle PDG
 422         const Double_t kCLight = 0.000299792458;
 423         //Double_t Bz = 5.;
 424         Double_t q  = 0.;
 425
 426         if( TMath::Abs(pdg)<10000 ){
 427            TParticlePDG *ppdg = TDatabasePDG::Instance()->GetParticle(pdg);
 428            if( ppdg ) q = ppdg->Charge()/3.;
 429         }
 430
 431         Double_t kappaMC = q/ptMC; // /Bz/kCLight;
 432         hnhit_ref->Fill(esdtrack->GetTPCNcls());
 433
 434         hPt->Fill(ptMC, (ptRec - ptMC)/(ptMC)*100.);
 435
 436         Float_t phiMC = TMath::ATan2(ref->Py(),ref->Px());
 437         hPhi->Fill(ptMC, (phiRec - phiMC)/phiMC*1000.);
 438
 439         Float_t lambdaMC = TMath::ATan2(ref->Pz(),ptMC);
 440         hLambda->Fill(ptMC, (lambdaRec - lambdaMC)/lambdaMC*1000.);
 441
 442         hY->Fill(ptMC, (local[1]-mclocal[1])*10.);
 443         hZ->Fill( fabs(mclocal[2]), (local[2]-mclocal[2] ) *10.); // Please check hY and hZ! KKK
 444
 445         if( tpctrack.GetSigmaY2()>0   && finite(tpctrack.GetSigmaY2())   )  hpullY   ->Fill( (local[1]-mclocal[1])/TMath::Sqrt(TMath::Abs(tpctrack.GetSigmaY2()))  );
 446         if( tpctrack.GetSigmaZ2()>0   && finite(tpctrack.GetSigmaZ2())   )  hpullZ   ->Fill( (local[2]-mclocal[2])/TMath::Sqrt(TMath::Abs(tpctrack.GetSigmaZ2()))  );
 447         if( tpctrack.GetSigmaSnp2()>0 && finite(tpctrack.GetSigmaSnp2()) )  hpullPhi ->Fill( (phiRec-phiMC)/TMath::Sqrt(TMath::Abs(tpctrack.GetSigmaSnp2()))       );
 448         if( tpctrack.GetSigmaTgl2()>0 && finite(tpctrack.GetSigmaTgl2()) )  hpullDzds->Fill( (lambdaRec-lambdaMC)/TMath::Sqrt(TMath::Abs(tpctrack.GetSigmaTgl2())) );
 449
 450         if( tpctrack.GetSigma1Pt2()>0 && finite(tpctrack.GetSigma1Pt2()) )
 451         if( q!=0 ) hpullK->Fill((tpctrack.GetSigned1Pt()-kappaMC)/TMath::Sqrt(TMath::Abs(tpctrack.GetSigma1Pt2())));
 452
 453         Float_t dedx = esdtrack->GetTPCsignal();
 454         hDedx->Fill(pRec,dedx,1.);
 455
 456         if(TMath::Abs(pdg)==211) // pions
 457            if(pRec>0.4 && pRec<0.5){
 458               he->Fill(dedx,1.);
 459            }
 460       } // end of loop over "good" tracks
 461
 462       /*
 463       cout<<"\nList of Not found tracks :\n";
 464       for ( i = 0; i< itrack_notfound; i++){
 465         cout<<track_notfound[i]<<"\t";
 466         if ((i%5)==4) cout<<"\n";
 467       }
 468       cout<<"\nList of fake  tracks :\n";
 469       for ( i = 0; i< itrack_fake; i++){
 470         cout<<track_fake[i]<<"\t";
 471         if ((i%5)==4) cout<<"\n";
 472       }
 473       cout<<"\nList of multiple found tracks :\n";
 474       for ( i=0; i<itrack_multifound; i++) {
 475           cout<<"id.   "<<track_multifound[i]
 476               <<"     found - "<<track_multifound_n[i]<<"times\n";
 477       }
 478
 479       cout<<"Number of found tracks ="<<nEntries<<endl;
 480       cout<<"Number of \"good\" tracks ="<<ngood<<endl;
 481       */
 482
 483   refs->Clear();
 484   } // end of the loop over events
 485
 486    delete event;
 487    delete esdTree;
 488    ef->Close();
 489
 490    delete tpcTree;
 491    refFile->Close();
 492
 493    Stat_t nGoodMC  = hGood->GetEntries();
 494    Stat_t nFakes   = hFake->GetEntries();
 495    Stat_t nClones  = hClone->GetEntries();
 496    Stat_t nRec     = hFound->GetEntries();
 497    Stat_t nRecMult = hFoundMult->GetEntries();
 498
 499
 500     if(nGoodMC!=0)         ::Info("\n\nAliHLTTPCTrackerEvaluation","Integral efficiency in (%) is about : %f\n", nRec/nGoodMC*100.);
 501     if(nRecMult+nFakes!=0) ::Info("AliHLTTPCTrackerEvaluation","Integral fake rate in (%) is about  : %f\n", nFakes/(nRecMult+nFakes)*100.);
 502     if(nRecMult+nFakes!=0) ::Info("AliHLTTPCTrackerEvaluation","Integral clone rate in (%) is about : %f\n", nClones/(nRecMult+nFakes)*100.);
 503
 504     ::Info("AliHLTTPCTrackerEvaluation", "Total number of selected MC tracks                    : %d\n", allgood);
 505     ::Info("AliHLTTPCTrackerEvaluation", "Total number of selected MC tracks with momentum cuts : %d\n", allselected);
 506     ::Info("AliHLTTPCTrackerEvaluation", "Total number of reconstructed tracks                  : %d\n", allfound);
 507
 508
 509    // ---------- Plotting the histograms ----------------------------------------------------------------------------
 510
 511    hFound    ->Sumw2();
 512    hFoundMult->Sumw2();
 513    hGood     ->Sumw2();
 514    hFake     ->Sumw2();
 515    hClone    ->Sumw2();
 516
 517    hEff     ->Divide(hFound,hGood,     1,1.,"b");
 518    hFakeEff ->Divide(hFake, hGood,     1,1.,"b");
 519    hCloneEff->Divide(hClone,hFoundMult,1,1.,"b");
 520    // KKK are you sure this is the correct definition? When plotted, it looks like we have 70% clones in the lower momenta...
 521    // perhaps fakes and clones should be displayed in another histogram??
 522
 523    gROOT->SetStyle("Plain");
 524    gStyle->SetOptStat(""); //gStyle->SetOptStat(111110);
 525    gStyle->SetOptFit(0);  //gStyle->SetOptFit(1);
 526    gStyle->SetPalette(1);
 527
 528
 529    //--------------- canvas 0 for EFFICIENCY -------------------
 530
 531    TCanvas *c0 = new TCanvas("c0","reconstruction efficiency",500,450);
 532
 533    c0->cd();
 534    hEff     ->Draw();
 535    hCloneEff->Draw("histsame");
 536    hFakeEff ->Draw("histsame");
 537
 538    TLegend *leg = new TLegend(0.3,0.5,0.85,0.8);
 539    leg->SetFillColor(10);
 540    leg->SetLineColor(10);
 541    leg->AddEntry(hEff,      "reconstruction eff.", "l");
 542    leg->AddEntry(hCloneEff, "clone contribution", "l");
 543    leg->AddEntry(hFakeEff,  "fake contribution",  "l");
 544    leg->Draw("same");
 545
 546    TLine *line1 = new TLine(0,1,6,1);
 547    line1->SetLineStyle(4);
 548    line1->Draw("same");
 549
 550    Float_t fakeData  = nFakes/(nRecMult+nFakes)*100.;
 551    Float_t cloneData = nClones/(nRecMult+nFakes)*100.;
 552    Char_t  fakeBuf[100];
 553    Char_t  cloneBuf[100];
 554
 555    sprintf(fakeBuf,"%d", (Int_t)fakeData);
 556    TString fakeStr = fakeBuf;
 557    fakeStr.Append(" %");
 558
 559    sprintf(cloneBuf,"%d", (Int_t)cloneData);
 560    TString cloneStr = cloneBuf;
 561    cloneStr.Append(" %");
 562
 563    TText *text = new TText(1.4,0.1,fakeStr.Data());
 564    text->SetTextSize(0.05);
 565    text->SetTextColor(kRed);
 566    text->Draw();
 567    text = new TText(0.6,0.3,cloneStr.Data());
 568    text->SetTextSize(0.05);
 569    text->SetTextColor(kBlue);
 570    text->Draw();
 571
 572
 573    //--------------- canvas 1 for RESOLUTIONS ----------------
 574
 575
 576    gStyle->SetOptStat(0);
 577    TCanvas *c1 = new TCanvas("c1","resolutions",1100,900); // please add the Y and Z projections
 578
 579    TF1 *fGauss = new TF1("gauss","gaus");
 580    for(Int_t i=0; i<15; i++){
 581      Float_t pMin = i*6/15.;// -.5;
 582      Float_t pMax = (i+1)*6/15.;// +.5;
 583      Int_t binx1 = hPt->GetXaxis()->FindBin(pMin);
 584      Int_t binx2 = hPt->GetXaxis()->FindBin(pMax);
 585      Float_t zMin = i*250./15.;
 586      Float_t zMax = (i+1)*250./15.;
 587      Int_t binz1 = hZ->GetXaxis()->FindBin(zMin);
 588      Int_t binz2 = hZ->GetXaxis()->FindBin(zMax);
 589      //if(i>0) binx1++;
 590
 591      TH1D *p;
 592      p = (hPt    ->ProjectionY("", binx1, binx2));
 593      //cout<<i<<" "<<pMin<<" "<<pMax<<" "<<binx1<<" "<<binx2<<": "<<p->GetEntries()<<endl;
 594      if(p->GetEntries()>50){
 595        fGauss->SetParameter(2,0);
 596        p->Fit("gauss","Q");
 597        hResPt->Fill(pMin, fGauss->GetParameter(2));fGauss->GetParError(2);
 598      }
 599
 600      //KKK I am resetting only the sigma and its error, perhaps we should do this with all the fit parameters ???
 601      // SG I don't know;
 602
 603      p    = (hPhi   ->ProjectionY("", binx1, binx2));
 604      if(p->GetEntries()>50){
 605        fGauss->SetParameter(2,0);
 606        p->Fit("gauss","Q");
 607        hResPhi->Fill(pMin, fGauss->GetParameter(2)); fGauss->GetParError(2);
 608      }
 609
 610      p = (hLambda->ProjectionY("", binx1, binx2));
 611      if(p->GetEntries()>50){
 612        fGauss->SetParameter(2,0);
 613        p->Fit("gauss","Q");
 614        hResLambda->Fill(pMin, fGauss->GetParameter(2)); fGauss->GetParError(2);
 615      }
 616
 617      p     = (hY     ->ProjectionY("", binx1, binx2));
 618      if(p->GetEntries()>50){
 619        fGauss->SetParameter(2,0);
 620        p->Fit("gauss","Q");
 621        hResY->Fill(pMin, fGauss->GetParameter(2)); fGauss->GetParError(2);
 622      }
 623
 624      p      = (hZ     ->ProjectionY("", binz1, binz2));
 625      //cout<<i<<" "<<zMin<<" "<<zMax<<" "<<binz1<<" "<<binz2<<endl;
 626      //cout<<p->GetEntries()<<endl;
 627      if(p->GetEntries()>50){
 628        fGauss->SetParameter(2,0);
 629        p->Fit("gauss","Q");
 630        hResZ->Fill(zMin, fGauss->GetParameter(2)); fGauss->GetParError(2);
 631      }
 632    }
 633
 634    c1->Clear(); c1->Divide(3,2);
 635
 636    TVirtualPad *p = c1->cd(1);
 637    p->SetGridy();
 638    hResPt->SetMarkerStyle(20);
 639    hResPt->Draw("P"); // KKK I haven't filled the errors for the sigma, perhaps a TGraph would make them easier to plot
 640
 641    p = c1->cd(2);
 642    p->SetGridy();
 643    hResPhi->SetMarkerStyle(20);
 644    hResPhi->Draw("P");
 645
 646    p=c1->cd(3);
 647    p->SetGridy();
 648    hResLambda->SetMarkerStyle(20);
 649    hResLambda->Draw("P");
 650
 651    p=c1->cd(4);
 652    p->SetGridy();
 653    hResY->SetMarkerStyle(20);
 654    hResY->Draw("P");
 655
 656    p=c1->cd(5);
 657    p->SetGridy();
 658    hResZ->SetMarkerStyle(20);
 659    hResZ->Draw("P");
 660    c1->Update();
 661 // pad 6 stays empty
 662
 663
 664    //----------------- optional canvases 2 and 3 for dE/dx and clusters ----------------
 665
 666    if(bDedxAndClus){
 667       TCanvas *c2 = new TCanvas("c2","dE/dx",500,450);
 668       c2->cd();
 669       hDedx->Draw();
 670
 671       TCanvas *c3 = new TCanvas("c","clusters",500,450);
 672       c3->cd();
 673       hClus->Draw();
 674       hnhit_ref->SetLineColor(kRed);
 675       hnhit_ref->Draw("same");
 676    }
 677
 678
 679    //----------------- canvas 3 for PULL VARIABLES --------------------
 680
 681    TCanvas *c4 = new TCanvas("c4","pull variables",1100,900);
 682
 683    gStyle->SetOptFit(2);
 684    gStyle->SetOptStat("e");
 685
 686    c4->Divide(3,2);
 687
 688    c4->cd(1);
 689    hpullY->Draw();
 690    hpullY->Fit("gaus","Q");
 691
 692    c4->cd(2);
 693    hpullZ->Draw();
 694    hpullZ->Fit("gaus","Q");
 695
 696    c4->cd(3);
 697    hpullPhi->Draw();
 698    hpullPhi->Fit("gaus","Q");
 699
 700    c4->cd(4);
 701    hpullDzds->Draw();
 702    hpullDzds->Fit("gaus","Q");
 703
 704    c4->cd(5);
 705    hpullK->Draw();
 706    hpullK->Fit("gaus","Q");
 707
 708
 709
 710 //-------------------- OUTPUT FILE ----------------------------------------
 711
 712    //TFile fc("CAtrackerEvaluation.root","RECREATE");
 713    //c1->Write();
 714    //c2->Write();
 715    //c3->Write();
 716    //fc.Close(); // KKK commented out temporarily
 717
 718    gBenchmark->Stop("AliHLTTPCTrackerEvaluation");
 719    gBenchmark->Show("AliHLTTPCTrackerEvaluation");
 720
 721    return 0;
 722 }
 723
 724
 725 //============================================================================
 726
 727 Int_t GoodTPCTracks(const Char_t *dir){
 728
 729   Char_t fname[100];
 730   sprintf(fname,"%s/galice.root",dir);
 731
 732   AliRunLoader *runLoader = AliRunLoader::Open(fname,"COMPARISON");
 733   if (!runLoader) {
 734      ::Error("GoodTPCTracks","Cannot start session!");
 735      return 1;
 736   }
 737
 738   // load MC information
 739
 740   runLoader->LoadgAlice();
 741   runLoader->LoadHeader();
 742   runLoader->LoadKinematics();
 743   runLoader->LoadTrackRefs();
 744
 745   AliTPCLoader *tpcLoader = (AliTPCLoader*)runLoader->GetLoader("TPCLoader");
 746   if (tpcLoader == 0x0) {
 747      ::Error("GoodTPCTracks","Cannot find TPCLoader!");
 748      delete runLoader;
 749      return 2;
 750   }
 751
 752   AliTPC *TPC = (AliTPC*)runLoader->GetAliRun()->GetDetector("TPC");
 753   Int_t tpcVersion = TPC->IsVersion();
 754
 755   //cout << "TPC version "<< tpcVersion << " has been found!\n";
 756
 757   if      (tpcVersion==1) tpcLoader->LoadRecPoints();
 758   else if (tpcVersion==2) tpcLoader->LoadDigits();
 759   else {
 760      ::Error("GoodTPCTracks","Wrong TPC version!");
 761      delete runLoader;
 762      return 3;
 763   }
 764
 765   AliCDBManager *manager = AliCDBManager::Instance();
 766   manager->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
 767   manager->SetRun(0);
 768
 769   // loading the TPC parameters
 770   AliTPCParamSR *digp = (AliTPCParamSR*)(AliTPCcalibDB::Instance()->GetParameters());
 771   if (!digp) {
 772     ::Error("AliHLTTPCTrackerEvaluation.C","TPC parameters have not been found!");
 773     delete runLoader;
 774     return 4;
 775   }
 776
 777   Int_t nrow_up     = digp->GetNRowUp();                           //get the number of pad rows in up sector
 778   Int_t nrows       = digp->GetNRowLow()+nrow_up;                  //get the number of pad rows in low sector and add
 779   Int_t zeroSup     = digp->GetZeroSup();
 780   Int_t gap         = Int_t(0.125*nrows), shift = Int_t(0.5*gap);
 781   Int_t good_number = Int_t(0.4*nrows);                            // will be used for selecting tracks with hits in 40% of the rows
 782
 783   Int_t nEvents = runLoader->GetNumberOfEvents();
 784   ::Info("GoodTPCTracks","Number of events : %d\n", nEvents);
 785
 786   sprintf(fname,"%s/GoodTPCTracks.root",dir);
 787   TFile *file = TFile::Open(fname,"recreate");
 788
 789   TClonesArray dummy("AliTrackReference",1000), *refs = &dummy;
 790   TTree tpcTree("tpcTree","Tree with information about the reconstructable TPC tracks");
 791   tpcTree.Branch("TPC",&refs);
 792
 793
 794
 795   //-------------- Loop over generated events ------------------------------
 796
 797   for(Int_t iEvent=0; iEvent<nEvents; iEvent++){
 798
 799       Int_t nt = 0;
 800       refs->Clear();
 801
 802       Int_t i;
 803
 804       runLoader->GetEvent(iEvent);
 805       file->cd();
 806
 807       Int_t nParticles = runLoader->GetHeader()->GetNtrack();
 808       //cout << "Event " << iEvent << ", Number of particles: " << nParticles << endl;
 809
 810       // ------- Set the selection criteria for the MC sample that will be used for the efficiency -------
 811
 812       Int_t *good = new Int_t[nParticles]; for(i=0; i<nParticles; i++) good[i] = 0;
 813
 814       switch (tpcVersion){
 815       case 1: // not used any more
 816          {
 817            AliTPCClustersArray *pca = new AliTPCClustersArray, &ca=*pca;
 818            ca.Setup(digp);
 819            ca.SetClusterType("AliTPCcluster");
 820            ca.ConnectTree(tpcLoader->TreeR());
 821            Int_t nrows=Int_t(ca.GetTree()->GetEntries());
 822            for (Int_t n=0; n<nrows; n++) {
 823              AliSegmentID *s=ca.LoadEntry(n);
 824              Int_t sec,row;
 825              digp->AdjustSectorRow(s->GetID(),sec,row);
 826              AliTPCClustersRow &clrow = *ca.GetRow(sec,row);
 827              Int_t ncl=clrow.GetArray()->GetEntriesFast();
 828              while (ncl--) {
 829                  AliTPCcluster *c=(AliTPCcluster*)clrow[ncl];
 830                  Int_t label=c->GetLabel(0);
 831                  if (label<0) continue; //noise cluster
 832                  label=TMath::Abs(label);
 833
 834                  if (sec>=digp->GetNInnerSector())
 835                     if (row==nrow_up-1) good[label]|=0x4000;
 836                  if (sec>=digp->GetNInnerSector())
 837                     if (row==nrow_up-1-gap) good[label]|=0x1000;
 838
 839                  if (sec>=digp->GetNInnerSector())
 840                     if (row==nrow_up-1-shift) good[label]|=0x2000;
 841                  if (sec>=digp->GetNInnerSector())
 842                     if (row==nrow_up-1-gap-shift) good[label]|=0x800;
 843
 844                  good[label]++;
 845              }
 846              ca.ClearRow(sec,row);
 847            }
 848            delete pca;
 849            }
 850          break;
 851
 852       case 2:
 853          {
 854            TTree *TD = tpcLoader->TreeD(); // get the digits tree
 855
 856            AliSimDigits da, *digits = &da;
 857            TD->GetBranch("Segment")->SetAddress(&digits);
 858
 859            Int_t *count = new Int_t[nParticles];
 860            Int_t i;
 861            for (i=0; i<nParticles; i++) count[i] = 0;
 862
 863            Int_t sectors_by_rows = (Int_t)TD->GetEntries();
 864
 865            for(i=0; i<sectors_by_rows; i++){
 866
 867              if (!TD->GetEvent(i)) continue;
 868              Int_t sec,row;
 869
 870              digp->AdjustSectorRow(digits->GetID(),sec,row);
 871              if(digits->First()){
 872                do {
 873                  Int_t it    = digits->CurrentRow();
 874                  Int_t ip    = digits->CurrentColumn();
 875                  Short_t dig = digits->GetDigit(it,ip);
 876                  Int_t idx0  = digits->GetTrackID(it,ip,0);
 877                  Int_t idx1  = digits->GetTrackID(it,ip,1);
 878                  Int_t idx2  = digits->GetTrackID(it,ip,2);
 879                  if(idx0>=0 && dig>=zeroSup && idx0<nParticles) count[idx0]+=1;
 880                  if(idx1>=0 && dig>=zeroSup && idx1<nParticles) count[idx1]+=1;
 881                  if(idx2>=0 && dig>=zeroSup && idx2<nParticles) count[idx2]+=1;
 882                } while (digits->Next());
 883              }
 884
 885              for(Int_t j=0; j<nParticles; j++){
 886
 887                  if(count[j]>1){
 888                     if(sec>=digp->GetNInnerSector())
 889                     if(row==nrow_up-1)               good[j]|=0x4000;
 890                     if(sec>=digp->GetNInnerSector())
 891                     if(row==nrow_up-1-gap)           good[j]|=0x1000;
 892
 893                     if(sec>=digp->GetNInnerSector())
 894                     if(row==nrow_up-1-shift)         good[j]|=0x2000;
 895                     if(sec>=digp->GetNInnerSector())
 896                     if(row==nrow_up-1-gap-shift)     good[j]|=0x800;
 897
 898                     good[j]++;
 899                  }
 900                  count[j] = 0;
 901              } // end of loop over particles
 902            } // end of loop over sectors_by_rows
 903            delete [] count;
 904          } // end of case 2
 905          break;
 906       } // end of switch
 907
 908
 909
 910
 911     // ---------------- Select sample of MC particles that will be used for forming the efficiency ------------------
 912
 913     AliStack *stack = runLoader->Stack();
 914
 915     for(i=0; i<nParticles; i++){ // loop over stack
 916
 917         if ((good[i]&0x5000) != 0x5000)//SG!!!
 918         if ((good[i]&0x2800) != 0x2800)     continue;
 919         if ((good[i]&0x7FF ) < good_number) continue;
 920         // N TPC rows with digits >= 64 => common checks of cluster finder & tracker
 921         // There are digits in rows (159&&139) || (149&&129)
 922
 923
 924         TParticle *part = (TParticle*)stack->Particle(i);
 925         if(part == 0x0){
 926            cerr << "Cannot get particle "<< i << endl;
 927            continue;
 928         }
 929
 930         if(part->Pt()<0.100) continue;                        // skip particles with pt below 0.1 GeV/c
 931         if(TMath::Abs(part->Pz()/part->Pt())>0.999) continue; // reject tracks outside the TPC acceptance, below 45 degrees
 932
 933         Double_t vx = part->Vx();
 934         Double_t vy = part->Vy();
 935         Double_t vz = part->Vz();
 936         //if (TMath::Sqrt(vx*vx+vy*vy)>3.5) continue; // vertex cuts
 937         //if (TMath::Abs(vz) > 50.) continue;
 938
 939         AliTrackReference *ref = new((*refs)[nt])AliTrackReference();
 940
 941         ref->SetLabel(i);
 942         Int_t pdg = part->GetPdgCode();
 943         ref->SetLength(pdg);  //This will the particle's PDG !
 944         ref->SetMomentum(0.,0.,0.);
 945         ref->SetPosition(0.,0.,0.);
 946
 947         nt++;
 948     } // end of loop over stack
 949
 950     //----------- Check whether there is also information at the entrance of the TPC
 951
 952     TTree   *TR     = runLoader->TreeTR();
 953     TBranch *branch = TR->GetBranch("TrackReferences");
 954     if(branch==0){
 955        ::Error("GoodTPCTracks","No track references!");
 956        delete runLoader;
 957        return 5;
 958     }
 959
 960     TClonesArray tpcdummy("AliTrackReference",1000), *tpcRefs = &tpcdummy;
 961     branch->SetAddress(&tpcRefs);
 962
 963     for(Int_t r=0; r<(Int_t)TR->GetEntries(); r++){
 964
 965         //cerr<<r<<' '<<(Int_t)TR->GetEntries()<<'\r';
 966         TR->GetEvent(r);
 967
 968         if(!tpcRefs->GetEntriesFast()) continue;
 969
 970         AliTrackReference *tpcRef = 0x0;
 971         for(Int_t iref=0; iref<tpcRefs->GetEntriesFast(); ++iref){
 972
 973             tpcRef = (AliTrackReference*)tpcRefs->UncheckedAt(iref);
 974             if(tpcRef->DetectorId() == AliTrackReference::kTPC) break;
 975             tpcRef = 0x0;
 976         }
 977
 978         if(!tpcRef) continue;
 979
 980         Int_t j;
 981         AliTrackReference *ref = 0;
 982
 983         for(j=0; j<nt; j++){
 984
 985           ref = (AliTrackReference*)refs->UncheckedAt(j);
 986           if(ref->Label()==tpcRef->Label()) break;
 987           ref = 0;
 988         }
 989         if(ref==0) continue;
 990
 991         ref->SetMomentum(tpcRef->Px(),tpcRef->Py(),tpcRef->Pz());
 992         ref->SetPosition(tpcRef->X(), tpcRef->Y(), tpcRef->Z());
 993
 994         tpcRefs->Clear();
 995     }
 996
 997     tpcTree.Fill();
 998
 999     delete [] good;
1000   } //---------- end of the loop over generated events
1001
1002   tpcTree.Write();
1003   file->Close();
1004
1005   delete runLoader;
1006   return 0;
1007 }
1008 #endif