Main Updates:
authormfloris <mfloris@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 5 Nov 2010 09:07:04 +0000 (09:07 +0000)
committermfloris <mfloris@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 5 Nov 2010 09:07:04 +0000 (09:07 +0000)
- AliAnalysisMultPbCentralitySelector: when running on data, uses
  AliESDcentrality by default, when running on MC a cut on the
  selected tracks

- AliAnalysisMultPbTrackHistoManager&AliAnalysisTaskMultPbTracks.cxx:
  added histograms for particle species, multiplicity histogram,
  changed pt binning (same as 900 GeV pp now)

- AliAnalysisTaskTriggerStudy: added venn-like histo

- AddTaskMultPbPbTracks.C: forgot in the initial commit

- run.C&run.sh: changes to run on grid (not yet working), loading
  extra classes through a TList, possibility to give a custom prefix
  to the output folder.

- correct.C: minor chages for specific checks

PWG0/multPbPb/AddTaskMultPbPbTracks.C [new file with mode: 0644]
PWG0/multPbPb/AliAnalysisMultPbCentralitySelector.cxx
PWG0/multPbPb/AliAnalysisMultPbCentralitySelector.h
PWG0/multPbPb/AliAnalysisMultPbTrackHistoManager.cxx
PWG0/multPbPb/AliAnalysisMultPbTrackHistoManager.h
PWG0/multPbPb/AliAnalysisTaskMultPbTracks.cxx
PWG0/multPbPb/AliAnalysisTaskMultPbTracks.h
PWG0/multPbPb/AliAnalysisTaskTriggerStudy.cxx
PWG0/multPbPb/correct.C
PWG0/multPbPb/run.C
PWG0/multPbPb/run.sh

diff --git a/PWG0/multPbPb/AddTaskMultPbPbTracks.C b/PWG0/multPbPb/AddTaskMultPbPbTracks.C
new file mode 100644 (file)
index 0000000..a99b2b5
--- /dev/null
@@ -0,0 +1,66 @@
+AliAnalysisTaskMultPbTracks * AddTaskMultPbPbTracks(const char * outfilename, AliESDtrackCuts * esdTrackCuts = 0)
+{
+  // TODO: add some parameters to set the centrality for this task, and maybe the name of the task
+  // TODO: shall I use the same file and different dirs for the different centralities?
+
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) {
+    ::Error("AddTaskPhysicsSelection", "No analysis manager to connect to.");
+    return NULL;
+  }  
+  
+  // Check the analysis type using the event handlers connected to the analysis manager.
+  //==============================================================================
+  if (!mgr->GetInputEventHandler()) {
+    ::Error("AddTaskPhysicsSelection", "This task requires an input event handler");
+    return NULL;
+  }
+  TString inputDataType = mgr->GetInputEventHandler()->GetDataType(); // can be "ESD" or "AOD"
+  
+  if (inputDataType != "ESD") {
+    Printf("ERROR! This task can only run on ESDs!");
+  }
+
+  // Configure analysis
+  //===========================================================================
+    
+    
+
+  AliAnalysisTaskMultPbTracks *task = new AliAnalysisTaskMultPbTracks("TaskMultPbTracks");
+  mgr->AddTask(task);
+  
+  // Set Cuts
+  if (!esdTrackCuts)
+    {
+      printf("ERROR: esdTrackCuts could not be created\n");
+      return;
+    }  
+  task->SetTrackCuts(esdTrackCuts);
+  
+  // TODO:
+  // IO into folders in a file?
+  // FIXME:
+  // Add centrality selection configuration (included arguments)
+
+  // Set I/O
+  AliAnalysisDataContainer *cinput0 = mgr->GetCommonInputContainer();
+  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("cmultPbTracksOutHM",
+                                                           AliAnalysisMultPbTrackHistoManager::Class(),
+                                                           AliAnalysisManager::kOutputContainer,
+                                                           outfilename);
+  AliAnalysisDataContainer *coutput2 = mgr->CreateContainer("cmultPbTracksOutCT",
+                                                           AliESDtrackCuts::Class(),
+                                                           AliAnalysisManager::kOutputContainer,
+                                                           outfilename);
+  // AliAnalysisDataContainer *coutput3 = mgr->CreateContainer("cmultPbTracksOutCM",
+  //                                                       AliAnalysisMultPbCentralitySelector::Class(),
+  //                                                       AliAnalysisManager::kOutputContainer,
+  //                                                       outfilename);
+
+  mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
+  mgr->ConnectOutput(task,1,coutput1);
+  mgr->ConnectOutput(task,2,coutput2);
+  //  mgr->ConnectOutput(task,3,coutput3);
+
+  return task;
+}   
index 9a3bd65..6029b87 100644 (file)
@@ -1,4 +1,45 @@
 #include "AliAnalysisMultPbCentralitySelector.h"
+#include "AliESDtrackCuts.h"
+#include "AliESDCentrality.h"
+#include "AliESDEvent.h"
+#include "AliLog.h"
 
 ClassImp(AliAnalysisMultPbCentralitySelector)
 
+Bool_t AliAnalysisMultPbCentralitySelector::IsCentralityBinSelected(AliESDEvent* aEsd, AliESDtrackCuts * trackCuts) {
+
+  // Centrality selection
+  // On MC cuts on the number of good tracks,
+  // On data cuts using AliESDCentrality and the cut requested in ntracks
+  if(fIsMC) {
+    if(!trackCuts){
+      AliFatal("Track cuts object is invalid");
+    }
+    if (trackCuts->CountAcceptedTracks(aEsd) < fMultMin) return kFALSE;
+    if (trackCuts->CountAcceptedTracks(aEsd) > fMultMax) return kFALSE;                                                       
+  }
+
+  AliESDCentrality *centrality = aEsd->GetCentrality();
+  if(!centrality) {
+    AliFatal("Centrality object not available"); 
+  }
+  else {
+    Int_t centrBin = centrality->GetCentralityClass5(fCentrEstimator.Data()) ;    
+    if (centrBin != fCentrBin && fCentrBin != -1) return kFALSE;
+  }
+
+  return kTRUE;
+
+}
+
+void AliAnalysisMultPbCentralitySelector::Print(Option_t* option ) const {
+  // Print some information
+
+  Printf("AliAnalysisMultPbCentralitySelector");
+  Printf(" - Centrality estimator [%s]",fCentrEstimator.Data());
+
+  if ( fIsMC ) {    
+    Printf("Running on Monte Carlo, actual cut was on tracks multiplicity [%d - %d]",fMultMin,fMultMax);    
+  } 
+  
+}
index eb07d7c..0053722 100644 (file)
@@ -24,30 +24,34 @@ class TH1F;
 class TCollection;
 class AliTriggerAnalysis;
 class AliAnalysisTaskSE;
