[u/mrichter/AliRoot.git] / PWGLF / 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(),
    fAODEP(),
    fMCESDFMD(),
    fMCHistos(),
    fMCAODFMD(),
    fRingSums(),
    fMCRingSums(),
    fPrimary(0),
    fEventInspector(),
    fSharingFilter(),
    fDensityCalculator(),
    fCorrections(),
    fHistCollector(),
    fEventPlaneFinder(),
    fList(0)
{
  // 
  // Constructor
  //
}

//____________________________________________________________________
AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const char* name)
  : AliForwardMultiplicityBase(name), 
    fHData(0),
    fESDFMD(),
    fHistos(),
    fAODFMD(kFALSE),
    fAODEP(kFALSE),
    fMCESDFMD(),
    fMCHistos(),
    fMCAODFMD(kTRUE),
    fRingSums(),
    fMCRingSums(),
    fPrimary(0),
    fEventInspector("event"),
    fSharingFilter("sharing"), 
    fDensityCalculator("density"),
    fCorrections("corrections"),
    fHistCollector("collector"),
    fEventPlaneFinder("eventplane"),
    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),
    fAODEP(o.fAODEP),
    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),
    fEventPlaneFinder(o.fEventPlaneFinder),
    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 
  //
  if (&o == this) return *this;
  AliForwardMultiplicityBase::operator=(o);

  fHData             = o.fHData;
  fEventInspector    = o.fEventInspector;
  fSharingFilter     = o.fSharingFilter;
  fDensityCalculator = o.fDensityCalculator;
  fCorrections       = o.fCorrections;
  fHistCollector     = o.fHistCollector;
  fEventPlaneFinder  = o.fEventPlaneFinder;
  fHistos            = o.fHistos;
  fAODFMD            = o.fAODFMD;
  fAODEP             = o.fAODEP;
  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);
  fEventPlaneFinder.SetDebug(dbg);
}
//____________________________________________________________________
void
AliForwardMCMultiplicityTask::SetOnlyPrimary(Bool_t use)
{
  fSharingFilter.GetTrackDensity().SetUseOnlyPrimary(use);
  fCorrections.SetSecondaryForMC(!use);
}

//____________________________________________________________________
Bool_t
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 false;

  fHistos.Init(*pe);
  fAODFMD.Init(*pe);
  fAODEP.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);
  fEventPlaneFinder.Init(*pe);

  this->Print();

  return true;
}

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

  TObject* epobj = &fAODFMD;
  ah->AddBranch("AliAODForwardEP", &epobj);
  
  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);
  fEventPlaneFinder.DefineOutput(fList);

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

  // Read production details 
  if (fFirstEvent) 
    fEventInspector.ReadProductionDetails(MCEvent());
    
  // Get the input data 
  AliESDEvent* esd     = GetESDEvent();
  AliMCEvent*  mcEvent = MCEvent();
  if (!esd || !mcEvent) return;

  // Clear stuff 
  fHistos.Clear();
  fESDFMD.Clear();
  fAODFMD.Clear();
  fAODEP.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      = -1;
  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, vz)) { 
    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, cent, vz)) { 
    AliWarning("Density calculator failed!");
    return;
  }
  if (!fDensityCalculator.CalculateMC(fMCESDFMD, fMCHistos)) { 
    AliWarning("MC Density calculator failed!");
    return;
  }
  fDensityCalculator.CompareResults(fHistos, fMCHistos);
  
  if (fEventInspector.GetCollisionSystem() == AliFMDEventInspector::kPbPb) {
    if (!fEventPlaneFinder.FindEventplane(esd, fAODEP, &(fAODFMD.GetHistogram()) , &fHistos))
      AliWarning("Eventplane finder failed!");
  }

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