]>
Commit | Line | Data |
---|---|---|
2b556440 | 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 "TProfile2D.h" | |
16 | #include "AliAODEvent.h" | |
17 | #include "AliAODForwardMult.h" | |
18 | #include "AliAODCentralMult.h" | |
008eda2b | 19 | #include "AliGenEventHeader.h" |
2b556440 | 20 | |
21 | ClassImp(AliForwardMCFlowTaskQC) | |
936b0a6c | 22 | #if 0 |
23 | ; | |
24 | #endif | |
2b556440 | 25 | //_____________________________________________________________________ |
936b0a6c | 26 | AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC() |
27 | : AliForwardFlowTaskQC(), | |
28 | fBinsFMDTR(), // List of FMDTR analysis objects | |
29 | fBinsSPDTR(), // List of SPDTR analysis objects | |
30 | fBinsMC(), // List of MC truth analysis objects | |
31 | fdNdedpMC(), // MC truth d^2N/detadphi histogram | |
32 | fWeights(), // Flow weights | |
33 | fImpactParToCent(), // Impact parameter to centrality graph | |
34 | fUseImpactPar(0), // Use impact par for centrality | |
008eda2b | 35 | fFMDMinEta(-6), // FMD min eta coverage for this vtx |
36 | fFMDMaxEta(6), // FMD max eta coverage for this vtx | |
936b0a6c | 37 | fAddFlow(0), // Add flow to MC truth |
38 | fAddType(0), // Add type of flow to MC truth | |
39 | fAddOrder(0) // Add order of flow to MC truth | |
40 | {} | |
2b556440 | 41 | // |
42 | // Default Constructor | |
43 | // | |
44 | //_____________________________________________________________________ | |
936b0a6c | 45 | AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC(const char* name) |
46 | : AliForwardFlowTaskQC(name), | |
47 | fBinsFMDTR(), // List of FMDTR analysis objects | |
48 | fBinsSPDTR(), // List of SPDTR analysis objects | |
49 | fBinsMC(), // List of MC truth analysis objects | |
50 | fdNdedpMC(), // MC truth d^2N/detadphi histogram | |
51 | fWeights(), // Flow weights | |
52 | fImpactParToCent(), // Impact parameter to centrality graph | |
53 | fUseImpactPar(0), // Use impact par for centrality | |
008eda2b | 54 | fFMDMinEta(-6), // FMD min eta coverage for this vtx |
55 | fFMDMaxEta(6), // FMD max eta coverage for this vtx | |
936b0a6c | 56 | fAddFlow(0), // Add flow to MC truth |
57 | fAddType(0), // Add type of flow to MC truth | |
58 | fAddOrder(0) // Add order of flow to MC truth | |
2b556440 | 59 | { |
60 | // | |
61 | // Constructor | |
62 | // Parameters: | |
63 | // name: Name of task | |
64 | // | |
2b556440 | 65 | |
936b0a6c | 66 | // Double_t impactParam[] = {0.,1.75,4.225,5.965,7.765,9.215,10.46, |
67 | // 11.565,12.575,13.515,16.679}; | |
68 | // Double_t centrality[] = {0.,2.5,7.5,15,25,35,45,55,65,75,90}; | |
69 | Double_t impactParam[] = { 0.00, 3.72, 5.23, 7.31, 8.88, 10.20, | |
70 | 11.38, 12.47, 13.50, 14.51, 16.679}; | |
71 | Double_t centrality[] = { 0., 5., 10., 20., 30., 40., | |
72 | 50., 60., 70., 80., 100.}; | |
2b556440 | 73 | |
74 | Int_t nPoints = sizeof(impactParam)/sizeof(Double_t); | |
75 | fImpactParToCent = new TGraph(nPoints, impactParam, centrality); | |
76 | ||
77 | } | |
78 | //_____________________________________________________________________ | |
936b0a6c | 79 | AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC(const AliForwardMCFlowTaskQC& o) |
80 | : AliForwardFlowTaskQC(o), | |
81 | fBinsFMDTR(), // List of FMDTR analysis objects | |
82 | fBinsSPDTR(), // List of SPDTR analysis objects | |
83 | fBinsMC(), // List of MC truth analysis objects | |
84 | fdNdedpMC(o.fdNdedpMC), // MC truth d^2N/detadphi histogram | |
85 | fWeights(o.fWeights), // Flow weights | |
86 | fImpactParToCent(o.fImpactParToCent), // Impact parameter to centrality | |
87 | fUseImpactPar(o.fUseImpactPar), // Use impact par for centrality | |
008eda2b | 88 | fFMDMinEta(o.fFMDMinEta), // FMD min eta coverage for this vtx |
89 | fFMDMaxEta(o.fFMDMaxEta), // FMD max eta coverage for this vtx | |
936b0a6c | 90 | fAddFlow(o.fAddFlow), // Add flow to MC truth |
91 | fAddType(o.fAddType), // Add type of flow to MC truth | |
92 | fAddOrder(o.fAddOrder) // Add order of flow to MC truth | |
93 | { | |
2b556440 | 94 | // |
95 | // Copy Constructor | |
96 | // | |
936b0a6c | 97 | } |
2b556440 | 98 | //_____________________________________________________________________ |
99 | AliForwardMCFlowTaskQC& | |
100 | AliForwardMCFlowTaskQC::operator=(const AliForwardMCFlowTaskQC& o) | |
101 | { | |
102 | // | |
103 | // Assignment operator | |
104 | // Parameters: | |
105 | // o Object to copy from | |
106 | // | |
eeb6c967 | 107 | if (&o == this) return *this; |
2b556440 | 108 | fdNdedpMC = o.fdNdedpMC; |
936b0a6c | 109 | fWeights = o.fWeights; |
2b556440 | 110 | fImpactParToCent = o.fImpactParToCent; |
4b5b52b7 | 111 | fUseImpactPar = o.fUseImpactPar; |
008eda2b | 112 | fFMDMinEta = o.fFMDMinEta; |
113 | fFMDMaxEta = o.fFMDMaxEta; | |
2b556440 | 114 | fAddFlow = o.fAddFlow; |
115 | fAddType = o.fAddType; | |
116 | fAddOrder = o.fAddOrder; | |
117 | return *this; | |
118 | } | |
119 | //_____________________________________________________________________ | |
120 | void AliForwardMCFlowTaskQC::InitVertexBins() | |
121 | { | |
122 | // | |
123 | // Initiate VertexBin lists | |
124 | // | |
125 | AliForwardFlowTaskQC::InitVertexBins(); | |
126 | ||
68589651 | 127 | Int_t moment = 0; |
128 | for(UShort_t n = 0; n < fV.GetSize(); n++) { | |
129 | moment = fV.At(n); | |
130 | for (Int_t v = 1; v <= fVtxAxis->GetNbins(); v++) { | |
a7989168 | 131 | Int_t vL = Int_t(fVtxAxis->GetBinLowEdge(v)); |
132 | Int_t vH = Int_t(fVtxAxis->GetBinUpEdge(v)); | |
68589651 | 133 | fBinsFMDTR.Add(new VertexBin(vL, vH, moment, "FMDTR", fFlowFlags, fFMDCut, fEtaGap)); |
134 | fBinsSPDTR.Add(new VertexBin(vL, vH, moment, "SPDTR", fFlowFlags, fSPDCut, fEtaGap)); | |
135 | fBinsMC.Add(new VertexBin(vL, vH, moment, "MC", fFlowFlags, -1, fEtaGap)); | |
2b556440 | 136 | } |
137 | } | |
2b556440 | 138 | } |
139 | //_____________________________________________________________________ | |
140 | void AliForwardMCFlowTaskQC::InitHists() | |
141 | { | |
142 | // | |
143 | // Initiate diagnostics hists and add to outputlist | |
144 | // | |
145 | AliForwardFlowTaskQC::InitHists(); | |
146 | ||
68589651 | 147 | fdNdedpMC = TH2D(Form("fdNdedpMC%s", ((fFlowFlags & kEtaGap) ? "_etaGap" : "")), |
148 | Form("fdNdedpMC%s", ((fFlowFlags & kEtaGap) ? "_etaGap" : "")), | |
149 | 48, -6., 6., 200, 0., 2.*TMath::Pi()); | |
150 | fdNdedpMC.Sumw2(); | |
151 | ||
2b556440 | 152 | TIter nextFMDTR(&fBinsFMDTR); |
153 | VertexBin* bin = 0; | |
154 | while ((bin = static_cast<VertexBin*>(nextFMDTR()))) { | |
155 | bin->AddOutput(fSumList); | |
156 | } | |
157 | TIter nextSPDTR(&fBinsSPDTR); | |
158 | while ((bin = static_cast<VertexBin*>(nextSPDTR()))) { | |
159 | bin->AddOutput(fSumList); | |
160 | } | |
161 | TIter nextMC(&fBinsMC); | |
162 | while ((bin = static_cast<VertexBin*>(nextMC()))) { | |
163 | bin->AddOutput(fSumList); | |
164 | } | |
165 | ||
3112bd03 | 166 | TList* wList = new TList(); |
167 | wList->SetName("FlowWeights"); | |
168 | fWeights.Init(wList); | |
169 | fSumList->Add(wList); | |
170 | ||
2b556440 | 171 | } |
172 | //_____________________________________________________________________ | |
173 | Bool_t AliForwardMCFlowTaskQC::Analyze() | |
174 | { | |
175 | // | |
176 | // Load FMD and SPD MC objects from aod tree and call the cumulants | |
177 | // calculation for the correct vertexbin | |
178 | // | |
179 | if (!AliForwardFlowTaskQC::Analyze()) return kFALSE; | |
180 | ||
181 | // Run analysis on trackrefs from FMD and SPD | |
936b0a6c | 182 | const AliAODForwardMult* aodfmult = |
183 | static_cast<AliAODForwardMult*>(fAOD->FindListObject("ForwardMC")); | |
184 | const AliAODCentralMult* aodcmult = | |
185 | static_cast<AliAODCentralMult*>(fAOD->FindListObject("CentralClustersMC")); | |
d420e249 | 186 | Int_t vtx = fVtxAxis->FindBin(fVtx)-1; |
2b556440 | 187 | |
4b5b52b7 | 188 | // if objects are present, get histograms |
d420e249 | 189 | if (aodfmult) { |
190 | const TH2D& fmdTRdNdetadphi = aodfmult->GetHistogram(); | |
68589651 | 191 | if ((fFlowFlags & kEtaGap)) { |
192 | FillVtxBinListEtaGap(fBinsFMDTR, fmdTRdNdetadphi, fmdTRdNdetadphi, vtx); | |
193 | } else { | |
194 | FillVtxBinList(fBinsFMDTR, fmdTRdNdetadphi, vtx); | |
195 | } | |
008eda2b | 196 | |
68589651 | 197 | if (aodcmult) { |
198 | const TH2D& spdTRdNdetadphi = aodcmult->GetHistogram(); | |
199 | if ((fFlowFlags & kEtaGap)) { | |
200 | FillVtxBinListEtaGap(fBinsSPDTR, fmdTRdNdetadphi, spdTRdNdetadphi, vtx); | |
201 | } else { | |
202 | FillVtxBinList(fBinsSPDTR, spdTRdNdetadphi, vtx); | |
203 | } | |
204 | } | |
d420e249 | 205 | } |
2b556440 | 206 | // Run analysis on MC branch |
207 | if (!LoopAODMC()) return kFALSE; | |
68589651 | 208 | if ((fFlowFlags & kEtaGap)) { |
209 | FillVtxBinListEtaGap(fBinsMC, fdNdedpMC, fdNdedpMC, vtx); | |
210 | } else { | |
211 | FillVtxBinList(fBinsMC, fdNdedpMC, vtx); | |
212 | } | |
2b556440 | 213 | |
214 | return kTRUE; | |
215 | } | |
216 | //_____________________________________________________________________ | |
217 | void AliForwardMCFlowTaskQC::Finalize() | |
218 | { | |
219 | // | |
220 | // Finalize command, called by Terminate() | |
221 | // | |
222 | AliForwardFlowTaskQC::Finalize(); | |
223 | ||
4b5b52b7 | 224 | EndVtxBinList(fBinsFMDTR); |
225 | EndVtxBinList(fBinsSPDTR); | |
226 | EndVtxBinList(fBinsMC); | |
2b556440 | 227 | |
4b5b52b7 | 228 | return; |
229 | } | |
230 | //_____________________________________________________________________ | |
d420e249 | 231 | Bool_t AliForwardMCFlowTaskQC::CheckTrigger(const AliAODForwardMult* aodfm) const |
232 | { | |
233 | // | |
234 | // Function to look for a trigger string in the event. | |
235 | // | |
236 | // Parameters: | |
237 | // AliAODForwardMult: forward mult object with trigger and vertex info | |
238 | // | |
239 | // Returns true if B trigger is present - for some reason this is the one we use in MC | |
240 | // | |
241 | return aodfm->IsTriggerBits(AliAODForwardMult::kB); | |
242 | } | |
243 | // _____________________________________________________________________ | |
4b5b52b7 | 244 | Bool_t AliForwardMCFlowTaskQC::GetCentrality(const AliAODForwardMult* aodfm) |
245 | { | |
246 | // | |
247 | // Function to use centrality parametrization from impact parameter | |
248 | // if flag is not set call AliForwardFlowTaskQC::GetCentrality | |
249 | // | |
250 | // Parameters: | |
251 | // AliAODForwardMult: forward mult object with trigger and vertex info | |
252 | // | |
253 | // Returns true when centrality is set. | |
254 | // | |
255 | if (fUseImpactPar) { | |
256 | fCent = GetCentFromB(); | |
d420e249 | 257 | if (fCent != -1) return kTRUE; |
2b556440 | 258 | } |
d420e249 | 259 | return AliForwardFlowTaskQC::GetCentrality(aodfm); |
2b556440 | 260 | } |
261 | //_____________________________________________________________________ | |
008eda2b | 262 | void AliForwardMCFlowTaskQC::GetFMDLimits() |
263 | { | |
264 | ||
265 | const AliAODForwardMult* aodfmult = | |
266 | static_cast<AliAODForwardMult*>(fAOD->FindListObject("Forward")); | |
267 | const TH2D& h = aodfmult->GetHistogram(); | |
268 | ||
269 | for (Int_t e = 1; ; e++) { | |
270 | if (h.GetBinContent(e, 0) != 0) { | |
271 | fFMDMinEta = h.GetXaxis()->GetBinLowEdge(e); | |
272 | break; | |
273 | } | |
274 | } | |
275 | for (Int_t e = h.GetNbinsX(); ; e--) { | |
276 | if (h.GetBinContent(e, 0) != 0) { | |
277 | fFMDMaxEta = h.GetXaxis()->GetBinLowEdge(e); | |
278 | break; | |
279 | } | |
280 | } | |
281 | ||
282 | return; | |
283 | } | |
284 | //_____________________________________________________________________ | |
2b556440 | 285 | Bool_t AliForwardMCFlowTaskQC::LoopAODMC() |
286 | { | |
287 | // | |
288 | // Loop over AliAODParticle branch and fill d^2N/detadphi-histograms. | |
289 | // Add flow if set to do so in AddTask function | |
68589651 | 290 | // |
2b556440 | 291 | fdNdedpMC.Reset(); |
008eda2b | 292 | GetFMDLimits(); |
2b556440 | 293 | |
294 | //retreive MC particles from event | |
936b0a6c | 295 | TClonesArray* mcArray = |
d420e249 | 296 | static_cast<TClonesArray*>(fAOD->FindListObject( |
297 | AliAODMCParticle::StdBranchName())); | |
2b556440 | 298 | if(!mcArray){ |
d420e249 | 299 | AliWarning("No MC array found in AOD. Try making it again."); |
2b556440 | 300 | return kFALSE; |
301 | } | |
302 | ||
936b0a6c | 303 | AliAODMCHeader* header = |
d420e249 | 304 | dynamic_cast<AliAODMCHeader*>(fAOD->FindListObject( |
305 | AliAODMCHeader::StdBranchName())); | |
936b0a6c | 306 | if (!header) |
2b556440 | 307 | AliWarning("No header file found."); |
008eda2b | 308 | |
9c686cb4 | 309 | // Double_t rp = header->GetReactionPlaneAngle(); |
2b556440 | 310 | |
311 | Int_t ntracks = mcArray->GetEntriesFast(); | |
008eda2b | 312 | // TODO: Make this bit smarter... |
313 | if (header->GetNCocktailHeaders() > 1) { | |
314 | ntracks = header->GetCocktailHeader(0)->NProduced(); | |
315 | } | |
2b556440 | 316 | |
936b0a6c | 317 | UShort_t flowFlags = 0; |
318 | if (fAddFlow.Length() > 1) { | |
319 | if (fAddFlow.Contains("pt")) flowFlags |= AliForwardFlowWeights::kPt; | |
320 | if (fAddFlow.Contains("pid")) flowFlags |= AliForwardFlowWeights::kPt; | |
321 | if (fAddFlow.Contains("eta")) flowFlags |= AliForwardFlowWeights::kEta; | |
322 | if (fAddFlow.Contains("cent")) flowFlags |= AliForwardFlowWeights::kCent; | |
323 | } | |
9c686cb4 | 324 | // Double_t b = header->GetImpactParameter(); |
936b0a6c | 325 | |
2b556440 | 326 | // Track loop: chek how many particles will be accepted |
327 | Double_t weight = 0; | |
328 | for (Int_t it = 0; it < ntracks; it++) { | |
329 | weight = 1; | |
330 | AliAODMCParticle* particle = (AliAODMCParticle*)mcArray->At(it); | |
331 | if (!particle) { | |
332 | AliError(Form("Could not receive track %d", it)); | |
333 | continue; | |
334 | } | |
008eda2b | 335 | if (!particle->IsPrimary()) continue; |
2b556440 | 336 | if (particle->Charge() == 0) continue; |
68589651 | 337 | //Double_t pT = particle->Pt(); |
2b556440 | 338 | Double_t eta = particle->Eta(); |
339 | Double_t phi = particle->Phi(); | |
008eda2b | 340 | if (eta >= fFMDMinEta && eta <= fFMDMaxEta) { |
2b556440 | 341 | // Add flow if it is in the argument |
68589651 | 342 | /* FLOW WEIGHTS DISABLED IN THE VERSION - COMING BACK SOON |
343 | if (flowFlags != 0) { | |
344 | // weight = fWeights.CalcWeight(eta, pT, phi, particle->PdgCode(), | |
345 | // rp, fCent, fAddType, fAddOrder, | |
346 | // flowFlags) + 1; | |
03f52e9d | 347 | weight = fWeights.CalcWeight(eta, pT, phi, particle->PdgCode(), |
348 | rp, b); | |
68589651 | 349 | // Printf("%f", weight); |
350 | }*/ | |
03f52e9d | 351 | fdNdedpMC.Fill(eta, phi, weight); |
2b556440 | 352 | } |
353 | } | |
008eda2b | 354 | Int_t sBin = fdNdedpMC.GetXaxis()->FindBin(fFMDMinEta); |
355 | Int_t eBin = fdNdedpMC.GetXaxis()->FindBin(fFMDMaxEta); | |
d420e249 | 356 | for ( ; sBin < eBin; sBin++) fdNdedpMC.SetBinContent(sBin, 0, 1.); |
2b556440 | 357 | |
358 | return kTRUE; | |
359 | } | |
360 | //_____________________________________________________________________ | |
2b556440 | 361 | Double_t AliForwardMCFlowTaskQC::GetCentFromB() const |
362 | { | |
363 | // | |
364 | // Get centrality from MC impact parameter. | |
2b556440 | 365 | // |
366 | Double_t cent = -1.; | |
936b0a6c | 367 | AliAODMCHeader* header = |
368 | static_cast<AliAODMCHeader*>(fAOD->FindListObject(AliAODMCHeader:: | |
369 | StdBranchName())); | |
2b556440 | 370 | if (!header) return cent; |
9c686cb4 | 371 | Double_t b = header->GetImpactParameter(); |
2b556440 | 372 | |
373 | cent = fImpactParToCent->Eval(b); | |
374 | ||
375 | return cent; | |
376 | } | |
377 | //_____________________________________________________________________ | |
d420e249 | 378 | void AliForwardMCFlowTaskQC::PrintFlowSetup() const |
379 | { | |
380 | // | |
381 | // Print the setup of the flow task | |
382 | // | |
383 | Printf("AliForwardMCFlowTaskQC::Print"); | |
384 | Printf("Number of bins in vertex axis:\t%d", fVtxAxis->GetNbins()); | |
385 | Printf("Range of vertex axis :\t[%3.1f,%3.1f]", | |
386 | fVtxAxis->GetXmin(), fVtxAxis->GetXmax()); | |
387 | printf("Doing flow analysis for :\t"); | |
68589651 | 388 | for (Int_t n = 0; n < fV.GetSize(); n++) printf("v%d ", fV.At(n)); |
d420e249 | 389 | printf("\n"); |
9c686cb4 | 390 | Printf("Satellite vertex flag :\t%s", ((fFlowFlags & kSatVtx) ? |
391 | "true" : "false")); | |
392 | Printf("Symmetrize ref. flow wrt eta = 0:\t%s", ((fFlowFlags & kSymEta) ? | |
393 | "true" : "false")); | |
394 | Printf("Use an eta-gap for ref. flow :\t%s", ((fFlowFlags & kEtaGap) ? | |
395 | "true" : "false")); | |
d420e249 | 396 | Printf("FMD sigma cut: :\t%f", fFMDCut); |
397 | Printf("SPD sigma cut: :\t%f", fSPDCut); | |
68589651 | 398 | if ((fFlowFlags & kEtaGap)) |
399 | Printf("Eta gap: :\t%f", fEtaGap); | |
d420e249 | 400 | } |
401 | //_____________________________________________________________________ | |
2b556440 | 402 | // |
403 | // | |
404 | // EOF | |
405 |