-
+class AliESDtrackCuts;
 
 class AliAnalysisMultPbCentralitySelector : public AliAnalysisCuts
 {
 public:
 
-  AliAnalysisMultPbCentralitySelector() : fIsMC (0) {;}
+  AliAnalysisMultPbCentralitySelector() : fIsMC (0), fCentrEstimator(""), fCentrBin(-1), fMultMin(0), fMultMax(1000000) {;}
   virtual ~AliAnalysisMultPbCentralitySelector(){}
     
   // AliAnalysisCuts interface
-  virtual UInt_t GetSelectionMask(const TObject* obj) { return IsCentralityBinSelected((const AliESDEvent*) obj); }
+  virtual UInt_t GetSelectionMask(const TObject* obj) { return (UInt_t) IsCentralityBinSelected((AliESDEvent*) obj, NULL); }
   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;}
+  virtual Bool_t IsSelected(TObject* obj)  {return (UInt_t) IsCentralityBinSelected ( (AliESDEvent*) obj, NULL);}
     
-  void SetAnalyzeMC(Bool_t flag = kTRUE) { fIsMC = flag; }
+  Bool_t IsCentralityBinSelected(AliESDEvent* aEsd, AliESDtrackCuts * trackCuts);
     
-  virtual void Print(Option_t* option = "") const { Printf ("Multiplitity Selector [AliAnalysisMultPbCentralitySelector] [%s]", option);}
+  void SetAnalyzeMC(Bool_t flag = kTRUE, Double_t multMin = 0, Double_t multMax=10000) { fIsMC = flag; fMultMin = multMin; fMultMax = multMax; }
+  void SetCentralityEstimator(const char * estimator) { fCentrEstimator = estimator; }
+
+  virtual void Print(Option_t* option = "") const ;
   virtual Long64_t Merge(TCollection* list){list->GetEntries();return 0;}
   
 protected:
   Bool_t fIsMC;             // flag if MC is analyzed
-
+  TString fCentrEstimator;  // Centrality estimator for AliESDCentrality
+  TString fCentrBin; // centrality bin to be selected
+  Int_t fMultMin ; // Minimum multiplicity, because on MC we cut on tracks rather than on the estimator  
+  Int_t fMultMax ; // Maximum multiplicity, because on MC we cut on tracks rather than on the estimator  
   ClassDef(AliAnalysisMultPbCentralitySelector, 1)
     
   private:
index dda337c..ca2aa7a 100644 (file)
@@ -4,6 +4,8 @@
 #include "TH3D.h"
 #include "TH1I.h"
 #include "TROOT.h"
+#include "TMCProcess.h"
+
 #include <iostream>
 using namespace std;
 ClassImp(AliAnalysisMultPbTrackHistoManager)
@@ -12,6 +14,7 @@ const char * AliAnalysisMultPbTrackHistoManager::kStatStepNames[]     = { "All E
 const char * AliAnalysisMultPbTrackHistoManager::kHistoPtEtaVzNames[] = { "hGenPtEtaVz", "hRecPtEtaVz", "hRecPtEtaVzPrim", 
                                                                          "hRecPtEtaVzSecWeak", "hRecPtEtaVzSecMaterial", "hRecPtEtaVzFake"};
 const char * AliAnalysisMultPbTrackHistoManager::kHistoDCANames[]     = { "hGenDCA", "hRecDCA", "hRecDCAPrim", "hRecDCASecWeak","hRecDCASecMaterial", "hRecDCAFake"};
+const char * AliAnalysisMultPbTrackHistoManager::kHistoPrefix[]     = { "hGen", "hRec", "hRecPrim", "hRecSecWeak","hRecSecMaterial", "hRecFake"};
 
 
 
@@ -58,6 +61,20 @@ TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoDCA(Histo_t id) {
 
 }
 
+TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoMult(Histo_t id) {
+  // Returns a 3D histo of Pt/eta/vtx. It it does not exist, books it.
+
+  TString name = kHistoPrefix[id];
+  name += "Mult";
+  TH1D * h = (TH1D*) GetHisto(name.Data());
+  if (!h) {
+    h = BookHistoMult(name.Data(), Form("Multiplicity distribution (%s)",kHistoPrefix[id]));
+  }
+
+  return h;
+
+}
+
 TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoPt (Histo_t id, 
                                                       Float_t minEta, Float_t maxEta, 
                                                       Float_t minVz, Float_t maxVz, 
@@ -84,7 +101,8 @@ TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoPt (Histo_t id,
 
   
   if (gROOT->FindObjectAny(hname.Data())){
-    AliError(Form("An object called %s already exists",hname.Data()));
+    AliError(Form("An object called %s already exists,adding suffix",hname.Data()));
+    hname += "_2";
   }
 
   TH1D * h = h3D->ProjectionX(hname.Data(), min1, max1, min2, max2, "E");
@@ -119,7 +137,8 @@ TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoVz (Histo_t id,
   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()));
+    AliError(Form("An object called %s already exists, adding suffix",hname.Data()));
+    hname+="_2";
   }
 
   TH1D * h = h3D->ProjectionZ(hname.Data(), min1, max1, min2, max2, "E");
@@ -153,7 +172,8 @@ TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoEta (Histo_t id,
   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()));
+    AliError(Form("An object called %s already exists, adding suffix",hname.Data()));
+    hname+="_2";
   }
 
   TH1D * h = h3D->ProjectionY(hname.Data(), min1, max1, min2, max2, "E");
@@ -161,6 +181,36 @@ TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoEta (Histo_t id,
   return h;
 }
 
+TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoSpecies(Histo_t id) {
+
+  // Returns histogram with particle specties
+
+  TString name = TString(kHistoPrefix[id])+"_Species";
+
+  TH1D * h =  (TH1D*) GetHisto(name);
+  if (!h) {
+    name+=fHNameSuffix;
+    Bool_t oldStatus = TH1::AddDirectoryStatus();
+    TH1::AddDirectory(kFALSE);
+
+    AliInfo(Form("Booking histo %s",name.Data()));
+
+    h = new TH1D (name.Data(), Form("Particle species (%s)",kHistoPrefix[id]), kPNoProcess+1, -0.5, kPNoProcess+1-0.5);                         
+    Int_t nbin = kPNoProcess+1;
+    for(Int_t ibin = 0; ibin < nbin; ibin++){
+      h->GetXaxis()->SetBinLabel(ibin+1,TMCProcessName[ibin]);      
+    }
+    TH1::AddDirectory(oldStatus);
+    fList->Add(h);
+
+
+  }
+  return h;
+  
+
+}
+
+
 
 TH1I * AliAnalysisMultPbTrackHistoManager::GetHistoStats() {
   // Returns histogram with event statistiscs (processed events at each step)
@@ -189,7 +239,7 @@ void AliAnalysisMultPbTrackHistoManager::ScaleHistos(Double_t nev, Option_t * op
     if (!h->InheritsFrom("TH1")) {
       AliFatal (Form("%s does not inherits from TH1, cannot scale",h->GetName()));
     }
-    cout << "Scaling " << h->GetName() << " " << nev << endl;
+    AliInfo(Form("Scaling %s, nev %2.2f", h->GetName(), nev));
     
     h->Scale(1./nev,option);
   }
@@ -208,15 +258,39 @@ TH3D * AliAnalysisMultPbTrackHistoManager::BookHistoPtEtaVz(const char * name, c
 
   AliInfo(Form("Booking %s",hname.Data()));
   
+  // Binning from Jacek task
+  const Int_t nptbins = 68;
+  const Double_t binsPt[] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,4.5,5.0,5.5,6.0,6.5,7.0,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0,16.0,18.0,20.0,22.0,24.0,26.0,28.0,30.0,32.0,34.0,36.0,40.0,45.0,50.0};
+  
+  const Int_t netabins=20;
+  Double_t binsEta[netabins+1];
+  Float_t minEta = -1;
+  Float_t maxEta =  1;
+  Float_t etaStep = (maxEta-minEta)/netabins;
+  for(Int_t ibin = 0; ibin < netabins; ibin++){    
+    binsEta[ibin]   = minEta + ibin*etaStep;
+    binsEta[ibin+1] = minEta + ibin*etaStep + etaStep;
+  }
+
+  const Int_t nvzbins=10;
+  Double_t binsVz[nvzbins+1];
+  Float_t minVz = -10;
+  Float_t maxVz =  10;
+  Float_t vzStep = (maxVz-minVz)/nvzbins;
+  for(Int_t ibin = 0; ibin < nvzbins; ibin++){    
+    binsVz[ibin]   = minVz + ibin*vzStep;
+    binsVz[ibin+1] = minVz + ibin*vzStep + vzStep;
+  }
 
   TH3D * h = new TH3D (hname,title, 
-                      50,0,10, // Pt
-                      20,-1, 1,   // Eta
-                      10,-10,10 // Vz
+                      nptbins,  binsPt,
+                      netabins, binsEta,
+                      nvzbins,  binsVz
                       );
 
-  h->SetYTitle("#eta");
   h->SetXTitle("p_{T}");
+  h->SetYTitle("#eta");
   h->SetZTitle("V_{z} (cm)");
   h->Sumw2();
   
@@ -248,6 +322,28 @@ TH1D * AliAnalysisMultPbTrackHistoManager::BookHistoDCA(const char * name, const
   TH1::AddDirectory(oldStatus);
   return h;
 }
+TH1D * AliAnalysisMultPbTrackHistoManager::BookHistoMult(const char * name, const char * title) {
+  // Books a multiplicity 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, 600, 0,6000);
+
+  h->SetXTitle("N tracks");
+  h->Sumw2();
+  
+  fList->Add(h);
+
+  TH1::AddDirectory(oldStatus);
+  return h;
+}
 
 TH1I * AliAnalysisMultPbTrackHistoManager::BookHistoStats() {
   // Books histogram with event statistiscs (processed events at each step)
index 4e74ce6..6f0591b 100644 (file)
@@ -42,6 +42,9 @@ public:
 
   TH1I * GetHistoStats();
   TH1D * GetHistoDCA(Histo_t id);
+  TH1D * GetHistoMult(Histo_t id);
+
+  TH1D * GetHistoSpecies(Histo_t id);
 
   // Misch utils
   void ScaleHistos (Double_t nev, Option_t * option="");
@@ -52,6 +55,7 @@ public:
   TH3D * BookHistoPtEtaVz(const char * name, const char * title);
   TH1D * BookHistoDCA(const char * name, const char * title);
   TH1I * BookHistoStats();
+  TH1D * BookHistoMult(const char * name, const char * title);
 
   // 
   TH1 * GetHisto(const char * name);
@@ -61,6 +65,7 @@ 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 
+  static const char * kHistoPrefix[];   // prefix for histo names // FIXME: remove the others and keep only this 
   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);
index 345417b..11532c8 100644 (file)
@@ -3,7 +3,8 @@
 // Author: Michele Floris, CERN
 // TODO:
 // - Add chi2/cluster plot for primary, secondaries and fakes
-
+// FIXME:
+// - re-implement centrality estimator to cut on tracks and multiplicity
 
 #include "AliAnalysisTaskMultPbTracks.h"
 #include "AliESDInputHandler.h"
@@ -32,7 +33,7 @@ AliAnalysisTaskMultPbTracks::AliAnalysisTaskMultPbTracks()
 
   DefineOutput(1, AliAnalysisMultPbTrackHistoManager::Class());
   DefineOutput(2, AliESDtrackCuts::Class());
-  //  DefineOutput(2, TH1I::Class());
+  //  DefineOutput(3, AliAnalysisMultPbCentralitySelector::Class());
 
   fHistoManager = new AliAnalysisMultPbTrackHistoManager("histoManager","Hitograms, Multiplicity, Track analysis");
   if(fIsMC) fHistoManager->SetSuffix("MC");
@@ -48,6 +49,7 @@ AliAnalysisTaskMultPbTracks::AliAnalysisTaskMultPbTracks(const char * name)
 
   DefineOutput(1, AliAnalysisMultPbTrackHistoManager::Class());
   DefineOutput(2, AliESDtrackCuts::Class());
+  //  DefineOutput(3, AliAnalysisMultPbCentralitySelector::Class());
 
   fHistoManager = new AliAnalysisMultPbTrackHistoManager("histoManager","Hitograms, Multiplicity, Track analysis");
   if(fIsMC) fHistoManager->SetSuffix("MC");
@@ -102,10 +104,12 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
   /* PostData(0) is taken care of by AliAnalysisTaskSE */
   PostData(1,fHistoManager);
   PostData(2,fTrackCuts);
+  //  PostData(3,fCentralityEstimator);
 
   // Cache histogram pointers
-  static TH3D * hTracks[AliAnalysisMultPbTrackHistoManager::kNHistos];
-  static TH1D * hDCA   [AliAnalysisMultPbTrackHistoManager::kNHistos];
+  static TH3D * hTracks  [AliAnalysisMultPbTrackHistoManager::kNHistos];
+  static TH1D * hDCA     [AliAnalysisMultPbTrackHistoManager::kNHistos];
+  static TH1D * hNTracks [AliAnalysisMultPbTrackHistoManager::kNHistos];
   hTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]        = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen       );
   hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec]        = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec       );
   hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]    = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim   );
