PWG0/ptSpectra/makeCorrectionPtEta.C

   1 /* $Id$ */
   2
   3 //
   4 // Script to calculate PtvsEta correction map using the CorrectionMatrix2D class.
   5 //
   6
   7 makeCorrectionPtEta(Char_t* dataDir, Int_t nRuns=20) {
   8
   9   Char_t str[256];
  10
  11   //gSystem->Load("/home/poghos/CERN2006/pp/esdTrackCuts/libESDtrackQuality.so");
  12   gSystem->Load("libPWG0base");
  13   gSystem->Load("libCorrectionMatrix2D.so");
  14   // ########################################################
  15   // selection of esd tracks
  16   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts();
  17   esdTrackCuts->DefineHistograms(1);
  18   esdTrackCuts->SetMinNClustersTPC(50);
  19   esdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
  20   esdTrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
  21   esdTrackCuts->SetRequireTPCRefit(kTRUE);
  22   esdTrackCuts->SetMinNsigmaToVertex(3);
  23   esdTrackCuts->SetAcceptKingDaughters(kFALSE);
  24   AliLog::SetClassDebugLevel("AliESDtrackCuts",1);
  25   AliLog::SetClassDebugLevel("ESDtrackQualityCuts",1);
  26
  27   // ########################################################
  28   // definition of PtEta correction object
  29   CorrectionMatrix2D* PtEtaMap = new CorrectionMatrix2D("CorrectionMatrix2");
  30   PtEtaMap->SetHist("PtEta",80,0.,10.,120,-2.,2.);
  31   //Float_t x[]={-20., -15., -5., 0., 6., 20.};
  32   //Float_t y[]={-6. , -3.,  -1., 2., 3., 6. };
  33   //PtEtaMap->SetHist("PtEta",5,x,5,y);
  34   PtEtaMap->SetHistTitle("p_{t}","#eta");
  35
  36
  37   // ########################################################
  38   // get the data dir
  39   Char_t execDir[256];
  40   sprintf(execDir,gSystem->pwd());
  41   TSystemDirectory* baseDir = new TSystemDirectory(".",dataDir);
  42   TList* dirList            = baseDir->GetListOfFiles();
  43   Int_t nDirs               = dirList->GetEntries();
  44   // go back to the dir where this script is executed
  45   gSystem->cd(execDir);
  46
  47   // ########################################################
  48   // definition of used pointers
  49   TFile* esdFile;
  50   TTree* esdTree;
  51   TBranch* esdBranch;
  52
  53   AliESD* esd =0;
  54   // ########################################################
  55   // loop over runs
  56   Int_t nRunCounter = 0;
  57   for (Int_t r=0; r<nDirs; r++) {
  58
  59     TSystemFile* presentDir = (TSystemFile*)dirList->At(r);
  60     if (!presentDir || !presentDir->IsDirectory())
  61       continue;
  62     // check that the files are there
  63     TString currentDataDir;
  64     currentDataDir.Form("%s/%s",dataDir, presentDir->GetName());
  65     cout << "Processing directory " << currentDataDir.Data() << endl;
  66     if ((!gSystem->Which(currentDataDir,"galice.root")) ||
  67           (!gSystem->Which(currentDataDir,"AliESDs.root")))
  68       continue;
  69
  70     if (nRunCounter++ >= nRuns)
  71       break;
  72
  73     cout << "run #" << nRunCounter << endl;
  74
  75     // #########################################################
  76     // setup galice and runloader
  77     if (gAlice) {
  78       delete gAlice->GetRunLoader();
  79       delete gAlice;
  80       gAlice=0;
  81     }
  82
  83     sprintf(str,"%s/galice.root",currentDataDir.Data());
  84     AliRunLoader* runLoader = AliRunLoader::Open(str);
  85
  86     runLoader->LoadgAlice();
  87     gAlice = runLoader->GetAliRun();
  88     runLoader->LoadKinematics();
  89     runLoader->LoadHeader();
  90
  91
  92     // #########################################################
  93     // open esd file and get the tree
  94
  95     sprintf(str,"%s/AliESDs.root",currentDataDir.Data());
  96     // close it first to avoid memory leak
  97     if (esdFile)
  98       if (esdFile->IsOpen())
  99         esdFile->Close();
 100
 101     esdFile = TFile::Open(str);
 102     esdTree = (TTree*)esdFile->Get("esdTree");
 103     if (!esdTree)
 104       continue;
 105     esdBranch = esdTree->GetBranch("ESD");
 106     esdBranch->SetAddress(&esd);
 107     if (!esdBranch)
 108       continue;
 109
 110     // ########################################################
 111     // Magnetic field
 112     AliTracker::SetFieldMap(gAlice->Field(),kTRUE); // kTRUE means uniform magnetic field
 113     //AliKalmanTrack::SetConvConst(1000/0.299792458/5.);
 114
 115     // ########################################################
 116     // getting number of events
 117     Int_t nEvents = (Int_t)runLoader->GetNumberOfEvents();
 118     Int_t nESDEvents = esdBranch->GetEntries();
 119
 120     if (nEvents!=nESDEvents) {
 121       cout << " Different number of events from runloader and esdtree!!!" << nEvents << " / " << nESDEvents << endl;
 122       return;
 123     }
 124     // ########################################################
 125     // loop over number of events
 126     cout << " looping over events..." << endl;
 127     for(Int_t i=0; i<nEvents; i++) {
 128
 129       esdBranch->GetEntry(i);
 130       runLoader->GetEvent(i);
 131
 132       // ########################################################
 133       // get the EDS vertex
 134       AliESDVertex* vtxESD = esd->GetVertex();
 135
 136       Double_t vtx[3];
 137       Double_t vtx_res[3];
 138       vtxESD->GetXYZ(vtx);
 139
 140       vtx_res[0] = vtxESD->GetXRes();
 141       vtx_res[1] = vtxESD->GetYRes();
 142       vtx_res[2] = vtxESD->GetZRes();
 143
 144       // the vertex should be reconstructed
 145       if (strcmp(vtxESD->GetName(),"default")==0)
 146         continue;
 147
 148       // the resolution should be reasonable???
 149       if (vtx_res[2]==0 || vtx_res[2]>0.1)
 150         continue;
 151
 152       // ########################################################
 153       // loop over mc particles
 154       AliStack* particleStack = runLoader->Stack();
 155       Int_t nPrim    = particleStack->GetNprimary();
 156
 157       for (Int_t i_mc=0; i_mc<nPrim; i_mc++) {
 158
 159         TParticle* part = particleStack->Particle(i_mc);
 160         if (!part || strcmp(part->GetName(),"XXX")==0)
 161           continue;
 162
 163         TParticlePDG* pdgPart = part->GetPDG();
 164
 165         Bool_t prim = kFALSE;
 166         // check if it's a primary particle - is there a better way ???
 167         if ((part->GetFirstDaughter() >= nPrim) || (part->GetFirstDaughter()==-1)) {
 168           if (TMath::Abs(pdgPart->PdgCode())>10 && pdgPart->PdgCode()!=21 && strcmp(pdgPart->ParticleClass(),"Unknown")!=0)
 169             prim = kTRUE;
 170         }
 171         if (!prim)
 172           continue;
 173
 174         if (pdgPart->Charge()==0)
 175           continue;
 176
 177         if (prim)
 178           PtEtaMap->FillGene(part->Pt(), part->Eta());
 179
 180       }// end of mc particle
 181
 182       // ########################################################
 183       // loop over esd tracks
 184       Int_t nTracks = esd->GetNumberOfTracks();
 185       for (Int_t t=0; t<nTracks; t++) {
 186         AliESDtrack* esdTrack = esd->GetTrack(t);
 187
 188         // cut the esd track?
 189         if (!esdTrackCuts->AcceptTrack(esdTrack))
 190           continue;
 191
 192       Double_t Ptr0[3];
 193       if(!esdTrack->GetPxPyPz(Ptr0)) continue;
 194       fTrP= new TVector3(Ptr0);
 195
 196         PtEtaMap->FillMeas(fTrP->Pt(), fTrP->Eta());
 197
 198       } // end of track loop
 199     } // end  of event loop
 200   } // end of run loop
 201
 202   PtEtaMap->Finish();
 203
 204   TFile* fout = new TFile("PtEtaCorrectionMap.root","RECREATE");
 205
 206   esdTrackCuts->SaveHistograms("EsdTrackCuts");
 207   PtEtaMap->SaveHistograms();
 208
 209   fout->Write();
 210   fout->Close();
 211
 212 }