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

/**
 * @file   AliForwardCreateResponseMatrices.cxx
 * @author Christian Holm Christensen <cholm@master.hehi.nbi.dk>
 * @date   Thu Feb  7 01:01:24 2013
 * 
 * @brief  
 * 
 * 
 * @ingroup pwglf_forward_multdist
 */

#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"


ClassImp(AliForwardCreateResponseMatrices)
#if 0
; // This is for Emacs - do not delete
#endif
//_____________________________________________________________________
AliForwardCreateResponseMatrices::AliForwardCreateResponseMatrices() 
  : AliBaseAODTask(),
    fBins(), 
    fIsSelected(false)
{
  //
  // Default Constructor
  //
}
 
//_____________________________________________________________________
AliForwardCreateResponseMatrices::AliForwardCreateResponseMatrices(const char* name) 
  : AliBaseAODTask(name,"AliForwardCreateResponseMatrices"),
    fBins(), 
    fIsSelected(false)
{
  //
  // Constructor
  //
}

//_____________________________________________________________________
Bool_t 
AliForwardCreateResponseMatrices::CheckEvent(const AliAODForwardMult& fwd)
{
  fIsSelected = AliBaseAODTask::CheckEvent(fwd);
  // We always return true, so that we can process MC truth;
  return true;
}

//_____________________________________________________________________
Bool_t AliForwardCreateResponseMatrices::Book()
{
  //
  // Create Output Objects
  //
  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);
 
  fSums->Add(dndetaSumForward);
  fSums->Add(dndetaSumCentral);
  fSums->Add(dndetaSumMC);

  TH1D* dndetaEventForward = new TH1D();
  TH1D* dndetaEventCentral = new TH1D();
  TH1D* dndetaEventMC = new TH1D();
  
  fSums->Add(dndetaEventForward);
  fSums->Add(dndetaEventCentral);
  fSums->Add(dndetaEventMC);

  //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, 400);
  }
 
  return true;
  
}


//_____________________________________________________________________
Bool_t AliForwardCreateResponseMatrices::Event(AliAODEvent& aod)
{
  //
  // User Exec
  //
  // Here, we used to get ForwardMC!
  AliAODForwardMult* AODforward = GetForward(aod);
  if (!AODforward) return false;

  // Here, we used to get CentralClusters!
  AliAODCentralMult* AODcentral = GetCentral(aod);
  if (!AODcentral) return false;

  TH2D& forward = AODforward->GetHistogram();
  TH2D& central = AODcentral->GetHistogram();
  
  //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   = GetPrimary(aod);
  
  TH1D* dndetaSumForward    = (TH1D*)fSums->FindObject("dndetaSumForward");
  TH1D* dndetaSumCentral    = (TH1D*)fSums->FindObject("dndetaSumCentral");
  TH1D* dndetaSumMC         = (TH1D*)fSums->FindObject("dndetaSumMC");
  TH1D* dndetaEventForward  = forward.ProjectionX("dndetaForward",1,
						  forward.GetNbinsY(),"");
  TH1D* dndetaEventCentral  = central.ProjectionX("dndetaCentral",1,
						  central.GetNbinsY(),"");
  TH1D* 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;
  normEventForward       = forward.ProjectionX("dndetaEventForwardNorm",0,0,"");
  normEventCentral       = central.ProjectionX("dndetaEventCentralNorm",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, fIsSelected,
		 isMCNSD, isESDNSD,  aod);
  }
    
  return true;  
}


//=====================================================================
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::CreateOutputObjectss(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, 
					       const 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;
  }
  AliAODMCHeader* header = 
    dynamic_cast<AliAODMCHeader*>(aodevent.
				  FindListObject(AliAODMCHeader::StdBranchName()));
  if (!header) {
    AliWarning("No header file found.");
    return;
  }
  
  
  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