@@ -120,6 +124,7 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
   hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]  = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat );
   hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak] = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak);
 
+  hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]        = fHistoManager->GetHistoMult(AliAnalysisMultPbTrackHistoManager::kHistoGen       );
 
   fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
   if (strcmp(fESD->ClassName(),"AliESDEvent")) {
@@ -130,13 +135,15 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
   fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatAll);
 
   // Centrality selection
+  // Bool_t isCentralitySelected = fCentrSelector->IsSelected(fESD);  
+  // if(!isCentralitySelected) return;
+
   AliESDCentrality *centrality = fESD->GetCentrality();
   if(!centrality) {
-    AliError("Centrality object not available"); // FIXME AliFatal here after debug
+    AliFatal("Centrality object not available"); 
   }
   else {
     Int_t centrBin = centrality->GetCentralityClass5(fCentralityEstimator.Data()) ;
-    cout << fCentralityEstimator.Data() << " BIN: " << centrBin << endl;
     
     if (centrBin != fCentrBin && fCentrBin != -1) return;
   }
@@ -150,27 +157,25 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
       AliError("No MC info found");
     } else {
       
-      //loop on the MC event
-      //      Int_t nMCTracks = fMCEvent->GetNumberOfTracks();
-      Int_t offset    = fMCEvent->GetPrimaryOffset();
-      Int_t nMCTracks = fMCEvent->GetNumberOfPrimaries()+offset;
-      for (Int_t ipart=offset; ipart<nMCTracks; ipart++) { 
+      //loop on the MC event, only over primaries, which are always
+      //      the first in stack.
+      Int_t nMCTracks = fMCEvent->GetNumberOfPrimaries();
+      Int_t nPhysicalPrimaries = 0;
+      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 (uncomment below)
-       if(!fMCEvent->Stack()->IsPhysicalPrimary(ipart)) continue;
-
        //check if current particle is a physical primary
-       // Bool_t physprim=fMCEvent->IsPhysicalPrimary(label);
-       // if (!physprim) continue;
-       // if (!track) return kFALSE;
-       // Bool_t transported = mcPart->Particle()->TestBit(kTransportBit);
-       // if(!transported) return kFALSE;
+       if(!IsPhysicalPrimaryAndTransportBit(ipart)) continue;
  
+       nPhysicalPrimaries++;
+       // Fill species histo
+       fHistoManager->GetHistoSpecies(AliAnalysisMultPbTrackHistoManager::kHistoGen)->Fill(mcPart->Particle()->GetUniqueID());
+
+       
        // Get MC vertex
        //FIXME: which vertex do I take for MC?
        TArrayF   vertex;
@@ -181,9 +186,12 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
        hTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]->Fill(mcPart->Pt(),mcPart->Eta(),zv);
        
       }
+      hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]->Fill(nPhysicalPrimaries);
+
     }
   }
   
+
   // FIXME: shall I take the primary vertex?
   const AliESDVertex* vtxESD = fESD->GetPrimaryVertex();
   if(!vtxESD) return;
@@ -242,7 +250,7 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
     if (fIsMC) {
       //      Int_t label = TMath::Abs(esdTrack->GetLabel()); // no fakes!!!
       Int_t label = esdTrack->GetLabel(); // 
-      AliMCParticle *mcPart  = (AliMCParticle*)fMCEvent->GetTrack(label);
+      AliMCParticle *mcPart  = label < 0 ? 0 : (AliMCParticle*)fMCEvent->GetTrack(label);
       if (!mcPart)  {
        if(accepted)
          hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
@@ -250,8 +258,9 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
          hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]->Fill(weightedDCA);
       }
       else {
-       if(fMCEvent->Stack()->IsPhysicalPrimary(label)) {
-         // FIXME add kTransportBit
+       if(IsPhysicalPrimaryAndTransportBit(label)) {
+         // Fill species histo
+         fHistoManager->GetHistoSpecies(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->Fill(mcPart->Particle()->GetUniqueID());
          if(accepted)
            hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
          if(acceptedNoDCA)
@@ -266,11 +275,13 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
            mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));
          }
          if(mfl==3){ // strangeness
+           fHistoManager->GetHistoSpecies(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak)->Fill(mcPart->Particle()->GetUniqueID());
            if(accepted)
              hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
            if(acceptedNoDCA)
              hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak]->Fill(weightedDCA);      
          }else{ // material
+           fHistoManager->GetHistoSpecies(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat)->Fill(mcPart->Particle()->GetUniqueID());
            if(accepted)
              hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
            if(acceptedNoDCA)
@@ -295,3 +306,13 @@ void   AliAnalysisTaskMultPbTracks::Terminate(Option_t *){
 }
 
 
+Bool_t AliAnalysisTaskMultPbTracks::IsPhysicalPrimaryAndTransportBit(Int_t ipart) {
+
+  Bool_t physprim=fMCEvent->IsPhysicalPrimary(ipart);
+  if (!physprim) return kFALSE;
+  Bool_t transported = ((AliMCParticle*)fMCEvent->GetTrack(ipart))->Particle()->TestBit(kTransportBit);
+  if(!transported) return kFALSE;
+
+  return kTRUE;
+
+}
index 8999008..3bd7979 100644 (file)
@@ -29,12 +29,14 @@ public:
   AliAnalysisTaskMultPbTracks(const char * name);
   AliAnalysisTaskMultPbTracks(const AliAnalysisTaskMultPbTracks& obj) ;
   ~AliAnalysisTaskMultPbTracks();
+  //void SetCentralitySelector(AliAnalysisMultPbCentralitySelector * centr) { fCentrSelector=centr;}
   void SetTrackCuts(AliESDtrackCuts * cuts) { fTrackCuts = cuts;}
   void SetCentralityBin(Int_t bin = 0) { fCentrBin = bin; }
   void SetCentralityEstimator(const char * centr) { fCentralityEstimator = centr; }
 
   void SetIsMC(Bool_t flag=kTRUE) { fIsMC = flag;}
   AliAnalysisMultPbTrackHistoManager * GetHistoManager() { return fHistoManager;}
+  Bool_t IsPhysicalPrimaryAndTransportBit(Int_t ipart) ;
 
   virtual void   UserCreateOutputObjects();
   virtual void   UserExec(Option_t *option);
@@ -48,6 +50,7 @@ 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
   Int_t fCentrBin; // centrality bin selected (5% XS percentiles)
   TString fCentralityEstimator; // Name of the centrality estimator, for AliESDCentrality
   AliESDtrackCuts * fTrackCuts; // track cuts
index 4b88d3a..ad5f5fb 100644 (file)
@@ -94,10 +94,6 @@ void AliAnalysisTaskTriggerStudy::UserExec(Option_t *)
     AliFatal("Not processing ESDs");
   }
 
