improved handling of tpc only tracks
authorjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 23 Jul 2008 10:15:09 +0000 (10:15 +0000)
committerjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 23 Jul 2008 10:15:09 +0000 (10:15 +0000)
resolution histograms
clean up of unused selectors
! change of nsigma in standard cuts from 3 to 4 !

22 files changed:
PWG0/AliCorrection.cxx
PWG0/AliCorrection.h
PWG0/AliPWG0Helper.cxx
PWG0/AliPWG0Helper.h
PWG0/CreateStandardCuts.C
PWG0/dNdEta/AliMultiplicityESDSelector.cxx [deleted file]
PWG0/dNdEta/AliMultiplicityESDSelector.h [deleted file]
PWG0/dNdEta/AliMultiplicityMCSelector.cxx [deleted file]
PWG0/dNdEta/AliMultiplicityMCSelector.h [deleted file]
PWG0/dNdEta/AliVertexSelector.cxx [deleted file]
PWG0/dNdEta/AliVertexSelector.h [deleted file]
PWG0/dNdEta/AlidNdEtaCorrectionTask.cxx
PWG0/dNdEta/AlidNdEtaCorrectionTask.h
PWG0/dNdEta/AlidNdEtaTask.cxx
PWG0/dNdEta/dNdEtaAnalysis.cxx
PWG0/dNdEta/drawPlots.C
PWG0/dNdEta/makeCorrection2.C [deleted file]
PWG0/dNdEta/makeSystematics.C [deleted file]
PWG0/dNdEta/plotsMultiplicity.C [deleted file]
PWG0/dNdEta/run.C
PWG0/dNdEta/runVertexSelector.C [deleted file]
PWG0/dNdEta/rundNdEtaAnalysis.C [deleted file]

index 3b8d455..98db421 100644 (file)
@@ -369,9 +369,15 @@ void AliCorrection::Scale(Double_t factor)
 }
 
 //____________________________________________________________________
-void AliCorrection::PrintInfo(Float_t ptCut)
+void AliCorrection::PrintStats(Float_t zRange, Float_t etaRange, Float_t ptCut)
 {
-  // prints some stats
+  // prints statistics and effective correction factors
+
+  Printf("AliCorrection::PrintInfo: Values in |eta| < %.2f, |vtx-z| < %.2f and pt > %.2f:", etaRange, zRange, ptCut);
+
+  // prevent to be on bin border
+  zRange -= 0.1;
+  etaRange -= 0.1;
 
   TH3* measured = GetTrackCorrection()->GetMeasuredHistogram();
   TH3* generated = GetTrackCorrection()->GetGeneratedHistogram();
@@ -379,24 +385,18 @@ void AliCorrection::PrintInfo(Float_t ptCut)
   TH2* measuredEvents = GetEventCorrection()->GetMeasuredHistogram();
   TH2* generatedEvents = GetEventCorrection()->GetGeneratedHistogram();
 
-  Printf("AliCorrection::PrintInfo: Whole phasespace:");
-
-  Printf("tracks measured: %.1f tracks generated: %.1f, events measured: %.1f, events generated %.1f", measured->Integral(), generated->Integral(), measuredEvents->Integral(), generatedEvents->Integral());
+  Float_t tracksM = measured->Integral(measured->GetXaxis()->FindBin(-zRange), measured->GetXaxis()->FindBin(zRange), measured->GetYaxis()->FindBin(-etaRange), measured->GetYaxis()->FindBin(etaRange), measured->GetZaxis()->FindBin(ptCut), measured->GetZaxis()->GetNbins());
+  Float_t tracksG = generated->Integral(generated->GetXaxis()->FindBin(-zRange), generated->GetXaxis()->FindBin(zRange), generated->GetYaxis()->FindBin(-etaRange), generated->GetYaxis()->FindBin(etaRange), generated->GetZaxis()->FindBin(ptCut), generated->GetZaxis()->GetNbins());
 
-  Printf("AliCorrection::PrintInfo: Values in |eta| < 0.8, |vtx-z| < 10 and pt > %.2f:", ptCut);
-
-  Float_t tracksM = measured->Integral(measured->GetXaxis()->FindBin(-9.9), measured->GetXaxis()->FindBin(9.9), measured->GetYaxis()->FindBin(-0.79), measured->GetYaxis()->FindBin(0.79), measured->GetZaxis()->FindBin(ptCut), measured->GetZaxis()->GetNbins());
-  Float_t tracksG = generated->Integral(generated->GetXaxis()->FindBin(-9.9), generated->GetXaxis()->FindBin(9.9), generated->GetYaxis()->FindBin(-0.79), generated->GetYaxis()->FindBin(0.79), generated->GetZaxis()->FindBin(ptCut), generated->GetZaxis()->GetNbins());
-
-  Float_t eventsM = measuredEvents->Integral(measuredEvents->GetXaxis()->FindBin(-9.9), measuredEvents->GetXaxis()->FindBin(9.9), 1, measuredEvents->GetNbinsY());
-  Float_t eventsG = generatedEvents->Integral(generatedEvents->GetXaxis()->FindBin(-9.9), generatedEvents->GetXaxis()->FindBin(9.9), 1, generatedEvents->GetNbinsY());
+  Float_t eventsM = measuredEvents->Integral(measuredEvents->GetXaxis()->FindBin(-zRange), measuredEvents->GetXaxis()->FindBin(zRange), 1, measuredEvents->GetNbinsY());
+  Float_t eventsG = generatedEvents->Integral(generatedEvents->GetXaxis()->FindBin(-zRange), generatedEvents->GetXaxis()->FindBin(zRange), 1, generatedEvents->GetNbinsY());
 
   Printf("tracks measured: %.1f tracks generated: %.1f, events measured: %.1f, events generated %.1f", tracksM, tracksG, eventsM, eventsG);
 
   if (tracksM > 0)
     Printf("Effective average correction factor for TRACKS: %.3f", tracksG / tracksM);
   if (eventsM > 0)
-  Printf("Effective average correction factor for EVENTS: %.3f", eventsG / eventsM);
+    Printf("Effective average correction factor for EVENTS: %.3f", eventsG / eventsM);
 
   if (eventsM > 0 && eventsG > 0)
   {
@@ -407,3 +407,22 @@ void AliCorrection::PrintInfo(Float_t ptCut)
     Printf("%.2f tracks/event measured, %.2f tracks/event after correction --> effective average correction factor is %.3f (pt cut %.2f GeV/c)", tracksM, tracksG, tracksG / tracksM, ptCut);
   }
 }
+
+//____________________________________________________________________
+void AliCorrection::PrintInfo(Float_t ptCut)
+{
+  // prints some stats
+
+  TH3* measured = GetTrackCorrection()->GetMeasuredHistogram();
+  TH3* generated = GetTrackCorrection()->GetGeneratedHistogram();
+
+  TH2* measuredEvents = GetEventCorrection()->GetMeasuredHistogram();
+  TH2* generatedEvents = GetEventCorrection()->GetGeneratedHistogram();
+
+  Printf("AliCorrection::PrintInfo: Whole phasespace:");
+
+  Printf("tracks measured: %.1f tracks generated: %.1f, events measured: %.1f, events generated %.1f", measured->Integral(), generated->Integral(), measuredEvents->Integral(), generatedEvents->Integral());
+
+  PrintStats(10, 1.0, ptCut);
+  PrintStats(10, 1.5, ptCut);
+}
index 845c48e..c1d8e1f 100644 (file)
@@ -53,6 +53,7 @@ public:
   virtual void ReduceInformation();
 
   virtual void Reset(Option_t* option = "");
+  void PrintStats(Float_t zRange, Float_t etaRange, Float_t ptCut);
   void PrintInfo(Float_t ptCut);
 
 protected:
index d88376e..ba01ff7 100644 (file)
@@ -47,7 +47,6 @@ Bool_t AliPWG0Helper::IsEventTriggered(ULong64_t triggerMask, Trigger trigger)
   //
   // this function needs the branch fTriggerMask
   
-
   // definitions from p-p.cfg
   ULong64_t spdFO = (1 << 14);
   ULong64_t v0left = (1 << 11);
@@ -67,49 +66,17 @@ Bool_t AliPWG0Helper::IsEventTriggered(ULong64_t triggerMask, Trigger trigger)
         return kTRUE;
       break;
     }
+    case kSPDFASTOR:
+    {
+      if (triggerMask & spdFO)
+        return kTRUE;
+      break;
+    }
   }
 
   return kFALSE;
 }
 
-//____________________________________________________________________
-Bool_t AliPWG0Helper::IsVertexReconstructed(const AliESD* aEsd)
-{
-  // see function with AliESDVertex argument
-
-  const AliESDVertex* vtxESD = aEsd->GetVertex();
-  return IsVertexReconstructed(vtxESD);
-}
-
-//____________________________________________________________________
-Bool_t AliPWG0Helper::IsVertexReconstructed(const AliESDVertex* vtxESD)
-{
-  // checks if the vertex is reasonable
-  //
-  // this function needs the branches fSPDVertex*
-
-  if (!vtxESD)
-    return kFALSE;
-
-  // the vertex should be reconstructed
-  if (strcmp(vtxESD->GetName(), "default")==0)
-    return kFALSE;
-
-  // check Ncontributors
-  if (vtxESD->GetNContributors() <= 0)
-    return kFALSE;
-
-  Double_t vtx_res[3];
-  vtx_res[0] = vtxESD->GetXRes();
-  vtx_res[1] = vtxESD->GetYRes();
-  vtx_res[2] = vtxESD->GetZRes();
-
-  if (vtx_res[2]==0 || vtx_res[2]>0.1)
-    return kFALSE;
-
-  return kTRUE;
-}
-
 //____________________________________________________________________
 const AliESDVertex* AliPWG0Helper::GetVertex(const AliESDEvent* aEsd, AnalysisMode analysisMode, Bool_t debug)
 {
@@ -197,8 +164,7 @@ Bool_t AliPWG0Helper::IsPrimaryCharged(TParticle* aParticle, Int_t aTotalPrimari
 
   if (strcmp(aParticle->GetName(),"XXX") == 0)
   {
-    if (adebug)
-      printf("WARNING: There is a particle named XXX.\n");
+    Printf("WARNING: There is a particle named XXX (pdg code %d).", pdgCode);
     return kFALSE;
   }
 
@@ -206,8 +172,7 @@ Bool_t AliPWG0Helper::IsPrimaryCharged(TParticle* aParticle, Int_t aTotalPrimari
 
   if (strcmp(pdgPart->ParticleClass(),"Unknown") == 0)
   {
-    if (adebug)
-      printf("WARNING: There is a particle with an unknown particle class (pdg code %d).\n", pdgCode);
+    Printf("WARNING: There is a particle with an unknown particle class (pdg code %d).", pdgCode);
     return kFALSE;
   }
 
@@ -332,7 +297,7 @@ const char* AliPWG0Helper::GetAxisTitle(TH3* hist, const char axis)
 }
 
 
-Int_t AliPWG0Helper::GetPythiaEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) {
+AliPWG0Helper::MCProcessType AliPWG0Helper::GetPythiaEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) {
 
   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader);
 
@@ -366,7 +331,7 @@ Int_t AliPWG0Helper::GetPythiaEventProcessType(AliGenEventHeader* aHeader, Bool_
 }
 
 
-Int_t AliPWG0Helper::GetDPMjetEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) {
+AliPWG0Helper::MCProcessType AliPWG0Helper::GetDPMjetEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) {
   //
   // get the process type of the event.
   //
@@ -402,7 +367,7 @@ Int_t AliPWG0Helper::GetDPMjetEventProcessType(AliGenEventHeader* aHeader, Bool_
 }
 
 
-Int_t AliPWG0Helper::GetEventProcessType(AliHeader* aHeader, Bool_t adebug) {
+AliPWG0Helper::MCProcessType AliPWG0Helper::GetEventProcessType(AliHeader* aHeader, Bool_t adebug) {
   //
   // get the process type of the event.
   //
@@ -502,45 +467,6 @@ Int_t AliPWG0Helper::FindPrimaryMotherLabel(AliStack* stack, Int_t label)
   return label;
 }
 
-void AliPWG0Helper::SetBranchStatusRecursive(TTree* tree, char *bname, Bool_t status, Bool_t debug)
-{
-  // Function  to switch on/off all data members of a top level branch
-  // this is needed for branches without a trailing dot ".", for those
-  // the root functionality with regular expressions does not work.
-  // Usage e.g.
-  // chain->SetBranchStatus("*", 0); 
-  // SetBranchStatusRecursive(chain,"SPDVertex",1);
-  // You need to give the full name of the top level branch zou want to access
-  //==========================================================================
-  // Author Christian.Klein-Boesing@cern.ch
-
-  if (!tree)
-    return;
-
-  TBranch *br = tree->GetBranch(bname);
-  if(!br) {
-    Printf("AliPWG0Helper::SetBranchStatusRecursive: Branch %s not found", bname);
-  }
-
-  TObjArray *leaves = tree->GetListOfLeaves();
-  Int_t nleaves = leaves->GetEntries();
-  TLeaf *leaf = 0;
-  TBranch *branch = 0;
-  TBranch *mother = 0;
-  for (Int_t i=0;i<nleaves;i++)  {
-    // the matched entry e.g. SPDVertex is its own Mother
-    leaf = (TLeaf*)leaves->UncheckedAt(i);
-    branch = (TBranch*)leaf->GetBranch();
-    mother = branch->GetMother();
-    if (mother==br){
-      if (debug)
-        Printf(">>>> AliPWG0Helper::SetBranchStatusRecursive: Setting status %s %s to %d", mother->GetName(), leaf->GetName(), status);
-      if (status) branch->ResetBit(kDoNotProcess);
-      else        branch->SetBit(kDoNotProcess);
-    }
-  }
-}
-
 //____________________________________________________________________
 void AliPWG0Helper::NormalizeToBinWidth(TH1* hist)
 {
@@ -594,6 +520,7 @@ void AliPWG0Helper::PrintConf(AnalysisMode analysisMode, Trigger trigger)
   {
     case kMB1 : str += "MB1"; break;
     case kMB2 : str += "MB2"; break;
+    case kSPDFASTOR : str += "SPD FASTOR"; break;
   }
 
   str += " <<<<";
index b1a1e1f..4439ee6 100644 (file)
@@ -22,22 +22,20 @@ class TTree;
 class AliPWG0Helper : public TObject
 {
   public:
-    enum Trigger { kMB1 = 0, kMB2 }; // definition from ALICE-INT-2005-025
+    enum Trigger { kMB1 = 0, kMB2, kSPDFASTOR }; // definition from ALICE-INT-2005-025
     enum AnalysisMode { kInvalid = -1, kSPD = 0, kTPC, kTPCITS };
     // in case we want to use bitmaps...
-    enum MCProcessType { kInvalidProcess = -1, kND = 0x1, kDD = 0x2, kSD = 0x4 }; 
+    enum MCProcessType { kInvalidProcess = -1, kND = 0x1, kDD = 0x2, kSD = 0x4 };
 
     static Bool_t IsEventTriggered(const AliESD* aEsd, Trigger trigger = kMB2);
     static Bool_t IsEventTriggered(ULong64_t triggerMask, Trigger trigger = kMB2);
-    static Bool_t IsVertexReconstructed(const AliESD* aEsd);
-    static Bool_t IsVertexReconstructed(const AliESDVertex* vtxESD);
     static const AliESDVertex* GetVertex(const AliESDEvent* aEsd, AnalysisMode analysisMethod, Bool_t debug = kFALSE);
 
     static Bool_t IsPrimaryCharged(TParticle* aParticle, Int_t aTotalPrimaries, Bool_t adebug = kFALSE);
 
-    static Int_t GetEventProcessType(AliHeader* aHeader, Bool_t adebug = kFALSE);
-    static Int_t GetPythiaEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug = kFALSE);
-    static Int_t GetDPMjetEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug = kFALSE);
+    static AliPWG0Helper::MCProcessType GetEventProcessType(AliHeader* aHeader, Bool_t adebug = kFALSE);
+    static AliPWG0Helper::MCProcessType GetPythiaEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug = kFALSE);
+    static AliPWG0Helper::MCProcessType GetDPMjetEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug = kFALSE);
     static Int_t GetLastProcessType() { return fgLastProcessType; }
 
     static TParticle* FindPrimaryMother(AliStack* stack, Int_t label);
@@ -47,8 +45,6 @@ class AliPWG0Helper : public TObject
     static void CreateDividedProjections(TH3* hist, TH3* hist2, const char* axis = 0, Bool_t putErrors = kFALSE, Bool_t save = kFALSE);
     static const char* GetAxisTitle(TH3* hist, const char axis);
 
-    static void SetBranchStatusRecursive(TTree* tree, char *bname, Bool_t status, Bool_t debug = kFALSE);
-
     static void NormalizeToBinWidth(TH1* hist);
     static void NormalizeToBinWidth(TH2* hist);
 
