]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG0/dNdEta/AlidNdEtaAnalysisMCSelector.cxx
EMCAL geometry can be created independently form anything now
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / AlidNdEtaAnalysisMCSelector.cxx
index 2e670eb6f0c098fe37bf5f08f3a947aca6b2cf05..4b57713ffbfbfac95bce7406ad9e1866565bb1bf 100644 (file)
 #include <AliGenEventHeader.h>
 #include <AliHeader.h>
 #include <AliStack.h>
+#include <AliESDtrack.h>
 
 #include "dNdEta/dNdEtaAnalysis.h"
-#include "AliPWG0Helper.h"
+#include <AliPWG0Helper.h>
+#include <AliCorrection.h>
+#include <AliCorrectionMatrix2D.h>
+#include "esdTrackCuts/AliESDtrackCuts.h"
 
 
 ClassImp(AlidNdEtaAnalysisMCSelector)
 
 AlidNdEtaAnalysisMCSelector::AlidNdEtaAnalysisMCSelector() :
   AliSelectorRL(),
+  fEsdTrackCuts(0),
   fdNdEtaAnalysis(0),
   fdNdEtaAnalysisTr(0),
   fdNdEtaAnalysisTrVtx(0),
+  fdNdEtaAnalysisTracks(0),
   fVertex(0),
-  fPartPt(0)
+  fPartPt(0),
+  fEvents(0)
 {
   //
   // Constructor. Initialization of pointers
@@ -44,6 +51,27 @@ AlidNdEtaAnalysisMCSelector::~AlidNdEtaAnalysisMCSelector()
   //
 }
 
+void AlidNdEtaAnalysisMCSelector::Begin(TTree* tree)
+{
+  // Begin function
+
+  ReadUserObjects(tree);
+}
+
+void AlidNdEtaAnalysisMCSelector::ReadUserObjects(TTree* tree)
+{
+  // read the user objects, called from slavebegin and begin
+
+  if (!fEsdTrackCuts && fInput)
+    fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fInput->FindObject("AliESDtrackCuts"));
+
+  if (!fEsdTrackCuts && tree)
+    fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (tree->GetUserInfo()->FindObject("AliESDtrackCuts"));
+
+  if (!fEsdTrackCuts)
+     AliDebug(AliLog::kError, "ERROR: Could not read EsdTrackCuts from input list.");
+}
+
 void AlidNdEtaAnalysisMCSelector::SlaveBegin(TTree * tree)
 {
   // The SlaveBegin() function is called after the Begin() function.
@@ -52,26 +80,39 @@ void AlidNdEtaAnalysisMCSelector::SlaveBegin(TTree * tree)
 
   AliSelectorRL::SlaveBegin(tree);
 
+  ReadUserObjects(tree);
+
   fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
   fdNdEtaAnalysisTr = new dNdEtaAnalysis("dndetaTr", "dndetaTr");
   fdNdEtaAnalysisTrVtx = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
+  fdNdEtaAnalysisTracks = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks");
   fVertex = new TH3F("vertex_check", "vertex_check", 50, -50, 50, 50, -50, 50, 50, -50, 50);
   for (Int_t i=0; i<3; ++i)
   {
-    fPartEta[i] = new TH1F(Form("dndeta_check_%d", i), Form("dndeta_check_%d", i), 120, -6, 6);
+    fPartEta[i] = new TH1F(Form("dndeta_check_%d", i), Form("dndeta_check_%d", i), 60, -6, 6);
     fPartEta[i]->Sumw2();
   }
   fPartPt =  new TH1F("dndeta_check_pt", "dndeta_check_pt", 1000, 0, 10);
   fPartPt->Sumw2();
+  fEvents = new TH1F("dndeta_check_vertex", "dndeta_check_vertex", 40, -20, 20);
 }
 
 void AlidNdEtaAnalysisMCSelector::Init(TTree *tree)
 {
   AliSelectorRL::Init(tree);
 
+  if (!tree)
+  {
+    AliDebug(AliLog::kError, "tree not available");
+    return;
+  }
+
   tree->SetBranchStatus("*", 0);
   tree->SetBranchStatus("fTriggerMask", 1);
   tree->SetBranchStatus("fSPDVertex*", 1);
+  tree->SetBranchStatus("fTracks.fLabel", 1);
+
+  AliESDtrackCuts::EnableNeededBranches(tree);
 }
 
 Bool_t AlidNdEtaAnalysisMCSelector::Process(Long64_t entry)
