adapting dndeta analysis to new AliCorrection schema
authorjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 1 Dec 2006 15:26:52 +0000 (15:26 +0000)
committerjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 1 Dec 2006 15:26:52 +0000 (15:26 +0000)
introducing vertexreco and trigger correction on track-level
adapting the corresponding selectors

12 files changed:
PWG0/dNdEta/AlidNdEtaAnalysisESDSelector.cxx
PWG0/dNdEta/AlidNdEtaAnalysisESDSelector.h
PWG0/dNdEta/AlidNdEtaAnalysisMCSelector.cxx
PWG0/dNdEta/AlidNdEtaAnalysisMCSelector.h
PWG0/dNdEta/AlidNdEtaCorrection.cxx
PWG0/dNdEta/AlidNdEtaCorrection.h
PWG0/dNdEta/AlidNdEtaCorrectionSelector.cxx
PWG0/dNdEta/AlidNdEtaCorrectionSelector.h
PWG0/dNdEta/dNdEtaAnalysis.cxx
PWG0/dNdEta/dNdEtaAnalysis.h
PWG0/dNdEta/drawPlots.C
PWG0/dNdEta/testAnalysis2.C

index 1ac813852175d17037faa8eb2df83bdbd56d7963..6f5b97e133a09b16c1b8522b9e7edb93051ae0ad 100644 (file)
 
 #include "esdTrackCuts/AliESDtrackCuts.h"
 #include "dNdEta/dNdEtaAnalysis.h"
-#include "dNdEta/AlidNdEtaCorrection.h"
 #include "AliPWG0Helper.h"
 
+#include "AliGenEventHeader.h"
+#include "AliHeader.h"
+#include "AliStack.h"
+#include "TParticle.h"
+
 ClassImp(AlidNdEtaAnalysisESDSelector)
 
 AlidNdEtaAnalysisESDSelector::AlidNdEtaAnalysisESDSelector() :
-  AliSelector(),
-  fdNdEtaAnalysisMBVtx(0),
-  fdNdEtaAnalysisMB(0),
+  AliSelectorRL(),
   fdNdEtaAnalysis(0),
-  fEsdTrackCuts(0),
-  fdNdEtaCorrection(0)
-  {
+  fEsdTrackCuts(0)
+{
   //
   // Constructor. Initialization of pointers
   //
@@ -64,16 +65,6 @@ void AlidNdEtaAnalysisESDSelector::ReadUserObjects(TTree* tree)
 
   if (!fEsdTrackCuts)
      AliDebug(AliLog::kError, "ERROR: Could not read EsdTrackCuts from input list.");
-
-
-  if (!fdNdEtaCorrection && fInput)
-    fdNdEtaCorrection = dynamic_cast<AlidNdEtaCorrection*> (fInput->FindObject("dndeta_correction"));
-
-  if (!fdNdEtaCorrection && tree)
-    fdNdEtaCorrection = dynamic_cast<AlidNdEtaCorrection*> (tree->GetUserInfo()->FindObject("dndeta_correction"));
-
-  if (!fdNdEtaCorrection)
-     AliDebug(AliLog::kError, "ERROR: Could not read dndeta_correction from input list.");
 }
 
 void AlidNdEtaAnalysisESDSelector::SlaveBegin(TTree* tree)
@@ -82,12 +73,10 @@ void AlidNdEtaAnalysisESDSelector::SlaveBegin(TTree* tree)
   // When running with PROOF SlaveBegin() is called on each slave server.
   // The tree argument is deprecated (on PROOF 0 is passed).
 
-  AliSelector::SlaveBegin(tree);
+  AliSelectorRL::SlaveBegin(tree);
 
   ReadUserObjects(tree);
 
-  fdNdEtaAnalysisMBVtx = new dNdEtaAnalysis("dndeta_mbvtx", "dndeta_mbvtx");
-  fdNdEtaAnalysisMB = new dNdEtaAnalysis("dndeta_mb", "dndeta_mb");
   fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
 }
 
@@ -95,7 +84,7 @@ void AlidNdEtaAnalysisESDSelector::Init(TTree* tree)
 {
   // read the user objects
 
-  AliSelector::Init(tree);
+  AliSelectorRL::Init(tree);
 
   // Enable only the needed branches
   if (tree)
@@ -103,6 +92,7 @@ void AlidNdEtaAnalysisESDSelector::Init(TTree* tree)
     tree->SetBranchStatus("*", 0);
     tree->SetBranchStatus("fTriggerMask", 1);
     tree->SetBranchStatus("fSPDVertex*", 1);
+    tree->SetBranchStatus("fTracks.fLabel", 1);
 
     AliESDtrackCuts::EnableNeededBranches(tree);
   }
@@ -110,25 +100,9 @@ void AlidNdEtaAnalysisESDSelector::Init(TTree* tree)
 
 Bool_t AlidNdEtaAnalysisESDSelector::Process(Long64_t entry)
 {
-  // The Process() function is called for each entry in the tree (or possibly
-  // keyed object in the case of PROOF) to be processed. The entry argument
-  // specifies which entry in the currently loaded tree is to be processed.
-  // It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
-  // to read either all or the required parts of the data. When processing
-  // keyed objects with PROOF, the object is already loaded and is available
-  // via the fObject pointer.
-  //
-  // This function should contain the "body" of the analysis. It can contain
-  // simple or elaborate selection criteria, run algorithms on the data
-  // of the event and typically fill histograms.
-
-  // WARNING when a selector is used with a TChain, you must use
-  //  the pointer to the current TTree to call GetEntry(entry).
-  //  The entry is always the local entry number in the current tree.
-  //  Assuming that fTree is the pointer to the TChain being processed,
-  //  use fTree->GetTree()->GetEntry(entry).
+  // loop over all events
 
-  if (AliSelector::Process(entry) == kFALSE)
+  if (AliSelectorRL::Process(entry) == kFALSE)
     return kFALSE;
 
   // Check prerequisites
@@ -144,12 +118,6 @@ Bool_t AlidNdEtaAnalysisESDSelector::Process(Long64_t entry)
     return kFALSE;
   }
 
-  if (!fdNdEtaCorrection)
-  {
-    AliDebug(AliLog::kError, "fdNdEtaCorrection not available");
-    return kFALSE;
-  }
-
   if (AliPWG0Helper::IsEventTriggered(fESD) == kFALSE)
   {
     AliDebug(AliLog::kDebug+1, Form("Skipping event %d because it was not triggered", (Int_t) entry));
@@ -162,31 +130,40 @@ Bool_t AlidNdEtaAnalysisESDSelector::Process(Long64_t entry)
     return kTRUE;
   }
 
+  AliHeader* header = GetHeader();
+  if (!header)
+  {
+    AliDebug(AliLog::kError, "Header not available");
+    return kFALSE;
+  }
+
+  AliStack* stack = GetStack();
+  if (!stack)
+  {
+    AliDebug(AliLog::kError, "Stack not available");
+    return kFALSE;
+  }
+
+  // get the MC vertex
+  AliGenEventHeader* genHeader = header->GenEventHeader();
+
+  TArrayF vtxMC(3);
+  genHeader->PrimaryVertex(vtxMC);
+
   // ########################################################
-  // get the EDS vertex
+  // get the ESD vertex
   const AliESDVertex* vtxESD = fESD->GetVertex();
   Double_t vtx[3];
   vtxESD->GetXYZ(vtx);
 
+  vtx[2] = vtxMC[2];
+
   // get number of "good" tracks
   TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
   Int_t nGoodTracks = list->GetEntries();
 
-  Float_t vertexRecoCorr = fdNdEtaCorrection->GetVertexRecoCorrection(vtx[2], nGoodTracks);
-  if (vertexRecoCorr <= 0)
-  {
-    AliDebug(AliLog::kError, Form("INFO: Skipping event because vertexRecoCorr is <= 0 (%f)", vertexRecoCorr));
-    delete list;
-    return kTRUE;
-  }
-
-  Float_t triggerCorr = fdNdEtaCorrection->GetTriggerBiasCorrection(vtx[2], nGoodTracks);
-  if (triggerCorr <= 0)
-  {
-    AliDebug(AliLog::kError, Form("INFO: Skipping event because triggerCorr is <= 0 (%f)", triggerCorr));
-    delete list;
-    return kTRUE;
-  }
+  // FAKE test!
+  //Int_t nContributors = vtxESD->GetNContributors();
 
   // loop over esd tracks
   for (Int_t t=0; t<nGoodTracks; t++)
@@ -198,6 +175,21 @@ Bool_t AlidNdEtaAnalysisESDSelector::Process(Long64_t entry)
       continue;
     }
 
+    Int_t label = TMath::Abs(esdTrack->GetLabel());
+
+    if (label == 0)
+    {
+      AliDebug(AliLog::kError, Form("Label is 0. Skipping! Track %d.", t));
+      continue;
+    }
+
+    TParticle* particle = stack->Particle(label);
+    if (!particle)
+    {
+      AliDebug(AliLog::kError, Form("ERROR: Could not retrieve particle %d.", esdTrack->GetLabel()));
+      continue;
+    }
+
     Double_t p[3];
     esdTrack->GetConstrainedPxPyPz(p); // ### TODO should be okay because we have a vertex, however GetInnerPxPyPy / GetOuterPxPyPy also exist
     TVector3 vector(p);
@@ -206,27 +198,17 @@ Bool_t AlidNdEtaAnalysisESDSelector::Process(Long64_t entry)
     Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
     Float_t pt = vector.Pt();
 
-    Float_t track2particleCorr = fdNdEtaCorrection->GetTrack2ParticleCorrection(vtx[2], eta, pt);
+    //eta = particle->Eta();
+    //pt = particle->Pt();
 
-    Float_t weight = track2particleCorr * vertexRecoCorr * triggerCorr;
-    if (weight <= 0)
-    {
-      AliDebug(AliLog::kError, Form("INFO: Skipping track because weight is <= 0 (track %d, weight %f) (vtx %f, eta %f, pt %f)", t, weight, vtx[2], eta, pt));
-      continue;
-    }
-
-    fdNdEtaAnalysisMBVtx->FillTrack(vtx[2], eta, pt, track2particleCorr);
-    fdNdEtaAnalysisMB->FillTrack(vtx[2], eta, pt, track2particleCorr * vertexRecoCorr);
-    fdNdEtaAnalysis->FillTrack(vtx[2], eta, pt, weight);
+    fdNdEtaAnalysis->FillTrack(vtx[2], eta, pt);
   } // end of track loop
 
   delete list;
   list = 0;
 
   // for event count per vertex
-  fdNdEtaAnalysisMBVtx->FillEvent(vtx[2], 1);
-  fdNdEtaAnalysisMB->FillEvent(vtx[2], vertexRecoCorr);
-  fdNdEtaAnalysis->FillEvent(vtx[2], vertexRecoCorr * triggerCorr);
+  fdNdEtaAnalysis->FillEvent(vtx[2], nGoodTracks);
 
   return kTRUE;
 }
@@ -237,7 +219,7 @@ void AlidNdEtaAnalysisESDSelector::SlaveTerminate()
   // have been processed. When running with PROOF SlaveTerminate() is called
   // on each slave server.
 
-  AliSelector::SlaveTerminate();
+  AliSelectorRL::SlaveTerminate();
 
   // Add the histograms to the output on each slave server
   if (!fOutput)
@@ -250,12 +232,6 @@ void AlidNdEtaAnalysisESDSelector::SlaveTerminate()
 
   fOutput->Add(fdNdEtaAnalysis);
   fdNdEtaAnalysis = 0;
-
-  fOutput->Add(fdNdEtaAnalysisMB);
-  fdNdEtaAnalysisMB = 0;
-
-  fOutput->Add(fdNdEtaAnalysisMBVtx);
-  fdNdEtaAnalysisMBVtx = 0;
 }
 
 void AlidNdEtaAnalysisESDSelector::Terminate()
@@ -264,44 +240,24 @@ void AlidNdEtaAnalysisESDSelector::Terminate()
   // a query. It always runs on the client, it can be used to present
   // the results graphically or save the results to file.
 
-  AliSelector::Terminate();
+  AliSelectorRL::Terminate();
 
   fdNdEtaAnalysis = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta"));
-  fdNdEtaAnalysisMB = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta_mb"));
-  fdNdEtaAnalysisMBVtx = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta_mbvtx"));
 
