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

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

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

  // 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	//____________________________________________________________________
	182	void
	183	AliForwardMCMultiplicityTask::InitializeSubs()
	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
19abe41d	192	if (!ReadCorrections(pe,pv,true)) return;
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
	249	this->Print();
	250	}
	251
	252	//____________________________________________________________________
	253	void
	254	AliForwardMCMultiplicityTask::UserCreateOutputObjects()
	255	{
7984e5f7	256	//
	257	// Create output objects
	258	//
	259	//
285e7b27	260	fList = new TList;
e1f47419	261	fList->SetOwner();
285e7b27	262
	263	AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
	264	AliAODHandler* ah =
	265	dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
	266	if (!ah) AliFatal("No AOD output handler set in analysis manager");
	267
	268
	269	TObject* obj = &fAODFMD;
	270	ah->AddBranch("AliAODForwardMult", &obj);
	271
	272	TObject* mcobj = &fMCAODFMD;
	273	ah->AddBranch("AliAODForwardMult", &mcobj);
	274
2b556440	275	TObject* epobj = &fAODFMD;
	276	ah->AddBranch("AliAODForwardEP", &epobj);
	277
4cbdf467	278	fPrimary = new TH2D("primary", "MC Primaries",
	279	200, -4, 6, 20, 0, 2*TMath::Pi());
	280	fPrimary->SetXTitle("#eta");
	281	fPrimary->SetYTitle("#varphi [radians]");
	282	fPrimary->SetZTitle("d^{2}N_{ch}/d#etad#phi");
	283	fPrimary->Sumw2();
	284	fPrimary->SetStats(0);
	285	fPrimary->SetDirectory(0);
	286	ah->AddBranch("TH2D", &fPrimary);
	287
285e7b27	288	fEventInspector.DefineOutput(fList);
285e7b27	289	fSharingFilter.DefineOutput(fList);
	290	fDensityCalculator.DefineOutput(fList);
	291	fCorrections.DefineOutput(fList);
12fffad7	292	fHistCollector.DefineOutput(fList);
2b556440	293	fEventPlaneFinder.DefineOutput(fList);
285e7b27	294
	295	PostData(1, fList);
	296	}
	297	//____________________________________________________________________
	298	void
	299	AliForwardMCMultiplicityTask::UserExec(Option_t*)
	300	{
7984e5f7	301	//
	302	// Process each event
	303	//
	304	// Parameters:
	305	// option Not used
	306	//
	307
f7cfc454	308	// Read production details
	309	if (fFirstEvent)
	310	fEventInspector.ReadProductionDetails(MCEvent());
	311
285e7b27	312	// Get the input data
e1f47419	313	AliESDEvent* esd = GetESDEvent();
e1f47419	314	AliMCEvent* mcEvent = MCEvent();
285e7b27	315
285e7b27	316	// Clear stuff
	317	fHistos.Clear();
	318	fESDFMD.Clear();
	319	fAODFMD.Clear();
2b556440	320	fAODEP.Clear();
285e7b27	321	fMCHistos.Clear();
	322	fMCESDFMD.Clear();
	323	fMCAODFMD.Clear();
4cbdf467	324	fPrimary->Reset();
285e7b27	325
5bb5d1f6	326	Bool_t lowFlux = kFALSE;
	327	UInt_t triggers = 0;
	328	UShort_t ivz = 0;
	329	Double_t vz = 0;
21d778b1	330	Double_t cent = -1;
5bb5d1f6	331	UShort_t nClusters = 0;
	332	UInt_t found = fEventInspector.Process(esd, triggers, lowFlux,
	333	ivz, vz, cent, nClusters);
e1f47419	334	UShort_t ivzMC = 0;
	335	Double_t vzMC = 0;
	336	Double_t phiR = 0;
	337	Double_t b = 0;
e308a636	338	Int_t npart = 0;
e308a636	339	Int_t nbin = 0;
e1f47419	340	// UInt_t foundMC =
e308a636	341	fEventInspector.ProcessMC(mcEvent, triggers, ivzMC, vzMC, b,
	342	npart, nbin, phiR);
	343	fEventInspector.CompareResults(vz, vzMC, cent, b, npart, nbin);
1dbfc345	344
e333578d	345	//Store all events
148c39de	346	MarkEventForStore();
e333578d	347
	348	Bool_t isAccepted = true;
	349	if (found & AliFMDEventInspector::kNoEvent) isAccepted = false; // return;
	350	if (found & AliFMDEventInspector::kNoTriggers) isAccepted = false; // return;
148c39de	351	//MarkEventForStore();
b30dee70	352	// Always set the B trigger - each MC event _is_ a collision
b30dee70	353	triggers \|= AliAODForwardMult::kB;
285e7b27	354	// Set trigger bits, and mark this event for storage
285e7b27	355	fAODFMD.SetTriggerBits(triggers);
e308a636	356	fAODFMD.SetSNN(fEventInspector.GetEnergy());
e308a636	357	fAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
e58000b7	358	fAODFMD.SetCentrality(cent);
5bb5d1f6	359	fAODFMD.SetNClusters(nClusters);
e308a636	360
	361	fMCAODFMD.SetTriggerBits(triggers);
	362	fMCAODFMD.SetSNN(fEventInspector.GetEnergy());
	363	fMCAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
	364	fMCAODFMD.SetCentrality(cent);
5bb5d1f6	365	fMCAODFMD.SetNClusters(nClusters);
e58000b7	366
e333578d	367	//All events should be stored - HHD
148c39de	368	//if (isAccepted) MarkEventForStore();
285e7b27	369
e1f47419	370	// Disable this check on SPD - will bias data
e1f47419	371	// if (found & AliFMDEventInspector::kNoSPD) isAccepted = false; // return;
e333578d	372	if (found & AliFMDEventInspector::kNoFMD) isAccepted = false; // return;
	373	if (found & AliFMDEventInspector::kNoVertex) isAccepted = false; // return;
	374
	375	if (isAccepted) {
	376	fAODFMD.SetIpZ(vz);
	377	fMCAODFMD.SetIpZ(vz);
	378	}
	379	if (found & AliFMDEventInspector::kBadVertex) isAccepted = false; // return;
285e7b27	380
	381	// We we do not want to use low flux specific code, we disable it here.
	382	if (!fEnableLowFlux) lowFlux = false;
148c39de	383
148c39de	384
285e7b27	385
	386	// Get FMD data
	387	AliESDFMD* esdFMD = esd->GetFMDData();
e1f47419	388
285e7b27	389	// Apply the sharing filter (or hit merging or clustering if you like)
6f4a5c0d	390	if (isAccepted && !fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, vz)) {
285e7b27	391	AliWarning("Sharing filter failed!");
	392	return;
	393	}
