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

/**
 * @file   AliForwardMultiplicityDistribution.cxx
 * @author Christian Holm Christensen <cholm@master.hehi.nbi.dk>
 * @date   Thu Feb  7 01:03:52 2013
 * 
 * @brief  
 * 
 * @ingroup pwglf_forward_multdist
 * 
 */
#include <TH1D.h>
#include "AliForwardMultiplicityDistribution.h"
#include "AliAODForwardMult.h"
#include "AliAODCentralMult.h"
#include "AliAODEvent.h"
#include "TString.h"

ClassImp(AliForwardMultiplicityDistribution)
#if 0
; // This is for Emacs - do not delete
#endif
//______________________________________________________________________
AliForwardMultiplicityDistribution::AliForwardMultiplicityDistribution() 
  : AliBaseAODTask(), 
    fBins(),       // eta bin list
    fNBins(-1)     // multiplicity axis' runs from 0 to fNbins
{
  //
  // Default Constructor
  //
}
 
//_____________________________________________________________________
AliForwardMultiplicityDistribution::AliForwardMultiplicityDistribution(const char* name) 
  : AliBaseAODTask(name),
    fBins(),       // eta bin list
    fNBins(-1)    // multiplicity axis' runs from 0 to fNbins
{
  //
  // Constructor
  //
}
//_____________________________________________________________________
Bool_t AliForwardMultiplicityDistribution::Book()
{
  TH1D* dndetaSumForward = new TH1D("dndetaSumForward","dndetaSumForward", 
				    200,-4,6);
  TH1D* dndetaSumCentral = new TH1D("dndetaSumCentral","dndetaSumCentral", 
				    200,-4,6);
  fSums->Add(dndetaSumForward);
  fSums->Add(dndetaSumCentral);
  
  
  TH1D* dndetaEventForward = new TH1D();
  TH1D* dndetaEventCentral = new TH1D();
  fSums->Add(dndetaEventForward);
  fSums->Add(dndetaEventCentral);
  
  //Loop over all individual eta bins, and define their hisograms 
  TIter next(&fBins);
  Bin * bin = 0;
  while ((bin = static_cast<Bin*>(next()))) { 
    bin->CreateOutputObjectss(fSums, fNBins);
  }
  return true;
}


//_____________________________________________________________________
Bool_t AliForwardMultiplicityDistribution::Event(AliAODEvent& aod)
{
  //
  // User Exec
  //
  AliAODForwardMult* AODforward = GetForward(aod);
  if (!AODforward) return false;

  AliAODCentralMult* AODcentral = GetCentral(aod);
  if (!AODcentral) return false;

  TH2D& forward = AODforward->GetHistogram();
  TH2D& central = AODcentral->GetHistogram();
  
  TH1D* dndetaSumForward   = (TH1D*)fSums->FindObject("dndetaSumForward");
  TH1D* dndetaSumCentral   = (TH1D*)fSums->FindObject("dndetaSumCentral");
  TH1D* dndetaEventForward = 0;//(TH1D*)fList->FindObject("dndetaEventForward");
  TH1D* dndetaEventCentral = 0;//(TH1D*)fList->FindObject("dndetaEventCentral");
  TH1D* normEventForward   = 0;
  TH1D* normEventCentral   = 0;

  // ACHTUNG ACHTUNG ACHTUNG! Serious memory leak here!
  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
  // 
  // ACHTUNG ACHTUNG ACHTUNG! Serious memory leak here!
  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);
  }
  
  // ACHTUNG ACHTUNG ACHTUNG! Serious memory leak here!
  // We should do 
  // 
  // delete dndetaEventForward;
  // delete dndetaEventCentral;
  // delete normEventForward;
  // delete normEventCentral;
  // 
  // and remove them from the output list - to prevent the memory leak
  
  return true;
  
}

//______________________________________________________________________
#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
33438b4c	1	/**
	2	* @file AliForwardMultiplicityDistribution.cxx
	3	* @author Christian Holm Christensen <cholm@master.hehi.nbi.dk>
	4	* @date Thu Feb 7 01:03:52 2013
	5	*
	6	* @brief
	7	*
	8	* @ingroup pwglf_forward_multdist
	9	*
	10	*/
