The code and scripts to do multiplicity distributions with FMD+SPD from C.Nygaard
authorhdalsgaa <hdalsgaa@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 31 Aug 2012 09:01:44 +0000 (09:01 +0000)
committerhdalsgaa <hdalsgaa@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 31 Aug 2012 09:01:44 +0000 (09:01 +0000)
15 files changed:
PWGLF/CMakelibPWGLFforward2.pkg
PWGLF/FORWARD/analysis2/AliForwardCreateResponseMatrices.cxx [new file with mode: 0644]
PWGLF/FORWARD/analysis2/AliForwardCreateResponseMatrices.h [new file with mode: 0644]
PWGLF/FORWARD/analysis2/AliForwardMultiplicityDistribution.cxx [new file with mode: 0644]
PWGLF/FORWARD/analysis2/AliForwardMultiplicityDistribution.h [new file with mode: 0644]
PWGLF/FORWARD/analysis2/scripts/multdist/AddTaskCreateResponseMatrices.C [new file with mode: 0644]
PWGLF/FORWARD/analysis2/scripts/multdist/AddTaskMultiplicity.C [new file with mode: 0644]
PWGLF/FORWARD/analysis2/scripts/multdist/CreateResponseMatrices.C [new file with mode: 0644]
PWGLF/FORWARD/analysis2/scripts/multdist/MakeMultAOD.C [new file with mode: 0644]
PWGLF/FORWARD/analysis2/scripts/multdist/MakeMultiplicityDistributions.C [new file with mode: 0644]
PWGLF/FORWARD/analysis2/scripts/multdist/doUnfolding.C [new file with mode: 0644]
PWGLF/FORWARD/analysis2/scripts/multdist/runMultiplicityDistributionAnalysis_README.txt [new file with mode: 0644]
PWGLF/FORWARD/analysis2/scripts/multdist/unfoldBase.C [new file with mode: 0644]
PWGLF/FORWARD/analysis2/scripts/multdist/unfoldChi2Method.C [new file with mode: 0644]
PWGLF/PWGLFforward2LinkDef.h

index a3715a67b6b74bc39ad7efa82e2dae123563cb48..b43b4b2ccb4765248420decbe5ef65269b240fa3 100644 (file)
@@ -76,6 +76,9 @@ set ( SRCS
   FORWARD/analysis2/AliPoissonCalculator.cxx
   FORWARD/analysis2/AliSPDMCTrackDensity.cxx
   FORWARD/GEO/AliAnalysisTaskZDCPbPb.cxx
+  FORWARD/analysis2/AliDisplacedVerticesESDFilterTask.cxx
+  FORWARD/analysis2/AliForwardMultiplicityDistribution.cxx
+  FORWARD/analysis2/AliForwardCreateResponseMatrices.cxx
 )
 
 string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" )
