[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliForwardMultiplicityTask.cxx

#include "AliForwardMultiplicityTask.h"
#include "AliTriggerAnalysis.h"
#include "AliPhysicsSelection.h"
#include "AliLog.h"
// #include "AliFMDAnaParameters.h"
#include "AliESDEvent.h"
#include "AliAODHandler.h"
#include "AliMultiplicity.h"
#include "AliInputEventHandler.h"
#include "AliForwardCorrectionManager.h"
#include "AliAnalysisManager.h"
#include <TH1.h>
#include <TDirectory.h>
#include <TTree.h>
#include <TROOT.h>
#include <iostream>
#include <iomanip>

//====================================================================
AliForwardMultiplicityTask::AliForwardMultiplicityTask()
  : AliAnalysisTaskSE(),
    fEnableLowFlux(true),
    fHData(0),
    fFirstEvent(true),
    fESDFMD(),
    fHistos(),
    fAODFMD(),
    fEventInspector(),
    fEnergyFitter(),
    fSharingFilter(),
    fDensityCalculator(),
    fCorrections(),
    fHistCollector(),
    fList(0)
{
}

//____________________________________________________________________
AliForwardMultiplicityTask::AliForwardMultiplicityTask(const char* name)
  : AliAnalysisTaskSE(name), 
    fEnableLowFlux(true),
    fHData(0),
    fFirstEvent(true),
    fESDFMD(),
    fHistos(),
    fAODFMD(kTRUE),
    fEventInspector("event"),
    fEnergyFitter("energy"),
    fSharingFilter("sharing"), 
    fDensityCalculator("density"),
    fCorrections("corrections"),
    fHistCollector("collector"),
    fList(0)
{
  DefineOutput(1, TList::Class());
}

//____________________________________________________________________
AliForwardMultiplicityTask::AliForwardMultiplicityTask(const AliForwardMultiplicityTask& o)
  : AliAnalysisTaskSE(o),
    fEnableLowFlux(o.fEnableLowFlux),
    fHData(o.fHData),
    fFirstEvent(true),
    fESDFMD(o.fESDFMD),
    fHistos(o.fHistos),
    fAODFMD(o.fAODFMD),
    fEventInspector(o.fEventInspector),
    fEnergyFitter(o.fEnergyFitter),
    fSharingFilter(o.fSharingFilter),
    fDensityCalculator(o.fDensityCalculator),
    fCorrections(o.fCorrections),
    fHistCollector(o.fHistCollector),
    fList(o.fList) 
{
  DefineOutput(1, TList::Class());
}

//____________________________________________________________________
AliForwardMultiplicityTask&
AliForwardMultiplicityTask::operator=(const AliForwardMultiplicityTask& o)
{
  AliAnalysisTaskSE::operator=(o);

  fEnableLowFlux     = o.fEnableLowFlux;
  fHData             = o.fHData;
  fFirstEvent        = o.fFirstEvent;
  fEventInspector    = o.fEventInspector;
  fEnergyFitter      = o.fEnergyFitter;
  fSharingFilter     = o.fSharingFilter;
  fDensityCalculator = o.fDensityCalculator;
  fCorrections       = o.fCorrections;
  fHistCollector     = o.fHistCollector;
  fHistos            = o.fHistos;
  fAODFMD            = o.fAODFMD;
  fList              = o.fList;

  return *this;
}

//____________________________________________________________________
void
AliForwardMultiplicityTask::SetDebug(Int_t dbg)
{
  fEventInspector.SetDebug(dbg);
  fEnergyFitter.SetDebug(dbg);
  fSharingFilter.SetDebug(dbg);
  fDensityCalculator.SetDebug(dbg);
  fCorrections.SetDebug(dbg);
  fHistCollector.SetDebug(dbg);
}

//____________________________________________________________________
void
AliForwardMultiplicityTask::Init()
{
  fFirstEvent = true;
}

//____________________________________________________________________
void
AliForwardMultiplicityTask::InitializeSubs()
{

  // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  // pars->Init(kTRUE);

  AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
  fcm.Init(fEventInspector.GetCollisionSystem(), 
	   fEventInspector.GetEnergy(),
	   fEventInspector.GetField());
  // Check that we have the energy loss fits, needed by 
  //   AliFMDSharingFilter 
  //   AliFMDDensityCalculator 
  if (!fcm.GetELossFit()) { 
    AliFatal(Form("No energy loss fits"));
    return;
  }
  // Check that we have the double hit correction - (optionally) used by 
  //  AliFMDDensityCalculator 
  if (!fcm.GetDoubleHit()) {
    AliWarning("No double hit corrections"); 
  }
  // Check that we have the secondary maps, needed by 
  //   AliFMDCorrections 
  //   AliFMDHistCollector
  if (!fcm.GetSecondaryMap()) {
    AliFatal("No secondary corrections");
    return;
  }
  // Check that we have the vertex bias correction, needed by 
  //   AliFMDCorrections 
  if (!fcm.GetVertexBias()) { 
    AliFatal("No event vertex bias corrections");
    return;
  }
  // Check that we have the merging efficiencies, optionally used by 
  //   AliFMDCorrections 
  if (!fcm.GetMergingEfficiency()) {
    AliWarning("No merging efficiencies");
  }


  const TAxis* pe = fcm.GetEtaAxis();
  const TAxis* pv = fcm.GetVertexAxis();
  if (!pe) AliFatal("No eta axis defined");
  if (!pv) AliFatal("No vertex axis defined");

  // TAxis e(pars->GetNetaBins(),  pars->GetEtaMin(),  pars->GetEtaMax());
  // TAxis v(pars->GetNvtxBins(), -pars->GetVtxCutZ(), pars->GetVtxCutZ());
  // pv->Dump();

  fHistos.Init(*pe);
  fAODFMD.Init(*pe);

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

  fEnergyFitter.Init(*pe);
  fEventInspector.Init(*pv);
  fHistCollector.Init(*pv);

  this->Print();
}

//____________________________________________________________________
void
AliForwardMultiplicityTask::UserCreateOutputObjects()
{
  fList = new TList;

  AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
  AliAODHandler*      ah = 
    dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
  if (!ah) AliFatal("No AOD output handler set in analysis manager");
    
    
  TObject* obj = &fAODFMD;
  ah->AddBranch("AliAODForwardMult", &obj);

  fEventInspector.DefineOutput(fList);
  fEnergyFitter.DefineOutput(fList);
  fSharingFilter.DefineOutput(fList);
  fDensityCalculator.DefineOutput(fList);
  fCorrections.DefineOutput(fList);

  PostData(1, fList);
}
//____________________________________________________________________
void
AliForwardMultiplicityTask::UserExec(Option_t*)
{
  // static Int_t cnt = 0;
  // cnt++;
  // Get the input data 
  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
  // AliInfo(Form("Event # %6d (esd=%p)", cnt, esd));
  if (!esd) { 
    AliWarning("No ESD event found for input event");
    return;
  }

  // On the first event, initialize the parameters 
  if (fFirstEvent && esd->GetESDRun()) { 
    fEventInspector.ReadRunDetails(esd);
    
    AliInfo(Form("Initializing with parameters from the ESD:\n"
		 "         AliESDEvent::GetBeamEnergy()   ->%f\n"
		 "         AliESDEvent::GetBeamType()     ->%s\n"
		 "         AliESDEvent::GetCurrentL3()    ->%f\n"
		 "         AliESDEvent::GetMagneticField()->%f\n"
		 "         AliESDEvent::GetRunNumber()    ->%d\n",
		 esd->GetBeamEnergy(), 
		 esd->GetBeamType(),
		 esd->GetCurrentL3(), 
		 esd->GetMagneticField(),
		 esd->GetRunNumber()));

	      
    // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
    // pars->SetParametersFromESD(esd);
    // pars->PrintStatus();
    fFirstEvent = false;

    InitializeSubs();
  }
  // Clear stuff 
  fHistos.Clear();
  fESDFMD.Clear();
  fAODFMD.Clear();

  Bool_t   lowFlux  = kFALSE;
  UInt_t   triggers = 0;
  UShort_t ivz      = 0;
  Double_t vz       = 0;
  UInt_t   found    = fEventInspector.Process(esd, triggers, lowFlux, ivz, vz);
  if (found & AliFMDEventInspector::kNoEvent)    return;
  if (found & AliFMDEventInspector::kNoTriggers) return;
 
  // Set trigger bits, and mark this event for storage 
  fAODFMD.SetTriggerBits(triggers);
  MarkEventForStore();

  if (found & AliFMDEventInspector::kNoSPD)     return;
  if (found & AliFMDEventInspector::kNoFMD)     return;
  if (found & AliFMDEventInspector::kNoVertex)  return;
  fAODFMD.SetIpZ(vz);

  if (found & AliFMDEventInspector::kBadVertex) return;

  // We we do not want to use low flux specific code, we disable it here. 
  if (!fEnableLowFlux) lowFlux = false;

  // Get FMD data 
  AliESDFMD* esdFMD = esd->GetFMDData();
  // Apply the sharing filter (or hit merging or clustering if you like)
  if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD)) { 
    AliWarning("Sharing filter failed!");
    return;
  }

  // Do the energy stuff 
  if (!fEnergyFitter.Accumulate(*esdFMD, triggers & AliAODForwardMult::kEmpty)){
    AliWarning("Energy fitter failed");
    return;
  }

  // Calculate the inclusive charged particle density 
  if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux)) { 
    AliWarning("Density calculator failed!");
    return;
  }
  
  // Do the secondary and other corrections. 
  if (!fCorrections.Correct(fHistos, ivz)) { 
    AliWarning("Corrections failed");
    return;
  }

  if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) {
    AliWarning("Histogram collector failed");
    return;
  }

  if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
    fHData->Add(&(fAODFMD.GetHistogram()));

  PostData(1, fList);
}