-  // FIXME: two options here: either we add a loop setting the name of
-  // the histos with the trigger class in them (there may be a global
-  // SetHistoNamePrefix) or I put a cut to select only collision
-  // classes
   
   // get the multiplicity object
   const AliMultiplicity* mult = fESD->GetMultiplicity();
@@ -209,19 +205,9 @@ void AliAnalysisTaskTriggerStudy::UserExec(Option_t *)
   //   // We don't care about neutrals and non-physical primaries
   //   if(mcPart->Charge() == 0) continue;
 
-  //   // FIXME: add kTransportBit (uncomment below)
-  //   if(!fMCEvent->Stack()->IsPhysicalPrimary(ipart)) continue;
-
-  //   //check if current particle is a physical primary
-  //   // Bool_t physprim=fMCEvent->IsPhysicalPrimary(label);
-  //   // if (!physprim) continue;
-  //   // if (!track) return kFALSE;
-  //   // Bool_t transported = mcPart->Particle()->TestBit(kTransportBit);
-  //   // if(!transported) return kFALSE;
+  // PHYSICAL PRIMARY
   //   // Get MC vertex
-  //   //FIXME: which vertex do I take for MC?
-  //   TArrayF   vertex;
+    //         TArrayF   vertex;
   //   fMCEvent->GenEventHeader()->PrimaryVertex(vertex);
   //   Float_t zv = vertex[2];
   //   //      Float_t zv = vtxESD->GetZ();
@@ -313,19 +299,19 @@ void AliAnalysisTaskTriggerStudy::FillTriggerOverlaps (const char * name, const
   
   Bool_t fo2 = nFastOrOffline>=2;
 
-  if(fo2 && v0A && v0C && OM2)     {cout << "Bin5: " << ibin << endl;  h->Fill(1);}
-  if(fo2 && !v0A && v0C && OM2)    {cout << "Bin6: " << ibin << endl;  h->Fill(2);}
-  if(fo2 && v0A && !v0C && OM2)    {cout << "Bin7: " << ibin << endl;  h->Fill(3);}
-  if(fo2 && v0A && v0C && !OM2)    {cout << "Bin8: " << ibin << endl;  h->Fill(4);}
-  if(fo2 && v0A && !v0C && !OM2)   {cout << "Bin9: " << ibin << endl;  h->Fill(5);}
-  if(fo2 && !v0A && v0C && !OM2)   {cout << "Bin10: " << ibin << endl; h->Fill(6);}
-  if(fo2 && !v0A && !v0C && OM2)   {cout << "Bin11: " << ibin << endl; h->Fill(7);}
-  if(!fo2 && v0A && !v0C && OM2)   {cout << "Bin12: " << ibin << endl; h->Fill(8);}
-  if(!fo2 && !v0A && v0C && OM2)   {cout << "Bin13: " << ibin << endl; h->Fill(9);}
-  if(!fo2 && v0A && v0C && !OM2)   {cout << "Bin14: " << ibin << endl; h->Fill(10);}
-  if(fo2 && !v0A && !v0C && !OM2)  {cout << "Bin1: " << ibin << endl;  h->Fill(11);}
-  if(!fo2 && v0A && !v0C && !OM2)  {cout << "Bin2: " << ibin << endl;  h->Fill(12);}
-  if(!fo2 && !v0A && v0C && !OM2)  {cout << "Bin3: " << ibin << endl;  h->Fill(13);}
-  if(!fo2 && !v0A && !v0C && OM2)  {cout << "Bin4: " << ibin << endl;  h->Fill(14);}
+  if(fo2 && v0A && v0C && OM2)     { h->Fill(1);}
+  if(fo2 && !v0A && v0C && OM2)    { h->Fill(2);}
+  if(fo2 && v0A && !v0C && OM2)    { h->Fill(3);}
+  if(fo2 && v0A && v0C && !OM2)    { h->Fill(4);}
+  if(fo2 && v0A && !v0C && !OM2)   { h->Fill(5);}
+  if(fo2 && !v0A && v0C && !OM2)   { h->Fill(6);}
+  if(fo2 && !v0A && !v0C && OM2)   { h->Fill(7);}
+  if(!fo2 && v0A && !v0C && OM2)   { h->Fill(8);}
+  if(!fo2 && !v0A && v0C && OM2)   { h->Fill(9);}
+  if(!fo2 && v0A && v0C && !OM2)   { h->Fill(10);}
+  if(fo2 && !v0A && !v0C && !OM2)  { h->Fill(11);}
+  if(!fo2 && v0A && !v0C && !OM2)  { h->Fill(12);}
+  if(!fo2 && !v0A && v0C && !OM2)  { h->Fill(13);}
+  if(!fo2 && !v0A && !v0C && OM2)  { h->Fill(14);}
 
 }
index 532dda7..244dfab 100644 (file)
@@ -24,50 +24,57 @@ TH1D * gHistoCompoments[kHistoFitCompoments];
 void LoadLibs();
 void LoadData(TString dataFolder, TString correctionFolder);
 void SetStyle();
-Double_t CheckSecondaries();
+void CheckSecondaries(Double_t & fracWeak, Double_t &fracMaterial);
+void CheckVz(); 
 void ShowAcceptanceInVzSlices() ;
 TH1D * GetRatioIntegratedFractions (TH1 * hNum, TH1 * hDenum) ;
