AliPhysics  c7b8e89 (c7b8e89)
AliForwardMCMultiplicityTask.cxx
Go to the documentation of this file.
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 //
17 #include "AliTriggerAnalysis.h"
18 #include "AliPhysicsSelection.h"
19 #include "AliLog.h"
20 #include "AliESDEvent.h"
21 #include "AliMCEvent.h"
22 #include "AliAODHandler.h"
23 #include "AliMultiplicity.h"
24 #include "AliInputEventHandler.h"
26 #include "AliAnalysisManager.h"
27 #include <TH1.h>
28 #include <TDirectory.h>
29 #include <TTree.h>
30 #include <TROOT.h>
31 #define MCAOD_SLOT 4
32 #define PRIMARY_SLOT 5
33 #ifdef POST_AOD
34 # define DEFINE(N,C) DefineOutput(N,C)
35 # define POST(N,O) PostData(N,O)
36 #else
37 # define DEFINE(N,C) do { } while(false)
38 # define POST(N,O) do { } while(false)
39 #endif
40 #ifndef ENABLE_TIMING
41 # define MAKE_SW(NAME) do {} while(false)
42 # define START_SW(NAME) do {} while(false)
43 # define FILL_SW(NAME,WHICH) do {} while(false)
44 #else
45 # define MAKE_SW(NAME) TStopwatch NAME
46 # define START_SW(NAME) if (fDoTiming) NAME.Start(true)
47 # define FILL_SW(NAME,WHICH) \
48  if (fDoTiming) fHTiming->Fill(WHICH,NAME.CpuTime())
49 #endif
50 
51 //====================================================================
54  fESDFMD(),
55  fMCESDFMD(),
56  fMCHistos(),
57  fMCAODFMD(),
58  fMCRingSums(),
59  fPrimary(0),
60  fEventInspector(),
61  // fMultEventClassifier(),
62  fESDFixer(),
63  fSharingFilter(),
64  fDensityCalculator(),
65  fCorrections(),
66  fHistCollector(),
67  fEventPlaneFinder()
68 {
69  //
70  // Constructor
71  //
72 }
73 
74 //____________________________________________________________________
77  fESDFMD(),
78  fMCESDFMD(),
79  fMCHistos(),
80  fMCAODFMD(kTRUE),
81  fMCRingSums(),
82  fPrimary(0),
83  fEventInspector("event"),
84  // fMultEventClassifier("multClass"),
85  fESDFixer("esdFizer"),
86  fSharingFilter("sharing"),
87  fDensityCalculator("density"),
88  fCorrections("corrections"),
89  fHistCollector("collector"),
90  fEventPlaneFinder("eventplane")
91 {
92  //
93  // Constructor
94  //
95  // Parameters:
96  // name Name of task
97  //
98  fPrimary = new TH2D("primary", "MC Primaries", 1,0,1,20,0,TMath::TwoPi());
99  fPrimary->SetXTitle("#eta");
100  fPrimary->SetYTitle("#varphi [radians]");
101  fPrimary->SetZTitle("d^{2}N_{ch}/d#etad#phi");
102  fPrimary->Sumw2();
103  fPrimary->SetStats(0);
104  fPrimary->SetDirectory(0);
105  DEFINE(MCAOD_SLOT,AliAODForwardMult::Class());
106  DEFINE(PRIM_SLOT, TH2D::Class());
107 }
108 
109 //____________________________________________________________________
110 void
112 {
115 }
116 
117 //____________________________________________________________________
118 void
120 {
121  //
122  // Create output objects
123  //
124  //
126 
127  TObject* mcobj = &fMCAODFMD;
128  ah->AddBranch("AliAODForwardMult", &mcobj);
129  ah->AddBranch("TH2D", &fPrimary);
130 }
131 
132 //____________________________________________________________________
133 Bool_t
135 {
136  // We do this to explicitly disable the noise corrector for MC
138 
141  POST(PRIM_SLOT, fPrimary);
142  return ret;
143 }
144 
145 //____________________________________________________________________
146 void
148 {
149  //
150  // Initialise the sub objects and stuff. Called on first event
151  //
152  //
154 
155  fMCHistos.Init(eta);
156  fMCAODFMD.Init(eta);
157  fMCRingSums.Init(eta);
158 
160  DMSG(fDebug,0,"Primary histogram rebinned to %d,%f,%f eta axis %d,%f,%f",
161  fPrimary->GetXaxis()->GetNbins(),
162  fPrimary->GetXaxis()->GetXmin(),
163  fPrimary->GetXaxis()->GetXmax(),
164  eta.GetNbins(),
165  eta.GetXmin(),
166  eta.GetXmax());
167 
168 
169  TList* mcRings = new TList;
170  mcRings->SetName("mcRingSums");
171  mcRings->SetOwner();
172  fList->Add(mcRings);
173 
174  mcRings->Add(fMCRingSums.Get(1, 'I'));
175  mcRings->Add(fMCRingSums.Get(2, 'I'));
176  mcRings->Add(fMCRingSums.Get(2, 'O'));
177  mcRings->Add(fMCRingSums.Get(3, 'I'));
178  mcRings->Add(fMCRingSums.Get(3, 'O'));
179  fMCRingSums.Get(1, 'I')->SetMarkerColor(AliForwardUtil::RingColor(1, 'I'));
180  fMCRingSums.Get(2, 'I')->SetMarkerColor(AliForwardUtil::RingColor(2, 'I'));
181  fMCRingSums.Get(2, 'O')->SetMarkerColor(AliForwardUtil::RingColor(2, 'O'));
182  fMCRingSums.Get(3, 'I')->SetMarkerColor(AliForwardUtil::RingColor(3, 'I'));
183  fMCRingSums.Get(3, 'O')->SetMarkerColor(AliForwardUtil::RingColor(3, 'O'));
184 }
185 
186 //____________________________________________________________________
187 Bool_t
189 {
190  if (fFirstEvent)
192  // Clear stuff
193  fHistos.Clear();
194  fESDFMD.Clear();
195  fAODFMD.Clear();
196  fAODEP.Clear();
197  fMCHistos.Clear();
198  fMCESDFMD.Clear();
199  fMCAODFMD.Clear();
200  fPrimary->Reset();
201  return true;
202 }
203 //____________________________________________________________________
204 Bool_t
206 {
207  //
208  // Process each event
209  //
210  // Parameters:
211  // option Not used
212  //
213  MAKE_SW(total);
214  MAKE_SW(individual);
215  START_SW(total);
216  START_SW(individual);
217 
218  // Read production details
219 
220  // Get the input data
221  AliMCEvent* mcEvent = MCEvent();
222  if (!mcEvent || !mcEvent->Stack()) return false;
223 
224  Bool_t lowFlux = kFALSE;
225  UInt_t triggers = 0;
226  UShort_t ivz = 0;
227  TVector3 ip(1024, 1024, 0);
228  Double_t cent = -1;
229  UShort_t nClusters = 0;
230  UInt_t found = fEventInspector.Process(&esd, triggers, lowFlux,
231  ivz, ip, cent, nClusters);
232  UShort_t ivzMC = 0;
233  TVector3 ipMC(1024, 1024, 0);
234  Double_t phiR = 0;
235  Double_t b = 0;
236  Double_t cMC = 0;
237  Int_t npart = 0;
238  Int_t nbin = 0;
239  // UInt_t foundMC =
240  fEventInspector.ProcessMC(mcEvent, triggers, ivzMC, ipMC, b, cMC,
241  npart, nbin, phiR);
242  fEventInspector.CompareResults(ip.Z(), ipMC.Z(), cent, cMC, b, npart, nbin);
243  // fMultEventClassifier.Process(&esd,&fAODRef);
244  FILL_SW(individual,kTimingEventInspector);
245 
246  //Store all events
248 
249  Bool_t isAccepted = true;
250  if (found & AliFMDEventInspector::kNoEvent) {
251  fHStatus->Fill(kStatusNoEvent);
252  isAccepted = false;
253  // return;
254  }
255  if (found & AliFMDEventInspector::kNoTriggers) {
256  fHStatus->Fill(kStatusNoTrigger);
257  isAccepted = false;
258  // return;
259  }
260  //MarkEventForStore();
261  // Always set the B trigger - each MC event _is_ a collision
262  triggers |= AliAODForwardMult::kB;
263  // Set trigger bits, and mark this event for storage
264  fAODFMD.SetTriggerBits(triggers);
267  fAODFMD.SetCentrality(cent);
268  fAODFMD.SetNClusters(nClusters);
269 
270  fMCAODFMD.SetTriggerBits(triggers);
273  fMCAODFMD.SetCentrality(cent);
274  fMCAODFMD.SetNClusters(nClusters);
275 
276  // Disable this check on SPD - will bias data
277  // if (found & AliFMDEventInspector::kNoSPD) isAccepted = false; // return;
278  if (found & AliFMDEventInspector::kNoFMD) {
279  fHStatus->Fill(kStatusNoFMD);
280  isAccepted = false;
281  // return;
282  }
283  if (found & AliFMDEventInspector::kNoVertex) {
284  fHStatus->Fill(kStatusNoVertex);
285  isAccepted = false;
286  // return;
287  }
288 
289  if (isAccepted) {
290  fAODFMD.SetIpZ(ip.Z());
291  fMCAODFMD.SetIpZ(ip.Z());
292  }
293  if (found & AliFMDEventInspector::kBadVertex) isAccepted = false; // return;
294 
295  // We we do not want to use low flux specific code, we disable it here.
296  if (!fEnableLowFlux) lowFlux = false;
297 
298  // Get FMD data
299  AliESDFMD* esdFMD = esd.GetFMDData();
300 
301  // Fix up the the ESD
302  GetESDFixer().Fix(*esdFMD, ip);
303 
304  // Apply the sharing filter (or hit merging or clustering if you like)
305  START_SW(individual);
306  if (isAccepted && !fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD,ip.Z())){
307  AliWarning("Sharing filter failed!");
309  return false;
310  }
311  if (!fSharingFilter.FilterMC(*esdFMD, *mcEvent, ip,fMCESDFMD,fPrimary)){
312  AliWarning("MC Sharing filter failed!");
314  return false;
315  }
316 
317  // Store some MC parameters in corners of histogram :-)
318  fPrimary->SetBinContent(0, 0, ipMC.Z());
319  fPrimary->SetBinContent(fPrimary->GetNbinsX()+1,0, phiR);
320  fPrimary->SetBinContent(fPrimary->GetNbinsX()+1,fPrimary->GetNbinsY(),cMC);
321 
322 
323  if (!isAccepted) {
324  // Exit on MC event w/o trigger, vertex, data - since there's no more
325  // to be done for MC
326  FILL_SW(individual,kTimingSharingFilter);
327  return false;
328  }
329 
330  //MarkEventForStore();
332  FILL_SW(individual,kTimingSharingFilter);
333 
334 
335  // Calculate the inclusive charged particle density
336  START_SW(individual);
337  if (!fDensityCalculator.Calculate(fESDFMD, fHistos, lowFlux, cent, ip)) {
338  AliWarning("Density calculator failed!");
340  return false;
341  }
343  AliWarning("MC Density calculator failed!");
345  return false;
346  }
348  FILL_SW(individual,kTimingDensityCalculator);
349 
351  START_SW(individual);
353  &(fAODFMD.GetHistogram()), &fHistos)){
354  AliWarning("Eventplane finder failed!");
356  }
357  FILL_SW(individual,kTimingEventPlaneFinder);
358  }
359 
360  // Do the secondary and other corrections.
361  START_SW(individual);
362  if (!fCorrections.Correct(fHistos, ivz)) {
363  AliWarning("Corrections failed");
365  return false;
366  }
367  if (!fCorrections.CorrectMC(fMCHistos, ivz)) {
368  AliWarning("MC Corrections failed");
370  return false;
371  }
373  FILL_SW(individual,kTimingCorrections);
374 
375  // Collect our 'super' histogram
377  START_SW(individual);
379  fRingSums,
380  ivz,
383  false,
384  add)) {
385  AliWarning("Histogram collector failed");
387  return false;
388  }
390  fMCRingSums,
391  ivz,
393  -1,
394  true,
395  add)) {
396  AliWarning("MC Histogram collector failed");
398  return false;
399  }
400  FILL_SW(individual,kTimingHistCollector);
401 #if 0
402  // Copy underflow bins to overflow bins - always full phi coverage
403  TH2& hMC = fMCAODFMD.GetHistogram();
404  Int_t nEta = hMC.GetNbinsX();
405  Int_t nY = hMC.GetNbinsY();
406  for (Int_t iEta = 1; iEta <= nEta; iEta++) {
407  hMC.SetBinContent(iEta, nY+1, hMC.GetBinContent(iEta, 0));
408  }
409 #endif
410 
411  if (add) {
412  fHData->Add(&(fAODFMD.GetHistogram()));
414  }
415  else {
416  fHStatus->Fill(kStatusNotAdded);
417  }
418  FILL_SW(total,kTimingTotal);
419 
420  return true;
421 }
422 
423 //____________________________________________________________________
424 Bool_t
426 {
430  return ret;
431 }
432 
433 //____________________________________________________________________
434 void
436  TList* output) const
437 {
439  MakeRingdNdeta(input, "mcRingSums", output, "mcRingResults", 24);
440 }
441 
442 
443 
444 //
445 // EOF
446 //
Bool_t FilterMC(const AliESDFMD &input, const AliMCEvent &event, const TVector3 &ip, AliESDFMD &output, TH2D *primary)
void SetIpZ(Float_t ipZ)
virtual void CreateBranches(AliAODHandler *ah)
void Clear(Option_t *option="")
#define POST(N, O)
double Double_t
Definition: External.C:58
#define DEFINE(N, C)
void SetTriggerBits(UInt_t bits)
void SetSecondaryForMC(Bool_t use)
virtual Bool_t CompareResults(Double_t vz, Double_t trueVz, Double_t cent, Double_t mcC, Double_t b, Int_t npart, Int_t nbin)
void SetNClusters(UShort_t n)
void SetSNN(UShort_t sNN)
void Clear(Option_t *option="")
virtual Bool_t Event(AliESDEvent &esd)
const AliFMDMCTrackDensity & GetTrackDensity() const
AliFMDMCDensityCalculator fDensityCalculator
void SetCentrality(Float_t c)
#define MAKE_SW(NAME)
#define DMSG(L, N, F,...)
virtual void MarkEventForStore() const
void Clear(Option_t *option="")
virtual Bool_t CompareResults(AliForwardUtil::Histos &esd, AliForwardUtil::Histos &mc)
virtual void CreateBranches(AliAODHandler *ah)
#define MCAOD_SLOT
virtual Bool_t Correct(AliForwardUtil::Histos &hists, UShort_t vtxBin)
UInt_t Process(const AliESDEvent *event, UInt_t &triggers, Bool_t &lowFlux, UShort_t &ivz, TVector3 &ip, Double_t &cent, UShort_t &nClusters)
int Int_t
Definition: External.C:63
static Color_t RingColor(UShort_t d, Char_t r)
void Init(const TAxis &etaAxis)
UShort_t GetCollisionSystem() const
virtual void EstimatedNdeta(const TList *input, TList *output) const
unsigned int UInt_t
Definition: External.C:33
static Bool_t IsTriggerBits(UInt_t bits, UInt_t trg)
virtual Bool_t CalculateMC(const AliESDFMD &fmd, AliForwardUtil::Histos &hists)
#define START_SW(NAME)
Bool_t Filter(const AliESDFMD &input, Bool_t lowFlux, AliESDFMD &output, Double_t zvtx)
#define FILL_SW(NAME, WHICH)
virtual Bool_t Collect(const AliForwardUtil::Histos &hists, AliForwardUtil::Histos &sums, UShort_t vtxBin, TH2D &out, Double_t cent=-1.0, Bool_t eta2phi=false, Bool_t add=true)
Definition: External.C:228
virtual Bool_t CompareResults(AliForwardUtil::Histos &esd, AliForwardUtil::Histos &mc)
void SetSystem(UShort_t sys)
void SetUseOnlyPrimary(Bool_t use)
void SetRecoNoiseFactor(Int_t f)
virtual void ReadProductionDetails(AliMCEvent *event)
virtual void EstimatedNdeta(const TList *input, TList *output) const
virtual Bool_t Calculate(const AliESDFMD &fmd, AliForwardUtil::Histos &hists, Bool_t lowFlux, Double_t cent=-1, const TVector3 &ip=TVector3(1024, 1024, 0))
virtual Bool_t CorrectMC(AliForwardUtil::Histos &hists, UShort_t vtxBin)
TH2D * Get(UShort_t d, Char_t r) const
virtual void InitMembers(const TAxis &pe, const TAxis &pv)
virtual void InitMembers(const TAxis &pe, const TAxis &pv)
Definition: External.C:220
void Init(const TAxis &etaAxis)
virtual void MakeRingdNdeta(const TList *input, const char *inName, TList *output, const char *outName, Int_t style=20) const
UInt_t ProcessMC(AliMCEvent *event, UInt_t &triggers, UShort_t &ivz, TVector3 &ip, Double_t &b, Double_t &c, Int_t &npart, Int_t &nbin, Double_t &phiR)
#define PRIMARY_SLOT
UShort_t GetEnergy() const
unsigned short UShort_t
Definition: External.C:28
Float_t GetCentrality() const
bool Bool_t
Definition: External.C:53
void CompareResults(const AliESDFMD &esd, const AliESDFMD &mc)
Bool_t FindEventplane(AliVEvent *esd, AliAODForwardEP &aodEp, TH2D *h, AliForwardUtil::Histos *hists)
void Fix(AliESDFMD &esd, const TVector3 &ip)
static void RebinEta(TH2D *hist, const TAxis &etaAxis)
const TH2D & GetHistogram() const