[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(),
    fRingSums(),
    fMCRingSums(),
    fPrimary(0),
    fEventInspector(),
    fSharingFilter(),
    fDensityCalculator(),
    fCorrections(),
    fHistCollector(),
    fList(0)
{
  // 
  // Constructor
  //
}

//____________________________________________________________________
AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const char* name)
  : AliForwardMultiplicityBase(name), 
    fHData(0),
    fESDFMD(),
    fHistos(),
    fAODFMD(kFALSE),
    fMCESDFMD(),
    fMCHistos(),
    fMCAODFMD(kTRUE),
    fRingSums(),
    fMCRingSums(),
    fPrimary(0),
    fEventInspector("event"),
    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),
    fRingSums(o.fRingSums),
    fMCRingSums(o.fMCRingSums),
    fPrimary(o.fPrimary),
    fEventInspector(o.fEventInspector),
    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;
  fSharingFilter     = o.fSharingFilter;
  fDensityCalculator = o.fDensityCalculator;
  fCorrections       = o.fCorrections;
  fHistCollector     = o.fHistCollector;
  fHistos            = o.fHistos;
  fAODFMD            = o.fAODFMD;
  fMCHistos          = o.fMCHistos;
  fMCAODFMD          = o.fMCAODFMD;
  fRingSums          = o.fRingSums;
  fMCRingSums        = o.fMCRingSums;
  fPrimary           = o.fPrimary;
  fList              = o.fList;

  return *this;
}

//____________________________________________________________________
void
AliForwardMCMultiplicityTask::SetDebug(Int_t dbg)
{
  // 
  // Set debug level 
  // 
  // Parameters:
  //    dbg debug level
  //
  fEventInspector.SetDebug(dbg);
  fSharingFilter.SetDebug(dbg);
  fDensityCalculator.SetDebug(dbg);
  fCorrections.SetDebug(dbg);
  fHistCollector.SetDebug(dbg);
}
//____________________________________________________________________
void
AliForwardMCMultiplicityTask::SetOnlyPrimary(Bool_t use)
{
  fSharingFilter.GetTrackDensity().SetUseOnlyPrimary(use);
  fCorrections.SetSecondaryForMC(!use);
}

//____________________________________________________________________
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);
  fRingSums.Init(*pe);
  fMCRingSums.Init(*pe);

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

  TList* mcRings = new TList;
  mcRings->SetName("mcRingSums");
  mcRings->SetOwner();
  fList->Add(mcRings);

  mcRings->Add(fMCRingSums.Get(1, 'I'));
  mcRings->Add(fMCRingSums.Get(2, 'I'));
  mcRings->Add(fMCRingSums.Get(2, 'O'));
  mcRings->Add(fMCRingSums.Get(3, 'I'));
  mcRings->Add(fMCRingSums.Get(3, 'O'));
  fMCRingSums.Get(1, 'I')->SetMarkerColor(AliForwardUtil::RingColor(1, 'I'));
  fMCRingSums.Get(2, 'I')->SetMarkerColor(AliForwardUtil::RingColor(2, 'I'));
  fMCRingSums.Get(2, 'O')->SetMarkerColor(AliForwardUtil::RingColor(2, 'O'));
  fMCRingSums.Get(3, 'I')->SetMarkerColor(AliForwardUtil::RingColor(3, 'I'));
  fMCRingSums.Get(3, 'O')->SetMarkerColor(AliForwardUtil::RingColor(3, 'O'));


  fEventInspector.Init(*pv);
  fSharingFilter.Init();
  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);
  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;
  UShort_t nClusters = 0;
  UInt_t   found     = fEventInspector.Process(esd, triggers, lowFlux, 
					       ivz, vz, cent, nClusters);
  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);
  fAODFMD.SetNClusters(nClusters);

  fMCAODFMD.SetTriggerBits(triggers);
  fMCAODFMD.SetSNN(fEventInspector.GetEnergy());
  fMCAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
  fMCAODFMD.SetCentrality(cent);
  fMCAODFMD.SetNClusters(nClusters);
  
  //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);

  // 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, fRingSums, 
			      ivz, fAODFMD.GetHistogram())) {
    AliWarning("Histogram collector failed");
    return;
  }
  if (!fHistCollector.Collect(fMCHistos, fMCRingSums, 
			      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);

  MakeRingdNdeta(list, "ringSums", list, "ringResults");
  MakeRingdNdeta(list, "mcRingSums", list, "mcRingResults", 24);

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