Commit	Line	Data
33438b4c	1	/**
	2	* @file AliForwardCreateResponseMatrices.cxx
	3	* @author Christian Holm Christensen <cholm@master.hehi.nbi.dk>
	4	* @date Thu Feb 7 01:01:24 2013
	5	*
	6	* @brief
	7	*
	8	*
	9	* @ingroup pwglf_forward_multdist
	10	*/
ad7be237	11
	12	#include <TH1D.h>
	13	#include "AliForwardCreateResponseMatrices.h"
	14	#include "AliForwardMultiplicityDistribution.h"
	15	#include "AliAODForwardMult.h"
	16	#include "AliAODCentralMult.h"
	17	#include "AliAODEvent.h"
	18	#include "AliFMDMCEventInspector.h"
	19	#include "AliAODMCHeader.h"
	20	#include "AliAODMCParticle.h"
	21
ad7be237	22
	23	ClassImp(AliForwardCreateResponseMatrices)
	24	#if 0
	25	; // This is for Emacs - do not delete
	26	#endif
	27	//_____________________________________________________________________
	28	AliForwardCreateResponseMatrices::AliForwardCreateResponseMatrices()
c8b1a7db	29	: AliBaseAODTask(),
	30	fBins(),
	31	fIsSelected(false)
ad7be237	32	{
	33	//
	34	// Default Constructor
	35	//
	36	}
	37
	38	//_____________________________________________________________________
	39	AliForwardCreateResponseMatrices::AliForwardCreateResponseMatrices(const char* name)
bb64e168	40	: AliBaseAODTask(name,"AliForwardCreateResponseMatrices"),
c8b1a7db	41	fBins(),
c8b1a7db	42	fIsSelected(false)
ad7be237	43	{
	44	//
	45	// Constructor
	46	//
ad7be237	47	}
	48
	49	//_____________________________________________________________________
c8b1a7db	50	Bool_t
	51	AliForwardCreateResponseMatrices::CheckEvent(const AliAODForwardMult& fwd)
	52	{
	53	fIsSelected = AliBaseAODTask::CheckEvent(fwd);
	54	// We always return true, so that we can process MC truth;
	55	return true;
	56	}
	57
	58	//_____________________________________________________________________
	59	Bool_t AliForwardCreateResponseMatrices::Book()
ad7be237	60	{
	61	//
	62	// Create Output Objects
	63	//
c8b1a7db	64	TH1D* dndetaSumForward = new TH1D("dndetaSumForward","dndetaSumForward",
	65	200,-4,6);
	66	TH1D* dndetaSumCentral = new TH1D("dndetaSumCentral","dndetaSumCentral",
	67	200,-4,6);
ad7be237	68	TH1D* dndetaSumMC = new TH1D("dndetaSumMC","dndetaSumMC", 200,-4,6);
ad7be237	69
c8b1a7db	70	fSums->Add(dndetaSumForward);
	71	fSums->Add(dndetaSumCentral);
	72	fSums->Add(dndetaSumMC);
ad7be237	73
	74	TH1D* dndetaEventForward = new TH1D();
	75	TH1D* dndetaEventCentral = new TH1D();
	76	TH1D* dndetaEventMC = new TH1D();
	77
c8b1a7db	78	fSums->Add(dndetaEventForward);
	79	fSums->Add(dndetaEventCentral);
	80	fSums->Add(dndetaEventMC);
ad7be237	81
	82	//Loop over all individual eta bins, and define their hisograms
	83	TIter next(&fBins);
	84	Bin * bin = 0;
	85	while ((bin = static_cast<Bin*>(next()))) {
c8b1a7db	86	bin->CreateOutputObjectss(fSums, 400);
ad7be237	87	}
ad7be237	88
c8b1a7db	89	return true;
ad7be237	90
	91	}
	92
	93
	94	//_____________________________________________________________________
