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