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

//====================================================================
// 
// Base class for classes that calculate the multiplicity in the
// forward regions event-by-event
// 
// Inputs: 
//   - AliESDEvent 
//
// Outputs: 
//   - AliAODForwardMult 
// 
// Histograms 
//   
// Corrections used 
#include "AliForwardMultiplicityBase.h"
#include "AliForwardCorrectionManager.h"
#include "AliForwardUtil.h"
#include "AliFMDCorrELossFit.h"
#include "AliLog.h"
#include "AliAODHandler.h"
#include "AliInputEventHandler.h"
#include "AliAnalysisManager.h"
#include "AliFMDEventInspector.h"
#include "AliFMDSharingFilter.h"
#include "AliFMDDensityCalculator.h"
#include "AliFMDCorrector.h"
#include "AliFMDHistCollector.h"
#include "AliFMDEventPlaneFinder.h"
#include "AliESDEvent.h"
#include <TROOT.h>
#include <TSystem.h>
#include <TAxis.h>
#include <TProfile.h>
#include <THStack.h>
#include <iostream>
#include <iomanip>

//====================================================================
AliForwardMultiplicityBase::AliForwardMultiplicityBase(const char* name) 
  : AliBaseESDTask(name, "AliForwardMultiplicityBase", 
		   &(AliForwardCorrectionManager::Instance())),
    fEnableLowFlux(false), 
    fStorePerRing(false),
    fHData(0),
    fHistos(),
    fAODFMD(false),
    fAODEP(false),
    fRingSums(),
    fDoTiming(false),
    fHTiming(0)
{
  DGUARD(fDebug, 3,"Named CTOR of AliForwardMultiplicityBase %s",name);
}


//____________________________________________________________________
void
AliForwardMultiplicityBase::SetDebug(Int_t dbg)
{
  // 
  // Set debug level 
  // 
  // Parameters:
  //    dbg debug level
  //
  AliBaseESDTask::        SetDebug(dbg);
  GetSharingFilter()	.SetDebug(dbg);
  GetDensityCalculator().SetDebug(dbg);
  GetCorrections()	.SetDebug(dbg);
  GetHistCollector()	.SetDebug(dbg);
  GetEventPlaneFinder()	.SetDebug(dbg);
}

//____________________________________________________________________
Bool_t
AliForwardMultiplicityBase::Book()
{
  // 
  // Create output objects 
  // 
  //
  DGUARD(fDebug,1,"Create user ouput");
  UInt_t what = AliForwardCorrectionManager::kAll;
  if (!fEnableLowFlux)
    what ^= AliForwardCorrectionManager::kDoubleHit;
  if (!GetCorrections().IsUseVertexBias())
    what ^= AliForwardCorrectionManager::kVertexBias;
  if (!GetCorrections().IsUseAcceptance())
    what ^= AliForwardCorrectionManager::kAcceptance;
  if (!GetCorrections().IsUseMergingEfficiency())
    what ^= AliForwardCorrectionManager::kMergingEfficiency;
  fNeededCorrections = what;

  GetSharingFilter()	.CreateOutputObjects(fList);
  GetDensityCalculator().CreateOutputObjects(fList);
  GetCorrections()	.CreateOutputObjects(fList);
  GetHistCollector()	.CreateOutputObjects(fList);
  GetEventPlaneFinder()	.CreateOutputObjects(fList);

  TAxis tmp(1, 0, 1);
  fHistos.Init(tmp);

  if (fDebug > 1) fDoTiming = true;
  if (fDoTiming) { 
    fHTiming = new TProfile("timing", "Timing of task", 
			    kTimingTotal, 0.5, kTimingTotal+.5);
    fHTiming->SetDirectory(0);
    fHTiming->SetFillColor(kRed+1);
    fHTiming->SetFillStyle(3001);
    fHTiming->SetMarkerStyle(20);
    fHTiming->SetMarkerColor(kBlack);
    fHTiming->SetLineColor(kBlack);
    fHTiming->SetXTitle("Part");
    fHTiming->SetYTitle("#LTt_{part}#GT [CPU]");
    fHTiming->SetStats(0);
    TAxis* xaxis = fHTiming->GetXaxis();
    xaxis->SetBinLabel(kTimingEventInspector,	
		       GetEventInspector()   .GetName());
    xaxis->SetBinLabel(kTimingSharingFilter,	
		       GetSharingFilter()    .GetName());
    xaxis->SetBinLabel(kTimingDensityCalculator,	
		       GetDensityCalculator().GetName());
    xaxis->SetBinLabel(kTimingCorrections,	
		       GetCorrections()      .GetName());
    xaxis->SetBinLabel(kTimingHistCollector,	
		       GetHistCollector()    .GetName());
    xaxis->SetBinLabel(kTimingEventPlaneFinder,	
		       GetEventPlaneFinder() .GetName());
    xaxis->SetBinLabel(kTimingTotal, "Total");
    fList->Add(fHTiming);
  }
  return true;
}
//____________________________________________________________________
void
AliForwardMultiplicityBase::CreateBranches(AliAODHandler* ah)
{
  TObject* obj   = &fAODFMD; ah->AddBranch("AliAODForwardMult", &obj);
  TObject* epobj = &fAODEP;  ah->AddBranch("AliAODForwardEP", &epobj);

  if (!fStorePerRing) return;
  
  AliWarning("Per-ring histograms in AOD\n"
	     "*********************************************************\n"
	     "* For each event 5 additional 2D histogram are stored   *\n"
	     "* in separate branches of the AODs.  This will increase *\n"
	     "* the size of the AODs - proceed with caution           *\n"
	     "*********************************************************");
  TObject* hists[] = { fHistos.fFMD1i, 
		       fHistos.fFMD2i, fHistos.fFMD2o, 
		       fHistos.fFMD3i, fHistos.fFMD3o };
  for (Int_t i = 0; i < 5; i++) { 
    ah->AddBranch("TH2D", &(hists[i]));
  }
}

