[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliForwardMultiplicity.cxx

#include "AliForwardMultiplicity.h"
#include "AliTriggerAnalysis.h"
#include "AliPhysicsSelection.h"
#include "AliLog.h"
#include "AliFMDAnaParameters.h"
#include "AliESDEvent.h"
#include "AliAODHandler.h"
#include "AliMultiplicity.h"
#include "AliInputEventHandler.h"
#include <TH1.h>
#include <TDirectory.h>
#include <TTree.h>

//====================================================================
AliForwardMultiplicity::AliForwardMultiplicity()
  : AliAnalysisTaskSE(),
    fHEventsTr(0), 
    fHEventsTrVtx(0),
    fHTriggers(0),
    fHData(0),
    fFirstEvent(true),
    fLowFluxCut(1000),
    fESDFMD(),
    fHistos(),
    fAODFMD(),
    fSharingFilter(),
    fDensityCalculator(),
    fCorrections(),
    fHistCollector(),
    fList(0), 
    fTree(0)
{
}

//____________________________________________________________________
AliForwardMultiplicity::AliForwardMultiplicity(const char* name)
  : AliAnalysisTaskSE(name), 
    fHEventsTr(0), 
    fHEventsTrVtx(0), 
    fHTriggers(0),
    fHData(0),
    fFirstEvent(true),
    fLowFluxCut(1000),
    fESDFMD(),
    fHistos(),
    fAODFMD(kTRUE),
    fSharingFilter("sharing"), 
    fDensityCalculator("density"),
    fCorrections("corrections"),
    fHistCollector("collector"),
    fList(0), 
    fTree(0)
{
  DefineOutput(1, TList::Class());
  // DefineOutput(2, TTree::Class());
}

//____________________________________________________________________
AliForwardMultiplicity::AliForwardMultiplicity(const AliForwardMultiplicity& o)
  : AliAnalysisTaskSE(o),
    fHEventsTr(o.fHEventsTr), 
    fHEventsTrVtx(o.fHEventsTrVtx), 
    fHTriggers(o.fHTriggers),
    fHData(o.fHData),
    fFirstEvent(true),
    fLowFluxCut(1000),
    fESDFMD(o.fESDFMD),
    fHistos(o.fHistos),
    fAODFMD(o.fAODFMD),
    fSharingFilter(o.fSharingFilter),
    fDensityCalculator(o.fDensityCalculator),
    fCorrections(o.fCorrections),
    fHistCollector(o.fHistCollector),
    fList(o.fList), 
    fTree(o.fTree)
{
}

//____________________________________________________________________
AliForwardMultiplicity&
AliForwardMultiplicity::operator=(const AliForwardMultiplicity& o)
{
  fHEventsTr         = o.fHEventsTr;
  fHEventsTrVtx      = o.fHEventsTrVtx;
  fHTriggers         = o.fHTriggers;
  fHData             = o.fHData;
  fFirstEvent        = o.fFirstEvent;
  fSharingFilter     = o.fSharingFilter;
  fDensityCalculator = o.fDensityCalculator;
  fCorrections       = o.fCorrections;
  fHistCollector     = o.fHistCollector;
  fHistos            = o.fHistos;
  fAODFMD            = o.fAODFMD;
  fList              = o.fList;
  fTree              = o.fTree;

  return *this;
}

//____________________________________________________________________
void
AliForwardMultiplicity::Init()
{
  fFirstEvent = true;
}

//____________________________________________________________________
void
AliForwardMultiplicity::InitializeSubs()
{
  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  pars->Init(kTRUE);

  fHEventsTr = new TH1I("nEvents", "Number of events w/trigger", 
		      pars->GetNvtxBins(), 
		      -pars->GetVtxCutZ(), 
		      pars->GetVtxCutZ());
  fHEventsTr->SetXTitle("v_{z} [cm]");
  fHEventsTr->SetYTitle("# of events");
  fHEventsTr->SetFillColor(kRed+1);
  fHEventsTr->SetFillStyle(3001);
  fHEventsTr->SetDirectory(0);
  // fHEventsTr->Sumw2();

  fHEventsTrVtx = new TH1I("nEventsTrVtx", 
			   "Number of events w/trigger and vertex", 
			   pars->GetNvtxBins(), 
			   -pars->GetVtxCutZ(), 
			   pars->GetVtxCutZ());
  fHEventsTrVtx->SetXTitle("v_{z} [cm]");
  fHEventsTrVtx->SetYTitle("# of events");
  fHEventsTrVtx->SetFillColor(kBlue+1);
  fHEventsTrVtx->SetFillStyle(3001);
  fHEventsTrVtx->SetDirectory(0);
  // fHEventsTrVtx->Sumw2();

      
  fHTriggers = new TH1I("triggers", "Triggers", 10, 0, 10);
  fHTriggers->SetFillColor(kRed+1);
  fHTriggers->SetFillStyle(3001);
  fHTriggers->SetStats(0);
  fHTriggers->SetDirectory(0);
  fHTriggers->GetXaxis()->SetBinLabel(1,"INEL");
  fHTriggers->GetXaxis()->SetBinLabel(2,"INEL>0");
  fHTriggers->GetXaxis()->SetBinLabel(3,"NSD");
  fHTriggers->GetXaxis()->SetBinLabel(4,"Empty");
  fHTriggers->GetXaxis()->SetBinLabel(5,"A");
  fHTriggers->GetXaxis()->SetBinLabel(6,"B");
  fHTriggers->GetXaxis()->SetBinLabel(7,"C");
  fHTriggers->GetXaxis()->SetBinLabel(8,"E");

  TAxis e(pars->GetNetaBins(), pars->GetEtaMin(), pars->GetEtaMax());
  fHistos.Init(e);
  fAODFMD.Init(e);

  fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
  fHData->SetStats(0);
  fHData->SetDirectory(0);
  fSharingFilter.Init();
  fHistCollector.Init(*(fHEventsTr->GetXaxis()));
}

//____________________________________________________________________
void
AliForwardMultiplicity::UserCreateOutputObjects()
{
  fList = new TList;

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

  // fTree = new TTree("T", "T");
  // fTree->Branch("forward", &fAODFMD);

  PostData(1, fList);
  // PostData(2, fTree);
}
//____________________________________________________________________
void
AliForwardMultiplicity::UserExec(Option_t*)
{
  // Get the input data 
  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
  if (!esd) { 
    AliWarning("No ESD event found for input event");
    return;
  }

#if 0
  static Int_t nEvents = 0;
  nEvents++;
  if (nEvents % 100 == 0) AliInfo(Form("Event # %6d", nEvents));
#endif

  // On the first event, initialize the parameters 
  if (fFirstEvent) { 
    AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
    pars->SetParametersFromESD(esd);
    pars->PrintStatus();
    fFirstEvent = false;

    InitializeSubs();
  }
  // Clear stuff 
  fHistos.Clear();
  fESDFMD.Clear();
  fAODFMD.Clear();

  // Read trigger information from the ESD and store in AOD object
  if (!ReadTriggers(esd)) { 
#ifdef VERBOSE
    AliWarning("Failed to read triggers from ESD");
#endif
    return;
  }

  // Mark this event for storage 
  MarkEventForStore();

  // Check if this is a high-flux event 
  const AliMultiplicity* testmult = esd->GetMultiplicity();
  if (!testmult) { 
#ifdef VERBOSE
    AliWarning("No central multiplicity object found");
#endif
    return;
  }
  Bool_t lowFlux = testmult->GetNumberOfTracklets() < fLowFluxCut;

  // Get the FMD ESD data 
  AliESDFMD* esdFMD = esd->GetFMDData();
  if (!esdFMD) { 
#ifdef VERBOSE
    AliWarning("No FMD data found in ESD");
#endif
    return;
  }

  // Get the vertex information 
  Double_t vz   = 0;
  Bool_t   vzOk = ReadVertex(esd, vz);

  fHEventsTr->Fill(vz);
  if (!vzOk) { 
#ifdef VERBOSE
    AliWarning("Failed to read vertex from ESD");
#endif
    return;
  }
  fHEventsTrVtx->Fill(vz);

  // Get the vertex bin 
  Int_t ivz = fHEventsTr->GetXaxis()->FindBin(vz)-1;
  fAODFMD.SetIpZ(vz);
  if (ivz < 0 || ivz >= fHEventsTr->GetXaxis()->GetNbins()) { 
#if 0
    AliWarning(Form("Vertex @ %f outside of range [%f,%f]", 
		    vz, fHEventsTr->GetXaxis()->GetXmin(), 
		    fHEventsTr->GetXaxis()->GetXmax()));
#endif
    return;
  }

  // Apply the sharing filter (or hit merging or clustering if you like)
  if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, vz)) { 
#ifdef VERBOSE
    AliWarning("Sharing filter failed!");
#endif
    return;
  }

  // Calculate the inclusive charged particle density 
  if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux)) { 
    AliWarning("Density calculator failed!");
    return;
  }
  
  // Do the secondary and other corrections. 
  if (!fCorrections.Correct(fHistos, ivz)) { 
    AliWarning("Corrections failed");
    return;
  }

  if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) {
    AliWarning("Histogram collector failed");
    return;
  }

  if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
    fHData->Add(&(fAODFMD.GetHistogram()));
}