//____________________________________________________________________
void
AliForwardMultiplicityTask::Terminate(Option_t*)
{
  TList* list = dynamic_cast<TList*>(GetOutputData(1));
  if (!list) {
    AliError(Form("No output list defined (%p)", GetOutputData(1)));
    if (GetOutputData(1)) GetOutputData(1)->Print();
    return;
  }
  
  // Get our histograms from the container 
  TH1I* hEventsTr    = 0;//static_cast<TH1I*>(list->FindObject("nEventsTr"));
  TH1I* hEventsTrVtx = 0;//static_cast<TH1I*>(list->FindObject("nEventsTrVtx"));
  TH1I* hTriggers    = 0;
  if (!fEventInspector.FetchHistograms(list, hEventsTr, 
				       hEventsTrVtx, hTriggers)) { 
    AliError(Form("Didn't get histograms from event selector "
		  "(hEventsTr=%p,hEventsTrVtx=%p)", 
		  hEventsTr, hEventsTrVtx));
    list->ls();
    return;
  }

  TH2D* hData        = static_cast<TH2D*>(list->FindObject("d2Ndetadphi"));
  if (!hData) { 
    AliError(Form("Couldn't get our summed histogram from output "
		  "list %s (d2Ndetadphi=%p)", list->GetName(), hData));
    list->ls();
    return;
  }
  
  // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
  TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
  TH1D* norm   = hData->ProjectionX("norm",   0,  1,  "");
  dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
  dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
  dNdeta->Divide(norm);
  dNdeta->SetStats(0);
  dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
		"width");
  list->Add(dNdeta);
  list->Add(norm);

  fEnergyFitter.Fit(list);
  fSharingFilter.ScaleHistograms(list,hEventsTr->Integral());
  fDensityCalculator.ScaleHistograms(list,hEventsTrVtx->Integral());
  fCorrections.ScaleHistograms(list,hEventsTrVtx->Integral());
}

