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

// 
// Calculate the multiplicity in the forward regions event-by-event 
// 
// Inputs: 
//   - AliESDEvent 
//   - Kinematics
//   - Track references
//
// Outputs: 
//   - AliAODForwardMult 
// 
// Histograms 
//   
// Corrections used 
// 
#include "AliForwardMCMultiplicityTask.h"
#include "AliTriggerAnalysis.h"
#include "AliPhysicsSelection.h"
#include "AliLog.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>

//====================================================================
AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask()
  : AliForwardMultiplicityBase(),
    fHData(0),
    fESDFMD(),
    fHistos(),
    fAODFMD(),
    fMCESDFMD(),
    fMCHistos(),
    fMCAODFMD(),
    fPrimary(0),
    fEventInspector(),
    fEnergyFitter(),
    fSharingFilter(),
    fDensityCalculator(),
    fCorrections(),
    fHistCollector(),
    fList(0)
{
  // 
  // Constructor
  //
}

//____________________________________________________________________
AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const char* name)
  : AliForwardMultiplicityBase(name), 
    fHData(0),
    fESDFMD(),
    fHistos(),
    fAODFMD(kFALSE),
    fMCESDFMD(),
    fMCHistos(),
    fMCAODFMD(kTRUE),
    fPrimary(0),
    fEventInspector("event"),
    fEnergyFitter("energy"),
    fSharingFilter("sharing"), 
    fDensityCalculator("density"),
    fCorrections("corrections"),
    fHistCollector("collector"),
    fList(0)
{
  // 
  // Constructor 
  // 
  // Parameters:
  //    name Name of task 
  //
  DefineOutput(1, TList::Class());
}

//____________________________________________________________________
AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const AliForwardMCMultiplicityTask& o)
  : AliForwardMultiplicityBase(o),
    fHData(o.fHData),
    fESDFMD(o.fESDFMD),
    fHistos(o.fHistos),
    fAODFMD(o.fAODFMD),
    fMCESDFMD(o.fMCESDFMD),
    fMCHistos(o.fMCHistos),
    fMCAODFMD(o.fMCAODFMD),
    fPrimary(o.fPrimary),
    fEventInspector(o.fEventInspector),
    fEnergyFitter(o.fEnergyFitter),
    fSharingFilter(o.fSharingFilter),
    fDensityCalculator(o.fDensityCalculator),
    fCorrections(o.fCorrections),
    fHistCollector(o.fHistCollector),
    fList(o.fList) 
{
  // 
  // Copy constructor 
  // 
  // Parameters:
  //    o Object to copy from 
  //
  DefineOutput(1, TList::Class());
}

//____________________________________________________________________
AliForwardMCMultiplicityTask&
AliForwardMCMultiplicityTask::operator=(const AliForwardMCMultiplicityTask& o)
{
  // 
  // Assignment operator 
  // 
  // Parameters:
  //    o Object to assign from 
  // 
  // Return:
  //    Reference to this object 
  //
  AliForwardMultiplicityBase::operator=(o);

  fHData             = o.fHData;
  fEventInspector    = o.fEventInspector;
  fEnergyFitter      = o.fEnergyFitter;
  fSharingFilter     = o.fSharingFilter;
  fDensityCalculator = o.fDensityCalculator;
  fCorrections       = o.fCorrections;
  fHistCollector     = o.fHistCollector;
  fHistos            = o.fHistos;
  fAODFMD            = o.fAODFMD;
  fMCHistos          = o.fMCHistos;
  fMCAODFMD          = o.fMCAODFMD;
  fPrimary           = o.fPrimary;
  fList              = o.fList;

  return *this;
}