4cbdf467	394	if (!fSharingFilter.FilterMC(esdFMD, mcEvent, vz, fMCESDFMD, fPrimary)) {
285e7b27	395	AliWarning("MC Sharing filter failed!");
	396	return;
	397	}
e333578d	398	if (!isAccepted) return; // Exit on MC event w/o trigger, vertex, data
148c39de	399	// HHD if (!isAccepted) return; // Exit on MC event w/o trigger, vertex, data
1dbfc345	400
148c39de	401	//MarkEventForStore();
285e7b27	402	fSharingFilter.CompareResults(fESDFMD, fMCESDFMD);
285e7b27	403
285e7b27	404	// Calculate the inclusive charged particle density
2d2513f0	405	if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux, cent, vz)) {
285e7b27	406	AliWarning("Density calculator failed!");
	407	return;
	408	}
	409	if (!fDensityCalculator.CalculateMC(fMCESDFMD, fMCHistos)) {
	410	AliWarning("MC Density calculator failed!");
	411	return;
	412	}
	413	fDensityCalculator.CompareResults(fHistos, fMCHistos);
	414
2b556440	415	if (fEventInspector.GetCollisionSystem() == AliFMDEventInspector::kPbPb) {
	416	if (!fEventPlaneFinder.FindEventplane(esd, fAODEP, &(fAODFMD.GetHistogram()) , &fHistos))
	417	AliWarning("Eventplane finder failed!");
	418	}
	419