index c801069..4231a49 100644 (file)
@@ -15,17 +15,17 @@ AliESDtrackCuts* CreateTrackCuts(AliPWG0Helper::AnalysisMode analysisMode, Bool_
   Double_t cov3 = 0.5;
   Double_t cov4 = 0.5;
   Double_t cov5 = 2;
-  Double_t nSigma = 3;
+  Double_t nSigma = 4;
+
+  Bool_t tpcRefit = kTRUE;
 
   TString tag("Global tracking");
 
   // TPC-only cuts
-  if (analysisMode == AliPWG0Helper::kTPC) 
+  if (analysisMode == AliPWG0Helper::kTPC)
   {
-    // beta cuts (still under investigation)
-    cov1 = 4;
-    cov2 = 4;
-    nSigma = 4;
+    // eventually replace by kTPCin?
+    tpcRefit = kFALSE;
     
     tag = "TPC-only tracking";
   }
@@ -41,8 +41,8 @@ AliESDtrackCuts* CreateTrackCuts(AliPWG0Helper::AnalysisMode analysisMode, Bool_
 
   esdTrackCuts->SetMinNsigmaToVertex(nSigma);
   esdTrackCuts->SetRequireSigmaToVertex(kTRUE);
-
-  esdTrackCuts->SetRequireTPCRefit(kTRUE);
+  
+  esdTrackCuts->SetRequireTPCRefit(tpcRefit);
   esdTrackCuts->SetAcceptKingDaughters(kFALSE);
 
   esdTrackCuts->SetMinNClustersTPC(50);
diff --git a/PWG0/dNdEta/AliMultiplicityESDSelector.cxx b/PWG0/dNdEta/AliMultiplicityESDSelector.cxx
deleted file mode 100644 (file)
index 7f26a66..0000000
+++ /dev/null
@@ -1,275 +0,0 @@
-/* $Id$ */
-
-#include "AliMultiplicityESDSelector.h"
-
-#include <TVector3.h>
-#include <TFile.h>
-#include <TH2F.h>
-#include <TH3F.h>
-#include <TTree.h>
-
-#include <AliLog.h>
-#include <AliESD.h>
-#include <AliMultiplicity.h>
-
-#include "esdTrackCuts/AliESDtrackCuts.h"
-#include "AliPWG0Helper.h"
-#include "dNdEta/AliMultiplicityCorrection.h"
-
-//#define TPCMEASUREMENT
-#define ITSMEASUREMENT
-
-ClassImp(AliMultiplicityESDSelector)
-
-AliMultiplicityESDSelector::AliMultiplicityESDSelector() :
-  AliSelector(),
-  fMultiplicity(0),
-  fEsdTrackCuts(0)
-{
-  //
-  // Constructor. Initialization of pointers
-  //
-}
-
-AliMultiplicityESDSelector::~AliMultiplicityESDSelector()
-{
-  //
-  // Destructor
-  //
-
-  // histograms are in the output list and deleted when the output
-  // list is deleted by the TSelector dtor
-}
-
-void AliMultiplicityESDSelector::Begin(TTree* tree)
-{
-  // Begin function
-
-  AliSelector::Begin(tree);
-
-  ReadUserObjects(tree);
-}
-
-void AliMultiplicityESDSelector::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 AliMultiplicityESDSelector::SlaveBegin(TTree* tree)
-{
-  // The SlaveBegin() function is called after the Begin() function.
-  // When running with PROOF SlaveBegin() is called on each slave server.
-  // The tree argument is deprecated (on PROOF 0 is passed).
-
-  AliSelector::SlaveBegin(tree);
-
-  ReadUserObjects(tree);
-
-  fMultiplicity = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
-}
-
-void AliMultiplicityESDSelector::Init(TTree* tree)
-{
-  // read the user objects
-
-  AliSelector::Init(tree);
-
-  // enable only the needed branches
-  if (tree)
-  {
-    tree->SetBranchStatus("*", 0);
-    tree->SetBranchStatus("fTriggerMask", 1);
-    tree->SetBranchStatus("fSPDVertex*", 1);
-
-    #ifdef ITSMEASUREMENT
-      tree->SetBranchStatus("fSPDMult*", 1);
-    #endif
-
-    #ifdef TPCMEASUREMENT
-      AliESDtrackCuts::EnableNeededBranches(tree);
-    #endif
-  }
-}
-
-Bool_t AliMultiplicityESDSelector::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).
-
-  if (AliSelector::Process(entry) == kFALSE)
-    return kFALSE;
-
-  // Check prerequisites
-  if (!fESD)
-  {
-    AliDebug(AliLog::kError, "ESD branch not available");
-    return kFALSE;
-  }
-
-  Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD);
-  Bool_t eventVertex = AliPWG0Helper::IsVertexReconstructed(fESD);
-
-  if (!eventTriggered || !eventVertex)
-    return kTRUE;
-
-  // get the ESD vertex
-  const AliESDVertex* vtxESD = fESD->GetVertex();
-  Double_t vtx[3];
-  vtxESD->GetXYZ(vtx);
-
-  Int_t nESDTracks05 = 0;
-  Int_t nESDTracks10 = 0;
-  Int_t nESDTracks15 = 0;
-  Int_t nESDTracks20 = 0;
-
-#ifdef ITSMEASUREMENT
-  // get tracklets
-  const AliMultiplicity* mult = fESD->GetMultiplicity();
-  if (!mult)
-  {
-    AliDebug(AliLog::kError, "AliMultiplicity not available");
-    return kFALSE;
-  }
-
-  // get multiplicity from ITS tracklets
-  for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
-  {
-    //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
-
-    // this removes non-tracklets. Very bad solution. SPD guys are working on better solution...
-    if (mult->GetDeltaPhi(i) < -1000)
-      continue;
-
-    Float_t theta = mult->GetTheta(i);
-    Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
-
-    if (TMath::Abs(eta) < 0.5)
-      nESDTracks05++;
-
-    if (TMath::Abs(eta) < 1.0)
-      nESDTracks10++;
-
-    if (TMath::Abs(eta) < 1.5)
-      nESDTracks15++;
-
-    if (TMath::Abs(eta) < 2.0)
-      nESDTracks20++;
-  }
-#endif
-
-#ifdef TPCMEASUREMENT
-  if (!fEsdTrackCuts)
-  {
-    AliDebug(AliLog::kError, "fESDTrackCuts not available");
-    return kFALSE;
-  }
-
-  // get multiplicity from ESD tracks
-  TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
-  Int_t nGoodTracks = list->GetEntries();
-  // loop over esd tracks
-  for (Int_t i=0; i<nGoodTracks; i++)
-  {
-    AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
-    if (!esdTrack)
-    {
-      AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
-      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);
-
-    Float_t theta = vector.Theta();
-    Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
-    Float_t pt = vector.Pt();
-
-    //if (pt < kPtCut)
-    //  continue;
-
-    if (TMath::Abs(eta) < 0.5)
-      nESDTracks05++;
-
-    if (TMath::Abs(eta) < 1.0)
-      nESDTracks10++;
-
-    if (TMath::Abs(eta) < 1.5)
-      nESDTracks15++;
-
-    if (TMath::Abs(eta) < 2.0)
-      nESDTracks20++;
-  }
-#endif
-
-  fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks10, nESDTracks15, nESDTracks20);
-
-  return kTRUE;
-}
-
-void AliMultiplicityESDSelector::SlaveTerminate()
-{
-  // The SlaveTerminate() function is called after all entries or objects
-  // have been processed. When running with PROOF SlaveTerminate() is called
-  // on each slave server.
-
-  AliSelector::SlaveTerminate();
-
-  // Add the histograms to the output on each slave server
-  if (!fOutput)
-  {
-    AliDebug(AliLog::kError, Form("ERROR: Output list not initialized."));
-    return;
-  }
-
-  fOutput->Add(fMultiplicity);
-}
-
-void AliMultiplicityESDSelector::Terminate()
-{
-  // The Terminate() function is the last function to be called during
-  // 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();
-
-  fMultiplicity = dynamic_cast<AliMultiplicityCorrection*> (fOutput->FindObject("Multiplicity"));
-
-  if (!fMultiplicity)
-  {
-    AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p", (void*) fMultiplicity));
-    return;
-  }
-
-  TFile* file = TFile::Open("multiplicityESD.root", "RECREATE");
-
-  fMultiplicity->SaveHistograms();
-
-  file->Close();
-
-  fMultiplicity->DrawHistograms();
-}
diff --git a/PWG0/dNdEta/AliMultiplicityESDSelector.h b/PWG0/dNdEta/AliMultiplicityESDSelector.h
deleted file mode 100644 (file)
index 6487fef..0000000
+++ /dev/null
@@ -1,36 +0,0 @@
-/* $Id$ */
-
-#ifndef ALIMULTIPLICITYESDSELECTOR_H
-#define ALIMULTIPLICITYESDSELECTOR_H
-
-#include "AliSelector.h"
-
-class AliESDtrackCuts;
-class AliMultiplicityCorrection;
-
-class AliMultiplicityESDSelector : public AliSelector {
-  public:
-    AliMultiplicityESDSelector();
-    virtual ~AliMultiplicityESDSelector();
-
-    virtual void    Begin(TTree* tree);
-    virtual void    SlaveBegin(TTree *tree);
-    virtual void    Init(TTree *tree);
-    virtual Bool_t  Process(Long64_t entry);
-    virtual void    SlaveTerminate();
-    virtual void    Terminate();
-
- protected:
-    void ReadUserObjects(TTree* tree);
-
-    AliMultiplicityCorrection* fMultiplicity; // object containing the extracted data
-    AliESDtrackCuts* fEsdTrackCuts;           // Object containing the parameters of the esd track cuts
-
- private:
-    AliMultiplicityESDSelector(const AliMultiplicityESDSelector&);
-    AliMultiplicityESDSelector& operator=(const AliMultiplicityESDSelector&);
-
-  ClassDef(AliMultiplicityESDSelector, 0);
-};
-
-#endif
diff --git a/PWG0/dNdEta/AliMultiplicityMCSelector.cxx b/PWG0/dNdEta/AliMultiplicityMCSelector.cxx
deleted file mode 100644 (file)
index c5b8a8d..0000000
+++ /dev/null
@@ -1,667 +0,0 @@
-/* $Id$ */
-
-#include "AliMultiplicityMCSelector.h"
-
-#include <TStyle.h>
-#include <TSystem.h>
-#include <TCanvas.h>
-#include <TVector3.h>
-#include <TChain.h>
-#include <TFile.h>
-#include <TH2F.h>
-#include <TH3F.h>
-#include <TParticle.h>
-#include <TRandom.h>
-#include <TNtuple.h>
-#include <TObjString.h>
-#include <TF1.h>
-
-#include <AliLog.h>
-#include <AliESD.h>
-#include <AliStack.h>
-#include <AliHeader.h>
-#include <AliGenEventHeader.h>
-#include <AliMultiplicity.h>
-
-#include "esdTrackCuts/AliESDtrackCuts.h"
-#include "AliPWG0Helper.h"
-#include "dNdEta/AliMultiplicityCorrection.h"
-#include "AliCorrection.h"
-#include "AliCorrectionMatrix3D.h"
-
-#define TPCMEASUREMENT
-//#define ITSMEASUREMENT
-
-ClassImp(AliMultiplicityMCSelector)
-
-AliMultiplicityMCSelector::AliMultiplicityMCSelector() :
-  AliSelectorRL(),
-  fMultiplicity(0),
-  fEsdTrackCuts(0),
-  fSystSkipParticles(kFALSE),
-  fSelectProcessType(0),
-  fParticleSpecies(0),
-  fPtSpectrum(0)
-{
-  //
-  // Constructor. Initialization of pointers
-  //
-
-  for (Int_t i = 0; i<4; i++)
-    fParticleCorrection[i] = 0;
-}
-
-AliMultiplicityMCSelector::~AliMultiplicityMCSelector()
-{
-  //
-  // Destructor
-  //
-
-  // histograms are in the output list and deleted when the output
-  // list is deleted by the TSelector dtor
-}
-
-void AliMultiplicityMCSelector::Begin(TTree* tree)
-{
-  // Begin function
-
-  AliSelectorRL::Begin(tree);
-
-  ReadUserObjects(tree);
-}
-
-void AliMultiplicityMCSelector::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.");
-
-  if (!fPtSpectrum && fInput)
-    fPtSpectrum = dynamic_cast<TH1*> (fInput->FindObject("pt-spectrum"));
-
-  if (!fPtSpectrum && tree)
-    fPtSpectrum = dynamic_cast<TH1*> (tree->GetUserInfo()->FindObject("pt-spectrum"));
-}
-
-void AliMultiplicityMCSelector::SlaveBegin(TTree* tree)
-{
-  // The SlaveBegin() function is called after the Begin() function.
-  // When running with PROOF SlaveBegin() is called on each slave server.
-  // The tree argument is deprecated (on PROOF 0 is passed).
-
-  AliSelector::SlaveBegin(tree);
-
-  ReadUserObjects(tree);
-
-  fMultiplicity = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
-
-  TString option(GetOption());
-  if (option.Contains("skip-particles"))
-  {
-    fSystSkipParticles = kTRUE;
-    AliInfo("WARNING: Systematic study enabled. Particles will be skipped.");
-  }
-
-  if (option.Contains("particle-efficiency"))
-    for (Int_t i = 0; i<4; i++)
-      fParticleCorrection[i] = new AliCorrection(Form("correction_%d", i), Form("correction_%d", i));
-
-  if (option.Contains("only-process-type-nd"))
-    fSelectProcessType = 1;
-
-  if (option.Contains("only-process-type-sd"))
-    fSelectProcessType = 2;
-
-  if (option.Contains("only-process-type-dd"))
-    fSelectProcessType = 3;
-
-  if (fSelectProcessType != 0)
-    AliInfo(Form("WARNING: Systematic study enabled. Only considering process type %d", fSelectProcessType));
-
-  if (option.Contains("pt-spectrum-hist"))
-  {
-    TFile* file = TFile::Open("ptspectrum_fit.root");
-    if (file)
-    {
-      TString subStr(option(option.Index("pt-spectrum")+17, 3));
-      TString histName(Form("ptspectrum_%s", subStr.Data()));
-      AliInfo(Form("Pt-Spectrum modification. Using %s.", histName.Data()));
-      fPtSpectrum = (TH1*) file->Get(histName);
-      if (!fPtSpectrum)
-       AliError("Histogram not found");
-    }
-    else
-      AliError("Could not open ptspectrum_fit.root. Pt Spectrum study could not be enabled.");
-
-    if (fPtSpectrum)
-      AliInfo("WARNING: Systematic study enabled. Pt spectrum will be modified");
-  }
-
-  if (option.Contains("pt-spectrum-func"))
-  {
-    if (fPtSpectrum)
-    {
-      Printf("Using function from input list for systematic p_t study");
-    }
-    else
-    {
-      fPtSpectrum = new TH1F("fPtSpectrum", "fPtSpectrum", 1, 0, 100);
-      fPtSpectrum->SetBinContent(1, 1);
-    }
-
-    if (fPtSpectrum)
-      AliInfo("WARNING: Systematic study enabled. Pt spectrum will be modified");
-  }
-
-  if (option.Contains("particle-species"))
-    fParticleSpecies = new TNtuple("fParticleSpecies", "fParticleSpecies", "vtx:Pi_True:K_True:p_True:o_True:Pi_Rec:K_Rec:p_Rec:o_Rec:nolabel:doublePrim:doubleCount");
-
-  // TODO set seed for random generator
-}
-
-void AliMultiplicityMCSelector::Init(TTree* tree)
-{
-  // read the user objects
-
-  AliSelectorRL::Init(tree);
-
-  // enable only the needed branches
-  if (tree)
-  {
-    tree->SetBranchStatus("*", 0);
-    tree->SetBranchStatus("fTriggerMask", 1);
-    tree->SetBranchStatus("fSPDVertex*", 1);
-
-    #ifdef ITSMEASUREMENT
-      tree->SetBranchStatus("fSPDMult*", 1);
-    #endif
-
-    #ifdef TPCMEASUREMENT
-      AliESDtrackCuts::EnableNeededBranches(tree);
-      tree->SetBranchStatus("fTracks.fLabel", 1);
-    #endif
-
-    tree->SetCacheSize(0);
-  }
-}
-
-Bool_t AliMultiplicityMCSelector::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).
-
-  if (AliSelectorRL::Process(entry) == kFALSE)
-    return kFALSE;
-
-  // Check prerequisites
-  if (!fESD)
-  {
-    AliDebug(AliLog::kError, "ESD branch not available");
-    return kFALSE;
-  }
-
-  if (!fEsdTrackCuts)
-  {
-    AliDebug(AliLog::kError, "fESDTrackCuts not available");
-    return kFALSE;
-  }
-
-  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;
-  }
-
-  if (fSelectProcessType > 0)
-  {
-    // getting process information
-    Int_t processtype = AliPWG0Helper::GetEventProcessType(header);
-
-    Bool_t processEvent = kFALSE;
-
-    // non diffractive
-    if (fSelectProcessType == 1 && processtype == AliPWG0Helper::kND )
-      processEvent = kTRUE;
-
-    // single diffractive
-    if (fSelectProcessType == 2 && processtype == AliPWG0Helper::kSD )
-      processEvent = kTRUE;
-
-    // double diffractive
-    if (fSelectProcessType == 3 && processtype == AliPWG0Helper::kDD )
-      processEvent = kTRUE;
-
-    if (!processEvent)
-    {
-      AliDebug(AliLog::kDebug, Form("Skipping this event, because it is not of the requested process type (%d)", processtype));
-      return kTRUE;
-    }
-  }
-
-  Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD);
-  eventTriggered = kTRUE; // no generated trigger for the simulated events
-  Bool_t eventVertex = AliPWG0Helper::IsVertexReconstructed(fESD);
-
-  // get the ESD vertex
-  const AliESDVertex* vtxESD = fESD->GetVertex();
-  Double_t vtx[3];
-  vtxESD->GetXYZ(vtx);
-
-  // get the MC vertex
-  AliGenEventHeader* genHeader = header->GenEventHeader();
-  TArrayF vtxMC(3);
-  genHeader->PrimaryVertex(vtxMC);
-
-  // get tracklets
-  const AliMultiplicity* mult = fESD->GetMultiplicity();
-  if (!mult)
-  {
-    AliDebug(AliLog::kError, "AliMultiplicity not available");
-    return kFALSE;
-  }
-
-  //const Float_t kPtCut = 0.3;
-
-  // get number of tracks from MC
-  // loop over mc particles
-  Int_t nPrim  = stack->GetNprimary();
-  Int_t nMCPart = stack->GetNtrack();
-
-  // tracks in different eta ranges
-  Int_t nMCTracks05 = 0;
-  Int_t nMCTracks09 = 0;
-  Int_t nMCTracks15 = 0;
-  Int_t nMCTracks20 = 0;
-  Int_t nMCTracksAll = 0;
-
-  // tracks per particle species, in |eta| < 2 (systematic study)
-  Int_t nMCTracksSpecies[4]; // (pi, K, p, other)
-  for (Int_t i = 0; i<4; ++i)
-    nMCTracksSpecies[i] = 0;
-  // eta range for nMCTracksSpecies, nESDTracksSpecies
-#ifdef ITSMEASUREMENT
-  const Float_t etaRange = 2.0;
-#else
-  const Float_t etaRange = 0.9;
-#endif
-
-  for (Int_t iMc = 0; iMc < nPrim; ++iMc)
-  {
-    AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc));
-
-    TParticle* particle = stack->Particle(iMc);
-
-    if (!particle)
-    {
-      AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
-      continue;
-    }
-
-    Bool_t debug = kFALSE;
-    if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim, debug) == kFALSE)
-    {
-      //printf("%d) DROPPED ", iMC);
-      //particle->Print();
-      continue;
-    }
-
-    //printf("%d) OK      ", iMC);
-    //particle->Print();
-
-    //if (particle->Pt() < kPtCut)
-    //  continue;
-
-    Int_t particleWeight = 1;
-
-    Float_t pt = particle->Pt();
-
-    // in case of systematic study, weight according to the change of the pt spectrum
-    // it cannot be just multiplied because we cannot count "half of a particle"
-    // instead a random generator decides if the particle is counted twice (if value > 1)
-    // or not (if value < 0)
-    if (fPtSpectrum)
-    {
-      Int_t bin = fPtSpectrum->FindBin(pt);
-      if (bin > 0 && bin <= fPtSpectrum->GetNbinsX())
-      {
-        Float_t factor = fPtSpectrum->GetBinContent(bin);
-        if (factor > 0)
-        {
-          Float_t random = gRandom->Uniform();
-          if (factor > 1 && random < factor - 1)
-          {
-            particleWeight = 2;
-          }
-          else if (factor < 1 && random < 1 - factor)
-            particleWeight = 0;
-        }
-      }
-    }
-
-    //Printf("MC weight is: %d", particleWeight);
-
-    if (TMath::Abs(particle->Eta()) < 0.5)
-      nMCTracks05 += particleWeight;
-
-    if (TMath::Abs(particle->Eta()) < 0.9)
-      nMCTracks09 += particleWeight;
-
-    if (TMath::Abs(particle->Eta()) < 1.5)
-      nMCTracks15 += particleWeight;
-
-    if (TMath::Abs(particle->Eta()) < 2.0)
-      nMCTracks20 += particleWeight;
-
-    nMCTracksAll += particleWeight;
-
-    if (fParticleCorrection[0] || fParticleSpecies)
-    {
-      Int_t id = -1;
-      switch (TMath::Abs(particle->GetPdgCode()))
-      {
-        case 211: id = 0; break;
-        case 321: id = 1; break;
-        case 2212: id = 2; break;
-        default: id = 3; break;
-      }
-      if (TMath::Abs(particle->Eta()) < etaRange)
-        nMCTracksSpecies[id]++;
-      if (fParticleCorrection[id])
-        fParticleCorrection[id]->GetTrackCorrection()->FillGene(vtxMC[2], particle->Eta(), particle->Pt());
-    }
-  }// end of mc particle
-
-  fMultiplicity->FillGenerated(vtxMC[2], eventTriggered, eventVertex, (Int_t) nMCTracks05, (Int_t) nMCTracks09, (Int_t) nMCTracks15, (Int_t) nMCTracks20, (Int_t) nMCTracksAll);
-
-  if (!eventTriggered || !eventVertex)
-    return kTRUE;
-
-  Int_t nESDTracks05 = 0;
-  Int_t nESDTracks09 = 0;
-  Int_t nESDTracks15 = 0;
-  Int_t nESDTracks20 = 0;
-
-  // tracks per particle species, in |eta| < 2 (systematic study)
-  Int_t nESDTracksSpecies[7]; // (pi, K, p, other, nolabel, doublecount_prim, doublecount_all)
-  for (Int_t i = 0; i<7; ++i)
-    nESDTracksSpecies[i] = 0;
-
-  Bool_t* foundPrimaries = new Bool_t[nPrim];   // to prevent double counting
-  for (Int_t i=0; i<nPrim; i++)
-    foundPrimaries[i] = kFALSE;
-
-  Bool_t* foundTracks = new Bool_t[nMCPart];    // to prevent double counting
-  for (Int_t i=0; i<nMCPart; i++)
-    foundTracks[i] = kFALSE;
-
-#ifdef ITSMEASUREMENT
-  // get multiplicity from ITS tracklets
-  for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
-  {
-    //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
-
-    // This removes non-tracklets in PDC06 data. Very bad solution. New solution is implemented for newer data. Keeping this for compatibility.
-    if (mult->GetDeltaPhi(i) < -1000)
-      continue;
-
-    // systematic study: 10% lower efficiency
-    if (fSystSkipParticles && (gRandom->Uniform() < 0.1))
-      continue;
-
-    Float_t theta = mult->GetTheta(i);
-    Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
-
-    // TODO define needed, because this only works with new AliRoot
-    Int_t label = mult->GetLabel(i);
-#endif
-#ifdef TPCMEASUREMENT
-  // get multiplicity from ESD tracks
-  TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
-  Int_t nGoodTracks = list->GetEntries();
-  // loop over esd tracks
-  for (Int_t i=0; i<nGoodTracks; i++)
-  {
-    AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
-    if (!esdTrack)
-    {
-      AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
-      continue;
-    }
-
-    Double_t p[3];
-    esdTrack->GetConstrainedPxPyPz(p); 
-    TVector3 vector(p);
-
-    Float_t theta = vector.Theta();
-    Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
-
-    Int_t label = esdTrack->GetLabel();
-
-    //Float_t pt = vector.Pt();
-    //if (pt < kPtCut)
-    //  continue;
-#endif
-
-    Int_t particleWeight = 1;
-
-    // in case of systematic study, weight according to the change of the pt spectrum
-    if (fPtSpectrum)
-    {
-      TParticle* mother = 0;
-
-      // preserve label for later
-      Int_t labelCopy = label;
-      if (labelCopy >= 0)
-        labelCopy = AliPWG0Helper::FindPrimaryMotherLabel(stack, labelCopy);
-      if (labelCopy >= 0)
-        mother = stack->Particle(labelCopy);
-
-      // in case of pt study we do not count particles w/o label, because they cannot be scaled
-      if (!mother)
-        continue;
-
-      // it cannot be just multiplied because we cannot count "half of a particle"
-      // instead a random generator decides if the particle is counted twice (if value > 1) 
-      // or not (if value < 0)
-      Int_t bin = fPtSpectrum->FindBin(mother->Pt());
-      if (bin > 0 && bin <= fPtSpectrum->GetNbinsX())
-      {
-        Float_t factor = fPtSpectrum->GetBinContent(bin);
-        if (factor > 0)
-        {
-          Float_t random = gRandom->Uniform();
-          if (factor > 1 && random < factor - 1)
-          {
-            particleWeight = 2;
-          }
-          else if (factor < 1 && random < 1 - factor)
-            particleWeight = 0;
-        }
-      }
-    }
-
-    //Printf("ESD weight is: %d", particleWeight);
-
-    if (TMath::Abs(eta) < 0.5)
-      nESDTracks05 += particleWeight;
-
-    if (TMath::Abs(eta) < 0.9)
-      nESDTracks09 += particleWeight;
-
-    if (TMath::Abs(eta) < 1.5)
-      nESDTracks15 += particleWeight;
-
-    if (TMath::Abs(eta) < 2.0)
-      nESDTracks20 += particleWeight;
-
-
-    if (fParticleCorrection[0] || fParticleSpecies)
-    {
-      Int_t motherLabel = -1;
-      TParticle* mother = 0;
-
-      // find mother
-      if (label >= 0)
-        motherLabel = AliPWG0Helper::FindPrimaryMotherLabel(stack, label);
-      if (motherLabel >= 0)
-        mother = stack->Particle(motherLabel);
-
-      if (!mother)
-      {
-        // count tracks that did not have a label
-        if (TMath::Abs(eta) < etaRange)
-          nESDTracksSpecies[4]++;
-        continue;
-      }
-
-      // get particle type (pion, proton, kaon, other)
-      Int_t id = -1;
-      switch (TMath::Abs(mother->GetPdgCode()))
-      {
-        case 211: id = 0; break;
-        case 321: id = 1; break;
-        case 2212: id = 2; break;
-        default: id = 3; break;
-      }
-
-      // double counting is ok for particle ratio study
-      if (TMath::Abs(eta) < etaRange)
-        nESDTracksSpecies[id]++;
-
-      // double counting is not ok for efficiency study
-
-      // check if we already counted this particle, this way distinguishes double counted particles (bug/artefact in tracking) or double counted primaries due to secondaries (physics)
-      if (foundTracks[label])
-      {
-        if (TMath::Abs(eta) < etaRange)
-          nESDTracksSpecies[6]++;
-        continue;
-      }
-      foundTracks[label] = kTRUE;
-
-      // particle (primary) already counted?
-      if (foundPrimaries[motherLabel])
-      {
-        if (TMath::Abs(eta) < etaRange)
-          nESDTracksSpecies[5]++;
-        continue;
-      }
-      foundPrimaries[motherLabel] = kTRUE;
-
-      if (fParticleCorrection[id])
-        fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], mother->Eta(), mother->Pt());
-    }
-  }
-
-  delete[] foundTracks;
-  delete[] foundPrimaries;
-
-  if ((Int_t) nMCTracks20 == 0 && nESDTracks20 > 3)
-    printf("WARNING: Event %lld has %d generated and %d reconstructed...\n", entry, nMCTracks20, nESDTracks20);
-
-  //printf("%.2f generated and %.2f reconstructed\n", nMCTracks20, nESDTracks20);
-
-  //Printf("%f %f %f %f %f %f %f %f %f %f", vtxMC[2], nMCTracks05, nMCTracks09, nMCTracks15, nMCTracks20, nMCTracksAll, nESDTracks05, nESDTracks09, nESDTracks15, nESDTracks20);
-
-  //Printf("%f %f %f %f %f %f %f %f %f %f", vtxMC[2], nMCTracks05, nMCTracks09, nMCTracks15, nMCTracks20, nMCTracksAll, nESDTracks05, nESDTracks09, nESDTracks15, nESDTracks20);
-
-  fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks09, nESDTracks15, nESDTracks20);
-
-  // fill response matrix using vtxMC (best guess)
-  fMultiplicity->FillCorrection(vtxMC[2],  nMCTracks05,  nMCTracks09,  nMCTracks15,  nMCTracks20,  nMCTracksAll,  nESDTracks05,  nESDTracks09,  nESDTracks15,  nESDTracks20);
-
-  if (fParticleSpecies)
-    fParticleSpecies->Fill(vtxMC[2], nMCTracksSpecies[0], nMCTracksSpecies[1], nMCTracksSpecies[2], nMCTracksSpecies[3], nESDTracksSpecies[0], nESDTracksSpecies[1], nESDTracksSpecies[2], nESDTracksSpecies[3], nESDTracksSpecies[4], nESDTracksSpecies[5], nESDTracksSpecies[6]);
-
-  //Printf("%d = %d %d %d %d %d %d %d", (Int_t) nESDTracks20, nESDTracksSpecies[0], nESDTracksSpecies[1], nESDTracksSpecies[2], nESDTracksSpecies[3], nESDTracksSpecies[4], nESDTracksSpecies[5], nESDTracksSpecies[6]);
-
-  return kTRUE;
-}
-
-void AliMultiplicityMCSelector::SlaveTerminate()
-{
-  // The SlaveTerminate() function is called after all entries or objects
-  // have been processed. When running with PROOF SlaveTerminate() is called
-  // on each slave server.
-
-  AliSelector::SlaveTerminate();
-
-  // Add the histograms to the output on each slave server
-  if (!fOutput)
-  {
-    AliDebug(AliLog::kError, Form("ERROR: Output list not initialized."));
-    return;
-  }
-
-  fOutput->Add(fMultiplicity);
-  for (Int_t i = 0; i < 4; ++i)
-    fOutput->Add(fParticleCorrection[i]);
-
-  fOutput->Add(fParticleSpecies);
-}
-
-void AliMultiplicityMCSelector::Terminate()
-{
-  // The Terminate() function is the last function to be called during
-  // 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();
-
-  fMultiplicity = dynamic_cast<AliMultiplicityCorrection*> (fOutput->FindObject("Multiplicity"));
-  for (Int_t i = 0; i < 4; ++i)
-    fParticleCorrection[i] = dynamic_cast<AliCorrection*> (fOutput->FindObject(Form("correction_%d", i)));
-  fParticleSpecies = dynamic_cast<TNtuple*> (fOutput->FindObject("fParticleSpecies"));
-
-  if (!fMultiplicity)
-  {
-    AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p", (void*) fMultiplicity));
-    return;
-  }
-
-  TFile* file = TFile::Open("multiplicityMC.root", "RECREATE");
-
-  fMultiplicity->SaveHistograms();
-  for (Int_t i = 0; i < 4; ++i)
-    if (fParticleCorrection[i])
-      fParticleCorrection[i]->SaveHistograms();
-  if (fParticleSpecies)
-    fParticleSpecies->Write();
-
-  TObjString option(GetOption());
-  option.Write();
-
-  file->Close();
-
-  //fMultiplicity->DrawHistograms();
-}
diff --git a/PWG0/dNdEta/AliMultiplicityMCSelector.h b/PWG0/dNdEta/AliMultiplicityMCSelector.h
deleted file mode 100644 (file)
index 79c552f..0000000
+++ /dev/null
@@ -1,49 +0,0 @@
-/* $Id$ */
-
-#ifndef ALIMULTIPLICITYMCSELECTOR_H
-#define ALIMULTIPLICITYMCSELECTOR_H
-
-#include "AliSelectorRL.h"
-
-class AliESDtrackCuts;
-class AliMultiplicityCorrection;
-class AliCorrection;
-class TNtuple;
-class TH1;
-
-class AliMultiplicityMCSelector : public AliSelectorRL {
-  public:
-    AliMultiplicityMCSelector();
-    virtual ~AliMultiplicityMCSelector();
-
-    virtual void    Begin(TTree* tree);
-    virtual void    SlaveBegin(TTree *tree);
-    virtual void    Init(TTree *tree);
-    virtual Bool_t  Process(Long64_t entry);
-    virtual void    SlaveTerminate();
-    virtual void    Terminate();
-
- protected:
-    void ReadUserObjects(TTree* tree);
-
-    AliMultiplicityCorrection* fMultiplicity; // object containing the extracted data
-    AliESDtrackCuts* fEsdTrackCuts;           // Object containing the parameters of the esd track cuts
-
-    Bool_t fSystSkipParticles;     // if true skips particles (systematic study)
-    AliCorrection* fParticleCorrection[4]; // correction from measured to generated particles for trigger, vertex sample in |eta| < 2;
-                                           // for each of the species: pi, k, p, other; for systematic study of pt cut off
-    Int_t fSelectProcessType;        // 0 = all (default), 1 = ND, 2 = SD, 3 = DD (for systematic study)
-    TNtuple *fParticleSpecies;       // per event: vtx_mc, (pi, k, p, rest (in |eta| < 2)) X (true, recon) + (nolabel,
-                                     // doubleTracks, doublePrimaries) [doubleTracks + doublePrimaries are already part of
-                                     // rec. particles!)
-
-    TH1* fPtSpectrum;                // function that modifies the pt spectrum (syst. study)
-
- private:
-    AliMultiplicityMCSelector(const AliMultiplicityMCSelector&);
-    AliMultiplicityMCSelector& operator=(const AliMultiplicityMCSelector&);
-
-  ClassDef(AliMultiplicityMCSelector, 0);
-};
-
-#endif
diff --git a/PWG0/dNdEta/AliVertexSelector.cxx b/PWG0/dNdEta/AliVertexSelector.cxx
deleted file mode 100644 (file)
index ab16e5f..0000000
+++ /dev/null
@@ -1,277 +0,0 @@
-/**************************************************************************
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- *                                                                        *
- * Author: The ALICE Off-line Project.                                    *
- * Contributors are mentioned in the code where appropriate.              *
- *                                                                        *
- * Permission to use, copy, modify and distribute this software and its   *
- * documentation strictly for non-commercial purposes is hereby granted   *
- * without fee, provided that the above copyright notice appears in all   *
- * copies and that both the copyright notice and this permission notice   *
- * appear in the supporting documentation. The authors make no claims     *
- * about the suitability of this software for any purpose. It is          *
- * provided "as is" without express or implied warranty.                  *
- **************************************************************************/
-
-/* $Id$ */
-
-// The ESD is available as member fESD
-//
-// The Process function is nearly empty. Implement your analysis there and look at the other listed below functions you
-// might need.
-//
-// The following methods can be overrriden. Please do not forgot to call the base class function.
-//
-//    Begin():        called everytime a loop on the tree starts,
-//                    a convenient place to create your histograms.
-//    SlaveBegin():   called after Begin(), when on PROOF called only on the
-//                    slave servers.
-//    Init():         called for each new tree. Enable/Disable branches here.
-//    Process():      called for each event, in this function you decide what
-//                    to read and fill your histograms.
-//    SlaveTerminate: called at the end of the loop on the tree, when on PROOF
-//                    called only on the slave servers.
-//    Terminate():    called at the end of the loop on the tree,
-//                    a convenient place to draw/fit your histograms.
-//
-//  Author: Jan.Fiete.Grosse-Oetringhaus@cern.ch
-
-#include "AliVertexSelector.h"
-
-#include <TH1F.h>
-#include <TH2F.h>
-#include <TFile.h>
-#include <TTree.h>
-#include <TProfile.h>
-#include <TCanvas.h>
-#include <TAxis.h>
-
-#include <AliLog.h>
-#include <AliESD.h>
-#include <AliHeader.h>
-#include <AliGenEventHeader.h>
-
-#include <AliPWG0Helper.h>
-#include "esdTrackCuts/AliESDtrackCuts.h"
-
-ClassImp(AliVertexSelector)
-
-AliVertexSelector::AliVertexSelector() :
-  AliSelectorRL(),
-  fEsdTrackCuts(0),
-  fVertexMC(0),
-  fVertexESD(0),
-  fVertexCorr(0),
-  fVertexCorr2(0),
-  fFakes(0)
-{
-  //
-  // Constructor. Initialization of pointers
-  //
-}
-
-AliVertexSelector::~AliVertexSelector()
-{
-  //
-  // Destructor
-  //
-}
-
-void AliVertexSelector::SlaveBegin(TTree* tree)
-{
-  AliSelectorRL::SlaveBegin(tree);
-
-  fVertexMC = new TH1F("fVertexMC", "fVertexMC;MC vtx-z;Count", 501, -10, 10);
-  fVertexESD = new TH1F("fVertexESD", "fVertexESD;ESD vtx-z;Count", 501, -10, 10);
-  fVertexCorr = new TH2F("fVertexCorr", "fVertexCorr;MC vtx-z;ESD vtx-z - MC vtx-z", 40, -10, 10, 200, -1, 1);
-  fVertexCorr2 = new TH2F("fVertexCorr2", "fVertexCorr2;MC vtx-z;ESD vtx-z - MC vtx-z", 40, -10, 10, 200, -1, 1);
-  fFakes = new TH2F("fFakes", "fFakes;type;label", 5, 0, 5, 1001, -500.5, 500.5);
-
-  if (!fEsdTrackCuts && fInput)
-    fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fInput->FindObject("AliESDtrackCuts"));
-}
-
-void AliVertexSelector::Init(TTree* tree)
-{
-  // read the user objects
-
-  AliSelectorRL::Init(tree);
-
-  // Enable only the needed branches
-  if (tree)
-  {
-    tree->SetBranchStatus("*", 0);
-    tree->SetBranchStatus("fTriggerMask", 1);
-    tree->SetBranchStatus("fSPDVertex*", 1);
-    tree->SetBranchStatus("fTracks.fLabel", 1);
-
-    AliESDtrackCuts::EnableNeededBranches(tree);
-  }
-}
-
-Bool_t AliVertexSelector::Process(Long64_t entry)
-{
-  //
-  // Implement your analysis here. Do not forget to call the parent class Process by
-  // if (AliSelector::Process(entry) == kFALSE)
-  //   return kFALSE;
-  //
-
-  if (AliSelectorRL::Process(entry) == kFALSE)
-    return kFALSE;
-
-  // Check prerequisites
-  if (!fESD)
-  {
-    AliDebug(AliLog::kError, "ESD branch not available");
-    return kFALSE;
-  }
-
-  if (!fEsdTrackCuts)
-  {
-    AliDebug(AliLog::kError, "fESDTrackCuts 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));
-    return kTRUE;
-  }
-
-  if (AliPWG0Helper::IsVertexReconstructed(fESD) == kFALSE)
-  {
-    AliDebug(AliLog::kDebug+1, Form("Skipping event %d because its vertex was not reconstructed", (Int_t) entry));
-    return kTRUE;
-  }
-
-  AliHeader* header = GetHeader();
-  if (!header)
-  {
-    AliDebug(AliLog::kError, "Header not available");
-    return kFALSE;
-  }
-
-  // get the MC vertex
-  AliGenEventHeader* genHeader = header->GenEventHeader();
-
-  TArrayF vtxMC(3);
-  genHeader->PrimaryVertex(vtxMC);
-
-  // ########################################################
-  // get the ESD vertex
-  const AliESDVertex* vtxESD = fESD->GetVertex();
-  Double_t vtx[3];
-  vtxESD->GetXYZ(vtx);
-
-  fVertexMC->Fill(vtxMC[2]);
-  fVertexESD->Fill(vtx[2]);
-  fVertexCorr->Fill(vtxMC[2], vtx[2] - vtxMC[2]);
-
-  Float_t correctedVertex = vtx[2] + 2.951034e-03 + 6.859620e-04 * vtxMC[2];
-
-  fVertexCorr2->Fill(vtxMC[2], correctedVertex - vtxMC[2]);
-
-  const Int_t max = 10000;
-
-  Bool_t used[max];
-  for (Int_t i=0; i<max; ++i)
-    used[i] = kFALSE;
-
-  Int_t fakeTracks = 0;
-
-  Int_t nTracks = fESD->GetNumberOfTracks();
-  for (Int_t t=0; t<nTracks; t++)
-  {
-    AliESDtrack* esdTrack = fESD->GetTrack(t);
-
-    if (!esdTrack)
-      continue;
-
-
-    if (!fEsdTrackCuts->AcceptTrack(esdTrack))
-      continue;
-
-    UInt_t status = esdTrack->GetStatus();
-    if ((!status & AliESDtrack::kTPCrefit))
-      continue;
-
-    Int_t label = TMath::Abs(esdTrack->GetLabel());
-    if (label < max)
-    {
-      if (used[label])
-      {
-        fFakes->Fill(1.5, esdTrack->GetLabel());
-        fakeTracks++;
-      }
-
-      used[label] = kTRUE;
-    }
-
-    fFakes->Fill(0.5, esdTrack->GetLabel());
-  }
-
-  if (fakeTracks > 10)
-    printf("In event %lld we have %d fake tracks.\n", entry, fakeTracks);
-
-  return kTRUE;
-}
-
-void AliVertexSelector::SlaveTerminate()
-{
-  // The SlaveTerminate() function is called after all entries or objects
-  // have been processed. When running with PROOF SlaveTerminate() is called
-  // on each slave server.
-
-  AliSelectorRL::SlaveTerminate();
-
-  fOutput->Add(fVertexMC);
-  fOutput->Add(fVertexESD);
-  fOutput->Add(fVertexCorr);
-  fOutput->Add(fVertexCorr2);
-  fOutput->Add(fFakes);
-}
-
-void AliVertexSelector::Terminate()
-{
-  // The Terminate() function is the last function to be called during
-  // a query. It always runs on the client, it can be used to present
-  // the results graphically or save the results to file.
-
-  AliSelectorRL::Terminate();
-
-  fVertexMC = dynamic_cast<TH1F*> (fOutput->FindObject("fVertexMC"));
-  fVertexESD = dynamic_cast<TH1F*> (fOutput->FindObject("fVertexESD"));
-  fVertexCorr = dynamic_cast<TH2F*> (fOutput->FindObject("fVertexCorr"));
-  fVertexCorr2 = dynamic_cast<TH2F*> (fOutput->FindObject("fVertexCorr2"));
-  fFakes = dynamic_cast<TH2F*> (fOutput->FindObject("fFakes"));
-
-  if (!fVertexMC || !fVertexESD || !fVertexCorr || !fVertexCorr2 || !fFakes)
-    return;
-
-  TFile* file = TFile::Open("vertex.root", "RECREATE");
-  fVertexMC->Write();
-  fVertexESD->Write();
-  fVertexCorr->Write();
-  fVertexCorr2->Write();
-  fFakes->Write();
-  file->Close();
-
-  TCanvas* canvas = new TCanvas("canvas", "canvas", 1000, 1000);
-  canvas->Divide(2, 2);
-
-  canvas->cd(1);
-  fVertexCorr->DrawCopy("COLZ");
-
-  canvas->cd(2);
-  fVertexCorr->ProfileX()->DrawCopy();
-
-  canvas->cd(3);
-  fVertexCorr2->DrawCopy("COLZ");
-
-  canvas->cd(4);
-  fFakes->DrawCopy();
-
-  printf("%f tracks, %f with label > 0, %f with label == 0, %f with label < 0\n", fFakes->Integral(1, 1, 0, fFakes->GetYaxis()->GetNbins()), fFakes->Integral(1, 1, fFakes->GetYaxis()->FindBin(0.0), fFakes->GetYaxis()->GetNbins()), fFakes->GetBinContent(1, fFakes->GetYaxis()->FindBin(0.0)), fFakes->Integral(1, 1, 0, fFakes->GetYaxis()->FindBin(0.0)));
-  printf("%f tracks, %f with label > 0, %f with label == 0, %f with label < 0\n", fFakes->Integral(2, 2, 0, fFakes->GetYaxis()->GetNbins()), fFakes->Integral(2, 2, fFakes->GetYaxis()->FindBin(0.0), fFakes->GetYaxis()->GetNbins()), fFakes->GetBinContent(2, fFakes->GetYaxis()->FindBin(0.0)), fFakes->Integral(2, 2, 0, fFakes->GetYaxis()->FindBin(0.0)));
-}
diff --git a/PWG0/dNdEta/AliVertexSelector.h b/PWG0/dNdEta/AliVertexSelector.h
deleted file mode 100644 (file)
index e2cd3c9..0000000
+++ /dev/null
@@ -1,39 +0,0 @@
-/* $Id$ */
-
-#ifndef AliVertexSelector_H
-#define AliVertexSelector_H
-
-#include "AliSelectorRL.h"
-
-class TH1F;
-class TH2F;
-class AliESDtrackCuts;
-
-class AliVertexSelector : public AliSelectorRL {
-  public:
-    AliVertexSelector();
-    virtual ~AliVertexSelector();
-
-    virtual void    SlaveBegin(TTree *tree);
-    virtual void    Init(TTree *tree);
-    virtual Bool_t  Process(Long64_t entry);
-    virtual void    SlaveTerminate();
-    virtual void    Terminate();
-
- private:
-    AliVertexSelector(const AliVertexSelector&);
-    AliVertexSelector& operator=(const AliVertexSelector&);
-
-    AliESDtrackCuts*  fEsdTrackCuts;     // Object containing the parameters of the esd track cuts
-
-    TH1F* fVertexMC;   // MC vertex distribution
-    TH1F* fVertexESD;  // ESD vertex distribution
-    TH2F* fVertexCorr; // vertex correlation (esd vtx - mc vtx vs. mc vtx)
-    TH2F* fVertexCorr2; // vertex correlation (esd vtx - mc vtx vs. mc vtx)
-
-    TH2F* fFakes; // counts the fakes
-
-  ClassDef(AliVertexSelector, 0);
-};
-
-#endif
index 9966f14..e2b242f 100644 (file)
@@ -8,8 +8,10 @@
 #include <TFile.h>
 #include <TH1F.h>
 #include <TH2F.h>
+#include <TH3F.h>
 #include <TProfile.h>
 #include <TParticle.h>
+#include <TParticlePDG.h>
 
 #include <AliLog.h>
 #include <AliESDVertex.h>
 
 ClassImp(AlidNdEtaCorrectionTask)
 
+AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask() :
+  AliAnalysisTask(),
+  fESD(0),
+  fOutput(0),
+  fOption(),
+  fAnalysisMode(AliPWG0Helper::kTPC),
+  fTrigger(AliPWG0Helper::kMB1),
+  fSignMode(0),
+  fOnlyPrimaries(kFALSE),
+  fEsdTrackCuts(0),
+  fdNdEtaCorrection(0),
+  fdNdEtaAnalysisMC(0),
+  fdNdEtaAnalysisESD(0),
+  fPIDParticles(0),
+  fPIDTracks(0),
+  fVertexCorrelation(0),
+  fVertexProfile(0),
+  fVertexShift(0),
+  fVertexShiftNorm(0),
+  fEtaCorrelation(0),
+  fEtaProfile(0),
+  fEtaResolution(0),
+  fpTResolution(0),
+  fEsdTrackCutsPrim(0),
+  fEsdTrackCutsSec(0),
+  fTemp1(0),
+  fTemp2(0),
+  fMultAll(0),
+  fMultTr(0),
+  fMultVtx(0),
+  fEventStats(0)
+{
+  //
+  // Constructor. Initialization of pointers
+  //
+
+  for (Int_t i=0; i<2; i++)
+    fdNdEtaCorrectionProcessType[i] = 0;
+  
+  for (Int_t i=0; i<8; i++)
+    fDeltaPhi[i] = 0;
+}
+
 AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask(const char* opt) :
   AliAnalysisTask("AlidNdEtaCorrectionTask", ""),
   fESD(0),
@@ -50,11 +95,16 @@ AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask(const char* opt) :
   fPIDTracks(0),
   fVertexCorrelation(0),
   fVertexProfile(0),
+  fVertexShift(0),
   fVertexShiftNorm(0),
   fEtaCorrelation(0),
   fEtaProfile(0),
-  fSigmaVertexTracks(0),
-  fSigmaVertexPrim(0),
+  fEtaResolution(0),
+  fpTResolution(0),
+  fEsdTrackCutsPrim(0),
+  fEsdTrackCutsSec(0),
+  fTemp1(0),
+  fTemp2(0),
   fMultAll(0),
   fMultTr(0),
   fMultVtx(0),
@@ -149,10 +199,10 @@ void AlidNdEtaCorrectionTask::CreateOutputObjects()
   fdNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction", fAnalysisMode);
   fOutput->Add(fdNdEtaCorrection);
 
-  fPIDParticles = new TH1F("pid_particles", "PID of generated primary particles", 10001, -5000.5, 5000.5);
+  fPIDParticles = new TH1F("fPIDParticles", "PID of generated primary particles", 10001, -5000.5, 5000.5);
   fOutput->Add(fPIDParticles);
 
-  fPIDTracks = new TH1F("pid_tracks", "MC PID of reconstructed tracks", 10001, -5000.5, 5000.5);
+  fPIDTracks = new TH1F("fPIDTracks", "MC PID of reconstructed tracks", 10001, -5000.5, 5000.5);
   fOutput->Add(fPIDTracks);
 
   fdNdEtaAnalysisMC = new dNdEtaAnalysis("dndetaMC", "dndetaMC", fAnalysisMode);
@@ -161,6 +211,14 @@ void AlidNdEtaCorrectionTask::CreateOutputObjects()
   fdNdEtaAnalysisESD = new dNdEtaAnalysis("dndetaESD", "dndetaESD", fAnalysisMode);
   fOutput->Add(fdNdEtaAnalysisESD);
 
+  if (fEsdTrackCuts)
+  {
+    fEsdTrackCutsPrim = dynamic_cast<AliESDtrackCuts*> (fEsdTrackCuts->Clone("fEsdTrackCutsPrim"));
+    fEsdTrackCutsSec  = dynamic_cast<AliESDtrackCuts*> (fEsdTrackCuts->Clone("fEsdTrackCutsSec"));
+    fOutput->Add(fEsdTrackCutsPrim);
+    fOutput->Add(fEsdTrackCutsSec);
+  }
+
   if (fOption.Contains("process-types")) {
     fdNdEtaCorrectionProcessType[0] = new AlidNdEtaCorrection("dndeta_correction_ND", "dndeta_correction_ND", fAnalysisMode);
     fdNdEtaCorrectionProcessType[1] = new AlidNdEtaCorrection("dndeta_correction_SD", "dndeta_correction_SD", fAnalysisMode);
@@ -171,19 +229,19 @@ void AlidNdEtaCorrectionTask::CreateOutputObjects()
     fOutput->Add(fdNdEtaCorrectionProcessType[2]);
   }
 
-  if (fOption.Contains("sigma-vertex"))
-  {
-    fSigmaVertexTracks = new TH1F("fSigmaVertexTracks", "fSigmaVertexTracks;Nsigma2vertex;NacceptedTracks", 60, 0.05, 6.05);
-    fSigmaVertexPrim = new TH1F("fSigmaVertexPrim", "fSigmaVertexPrim;Nsigma2vertex;NacceptedPrimaries", 60, 0.05, 6.05);
-    fOutput->Add(fSigmaVertexTracks);
-    fOutput->Add(fSigmaVertexPrim);
-    Printf("WARNING: sigma-vertex analysis enabled. This will produce weird results in the AliESDtrackCuts histograms");
-  }
+  /*
+  fTemp1 = new TH3F("fTemp1", "fTemp1 prim;b0;b1;1/pT", 200, 0, 10, 200, 0, 10, 200, 0, 10);
+  fOutput->Add(fTemp1);
+  fTemp2 = new TH3F("fTemp2", "fTemp2 sec;b0;b1;1/pT", 200, 0, 10, 200, 0, 10, 200, 0, 10);
+  fOutput->Add(fTemp2);
+  */
 
   fVertexCorrelation = new TH2F("fVertexCorrelation", "fVertexCorrelation;MC z-vtx;ESD z-vtx", 120, -30, 30, 120, -30, 30);
   fOutput->Add(fVertexCorrelation);
   fVertexProfile = new TProfile("fVertexProfile", "fVertexProfile;MC z-vtx;MC z-vtx - ESD z-vtx", 40, -20, 20);
   fOutput->Add(fVertexProfile);
+  fVertexShift = new TH1F("fVertexShift", "fVertexShift;(MC z-vtx - ESD z-vtx);Entries", 201, -2, 2);
+  fOutput->Add(fVertexShift);
   fVertexShiftNorm = new TH1F("fVertexShiftNorm", "fVertexShiftNorm;(MC z-vtx - ESD z-vtx) / #sigma_{ESD z-vtx};Entries", 200, -100, 100);
   fOutput->Add(fVertexShiftNorm);
 
@@ -191,6 +249,11 @@ void AlidNdEtaCorrectionTask::CreateOutputObjects()
   fOutput->Add(fEtaCorrelation);
   fEtaProfile = new TProfile("fEtaProfile", "fEtaProfile;MC #eta;MC #eta - ESD #eta", 120, -3, 3);
   fOutput->Add(fEtaProfile);
+  fEtaResolution = new TH1F("fEtaResolution", "fEtaResolution;MC #eta - ESD #eta", 201, -0.2, 0.2);
+  fOutput->Add(fEtaResolution);
+
+  fpTResolution = new TH1F("fpTResolution", "fpTResolution;MC p_{T} - ESD p_{T}", 201, -0.2, 0.2);
+  fOutput->Add(fpTResolution);
 
   fMultAll = new TH1F("fMultAll", "fMultAll", 500, -0.5, 499.5);
   fOutput->Add(fMultAll);
@@ -216,6 +279,12 @@ void AlidNdEtaCorrectionTask::CreateOutputObjects()
   fEventStats->GetYaxis()->SetBinLabel(2, "trg");
   fEventStats->GetYaxis()->SetBinLabel(3, "vtx");
   fEventStats->GetYaxis()->SetBinLabel(4, "trgvtx");
+
+  if (fEsdTrackCuts)
+  {
+    fEsdTrackCuts->SetName("fEsdTrackCuts");
+    fOutput->Add(fEsdTrackCuts);
+  }
 }
 
 void AlidNdEtaCorrectionTask::Exec(Option_t*)