-  if (!fdNdEtaAnalysis || !fdNdEtaAnalysisMB || !fdNdEtaAnalysisMBVtx)
+  if (!fdNdEtaAnalysis)
   {
-    AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p %p %p", (void*) fdNdEtaAnalysis, (void*) fdNdEtaAnalysisMB, (void*) fdNdEtaAnalysisMBVtx));
+    AliDebug(AliLog::kError, "ERROR: Histograms not available");
     return;
   }
 
-  if (fdNdEtaAnalysis)
-    fdNdEtaAnalysis->Finish(fdNdEtaCorrection, 0.3);
-
-  if (fdNdEtaAnalysisMB)
-    fdNdEtaAnalysisMB->Finish(fdNdEtaCorrection, 0.3);
-
-  if (fdNdEtaAnalysisMBVtx)
-    fdNdEtaAnalysisMBVtx->Finish(fdNdEtaCorrection, 0.3);
-
-  TFile* fout = new TFile("analysis_esd.root","RECREATE");
+  TFile* fout = new TFile("analysis_esd_raw.root", "RECREATE");
 
   if (fdNdEtaAnalysis)
     fdNdEtaAnalysis->SaveHistograms();
 
-  if (fdNdEtaAnalysisMB)
-    fdNdEtaAnalysisMB->SaveHistograms();
-
-  if (fdNdEtaAnalysisMBVtx)
-    fdNdEtaAnalysisMBVtx->SaveHistograms();
-
   if (fEsdTrackCuts)
     fEsdTrackCuts->SaveHistograms("esd_tracks_cuts");
 
-  if (fdNdEtaCorrection)
-    fdNdEtaCorrection->SaveHistograms();
-
   fout->Write();
   fout->Close();
 }
index c716c05f625b86c8fafd01896bada9c59c209cd8..8ea364a8e9ecfd951d89594cfe6857f047e1b5b3 100644 (file)
@@ -3,13 +3,13 @@
 #ifndef ALIDNDETAANALYSISESDSELECTOR_H
 #define ALIDNDETAANALYSISESDSELECTOR_H
 
-#include "AliSelector.h"
+#include "AliSelectorRL.h"
 
 class AliESDtrackCuts;
 class dNdEtaAnalysis;
 class AlidNdEtaCorrection;
 
-class AlidNdEtaAnalysisESDSelector : public AliSelector {
+class AlidNdEtaAnalysisESDSelector : public AliSelectorRL {
   public:
     AlidNdEtaAnalysisESDSelector();
     virtual ~AlidNdEtaAnalysisESDSelector();
@@ -24,14 +24,10 @@ class AlidNdEtaAnalysisESDSelector : public AliSelector {
  protected:
     void ReadUserObjects(TTree* tree);
 
-    dNdEtaAnalysis* fdNdEtaAnalysisMBVtx;   // contains the histograms for the triggered events with vertex
-    dNdEtaAnalysis* fdNdEtaAnalysisMB;      // contains the histograms corrected with vtx recon eff
-    dNdEtaAnalysis* fdNdEtaAnalysis;        // contains the histograms corrected with vtx recon eff and trigger bias eff
+    dNdEtaAnalysis* fdNdEtaAnalysis;        // contains the uncorrected histograms
 
     AliESDtrackCuts*  fEsdTrackCuts;     // Object containing the parameters of the esd track cuts
 
-    AlidNdEtaCorrection* fdNdEtaCorrection; // correction maps
-
  private:
     AlidNdEtaAnalysisESDSelector(const AlidNdEtaAnalysisESDSelector&);
     AlidNdEtaAnalysisESDSelector& operator=(const AlidNdEtaAnalysisESDSelector&);
index 2e670eb6f0c098fe37bf5f08f3a947aca6b2cf05..1f34d1c890f04d563a30fd21e3a5ecee757b72e9 100644 (file)
@@ -19,7 +19,9 @@
 #include <AliStack.h>
 
 #include "dNdEta/dNdEtaAnalysis.h"
-#include "AliPWG0Helper.h"
+#include <AliPWG0Helper.h>
+#include <AliCorrection.h>
+#include <AliCorrectionMatrix2D.h>
 
 
 ClassImp(AlidNdEtaAnalysisMCSelector)
@@ -30,7 +32,8 @@ AlidNdEtaAnalysisMCSelector::AlidNdEtaAnalysisMCSelector() :
   fdNdEtaAnalysisTr(0),
   fdNdEtaAnalysisTrVtx(0),
   fVertex(0),
-  fPartPt(0)
+  fPartPt(0),
+  fEvents(0)
 {
   //
   // Constructor. Initialization of pointers
@@ -58,11 +61,12 @@ void AlidNdEtaAnalysisMCSelector::SlaveBegin(TTree * tree)
   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)
@@ -114,6 +118,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 +131,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,12 +156,15 @@ 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);
   }
 
   return kTRUE;
@@ -179,6 +189,7 @@ void AlidNdEtaAnalysisMCSelector::SlaveTerminate()
   fOutput->Add(fdNdEtaAnalysisTr);
   fOutput->Add(fdNdEtaAnalysisTrVtx);
   fOutput->Add(fPartPt);
+  fOutput->Add(fEvents);
   for (Int_t i=0; i<3; ++i)
     fOutput->Add(fPartEta[i]);
 }
@@ -193,10 +204,11 @@ void AlidNdEtaAnalysisMCSelector::Terminate()
   fdNdEtaAnalysisTr = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTr"));
   fdNdEtaAnalysisTrVtx = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTrVtx"));
   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;
@@ -206,15 +218,14 @@ void AlidNdEtaAnalysisMCSelector::Terminate()
   fdNdEtaAnalysisTr->Finish(0, -1);
   fdNdEtaAnalysisTrVtx->Finish(0, -1);
 
-  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);
@@ -246,4 +257,10 @@ void AlidNdEtaAnalysisMCSelector::Terminate()
     new TCanvas("control2", "control2", 500, 500);
     fPartPt->DrawCopy();
   }
+
+  if (fEvents)
+  {
+    new TCanvas("control3", "control3", 500, 500);
+    fEvents->DrawCopy();
+  }
 }