ad7be237	11	#include <TH1D.h>
	12	#include "AliForwardMultiplicityDistribution.h"
	13	#include "AliAODForwardMult.h"
	14	#include "AliAODCentralMult.h"
	15	#include "AliAODEvent.h"
	16	#include "TString.h"
ad7be237	17
	18	ClassImp(AliForwardMultiplicityDistribution)
	19	#if 0
	20	; // This is for Emacs - do not delete
	21	#endif
	22	//______________________________________________________________________
	23	AliForwardMultiplicityDistribution::AliForwardMultiplicityDistribution()
c8b1a7db	24	: AliBaseAODTask(),
ad7be237	25	fBins(), // eta bin list
c8b1a7db	26	fNBins(-1) // multiplicity axis' runs from 0 to fNbins
ad7be237	27	{
	28	//
	29	// Default Constructor
	30	//
	31	}
	32
	33	//_____________________________________________________________________
	34	AliForwardMultiplicityDistribution::AliForwardMultiplicityDistribution(const char* name)
c8b1a7db	35	: AliBaseAODTask(name),
ad7be237	36	fBins(), // eta bin list
c8b1a7db	37	fNBins(-1) // multiplicity axis' runs from 0 to fNbins
ad7be237	38	{
	39	//
	40	// Constructor
	41	//
ad7be237	42	}
ad7be237	43	//_____________________________________________________________________
c8b1a7db	44	Bool_t AliForwardMultiplicityDistribution::Book()
ad7be237	45	{
c8b1a7db	46	TH1D* dndetaSumForward = new TH1D("dndetaSumForward","dndetaSumForward",
	47	200,-4,6);
	48	TH1D* dndetaSumCentral = new TH1D("dndetaSumCentral","dndetaSumCentral",
	49	200,-4,6);
	50	fSums->Add(dndetaSumForward);
	51	fSums->Add(dndetaSumCentral);
ad7be237	52
ad7be237	53
	54	TH1D* dndetaEventForward = new TH1D();
	55	TH1D* dndetaEventCentral = new TH1D();
c8b1a7db	56	fSums->Add(dndetaEventForward);
c8b1a7db	57	fSums->Add(dndetaEventCentral);
ad7be237	58
	59	//Loop over all individual eta bins, and define their hisograms
	60	TIter next(&fBins);
	61	Bin * bin = 0;
	62	while ((bin = static_cast<Bin*>(next()))) {
c8b1a7db	63	bin->CreateOutputObjectss(fSums, fNBins);
ad7be237	64	}
c8b1a7db	65	return true;
ad7be237	66	}
	67
	68
	69	//_____________________________________________________________________