//____________________________________________________________________
Bool_t
AliForwardMultiplicityBase::PreData(const TAxis& vertex, const TAxis& eta)
{
  // 
  // Initialise the sub objects and stuff.  Called on first event 
  // 
  //
  DGUARD(fDebug,1,"Initialize sub-algorithms");

  // Force this here so we select the proper quality 
  AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();

  UInt_t what = fNeededCorrections;
  // Check that we have the energy loss fits, needed by 
  //   AliFMDSharingFilter 
  //   AliFMDDensityCalculator 
  if (what & AliForwardCorrectionManager::kELossFits && !fcm.GetELossFit()) 
    AliFatal(Form("No energy loss fits"));
  
  // Check that we have the double hit correction - (optionally) used by 
  //  AliFMDDensityCalculator 
  if (what & AliForwardCorrectionManager::kDoubleHit && !fcm.GetDoubleHit()) 
    AliFatal("No double hit corrections"); 

  // Check that we have the secondary maps, needed by 
  //   AliFMDCorrector 
  //   AliFMDHistCollector
  if (what & AliForwardCorrectionManager::kSecondaryMap && 
      !fcm.GetSecondaryMap()) 
    AliFatal("No secondary corrections");

  // Check that we have the vertex bias correction, needed by 
  //   AliFMDCorrector 
  if (what & AliForwardCorrectionManager::kVertexBias && !fcm.GetVertexBias())
    AliFatal("No event vertex bias corrections");

  // Check that we have the merging efficiencies, optionally used by 
  //   AliFMDCorrector 
  if (what & AliForwardCorrectionManager::kMergingEfficiency && 
      !fcm.GetMergingEfficiency()) 
    AliFatal("No merging efficiencies");

  // Check that we have the acceptance correction, needed by 
  //   AliFMDCorrector 
  if (what & AliForwardCorrectionManager::kAcceptance && !fcm.GetAcceptance())
    AliFatal("No acceptance corrections");

  // const AliFMDCorrELossFit* fits = fcm.GetELossFit();
  // fits->CacheBins(GetDensityCalculator().GetMinQuality());

  InitMembers(eta,vertex);
  
  GetDensityCalculator().SetupForData(eta);
  GetSharingFilter()	.SetupForData(eta);
  GetCorrections()	.SetupForData(eta);
  GetHistCollector()	.SetupForData(vertex,eta);

  GetEventPlaneFinder() .SetRunNumber(GetEventInspector().GetRunNumber());
  GetEventPlaneFinder()	.SetupForData(eta);
  
  fAODFMD.SetBit(AliAODForwardMult::kSecondary, 
		 GetCorrections().IsUseSecondaryMap());
  fAODFMD.SetBit(AliAODForwardMult::kVertexBias, 
		 GetCorrections().IsUseVertexBias());
  fAODFMD.SetBit(AliAODForwardMult::kAcceptance, 
		 GetCorrections().IsUseAcceptance());
  fAODFMD.SetBit(AliAODForwardMult::kMergingEfficiency, 
		 GetCorrections().IsUseMergingEfficiency());
  fAODFMD.SetBit(AliAODForwardMult::kSum, 
		 GetHistCollector().GetMergeMethod() == 
		 AliFMDHistCollector::kSum);
  return true;
}

//____________________________________________________________________
void
AliForwardMultiplicityBase::InitMembers(const TAxis& eta, const TAxis& /*pv*/)
{
  fHistos.ReInit(eta);
  fAODFMD.Init(eta);
  fAODEP.Init(eta);
  fRingSums.Init(eta);

  fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
  fHData->SetStats(0);
  fHData->SetDirectory(0);
  fList->Add(fHData);

  TList* rings = new TList;
  rings->SetName("ringSums");
  rings->SetOwner();
  fList->Add(rings);

  rings->Add(fRingSums.Get(1, 'I'));
  rings->Add(fRingSums.Get(2, 'I'));
  rings->Add(fRingSums.Get(2, 'O'));
  rings->Add(fRingSums.Get(3, 'I'));
  rings->Add(fRingSums.Get(3, 'O'));
  fRingSums.Get(1, 'I')->SetMarkerColor(AliForwardUtil::RingColor(1, 'I'));
  fRingSums.Get(2, 'I')->SetMarkerColor(AliForwardUtil::RingColor(2, 'I'));
  fRingSums.Get(2, 'O')->SetMarkerColor(AliForwardUtil::RingColor(2, 'O'));
  fRingSums.Get(3, 'I')->SetMarkerColor(AliForwardUtil::RingColor(3, 'I'));
  fRingSums.Get(3, 'O')->SetMarkerColor(AliForwardUtil::RingColor(3, 'O'));
}
//____________________________________________________________________
Bool_t
AliForwardMultiplicityBase::Finalize()
{
  // 
  // End of job
  // 
  // Parameters:
  //    option Not used 
  //
  DGUARD(fDebug,1,"Processing the merged results");

  TList* list   = fList;
  TList* output = fResults;
  Double_t nTr = 0, nTrVtx = 0, nAcc = 0;
  MakeSimpledNdeta(list, output, nTr, nTrVtx, nAcc);

  EstimatedNdeta(list, output);

  GetSharingFilter()	.Terminate(list,output,Int_t(nTr));
  GetDensityCalculator().Terminate(list,output,Int_t(nTrVtx));
  GetCorrections()	.Terminate(list,output,Int_t(nTrVtx));

  TProfile* timing = static_cast<TProfile*>(list->FindObject("timing"));
  if (timing) { 
    TProfile* p = static_cast<TProfile*>(timing->Clone());
    p->SetDirectory(0);
    p->Scale(100. / p->GetBinContent(p->GetNbinsX()));
    p->SetYTitle("#LTt_{part}#GT/#LTt_{total}#GT [%]");
    p->SetTitle("Relative timing of task");
    output->Add(p);
  }
  return true;
}