-  TH1D * GetCumulativeHisto (TH1 * h) ;
+TH1D * GetCumulativeHisto (TH1 * h) ;
 static Double_t HistoSum(const double * x, const double* p);
 TF1 * GetFunctionHistoSum() ;
 TF1 * GetMTExp(Float_t norm=68, Float_t t=25) ;
 TF1 * GetHagedorn(Float_t norm=68, Float_t pt0=25, Float_t n=13) ;
 TF1 * GetLevy(Double_t temp=0.1, Double_t n=7, Double_t norm=10, const char * name="fLevy") ;
+void PrintCanvas(TCanvas* c,const TString formats) ;
 
+// global switches
+Bool_t doPrint=kFALSE;// disable PrintCanvas
 
 
-
-void correct(TString dataFolder = "output/LHC09d_000104892_p4/", TString correctionFolder = "output/LHC10a8_104867/") {
+void correct(TString dataFolder = "./output/LHC10g2d_130844_V0M_bin_10/", TString correctionFolder = "./output/LHC10g2a_130844_V0M_bin_10/") {
 
   // Load stuff and set some styles
   LoadLibs();
   LoadData(dataFolder,correctionFolder);
   SetStyle();
-  // ShowAcceptanceInVzSlices();
-  // return;
+  ShowAcceptanceInVzSlices();
+  return;
 
   // TODO add some cool printout for cuts and centrality selection
   
-  Double_t fractionWeak = CheckSecondaries();
-  cout << "Rescaling weak correction: " << fractionWeak << endl;
+  CheckVz();
+
+  Double_t fractionWeak = 1, fractionMaterial=1; 
+  //  CheckSecondaries(fractionWeak, fractionMaterial);
+  cout << "Rescaling secondaries correction, weak: " << fractionWeak << ", material: " << fractionMaterial <<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);
+  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");  
-
+  cdata->SetLogy();
   hDataPt->Draw();
   //  hMCPtRec->Draw("same");
 
   TCanvas * cMC = new TCanvas ("cMC", "Monte Carlo");
+  cMC->SetLogy();
   cMC->cd();
   hMCPtGen ->Draw();
   hMCPtRec ->Draw("same");
@@ -90,6 +97,7 @@ void correct(TString dataFolder = "output/LHC09d_000104892_p4/", TString correct
 
   TH1D * hCorSeM = (TH1D*) hMCPtSeM->Clone("hEffPt");
   hCorSeM->Divide(hMCPtSeM,hMCPtRec,1,1,"B");
+  hCorSeM->Scale(fractionMaterial);// rescale material correction
   hCorSeM->Multiply(hDataPt);
 
   TH1D * hCorSeW = (TH1D*) hMCPtSeW->Clone("hEffPt");
@@ -125,7 +133,7 @@ void correct(TString dataFolder = "output/LHC09d_000104892_p4/", TString correct
   hDataGen->Draw("same");
 }
 
-Double_t CheckSecondaries() {
+void CheckSecondaries(Double_t &fracWeak, Double_t &fracMaterial) {
   // Returns the fraction you need to rescale the secondaries from weak decays for.
 
   // Some shorthands
@@ -149,35 +157,54 @@ Double_t CheckSecondaries() {
   TCanvas * c1 = new TCanvas("cDCAFit", "Fit to the DCA distributions");  
   c1->SetLogy();
   // Draw all
-   hDataDCA->Draw();
-   hMCDCARec ->Draw("same");
-  hMCDCAPri ->Draw("same");
-  hMCDCASW ->Draw("same");
-  hMCDCASM ->Draw("same");
-  hMCDCAFak ->Draw("same");
-   return 1;
+  //  hDataDCA->Draw();
+  //  hMCDCARec ->Draw("same");
+  // hMCDCAPri ->Draw("same");
+  // hMCDCASW ->Draw("same");
+  // hMCDCASM ->Draw("same");
+  // hMCDCAFak ->Draw("same");
+  // return;
   
+  // Fit the DCA distribution. Uses a TF1 made by summing histograms
   TH1D * hMCPrimSMFak = (TH1D*) hMCDCAPri->Clone("hMCPrimSMFak");
   //  hMCPrimSMFak->Add(hMCDCASM);
   hMCPrimSMFak->Add(hMCDCAFak);
 
+  // Set the components which are used in HistoSum, the static
+  // function for GetFunctionHistoSum
   gHistoCompoments[0] = (TH1D*) hMCPrimSMFak->Clone();
   gHistoCompoments[1] = (TH1D*) hMCDCASW->Clone();
   gHistoCompoments[2] = (TH1D*) hMCDCASM->Clone();
   TF1 * fHistos = GetFunctionHistoSum();
-  fHistos->SetParameters(1,1);
+  fHistos->SetParameters(1,1,1);
+  fHistos->SetLineColor(kRed);
+  // Fit!
   hDataDCA->Fit(fHistos,"","",0,200);
+  // Rescale the components and draw to see how it looks like
   hMCPrimSMFak->Scale(fHistos->GetParameter(0));
   hMCDCASW    ->Scale(fHistos->GetParameter(1));
   hMCDCASM    ->Scale(fHistos->GetParameter(2));
   hMCPrimSMFak->Draw("same");
   hMCDCASW    ->Draw("same");
   hMCDCASM    ->Draw("same");
-  return fHistos->GetParameter(1)/fHistos->GetParameter(0);
+  // compute scaling factors
+  fracWeak     = hMCDCASW->Integral()/(hMCPrimSMFak->Integral()+hMCDCASW->Integral()+hMCDCASM->Integral());
+  fracMaterial = hMCDCASM->Integral()/(hMCPrimSMFak->Integral()+hMCDCASW->Integral()+hMCDCASM->Integral());
 
 
 }
 
+void CheckVz() {
+  // compares the Vz distribution in data and in MC
+  TCanvas * c1 = new TCanvas("cVz", "Vertex Z distribution");
+  c1->cd();
+  TH1D * hData  = hManData->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoRec       );
+  TH1D * hCorr  = hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoRec       );
+  hCorr->Draw("");
+  hData->Draw("same");
+
+}
+
 void LoadLibs() {
 
   gSystem->Load("libVMC");
@@ -297,8 +324,10 @@ void LoadData(TString dataFolder, TString correctionFolder){
   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));
+    // hManData->ScaleHistos(hEvStatData->GetBinContent(AliPhysicsSelection::kStatAccepted,irowGoodTrigger));
+    // hManCorr->ScaleHistos(hEvStatCorr->GetBinContent(AliPhysicsSelection::kStatAccepted,irowGoodTrigger));
+    hManData->ScaleHistos(hManData->GetHistoStats()->GetBinContent(2));
+    hManCorr->ScaleHistos(hManCorr->GetHistoStats()->GetBinContent(2));
   } else {
     cout << "WARNING!!! ARBITRARY SCALING" << endl;
     hManData->ScaleHistos(1000);
@@ -384,22 +413,32 @@ TH1D * GetRatioIntegratedFractions (TH1 * hNum, TH1 * hDenum) {
 
 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);
+  for(Int_t ivz = -10; ivz < -6; ivz+=2){
+    ivz=0;//FIXME  
+    Bool_t first = kTRUE;
+    TH1D * hMCPtPri  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim ,   -0.5,0.5,ivz,ivz+2);
+    TH1D * hMCPtGen  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen,        -0.5,0.5,ivz,ivz+2);
+    TH1D * hMCPtPriD  = hManData->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim ,   -0.5,0.5,ivz,ivz+2);
+    TH1D * hMCPtGenD  = hManData->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen,        -0.5,0.5,ivz,ivz+2);
     //    hEff= hMCPtGen;
