PWG0/esdTrackCuts/testESDtrackCuts.C

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