@@ -114,6 +155,8 @@ Bool_t AlidNdEtaAnalysisMCSelector::Process(Long64_t entry)
   // loop over mc particles
   Int_t nPrim  = stack->GetNprimary();
 
+  Int_t nAcceptedParticles = 0;
+
   for (Int_t iMc = 0; iMc < nPrim; ++iMc)
   {
     TParticle* particle = stack->Particle(iMc);
@@ -125,18 +168,19 @@ Bool_t AlidNdEtaAnalysisMCSelector::Process(Long64_t entry)
       continue;
 
     AliDebug(AliLog::kDebug+1, Form("Accepted primary %d, unique ID: %d", iMc, particle->GetUniqueID()));
+    ++nAcceptedParticles;
 
-    fdNdEtaAnalysis->FillTrack(vtxMC[2], particle->Eta(), particle->Pt(), 1);
+    fdNdEtaAnalysis->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
     fVertex->Fill(particle->Vx(), particle->Vy(), particle->Vz());
 
     if (eventTriggered)
     {
-      fdNdEtaAnalysisTr->FillTrack(vtxMC[2], particle->Eta(), particle->Pt(), 1);
+      fdNdEtaAnalysisTr->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
       if (vertexReconstructed)
-        fdNdEtaAnalysisTrVtx->FillTrack(vtxMC[2], particle->Eta(), particle->Pt(), 1);
+        fdNdEtaAnalysisTrVtx->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
     }
 
-    if (TMath::Abs(vtxMC[2]) < 10)
+    if (TMath::Abs(vtxMC[2]) < 20)
     {
       fPartEta[0]->Fill(particle->Eta());
 
@@ -149,14 +193,54 @@ Bool_t AlidNdEtaAnalysisMCSelector::Process(Long64_t entry)
     if (TMath::Abs(particle->Eta()) < 0.8)
       fPartPt->Fill(particle->Pt());
   }
-  fdNdEtaAnalysis->FillEvent(vtxMC[2], 1);
+
+  fEvents->Fill(vtxMC[2]);
+
+  fdNdEtaAnalysis->FillEvent(vtxMC[2], nAcceptedParticles);
   if (eventTriggered)
   {
-    fdNdEtaAnalysisTr->FillEvent(vtxMC[2], 1);
+    fdNdEtaAnalysisTr->FillEvent(vtxMC[2], nAcceptedParticles);
     if (vertexReconstructed)
-      fdNdEtaAnalysisTrVtx->FillEvent(vtxMC[2], 1);
+      fdNdEtaAnalysisTrVtx->FillEvent(vtxMC[2], nAcceptedParticles);
   }
 
+  if (!eventTriggered || !vertexReconstructed)
+    return kTRUE;
+
+  // from tracks is only done for triggered and vertex reconstructed events
+
+  // get number of "good" tracks
+  TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
+  Int_t nGoodTracks = list->GetEntries();
+
+  // loop over esd tracks
+  for (Int_t t=0; t<nGoodTracks; t++)
+  {
+    AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(t));
+    if (!esdTrack)
+    {
+      AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", t));
+      continue;
+    }
+
+    Int_t label = TMath::Abs(esdTrack->GetLabel());
+
+    TParticle* particle = stack->Particle(label);
+    if (!particle)
+    {
+      AliDebug(AliLog::kError, Form("ERROR: Could not retrieve particle %d.", esdTrack->GetLabel()));
+      continue;
+    }
+
+    fdNdEtaAnalysisTracks->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
+  } // end of track loop
+
+  delete list;
+  list = 0;
+
+  // for event count per vertex
+  fdNdEtaAnalysisTracks->FillEvent(vtxMC[2], nGoodTracks);
+
   return kTRUE;
 }
 