//____________________________________________________________________
void
AliForwardMultiplicityBase::EstimatedNdeta(const TList* input, 
					   TList*       output) const
{
  MakeRingdNdeta(input, "ringSums", output, "ringResults");
}

//____________________________________________________________________
Bool_t
AliForwardMultiplicityBase::MakeSimpledNdeta(const TList* input, 
					     TList*       output,
					     Double_t&    nTr, 
					     Double_t&    nTrVtx, 
					     Double_t&    nAcc)
{
  // Get our histograms from the container 
  TH1I* hEventsTr    = 0;
  TH1I* hEventsTrVtx = 0;
  TH1I* hEventsAcc   = 0;
  TH1I* hTriggers    = 0;
  if (!GetEventInspector().FetchHistograms(input, 
					   hEventsTr, 
					   hEventsTrVtx, 
					   hEventsAcc,
					   hTriggers)) { 
    AliError(Form("Didn't get histograms from event selector "
		  "(hEventsTr=%p,hEventsTrVtx=%p,hEventsAcc=%p,hTriggers=%p)", 
		  hEventsTr, hEventsTrVtx, hEventsAcc, hTriggers));
    input->ls();
    return false;
  }
  nTr             = hEventsTr->Integral();
  nTrVtx          = hEventsTrVtx->Integral();
  nAcc            = hEventsAcc->Integral();
  Double_t vtxEff = nTrVtx / nTr;
  TH2D*   hData   = static_cast<TH2D*>(input->FindObject("d2Ndetadphi"));
  if (!hData) { 
    AliError(Form("Couldn't get our summed histogram from output "
		  "list %s (d2Ndetadphi=%p)", input->GetName(), hData));
    input->ls();
    return false;
  }

  Int_t nY      = hData->GetNbinsY();
  TH1D* dNdeta  = hData->ProjectionX("dNdeta",  1,     nY, "e");
  TH1D* dNdeta_ = hData->ProjectionX("dNdeta_", 1,     nY, "e");
  TH1D* norm    = hData->ProjectionX("norm",    0,     0,  "");
  TH1D* phi     = hData->ProjectionX("phi",     nY+1,  nY+1,  "");
  dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
  dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
  dNdeta->SetMarkerColor(kRed+1);
  dNdeta->SetMarkerStyle(20);
  dNdeta->SetDirectory(0);

  dNdeta_->SetTitle("dN_{ch}/d#eta in the forward regions");
  dNdeta_->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
  dNdeta_->SetMarkerColor(kMagenta+1);
  dNdeta_->SetMarkerStyle(21);
  dNdeta_->SetDirectory(0);

  norm->SetTitle("Normalization to  #eta coverage");
  norm->SetYTitle("#eta coverage");
  norm->SetLineColor(kBlue+1);
  norm->SetMarkerColor(kBlue+1);
  norm->SetMarkerStyle(21);
  norm->SetFillColor(kBlue+1);
  norm->SetFillStyle(3005);
  norm->SetDirectory(0);

  phi->SetTitle("Normalization to  #phi acceptance");
  phi->SetYTitle("#phi acceptance");
  phi->SetLineColor(kGreen+1);
  phi->SetMarkerColor(kGreen+1);
  phi->SetMarkerStyle(20);
  phi->SetFillColor(kGreen+1);
  phi->SetFillStyle(3004);
  // phi->Scale(1. / nAcc);
  phi->SetDirectory(0);

  // dNdeta->Divide(norm);
  dNdeta->Divide(phi);
  dNdeta->SetStats(0);
  dNdeta->Scale(vtxEff,	"width");

  dNdeta_->Divide(norm);
  dNdeta_->SetStats(0);
  dNdeta_->Scale(vtxEff, "width");

  output->Add(dNdeta);
  output->Add(dNdeta_);
  output->Add(norm);
  output->Add(phi);

  return true;
}

					     
//____________________________________________________________________
void
AliForwardMultiplicityBase::MakeRingdNdeta(const TList* input, 
					   const char*  inName,
					   TList*       output,
					   const char*  outName,
					   Int_t        style) const
{
  // Make dN/deta for each ring found in the input list.  
  // 
  // A stack of all the dN/deta is also made for easy drawing. 
  // 
  // Note, that the distributions are normalised to the number of
  // observed events only - they should be corrected for 
  DGUARD(fDebug,3,"Make first-shot ring dN/deta");

  if (!input) return;
  TList* list = static_cast<TList*>(input->FindObject(inName));
  if (!list) { 
    AliWarning(Form("No list %s found in %s", inName, input->GetName()));
    return;
  }
  
  TList* out = new TList;
  out->SetName(outName);
  out->SetOwner();
  output->Add(out);

  THStack*     dndetaRings = new THStack("all", "dN/d#eta per ring");
  const char*  names[]     = { "FMD1I", "FMD2I", "FMD2O", "FMD3I", "FMD3O", 0 };
  const char** ptr         = names;
  
  while (*ptr) { 
    TList* thisList = new TList;
    thisList->SetOwner();
    thisList->SetName(*ptr);
    out->Add(thisList);

    TH2D* h = static_cast<TH2D*>(list->FindObject(Form("%s_cache", *ptr)));
    if (!h) { 
      AliWarning(Form("Didn't find %s_cache in %s", *ptr, list->GetName()));
      ptr++;
      continue;
    }
    TH2D* sumPhi = static_cast<TH2D*>(h->Clone("sum_phi"));
    sumPhi->SetDirectory(0);
    thisList->Add(sumPhi);

    TH2D* sumEta = static_cast<TH2D*>(h->Clone("sum_eta"));
    sumEta->SetDirectory(0);
    thisList->Add(sumEta);
    
    Int_t nY   = sumEta->GetNbinsY();
    TH1D* etaCov =static_cast<TH1D*>(h->ProjectionX("etaCov", 0,    0,    ""));
    TH1D* phiAcc =static_cast<TH1D*>(h->ProjectionX("phiAcc", nY+1, nY+1, ""));

    etaCov->SetTitle("Normalization to #eta coverage");
    etaCov->SetYTitle("#eta coverage");
    etaCov->SetMarkerColor(kBlue+1);
    etaCov->SetFillColor(kBlue+1);
    etaCov->SetFillStyle(3005);
    etaCov->SetDirectory(0);
    
    phiAcc->SetTitle("Normalization to #phi acceptance");
    phiAcc->SetYTitle("#phi acceptance");
    phiAcc->SetMarkerColor(kGreen+1);
    phiAcc->SetFillColor(kGreen+1);
    phiAcc->SetFillStyle(3004);
    // phiAcc->Scale(1. / nAcc);
    phiAcc->SetDirectory(0);

    // Double_t s = (etaCov->GetMaximum() > 0 ? 1. / etaCov->GetMaximum() : 1);
    for (Int_t i = 1; i <= sumEta->GetNbinsX(); i++) { 
      for (Int_t j = 1; j <= nY; j++) { 
	Double_t c = sumEta->GetBinContent(i, j);
	Double_t e = sumEta->GetBinError(i, j);
	Double_t a = etaCov->GetBinContent(i);
	Double_t p = phiAcc->GetBinContent(i);
	// Double_t t = p; // * a
	sumEta->SetBinContent(i, j, a <= 0 ? 0 : c / a);
	sumEta->SetBinError(  i, j, a <= 0 ? 0 : e / a);
	sumPhi->SetBinContent(i, j, p <= 0 ? 0 : c / p);
	sumPhi->SetBinError(  i, j, p <= 0 ? 0 : e / p);
      }
    }
    // etaCov->Scale(s);
    // phiAcc->Scale(s);

    TH1D* resPhi  =static_cast<TH1D*>(sumPhi->ProjectionX("dndeta_phi",
							  1,nY,"e"));
    resPhi->SetMarkerStyle(style);
    resPhi->SetDirectory(0);
    resPhi->Scale(1, "width");

    TH1D* resEta  =static_cast<TH1D*>(sumEta->ProjectionX("dndeta_eta",
							  1,nY,"e"));
    resEta->SetMarkerStyle(style);
    resEta->SetDirectory(0);
    resEta->Scale(1, "width");

    thisList->Add(resEta);
    thisList->Add(etaCov);
    thisList->Add(resPhi);
    thisList->Add(phiAcc);
    dndetaRings->Add(resPhi);
    ptr++;
  }
  out->Add(dndetaRings);
}
#define PFB(N,FLAG)				\
  do {									\
    AliForwardUtil::PrintName(N);					\
    std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
  } while(false)
