[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());
  DefineOutput(2, 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());
  DefineOutput(2, 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::SetupForData()
{
  // 
  // 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.SetupForData(*pv);
  fSharingFilter.SetupForData(*pe);
  fDensityCalculator.SetupForData(*pe);
  fCorrections.SetupForData(*pe);
  fHistCollector.SetupForData(*pv,*pe);
  fEventPlaneFinder.SetupForData(*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.CreateOutputObjects(fList);
  fSharingFilter.CreateOutputObjects(fList);
  fDensityCalculator.CreateOutputObjects(fList);
  fCorrections.CreateOutputObjects(fList);
  fHistCollector.CreateOutputObjects(fList);
  fEventPlaneFinder.CreateOutputObjects(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;
  TVector3 ip(1024, 1024, 0);
  Double_t cent      = -1;
  UShort_t nClusters = 0;
  UInt_t   found     = fEventInspector.Process(esd, triggers, lowFlux, 
					       ivz, ip, 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(ip.Z(), 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(ip.Z());
    fMCAODFMD.SetIpZ(ip.Z());
  }
  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,ip.Z())) { 
    AliWarning("Sharing filter failed!");
    return;
  }
  if (!fSharingFilter.FilterMC(*esdFMD, *mcEvent, ip.Z(),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, lowFlux, cent, ip)) { 
    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;
  }

  // Output list 
  TList* output = new TList;
  output->SetOwner();
  output->SetName(Form("%sResults", GetName()));

  Double_t nTr = 0, nTrVtx = 0, nAcc = 0;
  MakeSimpledNdeta(list, output, nTr, nTrVtx, nAcc);
  MakeRingdNdeta(list, "ringSums", output, "ringResults");
  MakeRingdNdeta(list, "mcRingSums", output, "mcRingResults", 24);

  fSharingFilter.Terminate(list,output,Int_t(nTr));
  fDensityCalculator.Terminate(list,output,Int_t(nTrVtx));
  fCorrections.Terminate(list,output,Int_t(nTrVtx));

  PostData(2, output);
}


//
// 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());
5934a3e3	87	DefineOutput(2, TList::Class());
285e7b27	88	}
	89
	90	//____________________________________________________________________
	91	AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const AliForwardMCMultiplicityTask& o)
	92	: AliForwardMultiplicityBase(o),
	93	fHData(o.fHData),
	94	fESDFMD(o.fESDFMD),
	95	fHistos(o.fHistos),
	96	fAODFMD(o.fAODFMD),
2b556440	97	fAODEP(o.fAODEP),
285e7b27	98	fMCESDFMD(o.fMCESDFMD),
	99	fMCHistos(o.fMCHistos),
	100	fMCAODFMD(o.fMCAODFMD),
5bb5d1f6	101	fRingSums(o.fRingSums),
5bb5d1f6	102	fMCRingSums(o.fMCRingSums),
4cbdf467	103	fPrimary(o.fPrimary),
285e7b27	104	fEventInspector(o.fEventInspector),
4cbdf467	105	fSharingFilter(o.fSharingFilter),
285e7b27	106	fDensityCalculator(o.fDensityCalculator),
	107	fCorrections(o.fCorrections),
	108	fHistCollector(o.fHistCollector),
