PWG0/dNdEta/testAnalysis.C

   1 //
   2 // Script to test the dN/dEta analysis using the dNdEtaAnalysis and
   3 // dNdEtaCorrection classes. Note that there is a cut on the events,
   4 // so the measurement will be biassed.
   5 //
   6
   7 testAnalysis(Char_t* dataDir, Int_t nRuns=20) {
   8
   9   Char_t str[256];
  10
  11   gSystem->Load("../esdTrackCuts/libESDtrackQuality.so");
  12   gSystem->Load("libdNdEta.so");
  13
  14   // ########################################################
  15   // selection of esd tracks
  16   ESDtrackQualityCuts* esdTrackCuts = new ESDtrackQualityCuts();
  17   esdTrackCuts->DefineHistograms(1);
  18
  19   esdTrackCuts->SetMinNClustersTPC(50);
  20   esdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
  21   esdTrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
  22   esdTrackCuts->SetRequireTPCRefit(kTRUE);
  23
  24   esdTrackCuts->SetMinNsigmaToVertex(3);
  25   esdTrackCuts->SetAcceptKingDaughters(kFALSE);
  26
  27   AliLog::SetClassDebugLevel("ESDtrackQualityCuts",1);
  28
  29   // ########################################################
  30   // definition of dNdEta objects
  31
  32   dNdEtaCorrection* dNdEtaMap = new dNdEtaCorrection();
  33
  34   dNdEtaMap->LoadHistograms("correction_map.root","dndeta_correction");
  35   dNdEtaMap->RemoveEdges(2,0,2);
  36
  37   dNdEtaAnalysis* analyse = new dNdEtaAnalysis("dndeta");
  38
  39   // ########################################################
  40   // get the data dir
  41   Char_t execDir[256];
  42   sprintf(execDir,gSystem->pwd());
  43   TSystemDirectory* baseDir = new TSystemDirectory(".",dataDir);
  44   TList* dirList            = baseDir->GetListOfFiles();
  45   Int_t nDirs               = dirList->GetEntries();
  46   // go back to the dir where this script is executed
  47   gSystem->cd(execDir);
  48
  49
  50   // ########################################################
  51   // definition of used pointers
  52   TFile* esdFile;
  53   TTree* esdTree;
  54   TBranch* esdBranch;
  55
  56   AliESD* esd =0;
  57
  58   // ########################################################
  59   // loop over runs
  60   Int_t nRunCounter = 0;
  61   for (Int_t r=1; r<=nDirs; r++) {
  62
  63     TSystemFile* presentDir = (TSystemFile*)dirList->At(r);
  64     if (!presentDir->IsDirectory())
  65       continue;
  66     // check that the files are there
  67     TString currentDataDir;
  68     currentDataDir.Form("%s/%s",dataDir, presentDir->GetName());
  69     if ((!gSystem->Which(currentDataDir.Data(),"galice.root")) ||
  70         (!gSystem->Which(currentDataDir.Data(),"AliESDs.root")))
  71       continue;
  72
  73     if (nRunCounter++ >= nRuns)
  74       break;
  75
  76     cout << "run #" << nRunCounter << endl;
  77
  78     // #########################################################
  79     // setup galice and runloader
  80     if (gAlice) {
  81       delete gAlice->GetRunLoader();
  82       delete gAlice;
  83       gAlice=0;
  84     }
  85
  86     sprintf(str,"%s/galice.root",currentDataDir.Data());
  87     AliRunLoader* runLoader = AliRunLoader::Open(str);
  88
  89     runLoader->LoadgAlice();
  90     gAlice = runLoader->GetAliRun();
  91     runLoader->LoadKinematics();
  92     runLoader->LoadHeader();
  93
  94     // #########################################################
  95     // open esd file and get the tree
  96
  97     sprintf(str,"%s/AliESDs.root",currentDataDir.Data());
  98     // close it first to avoid memory leak
  99     if (esdFile)
 100       if (esdFile->IsOpen())
 101         esdFile->Close();
 102
 103     esdFile = TFile::Open(str);
 104     esdTree = (TTree*)esdFile->Get("esdTree");
 105     if (!esdTree)
 106       continue;
 107     esdBranch = esdTree->GetBranch("ESD");
 108     esdBranch->SetAddress(&esd);
 109     if (!esdBranch)
 110       continue;
 111
 112     // ########################################################
 113     // Magnetic field
 114     AliTracker::SetFieldMap(gAlice->Field(),kTRUE); // kTRUE means uniform magnetic field
 115     //AliKalmanTrack::SetConvConst(1000/0.299792458/5.);
 116
 117     // ########################################################
 118     // getting number of events
 119     Int_t nEvents = (Int_t)runLoader->GetNumberOfEvents();
 120     Int_t nESDEvents = esdBranch->GetEntries();
 121
 122     if (nEvents!=nESDEvents) {
 123       cout << " Different number of events from runloader and esdtree!!!" << nEvents << " / " << nESDEvents << endl;
 124       return;
 125     }
 126     // ########################################################
 127     // loop over number of events
 128     cout << " looping over events..." << endl;
 129     for(Int_t i=1; i<nEvents; i++) {
 130
 131       esdBranch->GetEntry(i);
 132       runLoader->GetEvent(i);
 133
 134
 135       // ########################################################
 136       // get the EDS vertex
 137       AliESDVertex* vtxESD = esd->GetVertex();
 138
 139       Double_t vtx[3];
 140       Double_t vtx_res[3];
 141       vtxESD->GetXYZ(vtx);
 142
 143       vtx_res[0] = vtxESD->GetXRes();
 144       vtx_res[1] = vtxESD->GetYRes();
 145       vtx_res[2] = vtxESD->GetZRes();
 146
 147       // we need a good vertex
 148       //  => there will be a bias on the measurement, since this cuts away some events
 149       if (strcmp(vtxESD->GetName(),"default")==0)
 150         continue;
 151       if (vtx_res[2]==0 || vtx_res[2]>0.1)
 152         continue;
 153
 154       // ########################################################
 155       // loop over esd tracks
 156       Int_t nTracks = esd->GetNumberOfTracks();
 157       for (Int_t t=0; t<nTracks; t++) {
 158         AliESDtrack* esdTrack = esd->GetTrack(t);
 159
 160         if (!esdTrackCuts->AcceptTrack(esdTrack))
 161           continue;
 162
 163         AliTPCtrack* tpcTrack = new AliTPCtrack(*esdTrack);
 164
 165         if (tpcTrack->GetAlpha()==0) {
 166           cout << " Warning esd track alpha = 0" << endl;
 167           continue;
 168         }
 169
 170         Float_t theta = tpcTrack->Theta();
 171         Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
 172         Float_t correction = dNdEtaMap->GetCorrection(vtx[2], eta);
 173
 174         dNdEtaMap->FillMeas(vtx[2], eta);
 175
 176         analyse   ->FillTrack(vtx[2], eta, correction);
 177
 178       } // end of track loop
 179       analyse->FillEvent(vtx[2]);
 180
 181     } // end  of event loop
 182   } // end of run loop
 183
 184   analyse->Finish();
 185
 186   TFile* fout = new TFile("out.root","RECREATE");
 187
 188   esdTrackCuts->SaveHistograms("esd_tracks_cuts");
 189   dNdEtaMap->SaveHistograms();
 190   analyse   ->SaveHistograms();
 191
 192   fout->Write();
 193   fout->Close();
 194
 195 }