//____________________________________________________________________
void
AliForwardMultiplicityBase::Print(Option_t* option) const
{
  // 
  // Print information 
  // 
  // Parameters:
  //    option Not used
  //
  AliBaseESDTask::Print(option);
  gROOT->IncreaseDirLevel();
  PFB("Enable low flux code", fEnableLowFlux);
  PFB("Store per-ring hists", fStorePerRing);
  PFB("Make timing histogram", fDoTiming);
  // gROOT->IncreaseDirLevel();
  GetSharingFilter()    .Print(option);
  GetDensityCalculator().Print(option);
  GetCorrections()      .Print(option);
  GetHistCollector()    .Print(option);
  GetEventPlaneFinder() .Print(option);
  // gROOT->DecreaseDirLevel();
  gROOT->DecreaseDirLevel();
}


//
// EOF
//
Commit	Line	Data
1a26066e	1	//====================================================================
7984e5f7	2	//
	3	// Base class for classes that calculate the multiplicity in the
	4	// forward regions event-by-event
	5	//
	6	// Inputs:
	7	// - AliESDEvent
	8	//
	9	// Outputs:
	10	// - AliAODForwardMult
	11	//
	12	// Histograms
	13	//
	14	// Corrections used
1a26066e	15	#include "AliForwardMultiplicityBase.h"
1174780f	16	#include "AliForwardCorrectionManager.h"
fb3430ac	17	#include "AliForwardUtil.h"
8449e3e0	18	#include "AliFMDCorrELossFit.h"
1a26066e	19	#include "AliLog.h"
	20	#include "AliAODHandler.h"
	21	#include "AliInputEventHandler.h"
	22	#include "AliAnalysisManager.h"
1174780f	23	#include "AliFMDEventInspector.h"
1174780f	24	#include "AliFMDSharingFilter.h"
1174780f	25	#include "AliFMDDensityCalculator.h"
72cc12cd	26	#include "AliFMDCorrector.h"
1174780f	27	#include "AliFMDHistCollector.h"
2b556440	28	#include "AliFMDEventPlaneFinder.h"
7ec4d843	29	#include "AliESDEvent.h"
1a26066e	30	#include <TROOT.h>
2a276c75	31	#include <TSystem.h>
fb3430ac	32	#include <TAxis.h>
1f7aa5c7	33	#include <TProfile.h>
5bb5d1f6	34	#include <THStack.h>
3c5497d0	35	#include <iostream>
3c5497d0	36	#include <iomanip>
1a26066e	37
1a26066e	38	//====================================================================
19abe41d	39	AliForwardMultiplicityBase::AliForwardMultiplicityBase(const char* name)
c8b1a7db	40	: AliBaseESDTask(name, "AliForwardMultiplicityBase",
c8b1a7db	41	&(AliForwardCorrectionManager::Instance())),
9d05ffeb	42	fEnableLowFlux(false),
8449e3e0	43	fStorePerRing(false),
8449e3e0	44	fHData(0),
	45	fHistos(),
	46	fAODFMD(false),
	47	fAODEP(false),
