PWG0/esdTrackCuts/testESDtrackCuts.C

   1 testESDtrackCuts(Char_t* dataDir=, Int_t nRuns=10) {
   2
   3   Char_t str[256];
   4
   5   gSystem->Load("libESDtrackQuality.so");
   6
   7   // ########################################################
   8   // definition of ESD track cuts
   9
  10   ESDtrackQualityCuts* trackCuts = new ESDtrackQualityCuts();
  11   trackCuts->DefineHistograms(4);
  12
  13   trackCuts->SetMinNClustersTPC(50);
  14   trackCuts->SetMaxChi2PerClusterTPC(3.5);
  15   trackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
  16   trackCuts->SetRequireTPCRefit(kTRUE);
  17
  18   trackCuts->SetMinNsigmaToVertex(3);
  19   trackCuts->SetAcceptKingDaughters(kFALSE);
  20
  21   AliLog::SetClassDebugLevel("ESDtrackQualityCuts",1);
  22
  23   // ########################################################
  24   // definition of used pointers
  25   TFile* esdFile;
  26   TTree* esdTree;
  27   TBranch* esdBranch;
  28
  29   AliESD* esd = 0;
  30
  31   // ########################################################
  32   // get the data dir
  33   Char_t execDir[256];
  34   sprintf(execDir,gSystem->pwd());
  35   TSystemDirectory* baseDir = new TSystemDirectory(".",dataDir);
  36   TList* dirList            = baseDir->GetListOfFiles();
  37   Int_t nDirs               = dirList->GetEntries();
  38   // go back to the dir where this script is executed
  39   gSystem->cd(execDir);
  40
  41   // ########################################################
  42   // loop over runs
  43   Int_t nRunCounter = 0;
  44   for (Int_t r=1; r<=nDirs; r++) {
  45
  46     TSystemFile* presentDir = (TSystemFile*)dirList->At(r);
  47     if (!presentDir->IsDirectory())
  48       continue;
  49     // first check that the files are there
  50     sprintf(str,"%s/%s",dataDir, presentDir->GetName());
  51     if ((!gSystem->Which(str,"galice.root")) ||
  52         (!gSystem->Which(str,"AliESDs.root")))
  53       continue;
  54
  55     if (nRunCounter++ >= nRuns)
  56       break;
  57
  58     cout << "run #" << nRunCounter << endl;
  59
  60     // #########################################################
  61     // setup galice and runloader
  62     if (gAlice) {
  63       delete gAlice->GetRunLoader();
  64       delete gAlice;
  65       gAlice=0;
  66     }
  67
  68     sprintf(str,"%s/run%d/galice.root",dataDir,r);
  69     AliRunLoader* runLoader = AliRunLoader::Open(str);
  70
  71     runLoader->LoadgAlice();
  72     gAlice = runLoader->GetAliRun();
  73     runLoader->LoadHeader();
  74
  75     // #########################################################
  76     // open esd file and get the tree
  77
  78     sprintf(str,"%s/run%d/AliESDs.root",dataDir,r);
  79     // close it first to avoid memory leak
  80     if (esdFile)
  81       if (esdFile->IsOpen())
  82         esdFile->Close();
  83
  84     esdFile = TFile::Open(str);
  85     esdTree = (TTree*)esdFile->Get("esdTree");
  86     if (!esdTree)
  87       continue;
  88     esdBranch = esdTree->GetBranch("ESD");
  89     esdBranch->SetAddress(&esd);
  90     if (!esdBranch)
  91       continue;
  92
  93     // ########################################################
  94     // Magnetic field
  95     AliTracker::SetFieldMap(gAlice->Field(),kTRUE); // kTRUE means uniform magnetic field
  96
  97     // ########################################################
  98     // getting number of events
  99     Int_t nEvents    = (Int_t)runLoader->GetNumberOfEvents();
 100     Int_t nESDEvents = esdBranch->GetEntries();
 101
 102     if (nEvents!=nESDEvents)
 103       cout << " Warning: Different number of events from runloader and esdtree!!!" << nEvents << " / " << nESDEvents << endl;
 104
 105     // ########################################################
 106     // loop over number of events
 107     cout << " looping over events..." << endl;
 108     for(Int_t i=1; i<nEvents; i++) {
 109
 110       esdBranch->GetEntry(i);
 111       runLoader->GetEvent(i);
 112
 113       // ########################################################
 114       // get the EDS vertex
 115       AliESDVertex* vtxESD = esd->GetVertex();
 116
 117       Double_t vtxSigma[3];
 118       vtxESD->GetSigmaXYZ(vtxSigma);
 119
 120       // ########################################################
 121       // loop over esd tracks
 122       Int_t nTracks = esd->GetNumberOfTracks();
 123
 124       for (Int_t t=0; t<nTracks; t++) {
 125         AliESDtrack* esdTrack = esd->GetTrack(t);
 126
 127         //trackCuts->AcceptTrack(esdTrack, vtxESD, esd->GetMagneticField());
 128         trackCuts->AcceptTrack(esdTrack);
 129
 130       } // end of track loop
 131     } // end  of event loop
 132   } // end of run loop
 133
 134   TFile* fout = new TFile("out.root","RECREATE");
 135
 136   trackCuts->SaveHistograms("esd_track_cuts");
 137
 138   fout->Write();
 139   fout->Close();
 140
 141 }