//____________________________________________________________________
void
AliForwardMCMultiplicityTask::SetDebug(Int_t dbg)
{
  // 
  // Set debug level 
  // 
  // Parameters:
  //    dbg debug level
  //
  fEventInspector.SetDebug(dbg);
  fEnergyFitter.SetDebug(dbg);
  fSharingFilter.SetDebug(dbg);
  fDensityCalculator.SetDebug(dbg);
  fCorrections.SetDebug(dbg);
  fHistCollector.SetDebug(dbg);
}
//____________________________________________________________________
void
AliForwardMCMultiplicityTask::InitializeSubs()
{
  // 
  // Initialise the sub objects and stuff.  Called on first event 
  // 
  //
  const TAxis* pe = 0;
  const TAxis* pv = 0;

  if (!ReadCorrections(pe,pv,true)) return;

  fHistos.Init(*pe);
  fAODFMD.Init(*pe);
  fMCHistos.Init(*pe);
  fMCAODFMD.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);
  fDensityCalculator.Init(*pe);
  fCorrections.Init(*pe);
  fHistCollector.Init(*pv,*pe);

  this->Print();
}

//____________________________________________________________________
void
AliForwardMCMultiplicityTask::UserCreateOutputObjects()
{
  // 
  // Create output objects 
  // 
  //
  fList = new TList;
  fList->SetOwner();

  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);

  TObject* mcobj = &fMCAODFMD;
  ah->AddBranch("AliAODForwardMult", &mcobj);

  fPrimary = new TH2D("primary", "MC Primaries", 
		      200, -4, 6, 20, 0, 2*TMath::Pi());
  fPrimary->SetXTitle("#eta");
  fPrimary->SetYTitle("#varphi [radians]");
  fPrimary->SetZTitle("d^{2}N_{ch}/d#etad#phi");
  fPrimary->Sumw2();
  fPrimary->SetStats(0);
  fPrimary->SetDirectory(0);
  ah->AddBranch("TH2D", &fPrimary);

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

  PostData(1, fList);
}
//____________________________________________________________________
void
AliForwardMCMultiplicityTask::UserExec(Option_t*)
{
  // 
  // Process each event 
  // 
  // Parameters:
  //    option Not used
  //  

  // Get the input data 
  AliESDEvent* esd     = GetESDEvent();
  AliMCEvent*  mcEvent = MCEvent();

  // Clear stuff 
  fHistos.Clear();
  fESDFMD.Clear();
  fAODFMD.Clear();
  fMCHistos.Clear();
  fMCESDFMD.Clear();
  fMCAODFMD.Clear();
  fPrimary->Reset();

  Bool_t   lowFlux  = kFALSE;
  UInt_t   triggers = 0;
  UShort_t ivz      = 0;
  Double_t vz       = 0;
  Double_t cent     = 0;
  UInt_t   found    = fEventInspector.Process(esd, triggers, lowFlux, 
					      ivz, vz, cent);
  UShort_t ivzMC    = 0;
  Double_t vzMC     = 0;
  Double_t phiR     = 0;
  Double_t b        = 0;
  Int_t    npart    = 0;
  Int_t    nbin     = 0;
  // UInt_t   foundMC  = 
  fEventInspector.ProcessMC(mcEvent, triggers, ivzMC, vzMC, b, 
			    npart, nbin, phiR);
  fEventInspector.CompareResults(vz, vzMC, cent, b, npart, nbin);
  
  //Store all events
  MarkEventForStore();
  
  Bool_t isAccepted = true;
  if (found & AliFMDEventInspector::kNoEvent)    isAccepted = false; // return;
  if (found & AliFMDEventInspector::kNoTriggers) isAccepted = false; // return;
  //MarkEventForStore();
  // Always set the B trigger - each MC event _is_ a collision 
  triggers |= AliAODForwardMult::kB;
  // Set trigger bits, and mark this event for storage 
  fAODFMD.SetTriggerBits(triggers);
  fAODFMD.SetSNN(fEventInspector.GetEnergy());
  fAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
  fAODFMD.SetCentrality(cent);

  fMCAODFMD.SetTriggerBits(triggers);
  fMCAODFMD.SetSNN(fEventInspector.GetEnergy());
  fMCAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
  fMCAODFMD.SetCentrality(cent);
  
  //All events should be stored - HHD
  //if (isAccepted) MarkEventForStore();

  // Disable this check on SPD - will bias data 
  // if (found & AliFMDEventInspector::kNoSPD)  isAccepted = false; // return;
  if (found & AliFMDEventInspector::kNoFMD)     isAccepted = false; // return;
  if (found & AliFMDEventInspector::kNoVertex)  isAccepted = false; // return;

  if (isAccepted) {
    fAODFMD.SetIpZ(vz);
    fMCAODFMD.SetIpZ(vz);
  }
  if (found & AliFMDEventInspector::kBadVertex) isAccepted = false; // 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 (isAccepted && !fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD)) { 
    AliWarning("Sharing filter failed!");
    return;
  }
  if (!fSharingFilter.FilterMC(*esdFMD, *mcEvent, vz, fMCESDFMD, fPrimary)) { 
    AliWarning("MC Sharing filter failed!");
    return;
  }
  if (!isAccepted) return; // Exit on MC event w/o trigger, vertex, data
  // HHD if (!isAccepted) return; // Exit on MC event w/o trigger, vertex, data
  
  //MarkEventForStore();
  fSharingFilter.CompareResults(fESDFMD, fMCESDFMD);

  // Do the energy stuff 
  if (!fEnergyFitter.Accumulate(*esdFMD, cent, 
				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;
  }
  if (!fDensityCalculator.CalculateMC(fMCESDFMD, fMCHistos)) { 
    AliWarning("MC Density calculator failed!");
    return;
  }
  fDensityCalculator.CompareResults(fHistos, fMCHistos);
  
  // Do the secondary and other corrections. 
  if (!fCorrections.Correct(fHistos, ivz)) { 
    AliWarning("Corrections failed");
    return;
  }
  if (!fCorrections.CorrectMC(fMCHistos, ivz)) { 
    AliWarning("MC Corrections failed");
    return;
  }
  fCorrections.CompareResults(fHistos, fMCHistos);
    
  if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) {
    AliWarning("Histogram collector failed");
    return;
  }
  if (!fHistCollector.Collect(fMCHistos, ivz, fMCAODFMD.GetHistogram())) {
    AliWarning("MC Histogram collector failed");
    return;
  }

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

  PostData(1, fList);
}