1f7aa5c7	48	fRingSums(),
1f7aa5c7	49	fDoTiming(false),
c8b1a7db	50	fHTiming(0)
19abe41d	51	{
5ca83fee	52	DGUARD(fDebug, 3,"Named CTOR of AliForwardMultiplicityBase %s",name);
19abe41d	53	}
19abe41d	54
8449e3e0	55
	56	//____________________________________________________________________
	57	void
	58	AliForwardMultiplicityBase::SetDebug(Int_t dbg)
	59	{
	60	//
	61	// Set debug level
	62	//
	63	// Parameters:
	64	// dbg debug level
	65	//
c8b1a7db	66	AliBaseESDTask:: SetDebug(dbg);
8449e3e0	67	GetSharingFilter() .SetDebug(dbg);
	68	GetDensityCalculator().SetDebug(dbg);
	69	GetCorrections() .SetDebug(dbg);
	70	GetHistCollector() .SetDebug(dbg);
	71	GetEventPlaneFinder() .SetDebug(dbg);
	72	}
fb3430ac	73
8449e3e0	74	//____________________________________________________________________
c8b1a7db	75	Bool_t
c8b1a7db	76	AliForwardMultiplicityBase::Book()
8449e3e0	77	{
	78	//
	79	// Create output objects
	80	//
	81	//
	82	DGUARD(fDebug,1,"Create user ouput");
c8b1a7db	83	UInt_t what = AliForwardCorrectionManager::kAll;
	84	if (!fEnableLowFlux)
	85	what ^= AliForwardCorrectionManager::kDoubleHit;
	86	if (!GetCorrections().IsUseVertexBias())
	87	what ^= AliForwardCorrectionManager::kVertexBias;
	88	if (!GetCorrections().IsUseAcceptance())
	89	what ^= AliForwardCorrectionManager::kAcceptance;
	90	if (!GetCorrections().IsUseMergingEfficiency())
	91	what ^= AliForwardCorrectionManager::kMergingEfficiency;
	92	fNeededCorrections = what;
	93
8449e3e0	94	GetSharingFilter() .CreateOutputObjects(fList);
	95	GetDensityCalculator().CreateOutputObjects(fList);
	96	GetCorrections() .CreateOutputObjects(fList);
	97	GetHistCollector() .CreateOutputObjects(fList);
	98	GetEventPlaneFinder() .CreateOutputObjects(fList);
	99
0b7de667	100	TAxis tmp(1, 0, 1);
	101	fHistos.Init(tmp);
	102
1f7aa5c7	103	if (fDebug > 1) fDoTiming = true;
	104	if (fDoTiming) {
	105	fHTiming = new TProfile("timing", "Timing of task",
	106	kTimingTotal, 0.5, kTimingTotal+.5);
	107	fHTiming->SetDirectory(0);
	108	fHTiming->SetFillColor(kRed+1);
	109	fHTiming->SetFillStyle(3001);
	110	fHTiming->SetMarkerStyle(20);
	111	fHTiming->SetMarkerColor(kBlack);
	112	fHTiming->SetLineColor(kBlack);
	113	fHTiming->SetXTitle("Part");
	114	fHTiming->SetYTitle("#LTt_{part}#GT [CPU]");
	115	fHTiming->SetStats(0);
	116	TAxis* xaxis = fHTiming->GetXaxis();
	117	xaxis->SetBinLabel(kTimingEventInspector,
	118	GetEventInspector() .GetName());
	119	xaxis->SetBinLabel(kTimingSharingFilter,
	120	GetSharingFilter() .GetName());
	121	xaxis->SetBinLabel(kTimingDensityCalculator,
	122	GetDensityCalculator().GetName());
	123	xaxis->SetBinLabel(kTimingCorrections,
	124	GetCorrections() .GetName());
	125	xaxis->SetBinLabel(kTimingHistCollector,
	126	GetHistCollector() .GetName());
	127	xaxis->SetBinLabel(kTimingEventPlaneFinder,
	128	GetEventPlaneFinder() .GetName());
	129	xaxis->SetBinLabel(kTimingTotal, "Total");
	130	fList->Add(fHTiming);
	131	}
c8b1a7db	132	return true;
8449e3e0	133	}
	134	//____________________________________________________________________
	135	void
	136	AliForwardMultiplicityBase::CreateBranches(AliAODHandler* ah)
	137	{
c8b1a7db	138	TObject* obj = &fAODFMD; ah->AddBranch("AliAODForwardMult", &obj);
c8b1a7db	139	TObject* epobj = &fAODEP; ah->AddBranch("AliAODForwardEP", &epobj);
8449e3e0	140
8449e3e0	141	if (!fStorePerRing) return;
	142
	143	AliWarning("Per-ring histograms in AOD\n"
	144	"*********************************************************\n"
	145	"* For each event 5 additional 2D histogram are stored *\n"
	146	"* in separate branches of the AODs. This will increase *\n"
	147	"* the size of the AODs - proceed with caution *\n"
	148	"*********************************************************");
	149	TObject* hists[] = { fHistos.fFMD1i,
	150	fHistos.fFMD2i, fHistos.fFMD2o,
	151	fHistos.fFMD3i, fHistos.fFMD3o };
	152	for (Int_t i = 0; i < 5; i++) {
	153	ah->AddBranch("TH2D", &(hists[i]));
	154	}
	155	}
	156
	157	//____________________________________________________________________
	158	Bool_t