@@ -291,7 +360,8 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
     eventVertex = kTRUE;
     
     Double_t diff = vtxMC[2] - vtx[2];
-    if (vtxESD->GetZRes() > 0) 
+    fVertexShift->Fill(diff);
+    if (vtxESD->GetZRes() > 0)
       fVertexShiftNorm->Fill(diff / vtxESD->GetZRes());
   }
   else
@@ -328,7 +398,7 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
     ptArr = new Float_t[mult->GetNumberOfTracklets()];
     deltaPhiArr = new Float_t[mult->GetNumberOfTracklets()];
 
-    // get multiplicity from ITS tracklets
+    // get multiplicity from SPD tracklets
     for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
     {
       //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
@@ -364,7 +434,7 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
     }
 
     // get multiplicity from ESD tracks
-    TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
+    TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD, (fAnalysisMode == AliPWG0Helper::kTPC));
     Int_t nGoodTracks = list->GetEntries();
 
     Printf("Accepted %d tracks", nGoodTracks);
@@ -394,6 +464,78 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
     }
 
     delete list;
+
+    if (eventTriggered && vtxESD)
+    {
+      // collect values for primaries and secondaries
+      for (Int_t iTrack = 0; iTrack < fESD->GetNumberOfTracks(); iTrack++)
+      {
+        AliESDtrack* track = 0;
+
+        if (fAnalysisMode == AliPWG0Helper::kTPC)
+          track = AliESDtrackCuts::GetTPCOnlyTrack(fESD, iTrack);
+        else if (fAnalysisMode == AliPWG0Helper::kTPCITS)
+          track = fESD->GetTrack(iTrack);
+
+        if (!track)
+          continue;
+
+        Int_t label = TMath::Abs(track->GetLabel());
+        if (!stack->Particle(label) || !stack->Particle(label)->GetPDG())
+        {
+          Printf("WARNING: No particle for %d", label);
+          if (stack->Particle(label))
+            stack->Particle(label)->Print();
+          continue;
+        }
+
+        if (stack->Particle(label)->GetPDG()->Charge() == 0)
+          continue;
+
+        if (stack->IsPhysicalPrimary(label))
+        {
+          // primary
+          if (fEsdTrackCutsPrim->AcceptTrack(track)) {
+            if (track->Pt() > 0) {
+              Float_t b0, b1;
+              track->GetImpactParameters(b0, b1);
+              if (fTemp1)
+                ((TH3*) fTemp1)->Fill(b0, b1, 1.0 / track->Pt());
+              //fTemp1->Fill(TMath::Sqrt(b0*b0 + b1*b1), 1.0 / track->Pt());
+            }
+
+            if (AliESDtrackCuts::GetSigmaToVertex(track) > 900)
+            {
+              Printf("Track %d has nsigma of %f. Printing track and vertex...", iTrack, AliESDtrackCuts::GetSigmaToVertex(track));
+              Float_t b[2];
+              Float_t r[3];
+              track->GetImpactParameters(b, r);
+              Printf("Impact parameter %f %f and resolution: %f %f %f", b[0], b[1], r[0], r[1], r[2]);
+              track->Print("");
+              if (vtxESD)
+                vtxESD->Print();
+            }
+          }
+        }
+        else
+        {
+          // secondary
+          if (fEsdTrackCutsSec->AcceptTrack(track))  {
+            if (track->Pt() > 0) {
+              Float_t b0, b1;
+              track->GetImpactParameters(b0, b1);
+              if (fTemp2)
+                ((TH3*) fTemp2)->Fill(b0, b1, 1.0 / track->Pt());
+              //fTemp2->Fill(TMath::Sqrt(b0*b0 + b1*b1), 1.0 / track->Pt());
+            }
+          }
+        }
+
+        // TODO mem leak in the continue statements above
+        if (fAnalysisMode == AliPWG0Helper::kTPC)
+          delete track;
+      }
+    }
   }
   else
     return;
@@ -418,11 +560,18 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
       continue;
 
     if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
+    {
+      //if (TMath::Abs(particle->GetPdgCode()) > 3000 && TMath::Abs(particle->Eta()) < 1.0)
+      //  fPIDParticles->Fill(particle->GetPdgCode());
       continue;
+    }
 
     if (SignOK(particle->GetPDG()) == kFALSE)
       continue;
 
+    if (fPIDParticles && TMath::Abs(particle->Eta()) < 1.0)
+      fPIDParticles->Fill(particle->GetPdgCode());
+
     Float_t eta = particle->Eta();
     Float_t pt = particle->Pt();
 
@@ -481,6 +630,8 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
       continue;
     }
+
+    fPIDTracks->Fill(particle->GetPdgCode());
     
     // find particle that is filled in the correction map
     // this should be the particle that has been reconstructed
@@ -512,6 +663,10 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
 
       processed++;
 
+      // resolutions
+      fEtaResolution->Fill(particle->Eta() - etaArr[i]);
+      fpTResolution->Fill(particle->Pt() - ptArr[i]);
+
       Bool_t firstIsPrim = stack->IsPhysicalPrimary(label);
       // in case of primary the MC values are filled, otherwise (background) the reconstructed values
       if (label == label2 && firstIsPrim)
@@ -540,7 +695,7 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
           fdNdEtaCorrectionProcessType[0]->FillTrackedParticle(vtxMC[2], particle->Eta(), particle->Pt());
 
         // single diffractive
-        if (processType == AliPWG0Helper::kSD )
+        if (processType == AliPWG0Helper::kSD)
           fdNdEtaCorrectionProcessType[1]->FillTrackedParticle(vtxMC[2], particle->Eta(), particle->Pt());
 
         // double diffractive
@@ -637,54 +792,6 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
   delete[] labelArr2;
   delete[] ptArr;
   delete[] deltaPhiArr;
-
-  // fills the fSigmaVertex histogram (systematic study)
-  if (fSigmaVertexTracks)
-  {
-    // save the old value
-    Float_t oldSigmaVertex = fEsdTrackCuts->GetMinNsigmaToVertex();
-
-    // set to maximum
-    fEsdTrackCuts->SetMinNsigmaToVertex(6);
-
-    TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
-    Int_t nGoodTracks = list->GetEntries();
-
-    // loop over esd tracks
-    for (Int_t i=0; i<nGoodTracks; i++)
-    {
-      AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
-      if (!esdTrack)
-      {
-        AliError(Form("ERROR: Could not retrieve track %d.", i));
-        continue;
-      }
-
-      Float_t sigma2Vertex = fEsdTrackCuts->GetSigmaToVertex(esdTrack);
-
-      for (Double_t nSigma = 0.1; nSigma < 6.05; nSigma += 0.1)
-      {
-        if (sigma2Vertex < nSigma)
-        {
-          fSigmaVertexTracks->Fill(nSigma);
-
-          Int_t label = TMath::Abs(esdTrack->GetLabel());
-          TParticle* particle = stack->Particle(label);
-          if (!particle || label >= nPrim)
-            continue;
-
-          if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim))
-            fSigmaVertexPrim->Fill(nSigma);
-        }
-      }
-    }
-
-    delete list;
-    list = 0;
-
-    // set back the old value
-    fEsdTrackCuts->SetMinNsigmaToVertex(oldSigmaVertex);
-  }
 }
 
 void AlidNdEtaCorrectionTask::Terminate(Option_t *)
@@ -714,8 +821,18 @@ void AlidNdEtaCorrectionTask::Terminate(Option_t *)
   fileName.Form("correction_map%s.root", fOption.Data());
   TFile* fout = new TFile(fileName, "RECREATE");
 
+  fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCuts"));
   if (fEsdTrackCuts)
     fEsdTrackCuts->SaveHistograms("esd_track_cuts");
+
+  fEsdTrackCutsPrim = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCutsPrim"));
+  if (fEsdTrackCutsPrim)
+    fEsdTrackCutsPrim->SaveHistograms("esd_track_cuts_primaries");
+
+  fEsdTrackCutsSec = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCutsSec"));
+  if (fEsdTrackCutsSec)
+    fEsdTrackCutsSec->SaveHistograms("esd_track_cuts_secondaries");
+
   fdNdEtaCorrection->SaveHistograms();
   fdNdEtaAnalysisMC->SaveHistograms();
   fdNdEtaAnalysisESD->SaveHistograms();
@@ -727,12 +844,12 @@ void AlidNdEtaCorrectionTask::Terminate(Option_t *)
     if (fdNdEtaCorrectionProcessType[i])
       fdNdEtaCorrectionProcessType[i]->SaveHistograms();
 
-  fSigmaVertexTracks = dynamic_cast<TH1F*> (fOutput->FindObject("fSigmaVertexTracks"));
-  if (fSigmaVertexTracks)
-    fSigmaVertexTracks->Write();
-  fSigmaVertexPrim = dynamic_cast<TH1F*> (fOutput->FindObject("fSigmaVertexPrim"));
-  if (fSigmaVertexPrim)
-    fSigmaVertexPrim->Write();
+  fTemp1 = dynamic_cast<TH1*> (fOutput->FindObject("fTemp1"));
+  if (fTemp1)
+    fTemp1->Write();
+  fTemp2 = dynamic_cast<TH1*> (fOutput->FindObject("fTemp2"));
+  if (fTemp2)
+    fTemp2->Write();
 
   fVertexCorrelation = dynamic_cast<TH2F*> (fOutput->FindObject("fVertexCorrelation"));
   if (fVertexCorrelation)
@@ -740,6 +857,9 @@ void AlidNdEtaCorrectionTask::Terminate(Option_t *)
   fVertexProfile = dynamic_cast<TProfile*> (fOutput->FindObject("fVertexProfile"));
   if (fVertexProfile)
     fVertexProfile->Write();
+  fVertexShift = dynamic_cast<TH1F*> (fOutput->FindObject("fVertexShift"));
+  if (fVertexShift)
+    fVertexShift->Write();
   fVertexShiftNorm = dynamic_cast<TH1F*> (fOutput->FindObject("fVertexShiftNorm"));
   if (fVertexShiftNorm)
     fVertexShiftNorm->Write();
@@ -750,6 +870,12 @@ void AlidNdEtaCorrectionTask::Terminate(Option_t *)
   fEtaProfile = dynamic_cast<TProfile*> (fOutput->FindObject("fEtaProfile"));
   if (fEtaProfile)
     fEtaProfile->Write();
+  fEtaResolution = dynamic_cast<TH1F*> (fOutput->FindObject("fEtaResolution"));
+  if (fEtaResolution)
+    fEtaResolution->Write();
+  fpTResolution = dynamic_cast<TH1F*> (fOutput->FindObject("fpTResolution"));
+  if (fpTResolution)
+    fpTResolution->Write();
 
   fMultAll = dynamic_cast<TH1F*> (fOutput->FindObject("fMultAll"));
   if (fMultAll)
@@ -774,30 +900,35 @@ void AlidNdEtaCorrectionTask::Terminate(Option_t *)
   if (fEventStats)
     fEventStats->Write();
 
-  fout->Write();
-  fout->Close();
+  fPIDParticles = dynamic_cast<TH1F*> (fOutput->FindObject("fPIDParticles"));
+  if (fPIDParticles)
+    fPIDParticles->Write();
+
+  fPIDTracks = dynamic_cast<TH1F*> (fOutput->FindObject("fPIDTracks"));
+  if (fPIDTracks)
+    fPIDTracks->Write();
 
 //fdNdEtaCorrection->DrawHistograms();
+ //fdNdEtaCorrection->DrawHistograms();
 
   Printf("Writting result to %s", fileName.Data());
 
   if (fPIDParticles && fPIDTracks)
   {
-    new TCanvas("pidcanvas", "pidcanvas", 500, 500);
-
-    fPIDParticles->Draw();
-    fPIDTracks->SetLineColor(2);
-    fPIDTracks->Draw("SAME");
-
     TDatabasePDG* pdgDB = new TDatabasePDG;
 
     for (Int_t i=0; i <= fPIDParticles->GetNbinsX()+1; ++i)
       if (fPIDParticles->GetBinContent(i) > 0)
-        printf("PDG = %d (%s): generated: %d, reconstructed: %d, ratio: %f\n", (Int_t) fPIDParticles->GetBinCenter(i), pdgDB->GetParticle((Int_t) fPIDParticles->GetBinCenter(i))->GetName(), (Int_t) fPIDParticles->GetBinContent(i), (Int_t) fPIDTracks->GetBinContent(i), ((fPIDTracks->GetBinContent(i) > 0) ? fPIDParticles->GetBinContent(i) / fPIDTracks->GetBinContent(i) : -1));
+      {
+        TObject* pdgParticle = pdgDB->GetParticle((Int_t) fPIDParticles->GetBinCenter(i));
+        printf("PDG = %d (%s): generated: %d, reconstructed: %d, ratio: %f\n", (Int_t) fPIDParticles->GetBinCenter(i), (pdgParticle) ? pdgParticle->GetName() : "not found", (Int_t) fPIDParticles->GetBinContent(i), (Int_t) fPIDTracks->GetBinContent(i), ((fPIDTracks->GetBinContent(i) > 0) ? fPIDParticles->GetBinContent(i) / fPIDTracks->GetBinContent(i) : -1));
+      }
 
     delete pdgDB;
     pdgDB = 0;
   }
+
+  fout->Write();
+  fout->Close();
 }
 
 Bool_t AlidNdEtaCorrectionTask::SignOK(TParticlePDG* particle)
index 0c3b472..fe6a0e3 100644 (file)
@@ -10,6 +10,7 @@
 class AliESDtrackCuts;
 class dNdEtaAnalysis;
 class AlidNdEtaCorrection;
+class TH1;
 class TH1F;
 class AliESDEvent;
 class TParticlePDG;