//____________________________________________________________________
void
AliForwardMultiplicity::Terminate(Option_t*)
{
  TList* list = dynamic_cast<TList*>(GetOutputData(1));
  if (!list) {
    AliError("No output list defined");
    return;
  }
  // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
  TH1D* dNdeta = fHData->ProjectionX("dNdeta", 1, -1, "e");
  TH1D* norm   = fHData->ProjectionX("dNdeta", 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(fHEventsTrVtx->GetEntries())/fHEventsTr->GetEntries(),
		"width");

  list->Add(fHEventsTr);
  list->Add(fHEventsTrVtx);
  list->Add(fHTriggers);
  list->Add(fHData);
  list->Add(dNdeta);
  
  TList* last = new TList;
  last->SetName("LastEvent");
  list->Add(last);
  last->Add(&fAODFMD.GetHistogram());
  last->Add(fHistos.fFMD1i);
  last->Add(fHistos.fFMD2i);
  last->Add(fHistos.fFMD2o);
  last->Add(fHistos.fFMD3i);
  last->Add(fHistos.fFMD3o);


  fSharingFilter.ScaleHistograms(fHEventsTr->Integral());
  fSharingFilter.Output(list);

  fDensityCalculator.ScaleHistograms(fHEventsTrVtx->Integral());
  fDensityCalculator.Output(list);

  fCorrections.ScaleHistograms(fHEventsTrVtx->Integral());
  fCorrections.Output(list);
}