c8b1a7db	159	AliForwardMultiplicityBase::PreData(const TAxis& vertex, const TAxis& eta)
8449e3e0	160	{
	161	//
	162	// Initialise the sub objects and stuff. Called on first event
	163	//
	164	//
	165	DGUARD(fDebug,1,"Initialize sub-algorithms");
8449e3e0	166
c8b1a7db	167	// Force this here so we select the proper quality
	168	AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
	169
	170	UInt_t what = fNeededCorrections;
	171	// Check that we have the energy loss fits, needed by
	172	// AliFMDSharingFilter
	173	// AliFMDDensityCalculator
	174	if (what & AliForwardCorrectionManager::kELossFits && !fcm.GetELossFit())
	175	AliFatal(Form("No energy loss fits"));
8449e3e0	176
c8b1a7db	177	// Check that we have the double hit correction - (optionally) used by
	178	// AliFMDDensityCalculator
	179	if (what & AliForwardCorrectionManager::kDoubleHit && !fcm.GetDoubleHit())
	180	AliFatal("No double hit corrections");
	181
	182	// Check that we have the secondary maps, needed by
	183	// AliFMDCorrector
	184	// AliFMDHistCollector
	185	if (what & AliForwardCorrectionManager::kSecondaryMap &&
	186	!fcm.GetSecondaryMap())
	187	AliFatal("No secondary corrections");
	188
	189	// Check that we have the vertex bias correction, needed by
	190	// AliFMDCorrector
	191	if (what & AliForwardCorrectionManager::kVertexBias && !fcm.GetVertexBias())
	192	AliFatal("No event vertex bias corrections");
	193
	194	// Check that we have the merging efficiencies, optionally used by
	195	// AliFMDCorrector
	196	if (what & AliForwardCorrectionManager::kMergingEfficiency &&
	197	!fcm.GetMergingEfficiency())
	198	AliFatal("No merging efficiencies");
	199
	200	// Check that we have the acceptance correction, needed by
	201	// AliFMDCorrector
	202	if (what & AliForwardCorrectionManager::kAcceptance && !fcm.GetAcceptance())
	203	AliFatal("No acceptance corrections");
	204
7095962e CHC	205	// const AliFMDCorrELossFit* fits = fcm.GetELossFit();
7095962e CHC	206	// fits->CacheBins(GetDensityCalculator().GetMinQuality());
c8b1a7db	207
	208	InitMembers(eta,vertex);
	209
c8b1a7db	210	GetDensityCalculator().SetupForData(eta);
7095962e	211	GetSharingFilter() .SetupForData(eta);
c8b1a7db	212	GetCorrections() .SetupForData(eta);
	213	GetHistCollector() .SetupForData(vertex,eta);
	214
	215	GetEventPlaneFinder() .SetRunNumber(GetEventInspector().GetRunNumber());
	216	GetEventPlaneFinder() .SetupForData(eta);
8449e3e0	217
bfab35d9	218	fAODFMD.SetBit(AliAODForwardMult::kSecondary,
	219	GetCorrections().IsUseSecondaryMap());
	220	fAODFMD.SetBit(AliAODForwardMult::kVertexBias,
	221	GetCorrections().IsUseVertexBias());
	222	fAODFMD.SetBit(AliAODForwardMult::kAcceptance,
	223	GetCorrections().IsUseAcceptance());
	224	fAODFMD.SetBit(AliAODForwardMult::kMergingEfficiency,
	225	GetCorrections().IsUseMergingEfficiency());
	226	fAODFMD.SetBit(AliAODForwardMult::kSum,
	227	GetHistCollector().GetMergeMethod() ==
	228	AliFMDHistCollector::kSum);
8449e3e0	229	return true;
	230	}
	231
	232	//____________________________________________________________________
	233	void
c8b1a7db	234	AliForwardMultiplicityBase::InitMembers(const TAxis& eta, const TAxis& /pv/)
8449e3e0	235	{
c8b1a7db	236	fHistos.ReInit(eta);
	237	fAODFMD.Init(eta);
	238	fAODEP.Init(eta);
	239	fRingSums.Init(eta);
8449e3e0	240
	241	fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
	242	fHData->SetStats(0);
	243	fHData->SetDirectory(0);
	244	fList->Add(fHData);
	245
	246	TList* rings = new TList;
	247	rings->SetName("ringSums");
	248	rings->SetOwner();
	249	fList->Add(rings);
	250
	251	rings->Add(fRingSums.Get(1, 'I'));
	252	rings->Add(fRingSums.Get(2, 'I'));
	253	rings->Add(fRingSums.Get(2, 'O'));
	254	rings->Add(fRingSums.Get(3, 'I'));
	255	rings->Add(fRingSums.Get(3, 'O'));
	256	fRingSums.Get(1, 'I')->SetMarkerColor(AliForwardUtil::RingColor(1, 'I'));
	257	fRingSums.Get(2, 'I')->SetMarkerColor(AliForwardUtil::RingColor(2, 'I'));
	258	fRingSums.Get(2, 'O')->SetMarkerColor(AliForwardUtil::RingColor(2, 'O'));
	259	fRingSums.Get(3, 'I')->SetMarkerColor(AliForwardUtil::RingColor(3, 'I'));
	260	fRingSums.Get(3, 'O')->SetMarkerColor(AliForwardUtil::RingColor(3, 'O'));
	261	}