c8b1a7db	70	Bool_t AliForwardMultiplicityDistribution::Event(AliAODEvent& aod)
ad7be237	71	{
	72	//
	73	// User Exec
	74	//
c8b1a7db	75	AliAODForwardMult* AODforward = GetForward(aod);
c8b1a7db	76	if (!AODforward) return false;
ad7be237	77
c8b1a7db	78	AliAODCentralMult* AODcentral = GetCentral(aod);
	79	if (!AODcentral) return false;
	80
	81	TH2D& forward = AODforward->GetHistogram();
	82	TH2D& central = AODcentral->GetHistogram();
ad7be237	83
c8b1a7db	84	TH1D* dndetaSumForward = (TH1D*)fSums->FindObject("dndetaSumForward");
	85	TH1D* dndetaSumCentral = (TH1D*)fSums->FindObject("dndetaSumCentral");
	86	TH1D* dndetaEventForward = 0;//(TH1D*)fList->FindObject("dndetaEventForward");
	87	TH1D* dndetaEventCentral = 0;//(TH1D*)fList->FindObject("dndetaEventCentral");
73b32206	88	TH1D* normEventForward = 0;
73b32206	89	TH1D* normEventCentral = 0;
ad7be237	90
c8b1a7db	91	// ACHTUNG ACHTUNG ACHTUNG! Serious memory leak here!
ad7be237	92	dndetaEventForward = forward.ProjectionX("dndetaForward",1,forward.GetNbinsY(),"");
	93	dndetaEventCentral = central.ProjectionX("dndetaCentral",1,central.GetNbinsY(),"");
	94	dndetaSumForward->Add(dndetaEventForward);
	95	dndetaSumCentral->Add(dndetaEventCentral);
	96
c8b1a7db	97	// underflow eta bin of forward/central carry information on whether
	98	// there is acceptance. 1= acceptance, 0= no acceptance
	99	//
	100	// ACHTUNG ACHTUNG ACHTUNG! Serious memory leak here!
ad7be237	101	normEventForward = forward.ProjectionX("dndetaEventForwardNorm",0,0,"");
	102	normEventCentral = central.ProjectionX("dndetaEventCentralNorm",0,0,"");
	103
	104	Double_t VtxZ = AODforward->GetIpZ();
	105
	106
	107	// loop over all eta bins, and fill multiplicity histos
	108	TIter next(&fBins);
	109	Bin * bin = 0;
	110	while ((bin = static_cast<Bin*>(next()))) {
	111	bin->Process(dndetaEventForward, dndetaEventCentral,
	112	normEventForward, normEventCentral, VtxZ);
	113	}
	114
c8b1a7db	115	// ACHTUNG ACHTUNG ACHTUNG! Serious memory leak here!
	116	// We should do
	117	//
	118	// delete dndetaEventForward;
	119	// delete dndetaEventCentral;
	120	// delete normEventForward;
	121	// delete normEventCentral;
	122	//
	123	// and remove them from the output list - to prevent the memory leak
ad7be237	124
c8b1a7db	125	return true;
ad7be237	126
	127	}
	128
	129	//______________________________________________________________________
	130	#if 1
	131	const Char_t*
	132	#else
	133	TString
	134	#endif
	135	AliForwardMultiplicityDistribution::FormBinName(Double_t etaLow, Double_t etaHigh)
	136	{
	137	//
	138	// Form name of eta bin
	139	//
	140	TString sLow(TString::Format("%+03d",int(10*etaLow)));
	141	TString sHigh(TString::Format("%+03d",int(10*etaHigh)));
	142	sLow.ReplaceAll("+", "plus");
	143	sLow.ReplaceAll("-", "minus");
	144	sHigh.ReplaceAll("+", "plus");
	145	sHigh.ReplaceAll("-", "minus");
	146	#if 0
	147	TString* combined= new TString();
	148	combined->Append(TString::Format("%s_%s", sLow.Data(), sHigh.Data()));
	149	return combined->Data();
	150	#else
	151	static TString tmp;
	152	tmp = TString::Format("%s_%s", sLow.Data(), sHigh.Data());
	153	return tmp.Data();
	154	#endif
	155	}
	156
	157	//=====================================================================
	158	AliForwardMultiplicityDistribution::Bin::Bin()
	159	: TNamed(),
	160	fEtaLow(-1111), // low eta limit
	161	fEtaHigh(-1111), // high eta limit
	162	fHist(0), // multiplicity distribution hist
	163	fHistPlus05(0), // multiplicity distribution hist scaled up with 5%
	164	fHistPlus075(0), // multiplicity distribution hist scaled up with 7.5%
	165	fHistPlus10(0), // multiplicity distribution hist scaled up with 10%
	166	fHistMinus05(0), // multiplicity distribution hist scaled down with 5%
	167	fHistMinus075(0), // multiplicity distribution hist scaled down with 7.5%
	168	fHistMinus10(0), // multiplicity distribution hist scaled down with 10%
	169	fHistPlusSys(0), // multiplicity distribution hist scaled up with the event uncertainty
	170	fHistMinusSys(0), // multiplicity distribution hist scaled down with the event uncertainty
	171	fAcceptance(0), // histogram showing the 'holes' in acceptance.
	172	fVtxZvsNdataBins(0) // VtxZ vs. number of data acceptance bins (normalised to the eta range)
	173	{
	174	//
	175	// Default Constructor
	176	//
	177	}
	178	//_____________________________________________________________________
	179	AliForwardMultiplicityDistribution::Bin::Bin(Double_t etaLow, Double_t etaHigh)
	180	: TNamed("", ""),
	181	fEtaLow(etaLow), // low eta limit
	182	fEtaHigh(etaHigh), // high eta limit
	183	fHist(0), // multiplicity distribution hist
	184	fHistPlus05(0), // multiplicity distribution hist scaled up with 5%
	185	fHistPlus075(0), // multiplicity distribution hist scaled up with 7.5%
	186	fHistPlus10(0), // multiplicity distribution hist scaled up with 10%
	187	fHistMinus05(0), // multiplicity distribution hist scaled down with 5%
	188	fHistMinus075(0), // multiplicity distribution hist scaled down with 7.5%
	189	fHistMinus10(0), // multiplicity distribution hist scaled down with 10%
