AliPhysics  1811c8f (1811c8f)
AliForwardMCFlowTaskQC.cxx
Go to the documentation of this file.
1 //
2 // Calculate flow on MC data in the forward and central regions using the Q cumulants method.
3 //
4 // Inputs:
5 // - AliAODEvent
6 //
7 // Outputs:
8 // - AnalysisResults.root
9 //
10 #include "AliForwardMCFlowTaskQC.h"
11 #include "AliAODMCParticle.h"
12 #include "AliAODMCHeader.h"
13 #include "TGraph.h"
14 #include "TF1.h"
15 #include "AliAODEvent.h"
16 #include "AliAODForwardMult.h"
17 #include "AliAODCentralMult.h"
18 #include "AliGenEventHeader.h"
19 #include "AliAnalysisManager.h"
20 #include "AliInputEventHandler.h"
21 #include "AliGenEventHeaderTunedPbPb.h"
22 #include "AliCollisionGeometry.h"
23 
24 ClassImp(AliForwardMCFlowTaskQC)
25 #if 0
26 ;
27 #endif
28 //_____________________________________________________________________
31  fBinsForwardTR(), // List of FMDTR analysis objects
32  fBinsCentralTR(), // List of SPDTR analysis objects
33  fBinsMC(), // List of MC particle analysis objects
34  fAODMCHeader(), // MC Header
35  fHistdNdedpMC(), // MC particle d^2N/detadphi histogram
36  fHistFMDMCCorr(), // FMD MC correlation
37  fHistSPDMCCorr(), // SPD MC correlation
38  fWeights(0), // Flow weights
39  fImpactParToCent(), // Impact parameter to centrality graph
40  fUseImpactPar(0), // Use impact par for centrality
41  fUseMCVertex(0), // Get vertex from MC header?
42  fUseFlowWeights(0) // Add flow to MC particles
43 {}
44  //
45  // Default Constructor
46  //
47 //_____________________________________________________________________
49  : AliForwardFlowTaskQC(name),
50  fBinsForwardTR(), // List of FMDTR analysis objects
51  fBinsCentralTR(), // List of SPDTR analysis objects
52  fBinsMC(), // List of MC particles analysis objects
53  fAODMCHeader(0), // MC Header
54  fHistdNdedpMC(), // MC particles d^2N/detadphi histogram
55  fHistFMDMCCorr(), // FMD MC correlation
56  fHistSPDMCCorr(), // SPD MC correlation
57  fWeights(0), // Flow weights
58  fImpactParToCent(), // Impact parameter to centrality graph
59  fUseImpactPar(0), // Use impact par for centrality
60  fUseMCVertex(0), // Get vertex from MC header?
61  fUseFlowWeights(0) // Add flow to MC particles
62 {
63  //
64  // Constructor
65  // Parameters:
66  // name: Name of task
67  //
68 
69  // Double_t impactParam[] = {0.,1.75,4.225,5.965,7.765,9.215,10.46,
70  // 11.565,12.575,13.515,16.679};
71  // Double_t centrality[] = {0.,2.5,7.5,15,25,35,45,55,65,75,90};
72  Double_t impactParam[] = { 0.00, 3.72, 5.23, 7.31, 8.88, 10.20,
73  11.38, 12.47, 13.50, 14.51, 16.679};
74  Double_t centrality[] = { 0., 5., 10., 20., 30., 40.,
75  50., 60., 70., 80., 100.};
76 
77  Int_t nPoints = sizeof(impactParam)/sizeof(Double_t);
78  fImpactParToCent = new TGraph(nPoints, impactParam, centrality);
79 
80 }
81 //_____________________________________________________________________
84  fBinsForwardTR(), // List of FMDTR analysis objects
85  fBinsCentralTR(), // List of SPDTR analysis objects
86  fBinsMC(), // List of MC particles analysis objects
87  fAODMCHeader(o.fAODMCHeader), // MC Header
88  fHistdNdedpMC(o.fHistdNdedpMC), // MC particles d^2N/detadphi histogram
89  fHistFMDMCCorr(o.fHistFMDMCCorr), // FMD MC correlation
90  fHistSPDMCCorr(o.fHistSPDMCCorr), // SPD MC correlation
91  fWeights(o.fWeights), // Flow weights
92  fImpactParToCent(o.fImpactParToCent), // Impact parameter to centrality
93  fUseImpactPar(o.fUseImpactPar), // Use impact par for centrality
94  fUseMCVertex(o.fUseMCVertex), // Get vertex from MC header?
95  fUseFlowWeights(o.fUseFlowWeights) // Add flow to MC particles
96 {
97  //
98  // Copy Constructor
99  //
100 }
101 //_____________________________________________________________________
104 {
105  //
106  // Assignment operator
107  // Parameters:
108  // o Object to copy from
109  //
110  if (&o == this) return *this;
115  fWeights = o.fWeights;
120  return *this;
121 }
122 //_____________________________________________________________________
124 {
125  //
126  // Initiate VertexBin lists
127  //
129 
130  Bool_t isNUA = (fFlowFlags & kNUAcorr);
131  for (Int_t v = 1; v <= fVtxAxis->GetNbins(); v++) {
132  Int_t vL = Int_t(fVtxAxis->GetBinLowEdge(v));
133  Int_t vH = Int_t(fVtxAxis->GetBinUpEdge(v));
134  // FMD
135  if ((fFlowFlags & kFMD)) {
136  fBinsForwardTR.Add(new VertexBin(vL, vH, fMaxMoment, "FMDTR", fFlowFlags, fFMDCut, fEtaGap));
137  if (!(fFlowFlags & k3Cor))
138  fBinsCentralTR.Add(new VertexBin(vL, vH, fMaxMoment, "SPDTR", fFlowFlags|kSPD, fSPDCut, fEtaGap));
139  if (isNUA) fFlowFlags ^= kNUAcorr;
140  fBinsMC.Add(new VertexBin(vL, vH, fMaxMoment, "MC-FMD", fFlowFlags|kMC, -1, fEtaGap));
141  if ((fFlowFlags & kStdQC))
142  fBinsMC.Add(new VertexBin(vL, vH, fMaxMoment, "MC-SPD", fFlowFlags|kMC|kSPD, -1, fEtaGap));
143  if (isNUA) fFlowFlags ^= kNUAcorr;
144  }
145  // VZERO
146  else if ((fFlowFlags & kVZERO)) {
147  fBinsMC.Add(new VertexBin(vL, vH, fMaxMoment, "MC-VZERO", fFlowFlags|kMC, -1, fEtaGap));
148  }
149  }
150 }
151 //_____________________________________________________________________
153 {
154  //
155  // Initiate diagnostics hists and add to outputlist
156  //
158 
159  TString subDetName = ((fFlowFlags & kFMD) ? "FMD" : ((fFlowFlags & kVZERO) ? "VZERO" : "none"));
160  fHistdNdedpMC = TH2D(Form("fdNdedpMC%s%s", subDetName.Data(), GetQCType(fFlowFlags)),
161  Form("fdNdedpMC%s%s", subDetName.Data(), GetQCType(fFlowFlags)),
162  240, -6., 6., 200, 0., TMath::TwoPi());
163 
164  fHistFMDMCCorr = new TH2D("hFMDMCCorr", "hFMDMCCorr", 200, 0., 15000., 200, 0, 20000);
165  fHistSPDMCCorr = new TH2D("hSPDMCCorr", "hSPDMCCorr", 200, 0., 7500., 200, 0, 20000);
166  TList* dList = (TList*)fSumList->FindObject("Diagnostics");
167  if (!dList) {
168  dList = new TList();
169  dList->SetName("Diagnostics");
170  fSumList->Add(dList);
171  }
172  dList->Add(fHistFMDMCCorr);
173  dList->Add(fHistSPDMCCorr);
174 
175  TIter nextForwardTR(&fBinsForwardTR);
176  VertexBin* bin = 0;
177  while ((bin = static_cast<VertexBin*>(nextForwardTR()))) {
179  }
180  TIter nextCentralTR(&fBinsCentralTR);
181  while ((bin = static_cast<VertexBin*>(nextCentralTR()))) {
183  }
184  TIter nextMC(&fBinsMC);
185  while ((bin = static_cast<VertexBin*>(nextMC()))) {
187  }
188  if (fWeights) {
189  TList* wList = new TList();
190  wList->SetName("FlowWeights");
191  fWeights->Init(wList);
192  fSumList->Add(wList);
193  }
194 
195 }
196 //_____________________________________________________________________
198 {
199  //
200  // Load FMD and SPD MC objects from aod tree and call the cumulants
201  // calculation for the correct vertexbin
202  //
203  if (!AliForwardFlowTaskQC::Analyze()) return kFALSE;
204 
205  // Run analysis on trackrefs from FMD and SPD
206  AliAODForwardMult* aodfmult = static_cast<AliAODForwardMult*>(fAOD->FindListObject("ForwardMC"));
207  AliAODCentralMult* aodcmult = static_cast<AliAODCentralMult*>(fAOD->FindListObject("CentralClustersMC"));
208  Int_t vtx = fVtxAxis->FindBin(fVtx)-1;
209 
210  // if objects are present, get histograms
211  if (aodfmult) {
212  TH2D& fmdTRdNdetadphi = aodfmult->GetHistogram();
213  if ((fFlowFlags & kStdQC)) {
214  FillVtxBinList(fBinsForwardTR, fmdTRdNdetadphi, vtx);
215  } else if ((fFlowFlags & kEtaGap)) {
216  FillVtxBinListEtaGap(fBinsForwardTR, fmdTRdNdetadphi, fmdTRdNdetadphi, vtx/*, kDoVtxCut*/);
217  }
218  if (aodcmult) {
219  TH2D& spdTRdNdetadphi = aodcmult->GetHistogram();
220  if ((fFlowFlags & kStdQC)) {
221  FillVtxBinList(fBinsCentralTR, spdTRdNdetadphi, vtx);
222  } else if ((fFlowFlags & kEtaGap)) {
223  FillVtxBinListEtaGap(fBinsCentralTR, fmdTRdNdetadphi, spdTRdNdetadphi, vtx/*, kDoVtxCut*/);
224  } else if ((fFlowFlags & k3Cor)) {
225  FillVtxBinList3Cor(fBinsForwardTR, spdTRdNdetadphi, fmdTRdNdetadphi, vtx);
226  }
227  }
228  }
229  // Run analysis on MC branch
230  if (!FillMCHist()) return kFALSE;
231 
232  if ((fFlowFlags & kStdQC)) {
234  } else if ((fFlowFlags & kEtaGap)) {
236  } else if ((fFlowFlags & k3Cor)) {
238  }
239 
240  // Mult correlation diagnostics
241  if (aodfmult && aodcmult) {
242  AliAODForwardMult* fmult = static_cast<AliAODForwardMult*>(fAOD->FindListObject("Forward"));
243  AliAODCentralMult* cmult = static_cast<AliAODCentralMult*>(fAOD->FindListObject("CentralClusters"));
244  const TH2D& fhist = fmult->GetHistogram();
245  const TH2D& chist = cmult->GetHistogram();
246 
247  Double_t totForward = fhist.Integral(1, fhist.GetNbinsX(), 1, fhist.GetNbinsY());
248  Double_t totSPD = chist.Integral(1, chist.GetNbinsX(), 1, chist.GetNbinsY());
249  Double_t totMC = fHistdNdedpMC.Integral(1, fHistdNdedpMC.GetNbinsX(), 1, fHistdNdedpMC.GetNbinsY());
250  fHistFMDMCCorr->Fill(totForward, totMC);
251  fHistSPDMCCorr->Fill(totSPD, totMC);
252  }
253 
254  return kTRUE;
255 }
256 //_____________________________________________________________________
258 {
259  //
260  // Finalize command, called by Terminate()
261  //
263 
267 
268  return;
269 }
270 //_____________________________________________________________________
272 {
273  //
274  // Function to check that an AOD event meets the cuts
275  //
276  // Parameters:
277  // AliAODForwardMult: forward mult object with trigger and vertex info
278  //
279  // Return: false if there is no trigger or if the centrality or vertex
280  // is out of range. Otherwise true.
281  //
282  fAODMCHeader = static_cast<AliAODMCHeader*>(fAOD->FindListObject(AliAODMCHeader::StdBranchName()));
283  return AliForwardFlowTaskQC::CheckEvent(aodfm);
284 }
285 //_____________________________________________________________________
287 {
288  //
289  // Function to look for a trigger string in the event.
290  //
291  // Parameters:
292  // AliAODForwardMult: forward mult object with trigger and vertex info
293  //
294  // Returns true if B trigger is present - for some reason this is the one we use in MC
295  //
296  if (aodfm) return aodfm->IsTriggerBits(AliAODForwardMult::kB);
297  else return (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
298  ->IsEventSelected() & AliVEvent::kMB);
299 
300 }
301 // _____________________________________________________________________
303 {
304  //
305  // Function to use centrality parametrization from impact parameter
306  // if flag is not set call AliForwardFlowTaskQC::GetCentrality
307  //
308  // Parameters:
309  // AliAODForwardMult: forward mult object with trigger and vertex info
310  //
311  // Returns true when centrality is set.
312  //
313  AliGenEventHeaderTunedPbPb* header =
314  dynamic_cast<AliGenEventHeaderTunedPbPb*>(fAODMCHeader->GetCocktailHeader(0));
315  if (header) fCent = header->GetCentrality();
316  else if (fUseImpactPar) fCent = GetCentFromB();
317  else return AliForwardFlowTaskQC::GetCentrality(aodfm);
318 
319  if (fCentAxis->GetXmin() > fCent || fCent >= fCentAxis->GetXmax()) {
320  fHistEventSel->Fill(kInvCent);
321  return kFALSE;
322  }
323  if (fCent != 0) return kTRUE;
324  else {
325  fHistEventSel->Fill(kNoCent);
326  return kFALSE;
327  }
328 }
329 //_____________________________________________________________________
331 {
332  //
333  // Function to look for vertex determination in the event using the MC header.
334  //
335  // Parameters:
336  // AliAODForwardMult: Not used
337  //
338  // Returns true if vertex is determined
339  //
340  if (fUseMCVertex) {
341  if (fAODMCHeader) {
342  fVtx = fAODMCHeader->GetVtxZ();
343  if (fVtx < fVtxAxis->GetXmin() || fVtx > fVtxAxis->GetXmax()) {
344  fHistEventSel->Fill(kInvVtx);
345  return kFALSE;
346  }
347  return kTRUE;
348  } else {
349  fHistEventSel->Fill(kNoVtx);
350  return kFALSE;
351  }
352  } else
353  return AliForwardFlowTaskQC::GetVertex(aodfm);
354 }
355 //_____________________________________________________________________
357 {
358  //
359  // Loop over AliAODParticle branch and fill d^2N/detadphi-histograms.
360  // Add flow if set to do so in AddTask function
361  //
362  fHistdNdedpMC.Reset();
363  Double_t minEta = -3.75;
364  Double_t maxEta = 5.;
365 
366  //retreive MC particles from event
367  TClonesArray* mcArray =
368  static_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
369  if(!mcArray){
370  AliWarning("No MC array found in AOD. Try making it again.");
371  return kFALSE;
372  }
373 
374  if (!fAODMCHeader) AliWarning("No MC header found.");
375 
376  Int_t ntracks = mcArray->GetEntriesFast();
377  Double_t rp = -1, b = -1;
378  if (fAODMCHeader) {
379  rp = fAODMCHeader->GetReactionPlaneAngle();
380  b = fAODMCHeader->GetImpactParameter();
381  if (fAODMCHeader->GetNCocktailHeaders() > 1) {
382  ntracks = fAODMCHeader->GetCocktailHeader(0)->NProduced();
383  }
384  }
385 
386  for (Int_t it = 0; it < ntracks; it++) { // Track loop
387  Double_t weight = 1;
388  AliAODMCParticle* particle = (AliAODMCParticle*)mcArray->At(it);
389  if (!particle) {
390  AliError(Form("Could not receive track %d", it));
391  continue;
392  }
393  if (!particle->IsPrimary()) continue;
394  if (particle->Charge() == 0) continue;
395  Double_t pT = particle->Pt();
396  Double_t eta = particle->Eta();
397  Double_t phi = particle->Phi();
398  if (eta >= minEta && eta < maxEta) {
399  // Add flow if it is in the argument
400  if (fUseFlowWeights && fWeights) {
401  weight = fWeights->CalcWeight(eta, pT, phi, particle->PdgCode(), rp, b);
402 // Printf("%f", weight);
403  }
404  fHistdNdedpMC.Fill(eta, phi, weight);
405  }
406  }
407  // Set underflow bins for acceptance checks
408  Int_t sBin = fHistdNdedpMC.GetXaxis()->FindBin(minEta);
409  Int_t eBin = fHistdNdedpMC.GetXaxis()->FindBin(maxEta);
410  for ( ; sBin <= eBin; sBin++) {
411  fHistdNdedpMC.SetBinContent(sBin, 0, 1);
412  } // End of eta bin loop
413 
414  return kTRUE;
415 }
416 //_____________________________________________________________________
418 {
419  //
420  // Get centrality from MC impact parameter.
421  //
422  Double_t cent = -1.;
423  if (!fAODMCHeader) return cent;
424  Double_t b = fAODMCHeader->GetImpactParameter();
425  cent = fImpactParToCent->Eval(b);
426 
427  return cent;
428 }
429 //_____________________________________________________________________
430 //
431 //
432 // EOF
double Double_t
Definition: External.C:58
void AddOutput(TList *list, TAxis *centAxis)
void EndVtxBinList(const TList &list) const
virtual void Init(TList *l)
centrality
virtual Bool_t GetVertex(const AliAODForwardMult *aodfm)
virtual Bool_t GetVertex(const AliAODForwardMult *aodfm)
void FillVtxBinList3Cor(const TList &list, TH2D &hcent, TH2D &hfwd, Int_t vtx, UShort_t flags=0x0)
AliForwardFlowWeights * fWeights
Per-event per bin.
int Int_t
Definition: External.C:63
static Bool_t IsTriggerBits(UInt_t bits, UInt_t trg)
void FillVtxBinList(const TList &list, TH2D &h1, Int_t vtx, UShort_t flags=0x0) const
Definition: External.C:228
virtual Bool_t GetCentrality(const AliAODForwardMult *aodfm)
AliForwardMCFlowTaskQC & operator=(const AliForwardMCFlowTaskQC &o)
virtual Bool_t CheckEvent(const AliAODForwardMult *aodfm)
const TH2D & GetHistogram() const
virtual Bool_t GetCentrality(const AliAODForwardMult *aodfm)
virtual Bool_t CheckEvent(const AliAODForwardMult *aodfm)
Double_t CalcWeight(Double_t eta, Double_t pt, Double_t phi, Int_t id, Double_t phiR, Double_t bOrC, Int_t type, UShort_t order, UShort_t what) const
virtual Bool_t CheckTrigger(const AliAODForwardMult *aodfm) const
static const Char_t * GetQCType(UShort_t flags, Bool_t prependUS=kTRUE)
bool Bool_t
Definition: External.C:53
const TH2D & GetHistogram() const
void FillVtxBinListEtaGap(const TList &list, TH2D &href, TH2D &hdiff, Int_t vtx, UShort_t flags=0x0) const