7ec4d843	262	//____________________________________________________________________
7ec4d843	263	Bool_t
c8b1a7db	264	AliForwardMultiplicityBase::Finalize()
8449e3e0	265	{
	266	//
	267	// End of job
	268	//
	269	// Parameters:
	270	// option Not used
	271	//
	272	DGUARD(fDebug,1,"Processing the merged results");
	273
c8b1a7db	274	TList* list = fList;
c8b1a7db	275	TList* output = fResults;
8449e3e0	276	Double_t nTr = 0, nTrVtx = 0, nAcc = 0;
	277	MakeSimpledNdeta(list, output, nTr, nTrVtx, nAcc);
	278
	279	EstimatedNdeta(list, output);
	280
	281	GetSharingFilter() .Terminate(list,output,Int_t(nTr));
	282	GetDensityCalculator().Terminate(list,output,Int_t(nTrVtx));
	283	GetCorrections() .Terminate(list,output,Int_t(nTrVtx));
	284
1f7aa5c7	285	TProfile* timing = static_cast<TProfile*>(list->FindObject("timing"));
	286	if (timing) {
	287	TProfile* p = static_cast<TProfile*>(timing->Clone());
	288	p->SetDirectory(0);
	289	p->Scale(100. / p->GetBinContent(p->GetNbinsX()));
	290	p->SetYTitle("#LTt_{part}#GT/#LTt_{total}#GT [%]");
	291	p->SetTitle("Relative timing of task");
	292	output->Add(p);
	293	}
c8b1a7db	294	return true;
8449e3e0	295	}
	296
	297	//____________________________________________________________________
	298	void
	299	AliForwardMultiplicityBase::EstimatedNdeta(const TList* input,
	300	TList* output) const
	301	{
	302	MakeRingdNdeta(input, "ringSums", output, "ringResults");
	303	}
	304
5ca83fee	305	//____________________________________________________________________
	306	Bool_t
	307	AliForwardMultiplicityBase::MakeSimpledNdeta(const TList* input,
	308	TList* output,
	309	Double_t& nTr,
	310	Double_t& nTrVtx,
	311	Double_t& nAcc)
	312	{
	313	// Get our histograms from the container
	314	TH1I* hEventsTr = 0;
	315	TH1I* hEventsTrVtx = 0;
	316	TH1I* hEventsAcc = 0;
	317	TH1I* hTriggers = 0;
	318	if (!GetEventInspector().FetchHistograms(input,
	319	hEventsTr,
	320	hEventsTrVtx,
	321	hEventsAcc,
	322	hTriggers)) {
	323	AliError(Form("Didn't get histograms from event selector "
	324	"(hEventsTr=%p,hEventsTrVtx=%p,hEventsAcc=%p,hTriggers=%p)",
	325	hEventsTr, hEventsTrVtx, hEventsAcc, hTriggers));
	326	input->ls();
	327	return false;
	328	}
	329	nTr = hEventsTr->Integral();
	330	nTrVtx = hEventsTrVtx->Integral();
	331	nAcc = hEventsAcc->Integral();
	332	Double_t vtxEff = nTrVtx / nTr;
	333	TH2D* hData = static_cast<TH2D*>(input->FindObject("d2Ndetadphi"));
	334	if (!hData) {
	335	AliError(Form("Couldn't get our summed histogram from output "
	336	"list %s (d2Ndetadphi=%p)", input->GetName(), hData));
	337	input->ls();
	338	return false;
	339	}
	340
	341	Int_t nY = hData->GetNbinsY();
	342	TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, nY, "e");
	343	TH1D* dNdeta_ = hData->ProjectionX("dNdeta_", 1, nY, "e");
	344	TH1D* norm = hData->ProjectionX("norm", 0, 0, "");
	345	TH1D* phi = hData->ProjectionX("phi", nY+1, nY+1, "");
	346	dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
	347	dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
	348	dNdeta->SetMarkerColor(kRed+1);
	349	dNdeta->SetMarkerStyle(20);
	350	dNdeta->SetDirectory(0);
	351
	352	dNdeta_->SetTitle("dN_{ch}/d#eta in the forward regions");
	353	dNdeta_->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
	354	dNdeta_->SetMarkerColor(kMagenta+1);
	355	dNdeta_->SetMarkerStyle(21);
	356	dNdeta_->SetDirectory(0);
	357
bfab35d9	358	norm->SetTitle("Normalization to #eta coverage");
5ca83fee	359	norm->SetYTitle("#eta coverage");
bfab35d9	360	norm->SetLineColor(kBlue+1);
5ca83fee	361	norm->SetMarkerColor(kBlue+1);
	362	norm->SetMarkerStyle(21);
	363	norm->SetFillColor(kBlue+1);
	364	norm->SetFillStyle(3005);
	365	norm->SetDirectory(0);
	366
bfab35d9	367	phi->SetTitle("Normalization to #phi acceptance");
5ca83fee	368	phi->SetYTitle("#phi acceptance");
bfab35d9	369	phi->SetLineColor(kGreen+1);
5ca83fee	370	phi->SetMarkerColor(kGreen+1);
	371	phi->SetMarkerStyle(20);
	372	phi->SetFillColor(kGreen+1);
	373	phi->SetFillStyle(3004);
	374	// phi->Scale(1. / nAcc);
	375	phi->SetDirectory(0);
	376
	377	// dNdeta->Divide(norm);
	378	dNdeta->Divide(phi);
	379	dNdeta->SetStats(0);
	380	dNdeta->Scale(vtxEff, "width");
	381
	382	dNdeta_->Divide(norm);
	383	dNdeta_->SetStats(0);
	384	dNdeta_->Scale(vtxEff, "width");
	385
	386	output->Add(dNdeta);
	387	output->Add(dNdeta_);
	388	output->Add(norm);
	389	output->Add(phi);
	390
	391	return true;
	392	}
	393
	394
