PWG0/dNdEta/makeCorrection.C

   1 //
   2 // Script to make correction maps for dndeta measurements using the
   3 // dNdEtaCorrection class.
   4 //
   5
   6 makeCorrection(Char_t* dataDir, Int_t nRuns=20) {
   7
   8   Char_t str[256];
   9
  10   gSystem->Load("../esdTrackCuts/libESDtrackQuality.so");
  11   gSystem->Load("libdNdEta.so");
  12
  13   // ########################################################
  14   // selection of esd tracks
  15   ESDtrackQualityCuts* esdTrackCuts = new ESDtrackQualityCuts();
  16   esdTrackCuts->DefineHistograms(1);
  17
  18   esdTrackCuts->SetMinNClustersTPC(50);
  19   esdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
  20   esdTrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
  21   esdTrackCuts->SetRequireTPCRefit(kTRUE);
  22
  23   esdTrackCuts->SetMinNsigmaToVertex(3);
  24   esdTrackCuts->SetAcceptKingDaughters(kFALSE);
  25
  26   AliLog::SetClassDebugLevel("ESDtrackQualityCuts",1);
  27
  28   // ########################################################
  29   // definition of dNdEta correction object
  30
  31   dNdEtaCorrection* dNdEtaMap = new dNdEtaCorrection();
  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   // definition of used pointers
  45   TFile* esdFile;
  46   TTree* esdTree;
  47   TBranch* esdBranch;
  48
  49   AliESD* esd =0;
  50
  51   // ########################################################
  52   // loop over runs
  53   Int_t nRunCounter = 0;
  54   for (Int_t r=0; r<nDirs; r++) {
  55
  56     TSystemFile* presentDir = (TSystemFile*)dirList->At(r);
  57     if (!presentDir || !presentDir->IsDirectory())
  58       continue;
  59     // check that the files are there
  60     TString currentDataDir;
  61     currentDataDir.Form("%s/%s",dataDir, presentDir->GetName());
  62     cout << "Processing directory " << currentDataDir.Data() << endl;
  63     if ((!gSystem->Which(currentDataDir,"galice.root")) ||
  64           (!gSystem->Which(currentDataDir,"AliESDs.root")))
  65       continue;
  66
  67     if (nRunCounter++ >= nRuns)
  68       break;
  69
  70     cout << "run #" << nRunCounter << endl;
  71
  72     // #########################################################
  73     // setup galice and runloader
  74     if (gAlice) {
  75       delete gAlice->GetRunLoader();
  76       delete gAlice;
  77       gAlice=0;
  78     }
  79
  80     sprintf(str,"%s/galice.root",currentDataDir.Data());
  81     AliRunLoader* runLoader = AliRunLoader::Open(str);
  82
  83     runLoader->LoadgAlice();
  84     gAlice = runLoader->GetAliRun();
  85     runLoader->LoadKinematics();
  86     runLoader->LoadHeader();
  87
  88     // #########################################################
  89     // open esd file and get the tree
  90
  91     sprintf(str,"%s/AliESDs.root",currentDataDir.Data());
  92     // close it first to avoid memory leak
  93     if (esdFile)
  94       if (esdFile->IsOpen())
  95         esdFile->Close();
  96
  97     esdFile = TFile::Open(str);
  98     esdTree = (TTree*)esdFile->Get("esdTree");
  99     if (!esdTree)
 100       continue;
 101     esdBranch = esdTree->GetBranch("ESD");
 102     esdBranch->SetAddress(&esd);
 103     if (!esdBranch)
 104       continue;
 105
 106     // ########################################################
 107     // Magnetic field
 108     AliTracker::SetFieldMap(gAlice->Field(),kTRUE); // kTRUE means uniform magnetic field
 109     //AliKalmanTrack::SetConvConst(1000/0.299792458/5.);
 110
 111     // ########################################################
 112     // getting number of events
 113     Int_t nEvents = (Int_t)runLoader->GetNumberOfEvents();
 114     Int_t nESDEvents = esdBranch->GetEntries();
 115
 116     if (nEvents!=nESDEvents) {
 117       cout << " Different number of events from runloader and esdtree!!!" << nEvents << " / " << nESDEvents << endl;
 118       return;
 119     }
 120     // ########################################################
 121     // loop over number of events
 122     cout << " looping over events..." << endl;
 123     for(Int_t i=0; i<nEvents; i++) {
 124
 125       esdBranch->GetEntry(i);
 126       runLoader->GetEvent(i);
 127
 128       // ########################################################
 129       // get the EDS vertex
 130       AliESDVertex* vtxESD = esd->GetVertex();
 131
 132       Double_t vtx[3];
 133       Double_t vtx_res[3];
 134       vtxESD->GetXYZ(vtx);
 135
 136       vtx_res[0] = vtxESD->GetXRes();
 137       vtx_res[1] = vtxESD->GetYRes();
 138       vtx_res[2] = vtxESD->GetZRes();
 139
 140       // the vertex should be reconstructed
 141       if (strcmp(vtxESD->GetName(),"default")==0)
 142         continue;
 143
 144       // the resolution should be reasonable???
 145       if (vtx_res[2]==0 || vtx_res[2]>0.1)
 146         continue;
 147
 148
 149       // ########################################################
 150       // get the MC vertex
 151       AliGenPythiaEventHeader* genHeader =
 152         (AliGenPythiaEventHeader*)runLoader->GetHeader()->GenEventHeader();
 153
 154       TArrayF vtxMC(3);
 155       genHeader->PrimaryVertex(vtxMC);
 156       Double_t vtx[3];
 157       vtx[0] = vtxMC[0];
 158       vtx[1] = vtxMC[1];
 159       vtx[2] = vtxMC[2];
 160
 161
 162       // ########################################################
 163       // loop over mc particles
 164       AliStack* particleStack = runLoader->Stack();
 165       Int_t nPrim    = particleStack->GetNprimary();
 166
 167       for (Int_t i_mc=0; i_mc<nPrim; i_mc++) {
 168
 169         TParticle* part = particleStack->Particle(i_mc);
 170         if (!part || strcmp(part->GetName(),"XXX")==0)
 171           continue;
 172
 173         TParticlePDG* pdgPart = part->GetPDG();
 174
 175         Bool_t prim = kFALSE;
 176         // check if it's a primary particle - is there a better way ???
 177         if ((part->GetFirstDaughter() >= nPrim) || (part->GetFirstDaughter()==-1)) {
 178           if (TMath::Abs(pdgPart->PdgCode())>10 && pdgPart->PdgCode()!=21 && strcmp(pdgPart->ParticleClass(),"Unknown")!=0)
 179             prim = kTRUE;
 180         }
 181         if (!prim)
 182           continue;
 183
 184         if (pdgPart->Charge()==0)
 185           continue;
 186
 187         Float_t eta = part->Eta();
 188
 189         if (prim)
 190           dNdEtaMap->FillGene(vtx[2], eta);
 191
 192       }// end of mc particle
 193
 194       // ########################################################
 195       // loop over esd tracks
 196       Int_t nTracks = esd->GetNumberOfTracks();
 197       for (Int_t t=0; t<nTracks; t++) {
 198         AliESDtrack* esdTrack = esd->GetTrack(t);
 199
 200         // cut the esd track?
 201         if (!esdTrackCuts->AcceptTrack(esdTrack))
 202           continue;
 203
 204         AliTPCtrack* tpcTrack = new AliTPCtrack(*esdTrack);
 205         if (tpcTrack->GetAlpha()==0) {
 206           cout << " Warning esd track alpha = 0" << endl;
 207           continue;
 208         }
 209
 210         Float_t theta = tpcTrack->Theta();
 211         Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
 212
 213         // using the eta of the mc particle
 214         Int_t label = TMath::Abs(esdTrack->GetLabel());
 215         if (label<0) {
 216           cout << " cannot find corresponding mc part !!! " << label << endl;
 217           continue;
 218         }
 219         TParticle* mcPart = particleStack->Particle(label);
 220         eta = mcPart->Eta();
 221
 222         dNdEtaMap->FillMeas(vtx[2], eta);
 223
 224       } // end of track loop
 225     } // end  of event loop
 226   } // end of run loop
 227
 228   dNdEtaMap->Finish();
 229
 230   TFile* fout = new TFile("correction_map.root","RECREATE");
 231
 232   esdTrackCuts->SaveHistograms("esd_track_cuts");
 233   dNdEtaMap->SaveHistograms();
 234
 235   fout->Write();
 236   fout->Close();
 237
 238   dNdEtaMap->DrawHistograms();
 239
 240 }