@@ -18,7 +19,8 @@ class TProfile;
 
 class AlidNdEtaCorrectionTask : public AliAnalysisTask {
   public:
-    AlidNdEtaCorrectionTask(const char* opt = "");
+    AlidNdEtaCorrectionTask();
+    AlidNdEtaCorrectionTask(const char* opt);
     virtual ~AlidNdEtaCorrectionTask();
 
     virtual void   ConnectInputData(Option_t *);
@@ -55,19 +57,25 @@ class AlidNdEtaCorrectionTask : public AliAnalysisTask {
     // control histograms
     TH1F* fPIDParticles;                         //! pid of primary particles
     TH1F* fPIDTracks;                            //! pid of reconstructed tracks
+
     TH2F* fVertexCorrelation;                    //! ESD z-vtx vs MC z-vtx
     TProfile* fVertexProfile;                    //! Profile of MC z-vtx - ESD z-vtx vs. MC z-vtx
+    TH1F* fVertexShift;                          //! (MC z-vtx - ESD z-vtx) in +- 10 cm
     TH1F* fVertexShiftNorm;                      //! (MC z-vtx - ESD z-vtx) / (sigma_ESD-z-vtx) histogrammed
 
     TH2F* fEtaCorrelation;                       //! ESD eta vs MC eta
     TProfile* fEtaProfile;                       //! Profile of MC eta - ESD eta vs. MC eta
+    TH1F* fEtaResolution;                        //! MC eta - ESD eta in |eta| < 1
+
+    TH1F* fpTResolution;                         //! MC pT - ESD pT in |eta| < 1
+
+    AliESDtrackCuts*  fEsdTrackCutsPrim;         //! control histograms for primaries
+    AliESDtrackCuts*  fEsdTrackCutsSec;          //! control histograms for secondaries
 
     // histograms for systematic studies (must be enabled with option)
 
-    TH1F* fSigmaVertexTracks;                    //! (accepted tracks) vs (n of sigma to vertex cut)
-    TH1F* fSigmaVertexPrim;                      //! (accepted primaries) vs (n of sigma to vertex cut)
-                                                 // enable with option: sigma-vertex
+    TH1* fTemp1;                                 //! temp histogram for quick study of variables
+    TH1* fTemp2;                                 //! temp histogram for quick study of variables
 
     TH1F* fMultAll; //! primary particles  in |eta| < 1 and pT > 0.2 in all events
     TH1F* fMultTr; //! primary particles  in |eta| < 1 and pT > 0.2 in triggered events
index 46e1a02..2b8c665 100644 (file)
@@ -173,11 +173,11 @@ void AlidNdEtaTask::CreateOutputObjects()
   fPhi = new TH1F("fPhi", "fPhi;#phi in rad.;count", 720, 0, 2 * TMath::Pi());
   fOutput->Add(fPhi);
 
-  fEtaPhi = new TH2F("fEtaPhi", "fEtaPhi;#eta;#phi in rad.;count", 80, -4, 4, 80, 0, 2 * TMath::Pi());
+  fEtaPhi = new TH2F("fEtaPhi", "fEtaPhi;#eta;#phi in rad.;count", 80, -4, 4, 18*5, 0, 2 * TMath::Pi());
   fOutput->Add(fEtaPhi);
 
   if (fAnalysisMode == AliPWG0Helper::kSPD)
-    fDeltaPhi = new TH1F("fDeltaPhi", "fDeltaPhi;#Delta #phi;Entries", 1000, -3.14, 3.14);
+    fDeltaPhi = new TH1F("fDeltaPhi", "fDeltaPhi;#Delta #phi;Entries", 18*50, -3.14, 3.14);
 
   if (fReadMC)
   {
@@ -203,227 +203,269 @@ void AlidNdEtaTask::CreateOutputObjects()
     fPartPt->Sumw2();
     fOutput->Add(fPartPt);
   }
+
+  if (fEsdTrackCuts)
+  {
+    fEsdTrackCuts->SetName("fEsdTrackCuts");
+    fOutput->Add(fEsdTrackCuts);
+  }
 }
 
 void AlidNdEtaTask::Exec(Option_t*)
 {
   // process the event
 
-  // Check prerequisites
-  if (!fESD)
+  // these variables are also used in the MC section below; however, if ESD is off they stay with their default values
+  Bool_t eventTriggered = kFALSE;
+  const AliESDVertex* vtxESD = 0;
+
+  // ESD analysis
+  if (fESD)
   {
-    AliDebug(AliLog::kError, "ESD branch not available");
-    return;
-  }
+    // trigger definition
+    eventTriggered = AliPWG0Helper::IsEventTriggered(fESD->GetTriggerMask(), fTrigger);
 
-  // trigger definition
-  Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD->GetTriggerMask(), fTrigger);
+    // get the ESD vertex
+    vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
 
-  // get the ESD vertex
-  const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
-  
-  Double_t vtx[3];
-  if (vtxESD) {
-    vtxESD->GetXYZ(vtx);
-  }
-  else
-    Printf("No vertex");
-  
-  // fill z vertex resolution
-  if (vtxESD)
-    fVertexResolution->Fill(vtxESD->GetZRes());
-
-  // post the data already here
-  PostData(0, fOutput);
-
-  // needed for syst. studies
-  AliStack* stack = 0;
-  
-  if (fUseMCVertex || fUseMCKine) {
-    if (!fReadMC) {
-      Printf("ERROR: fUseMCVertex or fUseMCKine set without fReadMC set!");
-      return;
+    Double_t vtx[3];
+    if (vtxESD) {
+      vtxESD->GetXYZ(vtx);
     }
+    else
+      Printf("No vertex");
 
-    AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
-    if (!eventHandler) {
-      Printf("ERROR: Could not retrieve MC event handler");
-      return;
-    }
+    // fill z vertex resolution
+    if (vtxESD)
+      fVertexResolution->Fill(vtxESD->GetZRes());
 
-    AliMCEvent* mcEvent = eventHandler->MCEvent();
-    if (!mcEvent) {
-      Printf("ERROR: Could not retrieve MC event");
-      return;
-    }
+    // post the data already here
+    PostData(0, fOutput);
 
-    AliHeader* header = mcEvent->Header();
-    if (!header)
-    {
-      AliDebug(AliLog::kError, "Header not available");
-      return;
-    }
+    // needed for syst. studies
+    AliStack* stack = 0;
+    TArrayF vtxMC(3);
+
+    if (fUseMCVertex || fUseMCKine || fReadMC) {
+      if (!fReadMC) {
+        Printf("ERROR: fUseMCVertex or fUseMCKine set without fReadMC set!");
+        return;
+      }
+
+      AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
+      if (!eventHandler) {
+        Printf("ERROR: Could not retrieve MC event handler");
+        return;
+      }
+
+      AliMCEvent* mcEvent = eventHandler->MCEvent();
+      if (!mcEvent) {
+        Printf("ERROR: Could not retrieve MC event");
+        return;
+      }
+
+      AliHeader* header = mcEvent->Header();
+      if (!header)
+      {
+        AliDebug(AliLog::kError, "Header not available");
+        return;
+      }
 
-    if (fUseMCVertex)
-    {
-      Printf("WARNING: Replacing vertex by MC vertex. This is for systematical checks only.");
       // get the MC vertex
       AliGenEventHeader* genHeader = header->GenEventHeader();
-      TArrayF vtxMC(3);
+      if (!genHeader)
+      {
+        AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
+        return;
+      }
       genHeader->PrimaryVertex(vtxMC);
 
-      vtx[2] = vtxMC[2];
-    }
+      if (fUseMCVertex)
+      {
+        Printf("WARNING: Replacing vertex by MC vertex. This is for systematical checks only.");
+        vtx[2] = vtxMC[2];
+      }
 
-    if (fUseMCKine)
-    {
-      stack = mcEvent->Stack();
-      if (!stack)
+      if (fUseMCKine)
       {
-        AliDebug(AliLog::kError, "Stack not available");
-        return;
+        stack = mcEvent->Stack();
+        if (!stack)
+        {
+          AliDebug(AliLog::kError, "Stack not available");
+          return;
+        }
       }
     }
-  }
 
-  // create list of (label, eta, pt) tuples
-  Int_t inputCount = 0;
-  Int_t* labelArr = 0;
-  Float_t* etaArr = 0;
-  Float_t* ptArr = 0;
-  if (fAnalysisMode == AliPWG0Helper::kSPD)
-  {
-    // get tracklets
-    const AliMultiplicity* mult = fESD->GetMultiplicity();
-    if (!mult)
+    // create list of (label, eta, pt) tuples
+    Int_t inputCount = 0;
+    Int_t* labelArr = 0;
+    Float_t* etaArr = 0;
+    Float_t* ptArr = 0;
+    if (fAnalysisMode == AliPWG0Helper::kSPD)
     {
-      AliDebug(AliLog::kError, "AliMultiplicity not available");
-      return;
-    }
+      // get tracklets
+      const AliMultiplicity* mult = fESD->GetMultiplicity();
+      if (!mult)
+      {
+        AliDebug(AliLog::kError, "AliMultiplicity not available");
+        return;
+      }
 
-    labelArr = new Int_t[mult->GetNumberOfTracklets()];
-    etaArr = new Float_t[mult->GetNumberOfTracklets()];
-    ptArr = new Float_t[mult->GetNumberOfTracklets()];
+      labelArr = new Int_t[mult->GetNumberOfTracklets()];
+      etaArr = new Float_t[mult->GetNumberOfTracklets()];
+      ptArr = new Float_t[mult->GetNumberOfTracklets()];
 
-    if (fUseMCKine && stack)
-      Printf("Processing only primaries (MC information used). This is for systematical checks only.");
+      if (fUseMCKine && stack)
+        Printf("Processing only primaries (MC information used). This is for systematical checks only.");
 
-    // get multiplicity from ITS tracklets
-    for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
-    {
-      //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
+      // get multiplicity from ITS tracklets
+      for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
+      {
+        //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
+
+        if (fUseMCKine && stack)
+          if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
+            continue;
+
+        Float_t deltaPhi = mult->GetDeltaPhi(i);
+        // prevent values to be shifted by 2 Pi()
+        if (deltaPhi < -TMath::Pi())
+          deltaPhi += TMath::Pi() * 2;
+        if (deltaPhi > TMath::Pi())
+          deltaPhi -= TMath::Pi() * 2;
+
+        if (TMath::Abs(deltaPhi) > 1)
+          printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
+
+        Float_t phi = mult->GetPhi(i);
+        if (phi < 0)
+          phi += TMath::Pi() * 2;
+        fPhi->Fill(phi);
+        fEtaPhi->Fill(mult->GetEta(i), phi);
+
+        fDeltaPhi->Fill(deltaPhi);
+
+        etaArr[inputCount] = mult->GetEta(i);
+        labelArr[inputCount] = mult->GetLabel(i, 0);
+        ptArr[inputCount] = 0; // no pt for tracklets
+        ++inputCount;
+      }
 
-      if (fUseMCKine && stack)
-        if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
-          continue;
-      
-      Float_t deltaPhi = mult->GetDeltaPhi(i);
-      // prevent values to be shifted by 2 Pi()
-      if (deltaPhi < -TMath::Pi())
-        deltaPhi += TMath::Pi() * 2;
-      if (deltaPhi > TMath::Pi())
-        deltaPhi -= TMath::Pi() * 2;
-
-      if (TMath::Abs(deltaPhi) > 1)
-        printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
-
-      Float_t phi = mult->GetPhi(i);
-      if (phi < 0)
-        phi += TMath::Pi() * 2;
-      fPhi->Fill(phi);
-      fEtaPhi->Fill(mult->GetEta(i), phi);
-
-      fDeltaPhi->Fill(deltaPhi);
-
-      etaArr[inputCount] = mult->GetEta(i);
-      labelArr[inputCount] = mult->GetLabel(i, 0);
-      ptArr[inputCount] = 0; // no pt for tracklets
-      ++inputCount;
+      //Printf("Accepted %d tracklets", inputCount);
     }
-
-    //Printf("Accepted %d tracklets", inputCount);
-  }
-  else if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS)
-  {
-    if (!fEsdTrackCuts)
+    else if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS)
     {
-      AliDebug(AliLog::kError, "fESDTrackCuts not available");
-      return;
-    }
+      if (!fEsdTrackCuts)
+      {
+        AliDebug(AliLog::kError, "fESDTrackCuts not available");
+        return;
+      }
 
-    // get multiplicity from ESD tracks
-    TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
-    Int_t nGoodTracks = list->GetEntries();
+      // get multiplicity from ESD tracks
+      TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
+      Int_t nGoodTracks = list->GetEntries();
 
-    labelArr = new Int_t[nGoodTracks];
-    etaArr = new Float_t[nGoodTracks];
-    ptArr = new Float_t[nGoodTracks];
+      labelArr = new Int_t[nGoodTracks];
+      etaArr = new Float_t[nGoodTracks];
+      ptArr = new Float_t[nGoodTracks];
 
-    // loop over esd tracks
-    for (Int_t i=0; i<nGoodTracks; i++)
-    {
-      AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
-      if (!esdTrack)
+      // loop over esd tracks
+      for (Int_t i=0; i<nGoodTracks; i++)
       {
-        AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
-        continue;
-      }
+        AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
+        if (!esdTrack)
+        {
+          AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
+          continue;
+        }
 
-      Float_t phi = esdTrack->Phi();
-      if (phi < 0)
-        phi += TMath::Pi() * 2;
-      fPhi->Fill(phi);
-      fEtaPhi->Fill(esdTrack->Eta(), phi);
+        Float_t phi = esdTrack->Phi();
+        if (phi < 0)
+          phi += TMath::Pi() * 2;
+        fPhi->Fill(phi);
+        fEtaPhi->Fill(esdTrack->Eta(), phi);
 
-      etaArr[inputCount] = esdTrack->Eta();
-      labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
-      ptArr[inputCount] = esdTrack->Pt();
-      ++inputCount;
-    }
+        etaArr[inputCount] = esdTrack->Eta();
+        labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
+        ptArr[inputCount] = esdTrack->Pt();
+        ++inputCount;
+      }
 
-    Printf("Accepted %d tracks", nGoodTracks);
+      Printf("Accepted %d tracks", nGoodTracks);
 
-    delete list;
-  }
-  else
-    return;
-
-  // Processing of ESD information (always)
-  if (eventTriggered)
-  {
-    // control hist
-    fMult->Fill(inputCount);
-    fdNdEtaAnalysisESD->FillTriggeredEvent(inputCount);
+      delete list;
+    }
+    else
+      return;
 
-    if (vtxESD)
+    // Processing of ESD information (always)
+    if (eventTriggered)
     {
       // control hist
-      fMultVtx->Fill(inputCount);
+      fMult->Fill(inputCount);
+      fdNdEtaAnalysisESD->FillTriggeredEvent(inputCount);
 
-      for (Int_t i=0; i<inputCount; ++i)
+      if (vtxESD)
       {
-        Float_t eta = etaArr[i];
-        Float_t pt  = ptArr[i];
-
-        fdNdEtaAnalysisESD->FillTrack(vtx[2], eta, pt);
+        // control hist
+        fMultVtx->Fill(inputCount);
 
-        if (TMath::Abs(vtx[2]) < 20)
+        for (Int_t i=0; i<inputCount; ++i)
         {
-          fPartEta[0]->Fill(eta);
+          Float_t eta = etaArr[i];
+          Float_t pt  = ptArr[i];
 
-          if (vtx[2] < 0)
-            fPartEta[1]->Fill(eta);
-          else
-            fPartEta[2]->Fill(eta);
+          fdNdEtaAnalysisESD->FillTrack(vtx[2], eta, pt);
+
+          if (TMath::Abs(vtx[2]) < 20)
+          {
+            fPartEta[0]->Fill(eta);
+
+            if (vtx[2] < 0)
+              fPartEta[1]->Fill(eta);
+            else
+              fPartEta[2]->Fill(eta);
+          }
         }
-      }
 
-      // for event count per vertex
-      fdNdEtaAnalysisESD->FillEvent(vtx[2], inputCount);
+        // for event count per vertex
+        fdNdEtaAnalysisESD->FillEvent(vtx[2], inputCount);
 
-      // control hist
-      fEvents->Fill(vtx[2]);
+        // control hist
+        fEvents->Fill(vtx[2]);
+
+        if (fReadMC)
+        {
+          // from tracks is only done for triggered and vertex reconstructed events
+          for (Int_t i=0; i<inputCount; ++i)
+          {
+            Int_t label = labelArr[i];
+
+            if (label < 0)
+              continue;
+
+            //Printf("Getting particle of track %d", label);
+            TParticle* particle = stack->Particle(label);
+
+            if (!particle)
+            {
+              AliDebug(AliLog::kError, Form("ERROR: Could not retrieve particle %d.", label));
+              continue;
+            }
+
+            fdNdEtaAnalysisTracks->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
+          } // end of track loop
+
+          // for event count per vertex
+          fdNdEtaAnalysisTracks->FillEvent(vtxMC[2], inputCount);
+        }
+      }
     }
+
+    delete[] etaArr;
+    delete[] labelArr;
+    delete[] ptArr;
   }
 
   if (fReadMC)   // Processing of MC information (optional)
@@ -440,7 +482,7 @@ void AlidNdEtaTask::Exec(Option_t*)
       return;
     }
 
-    stack = mcEvent->Stack();
+    AliStack* stack = mcEvent->Stack();
     if (!stack)
     {
       AliDebug(AliLog::kError, "Stack not available");
@@ -479,16 +521,19 @@ void AlidNdEtaTask::Exec(Option_t*)
 
     for (Int_t iMc = 0; iMc < nPrim; ++iMc)
     {
+      if (!stack->IsPhysicalPrimary(iMc))
+        continue;
+
       //Printf("Getting particle %d", iMc);
       TParticle* particle = stack->Particle(iMc);
 
       if (!particle)
         continue;
 
-      if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
+      if (particle->GetPDG()->Charge() == 0)
         continue;
 
-      AliDebug(AliLog::kDebug+1, Form("Accepted primary %d, unique ID: %d", iMc, particle->GetUniqueID()));
+      //AliDebug(AliLog::kDebug+1, Form("Accepted primary %d, unique ID: %d", iMc, particle->GetUniqueID()));
       Float_t eta = particle->Eta();
       Float_t pt = particle->Pt();
 
@@ -523,38 +568,7 @@ void AlidNdEtaTask::Exec(Option_t*)
       if (vtxESD)
         fdNdEtaAnalysisTrVtx->FillEvent(vtxMC[2], nAcceptedParticles);
     }
-
-    if (eventTriggered && vtxESD)
-    {
-      // from tracks is only done for triggered and vertex reconstructed events
-
-      for (Int_t i=0; i<inputCount; ++i)
-      {
-        Int_t label = labelArr[i];
-
-        if (label < 0)
-          continue;
-
-        //Printf("Getting particle of track %d", label);
-        TParticle* particle = stack->Particle(label);
-
-        if (!particle)
-        {
-          AliDebug(AliLog::kError, Form("ERROR: Could not retrieve particle %d.", label));
-          continue;
-        }
-
-        fdNdEtaAnalysisTracks->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
-      } // end of track loop
-
-      // for event count per vertex
-      fdNdEtaAnalysisTracks->FillEvent(vtxMC[2], inputCount);
-    }
   }
-
-  delete[] etaArr;
-  delete[] labelArr;
-  delete[] ptArr;
 }
 
 void AlidNdEtaTask::Terminate(Option_t *)
@@ -564,22 +578,25 @@ void AlidNdEtaTask::Terminate(Option_t *)
   // the results graphically or save the results to file.
 
   fOutput = dynamic_cast<TList*> (GetOutputData(0));
-  if (!fOutput) {
+  if (!fOutput)
     Printf("ERROR: fOutput not available");
-    return;
-  }
 
-  fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("fdNdEtaAnalysisESD"));
-  fMult = dynamic_cast<TH1F*> (fOutput->FindObject("fMult"));
-  fMultVtx = dynamic_cast<TH1F*> (fOutput->FindObject("fMultVtx"));
-  for (Int_t i=0; i<3; ++i)
-    fPartEta[i] = dynamic_cast<TH1F*> (fOutput->FindObject(Form("dndeta_check_%d", i)));
-  fEvents = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_vertex"));
-  fVertexResolution = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_vertex_resolution_z"));
+  if (fOutput)
+  {
+    fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("fdNdEtaAnalysisESD"));
+    fMult = dynamic_cast<TH1F*> (fOutput->FindObject("fMult"));
+    fMultVtx = dynamic_cast<TH1F*> (fOutput->FindObject("fMultVtx"));
+    for (Int_t i=0; i<3; ++i)
+      fPartEta[i] = dynamic_cast<TH1F*> (fOutput->FindObject(Form("dndeta_check_%d", i)));
+    fEvents = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_vertex"));
+    fVertexResolution = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_vertex_resolution_z"));
 
-  fPhi = dynamic_cast<TH1F*> (fOutput->FindObject("fPhi"));
-  fEtaPhi = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaPhi"));
-  fDeltaPhi = dynamic_cast<TH1F*> (fOutput->FindObject("fDeltaPhi"));
+    fPhi = dynamic_cast<TH1F*> (fOutput->FindObject("fPhi"));
+    fEtaPhi = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaPhi"));
+    fDeltaPhi = dynamic_cast<TH1F*> (fOutput->FindObject("fDeltaPhi"));
+
+    fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCuts"));
+  }
 
   if (!fdNdEtaAnalysisESD)
   {
@@ -666,13 +683,16 @@ void AlidNdEtaTask::Terminate(Option_t *)
 
   if (fReadMC)
   {
-    fdNdEtaAnalysis = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta"));
-    fdNdEtaAnalysisNSD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaNSD"));
-    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"));
-    fVertex = dynamic_cast<TH3F*> (fOutput->FindObject("vertex_check"));
+    if (fOutput)
+    {
+      fdNdEtaAnalysis = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta"));
+      fdNdEtaAnalysisNSD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaNSD"));
+      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"));
+      fVertex = dynamic_cast<TH3F*> (fOutput->FindObject("vertex_check"));
+    }
 
     if (!fdNdEtaAnalysis || !fdNdEtaAnalysisTr || !fdNdEtaAnalysisTrVtx || !fPartPt)
     {
index e0956e3..5ad6e0c 100644 (file)
@@ -391,14 +391,18 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
       Int_t vertexBinEnd   = vertexHist->GetXaxis()->FindBin(vertexRangeEnd[vertexPos]);
 
       const Int_t *binBegin = 0;
-      // adjust acceptance range for SPD
+
+      // adjust acceptance range
       if (fAnalysisMode == AliPWG0Helper::kSPD)
       {
-        //const Int_t binBegin[30] = { 18, 16, 15, 14, 13, 13, 12, 11, 9, 7, 6, 5, 5, 5, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 1 };
+        //const Int_t binBeginSPD[30] = { 18, 16, 15, 14, 13, 13, 12, 11, 9, 7, 6, 5, 5, 5, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 1 }; // by eye
         const Int_t binBeginSPD[30] = { -1, -1, -1, -1, 16, 14, 12, 10, 9, 7, 6, 5, 5, 4, 4, 4, 4, 3, 3, 3, 3, 2, 2, 2, 2, 2, -1, -1, -1, -1}; // limit in correction map is 5
+
         //const Int_t binBegin[30] = { -1, -1, -1, 17, 15, 14, 12, 10, 8, 7, 6, 5, 4, 4, 4, 4, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, -1, -1, -1};  // limit in correction map is 10
-        binBegin = binBeginSPD;
 
+        //const Int_t binBeginSPD[30] = { -1, -1, -1, -1, 16, 15, 13, 11, 9, 8, 7, 6, 5, 5, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, -1, -1, -1, -1}; // limit 2
+
+        binBegin = binBeginSPD;
       }
       else if (fAnalysisMode == AliPWG0Helper::kTPC)
       {
@@ -418,14 +422,15 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
       if (vtxBegin == -1)
         continue;
       
-      //Float_t eta = vtxVsEta->GetYaxis()->GetBinCenter(iEta);
       //Printf("%d %d | %d %d", vtxBegin, vertexHist->GetXaxis()->FindBin(GetVtxMin(eta)), vtxEnd, vertexHist->GetXaxis()->FindBin(-GetVtxMin(-eta)));
       //vtxBegin = vertexHist->GetXaxis()->FindBin(GetVtxMin(eta));
       //vtxEnd = vertexHist->GetXaxis()->FindBin(-GetVtxMin(-eta));
-      //printf("Eta bin: %d (%f) Vertex range: %d; before: %d %d ", iEta, eta, vertexPos, vertexBinBegin, vertexBinEnd);
+
+      Float_t eta = vtxVsEta->GetYaxis()->GetBinCenter(iEta);
+      printf("Eta bin: %d (%f) Vertex range: %d; before: %d %d (range) %d %d (acceptance)", iEta, eta, vertexPos, vertexBinBegin, vertexBinEnd, vtxBegin, vtxEnd);
       vertexBinBegin = TMath::Max(vertexBinBegin, vtxBegin);
       vertexBinEnd =   TMath::Min(vertexBinEnd, vtxEnd);
-      //Printf("after:  %d %d", vertexBinBegin, vertexBinEnd);
+      Printf(" after:  %d %d", vertexBinBegin, vertexBinEnd);
 
       // no data for this bin
       if (vertexBinBegin > vertexBinEnd)
@@ -434,26 +439,37 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
         continue;
       }
 
-      Float_t totalEvents = vertexHist->Integral(vertexBinBegin, vertexBinEnd);
-      if (totalEvents == 0)
-      {
-        printf("WARNING: No events for hist %d %d %d\n", vertexPos, vertexBinBegin, vertexBinEnd);
-        continue;
-      }
-
+      Float_t totalEvents = 0;
       Float_t sum = 0;
       Float_t sumError2 = 0;
-      for (Int_t iVtx = vertexBinBegin; iVtx <= vertexBinEnd; iVtx++)
+      Float_t unusedTracks = 0;
+      Float_t unusedEvents = 0;
+      for (Int_t iVtx = 1; iVtx <= vtxVsEta->GetNbinsX(); iVtx++)
       {
-        if (vtxVsEta->GetBinContent(iVtx, iEta) != 0)
+        if (iVtx >= vertexBinBegin && iVtx <= vertexBinEnd)
         {
-          sum = sum + vtxVsEta->GetBinContent(iVtx, iEta);
+          if (vtxVsEta->GetBinContent(iVtx, iEta) != 0)
+          {
+            sum += vtxVsEta->GetBinContent(iVtx, iEta);
 
-          if (sumError2 > 10e30)
-            Printf("WARNING: sum of error2 is dangerously large - be prepared for crash... ");
+            if (sumError2 > 10e30)
+              Printf("WARNING: sum of error2 is dangerously large - be prepared for crash... ");
 
-          sumError2 = sumError2 + TMath::Power(vtxVsEta->GetBinError(iVtx, iEta),2);
+            sumError2 = sumError2 + TMath::Power(vtxVsEta->GetBinError(iVtx, iEta),2);
+          }
+          totalEvents += vertexHist->GetBinContent(iVtx);
         }
+        else
+        {
+          unusedTracks += vtxVsEta->GetBinContent(iVtx, iEta);
+          unusedEvents +=vertexHist->GetBinContent(iVtx);
+        }
+      }
+
+      if (totalEvents == 0)
+      {
+        printf("WARNING: No events for hist %d %d %d\n", vertexPos, vertexBinBegin, vertexBinEnd);
+        continue;
       }
 
       Float_t ptCutOffCorrection = 1;
@@ -486,7 +502,7 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
         fdNdEtaPtCutOffCorrected[vertexPos]->SetBinContent(bin, dndeta);
         fdNdEtaPtCutOffCorrected[vertexPos]->SetBinError(bin, error);
 
-        //Printf("Bin %d has dN/deta = %f +- %f; %f +- %f", bin, dndeta, error, fdNdEta[vertexPos]->GetBinContent(bin), fdNdEta[vertexPos]->GetBinError(bin));
+        Printf("Bin %d has dN/deta = %f +- %f; %.2f tracks %.2f events (outside acceptance: %.2f tracks, %.2f events)", bin, dndeta, error, sum, totalEvents, unusedTracks, unusedEvents);
       }
     }
   }
@@ -592,11 +608,11 @@ void dNdEtaAnalysis::LoadHistograms(const Char_t* dir)
 
   fPtDist = dynamic_cast<TH1F*> (gDirectory->Get(fPtDist->GetName()));
 
-  if (dynamic_cast<TNamed*> (gFile->Get("fTag")))
-    fTag = (dynamic_cast<TNamed*> (gFile->Get("fTag")))->GetTitle();
+  if (dynamic_cast<TNamed*> (gDirectory->Get("fTag")))
+    fTag = (dynamic_cast<TNamed*> (gDirectory->Get("fTag")))->GetTitle();
 
-  if (dynamic_cast<TParameter<Int_t>*> (gFile->Get("fAnalysisMode")))
-    fAnalysisMode = (AliPWG0Helper::AnalysisMode) (dynamic_cast<TParameter<Int_t>*> (gFile->Get("fAnalysisMode")))->GetVal();
+  if (dynamic_cast<TParameter<Int_t>*> (gDirectory->Get("fAnalysisMode")))
+    fAnalysisMode = (AliPWG0Helper::AnalysisMode) (dynamic_cast<TParameter<Int_t>*> (gDirectory->Get("fAnalysisMode")))->GetVal();
 
   gDirectory->cd("../");
 }
index 73a28c6..d3e0354 100644 (file)
@@ -1289,15 +1289,15 @@ void CompareTrack2Particle2D()
 
   Track2Particle2DCreatePlots("correction_maponly-positive.root");
 
-  TH2* posYX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_yx_div_meas_nTrackToNPart_yx")->Clone("pos_yx"));
-  TH2* posZX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zx_div_meas_nTrackToNPart_zx")->Clone("pos_zx"));
-  TH2* posZY = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zy_div_meas_nTrackToNPart_zy")->Clone("pos_zy"));
+  TH2* posYX = dynamic_cast<TH2*> (gROOT->FindObject("generated_yx_div_measured_yx")->Clone("pos_yx"));
+  TH2* posZX = dynamic_cast<TH2*> (gROOT->FindObject("generated_zx_div_measured_zx")->Clone("pos_zx"));
+  TH2* posZY = dynamic_cast<TH2*> (gROOT->FindObject("generated_zy_div_measured_zy")->Clone("pos_zy"));
 
   Track2Particle2DCreatePlots("correction_maponly-negative.root");
 
-  TH2* negYX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_yx_div_meas_nTrackToNPart_yx")->Clone("neg_yx"));
-  TH2* negZX = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zx_div_meas_nTrackToNPart_zx")->Clone("neg_zx"));
-  TH2* negZY = dynamic_cast<TH2*> (gROOT->FindObject("gene_nTrackToNPart_zy_div_meas_nTrackToNPart_zy")->Clone("neg_zy"));
+  TH2* negYX = dynamic_cast<TH2*> (gROOT->FindObject("generated_yx_div_measured_yx")->Clone("neg_yx"));
+  TH2* negZX = dynamic_cast<TH2*> (gROOT->FindObject("generated_zx_div_measured_zx")->Clone("neg_zx"));
+  TH2* negZY = dynamic_cast<TH2*> (gROOT->FindObject("generated_zy_div_measured_zy")->Clone("neg_zy"));
 
   posYX->Divide(negYX);
   posZX->Divide(negZX);
@@ -1729,8 +1729,10 @@ void DrawTrackletOrigin()
   }
 }
 