//____________________________________________________________________
void
AliForwardMultiplicity::MarkEventForStore() const
{
  // Make sure the AOD tree is filled 
  AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
  AliAODHandler*      ah = 
    dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
  if (!ah)  
    AliFatal("No AOD output handler set in analysis manager");

  ah->SetFillAOD(kTRUE);
}
//____________________________________________________________________
Bool_t
AliForwardMultiplicity::ReadTriggers(AliESDEvent* esd)
{
  // Get the analysis manager - should always be there 
  AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
  if (!am) { 
    AliWarning("No analysis manager defined!");
    return kFALSE;
  }

  // Get the input handler - should always be there 
  AliInputEventHandler* ih = 
    static_cast<AliInputEventHandler*>(am->GetInputEventHandler());
  if (!ih) { 
    AliWarning("No input handler");
    return kFALSE;
  }
  
  // Get the physics selection - add that by using the macro 
  // AddTaskPhysicsSelection.C 
  AliPhysicsSelection* ps = 
    static_cast<AliPhysicsSelection*>(ih->GetEventSelection());
  if (!ps) { 
    AliWarning("No physics selection");
    return kFALSE;
  }
  
  // Check if this is a collision candidate (INEL)
  Bool_t inel = ps->IsCollisionCandidate(esd);
  if (inel) { 
    fAODFMD.SetTriggerBits(AliAODForwardMult::kInel);
    fHTriggers->Fill(.5);
  }
  

  // IF this is inel, see if we have a tracklet 
  if (inel) { 
    const AliMultiplicity* spdmult = esd->GetMultiplicity();
    if (!spdmult) {
      AliWarning("No SPD multiplicity");
    }
    else { 
      Int_t n = spdmult->GetNumberOfTracklets();
      for (Int_t j = 0; j < n; j++) { 
	if(TMath::Abs(spdmult->GetEta(j)) < 1) { 
	  fAODFMD.SetTriggerBits(AliAODForwardMult::kInelGt0);
	  fHTriggers->Fill(1.5);
	  break;
	}
      }
    }
  }

  // Analyse some trigger stuff 
  AliTriggerAnalysis ta;
  if (ta.IsOfflineTriggerFired(esd, AliTriggerAnalysis::kNSD1)) {
    fAODFMD.SetTriggerBits(AliAODForwardMult::kNSD);
    fHTriggers->Fill(2.5);
  }

  // Get trigger stuff 
  TString triggers = esd->GetFiredTriggerClasses();
  if (triggers.Contains("CBEAMB-ABCE-NOPF-ALL")) {
    fAODFMD.SetTriggerBits(AliAODForwardMult::kEmpty);
    fHTriggers->Fill(3.5);
  }

  if (triggers.Contains("CINT1A-ABCE-NOPF-ALL")) {
    fAODFMD.SetTriggerBits(AliAODForwardMult::kA);
    fHTriggers->Fill(4.5);
  }

  if (triggers.Contains("CINT1B-ABCE-NOPF-ALL")) {
    fAODFMD.SetTriggerBits(AliAODForwardMult::kB);
    fHTriggers->Fill(5.5);
  }


  if (triggers.Contains("CINT1C-ABCE-NOPF-ALL")) {
    fAODFMD.SetTriggerBits(AliAODForwardMult::kC);
    fHTriggers->Fill(6.5);
  }

  if (triggers.Contains("CINT1-E-NOPF-ALL")) {
    fAODFMD.SetTriggerBits(AliAODForwardMult::kE);
    fHTriggers->Fill(7.5);
  }

  return kTRUE;
}
//____________________________________________________________________
Bool_t
AliForwardMultiplicity::ReadVertex(AliESDEvent* esd, Double_t& vz)
{
  // Get the vertex 
  const AliESDVertex* vertex = esd->GetPrimaryVertexSPD();
  if (!vertex) { 
#ifdef VERBOSE
    AliWarning("No SPD vertex found in ESD");
#endif
    return kFALSE;
  }

  // Check that enough tracklets contributed 
  if(vertex->GetNContributors() <= 0) {
#ifdef VERBOSE
    AliWarning(Form("Number of contributors to vertex is %d<=0",
		    vertex->GetNContributors()));
#endif
    return kFALSE;
  }

  // Check that the uncertainty isn't too large 
  if (vertex->GetZRes() > 0.1) { 
#ifdef VERBOSE
    AliWarning(Form("Uncertaintity in Z of vertex is too large %f > 0.1", 
		    vertex->GetZRes()));
#endif
    return kFALSE;
  }

  // Get the z coordiante 
  vz = vertex->GetZ();
	       
  return kTRUE;
}