285e7b27	420	// Do the secondary and other corrections.
	421	if (!fCorrections.Correct(fHistos, ivz)) {
	422	AliWarning("Corrections failed");
	423	return;
	424	}
	425	if (!fCorrections.CorrectMC(fMCHistos, ivz)) {
	426	AliWarning("MC Corrections failed");
	427	return;
	428	}
	429	fCorrections.CompareResults(fHistos, fMCHistos);
	430
5bb5d1f6	431	if (!fHistCollector.Collect(fHistos, fRingSums,
5bb5d1f6	432	ivz, fAODFMD.GetHistogram())) {
285e7b27	433	AliWarning("Histogram collector failed");
	434	return;
	435	}
5bb5d1f6	436	if (!fHistCollector.Collect(fMCHistos, fMCRingSums,
5bb5d1f6	437	ivz, fMCAODFMD.GetHistogram())) {
285e7b27	438	AliWarning("MC Histogram collector failed");
	439	return;
	440	}
	441
	442	if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
	443	fHData->Add(&(fAODFMD.GetHistogram()));
	444
	445	PostData(1, fList);
	446	}
	447
	448	//____________________________________________________________________
	449	void
	450	AliForwardMCMultiplicityTask::Terminate(Option_t*)
	451	{
7984e5f7	452	//
	453	// End of job
	454	//
	455	// Parameters:
	456	// option Not used
	457	//
285e7b27	458	TList* list = dynamic_cast<TList*>(GetOutputData(1));
	459	if (!list) {
	460	AliError(Form("No output list defined (%p)", GetOutputData(1)));
	461	if (GetOutputData(1)) GetOutputData(1)->Print();
	462	return;
	463	}
	464
	465	// Get our histograms from the container
	466	TH1I* hEventsTr = 0;
	467	TH1I* hEventsTrVtx = 0;
	468	TH1I* hTriggers = 0;
	469	if (!fEventInspector.FetchHistograms(list, hEventsTr,
	470	hEventsTrVtx, hTriggers)) {
	471	AliError(Form("Didn't get histograms from event selector "
	472	"(hEventsTr=%p,hEventsTrVtx=%p)",
	473	hEventsTr, hEventsTrVtx));
	474	list->ls();
	475	return;
	476	}
	477
	478	TH2D* hData = static_cast<TH2D*>(list->FindObject("d2Ndetadphi"));
	479	if (!hData) {
	480	AliError(Form("Couldn't get our summed histogram from output "
	481	"list %s (d2Ndetadphi=%p)", list->GetName(), hData));
	482	list->ls();
	483	return;
	484	}
	485
	486	// TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
	487	TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
	488	TH1D* norm = hData->ProjectionX("norm", 0, 1, "");
	489	dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
	490	dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
	491	dNdeta->Divide(norm);
	492	dNdeta->SetStats(0);
	493	dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
	494	"width");
	495	list->Add(dNdeta);
	496	list->Add(norm);
	497
5bb5d1f6	498	MakeRingdNdeta(list, "ringSums", list, "ringResults");
	499	MakeRingdNdeta(list, "mcRingSums", list, "mcRingResults", 24);
	500
6feacd76	501	fSharingFilter.ScaleHistograms(list,Int_t(hEventsTr->Integral()));
	502	fDensityCalculator.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
	503	fCorrections.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
285e7b27	504	}
285e7b27	505
285e7b27	506
	507	//
	508	// EOF
	509	//