c8b1a7db	95	Bool_t AliForwardCreateResponseMatrices::Event(AliAODEvent& aod)
ad7be237	96	{
	97	//
	98	// User Exec
	99	//
c8b1a7db	100	// Here, we used to get ForwardMC!
	101	AliAODForwardMult* AODforward = GetForward(aod);
	102	if (!AODforward) return false;
	103
	104	// Here, we used to get CentralClusters!
	105	AliAODCentralMult* AODcentral = GetCentral(aod);
	106	if (!AODcentral) return false;
	107
	108	TH2D& forward = AODforward->GetHistogram();
	109	TH2D& central = AODcentral->GetHistogram();
ad7be237	110
ad7be237	111	//check if event is NSD according to either MC truth or analysis (ESD)
c8b1a7db	112	Bool_t isMCNSD = kFALSE;
	113	Bool_t isESDNSD = kFALSE;
	114	if(AODforward->IsTriggerBits(AliAODForwardMult::kMCNSD)) isMCNSD=kTRUE;
	115	if(AODforward->IsTriggerBits(AliAODForwardMult::kNSD)) isESDNSD=kTRUE;
ad7be237	116
ad7be237	117	//primary dndeta dist
c8b1a7db	118	TH2D* primHist = GetPrimary(aod);
ad7be237	119
c8b1a7db	120	TH1D* dndetaSumForward = (TH1D*)fSums->FindObject("dndetaSumForward");
	121	TH1D* dndetaSumCentral = (TH1D*)fSums->FindObject("dndetaSumCentral");
	122	TH1D* dndetaSumMC = (TH1D*)fSums->FindObject("dndetaSumMC");
	123	TH1D* dndetaEventForward = forward.ProjectionX("dndetaForward",1,
	124	forward.GetNbinsY(),"");
	125	TH1D* dndetaEventCentral = central.ProjectionX("dndetaCentral",1,
	126	central.GetNbinsY(),"");
	127	TH1D* dndetaEventMC = primHist->ProjectionX("dndetaMC",1,
	128	primHist->GetNbinsY(),"");
ad7be237	129
ba192922	130	// underflow eta bin of forward/central carry information on whether
ba192922	131	// there is acceptance. 1= acceptance, 0= no acceptance
c8b1a7db	132	TH1D* normEventForward = 0;
	133	TH1D* normEventCentral = 0;
	134	normEventForward = forward.ProjectionX("dndetaEventForwardNorm",0,0,"");
	135	normEventCentral = central.ProjectionX("dndetaEventCentralNorm",0,0,"");
ad7be237	136
	137	dndetaSumForward->Add(dndetaEventForward);
	138	dndetaSumCentral->Add(dndetaEventCentral);
	139	dndetaSumMC->Add(dndetaEventMC);
	140
	141
	142	Double_t VtxZ = AODforward->GetIpZ();
	143
ba192922	144	// loop over all eta bins, and create response matrices,
ba192922	145	// trigger-vertex bias histograms etc
ad7be237	146	TIter next(&fBins);
	147	Bin * bin = 0;
	148	while ((bin = static_cast<Bin*>(next()))) {
	149	bin->Process(dndetaEventForward, dndetaEventCentral,
ba192922	150	normEventForward, normEventCentral,
c8b1a7db	151	dndetaEventMC, VtxZ, fIsSelected,
ba192922	152	isMCNSD, isESDNSD, aod);
ad7be237	153	}
ad7be237	154
c8b1a7db	155	return true;
ad7be237	156	}
	157
	158
	159	//=====================================================================
	160	AliForwardCreateResponseMatrices::Bin::Bin(Double_t etaLow, Double_t etaHigh)
	161	: TNamed("", ""),
	162	fEtaLow(etaLow), // low eta limit
	163	fEtaHigh(etaHigh), // high eta limit
	164	fHist(), // multiplicity histogram
	165	fHistMC(), // multiplicity histogram MC truth primaries
	166	fAcceptance(), // histogram showing the 'holes' in acceptance. BinContent of 1 shows a hole, and BinContent of 10 shows data coverage
	167	fVtxZvsNdataBins(), // VtxZ vs. number of data acceptance bins (normalised to the eta range)
	168	fResponseMatrix(), // Response matrix (MC truth vs. analysed multiplicity)
	169	fResponseMatrixPlus05(), // Response matrix with analysed multiplicity scaled up by 5%
	170	fResponseMatrixPlus075(), // Response matrix with analysed multiplicity scaled up by 7.5%
	171	fResponseMatrixPlus10(), // Response matrix with analysed multiplicity scaled up by 10%
	172	fResponseMatrixMinus05(), // Response matrix with analysed multiplicity scaled down by 5%
	173	fResponseMatrixMinus075(),// Response matrix with analysed multiplicity scaled down by 7.55%
	174	fResponseMatrixMinus10(), // Response matrix with analysed multiplicity scaled down by 10%
	175	fResponseMatrixMinusSys(),// Response matrix with analysed multiplicity scaled up by event mult uncertainty
	176	fResponseMatrixPlusSys(), // Response matrix with analysed multiplicity scaled down by event mult uncertainty
	177	fESDNSD(), // number of events found as NSD by the analysis vs. multiplicity
	178	fMCNSD(), // number of events found as NSD by the MC truth vs. multiplicity
	179	fMCESDNSD(), // number of events found as NSD by both analysis and MC truth vs. multiplicity
	180	fTriggerBias() // histogram for trigger vertex bias correction
	181	{
	182	//
	183	// Constructor
	184	//
	185	const char* name = AliForwardMultiplicityDistribution::FormBinName(etaLow,etaHigh);
	186
	187	SetName(name);
	188	SetTitle(Form("%+4.1f < #eta < %+4.1f", fEtaLow, fEtaHigh));
	189
	190	}
	191	//_____________________________________________________________________
	192	AliForwardCreateResponseMatrices::Bin::Bin()
	193	: TNamed(),
	194	fEtaLow(), // low eta limit
	195	fEtaHigh(), // high eta limit
	196	fHist(), // multiplicity histogram
	197	fHistMC(), // multiplicity histogram MC truth primaries
	198	fAcceptance(), // histogram showing the 'holes' in acceptance. BinContent of 1 shows a hole, and BinContent of 10 shows data coverage
	199	fVtxZvsNdataBins(), // VtxZ vs. number of data acceptance bins (normalised to the eta range)
	200	fResponseMatrix(), // Response matrix (MC truth vs. analysed multiplicity)
	201	fResponseMatrixPlus05(), // Response matrix with analysed multiplicity scaled up by 5%
	202	fResponseMatrixPlus075(), // Response matrix with analysed multiplicity scaled up by 7.5%
	203	fResponseMatrixPlus10(), // Response matrix with analysed multiplicity scaled up by 10%
	204	fResponseMatrixMinus05(), // Response matrix with analysed multiplicity scaled down by 5%
	205	fResponseMatrixMinus075(),// Response matrix with analysed multiplicity scaled down by 7.55%
	206	fResponseMatrixMinus10(), // Response matrix with analysed multiplicity scaled down by 10%
	207	fResponseMatrixMinusSys(),// Response matrix with analysed multiplicity scaled up by event mult uncertainty
	208	fResponseMatrixPlusSys(), // Response matrix with analysed multiplicity scaled down by event mult uncertainty
	209	fESDNSD(), // number of events found as NSD by the analysis vs. multiplicity
	210	fMCNSD(), // number of events found as NSD by the MC truth vs. multiplicity
	211	fMCESDNSD(), // number of events found as NSD by both analysis and MC truth vs. multiplicity
	212	fTriggerBias() // histogram for trigger vertex bias correction
	213	{
	214	//
	215	// Default Constructor
	216	//
	217	}
	218	//_____________________________________________________________________
