[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliForwardMultiplicityDistribution.cxx


#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->CreateOutputObjectss(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::CreateOutputObjectss(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
Commit	Line	Data
ad7be237	1
	2	#include <TH1D.h>
	3	#include "AliForwardMultiplicityDistribution.h"
	4	#include "AliAODForwardMult.h"
	5	#include "AliAODCentralMult.h"
	6	#include "AliAODEvent.h"
	7	#include "TString.h"
	8	using namespace std;
	9
	10	ClassImp(AliForwardMultiplicityDistribution)
	11	#if 0
	12	; // This is for Emacs - do not delete
	13	#endif
	14	//______________________________________________________________________
	15	AliForwardMultiplicityDistribution::AliForwardMultiplicityDistribution()
	16	: AliBasedNdetaTask(),
	17	fTrigger(0), // trigger histogram
	18	fBins(), // eta bin list
	19	fOutput(0), // output list
	20	fLowCent(-1), // lower centrality limit
	21	fHighCent(-1), // upper centrality limit
	22	fNBins(-1), // multiplicity axis' runs from 0 to fNbins
	23	fCent(0) // centrality
	24	{
	25	//
	26	// Default Constructor
	27	//
	28	}
	29
	30	//_____________________________________________________________________
	31	AliForwardMultiplicityDistribution::AliForwardMultiplicityDistribution(const char* name)
	32	: AliBasedNdetaTask(name),
	33	fTrigger(0), // trigger histogram
	34	fBins(), // eta bin list
	35	fOutput(0), // output list
	36	fLowCent(-1), // lower centrality limit
	37	fHighCent(-1), // upper centrality limit
	38	fNBins(-1), // multiplicity axis' runs from 0 to fNbins
	39	fCent(0) // centrality
	40	{
	41	//
	42	// Constructor
	43	//
	44	DefineOutput(1, TList::Class());
	45	}
	46
	47	//_____________________________________________________________________
	48	void AliForwardMultiplicityDistribution::UserCreateOutputObjects()
	49	{
	50	//
	51	// Create Output Objects
	52	//
	53	fOutput = new TList;
	54	fOutput->SetOwner();
	55	fOutput->SetName(GetName());
	56
	57	TH1D* dndetaSumForward = new TH1D("dndetaSumForward","dndetaSumForward", 200,-4,6);
	58	TH1D* dndetaSumCentral = new TH1D("dndetaSumCentral","dndetaSumCentral", 200,-4,6);
	59	fOutput->Add(dndetaSumForward);
	60	fOutput->Add(dndetaSumCentral);
	61
	62	fCent = new TH1D("cent", "Centrality", 100, 0, 100);
	63	fCent->SetDirectory(0);
	64	fCent->SetXTitle(0);
65	fOutput->Add(fCent);
66
67	TH1D* dndetaEventForward = new TH1D();
68	TH1D* dndetaEventCentral = new TH1D();
69	fOutput->Add(dndetaEventForward);
70	fOutput->Add(dndetaEventCentral);
71
72	//make trigger diagnostics histogram
73	fTrigger= new TH1I();
74	fTrigger= AliAODForwardMult::MakeTriggerHistogram("triggerHist");
75	fOutput->Add(fTrigger);
76
77	//Loop over all individual eta bins, and define their hisograms
78	TIter next(&fBins);
79	Bin * bin = 0;
80	while ((bin = static_cast<Bin*>(next()))) {
5934a3e3	81	bin->CreateOutputObjectss(fOutput, fNBins);
ad7be237	82	}
	83	//fOutput->ls();
	84	PostData(1, fOutput);
	85	}
	86
	87
	88	//_____________________________________________________________________
	89	void AliForwardMultiplicityDistribution::UserExec(Option_t /option*/)
	90	{
	91	//
	92	// User Exec
	93	//
	94	AliAODEvent* aod = dynamic_cast<AliAODEvent*>(InputEvent());
	95	if (!aod) {
	96	AliError("Cannot get the AOD event");
	97	return;
	98	}
	99	// Check AOD event
	100	AliAODForwardMult* AODforward = static_cast<AliAODForwardMult*>(aod->FindListObject("Forward"));
	101	if (!AODforward->CheckEvent(fTriggerMask, fVtxMin, fVtxMax, 0,0,fTrigger)){
	102	AliError("Invalid Event");
	103	return;
	104	}
	105
	106	// Fill centrality histogram, only useful for PbPb
	107	Float_t cent = AODforward->GetCentrality();
	108	if(fLowCent<fHighCent){
	109	if(cent<fLowCent\|\|cent>fHighCent) return;
	110	}
	111	fCent->Fill(cent);
	112
	113	TH2D forward;
	114	TH2D central;
	115
	116	// Get 2D eta-phi histograms for each event
	117	GetHistograms(aod, forward, central);
	118
	119	TH1D* dndetaSumForward = (TH1D*)fOutput->FindObject("dndetaSumForward");
	120	TH1D* dndetaSumCentral = (TH1D*)fOutput->FindObject("dndetaSumCentral");
	121	TH1D* dndetaEventForward = (TH1D*)fOutput->FindObject("dndetaEventForward");
	122	TH1D* dndetaEventCentral = (TH1D*)fOutput->FindObject("dndetaEventCentral");
	123	TH1D* normEventForward = 0;
	124	TH1D* normEventCentral = 0;
	125
	126	dndetaEventForward = forward.ProjectionX("dndetaForward",1,forward.GetNbinsY(),"");
	127	dndetaEventCentral = central.ProjectionX("dndetaCentral",1,central.GetNbinsY(),"");
	128	dndetaSumForward->Add(dndetaEventForward);
	129	dndetaSumCentral->Add(dndetaEventCentral);
	130
	131	// underflow eta bin of forward/central carry information on whether there is acceptance. 1= acceptance, 0= no acceptance
	132	normEventForward = forward.ProjectionX("dndetaEventForwardNorm",0,0,"");
	133	normEventCentral = central.ProjectionX("dndetaEventCentralNorm",0,0,"");
	134
	135	Double_t VtxZ = AODforward->GetIpZ();
	136
	137
	138	// loop over all eta bins, and fill multiplicity histos
	139	TIter next(&fBins);
	140	Bin * bin = 0;
	141	while ((bin = static_cast<Bin*>(next()))) {
	142	bin->Process(dndetaEventForward, dndetaEventCentral,
	143	normEventForward, normEventCentral, VtxZ);
	144	}
	145
146	PostData(1, fOutput);
147
148	}
149
150	//_____________________________________________________________________
151	void AliForwardMultiplicityDistribution::Terminate(Option_t /option*/)
152	{
153	//
154	// Terminate
155	//
156	}
157
158	//_________________________________________________________________-
159	TH2D* AliForwardMultiplicityDistribution::GetHistogram(const AliAODEvent*, Bool_t)
160	{
161	//
162	// implementation of pure virtual function, always returning 0
163	//
164	return 0;
165	}
166
167	void AliForwardMultiplicityDistribution::GetHistograms(const AliAODEvent* aod, TH2D& forward, TH2D& central, Bool_t /mc/)
168	{
169	//
170	// Get single event forward and central dN²/dEta dPhi histograms
171	//
172	TObject* forwardObj = 0;
173	TObject* centralObj = 0;
174
175
176	forwardObj = aod->FindListObject("Forward");
177	centralObj = aod->FindListObject("CentralClusters");
178
179	AliAODForwardMult* AODforward = static_cast<AliAODForwardMult*>(forwardObj);
180	AliAODCentralMult* AODcentral = static_cast<AliAODCentralMult*>(centralObj);
181
182	forward= AODforward->GetHistogram();
183	central= AODcentral->GetHistogram();
184
185	}
186
187	//______________________________________________________________________
188	#if 1
189	const Char_t*
190	#else
191	TString
192	#endif
193	AliForwardMultiplicityDistribution::FormBinName(Double_t etaLow, Double_t etaHigh)
194	{
195	//
196	// Form name of eta bin
197	//
198	TString sLow(TString::Format("%+03d",int(10*etaLow)));
199	TString sHigh(TString::Format("%+03d",int(10*etaHigh)));
200	sLow.ReplaceAll("+", "plus");
201	sLow.ReplaceAll("-", "minus");
202	sHigh.ReplaceAll("+", "plus");
203	sHigh.ReplaceAll("-", "minus");
204	#if 0
205	TString* combined= new TString();
206	combined->Append(TString::Format("%s_%s", sLow.Data(), sHigh.Data()));
207	return combined->Data();
208	#else
209	static TString tmp;
210	tmp = TString::Format("%s_%s", sLow.Data(), sHigh.Data());
211	return tmp.Data();
212	#endif
213	}
214
215	//=====================================================================
216	AliForwardMultiplicityDistribution::Bin::Bin()
217	: TNamed(),
218	fEtaLow(-1111), // low eta limit
219	fEtaHigh(-1111), // high eta limit
220	fHist(0), // multiplicity distribution hist
221	fHistPlus05(0), // multiplicity distribution hist scaled up with 5%
222	fHistPlus075(0), // multiplicity distribution hist scaled up with 7.5%
223	fHistPlus10(0), // multiplicity distribution hist scaled up with 10%
224	fHistMinus05(0), // multiplicity distribution hist scaled down with 5%
225	fHistMinus075(0), // multiplicity distribution hist scaled down with 7.5%
226	fHistMinus10(0), // multiplicity distribution hist scaled down with 10%
227	fHistPlusSys(0), // multiplicity distribution hist scaled up with the event uncertainty
228	fHistMinusSys(0), // multiplicity distribution hist scaled down with the event uncertainty
229	fAcceptance(0), // histogram showing the 'holes' in acceptance.
230	fVtxZvsNdataBins(0) // VtxZ vs. number of data acceptance bins (normalised to the eta range)
231	{
232	//
233	// Default Constructor
234	//
235	}
236	//_____________________________________________________________________
237	AliForwardMultiplicityDistribution::Bin::Bin(Double_t etaLow, Double_t etaHigh)
238	: TNamed("", ""),
239	fEtaLow(etaLow), // low eta limit
240	fEtaHigh(etaHigh), // high eta limit
241	fHist(0), // multiplicity distribution hist
242	fHistPlus05(0), // multiplicity distribution hist scaled up with 5%
243	fHistPlus075(0), // multiplicity distribution hist scaled up with 7.5%
244	fHistPlus10(0), // multiplicity distribution hist scaled up with 10%
245	fHistMinus05(0), // multiplicity distribution hist scaled down with 5%
246	fHistMinus075(0), // multiplicity distribution hist scaled down with 7.5%
247	fHistMinus10(0), // multiplicity distribution hist scaled down with 10%
248	fHistPlusSys(0), // multiplicity distribution hist scaled up with the event uncertainty
249	fHistMinusSys(0), // multiplicity distribution hist scaled down with the event uncertainty
250	fAcceptance(0), // histogram showing the 'holes' in acceptance.
251	fVtxZvsNdataBins(0) // VtxZ vs. number of data acceptance bins (normalised to the eta range)
252	{
253	//
254	// Constructor
255	//
256	const char* name = AliForwardMultiplicityDistribution::FormBinName(fEtaLow,fEtaHigh);
257
258	SetName(name);
259	SetTitle(TString::Format("%+4.1f < #eta < %+4.1f", fEtaLow, fEtaHigh));
260
261	}
262	//_____________________________________________________________________
5934a3e3	263	void AliForwardMultiplicityDistribution::Bin::CreateOutputObjectss(TList* cont, Int_t max)
ad7be237	264	{
	265	//
	266	// Define eta bin output histograms
	267	//
	268	TList* out = new TList;
	269	out->SetName(GetName());
	270	cont->Add(out);
	271
	272	fHist = new TH1D("mult", GetTitle(), max, -0.5, max-.5);
	273	fHistPlus05 = new TH1D("multPlus05", GetTitle(), max, -0.5, max-.5);
	274	fHistPlus075 = new TH1D("multPlus075", GetTitle(), max, -0.5, max-.5);
	275	fHistPlus10 = new TH1D("multPlus10", GetTitle(), max, -0.5, max-.5);
	276	fHistMinus05 = new TH1D("multMinus05", GetTitle(), max, -0.5, max-.5);
	277	fHistMinus075 = new TH1D("multMinus075", GetTitle(), max, -0.5, max-.5);
	278	fHistMinus10 = new TH1D("multMinus10", GetTitle(), max, -0.5, max-.5);
	279	fHistPlusSys = new TH1D("multPlusSys", GetTitle(), max, -0.5, max-.5);
	280	fHistMinusSys = new TH1D("multMinusSys", GetTitle(), max, -0.5, max-.5);
	281	fVtxZvsNdataBins = new TH2D("VtxZvsNdataBins", "VtxZ vs dataAcceptance/etaRange;z-vtz;dataAcceptance/etaRange", 20, -10,10, 130,0,1.3);
	282	fAcceptance = new TH2D("Acceptance","Acceptance;#eta;z-vtx",200,-4, 6 , 20,-10,10);
	283
	284	out->Add(fHist);
	285	out->Add(fHistPlus05);
	286	out->Add(fHistPlus075);
	287	out->Add(fHistPlus10);
	288	out->Add(fHistMinus05);
	289	out->Add(fHistMinus075);
	290	out->Add(fHistMinus10);
	291	out->Add(fHistPlusSys);
	292	out->Add(fHistMinusSys);
	293	out->Add(fVtxZvsNdataBins);
	294	out->Add(fAcceptance);
	295
	296	}
	297
	298	//_____________________________________________________________________
	299	void AliForwardMultiplicityDistribution::Bin::Process(TH1D* dndetaForward, TH1D* dndetaCentral,
	300	TH1D* normForward, TH1D* normCentral, Double_t VtxZ)
	301	{
	302	//
	303	// Process single eta bin
	304	//
	305
	306	// Diagnostics on event acceptance
	307	Int_t first = dndetaForward->GetXaxis()->FindBin(fEtaLow);
	308	Int_t last = dndetaForward->GetXaxis()->FindBin(fEtaHigh-.0001);
	309	Double_t acceptanceBins=0;
	310	for(Int_t n= first;n<=last;n++){
	311	if(normForward->GetBinContent(n)>0\|\|normCentral->GetBinContent(n)>0){
	312	acceptanceBins++;
	313	}
	314	fAcceptance->SetBinContent(n, fAcceptance->GetYaxis()->FindBin(VtxZ), 1);
	315	if(normForward->GetBinContent(n)>0\|\|normCentral->GetBinContent(n)>0)
	316	fAcceptance->SetBinContent(n, fAcceptance->GetYaxis()->FindBin(VtxZ),10);
	317	}
	318
	319	//std::cout << "% of bins covered in eta-range: " << acceptanceBins/(last-first+1) << std::endl;
	320	fVtxZvsNdataBins->Fill(VtxZ, (Double_t)acceptanceBins/(last-first+1));
	321
	322
	323	Double_t c = 0;
	324	Double_t e2 = 0;
	325	for (Int_t i = first; i <= last; i++) {
	326	Double_t cForward = 0;
	327	Double_t cCentral = 0;
328	Double_t e2Forward = 0;
329	Double_t e2Central = 0;
330	if (normForward->GetBinContent(i) > 0) {
331	cForward = dndetaForward->GetBinContent(i);
332	e2Forward = dndetaForward->GetBinError(i) * dndetaForward->GetBinError(i);
333	}
334	if (normCentral->GetBinContent(i) > 0) {
335	cCentral = dndetaCentral->GetBinContent(i);
336	e2Central = dndetaCentral->GetBinError(i) * dndetaCentral->GetBinError(i);
337	}
338	Double_t cc = 0;
339	Double_t ee2 = 0;
340	// if overlap between central/forward, use average of the two
341	if (cCentral > 0 && cForward > 0) {
342	cc = 0.5 * (cForward + cCentral);
343	ee2 = 0.5 * (e2Forward + e2Central);
344	}
345	else if (cCentral > 0) {
346	cc = cCentral;
347	ee2 = e2Central;
348	}
349	else if (cForward > 0) {
350	cc = cForward;
351	ee2 = e2Forward;
352	}
353	c += cc;
354	e2 += ee2;
355	}
356
357	//Systematic errors from here
358
359	Double_t fmd=0;
360	Double_t spd=0;
361	Double_t overlap=0;
362
363	// number of eta bins in fmd, spd and overlap respectively
364	for(Int_t i = first;i<=last;i++){
365	if(normForward->GetBinContent(i)>0&&normCentral->GetBinContent(i)<1)
366	fmd++;
367	if(normForward->GetBinContent(i)>0&&normCentral->GetBinContent(i)>0)
368	overlap++;
369	if(normCentral->GetBinContent(i)>0&&normForward->GetBinContent(i)<1)
370	spd++;
371	}
372
373	Double_t sysErrorSquared = 0;
374	// estimate of systematic uncertainty on the event multiplicity - hardcoded :(. estimates taken from Hans Hjersing Dalsgaard or Casper Nygaard phd theses.
375	Double_t fmdSysError= 0.08;
376	Double_t spdSysError= 0.067;
377
378	Double_t total = 0;
379	total= fmd+ spd + overlap;
380	// Combined systematc event uncertainty, by weighting with the number of eta-bins of fmd-only, spd-only and the overlap.
381	if(total){
382	sysErrorSquared= (fmdTMath::Power(fmdSysError,2)+ spdTMath::Power(spdSysError,2)+
383	0.5overlapTMath::Power(fmdSysError,2)+ 0.5overlapTMath::Power(spdSysError,2))/total;
384	}
385
386	//Strangeness correction, assumed to be 2.0 +- 1% (syst), taken from HHD and CN thesis
387	c=0.98*c;
388
389
390	// correction for missing material description. Correction is based on separate PbPb analysis of difference between nominel and displaced vertices
391	if(fEtaHigh<1.55&&fEtaHigh>1.45)
392	c=0.98*c;
393	if(fEtaHigh<2.05&&fEtaHigh>1.95)
394	c=0.963*c;
395	if(fEtaHigh<2.45&&fEtaHigh>2.35)
396	c=0.956*c;
397	if(fEtaHigh>2.95)
398	c=0.955*c;
399
400	fHist->Fill(c);
401	fHistPlusSys->Fill((1+TMath::Sqrt(sysErrorSquared))*c);
402	fHistMinusSys->Fill((1-TMath::Sqrt(sysErrorSquared))*c);
403	fHistPlus05->Fill(1.05*c);
404	fHistPlus075->Fill(1.075*c);
405	fHistPlus10->Fill(1.1*c);
406	fHistMinus05->Fill(0.95*c);
407	fHistMinus075->Fill(0.925*c);
408	fHistMinus10->Fill(0.9*c);
409
410
411	}
412
413
414
415
416	//_____________________________________________________________________
417	//
418	//
419	// EOF