-    TH1D * hEff = (TH1D*)hMCPtPri->Clone(Form("hEff_vz_%d_%d",ivz,ivz+4));
+    TH1D * hEff = (TH1D*)hMCPtPri->Clone(Form("hEff_vz_%d_%d",ivz,ivz+2));
     hEff->Divide(hMCPtPri,hMCPtGen,1,1,"B");
+    TH1D * hEffD = (TH1D*)hMCPtPriD->Clone(Form("hEffD_vz_%d_%d",ivz,ivz+2));
+    hEffD->Divide(hMCPtPriD,hMCPtGenD,1,1,"B");
+    hEffD->SetLineColor(kRed);
     cout << "ivz " << ivz << endl;
-    if(ivz < -9) {
+    if(first) {
+      first=kFALSE;
       cout << "First" << endl;
       hEff->Draw();
+      hEffD->Draw("same");
       // hMCPtGen->Draw();
       // hMCPtPri->Draw("same");
     }
     else {
       cout << "Same" << endl;
       hEff->Draw("same");
+      hEffD->Draw("same");
       // hMCPtGen->Draw("");
       // hMCPtPri->Draw("same");
     }
@@ -408,8 +447,8 @@ void ShowAcceptanceInVzSlices() {
   }
   
   //hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim )->Draw();
-   hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoRec )->Draw("");
-  hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoGen     )->Draw("same");      
+  //  hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoRec )->Draw("");
+  // hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoGen     )->Draw("same");      
 
 
 }
@@ -440,3 +479,21 @@ static Double_t HistoSum(const double * x, const double* p){
   return value;
 
 }
+
+void PrintCanvas(TCanvas* c,const TString formats) {
+  // print a canvas in every of the given comma-separated formats
+  // ensure the canvas is updated
+  if(!doPrint) return;
+  c->Update();
+  gSystem->ProcessEvents();
+  TObjArray * arr = formats.Tokenize(",");
+  TIterator * iter = arr->MakeIterator();
+  TObjString * element = 0;
+  TString name  =c ->GetName();
+  name.ReplaceAll(" ","_");
+  name.ReplaceAll("+","Plus");
+  name.ReplaceAll("-","");
+  while ((element = (TObjString*) iter->Next())) {
+    c->Print(name+ "."+element->GetString());
+  }
+}
index 97535da..906eaea 100644 (file)
@@ -1,12 +1,15 @@
 // TODO:
 // 1. Check cuts for 2010 (Jochen?)
 // 2. Run with many centrality bins at once
+#include <string.h>
 
-enum { kMyRunModeLocal = 0, kMyRunModeCAF};
+enum { kMyRunModeLocal = 0, kMyRunModeCAF, kMyRunModeGRID};
+
+TList * listToLoad = new TList();
 
 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, Int_t centrBin = 0, const char * centrEstimator = "VOM", const char* option = "",Int_t workers = -1)