//____________________________________________________________________
void
AliForwardMCMultiplicityTask::Terminate(Option_t*)
{
  // 
  // End of job
  // 
  // Parameters:
  //    option Not used 
  //
  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;
  TH1I* hEventsTrVtx = 0;
  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,Int_t(hEventsTr->Integral()));
  fDensityCalculator.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
  fCorrections.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
}


//
// EOF
//
Commit	Line	Data
7984e5f7	1	//
	2	// Calculate the multiplicity in the forward regions event-by-event
	3	//
	4	// Inputs:
	5	// - AliESDEvent
	6	// - Kinematics
	7	// - Track references
	8	//
	9	// Outputs:
	10	// - AliAODForwardMult
	11	//
	12	// Histograms
	13	//
	14	// Corrections used
	15	//
285e7b27	16	#include "AliForwardMCMultiplicityTask.h"
	17	#include "AliTriggerAnalysis.h"
	18	#include "AliPhysicsSelection.h"
	19	#include "AliLog.h"
285e7b27	20	#include "AliESDEvent.h"
	21	#include "AliAODHandler.h"
	22	#include "AliMultiplicity.h"
	23	#include "AliInputEventHandler.h"
	24	#include "AliForwardCorrectionManager.h"
	25	#include "AliAnalysisManager.h"
	26	#include <TH1.h>
	27	#include <TDirectory.h>
	28	#include <TTree.h>
	29	#include <TROOT.h>
285e7b27	30
	31	//====================================================================
	32	AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask()
	33	: AliForwardMultiplicityBase(),
	34	fHData(0),
	35	fESDFMD(),
	36	fHistos(),
	37	fAODFMD(),
	38	fMCESDFMD(),
	39	fMCHistos(),
	40	fMCAODFMD(),
