[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(*pe);
  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 = &fAODEP;
  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;
  Double_t cMC      = 0;
  Int_t    npart    = 0;
  Int_t    nbin     = 0;
  // UInt_t   foundMC  = 
  fEventInspector.ProcessMC(mcEvent, triggers, ivzMC, vzMC, b, cMC,
			    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
285e7b27	240	fEventInspector.Init(*pv);
e2ebf8c4	241	fSharingFilter.Init(*pe);
285e7b27	242	fDensityCalculator.Init(*pe);
285e7b27	243	fCorrections.Init(*pe);
12fffad7	244	fHistCollector.Init(pv,pe);
2b556440	245	fEventPlaneFinder.Init(*pe);
285e7b27	246
285e7b27	247	this->Print();
6ff251d8	248
6ff251d8	249	return true;
285e7b27	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
87ce9010	275	TObject* epobj = &fAODEP;
2b556440	276	ah->AddBranch("AliAODForwardEP", &epobj);
2b556440	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();
6ff251d8	315	if (!esd \|\| !mcEvent) return;
285e7b27	316
285e7b27	317	// Clear stuff
	318	fHistos.Clear();
	319	fESDFMD.Clear();
	320	fAODFMD.Clear();
2b556440	321	fAODEP.Clear();
285e7b27	322	fMCHistos.Clear();
	323	fMCESDFMD.Clear();
	324	fMCAODFMD.Clear();
4cbdf467	325	fPrimary->Reset();
285e7b27	326
5bb5d1f6	327	Bool_t lowFlux = kFALSE;
	328	UInt_t triggers = 0;
	329	UShort_t ivz = 0;
	330	Double_t vz = 0;
21d778b1	331	Double_t cent = -1;
5bb5d1f6	332	UShort_t nClusters = 0;
	333	UInt_t found = fEventInspector.Process(esd, triggers, lowFlux,
	334	ivz, vz, cent, nClusters);
e1f47419	335	UShort_t ivzMC = 0;
	336	Double_t vzMC = 0;
	337	Double_t phiR = 0;
	338	Double_t b = 0;
241cca4d	339	Double_t cMC = 0;
e308a636	340	Int_t npart = 0;
e308a636	341	Int_t nbin = 0;
e1f47419	342	// UInt_t foundMC =
241cca4d	343	fEventInspector.ProcessMC(mcEvent, triggers, ivzMC, vzMC, b, cMC,
e308a636	344	npart, nbin, phiR);
e308a636	345	fEventInspector.CompareResults(vz, vzMC, cent, b, npart, nbin);
1dbfc345	346
e333578d	347	//Store all events
148c39de	348	MarkEventForStore();
e333578d	349
	350	Bool_t isAccepted = true;
	351	if (found & AliFMDEventInspector::kNoEvent) isAccepted = false; // return;
	352	if (found & AliFMDEventInspector::kNoTriggers) isAccepted = false; // return;
148c39de	353	//MarkEventForStore();
b30dee70	354	// Always set the B trigger - each MC event _is_ a collision
b30dee70	355	triggers \|= AliAODForwardMult::kB;
285e7b27	356	// Set trigger bits, and mark this event for storage
285e7b27	357	fAODFMD.SetTriggerBits(triggers);
e308a636	358	fAODFMD.SetSNN(fEventInspector.GetEnergy());
e308a636	359	fAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
e58000b7	360	fAODFMD.SetCentrality(cent);
5bb5d1f6	361	fAODFMD.SetNClusters(nClusters);
e308a636	362
	363	fMCAODFMD.SetTriggerBits(triggers);
	364	fMCAODFMD.SetSNN(fEventInspector.GetEnergy());
	365	fMCAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
	366	fMCAODFMD.SetCentrality(cent);
5bb5d1f6	367	fMCAODFMD.SetNClusters(nClusters);
e58000b7	368
e333578d	369	//All events should be stored - HHD
148c39de	370	//if (isAccepted) MarkEventForStore();
285e7b27	371
e1f47419	372	// Disable this check on SPD - will bias data
e1f47419	373	// if (found & AliFMDEventInspector::kNoSPD) isAccepted = false; // return;
e333578d	374	if (found & AliFMDEventInspector::kNoFMD) isAccepted = false; // return;
	375	if (found & AliFMDEventInspector::kNoVertex) isAccepted = false; // return;
	376
	377	if (isAccepted) {
	378	fAODFMD.SetIpZ(vz);
	379	fMCAODFMD.SetIpZ(vz);
	380	}
	381	if (found & AliFMDEventInspector::kBadVertex) isAccepted = false; // return;
285e7b27	382
	383	// We we do not want to use low flux specific code, we disable it here.
	384	if (!fEnableLowFlux) lowFlux = false;
148c39de	385
148c39de	386
285e7b27	387
	388	// Get FMD data
	389	AliESDFMD* esdFMD = esd->GetFMDData();
e1f47419	390
285e7b27	391	// Apply the sharing filter (or hit merging or clustering if you like)
6f4a5c0d	392	if (isAccepted && !fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, vz)) {
285e7b27	393	AliWarning("Sharing filter failed!");
	394	return;
	395	}
4cbdf467	396	if (!fSharingFilter.FilterMC(esdFMD, mcEvent, vz, fMCESDFMD, fPrimary)) {
285e7b27	397	AliWarning("MC Sharing filter failed!");
	398	return;
	399	}
e333578d	400	if (!isAccepted) return; // Exit on MC event w/o trigger, vertex, data
148c39de	401	// HHD if (!isAccepted) return; // Exit on MC event w/o trigger, vertex, data
1dbfc345	402
148c39de	403	//MarkEventForStore();
285e7b27	404	fSharingFilter.CompareResults(fESDFMD, fMCESDFMD);
285e7b27	405
285e7b27	406	// Calculate the inclusive charged particle density
2d2513f0	407	if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux, cent, vz)) {
285e7b27	408	AliWarning("Density calculator failed!");
	409	return;
	410	}
	411	if (!fDensityCalculator.CalculateMC(fMCESDFMD, fMCHistos)) {
	412	AliWarning("MC Density calculator failed!");
	413	return;
	414	}
	415	fDensityCalculator.CompareResults(fHistos, fMCHistos);
	416