-void DetermineAcceptance(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction", Double_t ptmin=0.2)
+TH2* GetCorrection(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction", Double_t ptmin=0.2)
 {
+  // returns the correction factor with pt integrated out
+
   loadlibs();
 
   TFile::Open(fileName);
@@ -1768,6 +1770,13 @@ void DetermineAcceptance(const char* fileName = "correction_map.root", const cha
 //   proj = hist->Project3D("yx");
   proj->Draw("COLZ");
 
+  return proj;
+}
+
+void DetermineAcceptance(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction", Double_t ptmin=0.2)
+{
+  TH2* proj = GetCorrection(fileName, dirName, ptmin);
+
   const Float_t limit = 5;
 
   TString array = "{";
@@ -1821,3 +1830,171 @@ void AverageMultiplicity(const char* fileName = "correction_map.root", const cha
   Printf("%f %f --> %f", nEvents, nTracks, nTracks / nEvents);
 }
 
+void GetAverageCorrectionFactor(Float_t etaRange = 1.5, Float_t vertexRange = 9.9, const char* rawFile = "analysis_esd_raw.root", const char* mcFile = "analysis_mc.root")
+{
+  loadlibs();
+
+  TFile::Open(rawFile);
+  dNdEtaAnalysis* raw = new dNdEtaAnalysis("dndeta", "dndeta");
+  raw->LoadHistograms("fdNdEtaAnalysisESD");
+  raw->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->GetXaxis()->SetRangeUser(-vertexRange, vertexRange);
+  tracks = raw->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Project3D("y");
+  events = raw->GetData()->GetEventCorrection()->GetMeasuredHistogram()->ProjectionX("events", 0, raw->GetData()->GetEventCorrection()->GetMeasuredHistogram()->GetNbinsY() + 1);
+  Float_t nEvents = events->Integral(events->FindBin(-vertexRange), events->FindBin(vertexRange));
+  tracks->Scale(1.0 / nEvents / tracks->GetBinWidth(1));
+
+  TFile::Open(mcFile);
+  dNdEtaAnalysis* mc = new dNdEtaAnalysis("dndeta", "dndeta");
+  mc->LoadHistograms("dndetaTrVtx");
+  mcH = mc->GetdNdEtaPtCutOffCorrectedHistogram(0);
+
+  new TCanvas;
+  mcH->SetLineColor(2);
+  mcH->DrawCopy();
+  tracks->DrawCopy("SAME");
+
+  new TCanvas;
+  mcH->GetYaxis()->SetRangeUser(0, 5);
+  mcH->Divide(tracks);
+  mcH->DrawCopy();
+  mcH->Fit("pol0", "", "", -etaRange, etaRange);
+}
+
+void TrackCuts_Comparison(char* histName, Bool_t after = kFALSE, const char* fileName = "correction_map.root")
+{
+  // for the nsigmaplot it is needed to run with all cuts except the nsigmatovertex
+  //    --> manually disable it in the run.C
+
+  file = TFile::Open(fileName);
+
+  Int_t count = 0;
+  Int_t colors[] = { 1, 2, 3, 4, 5, 6 };
+
+  TLegend* legend = new TLegend(0.4, 0.6, 1, 1);
+  TLegend* legend2 = new TLegend(0.4, 0.6, 1, 1);
+  TLegend* legend3 = new TLegend(0.7, 0.5, 1, 0.7);
+
+  TCanvas* c1 = new TCanvas("c1", "c1", 800, 600);
+  //TCanvas* c2 = new TCanvas("c2", "c2", 800, 600);
+  TCanvas* c3 = new TCanvas("c3", "c3", 800, 600);
+
+  const char* folders2[] = { "before_cuts", "after_cuts" };
+  for (Int_t j = 0; j < ((after) ? 2 : 1); j++)
+  {
+    const char* folders1[] = { "esd_track_cuts", "esd_track_cuts_primaries", "esd_track_cuts_secondaries" };
+    TH1* base = 0;
+    TH1* prim = 0;
+    TH1* sec = 0;
+    for (Int_t i = 0; i < 3; i++)
+    {
+      TString folder;
+      folder.Form("%s/%s/%s", folders1[i], folders2[j], histName);
+      TH1* hist = (TH1*) file->Get(folder);
+      legend->AddEntry(hist, folder);
+
+      c1->cd();
+      hist->SetLineColor(colors[count]);
+      hist->DrawCopy((count == 0) ? "" : "SAME");
+
+      switch (i)
+      {
+        case 0: base = hist; break;
+        case 1: prim = hist; break;
+        case 2: sec = hist; break;
+      }
+
+      count++;
+    }
+
+    TH1* eff    = (TH1*) prim->Clone("eff"); eff->Reset();
+    TH1* purity = (TH1*) prim->Clone("purity"); purity->Reset();
+
+    for (Int_t bin = 1; bin <= prim->GetNbinsX(); bin++)
+    {
+      eff->SetBinContent(bin, prim->Integral(1, bin) / prim->Integral(1, prim->GetNbinsX() + 1));
+      purity->SetBinContent(bin, sec->Integral(1, bin) / (prim->Integral(1, bin) + sec->Integral(1, bin)));
+    }
+
+    eff->GetYaxis()->SetRangeUser(0, 1);
+    eff->SetLineColor(colors[0+j*2]);
+    eff->SetStats(kFALSE);
+    purity->SetLineColor(colors[1+j*2]);
+
+    legend3->AddEntry(eff, Form("%s: efficiency", folders2[j]));
+    legend3->AddEntry(purity, Form("%s: contamination", folders2[j]));
+
+    c3->cd();
+    eff->DrawCopy((j == 0) ? "" : "SAME");
+    purity->DrawCopy("SAME");
+  }
+
+  c1->cd();
+  c1->SetLogy();
+  c1->SetGridx();
+  c1->SetGridy();
+  legend->Draw();
+
+  //c2->cd();
+ // c2->SetGridx();
+ // c2->SetGridy();
+  //legend2->Draw();
+
+  c3->cd();
+  c3->SetGridx();
+  c3->SetGridy();
+  legend3->Draw();
+}
+
+void TrackCuts_DCA()
+{
+  file = TFile::Open("correction_map.root");
+  hist = (TH2*) file->Get("esd_track_cuts/before_cuts/dXYvsDZ");
+
+  TCanvas* c1 = new TCanvas("c1", "c1", 600, 600);
+  c1->SetLogz();
+  c1->SetRightMargin(0.12);
+  c1->SetBottomMargin(0.12);
+
+  hist->SetStats(kFALSE);
+  hist->Draw("COLZ");
+
+  ellipse = new TEllipse(0, 0, 4);
+  ellipse->SetLineWidth(2);
+  ellipse->SetLineStyle(2);
+  ellipse->SetFillStyle(0);
+  ellipse->Draw();
+
+  c1->SaveAs("trackcuts_dca_2d.eps");
+}
+
+void FindNSigma(TH2* hist, Int_t nSigma = 1)
+{
+  TH1* proj = hist->ProjectionY();
+  proj->Reset();
+
+  for (Int_t bin=1; bin<=proj->GetNbinsX(); bin++)
+  {
+    if (hist->Integral(1, hist->GetNbinsX(), bin, bin) == 0)
+      continue;
+
+    Int_t limit = -1;
+    for (limit = 1; limit<=hist->GetNbinsX(); limit++)
+  }
+}
+
+void ShowOnlyAccepted(TH2* input, const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction", Double_t ptmin=0.2)
+{
+  TH2* proj = GetCorrection(fileName, dirName, ptmin);
+
+  for (Int_t y=1; y<=proj->GetNbinsY(); ++y)
+    for (Int_t x=1; x<=proj->GetNbinsX(); ++x)
+      if (proj->GetBinContent(x, y) > 5 || proj->GetBinContent(x, y) == 0)
+      {
+        proj->SetBinContent(x, y, 0);
+      }
+      else
+        proj->SetBinContent(x, y, 1);
+
+
+  input->Multiply(proj);
+}
diff --git a/PWG0/dNdEta/makeCorrection2.C b/PWG0/dNdEta/makeCorrection2.C
deleted file mode 100644 (file)
index 21521ac..0000000
+++ /dev/null
@@ -1,89 +0,0 @@
-/* $Id$ */
-
-//
-// Script to make correction maps for dndeta measurements using the
-// dNdEtaCorrection class.
-//
-// implementation with TSelector
-//
-
-#include "../CreateESDChain.C"
-#include "../PWG0Helper.C"
-
-void makeCorrection2(Char_t* dataDir, Int_t nRuns=20, Int_t offset = 0, Bool_t debug = kFALSE, Bool_t aProof = kFALSE, const Char_t* option = "", const Char_t* proofConnect = "proof01@lxb6046")
-{
-  if (aProof)
-    connectProof(proofConnect);
-
-  TString libraries("libEG;libGeom;libESD;libPWG0base;libVMC;libMinuit;libSTEER;libPWG0dep;libEVGEN;libFASTSIM;libmicrocern;libpdf;libpythia6;libEGPythia6;libAliPythia6");
-  TString packages("PWG0base;PWG0dep");
-
-  if (!prepareQuery(libraries, packages, kTRUE))
-    return;
-
-  gROOT->ProcessLine(".L CreateCuts.C");
-
-  AliESDtrackCuts* esdTrackCuts = CreateTrackCuts();
-  if (!esdTrackCuts)
-  {
-    printf("ERROR: esdTrackCuts could not be created\n");
-    return;
-  }
-
-  TList inputList;
-  inputList.Add(esdTrackCuts);
-
-  TChain* chain = CreateESDChain(dataDir, nRuns, offset);
-
-  TString selector("AlidNdEtaCorrectionSelector.cxx++");
-  if (debug != kFALSE)
-    selector += "g";
-
-  Int_t result = executeQuery(chain, &inputList, selector, option);
-}
-
-void VerifyCorrection(Int_t draw = 0, const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction", const char* analysisDir = "dndetaESD", const char* checkDir = "dndetaMC")
-{
-  gSystem->Load("libPWG0base");
-
-  AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
-  dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
-  dNdEtaCorrection->SetName("dndeta_correction");
-  dNdEtaCorrection->SetTitle("dndeta_correction");
-
-  TFile::Open(correctionMapFile);
-
-  dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis(analysisDir, analysisDir);
-  fdNdEtaAnalysis->LoadHistograms();
-
-  // correct with track2particle
-  TH3F* hist = fdNdEtaAnalysis->GetHistogram();
-  for (Int_t x=1; x<=hist->GetNbinsX(); ++x)
-    for (Int_t y=1; y<=hist->GetNbinsY(); ++y)
-      for (Int_t z=1; z<=hist->GetNbinsZ(); ++z)
-      {
-        Float_t correction = dNdEtaCorrection->GetTrack2ParticleCorrection(hist->GetXaxis()->GetBinCenter(x), hist->GetYaxis()->GetBinCenter(y), hist->GetZaxis()->GetBinCenter(z));
-        Float_t value = hist->GetBinContent(x, y, z);
-
-        hist->SetBinContent(x, y, z, correction * value);
-      }
-
-  fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3);
-
-  dNdEtaAnalysis* fdNdEtaAnalysisMC = new dNdEtaAnalysis(checkDir, checkDir);
-  fdNdEtaAnalysisMC->LoadHistograms();
-
-  fdNdEtaAnalysisMC->Finish(0, 0.3);
-
-  if (draw == 0)
-    fdNdEtaAnalysis->DrawHistograms();
-  else if (draw == 1)
-    fdNdEtaAnalysisMC->DrawHistograms();
-  else
-  {
-    for (Int_t i=0; i<dNdEtaAnalysis::kVertexBinning; ++i)
-      fdNdEtaAnalysis->GetdNdEtaHistogram(i)->Divide(fdNdEtaAnalysisMC->GetdNdEtaHistogram(i));
-
-    fdNdEtaAnalysis->DrawHistograms();
-  }
-}
diff --git a/PWG0/dNdEta/makeSystematics.C b/PWG0/dNdEta/makeSystematics.C
deleted file mode 100644 (file)
index a3e92cb..0000000
+++ /dev/null
@@ -1,40 +0,0 @@
-/* $Id$ */
-
-//
-// Script to run the creation of input for systematics
-//
-
-#include "../CreateESDChain.C"
-#include "../PWG0Helper.C"
-
-void makeSystematics(Char_t* dataDir, Int_t nRuns=20, Int_t offset = 0, Bool_t debug = kFALSE, Bool_t aProof = kFALSE, const Char_t* option = "")
-{
-  if (aProof)
-    connectProof("proof01@lxb6046");
-
-  TString libraries("libEG;libGeom;libESD;libVMC;libMinuit;libSTEER;libEVGEN;libFASTSIM;libmicrocern;libpdf;libpythia6;libEGPythia6;libAliPythia6;libPWG0base;libPWG0dep");
-  TString packages("PWG0base;PWG0dep");
-
-  if (!prepareQuery(libraries, packages, kTRUE))
-    return;
-
-  gROOT->ProcessLine(".L CreateCuts.C");
-
-  AliESDtrackCuts* esdTrackCuts = CreateTrackCuts();
-  if (!esdTrackCuts)
-  {
-    printf("ERROR: esdTrackCuts could not be created\n");
-    return;
-  }
-
-  TList inputList;
-  inputList.Add(esdTrackCuts);
-
-  TChain* chain = CreateESDChain(dataDir, nRuns, offset);
-
-  TString selector("AlidNdEtaSystematicsSelector.cxx++");
-  if (debug != kFALSE)
-    selector += "g";
-
-  Int_t result = executeQuery(chain, &inputList, selector, option);
-}
diff --git a/PWG0/dNdEta/plotsMultiplicity.C b/PWG0/dNdEta/plotsMultiplicity.C
deleted file mode 100644 (file)
index 2f6a61b..0000000
+++ /dev/null
@@ -1,2875 +0,0 @@
-/* $Id$ */
-
-//
-// plots for the note about multplicity measurements
-//
-
-#if !defined(__CINT__) || defined(__MAKECINT__)
-
-#include <TCanvas.h>
-#include <TPad.h>
-#include <TH1F.h>
-#include <TH2F.h>
-#include <TH3F.h>
-#include <TLine.h>
-#include <TF1.h>
-#include <TSystem.h>
-#include <TFile.h>
-#include <TLegend.h>
-#include <TStopwatch.h>
-#include <TROOT.h>
-#include <TGraph.h>
-#include <TMath.h>
-#include <TPaveText.h>
-#include <TImage.h>
-#include <TLatex.h>
-
-#include "AliMultiplicityCorrection.h"
-#include "AliCorrection.h"
-#include "AliCorrectionMatrix3D.h"
-
-#endif
-
-const char* correctionFile = "multiplicityMC_2M.root";
-const char* measuredFile   = "multiplicityMC_1M_3.root";
-Int_t etaRange = 3;
-Int_t displayRange = 200; // axis range
-Int_t ratioRange = 151;   // range to calculate difference
-Int_t longDisplayRange = 200;
-
-const char* correctionFileTPC = "multiplicityMC_TPC_1.4M.root";
-const char* measuredFileTPC   = "multiplicityMC_TPC_0.6M.root";
-Int_t etaRangeTPC = 1;
-
-void SetTPC()
-{
-  correctionFile = correctionFileTPC;
-  measuredFile = measuredFileTPC;
-  etaRange = etaRangeTPC;
-  displayRange = 100;
-  ratioRange = 76;
-  longDisplayRange = 100;
-}
-
-void Smooth(TH1* hist, Int_t windowWidth = 20)
-{
-  TH1* clone = (TH1*) hist->Clone("clone");
-  for (Int_t bin=2; bin<=clone->GetNbinsX(); ++bin)
-  {
-    Int_t min = TMath::Max(2, bin-windowWidth);
-    Int_t max = TMath::Min(clone->GetNbinsX(), bin+windowWidth);
-    Float_t average = clone->Integral(min, max) / (max - min + 1);
-
-    hist->SetBinContent(bin, average);
-    hist->SetBinError(bin, 0);
-  }
-
-  delete clone;
-}
-
-void responseMatrixPlot()
-{
-  gSystem->Load("libPWG0base");
-
-  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
-
-  TFile::Open(correctionFile);
-  mult->LoadHistograms("Multiplicity");
-
-  TH1* hist = mult->GetCorrelation(etaRange)->Project3D("zy");
-  hist->SetStats(kFALSE);
-
-  hist->SetTitle(";true multiplicity;measured multiplicity;Entries");
-  hist->GetXaxis()->SetRangeUser(0, longDisplayRange);
-  hist->GetYaxis()->SetRangeUser(0, longDisplayRange);
-
-  TCanvas* canvas = new TCanvas("c1", "c1", 800, 600);
-  canvas->SetRightMargin(0.15);
-  canvas->SetTopMargin(0.05);
-
-  gPad->SetLogz();
-  hist->Draw("COLZ");
-
-  canvas->SaveAs("responsematrix.eps");
-}
-
-TCanvas* DrawResultRatio(TH1* mcHist, TH1* result, TString epsName)
-{
-  // normalize unfolded result to mc hist
-  result->Scale(1.0 / result->Integral(2, 200));
-  result->Scale(mcHist->Integral(2, 200));
-
-  TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
-  canvas->Range(0, 0, 1, 1);
-
-  TPad* pad1 = new TPad(Form("%s_pad1", epsName.Data()), "", 0, 0.5, 0.98, 0.98);
-  pad1->Draw();
-
-  TPad* pad2 = new TPad(Form("%s_pad2", epsName.Data()), "", 0, 0.02, 0.98, 0.5);
-  pad2->Draw();
-
-  pad1->SetRightMargin(0.05);
-  pad2->SetRightMargin(0.05);
-
-  // no border between them
-  pad1->SetBottomMargin(0);
-  pad2->SetTopMargin(0);
-
-  pad1->cd();
-
-  mcHist->GetXaxis()->SetLabelSize(0.06);
-  mcHist->GetYaxis()->SetLabelSize(0.06);
-  mcHist->GetXaxis()->SetTitleSize(0.06);
-  mcHist->GetYaxis()->SetTitleSize(0.06);
-  mcHist->GetYaxis()->SetTitleOffset(0.6);
-
-  mcHist->GetXaxis()->SetRangeUser(0, displayRange);
-
-  mcHist->SetTitle(";true multiplicity;Entries");
-  mcHist->SetStats(kFALSE);
-
-  mcHist->DrawCopy("HIST E");
-  gPad->SetLogy();
-
-  result->SetLineColor(2);
-  result->DrawCopy("SAME HISTE");
-
-  TLegend* legend = new TLegend(0.6, 0.65, 0.95, 0.9);
-  legend->AddEntry(mcHist, "true distribution");
-  legend->AddEntry(result, "unfolded distribution");
-  legend->SetFillColor(0);
-  legend->Draw();
-
-  pad2->cd();
-  pad2->SetBottomMargin(0.15);
-
-  // calculate ratio
-  mcHist->Sumw2();
-  TH1* ratio = (TH1*) mcHist->Clone("ratio");
-  result->Sumw2();
-  ratio->Divide(ratio, result, 1, 1, "");
-  ratio->GetYaxis()->SetTitle("Ratio (true / unfolded)");
-  ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
-
-  ratio->DrawCopy();
-
-  // get average of ratio
-  Float_t sum = 0;
-  for (Int_t i=2; i<=ratioRange; ++i)
-  {
-    sum += TMath::Abs(ratio->GetBinContent(i) - 1);
-  }
-  sum /= ratioRange-1;
-
-  printf("Average (2..%d) of |ratio - 1| is %f\n", ratioRange, sum);
-
-  TLine* line = new TLine(0, 1, displayRange, 1);
-  line->SetLineWidth(2);
-  line->Draw();
-
-  line = new TLine(0, 1.1, displayRange, 1.1);
-  line->SetLineWidth(2);
-  line->SetLineStyle(2);
-  line->Draw();
-  line = new TLine(0, 0.9, displayRange, 0.9);
-  line->SetLineWidth(2);
-  line->SetLineStyle(2);
-  line->Draw();
-
-  canvas->Modified();
-
-  canvas->SaveAs(epsName);
-
-  return canvas;
-}
-
-TCanvas* Draw2ResultRatio(TH1* mcHist, TH1* result1, TH1* result2, TString epsName)
-{
-  // draws the 3 plots in the upper plot
-  // draws the ratio between result1 and result2 in the lower plot
-
-  // normalize unfolded result to mc hist
-  result1->Scale(1.0 / result1->Integral(2, 200));
-  result1->Scale(mcHist->Integral(2, 200));
-  result2->Scale(1.0 / result2->Integral(2, 200));
-  result2->Scale(mcHist->Integral(2, 200));
-
-  TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
-  canvas->Range(0, 0, 1, 1);
-
-  TPad* pad1 = new TPad(Form("%s_pad1", epsName.Data()), "", 0, 0.5, 0.98, 0.98);
-  pad1->Draw();
-
-  TPad* pad2 = new TPad(Form("%s_pad2", epsName.Data()), "", 0, 0.02, 0.98, 0.5);
-  pad2->Draw();
-
-  pad1->SetRightMargin(0.05);
-  pad2->SetRightMargin(0.05);
-
-  // no border between them
-  pad1->SetBottomMargin(0);
-  pad2->SetTopMargin(0);
-
-  pad1->cd();
-
-  mcHist->GetXaxis()->SetLabelSize(0.06);
-  mcHist->GetYaxis()->SetLabelSize(0.06);
-  mcHist->GetXaxis()->SetTitleSize(0.06);
-  mcHist->GetYaxis()->SetTitleSize(0.06);
-  mcHist->GetYaxis()->SetTitleOffset(0.6);
-
-  mcHist->GetXaxis()->SetRangeUser(0, displayRange);
-
-  mcHist->SetTitle(";true multiplicity;Entries");
-  mcHist->SetStats(kFALSE);
-
-  mcHist->DrawCopy("HIST E");
-  gPad->SetLogy();
-
-  result1->SetLineColor(2);
-  result1->DrawCopy("SAME HISTE");
-
-  result2->SetLineColor(4);
-  result2->DrawCopy("SAME HISTE");
-
-  TLegend* legend = new TLegend(0.55, 0.6, 0.95, 0.9);
-  legend->AddEntry(mcHist, "true distribution");
-  legend->AddEntry(result1, "unfolded distribution (syst)");
-  legend->AddEntry(result2, "unfolded distribution (normal)");
-  legend->SetFillColor(0);
-  legend->Draw();
-
-  pad2->cd();
-  pad2->SetBottomMargin(0.15);
-
-  result1->GetXaxis()->SetLabelSize(0.06);
-  result1->GetYaxis()->SetLabelSize(0.06);
-  result1->GetXaxis()->SetTitleSize(0.06);
-  result1->GetYaxis()->SetTitleSize(0.06);
-  result1->GetYaxis()->SetTitleOffset(0.6);
-
-  result1->GetXaxis()->SetRangeUser(0, displayRange);
-
-  result1->SetTitle(";true multiplicity;Entries");
-  result1->SetStats(kFALSE);
-
-  // calculate ratio
-  result1->Sumw2();
-  TH1* ratio = (TH1*) result1->Clone("ratio");
-  result2->Sumw2();
-  ratio->Divide(ratio, result2, 1, 1, "");
-  ratio->GetYaxis()->SetTitle("Ratio (syst / normal)");
-  ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
-
-  ratio->DrawCopy();
-
-  // get average of ratio
-  Float_t sum = 0;
-  for (Int_t i=2; i<=ratioRange; ++i)
-  {
-    sum += TMath::Abs(ratio->GetBinContent(i) - 1);
-  }
-  sum /= ratioRange-1;
-
-  printf("Average (2..%d) of |ratio - 1| is %f\n", ratioRange, sum);
-
-  TLine* line = new TLine(0, 1, displayRange, 1);
-  line->SetLineWidth(2);
-  line->Draw();
-
-  line = new TLine(0, 1.1, displayRange, 1.1);
-  line->SetLineWidth(2);
-  line->SetLineStyle(2);
-  line->Draw();
-  line = new TLine(0, 0.9, displayRange, 0.9);
-  line->SetLineWidth(2);
-  line->SetLineStyle(2);
-  line->Draw();
-
-  canvas->Modified();
-
-  canvas->SaveAs(epsName);
-
-  return canvas;
-}
-
-TCanvas* DrawRatio(TH1* result, Int_t nResultSyst, TH1** resultSyst, TString epsName, Bool_t firstMarker = kFALSE, const char** legendStrings = 0, Bool_t errors = kFALSE)
-{
-  // compares n results with first results. E.g. one gained with the default response, another with a changed one to study
-  // a systematic effect
-
-  // normalize results
-  result->Scale(1.0 / result->Integral(2, 200));
-
-  TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
-  canvas->SetTopMargin(0.05);
-  canvas->SetRightMargin(0.05);
-
-  result->GetXaxis()->SetRangeUser(0, displayRange);
-  result->GetYaxis()->SetRangeUser(0.55, 1.45);
-  result->SetStats(kFALSE);
-
-  // to get the axis how we want it
-  TH1* dummy = (TH1*) result->Clone("dummy");
-  dummy->Reset();
-  dummy->SetTitle(";true multiplicity;Ratio");
-  dummy->DrawCopy();
-  delete dummy;
-
-  Int_t colors[] = {1, 2, 4, 6, 7, 8, 9, 10};
-
-  TLegend* legend = new TLegend(0.2, 0.75, 0.35, 0.95);
-  legend->SetFillColor(0);
-
-  for (Int_t n=0; n<nResultSyst; ++n)
-  {
-    resultSyst[n]->Scale(1.0 / resultSyst[n]->Integral(2, 200));
-
-    // calculate ratio
-    TH1* ratio = (TH1*) result->Clone("ratio");
-    ratio->Divide(ratio, resultSyst[n], 1, 1, "");
-    ratio->GetXaxis()->SetRangeUser(1, displayRange);
-
-    if (firstMarker)
-      ratio->SetMarkerStyle(5);
-
-    ratio->SetLineColor(colors[n / 2]);
-    if ((n % 2))
-      ratio->SetLineStyle(2);
-
-    TString drawStr("SAME HIST");
-    if (n == 0 && firstMarker)
-      drawStr = "SAME P";
-    if (errors)
-      drawStr += " E";
-
-    ratio->DrawCopy(drawStr);
-
-    if (legendStrings && legendStrings[n])
-      legend->AddEntry(ratio, legendStrings[n]);
-
-    // get average of ratio
-    Float_t sum = 0;
-    for (Int_t i=2; i<=ratioRange; ++i)
-      sum += TMath::Abs(ratio->GetBinContent(i) - 1);
-    sum /= ratioRange-1;
-
-    printf("%d) Average (2..%d) of |ratio - 1| is %f\n", n, ratioRange, sum);
-  }
-
-  if (legendStrings)
-    legend->Draw();
-
-  TLine* line = new TLine(0, 1, displayRange, 1);
-  line->SetLineWidth(2);
-  line->Draw();
-
-  line = new TLine(0, 1.1, displayRange, 1.1);
-  line->SetLineWidth(2);
-  line->SetLineStyle(2);
-  line->Draw();
-  line = new TLine(0, 0.9, displayRange, 0.9);
-  line->SetLineWidth(2);
-  line->SetLineStyle(2);
-  line->Draw();
-
-  canvas->SaveAs(epsName);
-  canvas->SaveAs(Form("%s.gif", epsName.Data()));
-
-  return canvas;
-}
-
-TCanvas* DrawRatio(Int_t nResultSyst, TH1** mc, TH1** result, TString epsName, Bool_t smooth = kFALSE, Bool_t dashed = kFALSE)
-{
-  // draws the ratios of each mc to the corresponding result
-
-  TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
-  canvas->SetRightMargin(0.05);
-  canvas->SetTopMargin(0.05);
-
-  for (Int_t n=0; n<nResultSyst; ++n)
-  {
-    // normalize
-    result[n]->Scale(1.0 / result[n]->Integral(2, 200));
-    mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
-
-    result[n]->GetXaxis()->SetRangeUser(0, displayRange);
-    result[n]->SetStats(kFALSE);
-
-    // calculate ratio
-    TH1* ratio = (TH1*) result[n]->Clone("ratio");
-    ratio->Divide(mc[n], ratio, 1, 1, "B");
-
-    // SetRangeUser(1, ...) would be the same, but the 0 should be still on the axis...
-    ratio->SetBinContent(1, 1); ratio->SetBinError(1, 0);
-
-    if (smooth)
-      Smooth(ratio);
-
-    ratio->SetTitle(Form(";true multiplicity;Ratio (true / unfolded)%s", ((smooth) ? " (smoothed)" : "")));
-    ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
-
-    if (dashed)
-    {
-      ratio->SetLineColor((n/2)+1);
-      ratio->SetLineStyle((n%2)+1);
-    }
-    else
-      ratio->SetLineColor(n+1);
-
-    ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
-
-    // get average of ratio
-    Float_t sum = 0;
-    for (Int_t i=2; i<=ratioRange; ++i)
-      sum += TMath::Abs(ratio->GetBinContent(i) - 1);
-    sum /= ratioRange-1;
-
-    printf("%d) Average (2..%d) of |ratio - 1| is %f\n", n, ratioRange, sum);
-  }
-
-  TLine* line = new TLine(0, 1, displayRange, 1);
-  line->SetLineWidth(2);
-  line->Draw();
-
-  line = new TLine(0, 1.1, displayRange, 1.1);
-  line->SetLineWidth(2);
-  line->SetLineStyle(2);
-  line->Draw();
-  line = new TLine(0, 0.9, displayRange, 0.9);
-  line->SetLineWidth(2);
-  line->SetLineStyle(2);
-  line->Draw();
-
-  canvas->Modified();
-
-  canvas->SaveAs(epsName);
-  canvas->SaveAs(Form("%s.gif", epsName.Data()));
-
-  return canvas;
-}
-
-TCanvas* DrawRatioDeduct(TH1* mcBase, TH1* resultBase, Int_t nResultSyst, TH1** mc, TH1** result, TString epsName)
-{
-  // draws the ratios of each mc to the corresponding result
-  // deducts from each ratio the ratio of mcBase / resultBase
-
-  // normalize
-  resultBase->Scale(1.0 / resultBase->Integral(2, 200));
-  mcBase->Scale(1.0 / mcBase->Integral(2, 200));
-
-  // calculate ratio
-  TH1* ratioBase = (TH1*) resultBase->Clone("ratioBase");
-  ratioBase->Divide(mcBase, ratioBase, 1, 1, "B");
-
-  TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
-  canvas->SetRightMargin(0.05);
-  canvas->SetTopMargin(0.05);
-
-  for (Int_t n=0; n<nResultSyst; ++n)
-  {
-    // normalize
-    result[n]->Scale(1.0 / result[n]->Integral(2, 200));
-    mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
-
-    result[n]->GetXaxis()->SetRangeUser(0, displayRange);
-    result[n]->SetStats(kFALSE);
-
-    // calculate ratio
-    TH1* ratio = (TH1*) result[n]->Clone("ratio");
-    ratio->Divide(mc[n], ratio, 1, 1, "B");
-    ratio->Add(ratioBase, -1);
-
-    ratio->SetTitle(";true multiplicity;Ratio_{syst} (t/u) - Ratio (t/u)");
-    ratio->GetYaxis()->SetRangeUser(-1, 1);
-    ratio->SetLineColor(n+1);
-    ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
-
-    // get average of ratio
-    Float_t sum = 0;
-    for (Int_t i=2; i<=ratioRange; ++i)
-      sum += TMath::Abs(ratio->GetBinContent(i));
-    sum /= ratioRange-1;
-
-    printf("%d) Average (2..%d) of |ratio - ratioBase| is %f\n", n, ratioRange, sum);
-  }
-
-  TLine* line = new TLine(0, 0, displayRange, 0);
-  line->SetLineWidth(2);
-  line->Draw();
-
-  line = new TLine(0, 0.1, displayRange, 0.1);
-  line->SetLineWidth(2);
-  line->SetLineStyle(2);
-  line->Draw();
-  line = new TLine(0, -0.1, displayRange, -0.1);
-  line->SetLineWidth(2);
-  line->SetLineStyle(2);
-  line->Draw();
-
-  canvas->Modified();
-
-  canvas->SaveAs(epsName);
-  canvas->SaveAs(Form("%s.gif", epsName.Data()));
-
-  return canvas;
-}
-
-TCanvas* DrawRatioDeductSmooth(TH1* mcBase, TH1* resultBase, Int_t nResultSyst, TH1** mc, TH1** result, TString epsName)
-{
-  // draws the ratios of each mc to the corresponding result
-  // deducts from each ratio the ratio of mcBase / resultBase
-  // smoothens the ratios by a sliding window
-
-  // normalize
-  resultBase->Scale(1.0 / resultBase->Integral(2, 200));
-  mcBase->Scale(1.0 / mcBase->Integral(2, 200));
-
-  // calculate ratio
-  TH1* ratioBase = (TH1*) resultBase->Clone("ratioBase");
-  ratioBase->Divide(mcBase, ratioBase, 1, 1, "B");
-
-  TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
-  canvas->SetRightMargin(0.05);
-  canvas->SetTopMargin(0.05);
-
-  for (Int_t n=0; n<nResultSyst; ++n)
-  {
-    // normalize
-    result[n]->Scale(1.0 / result[n]->Integral(2, 200));
-    mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
-
-    result[n]->GetXaxis()->SetRangeUser(0, displayRange);
-    result[n]->SetStats(kFALSE);
-
-    // calculate ratio
-    TH1* ratio = (TH1*) result[n]->Clone("ratio");
-    ratio->Divide(mc[n], ratio, 1, 1, "B");
-    ratio->Add(ratioBase, -1);
-
-    //new TCanvas; ratio->DrawCopy();
-    // clear 0 bin
-    ratio->SetBinContent(1, 0); ratio->SetBinError(1, 0);
-
-    Smooth(ratio);
-
-    //ratio->SetLineColor(1); ratio->DrawCopy("SAME");
-
-    canvas->cd();
-    ratio->SetTitle(";true multiplicity;Ratio_{syst} (t/u) - Ratio (t/u) (smoothed)");
-    ratio->GetYaxis()->SetRangeUser(-0.3, 0.3);
-    ratio->SetLineColor((n / 2)+1);
-    ratio->SetLineStyle((n % 2)+1);
-    ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
-
-    // get average of ratio
-    Float_t sum = 0;
-    for (Int_t i=2; i<=150; ++i)
-      sum += TMath::Abs(ratio->GetBinContent(i));
-    sum /= 149;
-
-    printf("%d) Average (2..150) of |ratio - ratioBase| is %f\n", n, sum);
-  }
-
-  TLine* line = new TLine(0, 0, displayRange, 0);
-  line->SetLineWidth(2);
-  line->Draw();
-
-  line = new TLine(0, 0.1, displayRange, 0.1);
-  line->SetLineWidth(2);
-  line->SetLineStyle(2);
-  line->Draw();
-  line = new TLine(0, -0.1, displayRange, -0.1);
-  line->SetLineWidth(2);
-  line->SetLineStyle(2);
-  line->Draw();
-
-  canvas->Modified();
-
-  canvas->SaveAs(epsName);
-  canvas->SaveAs(Form("%s.gif", epsName.Data()));
-
-  return canvas;
-}
-
-void DrawResiduals(TH1* measured, TH1* unfoldedFolded, const char* epsName)
-{
-  // normalize
-  unfoldedFolded->Scale(1.0 / unfoldedFolded->Integral(2, 200));
-  unfoldedFolded->Scale(measured->Integral(2, 200));
-
-  TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
-  canvas->Range(0, 0, 1, 1);
-
-  TPad* pad1 = new TPad(Form("%s_pad1", epsName), "", 0, 0.5, 0.98, 0.98);
-  pad1->Draw();
-  pad1->SetGridx();
-  pad1->SetGridy();
-
-  TPad* pad2 = new TPad(Form("%s_pad2", epsName), "", 0, 0.02, 0.98, 0.5);
-  pad2->Draw();
-  pad2->SetGridx();
-  pad2->SetGridy();
-
-  TPad* pad3 = new TPad(Form("%s_pad3", epsName), "", 0.15, 0.5, 0.35, 0.75);
-  pad3->SetGridx();
-  pad3->SetGridy();
-  pad3->SetRightMargin(0.05);
-  pad3->SetTopMargin(0.05);
-  pad3->Draw();
-
-  pad1->SetRightMargin(0.05);
-  pad2->SetRightMargin(0.05);
-
-  // no border between them
-  pad1->SetBottomMargin(0);
-  pad2->SetTopMargin(0);
-
-  pad1->cd();
-
-  measured->GetXaxis()->SetLabelSize(0.06);
-  measured->GetYaxis()->SetLabelSize(0.06);
-  measured->GetXaxis()->SetTitleSize(0.06);
-  measured->GetYaxis()->SetTitleSize(0.06);
-  measured->GetYaxis()->SetTitleOffset(0.6);
-
-  measured->GetXaxis()->SetRangeUser(0, 150);
-
-  measured->SetTitle(";measured multiplicity;Entries");
-  measured->SetStats(kFALSE);
-
-  measured->DrawCopy("HIST");
-  gPad->SetLogy();
-
-  unfoldedFolded->SetMarkerStyle(5);
-  unfoldedFolded->SetMarkerColor(2);
-  unfoldedFolded->SetLineColor(0);
-  unfoldedFolded->DrawCopy("SAME P");
-
-  TLegend* legend = new TLegend(0.6, 0.65, 0.95, 0.9);
-  legend->AddEntry(measured, "measured distribution");
-  legend->AddEntry(unfoldedFolded, "R #otimes unfolded distribution");
-  legend->SetFillColor(0);
-  legend->Draw();
-
-  pad2->cd();
-  pad2->SetBottomMargin(0.15);
-
-  // calculate ratio
-  measured->Sumw2();
-  TH1* residual = (TH1*) measured->Clone("residual");
-  unfoldedFolded->Sumw2();
-
-  residual->Add(unfoldedFolded, -1);
-
-  // projection
-  TH1* residualHist = new TH1F("residualHist", ";", 15, -3, 3);
-
-  for (Int_t i=1; i<=residual->GetNbinsX(); ++i)
-  {
-    if (measured->GetBinError(i) > 0)
-    {
-      residual->SetBinContent(i, residual->GetBinContent(i) / measured->GetBinError(i));
-      residual->SetBinError(i, 1);
-
-      residualHist->Fill(residual->GetBinContent(i));
-    }
-    else
-    {
-      residual->SetBinContent(i, 0);
-      residual->SetBinError(i, 0);
-    }
-  }
-
-  residual->GetYaxis()->SetTitle("Residuals   1/e (M - R #otimes U)");
-  residual->GetYaxis()->SetRangeUser(-4.5, 4.5);
-  residual->DrawCopy();
-
-  TLine* line = new TLine(-0.5, 0, 150.5, 0);
-  line->SetLineWidth(2);
-  line->Draw();
-
-  pad3->cd();
-  residualHist->SetStats(kFALSE);
-  residualHist->GetXaxis()->SetLabelSize(0.08);
-  residualHist->Fit("gaus");
-  residualHist->Draw();
-
-  canvas->Modified();
-  canvas->SaveAs(canvas->GetName());
-
-  //const char* epsName2 = "proj.eps";
-  //TCanvas* canvas = new TCanvas(epsName2, epsName2, 800, 600);
-  //canvas->SetGridx();
-  //canvas->SetGridy();
-
-  //canvas->SaveAs(canvas->GetName());
-}
-
-void bayesianExample()
-{
-  TStopwatch watch;
-  watch.Start();
-
-  gSystem->Load("libPWG0base");
-
-  TFile::Open(correctionFile);
-  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
-  mult->LoadHistograms("Multiplicity");
-
-  TFile::Open(measuredFile);
-  AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
-  mult2->LoadHistograms("Multiplicity");
-
-  mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
-
-  mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
-
-  TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
-  TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
-
-  //mult->DrawComparison("bayesianExample", etaRange, kFALSE, kTRUE, mcHist, kTRUE);
-  DrawResultRatio(mcHist, result, "bayesianExample.eps");
-
-  //Printf("KolmogorovTest says PROB = %f", mcHist->KolmogorovTest(result, "D"));
-  //Printf("Chi2Test says PROB = %f", mcHist->Chi2Test(result));
-
-  // draw residual plot
-
-  // TODO take out efficiency correction if other than AliMultiplicityCorrection::kTrVtx
-  TH2* convoluted = mult->CalculateMultiplicityESD(result, etaRange);
-  TH1* convolutedProj = convoluted->ProjectionY("convolutedProj", -1, -1, "e");
-
-  TH1* measured = mult2->GetMultiplicityESD(etaRange)->ProjectionY("measured");
-
-  DrawResiduals(measured, convolutedProj, "bayesianResiduals.eps");
-
-  watch.Stop();
-  watch.Print();
-}
-
-void chi2FluctuationResult()
-{
-  gSystem->Load("libPWG0base");
-
-  TFile::Open(correctionFile);
-  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
-  mult->LoadHistograms("Multiplicity");
-
-  TFile::Open(measuredFile);
-  AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
-  mult2->LoadHistograms("Multiplicity");
-
-  mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
-  mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0);
-  mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
-
-  TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
-  //TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
-
-  mult->DrawComparison("MinuitChi2", etaRange, kFALSE, kTRUE, mcHist, kTRUE);
-
-  TCanvas* canvas = (TCanvas*) gROOT->FindObject("MinuitChi2_DrawComparison_3");
-  canvas->SaveAs("chi2FluctuationResult.eps");
-}
-
-void chi2Example()
-{
-  TStopwatch watch;
-  watch.Start();
-
-  gSystem->Load("libPWG0base");
-
-  TFile::Open(correctionFile);
-  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
-  mult->LoadHistograms("Multiplicity");
-
-  TFile::Open(measuredFile);
-  AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
-  mult2->LoadHistograms("Multiplicity");
-
-  mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
-  mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
-  mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
-
-  TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
-  TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
-
-  DrawResultRatio(mcHist, result, "chi2Example.eps");
-
-  watch.Stop();
-  watch.Print();
-}
-
-void chi2ExampleTPC()
-{
-  TStopwatch watch;
-  watch.Start();
-
-  gSystem->Load("libPWG0base");
-
-  TFile::Open(correctionFileTPC);
-  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
-  mult->LoadHistograms("Multiplicity");
-
-  TFile::Open(measuredFileTPC);
-  AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
-  mult2->LoadHistograms("Multiplicity");
-
-  mult->SetMultiplicityESD(etaRangeTPC, mult2->GetMultiplicityESD(etaRangeTPC));
-  mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
-  mult->ApplyMinuitFit(etaRangeTPC, kFALSE, AliMultiplicityCorrection::kTrVtx);
-
-  TH1* mcHist = mult2->GetMultiplicityVtx(etaRangeTPC)->ProjectionY("mymc");
-  TH1* result = mult->GetMultiplicityESDCorrected(etaRangeTPC);
-
-  DrawResultRatio(mcHist, result, "chi2ExampleTPC.eps");
-
-  watch.Stop();
-  watch.Print();
-}
-
-void bayesianNBD()
-{
-  gSystem->Load("libPWG0base");
-  TFile::Open("multiplicityMC_3M.root");
-
-  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
-  mult->LoadHistograms("Multiplicity");
-
-  TFile::Open("multiplicityMC_3M_NBD.root");
-  AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
-  mult2->LoadHistograms("Multiplicity");
-
-  mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
-  mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1.0, 100);
-
-  TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
-
-  mcHist->Sumw2();
-  TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
-
-  //mult->DrawComparison("bayesianNBD", etaRange, kFALSE, kTRUE, mcHist);
-  DrawResultRatio(mcHist, result, "bayesianNBD.eps");
-}
-
-void minimizationNBD()
-{
-  gSystem->Load("libPWG0base");
-  TFile::Open("multiplicityMC_3M.root");
-
-  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
-  mult->LoadHistograms("Multiplicity");
-
-  TFile::Open("multiplicityMC_3M_NBD.root");
-  AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
-  mult2->LoadHistograms("Multiplicity");
-
-  mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
-  mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
-  mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
-
-  TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
-
-  mcHist->Sumw2();
-  TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
-
-  //mult->DrawComparison("minimizationNBD", etaRange, kFALSE, kTRUE, mcHist);
-  DrawResultRatio(mcHist, result, "minimizationNBD.eps");
-}
-
-void minimizationInfluenceAlpha()
-{
-  gSystem->Load("libPWG0base");
-
-  TFile::Open(measuredFile);
-  AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
-  mult2->LoadHistograms("Multiplicity");
-
-  TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
-  mcHist->Scale(1.0 / mcHist->Integral());
-  mcHist->GetXaxis()->SetRangeUser(0, 200);
-  mcHist->SetStats(kFALSE);
-  mcHist->SetTitle(";true multiplicity;P_{N}");
-
-  TCanvas* canvas = new TCanvas("minimizationInfluenceAlpha", "minimizationInfluenceAlpha", 1000, 300);
-  canvas->Divide(3, 1);
-
-  TFile::Open("eval-2M-1M/EvaluateChi2MethodDetail.root");
-
-  TH1* hist1 = (TH1*) gFile->Get("MinuitChi2_00_2_100.000000");
-  TH1* hist2 = (TH1*) gFile->Get("MinuitChi2_03_2_100000.000000");
-  TH1* hist3 = (TH1*) gFile->Get("MinuitChi2_06_2_100000000.000000");
-
-  mcHist->Rebin(2);  mcHist->Scale(0.5);
-  hist1->Rebin(2);   hist1->Scale(0.5);
-  hist2->Rebin(2);   hist2->Scale(0.5);
-  hist3->Rebin(2);   hist3->Scale(0.5);
-
-  mcHist->GetXaxis()->SetRangeUser(0, 200);
-
-  canvas->cd(1);
-  gPad->SetLogy();
-  mcHist->Draw();
-  hist1->SetMarkerStyle(5);
-  hist1->SetMarkerColor(2);
-  hist1->Draw("SAME PE");
-
-  canvas->cd(2);
-  gPad->SetLogy();
-  mcHist->Draw();
-  hist2->SetMarkerStyle(5);
-  hist2->SetMarkerColor(2);
-  hist2->Draw("SAME PE");
-
-  canvas->cd(3);
-  gPad->SetLogy();
-  mcHist->Draw();
-  hist3->SetMarkerStyle(5);
-  hist3->SetMarkerColor(2);
-  hist3->Draw("SAME PE");
-
-  canvas->SaveAs("minimizationInfluenceAlpha.eps");
-}
-
-void NBDFit()
-{
-  gSystem->Load("libPWG0base");
-
-  TFile::Open(correctionFile);
-  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
-  mult->LoadHistograms("Multiplicity");
-
-  TH1* fCurrentESD = mult->GetMultiplicityVtx(etaRange)->ProjectionY();
-  fCurrentESD->Sumw2();
-  fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
-
-  TF1* func = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])");
-  func->SetParNames("scaling", "averagen", "k");
-  func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
-  func->SetParLimits(1, 0.001, 1000);
-  func->SetParLimits(2, 0.001, 1000);
-  func->SetParameters(fCurrentESD->GetMaximum() * 100, 10, 2);
-
-  TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500);
-  lognormal->SetParNames("scaling", "mean", "sigma");
-  lognormal->SetParameters(1, 1, 1);
-  lognormal->SetParLimits(0, 0, 10);
-  lognormal->SetParLimits(1, 0, 100);
-  lognormal->SetParLimits(2, 1e-3, 10);
-
-  TCanvas* canvas = new TCanvas("c1", "c1", 700, 400);
-  fCurrentESD->SetStats(kFALSE);
-  fCurrentESD->GetYaxis()->SetTitleOffset(1.3);
-  fCurrentESD->SetTitle(";true multiplicity (N);P_{N}");
-  fCurrentESD->Draw("HIST");
-  fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
-  fCurrentESD->Fit(func, "W0", "", 0, 50);
-  func->SetRange(0, 100);
-  func->Draw("SAME");
-  printf("chi2 = %f\n", func->GetChisquare());
-
-  fCurrentESD->Fit(lognormal, "W0", "", 0.01, 100);
-  lognormal->SetLineColor(2);
-  lognormal->SetLineStyle(2);
-  lognormal->SetRange(0, 100);
-  lognormal->Draw("SAME");
-
-  canvas->SaveAs("NBDFit.eps");
-}
-
-void DifferentSamples()
-{
-  // data generated by runMultiplicitySelector.C DifferentSamples
-
-  const char* name = "DifferentSamples";
-
-  TFile* file = TFile::Open(Form("%s.root", name));
-
-  TCanvas* canvas = new TCanvas(name, name, 800, 600);
-  canvas->Divide(2, 2);
-
-  for (Int_t i=0; i<4; ++i)
-  {
-    canvas->cd(i+1);
-    gPad->SetTopMargin(0.05);
-    gPad->SetRightMargin(0.05);
-    TH1* chi2Result = (TH1*) file->Get(Form("chi2Result_%d", i));
-    TH1* bayesResult = (TH1*) file->Get(Form("bayesResult_%d", i));
-    TH1* mc = (TH1*) file->Get(Form("mc_%d", i));
-    mc->Sumw2();
-
-    chi2Result->Divide(chi2Result, mc, 1, 1, "");
-    bayesResult->Divide(bayesResult, mc, 1, 1, "");
-
-    chi2Result->SetTitle(";true multiplicity;unfolded measured/MC");
-    chi2Result->GetXaxis()->SetRangeUser(0, 150);
-    chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
-    chi2Result->GetYaxis()->SetTitleOffset(1.2);
-    chi2Result->SetLineColor(1);
-    chi2Result->SetStats(kFALSE);
-
-    bayesResult->SetStats(kFALSE);
-    bayesResult->SetLineColor(2);
-
-    chi2Result->DrawCopy("HIST");
-    bayesResult->DrawCopy("SAME HIST");
-
-    TLine* line = new TLine(0, 1, 150, 1);
-    line->Draw();
-  }
-
-  canvas->SaveAs(Form("%s.eps", canvas->GetName()));
-}
-
-void StartingConditions()
-{
-  // data generated by runMultiplicitySelector.C StartingConditions
-
-  const char* name = "StartingConditions";
-
-  TFile* file = TFile::Open(Form("%s.root", name));
-
-  TCanvas* canvas = new TCanvas(name, name, 800, 400);
-  canvas->Divide(2, 1);
-
-  TH1* mc = (TH1*) file->Get("mc");
-  mc->Sumw2();
-  mc->Scale(1.0 / mc->Integral());
-
-  //Int_t marker[] = {24, 25, 26, 27, 28, 2, 3, 4, 5};
-
-  TLegend* legend = new TLegend(0.6, 0.7, 0.95, 0.95);
-  legend->SetFillColor(0);
-
-  const char* names[] = { "True", "Measured 1", "Measured 2", "Measured 3", "NBD", "Flat" };
-
-  for (Int_t i=0; i<6; ++i)
-  {
-    Int_t id = i;
-    if (id > 2)
-      id += 2;
-
-    TH1* chi2Result = (TH1*) file->Get(Form("chi2Result_%d", id));
-    TH1* bayesResult = (TH1*) file->Get(Form("bayesResult_%d", id));
-
-    chi2Result->Divide(chi2Result, mc, 1, 1, "");
-    bayesResult->Divide(bayesResult, mc, 1, 1, "");
-
-    chi2Result->SetTitle("a) #chi^{2} minimization;true multiplicity;unfolded / MC");
-    chi2Result->GetXaxis()->SetRangeUser(0, 150);
-    chi2Result->GetYaxis()->SetRangeUser(0.8, 1.2);
-    chi2Result->GetYaxis()->SetTitleOffset(1.5);
-    //chi2Result->SetMarkerStyle(marker[i]);
-    chi2Result->SetLineColor(i+1);
-    chi2Result->SetMarkerColor(i+1);
-    chi2Result->SetStats(kFALSE);
-
-    bayesResult->SetTitle("b) Bayesian method;true multiplicity;unfolded / MC");
-    bayesResult->GetXaxis()->SetRangeUser(0, 150);
-    bayesResult->GetYaxis()->SetRangeUser(0.8, 1.2);
-    bayesResult->GetYaxis()->SetTitleOffset(1.5);
-    bayesResult->SetStats(kFALSE);
-    //bayesResult->SetLineColor(2);
-    bayesResult->SetLineColor(i+1);
-
-    canvas->cd(1);
-    gPad->SetLeftMargin(0.12);
-    chi2Result->DrawCopy((i == 0) ? "HIST" : "HIST SAME");
-
-    canvas->cd(2);
-    gPad->SetLeftMargin(0.12);
-    bayesResult->DrawCopy((i == 0) ? "HIST" : "HIST SAME");
-
-    //TLine* line = new TLine(0, 1, 150, 1);
-    //line->Draw();
-
-    legend->AddEntry(chi2Result, names[i]);
-  }
-
-  canvas->cd(1);
-  legend->Draw();
-  canvas->SaveAs(Form("%s.eps", canvas->GetName()));
-}
-
-void StatisticsPlot()
-{
-  const char* name = "StatisticsPlot";
-
-  TFile* file = TFile::Open(Form("%s.root", name));
-
-  TCanvas* canvas = new TCanvas(name, name, 600, 400);
-
-  TGraph* fitResultsChi2 = (TGraph*) file->Get("fitResultsChi2");
-  fitResultsChi2->SetTitle(";number of measured events;P_{1}");
-  fitResultsChi2->GetYaxis()->SetRangeUser(0, 2);
-  fitResultsChi2->Draw("AP");
-
-  TF1* f = new TF1("f", "[0]/x", 1, 1e4);
-  fitResultsChi2->Fit(f);
-
-  canvas->SaveAs(Form("%s.eps", canvas->GetName()));
-
-  TH1* mc[5];
-  TH1* result[5];
-
-  const char* plotname = "chi2Result";
-
-  name = "StatisticsPlotRatios";
-  canvas = new TCanvas(name, name, 600, 400);
-
-  for (Int_t i=0; i<5; ++i)
-  {
-    mc[i] = (TH1*) file->Get(Form("mc_%d", i));
-    result[i] = (TH1*) file->Get(Form("%s_%d", plotname, i));
-
-    result[i]->SetLineColor(i+1);
-    result[i]->Draw(((i == 0) ? "" : "SAME"));
-  }
-}
-
-void SystematicLowEfficiency()
-{
-  gSystem->Load("libPWG0base");
-
-  TFile::Open(correctionFile);
-  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
-  mult->LoadHistograms("Multiplicity");
-
-  // calculate result with systematic effect
-  TFile::Open("multiplicityMC_100k_1_loweff.root");
-  AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
-  mult2->LoadHistograms("Multiplicity");
-
-  mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
-  mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
-  mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
-
-  TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
-  TH1* result1 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1");
-
-  DrawResultRatio(mcHist, result1, "SystematicLowEfficiencyLow.eps");
-
-  // calculate normal result
-  TFile::Open("multiplicityMC_100k_1.root");
-  mult2->LoadHistograms("Multiplicity");
-
-  mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
-  mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
-
-  TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
-
-  DrawResultRatio(mcHist, result2, "SystematicLowEfficiencyOK.eps");
-
-  Draw2ResultRatio(mcHist, result1, result2, "SystematicLowEfficiency.eps");
-}
-
-void SystematicMisalignment()
-{
-  gSystem->Load("libPWG0base");
-
-  TFile::Open(correctionFile);
-  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
-  mult->LoadHistograms("Multiplicity");
-
-  TFile::Open("multiplicityMC_100k_fullmis.root");
-  AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
-  mult2->LoadHistograms("Multiplicity");
-
-  mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
-  mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
-  mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
-
-  TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
-  TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
-
-  DrawResultRatio(mcHist, result, "SystematicMisalignment.eps");
-}
-
-void SystematicMisalignmentTPC()
-{
-  gSystem->Load("libPWG0base");
-
-  SetTPC();
-
-  TFile::Open(correctionFile);
-  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
-  mult->LoadHistograms("Multiplicity");
-
-  TFile::Open("multiplicityMC_TPC_100k_fullmis.root");
-  AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
-  mult2->LoadHistograms("Multiplicity");
-
-  mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
-  mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
-  mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
-
-  TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
-  TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
-
-  DrawResultRatio(mcHist, result, "SystematicMisalignmentTPC.eps");
-}
-
-void EfficiencySpecies()
-{
-  gSystem->Load("libPWG0base");
-
-  Int_t marker[] = {24, 25, 26};
-  Int_t color[] = {1, 2, 4};
-
-  // SPD TPC
-  const char* fileName[] = { "multiplicityMC_400k_syst.root", "multiplicityMC_TPC_4kfiles_syst.root" };
-  Float_t etaRange[] = {0.49, 0.9};
-  const char* titles[] = { "SPD Tracklets", "TPC Tracks" };
-
-  TCanvas* canvas = new TCanvas("EfficiencySpecies", "EfficiencySpecies", 1000, 500);
-  canvas->Divide(2, 1);
-
-  for (Int_t loop=0; loop<1; ++loop)
-  {
-    Printf("%s", fileName[loop]);
-
-    AliCorrection* correction[4];
-
-    canvas->cd(loop+1);
-
-    gPad->SetGridx();
-    gPad->SetGridy();
-    gPad->SetRightMargin(0.05);
-    //gPad->SetTopMargin(0.05);
-
-    TLegend* legend = new TLegend(0.7, 0.4, 0.85, 0.6);
-    legend->SetFillColor(0);
-    legend->SetEntrySeparation(0.2);
-
-    Float_t below = 0;
-    Float_t total = 0;
-
-    TFile* file = TFile::Open(fileName[loop]);
-    if (!file)
-    {
-      Printf("Could not open %s", fileName[loop]);
-      return;
-    }
-
-    Float_t sumGen = 0;
-    Float_t sumMeas = 0;
-
-    for (Int_t i=0; i<3; ++i)
-    {
-      Printf("correction %d", i);
-
-      TString name; name.Form("correction_%d", i);
-      correction[i] = new AliCorrection(name, name);
-      correction[i]->LoadHistograms();
-
-      TH3* gene = correction[i]->GetTrackCorrection()->GetGeneratedHistogram();
-      TH3* meas = correction[i]->GetTrackCorrection()->GetMeasuredHistogram();
-
-      // limit vtx axis
-      gene->GetXaxis()->SetRangeUser(-3.9, 3.9);
-      meas->GetXaxis()->SetRangeUser(-3.9, 3.9);
-
-      // empty over/underflow bin in eta, setting range to +-2 is not enough because this is the maximum range, Project3D takes them into account then (might be a bug)
-      /*for (Int_t x = 1; x <= gene->GetNbinsX(); x++)
-        for (Int_t z = 1; z <= gene->GetNbinsZ(); z++)
-        {
-          gene->SetBinContent(x, 0, z, 0);
-          gene->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
-          meas->SetBinContent(x, 0, z, 0);
-          meas->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
-        }*/
-
-      // limit eta axis
-      gene->GetYaxis()->SetRangeUser(-etaRange[loop], etaRange[loop]);
-      meas->GetYaxis()->SetRangeUser(-etaRange[loop], etaRange[loop]);
-
-      TH1* genePt = gene->Project3D(Form("z_%d", i));
-      TH1* measPt = meas->Project3D(Form("z_%d", i));
-
-      genePt->Sumw2();
-      measPt->Sumw2();
-
-      sumGen += genePt->Integral();
-      sumMeas += measPt->Integral();
-
-      TH1* effPt = (TH1*) genePt->Clone(Form("effPt_%d", i));
-      effPt->Reset();
-      effPt->Divide(measPt, genePt, 1, 1, "B");
-
-      Int_t bin = 0;
-      for (bin=20; bin>=1; bin--)
-      {
-        if (effPt->GetBinContent(bin) < 0.5)
-          break;
-      }
-
-      Printf("Eff. below 50%% at bin %d, i.e. %.3f GeV/c", bin, effPt->GetXaxis()->GetBinUpEdge(bin));
-
-      Float_t fraction = genePt->Integral(1, bin) / genePt->Integral();
-      Printf("%.4f of the particles are below that momentum", fraction);
-
-      below += genePt->Integral(1, bin);
-      total += genePt->Integral();
-
-      effPt->SetLineColor(color[i]);
-      effPt->SetMarkerColor(color[i]);
-      effPt->SetMarkerStyle(marker[i]);
-
-      effPt->GetXaxis()->SetRangeUser(0.06, 1);
-      effPt->GetYaxis()->SetRangeUser(0, 1);
-
-      effPt->GetYaxis()->SetTitleOffset(1.2);
-
-      effPt->SetStats(kFALSE);
-      effPt->SetTitle(titles[loop]);
-      effPt->GetYaxis()->SetTitle("Efficiency");
-
-      effPt->DrawCopy((i == 0) ? "" : "SAME");
-
-      legend->AddEntry(effPt, ((i == 0) ? "#pi^{#pm}" : ((i == 1) ? "K^{#pm}" : "p,#bar{p}")));
-    }
-
-    Printf("In total %.4f of the particles are below their effective pt cut off", (Float_t) below / total);
-
-    Printf("%f measured, %f generated, effiency: %f", sumGen, sumMeas, sumMeas / sumGen);
-
-    legend->Draw();
-  }
-
-  canvas->SaveAs(Form("%s.eps", canvas->GetName()));
-}
-
-void ParticleSpeciesComparison1(Bool_t chi2 = kTRUE, const char* fileNameMC = "multiplicityMC_400k_syst_species.root", const char* fileNameESD = "multiplicityMC_100k_syst.root")
-{
-  gSystem->Load("libPWG0base");
-
-  TFile::Open(fileNameESD);
-  TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", etaRange));
-  TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", etaRange));
-
-  TH1* results[10];
-
-  // loop over cases (normal, enhanced/reduced ratios)
-  Int_t nMax = 7;
-  for (Int_t i = 0; i<nMax; ++i)
-  {
-    TString folder;
-    folder.Form("Multiplicity_%d", i);
-
-    AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folder, folder);
-
-    TFile::Open(fileNameMC);
-    mult->LoadHistograms();
-
-    mult->SetMultiplicityESD(etaRange, hist);
-
-    if (chi2)
-    {
-      mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
-      mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
-      //mult->DrawComparison(Form("ParticleSpeciesComparison_MinuitChi2_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist"));
-    }
-    else
-    {
-      mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
-      //mult->DrawComparison(Form("ParticleSpeciesComparison_Bayesian_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist2"));
-    }
-
-    //Float_t averageRatio = 0;
-    //mult->GetComparisonResults(0, 0, 0, &averageRatio);
-
-    results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
-
-    //Printf("Case %d. Average ratio is %f", i, averageRatio);
-  }
-
-  DrawResultRatio(hist2->ProjectionY("mymchist", -1, -1, "e"), results[0], "ParticleSpeciesComparison1_1.eps");
-
-  TH1* mc = hist2->ProjectionY("mymchist2", -1, -1, "e");
-
-  for (Int_t i=1; i<=results[0]->GetNbinsX(); i++)
-  {
-    results[0]->SetBinError(i, 0);
-    mc->SetBinError(i, 0);
-  }
-
-  const char* legendStrings[] = { "#pi^{#pm}", 0, "K^{#pm}", 0, "p,#bar{p}", 0 };
-
-  DrawRatio(results[0], nMax-1, results+1, "ParticleSpeciesComparison1_2.eps", kFALSE, legendStrings);
-
-  //not valid: draw chi2 uncertainty on top!
-  /*TFile::Open("bayesianUncertainty_400k_100k_syst.root");
-  TH1* errorHist = (TH1*) gFile->Get("errorBoth");
-  errorHist->SetLineColor(1);
-  errorHist->SetLineWidth(2);
-  TH1* errorHist2 = (TH1*) errorHist->Clone("errorHist2");
-  for (Int_t i=1; i<=errorHist->GetNbinsX(); i++)
-  {
-    errorHist->SetBinContent(i, errorHist->GetBinContent(i) + 1);
-    errorHist2->SetBinContent(i, 1 - errorHist2->GetBinContent(i));
-  }
-
-  errorHist->DrawCopy("SAME");
-  errorHist2->DrawCopy("SAME");*/
-
-  //canvas->SaveAs(canvas->GetName());
-
-  DrawRatio(mc, nMax, results, "ParticleSpeciesComparison1_3.eps", kTRUE, 0);
-
-  //errorHist->DrawCopy("SAME");
-  //errorHist2->DrawCopy("SAME");
-
-  //canvas2->SaveAs(canvas2->GetName());
-}
-
-/*void ParticleSpeciesComparison2()
-{
-  gSystem->Load("libPWG0base");
-
-  const char* fileNameMC = "multiplicityMC_400k_syst.root";
-  const char* fileNameESD = "out.root"; // based on multiplicityMC_100k_syst.root
-  Bool_t chi2 = 0;
-
-  TFile::Open(fileNameMC);
-  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
-  mult->LoadHistograms();
-
-  TH1* mc[10];
-  TH1* results[10];
-
-  // loop over cases (normal, enhanced/reduced ratios)
-  Int_t nMax = 7;
-  for (Int_t i = 0; i<nMax; ++i)
-  {
-    TString folder;
-    folder.Form("Multiplicity_%d", i);
-
-    AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection(folder, folder);
-
-    TFile::Open(fileNameESD);
-    mult2->LoadHistograms();
-
-    mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
-
-    if (chi2)
-    {
-      mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
-      mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
-      //mult->DrawComparison(Form("ParticleSpeciesComparison_MinuitChi2_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist"));
-    }
-    else
-    {
-      mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
-      //mult->DrawComparison(Form("ParticleSpeciesComparison_Bayesian_%d", i), etaRange, kFALSE, kTRUE, hist2->ProjectionY("mymchist2"));
-    }
-
-    //Float_t averageRatio = 0;
-    //mult->GetComparisonResults(0, 0, 0, &averageRatio);
-
-    results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
-
-    TH2F* hist2 = mult2->GetMultiplicityVtx(etaRange);
-    mc[i] = (TH1*) hist2->ProjectionY(Form("mymchist_%d", i), -1, -1, "e");
-
-    //TString fileName; fileName.Form("ParticleSpeciesComparison2_%d.eps", i);
-    //DrawResultRatio(hist2->ProjectionY("mymchist", -1, -1, "e"), results[i], fileName);
-
-    //Printf("Case %d. Average ratio is %f", i, averageRatio);
-  }
-
-  DrawRatio(nMax, mc, results, "ParticleSpeciesComparison2.eps");
-}*/
-
-TH1* Invert(TH1* eff)
-{
-  // calculate corr = 1 / eff
-
-  TH1* corr = (TH1*) eff->Clone(Form("%s_invert", eff->GetName()));
-  corr->Reset();
-
-  for (Int_t i=1; i<=eff->GetNbinsX(); i++)
-  {
-    if (eff->GetBinContent(i) > 0)
-    {
-      corr->SetBinContent(i, 1.0 / eff->GetBinContent(i));
-      corr->SetBinError(i, eff->GetBinError(i) / eff->GetBinContent(i) * corr->GetBinContent(i));
-    }
-  }
-
-  return corr;
-}
-
-void TriggerVertexCorrection()
-{
-  //
-  // plots the correction performed on the unfolded spectrum to gain the spectrum for the full inelastic sample
-  //
-
-  gSystem->Load("libPWG0base");
-
-  TFile::Open(correctionFile);
-  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
-  mult->LoadHistograms("Multiplicity");
-
-  TH1* corrINEL = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kINEL));
-  TH1* corrMB   = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kMB));
-
-  TCanvas* canvas = new TCanvas("TriggerVertexCorrection", "TriggerVertexCorrection", 800, 600);
-
-  corrINEL->SetStats(kFALSE);
-  corrINEL->GetXaxis()->SetRangeUser(0, 30);
-  corrINEL->GetYaxis()->SetRangeUser(0, 5);
-  corrINEL->SetTitle(";true multiplicity;correction factor");
-  corrINEL->SetMarkerStyle(22);
-  corrINEL->Draw("PE");
-
-  corrMB->SetStats(kFALSE);
-  corrMB->SetLineColor(2);
-  corrMB->SetMarkerStyle(25);
-  corrMB->SetMarkerColor(2);
-  corrMB->Draw("SAME PE");
-
-  TLegend* legend = new TLegend(0.3, 0.5, 0.85, 0.65);
-  legend->SetFillColor(0);
-  legend->AddEntry(corrINEL, "correction to inelastic sample");
-  legend->AddEntry(corrMB, "correction to minimum bias sample");
-
-  legend->Draw();
-
-  canvas->SaveAs(Form("%s.eps", canvas->GetName()));
-}
-
-void StatisticalUncertainty(Int_t methodType, Bool_t mc = kFALSE)
-{
-  gSystem->Load("libPWG0base");
-
-  TFile::Open(correctionFile);
-  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
-  mult->LoadHistograms("Multiplicity");
-
-  TFile::Open(measuredFile);
-  AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
-  mult2->LoadHistograms("Multiplicity");
-
-  mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
-
-  TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
-
-  TH1* errorResponse = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorResponse");
-
-  TH1* errorMeasured = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kFALSE, ((mc) ? mcHist : 0))->Clone("errorMeasured");
-  TH1* errorBoth = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorBoth");
-
-  if (!mc)
-  {
-    TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
-    DrawResultRatio(mcHist, result, "StatisticalUncertainty2.eps");
-  }
-
-  TCanvas* canvas = new TCanvas("StatisticalUncertainty", "StatisticalUncertainty", 600, 400);
-  canvas->SetGridx();
-  canvas->SetGridy();
-  canvas->SetRightMargin(0.05);
-  canvas->SetTopMargin(0.05);
-
-  errorResponse->SetLineColor(1);
-  errorResponse->GetXaxis()->SetRangeUser(0, 200);
-  errorResponse->GetYaxis()->SetRangeUser(0, 0.3);
-  errorResponse->SetStats(kFALSE);
-  errorResponse->SetTitle(";true multiplicity;Uncertainty");
-
-  errorResponse->Draw();
-
-  errorMeasured->SetLineColor(2);
-  errorMeasured->Draw("SAME");
-
-  errorBoth->SetLineColor(4);
-  errorBoth->Draw("SAME");
-
-  Printf("Average errorResponse: %f", errorResponse->Integral(2, 150) / 149);
-  Printf("Average errorMeasured: %f", errorMeasured->Integral(2, 150) / 149);
-  Printf("Average errorBoth: %f", errorBoth->Integral(2, 150) / 149);
-
-  canvas->SaveAs(Form("%s.eps", canvas->GetName()));
-
-  TFile* file = new TFile(Form("%s.root", canvas->GetName()), "RECREATE");
-  errorResponse->Write();
-  errorMeasured->Write();
-  errorBoth->Write();
-  file->Close();
-}
-
-void StatisticalUncertaintyCompare(const char* det = "SPD")
-{
-  TFile* file1 = TFile::Open(Form("StatisticalUncertainty%sBayesian.root", det));
-  TH1* errorResponse = (TH1*) file1->Get("errorResponse");
-  TH1* errorMeasured = (TH1*) file1->Get("errorMeasured");
-  TH1* errorBoth = (TH1*) file1->Get("errorBoth");
-
-  TString str;
-  str.Form("StatisticalUncertaintyCompare%s", det);
-
-  TCanvas* canvas = new TCanvas(str, str, 600, 400);
-  canvas->SetGridx();
-  canvas->SetGridy();
-  canvas->SetRightMargin(0.05);
-  canvas->SetTopMargin(0.05);
-
-  errorResponse->SetLineColor(1);
-  errorResponse->GetXaxis()->SetRangeUser(1, (strcmp(det, "TPC") ? 200 : 100));
-  errorResponse->GetYaxis()->SetRangeUser(0, 0.3);
-  errorResponse->SetStats(kFALSE);
-  errorResponse->GetYaxis()->SetTitleOffset(1.2);
-  errorResponse->SetTitle(";true multiplicity;#sigma(U-T)/T");
-
-  errorResponse->Draw();
-
-  errorMeasured->SetLineColor(2);
-  errorMeasured->Draw("SAME");
-
-  errorBoth->SetLineColor(4);
-  errorBoth->Draw("SAME");
-
-  TFile* file2 = TFile::Open(Form("StatisticalUncertainty%sChi2.root", det));
-  TH1* errorBoth2 = (TH1*) file2->Get("errorBoth");
-
-  errorBoth2->SetLineColor(4);
-  errorBoth2->SetLineStyle(2);
-  errorBoth2->Draw("SAME");
-
-  TLegend* legend = new TLegend(0.2, 0.6, 0.6, 0.9);
-  legend->SetFillColor(0);
-  legend->AddEntry(errorResponse, "response matrix (Bayesian)");
-  legend->AddEntry(errorMeasured, "measured (Bayesian)");
-  legend->AddEntry(errorBoth, "both (Bayesian)");
-  legend->AddEntry(errorBoth2, "both (#chi^{2} minimization)");
-  legend->Draw();
-
-  canvas->SaveAs(Form("%s.eps", canvas->GetName()));
-}
-
-void EfficiencyComparison(Int_t eventType = 2, Bool_t uncertainty = kTRUE)
-{
- const char* files[] = { "multiplicityMC_400k_syst_nd.root", "multiplicityMC_400k_syst_sd.root", "multiplicityMC_400k_syst_dd.root", "multiplicityMC_400k_syst_xsection.root" };
-
-  gSystem->Load("libPWG0base");
-
-  TCanvas* canvas = new TCanvas("EfficiencyComparison", "EfficiencyComparison", 800, 500);
-  canvas->SetGridx();
-  canvas->SetGridy();
-  canvas->SetRightMargin(0.05);
-  canvas->SetTopMargin(0.05);
-
-  AliMultiplicityCorrection* data[4];
-  TH1* effArray[4];
-
-  Int_t markers[] = { 24, 25, 26, 5 };
-  Int_t colors[] = { 1, 2, 3, 4 };
-
-  TLegend* legend = new TLegend(0.45, 0.45, 0.9, 0.7);
-  legend->SetFillColor(0);
-
-  TH1* effError = 0;
-
-  for (Int_t i=0; i<4; ++i)
-  {
-    TString name;
-    name.Form("Multiplicity_%d", i);
-
-    TFile::Open(files[i]);
-    data[i] = new AliMultiplicityCorrection(name, name);
-
-    if (i < 3)
-    {
-      data[i]->LoadHistograms("Multiplicity");
-    }
-    else
-      data[i]->LoadHistograms("Multiplicity_0");
-
-    TH1* eff = (TH1*) data[i]->GetEfficiency(etaRange, (AliMultiplicityCorrection::EventType) eventType)->Clone(Form("eff_%d", i));
-    effArray[i] = eff;
-
-    eff->GetXaxis()->SetRangeUser(0, 15);
-    eff->GetYaxis()->SetRangeUser(0, 1.1);
-    eff->SetStats(kFALSE);
-    eff->SetTitle(";true multiplicity;Efficiency");
-    eff->SetLineColor(colors[i]);
-    eff->SetMarkerColor(colors[i]);
-    eff->SetMarkerStyle(markers[i]);
-
-    if (i == 3)
-    {
-      for (Int_t bin=1; bin<=eff->GetNbinsX(); bin++)
-        eff->SetBinError(bin, 0);
-
-      // loop over cross section combinations
-      for (Int_t j=1; j<7; ++j)
-      {
-        AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multtmp", "Multtmp");
-        mult->LoadHistograms(Form("Multiplicity_%d", j));
-
-        TH1* eff2 = mult->GetEfficiency(etaRange, (AliMultiplicityCorrection::EventType) eventType);
-
-        for (Int_t bin=1; bin<=eff->GetNbinsX(); bin++)
-        {
-          // TODO we could also do asymmetric errors here
-          Float_t deviation = TMath::Abs(eff->GetBinContent(bin) - eff2->GetBinContent(bin));
-
-          eff->SetBinError(bin, TMath::Max(eff->GetBinError(bin), (Double_t) deviation));
-        }
-      }
-
-      for (Int_t bin=1; bin<=20; bin++)
-        if (eff->GetBinContent(bin) > 0)
-          Printf("Bin %d: Error: %.2f", bin, 100.0 * eff->GetBinError(bin) / eff->GetBinContent(bin));
-      
-      if (uncertainty) {
-       effError = (TH1*) eff->Clone("effError");
-       effError->Reset();
-
-       for (Int_t bin=2; bin<=eff->GetNbinsX(); bin++)
-         if (eff->GetBinContent(bin) > 0)
-           effError->SetBinContent(bin, 10.0 * eff->GetBinError(bin) / eff->GetBinContent(bin));
-
-       effError->SetLineColor(1);
-       effError->SetMarkerStyle(1);
-       effError->DrawCopy("SAME HIST");
-      }
-    }
-
-    eff->SetBinContent(1, 0);
-    eff->SetBinError(1, 0);
-
-    canvas->cd();
-    if (i == 0)
-    {
-      eff->DrawCopy("P");
-    }
-    else
-      eff->DrawCopy("SAME P");
-
-    legend->AddEntry(eff, (((i == 0) ? "non diffractive" : ((i == 1) ? "single diffractive" : ((i == 2) ? "double diffractive" : "Pythia combined")))));
-  }
-
-  if (uncertainty)
-    legend->AddEntry(effError, "relative syst. uncertainty #times 10");
-
-  legend->Draw();
-
-  canvas->SaveAs(Form("%s.eps", canvas->GetName()));
-}
-
-void ModelDependencyPlot()
-{
-  gSystem->Load("libPWG0base");
-
-  TFile::Open("multiplicityMC_3M.root");
-  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
-  mult->LoadHistograms("Multiplicity");
-
-  TH2* proj = (TH2*) mult->GetCorrelation(3)->Project3D("zy");
-
-  TCanvas* canvas = new TCanvas("ModelDependencyPlot", "ModelDependencyPlot", 800, 400);
-  canvas->SetGridx();
-  canvas->SetGridy();
-  //canvas->SetRightMargin(0.05);
-  //canvas->SetTopMargin(0.05);
-
-  canvas->Divide(2, 1);
-
-  canvas->cd(2);
-  gPad->SetLogy();
-  Int_t selectedMult = 30;
-  Int_t yMax = 200000;
-
-  TH1* full = proj->ProjectionX("full");
-  TH1* selected = proj->ProjectionY("selected", proj->GetXaxis()->FindBin(selectedMult), proj->GetXaxis()->FindBin(selectedMult)); 
-
-  full->SetStats(kFALSE);
-  full->GetXaxis()->SetRangeUser(0, 200);
-  full->GetYaxis()->SetRangeUser(5, yMax);
-  full->SetTitle(";multiplicity");
-
-  selected->SetLineColor(0);
-  selected->SetMarkerColor(2);
-  selected->SetMarkerStyle(7);
-
-  full->Draw();
-  selected->Draw("SAME P");
-
-  TLegend* legend = new TLegend(0.5, 0.65, 0.85, 0.85);
-  legend->SetFillColor(0);
-  legend->AddEntry(full, "true");
-  legend->AddEntry(selected, "measured");
-  legend->Draw();
-  TLine* line = new TLine(selectedMult, 5, selectedMult, yMax);
-  line->SetLineWidth(2);
-  line->Draw();
-
-  canvas->cd(1);
-  gPad->SetLogy();
-
-  full = proj->ProjectionY("full2");
-  selected = proj->ProjectionX("selected2", proj->GetYaxis()->FindBin(selectedMult), proj->GetYaxis()->FindBin(selectedMult));
-
-  full->SetStats(kFALSE);
-  full->GetXaxis()->SetRangeUser(0, 200);
-  full->GetYaxis()->SetRangeUser(5, yMax);
-  full->SetTitle(";multiplicity");
-
-  full->SetLineColor(0);
-  full->SetMarkerColor(2);
-  full->SetMarkerStyle(7);
-
-  full->Draw("P");
-  selected->Draw("SAME");
-
-  legend->Draw();
-
-  line = new TLine(selectedMult, 5, selectedMult, yMax);
-  line->SetLineWidth(2);
-  line->Draw();
-
-  canvas->SaveAs(Form("%s.eps", canvas->GetName()));
-}
-
-void SystematicpTSpectrum()
-{
-  gSystem->Load("libPWG0base");
-
-  TFile::Open("multiplicityMC_400k_syst_ptspectrum.root");
-  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
-  mult->LoadHistograms("Multiplicity");
-
-  TFile::Open("multiplicityMC_100k_syst.root");
-  AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
-  mult2->LoadHistograms("Multiplicity");
-
-  mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
-  mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
-  mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
-
-  TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
-  TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
-
-  DrawResultRatio(mcHist, result, "SystematicpTSpectrum.eps");
-}
-
-// to be deleted
-/*void covMatrix(Bool_t mc = kTRUE)
-{
-  gSystem->Load("libPWG0base");
-
-  TFile::Open(correctionFile);
-  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
-  mult->LoadHistograms("Multiplicity");
-
-  TFile::Open(measuredFile);
-  AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
-  mult2->LoadHistograms("Multiplicity");
-
-  mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
-
-  TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
-
-  mult->BayesianStatisticsEffect(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kTRUE, 1, 100, ((mc) ? mcHist : 0));
-}*/
-
-Double_t FitPtFunc(Double_t *x, Double_t *par)
-{
-  Double_t xx = x[0];
-
-  Float_t val1 = par[1] + par[2] * xx + par[3] * xx * xx;
-  Float_t val2 = TMath::Exp(par[4] + par[5] * xx);
-
-  const Float_t kTransitionWidth = 0;
-
-  // power law part
-  if (xx < par[0] - kTransitionWidth)
-  {
-    return val1;
-  }
-  /*else if (xx < par[0] + kTransitionWidth)
-  {
-    // smooth transition
-    Float_t factor = (xx - par[0] + kTransitionWidth) / kTransitionWidth / 2;
-    return (1 - factor) * val1 + factor * val2;
-  }*/
-  else
-  {
-    return val2;
-  }
-}
-
-void FitPtNew(const char* fileName = "TruePt14TeV.root")
-{
-  gSystem->Load("libANALYSIS");
-  gSystem->Load("libPWG0base");
-
-  TFile::Open(fileName);
-
-  TH1* genePt = (TH1*) gFile->Get("fHistPt");
-  genePt->Sumw2();
-
-  // normalize by bin width
-  for (Int_t x=1; x<genePt->GetNbinsX(); x++)
-    genePt->SetBinContent(x, genePt->GetBinContent(x) / genePt->GetBinWidth(x));
-
-  genePt->GetXaxis()->SetRangeUser(0.05, 2.0);
-
-  genePt->Scale(1.0 / genePt->Integral());
-
-  TF1* func = new TF1("func", "[0]*TMath::Exp([1]*x*x)", 0.001, 100);
-  //func->SetLineColor(2);
-  func->SetParameters(1, -1);
-
-  genePt->SetMarkerStyle(25);
-  genePt->SetTitle("");
-  genePt->SetStats(kFALSE);
-  genePt->GetYaxis()->SetRangeUser(1e-4, genePt->GetMaximum() * 1.2);
-  //func->Draw("SAME");
-
-  genePt->Fit(func, "0", "", 0.05, 1);
-
-  new TCanvas;
-  genePt->DrawCopy("P");
-  func->SetRange(0.02, 8);
-  func->DrawCopy("SAME");
-  gPad->SetLogy();
-}
-
-void FitPt(const char* fileName = "firstplots100k_truept.root")
-{
-  gSystem->Load("libPWG0base");
-
-  TFile::Open(fileName);
-
-  /*
-  // merge corrections
-  AliCorrection* correction[4];
-  TList list;
-
-  for (Int_t i=0; i<4; ++i)
-  {
-    Printf("correction %d", i);
-
-    TString name; name.Form("correction_%d", i);
-    correction[i] = new AliCorrection(name, name);
-    correction[i]->LoadHistograms();
-
-    if (i > 0)
-      list.Add(correction[i]);
-  }
-
-  correction[0]->Merge(&list);
-
-  TH3* gene = correction[0]->GetTrackCorrection()->GetGeneratedHistogram();
-
-  // limit vtx, eta axis
-  gene->GetXaxis()->SetRangeUser(-5.9, 5.9);
-  gene->GetYaxis()->SetRangeUser(-1.99, 0.99);
-
-  TH1* genePt = gene->Project3D("z");*/
-  TH1* genePt = (TH1*) gFile->Get("fdNdpTTrue");
-  if (!genePt)
-    genePt = (TH1*) gFile->Get("fHistPt");
-  genePt->Sumw2();
-
-  //genePt->Scale(1.0 / genePt->Integral());
-
-  // normalize by bin width
-  for (Int_t x=1; x<genePt->GetNbinsX(); x++)
-    genePt->SetBinContent(x, genePt->GetBinContent(x) / genePt->GetBinWidth(x));
-
-  /// genePt->GetXaxis()->GetBinCenter(x));
-
-  genePt->GetXaxis()->SetRangeUser(0, 7.9);
-  //genePt->GetYaxis()->SetTitle("a.u.");
-
-  //TF1* func = new TF1("func", "[0]*TMath::Exp([1]*x*x)", 0.001, 100);
-  TF1* func = new TF1("func", "[0]*TMath::Exp([1]*x)+[2]/(1+(x*[4])**[3])", 0.001, 100);
-  //func->SetLineColor(2);
-  func->SetParameters(1, -1, 1, 1, 1);
-  func->SetParLimits(3, 1, 10);
-  func->SetParLimits(4, 0, 10);
-
-  //TF1* func = new TF1("func", "[1]*x**[0]", 0.001, 100);
-
-  //TF1* func = new TF1("func", &FitPtFunc, 0, 2, 6);
-  //func->SetParameters(0.3, -2.34909e-01, 1.54394e+01, -3.04134e+01, 1.41912e+00, -2.79284e+00);
-  //func->FixParameter(0, 0.314);
-  //func->SetParLimits(0, 0.1, 0.3);
-
-  genePt->SetMarkerStyle(25);
-  genePt->SetTitle("");
-  genePt->SetStats(kFALSE);
-  genePt->GetYaxis()->SetRangeUser(1e-4, genePt->GetMaximum() * 1.2);
-  //func->Draw("SAME");
-
-  // fit only exp. part
-  func->SetParameters(1, -1);
-  func->FixParameter(2, 0);
-  func->FixParameter(3, 1);
-  func->FixParameter(4, 1);
-  genePt->Fit(func, "0", "", 0.2, 1);
-
-  new TCanvas;
-  genePt->DrawCopy("P");
-  func->SetRange(0.02, 8);
-  func->DrawCopy("SAME");
-  gPad->SetLogy();
-
-  // now fix exp. parameters and fit second part
-  Double_t param0 = func->GetParameter(0);
-  Double_t param1 = func->GetParameter(1);
-  func->SetParameters(0, -1, 1, 1, 1);
-  func->FixParameter(0, 0);
-  func->FixParameter(1, -1);
-  func->ReleaseParameter(2);
-  func->ReleaseParameter(3);
-  func->ReleaseParameter(4);
-  func->SetParLimits(3, 1, 10);
-  func->SetParLimits(4, 0, 10);
-
-  genePt->Fit(func, "0", "", 1.5, 4);
-
-  new TCanvas;
-  genePt->DrawCopy("P");
-  func->SetRange(0.02, 8);
-  func->DrawCopy("SAME");
-  gPad->SetLogy();
-
-  // fit both
-  func->SetParameter(0, param0);
-  func->SetParameter(1, param1);
-  func->ReleaseParameter(0);
-  func->ReleaseParameter(1);
-
-  new TCanvas;
-  genePt->DrawCopy("P");
-  func->SetRange(0.02, 5);
-  func->DrawCopy("SAME");
-  gPad->SetLogy();
-
-  genePt->Fit(func, "0", "", 0.2, 4);
-
-  TCanvas* canvas = new TCanvas("FitPt", "FitPt", 800, 400);
-  canvas->Divide(2, 1);
-  canvas->cd(1);
-
-  gPad->SetGridx();
-  gPad->SetGridy();
-  gPad->SetLeftMargin(0.13);
-  gPad->SetRightMargin(0.05);
-  gPad->SetTopMargin(0.05);
-
-  genePt->GetXaxis()->SetRangeUser(0, 4.9);
-  genePt->GetYaxis()->SetRangeUser(1e-2, 1e4);
-  genePt->GetYaxis()->SetTitleOffset(1.4);
-  genePt->GetXaxis()->SetTitleOffset(1.1);
-  genePt->DrawCopy("P");
-  func->SetRange(0.02, 5);
-  func->DrawCopy("SAME");
-  gPad->SetLogy();
-
-  canvas->cd(2);
-
-  TH1* genePtClone = (TH1*) genePt->Clone("genePtClone");
-  genePtClone->Reset();
-  genePtClone->DrawCopy("P");
-
-  gPad->SetGridx();
-  gPad->SetGridy();
-  gPad->SetLeftMargin(0.13);
-  gPad->SetRightMargin(0.05);
-  gPad->SetTopMargin(0.05);
-
-  func->DrawCopy("SAME");
-  gPad->SetLogy();
-
-  canvas->SaveAs(Form("%s.eps", canvas->GetName()));
-
-  TH1* first = (TH1*) func->GetHistogram()->Clone("first");
-
-  TCanvas* canvas2 = new TCanvas("FitPt2", "FitPt2", 600, 400);
-
-  TFile* file = TFile::Open("ptspectrum_fit.root", "RECREATE");
-
-  for (Int_t param=0; param<5; param++)
-  {
-    for (Int_t sign=0; sign<2; sign++)
-    {
-      TF1* func2 = (TF1*) func->Clone(Form("func_%d_%d", param, sign));  // new TF1(Form("func_%d_%d", param, sign), &FitPtFunc, 0, 2, 6);
-      func2->SetParameters(func->GetParameters());
-      //TF1* func2 = (TF1*) func->Clone(); // SetParameter after this does not work
-
-      Float_t factor = ((sign == 0) ? 0.9 : 1.1);
-      func2->SetParameter(param, func2->GetParameter(param) * factor);
-      //func2->Print();
-
-      canvas->cd(2);
-      func2->SetLineWidth(1);
-      func2->SetLineColor(2);
-      func2->DrawCopy("SAME");
-
-      canvas2->cd();
-      TH1* second = func2->GetHistogram();
-      second->Divide(first);
-      second->SetLineColor(param + 1);
-      second->GetYaxis()->SetRangeUser(0, 2);
-      second->DrawCopy((param == 0 && sign == 0) ? "" : "SAME");
-      second->Clone(Form("ptspectrum_%d_%d", param, sign))->Write();
-    }
-  }
-
-  canvas->SaveAs(Form("%s.eps", canvas->GetName()));
-  canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
-
-  file->Close();
-}
-
-void DrawSystematicpT()
-{
-  TFile* file = TFile::Open("SystematicpT.root");
-
-  TH1* mcHist2 = (TH1*) file->Get("mymc_unity");
-  TH1* result2 = (TH1*) file->Get("result_unity");
-
-  TH1* mcHist[12];
-  TH1* result[12];
-
-  Int_t nParams = 5;
-
-  for (Int_t id=0; id<nParams*2; ++id)
-  {
-    mcHist[id] = (TH1*) file->Get(Form("mymc_%d_%d.root", id / 2, id % 2));
-    result[id] = (TH1*) file->Get(Form("result_%d_%d.root", id / 2, id % 2));
-  }
-
-  DrawResultRatio(mcHist2, result2, "SystematicpT_OK.eps");
-
-  //DrawRatioDeduct(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
-
-  DrawRatio(nParams*2, mcHist, result, "SystematicpT_Ratios.eps", kTRUE, kTRUE);
-
-  //DrawRatioDeductSmooth(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
-
-  // does not make sense: mc is different
-  //Draw2ResultRatio(mcHist, result1, result2, "SystematicpT.eps");
-}
-
-void SystematicpT(Bool_t chi2 = 1)
-{
-  gSystem->Load("libPWG0base");
-
-  TFile::Open("ptspectrum900.root");
-  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
-  mult->LoadHistograms("Multiplicity");
-
-  AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
-
-  TH1* mcHist[12];
-  TH1* result[12];
-
-  Int_t nParams = 5;
-
-  for (Int_t param=0; param<nParams; param++)
-  {
-    for (Int_t sign=0; sign<2; sign++)
-    {
-      // calculate result with systematic effect
-      TFile::Open(Form("ptspectrum100_%d_%d.root", param, sign));
-      mult2->LoadHistograms("Multiplicity");
-
-      mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
-
-      if (chi2)
-      {
-        mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
-        mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
-      }
-      else
-        mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0);
-
-      Int_t id = param * 2 + sign;
-
-      mcHist[id] = mult2->GetMultiplicityVtx(etaRange)->ProjectionY(Form("mymc_%d_%d.root", param, sign));
-      result[id] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d_%d.root", param, sign));
-
-      TString tmp; tmp.Form("SystematicpT_%d_%d.eps", param, sign);
-      DrawResultRatio(mcHist[id], result[id], tmp);
-    }
-  }
-
-  // calculate normal result
-  TFile::Open("ptspectrum100_1.root");
-  mult2->LoadHistograms("Multiplicity");
-  TH1* mcHist2 = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc_unity");
-
-  mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
-
-  if (chi2)
-  {
-    mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
-  }
-  else
-    mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
-
-  TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result_unity");
-
-  TFile* file = TFile::Open("SystematicpT.root", "RECREATE");
-  mcHist2->Write();
-  result2->Write();
-  for (Int_t id=0; id<nParams*2; ++id)
-  {
-    mcHist[id]->Write();
-    result[id]->Write();
-  }
-  file->Close();
-
-  DrawSystematicpT();
-}
-
-void DrawSystematicpT2()
-{
-  //displayRange = 200;
-
-  // read from file
-  TFile* file = TFile::Open("SystematicpT2.root");
-  TH1* mcHist = (TH1*) file->Get("mymc");
-  TH1* result[12];
-  result[0] = (TH1*) file->Get("result_unity");
-  Int_t nParams = 5;
-  for (Int_t id=0; id<nParams*2; ++id)
-    result[id+1] = (TH1*) file->Get(Form("result_%d_%d", id / 2, id % 2));
-
-  DrawResultRatio((TH1*) mcHist->Clone(), (TH1*) result[0]->Clone(), "SystematicpT_OK.eps");
-  DrawRatio(mcHist, nParams*2+1, result, "SystematicpT_Ratios_MC.eps", kTRUE);
-  DrawRatio(result[0], nParams*2, result+1, "SystematicpT_Ratios.eps");
-}
-
-void SystematicpT2(Bool_t tpc = kTRUE, Bool_t chi2 = kTRUE)
-{
-  gSystem->Load("libPWG0base");
-
-  if (tpc)
-  {
-    SetTPC();
-    TFile::Open("multiplicityMC_TPC_0.6M_syst_pt_unity.root");
-  }
-  else
-    TFile::Open("ptspectrum100_1.root");
-
-  AliMultiplicityCorrection* measured = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
-  measured->LoadHistograms("Multiplicity");
-  TH1* mcHist = measured->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
-
-  TH1* result[12];
-
-  Int_t nParams = 5;
-
-  // -1 = unity change, 0...4 parameters
-  for (Int_t id=-1; id<nParams*2; id++)
-  {
-    Int_t param = id / 2;
-    Int_t sign = id % 2;
-
-    TString idStr;
-    if (id == -1)
-    {
-      idStr = "unity";
-    }
-    else
-      idStr.Form("%d_%d", param, sign);
-
-    // calculate result with systematic effect
-    if (tpc)
-    {
-      TFile::Open(Form("multiplicityMC_TPC_1.3M_syst_pt_%s.root", idStr.Data()));
-    }
-    else
-      TFile::Open(Form("ptspectrum900_%s.root", idStr.Data()));
-
-    AliMultiplicityCorrection* response = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
-    response->LoadHistograms("Multiplicity");
-
-    response->SetMultiplicityESD(etaRange, measured->GetMultiplicityESD(etaRange));
-
-    if (chi2)
-    {
-      response->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
-      response->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
-    }
-    else
-      response->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0);
-
-    result[id+1] = (TH1*) response->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%s", idStr.Data()));
-
-    TString tmp; tmp.Form("SystematicpT_%s.eps", idStr.Data());
-    DrawResultRatio(mcHist, result[id+1], tmp);
-  }
-
-  TFile* file = TFile::Open("SystematicpT2.root", "RECREATE");
-  mcHist->Write();
-  for (Int_t id=0; id<nParams*2+1; ++id)
-    result[id]->Write();
-  file->Close();
-
-  DrawSystematicpT2();
-}
-
-void SystematicpTCutOff(Bool_t chi2 = kTRUE)
-{
-  // only needed for TPC
-  SetTPC();
-
-  gSystem->Load("libPWG0base");
-
-  TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_unity.root");
-  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
-  mult->LoadHistograms("Multiplicity");
-
-  TFile::Open("multiplicityMC_TPC_0.6M_syst_pt_unity.root");
-  AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
-  mult2->LoadHistograms("Multiplicity");
-
-  // "normal" result
-  mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
-
-  if (chi2)
-  {
-    mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
-    mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
-  }
-  else
-    mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
-
-  TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
-  TH1* result1 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1");
-
-  TH1* syst[2];
-
-  // change of pt spectrum (down)
-  TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_red.root");
-  AliMultiplicityCorrection* mult3 = new AliMultiplicityCorrection("Multiplicity3", "Multiplicity3");
-  mult3->LoadHistograms("Multiplicity");
-  mult3->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
-  if (chi2)
-  {
-    mult3->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
-    mult3->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
-  }
-  else
-    mult3->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
-  syst[0] = (TH1*) mult3->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
-
-  // change of pt spectrum (up)
-  TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_inc.root");
-  AliMultiplicityCorrection* mult4 = new AliMultiplicityCorrection("Multiplicity4", "Multiplicity4");
-  mult4->LoadHistograms("Multiplicity");
-  mult4->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
-  if (chi2)
-  {
-    mult4->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
-    mult4->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
-  }
-  else
-    mult4->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
-  syst[1] = (TH1*) mult4->GetMultiplicityESDCorrected(etaRange)->Clone("result3");
-
-  DrawRatio(result1, 2, syst, "SystematicpTCutOff.eps", kFALSE, 0, kTRUE);
-
-  Draw2ResultRatio(mcHist, result1, syst[0], "SystematicpTCutOff1.eps");
-  Draw2ResultRatio(mcHist, result1, syst[1], "SystematicpTCutOff2.eps");
-}
-
-TH1* SystematicsSummary(Bool_t tpc = 1)
-{
-  Int_t nEffects = 7;
-
-  TH1* effects[10];
-  const char** names = 0;
-  Int_t colors[] = { 1, 2, 3, 4, 6, 7, 8 };
-  Int_t markers[] = { 20, 21, 22, 23, 24, 25, 26 };
-
-  for (Int_t i=0; i<nEffects; ++i)
-    effects[i] = new TH1F("SystematicsSummary", ";true multiplicity;Effect", 201, -0.5, 200.5);
-
-  if (tpc)
-  {
-    SetTPC();
-
-    const char* namesTPC[] = { "Unfolding method (#chi^{2})", "Rel. cross-section", "Particle composition", "p_{t} cut off", "Track selection", "Secondaries", "p_{t} spectrum" };
-    names = namesTPC;
-
-    // method
-    TFile* file = TFile::Open("StatisticalUncertaintyTPCChi2.root");
-    TH1* hist = (TH1*) file->Get("errorBoth");
-
-    // smooth a bit, but skip 0 bin
-    effects[0]->SetBinContent(2, hist->GetBinContent(2));
-    for (Int_t i=3; i<=200; ++i)
-      effects[0]->SetBinContent(i, (hist->GetBinContent(i) + hist->GetBinContent(i+1)) / 2);
-
-    // relative x-section
-    effects[1]->SetBinContent(2, 0.005);
-    effects[1]->SetBinContent(3, 0.0025);
-    effects[1]->SetBinContent(4, 0.0025);
-
-    // particle composition
-    for (Int_t i=2; i<=101; ++i)
-    {
-      if (i < 41)
-      {
-        effects[2]->SetBinContent(i, 0.01);
-      }
-      else if (i < 76)
-      {
-        effects[2]->SetBinContent(i, 0.02);
-      }
-      else
-        effects[2]->SetBinContent(i, 0.02 + 0.08 / 25 * (i - 76));
-    }
-
-    // pt cut off (only tpc)
-    for (Int_t i=2; i<=101; ++i)
-    {
-      if (i < 11)
-      {
-        effects[3]->SetBinContent(i, 0.05);
-      }
-      else if (i < 51)
-      {
-        effects[3]->SetBinContent(i, 0.01);
-      }
-      else
-        effects[3]->SetBinContent(i, 0.01 + 0.1 / 30 * (i - 51));
-    }
-
-    // track selection (only tpc)
-    for (Int_t i=2; i<=101; ++i)
-      effects[4]->SetBinContent(i, 0.03);
-
-    // secondaries
-    for (Int_t i=2; i<=101; ++i)
-      effects[5]->SetBinContent(i, 0.01);
-
-    // pt spectrum
-    for (Int_t i=2; i<=101; ++i)
-    {
-      if (i < 21)
-      {
-        effects[6]->SetBinContent(i, 0.05);
-      }
-      else if (i < 51)
-      {
-        effects[6]->SetBinContent(i, 0.02);
-      }
-      else
-        effects[6]->SetBinContent(i, 0.02 + 0.13 / 25 * (i - 51));
-    }
-
-  }
-  else
-  {
-    displayRange = 200;
-    nEffects = 5;
-
-    const char* namesSPD[] = { "Unfolding Method (#chi^{2})", "Rel. cross-section", "Particle composition", "Secondaries", "p_{t} spectrum"};
-    names = namesSPD;
-
-    // method
-    TFile* file = TFile::Open("StatisticalUncertaintySPDChi2.root");
-    TH1* hist = (TH1*) file->Get("errorBoth");
-
-    // smooth a bit, but skip 0 bin
-    effects[0]->SetBinContent(2, hist->GetBinContent(2));
-    for (Int_t i=3; i<=201; ++i)
-      effects[0]->SetBinContent(i, (hist->GetBinContent(i) + hist->GetBinContent(i+1)) / 2);
-
-    // relative x-section
-    effects[1]->SetBinContent(2, 0.01);
-    effects[1]->SetBinContent(3, 0.005);
-
-    // particle composition
-    for (Int_t i=2; i<=201; ++i)
-    {
-      if (i < 6)
-      {
-        effects[2]->SetBinContent(i, 0.3);
-      }
-      else if (i < 11)
-      {
-        effects[2]->SetBinContent(i, 0.05);
-      }
-      else if (i < 121)
-      {
-        effects[2]->SetBinContent(i, 0.02);
-      }
-      else if (i < 151)
-      {
-        effects[2]->SetBinContent(i, 0.02 + 0.04 / 30 * (i - 121));
-      }
-      else
-        effects[2]->SetBinContent(i, 0.06 + 0.1 / 30 * (i - 151));
-    }
-
-    // secondaries
-    for (Int_t i=2; i<=201; ++i)
-      effects[3]->SetBinContent(i, 0.01);
-
-    // pt spectrum
-    for (Int_t i=2; i<=201; ++i)
-    {
-      if (i < 6)
-      {
-        effects[4]->SetBinContent(i, 1);
-      }
-      else if (i < 121)
-      {
-        effects[4]->SetBinContent(i, 0.03);
-      }
-      else if (i < 151)
-      {
-        effects[4]->SetBinContent(i, 0.03 + 0.07 / 30 * (i - 121));
-      }
-      else
-        effects[4]->SetBinContent(i, 0.1);
-    }
-  }
-
-  TCanvas* canvas = new TCanvas("SystematicsSummary.eps", "SystematicsSummary.eps", 800, 400);
-  canvas->SetRightMargin(0.25);
-  canvas->SetTopMargin(0.05);
-  TLegend* legend = new TLegend(0.2, 0.4, 0.5, 0.4 + 0.5 * nEffects / 7);
-  legend->SetFillColor(0);
-
-  for (Int_t i=0; i<nEffects; ++i)
-  {
-    TH1* current = (TH1*) effects[i]->Clone(Form("current_%d", i));
-    /*current->Reset();
-    for (Int_t j=0; j<nEffects-i; ++j)
-      current->Add(effects[j]);*/
-
-    current->SetLineColor(colors[i]);
-    //current->SetFillColor(colors[i]);
-    current->SetMarkerColor(colors[i]);
-    //current->SetMarkerStyle(markers[i]);
-
-    current->SetStats(kFALSE);
-    current->GetYaxis()->SetRangeUser(0, 0.4);
-    current->GetXaxis()->SetRangeUser(0, displayRange);
-    current->DrawCopy(((i == 0) ? "" : "SAME"));
-    legend->AddEntry(current, names[i]);
-
-    TLatex* text = new TLatex(displayRange+5, current->GetBinContent(displayRange+1), names[i]);
-    text->SetTextColor(colors[i]);
-    text->Draw();
-  }
-
-  // add total in square
-  TH1* total = (TH1*) effects[0]->Clone("total");
-  total->Reset();
-
-  for (Int_t i=0; i<nEffects; ++i)
-  {
-    //Printf("%d %f", i, effects[i]->GetBinContent(20));
-    effects[i]->Multiply(effects[i]);
-    total->Add(effects[i]);
-  }
-
-  for (Int_t i=1; i<=total->GetNbinsX(); ++i)
-    if (total->GetBinContent(i) > 0)
-      total->SetBinContent(i, TMath::Min(sqrt(total->GetBinContent(i)), 1.0));
-
-  //Printf("%f", total->GetBinContent(20));
-
-  total->SetMarkerStyle(3);
-  total->SetMarkerColor(1);
-  legend->AddEntry(total, "total");
-  total->DrawCopy("SAME P");
-
-  legend->Draw();
-
-  canvas->SaveAs(canvas->GetName());
-
-  return total;
-}
-
-void finalPlot(Bool_t tpc = kTRUE, Bool_t chi2 = kTRUE, Bool_t small = kFALSE)
-{
-  gSystem->Load("libPWG0base");
-
-  if (tpc)
-    SetTPC();
-
-  if (!chi2)
-    Printf("WARNING: Bayesian set. This is only for test!");
-
-  // systematic error
-  TH1* error = SystematicsSummary(tpc);
-
-  TFile::Open(correctionFile);
-  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
-  mult->LoadHistograms("Multiplicity");
-
-  TFile::Open(measuredFile);
-  AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
-  mult2->LoadHistograms("Multiplicity");
-
-  mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
-
-  if (chi2)
-  {
-    mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
-    mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kINEL);
-  }
-  else
-    mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kINEL, 1, 100, 0, kFALSE);
-
-  TH1* mcHist = mult2->GetMultiplicityINEL(etaRange)->ProjectionY("mymc");
-  TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
-
-  DrawResultRatio(mcHist, result, "finalPlotCheck.eps");
-
-  // normalize result
-  result->Scale(1.0 / result->Integral(2, 200));
-
-  result->GetXaxis()->SetRangeUser(0, ((tpc) ? displayRange : 200));
-  result->SetBinContent(1, 0); result->SetBinError(1, 0);
-  result->SetTitle(";true multiplicity;Probability");
-  result->SetLineColor(1);
-  result->SetStats(kFALSE);
-
-  TH1* systError = (TH1*) result->Clone("systError");
-  for (Int_t i=2; i<=systError->GetNbinsX(); ++i)
-    systError->SetBinError(i, systError->GetBinContent(i) * error->GetBinContent(i));
-
-  // change error drawing style
-  systError->SetFillColor(15);
-
-  TCanvas* canvas = new TCanvas("finalPlot.eps", "finalPlot.eps", (small) ? 600 : 800, 400);
-  canvas->SetRightMargin(0.05);
-  canvas->SetTopMargin(0.05);
-
-  systError->Draw("E2 ][");
-  result->DrawCopy("SAME E ][");
-  canvas->SetLogy();
-
-  //TPaveText* text = new TPaveText(10, 1e-3, 50, 1e-4, "B");
-  TPaveText* text = new TPaveText(0.15, 0.2, 0.5, 0.4, "B NDC");
-  text->SetFillColor(0);
-  text->SetTextAlign(12);
-  text->AddText("Systematic errors summed quadratically");
-  text->AddText("0.6 million minimum bias events");
-  text->AddText("corrected to inelastic events");
-  text->Draw("B");
-
-  TPaveText* text2 = new TPaveText(0.4, 0.7, 0.6, 0.85, "B NDC");
-  text2->SetFillColor(0);
-  text2->SetTextAlign(12);
-  text2->AddText("#sqrt{s} = 14 TeV");
-  if (tpc)
-  {
-    text2->AddText("|#eta| < 0.9");
-  }
-  else
-    text2->AddText("|#eta| < 2.0");
-  text2->AddText("simulated data (PYTHIA)");
-  text2->Draw("B");
-
-  if (tpc)
-  {
-    TText* text3 = new TText(0.75, 0.6, "TPC - full tracking");
-    text3->SetNDC();
-    text3->Draw();
-  }
-  else
-  {
-    TText* text3 = new TText(0.75, 0.6, "SPD - Tracklets");
-    text3->SetNDC();
-    text3->Draw();
-  }
-
-  // alice logo
-  TPad* pad = new TPad("pad", "pad", 0.8, 0.7, 0.9, 0.9);
-  pad->Draw();
-  pad->cd();
-  TImage* img = TImage::Open("$HOME/alice.png");
-  img->SetImageQuality(TAttImage::kImgBest);
-  img->Draw();
-
-  canvas->Modified();
-
-/*  TText* text = new TText(10, 1e-4, "Systematic errors summed quadratically");
-  text->SetTextSize(0.04);
-  text->DrawText(10, 5e-5, "0.6 #cdot 10^{6} minimum bias events");
-  text->DrawText(10, 3e-5, "TPC tracks in |#eta| < 0.9");
-  text->DrawText(10, 1e-5, "corrected to ineleastic events in |#eta| < 0.9");
-  text->Draw();*/
-
-
-  canvas->SaveAs(canvas->GetName());
-}
-
-void BlobelUnfoldingExample()
-{
-  const Int_t kSize = 20;
-
-  TMatrixD matrix(kSize, kSize);
-  for (Int_t x=0; x<kSize; x++)
-  {
-    for (Int_t y=0; y<kSize; y++)
-    {
-      if (x == y)
-      {
-        if (x == 0 || x == kSize -1)
-        {
-          matrix(x, y) = 0.75;
-        }
-        else
-          matrix(x, y) = 0.5;
-      }
-      else if (TMath::Abs(x - y) == 1)
-      {
-        matrix(x, y) = 0.25;
-      }
-    }
-  }
-
-  //matrix.Print();
-
-  TMatrixD inverted(matrix);
-  inverted.Invert();
-
-  //inverted.Print();
-
-  TH1F* inputDist = new TH1F("inputDist", ";t;#tilde{T}(t)", kSize, -0.5, (Float_t) kSize - 0.5);
-  TVectorD inputDistVector(kSize);
-  TH1F* unfolded = inputDist->Clone("unfolded");
-  TH1F* measuredIdealDist = inputDist->Clone("measuredIdealDist");
-  measuredIdealDist->SetTitle(";m;#tilde{M}(m)");
-  TH1F* measuredDist = measuredIdealDist->Clone("measuredDist");
-
-  TF1* gaus = new TF1("func", "gaus(0)", -0.5, kSize);
-  // norm: 1/(sqrt(2pi)sigma)
-  gaus->SetParameters(10000 / sqrt(2 * TMath::Pi()) / ((Float_t) kSize / 8), (Float_t) kSize / 2, (Float_t) kSize / 8);
-  //gaus->Print();
-
-  for (Int_t x=1; x<=inputDist->GetNbinsX(); x++)
-  {
-    Float_t value = gaus->Eval(inputDist->GetBinCenter(x));
-    inputDist->SetBinContent(x, value);
-    inputDistVector(x-1) = value;
-  }
-
-  TVectorD measuredDistIdealVector = matrix * inputDistVector;
-  
-  for (Int_t x=1; x<=measuredIdealDist->GetNbinsX(); x++)
-    measuredIdealDist->SetBinContent(x, measuredDistIdealVector(x-1));
-
-  measuredDist->FillRandom(measuredIdealDist, 10000);
-
-  // fill error matrix before scaling
-  TMatrixD covarianceMatrixMeasured(kSize, kSize);
-  for (Int_t x=1; x<=unfolded->GetNbinsX(); x++)
-    covarianceMatrixMeasured(x-1, x-1) = TMath::Sqrt(measuredDist->GetBinContent(x));
-
-  TMatrixD covarianceMatrix = inverted * covarianceMatrixMeasured * inverted;
-  //covarianceMatrix.Print();
-
-  TVectorD measuredDistVector(kSize);
-  for (Int_t x=1; x<=measuredDist->GetNbinsX(); x++)
-    measuredDistVector(x-1) = measuredDist->GetBinContent(x);
-
-  TVectorD unfoldedVector = inverted * measuredDistVector;
-  for (Int_t x=1; x<=unfolded->GetNbinsX(); x++)
-    unfolded->SetBinContent(x, unfoldedVector(x-1));
-
-  TCanvas* canvas = new TCanvas("BlobelUnfoldingExample", "BlobelUnfoldingExample", 1000, 500);
-  canvas->SetTopMargin(0.05);
-  canvas->Divide(2, 1);
-
-  canvas->cd(1);
-  canvas->cd(1)->SetLeftMargin(0.15);
-  canvas->cd(1)->SetRightMargin(0.05);
-  measuredDist->GetYaxis()->SetTitleOffset(1.7);
-  measuredDist->SetStats(0);
-  measuredDist->DrawCopy();
-  gaus->Draw("SAME");
-
-  canvas->cd(2);
-  canvas->cd(2)->SetLeftMargin(0.15);
-  canvas->cd(2)->SetRightMargin(0.05);
-  unfolded->GetYaxis()->SetTitleOffset(1.7);
-  unfolded->SetStats(0);
-  unfolded->DrawCopy();
-  gaus->Draw("SAME");
-
-  canvas->SaveAs("BlobelUnfoldingExample.eps");
-}
-
-void E735Fit()
-{
-  TH1* fCurrentESD = new TH1F("mult", "mult", 501, -0.5, 500.5);
-  fCurrentESD->Sumw2();
-
-  // Open the input stream
-  ifstream in;
-  in.open("e735data.txt");
-
-  while(in.good())
-  {
-    Float_t x, y, ye;
-    in >> x >> y >> ye;
-
-    //Printf("%f %f %f", x, y, ye);
-    fCurrentESD->SetBinContent(fCurrentESD->FindBin(x), y);
-    fCurrentESD->SetBinError(fCurrentESD->FindBin(x), ye);
-  }
-
-  in.close();
-
-  //new TCanvas; fCurrentESD->DrawCopy(); gPad->SetLogy();
-
-  fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
-
-  TF1* func = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])");
-  func->SetParNames("scaling", "averagen", "k");
-  func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
-  func->SetParLimits(1, 0.001, 1000);
-  func->SetParLimits(2, 0.001, 1000);
-  func->SetParameters(fCurrentESD->GetMaximum() * 100, 10, 2);
-
-  TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500);
-  lognormal->SetParNames("scaling", "mean", "sigma");
-  lognormal->SetParameters(1, 1, 1);
-  lognormal->SetParLimits(0, 0, 10);
-  lognormal->SetParLimits(1, 0, 100);
-  lognormal->SetParLimits(2, 1e-3, 10);
-
-  TCanvas* canvas = new TCanvas("c1", "c1", 700, 400);
-  fCurrentESD->SetStats(kFALSE);
-  fCurrentESD->GetYaxis()->SetTitleOffset(1.3);
-  fCurrentESD->SetTitle(";true multiplicity (N);P_{N}");
-  fCurrentESD->Draw("");
-  fCurrentESD->GetXaxis()->SetRangeUser(0, 250);
-  fCurrentESD->Fit(func, "0", "", 0, 150);
-  func->SetRange(0, 250);
-  func->Draw("SAME");
-  printf("chi2 = %f\n", func->GetChisquare());
-
-  fCurrentESD->Fit(lognormal, "0", "", 0.01, 150);
-  lognormal->SetLineColor(2);
-  lognormal->SetLineStyle(2);
-  lognormal->SetRange(0, 250);
-  lognormal->Draw("SAME");
-
-  gPad->SetLogy();
-
-  canvas->SaveAs("E735Fit.eps");
-}
index 7676fe7..1aeca05 100644 (file)
@@ -53,8 +53,8 @@ void run(Int_t runWhat, const Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool
     }
     else
     {
-      gProof->UploadPackage("$ALICE_ROOT/AF-v4-12");
-      gProof->EnablePackage("$ALICE_ROOT/AF-v4-12");
+      gProof->UploadPackage("/afs/cern.ch/alice/caf/sw/ALICE/PARs/v4-14-Release/AF-v4-14");
+      gProof->EnablePackage("/afs/cern.ch/alice/caf/sw/ALICE/PARs/v4-14-Release/AF-v4-14");
     }
 
     gProof->UploadPackage("$ALICE_ROOT/PWG0base");