4cbdf467	41	fPrimary(0),
285e7b27	42	fEventInspector(),
	43	fEnergyFitter(),
	44	fSharingFilter(),
	45	fDensityCalculator(),
	46	fCorrections(),
	47	fHistCollector(),
	48	fList(0)
	49	{
7984e5f7	50	//
	51	// Constructor
	52	//
285e7b27	53	}
	54
	55	//____________________________________________________________________
	56	AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const char* name)
	57	: AliForwardMultiplicityBase(name),
	58	fHData(0),
	59	fESDFMD(),
	60	fHistos(),
	61	fAODFMD(kFALSE),
	62	fMCESDFMD(),
	63	fMCHistos(),
	64	fMCAODFMD(kTRUE),
4cbdf467	65	fPrimary(0),
285e7b27	66	fEventInspector("event"),
	67	fEnergyFitter("energy"),
	68	fSharingFilter("sharing"),
	69	fDensityCalculator("density"),
	70	fCorrections("corrections"),
	71	fHistCollector("collector"),
	72	fList(0)
	73	{
7984e5f7	74	//
	75	// Constructor
	76	//
	77	// Parameters:
	78	// name Name of task
	79	//
285e7b27	80	DefineOutput(1, TList::Class());
	81	}
	82
	83	//____________________________________________________________________
	84	AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const AliForwardMCMultiplicityTask& o)
	85	: AliForwardMultiplicityBase(o),
	86	fHData(o.fHData),
	87	fESDFMD(o.fESDFMD),
	88	fHistos(o.fHistos),
	89	fAODFMD(o.fAODFMD),
	90	fMCESDFMD(o.fMCESDFMD),
	91	fMCHistos(o.fMCHistos),
	92	fMCAODFMD(o.fMCAODFMD),