//____________________________________________________________________
void
AliForwardMultiplicityTask::MarkEventForStore() const
{
  // Make sure the AOD tree is filled 
  AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
  AliAODHandler*      ah = 
    dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
  if (!ah)  
    AliFatal("No AOD output handler set in analysis manager");

  ah->SetFillAOD(kTRUE);
}

//____________________________________________________________________
void
AliForwardMultiplicityTask::Print(Option_t* option) const
{
  std::cout << "AliForwardMultiplicityTask: " << GetName() << "\n" 
	    << "  Enable low flux code:   " << (fEnableLowFlux ? "yes" : "no")
	    << std::endl;
  gROOT->IncreaseDirLevel();
  fEventInspector   .Print(option);
  fEnergyFitter     .Print(option);    
  fSharingFilter    .Print(option);
  fDensityCalculator.Print(option);
  fCorrections      .Print(option);
  fHistCollector    .Print(option);
  gROOT->DecreaseDirLevel();
}

//
// EOF
//
Commit	Line	Data
0bd4b00f	1	#include "AliForwardMultiplicityTask.h"
7e4038b5	2	#include "AliTriggerAnalysis.h"
	3	#include "AliPhysicsSelection.h"
	4	#include "AliLog.h"
0bd4b00f	5	// #include "AliFMDAnaParameters.h"
7e4038b5	6	#include "AliESDEvent.h"
	7	#include "AliAODHandler.h"
	8	#include "AliMultiplicity.h"
	9	#include "AliInputEventHandler.h"