//____________________________________________________________________
void
AliForwardMultiplicity::Print(Option_t*) const
{
}

//
// EOF
//
Commit	Line	Data
7e4038b5	1	#include "AliForwardMultiplicity.h"
	2	#include "AliTriggerAnalysis.h"
	3	#include "AliPhysicsSelection.h"
	4	#include "AliLog.h"
	5	#include "AliFMDAnaParameters.h"
	6	#include "AliESDEvent.h"
	7	#include "AliAODHandler.h"
	8	#include "AliMultiplicity.h"
	9	#include "AliInputEventHandler.h"
	10	#include <TH1.h>
	11	#include <TDirectory.h>
	12	#include <TTree.h>
	13
	14	//====================================================================
	15	AliForwardMultiplicity::AliForwardMultiplicity()
	16	: AliAnalysisTaskSE(),
	17	fHEventsTr(0),
	18	fHEventsTrVtx(0),
	19	fHTriggers(0),
	20	fHData(0),
	21	fFirstEvent(true),
	22	fLowFluxCut(1000),
	23	fESDFMD(),
	24	fHistos(),
	25	fAODFMD(),
	26	fSharingFilter(),
	27	fDensityCalculator(),
	28	fCorrections(),
	29	fHistCollector(),
	30	fList(0),
	31	fTree(0)
	32	{
	33	}
	34
	35	//____________________________________________________________________
	36	AliForwardMultiplicity::AliForwardMultiplicity(const char* name)
	37	: AliAnalysisTaskSE(name),
	38	fHEventsTr(0),
	39	fHEventsTrVtx(0),
	40	fHTriggers(0),
	41	fHData(0),
	42	fFirstEvent(true),
	43	fLowFluxCut(1000),
	44	fESDFMD(),
	45	fHistos(),
	46	fAODFMD(kTRUE),
	47	fSharingFilter("sharing"),
	48	fDensityCalculator("density"),
	49	fCorrections("corrections"),
	50	fHistCollector("collector"),
	51	fList(0),
	52	fTree(0)
	53	{
	54	DefineOutput(1, TList::Class());
	55	// DefineOutput(2, TTree::Class());
	56	}
	57
	58	//____________________________________________________________________
	59	AliForwardMultiplicity::AliForwardMultiplicity(const AliForwardMultiplicity& o)
	60	: AliAnalysisTaskSE(o),
	61	fHEventsTr(o.fHEventsTr),
	62	fHEventsTrVtx(o.fHEventsTrVtx),
	63	fHTriggers(o.fHTriggers),
	64	fHData(o.fHData),