4cbdf467	93	fPrimary(o.fPrimary),
285e7b27	94	fEventInspector(o.fEventInspector),
285e7b27	95	fEnergyFitter(o.fEnergyFitter),
4cbdf467	96	fSharingFilter(o.fSharingFilter),
285e7b27	97	fDensityCalculator(o.fDensityCalculator),
	98	fCorrections(o.fCorrections),
	99	fHistCollector(o.fHistCollector),
	100	fList(o.fList)
	101	{
7984e5f7	102	//
	103	// Copy constructor
	104	//
	105	// Parameters:
	106	// o Object to copy from
	107	//
285e7b27	108	DefineOutput(1, TList::Class());
	109	}
	110
	111	//____________________________________________________________________
	112	AliForwardMCMultiplicityTask&
	113	AliForwardMCMultiplicityTask::operator=(const AliForwardMCMultiplicityTask& o)
	114	{
7984e5f7	115	//
	116	// Assignment operator
	117	//
	118	// Parameters:
	119	// o Object to assign from
	120	//
	121	// Return:
	122	// Reference to this object
	123	//
285e7b27	124	AliForwardMultiplicityBase::operator=(o);
	125
	126	fHData = o.fHData;
	127	fEventInspector = o.fEventInspector;
	128	fEnergyFitter = o.fEnergyFitter;
	129	fSharingFilter = o.fSharingFilter;
	130	fDensityCalculator = o.fDensityCalculator;
	131	fCorrections = o.fCorrections;
	132	fHistCollector = o.fHistCollector;
	133	fHistos = o.fHistos;
	134	fAODFMD = o.fAODFMD;
	135	fMCHistos = o.fMCHistos;
	136	fMCAODFMD = o.fMCAODFMD;
4cbdf467	137	fPrimary = o.fPrimary;
285e7b27	138	fList = o.fList;
	139
	140	return *this;
	141	}
	142
	143	//____________________________________________________________________
	144	void
	145	AliForwardMCMultiplicityTask::SetDebug(Int_t dbg)
	146	{
7984e5f7	147	//
	148	// Set debug level
	149	//
	150	// Parameters:
	151	// dbg debug level
	152	//
285e7b27	153	fEventInspector.SetDebug(dbg);
	154	fEnergyFitter.SetDebug(dbg);
	155	fSharingFilter.SetDebug(dbg);
	156	fDensityCalculator.SetDebug(dbg);
	157	fCorrections.SetDebug(dbg);
	158	fHistCollector.SetDebug(dbg);
	159	}
	160	//____________________________________________________________________
	161	void
	162	AliForwardMCMultiplicityTask::InitializeSubs()
	163	{
7984e5f7	164	//
	165	// Initialise the sub objects and stuff. Called on first event
	166	//
	167	//
7ec4d843	168	const TAxis* pe = 0;
7ec4d843	169	const TAxis* pv = 0;
1174780f	170
19abe41d	171	if (!ReadCorrections(pe,pv,true)) return;
285e7b27	172
	173	fHistos.Init(*pe);
	174	fAODFMD.Init(*pe);
	175	fMCHistos.Init(*pe);
	176	fMCAODFMD.Init(*pe);
	177
	178	fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
	179	fHData->SetStats(0);
	180	fHData->SetDirectory(0);
4cbdf467	181
285e7b27	182	fList->Add(fHData);
	183
	184	fEnergyFitter.Init(*pe);
	185	fEventInspector.Init(*pv);
	186	fDensityCalculator.Init(*pe);
	187	fCorrections.Init(*pe);
12fffad7	188	fHistCollector.Init(pv,pe);
285e7b27	189
	190	this->Print();
	191	}
	192
	193	//____________________________________________________________________
	194	void
	195	AliForwardMCMultiplicityTask::UserCreateOutputObjects()
	196	{
7984e5f7	197	//
	198	// Create output objects
	199	//
	200	//
285e7b27	201	fList = new TList;
e1f47419	202	fList->SetOwner();
285e7b27	203
	204	AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
	205	AliAODHandler* ah =
	206	dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
	207	if (!ah) AliFatal("No AOD output handler set in analysis manager");
	208
	209
	210	TObject* obj = &fAODFMD;
	211	ah->AddBranch("AliAODForwardMult", &obj);
	212
	213	TObject* mcobj = &fMCAODFMD;
	214	ah->AddBranch("AliAODForwardMult", &mcobj);
	215
4cbdf467	216	fPrimary = new TH2D("primary", "MC Primaries",
	217	200, -4, 6, 20, 0, 2*TMath::Pi());
	218	fPrimary->SetXTitle("#eta");
	219	fPrimary->SetYTitle("#varphi [radians]");
	220	fPrimary->SetZTitle("d^{2}N_{ch}/d#etad#phi");
	221	fPrimary->Sumw2();
	222	fPrimary->SetStats(0);
	223	fPrimary->SetDirectory(0);
	224	ah->AddBranch("TH2D", &fPrimary);
	225
285e7b27	226	fEventInspector.DefineOutput(fList);
	227	fEnergyFitter.DefineOutput(fList);
	228	fSharingFilter.DefineOutput(fList);
	229	fDensityCalculator.DefineOutput(fList);
	230	fCorrections.DefineOutput(fList);
12fffad7	231	fHistCollector.DefineOutput(fList);
285e7b27	232
	233	PostData(1, fList);
	234	}
	235	//____________________________________________________________________
	236	void
	237	AliForwardMCMultiplicityTask::UserExec(Option_t*)
	238	{
7984e5f7	239	//
	240	// Process each event
	241	//
	242	// Parameters:
	243	// option Not used
	244	//
	245
285e7b27	246	// Get the input data
e1f47419	247	AliESDEvent* esd = GetESDEvent();
e1f47419	248	AliMCEvent* mcEvent = MCEvent();
285e7b27	249
285e7b27	250	// Clear stuff
	251	fHistos.Clear();
	252	fESDFMD.Clear();
	253	fAODFMD.Clear();
	254	fMCHistos.Clear();
	255	fMCESDFMD.Clear();
	256	fMCAODFMD.Clear();
4cbdf467	257	fPrimary->Reset();
285e7b27	258
	259	Bool_t lowFlux = kFALSE;
	260	UInt_t triggers = 0;
	261	UShort_t ivz = 0;
	262	Double_t vz = 0;
5e4d905e	263	Double_t cent = 0;
	264	UInt_t found = fEventInspector.Process(esd, triggers, lowFlux,
	265	ivz, vz, cent);
e1f47419	266	UShort_t ivzMC = 0;
	267	Double_t vzMC = 0;
	268	Double_t phiR = 0;
	269	Double_t b = 0;
e308a636	270	Int_t npart = 0;
e308a636	271	Int_t nbin = 0;
e1f47419	272	// UInt_t foundMC =
e308a636	273	fEventInspector.ProcessMC(mcEvent, triggers, ivzMC, vzMC, b,
	274	npart, nbin, phiR);
	275	fEventInspector.CompareResults(vz, vzMC, cent, b, npart, nbin);
1dbfc345	276
e333578d	277	//Store all events
148c39de	278	MarkEventForStore();
e333578d	279
	280	Bool_t isAccepted = true;
	281	if (found & AliFMDEventInspector::kNoEvent) isAccepted = false; // return;
	282	if (found & AliFMDEventInspector::kNoTriggers) isAccepted = false; // return;
148c39de	283	//MarkEventForStore();
b30dee70	284	// Always set the B trigger - each MC event _is_ a collision
b30dee70	285	triggers \|= AliAODForwardMult::kB;
285e7b27	286	// Set trigger bits, and mark this event for storage
285e7b27	287	fAODFMD.SetTriggerBits(triggers);
e308a636	288	fAODFMD.SetSNN(fEventInspector.GetEnergy());
e308a636	289	fAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
e58000b7	290	fAODFMD.SetCentrality(cent);
e308a636	291
	292	fMCAODFMD.SetTriggerBits(triggers);
	293	fMCAODFMD.SetSNN(fEventInspector.GetEnergy());
	294	fMCAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
	295	fMCAODFMD.SetCentrality(cent);
e58000b7	296
e333578d	297	//All events should be stored - HHD
148c39de	298	//if (isAccepted) MarkEventForStore();
285e7b27	299
e1f47419	300	// Disable this check on SPD - will bias data
e1f47419	301	// if (found & AliFMDEventInspector::kNoSPD) isAccepted = false; // return;
e333578d	302	if (found & AliFMDEventInspector::kNoFMD) isAccepted = false; // return;
	303	if (found & AliFMDEventInspector::kNoVertex) isAccepted = false; // return;
	304
	305	if (isAccepted) {
	306	fAODFMD.SetIpZ(vz);
	307	fMCAODFMD.SetIpZ(vz);
	308	}
	309	if (found & AliFMDEventInspector::kBadVertex) isAccepted = false; // return;
285e7b27	310
	311	// We we do not want to use low flux specific code, we disable it here.
	312	if (!fEnableLowFlux) lowFlux = false;
148c39de	313
148c39de	314
285e7b27	315
	316	// Get FMD data
	317	AliESDFMD* esdFMD = esd->GetFMDData();
e1f47419	318
285e7b27	319	// Apply the sharing filter (or hit merging or clustering if you like)
e333578d	320	if (isAccepted && !fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD)) {
285e7b27	321	AliWarning("Sharing filter failed!");
	322	return;
	323	}