@@ -66,6 +66,7 @@ void run(Int_t runWhat, const Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool
     gSystem->Load("libTree");
     gSystem->Load("libSTEERBase");
     gSystem->Load("libESD");
+    gSystem->Load("libAOD");
     gSystem->Load("libANALYSIS");
     gSystem->Load("libANALYSISalice");
     gSystem->Load("libPWG0base");
@@ -101,6 +102,7 @@ void run(Int_t runWhat, const Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool
       printf("ERROR: esdTrackCuts could not be created\n");
       return;
     }
+    esdTrackCuts->SetHistogramsOn(kTRUE);
   }
 
   if (runWhat == 0)
@@ -170,6 +172,7 @@ void run(Int_t runWhat, const Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool
     // Create chain of input files
     gROOT->LoadMacro("../CreateESDChain.C");
     chain = CreateESDChain(data, nRuns, offset);
+    //chain = CreateChain("TE", data, nRuns, offset);
 
     mgr->StartAnalysis((aProof > 0) ? "proof" : "local", chain);
   }
@@ -236,7 +239,15 @@ void FinishAnalysisAll(const char* dataInput = "analysis_esd_raw.root", const ch
   file->cd();
   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks");
   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
-  fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone, "ESD raw");
+  fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone, "ESD raw with pt cut");
+  //fdNdEtaAnalysis->DrawHistograms(kTRUE);
+  file2->cd();
+  fdNdEtaAnalysis->SaveHistograms();
+
+  file->cd();
+  fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTracksAll", "dndetaTracksAll");
+  fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
+  fdNdEtaAnalysis->Finish(0, -1, AlidNdEtaCorrection::kNone, "ESD raw");
   //fdNdEtaAnalysis->DrawHistograms(kTRUE);
   file2->cd();
   fdNdEtaAnalysis->SaveHistograms();
