PWG0/esdTrackCuts/testESDtrackCuts.C

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