5bb5d1f6	395	//____________________________________________________________________
	396	void
	397	AliForwardMultiplicityBase::MakeRingdNdeta(const TList* input,
	398	const char* inName,
	399	TList* output,
	400	const char* outName,
	401	Int_t style) const
	402	{
	403	// Make dN/deta for each ring found in the input list.
	404	//
	405	// A stack of all the dN/deta is also made for easy drawing.
	406	//
	407	// Note, that the distributions are normalised to the number of
	408	// observed events only - they should be corrected for
f53fb4f6	409	DGUARD(fDebug,3,"Make first-shot ring dN/deta");
f53fb4f6	410
5bb5d1f6	411	if (!input) return;
	412	TList* list = static_cast<TList*>(input->FindObject(inName));
	413	if (!list) {
	414	AliWarning(Form("No list %s found in %s", inName, input->GetName()));
	415	return;
	416	}
	417
	418	TList* out = new TList;
	419	out->SetName(outName);
	420	out->SetOwner();
	421	output->Add(out);
	422
	423	THStack* dndetaRings = new THStack("all", "dN/d#eta per ring");
	424	const char* names[] = { "FMD1I", "FMD2I", "FMD2O", "FMD3I", "FMD3O", 0 };
	425	const char** ptr = names;
	426
	427	while (*ptr) {
	428	TList* thisList = new TList;
	429	thisList->SetOwner();
	430	thisList->SetName(*ptr);
	431	out->Add(thisList);
	432
	433	TH2D* h = static_cast<TH2D>(list->FindObject(Form("%s_cache", ptr)));
	434	if (!h) {
	435	AliWarning(Form("Didn't find %s_cache in %s", *ptr, list->GetName()));
	436	ptr++;
	437	continue;
	438	}
5ca83fee	439	TH2D* sumPhi = static_cast<TH2D*>(h->Clone("sum_phi"));
	440	sumPhi->SetDirectory(0);
	441	thisList->Add(sumPhi);
	442
	443	TH2D* sumEta = static_cast<TH2D*>(h->Clone("sum_eta"));
	444	sumEta->SetDirectory(0);
	445	thisList->Add(sumEta);
5bb5d1f6	446
5ca83fee	447	Int_t nY = sumEta->GetNbinsY();
	448	TH1D* etaCov =static_cast<TH1D*>(h->ProjectionX("etaCov", 0, 0, ""));
	449	TH1D* phiAcc =static_cast<TH1D*>(h->ProjectionX("phiAcc", nY+1, nY+1, ""));
5bb5d1f6	450
5ca83fee	451	etaCov->SetTitle("Normalization to #eta coverage");
	452	etaCov->SetYTitle("#eta coverage");
	453	etaCov->SetMarkerColor(kBlue+1);
	454	etaCov->SetFillColor(kBlue+1);
	455	etaCov->SetFillStyle(3005);
	456	etaCov->SetDirectory(0);
65abd48b	457
5ca83fee	458	phiAcc->SetTitle("Normalization to #phi acceptance");
	459	phiAcc->SetYTitle("#phi acceptance");
	460	phiAcc->SetMarkerColor(kGreen+1);
	461	phiAcc->SetFillColor(kGreen+1);
	462	phiAcc->SetFillStyle(3004);
	463	// phiAcc->Scale(1. / nAcc);
	464	phiAcc->SetDirectory(0);
	465
	466	// Double_t s = (etaCov->GetMaximum() > 0 ? 1. / etaCov->GetMaximum() : 1);
	467	for (Int_t i = 1; i <= sumEta->GetNbinsX(); i++) {
	468	for (Int_t j = 1; j <= nY; j++) {
	469	Double_t c = sumEta->GetBinContent(i, j);
	470	Double_t e = sumEta->GetBinError(i, j);
	471	Double_t a = etaCov->GetBinContent(i);
	472	Double_t p = phiAcc->GetBinContent(i);
	473	// Double_t t = p; // * a
	474	sumEta->SetBinContent(i, j, a <= 0 ? 0 : c / a);
	475	sumEta->SetBinError( i, j, a <= 0 ? 0 : e / a);
	476	sumPhi->SetBinContent(i, j, p <= 0 ? 0 : c / p);
	477	sumPhi->SetBinError( i, j, p <= 0 ? 0 : e / p);
	478	}
65abd48b	479	}
5ca83fee	480	// etaCov->Scale(s);
	481	// phiAcc->Scale(s);
	482
	483	TH1D* resPhi =static_cast<TH1D*>(sumPhi->ProjectionX("dndeta_phi",
	484	1,nY,"e"));
	485	resPhi->SetMarkerStyle(style);
	486	resPhi->SetDirectory(0);
	487	resPhi->Scale(1, "width");
	488
	489	TH1D* resEta =static_cast<TH1D*>(sumEta->ProjectionX("dndeta_eta",
	490	1,nY,"e"));
	491	resEta->SetMarkerStyle(style);
	492	resEta->SetDirectory(0);
	493	resEta->Scale(1, "width");
	494
	495	thisList->Add(resEta);
	496	thisList->Add(etaCov);
	497	thisList->Add(resPhi);
	498	thisList->Add(phiAcc);
	499	dndetaRings->Add(resPhi);
5bb5d1f6	500	ptr++;
	501	}
	502	out->Add(dndetaRings);
	503	}
c8b1a7db	504	#define PFB(N,FLAG) \
	505	do { \
	506	AliForwardUtil::PrintName(N); \
	507	std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
	508	} while(false)
1a26066e	509	//____________________________________________________________________
	510	void
	511	AliForwardMultiplicityBase::Print(Option_t* option) const
	512	{
7984e5f7	513	//
	514	// Print information
	515	//
	516	// Parameters:
	517	// option Not used
	518	//
c8b1a7db	519	AliBaseESDTask::Print(option);
1174780f	520	gROOT->IncreaseDirLevel();
c8b1a7db	521	PFB("Enable low flux code", fEnableLowFlux);
	522	PFB("Store per-ring hists", fStorePerRing);
	523	PFB("Make timing histogram", fDoTiming);
	524	// gROOT->IncreaseDirLevel();
1174780f	525	GetSharingFilter() .Print(option);
	526	GetDensityCalculator().Print(option);
	527	GetCorrections() .Print(option);
	528	GetHistCollector() .Print(option);
2b556440	529	GetEventPlaneFinder() .Print(option);
c8b1a7db	530	// gROOT->DecreaseDirLevel();
1174780f	531	gROOT->DecreaseDirLevel();
1a26066e	532	}
1a26066e	533
1174780f	534
1a26066e	535	//
	536	// EOF
	537	//