Classes for track-based multiplicity analysis in PbPb. Initial commit.
authormfloris <mfloris@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 29 Oct 2010 19:35:39 +0000 (19:35 +0000)
committermfloris <mfloris@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 29 Oct 2010 19:35:39 +0000 (19:35 +0000)
PWG0/multPbPb/AliAnalysisMultPbCentralitySelector.cxx [new file with mode: 0644]
PWG0/multPbPb/AliAnalysisMultPbCentralitySelector.h [new file with mode: 0644]
PWG0/multPbPb/AliAnalysisMultPbTrackHistoManager.cxx [new file with mode: 0644]
PWG0/multPbPb/AliAnalysisMultPbTrackHistoManager.h [new file with mode: 0644]
PWG0/multPbPb/AliAnalysisTaskMultPbTracks.cxx [new file with mode: 0644]
PWG0/multPbPb/AliAnalysisTaskMultPbTracks.h [new file with mode: 0644]
PWG0/multPbPb/correct.C [new file with mode: 0644]
PWG0/multPbPb/run.C [new file with mode: 0644]
PWG0/multPbPb/run.sh [new file with mode: 0755]

diff --git a/PWG0/multPbPb/AliAnalysisMultPbCentralitySelector.cxx b/PWG0/multPbPb/AliAnalysisMultPbCentralitySelector.cxx
new file mode 100644 (file)
index 0000000..9a3bd65
--- /dev/null
@@ -0,0 +1,4 @@
+#include "AliAnalysisMultPbCentralitySelector.h"
+
+ClassImp(AliAnalysisMultPbCentralitySelector)
+
diff --git a/PWG0/multPbPb/AliAnalysisMultPbCentralitySelector.h b/PWG0/multPbPb/AliAnalysisMultPbCentralitySelector.h
new file mode 100644 (file)
index 0000000..eb07d7c
--- /dev/null
@@ -0,0 +1,58 @@
+#ifndef ALIANALYSISMULTPBCENTRALITYSELECTOR_H
+#define ALIANALYSISMULTPBCENTRALITYSELECTOR_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//-------------------------------------------------------------------------
+//                    AliAnalysisMultPbCentralitySelector
+// 
+// This class selects collision candidates from data runs, applying selection cuts on triggers 
+// and background rejection based on the content of the ESD
+//
+// Author: Michele Floris, CERN
+//-------------------------------------------------------------------------
+
+#include <AliAnalysisCuts.h>
+#include <AliLog.h>
+
+#define VERBOSE_STAT
+
+class AliESDEvent;
+class TH2F;
+class TH1F;
+class TCollection;
+class AliTriggerAnalysis;
+class AliAnalysisTaskSE;
+
+
+class AliAnalysisMultPbCentralitySelector : public AliAnalysisCuts
+{
+public:
+
+  AliAnalysisMultPbCentralitySelector() : fIsMC (0) {;}
+  virtual ~AliAnalysisMultPbCentralitySelector(){}
+    
+  // AliAnalysisCuts interface
+  virtual UInt_t GetSelectionMask(const TObject* obj) { return IsCentralityBinSelected((const AliESDEvent*) obj); }
+  virtual Bool_t IsSelected(TList*) { AliFatal("Not implemented"); return kFALSE; }
+  virtual Bool_t IsSelected(TObject* obj)  {return IsCentralityBinSelected ( (AliESDEvent*) obj);}
+    
+  UInt_t IsCentralityBinSelected(const AliESDEvent* aEsd){ return kTRUE;}
+    
+  void SetAnalyzeMC(Bool_t flag = kTRUE) { fIsMC = flag; }
+    
+  virtual void Print(Option_t* option = "") const { Printf ("Multiplitity Selector [AliAnalysisMultPbCentralitySelector] [%s]", option);}
+  virtual Long64_t Merge(TCollection* list){list->GetEntries();return 0;}
+  
+protected:
+  Bool_t fIsMC;             // flag if MC is analyzed
+
+  ClassDef(AliAnalysisMultPbCentralitySelector, 1)
+    
+  private:
+  AliAnalysisMultPbCentralitySelector(const AliAnalysisMultPbCentralitySelector&); // not implemented
+  AliAnalysisMultPbCentralitySelector& operator=(const AliAnalysisMultPbCentralitySelector&); // not implemented
+};
+
+#endif
diff --git a/PWG0/multPbPb/AliAnalysisMultPbTrackHistoManager.cxx b/PWG0/multPbPb/AliAnalysisMultPbTrackHistoManager.cxx
new file mode 100644 (file)
index 0000000..a8b7d8b
--- /dev/null
@@ -0,0 +1,272 @@
+#include "AliAnalysisMultPbTrackHistoManager.h"
+#include "AliLog.h"
+#include "TH1.h"
+#include "TH3D.h"
+#include "TH1I.h"
+#include "TROOT.h"
+#include <iostream>
+using namespace std;
+ClassImp(AliAnalysisMultPbTrackHistoManager)
+
+const char * AliAnalysisMultPbTrackHistoManager::kStatStepNames[]     = { "All Events", "After centrality selection", "With Vertex" };
+const char * AliAnalysisMultPbTrackHistoManager::kHistoPtEtaVzNames[] = { "hGenPtEtaVz", "hRecPtEtaVz", "hRecPtEtaVzPrim", 
+                                                                         "hRecPtEtaVzSecWeak", "hRecPtEtaVzSecMaterial", "hRecPtEtaVzFake"};
+const char * AliAnalysisMultPbTrackHistoManager::kHistoDCANames[]     = { "hGenDCA", "hRecDCA", "hRecDCAPrim", "hRecDCASecWeak","hRecDCASecMaterial", "hRecDCAFake"};
+
+
+
+AliAnalysisMultPbTrackHistoManager::AliAnalysisMultPbTrackHistoManager() : AliHistoListWrapper(), fHNameSuffix(""){ 
+  // standard ctor
+
+}
+
+AliAnalysisMultPbTrackHistoManager::AliAnalysisMultPbTrackHistoManager(const char * name, const char * title): AliHistoListWrapper(name,title), fHNameSuffix("")  {
+  //named ctor
+};
+
+AliAnalysisMultPbTrackHistoManager::AliAnalysisMultPbTrackHistoManager(const AliAnalysisMultPbTrackHistoManager& obj) : AliHistoListWrapper (obj) {
+  // copy ctor
+  AliError("Not Implemented");
+};
+
+AliAnalysisMultPbTrackHistoManager::~AliAnalysisMultPbTrackHistoManager() {
+  // dtor
+
+}
+
+TH3D * AliAnalysisMultPbTrackHistoManager::GetHistoPtEtaVz(Histo_t id) {
+  // Returns a 3D histo of Pt/eta/vtx. It it does not exist, books it.
+
+  TH3D * h = (TH3D*) GetHisto(kHistoPtEtaVzNames[id]);
+  if (!h) {
+    h = BookHistoPtEtaVz(kHistoPtEtaVzNames[id], Form("Pt Eta Vz distribution (%s)",kHistoPtEtaVzNames[id]));
+  }
+
+  return h;
+
+}
+
+TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoDCA(Histo_t id) {
+  // Returns a 3D histo of Pt/eta/vtx. It it does not exist, books it.
+
+  TH1D * h = (TH1D*) GetHisto(kHistoDCANames[id]);
+  if (!h) {
+    h = BookHistoDCA(kHistoDCANames[id], Form("Pt Eta Vz distribution (%s)",kHistoDCANames[id]));
+  }
+
+  return h;
+
+}
+
+TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoPt (Histo_t id, 
+                                                      Float_t minEta, Float_t maxEta, 
+                                                      Float_t minVz, Float_t maxVz, 
+                                                      Bool_t scaleWidth) {
+  // Returns a projection of the 3D histo pt/eta/vz.
+  // WARNING: since that is a histo, the requested range will be discretized to the binning.
+  // Always avoids under (over) flows
+  // If scaleWidth = kTRUE, the projection is scaled for the bin width (default)
+
+  TH3D * h3D = GetHistoPtEtaVz(id);
+
+  // Get range in terms of bin numners.  If the float range is
+  // less than -11111 take the range from the first to the last bin (i.e. no
+  // under/over-flows)
+  Int_t min1 = minEta  < -11111 ? 1 : h3D ->GetYaxis()->FindBin(minEta);
+  Int_t min2  = minVz  < -11111 ? 1 : h3D ->GetZaxis()->FindBin(minVz) ;
+
+  Int_t max1 = maxEta  < -11111 ? h3D->GetNbinsY() : h3D ->GetYaxis()->FindBin(maxEta-0.00001);
+  Int_t max2  = maxVz  < -11111 ? h3D->GetNbinsZ() : h3D ->GetZaxis()->FindBin(maxVz -0.00001);
+
+
+  TString hname = h3D->GetName();
+  hname = hname +  "_pt_" + long (min1)  + "_" + long(max1) + "_" + long (min2)  + "_" + long(max2);
+
+  
+  if (gROOT->FindObjectAny(hname.Data())){
+    AliError(Form("An object called %s already exists",hname.Data()));
+  }
+
+  TH1D * h = h3D->ProjectionX(hname.Data(), min1, max1, min2, max2, "E");
+  if(scaleWidth) h ->Scale(1.,"width");
+
+  return h;
+
+}
+
+TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoVz (Histo_t id, 
+                                                      Float_t minPt, Float_t maxPt,
+                                                      Float_t minEta, Float_t maxEta,
+                                                      Bool_t scaleWidth) { 
+  // Returns a projection of the 3D histo pt/eta/vz.
+  // WARNING: since that is a histo, the requested range will be discretized to the binning.
+  // Always avoids under (over) flows
+  // If scaleWidth = kTRUE, the projection is scaled for the bin width (default)
+
+  TH3D * h3D = GetHistoPtEtaVz(id);
+
+  // Get range in terms of bin numners.  If the float range is
+  // less than -11111 take the range from the first to the last bin (i.e. no
+  // under/over-flows)
+  Int_t min1  = minPt  < -11111 ? 1 : h3D ->GetXaxis()->FindBin(minPt) ;
+  Int_t min2  = minEta < -11111 ? 1 : h3D ->GetYaxis()->FindBin(minEta);
+
+  Int_t max1  = maxPt  < -11111 ? h3D->GetNbinsX() : h3D ->GetXaxis()->FindBin(maxPt -0.00001);
+  Int_t max2  = maxEta < -11111 ? h3D->GetNbinsY() : h3D ->GetYaxis()->FindBin(maxEta-0.00001);
+
+
+  TString hname = h3D->GetName();
+  hname = hname +  "_Vz_" + long (min1)  + "_" + long(max1) + "_" + long (min2)  + "_" + long(max2);
+
+  if (gROOT->FindObjectAny(hname.Data())){
+    AliError(Form("An object called %s already exists",hname.Data()));
+  }
+
+  TH1D * h = h3D->ProjectionZ(hname.Data(), min1, max1, min2, max2, "E");
+  if(scaleWidth) h ->Scale(1.,"width");
+  return h;
+
+
+}
+
+TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoEta (Histo_t id, 
+                                                       Float_t minPt, Float_t maxPt, 
+                                                       Float_t minVz, Float_t maxVz,
+                                                       Bool_t scaleWidth) {
+  // Returns a projection of the 3D histo pt/eta/vz.
+  // WARNING: since that is a histo, the requested range will be discretized to the binning.
+  // Always avoids under (over) flows
+  // If scaleWidth = kTRUE, the projection is scaled for the bin width (default)
+
+  TH3D * h3D = GetHistoPtEtaVz(id);
+
+  // Get range in terms of bin numners.  If the float range is
+  // less than -11111 take the range from the first to the last bin (i.e. no
+  // under/over-flows)
+  Int_t min1 = minPt < -11111 ? 1 : h3D ->GetXaxis()->FindBin(minPt) ;
+  Int_t min2 = minVz < -11111 ? 1 : h3D ->GetYaxis()->FindBin(minVz);
+
+  Int_t max1 = maxPt < -11111 ? h3D->GetNbinsX() : h3D ->GetXaxis()->FindBin(maxPt -0.00001);
+  Int_t max2 = maxVz < -11111 ? h3D->GetNbinsY() : h3D ->GetYaxis()->FindBin(maxVz-0.00001);
+
+  TString hname = h3D->GetName();
+  hname = hname +  "_Eta_" + long (min1)  + "_" + long(max1) + "_" + long (min2)  + "_" + long(max2);
+
+  if (gROOT->FindObjectAny(hname.Data())){
+    AliError(Form("An object called %s already exists",hname.Data()));
+  }
+
+  TH1D * h = h3D->ProjectionY(hname.Data(), min1, max1, min2, max2, "E");
+  if(scaleWidth) h ->Scale(1.,"width");
+  return h;
+}
+
+
+TH1I * AliAnalysisMultPbTrackHistoManager::GetHistoStats() {
+  // Returns histogram with event statistiscs (processed events at each step)
+
+  TH1I * h =  (TH1I*) GetHisto("hStats");
+  if (!h) h = BookHistoStats();
+  return h;
+  
+} 
+
+
+
+TH1 * AliAnalysisMultPbTrackHistoManager::GetHisto(const char * name) {
+  //Search list for histo
+  // TODO: keep track of histo id rather than searching by name?
+  return (TH1*) fList->FindObject(TString(name)+fHNameSuffix);
+
+}
+
+void AliAnalysisMultPbTrackHistoManager::ScaleHistos(Double_t nev, Option_t * option) {
+  // Scales all histos in the list for nev
+  // option can be used to pass further options to TH1::Scale
+  TH1 * h = 0;
+  TIter iter = fList->MakeIterator();
+  while ((h = (TH1*) iter.Next())) {
+    if (!h->InheritsFrom("TH1")) {
+      AliFatal (Form("%s does not inherits from TH1, cannot scale",h->GetName()));
+    }
+    cout << "Scaling " << h->GetName() << " " << nev << endl;
+    
+    h->Scale(1./nev,option);
+  }
+
+}
+
+TH3D * AliAnalysisMultPbTrackHistoManager::BookHistoPtEtaVz(const char * name, const char * title) {
+  // Books a 3D histo of Pt/eta/vtx
+  // TODO: make the binning settable, variable binning?
+
+  Bool_t oldStatus = TH1::AddDirectoryStatus();
+  TH1::AddDirectory(kFALSE);
+
+  TString hname = name;
+  hname+=fHNameSuffix;
+
+  AliInfo(Form("Booking %s",hname.Data()));
+  
+
+  TH3D * h = new TH3D (hname,title, 
+                      50,0,10, // Pt
+                      20,-1, 1,   // Eta
+                      10,-10,10 // Vz
+                      );
+
+  h->SetYTitle("#eta");
+  h->SetXTitle("p_{T}");
+  h->SetZTitle("V_{z} (cm)");
+  h->Sumw2();
+  
+  fList->Add(h);
+
+  TH1::AddDirectory(oldStatus);
+  return h;
+}
+
+TH1D * AliAnalysisMultPbTrackHistoManager::BookHistoDCA(const char * name, const char * title) {
+  // Books a DCA histo 
+
+  Bool_t oldStatus = TH1::AddDirectoryStatus();
+  TH1::AddDirectory(kFALSE);
+
+  TString hname = name;
+  hname+=fHNameSuffix;
+
+  AliInfo(Form("Booking %s",hname.Data()));
+  
+
+  TH1D * h = new TH1D (hname,title, 100,0,50);
+
+  h->SetXTitle("#Delta DCA");
+  h->Sumw2();
+  
+  fList->Add(h);
+
+  TH1::AddDirectory(oldStatus);
+  return h;
+}
+
+TH1I * AliAnalysisMultPbTrackHistoManager::BookHistoStats() {
+  // Books histogram with event statistiscs (processed events at each step)
+
+  Bool_t oldStatus = TH1::AddDirectoryStatus();
+  TH1::AddDirectory(kFALSE);
+
+  AliInfo(Form("Booking Stat histo"));
+
+  TH1I * h = new TH1I (TString("hStats")+fHNameSuffix, "Number of processed events", kNStatBins, -0.5, kNStatBins-0.5);
+  for(Int_t istatbin = 0; istatbin < kNStatBins; istatbin++){
+    h->GetXaxis()->SetBinLabel(istatbin+1,kStatStepNames[istatbin]);
+  }
+  TH1::AddDirectory(oldStatus);
+  fList->Add(h);
+  return h;
+}
+
+
+
+
diff --git a/PWG0/multPbPb/AliAnalysisMultPbTrackHistoManager.h b/PWG0/multPbPb/AliAnalysisMultPbTrackHistoManager.h
new file mode 100644 (file)
index 0000000..4e74ce6
--- /dev/null
@@ -0,0 +1,73 @@
+#ifndef ALIANALYSISMULTPBTRACKHISTOMANAGER_H
+#define ALIANALYSISMULTPBTRACKHISTOMANAGER_H
+
+#include "AliHistoListWrapper.h"
+#include "TH3D.h"
+
+class TH3D;
+class TH1D;
+class TH1I;
+
+//-------------------------------------------------------------------------
+//                      AliAnalysisMultPbTrackHistoManager
+// 
+// 
+//
+//
+// Author: Michele Floris, CERN
+//-------------------------------------------------------------------------
+
+
+
+class AliAnalysisMultPbTrackHistoManager : public AliHistoListWrapper {
+
+public:
+
+  typedef enum {kHistoGen, kHistoRec, kHistoRecPrim, kHistoRecSecWeak, kHistoRecSecMat, kHistoRecFake, kNHistos} Histo_t;
+  typedef enum {kStatAll, kStatCentr, kStatVtx, kNStatBins} Stat_t;
+
+  AliAnalysisMultPbTrackHistoManager();
+  AliAnalysisMultPbTrackHistoManager(const char * name,const char * title);
+  AliAnalysisMultPbTrackHistoManager(const AliAnalysisMultPbTrackHistoManager& obj) ;
+  ~AliAnalysisMultPbTrackHistoManager();
+  
+  // Setters
+  void SetSuffix(const char * suffix) { fHNameSuffix = suffix;}
+
+  // Histo getters
+  TH3D * GetHistoPtEtaVz(Histo_t id);
+  TH1D * GetHistoPt (Histo_t id, Float_t minEta = -22222, Float_t maxEta = -22222, Float_t minVz  = -22222, Float_t maxVz  = -22222, Bool_t scaleWidth = kTRUE);
+  TH1D * GetHistoEta(Histo_t id, Float_t minPt  = -22222, Float_t maxPt  = -22222, Float_t minVz  = -22222, Float_t maxVz  = -22222, Bool_t scaleWidth = kTRUE);
+  TH1D * GetHistoVz (Histo_t id, Float_t minPt  = -22222, Float_t maxPt  = -22222, Float_t minEta = -22222, Float_t maxEta = -22222, Bool_t scaleWidth = kTRUE);
+
+  TH1I * GetHistoStats();
+  TH1D * GetHistoDCA(Histo_t id);
+
+  // Misch utils
+  void ScaleHistos (Double_t nev, Option_t * option="");
+  
+
+
+  // Histo bookers
+  TH3D * BookHistoPtEtaVz(const char * name, const char * title);
+  TH1D * BookHistoDCA(const char * name, const char * title);
+  TH1I * BookHistoStats();
+
+  // 
+  TH1 * GetHisto(const char * name);
+
+private:
+
+  static const char * kStatStepNames[];       // names of the step hist
+  static const char * kHistoPtEtaVzNames[];   // names of the 3D histograms pt/eta/vz
+  static const char * kHistoDCANames[];   // names of the DCA histograms 
+  TString fHNameSuffix; // Suffix added to all histo names. Useful if you have in the same session e.g. MC and data.
+
+  AliAnalysisMultPbTrackHistoManager& operator=(const AliAnalysisMultPbTrackHistoManager& task);
+  
+  ClassDef(AliAnalysisMultPbTrackHistoManager, 2)
+
+
+};
+
+#endif
diff --git a/PWG0/multPbPb/AliAnalysisTaskMultPbTracks.cxx b/PWG0/multPbPb/AliAnalysisTaskMultPbTracks.cxx
new file mode 100644 (file)
index 0000000..bd33b8e
--- /dev/null
@@ -0,0 +1,277 @@
+// AliAnalysisTaskMultPbTracks
+
+// Author: Michele Floris, CERN
+// TODO:
+// - Add chi2/cluster plot for primary, secondaries and fakes
+
+
+#include "AliAnalysisTaskMultPbTracks.h"
+#include "AliESDInputHandler.h"
+#include "AliAnalysisMultPbTrackHistoManager.h"
+#include "AliAnalysisManager.h"
+#include "AliESDtrackCuts.h"
+#include "AliAnalysisMultPbCentralitySelector.h"
+#include "AliMCEvent.h"
+#include "AliStack.h"
+#include "TH1I.h"
+#include "TH3D.h"
+#include "AliMCParticle.h"
+#include "AliGenEventHeader.h"
+
+#include <iostream>
+
+using namespace std;
+
+ClassImp(AliAnalysisTaskMultPbTracks)
+
+AliAnalysisTaskMultPbTracks::AliAnalysisTaskMultPbTracks()
+: AliAnalysisTaskSE("TaskMultPbTracks"),fESD(0),fHistoManager(0),fCentrSelector(0),fTrackCuts(0),fTrackCutsNoDCA(0),fIsMC(0)
+{
+  // constructor
+
+  DefineOutput(1, AliAnalysisMultPbTrackHistoManager::Class());
+  DefineOutput(2, AliAnalysisMultPbCentralitySelector::Class());
+  DefineOutput(3, AliESDtrackCuts::Class());
+  //  DefineOutput(2, TH1I::Class());
+
+  fHistoManager = new AliAnalysisMultPbTrackHistoManager("histoManager","Hitograms, Multiplicity, Track analysis");
+  if(fIsMC) fHistoManager->SetSuffix("MC");
+
+}
+AliAnalysisTaskMultPbTracks::AliAnalysisTaskMultPbTracks(const char * name)
+  : AliAnalysisTaskSE(name),fESD(0),fHistoManager(0),fCentrSelector(0),fTrackCuts(0),fTrackCutsNoDCA(0),fIsMC(0)
+{
+  //
+  // Standard constructur which should be used
+  //
+
+  DefineOutput(1, AliAnalysisMultPbTrackHistoManager::Class());
+  DefineOutput(2, AliAnalysisMultPbCentralitySelector::Class());
+  DefineOutput(3, AliESDtrackCuts::Class());
+
+  fHistoManager = new AliAnalysisMultPbTrackHistoManager("histoManager","Hitograms, Multiplicity, Track analysis");
+  if(fIsMC) fHistoManager->SetSuffix("MC");
+
+}
+
+AliAnalysisTaskMultPbTracks::AliAnalysisTaskMultPbTracks(const AliAnalysisTaskMultPbTracks& obj) : 
+  AliAnalysisTaskSE(obj) ,fESD (0), fHistoManager(0), fCentrSelector(0), fTrackCuts(0),fTrackCutsNoDCA(0),fIsMC(0)
+{
+  //copy ctor
+  fESD = obj.fESD ;
+  fHistoManager= obj.fHistoManager;
+  fCentrSelector = obj.fCentrSelector;
+  fTrackCuts  = obj.fTrackCuts;
+  fTrackCutsNoDCA  = obj.fTrackCutsNoDCA;
+
+}
+
+AliAnalysisTaskMultPbTracks::~AliAnalysisTaskMultPbTracks(){
+  // destructor
+
+  if(!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
+    if(fHistoManager) {
+      delete fHistoManager;
+      fHistoManager = 0;
+    }
+  }
+  // Histo list should not be destroyed: fListWrapper is owner!
+
+}
+void AliAnalysisTaskMultPbTracks::UserCreateOutputObjects()
+{
+  // Called once
+
+  // For the DCA distributions, we want to use exactly the same cuts
+  // as for the other distributions, with the exception of the DCA cut
+  fTrackCutsNoDCA = new AliESDtrackCuts(*fTrackCuts); // clone cuts
+  // Reset all DCA cuts; FIXME: is this all?
+  fTrackCutsNoDCA->SetMaxDCAToVertexXY();
+  fTrackCutsNoDCA->SetMaxDCAToVertexZ ();
+  fTrackCutsNoDCA->SetMaxDCAToVertexXYPtDep();
+  fTrackCutsNoDCA->SetMaxDCAToVertexZPtDep();
+
+}
+
+
+void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
+{
+  // User code
+
+  /* PostData(0) is taken care of by AliAnalysisTaskSE */
+  PostData(1,fHistoManager);
+  PostData(2,fCentrSelector);
+  PostData(3,fTrackCuts);
+
+  TH3D * hTracks[AliAnalysisMultPbTrackHistoManager::kNHistos];
+  TH1D * hDCA   [AliAnalysisMultPbTrackHistoManager::kNHistos];
+  hTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]        = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen       );
+  hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec]        = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec       );
+  hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]    = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim   );
+  hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]    = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecFake   );
+  hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]  = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat );
+  hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak] = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak);
+
+  hDCA[AliAnalysisMultPbTrackHistoManager::kHistoGen]        = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen       );
+  hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRec]        = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec       );
+  hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]    = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim   );
+  hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]    = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake   );
+  hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]  = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat );
+  hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak] = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak);
+
+
+  fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
+  if (strcmp(fESD->ClassName(),"AliESDEvent")) {
+    AliFatal("Not processing ESDs");
+  }
+  
+  // FIXME: use physics selection here to keep track of events lost?
+  fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatAll);
+
+  Bool_t isCentralitySelected = fCentrSelector->IsSelected(fESD);  
+  if(!isCentralitySelected) return;
+
+  fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatCentr);
+
+  if (fIsMC) {
+    
+
+    if (!fMCEvent) {
+      AliError("No MC info found");
+    } else {
+      
+      //loop on the MC event
+      Int_t nMCTracks = fMCEvent->GetNumberOfTracks();
+      for (Int_t ipart=0; ipart<nMCTracks; ipart++) { 
+       
+       AliMCParticle *mcPart  = (AliMCParticle*)fMCEvent->GetTrack(ipart);
+       
+       // We don't care about neutrals and non-physical primaries
+       if(mcPart->Charge() == 0) continue;
+
+       // FIXME: add kTransportBit
+       if(!fMCEvent->Stack()->IsPhysicalPrimary(ipart)) continue;
+       
+       // Get MC vertex
+       //FIXME: which vertex do I take for MC?
+       TArrayF   vertex;
+       fMCEvent->GenEventHeader()->PrimaryVertex(vertex);
+       Float_t zv = vertex[2];
+       //      Float_t zv = vtxESD->GetZ();
+       // Fill generated histo
+       hTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]->Fill(mcPart->Pt(),mcPart->Eta(),zv);
+       
+      }
+    }
+  }
+  
+  // FIXME: shall I take the primary vertex?
+  const AliESDVertex* vtxESD = fESD->GetPrimaryVertex();
+  if(!vtxESD) return;
+  fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatVtx);
+
+
+  // loop on tracks
+  Int_t ntracks = fESD->GetNumberOfTracks();
+
+
+  for (Int_t itrack = 0; itrack<ntracks; itrack++) {    
+    AliESDtrack * esdTrack = fESD->GetTrack(itrack);
+
+    // Fill DCA distibutions, without the DCA cut, Fill the other stuff, with all cuts!
+    Bool_t accepted = fTrackCuts->AcceptTrack(esdTrack);
+    Bool_t acceptedNoDCA = fTrackCutsNoDCA->AcceptTrack(esdTrack);
+
+    // Compute weighted offset
+    // TODO: other choiches for the DCA?
+    Double_t b = fESD->GetMagneticField();
+    Double_t dca[2];
+    Double_t cov[3];
+    Double_t weightedDCA = 10;
+    
+
+
+    if (esdTrack->PropagateToDCA(vtxESD, b,3., dca, cov)) {
+      Double_t det = cov[0]*cov[2]-cov[1]*cov[1]; 
+      if (det<=0) {
+       AliError("DCA Covariance matrix is not positive definite!");
+      }
+      else {
+       weightedDCA = (dca[0]*dca[0]*cov[2]+dca[1]*dca[1]*cov[0]-2*dca[0]*dca[1]*cov[1])/det; 
+       weightedDCA = weightedDCA>0 ? TMath::Sqrt(weightedDCA) : 10;
+      }
+      //      printf("dR:%e dZ%e  sigR2:%e sigRZ:%e sigZ2:%e\n",dca[0],dca[1],cov[0],cov[1],cov[2]);
+    }
+    
+
+    
+
+    // Alternative: p*DCA
+    // Float_t xz[2]; 
+    // esdTrack->GetDZ(vtxESD->GetX(),vtxESD->GetY(),vtxESD->GetZ(), fESD->GetMagneticField(), xz); 
+    // Float_t dca = xz[0]*esdTrack->P();
+
+    // for each track
+    if(accepted) 
+      hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
+    if(acceptedNoDCA)
+      hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRec]->Fill(weightedDCA);
+
+    // Fill separately primaries and secondaries
+    // FIXME: fakes? Is this the correct way to account for them?
+    // Get label and corresponding mcPart;
+    if (fIsMC) {
+      //      Int_t label = TMath::Abs(esdTrack->GetLabel()); // no fakes!!!
+      Int_t label = esdTrack->GetLabel(); // 
+      AliMCParticle *mcPart  = (AliMCParticle*)fMCEvent->GetTrack(label);
+      if (!mcPart)  {
+       if(accepted)
+         hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
+       if(acceptedNoDCA)
+         hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]->Fill(weightedDCA);
+      }
+      else {
+       if(fMCEvent->Stack()->IsPhysicalPrimary(label)) {
+         if(accepted)
+           hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
+         if(acceptedNoDCA)
+           hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]->Fill(weightedDCA);
+       } 
+       else {
+         Int_t mfl=0;
+         Int_t indexMoth=mcPart->Particle()->GetFirstMother();
+         if(indexMoth>=0){
+           TParticle* moth = fMCEvent->Stack()->Particle(indexMoth);
+           Float_t codemoth = TMath::Abs(moth->GetPdgCode());
+           mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));
+         }
+         if(mfl==3){ // strangeness
+           if(accepted)
+             hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
+           if(acceptedNoDCA)
+             hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak]->Fill(weightedDCA);      
+         }else{ // material
+           if(accepted)
+             hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
+           if(acceptedNoDCA)
+             hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]->Fill(weightedDCA);       
+
+         }
+
+
+       }
+      }
+    }
+
+
+  }
+
+
+}
+
+void   AliAnalysisTaskMultPbTracks::Terminate(Option_t *){
+  // terminate
+
+}
+
+
diff --git a/PWG0/multPbPb/AliAnalysisTaskMultPbTracks.h b/PWG0/multPbPb/AliAnalysisTaskMultPbTracks.h
new file mode 100644 (file)
index 0000000..7f8a39e
--- /dev/null
@@ -0,0 +1,63 @@
+#ifndef ALIANALYSISTASKMULTPBTRACKS_H
+#define ALIANALYSISTASKMULTPBTRACKS_H
+
+#include "AliAnalysisTaskSE.h"
+#include "AliESDtrackCuts.h" // if I don't include this, nothing compiles
+
+//-------------------------------------------------------------------------
+//                      AliAnalysisTaskMultPbTracks
+// 
+// 
+//
+//
+// Author: Michele Floris, CERN
+//-------------------------------------------------------------------------
+
+
+class AliESDEvent;
+class AliESDtrackCuts;
+class AliAnalysisMultPbCentralitySelector;
+class AliAnalysisMultPbTrackHistoManager;
+
+
+
+class AliAnalysisTaskMultPbTracks : public AliAnalysisTaskSE {
+
+public:
+
+  AliAnalysisTaskMultPbTracks();
+  AliAnalysisTaskMultPbTracks(const char * name);
+  AliAnalysisTaskMultPbTracks(const AliAnalysisTaskMultPbTracks& obj) ;
+  ~AliAnalysisTaskMultPbTracks();
+  void SetTrackCuts(AliESDtrackCuts * cuts) { fTrackCuts = cuts;}
+  void SetCentralitySelector(AliAnalysisMultPbCentralitySelector * centr) { fCentrSelector=centr;}
+
+  void SetIsMC(Bool_t flag=kTRUE) { fIsMC = flag;}
+  AliAnalysisMultPbTrackHistoManager * GetHistoManager() { return fHistoManager;}
+
+  virtual void   UserCreateOutputObjects();
+  virtual void   UserExec(Option_t *option);
+  virtual void   Terminate(Option_t *);
+  
+  
+
+private:
+
+  //
+  AliESDEvent *  fESD;    //! ESD object  AliVEvent*     fEvent;
+  //  TList * fListHisto;     // list of output object
+  AliAnalysisMultPbTrackHistoManager  * fHistoManager; // wrapper for the list, takes care of merging + histo booking and getters
+  AliAnalysisMultPbCentralitySelector * fCentrSelector; // centrality selector
+  AliESDtrackCuts * fTrackCuts; // track cuts
+  AliESDtrackCuts * fTrackCutsNoDCA; // copy of the previous one, but with no DCA cuts
+  
+  Bool_t fIsMC; // true if processing montecarlo
+
+  AliAnalysisTaskMultPbTracks& operator=(const AliAnalysisTaskMultPbTracks& task);
+  
+  ClassDef(AliAnalysisTaskMultPbTracks, 2)
+
+
+};
+
+#endif
diff --git a/PWG0/multPbPb/correct.C b/PWG0/multPbPb/correct.C
new file mode 100644 (file)
index 0000000..7bbf02d
--- /dev/null
@@ -0,0 +1,399 @@
+class AliAnalysisMultPbTrackHistoManager;
+
+AliAnalysisMultPbTrackHistoManager * hManData = 0;
+AliAnalysisMultPbTrackHistoManager * hManCorr = 0;
+TH2D * hEvStatData = 0;
+TH2D * hEvStatCorr = 0;
+
+void LoadLibs();
+void LoadData(TString dataFolder, TString correctionFolder);
+void SetStyle();
+Double_t CheckSecondaries();
+void ShowAcceptanceInVzSlices() ;
+
+void correct(TString dataFolder = "output/LHC09d_000104892_p4/", TString correctionFolder = "output/LHC10a8_104867/") {
+
+  // Load stuff and set some styles
+  LoadLibs();
+  LoadData(dataFolder,correctionFolder);
+  SetStyle();
+  // ShowAcceptanceInVzSlices();
+  // return;
+
+  // TODO add some cool printout for cuts and centrality selection
+  
+  Double_t fractionWeak = CheckSecondaries();
+  cout << "Rescaling weak correction: " << fractionWeak << endl;
+  
+
+  // Some shorthands
+  // FIXME: Gen should be projected including overflow in z?
+  TH1D * hDataPt   = (TH1D*) hManData->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRec,        -0.5,0.5,-10,10)->Clone("hDataPt");
+  TH1D * hMCPtGen  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen,        -0.5,0.5,-10,10);
+  TH1D * hMCPtRec  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRec,        -0.5,0.5,-10,10);
+  TH1D * hMCPtPri  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim,    -0.5,0.5,-10,10);
+  TH1D * hMCPtSeM  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat,  -0.5,0.5,-10,10);
+  TH1D * hMCPtSeW  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak, -0.5,0.5,-10,10);
+  TH1D * hMCPtFak  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecFake,    -0.5,0.5,-10,10);
+  TCanvas * cdata = new TCanvas ("cData", "Data");  
+
+  hDataPt->Draw();
+  //  hMCPtRec->Draw("same");
+
+  TCanvas * cMC = new TCanvas ("cMC", "Monte Carlo");
+  hMCPtGen ->Draw();
+  hMCPtRec ->Draw("same");
+  hMCPtPri ->Draw("same");
+  hMCPtSeM ->Draw("same");
+  hMCPtSeW ->Draw("same");
+  hMCPtFak ->Draw("same");
+
+  cout << "Fake/All Rec  = " << hMCPtFak->Integral()/hMCPtRec->Integral()  << endl;
+  cout << "SM/All   Rec  = " << hMCPtSeM->Integral()/hMCPtRec->Integral()  << endl;
+  cout << "SW/All   Rec  = " << hMCPtSeW->Integral()/hMCPtRec->Integral()  << endl;
+
+
+  // Compute efficiency and subtract secondaries and fakes after rescaling to data
+  // PRIM_DATA = ALL_DATA - SEC_MC/ALL_MC*ALL_DATA - FAK_MC/ALL_MC*ALL_DATA
+  // TRUE_DATA = PRIM_DATA * GEN_MC/PRIM_MC
+
+  TH1D * hEffPt = (TH1D*) hMCPtPri->Clone("hEffPt");
+  hEffPt->Divide(hMCPtPri,hMCPtGen,1,1,"B");
+
+  TH1D * hCorSeM = (TH1D*) hMCPtSeM->Clone("hEffPt");
+  hCorSeM->Divide(hMCPtSeM,hMCPtRec,1,1,"B");
+  hCorSeM->Multiply(hDataPt);
+
+  TH1D * hCorSeW = (TH1D*) hMCPtSeW->Clone("hEffPt");
+  hCorSeW->Divide(hMCPtSeW,hMCPtRec,1,1,"B");
+  hCorSeW->Scale(fractionWeak);// rescale weak correction
+  hCorSeW->Multiply(hDataPt); 
+
+  TH1D * hCorFak = (TH1D*) hMCPtFak->Clone("hEffPt");
+  hCorFak->Divide(hMCPtFak,hMCPtRec,1,1,"B");
+  hCorFak->Multiply(hDataPt);
+
+  TH1D * hCorrected = (TH1D*) hDataPt->Clone("hCorrected");
+  hCorrected->Add(hCorSeM,-1);
+  hCorrected->Add(hCorSeW,-1);
+  hCorrected->Add(hCorFak,-1);
+  hCorrected->Divide(hEffPt);
+  hCorrected->SetMarkerStyle(kOpenStar);
+
+  cdata->cd();
+  hDataPt->Draw();
+  hCorrected->SetLineColor(kBlack);
+  hCorrected->SetMarkerColor(kBlack);
+  hCorrected->Draw("same");
+  hMCPtGen->DrawClone("same");
+  TF1 * f = GetLevy();
+  hCorrected->Fit(f,"", "same");
+  hCorrected->Fit(f,"IME", "same");
+  cout << "dN/deta = " << f->Integral(0,100) << " +- " << f->IntegralError(0,100) << endl;
+  cout << "Generated dN/deta (correction) = " << hMCPtGen->Integral() << endl;
+  // FIXME: comment this out
+  TH1D * hDataGen  = hManData->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen,        -0.5,0.5,-10,10);
+  cout << "Generated dN/deta (data) =       " << hDataGen->Integral() << endl;
+  hDataGen->Draw("same");
+}
+
+Double_t CheckSecondaries() {
+  // Returns the fraction you need to rescale the secondaries from weak decays for.
+
+  // Some shorthands
+  TH1D * hDataDCA   = hManData->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec       );
+  TH1D * hMCDCAGen  = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen       );
+  TH1D * hMCDCARec  = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec       );
+  TH1D * hMCDCAPri  = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim   );
+  TH1D * hMCDCASW  = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak    );
+  TH1D * hMCDCASM  = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat    );
+  TH1D * hMCDCAFak  = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake   );
+
+  TCanvas * cCumulative = new TCanvas("cDCAculumative","DCA cumulative distributions");
+  GetCumulativeHisto(hMCDCAPri )->Draw();
+  GetRatioIntegratedFractions(hMCDCASW, hMCDCARec  )->Draw("same");
+  GetRatioIntegratedFractions(hMCDCASM, hMCDCARec  )->Draw("same");
+  GetRatioIntegratedFractions(hMCDCAPri,hMCDCARec  )->Draw("same");
+
+
+  TCanvas * c1 = new TCanvas("cDCAFit", "Fit to the DCA distributions");  
+  c1->SetLogy();
+  // Draw all
+  //  hDataDCA->Draw();
+  // //  hMCDCAGen ->Draw("same");
+  //  hMCDCARec ->Draw("same");
+  // hMCDCAPri ->Draw("same");
+  // hMCDCASW ->Draw("same");
+  // hMCDCASM ->Draw("same");
+  // hMCDCAFak ->Draw("same");
+  //  return;
+  
+  TH1D * hMCPrimSMFak = (TH1D*) hMCDCAPri->Clone("hMCPrimSMFak");
+  hMCPrimSMFak->Add(hMCDCASM);
+  hMCPrimSMFak->Add(hMCDCAFak);
+  Int_t ncomp = 2;
+  MyHistoFitter * fitter = new MyHistoFitter (ncomp);
+  //  fitter->SetFCN(fcnhistoNoMCErr);
+  fitter->SetHistoToFit(hDataDCA);
+  fitter->SetMaxXValue(50);
+  fitter->SetHistoComponent(0,hMCPrimSMFak);
+  //fitter->SetHistoComponent(0,hMCPrimSM);
+  fitter->SetHistoComponent(1,hMCDCASW);
+  //  fitter->SetHistoComponent(2,hMCDCAFak);
+  fitter->SetParName(0,"Primaries + Secondaries Material + Fakes");
+  fitter->SetParName(1,"Secondaries (Weak)");
+  //  fitter->SetParName(2,"Fakes");
+  fitter->SetFactor(0,1);
+  fitter->SetFactor(1,1);
+  //  fitter->SetFactor(2,1);
+  fitter->SetFactorLimits(0,0.5,3);
+  fitter->SetFactorLimits(1,0.5,3);
+  //  fitter->SetFactorLimits(2,0.5,3);
+  //  fitter->SetFactorLimits(2,0.5,30);
+  //  fitter->FixFactor(2);// FIXME: THIS IS BUGGY!
+  //  fitter->SetFactorLimits(3,0.9,1.1);
+  fitter->SetScaleHistos();
+  fitter->Fit();
+  fitter->PrintFitResults();
+  fitter->GetHistoToFit()->Draw("");
+  fitter->GetHistoFitted()->Draw("chist,same");
+  for(Int_t icomp = 0; icomp < ncomp; icomp++){
+    fitter->GetHistoComponent(icomp)->Draw("same");
+  }
+  
+  
+  return fitter->GetFactor(1)/fitter->GetFactor(0);
+
+}
+
+void LoadLibs() {
+
+  gSystem->Load("libVMC");
+  gSystem->Load("libTree");
+  gSystem->Load("libSTEERBase");
+  gSystem->Load("libESD");
+  gSystem->Load("libAOD");
+  gSystem->Load("libANALYSIS");
+  gSystem->Load("libANALYSISalice");
+  gSystem->Load("libCORRFW");
+  gSystem->Load("libMinuit");
+  gSystem->Load("libPWG2Spectra");
+  gSystem->Load("libPWG0base"); 
+   
+  gROOT->ProcessLine(gSystem->ExpandPathName(".include $ALICE_ROOT/PWG0"));
+  // Load helper classes
+  // TODO: replace this by a list of TOBJStrings
+  TString taskName("AliAnalysisTaskMultPbTracks.cxx+");
+  TString histoManName("AliAnalysisMultPbTrackHistoManager.cxx+");
+  TString centrName("AliAnalysisMultPbCentralitySelector.cxx+");
+  TString listName("$ALICE_ROOT/PWG1/background/AliHistoListWrapper.cxx+");
+
+  Bool_t debug=0;
+  gROOT->LoadMacro(listName    +(debug?"+g":""));   
+  gROOT->LoadMacro(histoManName+(debug?"+g":""));
+  gROOT->LoadMacro(centrName   +(debug?"+g":""));   
+  gROOT->LoadMacro(taskName    +(debug?"+g":""));   
+
+  // Histo fitter
+  gROOT->LoadMacro("/Users/mfloris/Work/ALICE/ANALYSIS/HistoFitter/fcn.cxx+g");
+  gROOT->LoadMacro("/Users/mfloris/Work/ALICE/ANALYSIS/HistoFitter/MyHistoFitter.cxx+g");
+
+
+}
+
+
+void SetStyle() {
+
+  hManData->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetLineColor(kBlack);    
+  hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen)    ->SetLineColor(kRed  );
+  hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetLineColor(kRed  );
+  hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetLineColor(kGreen);
+  hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetLineColor(kBlue );
+  hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetLineColor(kBlue );
+  hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetLineColor(kCyan );
+
+  hManData->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerColor(kBlack);    
+  hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen)    ->SetMarkerColor(kRed  );
+  hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerColor(kRed  );
+  hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetMarkerColor(kGreen);
+  hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetMarkerColor(kBlue );
+  hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetMarkerColor(kBlue );
+  hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetMarkerColor(kCyan );
+
+  hManData->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerStyle(kFullCircle);    
+  hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen)    ->SetMarkerStyle(kFullSquare);
+  hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerStyle(kOpenSquare);
+  hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetMarkerStyle(kOpenSquare);
+  hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetMarkerStyle(kOpenSquare);
+  hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetMarkerStyle(kOpenCircle);
+  hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetMarkerStyle(kOpenSquare);
+
+ hManData->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetLineColor(kBlack);    
+  hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen)    ->SetLineColor(kRed  );
+  hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetLineColor(kRed  );
+  hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetLineColor(kGreen);
+  hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetLineColor(kBlue );
+  hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetLineColor(kBlue );
+  hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetLineColor(kCyan );
+
+  hManData->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerColor(kBlack);    
+  hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen)    ->SetMarkerColor(kRed  );
+  hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerColor(kRed  );
+  hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetMarkerColor(kGreen);
+  hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetMarkerColor(kBlue );
+  hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetMarkerColor(kBlue );
+  hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetMarkerColor(kCyan );
+
+  hManData->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerStyle(kFullCircle);    
+  hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen)    ->SetMarkerStyle(kFullSquare);
+  hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerStyle(kOpenSquare);
+  hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetMarkerStyle(kOpenSquare);
+  hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetMarkerStyle(kOpenSquare);
+  hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetMarkerStyle(kOpenSquare);
+  hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetMarkerStyle(kOpenSquare);
+
+
+}
+
+void LoadData(TString dataFolder, TString correctionFolder){
+
+  // Get histo manager for data and MC + stat histos
+  TFile * fData = new TFile(dataFolder+"multPbPbtracks.root");
+  TFile * fCorr = new TFile(correctionFolder+"multPbPbtracks.root");
+  TFile * fStatData = new TFile(dataFolder+"event_stat.root");
+  TFile * fStatCorr = new TFile(correctionFolder+"event_stat.root");
+
+  hManData = (AliAnalysisMultPbTrackHistoManager*) fData->Get("histoManager");
+  hManCorr = (AliAnalysisMultPbTrackHistoManager*) fCorr->Get("histoManager");
+  AliESDtrackCuts * cutsData = (AliESDtrackCuts*) fData->Get("AliESDtrackCuts");
+  AliESDtrackCuts * cutsCorr = (AliESDtrackCuts*) fCorr->Get("AliESDtrackCuts");
+  if (cutsData != cutsCorr) {
+    cout << "ERROR: MC and data do not have the same cuts" << endl;
+    // FIXME: exit here
+  }
+  cutsData->Print();
+  hEvStatData = (TH2D*) fStatData->Get("fHistStatistics");
+  hEvStatCorr = (TH2D*) fStatCorr->Get("fHistStatistics");
+
+  // Normalize
+  Int_t irowGoodTrigger = 1;
+  if (hEvStatCorr && hEvStatData) {
+    //  hManData->ScaleHistos(75351.36/1.015);// Nint for run 104892 estimated correcting for the trigger efficiency, multiplied for the physics selection efficiency which I'm not correcting for the time being
+    hManData->ScaleHistos(hEvStatData->GetBinContent(AliPhysicsSelection::kStatAccepted,irowGoodTrigger));
+    hManCorr->ScaleHistos(hEvStatCorr->GetBinContent(AliPhysicsSelection::kStatAccepted,irowGoodTrigger));
+  } else {
+    cout << "WARNING!!! ARBITRARY SCALING" << endl;
+    hManData->ScaleHistos(1000);
+    hManCorr->ScaleHistos(1000);    
+  }
+}
+
+TF1 * GetHagedorn(Float_t norm=68, Float_t pt0=25, Float_t n=13) {
+
+  TF1 * f =0;
+  Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
+  
+  f=new TF1("fHagedorn",Form("(x/sqrt(x*x+%f*%f))*x*[0]*( 1 + x/[1] )^(-[2])",mass,mass),0,10);
+  f->SetParameters(norm, pt0, n);
+  f->SetParLimits(1, 0.01, 10);
+  f->SetParNames("norm", "pt0", "n");
+  f->SetLineWidth(1);
+  return f;
+
+
+}
+
+TF1 * GetMTExp(Float_t norm=68, Float_t t=25) {
+
+  TF1 * f =0;
+  Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
+  
+  f=new TF1("fMTExp",Form("(x/sqrt(x*x+%f*%f))*x*[0]*exp(-sqrt(x*x+%f*%f)/[1])",mass,mass,mass,mass),0,10);
+  f->SetParameters(norm, t);
+  //  f->SetParLimits(1, 0.01);
+  f->SetParNames("norm", "T");
+  f->SetLineWidth(1);
+  return f;
+
+
+}
+
+TF1 * GetLevy(Double_t temp=0.1, Double_t n=7, Double_t norm=10, const char * name="fLevy") {
+
+  char formula[500];
+  Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
+  
+
+  sprintf(formula,"(x/sqrt(x*x+[3]*[3]))*x*( [0]*([1]-1)*([1]-2)  )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) )  ) * ( 1 + (sqrt([3]*[3]+x*x) -[3])/([1]*[2])  )^(-[1])");
+  TF1* f=new TF1(name,formula,0,10);
+  f->SetParameters(norm, n, temp,mass);
+  f->SetParLimits(2, 0.01, 10);
+  f->SetParNames("norm (dN/dy)", "n", "T", "mass");
+  f->FixParameter(3,mass);
+  f->SetLineWidth(1);
+  return f;
+
+
+}
+
+
+TH1D * GetCumulativeHisto (TH1 * h) { 
+  // Returns a cumulative histogram
+  // FIXME: put this method in a tools class
+  TH1D * hInt = h->Clone(TString(h->GetName())+"_Int");
+  hInt->Reset();
+  Float_t integral = h->Integral(-1,-1); // Consider under/over flow!
+  Int_t nbin = h->GetNbinsX();
+  for(Int_t ibin = 1; ibin <= nbin; ibin++){
+    Double_t content = h->Integral(-1,ibin) / integral;
+    hInt->SetBinContent(ibin, content);
+  }
+  return hInt;
+}
+
+TH1D * GetRatioIntegratedFractions (TH1 * hNum, TH1 * hDenum) { 
+  // Returns the the ratio of integrated histograms 
+  // FIXME: put this method in a tools class
+  TH1D * hRatio = hNum->Clone(TString(hNum->GetName())+hDenum->GetName()+"_RatioIntegrated");
+  hRatio->Reset();
+  Int_t nbin = hNum->GetNbinsX();
+  for(Int_t ibin = 1; ibin <= nbin; ibin++){
+    Double_t content = hNum->Integral(-1,ibin) / hDenum->Integral(-1,ibin);// consider underflow
+    hRatio->SetBinContent(ibin, content);
+  }
+  return hRatio;
+}
+
+void ShowAcceptanceInVzSlices() {
+  TCanvas * cvz = new TCanvas("cvz","acc #times eff vs vz");
+  for(Int_t ivz = -10; ivz < 10; ivz+=4){
+    TH1D * hMCPtPri  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim ,   -0.5,0.5,ivz,ivz+4);
+    TH1D * hMCPtGen  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen,        -0.5,0.5,ivz,ivz+4);
+    //    hEff= hMCPtGen;
+    TH1D * hEff = hMCPtPri->Clone(Form("hEff_vz_%d_%d",ivz,ivz+4));
+    hEff->Divide(hMCPtPri,hMCPtGen,1,1,"B");
+    cout << "ivz " << ivz << endl;
+    if(ivz < -9) {
+      cout << "First" << endl;
+      hEff->Draw();
+      // hMCPtGen->Draw();
+      // hMCPtPri->Draw("same");
+    }
+    else {
+      cout << "Same" << endl;
+      hEff->Draw("same");
+      // hMCPtGen->Draw("");
+      // hMCPtPri->Draw("same");
+    }
+    cvz->Update();
+    //    cvz->WaitPrimitive();
+  }
+  
+  //hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim )->Draw();
+   hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoRec )->Draw("");
+  hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoGen     )->Draw("same");      
+
+
+}
diff --git a/PWG0/multPbPb/run.C b/PWG0/multPbPb/run.C
new file mode 100644 (file)
index 0000000..b1a8672
--- /dev/null
@@ -0,0 +1,205 @@
+// TODO:
+// 1. Check cuts for 2010 (Jochen?)
+
+
+enum { kMyRunModeLocal = 0, kMyRunModeCAF};
+
+TChain * GetAnalysisChain(const char * incollection);
+
+void run(Char_t* data, Long64_t nev = -1, Long64_t offset = 0, Bool_t debug = kFALSE, Int_t runMode = 0, Bool_t isMC = 0, const char* option = "",Int_t workers = -1)
+{
+  // runMode:
+  //
+  // 0 local 
+  // 1 proof
+
+  if (nev < 0)
+    nev = 1234567890;
+
+  InitAndLoadLibs(runMode,workers,debug);
+
+  // Create the analysis manager
+  mgr = new AliAnalysisManager;
+
+  // Add ESD handler
+  AliESDInputHandler* esdH = new AliESDInputHandler;
+  // Do I need any of this? 
+  esdH->SetInactiveBranches("AliESDACORDE FMD ALIESDTZERO ALIESDZDC AliRawDataErrorLogs CaloClusters Cascades EMCALCells EMCALTrigger ESDfriend Kinks AliESDTZERO ALIESDACORDE MuonTracks TrdTracks");
+  mgr->SetInputEventHandler(esdH);
+
+  if(isMC) {
+    AliMCEventHandler* handler = new AliMCEventHandler;
+    mgr->SetMCtruthEventHandler(handler);
+  }
+
+  // physics selection
+  gROOT->ProcessLine(".L $ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
+  physicsSelectionTask = AddTaskPhysicsSelection(isMC);
+
+
+
+  // Parse option strings
+  TString optionStr(option);
+  
+  // remove SAVE option if set
+  // This  is copied from a macro by Jan. The reason I kept it is that I may want to pass textual options to the new task at some point
+  Bool_t doSave = kFALSE;
+  TString optionStr(option);
+  if (optionStr.Contains("SAVE"))
+  {
+    optionStr = optionStr(0,optionStr.Index("SAVE")) + optionStr(optionStr.Index("SAVE")+4, optionStr.Length());
+    doSave = kTRUE;
+  }
+
+  AliESDtrackCuts * cuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(kTRUE);
+  TString pathsuffix = "";
+  // cuts->SetPtRange(0.15,0.2);// FIXME pt cut
+  // const char * pathsuffix = "_pt_015_020_nofakes";
+
+  if (optionStr.Contains("ITSsa")) {
+    delete cuts;
+    cuts = AliESDtrackCuts::GetStandardITSPureSATrackCuts2009();
+    cout << ">>>> USING ITS sa tracks" << endl;
+    pathsuffix="ITSsa";
+  }
+
+  if (optionStr.Contains("TPC")) {
+    delete cuts;
+    cuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
+    cout << ">>>> USING TPC only tracks" << endl;
+    pathsuffix="TPC";
+  }
+  
+  
+  // load my task
+  gROOT->ProcessLine(".L $ALICE_ROOT/PWG0/multPbPb/AddTaskMultPbPbTracks.C");
+  AliAnalysisTaskMultPbTracks * task = AddTaskMultPbPbTracks("multPbPbtracks.root", cuts); // kTRUE enables DCA cut
+  task->SetIsMC(isMC);
+   if(isMC) task->GetHistoManager()->SetSuffix("MC");
+  
+  if (!mgr->InitAnalysis()) return;
+       
+  mgr->PrintStatus();
+  
+  if (runMode == kMyRunModeLocal ) {
+    // If running in local mode, create chain of ESD files
+    TChain * chain = GetAnalysisChain(data);
+    chain->Print();
+    mgr->StartAnalysis("local",chain,nev);
+  } else if (runMode == kMyRunModeCAF) {
+    mgr->StartAnalysis("proof",TString(data)+"#esdTree",nev);
+  } else {
+    cout << "ERROR: unknown run mode" << endl;        
+  }
+
+  if (doSave) MoveOutput(data, pathsuffix.Data());
+
+  
+}
+
+
+void MoveOutput(const char * data, const char * suffix = ""){
+
+  TString path("output/");
+  path = path + TString(data).Tokenize("/")->Last()->GetName() + suffix;
+  
+  TString fileName = "multPbPbtracks.root";
+  gSystem->mkdir(path, kTRUE);
+  gSystem->Rename(fileName, path + "/" + fileName);
+  gSystem->Rename("event_stat.root", path + "/event_stat.root");      
+  Printf(">>>>> Moved files to %s", path.Data());
+}  
+
+
+
+TChain * GetAnalysisChain(const char * incollection){
+  // Builds a chain of esd files
+  // incollection can be
+  // - a single root file
+  // - an xml collection of files on alien
+  // - a ASCII containing a list of local root files
+  TChain* analysisChain = 0;
+  // chain
+  analysisChain = new TChain("esdTree");
+  if (TString(incollection).Contains(".root")){
+    analysisChain->Add(incollection);
+  }
+  else if (TString(incollection).Contains("xml")){
+    TGrid::Connect("alien://");
+    TAlienCollection * coll = TAlienCollection::Open (incollection);
+    while(coll->Next()){
+      analysisChain->Add(TString("alien://")+coll->GetLFN());
+    }
+  } else {
+    ifstream file_collect(incollection);
+    TString line;
+    while (line.ReadLine(file_collect) ) {
+      analysisChain->Add(line.Data());
+    }
+  }
+  analysisChain->GetListOfFiles()->Print();
+
+  return analysisChain;
+}
+
+
+void InitAndLoadLibs(Int_t runMode=kMyRunModeLocal, Int_t workers=0,Bool_t debug=0) {
+
+  if (runMode == kMyRunModeCAF)
+  {
+    cout << "Init in CAF mode" << endl;
+    
+    gEnv->SetValue("XSec.GSI.DelegProxy", "2");
+    TProof::Open("alice-caf.cern.ch", workers>0 ? Form("workers=%d",workers) : "");
+    
+    // Enable the needed package
+    gProof->UploadPackage("$ALICE_ROOT/STEERBase");
+    gProof->EnablePackage("$ALICE_ROOT/STEERBase");
+    gProof->UploadPackage("$ALICE_ROOT/ESD");
+    gProof->EnablePackage("$ALICE_ROOT/ESD");
+    gProof->UploadPackage("$ALICE_ROOT/AOD");
+    gProof->EnablePackage("$ALICE_ROOT/AOD");
+    gProof->UploadPackage("$ALICE_ROOT/ANALYSIS");
+    gProof->EnablePackage("$ALICE_ROOT/ANALYSIS");
+    gProof->UploadPackage("$ALICE_ROOT/ANALYSISalice");
+    gProof->EnablePackage("$ALICE_ROOT/ANALYSISalice");
+    gProof->UploadPackage("$ALICE_ROOT/PWG0base");
+    gProof->EnablePackage("$ALICE_ROOT/PWG0base");
+  }
+  else
+  {
+    cout << "Init in Local mode" << endl;
+
+    gSystem->Load("libVMC");
+    gSystem->Load("libTree");
+    gSystem->Load("libSTEERBase");
+    gSystem->Load("libESD");
+    gSystem->Load("libAOD");
+    gSystem->Load("libANALYSIS");
+    gSystem->Load("libANALYSISalice");
+    gSystem->Load("libPWG0base");
+    
+    gROOT->ProcessLine(gSystem->ExpandPathName(".include $ALICE_ROOT/PWG0"));
+  }
+  // Load helper classes
+  // TODO: replace this by a list of TOBJStrings
+  TString taskName("AliAnalysisTaskMultPbTracks.cxx+");
+  TString histoManName("AliAnalysisMultPbTrackHistoManager.cxx+");
+  TString centrName("AliAnalysisMultPbCentralitySelector.cxx+");
+  TString listName("AliHistoListWrapper.cxx+");
+
+  // Create, add task
+  if (runMode == kMyRunModeCAF) {
+    gProof->Load(listName+(debug?"+g":""));   
+    gProof->Load(histoManName+(debug?"+g":""));
+    gProof->Load(centrName+(debug?"+g":""));   
+    gProof->Load(taskName+(debug?"+g":""));
+  } else {
+    gROOT->LoadMacro(listName+(debug?"+g":""));   
+    gROOT->LoadMacro(histoManName+(debug?"+g":""));
+    gROOT->LoadMacro(centrName+(debug?"+g":""));   
+    gROOT->LoadMacro(taskName+(debug?"+g":""));    
+  }
+
+
+}
diff --git a/PWG0/multPbPb/run.sh b/PWG0/multPbPb/run.sh
new file mode 100755 (executable)
index 0000000..caf14e5
--- /dev/null
@@ -0,0 +1,104 @@
+#!/bin/bash
+
+# defaults
+isMC=kFALSE
+run=no
+correct=no
+nev=-1
+offset=0
+debug=kFALSE
+runmode=1
+dataset=/PWG3/zampolli/run000104867_90_92_pass4
+ropt="-l"
+option="SAVE"
+workers=26
+analysismode=9; #SPD + field on
+
+give_help() {
+
+cat <<ENDOFGUIDE
+This scripts runs the mupliplicity analysis.
+
+Available options:
+ Mode control, at least one of the following options should be used
+  -r <mode>                    Run the task
+                               Modes [default=$runmode]:
+                                  0 local
+                                  1 caf    
+  -c                           Run the correction
+ Proof settings
+  -w nworkers                  Set the number of worker nodes
+  -n <nev>                     Number of events to be analized 
+ Misc
+  -d <dataset>                 Dataset or data collection (according to run mode) [default=$dataset]
+  -o <option>                  Misc option [default=$option]
+                               Available options: 
+                                - SAVE:  move results to a different output folder
+                                - ITSsa: Use ITSsa tracks
+                                - TPC:   Use TPC only tracks
+  -t <option>                  Command line option for root [defaul=$ropt]
+  -m                           Use this to run on Monte Carlo
+  -g                           Debug mode
+  -h                           This help
+ENDOFGUIDE
+
+}
+
+while getopts "r:cgmd:o:w:n:" opt; do
+  case $opt in
+    r)
+      run=yes
+      runmode=$OPTARG
+      ;;      
+    c)
+      correct=yes
+      ;;      
+    g)
+      debug=kTRUE;
+      ;;
+    m)
+      isMC=kTRUE
+      ;;
+    d)
+      dataset=$OPTARG
+     ;;
+    o)
+      option=$OPTARG
+      ;;
+    t)
+      ropt=$OPTARG
+      ;;
+    w) 
+      workers=$OPTARG
+      ;;
+    n) 
+      nev=$OPTARG
+      ;;
+    h)
+      give_help
+      exit 1
+      ;;
+    \?)
+      echo "Invalid option: -$OPTARG" >&2
+      give_help
+      exit 1
+      ;;
+    :)
+      echo "Option -$OPTARG requires an argument." >&2
+      give_help
+      exit 1
+      ;;
+  esac
+done
+
+if [ "$run" = "$correct" ]
+    then 
+    echo "One and only one option between -r and -c must be selected"
+    give_help
+    exit 1
+fi
+
+if [ "$run" = "yes" ]
+    then
+    root $ropt run.C\(\"$dataset\",$nev,$offset,$debug,$runmode,$isMC,\"$option\",$workers\)
+fi