First version of a class to compute the number of collisions corresponding to a data...
authormfloris <mfloris@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 3 Jun 2010 13:29:56 +0000 (13:29 +0000)
committermfloris <mfloris@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 3 Jun 2010 13:29:56 +0000 (13:29 +0000)
(AliCollisionNormalization) and of a corresponding analysis task, demonstrating usage
(AliCollisionNormalizationTask). Macros to run the task are also provided.

Minor modification fo AliPhysiscsSelection: the enum of stat table is now public.

ANALYSIS/ANALYSISaliceLinkDef.h
ANALYSIS/AliCollisionNormalization.cxx [new file with mode: 0644]
ANALYSIS/AliCollisionNormalization.h [new file with mode: 0644]
ANALYSIS/AliCollisionNormalizationTask.cxx [new file with mode: 0644]
ANALYSIS/AliCollisionNormalizationTask.h [new file with mode: 0644]
ANALYSIS/AliPhysicsSelection.h
ANALYSIS/libANALYSISalice.pkg
ANALYSIS/macros/runProofNormalization.C [new file with mode: 0644]
ANALYSIS/macros/runTaskNormalization.C [new file with mode: 0644]

index 098dd1e..1878e3d 100644 (file)
@@ -31,6 +31,9 @@
 #pragma link C++ class  AliAnalysisFilter+;
 #pragma link C++ class  AliAnalysisCuts+;
 
+#pragma link C++ class AliCollisionNormalization+;
+#pragma link C++ class AliCollisionNormalizationTask+;
+
 #ifdef WITHXML
 #pragma link C++ class AliTagAnalysis+;
 #pragma link C++ class AliXMLCollection+;