0bd4b00f	10	#include "AliForwardCorrectionManager.h"
0bd4b00f	11	#include "AliAnalysisManager.h"
7e4038b5	12	#include <TH1.h>
	13	#include <TDirectory.h>
	14	#include <TTree.h>
0bd4b00f	15	#include <TROOT.h>
	16	#include <iostream>
	17	#include <iomanip>
7e4038b5	18
7e4038b5	19	//====================================================================
0bd4b00f	20	AliForwardMultiplicityTask::AliForwardMultiplicityTask()
7e4038b5	21	: AliAnalysisTaskSE(),
0bd4b00f	22	fEnableLowFlux(true),
7e4038b5	23	fHData(0),
7e4038b5	24	fFirstEvent(true),
7e4038b5	25	fESDFMD(),
	26	fHistos(),
	27	fAODFMD(),
fffea31d	28	fEventInspector(),
fffea31d	29	fEnergyFitter(),
7e4038b5	30	fSharingFilter(),
	31	fDensityCalculator(),
	32	fCorrections(),
	33	fHistCollector(),
fffea31d	34	fList(0)
7e4038b5	35	{
	36	}
	37
	38	//____________________________________________________________________
0bd4b00f	39	AliForwardMultiplicityTask::AliForwardMultiplicityTask(const char* name)
7e4038b5	40	: AliAnalysisTaskSE(name),
0bd4b00f	41	fEnableLowFlux(true),
7e4038b5	42	fHData(0),
7e4038b5	43	fFirstEvent(true),
7e4038b5	44	fESDFMD(),
	45	fHistos(),
	46	fAODFMD(kTRUE),
fffea31d	47	fEventInspector("event"),
fffea31d	48	fEnergyFitter("energy"),
7e4038b5	49	fSharingFilter("sharing"),
	50	fDensityCalculator("density"),
	51	fCorrections("corrections"),
	52	fHistCollector("collector"),
fffea31d	53	fList(0)
7e4038b5	54	{
7e4038b5	55	DefineOutput(1, TList::Class());
7e4038b5	56	}
	57
	58	//____________________________________________________________________
0bd4b00f	59	AliForwardMultiplicityTask::AliForwardMultiplicityTask(const AliForwardMultiplicityTask& o)
7e4038b5	60	: AliAnalysisTaskSE(o),
0bd4b00f	61	fEnableLowFlux(o.fEnableLowFlux),
7e4038b5	62	fHData(o.fHData),
7e4038b5	63	fFirstEvent(true),
7e4038b5	64	fESDFMD(o.fESDFMD),
	65	fHistos(o.fHistos),
	66	fAODFMD(o.fAODFMD),
fffea31d	67	fEventInspector(o.fEventInspector),
fffea31d	68	fEnergyFitter(o.fEnergyFitter),
7e4038b5	69	fSharingFilter(o.fSharingFilter),
	70	fDensityCalculator(o.fDensityCalculator),
	71	fCorrections(o.fCorrections),
	72	fHistCollector(o.fHistCollector),
fffea31d	73	fList(o.fList)
7e4038b5	74	{
0bd4b00f	75	DefineOutput(1, TList::Class());
7e4038b5	76	}
	77
	78	//____________________________________________________________________
0bd4b00f	79	AliForwardMultiplicityTask&
0bd4b00f	80	AliForwardMultiplicityTask::operator=(const AliForwardMultiplicityTask& o)
7e4038b5	81	{
fffea31d	82	AliAnalysisTaskSE::operator=(o);
fffea31d	83
0bd4b00f	84	fEnableLowFlux = o.fEnableLowFlux;
7e4038b5	85	fHData = o.fHData;
7e4038b5	86	fFirstEvent = o.fFirstEvent;
fffea31d	87	fEventInspector = o.fEventInspector;
fffea31d	88	fEnergyFitter = o.fEnergyFitter;
7e4038b5	89	fSharingFilter = o.fSharingFilter;
	90	fDensityCalculator = o.fDensityCalculator;
	91	fCorrections = o.fCorrections;
	92	fHistCollector = o.fHistCollector;
	93	fHistos = o.fHistos;
	94	fAODFMD = o.fAODFMD;
	95	fList = o.fList;
7e4038b5	96
	97	return *this;
	98	}
	99
