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("../libPWG0base.so");
  11
  12   // ########################################################
  13   // selection of esd tracks
  14   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts();
  15   esdTrackCuts->DefineHistograms(1);
  16
  17   esdTrackCuts->SetMinNClustersTPC(50);
  18   esdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
  19   esdTrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
  20   esdTrackCuts->SetRequireTPCRefit(kTRUE);
  21
  22   esdTrackCuts->SetMinNsigmaToVertex(3);
  23   esdTrackCuts->SetAcceptKingDaughters(kFALSE);
  24
  25   AliLog::SetClassDebugLevel("AliESDtrackCuts",1);
  26
  27   // ########################################################
  28   // definition of dNdEta correction object
  29
  30   AlidNdEtaCorrection* dNdEtaMap = new AlidNdEtaCorrection();
  31
  32   // ########################################################
  33   // get the data dir
  34   Char_t execDir[256];
  35   sprintf(execDir,gSystem->pwd());
  36   TSystemDirectory* baseDir = new TSystemDirectory(".",dataDir);
  37   TList* dirList            = baseDir->GetListOfFiles();
  38   Int_t nDirs               = dirList->GetEntries();
  39   // go back to the dir where this script is executed
  40   gSystem->cd(execDir);
  41
  42   // ########################################################
  43   // definition of used pointers
  44   TFile* esdFile;
  45   TTree* esdTree;
  46   TBranch* esdBranch;
  47
  48   AliESD* esd =0;
  49
  50   // ########################################################
  51   // loop over runs
  52   Int_t nRunCounter = 0;
  53   for (Int_t r=0; r<nDirs; r++) {
  54
  55     TSystemFile* presentDir = (TSystemFile*)dirList->At(r);
  56     if (!presentDir || !presentDir->IsDirectory())
  57       continue;
  58     // check that the files are there
  59     TString currentDataDir;
  60     currentDataDir.Form("%s/%s",dataDir, presentDir->GetName());
  61     cout << "Processing directory " << currentDataDir.Data() << endl;
  62     if ((!gSystem->Which(currentDataDir,"galice.root")) ||
  63           (!gSystem->Which(currentDataDir,"AliESDs.root")))
  64       continue;
  65
  66     if (nRunCounter++ >= nRuns)
  67       break;
  68
  69     cout << "run #" << nRunCounter << endl;
  70
  71     // #########################################################
  72     // setup galice and runloader
  73     if (gAlice) {
  74       delete gAlice->GetRunLoader();
  75       delete gAlice;
  76       gAlice=0;
  77     }
  78
  79     sprintf(str,"%s/galice.root",currentDataDir.Data());
  80     AliRunLoader* runLoader = AliRunLoader::Open(str);
  81
  82     runLoader->LoadgAlice();
  83     gAlice = runLoader->GetAliRun();
  84     runLoader->LoadKinematics();
  85     runLoader->LoadHeader();
  86
  87     // #########################################################
  88     // open esd file and get the tree
  89
  90     sprintf(str,"%s/AliESDs.root",currentDataDir.Data());
  91     // close it first to avoid memory leak
  92     if (esdFile)
  93       if (esdFile->IsOpen())
  94         esdFile->Close();
  95
  96     esdFile = TFile::Open(str);
  97     esdTree = (TTree*)esdFile->Get("esdTree");
  98     if (!esdTree)
  99       continue;
 100     esdBranch = esdTree->GetBranch("ESD");
 101     esdBranch->SetAddress(&esd);
 102     if (!esdBranch)
 103       continue;
 104
 105     // ########################################################
 106     // Magnetic field
 107     AliTracker::SetFieldMap(gAlice->Field(),kTRUE); // kTRUE means uniform magnetic field
 108     //AliKalmanTrack::SetConvConst(1000/0.299792458/5.);
 109
 110     // ########################################################
 111     // getting number of events
 112     Int_t nEvents = (Int_t)runLoader->GetNumberOfEvents();
 113     Int_t nESDEvents = esdBranch->GetEntries();
 114
 115     Int_t nEventsTriggered = 0;
 116     Int_t nEventsAll       = 0;
 117
 118     if (nEvents!=nESDEvents) {
 119       cout << " Different number of events from runloader and esdtree!!!" << nEvents << " / " << nESDEvents << endl;
 120       return;
 121     }
 122    // ########################################################
 123     // loop over number of events
 124     cout << " looping over events..." << endl;
 125     for(Int_t i=0; i<nEvents; i++) {
 126
 127       esdBranch->GetEntry(i);
 128       runLoader->GetEvent(i);
 129
 130       // ########################################################
 131       // get the EDS vertex
 132       AliESDVertex* vtxESD = esd->GetVertex();
 133
 134       Double_t vtx[3];
 135       Double_t vtx_res[3];
 136       vtxESD->GetXYZ(vtx);
 137
 138       vtx_res[0] = vtxESD->GetXRes();
 139       vtx_res[1] = vtxESD->GetYRes();
 140       vtx_res[2] = vtxESD->GetZRes();
 141
 142       Bool_t vertexReconstructed = kTRUE;
 143
 144       // the vertex should be reconstructed
 145       if (strcmp(vtxESD->GetName(),"default")==0)
 146         vertexReconstructed = kFALSE;
 147
 148       // the resolution should be reasonable???
 149       if (vtx_res[2]==0 || vtx_res[2]>0.01)
 150         vertexReconstructed = kFALSE;
 151
 152       // ########################################################
 153       // get the trigger info
 154
 155       Bool_t triggered = kFALSE;
 156
 157       // MB should be
 158       // ITS_SPD_GFO_L0  : 32
 159       // VZERO_OR_LEFT   : 1
 160       // VZERO_OR_RIGHT  : 2
 161
 162       ULong64_t triggerMask = esd->GetTriggerMask();
 163
 164       nEventsAll++;
 165
 166       if (triggerMask&32 && ((triggerMask&1) || (triggerMask&2))) {
 167         triggered = kTRUE;
 168         nEventsTriggered++;
 169       }
 170
 171       // ########################################################
 172       // get the MC vertex
 173       AliGenPythiaEventHeader* genHeader =
 174         (AliGenPythiaEventHeader*)runLoader->GetHeader()->GenEventHeader();
 175
 176       TArrayF vtxMC(3);
 177       genHeader->PrimaryVertex(vtxMC);
 178       Double_t vtx[3];
 179       vtx[0] = vtxMC[0];
 180       vtx[1] = vtxMC[1];
 181       vtx[2] = vtxMC[2];
 182
 183       // ########################################################
 184       // loop over mc particles
 185       AliStack* particleStack = runLoader->Stack();
 186       Int_t nPrim    = particleStack->GetNprimary();
 187
 188       for (Int_t i_mc=0; i_mc<nPrim; i_mc++) {
 189
 190         TParticle* part = particleStack->Particle(i_mc);
 191         if (!part || strcmp(part->GetName(),"XXX")==0)
 192           continue;
 193
 194         TParticlePDG* pdgPart = part->GetPDG();
 195
 196         Bool_t prim = kFALSE;
 197         // check if it's a primary particle - is there a better way ???
 198         if ((part->GetFirstDaughter() >= nPrim) || (part->GetFirstDaughter()==-1)) {
 199           if (TMath::Abs(pdgPart->PdgCode())>10 && pdgPart->PdgCode()!=21 && strcmp(pdgPart->ParticleClass(),"Unknown")!=0)
 200             prim = kTRUE;
 201         }
 202         if (!prim)
 203           continue;
 204
 205         if (pdgPart->Charge()==0)
 206           continue;
 207
 208         Float_t eta = part->Eta();
 209         Float_t pt  = part->Pt();
 210
 211         if (prim) {
 212
 213           dNdEtaMap->FillParticleAllEvents(eta, pt);
 214
 215           if (triggered)
 216             dNdEtaMap->FillParticleWhenEventTriggered(eta, pt);
 217
 218           if (vertexReconstructed)
 219             dNdEtaMap->FillParticle(vtx[2], eta, 1.);
 220         }
 221
 222       }// end of mc particle
 223
 224       // ########################################################
 225       // loop over esd tracks
 226       Int_t nGoodTracks = 0;
 227
 228       Int_t nTracks = esd->GetNumberOfTracks();
 229       for (Int_t t=0; t<nTracks; t++) {
 230         AliESDtrack* esdTrack = esd->GetTrack(t);
 231
 232         // cut the esd track?
 233         if (!esdTrackCuts->AcceptTrack(esdTrack))
 234           continue;
 235
 236         nGoodTracks++;
 237
 238         Double_t p[3];
 239         esdTrack->GetPxPyPz(p);
 240         Float_t momentum = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2) + TMath::Power(p[2],2));
 241
 242         Float_t eta = -100.;
 243         if((momentum != TMath::Abs(p[2]))&&(momentum != 0))
 244           eta = 0.5*TMath::Log((momentum + p[2])/(momentum - p[2]));
 245
 246         // using the eta of the mc particle
 247         Int_t label = TMath::Abs(esdTrack->GetLabel());
 248         if (label<0) {
 249           cout << " cannot find corresponding mc part !!! " << label << endl;
 250           continue;
 251         }
 252         TParticle* mcPart = particleStack->Particle(label);
 253         eta = mcPart->Eta();
 254         Float_t pt = mcPart->Pt();
 255
 256         if (vertexReconstructed)
 257           dNdEtaMap->FillParticleWhenMeasuredTrack(vtx[2], eta, pt);
 258
 259       } // end of track loop
 260
 261       dNdEtaMap->FillEvent(vtxMC[2], nGoodTracks);
 262
 263       if (vertexReconstructed)
 264         dNdEtaMap->FillEventWithReconstructedVertex(vtxMC[2], nGoodTracks);
 265
 266     } // end  of event loop
 267   } // end of run loop
 268
 269   dNdEtaMap->Finish(nEventsAll, nEventsTriggered);
 270
 271   TFile* fout = new TFile("correction_map.root","RECREATE");
 272
 273   esdTrackCuts->SaveHistograms("esd_track_cuts");
 274   dNdEtaMap->SaveHistograms();
 275
 276   fout->Write();
 277   fout->Close();
 278
 279   dNdEtaMap->DrawHistograms();
 280
 281 }