+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, Int_t centrBin = 0, const char * centrEstimator = "VOM", const char* option = "",TString customSuffix = "", Int_t workers = -1)
 {
   // runMode:
   //
@@ -33,6 +36,21 @@ void run(Char_t* data, Long64_t nev = -1, Long64_t offset = 0, Bool_t debug = kF
     mgr->SetMCtruthEventHandler(handler);
   }
 
+
+  // If we are running on grid, we need the alien handler
+  if (runMode == kMyRunModeGRID) {
+    // Create and configure the alien handler plugin
+    gROOT->LoadMacro("CreateAlienHandler.C");
+    AliAnalysisGrid *alienHandler = CreateAlienHandler(data, listToLoad, "test", isMC);  
+    if (!alienHandler) {
+      cout << "Cannot create alien handler" << endl;    
+      exit(1);
+    }
+    mgr->SetGridHandler(alienHandler);  
+  }
+
+
+
   // physics selection
   gROOT->ProcessLine(".L $ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
   physicsSelectionTask = AddTaskPhysicsSelection(isMC);
@@ -58,7 +76,7 @@ void run(Char_t* data, Long64_t nev = -1, Long64_t offset = 0, Bool_t debug = kF
     doSave = kTRUE;
   }
 
-  AliESDtrackCuts * cuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(kTRUE);
+  AliESDtrackCuts * cuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
   TString pathsuffix = "";
   // cuts->SetPtRange(0.15,0.2);// FIXME pt cut
   // const char * pathsuffix = "_pt_015_020_nofakes";
@@ -88,6 +106,11 @@ void run(Char_t* data, Long64_t nev = -1, Long64_t offset = 0, Bool_t debug = kF
   gROOT->ProcessLine(".L $ALICE_ROOT/PWG0/multPbPb/AddTaskMultPbPbTracks.C");
   AliAnalysisTaskMultPbTracks * task = AddTaskMultPbPbTracks("multPbPbtracks.root", cuts); // kTRUE enables DCA cut
   task->SetIsMC(useMCKinematics);
+  if(useMCKinematics) task->GetHistoManager()->SetSuffix("MC");
+  if(customSuffix!=""){
+    cout << "Setting custom suffix: " << customSuffix << endl;    
+    task->GetHistoManager()->SetSuffix(customSuffix);
+  }
   task->SetCentralityBin(centrBin);
   task->SetCentralityEstimator(centrEstimator);
   
@@ -103,10 +126,13 @@ void run(Char_t* data, Long64_t nev = -1, Long64_t offset = 0, Bool_t debug = kF
     mgr->StartAnalysis("local",chain,nev);
   } else if (runMode == kMyRunModeCAF) {
     mgr->StartAnalysis("proof",TString(data)+"#esdTree",nev);
+  } else if (runMode == kMyRunModeGRID) {
+    mgr->StartAnalysis("grid");
   } else {
     cout << "ERROR: unknown run mode" << endl;        
   }
 
+  pathsuffix = pathsuffix + "_" + centrEstimator + "_bin_"+long(centrBin);
   if (doSave) MoveOutput(data, pathsuffix.Data());
 
   
@@ -159,6 +185,14 @@ TChain * GetAnalysisChain(const char * incollection){
 
 
 void InitAndLoadLibs(Int_t runMode=kMyRunModeLocal, Int_t workers=0,Bool_t debug=0) {
+  // Loads libs and par files + custom task and classes
+
+  // Custom stuff to be loaded
+  listToLoad->Add(new TObjString("$ALICE_ROOT/ANALYSIS/AliCentralitySelectionTask.cxx+"));
+  listToLoad->Add(new TObjString("$ALICE_ROOT/PWG1/background/AliHistoListWrapper.cxx+"));
+  listToLoad->Add(new TObjString("$ALICE_ROOT/PWG0/multPbPb/AliAnalysisMultPbTrackHistoManager.cxx+"));
+  listToLoad->Add(new TObjString("$ALICE_ROOT/PWG0/multPbPb/AliAnalysisTaskMultPbTracks.cxx+"));
+
 
   if (runMode == kMyRunModeCAF)
   {
@@ -185,7 +219,7 @@ void InitAndLoadLibs(Int_t runMode=kMyRunModeLocal, Int_t workers=0,Bool_t debug
   }
   else
   {
-    cout << "Init in Local mode" << endl;
+    cout << "Init in Local or Grid mode" << endl;
 
     gSystem->Load("libVMC");
     gSystem->Load("libTree");
@@ -201,30 +235,16 @@ void InitAndLoadLibs(Int_t runMode=kMyRunModeLocal, Int_t workers=0,Bool_t debug
     //    gROOT->ProcessLine(gSystem->ExpandPathName(".include $ALICE_ROOT/PWG1/background/"));
   }
   // Load helper classes
-  // TODO: replace this by a list of TOBJStrings
-  TString taskName("$ALICE_ROOT/PWG0/multPbPb/AliAnalysisTaskMultPbTracks.cxx+");
-  TString histoManName("$ALICE_ROOT/PWG0/multPbPb/AliAnalysisMultPbTrackHistoManager.cxx+");
-  TString listName("$ALICE_ROOT/PWG1/background/AliHistoListWrapper.cxx+");
-
-  gSystem->ExpandPathName(taskName);
-  gSystem->ExpandPathName(histoManName);
-  gSystem->ExpandPathName(listName);
-
-
-
-  // Create, add task
-  if (runMode == kMyRunModeCAF) {
-    gProof->Load(listName+(debug?"+g":""));   
-    gProof->Load(histoManName+(debug?"+g":""));
-    gProof->Load(taskName+(debug?"+g":""));
-    gProof->Load("$ALICE_ROOT/ANALYSIS/AliCentralitySelectionTask.cxx++g");      
-  } else {
-    gROOT->LoadMacro(listName+(debug?"+g":""));   
-    gROOT->LoadMacro(histoManName+(debug?"+g":""));
-    gROOT->LoadMacro(taskName+(debug?"+g":""));    
-    gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/AliCentralitySelectionTask.cxx++g");
-
+  TIterator * iter = listToLoad->MakeIterator();
+  TObjString * name = 0;
+  while (name = (TObjString *)iter->Next()) {
+    gSystem->ExpandPathName(name->String());
+    cout << name->String().Data();
+    if (runMode == kMyRunModeCAF) {
+      gProof->Load(name->String()+(debug?"+g":""));   
+    } else {
+      gROOT->LoadMacro(name->String()+(debug?"+g":""));   
+    }
   }
 
-
 }
index 94d78de..004fef6 100755 (executable)
@@ -16,11 +16,12 @@ analysismode=9; #SPD + field on
 centrBin=0
 centrEstimator="V0M"
 runTriggerStudy=no
+customSuffix=""
 
 give_help() {
 
 cat <<ENDOFGUIDE
-This scripts runs the mupliplicity analysis and the trigger study class
+This scripts runs the mupliplicity analysis and the trigger study task
 
 Available options:
  Mode control, at least one of the following options should be used
@@ -28,6 +29,7 @@ Available options:
                                Modes [default=$runmode]:
                                   0 local
                                   1 caf    
+                                  2 grid    
   -c                           Run the correction
   -s                           Run the trigger study task (by default it runs the multiplicity analysis)
  Proof settings
@@ -35,6 +37,10 @@ Available options:
   -n <nev>                     Number of events to be analized 
  Misc
   -d <dataset>                 Dataset or data collection (according to run mode) [default=$dataset]
+                                - local mode: a single ESD file, an xml collection of files on 
+                                  grid or a text file with a ESD per line
+                                - caf mode: a dataset
+                                - grid mode: a directory on alien
  Options specific to the multiplicity analysis
   -b <bin>                     Set centrality bin [default=$centrBin]
   -e <estimator>               Set centrality estimator [default=$centrEstimator]
@@ -56,18 +62,22 @@ Available options:
                                 * == can be used in trigger studies task
   -t <option>                  Command line option for root [defaul=$ropt]
   -m                           Use this to run on Monte Carlo
+  -x  <suffix>                 Set a custom suffix in the histo manager        
   -g                           Debug mode
   -h                           This help
 ENDOFGUIDE
 
 }
 
-while getopts "sr:cgmd:o:w:n:e:b:" opt; do
+while getopts "x:sr:cgmd:o:w:n:e:b:t:" opt; do
   case $opt in
     r)
       run=yes
       runmode=$OPTARG
       ;;      
+    x)
+      customSuffix=$OPTARG
+      ;;      
     s)
       runTriggerStudy=yes
       ;;      
@@ -132,7 +142,7 @@ if [ "$run" = "yes" ]
        then
        root $ropt runTriggerStudy.C\(\"$dataset\",$nev,$offset,$debug,$runmode,$isMC,\"$option\",$workers\)
     else
-       root $ropt run.C\(\"$dataset\",$nev,$offset,$debug,$runmode,$isMC,$centrBin,\"$centrEstimator\",\"$option\",$workers\)
+       root $ropt run.C\(\"$dataset\",$nev,$offset,$debug,$runmode,$isMC,$centrBin,\"$centrEstimator\",\"$option\",\"$customSuffix\",$workers\)
     fi
 fi