diff --git a/PWG0/dNdEta/runVertexSelector.C b/PWG0/dNdEta/runVertexSelector.C
deleted file mode 100644 (file)
index 1db53d6..0000000
+++ /dev/null
@@ -1,40 +0,0 @@
-/* $Id$ */
-
-//
-// script to run the AliVertexSelector
-//
-
-#include "../CreateESDChain.C"
-#include "../PWG0Helper.C"
-
-void runVertexSelector(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aProof = kFALSE)
-{
-  if (aProof)
-    connectProof("jgrosseo@lxb6046");
-
-  TString libraries("libEG;libGeom;libESD;libPWG0base");
-  TString packages("PWG0base");
-
-  //libraries += ";libVMC;libMinuit;libSTEER;libPWG0dep;libEVGEN;libFASTSIM;libmicrocern;libpdf;libpythia6;libEGPythia6;libAliPythia6";
-
-  if (!prepareQuery(libraries, packages, 1))
-    return;
-
-  TChain* chain = CreateESDChain(data, nRuns, offset);
-
-  TList inputList;
-
-  gROOT->ProcessLine(".L CreateCuts.C");
-
-  AliESDtrackCuts* esdTrackCuts = CreateTrackCuts();
-  if (!esdTrackCuts)
-  {
-    printf("ERROR: esdTrackCuts could not be created\n");
-    return;
-  }
-
-  inputList.Add(esdTrackCuts);
-
-  executeQuery(chain, &inputList, "AliVertexSelector.cxx+");
-}
-
diff --git a/PWG0/dNdEta/rundNdEtaAnalysis.C b/PWG0/dNdEta/rundNdEtaAnalysis.C
deleted file mode 100644 (file)
index 7c97279..0000000
+++ /dev/null
@@ -1,72 +0,0 @@
-/* $Id$ */
-
-//
-// Script to test the dN/dEta analysis using the dNdEtaAnalysis and
-// dNdEtaCorrection classes. Note that there is a cut on the events,
-// so the measurement will be biassed.
-//
-// implementation with TSelector
-//
-
-#include "../CreateESDChain.C"
-#include "../PWG0Helper.C"
-
-void rundNdEtaAnalysis(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 = "lxb6046")
-{
-  if (aProof)
-    connectProof(proofServer);
-
-  TString libraries("libEG;libGeom;libESD;libPWG0base");
-  TString packages("PWG0base");
-
-  if (!prepareQuery(libraries, packages, kTRUE))
-    return;
-
-  gROOT->ProcessLine(".L CreateCuts.C");
-  gROOT->ProcessLine(".L drawPlots.C");
-
-  TChain* chain = CreateESDChain(data, nRuns, offset);
-
-  TList inputList;
-
-  // selection of esd tracks
-  AliESDtrackCuts* esdTrackCuts = CreateTrackCuts();
-  if (!esdTrackCuts)
-  {
-    printf("ERROR: esdTrackCuts could not be created\n");
-    return;
-  }
-
-  inputList.Add(esdTrackCuts);
-
-  TString selectorName = ((aMC == kFALSE) ? "AlidNdEtaAnalysisESDSelector" : "AlidNdEtaAnalysisMCSelector");
-  AliLog::SetClassDebugLevel(selectorName, AliLog::kInfo);
-
-  selectorName += ".cxx+";
-
-  if (aDebug != kFALSE)
-    selectorName += "+g";
-
-  Int_t result = executeQuery(chain, &inputList, selectorName, option);
-
-  if (result >= 0)
-  {
-    if (aMC)
-    {
-      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(kTRUE);
-    }
-    else
-      FinishAnalysisAll("analysis_esd_raw.root", "analysis_esd.root", correctionMapFile, correctionMapFolder);
-  }
-}
-