2b556440	109	fEventPlaneFinder(o.fEventPlaneFinder),
285e7b27	110	fList(o.fList)
285e7b27	111	{
7984e5f7	112	//
	113	// Copy constructor
	114	//
	115	// Parameters:
	116	// o Object to copy from
	117	//
285e7b27	118	DefineOutput(1, TList::Class());
5934a3e3	119	DefineOutput(2, TList::Class());
285e7b27	120	}
	121
	122	//____________________________________________________________________
	123	AliForwardMCMultiplicityTask&
	124	AliForwardMCMultiplicityTask::operator=(const AliForwardMCMultiplicityTask& o)
	125	{
7984e5f7	126	//
	127	// Assignment operator
	128	//
	129	// Parameters:
	130	// o Object to assign from
	131	//
	132	// Return:
	133	// Reference to this object
	134	//
d015ecfe	135	if (&o == this) return *this;
285e7b27	136	AliForwardMultiplicityBase::operator=(o);
	137
	138	fHData = o.fHData;
	139	fEventInspector = o.fEventInspector;
285e7b27	140	fSharingFilter = o.fSharingFilter;
	141	fDensityCalculator = o.fDensityCalculator;
	142	fCorrections = o.fCorrections;
	143	fHistCollector = o.fHistCollector;
2b556440	144	fEventPlaneFinder = o.fEventPlaneFinder;
285e7b27	145	fHistos = o.fHistos;
285e7b27	146	fAODFMD = o.fAODFMD;
2b556440	147	fAODEP = o.fAODEP;
285e7b27	148	fMCHistos = o.fMCHistos;
285e7b27	149	fMCAODFMD = o.fMCAODFMD;
5bb5d1f6	150	fRingSums = o.fRingSums;
5bb5d1f6	151	fMCRingSums = o.fMCRingSums;
4cbdf467	152	fPrimary = o.fPrimary;
285e7b27	153	fList = o.fList;
	154
	155	return *this;
	156	}
	157
	158	//____________________________________________________________________
	159	void
	160	AliForwardMCMultiplicityTask::SetDebug(Int_t dbg)
	161	{
7984e5f7	162	//
	163	// Set debug level
	164	//
	165	// Parameters:
	166	// dbg debug level
	167	//
285e7b27	168	fEventInspector.SetDebug(dbg);
285e7b27	169	fSharingFilter.SetDebug(dbg);
	170	fDensityCalculator.SetDebug(dbg);
	171	fCorrections.SetDebug(dbg);
	172	fHistCollector.SetDebug(dbg);
2b556440	173	fEventPlaneFinder.SetDebug(dbg);
285e7b27	174	}
5bb5d1f6	175	//____________________________________________________________________
	176	void
	177	AliForwardMCMultiplicityTask::SetOnlyPrimary(Bool_t use)
	178	{
42856690	179	fSharingFilter.GetTrackDensity().SetUseOnlyPrimary(use);
5bb5d1f6	180	fCorrections.SetSecondaryForMC(!use);
	181	}
	182
