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


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