From: tpogos <tpogos@f7af4fe6-9843-0410-8265-dc069ae4e863>
Date: Mon, 22 May 2006 09:10:18 +0000 (+0000)
Subject: (martin) pt vs eta correction matrix calculation macro using CorrectionMatrix2D class.
X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=commitdiff_plain;h=559d0c87ae779c3681ae7486db74297306ad8138

(martin) pt vs eta correction matrix calculation macro using CorrectionMatrix2D class.
---

diff --git a/PWG0/ptSpectra/makeCorrectionPtEta.C b/PWG0/ptSpectra/makeCorrectionPtEta.C
new file mode 100644
index 00000000000..644c73c9da1
--- /dev/null
+++ b/PWG0/ptSpectra/makeCorrectionPtEta.C
@@ -0,0 +1,210 @@
+//
+// Script to calculate PtvsEta correction map using the CorrectionMatrix2D class.
+// 
+
+makeCorrectionPtEta(Char_t* dataDir, Int_t nRuns=20) {
+
+  Char_t str[256];
+
+  //gSystem->Load("/home/poghos/CERN2006/pp/esdTrackCuts/libESDtrackQuality.so");
+  gSystem->Load("libPWG0base");
+  gSystem->Load("libCorrectionMatrix2D.so");
+  // ########################################################
+  // selection of esd tracks
+  AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts();    
+  esdTrackCuts->DefineHistograms(1);
+  esdTrackCuts->SetMinNClustersTPC(50);
+  esdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
+  esdTrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
+  esdTrackCuts->SetRequireTPCRefit(kTRUE);
+  esdTrackCuts->SetMinNsigmaToVertex(3);
+  esdTrackCuts->SetAcceptKingDaughters(kFALSE);
+  AliLog::SetClassDebugLevel("AliESDtrackCuts",1);
+  AliLog::SetClassDebugLevel("ESDtrackQualityCuts",1);
+
+  // ########################################################
+  // definition of PtEta correction object
+  CorrectionMatrix2D* PtEtaMap = new CorrectionMatrix2D("CorrectionMatrix2");
+  PtEtaMap->SetHist("PtEta",80,0.,10.,120,-2.,2.);
+  //Float_t x[]={-20., -15., -5., 0., 6., 20.};
+  //Float_t y[]={-6. , -3.,  -1., 2., 3., 6. };
+  //PtEtaMap->SetHist("PtEta",5,x,5,y);
+  PtEtaMap->SetHistTitle("p_{t}","#eta");
+
+
+  // ########################################################
+  // get the data dir  
+  Char_t execDir[256];
+  sprintf(execDir,gSystem->pwd());
+  TSystemDirectory* baseDir = new TSystemDirectory(".",dataDir);
+  TList* dirList            = baseDir->GetListOfFiles();
+  Int_t nDirs               = dirList->GetEntries();
+  // go back to the dir where this script is executed
+  gSystem->cd(execDir);
+
+  // ########################################################
+  // definition of used pointers
+  TFile* esdFile;
+  TTree* esdTree;
+  TBranch* esdBranch;
+
+  AliESD* esd =0;
+  // ########################################################
+  // loop over runs
+  Int_t nRunCounter = 0;
+  for (Int_t r=0; r<nDirs; r++) {
+
+    TSystemFile* presentDir = (TSystemFile*)dirList->At(r);
+    if (!presentDir || !presentDir->IsDirectory())
+      continue;
+    // check that the files are there
+    TString currentDataDir;
+    currentDataDir.Form("%s/%s",dataDir, presentDir->GetName());
+    cout << "Processing directory " << currentDataDir.Data() << endl;
+    if ((!gSystem->Which(currentDataDir,"galice.root")) ||
+          (!gSystem->Which(currentDataDir,"AliESDs.root")))
+      continue;
+
+    if (nRunCounter++ >= nRuns)
+      break;
+
+    cout << "run #" << nRunCounter << endl;
+
+    // #########################################################
+    // setup galice and runloader
+    if (gAlice) {
+      delete gAlice->GetRunLoader();
+      delete gAlice;
+      gAlice=0;
+    }
+
+    sprintf(str,"%s/galice.root",currentDataDir.Data());
+    AliRunLoader* runLoader = AliRunLoader::Open(str);
+
+    runLoader->LoadgAlice();
+    gAlice = runLoader->GetAliRun();
+    runLoader->LoadKinematics();
+    runLoader->LoadHeader();
+
+   
+    // #########################################################
+    // open esd file and get the tree
+
+    sprintf(str,"%s/AliESDs.root",currentDataDir.Data());
+    // close it first to avoid memory leak
+    if (esdFile)
+      if (esdFile->IsOpen())
+        esdFile->Close();
+
+    esdFile = TFile::Open(str);
+    esdTree = (TTree*)esdFile->Get("esdTree");
+    if (!esdTree)
+      continue;
+    esdBranch = esdTree->GetBranch("ESD");
+    esdBranch->SetAddress(&esd);
+    if (!esdBranch)
+      continue;
+
+    // ########################################################
+    // Magnetic field
+    AliTracker::SetFieldMap(gAlice->Field(),kTRUE); // kTRUE means uniform magnetic field
+    //AliKalmanTrack::SetConvConst(1000/0.299792458/5.);
+
+    // ########################################################
+    // getting number of events
+    Int_t nEvents = (Int_t)runLoader->GetNumberOfEvents();
+    Int_t nESDEvents = esdBranch->GetEntries();
+    
+    if (nEvents!=nESDEvents) {
+      cout << " Different number of events from runloader and esdtree!!!" << nEvents << " / " << nESDEvents << endl;
+      return;
+    }
+    // ########################################################
+    // loop over number of events
+    cout << " looping over events..." << endl;
+    for(Int_t i=0; i<nEvents; i++) {
+
+      esdBranch->GetEntry(i);
+      runLoader->GetEvent(i);              
+
+      // ########################################################
+      // get the EDS vertex
+      AliESDVertex* vtxESD = esd->GetVertex();
+
+      Double_t vtx[3];
+      Double_t vtx_res[3];
+      vtxESD->GetXYZ(vtx);          
+    
+      vtx_res[0] = vtxESD->GetXRes();
+      vtx_res[1] = vtxESD->GetYRes();
+      vtx_res[2] = vtxESD->GetZRes();
+
+      // the vertex should be reconstructed
+      if (strcmp(vtxESD->GetName(),"default")==0) 
+	continue;
+
+      // the resolution should be reasonable???
+      if (vtx_res[2]==0 || vtx_res[2]>0.1)
+	continue;
+
+      // ########################################################
+      // loop over mc particles
+      AliStack* particleStack = runLoader->Stack();
+      Int_t nPrim    = particleStack->GetNprimary();
+
+      for (Int_t i_mc=0; i_mc<nPrim; i_mc++) {
+      
+	TParticle* part = particleStack->Particle(i_mc);      
+	if (!part || strcmp(part->GetName(),"XXX")==0) 
+	  continue;
+      
+	TParticlePDG* pdgPart = part->GetPDG();
+
+	Bool_t prim = kFALSE;
+	// check if it's a primary particle - is there a better way ???
+	if ((part->GetFirstDaughter() >= nPrim) || (part->GetFirstDaughter()==-1)) {
+	  if (TMath::Abs(pdgPart->PdgCode())>10 && pdgPart->PdgCode()!=21 && strcmp(pdgPart->ParticleClass(),"Unknown")!=0)
+	    prim = kTRUE;
+	}
+	if (!prim)
+	  continue;
+
+	if (pdgPart->Charge()==0)
+	  continue;
+
+	if (prim)
+	  PtEtaMap->FillGene(part->Pt(), part->Eta());	
+	
+      }// end of mc particle
+
+      // ########################################################
+      // loop over esd tracks      
+      Int_t nTracks = esd->GetNumberOfTracks();            
+      for (Int_t t=0; t<nTracks; t++) {
+	AliESDtrack* esdTrack = esd->GetTrack(t);      
+	
+	// cut the esd track?
+	if (!esdTrackCuts->AcceptTrack(esdTrack))
+	  continue;
+	
+      Double_t Ptr0[3];
+      if(!esdTrack->GetPxPyPz(Ptr0)) continue;
+      fTrP= new TVector3(Ptr0);
+	
+	PtEtaMap->FillMeas(fTrP->Pt(), fTrP->Eta());	
+
+      } // end of track loop
+    } // end  of event loop
+  } // end of run loop
+
+  PtEtaMap->Finish();  
+
+  TFile* fout = new TFile("PtEtaCorrectionMap.root","RECREATE");
+  
+  esdTrackCuts->SaveHistograms("EsdTrackCuts");
+  PtEtaMap->SaveHistograms();
+  
+  fout->Write();
+  fout->Close();
+
+}