index 72a6e893509c9132e00283f6146f97d270e8eebd..58b10d3aa20247bafbd94f07af8fa734117fb793 100644 (file)
@@ -20,16 +20,16 @@ class AlidNdEtaAnalysisMCSelector : public AliSelectorRL {
     virtual Bool_t  Process(Long64_t entry);
     virtual void    Terminate();
 
- protected:
-
  private:
-    dNdEtaAnalysis* fdNdEtaAnalysis;      // contains the dndeta from the full sample
+    dNdEtaAnalysis* fdNdEtaAnalysis;        // contains the dndeta from the full sample
     dNdEtaAnalysis* fdNdEtaAnalysisTr;      // contains the dndeta from the triggered events
-    dNdEtaAnalysis* fdNdEtaAnalysisTrVtx;      // contains the dndeta from the triggered events with vertex
+    dNdEtaAnalysis* fdNdEtaAnalysisTrVtx;   // contains the dndeta from the triggered events with vertex
 
-    TH3F* fVertex;  //! vertex of counted particles
-    TH1F* fPartEta[3]; //! counted particles as function of eta (full vertex range, below 0 range, above 0 range)
-    TH1F* fPartPt; //! counted particles as function of pt
+    // the following are control histograms to check the dNdEtaAnalysis class
+    TH3F* fVertex;     // vertex of counted particles
+    TH1F* fPartEta[3]; // counted particles as function of eta (full vertex range, below 0 range, above 0 range)
+    TH1F* fPartPt;     // counted particles as function of pt
+    TH1F* fEvents;     // events counted as function of vtx
 
     AlidNdEtaAnalysisMCSelector(const AlidNdEtaAnalysisMCSelector&);
     AlidNdEtaAnalysisMCSelector& operator=(const AlidNdEtaAnalysisMCSelector&);
index ac466918f4ea7c246ee32fbf8548a87862d68952..248f34e3c2bae3f05249ec8abb0e4d573e917edf 100644 (file)
@@ -2,9 +2,13 @@
 
 #include "AlidNdEtaCorrection.h"
 
+#include <AliLog.h>
 #include <TCanvas.h>
 #include <TH3F.h>
 #include <TH1D.h>
+#include <AliCorrection.h>
+#include <AliCorrectionMatrix2D.h>
+#include <AliCorrectionMatrix3D.h>
 
 //____________________________________________________________________
 ClassImp(AlidNdEtaCorrection)
@@ -30,35 +34,16 @@ AlidNdEtaCorrection::AlidNdEtaCorrection(const Char_t* name, const Char_t* title
   fTriggerBiasCorrectionMBToNSD(0),
   fTriggerBiasCorrectionMBToND(0)
 {
+  //
   // constructor
   //
 
-  Float_t binLimitsPt[] = {0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.5, 2.0, 5.0, 10.0, 100.0};
-
-  TString matrixName;
-  matrixName.Form("%s_nTrackToNPart", name);
-
-  fTrack2ParticleCorrection = new AliCorrectionMatrix3D(matrixName, matrixName, 40, -20, 20, 20, -2, 2, 15, binLimitsPt);
+  fTrack2ParticleCorrection = new AliCorrection("Track2Particle", "Track2Particle");
+  fVertexRecoCorrection     = new AliCorrection("VertexReconstruction", "VertexReconstruction");
 
-  Float_t binLimitsN[]   = {-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
-                           10.5, 12.5, 14.5, 16.5, 18.5, 20.5, 25.5, 30.5, 40.5, 50.5, 100.5, 300.5};
-  Float_t binLimitsVtx[] = {-20,-15,-10,-6,-3,0,3,6,10,15,20};
-
-  matrixName.Form("%s_vtxReco", name);
-  fVertexRecoCorrection        = new AliCorrectionMatrix2D(matrixName, matrixName, 10,binLimitsVtx ,22,binLimitsN);
-
-  matrixName.Form("%s_triggerBias_MBToINEL", name);
-  fTriggerBiasCorrectionMBToINEL       = new AliCorrectionMatrix2D(matrixName, matrixName, 10,binLimitsVtx ,22,binLimitsN);
-  matrixName.Form("%s_triggerBias_MBToNSD", name);
-  fTriggerBiasCorrectionMBToNSD        = new AliCorrectionMatrix2D(matrixName, matrixName, 10,binLimitsVtx ,22,binLimitsN);
-  matrixName.Form("%s_triggerBias_MBToND", name);
-  fTriggerBiasCorrectionMBToND         = new AliCorrectionMatrix2D(matrixName, matrixName, 10,binLimitsVtx ,22,binLimitsN);
-  
-  fTrack2ParticleCorrection      ->SetAxisTitles("vtx z [cm]", "#eta", "p_{T} [GeV/c]");
-  fVertexRecoCorrection          ->SetAxisTitles("vtx z [cm]", "Ntracks");
-  fTriggerBiasCorrectionMBToINEL ->SetAxisTitles("vtx z [cm]", "Ntracks");
-  fTriggerBiasCorrectionMBToNSD  ->SetAxisTitles("vtx z [cm]", "Ntracks");
-  fTriggerBiasCorrectionMBToND   ->SetAxisTitles("vtx z [cm]", "Ntracks");
+  fTriggerBiasCorrectionMBToINEL = new AliCorrection("TriggerBias_MBToINEL", "TriggerBias_MBToINEL");
+  fTriggerBiasCorrectionMBToNSD  = new AliCorrection("TriggerBias_MBToNSD", "TriggerBias_MBToNSD");
+  fTriggerBiasCorrectionMBToND   = new AliCorrection("TriggerBias_MBToND", "TriggerBias_MBToND");
 }
 
 //____________________________________________________________________
@@ -101,20 +86,6 @@ AlidNdEtaCorrection::Finish() {
   // divide the histograms in the AliCorrectionMatrix2D objects to get the corrections
 
   fTrack2ParticleCorrection->Divide();
-
-  TH3F* hist = fTrack2ParticleCorrection->GetCorrectionHistogram();
-  Int_t emptyBins = 0;
-  for (Int_t x=hist->GetXaxis()->FindBin(-10); x<=hist->GetXaxis()->FindBin(10); ++x)
-    for (Int_t y=hist->GetYaxis()->FindBin(-0.8); y<=hist->GetYaxis()->FindBin(0.8); ++y)
-      for (Int_t z=hist->GetZaxis()->FindBin(0.3); z<=hist->GetZaxis()->FindBin(9.9); ++z)
-        if (hist->GetBinContent(x, y, z) == 0)
-        {
-          printf("Empty bin in fTrack2ParticleCorrection at vtx = %f, eta = %f, pt = %f\n", hist->GetXaxis()->GetBinCenter(x), hist->GetYaxis()->GetBinCenter(y), hist->GetZaxis()->GetBinCenter(z));
-          ++emptyBins;
-        }
-
-  printf("INFO: In the central region fTrack2ParticleCorrection has %d empty bins\n", emptyBins);
-
   fVertexRecoCorrection->Divide();
   fTriggerBiasCorrectionMBToINEL->Divide();
   fTriggerBiasCorrectionMBToNSD->Divide();
@@ -122,8 +93,8 @@ AlidNdEtaCorrection::Finish() {
 }
 
 //____________________________________________________________________
-Long64_t
-AlidNdEtaCorrection::Merge(TCollection* list) {
+Long64_t AlidNdEtaCorrection::Merge(TCollection* list)
+{
   // Merge a list of dNdEtaCorrection objects with this (needed for
   // PROOF).
   // Returns the number of merged objects (including this).
@@ -151,16 +122,13 @@ AlidNdEtaCorrection::Merge(TCollection* list) {
     if (entry == 0)
       continue;
 
-    collectionNtrackToNparticle  ->Add(entry->GetTrack2ParticleCorrection());
-    collectionVertexReco         ->Add(entry->GetVertexRecoCorrection());
-    collectionTriggerBiasMBToINEL->Add(entry->GetTriggerBiasCorrection("INEL"));
-    collectionTriggerBiasMBToNSD ->Add(entry->GetTriggerBiasCorrection("NSD"));
-    collectionTriggerBiasMBToND  ->Add(entry->GetTriggerBiasCorrection("ND"));
+    collectionNtrackToNparticle  ->Add(entry->fTrack2ParticleCorrection);
+    collectionVertexReco         ->Add(entry->fVertexRecoCorrection);
+    collectionTriggerBiasMBToINEL->Add(entry->fTriggerBiasCorrectionMBToINEL);
+    collectionTriggerBiasMBToNSD ->Add(entry->fTriggerBiasCorrectionMBToNSD);
+    collectionTriggerBiasMBToND  ->Add(entry->fTriggerBiasCorrectionMBToND);
 
     count++;
-
-    //fNEvents += entry->fNEvents;
-    //fNTriggeredEvents += entry->fNTriggeredEvents;
   }
   fTrack2ParticleCorrection      ->Merge(collectionNtrackToNparticle);
   fVertexRecoCorrection          ->Merge(collectionVertexReco);
@@ -177,29 +145,34 @@ AlidNdEtaCorrection::Merge(TCollection* list) {
   return count+1;
 }
 
-
-
 //____________________________________________________________________
-Bool_t
-AlidNdEtaCorrection::LoadHistograms(const Char_t* fileName, const Char_t* dir) {
+Bool_t AlidNdEtaCorrection::LoadHistograms(const Char_t* dir)
+{
   //
   // loads the histograms
+  // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
   //
 
-  fTrack2ParticleCorrection      ->LoadHistograms(fileName, dir);
-  fVertexRecoCorrection          ->LoadHistograms(fileName, dir);
-  fTriggerBiasCorrectionMBToINEL ->LoadHistograms(fileName, dir);
-  fTriggerBiasCorrectionMBToNSD  ->LoadHistograms(fileName, dir);
-  fTriggerBiasCorrectionMBToND   ->LoadHistograms(fileName, dir);
+  if (!dir)
+    dir = GetName();
 
-  
+  if (!gDirectory->cd(dir))
+    return kFALSE;
+
+  fTrack2ParticleCorrection      ->LoadHistograms();
+  fVertexRecoCorrection          ->LoadHistograms();
+  fTriggerBiasCorrectionMBToINEL ->LoadHistograms();
+  fTriggerBiasCorrectionMBToNSD  ->LoadHistograms();
+  fTriggerBiasCorrectionMBToND   ->LoadHistograms();
+
+  gDirectory->cd("..");
 
   return kTRUE;
 }
 
 //____________________________________________________________________
-void
-AlidNdEtaCorrection::SaveHistograms() {
+void AlidNdEtaCorrection::SaveHistograms()
+{
   //
   // save the histograms
   //
@@ -213,96 +186,107 @@ AlidNdEtaCorrection::SaveHistograms() {
   fTriggerBiasCorrectionMBToNSD ->SaveHistograms();
   fTriggerBiasCorrectionMBToND  ->SaveHistograms();
 
-  gDirectory->cd("../");
+  gDirectory->cd("..");
 }
 
 //____________________________________________________________________
 void AlidNdEtaCorrection::DrawHistograms()
 {
   //
-  // call the draw histogram method of the two AliCorrectionMatrix2D objects
+  // call the draw histogram method of the correction
+  //
 
   fTrack2ParticleCorrection     ->DrawHistograms();
   fVertexRecoCorrection         ->DrawHistograms();
   fTriggerBiasCorrectionMBToINEL->DrawHistograms();
   fTriggerBiasCorrectionMBToNSD ->DrawHistograms();
   fTriggerBiasCorrectionMBToND  ->DrawHistograms();
-
 }
 
 //____________________________________________________________________
-void AlidNdEtaCorrection::FillEventWithTrigger(Float_t vtx, Float_t n) 
+void AlidNdEtaCorrection::FillMCParticle(Float_t vtx, Float_t eta, Float_t pt, Bool_t trigger, Bool_t vertex, Int_t processType)
 {
-  // fill events with trigger.
-  // used to calculate vertex reco and trigger bias corrections
+  // fills a particle in the corrections
+  // it is filled in generated or measured depending of the flags
 
-  fVertexRecoCorrection->FillGene(vtx, n);
-  fTriggerBiasCorrectionMBToINEL ->FillMeas(vtx, n);
-  fTriggerBiasCorrectionMBToNSD  ->FillMeas(vtx, n);
-  fTriggerBiasCorrectionMBToND   ->FillMeas(vtx, n);  
-}
+  fTriggerBiasCorrectionMBToINEL->GetTrackCorrection()->FillGene(vtx, eta, pt);
 
-//____________________________________________________________________
-void AlidNdEtaCorrection::FillEventAll(Float_t vtx, Float_t n, Char_t* opt)  
-{
-  // fill all events 
-  // used to calculate trigger bias corrections
-
-  if (strcmp(opt,"INEL")==0) 
-    fTriggerBiasCorrectionMBToINEL->FillGene(vtx, n);
-  else if (strcmp(opt,"NSD")==0)  
-    fTriggerBiasCorrectionMBToNSD->FillGene(vtx, n);
-  else if (strcmp(opt,"ND")==0)   
-    fTriggerBiasCorrectionMBToND->FillGene(vtx, n);
-  else 
-    AliDebug(AliLog::kWarning, Form(" event type %s unknown (use INEL, NSD or ND)",opt));
+  if (processType != 92 && processType != 93)
+    fTriggerBiasCorrectionMBToNSD->GetTrackCorrection()->FillGene(vtx, eta, pt);
+
+  if (processType!=92 && processType!=93 && processType!=94)
+    fTriggerBiasCorrectionMBToND->GetTrackCorrection()->FillGene(vtx, eta, pt);
+
+  if (!trigger)
+    return;
+
+  fTriggerBiasCorrectionMBToINEL->GetTrackCorrection()->FillMeas(vtx, eta, pt);
+  fTriggerBiasCorrectionMBToNSD->GetTrackCorrection()->FillMeas(vtx, eta, pt);
+  fTriggerBiasCorrectionMBToND->GetTrackCorrection()->FillMeas(vtx, eta, pt);
+  fVertexRecoCorrection->GetTrackCorrection()->FillGene(vtx, eta, pt);
+
+  if (!vertex)
+    return;
+
+  fVertexRecoCorrection->GetTrackCorrection()->FillMeas(vtx, eta, pt);
+  fTrack2ParticleCorrection->GetTrackCorrection()->FillGene(vtx, eta, pt);
 }
 
 //____________________________________________________________________
-AliCorrectionMatrix2D* AlidNdEtaCorrection::GetTriggerBiasCorrection(Char_t* opt)  
+void AlidNdEtaCorrection::FillTrackedParticle(Float_t vtx, Float_t eta, Float_t pt)
 {
-  // returning the trigger bias correction matrix 
-  // option can be used to specify to which collision process (INEL, NSD or ND)  
-  if (strcmp(opt,"INEL")==0) 
-    return fTriggerBiasCorrectionMBToINEL;
-  else if (strcmp(opt,"NSD")==0)  
-    return fTriggerBiasCorrectionMBToNSD;
-  else if (strcmp(opt,"ND")==0)   
-    return fTriggerBiasCorrectionMBToND;
-  else 
-    {
-      AliDebug(AliLog::kWarning, Form(" %s is unknown (use INEL, NSD or ND). returning INEL ",opt)); 
-      return fTriggerBiasCorrectionMBToINEL;
-    }
+  // fills a tracked particle in the corrections
+
+  fTrack2ParticleCorrection->GetTrackCorrection()->FillMeas(vtx, eta, pt);
 }
 
 //____________________________________________________________________
-Float_t AlidNdEtaCorrection::GetTriggerBiasCorrection(Float_t vtx, Float_t n, Char_t* opt) 
+void AlidNdEtaCorrection::FillEvent(Float_t vtx, Float_t n, Bool_t trigger, Bool_t vertex, Int_t processType)
 {
-  // returning the trigger bias correction matrix 
-  // 3rd option can be used to specify to which collision process (INEL, NSD or ND)  
-
-  if (strcmp(opt,"INEL")==0) 
-    return fTriggerBiasCorrectionMBToINEL->GetCorrection(vtx, n);
-  else if (strcmp(opt,"NSD")==0)  
-    return fTriggerBiasCorrectionMBToNSD->GetCorrection(vtx, n);
-  else if (strcmp(opt,"ND")==0)   
-    return fTriggerBiasCorrectionMBToND->GetCorrection(vtx, n);
-  else {
-    AliDebug(AliLog::kWarning, Form(" %s unknown (use INEL, NSD or ND). returning corr for INEL",opt)); 
-    return fTriggerBiasCorrectionMBToINEL->GetCorrection(vtx, n);
-  }
-}
+  // fills an event int he correction
+  // it is filled in generated or measured depending of the flags
+
+  fTriggerBiasCorrectionMBToINEL->GetEventCorrection()->FillGene(vtx, n);
+
+  if (processType != 92 && processType != 93)
+    fTriggerBiasCorrectionMBToNSD->GetEventCorrection()->FillGene(vtx, n);
 
+  if (processType!=92 && processType!=93 && processType!=94)
+    fTriggerBiasCorrectionMBToND->GetEventCorrection()->FillGene(vtx, n);
+
+  if (!trigger)
+    return;
+
+  fTriggerBiasCorrectionMBToINEL->GetEventCorrection()->FillMeas(vtx, n);
+  fTriggerBiasCorrectionMBToNSD->GetEventCorrection()->FillMeas(vtx, n);
+  fTriggerBiasCorrectionMBToND->GetEventCorrection()->FillMeas(vtx, n);
+  fVertexRecoCorrection->GetEventCorrection()->FillGene(vtx, n);
+
+  if (!vertex)
+    return;
+
+  fVertexRecoCorrection->GetEventCorrection()->FillMeas(vtx, n);
+}
 
 //____________________________________________________________________
-Float_t AlidNdEtaCorrection::GetMeasuredFraction(Float_t ptCutOff, Float_t eta, Bool_t debug)
+Float_t AlidNdEtaCorrection::GetMeasuredFraction(CorrectionType correctionType, Float_t ptCutOff, Float_t eta, Bool_t debug)
 {
   // calculates the fraction of particles measured (some are missed due to the pt cut off)
-  // uses the generated particle histogram from fTrack2ParticleCorrection
+  //
+  // uses the generated particle histogram from the correction passed, e.g. pass GetTrack2ParticleCorrection()
 
-  const TH3F* generated = fTrack2ParticleCorrection->GetGeneratedHistogram();
+  const TH3F* generated = 0;
 
+  switch (correctionType)
+  {
+    case kNone : return -1;
+    case kTrack2Particle : generated = fTrack2ParticleCorrection->GetTrackCorrection()->GetGeneratedHistogram(); break;
+    case kVertexReco : generated = fVertexRecoCorrection->GetTrackCorrection()->GetGeneratedHistogram(); break;
+    case kINEL : generated = fTriggerBiasCorrectionMBToINEL->GetTrackCorrection()->GetGeneratedHistogram(); break;
+    case kNSD: generated = fTriggerBiasCorrectionMBToNSD->GetTrackCorrection()->GetGeneratedHistogram(); break;
+    case kND: generated = fTriggerBiasCorrectionMBToND->GetTrackCorrection()->GetGeneratedHistogram(); break;
+  }
+  
   // find eta borders, if eta is negative assume -0.8 ... 0.8
   Int_t etaBegin = 0;
   Int_t etaEnd = 0;
@@ -317,13 +301,15 @@ Float_t AlidNdEtaCorrection::GetMeasuredFraction(Float_t ptCutOff, Float_t eta,
     etaEnd = etaBegin;
   }
 
-  Int_t vertexBegin = generated->GetXaxis()->FindBin(-10);
-  Int_t vertexEnd = generated->GetXaxis()->FindBin(10);
+  Int_t vertexBegin = generated->GetXaxis()->FindBin(-19.99);
+  Int_t vertexEnd = generated->GetXaxis()->FindBin(19.99);
 
   TH1D* ptProj = dynamic_cast<TH1D*> (generated->ProjectionZ(Form("%s_pt", generated->GetName()), vertexBegin, vertexEnd, etaBegin, etaEnd));
+  printf("GetMeasuredFraction: bin range %d %d %d %d\n", vertexBegin, vertexEnd, etaBegin, etaEnd);
   ptProj->GetXaxis()->SetTitle(generated->GetZaxis()->GetTitle());
 
   Int_t ptBin = ptProj->FindBin(ptCutOff);
+  printf("GetMeasuredFraction: bin range %d %d\n", ptBin, ptProj->GetNbinsX());
   Float_t abovePtCut = ptProj->Integral(ptBin, ptProj->GetNbinsX());
   Float_t all = ptProj->Integral();
 
@@ -332,6 +318,8 @@ Float_t AlidNdEtaCorrection::GetMeasuredFraction(Float_t ptCutOff, Float_t eta,
 
   Float_t fraction = abovePtCut / all;
 
+  printf("GetMeasuredFraction: all %f above %f fraction %f\n", all, abovePtCut, fraction);
+
   if (debug)
   {
     new TCanvas;
@@ -346,6 +334,7 @@ Float_t AlidNdEtaCorrection::GetMeasuredFraction(Float_t ptCutOff, Float_t eta,
   return fraction;
 }
 
+//____________________________________________________________________
 void AlidNdEtaCorrection::ReduceInformation()
 {
   // this function deletes the measured and generated histograms from the corrections to reduce the amount of data
index 1cba5b0d0dc0bad0a3702b6c9aaba9032db81052..71c740cd4effdcc8e7c141ad9a97f08b26cc3509 100644 (file)
 
 #include <TNamed.h>
 
-#include <AliCorrectionMatrix2D.h>
-#include <AliCorrectionMatrix3D.h>
-#include <AliLog.h>
+class AliCorrection;
 
 class AlidNdEtaCorrection : public TNamed
 {
 public:
+  enum CorrectionType {
+    kNone,
+    kTrack2Particle,  // measured events
+    kVertexReco,      // MB sample
+    kINEL,
+    kNSD,
+    kND
+  };
+
   AlidNdEtaCorrection();
   AlidNdEtaCorrection(const Char_t* name, const Char_t* title);
 
+  virtual Long64_t Merge(TCollection* list);
+
   ~AlidNdEtaCorrection();
 
-  // fTrack2ParticleCorrection
-  void FillParticle(Float_t vtx, Float_t eta, Float_t pt)                  {fTrack2ParticleCorrection->FillGene(vtx, eta, pt);}
-  void FillParticleWhenMeasuredTrack(Float_t vtx, Float_t eta, Float_t pt) {fTrack2ParticleCorrection->FillMeas(vtx, eta, pt);}
+  void FillMCParticle(Float_t vtx, Float_t eta, Float_t pt, Bool_t trigger, Bool_t vertex, Int_t processType);
+  void FillTrackedParticle(Float_t vtx, Float_t eta, Float_t pt);
+  void FillEvent(Float_t vtx, Float_t n, Bool_t trigger, Bool_t vertex, Int_t processType);
 
-  // fVertexRecoCorrection & fTriggerBiasCorrection
-  void FillEventWithTriggerWithReconstructedVertex(Float_t vtx, Float_t n) {fVertexRecoCorrection->FillMeas(vtx, n);}
-  void FillEventWithTrigger(Float_t vtx, Float_t n);                         
-  void FillEventAll(Float_t vtx, Float_t n, Char_t* opt="INEL"); 
-    
   void Finish();
 
-  AliCorrectionMatrix3D* GetTrack2ParticleCorrection()  {return fTrack2ParticleCorrection;}
-  AliCorrectionMatrix2D* GetVertexRecoCorrection()      {return fVertexRecoCorrection;}
-  AliCorrectionMatrix2D* GetTriggerBiasCorrectionINEL() {return fTriggerBiasCorrectionMBToINEL;}
-  AliCorrectionMatrix2D* GetTriggerBiasCorrectionNSD()  {return fTriggerBiasCorrectionMBToNSD;}
-  AliCorrectionMatrix2D* GetTriggerBiasCorrectionND()   {return fTriggerBiasCorrectionMBToND;}
-  AliCorrectionMatrix2D* GetTriggerBiasCorrection(Char_t* opt="INEL");
-
-  virtual Long64_t Merge(TCollection* list);
+  AliCorrection* GetTrack2ParticleCorrection()  {return fTrack2ParticleCorrection;}
+  AliCorrection* GetVertexRecoCorrection()      {return fVertexRecoCorrection;}
+  AliCorrection* GetTriggerBiasCorrectionINEL() {return fTriggerBiasCorrectionMBToINEL;}
+  AliCorrection* GetTriggerBiasCorrectionNSD()  {return fTriggerBiasCorrectionMBToNSD;}
+  AliCorrection* GetTriggerBiasCorrectionND()   {return fTriggerBiasCorrectionMBToND;}
 
   void    SaveHistograms();
-  Bool_t  LoadHistograms(const Char_t* fileName, const Char_t* dir = "dndeta_correction");
-  Bool_t  LoadCorrection(const Char_t* fileName, const Char_t* dir = "dndeta_correction") {return LoadHistograms(fileName, dir);}
-  
+  Bool_t  LoadHistograms(const Char_t* dir = 0);
   void DrawHistograms();
-  
-  //  void RemoveEdges(Float_t cut=2, Int_t nBinsVtx=0, Int_t nBinsEta=0);
-  
-  Float_t GetTrack2ParticleCorrection(Float_t vtx, Float_t eta, Float_t pt)
-    {return fTrack2ParticleCorrection->GetCorrection(vtx, eta, pt);}
-  
-  Float_t GetVertexRecoCorrection     (Float_t vtx, Float_t n) {return fVertexRecoCorrection->GetCorrection(vtx, n);}
-  Float_t GetTriggerBiasCorrectionINEL(Float_t vtx, Float_t n) {return fTriggerBiasCorrectionMBToINEL->GetCorrection(vtx, n);} 
-  Float_t GetTriggerBiasCorrectionNSD (Float_t vtx, Float_t n) {return fTriggerBiasCorrectionMBToNSD->GetCorrection(vtx, n);} 
-  Float_t GetTriggerBiasCorrectionND  (Float_t vtx, Float_t n) {return fTriggerBiasCorrectionMBToND->GetCorrection(vtx, n);} 
-  Float_t GetTriggerBiasCorrection    (Float_t vtx, Float_t n, Char_t* opt="INEL");
-
-  Float_t GetMeasuredFraction(Float_t ptCutOff, Float_t eta = -100, Bool_t debug = kFALSE);
 
+  Float_t GetMeasuredFraction(CorrectionType correctionType, Float_t ptCutOff, Float_t eta = -100, Bool_t debug = kFALSE);
 
   void ReduceInformation();
 
 protected:
-  AliCorrectionMatrix3D* fTrack2ParticleCorrection;       //-> handles the track-to-particle correction, function of vtx_z, eta, pt
-  AliCorrectionMatrix2D* fVertexRecoCorrection;           //-> handles the vertex reconstruction efficiency, function of n_clustersITS and vtx_z
-  AliCorrectionMatrix2D* fTriggerBiasCorrectionMBToINEL;  //-> handles the trigger bias MB->INEL, function of n and vtx_z
-  AliCorrectionMatrix2D* fTriggerBiasCorrectionMBToNSD;   //-> handles the trigger bias MB->NSD,  function of n and vtx_z
-  AliCorrectionMatrix2D* fTriggerBiasCorrectionMBToND;    //-> handles the trigger bias MB->ND,   function of n and vtx_z
-
-  //Long64_t fNEvents;
-  //Long64_t fNTriggeredEvents;
+  AliCorrection* fTrack2ParticleCorrection;       //-> handles the track-to-particle correction (only track level (vtx_z, eta, pt))
+  AliCorrection* fVertexRecoCorrection;           //-> handles the vertex reconstruction efficiency, (n, vtx_z)
+  AliCorrection* fTriggerBiasCorrectionMBToINEL;  //-> handles the trigger bias MB->INEL, function of n and vtx_z
+  AliCorrection* fTriggerBiasCorrectionMBToNSD;   //-> handles the trigger bias MB->NSD,  function of n and vtx_z
+  AliCorrection* fTriggerBiasCorrectionMBToND;    //-> handles the trigger bias MB->ND,   function of n and vtx_z
 
 private:
   AlidNdEtaCorrection(const AlidNdEtaCorrection&);
index 74e018cf1e67014461b335f4e7ec46a2d52bc235..a0e1079d22f7326dea11f907367f85c33f8d609d 100644 (file)
 #include "AliPWG0Helper.h"
 #include "AliPWG0depHelper.h"
 
+#include "dNdEta/dNdEtaAnalysis.h"
+
 ClassImp(AlidNdEtaCorrectionSelector)
 
 AlidNdEtaCorrectionSelector::AlidNdEtaCorrectionSelector() :
   AliSelectorRL(),
   fEsdTrackCuts(0),
   fdNdEtaCorrection(0),
+  fdNdEtaAnalysisMC(0),
+  fdNdEtaAnalysisESD(0),
   fPIDParticles(0),
   fPIDTracks(0),
   fClustersITSPos(0),
@@ -148,6 +152,9 @@ void AlidNdEtaCorrectionSelector::SlaveBegin(TTree * tree)
 
   fClustersITSNeg = new TH1F("clusters_its_neg", "clusters_its_neg", 7, -0.5, 6.5);
   fClustersTPCNeg = new TH1F("clusters_tpc_neg", "clusters_tpc_neg", 160, -0.5, 159.5);
+
+  fdNdEtaAnalysisMC = new dNdEtaAnalysis("dndetaMC", "dndetaMC");
+  fdNdEtaAnalysisESD = new dNdEtaAnalysis("dndetaESD", "dndetaESD");
 }
 
 void AlidNdEtaCorrectionSelector::Init(TTree* tree)
@@ -230,9 +237,17 @@ Bool_t AlidNdEtaCorrectionSelector::Process(Long64_t entry)
   // get the MC vertex
   AliGenEventHeader* genHeader = header->GenEventHeader();
 
+  // primary vertex (from MC)
   TArrayF vtxMC(3);
   genHeader->PrimaryVertex(vtxMC);
 
+  // get process type
+  Int_t processType = AliPWG0depHelper::GetPythiaEventProcessType(header);
+  AliDebug(AliLog::kDebug+1, Form("Found pythia procces type %d", processType));
+
+  if (processType<0)
+    AliDebug(AliLog::kError, Form("Unknown Pythia process type %d.", processType));
+
   // loop over mc particles
   Int_t nPrim  = stack->GetNprimary();
 
@@ -257,13 +272,16 @@ Bool_t AlidNdEtaCorrectionSelector::Process(Long64_t entry)
     Float_t eta = particle->Eta();
     Float_t pt = particle->Pt();
 
-    if (vertexReconstructed) {
-      fdNdEtaCorrection->FillParticle(vtxMC[2], eta, pt);
+    fdNdEtaCorrection->FillMCParticle(vtxMC[2], eta, pt, eventTriggered, vertexReconstructed, processType);
 
-      if (pt > 0.1 && pt < 0.2)
-       fPIDParticles->Fill(particle->GetPdgCode());
+    if (eventTriggered)
+    {
+      if (vertexReconstructed)
+      {
+        fdNdEtaAnalysisMC->FillTrack(vtxMC[2], eta, pt);
+      }
     }
-  }// end of mc particle
+  } // end of mc particle
 
   // ########################################################
   // loop over esd tracks
@@ -301,12 +319,13 @@ Bool_t AlidNdEtaCorrectionSelector::Process(Long64_t entry)
     if (SignOK(particle->GetPDG()) == kFALSE)
         continue;
 
-    if (vertexReconstructed)
+    if (eventTriggered && vertexReconstructed)
     {
-      fdNdEtaCorrection->FillParticleWhenMeasuredTrack(vtxMC[2], particle->Eta(), particle->Pt());
+      fdNdEtaCorrection->FillTrackedParticle(vtxMC[2], particle->Eta(), particle->Pt());
+      fdNdEtaAnalysisESD->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
       if (particle->Pt() > 0.1 && particle->Pt() < 0.2)
-       {
-         fPIDTracks->Fill(particle->GetPdgCode());
+      {
+        fPIDTracks->Fill(particle->GetPdgCode());
         if (particle->GetPDG()->Charge() > 0)
         {
           fClustersITSPos->Fill(esdTrack->GetITSclusters(0));
@@ -321,28 +340,23 @@ Bool_t AlidNdEtaCorrectionSelector::Process(Long64_t entry)
     }
   } // end of track loop
 
+  // get the ESD vertex
+  const AliESDVertex* vtxESD = fESD->GetVertex();
+  // FAKE test!
+  //nGoodTracks = vtxESD->GetNContributors();
+
+  if (eventTriggered && vertexReconstructed)
+    fdNdEtaAnalysisMC->FillEvent(vtxMC[2], nGoodTracks);
+
   // stuff regarding the vertex reco correction and trigger bias correction
+  fdNdEtaCorrection->FillEvent(vtxMC[2], nGoodTracks, eventTriggered, vertexReconstructed, processType);
   if (eventTriggered) {
-    fdNdEtaCorrection->FillEventWithTrigger(vtxMC[2], nGoodTracks);
     if (vertexReconstructed)
-      fdNdEtaCorrection->FillEventWithTriggerWithReconstructedVertex(vtxMC[2], nGoodTracks);
+    {
+      fdNdEtaAnalysisESD->FillEvent(vtxMC[2], nGoodTracks);
+    }
   }
 
-  // getting process information
-  Int_t processtype = AliPWG0depHelper::GetPythiaEventProcessType(header);
-  AliDebug(AliLog::kDebug+1,Form(" Found pythia procces type %d", processtype));
-
-  if (processtype<0)
-    AliDebug(AliLog::kError, Form("Unkown Pythia process type %d.", processtype));
-  
-  fdNdEtaCorrection->FillEventAll(vtxMC[2], nGoodTracks, "INEL");
-  
-  if (processtype!=92 && processtype!=93)
-    fdNdEtaCorrection->FillEventAll(vtxMC[2], nGoodTracks, "NSD");
-  
-  if (processtype!=92 && processtype!=93 && processtype!=94)
-    fdNdEtaCorrection->FillEventAll(vtxMC[2], nGoodTracks, "ND");
-
   return kTRUE;
 }
 
@@ -362,6 +376,8 @@ void AlidNdEtaCorrectionSelector::SlaveTerminate()
   }
 
   fOutput->Add(fdNdEtaCorrection);
+  fOutput->Add(fdNdEtaAnalysisMC);
+  fOutput->Add(fdNdEtaAnalysisESD);
 }
 
 void AlidNdEtaCorrectionSelector::Terminate()
@@ -373,7 +389,9 @@ void AlidNdEtaCorrectionSelector::Terminate()
   AliSelectorRL::Terminate();
 
   fdNdEtaCorrection = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction"));
-  if (!fdNdEtaCorrection)
+  fdNdEtaAnalysisMC = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaMC"));
+  fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaESD"));
+  if (!fdNdEtaCorrection || !fdNdEtaAnalysisMC || !fdNdEtaAnalysisESD)
   {
     AliDebug(AliLog::kError, "Could not read object from output list");
     return;
@@ -386,6 +404,8 @@ void AlidNdEtaCorrectionSelector::Terminate()
   if (fEsdTrackCuts)
     fEsdTrackCuts->SaveHistograms("esd_track_cuts");
   fdNdEtaCorrection->SaveHistograms();
+  fdNdEtaAnalysisMC->SaveHistograms();
+  fdNdEtaAnalysisESD->SaveHistograms();
 
   fout->Write();
   fout->Close();
index 3a5e3a1f5d4b55ae0d4cc2baaed24df38f878ef2..da6b4948658919b48d6fb424773d297f0daa9c3f 100644 (file)
@@ -9,6 +9,7 @@ class AliESDtrackCuts;
 class AlidNdEtaCorrection;
 class TH1F;
 class TParticlePDG;
+class dNdEtaAnalysis;
 
 class AlidNdEtaCorrectionSelector : public AliSelectorRL {
   public:
@@ -30,6 +31,9 @@ class AlidNdEtaCorrectionSelector : public AliSelectorRL {
 
     AlidNdEtaCorrection* fdNdEtaCorrection;      // contains the intermediate histograms (on each slave)
 
+    dNdEtaAnalysis* fdNdEtaAnalysisMC; // analysis from MC (only triggered, vertex events)
+    dNdEtaAnalysis* fdNdEtaAnalysisESD; // analysis from ESD (not yet corrected!)
+
     TH1F* fPIDParticles; // pid of primary particles
     TH1F* fPIDTracks; // pid of reconstructed tracks
 
index b6aa85a6e044611885842b65ab75934f108af559..899f22eb44a80bef49554b431a6a7adb57c51b8e 100644 (file)
@@ -4,17 +4,21 @@
 
 #include <TFile.h>
 #include <TH3F.h>
-#include <TH2D.h>
-#include <TH1D.h>
+#include <TH2F.h>
+#include <TH1F.h>
 #include <TMath.h>
 #include <TCanvas.h>
 #include <TCollection.h>
 #include <TIterator.h>
 #include <TList.h>
 #include <TLegend.h>
+#include <TLine.h>
 
 #include "AlidNdEtaCorrection.h"
-#include "AliPWG0Helper.h"
+#include <AliCorrection.h>
+#include <AliPWG0Helper.h>
+#include <AliCorrectionMatrix2D.h>
+#include <AliCorrectionMatrix3D.h>
 
 //____________________________________________________________________
 ClassImp(dNdEtaAnalysis)
@@ -23,8 +27,6 @@ ClassImp(dNdEtaAnalysis)
 dNdEtaAnalysis::dNdEtaAnalysis() :
   TNamed(),
   fData(0),
-  fDataUncorrected(0),
-  fVtx(0),
   fPtDist(0)
 {
   // default constructor
@@ -40,45 +42,31 @@ dNdEtaAnalysis::dNdEtaAnalysis() :
 dNdEtaAnalysis::dNdEtaAnalysis(Char_t* name, Char_t* title) :
   TNamed(name, title),
   fData(0),
-  fDataUncorrected(0),
-  fVtx(0),
   fPtDist(0)
 {
   // constructor
 
-  fData  = new TH3F(Form("%s_analysis", name),"Input data",80,-20,20,40,-2,2,100, 0, 10);
-  fData->SetXTitle("vtx z [cm]");
-  fData->SetYTitle("#eta");
-  fData->SetZTitle("p_{T}");
-  fData->Sumw2();
+  Float_t binLimitsPt[] = {0.0, 0.05, 0.1, 0.125, 0.15, 0.175, 0.2, 0.225, 0.25, 0.275, 0.3, 0.325, 0.35, 0.375, 0.4, 0.425, 0.45, 0.475, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.5, 2.0, 5.0, 10.0, 100.0};
 
-  fDataUncorrected = dynamic_cast<TH3F*> (fData->Clone(Form("%s_analysis_uncorrected", name)));
-  fDataUncorrected->SetTitle(Form("%s uncorrected", fData->GetTitle()));
-  fVtx       = dynamic_cast<TH1D*> (fData->Project3D("x"));
-  fVtx->SetName(Form("%s_vtx", name));
-  fVtx->SetTitle("Vertex distribution");
-  fVtx->GetXaxis()->SetTitle(fData->GetXaxis()->GetTitle());
-  //fVtx->Sumw2();
+  fData = new AliCorrection("Analysis", Form("%s Analysis", title));
 
-  fdNdEta[0] = dynamic_cast<TH1D*> (fData->Project3D("y"));
-  fdNdEta[0]->SetName(Form("%s_dNdEta", name));
-  fdNdEta[0]->SetTitle("dN_{ch}/d#eta");
-  fdNdEta[0]->GetXaxis()->SetTitle(fData->GetYaxis()->GetTitle());
-  fdNdEta[0]->SetYTitle("dN_{ch}/d#eta");
+  // do not add this hists to the directory
+  Bool_t oldStatus = TH1::AddDirectoryStatus();
+  TH1::AddDirectory(kFALSE);
 
-  fdNdEtaPtCutOffCorrected[0] = dynamic_cast<TH1D*> (fdNdEta[0]->Clone(Form("%s_corrected", fdNdEta[0]->GetName())));
+  fdNdEta[0] = new TH1F("dNdEta", "dN_{ch}/d#eta;#eta;dN_{ch}/d#eta", 20, -2, 2);
+
+  fdNdEtaPtCutOffCorrected[0] = dynamic_cast<TH1F*> (fdNdEta[0]->Clone(Form("%s_corrected", fdNdEta[0]->GetName())));
 
   for (Int_t i=1; i<kVertexBinning; ++i)
   {
-    fdNdEta[i]    = dynamic_cast<TH1D*> (fdNdEta[0]->Clone(Form("%s_%d", fdNdEta[0]->GetName(), i)));
-    fdNdEtaPtCutOffCorrected[i]    = dynamic_cast<TH1D*> (fdNdEtaPtCutOffCorrected[0]->Clone(Form("%s_%d", fdNdEtaPtCutOffCorrected[0]->GetName(), i)));
+    fdNdEta[i]    = dynamic_cast<TH1F*> (fdNdEta[0]->Clone(Form("%s_%d", fdNdEta[0]->GetName(), i)));
+    fdNdEtaPtCutOffCorrected[i]    = dynamic_cast<TH1F*> (fdNdEtaPtCutOffCorrected[0]->Clone(Form("%s_%d", fdNdEtaPtCutOffCorrected[0]->GetName(), i)));
   }
 
-  fPtDist = dynamic_cast<TH1D*> (fData->Project3D("z"));
-  fPtDist->SetName(Form("%s_pt", name));
-  fPtDist->SetTitle("p_{T} [GeV/c]");
-  fPtDist->GetXaxis()->SetTitle(fData->GetZaxis()->GetTitle());
-  fPtDist->SetYTitle("#frac{dN}{d#eta dp_{T}} [c/GeV]");
+  fPtDist = new TH1F("Pt", "p_{T} distribution;p_{T} [GeV/c];#frac{dN}{d#eta dp_{T}} [c/GeV]", 28, binLimitsPt);
+
+  TH1::AddDirectory(oldStatus);
 }
 
 //____________________________________________________________________
@@ -92,18 +80,6 @@ dNdEtaAnalysis::~dNdEtaAnalysis()
     fData = 0;
   }
 
-  if (fDataUncorrected)
-  {
-    delete fDataUncorrected;
-    fDataUncorrected = 0;
-  }
-
-  if (fVtx)
-  {
-    delete fVtx;
-    fVtx = 0;
-  }
-
   for (Int_t i=0; i<kVertexBinning; ++i)
   {
     if (fdNdEta[i])
@@ -129,8 +105,6 @@ dNdEtaAnalysis::~dNdEtaAnalysis()
 dNdEtaAnalysis::dNdEtaAnalysis(const dNdEtaAnalysis &c) :
   TNamed(c),
   fData(0),
-  fDataUncorrected(0),
-  fVtx(0),
   fPtDist(0)
 {
   //
@@ -160,92 +134,147 @@ void dNdEtaAnalysis::Copy(TObject &c) const
 
   dNdEtaAnalysis& target = (dNdEtaAnalysis &) c;
 
-  target.fData = dynamic_cast<TH3F*> (fData->Clone());
-  target.fDataUncorrected = dynamic_cast<TH3F*> (fDataUncorrected->Clone());
-  target.fVtx = dynamic_cast<TH1D*> (fVtx->Clone());
+  target.fData = dynamic_cast<AliCorrection*> (fData->Clone());
 
   for (Int_t i=0; i<kVertexBinning; ++i)
   {
-    target.fdNdEta[i] = dynamic_cast<TH1D*> (fdNdEta[i]->Clone());
-    target.fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1D*> (fdNdEtaPtCutOffCorrected[i]->Clone());
+    target.fdNdEta[i] = dynamic_cast<TH1F*> (fdNdEta[i]->Clone());
+    target.fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (fdNdEtaPtCutOffCorrected[i]->Clone());
   }
 
-  target.fPtDist = dynamic_cast<TH1D*> (fPtDist->Clone());
+  target.fPtDist = dynamic_cast<TH1F*> (fPtDist->Clone());
 
   TNamed::Copy((TNamed &) c);
 }
 
 //____________________________________________________________________
-void dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t pt, Float_t weight)
+void dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t pt)
 {
   // fills a track into the histograms
 
-  fDataUncorrected->Fill(vtx, eta, pt);
-  fData->Fill(vtx, eta, pt, weight);
+  fData->GetTrackCorrection()->FillMeas(vtx, eta, pt);
 }
 
 //____________________________________________________________________
-void dNdEtaAnalysis::FillEvent(Float_t vtx, Float_t weight)
+void dNdEtaAnalysis::FillEvent(Float_t vtx, Float_t n)
 {
   // fills an event into the histograms
 
-  fVtx->Fill(vtx, weight);
+  fData->GetEventCorrection()->FillMeas(vtx, n);
 }
 
 //____________________________________________________________________
-void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut)
+void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, AlidNdEtaCorrection::CorrectionType correctionType)
 {
-  // correct with correction values if available
+  //
+  // correct with the given correction values and calculate dNdEta and pT distribution
+  // the corrections that are applied can be steered by the flag correctionType
+  //
+
+  // TODO put tag somewhere which corrections have been applied
+
+  // set corrections to 1
+  fData->SetCorrectionToUnity();
+
+  if (correction && correctionType != AlidNdEtaCorrection::kNone)
+  {
+    TH3F* trackCorr = fData->GetTrackCorrection()->GetCorrectionHistogram();
+    TH2F* eventCorr = fData->GetEventCorrection()->GetCorrectionHistogram();
+
+    if (correctionType >= AlidNdEtaCorrection::kTrack2Particle)
+      trackCorr->Multiply(correction->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetCorrectionHistogram());
 
-  // TODO what do we do with the error?
-  if (!correction)
+    if (correctionType >= AlidNdEtaCorrection::kVertexReco)
+    {
+      trackCorr->Multiply(correction->GetVertexRecoCorrection()->GetTrackCorrection()->GetCorrectionHistogram());
+      eventCorr->Multiply(correction->GetVertexRecoCorrection()->GetEventCorrection()->GetCorrectionHistogram());
+    }
+
+    switch (correctionType)
+    {
+      case AlidNdEtaCorrection::kINEL :
+      {
+        trackCorr->Multiply(correction->GetTriggerBiasCorrectionINEL()->GetTrackCorrection()->GetCorrectionHistogram());
+        eventCorr->Multiply(correction->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetCorrectionHistogram());
+        break;
+      }
+      case AlidNdEtaCorrection::kNSD :
+      {
+        trackCorr->Multiply(correction->GetTriggerBiasCorrectionNSD()->GetTrackCorrection()->GetCorrectionHistogram());
+        eventCorr->Multiply(correction->GetTriggerBiasCorrectionNSD()->GetEventCorrection()->GetCorrectionHistogram());
+        break;
+      }
+      case AlidNdEtaCorrection::kND :
+      {
+        trackCorr->Multiply(correction->GetTriggerBiasCorrectionND()->GetTrackCorrection()->GetCorrectionHistogram());
+        eventCorr->Multiply(correction->GetTriggerBiasCorrectionND()->GetEventCorrection()->GetCorrectionHistogram());
+        break;
+      }
+      default : break;
+    }
+  }
+  else
     printf("INFO: No correction applied\n");
 
-  // In fData we have the track2particle and vertex reconstruction efficiency correction already applied
+  fData->Multiply();
+
+  TH3F* dataHist = fData->GetTrackCorrection()->GetGeneratedHistogram();
+
+  // integrate multiplicity axis out (include under/overflow bins!!!)
+  TH2F* tmp = fData->GetEventCorrection()->GetGeneratedHistogram();
+  TH1D* vertexHist = tmp->ProjectionX("_px", 0, tmp->GetNbinsY() + 1, "e");
 
   // create pt hist
   {
     // reset all ranges
-    fData->GetXaxis()->SetRange(0, 0);
-    fData->GetYaxis()->SetRange(0, 0);
-    fData->GetZaxis()->SetRange(0, 0);
+    dataHist->GetXaxis()->SetRange(0, 0);
+    dataHist->GetYaxis()->SetRange(0, 0);
+    dataHist->GetZaxis()->SetRange(0, 0);
 
     // vtx cut
-    Int_t vertexBinBegin = fData->GetXaxis()->FindBin(-5);
-    Int_t vertexBinEnd = fData->GetXaxis()->FindBin(5);
-    fData->GetXaxis()->SetRange(vertexBinBegin, vertexBinEnd);
-    Float_t nEvents = fVtx->Integral(vertexBinBegin, vertexBinEnd);
+    Int_t vertexBinBegin = dataHist->GetXaxis()->FindBin(-5);
+    Int_t vertexBinEnd = dataHist->GetXaxis()->FindBin(5);
+    dataHist->GetXaxis()->SetRange(vertexBinBegin, vertexBinEnd);
+    Float_t nEvents = vertexHist->Integral(vertexBinBegin, vertexBinEnd);
 
-    // eta cut
-    fData->GetYaxis()->SetRange(fData->GetYaxis()->FindBin(-0.8), fData->GetYaxis()->FindBin(0.8));
-    Float_t etaWidth = 1.6;
+    if (nEvents > 0)
+    {
+      // eta cut
+      dataHist->GetYaxis()->SetRange(dataHist->GetYaxis()->FindBin(-0.8), dataHist->GetYaxis()->FindBin(0.8));
+      Float_t etaWidth = 1.6;
 
-    TH1D* ptHist = dynamic_cast<TH1D*> (fData->Project3D("ze"));
+      TH1D* ptHist = dynamic_cast<TH1D*> (dataHist->Project3D("ze"));
 
-    for (Int_t i=1; i<=fPtDist->GetNbinsX(); ++i)
-    {
-      Float_t binSize = fPtDist->GetBinWidth(i);
-      fPtDist->SetBinContent(i, ptHist->GetBinContent(i) / binSize / nEvents / etaWidth);
-      fPtDist->SetBinError(i, ptHist->GetBinError(i) / binSize / nEvents / etaWidth);
-    }
+      for (Int_t i=1; i<=fPtDist->GetNbinsX(); ++i)
+      {
+        Float_t binSize = fPtDist->GetBinWidth(i);
+        fPtDist->SetBinContent(i, ptHist->GetBinContent(i) / binSize / nEvents / etaWidth);
+        fPtDist->SetBinError(i, ptHist->GetBinError(i) / binSize / nEvents / etaWidth);
+      }
 
-    delete ptHist;
+      delete ptHist;
+    }
+    else
+      printf("ERROR: nEvents is 0!\n");
   }
 
   // reset all ranges
-  fData->GetXaxis()->SetRange(0, 0);
-  fData->GetYaxis()->SetRange(0, 0);
-  fData->GetZaxis()->SetRange(0, 0);
+  dataHist->GetXaxis()->SetRange(0, 0);
+  dataHist->GetYaxis()->SetRange(0, 0);
+  dataHist->GetZaxis()->SetRange(0, 0);
 
   // integrate over pt (with pt cut)
   Int_t ptLowBin = 1;
   if (ptCut > 0)
-    ptLowBin = fData->GetZaxis()->FindBin(ptCut);
+    ptLowBin = dataHist->GetZaxis()->FindBin(ptCut);
 
-  fData->GetZaxis()->SetRange(ptLowBin, fData->GetZaxis()->GetNbins());
-  TH2D* vtxVsEta = dynamic_cast<TH2D*> (fData->Project3D("yx2e"));
-  vtxVsEta->GetXaxis()->SetTitle(fData->GetXaxis()->GetTitle());
-  vtxVsEta->GetYaxis()->SetTitle(fData->GetYaxis()->GetTitle());
+  dataHist->GetZaxis()->SetRange(ptLowBin, dataHist->GetZaxis()->GetNbins());
+  printf("range %d %d\n", ptLowBin, dataHist->GetZaxis()->GetNbins());
+  TH2D* vtxVsEta = dynamic_cast<TH2D*> (dataHist->Project3D("yx2e"));
+
+  dataHist->GetZaxis()->SetRange(0, 0);
+  vtxVsEta->GetXaxis()->SetTitle(dataHist->GetXaxis()->GetTitle());
+  vtxVsEta->GetYaxis()->SetTitle(dataHist->GetYaxis()->GetTitle());
 
   if (vtxVsEta == 0)
   {
@@ -253,17 +282,14 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut)
     return;
   }
 
-  //new TCanvas;
-  //vtxVsEta->Draw("COLZ");
-
-  for (Int_t iEta=0; iEta<=vtxVsEta->GetNbinsY(); iEta++)
+  for (Int_t iEta=1; iEta<=vtxVsEta->GetNbinsY(); iEta++)
   {
     // do we have several histograms for different vertex positions?
-    Int_t vertexBinWidth = fVtx->GetNbinsX() / (kVertexBinning-1);
+    Int_t vertexBinWidth = vertexHist->GetNbinsX() / (kVertexBinning-1);
     for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
     {
       Int_t vertexBinBegin = 1;
-      Int_t vertexBinEnd = fVtx->GetNbinsX() + 1;
+      Int_t vertexBinEnd = vertexHist->GetNbinsX() + 1;
 
       // the first histogram is always for the whole vertex range
       if (vertexPos > 0)
@@ -272,7 +298,7 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut)
         vertexBinEnd = vertexBinBegin + vertexBinWidth;
       }
 
-      Float_t totalEvents = fVtx->Integral(vertexBinBegin, vertexBinEnd - 1);
+      Float_t totalEvents = vertexHist->Integral(vertexBinBegin, vertexBinEnd - 1);
       if (totalEvents == 0)
       {
         printf("WARNING: No events for hist %d %d %d\n", vertexPos, vertexBinBegin, vertexBinEnd);
@@ -291,8 +317,8 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut)
       }
 
       Float_t ptCutOffCorrection = 1;
-      if (correction)
-        ptCutOffCorrection = correction->GetMeasuredFraction(ptCut, vtxVsEta->GetYaxis()->GetBinCenter(iEta));
+      if (correction && ptCut > 0)
+        ptCutOffCorrection = correction->GetMeasuredFraction(correctionType, ptCut, vtxVsEta->GetYaxis()->GetBinCenter(iEta));
 
       if (ptCutOffCorrection <= 0)
       {
@@ -300,6 +326,8 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut)
         continue;
       }
 
+      printf("Eta: %d Vertex Range: %d %d, Event Count %f, Track Sum: %f, Track Sum corrected: %f\n", iEta, vertexBinBegin, vertexBinEnd, totalEvents, sum, sum / ptCutOffCorrection);
+
       Float_t dndeta = sum / totalEvents;
       Float_t error  = TMath::Sqrt(sumError2) / totalEvents;
 
@@ -328,25 +356,11 @@ void dNdEtaAnalysis::SaveHistograms()
 
   if (fData)
   {
-    fData->Write();
-    AliPWG0Helper::CreateProjections(fData);
+    fData->SaveHistograms();
   }
   else
     printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fData is 0\n");
 
-  if (fDataUncorrected)
-  {
-    fDataUncorrected->Write();
-    AliPWG0Helper::CreateProjections(fDataUncorrected);
-  }
-  else
-    printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fDataUncorrected is 0\n");
-
-  if (fVtx)
-    fVtx       ->Write();
-  else
-    printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fVtx is 0\n");
-
   if (fPtDist)
     fPtDist       ->Write();
   else
@@ -368,87 +382,123 @@ void dNdEtaAnalysis::SaveHistograms()
   gDirectory->cd("../");
 }
 
-void dNdEtaAnalysis::LoadHistograms()
+void dNdEtaAnalysis::LoadHistograms(const Char_t* dir)
 {
   // loads the histograms from a directory with the name of this class (retrieved from TNamed)
 
-  gDirectory->cd(GetName());
+  if (!dir)
+    dir = GetName();
 
-  fData = dynamic_cast<TH3F*> (gDirectory->Get(fData->GetName()));
-  fDataUncorrected = dynamic_cast<TH3F*> (gDirectory->Get(fDataUncorrected->GetName()));
+  gDirectory->cd(dir);
 
-  fVtx = dynamic_cast<TH1D*> (gDirectory->Get(fVtx->GetName()));
+  fData->LoadHistograms();
 
   for (Int_t i=0; i<kVertexBinning; ++i)
   {
-    fdNdEta[i] = dynamic_cast<TH1D*> (gDirectory->Get(fdNdEta[i]->GetName()));
-    fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1D*> (gDirectory->Get(fdNdEtaPtCutOffCorrected[i]->GetName()));
+    fdNdEta[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEta[i]->GetName()));
+    fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEtaPtCutOffCorrected[i]->GetName()));
   }
 
-  fPtDist = dynamic_cast<TH1D*> (gDirectory->Get(fPtDist->GetName()));
+  fPtDist = dynamic_cast<TH1F*> (gDirectory->Get(fPtDist->GetName()));
 
   gDirectory->cd("../");
 }
 
 //____________________________________________________________________
-void dNdEtaAnalysis::DrawHistograms()
+void dNdEtaAnalysis::DrawHistograms(Bool_t simple)
 {
   // draws the histograms
+  
+  if (!simple)
+  {
+    if (fData)
+      fData->DrawHistograms(GetName());
 
-  TCanvas* canvas = new TCanvas("dNdEtaAnalysis", "dNdEtaAnalysis", 800, 800);
-  canvas->Divide(2, 2);
-
-  canvas->cd(1);
-  if (fData)
-    fData->Draw("COLZ");
-
-  /*canvas->cd(2);
-  if (fDataUncorrected)
-    fDataUncorrected->Draw("COLZ");*/
+    TCanvas* canvas = new TCanvas(Form("%s_dNdEtaAnalysis", GetName()), Form("%s_dNdEtaAnalysis", GetName()), 800, 400);
+    canvas->Divide(2, 1);
 
-  canvas->cd(2);
-  if (fVtx)
-    fVtx->Draw();
+    canvas->cd(1);
+    if (fdNdEtaPtCutOffCorrected[0])
+      fdNdEtaPtCutOffCorrected[0]->Draw();
 
-  canvas->cd(3);
-  if (fdNdEtaPtCutOffCorrected[0])
-    fdNdEtaPtCutOffCorrected[0]->Draw();
+    if (fdNdEta[0])
+    {
+      fdNdEta[0]->SetLineColor(kRed);
+      fdNdEta[0]->Draw("SAME");
+    }
 
-  if (fdNdEta[0])
-  {
-    fdNdEta[0]->SetLineColor(kRed);
-    fdNdEta[0]->Draw("SAME");
+    canvas->cd(2);
+    if (fPtDist)
+      fPtDist->Draw();
   }
 
-  canvas->cd(4);
-  if (fPtDist)
-    fPtDist->Draw();
-
     // histograms for different vertices?
   if (kVertexBinning > 0)
   {
       // doesnt work, but i dont get it, giving up...
-    /*TCanvas* canvas2 =*/ new TCanvas("dNdEtaAnalysisVtx", "dNdEtaAnalysisVtx", 450, 450);
+    TCanvas* canvas2 = new TCanvas(Form("%s_dNdEtaAnalysisVtx", GetName()), Form("%s_dNdEtaAnalysisVtx", GetName()), 450, 450);
+    TCanvas* canvas3 = 0;
+    if (!simple)
+      canvas3 = new TCanvas(Form("%s_dNdEtaAnalysisVtx_noptcutoff", GetName()), Form("%s_dNdEtaAnalysisVtx_noptcutoff", GetName()), 450, 450);
+
     //Int_t yPads = (Int_t) TMath::Ceil(((Double_t) kVertexBinning - 1) / 2);
     //printf("%d\n", yPads);
     //canvas2->Divide(2, yPads);
 
-    TLegend* legend = new TLegend(0.7, 0.7, 0.9, 0.9);
+    TLegend* legend = new TLegend(0.4, 0.7, 0.6, 0.9);
 
     for (Int_t i=0; i<kVertexBinning; ++i)
     {
-      //canvas2->cd(i-1);
-      //printf("%d\n", i);
       if (fdNdEtaPtCutOffCorrected[i])
       {
+        canvas2->cd();
+
         fdNdEtaPtCutOffCorrected[i]->SetLineColor(i+1);
         fdNdEtaPtCutOffCorrected[i]->Draw((i == 0) ? "" : "SAME");
         legend->AddEntry(fdNdEtaPtCutOffCorrected[i], (i == 0) ? "Vtx All" : Form("Vtx Bin %d", i-1));
       }
+      if (canvas3 && fdNdEta[i])
+      {
+        canvas3->cd();
+
+        fdNdEta[i]->SetLineColor(i+1);
+        fdNdEta[i]->Draw((i == 0) ? "" : "SAME");
+      }
     }
 
+    canvas2->cd();
     legend->Draw();
+    canvas2->SaveAs(Form("%s_%s.gif", canvas2->GetName(), GetName()));
+
+    if (canvas3)
+    {
+      canvas3->cd();
+      legend->Draw();
+    }
   }
+  
+  if (kVertexBinning == 3)
+  {
+     TH1* clone = dynamic_cast<TH1*> (fdNdEtaPtCutOffCorrected[1]->Clone("clone"));
+     TH1* clone2 = dynamic_cast<TH1*> (fdNdEtaPtCutOffCorrected[2]->Clone("clone2"));
+     
+     if (clone && clone2)
+     {
+        TCanvas* canvas4 = new TCanvas(Form("%s_dNdEtaAnalysisVtxRatios", GetName()), Form("%s_dNdEtaAnalysisVtxRatios", GetName()), 450, 450);
+    
+        clone->Divide(fdNdEtaPtCutOffCorrected[0]);
+        clone->GetYaxis()->SetRangeUser(0.95, 1.05);
+        clone->Draw();
+        
+        clone2->Divide(fdNdEtaPtCutOffCorrected[0]);
+        clone2->Draw("SAME");
+
+        TLine* line = new TLine(-1, 1, 1, 1);
+        line->Draw();
+
+        canvas4->SaveAs(Form("%s_%s.gif", canvas4->GetName(), GetName()));
+     }
+   }
 }
 
 Long64_t dNdEtaAnalysis::Merge(TCollection* list)
@@ -467,7 +517,7 @@ Long64_t dNdEtaAnalysis::Merge(TCollection* list)
   TObject* obj;
 
   // sub collections
-  const Int_t nCollections = 2 * kVertexBinning + 4; // 4 standalone hists, two arrays of size kVertexBinning
+  const Int_t nCollections = 2 * kVertexBinning + 2; // 2 standalone hists, two arrays of size kVertexBinning
   TList* collections[nCollections];
   for (Int_t i=0; i<nCollections; ++i)
     collections[i] = new TList;
@@ -480,29 +530,23 @@ Long64_t dNdEtaAnalysis::Merge(TCollection* list)
       continue;
 
     collections[0]->Add(entry->fData);
-    if (fDataUncorrected)
-      collections[1]->Add(entry->fDataUncorrected);
-    collections[2]->Add(entry->fVtx);
-    collections[3]->Add(entry->fPtDist);
+    collections[1]->Add(entry->fPtDist);
 
     for (Int_t i=0; i<kVertexBinning; ++i)
     {
-      collections[4+i]->Add(entry->fdNdEta[i]);
-      collections[4+kVertexBinning+i]->Add(entry->fdNdEtaPtCutOffCorrected[i]);
+      collections[2+i]->Add(entry->fdNdEta[i]);
+      collections[2+kVertexBinning+i]->Add(entry->fdNdEtaPtCutOffCorrected[i]);
     }
 
     ++count;
   }
 
   fData->Merge(collections[0]);
-  if (fDataUncorrected)
-    fDataUncorrected->Merge(collections[1]);
-  fVtx->Merge(collections[2]);
-  fPtDist->Merge(collections[3]);
+  fPtDist->Merge(collections[1]);
   for (Int_t i=0; i<kVertexBinning; ++i)
   {
-    fdNdEta[i]->Merge(collections[4+i]);
-    fdNdEta[i]->Merge(collections[4+kVertexBinning+i]);
+    fdNdEta[i]->Merge(collections[2+i]);
+    fdNdEta[i]->Merge(collections[2+kVertexBinning+i]);
   }
 
   for (Int_t i=0; i<nCollections; ++i)
index 2452191ffdeda458b40efdededcc5e883dfdae79..c8568d56120ea8a43bd391f3798941e735d5cb8d 100644 (file)
 // - add functionality to make dn/deta for different mult classes?
 
 #include <TNamed.h>
+#include "AlidNdEtaCorrection.h"
 
-class TH3F;
-class TH1D;
+class TH1F;
 class TCollection;
+
 class AlidNdEtaCorrection;
+class AliCorrection;
 
 class dNdEtaAnalysis : public TNamed
 {
 public:
-  enum { kVertexBinning = 1+3 }; // the first is for the whole vertex range, the others divide the vertex range
+  enum { kVertexBinning = 1+2 }; // the first is for the whole vertex range, the others divide the vertex range
 
   dNdEtaAnalysis();
   dNdEtaAnalysis(Char_t* name, Char_t* title);
@@ -37,34 +39,31 @@ public:
   dNdEtaAnalysis &operator=(const dNdEtaAnalysis &c);
   virtual void Copy(TObject &c) const;
 
-  void FillTrack(Float_t vtx, Float_t eta, Float_t pt, Float_t weight);
-  void FillEvent(Float_t vtx, Float_t weight);
+  void FillTrack(Float_t vtx, Float_t eta, Float_t pt);
+  void FillEvent(Float_t vtx, Float_t n);
 
-  void Finish(AlidNdEtaCorrection* correction, Float_t ptCut);
+  void Finish(AlidNdEtaCorrection* correction, Float_t ptCut, AlidNdEtaCorrection::CorrectionType correctionType);
 
-  void DrawHistograms();
-  void LoadHistograms();
+  void DrawHistograms(Bool_t simple = kFALSE);
+  void LoadHistograms(const Char_t* dir = 0);
   void SaveHistograms();
 
   virtual Long64_t Merge(TCollection* list);
 
-  TH3F* GetHistogram() { return fData; }
-  TH3F* GetUncorrectedHistogram() { return fDataUncorrected; }
-  TH1D* GetVtxHistogram() { return fVtx; }
-  TH1D* GetPtHistogram() { return fPtDist; }
-  TH1D* GetdNdEtaHistogram(Int_t i = 0) { return fdNdEta[i]; }
-  TH1D* GetdNdEtaPtCutOffCorrectedHistogram(Int_t i = 0) { return fdNdEtaPtCutOffCorrected[i]; }
+  AliCorrection* GetData() { return fData; }
 
-protected:
-  TH3F* fData;              // histogram Eta vs vtx (track count)
-  TH3F* fDataUncorrected;   // uncorrected histograms Eta vs vtx (track count)
+  TH1F* GetPtHistogram() { return fPtDist; }
 
-  TH1D* fVtx;                   // vtx histogram (event count)
+  TH1F* GetdNdEtaHistogram(Int_t i = 0) { return fdNdEta[i]; }
+  TH1F* GetdNdEtaPtCutOffCorrectedHistogram(Int_t i = 0) { return fdNdEtaPtCutOffCorrected[i]; }
+
+protected:
+  AliCorrection* fData;     // we store the data in an AliCorrection
 
-  TH1D* fPtDist; // pt distribution
+  TH1F* fPtDist; // pt distribution
 
-  TH1D* fdNdEta[kVertexBinning]; // dndeta results for different vertex bins (0 = full range)
-  TH1D* fdNdEtaPtCutOffCorrected[kVertexBinning];  // dndeta results for different vertex bins (0 = full range), pt cut off corrected
+  TH1F* fdNdEta[kVertexBinning]; // dndeta results for different vertex bins (0 = full range)
+  TH1F* fdNdEtaPtCutOffCorrected[kVertexBinning];  // dndeta results for different vertex bins (0 = full range), pt cut off corrected
 
   ClassDef(dNdEtaAnalysis, 1)
 };
index 420370519611d77a4123ade5b63cd4c0a890d268..d32a048183a93db6772d7fa5f0bb6fa1c5c40bbe 100644 (file)
@@ -205,30 +205,47 @@ void CompareTrack2ParticleWithAnalysisData(const char* correctionMapFile = "corr
 void dNdEta(Bool_t onlyESD = kFALSE)
 {
   TFile* file = TFile::Open("analysis_esd.root");
-  TH1* histESD = (TH1*) file->Get("dndeta/dndeta_dNdEta_corrected_2");
-  TH1* histESDNoPt = (TH1*) file->Get("dndeta/dndeta_dNdEta_2");
-  TH1* histESDMB = (TH1*) file->Get("dndeta_mb/dndeta_mb_dNdEta_corrected_2");
-  TH1* histESDMBVtx = (TH1*) file->Get("dndeta_mbvtx/dndeta_mbvtx_dNdEta_corrected_2");
+  TH1* histESD = (TH1*) file->Get("dndeta/dNdEta_corrected");
+  TH1* histESDNoPt = (TH1*) file->Get("dndeta/dNdEta");
+  TH1* histESDMB = (TH1*) file->Get("dndetaTr/dNdEta_corrected");
+  TH1* histESDMBNoPt = (TH1*) file->Get("dndetaTr/dNdEta");
+  TH1* histESDMBVtx = (TH1*) file->Get("dndetaTrVtx/dNdEta_corrected");
+  TH1* histESDMBVtxNoPt = (TH1*) file->Get("dndetaTrVtx/dNdEta");
 
   TCanvas* canvas = new TCanvas("dNdEta1", "dNdEta1", 500, 500);
 
   Prepare1DPlot(histESD);
-  Prepare1DPlot(histESDNoPt);
   Prepare1DPlot(histESDMB);
   Prepare1DPlot(histESDMBVtx);
 
-  histESD->SetLineColor(0);
-  histESDMB->SetLineColor(0);
-  histESDMBVtx->SetLineColor(0);
+  Prepare1DPlot(histESDNoPt);
+  Prepare1DPlot(histESDMBNoPt);
+  Prepare1DPlot(histESDMBVtxNoPt);
+
+  histESD->SetLineWidth(0);
+  histESDMB->SetLineWidth(0);
+  histESDMBVtx->SetLineWidth(0);
+
+  histESDNoPt->SetLineWidth(0);
+  histESDMBNoPt->SetLineWidth(0);
+  histESDMBVtxNoPt->SetLineWidth(0);
 
   histESD->SetMarkerColor(kRed);
   histESDMB->SetMarkerColor(kBlue);
   histESDMBVtx->SetMarkerColor(103);
 
+  histESDNoPt->SetMarkerColor(kRed);
+  histESDMBNoPt->SetMarkerColor(kBlue);
+  histESDMBVtxNoPt->SetMarkerColor(103);
+
   histESD->SetMarkerStyle(20);
   histESDMB->SetMarkerStyle(21);
   histESDMBVtx->SetMarkerStyle(22);
 
+  histESDNoPt->SetMarkerStyle(20);
+  histESDMBNoPt->SetMarkerStyle(21);
+  histESDMBVtxNoPt->SetMarkerStyle(22);
+
   TH2F* dummy = new TH2F("dummy", "", 100, -1.5, 1.5, 1000, 0.1, histESDMBVtx->GetMaximum() * 1.1);
   Prepare1DPlot(dummy);
   dummy->SetStats(kFALSE);
@@ -241,7 +258,10 @@ void dNdEta(Bool_t onlyESD = kFALSE)
   histESDMBVtx->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
   histESDMB->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
   histESD->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
+
   histESDNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
+  histESDMBNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
+  histESDMBVtxNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 
   dummy->DrawCopy();
   histESDMBVtx->Draw("SAME");
@@ -254,34 +274,63 @@ void dNdEta(Bool_t onlyESD = kFALSE)
   if (onlyESD)
     return;
 
+  gSystem->Load("libPWG0base");
+
   TFile* file2 = TFile::Open("analysis_mc.root");
-  TH1* histMC = (TH1*) file2->Get("dndeta/dndeta_dNdEta_corrected_2")->Clone("cloned");
+  TH1* histMC = (TH1*) file2->Get("dndeta/dNdEta_corrected")->Clone("cloned");
+  TH1* histMCTr = (TH1*) file2->Get("dndetaTr/dNdEta_corrected")->Clone("cloned2");
+  TH1* histMCTrVtx = (TH1*) file2->Get("dndetaTrVtx/dNdEta_corrected")->Clone("cloned3");
 
-  gSystem->Load("libPWG0base");
   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
   fdNdEtaAnalysis->LoadHistograms();
-  fdNdEtaAnalysis->Finish(0, 0.3);
-  TH1* histMCPtCut = fdNdEtaAnalysis->GetdNdEtaHistogram(2);
+  fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
+  TH1* histMCPtCut = fdNdEtaAnalysis->GetdNdEtaHistogram(0);
+
+  fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTr", "dndetaTr");
+  fdNdEtaAnalysis->LoadHistograms();
+  fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
+  TH1* histMCTrPtCut = fdNdEtaAnalysis->GetdNdEtaHistogram(0);
+
+  fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
+  fdNdEtaAnalysis->LoadHistograms();
+  fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
+  TH1* histMCTrVtxPtCut = fdNdEtaAnalysis->GetdNdEtaHistogram(0);
 
   TCanvas* canvas2 = new TCanvas("dNdEta2", "dNdEta2", 500, 500);
 
   Prepare1DPlot(histMC);
+  Prepare1DPlot(histMCTr);
+  Prepare1DPlot(histMCTrVtx);
+
   Prepare1DPlot(histMCPtCut);
+  Prepare1DPlot(histMCTrPtCut);
+  Prepare1DPlot(histMCTrVtxPtCut);
 
   histMC->SetLineColor(kBlue);
-  histMCPtCut->SetLineColor(104);
-  histESDNoPt->SetLineColor(102);
+  histMCTr->SetLineColor(kBlue);
+  histMCTrVtx->SetLineColor(kBlue);
+
+  histMCPtCut->SetLineColor(kBlue);
+  histMCTrPtCut->SetLineColor(kBlue);
+  histMCTrVtxPtCut->SetLineColor(kBlue);
 
   TH2* dummy2 = (TH2F*) dummy->Clone("dummy2");
   Prepare1DPlot(dummy2);
-  dummy2->GetYaxis()->SetRangeUser(0, histESD->GetMaximum() * 1.1);
+  dummy2->GetYaxis()->SetRangeUser(0, histESDMBVtx->GetMaximum() * 1.1);
 
   dummy2->DrawCopy();
   histMC->Draw("SAME");
-//  histMC->Draw();
+  histMCTr->Draw("SAME");
+  histMCTrVtx->Draw("SAME");
   histESD->Draw("SAME");
+  histESDMB->Draw("SAME");
+  histESDMBVtx->Draw("SAME");
   histESDNoPt->Draw("SAME");
+  histESDMBNoPt->Draw("SAME");
+  histESDMBVtxNoPt->Draw("SAME");
   histMCPtCut->Draw("SAME");
+  histMCTrPtCut->Draw("SAME");
+  histMCTrVtxPtCut->Draw("SAME");
 
   canvas2->SaveAs("dNdEta2.gif");
   canvas2->SaveAs("dNdEta2.eps");
index d0d58c61c84cafba2c3b84decdf643a3c88d3b5a..f9c2cc95f5548125e22ce3df6f3dd5846cc7b387 100644 (file)
 void testAnalysis2(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aMC = kFALSE, Bool_t aDebug = kFALSE, Bool_t aProof = kFALSE, const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction", const char* option = "", const char* proofServer = "jgrosseo@lxb6046")
 {
   if (aProof)
+  {
     connectProof(proofServer);
+    gProof->AddInput(new TParameter<long>("PROOF_MaxSlavesPerNode", (long)2));
+    gProof->AddInput(new TNamed("PROOF_Packetizer", "TAdaptivePacketizer"));
+  }
 
   TString libraries("libEG;libGeom;libESD;libPWG0base");
-  TString packages("PWG0base");
+  TString packages("adaptivepacketizer;PWG0base");
   if (aMC != kFALSE)
   {
     libraries += ";libVMC;libMinuit;libSTEER;libEVGEN;libFASTSIM;libmicrocern;libpdf;libpythia6;libEGPythia6;libAliPythia6;libPWG0dep";
@@ -33,7 +37,6 @@ void testAnalysis2(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aMC = kF
 
   TChain* chain = CreateESDChain(data, nRuns, offset);
 
-
   TList inputList;
 
   if (aMC == kFALSE)
@@ -47,14 +50,6 @@ void testAnalysis2(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aMC = kF
     }
 
     inputList.Add(esdTrackCuts);
-
-    AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
-    dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
-    dNdEtaCorrection->ReduceInformation();
-    dNdEtaCorrection->SetName("dndeta_correction");
-    dNdEtaCorrection->SetTitle("dndeta_correction");
-
-    inputList.Add(dNdEtaCorrection);
   }
 
   TString selectorName = ((aMC == kFALSE) ? "AlidNdEtaAnalysisESDSelector" : "AlidNdEtaAnalysisMCSelector");
@@ -63,23 +58,98 @@ void testAnalysis2(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aMC = kF
   selectorName += ".cxx+";
 
   if (aDebug != kFALSE)
-    selectorName += "g";
+    selectorName += "+g";
 
   Int_t result = executeQuery(chain, &inputList, selectorName, option);
 
   if (result >= 0)
   {
-    dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
-
-    TFile* file = TFile::Open(aMC ? "analysis_mc.root" : "analysis_esd.root");
-    if (!file)
+    if (aMC)
     {
-      cout << "Error. File not found" << endl;
-      return;
+      dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
+
+      TFile* file = TFile::Open("analysis_mc.root");
+
+      if (!file)
+      {
+        cout << "Error. File not found" << endl;
+        return;
+      }
+      fdNdEtaAnalysis->LoadHistograms();
+      fdNdEtaAnalysis->DrawHistograms();
     }
-    fdNdEtaAnalysis->LoadHistograms();
-    fdNdEtaAnalysis->DrawHistograms();
+    else
+      FinishAnalysisAll(correctionMapFile, correctionMapFolder);
+  }
+}
+
+void FinishAnalysisAll(const char* dataInput = "analysis_esd_raw.root", const char* dataOutput = "analysis_esd.root", const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
+{
+  gSystem->Load("libPWG0base");
+
+  AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
+  TFile::Open(correctionMapFile);
+  dNdEtaCorrection->LoadHistograms();
 
-    dNdEta(kTRUE);
+  TFile* file = TFile::Open(dataInput);
+
+  if (!file)
+  {
+    cout << "Error. File not found" << endl;
+    return;
+  }
+
+  dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
+  fdNdEtaAnalysis->LoadHistograms("dndeta");
+  fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kINEL);
+  fdNdEtaAnalysis->DrawHistograms(kTRUE);
+  TFile* file2 = TFile::Open(dataOutput, "RECREATE");
+  fdNdEtaAnalysis->SaveHistograms();
+
+  file->cd();
+  fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTr", "dndetaTr");
+  fdNdEtaAnalysis->LoadHistograms("dndeta");
+  fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kVertexReco);
+  fdNdEtaAnalysis->DrawHistograms(kTRUE);
+  file2->cd();
+  fdNdEtaAnalysis->SaveHistograms();
+
+  file->cd();
+  fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
+  fdNdEtaAnalysis->LoadHistograms("dndeta");
+  fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kTrack2Particle);
+  fdNdEtaAnalysis->DrawHistograms(kTRUE);
+  file2->cd();
+  fdNdEtaAnalysis->SaveHistograms();
+}
+
+void FinishAnalysis(const char* analysisFile = "analysis_esd.root", const char* analysisDir = "dndeta", const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction", Bool_t useUncorrected = kFALSE, Bool_t simple = kFALSE)
+{
+  gSystem->Load("libPWG0base");
+
+  TFile* file = TFile::Open(analysisFile);
+
+  dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis(analysisDir, analysisDir);
+  fdNdEtaAnalysis->LoadHistograms();
+
+  if (correctionMapFile)
+  {
+    AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
+    TFile::Open(correctionMapFile);
+    dNdEtaCorrection->LoadHistograms();
+
+    fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3);
   }
+
+  fdNdEtaAnalysis->DrawHistograms(simple);
+
+  TH1* hist = fdNdEtaAnalysis->GetdNdEtaHistogram(1);
+  Int_t binLeft = hist->GetXaxis()->FindBin(-0.5);
+  Int_t binRight = hist->GetXaxis()->FindBin(0.5);
+  Float_t value1 = hist->Integral(binLeft, binRight);
+
+  hist = fdNdEtaAnalysis->GetdNdEtaHistogram(2);
+  Float_t value2 = hist->Integral(binLeft, binRight);
+
+  printf("Ratio is %f, values are %f %f\n", value1 / value2, value1, value2);
 }