fffea31d	100	//____________________________________________________________________
fffea31d	101	void
0bd4b00f	102	AliForwardMultiplicityTask::SetDebug(Int_t dbg)
fffea31d	103	{
	104	fEventInspector.SetDebug(dbg);
	105	fEnergyFitter.SetDebug(dbg);
	106	fSharingFilter.SetDebug(dbg);
	107	fDensityCalculator.SetDebug(dbg);
	108	fCorrections.SetDebug(dbg);
	109	fHistCollector.SetDebug(dbg);
	110	}
	111
7e4038b5	112	//____________________________________________________________________
7e4038b5	113	void
0bd4b00f	114	AliForwardMultiplicityTask::Init()
7e4038b5	115	{
	116	fFirstEvent = true;
	117	}
	118
	119	//____________________________________________________________________
	120	void
0bd4b00f	121	AliForwardMultiplicityTask::InitializeSubs()
7e4038b5	122	{
7e4038b5	123
0bd4b00f	124	// AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
	125	// pars->Init(kTRUE);
	126
	127	AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
	128	fcm.Init(fEventInspector.GetCollisionSystem(),
	129	fEventInspector.GetEnergy(),
	130	fEventInspector.GetField());
	131	// Check that we have the energy loss fits, needed by
	132	// AliFMDSharingFilter
	133	// AliFMDDensityCalculator
	134	if (!fcm.GetELossFit()) {
	135	AliFatal(Form("No energy loss fits"));
	136	return;
	137	}
	138	// Check that we have the double hit correction - (optionally) used by
	139	// AliFMDDensityCalculator
	140	if (!fcm.GetDoubleHit()) {
	141	AliWarning("No double hit corrections");
	142	}
	143	// Check that we have the secondary maps, needed by
	144	// AliFMDCorrections
	145	// AliFMDHistCollector
	146	if (!fcm.GetSecondaryMap()) {
	147	AliFatal("No secondary corrections");
	148	return;
	149	}
	150	// Check that we have the vertex bias correction, needed by
	151	// AliFMDCorrections
	152	if (!fcm.GetVertexBias()) {
	153	AliFatal("No event vertex bias corrections");
	154	return;
	155	}
	156	// Check that we have the merging efficiencies, optionally used by
	157	// AliFMDCorrections
	158	if (!fcm.GetMergingEfficiency()) {
	159	AliWarning("No merging efficiencies");
	160	}
	161
	162
	163	const TAxis* pe = fcm.GetEtaAxis();
	164	const TAxis* pv = fcm.GetVertexAxis();
	165	if (!pe) AliFatal("No eta axis defined");
	166	if (!pv) AliFatal("No vertex axis defined");
fffea31d	167
0bd4b00f	168	// TAxis e(pars->GetNetaBins(), pars->GetEtaMin(), pars->GetEtaMax());
	169	// TAxis v(pars->GetNvtxBins(), -pars->GetVtxCutZ(), pars->GetVtxCutZ());
	170	// pv->Dump();
	171
	172	fHistos.Init(*pe);
	173	fAODFMD.Init(*pe);
7e4038b5	174
	175	fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
	176	fHData->SetStats(0);
	177	fHData->SetDirectory(0);
9d99b0dd	178	fList->Add(fHData);
9d99b0dd	179
0bd4b00f	180	fEnergyFitter.Init(*pe);
	181	fEventInspector.Init(*pv);
	182	fHistCollector.Init(*pv);
	183
	184	this->Print();
7e4038b5	185	}
	186
	187	//____________________________________________________________________
	188	void