2b556440	417	if (fEventInspector.GetCollisionSystem() == AliFMDEventInspector::kPbPb) {
	418	if (!fEventPlaneFinder.FindEventplane(esd, fAODEP, &(fAODFMD.GetHistogram()) , &fHistos))
	419	AliWarning("Eventplane finder failed!");
	420	}
	421
285e7b27	422	// Do the secondary and other corrections.
	423	if (!fCorrections.Correct(fHistos, ivz)) {
	424	AliWarning("Corrections failed");
	425	return;
	426	}
	427	if (!fCorrections.CorrectMC(fMCHistos, ivz)) {
	428	AliWarning("MC Corrections failed");
	429	return;
	430	}
	431	fCorrections.CompareResults(fHistos, fMCHistos);
	432
5bb5d1f6	433	if (!fHistCollector.Collect(fHistos, fRingSums,
5bb5d1f6	434	ivz, fAODFMD.GetHistogram())) {
285e7b27	435	AliWarning("Histogram collector failed");
	436	return;
	437	}
5bb5d1f6	438	if (!fHistCollector.Collect(fMCHistos, fMCRingSums,
5bb5d1f6	439	ivz, fMCAODFMD.GetHistogram())) {
285e7b27	440	AliWarning("MC Histogram collector failed");
	441	return;
	442	}
	443
	444	if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
	445	fHData->Add(&(fAODFMD.GetHistogram()));
	446
	447	PostData(1, fList);
	448	}
	449
	450	//____________________________________________________________________
	451	void
	452	AliForwardMCMultiplicityTask::Terminate(Option_t*)
	453	{
7984e5f7	454	//
	455	// End of job
	456	//
	457	// Parameters:
	458	// option Not used
	459	//
285e7b27	460	TList* list = dynamic_cast<TList*>(GetOutputData(1));
	461	if (!list) {
	462	AliError(Form("No output list defined (%p)", GetOutputData(1)));
	463	if (GetOutputData(1)) GetOutputData(1)->Print();
	464	return;
	465	}
	466
	467	// Get our histograms from the container
	468	TH1I* hEventsTr = 0;
	469	TH1I* hEventsTrVtx = 0;
	470	TH1I* hTriggers = 0;
	471	if (!fEventInspector.FetchHistograms(list, hEventsTr,
	472	hEventsTrVtx, hTriggers)) {
	473	AliError(Form("Didn't get histograms from event selector "
	474	"(hEventsTr=%p,hEventsTrVtx=%p)",
	475	hEventsTr, hEventsTrVtx));
	476	list->ls();
	477	return;
	478	}
	479
	480	TH2D* hData = static_cast<TH2D*>(list->FindObject("d2Ndetadphi"));
	481	if (!hData) {
	482	AliError(Form("Couldn't get our summed histogram from output "
	483	"list %s (d2Ndetadphi=%p)", list->GetName(), hData));
	484	list->ls();
	485	return;
	486	}
	487
	488	// TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
	489	TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
	490	TH1D* norm = hData->ProjectionX("norm", 0, 1, "");
	491	dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
	492	dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
	493	dNdeta->Divide(norm);
	494	dNdeta->SetStats(0);
	495	dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
	496	"width");
	497	list->Add(dNdeta);
	498	list->Add(norm);
	499
5bb5d1f6	500	MakeRingdNdeta(list, "ringSums", list, "ringResults");
	501	MakeRingdNdeta(list, "mcRingSums", list, "mcRingResults", 24);
	502
6feacd76	503	fSharingFilter.ScaleHistograms(list,Int_t(hEventsTr->Integral()));
	504	fDensityCalculator.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
	505	fCorrections.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
285e7b27	506	}
285e7b27	507
285e7b27	508
	509	//
	510	// EOF
	511	//