65	fFirstEvent(true),
66	fLowFluxCut(1000),
67	fESDFMD(o.fESDFMD),
68	fHistos(o.fHistos),
69	fAODFMD(o.fAODFMD),
70	fSharingFilter(o.fSharingFilter),
71	fDensityCalculator(o.fDensityCalculator),
72	fCorrections(o.fCorrections),
73	fHistCollector(o.fHistCollector),
74	fList(o.fList),
75	fTree(o.fTree)
76	{
77	}
78
79	//____________________________________________________________________
80	AliForwardMultiplicity&
81	AliForwardMultiplicity::operator=(const AliForwardMultiplicity& o)
82	{
83	fHEventsTr = o.fHEventsTr;
84	fHEventsTrVtx = o.fHEventsTrVtx;
85	fHTriggers = o.fHTriggers;
86	fHData = o.fHData;
87	fFirstEvent = o.fFirstEvent;
88	fSharingFilter = o.fSharingFilter;
89	fDensityCalculator = o.fDensityCalculator;
90	fCorrections = o.fCorrections;
91	fHistCollector = o.fHistCollector;
92	fHistos = o.fHistos;
93	fAODFMD = o.fAODFMD;
94	fList = o.fList;
95	fTree = o.fTree;
96
97	return *this;
98	}
99
100	//____________________________________________________________________
101	void
102	AliForwardMultiplicity::Init()
103	{
104	fFirstEvent = true;
105	}
106
107	//____________________________________________________________________
108	void
109	AliForwardMultiplicity::InitializeSubs()
110	{
111	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
112	pars->Init(kTRUE);
113
114	fHEventsTr = new TH1I("nEvents", "Number of events w/trigger",
115	pars->GetNvtxBins(),
116	-pars->GetVtxCutZ(),
117	pars->GetVtxCutZ());
118	fHEventsTr->SetXTitle("v_{z} [cm]");
119	fHEventsTr->SetYTitle("# of events");
120	fHEventsTr->SetFillColor(kRed+1);
121	fHEventsTr->SetFillStyle(3001);
122	fHEventsTr->SetDirectory(0);
123	// fHEventsTr->Sumw2();
124
125	fHEventsTrVtx = new TH1I("nEventsTrVtx",
126	"Number of events w/trigger and vertex",
127	pars->GetNvtxBins(),
128	-pars->GetVtxCutZ(),
129	pars->GetVtxCutZ());
130	fHEventsTrVtx->SetXTitle("v_{z} [cm]");
131	fHEventsTrVtx->SetYTitle("# of events");
132	fHEventsTrVtx->SetFillColor(kBlue+1);
133	fHEventsTrVtx->SetFillStyle(3001);
134	fHEventsTrVtx->SetDirectory(0);
135	// fHEventsTrVtx->Sumw2();
136
137
138	fHTriggers = new TH1I("triggers", "Triggers", 10, 0, 10);
139	fHTriggers->SetFillColor(kRed+1);
140	fHTriggers->SetFillStyle(3001);
141	fHTriggers->SetStats(0);
142	fHTriggers->SetDirectory(0);
143	fHTriggers->GetXaxis()->SetBinLabel(1,"INEL");
144	fHTriggers->GetXaxis()->SetBinLabel(2,"INEL>0");
145	fHTriggers->GetXaxis()->SetBinLabel(3,"NSD");
146	fHTriggers->GetXaxis()->SetBinLabel(4,"Empty");
147	fHTriggers->GetXaxis()->SetBinLabel(5,"A");
148	fHTriggers->GetXaxis()->SetBinLabel(6,"B");
149	fHTriggers->GetXaxis()->SetBinLabel(7,"C");
150	fHTriggers->GetXaxis()->SetBinLabel(8,"E");
151
152	TAxis e(pars->GetNetaBins(), pars->GetEtaMin(), pars->GetEtaMax());
153	fHistos.Init(e);
154	fAODFMD.Init(e);
155
156	fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
157	fHData->SetStats(0);
158	fHData->SetDirectory(0);
159	fSharingFilter.Init();
160	fHistCollector.Init(*(fHEventsTr->GetXaxis()));
161	}
162
163	//____________________________________________________________________
164	void
165	AliForwardMultiplicity::UserCreateOutputObjects()
166	{
167	fList = new TList;
168
169	AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
170	AliAODHandler* ah =
171	dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
172	if (!ah)
173	AliFatal("No AOD output handler set in analysis manager");
174
175
176	TObject* obj = &fAODFMD;
177	ah->AddBranch("AliAODForwardMult", &obj);
178
179	// fTree = new TTree("T", "T");
180	// fTree->Branch("forward", &fAODFMD);
181
182	PostData(1, fList);
183	// PostData(2, fTree);
184	}
185	//____________________________________________________________________
186	void
187	AliForwardMultiplicity::UserExec(Option_t*)
188	{
189	// Get the input data
190	AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
191	if (!esd) {
192	AliWarning("No ESD event found for input event");
193	return;
194	}
195
196	#if 0
197	static Int_t nEvents = 0;
198	nEvents++;
199	if (nEvents % 100 == 0) AliInfo(Form("Event # %6d", nEvents));
200	#endif
201
202	// On the first event, initialize the parameters
203	if (fFirstEvent) {
204	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
205	pars->SetParametersFromESD(esd);
206	pars->PrintStatus();
207	fFirstEvent = false;
208
209	InitializeSubs();
210	}
211	// Clear stuff
212	fHistos.Clear();
213	fESDFMD.Clear();
214	fAODFMD.Clear();
215
216	// Read trigger information from the ESD and store in AOD object
217	if (!ReadTriggers(esd)) {
218	#ifdef VERBOSE
219	AliWarning("Failed to read triggers from ESD");
220	#endif
221	return;
222	}
223
224	// Mark this event for storage
225	MarkEventForStore();
226
227	// Check if this is a high-flux event
228	const AliMultiplicity* testmult = esd->GetMultiplicity();
229	if (!testmult) {
230	#ifdef VERBOSE
231	AliWarning("No central multiplicity object found");
232	#endif
233	return;
234	}
235	Bool_t lowFlux = testmult->GetNumberOfTracklets() < fLowFluxCut;
236
237	// Get the FMD ESD data
238	AliESDFMD* esdFMD = esd->GetFMDData();
239	if (!esdFMD) {
240	#ifdef VERBOSE
241	AliWarning("No FMD data found in ESD");
242	#endif
243	return;
244	}
245
246	// Get the vertex information
247	Double_t vz = 0;
248	Bool_t vzOk = ReadVertex(esd, vz);
249
250	fHEventsTr->Fill(vz);
251	if (!vzOk) {
252	#ifdef VERBOSE
253	AliWarning("Failed to read vertex from ESD");
254	#endif
255	return;
256	}
257	fHEventsTrVtx->Fill(vz);
258
259	// Get the vertex bin
260	Int_t ivz = fHEventsTr->GetXaxis()->FindBin(vz)-1;
261	fAODFMD.SetIpZ(vz);
262	if (ivz < 0 \|\| ivz >= fHEventsTr->GetXaxis()->GetNbins()) {
263	#if 0
264	AliWarning(Form("Vertex @ %f outside of range [%f,%f]",
265	vz, fHEventsTr->GetXaxis()->GetXmin(),
266	fHEventsTr->GetXaxis()->GetXmax()));
267	#endif
268	return;
269	}
270
271	// Apply the sharing filter (or hit merging or clustering if you like)
272	if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, vz)) {
273	#ifdef VERBOSE
274	AliWarning("Sharing filter failed!");
275	#endif
276	return;
277	}
278
279	// Calculate the inclusive charged particle density
280	if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux)) {
281	AliWarning("Density calculator failed!");
282	return;
283	}
284
285	// Do the secondary and other corrections.
286	if (!fCorrections.Correct(fHistos, ivz)) {
287	AliWarning("Corrections failed");
288	return;
289	}
290
291	if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) {
292	AliWarning("Histogram collector failed");
293	return;
294	}
295
296	if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
297	fHData->Add(&(fAODFMD.GetHistogram()));
298	}
299
300	//____________________________________________________________________
301	void
302	AliForwardMultiplicity::Terminate(Option_t*)
303	{
304	TList* list = dynamic_cast<TList*>(GetOutputData(1));
305	if (!list) {
306	AliError("No output list defined");
307	return;
308	}
309	// TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
310	TH1D* dNdeta = fHData->ProjectionX("dNdeta", 1, -1, "e");
311	TH1D* norm = fHData->ProjectionX("dNdeta", 0, 1, "");
312	dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
313	dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
314	dNdeta->Divide(norm);
315	dNdeta->SetStats(0);
316	dNdeta->Scale(Double_t(fHEventsTrVtx->GetEntries())/fHEventsTr->GetEntries(),
317	"width");
318
319	list->Add(fHEventsTr);
320	list->Add(fHEventsTrVtx);
321	list->Add(fHTriggers);
322	list->Add(fHData);
323	list->Add(dNdeta);
324
325	TList* last = new TList;
326	last->SetName("LastEvent");
327	list->Add(last);
328	last->Add(&fAODFMD.GetHistogram());
329	last->Add(fHistos.fFMD1i);
330	last->Add(fHistos.fFMD2i);
331	last->Add(fHistos.fFMD2o);
332	last->Add(fHistos.fFMD3i);
333	last->Add(fHistos.fFMD3o);
334
335
336	fSharingFilter.ScaleHistograms(fHEventsTr->Integral());
337	fSharingFilter.Output(list);
338
339	fDensityCalculator.ScaleHistograms(fHEventsTrVtx->Integral());
340	fDensityCalculator.Output(list);
341
342	fCorrections.ScaleHistograms(fHEventsTrVtx->Integral());
343	fCorrections.Output(list);
344	}
345
346	//____________________________________________________________________
347	void
348	AliForwardMultiplicity::MarkEventForStore() const
349	{
350	// Make sure the AOD tree is filled
351	AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
352	AliAODHandler* ah =
353	dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
354	if (!ah)
355	AliFatal("No AOD output handler set in analysis manager");
356
357	ah->SetFillAOD(kTRUE);
358	}
359	//____________________________________________________________________
360	Bool_t
361	AliForwardMultiplicity::ReadTriggers(AliESDEvent* esd)
362	{
363	// Get the analysis manager - should always be there
364	AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
365	if (!am) {
366	AliWarning("No analysis manager defined!");
367	return kFALSE;
368	}
369
370	// Get the input handler - should always be there
371	AliInputEventHandler* ih =
372	static_cast<AliInputEventHandler*>(am->GetInputEventHandler());
373	if (!ih) {
374	AliWarning("No input handler");
375	return kFALSE;
376	}
377
378	// Get the physics selection - add that by using the macro
379	// AddTaskPhysicsSelection.C
380	AliPhysicsSelection* ps =
381	static_cast<AliPhysicsSelection*>(ih->GetEventSelection());
382	if (!ps) {
383	AliWarning("No physics selection");
384	return kFALSE;
385	}
386
387	// Check if this is a collision candidate (INEL)
388	Bool_t inel = ps->IsCollisionCandidate(esd);
389	if (inel) {
390	fAODFMD.SetTriggerBits(AliAODForwardMult::kInel);
391	fHTriggers->Fill(.5);
392	}
393
394
395	// IF this is inel, see if we have a tracklet
396	if (inel) {
397	const AliMultiplicity* spdmult = esd->GetMultiplicity();
398	if (!spdmult) {
399	AliWarning("No SPD multiplicity");
400	}
401	else {
402	Int_t n = spdmult->GetNumberOfTracklets();
403	for (Int_t j = 0; j < n; j++) {
404	if(TMath::Abs(spdmult->GetEta(j)) < 1) {
405	fAODFMD.SetTriggerBits(AliAODForwardMult::kInelGt0);
406	fHTriggers->Fill(1.5);
407	break;
408	}
409	}
410	}
411	}
412
413	// Analyse some trigger stuff
414	AliTriggerAnalysis ta;
415	if (ta.IsOfflineTriggerFired(esd, AliTriggerAnalysis::kNSD1)) {
416	fAODFMD.SetTriggerBits(AliAODForwardMult::kNSD);
417	fHTriggers->Fill(2.5);
418	}
419
420	// Get trigger stuff
421	TString triggers = esd->GetFiredTriggerClasses();
422	if (triggers.Contains("CBEAMB-ABCE-NOPF-ALL")) {
423	fAODFMD.SetTriggerBits(AliAODForwardMult::kEmpty);
424	fHTriggers->Fill(3.5);
425	}
426
427	if (triggers.Contains("CINT1A-ABCE-NOPF-ALL")) {
428	fAODFMD.SetTriggerBits(AliAODForwardMult::kA);
429	fHTriggers->Fill(4.5);
430	}
431
432	if (triggers.Contains("CINT1B-ABCE-NOPF-ALL")) {
433	fAODFMD.SetTriggerBits(AliAODForwardMult::kB);
434	fHTriggers->Fill(5.5);
435	}
436
437
438	if (triggers.Contains("CINT1C-ABCE-NOPF-ALL")) {
439	fAODFMD.SetTriggerBits(AliAODForwardMult::kC);
440	fHTriggers->Fill(6.5);
441	}
442
443	if (triggers.Contains("CINT1-E-NOPF-ALL")) {
444	fAODFMD.SetTriggerBits(AliAODForwardMult::kE);
445	fHTriggers->Fill(7.5);
446	}
447
448	return kTRUE;
449	}
450	//____________________________________________________________________
451	Bool_t
452	AliForwardMultiplicity::ReadVertex(AliESDEvent* esd, Double_t& vz)
453	{
454	// Get the vertex
455	const AliESDVertex* vertex = esd->GetPrimaryVertexSPD();
456	if (!vertex) {
457	#ifdef VERBOSE
458	AliWarning("No SPD vertex found in ESD");
459	#endif
460	return kFALSE;
461	}
462
463	// Check that enough tracklets contributed
464	if(vertex->GetNContributors() <= 0) {
465	#ifdef VERBOSE
466	AliWarning(Form("Number of contributors to vertex is %d<=0",
467	vertex->GetNContributors()));
468	#endif
469	return kFALSE;
470	}
471
472	// Check that the uncertainty isn't too large
473	if (vertex->GetZRes() > 0.1) {
474	#ifdef VERBOSE
475	AliWarning(Form("Uncertaintity in Z of vertex is too large %f > 0.1",
476	vertex->GetZRes()));
477	#endif
478	return kFALSE;
479	}
480
481	// Get the z coordiante
482	vz = vertex->GetZ();
483
484	return kTRUE;
485	}
486
487	//____________________________________________________________________
488	void
489	AliForwardMultiplicity::Print(Option_t*) const
490	{
491	}
492
493	//
494	// EOF
495	//