@@ -178,7 +262,10 @@ void AlidNdEtaAnalysisMCSelector::SlaveTerminate()
   fOutput->Add(fdNdEtaAnalysis);
   fOutput->Add(fdNdEtaAnalysisTr);
   fOutput->Add(fdNdEtaAnalysisTrVtx);
+  fOutput->Add(fdNdEtaAnalysisTracks);
+
   fOutput->Add(fPartPt);
+  fOutput->Add(fEvents);
   for (Int_t i=0; i<3; ++i)
     fOutput->Add(fPartEta[i]);
 }
@@ -192,29 +279,31 @@ void AlidNdEtaAnalysisMCSelector::Terminate()
   fdNdEtaAnalysis = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta"));
   fdNdEtaAnalysisTr = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTr"));
   fdNdEtaAnalysisTrVtx = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTrVtx"));
+  fdNdEtaAnalysisTracks = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTracks"));
   fPartPt = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_pt"));
+  fEvents = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_vertex"));
   for (Int_t i=0; i<3; ++i)
     fPartEta[i] = dynamic_cast<TH1F*> (fOutput->FindObject(Form("dndeta_check_%d", i)));
 
-  if (!fdNdEtaAnalysis || !fdNdEtaAnalysisTr || !fdNdEtaAnalysisTrVtx || !fPartPt)
+  if (!fdNdEtaAnalysis || !fdNdEtaAnalysisTr || !fdNdEtaAnalysisTrVtx || !fPartPt || !fEvents)
   {
     AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p %p", (void*) fdNdEtaAnalysis, (void*) fPartPt));
     return;
   }
 
-  fdNdEtaAnalysis->Finish(0, -1);
-  fdNdEtaAnalysisTr->Finish(0, -1);
-  fdNdEtaAnalysisTrVtx->Finish(0, -1);
+  fdNdEtaAnalysis->Finish(0, -1, AlidNdEtaCorrection::kNone);
+  fdNdEtaAnalysisTr->Finish(0, -1, AlidNdEtaCorrection::kNone);
+  fdNdEtaAnalysisTrVtx->Finish(0, -1, AlidNdEtaCorrection::kNone);
+  fdNdEtaAnalysisTracks->Finish(0, -1, AlidNdEtaCorrection::kNone);
 
-  Int_t events = (Int_t) fdNdEtaAnalysis->GetVtxHistogram()->Integral();
+  Int_t events = (Int_t) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Integral();
   fPartPt->Scale(1.0/events);
   fPartPt->Scale(1.0/fPartPt->GetBinWidth(1));
 
   if (fPartEta[0])
   {
-    TH1* vtxHist = fdNdEtaAnalysis->GetVtxHistogram();
-    Int_t events1 = (Int_t) vtxHist->Integral(vtxHist->GetXaxis()->FindBin(-9.9), vtxHist->GetXaxis()->FindBin(-0.1));
-    Int_t events2 = (Int_t) vtxHist->Integral(vtxHist->GetXaxis()->FindBin(0.1), vtxHist->GetXaxis()->FindBin(9.9));
+    Int_t events1 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(-19.9), fEvents->GetXaxis()->FindBin(-0.001));
+    Int_t events2 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(0.001), fEvents->GetXaxis()->FindBin(19.9));
 
     fPartEta[0]->Scale(1.0 / (events1 + events2));
     fPartEta[1]->Scale(1.0 / events1);
@@ -236,6 +325,7 @@ void AlidNdEtaAnalysisMCSelector::Terminate()
   fdNdEtaAnalysis->SaveHistograms();
   fdNdEtaAnalysisTr->SaveHistograms();
   fdNdEtaAnalysisTrVtx->SaveHistograms();
+  fdNdEtaAnalysisTracks->SaveHistograms();
   fPartPt->Write();
 
   fout->Write();
@@ -246,4 +336,10 @@ void AlidNdEtaAnalysisMCSelector::Terminate()
     new TCanvas("control2", "control2", 500, 500);
     fPartPt->DrawCopy();
   }
+
+  if (fEvents)
+  {
+    new TCanvas("control3", "control3", 500, 500);
+    fEvents->DrawCopy();
+  }
 }