0bd4b00f	189	AliForwardMultiplicityTask::UserCreateOutputObjects()
7e4038b5	190	{
	191	fList = new TList;
	192
	193	AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
0bd4b00f	194	AliAODHandler* ah =
0bd4b00f	195	dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
fffea31d	196	if (!ah) AliFatal("No AOD output handler set in analysis manager");
7e4038b5	197
	198
	199	TObject* obj = &fAODFMD;
	200	ah->AddBranch("AliAODForwardMult", &obj);
	201
fffea31d	202	fEventInspector.DefineOutput(fList);
fffea31d	203	fEnergyFitter.DefineOutput(fList);
9d99b0dd	204	fSharingFilter.DefineOutput(fList);
	205	fDensityCalculator.DefineOutput(fList);
	206	fCorrections.DefineOutput(fList);
0bd4b00f	207
0bd4b00f	208	PostData(1, fList);
7e4038b5	209	}
	210	//____________________________________________________________________
	211	void
0bd4b00f	212	AliForwardMultiplicityTask::UserExec(Option_t*)
7e4038b5	213	{
0bd4b00f	214	// static Int_t cnt = 0;
0bd4b00f	215	// cnt++;
7e4038b5	216	// Get the input data
7e4038b5	217	AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
0bd4b00f	218	// AliInfo(Form("Event # %6d (esd=%p)", cnt, esd));
7e4038b5	219	if (!esd) {
	220	AliWarning("No ESD event found for input event");
	221	return;
	222	}
	223
7e4038b5	224	// On the first event, initialize the parameters
fffea31d	225	if (fFirstEvent && esd->GetESDRun()) {
0bd4b00f	226	fEventInspector.ReadRunDetails(esd);
0bd4b00f	227
fffea31d	228	AliInfo(Form("Initializing with parameters from the ESD:\n"
	229	" AliESDEvent::GetBeamEnergy() ->%f\n"
	230	" AliESDEvent::GetBeamType() ->%s\n"
	231	" AliESDEvent::GetCurrentL3() ->%f\n"
	232	" AliESDEvent::GetMagneticField()->%f\n"
	233	" AliESDEvent::GetRunNumber() ->%d\n",
	234	esd->GetBeamEnergy(),
	235	esd->GetBeamType(),
	236	esd->GetCurrentL3(),
	237	esd->GetMagneticField(),
	238	esd->GetRunNumber()));
0bd4b00f	239
	240
	241
	242	// AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
	243	// pars->SetParametersFromESD(esd);
	244	// pars->PrintStatus();
7e4038b5	245	fFirstEvent = false;
	246
	247	InitializeSubs();
	248	}
	249	// Clear stuff
	250	fHistos.Clear();
	251	fESDFMD.Clear();
	252	fAODFMD.Clear();
	253
fffea31d	254	Bool_t lowFlux = kFALSE;
fffea31d	255	UInt_t triggers = 0;
0bd4b00f	256	UShort_t ivz = 0;
fffea31d	257	Double_t vz = 0;
	258	UInt_t found = fEventInspector.Process(esd, triggers, lowFlux, ivz, vz);
	259	if (found & AliFMDEventInspector::kNoEvent) return;
	260	if (found & AliFMDEventInspector::kNoTriggers) return;
	261
	262	// Set trigger bits, and mark this event for storage
9d99b0dd	263	fAODFMD.SetTriggerBits(triggers);
7e4038b5	264	MarkEventForStore();
7e4038b5	265
fffea31d	266	if (found & AliFMDEventInspector::kNoSPD) return;
	267	if (found & AliFMDEventInspector::kNoFMD) return;
	268	if (found & AliFMDEventInspector::kNoVertex) return;
	269	fAODFMD.SetIpZ(vz);
7e4038b5	270
fffea31d	271	if (found & AliFMDEventInspector::kBadVertex) return;
7e4038b5	272
0bd4b00f	273	// We we do not want to use low flux specific code, we disable it here.
	274	if (!fEnableLowFlux) lowFlux = false;
	275
fffea31d	276	// Get FMD data
	277	AliESDFMD* esdFMD = esd->GetFMDData();
	278	// Apply the sharing filter (or hit merging or clustering if you like)
	279	if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD)) {
	280	AliWarning("Sharing filter failed!");
7e4038b5	281	return;
	282	}
	283
