PWG0/ptSpectra/makeCorrectionPtEtaVSz.C

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