diff --git a/PWGLF/FORWARD/analysis2/AliForwardCreateResponseMatrices.cxx b/PWGLF/FORWARD/analysis2/AliForwardCreateResponseMatrices.cxx
new file mode 100644 (file)
index 0000000..aee0ae5
--- /dev/null
@@ -0,0 +1,481 @@
+
+#include <TH1D.h>
+#include "AliForwardCreateResponseMatrices.h"
+#include "AliForwardMultiplicityDistribution.h"
+#include "AliAODForwardMult.h"
+#include "AliAODCentralMult.h"
+#include "AliAODEvent.h"
+#include "AliFMDMCEventInspector.h"
+#include "AliAODMCHeader.h"
+#include "AliAODMCParticle.h"
+
+using namespace std;
+
+ClassImp(AliForwardCreateResponseMatrices)
+#if 0
+; // This is for Emacs - do not delete
+#endif
+//_____________________________________________________________________
+AliForwardCreateResponseMatrices::AliForwardCreateResponseMatrices() 
+  : AliBasedNdetaTask(),
+    fTrigger(),
+    fBins(),
+    fOutput()
+{
+  //
+  // Default Constructor
+  //
+}
+//_____________________________________________________________________
+AliForwardCreateResponseMatrices::AliForwardCreateResponseMatrices(const char* name) 
+  : AliBasedNdetaTask(name),
+    fTrigger(),
+    fBins(),
+    fOutput()
+{
+  //
+  // Constructor
+  //
+  DefineOutput(1, TList::Class());
+}
+
+//_____________________________________________________________________
+void AliForwardCreateResponseMatrices::UserCreateOutputObjects()
+{
+  //
+  // Create Output Objects
+  //
+  fOutput = new TList;
+  fOutput->SetOwner();
+  fOutput->SetName(GetName());
+  
+  TH1D* dndetaSumForward = new TH1D("dndetaSumForward","dndetaSumForward", 200,-4,6);
+  TH1D* dndetaSumCentral = new TH1D("dndetaSumCentral","dndetaSumCentral", 200,-4,6);
+  TH1D* dndetaSumMC = new TH1D("dndetaSumMC","dndetaSumMC", 200,-4,6);
+  fOutput->Add(dndetaSumForward);
+  fOutput->Add(dndetaSumCentral);
+  fOutput->Add(dndetaSumMC);
+
+  TH1D* dndetaEventForward = new TH1D();
+  TH1D* dndetaEventCentral = new TH1D();
+  TH1D* dndetaEventMC = new TH1D();
+  
+  fOutput->Add(dndetaEventForward);
+  fOutput->Add(dndetaEventCentral);
+  fOutput->Add(dndetaEventMC);
+
+  // make trigger diagnostics histogram
+  fTrigger= new TH1I();
+  fTrigger= AliAODForwardMult::MakeTriggerHistogram("triggerHist");
+  fOutput->Add(fTrigger);
+  
+
+  //Loop over all individual eta bins, and define their hisograms 
+  TIter next(&fBins);
+  Bin * bin = 0;
+  while ((bin = static_cast<Bin*>(next()))) { 
+    bin->DefineOutputs(fOutput, 400);
+  }
+  PostData(1, fOutput);
+  
+}
+
+
+//_____________________________________________________________________
+void AliForwardCreateResponseMatrices::UserExec(Option_t */*option*/)
+{
+  //
+  // User Exec
+  //
+  AliAODEvent* aod = dynamic_cast<AliAODEvent*>(InputEvent());
+  if (!aod) {
+    AliError("Cannot get the AOD event");
+    return;  } 
+  AliAODForwardMult* AODforward = static_cast<AliAODForwardMult*>(aod->FindListObject("Forward"));
+  if (!AODforward) {
+    AliError("Cannot get the AODforwardMult object");
+    return;  } 
+  
+  // Check the AOD event
+  Bool_t selectedTrigger = AODforward->CheckEvent(fTriggerMask, fVtxMin, fVtxMax,0,0,fTrigger);
+  
+  TH2D forward;
+  TH2D central;
+  
+  
+  // Get 2D eta-phi histograms for each event
+  GetHistograms(aod, forward, central);
+  
+  //check if event is NSD according to either MC truth or analysis (ESD)
+  Bool_t isMCNSD=kFALSE;
+  Bool_t isESDNSD=kFALSE;
+  if(AODforward->IsTriggerBits(AliAODForwardMult::kMCNSD))
+    isMCNSD=kTRUE; 
+  if(AODforward->IsTriggerBits(AliAODForwardMult::kNSD))
+    isESDNSD=kTRUE;
+  
+  //primary dndeta dist
+  TH2D* primHist;
+  TObject* oPrimary   = aod->FindListObject("primary");
+  primHist   = static_cast<TH2D*>(oPrimary);
+  
+  TH1D* dndetaSumForward    = (TH1D*)fOutput->FindObject("dndetaSumForward");
+  TH1D* dndetaSumCentral    = (TH1D*)fOutput->FindObject("dndetaSumCentral");
+  TH1D* dndetaSumMC         = (TH1D*)fOutput->FindObject("dndetaSumMC");
+  TH1D* dndetaEventForward  = (TH1D*)fOutput->FindObject("dndetaEventForward");
+  TH1D* dndetaEventCentral  = (TH1D*)fOutput->FindObject("dndetaEventCentral");
+  TH1D* dndetaEventMC       = (TH1D*)fOutput->FindObject("dndetaEventMC");
+
+  dndetaEventForward = forward.ProjectionX("dndetaForward",1,forward.GetNbinsY(),"");
+  dndetaEventCentral = central.ProjectionX("dndetaCentral",1,central.GetNbinsY(),"");
+  dndetaEventMC      = primHist->ProjectionX("dndetaMC",1,primHist->GetNbinsY(),"");
+   
+  // underflow eta bin of forward/central carry information on whether there is acceptance. 1= acceptance, 0= no acceptance
+  TH1D* normEventForward    = 0;
+  TH1D* normEventCentral    = 0;
+  TH1D* normEventMC         = 0;
+  normEventForward   = forward.ProjectionX("dndetaEventForwardNorm",0,0,"");
+  normEventCentral   = central.ProjectionX("dndetaEventCentralNorm",0,0,"");
+  normEventMC        = primHist->ProjectionX("dndetaEventNormMC",0,0,"");
+
+  dndetaSumForward->Add(dndetaEventForward);
+  dndetaSumCentral->Add(dndetaEventCentral);
+  dndetaSumMC->Add(dndetaEventMC);
+     
+  
+  Double_t VtxZ = AODforward->GetIpZ();
+
+  // loop over all eta bins, and create response matrices, trigger-vertex bias histograms etc
+  TIter next(&fBins);
+  Bin * bin = 0;
+  while ((bin = static_cast<Bin*>(next()))) { 
+    bin->Process(dndetaEventForward, dndetaEventCentral, 
+                normEventForward,   normEventCentral, dndetaEventMC, VtxZ, selectedTrigger,isMCNSD, isESDNSD,  aod);
+  }
+    
+  PostData(1, fOutput);
+  
+  
+}
+
+//_____________________________________________________________________
+void AliForwardCreateResponseMatrices::Terminate(Option_t */*option*/)
+{
+  //
+  // Terminate
+  //
+}
+
+//_________________________________________________________________-
+TH2D* AliForwardCreateResponseMatrices::GetHistogram(const AliAODEvent*, Bool_t )
+{
+  //
+  // implementation of pure virtual function, always returning 0
+  //
+  return 0;
+}
+
+void AliForwardCreateResponseMatrices::GetHistograms(const AliAODEvent* aod, TH2D& forward, TH2D& central)
+{
+  //
+  // Get single event forward and central d²N/dEta dPhi histograms 
+  //
+  TObject* forwardObj = 0;
+  TObject* centralObj = 0;
+  
+  forwardObj = aod->FindListObject("ForwardMC");
+  centralObj = aod->FindListObject("CentralClusters");
+   
+  if (!forwardObj) {
+    AliWarning("No MC forward object found in AOD");
+  }
+  if (!centralObj) {
+    AliWarning("No MC central object found in AOD");
+  }
+
+  AliAODForwardMult* AODforward = static_cast<AliAODForwardMult*>(forwardObj);
+  AliAODCentralMult* AODcentral = static_cast<AliAODCentralMult*>(centralObj);
+  
+  forward= AODforward->GetHistogram();
+  central= AODcentral->GetHistogram();
+  
+  
+}
+
+
+//=====================================================================
+AliForwardCreateResponseMatrices::Bin::Bin(Double_t etaLow, Double_t etaHigh) 
+  : TNamed("", ""), 
+    fEtaLow(etaLow),          // low eta limit 
+    fEtaHigh(etaHigh),        // high eta limit 
+    fHist(),                  // multiplicity histogram 
+    fHistMC(),                // multiplicity histogram MC truth primaries
+    fAcceptance(),            // histogram showing the 'holes' in acceptance. BinContent of 1 shows a hole, and BinContent of 10 shows data coverage
+    fVtxZvsNdataBins(),       // VtxZ vs. number of data acceptance bins (normalised to the eta range)
+    fResponseMatrix(),        // Response matrix (MC truth vs. analysed multiplicity)
+    fResponseMatrixPlus05(),  // Response matrix with analysed multiplicity scaled up by 5%
+    fResponseMatrixPlus075(), // Response matrix  with analysed multiplicity scaled up by 7.5%
+    fResponseMatrixPlus10(),  // Response matrix with analysed multiplicity scaled up by 10%
+    fResponseMatrixMinus05(), // Response matrix with analysed multiplicity scaled down by 5%
+    fResponseMatrixMinus075(),// Response matrix with analysed multiplicity scaled down by 7.55%
+    fResponseMatrixMinus10(), // Response matrix with analysed multiplicity scaled down by 10%
+    fResponseMatrixMinusSys(),// Response matrix with analysed multiplicity scaled up by event mult uncertainty
+    fResponseMatrixPlusSys(), // Response matrix with analysed multiplicity scaled down by event mult uncertainty
+    fESDNSD(),                // number of events found as NSD by the analysis vs. multiplicity
+    fMCNSD(),                 // number of events found as NSD by the MC truth vs. multiplicity
+    fMCESDNSD(),              // number of events found as NSD by both analysis and MC truth vs. multiplicity   
+    fTriggerBias()             // histogram for trigger vertex bias correction
+{
+  //
+  // Constructor
+  //
+  const char* name = AliForwardMultiplicityDistribution::FormBinName(etaLow,etaHigh);
+
+  SetName(name);
+  SetTitle(Form("%+4.1f < #eta < %+4.1f", fEtaLow, fEtaHigh));
+  
+}
+//_____________________________________________________________________
+AliForwardCreateResponseMatrices::Bin::Bin() 
+  : TNamed(), 
+    fEtaLow(),          // low eta limit 
+    fEtaHigh(),        // high eta limit 
+    fHist(),                  // multiplicity histogram 
+    fHistMC(),                // multiplicity histogram MC truth primaries
+    fAcceptance(),            // histogram showing the 'holes' in acceptance. BinContent of 1 shows a hole, and BinContent of 10 shows data coverage
+    fVtxZvsNdataBins(),       // VtxZ vs. number of data acceptance bins (normalised to the eta range)
+    fResponseMatrix(),        // Response matrix (MC truth vs. analysed multiplicity)
+    fResponseMatrixPlus05(),  // Response matrix with analysed multiplicity scaled up by 5%
+    fResponseMatrixPlus075(), // Response matrix  with analysed multiplicity scaled up by 7.5%
+    fResponseMatrixPlus10(),  // Response matrix with analysed multiplicity scaled up by 10%
+    fResponseMatrixMinus05(), // Response matrix with analysed multiplicity scaled down by 5%
+    fResponseMatrixMinus075(),// Response matrix with analysed multiplicity scaled down by 7.55%
+    fResponseMatrixMinus10(), // Response matrix with analysed multiplicity scaled down by 10%
+    fResponseMatrixMinusSys(),// Response matrix with analysed multiplicity scaled up by event mult uncertainty
+    fResponseMatrixPlusSys(), // Response matrix with analysed multiplicity scaled down by event mult uncertainty
+    fESDNSD(),                // number of events found as NSD by the analysis vs. multiplicity
+    fMCNSD(),                 // number of events found as NSD by the MC truth vs. multiplicity
+    fMCESDNSD(),              // number of events found as NSD by both analysis and MC truth vs. multiplicity        
+    fTriggerBias()             // histogram for trigger vertex bias correction
+{
+  //
+  // Default Constructor
+  //
+}
+//_____________________________________________________________________
+void AliForwardCreateResponseMatrices::Bin::DefineOutputs(TList* cont,  Int_t max)
+{
+  //
+  // Define eta bin output histos
+  //
+  TList* out = new TList;
+  out->SetName(GetName());
+  cont->Add(out);
+  
+  fHist                    = new TH1D("mult", GetTitle(), max, -0.5, max-.5);
+  fHistMC                  = new TH1D("multTruth", GetTitle(), max, -0.5, max-.5);
+  fVtxZvsNdataBins         = new TH2D("VtxZvsNdataBins", "VtxZ vs dataAcceptance/etaRange;z-vtz;dataAcceptance/etaRange", 20, -10,10, 130,0,1.3);
+  fAcceptance              = new TH2D("Acceptance","Acceptance;#eta;z-vtx",200,-4, 6 , 20,-10,10);
+  fResponseMatrix          = new TH2D("responseMatrix","Response Matrix;MC_{truth};MC_{measured}", max, -0.5, max-.5, max, -0.5, max-.5);
+  fResponseMatrixPlus05    = new TH2D("responseMatrixPlus05","Response Matrix;MC_{truth};MC_{measured}", max, -0.5, max-.5, max, -0.5, max-.5);
+  fResponseMatrixPlus075   = new TH2D("responseMatrixPlus075","Response Matrix;MC_{truth};MC_{measured}", max, -0.5, max-.5, max, -0.5, max-.5);
+  fResponseMatrixPlus10    = new TH2D("responseMatrixPlus10","Response Matrix;MC_{truth};MC_{measured}", max, -0.5, max-.5, max, -0.5, max-.5);
+  fResponseMatrixMinus05   = new TH2D("responseMatrixMinus05","Response Matrix;MC_{truth};MC_{measured}", max, -0.5, max-.5, max, -0.5, max-.5);
+  fResponseMatrixMinus075  = new TH2D("responseMatrixMinus075","Response Matrix;MC_{truth};MC_{measured}", max, -0.5, max-.5, max, -0.5, max-.5);
+  fResponseMatrixMinus10   = new TH2D("responseMatrixMinus10","Response Matrix;MC_{truth};MC_{measured}", max, -0.5, max-.5, max, -0.5, max-.5);
+  fResponseMatrixMinusSys  = new TH2D("responseMatrixMinusSys","Response Matrix;MC_{truth};MC_{measured}", max, -0.5, max-.5, max, -0.5, max-.5);
+  fResponseMatrixPlusSys   = new TH2D("responseMatrixPlusSys","Response Matrix;MC_{truth};MC_{measured}", max, -0.5, max-.5, max, -0.5, max-.5);
+  fTriggerBias             = new TH1D("triggerBias","triggerBias;MC_{truth};Correction}", max, -0.5, max-.5);
+  fMCNSD                   = new TH1D("fMCNSD","fMCNSD", 300,-0.5,299.5);
+  fESDNSD                  = new TH1D("fESDNSD","fESDNSD", 300,-0.5,299.5);
+  fMCESDNSD                = new TH1D("fMCESDNSD","fMCESDNSD", 300,-0.5,299.5);
+  
+  out->Add(fMCNSD);
+  out->Add(fESDNSD);
+  out->Add(fMCESDNSD);
+  out->Add(fHist);
+  out->Add(fHistMC);
+  out->Add(fVtxZvsNdataBins);
+  out->Add(fAcceptance);
+  out->Add(fResponseMatrix);
+  out->Add(fResponseMatrixPlus05);
+  out->Add(fResponseMatrixPlus075);
+  out->Add(fResponseMatrixPlus10);
+  out->Add(fResponseMatrixMinus05);
+  out->Add(fResponseMatrixMinus075);
+  out->Add(fResponseMatrixMinus10);
+  out->Add(fResponseMatrixPlusSys);
+  out->Add(fResponseMatrixMinusSys);
+  out->Add(fTriggerBias);
+  
+}
+
+//_____________________________________________________________________
+void AliForwardCreateResponseMatrices::Bin::Process(TH1D* dndetaForward, TH1D* dndetaCentral,
+                                                   TH1D* normForward,   TH1D* normCentral, TH1D* mc, Double_t VtxZ, Bool_t selectedTrigger, Bool_t isMCNSD, Bool_t isESDNSD, AliAODEvent* aodevent) 
+{
+  //
+  // Process a single eta bin
+  //  
+
+  // Diagnostics on event acceptance
+  Int_t    first = dndetaForward->GetXaxis()->FindBin(fEtaLow);
+  Int_t    last  = dndetaForward->GetXaxis()->FindBin(fEtaHigh-.0001);
+  Double_t acceptanceBins=0;
+  for(Int_t n= first;n<=last;n++){
+    if(normForward->GetBinContent(n)>0||normCentral->GetBinContent(n)>0){
+      acceptanceBins++;
+    }
+    fAcceptance->SetBinContent(n, fAcceptance->GetYaxis()->FindBin(VtxZ), 1);
+    if(normForward->GetBinContent(n)>0||normCentral->GetBinContent(n)>0)
+      fAcceptance->SetBinContent(n, fAcceptance->GetYaxis()->FindBin(VtxZ),10);
+  }
+  fVtxZvsNdataBins->Fill(VtxZ, (Double_t)acceptanceBins/(last-first+1));
+
+
+
+  Double_t c        = 0;
+  Double_t e2       = 0;
+  Double_t cPrimary = 0;
+  Double_t e2Primary= 0;
+    
+  for (Int_t i = first; i <= last; i++){ 
+    Double_t cForward  = 0;
+    Double_t cCentral  = 0;
+    Double_t e2Forward = 0;
+    Double_t e2Central = 0;
+    Double_t cMC  = 0;
+    Double_t e2MC = 0;
+    cMC= mc->GetBinContent(i);
+    e2MC= mc->GetBinError(i) * mc->GetBinError(i);
+    cPrimary+=cMC;
+    e2Primary+=e2MC;
+    if (normForward->GetBinContent(i) > 0) {
+      cForward  = dndetaForward->GetBinContent(i);
+      e2Forward = dndetaForward->GetBinError(i) * dndetaForward->GetBinError(i);
+    }
+    if (normCentral->GetBinContent(i) > 0) { 
+      cCentral  = dndetaCentral->GetBinContent(i);
+      e2Central = dndetaCentral->GetBinError(i) * dndetaCentral->GetBinError(i);
+    }
+    Double_t cc  = 0;
+    Double_t ee2 = 0;
+    if (cCentral > 0 && cForward > 0) { 
+      cc  = 0.5 * (cForward  + cCentral);
+      ee2 = 0.5 * (e2Forward + e2Central);
+    }
+    else if (cCentral > 0) { 
+      cc  = cCentral;
+      ee2 = e2Central;
+    }
+    else if (cForward > 0) { 
+      cc  = cForward;
+      ee2 = e2Forward;
+    }
+    c  += cc;
+    e2 += ee2;
+  }
+  
+  //retreive MC particles from event
+  TClonesArray* mcArray = (TClonesArray*)aodevent->FindListObject(AliAODMCParticle::StdBranchName());
+  if(!mcArray){
+        AliWarning("No MC array found in AOD. Try making it again.");
+    // return kFALSE;
+  }
+  AliAODMCHeader* header = dynamic_cast<AliAODMCHeader*>(aodevent->FindListObject(AliAODMCHeader::StdBranchName()));
+  if (!header) {
+    AliWarning("No header file found.");
+  }
+  
+  
+  Int_t ntracks = mcArray->GetEntriesFast();
+  // Track loop: find MC truth multiplicity in selected eta bin
+  Double_t mult = 0;
+  for (Int_t it = 0; it < ntracks; it++) {
+    AliAODMCParticle* particle = (AliAODMCParticle*)mcArray->At(it);
+    if (!particle) {
+      AliError(Form("Could not receive track %d", it));
+      continue;
+    }
+    if (!particle->IsPhysicalPrimary()) continue;
+    if (particle->Charge() == 0) continue;
+    if(particle->Eta()>fEtaLow&&particle->Eta()<fEtaHigh-0.0001)
+      mult++;
+  }
+  //fill fMCNSD with multiplicity of MC truth NSD events with MC truth |vtxz|<4
+  if(header->GetVtxZ()>-4&&header->GetVtxZ()<4){
+    if(isMCNSD){
+      fMCNSD->Fill(mult);
+    }
+  }
+  //fill fESDNSD with multiplicity from events with a reconstructed NSD trigger and reconstructed |vtxz|<4
+  if(VtxZ>-4&&VtxZ<4){
+    if(isESDNSD){
+      fESDNSD->Fill(mult);
+    }
+  }
+  // fullfilling both requirements of fMCNSD and fESDNSD
+  if(header->GetVtxZ()>-4&&header->GetVtxZ()<4&&VtxZ>-4&&VtxZ<4&&isMCNSD&&isESDNSD){
+    fMCESDNSD->Fill(mult);
+  }
+  
+
+
+//Systematic errors from here
+
+  Double_t fmd=0;
+  Double_t spd=0;
+  Double_t overlap=0;
+    
+  // number of eta bins in fmd, spd and overlap respectively
+  for(Int_t i = first;i<=last;i++){
+    if(normForward->GetBinContent(i)>0&&normCentral->GetBinContent(i)<1)
+      fmd++;
+    if(normForward->GetBinContent(i)>0&&normCentral->GetBinContent(i)>0)
+      overlap++;
+    if(normCentral->GetBinContent(i)>0&&normForward->GetBinContent(i)<1)
+      spd++;
+  }
+  
+  Double_t sysErrorSquared = 0;  
+
+  // estimate of systematic uncertainty on the event multiplicity - hardcoded :(. estimates taken from Hans Hjersing Dalsgaard or Casper Nygaard phd theses.
+  Double_t fmdSysError= 0.08;
+  Double_t spdSysError= 0.04;
+  Double_t total = 0;
+  total= fmd + spd + overlap; 
+  if(total){  
+    // Combined systematc event uncertainty, by weighting with the number of eta-bins of fmd-only, spd-only and the overlap.
+    sysErrorSquared= (fmd*TMath::Power(fmdSysError,2)+ spd*TMath::Power(spdSysError,2)+
+                     0.5*overlap*TMath::Power(fmdSysError,2)+ 0.5*overlap*TMath::Power(spdSysError,2))/total;
+  }
+  
+  
+  if(selectedTrigger){
+    fHist->Fill(c);
+    fHistMC->Fill(cPrimary);
+    fResponseMatrix->Fill(cPrimary, c);
+    fResponseMatrixPlusSys->Fill(cPrimary,(1+TMath::Sqrt(sysErrorSquared))*c);
+    fResponseMatrixMinusSys->Fill(cPrimary,(1-TMath::Sqrt(sysErrorSquared))*c);
+    fResponseMatrixPlus05->Fill(cPrimary, 1.05*c);
+    fResponseMatrixPlus075->Fill(cPrimary, 1.075*c);
+    fResponseMatrixPlus10->Fill(cPrimary, 1.1*c);
+    fResponseMatrixMinus05->Fill(cPrimary, 0.95*c);
+    fResponseMatrixMinus075->Fill(cPrimary, 0.925*c);
+    fResponseMatrixMinus10->Fill(cPrimary, 0.9*c);
+     
+  }
+  
+}
+
+
+
+
+//_____________________________________________________________________
+//
+//
+// EOF
diff --git a/PWGLF/FORWARD/analysis2/AliForwardCreateResponseMatrices.h b/PWGLF/FORWARD/analysis2/AliForwardCreateResponseMatrices.h
new file mode 100644 (file)
index 0000000..65a565e
--- /dev/null
@@ -0,0 +1,138 @@
+#ifndef ALIFORWARDCREATERESPONSEMATRICES_H
+#define ALIFORWARDCREATERESPONSEMATRICES_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+  
+#include "AliAnalysisTaskSE.h"
+#include "AliBasedNdetaTask.h"
+#include <TList.h>
+#include <iostream>
+
+class TH2D;
+
+/**
+ * Task to do the multiplicity distibution
+ * 
+ */
+class AliForwardCreateResponseMatrices : public AliBasedNdetaTask
+{
+public:
+  /** 
+   * 
+   * Default Constructor 
+   * 
+   */
+  AliForwardCreateResponseMatrices();
+  /**
+   *  Constructor
+   *
+   */
+  AliForwardCreateResponseMatrices(const char* name);
+  /**
+   * Copy Constructor
+   *
+   */
+  AliForwardCreateResponseMatrices(const AliForwardCreateResponseMatrices& o) : AliBasedNdetaTask(o), fTrigger(),fBins(), fOutput() { }
+  /**
+   * Assignment operator
+   *
+   */
+  AliForwardCreateResponseMatrices& operator=(const AliForwardCreateResponseMatrices&){return *this;}
+  /** 
+   * 
+   * Destructor
+   * 
+   */
+  virtual ~AliForwardCreateResponseMatrices(){}
+  /**
+   *  Embedded Class begins here
+   */
+  struct Bin : public TNamed
+  {
+    /**
+     * Default Constructor
+     */
+    Bin();
+    /**
+     * Constructor
+     */
+    Bin(Double_t etaLow, Double_t etaHigh);
+    /**
+     * Copy Constructor
+     */ 
+    Bin(const Bin&);
+    /**
+     * Assignment Operator
+     */
+    Bin&operator=(const Bin&){return*this;}
+    /**
+     * Destructor
+     */
+    ~Bin(){}
+    /**
+     * Define outputs of a single eta bin
+     */
+    virtual void DefineOutputs(TList* cont,  Int_t max);
+    /**
+     * Process a single eta bin
+     */    
+    virtual void Process(TH1D* dndetaForward, TH1D* dndetaCentral, TH1D* normForward,   TH1D* normCentral, TH1D* dndetaMC, Double_t VtxZ, Bool_t selectedTrigger,  Bool_t isMCNSDm, Bool_t isESDNSD, AliAODEvent* aodevent);
+    Double_t fEtaLow;                  // low eta limit 
+    Double_t fEtaHigh;                 // high eta limit 
+    TH1D*    fHist;                    // multiplicity histogram 
+    TH1D*    fHistMC;                  // multiplicity histogram MC truth primaries
+    TH2D*    fAcceptance;              // histogram showing the 'holes' in acceptance. BinContent of 1 shows a hole, and BinContent of 10 shows data coverage
+    TH2D*    fVtxZvsNdataBins;         // VtxZ vs. number of data acceptance bins (normalised to the eta range) 
+    TH2D*    fResponseMatrix;          //Response matrix (MC truth vs. analysed multiplicity)
+    TH2D*    fResponseMatrixPlus05;    //Response matrix with analysed multiplicity scaled up by 5%
+    TH2D*    fResponseMatrixPlus075;   //Response matrix  with analysed multiplicity scaled up by 7.5%
+    TH2D*    fResponseMatrixPlus10;    //Response matrix with analysed multiplicity scaled up by 10%
+    TH2D*    fResponseMatrixMinus05;   //Response matrix with analysed multiplicity scaled down by 5%
+    TH2D*    fResponseMatrixMinus075;  //Response matrix with analysed multiplicity scaled down by 7.55%
+    TH2D*    fResponseMatrixMinus10;   //Response matrix with analysed multiplicity scaled down by 10%
+    TH2D*    fResponseMatrixMinusSys;  //Response matrix with analysed multiplicity scaled up by event mult uncertainty
+    TH2D*    fResponseMatrixPlusSys;   //Response matrix with analysed multiplicity scaled down by event mult uncertainty
+    TH1D*    fESDNSD;                  //number of events found as NSD by the analysis vs. multiplicity
+    TH1D*    fMCNSD;                   //number of events found as NSD by the MC truth vs. multiplicity
+    TH1D*    fMCESDNSD;                //number of events found as NSD by both analysis and MC truth vs. multiplicity
+    TH1D*    fTriggerBias;             // histogram for trigger vertex bias correction
+   ClassDef(Bin,1); // Manager of data 
+  };
+  /**
+   * Create Output Objects
+   */
+  virtual void UserCreateOutputObjects();
+  /**
+   * User Exec
+   */
+  void UserExec(Option_t *option);
+  /**
+   * Terminate
+   */
+  void Terminate(Option_t *option);
+  /** 
+   * implementation of pure virtual function, always returning 0
+   */
+  virtual TH2D* GetHistogram(const AliAODEvent* aod, Bool_t mc);
+  /** 
+   * Get single event forward and central @f$d^2N/d\eta d\phi@f$
+   * histograms
+   * 
+   */  
+  virtual void GetHistograms(const AliAODEvent* aod, TH2D& forward, TH2D& central); 
+  /** 
+   * Add another eta bin to the task
+   */
+  void AddBin(Double_t etaLow, Double_t etaHigh){fBins.Add(new Bin(etaLow, etaHigh)); }
+ protected:
+  TH1I* fTrigger;  //Trigger histogram
+  TList  fBins;    // List of eta bins
+  TList*  fOutput; // Output list
+  ClassDef(AliForwardCreateResponseMatrices, 1); 
+};
+
+#endif
+// Local Variables:
+//   mode: C++
+// End:
diff --git a/PWGLF/FORWARD/analysis2/AliForwardMultiplicityDistribution.cxx b/PWGLF/FORWARD/analysis2/AliForwardMultiplicityDistribution.cxx
new file mode 100644 (file)
index 0000000..a7117dd
--- /dev/null
@@ -0,0 +1,419 @@
+
+#include <TH1D.h>
+#include "AliForwardMultiplicityDistribution.h"
+#include "AliAODForwardMult.h"
+#include "AliAODCentralMult.h"
+#include "AliAODEvent.h"
+#include "TString.h"
+using namespace std;
+
+ClassImp(AliForwardMultiplicityDistribution)
+#if 0
+; // This is for Emacs - do not delete
+#endif
+//______________________________________________________________________
+AliForwardMultiplicityDistribution::AliForwardMultiplicityDistribution() 
+  : AliBasedNdetaTask(), 
+    fTrigger(0),   // trigger histogram
+    fBins(),       // eta bin list
+    fOutput(0),    // output list
+    fLowCent(-1),  // lower centrality limit
+    fHighCent(-1), // upper centrality limit
+    fNBins(-1),    // multiplicity axis' runs from 0 to fNbins
+    fCent(0)       // centrality
+{
+  //
+  // Default Constructor
+  //
+}
+//_____________________________________________________________________
+AliForwardMultiplicityDistribution::AliForwardMultiplicityDistribution(const char* name) 
+  : AliBasedNdetaTask(name),
+    fTrigger(0),   // trigger histogram
+    fBins(),       // eta bin list
+    fOutput(0),    // output list
+    fLowCent(-1),  // lower centrality limit
+    fHighCent(-1), // upper centrality limit
+    fNBins(-1),    // multiplicity axis' runs from 0 to fNbins
+    fCent(0)       // centrality
+{
+  //
+  // Constructor
+  //
+  DefineOutput(1, TList::Class());
+}
+
+//_____________________________________________________________________
+void AliForwardMultiplicityDistribution::UserCreateOutputObjects()
+{
+  //
+  // Create Output Objects
+  //
+  fOutput = new TList;
+  fOutput->SetOwner();
+  fOutput->SetName(GetName());
+  
+  TH1D* dndetaSumForward = new TH1D("dndetaSumForward","dndetaSumForward", 200,-4,6);
+  TH1D* dndetaSumCentral = new TH1D("dndetaSumCentral","dndetaSumCentral", 200,-4,6);
+  fOutput->Add(dndetaSumForward);
+  fOutput->Add(dndetaSumCentral);
+  
+  fCent = new TH1D("cent", "Centrality", 100, 0, 100);
+  fCent->SetDirectory(0);
+  fCent->SetXTitle(0);
+  fOutput->Add(fCent);
+  
+  TH1D* dndetaEventForward = new TH1D();
+  TH1D* dndetaEventCentral = new TH1D();
+  fOutput->Add(dndetaEventForward);
+  fOutput->Add(dndetaEventCentral);
+  
+  //make trigger diagnostics histogram
+  fTrigger= new TH1I();
+  fTrigger= AliAODForwardMult::MakeTriggerHistogram("triggerHist");
+  fOutput->Add(fTrigger);
+  
+  //Loop over all individual eta bins, and define their hisograms 
+  TIter next(&fBins);
+  Bin * bin = 0;
+  while ((bin = static_cast<Bin*>(next()))) { 
+    bin->DefineOutputs(fOutput, fNBins);
+  }
+  //fOutput->ls();
+  PostData(1, fOutput);
+}
+
+
+//_____________________________________________________________________
+void AliForwardMultiplicityDistribution::UserExec(Option_t */*option*/)
+{
+  //
+  // User Exec
+  //
+  AliAODEvent* aod = dynamic_cast<AliAODEvent*>(InputEvent());
+  if (!aod) {
+    AliError("Cannot get the AOD event");
+    return;
+  } 
+  // Check AOD event
+  AliAODForwardMult* AODforward = static_cast<AliAODForwardMult*>(aod->FindListObject("Forward"));
+  if (!AODforward->CheckEvent(fTriggerMask, fVtxMin, fVtxMax, 0,0,fTrigger)){
+    AliError("Invalid Event");
+    return;
+  }
+  
+  // Fill centrality histogram, only useful for PbPb 
+  Float_t cent = AODforward->GetCentrality();
+  if(fLowCent<fHighCent){
+    if(cent<fLowCent||cent>fHighCent) return;
+  }
+  fCent->Fill(cent);
+
+  TH2D forward;
+  TH2D central;
+  
+  // Get 2D eta-phi histograms for each event
+  GetHistograms(aod, forward, central);
+  
+  TH1D* dndetaSumForward    = (TH1D*)fOutput->FindObject("dndetaSumForward");
+  TH1D* dndetaSumCentral    = (TH1D*)fOutput->FindObject("dndetaSumCentral");
+  TH1D* dndetaEventForward  = (TH1D*)fOutput->FindObject("dndetaEventForward");
+  TH1D* dndetaEventCentral  = (TH1D*)fOutput->FindObject("dndetaEventCentral");
+  TH1D* normEventForward    = 0;
+  TH1D* normEventCentral    = 0;
+
+  dndetaEventForward = forward.ProjectionX("dndetaForward",1,forward.GetNbinsY(),"");
+  dndetaEventCentral = central.ProjectionX("dndetaCentral",1,central.GetNbinsY(),"");
+  dndetaSumForward->Add(dndetaEventForward);
+  dndetaSumCentral->Add(dndetaEventCentral);
+  
+  // underflow eta bin of forward/central carry information on whether there is acceptance. 1= acceptance, 0= no acceptance
+  normEventForward   = forward.ProjectionX("dndetaEventForwardNorm",0,0,"");
+  normEventCentral   = central.ProjectionX("dndetaEventCentralNorm",0,0,"");
+  
+  Double_t VtxZ = AODforward->GetIpZ();
+
+  
+  // loop over all eta bins, and fill multiplicity histos
+  TIter next(&fBins);
+  Bin * bin = 0;
+  while ((bin = static_cast<Bin*>(next()))) { 
+    bin->Process(dndetaEventForward, dndetaEventCentral, 
+                normEventForward,   normEventCentral, VtxZ);
+  }
+  
+  PostData(1, fOutput);
+  
+}
+
+//_____________________________________________________________________
+void AliForwardMultiplicityDistribution::Terminate(Option_t */*option*/)
+{
+  //
+  // Terminate
+  //
+}
+
+//_________________________________________________________________-
+TH2D* AliForwardMultiplicityDistribution::GetHistogram(const AliAODEvent*, Bool_t)
+{
+  //
+  // implementation of pure virtual function, always returning 0
+  //
+  return 0;
+}
+
+void AliForwardMultiplicityDistribution::GetHistograms(const AliAODEvent* aod, TH2D& forward, TH2D& central, Bool_t /*mc*/)
+{
+  //
+  // Get single event forward and central dN²/dEta dPhi histograms 
+  //
+  TObject* forwardObj = 0;
+  TObject* centralObj = 0;
+
+    
+  forwardObj = aod->FindListObject("Forward");
+  centralObj = aod->FindListObject("CentralClusters");
+
+  AliAODForwardMult* AODforward = static_cast<AliAODForwardMult*>(forwardObj);
+  AliAODCentralMult* AODcentral = static_cast<AliAODCentralMult*>(centralObj);
+  
+  forward= AODforward->GetHistogram();
+  central= AODcentral->GetHistogram();
+  
+}
+
+//______________________________________________________________________
+#if 1
+const Char_t*
+#else
+TString
+#endif
+AliForwardMultiplicityDistribution::FormBinName(Double_t etaLow, Double_t etaHigh)
+{
+  //
+  //  Form name of eta bin
+  //
+  TString sLow(TString::Format("%+03d",int(10*etaLow)));
+  TString sHigh(TString::Format("%+03d",int(10*etaHigh)));
+  sLow.ReplaceAll("+", "plus");
+  sLow.ReplaceAll("-", "minus");
+  sHigh.ReplaceAll("+", "plus");
+  sHigh.ReplaceAll("-", "minus");
+ #if 0
+  TString* combined= new TString();
+  combined->Append(TString::Format("%s_%s", sLow.Data(), sHigh.Data()));
+  return combined->Data();
+#else 
+  static TString tmp;
+  tmp = TString::Format("%s_%s", sLow.Data(), sHigh.Data());
+  return tmp.Data();
+#endif
+}
+
+//=====================================================================
+AliForwardMultiplicityDistribution::Bin::Bin() 
+  : TNamed(), 
+    fEtaLow(-1111),     // low eta limit 
+    fEtaHigh(-1111),    // high eta limit 
+    fHist(0),           // multiplicity distribution hist 
+    fHistPlus05(0),     // multiplicity distribution hist scaled up with 5%
+    fHistPlus075(0),    // multiplicity distribution hist scaled up with 7.5%
+    fHistPlus10(0),     // multiplicity distribution hist scaled up with 10%
+    fHistMinus05(0),    // multiplicity distribution hist scaled down with 5%
+    fHistMinus075(0),   // multiplicity distribution hist scaled down with 7.5%
+    fHistMinus10(0),    // multiplicity distribution hist scaled down with 10% 
+    fHistPlusSys(0),    // multiplicity distribution hist scaled up with the event uncertainty
+    fHistMinusSys(0),   // multiplicity distribution hist scaled down with the event uncertainty
+    fAcceptance(0),     // histogram showing the 'holes' in acceptance. 
+    fVtxZvsNdataBins(0) // VtxZ vs. number of data acceptance bins (normalised to the eta range) 
+{
+  // 
+  // Default Constructor
+  //
+}
+//_____________________________________________________________________
+AliForwardMultiplicityDistribution::Bin::Bin(Double_t etaLow, Double_t etaHigh) 
+  : TNamed("", ""), 
+    fEtaLow(etaLow),    // low eta limit 
+    fEtaHigh(etaHigh),  // high eta limit 
+    fHist(0),           // multiplicity distribution hist 
+    fHistPlus05(0),     // multiplicity distribution hist scaled up with 5%
+    fHistPlus075(0),    // multiplicity distribution hist scaled up with 7.5%
+    fHistPlus10(0),     // multiplicity distribution hist scaled up with 10%
+    fHistMinus05(0),    // multiplicity distribution hist scaled down with 5%
+    fHistMinus075(0),   // multiplicity distribution hist scaled down with 7.5%
+    fHistMinus10(0),    // multiplicity distribution hist scaled down with 10% 
+    fHistPlusSys(0),    // multiplicity distribution hist scaled up with the event uncertainty
+    fHistMinusSys(0),   // multiplicity distribution hist scaled down with the event uncertainty
+    fAcceptance(0),     // histogram showing the 'holes' in acceptance. 
+    fVtxZvsNdataBins(0) // VtxZ vs. number of data acceptance bins (normalised to the eta range) 
+{
+  // 
+  // Constructor
+  //
+  const char* name = AliForwardMultiplicityDistribution::FormBinName(fEtaLow,fEtaHigh);
+
+  SetName(name);
+  SetTitle(TString::Format("%+4.1f < #eta < %+4.1f", fEtaLow, fEtaHigh));
+
+}
+//_____________________________________________________________________
+void AliForwardMultiplicityDistribution::Bin::DefineOutputs(TList* cont,  Int_t max)
+{
+  //
+  // Define eta bin output histograms
+  //
+  TList* out = new TList;
+  out->SetName(GetName());
+  cont->Add(out);
+  fHist            = new TH1D("mult", GetTitle(), max, -0.5, max-.5);
+  fHistPlus05      = new TH1D("multPlus05", GetTitle(), max, -0.5, max-.5);
+  fHistPlus075     = new TH1D("multPlus075", GetTitle(), max, -0.5, max-.5);
+  fHistPlus10      = new TH1D("multPlus10", GetTitle(), max, -0.5, max-.5);
+  fHistMinus05     = new TH1D("multMinus05", GetTitle(), max, -0.5, max-.5);
+  fHistMinus075    = new TH1D("multMinus075", GetTitle(), max, -0.5, max-.5);
+  fHistMinus10     = new TH1D("multMinus10", GetTitle(), max, -0.5, max-.5);
+  fHistPlusSys     = new TH1D("multPlusSys", GetTitle(), max, -0.5, max-.5);
+  fHistMinusSys    = new TH1D("multMinusSys", GetTitle(), max, -0.5, max-.5);
+  fVtxZvsNdataBins = new TH2D("VtxZvsNdataBins", "VtxZ vs dataAcceptance/etaRange;z-vtz;dataAcceptance/etaRange", 20, -10,10, 130,0,1.3);
+  fAcceptance      = new TH2D("Acceptance","Acceptance;#eta;z-vtx",200,-4, 6 , 20,-10,10);
+
+  out->Add(fHist);
+  out->Add(fHistPlus05);
+  out->Add(fHistPlus075);
+  out->Add(fHistPlus10);
+  out->Add(fHistMinus05);
+  out->Add(fHistMinus075);
+  out->Add(fHistMinus10);
+  out->Add(fHistPlusSys);
+  out->Add(fHistMinusSys);
+  out->Add(fVtxZvsNdataBins);
+  out->Add(fAcceptance);
+
+}
+//_____________________________________________________________________
+void AliForwardMultiplicityDistribution::Bin::Process(TH1D* dndetaForward, TH1D* dndetaCentral,
+                                                     TH1D* normForward,   TH1D* normCentral, Double_t VtxZ) 
+{
+  //
+  // Process single eta bin
+  //
+  
+  // Diagnostics on event acceptance
+  Int_t    first = dndetaForward->GetXaxis()->FindBin(fEtaLow);
+  Int_t    last  = dndetaForward->GetXaxis()->FindBin(fEtaHigh-.0001);
+  Double_t acceptanceBins=0;
+  for(Int_t n= first;n<=last;n++){
+    if(normForward->GetBinContent(n)>0||normCentral->GetBinContent(n)>0){
+      acceptanceBins++;
+    }
+    fAcceptance->SetBinContent(n, fAcceptance->GetYaxis()->FindBin(VtxZ), 1);
+    if(normForward->GetBinContent(n)>0||normCentral->GetBinContent(n)>0)
+      fAcceptance->SetBinContent(n, fAcceptance->GetYaxis()->FindBin(VtxZ),10);
+  }
+
+  //std::cout << "% of bins covered in eta-range: " << acceptanceBins/(last-first+1) << std::endl;
+  fVtxZvsNdataBins->Fill(VtxZ, (Double_t)acceptanceBins/(last-first+1));
+  
+
+  Double_t c     = 0;
+  Double_t e2    = 0;
+  for (Int_t i = first; i <= last; i++) { 
+    Double_t cForward  = 0;
+    Double_t cCentral  = 0;
+    Double_t e2Forward = 0;
+    Double_t e2Central = 0;
+    if (normForward->GetBinContent(i) > 0) {
+      cForward  = dndetaForward->GetBinContent(i);
+      e2Forward = dndetaForward->GetBinError(i) * dndetaForward->GetBinError(i);
+    }
+    if (normCentral->GetBinContent(i) > 0) { 
+      cCentral  = dndetaCentral->GetBinContent(i);
+      e2Central = dndetaCentral->GetBinError(i) * dndetaCentral->GetBinError(i);
+    }
+    Double_t cc  = 0;
+    Double_t ee2 = 0;
+    // if overlap between central/forward, use average of the two
+    if (cCentral > 0 && cForward > 0) { 
+      cc  = 0.5 * (cForward  + cCentral);
+      ee2 = 0.5 * (e2Forward + e2Central);
+    }
+    else if (cCentral > 0) { 
+      cc  = cCentral;
+      ee2 = e2Central;
+    }
+    else if (cForward > 0) { 
+      cc  = cForward;
+      ee2 = e2Forward;
+    }
+    c  += cc;
+    e2 += ee2;
+  }
+  
+  //Systematic errors from here
+  
+  Double_t fmd=0;
+  Double_t spd=0;
+  Double_t overlap=0;
+
+  // number of eta bins in fmd, spd and overlap respectively
+  for(Int_t i = first;i<=last;i++){
+    if(normForward->GetBinContent(i)>0&&normCentral->GetBinContent(i)<1)
+      fmd++;
+    if(normForward->GetBinContent(i)>0&&normCentral->GetBinContent(i)>0)
+      overlap++;
+    if(normCentral->GetBinContent(i)>0&&normForward->GetBinContent(i)<1)
+      spd++;
+  }
+  
+  Double_t sysErrorSquared = 0;
+  // estimate of systematic uncertainty on the event multiplicity - hardcoded :(. estimates taken from Hans Hjersing Dalsgaard or Casper Nygaard phd theses.
+  Double_t fmdSysError= 0.08;
+  Double_t spdSysError= 0.067;
+  
+  Double_t total = 0;
+  total= fmd+ spd + overlap; 
+  // Combined systematc event uncertainty, by weighting with the number of eta-bins of fmd-only, spd-only and the overlap.
+  if(total){
+    sysErrorSquared= (fmd*TMath::Power(fmdSysError,2)+ spd*TMath::Power(spdSysError,2)+
+                     0.5*overlap*TMath::Power(fmdSysError,2)+ 0.5*overlap*TMath::Power(spdSysError,2))/total;
+  }
+    
+  //Strangeness correction, assumed to be 2.0 +- 1% (syst),  taken from HHD and CN thesis
+  c=0.98*c;
+  
+
+  // correction for missing material description. Correction is based on separate PbPb analysis of difference between nominel and displaced vertices
+  if(fEtaHigh<1.55&&fEtaHigh>1.45)
+    c=0.98*c;
+  if(fEtaHigh<2.05&&fEtaHigh>1.95)
+    c=0.963*c;
+  if(fEtaHigh<2.45&&fEtaHigh>2.35)
+    c=0.956*c;
+  if(fEtaHigh>2.95)
+    c=0.955*c;
+  
+  fHist->Fill(c);
+  fHistPlusSys->Fill((1+TMath::Sqrt(sysErrorSquared))*c);
+  fHistMinusSys->Fill((1-TMath::Sqrt(sysErrorSquared))*c);
+  fHistPlus05->Fill(1.05*c);
+  fHistPlus075->Fill(1.075*c);
+  fHistPlus10->Fill(1.1*c);
+  fHistMinus05->Fill(0.95*c);
+  fHistMinus075->Fill(0.925*c);
+  fHistMinus10->Fill(0.9*c);
+
+  
+}
+
+
+
+//_____________________________________________________________________
+//
+//
+// EOF
diff --git a/PWGLF/FORWARD/analysis2/AliForwardMultiplicityDistribution.h b/PWGLF/FORWARD/analysis2/AliForwardMultiplicityDistribution.h
new file mode 100644 (file)
index 0000000..a916d9a
--- /dev/null
@@ -0,0 +1,141 @@
+#ifndef ALIFORWARDMULTIPLICITYDISTRIBUTION_H
+#define ALIFORWARDMULTIPLICITYDISTRIBUTION_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+  
+#include "AliAnalysisTaskSE.h"
+#include "AliBasedNdetaTask.h"
+#include <TList.h>
+#include <iostream>
+
+class TH2D;
+
+/**
+ * Task to do the multiplicity distibution
+ * 
+ */
+class AliForwardMultiplicityDistribution : public AliBasedNdetaTask
+{
+public:
+  /**
+   * Default Constructor
+   */
+  AliForwardMultiplicityDistribution();
+  /**
+   * Constructor
+   */
+  AliForwardMultiplicityDistribution(const char* name);
+  /**
+   * Copy Constructor
+   */ 
+  AliForwardMultiplicityDistribution(const AliForwardMultiplicityDistribution& o) : AliBasedNdetaTask(o), fTrigger(o.fTrigger),fBins(), fOutput(o.fOutput), fLowCent(o.fLowCent), fHighCent(o.fHighCent),fNBins(o.fNBins), fCent(o.fCent){ }
+  /**
+   * Assignment Operator
+   */
+  AliForwardMultiplicityDistribution& operator=(const AliForwardMultiplicityDistribution&){return *this;}
+  /**
+   * Destructor
+   */
+  virtual ~AliForwardMultiplicityDistribution(){}
+  /**
+   * Embedded Class begins here
+   */
+  struct Bin : public TNamed
+  {
+    /**
+     * Default Constructor
+     */
+    Bin();
+    /**
+     * Constructor
+     */
+    Bin(Double_t etaLow, Double_t etaHigh);
+    /**
+     * Copy Constructor
+     */    
+    Bin(const Bin&);
+    /**
+     * Assignment operator
+     */
+    Bin&operator=(const Bin&){return*this;}
+    /**
+     * Destructor
+     */    
+    ~Bin(){}
+    /**
+     *  Define outputs of a single eta bin
+     */
+    virtual void DefineOutputs(TList* cont,  Int_t max);
+    /**
+     * Process a single eta bin
+     */
+    virtual void Process(TH1D* dndetaForward, TH1D* dndetaCentral, TH1D* normForward,   TH1D* normCentral, Double_t VtxZ);
+    Double_t fEtaLow;           // low eta limit 
+    Double_t fEtaHigh;          // high eta limit 
+    TH1D*    fHist;             // multiplicity distribution hist 
+    TH1D*    fHistPlus05;       // multiplicity distribution hist scaled up with 5%
+    TH1D*    fHistPlus075;      // multiplicity distribution hist scaled up with 7.5%
+    TH1D*    fHistPlus10;       // multiplicity distribution hist scaled up with 10%
+    TH1D*    fHistMinus05;      // multiplicity distribution hist scaled down with 5%
+    TH1D*    fHistMinus075;     // multiplicity distribution hist scaled down with 7.5%
+    TH1D*    fHistMinus10;      // multiplicity distribution hist scaled down with 10%
+    TH1D*    fHistPlusSys;      // multiplicity distribution hist scaled up with the event uncertainty
+    TH1D*    fHistMinusSys;     // multiplicity distribution hist scaled down with the event uncertainty
+    TH2D*    fAcceptance;       // histogram showing the 'holes' in acceptance. 
+                                // BinContent of 1 shows a hole, and BinContent of 10 shows data coverage
+    TH2D*    fVtxZvsNdataBins;  // VtxZ vs. number of data acceptance bins (normalised to the eta range) 
+    
+    ClassDef(Bin,1);  // Manager of data 
+  };
+  /**
+   * Create Output Objects
+   */
+  virtual void UserCreateOutputObjects();
+  /**
+   * User Exec
+   */
+  void UserExec(Option_t *option);
+  /**
+   * Terminate
+   */
+  void Terminate(Option_t *option);
+  /**
+   * Set Centrality
+   */
+  void SetCentrality(Int_t lowCent, Int_t highCent){fLowCent= lowCent; fHighCent= highCent;}
+  /**
+   * Set fNBins, multiplicity histos run from 0 to fNBins
+   */
+  void SetNBins(Int_t n){fNBins= n;}
+  /** 
+   * implementation of pure virtual function, always returning 0
+   */
+  virtual TH2D* GetHistogram(const AliAODEvent* aod, Bool_t mc);
+  /**
+   * Get single event forward and central dN²/dEta dPhi histograms 
+   */
+  virtual void GetHistograms(const AliAODEvent* aod, TH2D& forward, TH2D& central , Bool_t mc=false);
+  /** 
+   * Add another eta bin to the task
+   */
+  void AddBin(Double_t etaLow, Double_t etaHigh){fBins.Add(new Bin(etaLow, etaHigh)); }
+  /**
+   *  Form name of eta bin
+   */
+  static const Char_t* FormBinName(Double_t etaLow, Double_t etaHigh);
+protected:
+  TH1I*  fTrigger;   // trigger histogram
+  TList  fBins;      // eta bin list
+  TList* fOutput;    // output list
+  Int_t  fLowCent;   // lower centrality limit
+  Int_t  fHighCent;  // upper centrality limit
+  Int_t  fNBins;     // multiplicity axis' runs from 0 to fNbins
+  TH1D*  fCent;      // centrality
+  ClassDef(AliForwardMultiplicityDistribution, 1); 
+};
+
+#endif
+// Local Variables:
+//   mode: C++
+// End:
diff --git a/PWGLF/FORWARD/analysis2/scripts/multdist/AddTaskCreateResponseMatrices.C b/PWGLF/FORWARD/analysis2/scripts/multdist/AddTaskCreateResponseMatrices.C
new file mode 100644 (file)
index 0000000..9b92528
--- /dev/null
@@ -0,0 +1,94 @@
+/**
+ * @file   AddTaskCentraldNdeta.C
+ * @author Christian Holm Christensen <cholm@nbi.dk>
+ * @date   Fri Jan 28 10:22:26 2011
+ * 
+ * @brief Script to add a multiplicity task for the central
+ *        @f$\eta@f$ region
+ * 
+ * 
+ */
+AliAnalysisTask*
+AddTaskCreateResponseMatrices(const char* trig="INEL", Double_t vzMin=-10, Double_t vzMax=10)
+{
+  // analysis manager
+  AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
+  
+  // Make our object.  2nd argumenent is absolute max Eta 
+  // 3rd argument is absolute max Vz
+  AliForwardCreateResponseMatrices* task = new AliForwardCreateResponseMatrices("ResponseMatrices");
+  //AliForwarddNdetaTask* task = new AliForwarddNdetaTask("Forward");
+  
+  task->SetVertexRange(vzMin, vzMax);
+  task->SetTriggerMask(trig);
+
+
+  //task->SetTriggerMask(AliAODForwardMult::kMCNSD); // trig);
+  
+  mgr->AddTask(task);
+  //Add Full eta-ranges
+  task->AddBin(-3.4,5.1);
+  
+  //Add Symmetric eta bins.
+  Double_t limits[] = { 3.4, 3.0, 2.5, 2.4, 2.0, 1.5, 1.4, 1.0, 0.5, 0. };
+  Double_t* limit = limits;
+  while ((*limit) > 0.1) { 
+    task->AddBin(-(*limit), +(*limit));
+    limit++;
+  }
+  
+  //Add 0-<eta> ranges
+  task->AddBin(0,5.1);
+  task->AddBin(0,5.0);
+  task->AddBin(0,4.5);
+  task->AddBin(0,4.0);
+  task->AddBin(0,3.5);
+
+ limit = limits;
+ while ((*limit) > 0.1) { 
+   task->AddBin(0, +(*limit));
+   limit++;
+ }
+
+ limit = limits;
+ while ((*limit) > 0.1) { 
+   task->AddBin(-(*limit),0);
+   limit++;
+ }
+
+
+ //Add 0.5 eta intervals
+ for (Double_t l = -3; l < 5; l += 0.5){ 
+   task->AddBin(l, l+.5);
+ }
+ //Add 0.20 eta intervals
+ for (Double_t l = -3; l < 5; l += 0.2){ 
+   task->AddBin(l, l+.2);
+ }
+  
+   
+  // create containers for input/output
+  AliAnalysisDataContainer *sums = 
+    mgr->CreateContainer("ResponseMatrices", TList::Class(), 
+                        AliAnalysisManager::kOutputContainer, 
+                       AliAnalysisManager::GetCommonFileName());
+ /*
+  AliAnalysisDataContainer *output = 
+  mgr->CreateContainer("CentralResults", TList::Class(), 
+  AliAnalysisManager::kParamContainer, 
+  AliAnalysisManager::GetCommonFileName());
+ */
+  // connect input/output
+ mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
+ mgr->ConnectOutput(task, 1, sums);
+ // mgr->ConnectOutput(task, 2, output);
+ return task;
+}
+
+  
+//________________________________________________________________________
+//
+// EOF
+// 
diff --git a/PWGLF/FORWARD/analysis2/scripts/multdist/AddTaskMultiplicity.C b/PWGLF/FORWARD/analysis2/scripts/multdist/AddTaskMultiplicity.C
new file mode 100644 (file)
index 0000000..3f0499d
--- /dev/null
@@ -0,0 +1,117 @@
+/**
+ * @file   AddTaskCentraldNdeta.C
+ * @author Christian Holm Christensen <cholm@nbi.dk>
+ * @date   Fri Jan 28 10:22:26 2011
+ * 
+ * @brief Script to add a multiplicity task for the central
+ *        @f$\eta@f$ region
+ * 
+ * 
+ */
+AliAnalysisTask*
+AddTaskMultiplicity(const char* trig="INEL", Double_t vzMin=-10, Double_t vzMax=10, Int_t lowCent, Int_t highCent, Int_t nBins)
+{
+  // analysis manager
+  AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
+  
+  // Make our object.  2nd argumenent is absolute max Eta 
+  // 3rd argument is absolute max Vz
+  AliForwardMultiplicityDistribution* task = new AliForwardMultiplicityDistribution("Mult");
+  //AliForwarddNdetaTask* task = new AliForwarddNdetaTask("Forward");
+
+  task->SetVertexRange(vzMin, vzMax);
+  task->SetTriggerMask(trig);
+  task->SetCentrality(lowCent, highCent);
+  task->SetNBins(nBins);
+  mgr->AddTask(task);
+  /*
+ if (useCent) {
+   Short_t bins[] = { 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 };
+   task->SetCentralityAxis(11, bins);
+   task->InitCentBins();
+ }
+  */
+  //Add Full eta-ranges
+  task->AddBin(-3.4,5.1);
+  
+  //Add Symmetric eta bins.
+  Double_t limits[] = { 3.4, 3.0, 2.5, 2.4, 2.0, 1.5, 1.4, 1.0, 0.5, 0. };
+  Double_t* limit = limits;
+  while ((*limit) > 0.1) { 
+    task->AddBin(-(*limit), +(*limit));
+    limit++;
+  }
+  
+  //Add 0-<eta> ranges
+  task->AddBin(0,5.1);
+  task->AddBin(0,5.0);
+  task->AddBin(0,4.5);
+  task->AddBin(0,4.0);
+  task->AddBin(0,3.5);
+
+ limit = limits;
+ while ((*limit) > 0.1) { 
+   task->AddBin(0, +(*limit));
+   limit++;
+ }
+
+ limit = limits;
+ while ((*limit) > 0.1) { 
+   task->AddBin(-(*limit),0);
+   limit++;
+ }
+
+
+ //Add 0.5 eta intervals
+ for (Double_t l = -3; l < 5; l += 0.5){ 
+   task->AddBin(l, l+.5);
+ }
+ //Add 0.20 eta intervals
+ for (Double_t l = -3; l < 5; l += 0.2){ 
+   task->AddBin(l, l+.2);
+ }
+ /*
+ task->AddBin(-3.0,-2.5);
+ task->AddBin(-2.5,-2.0);
+ task->AddBin(-2.0,-1.5);
+ task->AddBin(-1.5,-1.0);
+ task->AddBin(-1.0,-0.5);
+ task->AddBin(0.5,1.0);
+ task->AddBin(1.0,1.5);
+ task->AddBin(1.5,2.0);
+ task->AddBin(2.0,2.5);
+ task->AddBin(2.5,3.0);
+ task->AddBin(3.0,3.5);
+ task->AddBin(3.5,4.0);
+ task->AddBin(4.0,4.5);
+ task->AddBin(4.5,5.0);
+ */
+
+ // create containers for input/output
+ AliAnalysisDataContainer *sums = 
+   mgr->CreateContainer("CentralSums", TList::Class(), 
+                       AliAnalysisManager::kOutputContainer, 
+                       AliAnalysisManager::GetCommonFileName());
+ /*
+  AliAnalysisDataContainer *output = 
+  mgr->CreateContainer("CentralResults", TList::Class(), 
+  AliAnalysisManager::kParamContainer, 
+  AliAnalysisManager::GetCommonFileName());
+ */
+  // connect input/output
+ mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
+ mgr->ConnectOutput(task, 1, sums);
+ // mgr->ConnectOutput(task, 2, output);
+ return task;
+}
+
+  
+//________________________________________________________________________
+//
+// EOF
+// 
diff --git a/PWGLF/FORWARD/analysis2/scripts/multdist/CreateResponseMatrices.C b/PWGLF/FORWARD/analysis2/scripts/multdist/CreateResponseMatrices.C
new file mode 100644 (file)
index 0000000..3bb1d0c
--- /dev/null
@@ -0,0 +1,84 @@
+/**
+ * @file 
+ * 
+ * @ingroup pwg2_forward_scripts
+ */
+
+/** 
+ * Run first pass of the analysis - that is read in ESD and produce AOD
+ * 
+ * @param esddir    ESD input directory. Any file matching the pattern 
+ *                  *AliESDs*.root are added to the chain 
+ * @param nEvents   Number of events to process.  If 0 or less, then 
+ *                  all events are analysed
+ * @param proof     Run in proof mode 
+ * @param mc        Run over MC data
+ *
+ * If PROOF mode is selected, then Terminate will be run on the master node 
+ * in any case. 
+ * 
+ *
+ * @ingroup pwg2_forward_scripts
+ */
+void CreateResponseMatrices(const char* aoddir=".", 
+               Int_t       nEvents=-1, 
+               const char* trig="INEL",
+               Double_t    vzMin=-10,
+               Double_t    vzMax=10,
+               char* output="responseMatrices.root",
+               Int_t       proof=0)
+{
+  // --- Libraries to load -------------------------------------------
+  gROOT->Macro("$ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/LoadLibs.C");
+
+  // --- Check for proof mode, and possibly upload pars --------------
+  if (proof> 0) { 
+    gROOT->LoadMacro("$ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/LoadPars.C");
+    LoadPars(proof);
+  }
+  
+  // --- Our data chain ----------------------------------------------
+  gROOT->LoadMacro("$ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/MakeChain.C");
+  TChain* chain = MakeChain("AOD", aoddir,true);
+  // If 0 or less events is select, choose all 
+  if (nEvents <= 0) nEvents = chain->GetEntries();
+
+  // --- Creating the manager and handlers ---------------------------
+  AliAnalysisManager *mgr  = new AliAnalysisManager("Forward Train", 
+                                                   "Forward Multiplicity");
+  AliAnalysisManager::SetCommonFileName(output);
+
+  // --- ESD input handler -------------------------------------------
+  AliAODInputHandler *aodInputHandler = new AliAODInputHandler();
+  mgr->SetInputEventHandler(aodInputHandler);      
+       
+  // --- Add tasks ---------------------------------------------------
+  // Forward 
+  gROOT->LoadMacro("$ALICE_ROOT/PWGLF/FORWARD/analysis2/AddTaskCreateResponseMatrices.C");
+  AddTaskCreateResponseMatrices(trig, vzMin, vzMax);
+  
+  
+  // --- Run the analysis --------------------------------------------
+  TStopwatch t;
+  if (!mgr->InitAnalysis()) {
+    Error("MakeMult", "Failed to initialize analysis train!");
+    return;
+  }
+  // Skip terminate if we're so requested and not in Proof or full mode
+  mgr->SetSkipTerminate(false);
+  // Some informative output 
+  mgr->PrintStatus();
+  // mgr->SetDebugLevel(3);
+  if (mgr->GetDebugLevel() < 1 && !proof) 
+    mgr->SetUseProgressBar(kTRUE,100);
+
+  // Run the train 
+  t.Start();
+  Printf("=== RUNNING ANALYSIS ==================================");
+  mgr->StartAnalysis(proof ? "proof" : "local", chain, nEvents);
+  t.Stop();
+  t.Print();
+}
+//
+// EOF
+//
diff --git a/PWGLF/FORWARD/analysis2/scripts/multdist/MakeMultAOD.C b/PWGLF/FORWARD/analysis2/scripts/multdist/MakeMultAOD.C
new file mode 100644 (file)
index 0000000..cef2088
--- /dev/null
@@ -0,0 +1,196 @@
+/*
+void Setup900GeV(MakeAODTrain p);
+void Setup900GeVMC(MakeAODTrain p);
+void Add900GeVruns(MakeAODTrain p);
+*/
+class MakeAODTrain;
+
+void MakeMultAOD()
+{
+  bool        usePar  =  true;
+  bool        mc      =  true;
+  Int_t       nEvents = -1;
+  UShort_t    proof   = 0;
+
+  const char* name    = "test_flatMult_withbg";
+  UShort_t    type    = 1;  // pp==1, PbPb==2
+  UShort_t    cms     = 900; // 2750 for 2760GeV
+  Short_t     field   = 5; // -5 for 2760GeV!!!
+  Bool_t      useCent = false;
+  const char* builder = "$(ALICE_ROOT)/PWG2/FORWARD/analysis2/trains/BuildTrain.C";
+  gROOT->LoadMacro(builder);
+  BuildTrain("MakeAODTrain");
+
+  std::cout << "making AOD train " << std::endl;
+  MakeAODTrain t(name, type, cms, field, useCent);
+  t.SetDataSet("");
+  t.SetNReplica(2);
+  t.SetAllowOverwrite(true);
+  t.SetProofServer(Form("workers=%d",proof));
+  
+  t.SetROOTVersion("v5-28-00f");
+  t.SetAliROOTVersion("v4-21-33-AN");
+
+  //if pp coll
+  if(type==1){
+    if (cms == 900) { 
+      if (!mc) Setup900GeV(t);
+      else     Setup900GeVMC(t);
+    }
+    if (cms == 2750) {
+      if (!mc) Setup2760GeV(t);
+      else Setup2760GeVMC(t);
+    }
+    
+    if (cms == 7000) { 
+      if (!mc) Setup7000GeV(t);
+      else     Setup7000GeVMC(t);
+    }
+  }
+  else
+    SetupPbPb(t);
+   
+    
+  t.Run("GRID", "FULL", nEvents, mc, usePar);
+}
+//________________________________________________________
+
+void SetupPbPb(MakeAODTrain& p){
+  p.SetDataDir("/alice/data/2010/LHC10h") ;
+  p.SetESDPass(2); 
+
+
+  p.AddRun(138190);
+  p.AddRun(138534);
+  p.AddRun(138364);
+  p.AddRun(138442);
+  //p.AddRun(139465);
+  p.AddRun(138396);
+  //p.AddRun(137722);
+  //p.AddRun(139107);
+  //p.AddRun(139437);
+  p.AddRun(138653);
+  //p.AddRun(139038);
+  //p.AddRun(137608);
+}
+
+void Setup900GeV(MakeAODTrain& p){
+  p.SetDataDir("/alice/data/2010/LHC10c");
+  p.SetESDPass(3); 
+
+  Add900GeVruns(p);
+
+}
+
+void Setup900GeVMC(MakeAODTrain& p){
+  //Flat mult
+  p.SetDataDir("/alice/sim/LHC10f1");
+
+  //nomal MC (Pythia)
+  //p.SetDataDir("/alice/sim/LHC11b1a");
+
+  //900 Gev normal Phojet MC
+  //p.SetDataDir("/alice/sim/LHC11c1");
+  
+  Add900GeVruns(p);
+
+}
+
+void Add900GeVruns(MakeAODTrain& p){
+  p.AddRun(118506);
+  p.AddRun(118507);
+  p.AddRun(118512);
+  p.AddRun(118556);
+  p.AddRun(118558);
+  p.AddRun(118560);
+  p.AddRun(118561);
+  p.AddRun(121039);
+  p.AddRun(121040);
+  
+}
+void Setup2760GeV(MakeAODTrain& p){
+  p.SetDataDir("/alice/data/2011/LHC11a");
+  p.SetESDPass(2);
+  p.SetPassPostfix("_with_SDD");
+  Add2760GeVruns(p);
+}
+void Setup2760GeVMC(MakeAODTrain& p){
+  p.SetDataDir("/alice/sim/LHC11b10c");
+  Add2760GeVruns(p);
+}
+
+void Add2760GeVruns(MakeAODTrain& p){
+  /*
+  p.AddRun(146860); 
+  p.AddRun(146859);    
+  p.AddRun(146858);    
+  p.AddRun(146857);    
+  p.AddRun(146856);
+  */
+  p.AddRun(146824); 
+  p.AddRun(146817);
+  p.AddRun(146807);
+  p.AddRun(146806);
+  p.AddRun(146805);
+  p.AddRun(146804);
+  p.AddRun(146803);
+  p.AddRun(146802);
+  p.AddRun(146801);
+  p.AddRun(146748);
+  p.AddRun(146747);
+  p.AddRun(146746);
+  p.AddRun(146689);
+  p.AddRun(146688);
+  p.AddRun(146686);
+  
+}
+
+void Setup7000GeV(MakeAODTrain& p){
+  p.SetDataDir("/alice/data/2010/LHC10d");
+  p.AddRun(126424);
+  p.AddRun(126422);  
+  p.AddRun(126409);
+  p.AddRun(126408);
+  p.AddRun(126407);
+  p.AddRun(126406);
+  p.AddRun(126404);
+  p.AddRun(126359);
+  p.AddRun(126352);
+  p.AddRun(126351);
+  p.AddRun(126097);
+  p.AddRun(125855);
+  p.AddRun(125851);
+  p.AddRun(125850);
+  p.AddRun(125849);
+  p.AddRun(125848);
+  p.AddRun(125847);
+  p.AddRun(125101);
+  p.AddRun(125097);
+  p.AddRun(125085);
+  
+  p.SetESDPass(2); 
+}
+void Setup7000GeVMC(MakeAODTrain& p){
+  p.SetDataDir("/alice/sim/LHC10h16");
+  p.AddRun(125186);
+  p.AddRun(125633);
+  p.AddRun(125842);
+  p.AddRun(125851);
+  p.AddRun(126008);
+  p.AddRun(126081);
+  p.AddRun(126082);
+  p.AddRun(126090);
+  p.AddRun(126097);
+  p.AddRun(126167);
+  p.AddRun(126283);
+  p.AddRun(126359);
+  p.AddRun(126408);
+  p.AddRun(126409);
+  p.AddRun(126422);
+  p.AddRun(126425);
+  p.AddRun(126437);
+}
+
+
+//EOF
diff --git a/PWGLF/FORWARD/analysis2/scripts/multdist/MakeMultiplicityDistributions.C b/PWGLF/FORWARD/analysis2/scripts/multdist/MakeMultiplicityDistributions.C
new file mode 100644 (file)
index 0000000..74f2f71
--- /dev/null
@@ -0,0 +1,93 @@
+/**
+ * @file 
+ * 
+ * @ingroup pwg2_forward_scripts
+ */
+
+/** 
+ * Run first pass of the analysis - that is read in ESD and produce AOD
+ * 
+ * @param esddir    ESD input directory. Any file matching the pattern 
+ *                  *AliESDs*.root are added to the chain 
+ * @param nEvents   Number of events to process.  If 0 or less, then 
+ *                  all events are analysed
+ * @param proof     Run in proof mode 
+ * @param mc        Run over MC data
+ *
+ * If PROOF mode is selected, then Terminate will be run on the master node 
+ * in any case. 
+ * 
+ *
+ * @ingroup pwg2_forward_scripts
+ */
+void MakeMultiplicityDistributions(const char* aoddir=".", 
+                                  Int_t       nEvents=-1, 
+                                  const char* trig="INEL",
+                                  Double_t    vzMin=-10,
+                                  Double_t    vzMax=10,
+                                  Int_t    lowCent= 0,
+                                  Int_t    highCent= 0,
+                                  char*       output="forward_multiplicity.root",         
+                                  Int_t       nBins= 500,
+                                  Int_t       proof=0)
+{
+  // --- Libraries to load -------------------------------------------
+  gROOT->Macro("$ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/LoadLibs.C");
+
+  // --- Check for proof mode, and possibly upload pars --------------
+  if (proof> 0) { 
+    gROOT->LoadMacro("$ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/LoadPars.C");
+    LoadPars(proof);
+  }
+  
+  // --- Our data chain ----------------------------------------------
+  std::cout << "making chain of files" << std::endl;
+  gROOT->LoadMacro("$ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/MakeChain.C");
+  TChain* chain = MakeChain("AOD", aoddir,true);
+  // If 0 or less events is select, choose all 
+  if (nEvents <= 0) nEvents = chain->GetEntries();
+
+  // --- Creating the manager and handlers ---------------------------
+  AliAnalysisManager *mgr  = new AliAnalysisManager("Forward Train", 
+                                                   "Forward Multiplicity");
+
+  if(lowCent>=highCent)
+    AliAnalysisManager::SetCommonFileName(output);
+  else{
+    Char_t* name = Form("%s_cent_%d_%d.root",output,lowCent, highCent);
+    AliAnalysisManager::SetCommonFileName(name);
+  }
+  // --- ESD input handler -------------------------------------------
+  AliAODInputHandler *aodInputHandler = new AliAODInputHandler();
+  mgr->SetInputEventHandler(aodInputHandler);      
+       
+  // --- Add tasks ---------------------------------------------------
+  // Forward 
+  gROOT->LoadMacro("$ALICE_ROOT/PWGLF/FORWARD/analysis2/AddTaskMultiplicity.C");
+  AddTaskMultiplicity(trig, vzMin, vzMax, lowCent, highCent, nBins);
+  
+  
+  // --- Run the analysis --------------------------------------------
+  TStopwatch t;
+  if (!mgr->InitAnalysis()) {
+    Error("MakeMult", "Failed to initialize analysis train!");
+    return;
+  }
+  // Skip terminate if we're so requested and not in Proof or full mode
+  mgr->SetSkipTerminate(false);
+  // Some informative output 
+  mgr->PrintStatus();
+  // mgr->SetDebugLevel(3);
+  if (mgr->GetDebugLevel() < 1 && !proof) 
+    mgr->SetUseProgressBar(kTRUE,100);
+
+  // Run the train 
+  t.Start();
+  Printf("=== RUNNING ANALYSIS ==================================");
+  mgr->StartAnalysis(proof ? "proof" : "local", chain, nEvents);
+  t.Stop();
+  t.Print();
+}
+//
+// EOF
+//
diff --git a/PWGLF/FORWARD/analysis2/scripts/multdist/doUnfolding.C b/PWGLF/FORWARD/analysis2/scripts/multdist/doUnfolding.C
new file mode 100644 (file)
index 0000000..83629ae
--- /dev/null
@@ -0,0 +1,15 @@
+{
+  gROOT->Macro("$ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/LoadLibs.C");
+  gSystem->Load("libPWGUDbase");
+  //gSystem->Load("libPWGmuondep");
+  gSystem->Load("libPWGUDselectors");
+  gSystem->Load("$HOME/Desktop/RooUnfold-1.0.3/libRooUnfold.so");
+  gSystem->AddIncludePath("-I${ALICE_ROOT}/include -I${ALICE_ROOT}/PWGPP/");
+
+  gROOT->LoadMacro("unfoldBase.C++g");
+  
+  unfoldBase("900GeV_materialCorrection_unfolded.root",kBayes, "/home/caz/ALICE/AliRoot.old/PWG2/FORWARD/analysis2/finalResponse_900GeV.root", "/home/caz/ALICE/thesis/analysisFiles/900GeV/900GeV_materialCorrection.root");
+  
+  
+}
+
diff --git a/PWGLF/FORWARD/analysis2/scripts/multdist/runMultiplicityDistributionAnalysis_README.txt b/PWGLF/FORWARD/analysis2/scripts/multdist/runMultiplicityDistributionAnalysis_README.txt
new file mode 100644 (file)
index 0000000..fc940aa
--- /dev/null
@@ -0,0 +1,112 @@
+****************************************************************
+* Instructions for running multiplicity distribution analysis  *
+*                    ------------                              *
+* Casper Nygaard, nygaard.casper@gmail.com, July 21st, 2012    *
+****************************************************************
+
+In these instructions the abbrevation $FMD corresponds to the
+$ALICE_ROOT/PWGLF/FORWARD/analysis2/ directory.
+
+For running the multiplicity distribution analysis the following files
+are used (in addition to some general AliRoot FORWARD code):
+
+$FMD/MakeMultAOD.C
+$FMD/MakeMultiplicityDistribution.C
+$FMD/CreateResponseMatrices.C
+$FMD/AliForwardMultiplicityDistribution.cxx/.h
+$FMD/AliForwardCreateResponseMatrices.cxx/.h
+$FMD/AddTaskMultiplicity.C
+$FMD/AddTaskCreateResponseMatrices.C
+
+$FMD/scripts/doUnfolding.C
+$FMD/scripts/unfoldBase.C
+$FMD/scripts/unfoldChi2Method.C
+
+Additionally the unfolding ROOT framework RooUnfold should be
+installed (http://hepunx.rl.ac.uk/~adye/software/unfold/RooUnfold.html).
+
+
+1) CREATE AODs
+------------------------------------------------------------
+Create AODs on the GRID, similarly to the standard
+FORWARD/analysis2 method. 
+
+One can use $FMD/MakeMultAOD.C macro to create the AODs, specifying
+which runs/ESD pass/etc. to use.
+
+
+2) CREATE BASIC MULTIPLICITY FILES
+------------------------------------------------------------
+Run over the (local) AODs to create basic multiplicity distribution
+analysis files, and response matrices.
+
+The basic idea is to run over data AODs to create raw multiplicity
+distrbutions in a number of eta bins, and run over MC AODs to create
+response matrices over the same eta bins.
+
+The output file is a .root file with a folder for each eta-bin,
+containing 'raw' multiplicity distributions or response matrices.
+
+
+Multiplicity distributions:
+aliroot -l 
+.L MakeMultiplicityDistribution.C
+MakeMultiplicityDistributions(const char* aoddir, Int_t nEvents, 
+                             const char* trig, 
+                             Double_t vzMin , Double_t vzMax, 
+                             Int_t lowCent, Int_t highCent, 
+                             char* output, Int_t nBins)
+
+aoddir   = directory of the local AOD files
+nEvents  = number of events to run over, -1 for all
+trig     = trigger word to analyse, "NSD" for Non-Single-Diffractive
+           or "INEL" for inelastic
+vzMin    = lower bound for z-vertex
+vzMax    = upper bound for z-vertex. To have full continuous
+           eta-coverage between the SPD and FMD: -4<vz<4
+lowCent  = lower bound for centrality
+highCent = upper bound for centrality. Only relevant for PbPb. In pp
+           set both to zero.
+output   = output filename
+nBins    = max multiplicity. The multiplicity axis runs from 
+          [-0.5 ; nBins-0.5] in nBins bins. For pp 500 suffices, for PbPb the most
+           central collisions need nBins=30000. 
+
+response matrices:
+
+aliroot -l
+.L CreateResponseMatrices.C
+CreateResponseMatrices(const char* aoddir, Int_t nEvents, 
+                      const char* trig, 
+                      Double_t vzMin, Double_t vzMax, 
+                      char* output")
+
+Parameters are the same as for
+MakeMultiplicityDistribution(). lowCent, highCent, nBins are not
+included since at the moment unfolding is not attempted for PbPb.
+
+
+3) UNFOLDING
+-------------------------------------------------------------
+The 'raw' distributions must be unfolded to correct for detector
+effects. For pp the unfolding can be done in a number of ways. I
+haveused two methods; a Single Value Decomposition (SVD) method and a
+Bayesen Iterative method. 
+
+The SVD method is implemented in AliRoot by Jan Fiete, and the
+Bayesian method uses the ROOT framework RooUnfold.
+
+doUnfolding.C is a small steering routine, which calls unfoldBase.C,
+which handles the unfolding itself. 
+
+aliroot -l
+.L unfoldBase.C
+unfoldBase(const Char_t* outputFile, Method method, 
+          const Char_t* responseFileName, const Char_t* dataFileName)
+
+outputFile        = name of unfolded distributions file
+method            = unfolding method, possible values are kBayes or kSvd
+responseFilenName = response matrix filename
+dataFileName      = uncorrected multiplicity filename
+
+
diff --git a/PWGLF/FORWARD/analysis2/scripts/multdist/unfoldBase.C b/PWGLF/FORWARD/analysis2/scripts/multdist/unfoldBase.C
new file mode 100644 (file)
index 0000000..f54f29f
--- /dev/null
@@ -0,0 +1,446 @@
+#ifndef __CINT__
+#include <iostream>
+#include <TH1.h>
+#include <TH2.h>
+#include <TF1.h>
+#include <TList.h>
+#include <TDirectory.h>
+#include <TFile.h>
+#include <TCanvas.h>
+#include <TPad.h>
+#include <THistPainter.h>
+#include <TObject.h>
+#include <TMath.h>
+#include <TFile.h>
+#include <TGraphErrors.h>
+#include <TLatex.h>
+#include <TLine.h>
+#include <TMarker.h>
+#include <TStyle.h>
+#include <TVirtualFitter.h>
+//#include <AliUnfolding.h>
+//#include <AliUnfolding.cxx>
+//#include <AliPWG0Helper.h>
+//#include <AliMultiplicityCorrection.h>
+#include <TMinuit.h>
+#include <TBox.h>
+#include <TGaxis.h>
+#include "TRandom.h"
+#include "/home/caz/ALICE/AliRoot/PWGLF/FORWARD/analysis2/AliForwardMultiplicityDistribution.h"
+#include "/home/caz/ALICE/AliRoot/PWGLF/FORWARD/analysis2/AliForwardMultiplicityDistribution.cxx"
+#include <TROOT.h>
+#include <TSystem.h>
+#include "TH1D.h"
+#include "/home/caz/Desktop/RooUnfold-1.0.3/src/RooUnfoldResponse.h"
+#include "/home/caz/Desktop/RooUnfold-1.0.3/src/RooUnfoldBayes.h"
+#include "unfoldChi2Method.C"
+//#include "RooUnfoldSvd.h"
+//#include "RooUnfoldBinByBin.h"
+#else
+class RooUnfoldResponse;
+class RooUnfoldBayes;
+class RooUnfold;
+class TF1;
+class TFile;
+class TH1D;
+class TH1;
+class TH1F;
+class TH2D;
+class TGraphErrors;
+class TGaxis;
+class AliUnfolding;
+#endif
+using namespace std;
+enum Method {  kBayes, kSvd };
+enum { kBla = 0, 
+       kPol0, 
+       kPol1, 
+       kLog, 
+       kEntropy, 
+       kCurvature, 
+       kRatio };
+
+void    createOneUnfold(TList* list, TFile* dataFile, TFile* out, Double_t l, Double_t h, Method method);
+TList*  getList(TFile* f, const Char_t* parent, const Char_t* name);
+const char* getPostfix(const TH1* h);
+void    unfoldBayesSet(const TObjArray& a, TDirectory* dir, TH2D* response, TH1D* triggerBias);
+void    unfoldChi2minSet(const TObjArray& a, TDirectory* dir, TH2D* response, TH1D* triggerBias, Int_t limit);
+void    normaliseAndCorrect(TH1D* unfolded, TH1D* triggerBias);
+void    normaliseAndCorrect(TH1D* unfolded, TH1D* triggerBias, Int_t limit);
+
+TH1D*   unfoldBayes(TH1D* data, RooUnfoldResponse* response);
+TH1D*   unfoldChi2min(TH1F* data, TH2D* response, TH1F* eff, Int_t regFunc, Float_t regWeight);
+TH1D*   getTriggerBiasHistogram(TList* responseList);
+void    doTriggerBiasCorrection(TH1D*& hist, TH1D* triggerBias);
+
+// --- Steering routine ----------------------------------------------
+void unfoldBase(const Char_t* outputFile="unfoldOutput.root", Method method= kBayes,
+                 const Char_t* responseFileName="$ALICE_ROOT/PWGLF/FORWARD/analysis2/responseMatrices.root", 
+         const Char_t* dataFileName="$ALICE_ROOT/PWGLF/FORWARD/analysis2/forward_multiplicity.root"){
+
+  /*#ifdef __CINT__
+  gSystem->Load("libANALYSIS");
+  gSystem->Load("libANALYSISalice");
+  
+  gSystem->Load("libPWG0base");
+  gSystem->Load("libPWG0dep");
+  gSystem->Load("libPWG0selectors");
+  
+  gROOT->Macro("$ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/LoadLibs.C");
+  gROOT->LoadMacro("/home/caz/ppData/FirstAnalysis/multDists/Unfold.C");
+#endif
+  */
+  TFile* responseFile = TFile::Open(responseFileName,"READ");
+  if (!responseFile) { 
+    Error("createUnfoldingFile", "Couldn't open %s", responseFileName);
+    return;
+  }
+  TList* list = static_cast<TList*>(responseFile->Get("ResponseMatrices"));
+  TFile* dataFile = TFile::Open(dataFileName,"READ");
+  if (!dataFile) { 
+    Error("createUnfoldingFile", "Couldn't open %s", dataFileName);
+    return;
+  }
+  TFile* out = TFile::Open(outputFile, "RECREATE");
+  cout << " Loaded responseFile and dataFile" << endl;
+
+  
+  Double_t limits[] = { 2.4, 1, 0. };
+  // Double_t limits[] = {5.1, 3.4, 3.0, 2.4, 2.0, 1.5, 1.0, 0.5, 0. };
+  Double_t* limit = limits;
+  while ((*limit) > 0.1){
+    if((*limit) >5)
+      createOneUnfold(list, dataFile, out, -3.4,(*limit), method);
+    else{
+      createOneUnfold(list, dataFile, out, -(*limit),(*limit), method);
+    }
+    limit++; 
+  }
+  
+  /*
+  limit = limits;
+  while ((*limit) > 0.1) { 
+    createOneUnfold(list, dataFile, out, 0,(*limit),method);
+    createOneUnfold(list, dataFile, out, -(*limit), 0, method);
+    limit++; }
+  
+  
+  for (Double_t l = -3; l < -0.5; l += 0.5){ 
+    createOneUnfold(list, dataFile, out, l, l+0.5, method);  }
+  for (Double_t l = 0.5; l < 5; l += 0.5){ 
+    createOneUnfold(list, dataFile, out, l, l+0.5, method);  }
+  for (Double_t l = -3; l < 5; l += 0.2){ 
+    createOneUnfold(list, dataFile, out, l, l+0.2, method);  }
+  */
+    
+  out->Close();
+  responseFile->Close();
+  dataFile->Close();
+}
+
+
+// --- Unfold data for one eta bin -----------------------------------
+void createOneUnfold(TList* list, TFile* dataFile, TFile* out, Double_t l, Double_t h, Method method)
+{
+  const char* name= AliForwardMultiplicityDistribution::FormBinName(l, h);
+  std::cout << "Name=" << name << std::endl;
+  TDirectory* dir = out->mkdir(name);
+  dir->cd();
+  
+  TList* multList = getList(dataFile, "CentralSums", name);
+  if (!multList) return;
+  
+  //_____if unfolding MC response files____________________
+  //TList* multList = getList(dataFile, "ResponseMatrices", name);
+  //if (!multList) return;
+  //cout << multList << endl;
+    
+  Int_t maxMult;
+  TH1D* tmp= (TH1D*)multList->FindObject("mult");
+  for(Int_t n=5; n<=tmp->GetNbinsX();n++){
+    if(tmp->GetBinContent(n)<1e-9){
+      maxMult= n;
+      break;
+    }
+  }
+
+  
+  TObjArray a;
+  a.Add(multList->FindObject("mult"));
+  a.Add(multList->FindObject("multPlusSys"));
+  a.Add(multList->FindObject("multMinusSys"));
+  
+  TList* responseList = 0;
+  responseList = static_cast<TList*>(list->FindObject(name));
+    
+  TH2D* response= (TH2D*)responseList->FindObject("responseMatrix");
+  response->SetDirectory(0);
+  
+  TH1D* triggerBias = getTriggerBiasHistogram(responseList);
+  
+  switch (method) { 
+  case kBayes: 
+    unfoldBayesSet(a, dir, response, triggerBias);
+    break;
+  case kSvd: 
+     unfoldChi2minSet(a, dir, response, triggerBias, maxMult+30);
+     break;
+  }
+  
+  
+  TIter next(&a);
+  TH1D* data = 0;
+  while ((data = static_cast<TH1D*>(next()))) {
+    normaliseAndCorrect(data,triggerBias);
+  }
+  
+  dir->Add(multList->FindObject("multMC"));
+  
+  dir->Add(&a);
+  dir->Write();
+  
+  //delete list;
+  //delete responseList;
+  
+  //delete triggerBias;
+  delete multList;
+  delete dir;
+}
+
+// --- Get a list from a file ----------------------------------------
+TList*
+getList(TFile* f, const Char_t* parent, const Char_t* name)
+{
+  TList* list = static_cast<TList*>(f->Get(parent));
+  if (!list) { 
+    Error("getList", "Couldn't find list %s in %s", parent, f->GetName());
+    return 0;
+  }
+  if (!name) return list;
+  TList* ret = static_cast<TList*>(list->FindObject(name));
+  if (!ret) { 
+    Error("getList", "Couldn't find %s in %s", name, parent);
+    return 0;
+  }
+  
+  return ret;
+}
+
+// --- Get end of input histogram name -------------------------------
+const Char_t* 
+getPostfix(const TH1* h)
+{
+  static TString t;
+  t = h->GetName();
+  t.ReplaceAll("mult", "");
+  return t.Data();
+}
+
+// --- Unfold a single hist using Bayes ------------------------------
+void
+unfoldBayesSet(const TObjArray& a, TDirectory* dir, TH2D* response, TH1D* triggerBias)
+{
+  TH1* projY= (TH1D*) response->ProjectionY("projY",1,response->GetNbinsX(),"");
+  TH1* projX= (TH1D*) response->ProjectionX("projX",1,response->GetNbinsY(),"");
+  projX->SetDirectory(0);
+  projY->SetDirectory(0);
+  
+  TH2D* responseTranspose = (TH2D*)response->Clone("response");
+  for(Int_t i=1;i<=response->GetNbinsX();i++){
+    for(Int_t j=1;j<=response->GetNbinsY();j++){
+      responseTranspose->SetBinContent(j,i, response->GetBinContent(i,j));
+      responseTranspose->SetBinError(j,i, response->GetBinError(i,j));
+    }
+  }
+  RooUnfoldResponse* responseObject = new RooUnfoldResponse(projY, projX, responseTranspose,"test", "bla");
+
+  Int_t dC = gStyle->GetNumberOfColors() / a.GetEntriesFast(); 
+  Int_t iC = 0;
+  TIter next(&a);
+  TH1D* h = 0;
+  while ((h = static_cast<TH1D*>(next()))) {
+    TH1D* unfolded = unfoldBayes(h, responseObject);
+    normaliseAndCorrect(unfolded,triggerBias);
+    unfolded->SetDirectory(dir);
+    unfolded->SetMarkerColor(gStyle->GetColorPalette(iC += dC));
+    unfolded->SetMarkerStyle(20);
+  }
+  
+  delete projY;
+  delete projX;
+  delete response;
+  delete responseTranspose;
+  delete responseObject;
+}
+
+// --- Unfold a single hist using Jan-Fietes SVD (Chi2) method --------------------------
+void
+unfoldChi2minSet(const TObjArray& a, TDirectory* dir, TH2D* response, TH1D* triggerBias, Int_t limit)
+{
+
+  TIter next(&a);
+  TH1D* data = 0;
+  while ((data = static_cast<TH1D*>(next()))) {
+    Int_t max = (limit < 0 ? data->GetNbinsX() : limit);
+    TH1F* dataTmp     = new TH1F(data->GetName(),data->GetTitle(), max,-0.5, max-.5);
+    dataTmp->SetDirectory(0);
+    for(Int_t k=1;k<= max; k++){
+      dataTmp->SetBinContent(k,data->GetBinContent(k));
+      dataTmp->SetBinError(k, data->GetBinError(k));
+    }
+    
+    TH2D* responseTmp = response; 
+    TH1F* eff         = new TH1F("eff","unfolded",max, -.5, max-.5);
+    if (limit > 0) {
+      responseTmp = new TH2D("res","res", max, -.5, max-.5,max, -.5, max-.5);
+      for(Int_t k=1;k<= max; k++){
+       eff->SetBinContent(k,1);
+       if (limit > 0) {
+         for(Int_t j=1;j<=max; j++){
+           responseTmp->SetBinContent(k,j,response->GetBinContent(k,j));
+           responseTmp->SetBinError(k,j, response->GetBinError(k,j));
+         }
+       }
+      }    
+    }
+    Int_t nR = 8;
+    Int_t dC = gStyle->GetNumberOfColors() / (nR * a.GetEntriesFast()); 
+    Int_t iC = 0;
+    for(Int_t f = 3; f <= 3; f++){
+      iC = 0;
+      for(Float_t w = 1e-5; w <= 1e5 ;w *= 10){
+       TH1D* unfolded = unfoldChi2min(dataTmp,responseTmp, eff, f, w);
+       normaliseAndCorrect(unfolded,triggerBias, limit);
+       unfolded->SetDirectory(dir);
+       unfolded->SetMarkerColor(gStyle->GetColorPalette(iC += dC));
+       unfolded->SetMarkerStyle(20);
+      }
+    }
+    if (limit > 0) delete responseTmp;
+    delete dataTmp;
+    delete eff;
+  }
+}
+
+// --- Normalise to 1 and correct for trigger ------------------------
+void normaliseAndCorrect(TH1D* unfolded, TH1D* triggerBias)
+{
+  if (triggerBias) doTriggerBiasCorrection(unfolded, triggerBias);
+  unfolded->SetDirectory(0);
+  unfolded->Sumw2();
+  unfolded->Scale(1/unfolded->Integral());
+}
+
+void normaliseAndCorrect(TH1D* unfolded, TH1D* triggerBias, Int_t limit)
+{
+  if (triggerBias) doTriggerBiasCorrection(unfolded, triggerBias);
+  unfolded->SetDirectory(0);
+  unfolded->Sumw2();
+  unfolded->Scale(1/unfolded->Integral(1,limit));
+}
+
+// --- Use RooUnfold stuff -------------------------------------------
+TH1D* unfoldBayes(TH1D* data, RooUnfoldResponse* response)
+{
+  
+  
+  RooUnfold* unfold = new RooUnfoldBayes(response, data, 20);
+  TH1D* unfolded= (TH1D*) unfold->Hreco();
+  
+  TString t = TString::Format("unfolded_bayes_%s",getPostfix(data));
+  
+  unfolded->SetName(t.Data());
+  delete unfold;
+
+  return unfolded;
+}
+
+// --- Use Jan-Fiete routines ----------------------------------------
+TH1D* unfoldChi2min(TH1F* data, TH2D* response, TH1F* eff, 
+                   Int_t regFunc, Float_t regWeight)
+{  
+  TH1F* tmp = static_cast<TH1F*>(data->Clone("tmp"));
+  tmp->Reset();
+  UnfoldChi2Min(data, tmp, response, eff, regFunc, regWeight);  
+  
+  Int_t max = data->GetNbinsX();
+  TH1D* unfolded= new TH1D("unfolded","unfolded", max, -.5, max);  
+  for(Int_t n=1;n<=max; n++){
+    unfolded->SetBinContent(n,tmp->GetBinContent(n));
+    unfolded->SetBinError(n, tmp->GetBinError(n));
+  }
+
+  TString regString="0";
+  switch (regFunc) { 
+  case 1: regString="pol0";      break;
+  case 2: regString="pol1";      break;
+  case 3: regString="log";       break;
+  case 4: regString="entropy";   break;
+  case 5: regString="curvature"; break;
+  default:
+    Error("unfold", "Invalid regFunc type");
+    return 0;
+  }
+  delete tmp; 
+
+  Info("", "data name: %s, postfix: %s", data->GetName(), getPostfix(data));
+  TString n = TString::Format("unfolded_chi2_%s_%s_%1.0e", getPostfix(data),regString.Data(),regWeight);
+  n.ReplaceAll("-0", "m");
+  n.ReplaceAll("+0", "p");
+  //n.ReplaceAll(".", "d");
+  unfolded->SetName(n);
+  return unfolded;
+}
+
+
+// --- Find trigger bias histogram -----------------------------------
+TH1D* getTriggerBiasHistogram(TList* responseList){
+  TH1D* hMCNSD;
+  TH1D* hMCESDNSD;
+  TH1D* hESDNSD;
+  hMCNSD = (TH1D*)responseList->FindObject("fMCNSD");
+  hMCESDNSD = (TH1D*)responseList->FindObject("fMCESDNSD");
+  hESDNSD = (TH1D*)responseList->FindObject("fESDNSD");
+    
+  hMCNSD->Sumw2();
+  hESDNSD->Sumw2();
+  hMCESDNSD->Sumw2();
+  
+  TH1D* corr     = new TH1D("corr", "corr", 50,-0.5,49.5);
+  for(Int_t n=1;n<corr->GetNbinsX();n++){
+    Double_t errorSquared=0;
+    errorSquared= (1/hESDNSD->GetBinContent(n)+   1/hMCNSD->GetBinContent(n));   
+    corr->SetBinContent(n, hESDNSD->GetBinContent(n)/hMCNSD->GetBinContent(n));
+    corr->SetBinError(n,corr->GetBinContent(n)*TMath::Sqrt(errorSquared));
+  }
+  
+  hMCESDNSD->Divide(hMCNSD);
+  hESDNSD->Divide(hMCNSD);
+   
+  
+  delete hMCNSD;
+  delete hMCESDNSD;
+  
+  return corr;
+    
+}
+
+// --- Apply trigger bias histogram ----------------------------------
+void doTriggerBiasCorrection(TH1D*& hist, TH1D* triggerBias){
+  for(Int_t i = 1; i<=35;i++){
+    if(triggerBias->GetBinContent(i)>1e-5&&hist->GetBinContent(i)>0){
+      hist->SetBinContent(i, hist->GetBinContent(i)/triggerBias->GetBinContent(i));
+      hist->SetBinError(i,TMath::Sqrt(TMath::Power(hist->GetBinError(i)/hist->GetBinContent(i),2)+TMath::Power(triggerBias->GetBinError(i)/triggerBias->GetBinContent(i),2))*hist->GetBinContent(i));
+      //cout << hist->GetBinContent(i) << " +-  " << hist->GetBinError(i) << endl; 
+    }
+    //if(triggerBias->GetBinContent(i)<1e-5)
+    //  hist->SetBinContent(i, hist->GetBinContent(i)/triggerBias->GetBinContent(i+1));
+    
+  }
+}
+
+//
+// EOF
+//
diff --git a/PWGLF/FORWARD/analysis2/scripts/multdist/unfoldChi2Method.C b/PWGLF/FORWARD/analysis2/scripts/multdist/unfoldChi2Method.C
new file mode 100644 (file)
index 0000000..9ccb1db
--- /dev/null
@@ -0,0 +1,80 @@
+#ifndef __CINT__
+#include <iostream>
+#include <TH1.h>
+#include <TH2.h>
+#include <TF1.h>
+#include <TList.h>
+#include <TDirectory.h>
+#include <TFile.h>
+#include <TCanvas.h>
+#include <TPad.h>
+#include <THistPainter.h>
+#include <TObject.h>
+#include <TMath.h>
+#include <TFile.h>
+#include <TGraphErrors.h>
+#include <TLatex.h>
+#include <TLine.h>
+#include <TMarker.h>
+#include <TStyle.h>
+#include <TVirtualFitter.h>
+#include "/home/caz/ALICE/AliRoot/ANALYSIS/AliUnfolding.h"
+#include "/home/caz/ALICE/AliRoot/ANALYSIS/AliUnfolding.cxx"
+#include <AliPWG0Helper.h>
+// #include <AliMultiplicityCorrection.h>
+//#include <TMinuit.h>
+#include <TBox.h>
+#include <TGaxis.h>
+#else
+class TF1;
+class TH1D;
+class TH1F;
+class TGraphErrors;
+class TGaxis;
+// class AliUnfolding;
+#endif
+using namespace std;
+
+
+void UnfoldChi2Min(TH1F*   data, 
+                  TH1F*&  unfolded, 
+                  TH2D*   response=0 ,
+                  TH1F*   eff =0,
+                  Int_t   regularization=1,
+                  Float_t   regularizationWeight=1000,
+                  Float_t minInitialValue=0.0001, 
+                  Bool_t  skipBin0InChi2= kFALSE,
+                  Int_t   skipNbins=0 ){
+#ifdef __CINT__
+  gSystem->Load("libANALYSIS");
+  gSystem->Load("libANALYSISalice");
+  gSystem->Load("libPWG0base");
+  gSystem->Load("libPWG0dep");
+  gSystem->Load("libPWG0selectors");
+#endif
+
+  cout << data << "  " << unfolded << endl;
+  
+  AliUnfolding::RegularizationType type = AliUnfolding::kNone;
+  switch (regularization) {
+  case 0: type = AliUnfolding::kNone; break; 
+  case 1: type = AliUnfolding::kPol0; break; 
+  case 2: type = AliUnfolding::kPol1; break; 
+  case 3: type = AliUnfolding::kLog; break; 
+  case 4: type = AliUnfolding::kEntropy; break; 
+  case 5: type = AliUnfolding::kCurvature; break; 
+  case 6: type = AliUnfolding::kRatio; break; 
+  };
+
+  AliUnfolding::SetUnfoldingMethod(AliUnfolding::kChi2Minimization);
+  AliUnfolding::SetNbins(data->GetNbinsX(), unfolded->GetNbinsX()); 
+  AliUnfolding::SetMinimumInitialValue(kTRUE, minInitialValue);
+  AliUnfolding::SetSkip0BinInChi2(skipBin0InChi2);
+  AliUnfolding::SetSkipBinsBegin(skipNbins);
+  AliUnfolding::SetChi2Regularization(type, regularizationWeight);
+  //AliUnfolding::SetCreateOverflowBin(1e-6);
+    
+  AliUnfolding::Unfold(response, eff, data, data, unfolded, kFALSE);
+}
index 6ef6ba19bf56e94b6bbb6cef67fae45aa31ff22c..be56af2c1874755b6af94d1fbc002908a6d16cd9 100644 (file)
 #pragma link C++ class AliSPDMCTrackDensity+;
 
 
+#pragma link C++ class AliDisplacedVerticesESDFilterTask+;
+
+#pragma link C++ class AliForwardCreateResponseMatrices+;
+#pragma link C++ class AliForwardMultiplicityDistribution+;
+
 #else
 # error Not for compilation 
 #endif