5934a3e3	219	void AliForwardCreateResponseMatrices::Bin::CreateOutputObjectss(TList* cont, Int_t max)
ad7be237	220	{
	221	//
	222	// Define eta bin output histos
	223	//
	224	TList* out = new TList;
	225	out->SetName(GetName());
	226	cont->Add(out);
	227
	228	fHist = new TH1D("mult", GetTitle(), max, -0.5, max-.5);
	229	fHistMC = new TH1D("multTruth", GetTitle(), max, -0.5, max-.5);
	230	fVtxZvsNdataBins = new TH2D("VtxZvsNdataBins", "VtxZ vs dataAcceptance/etaRange;z-vtz;dataAcceptance/etaRange", 20, -10,10, 130,0,1.3);
	231	fAcceptance = new TH2D("Acceptance","Acceptance;#eta;z-vtx",200,-4, 6 , 20,-10,10);
	232	fResponseMatrix = new TH2D("responseMatrix","Response Matrix;MC_{truth};MC_{measured}", max, -0.5, max-.5, max, -0.5, max-.5);
	233	fResponseMatrixPlus05 = new TH2D("responseMatrixPlus05","Response Matrix;MC_{truth};MC_{measured}", max, -0.5, max-.5, max, -0.5, max-.5);
	234	fResponseMatrixPlus075 = new TH2D("responseMatrixPlus075","Response Matrix;MC_{truth};MC_{measured}", max, -0.5, max-.5, max, -0.5, max-.5);
	235	fResponseMatrixPlus10 = new TH2D("responseMatrixPlus10","Response Matrix;MC_{truth};MC_{measured}", max, -0.5, max-.5, max, -0.5, max-.5);
	236	fResponseMatrixMinus05 = new TH2D("responseMatrixMinus05","Response Matrix;MC_{truth};MC_{measured}", max, -0.5, max-.5, max, -0.5, max-.5);
	237	fResponseMatrixMinus075 = new TH2D("responseMatrixMinus075","Response Matrix;MC_{truth};MC_{measured}", max, -0.5, max-.5, max, -0.5, max-.5);
	238	fResponseMatrixMinus10 = new TH2D("responseMatrixMinus10","Response Matrix;MC_{truth};MC_{measured}", max, -0.5, max-.5, max, -0.5, max-.5);
	239	fResponseMatrixMinusSys = new TH2D("responseMatrixMinusSys","Response Matrix;MC_{truth};MC_{measured}", max, -0.5, max-.5, max, -0.5, max-.5);
	240	fResponseMatrixPlusSys = new TH2D("responseMatrixPlusSys","Response Matrix;MC_{truth};MC_{measured}", max, -0.5, max-.5, max, -0.5, max-.5);
	241	fTriggerBias = new TH1D("triggerBias","triggerBias;MC_{truth};Correction}", max, -0.5, max-.5);
	242	fMCNSD = new TH1D("fMCNSD","fMCNSD", 300,-0.5,299.5);
	243	fESDNSD = new TH1D("fESDNSD","fESDNSD", 300,-0.5,299.5);
	244	fMCESDNSD = new TH1D("fMCESDNSD","fMCESDNSD", 300,-0.5,299.5);
	245
	246	out->Add(fMCNSD);
	247	out->Add(fESDNSD);
	248	out->Add(fMCESDNSD);
	249	out->Add(fHist);
	250	out->Add(fHistMC);
	251	out->Add(fVtxZvsNdataBins);
	252	out->Add(fAcceptance);
	253	out->Add(fResponseMatrix);
	254	out->Add(fResponseMatrixPlus05);
	255	out->Add(fResponseMatrixPlus075);
	256	out->Add(fResponseMatrixPlus10);
	257	out->Add(fResponseMatrixMinus05);
	258	out->Add(fResponseMatrixMinus075);
	259	out->Add(fResponseMatrixMinus10);
	260	out->Add(fResponseMatrixPlusSys);
	261	out->Add(fResponseMatrixMinusSys);
	262	out->Add(fTriggerBias);
	263
	264	}
	265
	266
	267	//_____________________________________________________________________