diff --git a/ANALYSIS/AliCollisionNormalization.cxx b/ANALYSIS/AliCollisionNormalization.cxx
new file mode 100644 (file)
index 0000000..69928f4
--- /dev/null
@@ -0,0 +1,626 @@
+#include "AliCollisionNormalization.h"
+#include "AliPhysicsSelection.h"
+#include "AliLog.h"
+#include "TFile.h"
+#include "TCanvas.h"
+#include "AliGenPythiaEventHeader.h"
+#include "AliGenDPMjetEventHeader.h"
+#include "AliGenEventHeader.h"
+#include "AliMCEvent.h"
+
+ClassImp(AliCollisionNormalization)
+
+const char * AliCollisionNormalization::fProcLabel[] = {"SD","DD","ND", "Unknown"};
+
+AliCollisionNormalization::AliCollisionNormalization() :
+  TObject(),
+  fNbinsVz(0),
+  fMinVz(0),
+  fMaxVz(0),
+  fZRange(9.99),
+  fIsMC(0),
+  fReferenceXS(0),
+  fVerbose(0),
+  fHistVzData       (0),
+  fHistProcTypes    (0),
+  fHistStatBin0     (0)
+{
+
+  // ctor
+
+  fNbinsVz =  30;
+  fMinVz   = -15;
+  fMaxVz   = +15;
+
+  for(Int_t iproc = 0; iproc < kNProcs; iproc++){
+      fHistVzMCGen[iproc] = 0;
+      fHistVzMCRec[iproc] = 0;
+      fHistVzMCTrg[iproc] = 0;
+  }
+  
+  
+  BookAllHistos();
+}
+
+AliCollisionNormalization::AliCollisionNormalization(Int_t nbinz, Float_t minz, Float_t maxz):
+  TObject(),
+  fNbinsVz(0),
+  fMinVz(0),
+  fMaxVz(0),
+  fZRange(9.99),
+  fIsMC(0),
+  fReferenceXS(0),
+  fVerbose(0),
+  fHistVzData       (0),
+  fHistProcTypes    (0),
+  fHistStatBin0     (0)
+{
+
+  // ctor, allows setting binning
+  fNbinsVz = nbinz;
+  fMinVz   = minz;
+  fMaxVz   = maxz;
+
+  for(Int_t iproc = 0; iproc < kNProcs; iproc++){
+      fHistVzMCGen[iproc] = 0;
+      fHistVzMCRec[iproc] = 0;
+      fHistVzMCTrg[iproc] = 0;
+  }
+
+  BookAllHistos();
+}
+
+AliCollisionNormalization::AliCollisionNormalization(const char * dataFile, const char * dataListName, 
+                                                    const char * mcFile,   const char * mcListName,
+                                                    const char * eventStatFile) :
+  TObject(),
+  fNbinsVz(0),
+  fMinVz(0),
+  fMaxVz(0),
+  fZRange(9.99),
+  fIsMC(0),
+  fReferenceXS(0),
+  fVerbose(0),
+  fHistVzData       (0),
+  fHistProcTypes    (0),
+  fHistStatBin0     (0)
+{
+
+  // ctor, loads histograms from file
+  for(Int_t iproc = 0; iproc < kNProcs; iproc++){
+      fHistVzMCGen[iproc] = 0;
+      fHistVzMCRec[iproc] = 0;
+      fHistVzMCTrg[iproc] = 0;
+  }
+
+
+  TFile * fdata = new TFile (dataFile);
+  TFile * fmc   = new TFile (mcFile  );
+  
+  TList * ldata = (TList*) fdata->Get(dataListName);
+  TList * lmc   = (TList*) fmc  ->Get(mcListName  );
+
+  AliCollisionNormalization * cndata = (AliCollisionNormalization*) ldata->FindObject("AliCollisionNormalization");
+  AliCollisionNormalization * cnmc   = (AliCollisionNormalization*) lmc  ->FindObject("AliCollisionNormalization");
+
+
+  // Assign or book all histos
+  for(Int_t iproc = 0; iproc < kNProcs; iproc++){
+    fHistVzMCGen[iproc]= cnmc->GetVzMCGen(iproc);
+    fHistVzMCRec[iproc]= cnmc->GetVzMCRec(iproc);
+    fHistVzMCTrg[iproc]= cnmc->GetVzMCTrg(iproc);
+  }
+  fHistVzData        = cndata->GetVzData();
+  fHistProcTypes     = cnmc->GetHistProcTypes();
+    
+  TFile * fstat = new TFile(eventStatFile);
+  fHistStatBin0      =  (TH1F*) fstat->Get("fHistStatistics_Bin0");
+  
+}
+
+
+AliCollisionNormalization::~AliCollisionNormalization(){
+
+  // dtor
+  for(Int_t iproc = 0; iproc < kNProcs; iproc++){
+    if(fHistVzMCGen[iproc]) { delete fHistVzMCGen[iproc]      ; fHistVzMCGen[iproc]      =0;}
+    if(fHistVzMCRec[iproc]) { delete fHistVzMCRec[iproc]      ; fHistVzMCRec[iproc]      =0;}
+    if(fHistVzMCTrg[iproc]) { delete fHistVzMCTrg[iproc]      ; fHistVzMCTrg[iproc]      =0;}    
+  }
+  
+  if(fHistVzData       ) { delete fHistVzData       ; fHistVzData       =0;}
+  if(fHistStatBin0     ) { delete fHistStatBin0     ; fHistStatBin0     =0;}
+  if(fHistProcTypes    ) { delete fHistProcTypes    ; fHistProcTypes    =0;}
+
+}
+  
+void AliCollisionNormalization::BookAllHistos(){
+  // books all histos
+  // Book histos of vz distributions vs multiplicity
+  // if vzOnly == kTRUE, it returns a 1d histo with vz dist only
+
+  // Do not attach histos to the current directory
+  Bool_t oldStatus = TH1::AddDirectoryStatus();
+  TH1::AddDirectory(kFALSE);
+
+  for(Int_t iproc = 0; iproc < kNProcs; iproc++){
+    fHistVzMCGen [iproc]      = (TH2F*) BookVzHisto(TString("fHistVzMCGen")+ fProcLabel[iproc]      ,"Vz distribution of generated events vs rec multiplicity    ");
+    fHistVzMCRec [iproc]      = (TH2F*) BookVzHisto(TString("fHistVzMCRec")+ fProcLabel[iproc]      ,"Vz distribution of reconstructed events vs rec multiplicity");
+    fHistVzMCTrg [iproc]      = (TH2F*) BookVzHisto(TString("fHistVzMCTrg")+ fProcLabel[iproc]      ,"Vz distribution of triggered events vs rec multiplicity    ");    
+  }
+  fHistVzData        = (TH2F*) BookVzHisto("fHistVzData"       ,"Vz distribution of triggered events vs rec multiplicity    ");
+  fHistProcTypes     = new TH1F   ("fHistProcTypes", "Number of events in the different process classes", kNProcs, -0.5 , kNProcs-0.5);
+
+  fHistProcTypes->GetXaxis()->SetBinLabel(kProcSD+1,"SD");
+  fHistProcTypes->GetXaxis()->SetBinLabel(kProcND+1,"ND");
+  fHistProcTypes->GetXaxis()->SetBinLabel(kProcDD+1,"DD");
+  fHistProcTypes->GetXaxis()->SetBinLabel(kProcUnknown+1,"Unknown");
+
+  TH1::AddDirectory(oldStatus);
+
+}
+
+TH1 * AliCollisionNormalization::BookVzHisto(const char * name , const char * title, Bool_t vzOnly) {
+
+  // Book histos of vz distributions vs multiplicity
+  // if vzOnly == kTRUE, it returns a 1d histo with vz dist only
+
+  // Do not attach histos to the current directory
+  Bool_t oldStatus = TH1::AddDirectoryStatus();
+  TH1::AddDirectory(kFALSE);
+
+  TH1 * h;
+  Double_t binLimitsVtx[] = {-30,-25,-20,-15,-10,-7,-5.5,-4,-3,-2,-1,0,1,2,3,4,5.5,7,10,15,20,25,30};
+  if (vzOnly) {
+    h = new TH1F(name,title,22,binLimitsVtx);
+  } else {
+    h = new TH2F(name,title,22,binLimitsVtx,100,Double_t(-0.5),Double_t(99.5));
+  }
+  h->SetXTitle("V_{z} (cm)");
+  h->SetYTitle("n_{trk}");
+  h->Sumw2();
+
+  TH1::AddDirectory(oldStatus);
+
+  return h;
+
+}
+
+TH2F * AliCollisionNormalization::GetVzMCGen (Int_t procType) { 
+
+  //  returns MC gen histo. If proc type < 0 sums over all proc types, reweighting the XS
+
+  if(procType>=0)  return fHistVzMCGen[procType] ;
+
+  TH2F * sum = (TH2F*) fHistVzMCGen[kProcSD]->Clone();
+  sum->Reset();
+
+  for(Int_t iproc = 0; iproc < kProcUnknown; iproc++){
+    sum->Add(fHistVzMCGen[iproc],GetProcessWeight(iproc));    
+  }
+  
+  return sum;
+}
+TH2F *   AliCollisionNormalization::GetVzMCRec       (Int_t procType) { 
+
+  //  returns MC rec histo. If proc type < 0 sums over all proc types, reweighting the XS
+
+  if(procType>=0)  return fHistVzMCRec[procType] ;
+
+  TH2F * sum = (TH2F*) fHistVzMCRec[kProcSD]->Clone();
+  sum->Reset();
+
+  for(Int_t iproc = 0; iproc < kProcUnknown; iproc++){
+    sum->Add(fHistVzMCRec[iproc],GetProcessWeight(iproc));    
+  }
+  
+  return sum;
+
+}
+
+
+TH2F *   AliCollisionNormalization::GetVzMCTrg       (Int_t procType) { 
+
+  //  returns MC trg histo. If proc type < 0 sums over all proc types, reweighting the XS
+
+  if(procType>=0)  return fHistVzMCTrg[procType] ;
+
+  TH2F * sum = (TH2F*) fHistVzMCTrg[kProcSD]->Clone();
+  sum->Reset();
+
+  for(Int_t iproc = 0; iproc < kProcUnknown; iproc++){
+    sum->Add(fHistVzMCTrg[iproc],GetProcessWeight(iproc));    
+  }
+  
+  return sum;
+
+}
+
+Double_t AliCollisionNormalization::ComputeNint() {
+
+  // Compute the number of collisions applying all corrections
+  // TODO: check error propagation
+
+  TH1 * hVzData = fHistVzData->ProjectionX("_px",2,-1,"E"); // Skip zero bin 
+  Int_t allEventsWithVertex = (Int_t) fHistVzData->Integral(0, fHistVzData->GetNbinsX()+1,
+                                                           2, fHistVzData->GetNbinsY()+1); // include under/overflow!, skip 0 bin
+
+  // Assign histos reweighted with measured XS
+  TH2F * histVzMCGen = GetVzMCGen(-1);
+  TH2F * histVzMCTrg = GetVzMCTrg(-1);
+  TH2F * histVzMCRec = GetVzMCRec(-1);
+
+  // Get or compute BG. This assumes the CINT1B suite
+  Double_t triggeredEventsWith0MultWithBG = fHistVzData->Integral(0, fHistVzData->GetNbinsX()+1,1, 1);
+  Double_t bg = 0; // This will include beam gas + accidentals
+  if (fHistStatBin0->GetNbinsY() > 4) { // FIXME: we need a better criterion to decide...
+    AliInfo("Using BG computed by Physics Selection");
+    bg = fHistStatBin0->GetBinContent(fHistStatBin0->GetNbinsX(),AliPhysicsSelection::kStatRowBG);
+    bg += fHistStatBin0->GetBinContent(fHistStatBin0->GetNbinsX(),AliPhysicsSelection::kStatRowAcc);
+  } else {
+    AliInfo("Computing BG using CINT1A/B/C/E, ignoring intensities");
+    Int_t icol = fHistStatBin0->GetNbinsX();
+    Int_t cint1B = (Int_t) fHistStatBin0->GetBinContent(icol,1);       
+    Int_t cint1A = (Int_t) fHistStatBin0->GetBinContent(icol,2);       
+    Int_t cint1C = (Int_t) fHistStatBin0->GetBinContent(icol,3);       
+    Int_t cint1E = (Int_t) fHistStatBin0->GetBinContent(icol,4);      
+    bg   = cint1A + cint1C-2*cint1E ;
+    if (cint1B != triggeredEventsWith0MultWithBG) {
+      AliWarning(Form("Events in bin0 from physics selection and local counter not consistent: %d - %d", cint1B, triggeredEventsWith0MultWithBG));
+    }
+  }
+
+  Double_t triggeredEventsWith0Mult = triggeredEventsWith0MultWithBG - bg; 
+  if(fVerbose > 0) AliInfo(Form("Measured events with vertex: %d",allEventsWithVertex));
+  Double_t bin0 = fHistVzData->Integral(0, fHistVzData->GetNbinsX()+1,
+                                       1, 1);
+  if(fVerbose > 0) AliInfo(Form("Zero Bin, Meas: %2.2f, BG: %2.2f, Meas - BG: %2.2f", bin0, bg, triggeredEventsWith0Mult)); 
+
+  // This pointers are here so that I use the same names used in Jan Fiete's code (allows for faster/easier comparison)
+
+  TH2* eTrig =    histVzMCTrg; 
+  TH2* eTrigVtx = histVzMCRec; 
+  TH1* eTrigVtx_projx = eTrigVtx->ProjectionX("eTrigVtx_projx", 2, eTrigVtx->GetNbinsX()+1);
+
+  // compute trigger and vertex efficiency
+  TH2 * effVtxTrig = (TH2*) histVzMCRec->Clone("effVtxTrig");
+  effVtxTrig->Reset();
+  effVtxTrig->Divide(histVzMCRec,histVzMCGen,1,1,"B");
+  // Apply correction to data to get corrected events
+  TH2 * correctedEvents = (TH2*) fHistVzData->Clone("correctedEvents");
+  correctedEvents->Divide(effVtxTrig);
+
+  //  TH1 * correctedEvents = fHistVzData->ProjectionX("eTrigVtx_projx", 2, eTrigVtx->GetNbinsX()+1); 
+
+  // loop over vertex bins
+  for (Int_t i = 1; i <= fHistVzData->GetNbinsX(); i++)
+    {
+      Double_t alpha = (Double_t) hVzData->GetBinContent(i) / allEventsWithVertex;
+      Double_t events = alpha * triggeredEventsWith0Mult;
+      
+      if (histVzMCRec->GetBinContent(i,1) == 0)
+        continue;
+
+      Double_t fZ = eTrigVtx_projx->Integral(0, eTrigVtx_projx->GetNbinsX()+1) / eTrigVtx_projx->GetBinContent(i) *
+        eTrig->GetBinContent(i, 1) / eTrig->Integral(0, eTrig->GetNbinsX()+1, 1, 1);
+
+      events *= fZ;
+
+      // multiply with trigger correction
+      Double_t correction = histVzMCGen->GetBinContent(i,1)/histVzMCTrg->GetBinContent(i,1);
+      events *= correction;
+
+      if (fVerbose > 1) Printf("  Bin %d, alpha is %.2f%%, fZ is %.3f, number of events with 0 mult.: %.2f (MC comparison: %.2f)", i, alpha * 100., fZ, events, 
+                              histVzMCRec->GetBinContent(i,1));
+      correctedEvents->SetBinContent(i, 1, events);
+    }
+
+  
+  
+  // Integrate correctedEvents over full z range
+  Double_t allEvents       = correctedEvents->Integral(0, correctedEvents->GetNbinsX()+1,0, correctedEvents->GetNbinsY()+1);
+  // Integrate correctedEvents over needed z range
+  Double_t allEventsZrange = correctedEvents->Integral(correctedEvents->GetXaxis()->FindBin(-fZRange), correctedEvents->GetXaxis()->FindBin(fZRange),
+                                                      0, correctedEvents->GetNbinsY()+1);
+  
+  if(fVerbose > 1) AliInfo(Form("Results in |Vz| < %3.3f",fZRange));
+  if(fVerbose > 1) AliInfo(Form(" Events in Bin0: %2.2f, With > 1 track: %2.2f, All corrected: %2.2f",
+                               correctedEvents->Integral(correctedEvents->GetXaxis()->FindBin(-fZRange), correctedEvents->GetXaxis()->FindBin(fZRange),1,1),
+                               correctedEvents->Integral(correctedEvents->GetXaxis()->FindBin(-fZRange), correctedEvents->GetXaxis()->FindBin(fZRange), 
+                                                         2,correctedEvents->GetNbinsX()+1),
+                               allEventsZrange));
+  if(fVerbose > 1) {    
+    Int_t nbin = histVzMCRec->GetNbinsX();
+    AliInfo(Form("Efficiency in the zero bin: %3.3f", histVzMCRec->Integral(0,nbin+1,1,1)/histVzMCGen->Integral(0,nbin+1,1,1) ));
+  }
+  
+  AliInfo(Form("Number of collisions in full phase space: %2.2f", allEvents));
+
+  return allEvents;
+}
+
+Long64_t AliCollisionNormalization::Merge(TCollection* list)
+{
+  // Merge a list of AliCollisionNormalization objects with this
+  // (needed for PROOF).  
+  // Returns the number of merged objects (including this).
+
+  if (!list)
+    return 0;
+
+  if (list->IsEmpty())
+    return 1;
+
+  TIterator* iter = list->MakeIterator();
+  TObject* obj;
+  
+  // collections of all histograms
+  const Int_t nHists = kNProcs*3+4;
+  TList collections[nHists];
+
+  Int_t count = 0;
+  while ((obj = iter->Next())) {
+
+    AliCollisionNormalization* entry = dynamic_cast<AliCollisionNormalization*> (obj);
+    if (entry == 0) 
+      continue;
+    Int_t ihist = -1;
+
+    for(Int_t iproc = 0; iproc < kNProcs; iproc++){
+      if (entry->fHistVzMCGen[iproc]      ) collections[++ihist].Add(entry->fHistVzMCGen[iproc]      );
+      if (entry->fHistVzMCRec[iproc]      ) collections[++ihist].Add(entry->fHistVzMCRec[iproc]      );
+      if (entry->fHistVzMCTrg[iproc]      ) collections[++ihist].Add(entry->fHistVzMCTrg[iproc]      );
+    }
+    if (entry->fHistVzData       ) collections[++ihist].Add(entry->fHistVzData       );
+    if (entry->fHistProcTypes    ) collections[++ihist].Add(entry->fHistProcTypes    );
+    if (entry->fHistStatBin0     ) collections[++ihist].Add(entry->fHistStatBin0     );
+
+    count++;
+  }
+
+  Int_t ihist = -1;
+  for(Int_t iproc = 0; iproc < kNProcs; iproc++){
+    if (fHistVzMCGen[iproc]      ) fHistVzMCGen[iproc]      ->Merge(&collections[++ihist]);
+    if (fHistVzMCRec[iproc]      ) fHistVzMCRec[iproc]      ->Merge(&collections[++ihist]);
+    if (fHistVzMCTrg[iproc]      ) fHistVzMCTrg[iproc]      ->Merge(&collections[++ihist]);
+
+  }  
+  if (fHistVzData       ) fHistVzData       ->Merge(&collections[++ihist]);
+  if (fHistProcTypes    ) fHistProcTypes    ->Merge(&collections[++ihist]);
+  if (fHistStatBin0     ) fHistStatBin0     ->Merge(&collections[++ihist]);
+    
+  
+  delete iter;
+
+  return count+1;
+}
+
+void AliCollisionNormalization::FillVzMCGen(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt) {
+
+  // Fill MC gen histo and the process types statistics
+
+  Int_t evtype = GetProcessType(mcEvt);
+  // When I fill gen histos, I also fill statistics of process types (used to detemine ratios afterwards).
+  fHistProcTypes->Fill(evtype);
+  fHistVzMCGen[evtype]->Fill(vz,ntrk);
+}
+
+void AliCollisionNormalization::FillVzMCRec(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt){
+
+  // Fill MC rec histo
+  Int_t evtype = GetProcessType(mcEvt);
+  fHistVzMCRec[evtype]->Fill(vz,ntrk);
+
+}
+void AliCollisionNormalization::FillVzMCTrg(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt) {      
+
+  // Fill MC trg histo
+  Int_t evtype = GetProcessType(mcEvt);
+  fHistVzMCTrg[evtype]->Fill(vz,ntrk);
+}
+
+
+Int_t AliCollisionNormalization::GetProcessType(AliMCEvent * mcEvt) {
+
+  // Determine if the event was generated with pythia or phojet and return the process type
+
+  // Check if mcEvt is fine
+  if (!mcEvt) {
+    AliFatal("NULL mc event");
+  } 
+
+  // Determine if it was a pythia or phojet header, and return the correct process type
+  AliGenPythiaEventHeader * headPy  = 0;
+  AliGenDPMjetEventHeader * headPho = 0;
+  AliGenEventHeader * htmp = mcEvt->GenEventHeader();
+  if(!htmp) {
+    AliFatal("Cannot Get MC Header!!");
+  }
+  if( TString(htmp->IsA()->GetName()) == "AliGenPythiaEventHeader") {
+    headPy =  (AliGenPythiaEventHeader*) htmp;
+  } else if (TString(htmp->IsA()->GetName()) == "AliGenDPMjetEventHeader") {
+    headPho = (AliGenDPMjetEventHeader*) htmp;
+  } else {
+    AliError("Unknown header");
+  }
+
+  // Determine process type
+  if(headPy)   {
+    if(headPy->ProcessType() == 92 || headPy->ProcessType() == 93) {
+      // single difractive
+      return kProcSD;
+    } else if (headPy->ProcessType() == 94) {
+      // double diffractive
+      return kProcDD;
+    }
+    else if(headPy->ProcessType() != 92 && headPy->ProcessType() != 93 && headPy->ProcessType() != 94) {    
+      // non difractive
+      return kProcND; 
+    }
+  } else if (headPho) {
+    if(headPho->ProcessType() == 5 || headPho->ProcessType() == 6 ) {
+      // single difractive
+      return kProcSD;
+    } else if (headPho->ProcessType() == 7) { 
+      // double diffractive
+      return kProcDD;      
+    } else if(headPho->ProcessType() != 5 && headPho->ProcessType() != 6  && headPho->ProcessType() != 7 ) {
+      // non difractive
+      return kProcND; 
+    }       
+  }
+  
+
+  // no process type found?
+  AliError(Form("Unknown header: %s", htmp->IsA()->GetName()));
+  return kProcUnknown;
+}
+
+Double_t AliCollisionNormalization::GetProcessWeight(Int_t proc) {
+
+  // Return a weight to be used when combining the MC histos to
+  // compute efficiency, defined as the ratio XS measured / XS in
+  // generator
+
+  Float_t ref_SD,  ref_DD,  ref_ND,  error_SD,  error_DD,  error_ND;
+  GetRelativeFractions(fReferenceXS,ref_SD,  ref_DD,  ref_ND,  error_SD,  error_DD,  error_ND);
+
+  static Double_t total = fHistProcTypes->Integral();
+  if (fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcUnknown)) > 0) {
+    AliError("There were unknown processes!!!");
+  }
+  static Double_t SD = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcSD))/total;
+  static Double_t DD = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcDD))/total;
+  static Double_t ND = fHistProcTypes->GetBinContent(fHistProcTypes->FindBin(kProcND))/total;
+
+  if (fVerbose > 2) {
+    AliInfo(Form("Total MC evts: %f",total));
+    AliInfo(Form(" Frac SD %4.4f", SD));
+    AliInfo(Form(" Frac DD %4.4f", DD));
+    AliInfo(Form(" Frac ND %4.4f", ND));
+  }
+  
+  switch(proc) {
+  case kProcSD:
+    return ref_SD/SD;
+    break;
+  case kProcDD:
+    return ref_DD/DD;
+    break;
+  case kProcND:
+    return ref_ND/ND;
+    break;
+  default:
+    AliError("Unknown process");
+  }
+  
+  return 0;
+
+} 
+
+void AliCollisionNormalization::GetRelativeFractions(Int_t origin, Float_t& ref_SD, Float_t& ref_DD, Float_t& ref_ND, Float_t& error_SD, Float_t& error_DD, Float_t& error_ND)
+{
+  // Returns fraction of XS (SD, ND, DD) and corresponding error
+  // Stolen from Jan Fiete's drawSystematics macro
+
+  // origin: 
+  //   -1 = Pythia (test)
+  //   0 = UA5
+  //   1 = Data 1.8 TeV
+  //   2 = Tel-Aviv
+  //   3 = Durham
+  //
+
+  switch (origin)
+  {
+    case -10: // Pythia default at 7 GeV, 50% error
+      AliInfo("PYTHIA x-sections");
+      ref_SD = 0.192637; error_SD = ref_SD * 0.5;
+      ref_DD = 0.129877; error_DD = ref_DD * 0.5;
+      ref_ND = 0.677486; error_ND = 0;
+      break;
+
+    case -1: // Pythia default at 900 GeV, as test
+      AliInfo("PYTHIA x-sections");
+      ref_SD = 0.223788;
+      ref_DD = 0.123315;
+      ref_ND = 0.652897;
+      break;
+      
+    case 0: // UA5
+      AliInfo("UA5 x-sections a la first paper");
+      ref_SD = 0.153; error_SD = 0.05;
+      ref_DD = 0.080; error_DD = 0.05;
+      ref_ND = 0.767; error_ND = 0;
+      break;
+      
+    case 10: // UA5
+      AliInfo("UA5 x-sections hadron level definition for Pythia"); 
+      // Fractions in Pythia with UA5 cuts selection for SD
+      // ND: 0.688662
+      // SD: 0.188588 --> this should be 15.3
+      // DD: 0.122750
+      ref_SD = 0.224 * 0.153 / 0.189; error_SD = 0.023 * 0.224 / 0.189;
+      ref_DD = 0.095;                 error_DD = 0.06; 
+      ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
+      break;
+    
+    case 11: // UA5
+      AliInfo("UA5 x-sections hadron level definition for Phojet"); 
+      // Fractions in Phojet with UA5 cuts selection for SD
+      // ND: 0.783573
+      // SD: 0.151601 --> this should be 15.3
+      // DD: 0.064827
+      ref_SD = 0.191 * 0.153 / 0.152; error_SD = 0.023 * 0.191 / 0.152;
+      ref_DD = 0.095;                 error_DD = 0.06; 
+      ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
+      break;
+      
+    case 20: // E710, 1.8 TeV
+      AliInfo("E710 x-sections hadron level definition for Pythia");
+      // ND: 0.705709
+      // SD: 0.166590 --> this should be 15.9
+      // DD: 0.127701
+      ref_SD = 0.217 * 0.159 / 0.167; error_SD = 0.024 * 0.217 / 0.167;
+      ref_DD = 0.075 * 1.43;          error_DD = 0.02 * 1.43; 
+      ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
+      break;
+    
+    case 21: // E710, 1.8 TeV
+      AliInfo("E710 x-sections hadron level definition for Phojet"); 
+      // ND: 0.817462
+      // SD: 0.125506 --> this should be 15.9
+      // DD: 0.057032
+      ref_SD = 0.161 * 0.159 / 0.126; error_SD = 0.024 * 0.161 / 0.126;
+      ref_DD = 0.075 * 1.43;         error_DD = 0.02 * 1.43;
+      ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
+      break;
+    
+    case 1: // data 1.8 TeV
+      AliInfo("??? x-sections");
+      ref_SD = 0.152;
+      ref_DD = 0.092;
+      ref_ND = 1 - ref_SD - ref_DD;
+      break;
+      
+    case 2: // tel-aviv model
+      AliInfo("Tel-aviv model x-sections");
+      ref_SD = 0.171;
+      ref_DD = 0.094;
+      ref_ND = 1 - ref_SD - ref_DD;
+      break;
+    
+    case 3: // durham model
+      AliInfo("Durham model x-sections");
+      ref_SD = 0.190;
+      ref_DD = 0.125;
+      ref_ND = 1 - ref_SD - ref_DD;
+      break;
+    
+    default:
+      AliFatal(Form("Unknown origin %d", origin));
+  }
+}
+
diff --git a/ANALYSIS/AliCollisionNormalization.h b/ANALYSIS/AliCollisionNormalization.h
new file mode 100644 (file)
index 0000000..a53b643
--- /dev/null
@@ -0,0 +1,105 @@
+
+#ifndef ALICOLLISIONNORMALIZATION_H
+#define ALICOLLISIONNORMALIZATION_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//-------------------------------------------------------------------------
+//                      Implementation of   Class AliCollisionNormalization
+//
+//  This class is used to store the vertex ditributions in the data
+//  and in Monte Carlo, needed to compute the real number of
+//  collisions a given sample is corresponding to.
+//  The strategy matches what described in CERN-THESIS-2009-033 p 119.
+//
+//    Author:     Michele Floris, CERN
+//-------------------------------------------------------------------------
+
+#include "TH2F.h"
+
+class TH1F;
+class TH1I;
+class AliMCEvent;
+
+class AliCollisionNormalization : public TObject
+
+{
+
+
+public:
+  enum { kNevBin0, kNevCollisions, kNevNbin };
+  typedef enum { kProcSD, kProcDD, kProcND, kProcUnknown, kNProcs } ProcType_t; 
+  AliCollisionNormalization();
+  AliCollisionNormalization(Int_t nbinz, Float_t minz, Float_t maxz);
+  AliCollisionNormalization(const char * dataFile, const char * dataListName, 
+                           const char * mcFile,   const char * mcListName,
+                           const char * eventStatFile);
+
+  ~AliCollisionNormalization();
+  
+  void SetMC(Bool_t flag = kTRUE) { fIsMC = flag;}
+
+  void BookAllHistos();
+  TH1 * BookVzHisto(const char * name , const char * title, Bool_t vzOnly=kFALSE);
+
+  void FillVzMCGen(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt);      
+  void FillVzMCRec(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt);      
+  void FillVzMCTrg(Float_t vz, Int_t ntrk, AliMCEvent * mcEvt);      
+  void FillVzData (Float_t vz, Int_t ntrk)      {fHistVzData       ->Fill(vz,ntrk);}
+
+  TH2F *   GetVzMCGen       (Int_t procType) ;
+  TH2F *   GetVzMCRec       (Int_t procType) ;
+  TH2F *   GetVzMCTrg       (Int_t procType) ;
+  TH2F *   GetVzData        () { return fHistVzData       ; }
+  TH1F *   GetStatBin0      () { return fHistStatBin0     ; }
+  TH1F *   GetHistProcTypes () { return fHistProcTypes    ; }
+   
+
+  Int_t GetProcessType(AliMCEvent * mcEvt) ;
+  Double_t GetProcessWeight(Int_t proctype);
+
+  void SetReferencsXS(Int_t ref) { fReferenceXS = ref;}
+  
+  Double_t ComputeNint();
+  void SetZRange(Float_t zrange) { fZRange = zrange ;}
+
+  void SetReferenceXS(Int_t ref) { fReferenceXS = ref ;}
+  void GetRelativeFractions(Int_t origin, Float_t& ref_SD, Float_t& ref_DD, Float_t& ref_ND, Float_t& error_SD, Float_t& error_DD, Float_t& error_ND);
+
+  void SetVerbose(Int_t lev) { fVerbose = lev ;}
+
+  Long64_t Merge(TCollection* list);
+
+protected:
+
+  Int_t   fNbinsVz; // number of z bins in the vz histo
+  Float_t fMinVz  ; // lowest Z
+  Float_t fMaxVz  ; // highest Z
+
+  Float_t fZRange; // max |Z| vertex to be considered
+
+  Bool_t fIsMC; // True if processing MC
+  
+  Int_t fReferenceXS;                // index of reference cross section to be used to rescale process types in the calculation of the efficiency
+
+  Int_t fVerbose;                    // Determines the ammount of printout
+
+  TH2F * fHistVzMCGen[kNProcs]    ;    // Vz distribution of generated events vs rec multiplicity
+  TH2F * fHistVzMCRec[kNProcs]    ;    // Vz distribution of reconstructed events vs rec multiplicity
+  TH2F * fHistVzMCTrg[kNProcs]    ;    // Vz distribution of triggered events vs rec multiplicity
+  TH2F * fHistVzData              ;    // Vz distribution of triggered events vs rec multiplicity    
+  TH1F * fHistProcTypes           ;    // Number of evts for different Process types 
+
+  TH1F * fHistStatBin0     ; // event stat histogram, created by physiscs selection; used in ComputeNint;
+
+  static const char * fProcLabel[] ; // labels of the different process types
+  
+  ClassDef(AliCollisionNormalization, 1);
+    
+private:
+  AliCollisionNormalization(const AliCollisionNormalization&);
+  AliCollisionNormalization& operator=(const AliCollisionNormalization&);
+};
+
+#endif
diff --git a/ANALYSIS/AliCollisionNormalizationTask.cxx b/ANALYSIS/AliCollisionNormalizationTask.cxx
new file mode 100644 (file)
index 0000000..8344f9d
--- /dev/null
@@ -0,0 +1,220 @@
+/* $Id$ */
+
+
+// Simple task to test the collision normalization. Can be called for
+// both MC and data and shows how to fill the collision normalization
+// class
+//
+// Author: Michele Floris
+//         CERN
+
+
+#include "AliCollisionNormalizationTask.h"
+
+#include <TFile.h>
+#include <TH1F.h>
+#include <TH2F.h>
+
+#include <AliLog.h>
+#include <AliESDEvent.h>
+#include <AliHeader.h>
+
+#include "AliCollisionNormalization.h"
+#include "AliAnalysisManager.h"
+#include "AliInputEventHandler.h"
+
+//#include "AliBackgroundSelection.h"
+#include "AliMultiplicity.h"
+#include "AliMCEvent.h"
+
+ClassImp(AliCollisionNormalizationTask)
+
+AliCollisionNormalizationTask::AliCollisionNormalizationTask() :
+  AliAnalysisTaskSE("AliCollisionNormalizationTask"),
+  fOutput(0),
+  fIsMC(0),
+  fCollisionNormalization(0)
+{
+  //
+  // Default event handler
+  //
+
+  // Define input and output slots here
+  DefineOutput(1, TList::Class());
+  
+}
+
+AliCollisionNormalizationTask::AliCollisionNormalizationTask(const char* name) :
+  AliAnalysisTaskSE(name),
+  fOutput(0),
+  fIsMC(0),
+  fCollisionNormalization(new AliCollisionNormalization())
+{
+  //
+  // Constructor. Initialization of pointers
+  //
+
+  // Define input and output slots here
+  DefineOutput(1, TList::Class());
+  
+  //  AliLog::SetClassDebugLevel("AliCollisionNormalizationTask", AliLog::kWarning);
+}
+
+AliCollisionNormalizationTask::~AliCollisionNormalizationTask()
+{
+  //
+  // Destructor
+  //
+
+  // histograms are in the output list and deleted when the output
+  // list is deleted by the TSelector dtor
+
+  if (fOutput) {
+    delete fOutput;
+    fOutput = 0;
+  }
+}
+
+void AliCollisionNormalizationTask::UserCreateOutputObjects()
+{
+  // create result objects and add to output list
+
+  Printf("AliCollisionNormalizationTask::CreateOutputObjects");
+
+  fOutput = new TList;
+  fOutput->SetOwner();
+  
+  if (!fCollisionNormalization)
+    fCollisionNormalization = new AliCollisionNormalization;
+  
+  fOutput->Add(fCollisionNormalization);
+//   fOutput->Add(fCollisionNormalization->GetVzCorrZeroBin ());
+//   fOutput->Add(fCollisionNormalization->GetVzMCGen       ());
+//   fOutput->Add(fCollisionNormalization->GetVzMCRec       ());
+//   fOutput->Add(fCollisionNormalization->GetVzMCTrg       ());
+//   fOutput->Add(fCollisionNormalization->GetVzData        ());
+//   fOutput->Add(fCollisionNormalization->GetNEvents       ());
+//   fOutput->Add(fCollisionNormalization->GetStatBin0      ());
+
+}
+
+void AliCollisionNormalizationTask::UserExec(Option_t*)
+{
+  // process the event
+
+
+  PostData(1, fOutput);
+
+  // Get the ESD
+  AliESDEvent * aESD = dynamic_cast<AliESDEvent*>(fInputEvent);
+  if (strcmp(aESD->ClassName(),"AliESDEvent")) {
+    AliFatal("Not processing ESDs");
+  }
+  // Get MC event, if needed
+  AliMCEvent* mcEvent = fIsMC ? MCEvent() : 0;
+  if (!mcEvent && fIsMC){
+    AliFatal("Running on MC but no MC handler available");
+  }
+
+  // Physics selection. At least in the case of MC we cannot use
+  // yourTask->SelectCollisionCandidates();, because we also need to
+  // fill the "generated" histogram 
+  // NB never call IsEventSelected more than once per event
+  // (statistics histogram would be altered)
+
+  Bool_t isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
+
+  // Get the Multiplicity cut
+  const AliMultiplicity* mult = aESD->GetMultiplicity();
+  if (!mult){
+    AliError("Can't get mult object");
+    return;
+  }
+
+  // Check if the event is "reconstructed" according to some quality
+  // cuts. Tipically, we mean "it has a good-eonugh vertex"
+
+  Int_t ntracklet = mult->GetNumberOfTracklets();
+  const AliESDVertex * vtxESD = aESD->GetPrimaryVertexSPD();
+  if(vtxESD) {
+    // If there is a vertex from vertexer z with delta phi > 0.02 we
+    // don't consider it rec (we keep the event in bin0). If quality
+    // is good eneough we check the number of tracklets
+    // A similar selection is applied in IsEventInBinZero
+    // If the vertex is too far away, we don't consider it reconstructed either.
+    if ((vtxESD->IsFromVertexerZ() && vtxESD->GetDispersion() > 0.02) || TMath::Abs(vtxESD->GetZ()) > 15) {
+      vtxESD = 0;
+      ntracklet = 0; // Don't trust reconstructed tracklets if you don't trust clusters
+    } 
+  } 
+  
+  if (ntracklet > 0 && !vtxESD) {
+    AliError("No vertex but reconstructed tracklets?");
+  }
+
+  // assign vz. For MC we use generated vz
+  Float_t vz = 0;
+  if (!fIsMC) vz = vtxESD ? vtxESD->GetZ() : 0; // FIXME : is zv used anywhere in Gen?
+  else        vz = mcEvent->GetPrimaryVertex()->GetZ();
+
+  if (fIsMC) {
+    // Monte Carlo:  we fill 3 histos
+    fCollisionNormalization->FillVzMCGen(vz, ntracklet, mcEvent);      
+    // If triggered == passing the physics selection
+    if (isSelected) fCollisionNormalization->FillVzMCTrg(vz, ntracklet, mcEvent);
+    // If reconstructer == good enough vertex
+    if (vtxESD) fCollisionNormalization->FillVzMCRec(vz, ntracklet, mcEvent);    
+  } else {
+    if (isSelected) {
+      // Passing the trigger
+      fCollisionNormalization->FillVzData(vz,ntracklet);
+    }
+  }
+
+}
+
+void AliCollisionNormalizationTask::Terminate(Option_t *)
+{
+  // The Terminate() function is the last function to be called during
+  // a query. It always runs on the client, it can be used to present
+  // the results graphically or save the results to file.
+
+  fOutput = dynamic_cast<TList*> (GetOutputData(1));
+  if (!fOutput)
+    Printf("ERROR: fOutput not available");
+
+
+}
+
+Bool_t AliCollisionNormalizationTask::IsEventInBinZero() {
+
+  // Returns true if an event is to be assigned to the zero bin
+  //
+  // You should have your own version of this method in the class: in
+  // general, the definition of "reconstructed" event is subjective.
+
+  Bool_t isZeroBin = kTRUE;
+  const AliESDEvent* esd= dynamic_cast<AliESDEvent*>(fInputEvent);
+  const AliMultiplicity* mult = esd->GetMultiplicity();
+  if (!mult){
+    Printf("AliAnalysisTaskBGvsTime::IsBinZero: Can't get mult object");
+    return kFALSE;
+  }
+  Int_t ntracklet = mult->GetNumberOfTracklets();
+  const AliESDVertex * vtxESD = esd->GetPrimaryVertexSPD();
+  if(vtxESD) {
+    // If there is a vertex from vertexer z with delta phi > 0.02 we
+    // don't consider it rec (we keep the event in bin0). If quality
+    // is good eneough we check the number of tracklets
+    // if the vertex is more than 15 cm away, this is autamatically bin0
+    if( TMath::Abs(vtxESD->GetZ()) <= 15 ) {
+      if (vtxESD->IsFromVertexerZ()) {
+       if (vtxESD->GetDispersion()<=0.02 ) {
+         if(ntracklet>0) isZeroBin = kFALSE;
+       }
+      } else if(ntracklet>0) isZeroBin = kFALSE; // if the event is not from Vz we chek the n of tracklets
+    } 
+  }
+  return isZeroBin;
+
+}
diff --git a/ANALYSIS/AliCollisionNormalizationTask.h b/ANALYSIS/AliCollisionNormalizationTask.h
new file mode 100644 (file)
index 0000000..c88b1ce
--- /dev/null
@@ -0,0 +1,47 @@
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+
+/* $Id$ */
+
+#ifndef ALICOLLISIONNORMALIZATIONTASK_H
+#define ALICOLLISIONNORMALIZATIONTASK_H
+
+#include "AliAnalysisTaskSE.h"
+
+class AliCollisionNormalization;
+
+class AliCollisionNormalizationTask : public AliAnalysisTaskSE {
+  public:
+    AliCollisionNormalizationTask();
+    AliCollisionNormalizationTask(const char* name);
+
+    virtual ~AliCollisionNormalizationTask();
+
+    virtual void   UserCreateOutputObjects();
+    virtual void   UserExec(Option_t*);
+    virtual void   Terminate(Option_t*);
+
+  Bool_t IsEventInBinZero(); // returns true if the event has to be put in the bin0.
+  void SetMC(Bool_t flag = kTRUE) { fIsMC = flag; }
+  
+  //    void SetOption(const char* opt) { fOption = opt; }
+    
+  //    void SetCollisionNormalization(AliCollisionNormalization* physicsSelection) { fCollisionNormalization = physicsSelection; }
+  AliCollisionNormalization* GetCollisionNormalization() const { return fCollisionNormalization; }
+
+ protected:
+  TList* fOutput;                  //! list send on output slot 1
+  //    TString fOption;                 // option string  
+  Bool_t fIsMC;
+  
+  AliCollisionNormalization* fCollisionNormalization; // collision normalization class
+
+ private:
+    AliCollisionNormalizationTask(const AliCollisionNormalizationTask&);
+    AliCollisionNormalizationTask& operator=(const AliCollisionNormalizationTask&);
+
+  ClassDef(AliCollisionNormalizationTask, 1);
+};
+
+#endif
index 02fe485..6fbc606 100644 (file)
@@ -31,6 +31,7 @@ class AliAnalysisTaskSE;
 
 class AliPhysicsSelection : public AliAnalysisCuts
 {
+public:
 
   enum {kStatTriggerClass=1,kStatHWTrig,kStatV0ABG,kStatV0CBG,kStatMB1,kStatMB1Prime,kStatFMD,kStatFO1,kStatFO2,kStatV0A,kStatV0C,kStatSSD1,kStatFO1AndV0,kStatV0,kStatAny2Hits,kStatOffline,kStatBG,kStatAccepted};
 
@@ -43,7 +44,6 @@ class AliPhysicsSelection : public AliAnalysisCuts
   enum {kStatIdxAll=0,kStatIdxBin0=1};
 
 
-public:
   typedef Bool_t (*Bin0Callback_t)(const AliESDEvent *);
 
   AliPhysicsSelection();
index 2813d02..e7af22f 100644 (file)
@@ -9,7 +9,8 @@ SRCS =  AliAnalysisTaskSE.cxx AliAnalysisTaskME.cxx \
        AliTriggerAnalysis.cxx \
         AliPhysicsSelection.cxx AliBackgroundSelection.cxx \
         AliPhysicsSelectionTask.cxx \
-        AliAnalysisFilter.cxx AliAnalysisCuts.cxx 
+        AliAnalysisFilter.cxx AliAnalysisCuts.cxx \
+       AliCollisionNormalization.cxx AliCollisionNormalizationTask.cxx
 
 
 ifeq (yes,$(CHECKALIEN))
diff --git a/ANALYSIS/macros/runProofNormalization.C b/ANALYSIS/macros/runProofNormalization.C
new file mode 100644 (file)
index 0000000..c754bc6
--- /dev/null
@@ -0,0 +1,66 @@
+
+void runProofNormalization(const char * dataset = "LHC09b12_7TeV_0.5T", TString dataSetPath ="/PWG0/jgrosseo/",const char * filename = "LHC09b12_7TeV_0.5T_norm.root", Bool_t isMC = 1,Int_t nev =123456789) {
+
+  gEnv->SetValue("XSec.GSI.DelegProxy","2");
+  TProof::Open("alicecaf");
+  
+  //  gSystem->AddIncludePath("-I${ALICE_ROOT}/include/ -I${ALICE_ROOT}/PWG0/ -I${ALICE_ROOT}/PWG0/dNdEta/");
+  gSystem->AddIncludePath("-I${ALICE_ROOT}/include/");
+  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/CORRFW");
+  gProof->EnablePackage("$ALICE_ROOT/CORRFW");
+
+  // Make the analysis manager
+  AliAnalysisManager *mgr = new AliAnalysisManager("TestManager");
+  //  mgr->SetDebugLevel(3);
+  // Add ESD handler
+  AliESDInputHandler* esdH = new AliESDInputHandler; 
+
+  mgr->SetInputEventHandler(esdH);
+       
+  if(isMC) {
+    AliMCEventHandler *mc = new AliMCEventHandler();
+    mc->SetReadTR(kFALSE);
+    mgr->SetMCtruthEventHandler(mc);
+  }
+  // assign simple task
+//   gProof->Load("AliCollisionNormalization.cxx++g");   
+//   gProof->Load("AliCollisionNormalizationTask.cxx++g");   
+  //____________________________________________//
+  // assign simple task
+  AliCollisionNormalizationTask * task = new AliCollisionNormalizationTask("TaskNormalization");
+  //  task->SetMC();
+  task->SetMC(isMC);
+  mgr->AddTask(task);
+
+
+  gROOT->LoadMacro("$(ALICE_ROOT)/ANALYSIS/macros/AddTaskPhysicsSelection.C");
+  AliPhysicsSelectionTask* physSelTask = AddTaskPhysicsSelection(isMC,1,!isMC); // Use Physics Selection. Enable computation of BG if is not MC
+  //  task->SelectCollisionCandidates(); /// This should be disabled, at least for MC: we need all the events
+  physSelTask->GetPhysicsSelection()->SetBin0Callback("TaskNormalization");
+
+  AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();  
+  mgr->ConnectInput(task,0,cinput1);
+
+
+  
+  // Attach output
+  cOutput = mgr->CreateContainer("Norm", TList::Class(), AliAnalysisManager::kOutputContainer,filename);
+  mgr->ConnectOutput(task, 1, cOutput);      
+       
+  if (!mgr->InitAnalysis()) return;
+       
+  mgr->PrintStatus();
+  mgr->StartAnalysis("proof",dataSetPath+dataset+"#esdTree",nev);
+
+}
diff --git a/ANALYSIS/macros/runTaskNormalization.C b/ANALYSIS/macros/runTaskNormalization.C
new file mode 100644 (file)
index 0000000..179e3b2
--- /dev/null
@@ -0,0 +1,80 @@
+
+void runTaskNormalization(const char * incollection,const char * filename = "LHC09b12_7TeV_0.5T_norm.root", Bool_t isMC = 1,Int_t nev =123456789) {
+
+  
+  // Load libraries
+  gSystem->Load("libANALYSIS") ;
+  gSystem->Load("libANALYSISalice") ;
+  gSystem->Load("libCORRFW") ;
+  gSystem->Load("libITSbase") ;
+  gSystem->Load("libPWG0base") ;
+
+
+  // chain
+  TChain* analysisChain = 0;
+  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();
+
+
+  // Make the analysis manager
+  AliAnalysisManager *mgr = new AliAnalysisManager("TestManager");
+  //  mgr->SetDebugLevel(3);
+  // Add ESD handler
+  AliESDInputHandler* esdH = new AliESDInputHandler; 
+
+  mgr->SetInputEventHandler(esdH);
+       
+  if(isMC) {
+    AliMCEventHandler *mc = new AliMCEventHandler();
+    mc->SetReadTR(kFALSE);
+    mgr->SetMCtruthEventHandler(mc);
+  }
+  // assign simple task
+//   gROOT->LoadMacro("AliCollisionNormalization.cxx++g");   
+//   gROOT->LoadMacro("AliCollisionNormalizationTask.cxx++g");   
+  //____________________________________________//
+  // Physics selection
+  gROOT->LoadMacro("$(ALICE_ROOT)/ANALYSIS/macros/AddTaskPhysicsSelection.C");
+  AliPhysicsSelectionTask* physSelTask = AddTaskPhysicsSelection(isMC,1,!isMC); // Use Physics Selection. Enable computation of BG if is not MC
+  //  task->SelectCollisionCandidates(); /// This should be disabled, at least for MC: we need all the events
+  physSelTask->GetPhysicsSelection()->SetBin0Callback("TaskNormalization");
+
+  // assign simple task
+  AliCollisionNormalizationTask * task = new AliCollisionNormalizationTask("TaskNormalization");
+  //  task->SetMC();
+  task->SetMC(isMC);
+  mgr->AddTask(task);
+
+
+
+  AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();  
+  mgr->ConnectInput(task,0,cinput1);
+
+
+  
+  // Attach output
+  cOutput = mgr->CreateContainer("Norm", TList::Class(), AliAnalysisManager::kOutputContainer,filename);
+  mgr->ConnectOutput(task, 1, cOutput);      
+       
+  if (!mgr->InitAnalysis()) return;
+       
+  mgr->PrintStatus();
+  mgr->StartAnalysis("local",analysisChain,nev);
+
+}