fffea31d	284	// Do the energy stuff
	285	if (!fEnergyFitter.Accumulate(*esdFMD, triggers & AliAODForwardMult::kEmpty)){
	286	AliWarning("Energy fitter failed");
7e4038b5	287	return;
	288	}
	289
	290	// Calculate the inclusive charged particle density
	291	if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux)) {
	292	AliWarning("Density calculator failed!");
	293	return;
	294	}
	295
	296	// Do the secondary and other corrections.
	297	if (!fCorrections.Correct(fHistos, ivz)) {
	298	AliWarning("Corrections failed");
	299	return;
	300	}
	301
	302	if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) {
	303	AliWarning("Histogram collector failed");
	304	return;
	305	}
	306
	307	if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
	308	fHData->Add(&(fAODFMD.GetHistogram()));
2d68d438	309
2d68d438	310	PostData(1, fList);
7e4038b5	311	}
	312
	313	//____________________________________________________________________
	314	void
0bd4b00f	315	AliForwardMultiplicityTask::Terminate(Option_t*)
7e4038b5	316	{
	317	TList* list = dynamic_cast<TList*>(GetOutputData(1));
	318	if (!list) {
2d68d438	319	AliError(Form("No output list defined (%p)", GetOutputData(1)));
2d68d438	320	if (GetOutputData(1)) GetOutputData(1)->Print();
7e4038b5	321	return;
7e4038b5	322	}
9d99b0dd	323
9d99b0dd	324	// Get our histograms from the container
fffea31d	325	TH1I* hEventsTr = 0;//static_cast<TH1I*>(list->FindObject("nEventsTr"));
	326	TH1I* hEventsTrVtx = 0;//static_cast<TH1I*>(list->FindObject("nEventsTrVtx"));
	327	TH1I* hTriggers = 0;
	328	if (!fEventInspector.FetchHistograms(list, hEventsTr,
	329	hEventsTrVtx, hTriggers)) {
	330	AliError(Form("Didn't get histograms from event selector "
	331	"(hEventsTr=%p,hEventsTrVtx=%p)",
	332	hEventsTr, hEventsTrVtx));
	333	list->ls();
	334	return;
	335	}
	336
9d99b0dd	337	TH2D* hData = static_cast<TH2D*>(list->FindObject("d2Ndetadphi"));
fffea31d	338	if (!hData) {
	339	AliError(Form("Couldn't get our summed histogram from output "
	340	"list %s (d2Ndetadphi=%p)", list->GetName(), hData));
2d68d438	341	list->ls();
	342	return;
	343	}
9d99b0dd	344
7e4038b5	345	// TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
9d99b0dd	346	TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
0bd4b00f	347	TH1D* norm = hData->ProjectionX("norm", 0, 1, "");
7e4038b5	348	dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
	349	dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
	350	dNdeta->Divide(norm);
	351	dNdeta->SetStats(0);
9d99b0dd	352	dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
7e4038b5	353	"width");
7e4038b5	354	list->Add(dNdeta);
0bd4b00f	355	list->Add(norm);
fffea31d	356
fffea31d	357	fEnergyFitter.Fit(list);
9d99b0dd	358	fSharingFilter.ScaleHistograms(list,hEventsTr->Integral());
	359	fDensityCalculator.ScaleHistograms(list,hEventsTrVtx->Integral());
	360	fCorrections.ScaleHistograms(list,hEventsTrVtx->Integral());
7e4038b5	361	}
	362
	363	//____________________________________________________________________
	364	void
0bd4b00f	365	AliForwardMultiplicityTask::MarkEventForStore() const
7e4038b5	366	{
	367	// Make sure the AOD tree is filled
	368	AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
	369	AliAODHandler* ah =
	370	dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
	371	if (!ah)
	372	AliFatal("No AOD output handler set in analysis manager");
	373
	374	ah->SetFillAOD(kTRUE);
	375	}
7e4038b5	376
	377	//____________________________________________________________________
	378	void
0bd4b00f	379	AliForwardMultiplicityTask::Print(Option_t* option) const
7e4038b5	380	{
0bd4b00f	381	std::cout << "AliForwardMultiplicityTask: " << GetName() << "\n"
	382	<< " Enable low flux code: " << (fEnableLowFlux ? "yes" : "no")
	383	<< std::endl;
	384	gROOT->IncreaseDirLevel();
	385	fEventInspector .Print(option);
	386	fEnergyFitter .Print(option);
	387	fSharingFilter .Print(option);
	388	fDensityCalculator.Print(option);
	389	fCorrections .Print(option);
	390	fHistCollector .Print(option);
	391	gROOT->DecreaseDirLevel();
7e4038b5	392	}
	393
	394	//
	395	// EOF
	396	//