190	fHistPlusSys(0), // multiplicity distribution hist scaled up with the event uncertainty
191	fHistMinusSys(0), // multiplicity distribution hist scaled down with the event uncertainty
192	fAcceptance(0), // histogram showing the 'holes' in acceptance.
193	fVtxZvsNdataBins(0) // VtxZ vs. number of data acceptance bins (normalised to the eta range)
194	{
195	//
196	// Constructor
197	//
198	const char* name = AliForwardMultiplicityDistribution::FormBinName(fEtaLow,fEtaHigh);
199
200	SetName(name);
201	SetTitle(TString::Format("%+4.1f < #eta < %+4.1f", fEtaLow, fEtaHigh));
202
203	}
204	//_____________________________________________________________________
5934a3e3	205	void AliForwardMultiplicityDistribution::Bin::CreateOutputObjectss(TList* cont, Int_t max)
ad7be237	206	{
	207	//
	208	// Define eta bin output histograms
	209	//
	210	TList* out = new TList;
	211	out->SetName(GetName());
	212	cont->Add(out);
	213
	214	fHist = new TH1D("mult", GetTitle(), max, -0.5, max-.5);
	215	fHistPlus05 = new TH1D("multPlus05", GetTitle(), max, -0.5, max-.5);
	216	fHistPlus075 = new TH1D("multPlus075", GetTitle(), max, -0.5, max-.5);
	217	fHistPlus10 = new TH1D("multPlus10", GetTitle(), max, -0.5, max-.5);
	218	fHistMinus05 = new TH1D("multMinus05", GetTitle(), max, -0.5, max-.5);
	219	fHistMinus075 = new TH1D("multMinus075", GetTitle(), max, -0.5, max-.5);
	220	fHistMinus10 = new TH1D("multMinus10", GetTitle(), max, -0.5, max-.5);
	221	fHistPlusSys = new TH1D("multPlusSys", GetTitle(), max, -0.5, max-.5);
	222	fHistMinusSys = new TH1D("multMinusSys", GetTitle(), max, -0.5, max-.5);
	223	fVtxZvsNdataBins = new TH2D("VtxZvsNdataBins", "VtxZ vs dataAcceptance/etaRange;z-vtz;dataAcceptance/etaRange", 20, -10,10, 130,0,1.3);
	224	fAcceptance = new TH2D("Acceptance","Acceptance;#eta;z-vtx",200,-4, 6 , 20,-10,10);
	225
	226	out->Add(fHist);
	227	out->Add(fHistPlus05);
	228	out->Add(fHistPlus075);
	229	out->Add(fHistPlus10);
	230	out->Add(fHistMinus05);
	231	out->Add(fHistMinus075);
	232	out->Add(fHistMinus10);
	233	out->Add(fHistPlusSys);
	234	out->Add(fHistMinusSys);
	235	out->Add(fVtxZvsNdataBins);
	236	out->Add(fAcceptance);
	237
	238	}
	239
	240	//_____________________________________________________________________
	241	void AliForwardMultiplicityDistribution::Bin::Process(TH1D* dndetaForward, TH1D* dndetaCentral,
	242	TH1D* normForward, TH1D* normCentral, Double_t VtxZ)
	243	{
	244	//
	245	// Process single eta bin
	246	//
	247
	248	// Diagnostics on event acceptance
	249	Int_t first = dndetaForward->GetXaxis()->FindBin(fEtaLow);
	250	Int_t last = dndetaForward->GetXaxis()->FindBin(fEtaHigh-.0001);
	251	Double_t acceptanceBins=0;
	252	for(Int_t n= first;n<=last;n++){
	253	if(normForward->GetBinContent(n)>0\|\|normCentral->GetBinContent(n)>0){
	254	acceptanceBins++;
	255	}
	256	fAcceptance->SetBinContent(n, fAcceptance->GetYaxis()->FindBin(VtxZ), 1);
	257	if(normForward->GetBinContent(n)>0\|\|normCentral->GetBinContent(n)>0)
	258	fAcceptance->SetBinContent(n, fAcceptance->GetYaxis()->FindBin(VtxZ),10);
	259	}
	260
	261	//std::cout << "% of bins covered in eta-range: " << acceptanceBins/(last-first+1) << std::endl;
	262	fVtxZvsNdataBins->Fill(VtxZ, (Double_t)acceptanceBins/(last-first+1));
	263
	264
	265	Double_t c = 0;
	266	Double_t e2 = 0;
	267	for (Int_t i = first; i <= last; i++) {
	268	Double_t cForward = 0;
	269	Double_t cCentral = 0;
270	Double_t e2Forward = 0;
271	Double_t e2Central = 0;
272	if (normForward->GetBinContent(i) > 0) {
273	cForward = dndetaForward->GetBinContent(i);
274	e2Forward = dndetaForward->GetBinError(i) * dndetaForward->GetBinError(i);
275	}
276	if (normCentral->GetBinContent(i) > 0) {
277	cCentral = dndetaCentral->GetBinContent(i);
278	e2Central = dndetaCentral->GetBinError(i) * dndetaCentral->GetBinError(i);
279	}
280	Double_t cc = 0;
281	Double_t ee2 = 0;
282	// if overlap between central/forward, use average of the two
283	if (cCentral > 0 && cForward > 0) {
284	cc = 0.5 * (cForward + cCentral);
285	ee2 = 0.5 * (e2Forward + e2Central);
286	}
287	else if (cCentral > 0) {
288	cc = cCentral;
289	ee2 = e2Central;
290	}
291	else if (cForward > 0) {
292	cc = cForward;
293	ee2 = e2Forward;
294	}
295	c += cc;
296	e2 += ee2;
297	}
298
299	//Systematic errors from here
300
301	Double_t fmd=0;
302	Double_t spd=0;
303	Double_t overlap=0;
304
305	// number of eta bins in fmd, spd and overlap respectively
306	for(Int_t i = first;i<=last;i++){
307	if(normForward->GetBinContent(i)>0&&normCentral->GetBinContent(i)<1)
308	fmd++;
309	if(normForward->GetBinContent(i)>0&&normCentral->GetBinContent(i)>0)
310	overlap++;
311	if(normCentral->GetBinContent(i)>0&&normForward->GetBinContent(i)<1)
312	spd++;
313	}
314
315	Double_t sysErrorSquared = 0;
316	// estimate of systematic uncertainty on the event multiplicity - hardcoded :(. estimates taken from Hans Hjersing Dalsgaard or Casper Nygaard phd theses.
317	Double_t fmdSysError= 0.08;
318	Double_t spdSysError= 0.067;
319
320	Double_t total = 0;
321	total= fmd+ spd + overlap;
322	// Combined systematc event uncertainty, by weighting with the number of eta-bins of fmd-only, spd-only and the overlap.
323	if(total){
324	sysErrorSquared= (fmdTMath::Power(fmdSysError,2)+ spdTMath::Power(spdSysError,2)+
325	0.5overlapTMath::Power(fmdSysError,2)+ 0.5overlapTMath::Power(spdSysError,2))/total;
326	}
327
328	//Strangeness correction, assumed to be 2.0 +- 1% (syst), taken from HHD and CN thesis
329	c=0.98*c;
330
331
332	// correction for missing material description. Correction is based on separate PbPb analysis of difference between nominel and displaced vertices
333	if(fEtaHigh<1.55&&fEtaHigh>1.45)
334	c=0.98*c;
335	if(fEtaHigh<2.05&&fEtaHigh>1.95)
336	c=0.963*c;
337	if(fEtaHigh<2.45&&fEtaHigh>2.35)
338	c=0.956*c;
339	if(fEtaHigh>2.95)
340	c=0.955*c;
341
342	fHist->Fill(c);
343	fHistPlusSys->Fill((1+TMath::Sqrt(sysErrorSquared))*c);
344	fHistMinusSys->Fill((1-TMath::Sqrt(sysErrorSquared))*c);
345	fHistPlus05->Fill(1.05*c);
346	fHistPlus075->Fill(1.075*c);
347	fHistPlus10->Fill(1.1*c);
348	fHistMinus05->Fill(0.95*c);
349	fHistMinus075->Fill(0.925*c);
350	fHistMinus10->Fill(0.9*c);
351
352
353	}
354
355
356
357
358	//_____________________________________________________________________
359	//
360	//
361	// EOF