c8b1a7db	268	void
	269	AliForwardCreateResponseMatrices::Bin::Process(TH1D* dndetaForward,
	270	TH1D* dndetaCentral,
	271	TH1D* normForward,
	272	TH1D* normCentral,
	273	TH1D* mc,
	274	Double_t VtxZ,
	275	Bool_t selectedTrigger,
	276	Bool_t isMCNSD,
	277	Bool_t isESDNSD,
	278	const AliAODEvent& aodevent)
ad7be237	279	{
	280	//
	281	// Process a single eta bin
	282	//
	283
	284	// Diagnostics on event acceptance
	285	Int_t first = dndetaForward->GetXaxis()->FindBin(fEtaLow);
	286	Int_t last = dndetaForward->GetXaxis()->FindBin(fEtaHigh-.0001);
	287	Double_t acceptanceBins=0;
	288	for(Int_t n= first;n<=last;n++){
	289	if(normForward->GetBinContent(n)>0\|\|normCentral->GetBinContent(n)>0){
	290	acceptanceBins++;
	291	}
	292	fAcceptance->SetBinContent(n, fAcceptance->GetYaxis()->FindBin(VtxZ), 1);
	293	if(normForward->GetBinContent(n)>0\|\|normCentral->GetBinContent(n)>0)
	294	fAcceptance->SetBinContent(n, fAcceptance->GetYaxis()->FindBin(VtxZ),10);
	295	}
	296	fVtxZvsNdataBins->Fill(VtxZ, (Double_t)acceptanceBins/(last-first+1));
	297
	298
	299
	300	Double_t c = 0;
	301	Double_t e2 = 0;
	302	Double_t cPrimary = 0;
	303	Double_t e2Primary= 0;
	304
	305	for (Int_t i = first; i <= last; i++){
	306	Double_t cForward = 0;
	307	Double_t cCentral = 0;
	308	Double_t e2Forward = 0;
	309	Double_t e2Central = 0;
	310	Double_t cMC = 0;
	311	Double_t e2MC = 0;
	312	cMC= mc->GetBinContent(i);
	313	e2MC= mc->GetBinError(i) * mc->GetBinError(i);
	314	cPrimary+=cMC;
	315	e2Primary+=e2MC;
	316	if (normForward->GetBinContent(i) > 0) {
	317	cForward = dndetaForward->GetBinContent(i);
	318	e2Forward = dndetaForward->GetBinError(i) * dndetaForward->GetBinError(i);
	319	}
	320	if (normCentral->GetBinContent(i) > 0) {
	321	cCentral = dndetaCentral->GetBinContent(i);
	322	e2Central = dndetaCentral->GetBinError(i) * dndetaCentral->GetBinError(i);
	323	}
	324	Double_t cc = 0;
	325	Double_t ee2 = 0;
	326	if (cCentral > 0 && cForward > 0) {
	327	cc = 0.5 * (cForward + cCentral);
	328	ee2 = 0.5 * (e2Forward + e2Central);
	329	}
	330	else if (cCentral > 0) {
	331	cc = cCentral;
	332	ee2 = e2Central;
	333	}
	334	else if (cForward > 0) {
	335	cc = cForward;
	336	ee2 = e2Forward;
	337	}
	338	c += cc;
	339	e2 += ee2;
	340	}
	341
73b32206	342	// retreive MC particles from event
c8b1a7db	343	TClonesArray* mcArray = (TClonesArray*)aodevent.FindListObject(AliAODMCParticle::StdBranchName());
ad7be237	344	if(!mcArray){
73b32206	345	AliWarning("No MC array found in AOD. Try making it again.");
73b32206	346	return;
ad7be237	347	}
73b32206	348	AliAODMCHeader* header =
c8b1a7db	349	dynamic_cast<AliAODMCHeader*>(aodevent.
73b32206	350	FindListObject(AliAODMCHeader::StdBranchName()));
ad7be237	351	if (!header) {
ad7be237	352	AliWarning("No header file found.");
73b32206	353	return;
ad7be237	354	}
	355
	356
	357	Int_t ntracks = mcArray->GetEntriesFast();
	358	// Track loop: find MC truth multiplicity in selected eta bin
	359	Double_t mult = 0;
	360	for (Int_t it = 0; it < ntracks; it++) {
	361	AliAODMCParticle* particle = (AliAODMCParticle*)mcArray->At(it);
	362	if (!particle) {
	363	AliError(Form("Could not receive track %d", it));
	364	continue;
	365	}
	366	if (!particle->IsPhysicalPrimary()) continue;
	367	if (particle->Charge() == 0) continue;
	368	if(particle->Eta()>fEtaLow&&particle->Eta()<fEtaHigh-0.0001)
	369	mult++;
	370	}
	371	//fill fMCNSD with multiplicity of MC truth NSD events with MC truth \|vtxz\|<4
	372	if(header->GetVtxZ()>-4&&header->GetVtxZ()<4){
	373	if(isMCNSD){
	374	fMCNSD->Fill(mult);
	375	}
	376	}
	377	//fill fESDNSD with multiplicity from events with a reconstructed NSD trigger and reconstructed \|vtxz\|<4
	378	if(VtxZ>-4&&VtxZ<4){
	379	if(isESDNSD){
	380	fESDNSD->Fill(mult);
	381	}
	382	}
	383	// fullfilling both requirements of fMCNSD and fESDNSD
	384	if(header->GetVtxZ()>-4&&header->GetVtxZ()<4&&VtxZ>-4&&VtxZ<4&&isMCNSD&&isESDNSD){
	385	fMCESDNSD->Fill(mult);
	386	}
	387
	388
	389
	390	//Systematic errors from here
	391
	392	Double_t fmd=0;
	393	Double_t spd=0;
	394	Double_t overlap=0;
	395
	396	// number of eta bins in fmd, spd and overlap respectively
	397	for(Int_t i = first;i<=last;i++){
	398	if(normForward->GetBinContent(i)>0&&normCentral->GetBinContent(i)<1)
	399	fmd++;
	400	if(normForward->GetBinContent(i)>0&&normCentral->GetBinContent(i)>0)
	401	overlap++;
	402	if(normCentral->GetBinContent(i)>0&&normForward->GetBinContent(i)<1)
	403	spd++;
	404	}
	405
	406	Double_t sysErrorSquared = 0;
	407
	408	// estimate of systematic uncertainty on the event multiplicity - hardcoded :(. estimates taken from Hans Hjersing Dalsgaard or Casper Nygaard phd theses.
	409	Double_t fmdSysError= 0.08;
	410	Double_t spdSysError= 0.04;
	411	Double_t total = 0;
	412	total= fmd + spd + overlap;
	413	if(total){
	414	// Combined systematc event uncertainty, by weighting with the number of eta-bins of fmd-only, spd-only and the overlap.
	415	sysErrorSquared= (fmdTMath::Power(fmdSysError,2)+ spdTMath::Power(spdSysError,2)+
	416	0.5overlapTMath::Power(fmdSysError,2)+ 0.5overlapTMath::Power(spdSysError,2))/total;
	417	}
418
419
420	if(selectedTrigger){
421	fHist->Fill(c);
422	fHistMC->Fill(cPrimary);
423	fResponseMatrix->Fill(cPrimary, c);
424	fResponseMatrixPlusSys->Fill(cPrimary,(1+TMath::Sqrt(sysErrorSquared))*c);
425	fResponseMatrixMinusSys->Fill(cPrimary,(1-TMath::Sqrt(sysErrorSquared))*c);
426	fResponseMatrixPlus05->Fill(cPrimary, 1.05*c);
427	fResponseMatrixPlus075->Fill(cPrimary, 1.075*c);
428	fResponseMatrixPlus10->Fill(cPrimary, 1.1*c);
429	fResponseMatrixMinus05->Fill(cPrimary, 0.95*c);
430	fResponseMatrixMinus075->Fill(cPrimary, 0.925*c);
431	fResponseMatrixMinus10->Fill(cPrimary, 0.9*c);
432
433	}
434
435	}
436
437
438
439
440	//_____________________________________________________________________
441	//
442	//
443	// EOF