285e7b27	183	//____________________________________________________________________
6ff251d8	184	Bool_t
5934a3e3	185	AliForwardMCMultiplicityTask::SetupForData()
285e7b27	186	{
7984e5f7	187	//
	188	// Initialise the sub objects and stuff. Called on first event
	189	//
	190	//
7ec4d843	191	const TAxis* pe = 0;
7ec4d843	192	const TAxis* pv = 0;
1174780f	193
6ff251d8	194	if (!ReadCorrections(pe,pv,true)) return false;
285e7b27	195
	196	fHistos.Init(*pe);
	197	fAODFMD.Init(*pe);
2b556440	198	fAODEP.Init(*pe);
285e7b27	199	fMCHistos.Init(*pe);
285e7b27	200	fMCAODFMD.Init(*pe);
5bb5d1f6	201	fRingSums.Init(*pe);
5bb5d1f6	202	fMCRingSums.Init(*pe);
285e7b27	203
	204	fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
	205	fHData->SetStats(0);
	206	fHData->SetDirectory(0);
4cbdf467	207
285e7b27	208	fList->Add(fHData);
285e7b27	209
5bb5d1f6	210	TList* rings = new TList;
	211	rings->SetName("ringSums");
	212	rings->SetOwner();
	213	fList->Add(rings);
	214
	215	rings->Add(fRingSums.Get(1, 'I'));
	216	rings->Add(fRingSums.Get(2, 'I'));
	217	rings->Add(fRingSums.Get(2, 'O'));
	218	rings->Add(fRingSums.Get(3, 'I'));
	219	rings->Add(fRingSums.Get(3, 'O'));
	220	fRingSums.Get(1, 'I')->SetMarkerColor(AliForwardUtil::RingColor(1, 'I'));
	221	fRingSums.Get(2, 'I')->SetMarkerColor(AliForwardUtil::RingColor(2, 'I'));
	222	fRingSums.Get(2, 'O')->SetMarkerColor(AliForwardUtil::RingColor(2, 'O'));
	223	fRingSums.Get(3, 'I')->SetMarkerColor(AliForwardUtil::RingColor(3, 'I'));
	224	fRingSums.Get(3, 'O')->SetMarkerColor(AliForwardUtil::RingColor(3, 'O'));
	225
	226	TList* mcRings = new TList;
	227	mcRings->SetName("mcRingSums");
	228	mcRings->SetOwner();
	229	fList->Add(mcRings);
	230
	231	mcRings->Add(fMCRingSums.Get(1, 'I'));
	232	mcRings->Add(fMCRingSums.Get(2, 'I'));
	233	mcRings->Add(fMCRingSums.Get(2, 'O'));
	234	mcRings->Add(fMCRingSums.Get(3, 'I'));
	235	mcRings->Add(fMCRingSums.Get(3, 'O'));
	236	fMCRingSums.Get(1, 'I')->SetMarkerColor(AliForwardUtil::RingColor(1, 'I'));
	237	fMCRingSums.Get(2, 'I')->SetMarkerColor(AliForwardUtil::RingColor(2, 'I'));
	238	fMCRingSums.Get(2, 'O')->SetMarkerColor(AliForwardUtil::RingColor(2, 'O'));
	239	fMCRingSums.Get(3, 'I')->SetMarkerColor(AliForwardUtil::RingColor(3, 'I'));
	240	fMCRingSums.Get(3, 'O')->SetMarkerColor(AliForwardUtil::RingColor(3, 'O'));
	241
5934a3e3	242	fEventInspector.SetupForData(*pv);
	243	fSharingFilter.SetupForData(*pe);
	244	fDensityCalculator.SetupForData(*pe);
	245	fCorrections.SetupForData(*pe);
	246	fHistCollector.SetupForData(pv,pe);
	247	fEventPlaneFinder.SetupForData(*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
87ce9010	277	TObject* epobj = &fAODEP;
2b556440	278	ah->AddBranch("AliAODForwardEP", &epobj);
2b556440	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
5934a3e3	290	fEventInspector.CreateOutputObjects(fList);
	291	fSharingFilter.CreateOutputObjects(fList);
	292	fDensityCalculator.CreateOutputObjects(fList);
	293	fCorrections.CreateOutputObjects(fList);
	294	fHistCollector.CreateOutputObjects(fList);
	295	fEventPlaneFinder.CreateOutputObjects(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;
5ca83fee	332	TVector3 ip(1024, 1024, 0);
21d778b1	333	Double_t cent = -1;
5bb5d1f6	334	UShort_t nClusters = 0;
5bb5d1f6	335	UInt_t found = fEventInspector.Process(esd, triggers, lowFlux,
5ca83fee	336	ivz, ip, cent, nClusters);
e1f47419	337	UShort_t ivzMC = 0;
	338	Double_t vzMC = 0;
	339	Double_t phiR = 0;
	340	Double_t b = 0;
241cca4d	341	Double_t cMC = 0;
e308a636	342	Int_t npart = 0;
e308a636	343	Int_t nbin = 0;
e1f47419	344	// UInt_t foundMC =
241cca4d	345	fEventInspector.ProcessMC(mcEvent, triggers, ivzMC, vzMC, b, cMC,
e308a636	346	npart, nbin, phiR);
5ca83fee	347	fEventInspector.CompareResults(ip.Z(), vzMC, cent, b, npart, nbin);
1dbfc345	348
e333578d	349	//Store all events
148c39de	350	MarkEventForStore();
e333578d	351
	352	Bool_t isAccepted = true;
	353	if (found & AliFMDEventInspector::kNoEvent) isAccepted = false; // return;
	354	if (found & AliFMDEventInspector::kNoTriggers) isAccepted = false; // return;
148c39de	355	//MarkEventForStore();
b30dee70	356	// Always set the B trigger - each MC event _is_ a collision
b30dee70	357	triggers \|= AliAODForwardMult::kB;
285e7b27	358	// Set trigger bits, and mark this event for storage
285e7b27	359	fAODFMD.SetTriggerBits(triggers);
e308a636	360	fAODFMD.SetSNN(fEventInspector.GetEnergy());
e308a636	361	fAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
e58000b7	362	fAODFMD.SetCentrality(cent);
5bb5d1f6	363	fAODFMD.SetNClusters(nClusters);
e308a636	364
	365	fMCAODFMD.SetTriggerBits(triggers);
	366	fMCAODFMD.SetSNN(fEventInspector.GetEnergy());
	367	fMCAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
	368	fMCAODFMD.SetCentrality(cent);
5bb5d1f6	369	fMCAODFMD.SetNClusters(nClusters);
e58000b7	370
e333578d	371	//All events should be stored - HHD
148c39de	372	//if (isAccepted) MarkEventForStore();
285e7b27	373
e1f47419	374	// Disable this check on SPD - will bias data
e1f47419	375	// if (found & AliFMDEventInspector::kNoSPD) isAccepted = false; // return;
e333578d	376	if (found & AliFMDEventInspector::kNoFMD) isAccepted = false; // return;
	377	if (found & AliFMDEventInspector::kNoVertex) isAccepted = false; // return;
	378
	379	if (isAccepted) {
5ca83fee	380	fAODFMD.SetIpZ(ip.Z());
5ca83fee	381	fMCAODFMD.SetIpZ(ip.Z());
e333578d	382	}
e333578d	383	if (found & AliFMDEventInspector::kBadVertex) isAccepted = false; // return;
285e7b27	384
	385	// We we do not want to use low flux specific code, we disable it here.
	386	if (!fEnableLowFlux) lowFlux = false;
148c39de	387
148c39de	388
285e7b27	389
	390	// Get FMD data
	391	AliESDFMD* esdFMD = esd->GetFMDData();
e1f47419	392
285e7b27	393	// Apply the sharing filter (or hit merging or clustering if you like)
5ca83fee	394	if (isAccepted && !fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD,ip.Z())) {
285e7b27	395	AliWarning("Sharing filter failed!");
	396	return;
	397	}
5ca83fee	398	if (!fSharingFilter.FilterMC(esdFMD, mcEvent, ip.Z(),fMCESDFMD,fPrimary)) {
285e7b27	399	AliWarning("MC Sharing filter failed!");
	400	return;
	401	}
e333578d	402	if (!isAccepted) return; // Exit on MC event w/o trigger, vertex, data
148c39de	403	// HHD if (!isAccepted) return; // Exit on MC event w/o trigger, vertex, data
1dbfc345	404
148c39de	405	//MarkEventForStore();
285e7b27	406	fSharingFilter.CompareResults(fESDFMD, fMCESDFMD);
285e7b27	407
285e7b27	408	// Calculate the inclusive charged particle density
5934a3e3	409	if (!fDensityCalculator.Calculate(fESDFMD, fHistos, lowFlux, cent, ip)) {
285e7b27	410	AliWarning("Density calculator failed!");
	411	return;
	412	}
	413	if (!fDensityCalculator.CalculateMC(fMCESDFMD, fMCHistos)) {
	414	AliWarning("MC Density calculator failed!");
	415	return;
	416	}
	417	fDensityCalculator.CompareResults(fHistos, fMCHistos);
	418
2b556440	419	if (fEventInspector.GetCollisionSystem() == AliFMDEventInspector::kPbPb) {
	420	if (!fEventPlaneFinder.FindEventplane(esd, fAODEP, &(fAODFMD.GetHistogram()) , &fHistos))
	421	AliWarning("Eventplane finder failed!");
	422	}
	423
285e7b27	424	// Do the secondary and other corrections.
	425	if (!fCorrections.Correct(fHistos, ivz)) {
	426	AliWarning("Corrections failed");
	427	return;
	428	}
	429	if (!fCorrections.CorrectMC(fMCHistos, ivz)) {
	430	AliWarning("MC Corrections failed");
	431	return;
	432	}
	433	fCorrections.CompareResults(fHistos, fMCHistos);
	434
5bb5d1f6	435	if (!fHistCollector.Collect(fHistos, fRingSums,
5bb5d1f6	436	ivz, fAODFMD.GetHistogram())) {
285e7b27	437	AliWarning("Histogram collector failed");
	438	return;
	439	}
5bb5d1f6	440	if (!fHistCollector.Collect(fMCHistos, fMCRingSums,
5bb5d1f6	441	ivz, fMCAODFMD.GetHistogram())) {
285e7b27	442	AliWarning("MC Histogram collector failed");
	443	return;
	444	}
	445
	446	if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
	447	fHData->Add(&(fAODFMD.GetHistogram()));
	448
	449	PostData(1, fList);
	450	}
	451
	452	//____________________________________________________________________
	453	void
	454	AliForwardMCMultiplicityTask::Terminate(Option_t*)
	455	{
7984e5f7	456	//
	457	// End of job
	458	//
	459	// Parameters:
	460	// option Not used
	461	//
285e7b27	462	TList* list = dynamic_cast<TList*>(GetOutputData(1));
	463	if (!list) {
	464	AliError(Form("No output list defined (%p)", GetOutputData(1)));
	465	if (GetOutputData(1)) GetOutputData(1)->Print();
	466	return;
	467	}
285e7b27	468
5934a3e3	469	// Output list
	470	TList* output = new TList;
	471	output->SetOwner();
	472	output->SetName(Form("%sResults", GetName()));
	473
5ca83fee	474	Double_t nTr = 0, nTrVtx = 0, nAcc = 0;
5934a3e3	475	MakeSimpledNdeta(list, output, nTr, nTrVtx, nAcc);
	476	MakeRingdNdeta(list, "ringSums", output, "ringResults");
	477	MakeRingdNdeta(list, "mcRingSums", output, "mcRingResults", 24);
	478
	479	fSharingFilter.Terminate(list,output,Int_t(nTr));
	480	fDensityCalculator.Terminate(list,output,Int_t(nTrVtx));
	481	fCorrections.Terminate(list,output,Int_t(nTrVtx));
5bb5d1f6	482
5934a3e3	483	PostData(2, output);
285e7b27	484	}
285e7b27	485
285e7b27	486
	487	//
	488	// EOF
	489	//