4cbdf467	324	if (!fSharingFilter.FilterMC(esdFMD, mcEvent, vz, fMCESDFMD, fPrimary)) {
285e7b27	325	AliWarning("MC Sharing filter failed!");
	326	return;
	327	}
e333578d	328	if (!isAccepted) return; // Exit on MC event w/o trigger, vertex, data
148c39de	329	// HHD if (!isAccepted) return; // Exit on MC event w/o trigger, vertex, data
1dbfc345	330
148c39de	331	//MarkEventForStore();
285e7b27	332	fSharingFilter.CompareResults(fESDFMD, fMCESDFMD);
	333
	334	// Do the energy stuff
5e4d905e	335	if (!fEnergyFitter.Accumulate(*esdFMD, cent,
5e4d905e	336	triggers & AliAODForwardMult::kEmpty)){
285e7b27	337	AliWarning("Energy fitter failed");
	338	return;
	339	}
	340
	341	// Calculate the inclusive charged particle density
	342	if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux)) {
	343	AliWarning("Density calculator failed!");
	344	return;
	345	}
	346	if (!fDensityCalculator.CalculateMC(fMCESDFMD, fMCHistos)) {
	347	AliWarning("MC Density calculator failed!");
	348	return;
	349	}
	350	fDensityCalculator.CompareResults(fHistos, fMCHistos);
	351
	352	// Do the secondary and other corrections.
	353	if (!fCorrections.Correct(fHistos, ivz)) {
	354	AliWarning("Corrections failed");
	355	return;
	356	}
	357	if (!fCorrections.CorrectMC(fMCHistos, ivz)) {
	358	AliWarning("MC Corrections failed");
	359	return;
	360	}
	361	fCorrections.CompareResults(fHistos, fMCHistos);
	362
	363	if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) {
	364	AliWarning("Histogram collector failed");
	365	return;
	366	}
	367	if (!fHistCollector.Collect(fMCHistos, ivz, fMCAODFMD.GetHistogram())) {
	368	AliWarning("MC Histogram collector failed");
	369	return;
	370	}
	371
	372	if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
	373	fHData->Add(&(fAODFMD.GetHistogram()));
	374
	375	PostData(1, fList);
	376	}
	377
	378	//____________________________________________________________________
	379	void
	380	AliForwardMCMultiplicityTask::Terminate(Option_t*)
	381	{
7984e5f7	382	//
	383	// End of job
	384	//
	385	// Parameters:
	386	// option Not used
	387	//
285e7b27	388	TList* list = dynamic_cast<TList*>(GetOutputData(1));
	389	if (!list) {
	390	AliError(Form("No output list defined (%p)", GetOutputData(1)));
	391	if (GetOutputData(1)) GetOutputData(1)->Print();
	392	return;
	393	}
	394
	395	// Get our histograms from the container
	396	TH1I* hEventsTr = 0;
	397	TH1I* hEventsTrVtx = 0;
	398	TH1I* hTriggers = 0;
	399	if (!fEventInspector.FetchHistograms(list, hEventsTr,
	400	hEventsTrVtx, hTriggers)) {
	401	AliError(Form("Didn't get histograms from event selector "
	402	"(hEventsTr=%p,hEventsTrVtx=%p)",
	403	hEventsTr, hEventsTrVtx));
	404	list->ls();
	405	return;
	406	}
	407
	408	TH2D* hData = static_cast<TH2D*>(list->FindObject("d2Ndetadphi"));
	409	if (!hData) {
	410	AliError(Form("Couldn't get our summed histogram from output "
	411	"list %s (d2Ndetadphi=%p)", list->GetName(), hData));
	412	list->ls();
	413	return;
	414	}
	415
	416	// TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
	417	TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
	418	TH1D* norm = hData->ProjectionX("norm", 0, 1, "");
	419	dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
	420	dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
	421	dNdeta->Divide(norm);
	422	dNdeta->SetStats(0);
	423	dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
	424	"width");
	425	list->Add(dNdeta);
	426	list->Add(norm);
	427
	428	fEnergyFitter.Fit(list);
6feacd76	429	fSharingFilter.ScaleHistograms(list,Int_t(hEventsTr->Integral()));
	430	fDensityCalculator.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
	431	fCorrections.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
285e7b27	432	}
285e7b27	433
285e7b27	434
	435	//
	436	// EOF
	437	//