]>
Commit | Line | Data |
---|---|---|
d2bea14e | 1 | // |
2 | // Calculate flow in the forward regions using the Q cumulants method | |
3 | // | |
4 | // Inputs: | |
5 | // - AliAODEvent | |
6 | // | |
7 | // Outputs: | |
8 | // - AnalysisResults.root | |
9 | // | |
10 | // TODO! | |
11 | // - Add centrality stuff | |
12 | // END OF TODO! | |
13 | // | |
14 | #include <TROOT.h> | |
15 | #include <TSystem.h> | |
16 | #include <TInterpreter.h> | |
17 | #include <TChain.h> | |
18 | #include <TFile.h> | |
19 | #include <TList.h> | |
20 | #include <iostream> | |
21 | #include <TMath.h> | |
22 | #include "AliLog.h" | |
23 | #include "AliForwardFlowTaskQC.h" | |
24 | #include "AliAnalysisManager.h" | |
25 | #include "AliAODHandler.h" | |
26 | #include "AliAODInputHandler.h" | |
27 | #include "AliAODMCParticle.h" | |
28 | #include "AliForwardFlowBase.h" | |
29 | #include "AliAODForwardMult.h" | |
30 | ||
31 | ClassImp(AliForwardFlowTaskQC) | |
9d05ffeb | 32 | #if 0 |
33 | ; // For emacs | |
34 | #endif | |
d2bea14e | 35 | |
36 | AliForwardFlowTaskQC::AliForwardFlowTaskQC() | |
9d05ffeb | 37 | : fDebug(0), // Debug flag |
38 | fOutputList(0), // Output list | |
39 | fAOD(0), // AOD input event | |
40 | fMC(kFALSE), // MC flag | |
41 | fEtaBins(20) // # of eta bins in histograms | |
d2bea14e | 42 | { |
43 | // | |
44 | // Default constructor | |
45 | // | |
46 | } | |
47 | //_____________________________________________________________________ | |
48 | AliForwardFlowTaskQC::AliForwardFlowTaskQC(const char* name) : | |
49 | AliAnalysisTaskSE(name), | |
50 | fDebug(0), // Debug flag | |
51 | fOutputList(0), // Output list | |
52 | fAOD(0), // AOD input event | |
53 | fMC(kFALSE), // MC flag | |
54 | fEtaBins(20) // # of Eta bins | |
55 | { | |
56 | // | |
57 | // Constructor | |
58 | // | |
59 | // Parameters: | |
60 | // name: Name of task | |
61 | // | |
62 | for (Int_t n = 0; n <= 4; n++) fv[n] = kTRUE; | |
63 | DefineOutput(1, TList::Class()); | |
64 | } | |
65 | //_____________________________________________________________________ | |
66 | AliForwardFlowTaskQC::AliForwardFlowTaskQC(const AliForwardFlowTaskQC& o) : | |
67 | AliAnalysisTaskSE(o), | |
68 | fDebug(o.fDebug), // Debug flag | |
69 | fOutputList(o.fOutputList), // Output list | |
70 | fAOD(o.fAOD), // AOD input event | |
71 | fMC(o.fMC), // MC flag | |
72 | fEtaBins(o.fEtaBins) // # of Eta bins | |
73 | { | |
74 | // | |
75 | // Copy constructor | |
76 | // | |
77 | // Parameters: | |
78 | // o Object to copy from | |
79 | // | |
80 | for (Int_t n = 0; n <= 4; n++) fv[n] = o.fv[n]; | |
81 | DefineOutput(1, TList::Class()); | |
82 | } | |
83 | //_____________________________________________________________________ | |
84 | void AliForwardFlowTaskQC::CreateOutputObjects() | |
85 | { | |
86 | // | |
87 | // Create output objects | |
88 | // | |
89 | if (!fOutputList) | |
90 | fOutputList = new TList(); | |
91 | fOutputList->SetName("QCumulants"); | |
92 | ||
93 | if (fEtaBins % 20) fEtaBins = 20; | |
94 | ||
95 | // Histograms for cumulants analysis | |
96 | ||
97 | // We loop over flow histograms here to add different orders of harmonics | |
98 | for (Int_t n = 1; n <= 4; n++) { | |
99 | if (!fv[n]) continue; | |
100 | // Only one flow histogram is needed for each type of data; | |
101 | // x-axis is eta-bins with differential flow, integrated is in underflowbin | |
102 | // y-axis bin 1: (w_<2> * <2>).Re() | |
103 | // y-axis bin 2: (w_<2> * <2>).Im() | |
104 | // y-axis bin 3: w_<2> = M(M-1) | |
105 | // y-axis bin 4: (w_<2> * <2> * <2>).Re() | |
106 | // y-axis bin 5: (w_<2> * <2> * <2>).Im() | |
107 | // y-axis bin 6: (w_<2> * w_<2'> * <2> * <2'>).Re() | |
108 | // y-axis bin 7: (w_<2> * w_<2'> * <2> * <2'>).Im() | |
109 | // y-axis bin 8: w_<2> * w_<2'> | |
110 | // y-axis bin 9: (w_<4> * <4>).Re() | |
111 | // y-axis bin 10: (w_<4> * <4>).Im() | |
112 | // y-axis bin 11: w_<4> | |
113 | // y-axis bin 12: (w_<4> * <4> * <4>).Re() | |
114 | // y-axis bin 13: (w_<4> * <4> * <4>).Im() | |
115 | // y-axis bin 14: w_<2> * w_<4> | |
116 | // y-axis bin 15: (w_<2> * w_<4> * <2> * <4>).Re() | |
117 | // y-axis bin 16: (w_<2> * w_<4> * <2> * <4>).Im() | |
118 | // y-axis bin 17: w_<2> * w_<4'> | |
119 | // y-axis bin 18: (w_<2> * w_<4'> * <2> * <4'>).Re() | |
120 | // y-axis bin 19: (w_<2> * w_<4'> * <2> * <4'>).Im() | |
121 | // y-axis bin 20: w_<4> * w_<2'> | |
122 | // y-axis bin 21: (w_<4> * w_<2'> * <4> * <2'>).Re() | |
123 | // y-axis bin 22: (w_<4> * w_<2'> * <4> * <2'>).Im() | |
124 | // y-axis bin 23: w_<4> * w_<4'> | |
125 | // y-axis bin 24: (w_<4> * w_<4'> * <4> * <4'>).Re() | |
126 | // y-axis bin 25: (w_<4> * w_<4'> * <4> * <4'>).Im() | |
127 | // y-axis bin 26: Qn or pn.Re() = <<cos2phi or psi>> | |
128 | // y-axis bin 27: Qn or pn.Im() = <<sin2phi or psi>> | |
129 | // y-axis bin 28: M or mp | |
130 | // y-axis bin 29: (Qn*Qn-Q2n).Re() = <<cos(2(phi1 or psi1+phi2))>> | |
131 | // y-axis bin 30: (Qn*Qn-Q2n).Im() = <<sin(2(phi1 or psi1+phi2))>> | |
132 | // y-axis bin 31: <<cos(2(phi1 or psi1-phi2-phi3))>> | |
133 | // y-axis bin 32: <<sin(2(phi1 or psi1-phi2-phi3))>> | |
134 | // y-axis bin 33: M*(M-1)*(M-2) or similar for diff | |
135 | // y-axis bin 34: <<cos(2(psi1+phi2-phi3>> | |
136 | // y-axis bin 35: <<sin(2(psi1+phi2-phi3>> | |
137 | TH2D* hFlowHist = new TH2D(Form("hQ%dCumuHist", n), Form("hQ%dCumuHist", n), fEtaBins, -4, 6, 35, 0.5, 35.5); | |
138 | hFlowHist->Sumw2(); | |
139 | fOutputList->Add(hFlowHist); | |
140 | ||
141 | TH2D* hFlowHistMC = new TH2D(Form("hQ%dCumuHistMC", n), Form("hQ%dCumuHistMC", n), fEtaBins, -4, 6, 35, 0.5, 35.5); | |
142 | hFlowHistMC->Sumw2(); | |
143 | fOutputList->Add(hFlowHistMC); | |
144 | ||
145 | TH2D* hFlowHistTrRef = new TH2D(Form("hQ%dCumuHistTrRef", n), Form("hQ%dCumuHistTrRef", n), fEtaBins, -4, 6, 35, 0.5, 35.5); | |
146 | hFlowHistTrRef->Sumw2(); | |
147 | fOutputList->Add(hFlowHistTrRef); | |
148 | ||
149 | // Output histograms | |
150 | TH1D* hCumulant2Flow = new TH1D(Form("hQ%dCumulant2Flow", n), Form("hQ%dCumulant2Flow", n), fEtaBins, -4, 6); | |
151 | hCumulant2Flow->Sumw2(); | |
152 | fOutputList->Add(hCumulant2Flow); | |
153 | ||
154 | TH1D* hCumulant2FlowMC = new TH1D(Form("hQ%dCumulant2FlowMC", n),Form("hQ%dCumulant2FlowMC", n), fEtaBins, -4, 6); | |
155 | hCumulant2FlowMC->Sumw2(); | |
156 | fOutputList->Add(hCumulant2FlowMC); | |
157 | ||
158 | TH1D* hCumulant2FlowTrRef = new TH1D(Form("hQ%dCumulant2FlowTrRef", n), Form("hQ%dCumulant2FlowTrRef", n), fEtaBins, -4, 6); | |
159 | hCumulant2FlowTrRef->Sumw2(); | |
160 | fOutputList->Add(hCumulant2FlowTrRef); | |
161 | ||
162 | ||
163 | TH1D* hCumulant4Flow = new TH1D(Form("hQ%dCumulant4Flow", n), Form("hQ%dCumulant4Flow", n), fEtaBins, -4, 6); | |
164 | hCumulant4Flow->Sumw2(); | |
165 | fOutputList->Add(hCumulant4Flow); | |
166 | ||
167 | TH1D* hCumulant4FlowMC = new TH1D(Form("hQ%dCumulant4FlowMC", n), Form("hQ%dCumulant4FlowMC", n), fEtaBins, -4, 6); | |
168 | hCumulant4FlowMC->Sumw2(); | |
169 | fOutputList->Add(hCumulant4FlowMC); | |
170 | ||
171 | TH1D* hCumulant4FlowTrRef = new TH1D(Form("hQ%dCumulant4FlowTrRef", n), Form("hQ%dCumulant4FlowTrRef", n), fEtaBins, -4, 6); | |
172 | hCumulant4FlowTrRef->Sumw2(); | |
173 | fOutputList->Add(hCumulant4FlowTrRef); | |
174 | } | |
175 | ||
176 | // Single Event histograms | |
177 | TH1D* hdNdphiSE = new TH1D("hdNdphiSE","hdNdphiSE", 20, 0, 2*TMath::Pi()); | |
178 | hdNdphiSE->Sumw2(); | |
179 | fOutputList->Add(hdNdphiSE); | |
180 | ||
181 | TH2D* hdNdetadphiSE = new TH2D("hdNdetadphiSE", "hdNdetadphiSE", fEtaBins, -4, 6, 20, 0, 2*TMath::Pi()); | |
182 | hdNdetadphiSE->Sumw2(); | |
183 | fOutputList->Add(hdNdetadphiSE); | |
184 | ||
185 | TH1D* hdNdphiSEMC = new TH1D("hdNdphiSEMC","hdNdphiSEMC", 20, 0, 2*TMath::Pi()); | |
186 | hdNdphiSEMC->Sumw2(); | |
187 | fOutputList->Add(hdNdphiSEMC); | |
188 | ||
189 | TH2D* hdNdetadphiSEMC = new TH2D("hdNdetadphiSEMC", "hdNdetadphiSEMC", fEtaBins, -4, 6, 20, 0, 2*TMath::Pi()); | |
190 | hdNdetadphiSEMC->Sumw2(); | |
191 | fOutputList->Add(hdNdetadphiSEMC); | |
192 | ||
193 | TH1D* hdNdphiSETrRef = new TH1D("hdNdphiSETrRef","hdNdphiSETrRef", 20, 0, 2*TMath::Pi()); | |
194 | hdNdphiSETrRef->Sumw2(); | |
195 | fOutputList->Add(hdNdphiSETrRef); | |
196 | ||
197 | TH2D* hdNdetadphiSETrRef = new TH2D("hdNdetadphiSETrRef", "hdNdetadphiSETrRef", fEtaBins, -4, 6, 20, 0, 2*TMath::Pi()); | |
198 | hdNdetadphiSETrRef->Sumw2(); | |
199 | fOutputList->Add(hdNdetadphiSETrRef); | |
200 | ||
201 | PostData(1, fOutputList); | |
202 | } | |
203 | //_____________________________________________________________________ | |
204 | void AliForwardFlowTaskQC::UserExec(Option_t */*option*/) | |
205 | { | |
206 | // | |
207 | // Process each event | |
208 | // | |
209 | // Parameters: | |
210 | // option: Not used | |
211 | // | |
212 | ||
213 | // Get input event | |
214 | fAOD = dynamic_cast<AliAODEvent*>(InputEvent()); | |
215 | if (!fAOD) return; | |
216 | ||
217 | // Load histograms and reset from last event | |
218 | TH1D* dNdphi = (TH1D*)fOutputList->FindObject("hdNdphiSE"); | |
219 | TH2D* dNdetadphi = (TH2D*)fOutputList->FindObject("hdNdetadphiSE"); | |
220 | ||
221 | dNdphi->Reset(); | |
222 | dNdetadphi->Reset(); | |
223 | ||
224 | // Initiate FlowCommon and fill histograms | |
225 | AliForwardFlowBase* common = new AliForwardFlowBase(fOutputList); | |
226 | ||
227 | if (!common->LoopAODFMD(fAOD)) return; | |
9d05ffeb | 228 | // else if (!common->LoopAODSPD(fAOD)) return; |
229 | // if (!common->LoopAODFMDandSPD(fAOD)) return; | |
d2bea14e | 230 | |
231 | // Run analysis | |
232 | for (Int_t n = 1; n <= 4; n++) { | |
233 | if (fv[n]) | |
234 | CumulantsMethod("", n); | |
235 | } | |
236 | ||
237 | // Find out if there's any MC data present | |
238 | if (!fMC) { | |
239 | TClonesArray* mcArray = 0; | |
240 | mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName()); | |
241 | if (mcArray) fMC = kTRUE; | |
242 | } | |
243 | if (fMC) | |
244 | ProcessPrimary(); | |
245 | ||
246 | } | |
247 | //_____________________________________________________________________ | |
9d05ffeb | 248 | void AliForwardFlowTaskQC::CumulantsMethod(TString type = "", |
249 | Int_t harmonic = 2) | |
d2bea14e | 250 | { |
251 | // | |
252 | // Calculate the Q cumulant of order n | |
253 | // | |
254 | // Parameters: | |
255 | // type: Determines which histograms should be used | |
256 | // - "" = data histograms | |
257 | // - "TrRef" = track reference histograms | |
258 | // - "MC" = MC truth histograms | |
259 | // harmonic: Which harmonic to calculate | |
260 | // | |
261 | Double_t n = harmonic; | |
262 | ||
263 | // We get histograms depending on if it's real data or MC truth data | |
264 | TH2D* flowHist = (TH2D*)fOutputList->FindObject(Form("hQ%dCumuHist%s", harmonic, type.Data())); | |
265 | TH1D* dNdphi = (TH1D*)fOutputList->FindObject(Form("hdNdphiSE%s", type.Data())); | |
266 | TH2D* dNdetadphi = (TH2D*)fOutputList->FindObject(Form("hdNdetadphiSE%s", type.Data())); | |
267 | ||
268 | // We create the objects needed for the analysis | |
269 | Double_t Mult = dNdphi->GetBinContent(0); | |
270 | ||
271 | Double_t QnRe = 0, Q2nRe = 0, QnIm = 0, Q2nIm = 0; | |
272 | Double_t pnRe = 0, p2nRe = 0, qnRe = 0, qnnRe = 0, pnIm = 0, p2nIm = 0, qnIm = 0, qnnIm = 0; | |
273 | Double_t avg2 = 0, avg4 = 0, avg2p = 0, avg4p = 0; | |
274 | Double_t w2avg2sq = 0, w2pavg2psq = 0, w4avg4sq = 0, w4pavg4psq = 0; | |
275 | Double_t w2w2pavg2avg2p = 0, w2w4avg2avg4 = 0, w2pw4pavg2pavg4p = 0; | |
276 | Double_t w2w4pavg2avg4p = 0, w4w2pavg4avg2p = 0, w4w4pavg4avg4p = 0; | |
277 | Double_t CosPhi1Phi2 = 0, CosPhi1Phi2Phi3m = 0, CosPhi1Phi2Phi3p = 0; | |
278 | Double_t SinPhi1Phi2 = 0, SinPhi1Phi2Phi3m = 0, SinPhi1Phi2Phi3p = 0; | |
279 | Double_t Phii = 0; | |
280 | Double_t multi = 0, multp = 0, mp = 0, mq = 0; | |
281 | Double_t W2 = 0, W4 = 0, W2p = 0, W4p = 0; | |
282 | ||
283 | // We loop over the data 1 time! | |
284 | for (Int_t eta = 1; eta <= dNdetadphi->GetNbinsX(); eta++) { | |
285 | // The values for each individual eta bins are reset | |
286 | mp = 0; | |
287 | pnRe = 0; | |
288 | p2nRe = 0; | |
289 | pnIm = 0; | |
290 | p2nIm = 0; | |
291 | ||
292 | for (Int_t phii = 1; phii <= dNdphi->GetNbinsX()+1; phii++) { | |
293 | Phii = dNdphi->GetXaxis()->GetBinCenter(phii); | |
294 | multi = dNdphi->GetBinContent(phii); | |
295 | ||
296 | // In the phi loop on the first eta loop the integrated flow | |
297 | // is calculated from the dNdphi histogram | |
298 | if(eta == 1) { | |
299 | QnRe += multi*TMath::Cos(n*Phii); | |
300 | QnIm += multi*TMath::Sin(n*Phii); | |
301 | Q2nRe += multi*TMath::Cos(2.*n*Phii); | |
302 | Q2nIm += multi*TMath::Sin(2.*n*Phii); | |
303 | } | |
304 | ||
305 | // For each eta bin the necesarry values for differential flow | |
306 | // is calculated. Here is the loop over the phi's. | |
307 | multp = dNdetadphi->GetBinContent(eta, phii); | |
308 | mp += multp; | |
309 | pnRe += multp*TMath::Cos(n*Phii); | |
310 | pnIm += multp*TMath::Sin(n*Phii); | |
311 | p2nRe += multp*TMath::Cos(2.*n*Phii); | |
312 | p2nIm += multp*TMath::Sin(2.*n*Phii); | |
313 | } | |
314 | ||
315 | // The integrated flow is calculated | |
316 | if (eta == 1) { | |
317 | Double_t Eta = flowHist->GetXaxis()->GetBinCenter(0); | |
318 | ||
319 | // 2-particle | |
320 | W2 = Mult * (Mult - 1.); | |
321 | avg2 = QnRe*QnRe + QnIm*QnIm - Mult; | |
322 | avg2 /= W2; | |
323 | w2avg2sq = W2 * avg2 * avg2; | |
324 | ||
325 | flowHist->Fill(Eta, 1, W2 * avg2); | |
326 | flowHist->Fill(Eta, 3, W2); | |
327 | flowHist->Fill(Eta, 4, w2avg2sq); | |
328 | ||
329 | flowHist->Fill(Eta, 26, QnRe); | |
330 | flowHist->Fill(Eta, 27, QnIm); | |
331 | flowHist->Fill(Eta, 28, Mult); | |
332 | ||
333 | // 4-particle | |
334 | W4 = Mult * (Mult - 1.) * (Mult - 2.) * (Mult - 3.); | |
335 | Double_t real = Q2nRe*QnRe*QnRe - Q2nRe*QnIm*QnIm + 2.*Q2nIm*QnRe*QnIm; | |
336 | ||
337 | avg4 = TMath::Power(QnRe*QnRe + QnIm*QnIm, 2); | |
338 | avg4 += Q2nRe*Q2nRe + Q2nIm*Q2nIm - 2.*real; | |
339 | avg4 -= 4.*(Mult - 2.)*(QnRe*QnRe + QnIm*QnIm) - 2.*Mult*(Mult - 3.); | |
340 | ||
341 | avg4 /= W4; | |
342 | w4avg4sq = W4 * avg4 * avg4; | |
343 | w2w4avg2avg4 = W2 * W4 * avg2 * avg4; | |
344 | ||
345 | flowHist->Fill(Eta, 9, W4 * avg4); | |
346 | flowHist->Fill(Eta, 11, W4); | |
347 | flowHist->Fill(Eta, 12, w4avg4sq); | |
348 | flowHist->Fill(Eta, 14, W2 * W4); | |
349 | flowHist->Fill(Eta, 15, w2w4avg2avg4); | |
350 | ||
351 | CosPhi1Phi2 = QnRe*QnRe - QnIm*QnIm - Q2nRe; | |
352 | SinPhi1Phi2 = 2.*QnRe*QnIm - Q2nIm; | |
353 | ||
354 | CosPhi1Phi2Phi3m = TMath::Power(QnRe, 3) + QnRe*QnIm*QnIm; | |
355 | CosPhi1Phi2Phi3m -= QnRe*Q2nRe + QnIm*Q2nIm + 2.*(Mult - 1.)*QnRe; | |
356 | SinPhi1Phi2Phi3m = -TMath::Power(QnIm, 3) - QnRe*QnRe*QnIm; | |
357 | SinPhi1Phi2Phi3m -= QnIm*Q2nRe - QnRe*Q2nIm + 2.*(Mult - 1.)*QnIm; | |
358 | ||
359 | flowHist->Fill(Eta, 29, CosPhi1Phi2); | |
360 | flowHist->Fill(Eta, 30, SinPhi1Phi2); | |
361 | flowHist->Fill(Eta, 31, CosPhi1Phi2Phi3m); | |
362 | flowHist->Fill(Eta, 32, SinPhi1Phi2Phi3m); | |
363 | flowHist->Fill(Eta, 33, Mult*(Mult-1.)*(Mult-2.)); | |
364 | ||
365 | // Count number of events | |
366 | flowHist->Fill(Eta, 0., 1.); | |
367 | } // end of harmonics loop | |
368 | ||
369 | // Differential flow calculations for each eta bin is done: | |
370 | if (mp == 0) continue; | |
371 | Double_t Eta = dNdetadphi->GetXaxis()->GetBinCenter(eta); | |
372 | ||
373 | mq = mp; | |
374 | qnRe = pnRe; | |
375 | qnIm = pnIm; | |
376 | qnnRe = p2nRe; | |
377 | qnnIm = p2nIm; | |
378 | ||
379 | // 2-particle differential flow | |
380 | W2p = mp * Mult - mq; | |
381 | avg2p = pnRe*QnRe + pnIm*QnIm - mq; | |
382 | avg2p /= W2p; | |
383 | w2pavg2psq = W2p * avg2p * avg2p; | |
384 | w2w2pavg2avg2p = W2 * W2p * avg2 * avg2p; | |
385 | ||
386 | flowHist->Fill(Eta, 1, W2p * avg2p); | |
387 | flowHist->Fill(Eta, 3, W2p); | |
388 | flowHist->Fill(Eta, 4, w2pavg2psq); | |
389 | flowHist->Fill(Eta, 6, w2w2pavg2avg2p); | |
390 | flowHist->Fill(Eta, 8, W2 * W2p); | |
391 | ||
392 | flowHist->Fill(Eta, 26, pnRe); | |
393 | flowHist->Fill(Eta, 27, pnIm); | |
394 | flowHist->Fill(Eta, 28, mp); | |
395 | ||
396 | // 4-particle differential flow | |
397 | W4p = (mp * Mult - 3.*mq)*(Mult - 1.)*(Mult - 2.); | |
398 | ||
399 | avg4p = pnRe*QnRe*(QnRe*QnRe + QnIm*QnIm) + pnIm*QnIm*(QnRe*QnRe + QnIm*QnIm); | |
400 | avg4p -= qnnRe*QnRe*QnRe - qnnRe*QnIm*QnIm + 2.*qnnIm*QnRe*QnIm; | |
401 | avg4p -= pnRe*QnRe*Q2nRe - pnRe*QnIm*Q2nIm + pnIm*QnRe*Q2nIm + pnIm*QnIm*Q2nRe; | |
402 | avg4p -= 2.*Mult*(pnRe*QnRe + pnIm*QnIm); | |
403 | ||
404 | avg4p += - 2.*mq*(QnRe*QnRe + QnIm*QnIm) + 7.*(qnRe*QnRe + qnIm*QnIm); | |
405 | avg4p += - (QnRe*qnRe + QnIm*qnIm) + (qnnRe*Q2nRe + qnnIm*Q2nIm); | |
406 | avg4p += 2.*(pnRe*QnRe + pnIm*QnIm) + 2.*mq*Mult - 6.*mq; | |
407 | avg4p /= W4p; | |
408 | ||
409 | w4pavg4psq = W4p * avg4p * avg4p; | |
410 | w2w4pavg2avg4p = W2 * W4p * avg2 * avg4p; | |
411 | w4w2pavg4avg2p = W4 * W2p * avg4 * avg2p; | |
412 | w4w4pavg4avg4p = W4 * W4p * avg4 * avg4p; | |
413 | w2pw4pavg2pavg4p = W2p * W4p * avg2p * avg4p; | |
414 | ||
415 | flowHist->Fill(Eta, 9, W4p * avg4p); | |
416 | flowHist->Fill(Eta, 11, W4p); | |
417 | flowHist->Fill(Eta, 12, w4pavg4psq); | |
418 | flowHist->Fill(Eta, 14, W2p * W4p); | |
419 | flowHist->Fill(Eta, 15, w2pw4pavg2pavg4p); | |
420 | flowHist->Fill(Eta, 17, W2 * W4p); | |
421 | flowHist->Fill(Eta, 18, w2w4pavg2avg4p); | |
422 | flowHist->Fill(Eta, 20, W4 * W2p); | |
423 | flowHist->Fill(Eta, 21, w4w2pavg4avg2p); | |
424 | flowHist->Fill(Eta, 23, W4 * W4p); | |
425 | flowHist->Fill(Eta, 24, w4w4pavg4avg4p); | |
426 | ||
427 | CosPhi1Phi2 = pnRe*QnRe - pnIm*QnIm - qnnRe; | |
428 | SinPhi1Phi2 = pnRe*QnIm + pnIm*QnRe - qnnIm; | |
429 | ||
430 | CosPhi1Phi2Phi3p = pnRe*(QnRe*QnRe + QnIm*QnIm - Mult); | |
431 | CosPhi1Phi2Phi3p -= qnnRe*QnRe - qnnIm*QnIm + mq*QnRe - 2.*qnRe; | |
432 | SinPhi1Phi2Phi3p = pnIm*(QnRe*QnRe + QnIm*QnIm - Mult); | |
433 | SinPhi1Phi2Phi3p -= qnnIm*QnRe - qnnRe*QnIm + mq*QnIm - 2.*qnIm; | |
434 | ||
435 | CosPhi1Phi2Phi3m = pnRe*(QnRe*QnRe - QnIm*QnIm) + 2.*pnIm*QnRe*QnIm; | |
436 | CosPhi1Phi2Phi3m -= pnRe*Q2nRe + pnIm*Q2nIm + 2.*mq*QnRe - 2.*qnRe; | |
437 | SinPhi1Phi2Phi3m = pnIm*(QnRe*QnRe - QnIm*QnIm) - 2.*pnRe*QnRe*QnIm; | |
438 | SinPhi1Phi2Phi3m += - pnIm*Q2nRe + pnRe*Q2nIm + 2.*mq*QnIm - 2.*qnIm; | |
439 | ||
440 | flowHist->Fill(Eta, 29, CosPhi1Phi2); | |
441 | flowHist->Fill(Eta, 30, SinPhi1Phi2); | |
442 | flowHist->Fill(Eta, 31, CosPhi1Phi2Phi3m); | |
443 | flowHist->Fill(Eta, 32, SinPhi1Phi2Phi3m); | |
444 | flowHist->Fill(Eta, 33, (mp*Mult-2.*mq)*(Mult-1.)); | |
445 | flowHist->Fill(Eta, 34, CosPhi1Phi2Phi3p); | |
446 | flowHist->Fill(Eta, 35, SinPhi1Phi2Phi3p); | |
447 | ||
448 | } | |
449 | ||
450 | } | |
451 | //_____________________________________________________________________ | |
452 | void AliForwardFlowTaskQC::Terminate(Option_t */*option*/) | |
453 | { | |
454 | // | |
455 | // End of job | |
456 | // Finalizes the Q cumulant calculations | |
457 | // | |
458 | // Parameters: | |
459 | // option Not used | |
460 | // | |
461 | ||
462 | TH2D* cumulantsHist; | |
463 | TH1D* cumulant2Hist; | |
464 | TH1D* cumulant4Hist; | |
465 | ||
466 | Int_t nLoops = (fMC ? 3 : 1); | |
467 | ||
468 | // Do a loop over the difference analysis types, calculating flow | |
469 | // 1 loop for real data, 3 for MC data | |
470 | // inside each is a nested loop over each harmonic (1, 2, 3 and 4 at the moment) | |
471 | for (Int_t loop = 1; loop <= nLoops; loop++) { | |
472 | ||
473 | TString type; | |
474 | if (loop == 1) type = ""; | |
475 | if (loop == 2) type = "MC"; | |
476 | if (loop == 3) type = "TrRef"; | |
477 | ||
478 | for (Int_t n = 1; n <= 4; n++) { | |
479 | if (!fv[n]) continue; | |
480 | ||
481 | cumulantsHist = (TH2D*)fOutputList->FindObject(Form("hQ%dCumuHist%s", n, type.Data())); | |
482 | cumulant2Hist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant2Flow%s", n, type.Data())); | |
483 | cumulant4Hist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant4Flow%s", n, type.Data())); | |
484 | ||
485 | // For flow calculations | |
486 | Double_t Avg2 = 0, c2 = 0, vTwo2 = 0, Avg4 = 0, c4 = 0, vTwo4 = 0; | |
487 | Double_t Avg2p = 0, d2 = 0, vTwo2diff = 0, Avg4p = 0, d4 = 0, vTwo4diff = 0; | |
488 | Double_t W2 = 0, W4 = 0, W2p = 0, W4p = 0, sqrtW2sq = 0, sqrtW2psq = 0, W2W2p = 0; | |
489 | Double_t W2W4 = 0, W2W4p = 0, W4W2p = 0, W4W4p = 0, W2pW4p = 0; | |
490 | Double_t sqrtW4sq = 0, sqrtW4psq = 0; | |
491 | Double_t W2avg2 = 0, W2pavg2p = 0, W4avg4 = 0, W4pavg4p = 0; | |
492 | Double_t AvgCos2Phi = 0, AvgSin2Phi = 0, Mult = 0, AvgCos2Phi1Phi2 = 0, AvgSin2Phi1Phi2 = 0; | |
493 | Double_t AvgCos2Phi1Phi2Phi3m = 0, AvgSin2Phi1Phi2Phi3m = 0, Multm1m2 = 0; | |
494 | Double_t AvgCos2Psi = 0, AvgSin2Psi = 0, mp = 0, AvgCos2Psi1Phi2 = 0, AvgSin2Psi1Phi2 = 0; | |
495 | Double_t AvgCos2Psi1Phi2Phi3m = 0, AvgSin2Psi1Phi2Phi3m = 0, mpqMult = 0; | |
496 | Double_t AvgCos2Psi1Phi2Phi3p = 0, AvgSin2Psi1Phi2Phi3p = 0; | |
497 | ||
498 | // For error calculations | |
499 | Double_t W2avg2sq = 0, W2W2pavg2avg2p = 0, W2pavg2psq = 0; | |
500 | Double_t W4avg4sq = 0, W2W4avg2avg4 = 0, W4W4pavg4avg4p = 0, W4pavg4psq = 0; | |
501 | Double_t W2W4pavg2avg4p = 0, W4W2pavg4avg2p = 0; | |
502 | Double_t sAvg2sq = 0, sAvg2psq = 0, sAvg4sq = 0, sAvg4psq = 0; | |
503 | Double_t vTwo2err = 0, vTwo2diffErr = 0, vTwo4err = 0, vTwo4diffErr = 0; | |
504 | Double_t Cov22p = 0, Cov24 = 0, Cov24p = 0, Cov42p = 0, Cov44p = 0, Cov2p2np = 0; | |
505 | Double_t W2pW4pavg2pavg4p = 0; | |
506 | ||
507 | for (Int_t eta = 0; eta <= cumulantsHist->GetNbinsX(); eta++) { | |
508 | if (eta == 0) { | |
509 | // 2-particle reference flow | |
510 | W2avg2 = cumulantsHist->GetBinContent(eta, 1); | |
511 | if (!W2avg2) continue; | |
512 | W2 = cumulantsHist->GetBinContent(eta, 3); | |
513 | AvgCos2Phi = cumulantsHist->GetBinContent(eta, 26); | |
514 | AvgSin2Phi = cumulantsHist->GetBinContent(eta, 27); | |
515 | Mult = cumulantsHist->GetBinContent(eta, 28); | |
516 | AvgCos2Phi /= Mult; | |
517 | AvgSin2Phi /= Mult; | |
518 | Avg2 = W2avg2 / W2; | |
519 | c2 = Avg2 - TMath::Power(AvgCos2Phi, 2) - TMath::Power(AvgSin2Phi, 2); | |
520 | vTwo2 = TMath::Sqrt(c2); | |
521 | cumulant2Hist->SetBinContent(eta, vTwo2); | |
522 | ||
523 | // 4-particle reference flow | |
524 | W4avg4 = cumulantsHist->GetBinContent(eta, 9); | |
525 | W4 = cumulantsHist->GetBinContent(eta, 11); | |
526 | AvgCos2Phi1Phi2 = cumulantsHist->GetBinContent(eta, 29); | |
527 | AvgSin2Phi1Phi2 = cumulantsHist->GetBinContent(eta, 30); | |
528 | AvgCos2Phi1Phi2 /= W2; | |
529 | AvgSin2Phi1Phi2 /= W2; | |
530 | AvgCos2Phi1Phi2Phi3m = cumulantsHist->GetBinContent(eta, 31); | |
531 | AvgSin2Phi1Phi2Phi3m = cumulantsHist->GetBinContent(eta, 32); | |
532 | Multm1m2 = cumulantsHist->GetBinContent(eta, 33); | |
533 | AvgCos2Phi1Phi2Phi3m /= Multm1m2; | |
534 | AvgSin2Phi1Phi2Phi3m /= Multm1m2; | |
535 | Avg4 = W4avg4 / W4; | |
536 | c4 = Avg4 - 2. * Avg2 * Avg2; | |
537 | c4 -= 4.*AvgCos2Phi*AvgCos2Phi1Phi2Phi3m; | |
538 | c4 += 4.*AvgSin2Phi*AvgSin2Phi1Phi2Phi3m; | |
539 | c4 -= TMath::Power(AvgCos2Phi1Phi2, 2) + TMath::Power(AvgSin2Phi1Phi2 , 2); | |
540 | c4 += 4.*AvgCos2Phi1Phi2*(TMath::Power(AvgCos2Phi, 2) - TMath::Power(AvgSin2Phi, 2)); | |
541 | c4 += 8.*AvgSin2Phi1Phi2*AvgSin2Phi*AvgCos2Phi; | |
542 | c4 += 8.*Avg2*(TMath::Power(AvgCos2Phi, 2) + TMath::Power(AvgSin2Phi, 2)); | |
543 | c4 -= 6.*TMath::Power(TMath::Power(AvgCos2Phi, 2)+TMath::Power(AvgSin2Phi, 2), 2); | |
544 | ||
545 | vTwo4 = TMath::Power(-c4, 0.25); | |
546 | cumulant4Hist->SetBinContent(eta, vTwo4); | |
547 | ||
548 | // 2-particle reference flow error calculations | |
549 | W2avg2sq = cumulantsHist->GetBinContent(eta, 4); | |
550 | sqrtW2sq = cumulantsHist->GetBinError(eta, 3); | |
551 | ||
552 | sAvg2sq = VarSQ(W2avg2sq, Avg2, W2, W2avg2, sqrtW2sq); | |
553 | vTwo2err = sqrtW2sq * TMath::Sqrt(sAvg2sq) / (2. * TMath::Sqrt(Avg2) * W2); | |
554 | cumulant2Hist->SetBinError(eta, vTwo2err); | |
555 | ||
556 | // 4-particle reference flow error calculations | |
557 | W4avg4sq = cumulantsHist->GetBinContent(eta, 12); | |
558 | sqrtW4sq = cumulantsHist->GetBinError(eta, 11); | |
559 | W2W4 = cumulantsHist->GetBinContent(eta, 14); | |
560 | W2W4avg2avg4 = cumulantsHist->GetBinContent(eta, 15); | |
561 | ||
562 | sAvg4sq = VarSQ(W4avg4sq, Avg4, W4, W4avg4, sqrtW4sq); | |
563 | Cov24 = CovXY(W2W4avg2avg4, W2W4, Avg2*Avg4, W2, W4); | |
564 | ||
565 | vTwo4err = Avg2*Avg2 * TMath::Power(sqrtW2sq, 2) * sAvg2sq / (W2*W2); | |
566 | vTwo4err += TMath::Power(sqrtW4sq, 2) * sAvg4sq / (16. * W4*W4); | |
567 | vTwo4err -= Avg2 * W2W4 * Cov24 / (2. * W2 * W4); | |
568 | vTwo4err /= TMath::Power(2. * Avg2*Avg2 - Avg4, 1.5); | |
569 | vTwo4err = TMath::Sqrt(vTwo4err); | |
570 | cumulant4Hist->SetBinError(eta, vTwo4err); | |
571 | ||
572 | continue; | |
573 | } | |
574 | ||
575 | // 2-particle differential flow | |
576 | W2pavg2p = cumulantsHist->GetBinContent(eta, 1); | |
577 | if (!W2pavg2p) continue; | |
578 | W2p = cumulantsHist->GetBinContent(eta, 3); | |
579 | AvgCos2Psi = cumulantsHist->GetBinContent(eta, 26); | |
580 | AvgSin2Psi = cumulantsHist->GetBinContent(eta, 27); | |
581 | mp = cumulantsHist->GetBinContent(eta, 28); | |
582 | AvgCos2Psi /= mp; | |
583 | AvgSin2Psi /= mp; | |
584 | Avg2p = W2pavg2p / W2p; | |
585 | d2 = Avg2p - AvgCos2Psi*AvgCos2Phi - AvgSin2Psi*AvgSin2Phi; | |
586 | vTwo2diff = d2 / TMath::Sqrt(c2); | |
587 | cumulant2Hist->SetBinContent(eta, vTwo2diff); | |
588 | ||
589 | // 4-particle differential flow | |
590 | W4pavg4p = cumulantsHist->GetBinContent(eta, 9); | |
591 | W4p = cumulantsHist->GetBinContent(eta, 11); | |
592 | AvgCos2Psi1Phi2 = cumulantsHist->GetBinContent(eta, 29); | |
593 | AvgSin2Psi1Phi2 = cumulantsHist->GetBinContent(eta, 30); | |
594 | AvgCos2Psi1Phi2 /= W2p; | |
595 | AvgSin2Psi1Phi2 /= W2p; | |
596 | AvgCos2Psi1Phi2Phi3m = cumulantsHist->GetBinContent(eta, 31); | |
597 | AvgSin2Psi1Phi2Phi3m = cumulantsHist->GetBinContent(eta, 32); | |
598 | mpqMult = cumulantsHist->GetBinContent(eta, 33); | |
599 | AvgCos2Psi1Phi2Phi3m /= mpqMult; | |
600 | AvgSin2Psi1Phi2Phi3m /= mpqMult; | |
601 | AvgCos2Psi1Phi2Phi3p = cumulantsHist->GetBinContent(eta, 34); | |
602 | AvgSin2Psi1Phi2Phi3p = cumulantsHist->GetBinContent(eta, 35); | |
603 | AvgCos2Psi1Phi2Phi3p /= mpqMult; | |
604 | AvgSin2Psi1Phi2Phi3p /= mpqMult; | |
605 | ||
606 | Avg4p = W4pavg4p / W4p; | |
607 | d4 = Avg4p - 2. * Avg2p * Avg2; | |
608 | d4 -= AvgCos2Psi*AvgCos2Phi1Phi2Phi3m; | |
609 | d4 += AvgSin2Psi*AvgSin2Phi1Phi2Phi3m; | |
610 | d4 -= AvgCos2Phi*AvgCos2Psi1Phi2Phi3m; | |
611 | d4 += AvgSin2Phi*AvgSin2Psi1Phi2Phi3m; | |
612 | d4 -= 2.*AvgCos2Phi*AvgCos2Psi1Phi2Phi3p; | |
613 | d4 -= 2.*AvgSin2Phi*AvgSin2Psi1Phi2Phi3p; | |
614 | d4 -= AvgCos2Psi1Phi2*AvgCos2Phi1Phi2; | |
615 | d4 -= AvgSin2Psi1Phi2*AvgSin2Phi1Phi2; | |
616 | d4 += 2.*AvgCos2Phi1Phi2*(AvgCos2Psi*AvgCos2Phi - AvgSin2Psi*AvgSin2Phi); | |
617 | d4 += 2.*AvgSin2Phi1Phi2*(AvgCos2Psi*AvgSin2Phi + AvgSin2Psi*AvgCos2Phi); | |
618 | d4 += 4.*Avg2*(AvgCos2Psi*AvgCos2Phi + AvgSin2Psi*AvgSin2Phi); | |
619 | d4 += 2.*AvgCos2Psi1Phi2*(TMath::Power(AvgCos2Phi, 2) - TMath::Power(AvgSin2Phi, 2)); | |
620 | d4 += 4.*AvgSin2Psi1Phi2*AvgCos2Phi*AvgSin2Phi; | |
621 | d4 += 4.*Avg2p*(TMath::Power(AvgCos2Phi, 2) + TMath::Power(AvgSin2Phi, 2)); | |
622 | d4 -= 6.*(TMath::Power(AvgCos2Phi, 2) - TMath::Power(AvgSin2Phi, 2)) | |
9d05ffeb | 623 | *(AvgCos2Psi*AvgCos2Phi-AvgSin2Psi*AvgSin2Phi); |
d2bea14e | 624 | d4 -= 12.*AvgCos2Phi*AvgSin2Phi*(AvgSin2Psi*AvgCos2Phi+AvgCos2Psi*AvgSin2Phi); |
625 | ||
626 | vTwo4diff = - d4 / TMath::Power(-c4, 0.75); | |
627 | cumulant4Hist->SetBinContent(eta, vTwo4diff); | |
628 | ||
629 | // 2-particle differential flow error calculations | |
630 | W2pavg2psq = cumulantsHist->GetBinContent(eta, 4); | |
631 | sqrtW2psq = cumulantsHist->GetBinError(eta, 3); | |
632 | W2W2pavg2avg2p = cumulantsHist->GetBinContent(eta, 6); | |
633 | W2W2p = cumulantsHist->GetBinContent(eta, 8); | |
634 | ||
635 | Cov22p = CovXY(W2W2pavg2avg2p, W2W2p, Avg2*Avg2p, W2, W2p); | |
636 | sAvg2psq = VarSQ(W2pavg2psq, Avg2p, W2p, W2pavg2p, sqrtW2psq); | |
637 | ||
638 | vTwo2diffErr = Avg2p*Avg2p*TMath::Power(sqrtW2psq, 2)*sAvg2sq/(W2*W2); | |
639 | vTwo2diffErr += 4.*Avg2*Avg2*TMath::Power(sqrtW2psq, 2)*sAvg2psq/(W2p*W2p); | |
640 | vTwo2diffErr -= 4.*Avg2*Avg2p*W2W2p*Cov22p/(W2*W2p); | |
641 | vTwo2diffErr /= (4. * TMath::Power(Avg2, 3)); | |
642 | vTwo2diffErr = TMath::Sqrt(vTwo2diffErr); | |
643 | cumulant2Hist->SetBinError(eta, vTwo2diffErr); | |
644 | ||
645 | // 4-particle differential flow error calculations | |
646 | sqrtW4psq = cumulantsHist->GetBinError(eta, 11); | |
647 | W4pavg4psq = cumulantsHist->GetBinContent(eta, 12); | |
648 | W2pW4p = cumulantsHist->GetBinContent(eta, 14); | |
649 | W2pW4pavg2pavg4p = cumulantsHist->GetBinContent(eta, 15); | |
650 | W2W4p = cumulantsHist->GetBinContent(eta, 17); | |
651 | W2W4pavg2avg4p = cumulantsHist->GetBinContent(eta, 18); | |
652 | W4W2p = cumulantsHist->GetBinContent(eta, 20); | |
653 | W4W2pavg4avg2p = cumulantsHist->GetBinContent(eta, 21); | |
654 | W4W4p = cumulantsHist->GetBinContent(eta, 23); | |
655 | W4W4pavg4avg4p = cumulantsHist->GetBinContent(eta, 24); | |
656 | ||
657 | sAvg4psq = VarSQ(W4pavg4psq, Avg4p, W4p, W4pavg4p, sqrtW4psq); | |
658 | Cov24p = CovXY(W2W4pavg2avg4p, W2W4p, Avg2*Avg4p, W2, W4p); | |
659 | Cov42p = CovXY(W4W2pavg4avg2p, W4W2p, Avg4*Avg2p, W4, W2p); | |
660 | Cov44p = CovXY(W4W4pavg4avg4p, W4W4p, Avg4*Avg4p, W4, W4p); | |
661 | Cov2p2np = CovXY(W2pW4pavg2pavg4p, W2pW4p, Avg2p*Avg4p, W2p, W4p); | |
662 | ||
663 | // Numbers on the side reference term number in paper (cite needed) loosely | |
9d05ffeb | 664 | /*1*/ vTwo4diffErr = TMath::Power(2.*Avg2*Avg2*Avg2p - 3.*Avg2*Avg4p + 2.*Avg4*Avg2p, 2) |
665 | * TMath::Power(sqrtW2sq, 2) * sAvg2sq / (W2*W2); | |
666 | /*2*/ vTwo4diffErr += 9. * TMath::Power(2.*Avg2*Avg2p - Avg4p, 2) * TMath::Power(sqrtW4sq, 2) | |
667 | * sAvg4sq / (16. * W4*W4); | |
668 | /*3*/ vTwo4diffErr += 4. * Avg2*Avg2 * TMath::Power(2.*Avg2*Avg2 - Avg4, 2) * TMath::Power(sqrtW2psq, 2) | |
669 | * sAvg2psq / (W2p*W2p); | |
670 | /*4*/ vTwo4diffErr += TMath::Power(2.*Avg2*Avg2 - Avg4, 2) * TMath::Power(sqrtW4psq, 2) * sAvg4psq | |
671 | / (W4p*W4p); | |
672 | /*5*/ vTwo4diffErr -= 1.5 * (2.*Avg2*Avg2p - Avg4p) * (2.*Avg2*Avg2*Avg2p - 3.*Avg2*Avg4p + 2.*Avg4*Avg2p) | |
673 | * W2W4 * Cov24 / (W2*W4); | |
674 | /*6*/ vTwo4diffErr -= 4. * Avg2 * (2.*Avg2*Avg2 - Avg4) | |
675 | * (2.*Avg2*Avg2*Avg2p - 3.*Avg2*Avg4p + 2.*Avg4*Avg2p) | |
676 | * W2W2p * Cov22p / (W2 * W2p); | |
677 | /*7*/ vTwo4diffErr += 2. * (2.*Avg2*Avg2 - Avg4) | |
678 | * (2.*Avg2*Avg2*Avg2p - 3.*Avg2*Avg4p + 2.*Avg4*Avg2p) | |
679 | * W2W4p * Cov24p / (W2 * W4p); | |
680 | /*8*/ vTwo4diffErr += 3.*Avg2*(2.*Avg2*Avg2 - Avg4)*(2.*Avg2*Avg2p - Avg4p) | |
681 | * W4W2p * Cov42p / (W4*W2p); | |
682 | /*9*/ vTwo4diffErr -= 1.5 * (2.*Avg2*Avg2 - Avg4)*(2.*Avg2*Avg2p - Avg4p) | |
683 | * W4W4p * Cov44p / (W4 * W4p); | |
684 | /*10*/vTwo4diffErr -= 4.*Avg2*TMath::Power(2.*Avg2*Avg2 - Avg4, 2) | |
685 | * W2pW4p * Cov2p2np / (W2p * W4p); | |
686 | /*11*/vTwo4diffErr /= TMath::Power(2.*Avg2*Avg2 - Avg4, 3.5); | |
687 | vTwo4diffErr = TMath::Sqrt(vTwo4diffErr); | |
d2bea14e | 688 | |
689 | cumulant4Hist->SetBinError(eta, vTwo4diffErr); | |
690 | } // End of eta loop | |
691 | ||
692 | // Number of events: | |
a0a945d3 | 693 | Int_t nEv = (Int_t)cumulantsHist->GetBinContent(0,0); |
d2bea14e | 694 | cumulant2Hist->SetBinContent(cumulant2Hist->GetNbinsX() + 1, nEv); |
695 | cumulant4Hist->SetBinContent(cumulant4Hist->GetNbinsX() + 1, nEv); | |
696 | } // End of harmonics loop | |
697 | } // End of type loop | |
698 | } | |
699 | // _____________________________________________________________________ | |
700 | void AliForwardFlowTaskQC::ProcessPrimary() | |
701 | { | |
702 | // | |
703 | // If fMC == kTRUE this function takes care of organizing the input | |
704 | // Monte Carlo data and histograms so AliForwardFlowTaskQC::QCumulants | |
705 | // can be run on it. | |
706 | // | |
707 | ||
708 | // Histograms are loaded and reset | |
709 | TH1D* dNdphi = (TH1D*)fOutputList->FindObject("hdNdphiSEMC"); | |
710 | TH2D* dNdetadphi = (TH2D*)fOutputList->FindObject("hdNdetadphiSEMC"); | |
711 | TH1D* dNdphiTrRef = (TH1D*)fOutputList->FindObject("hdNdphiSETrRef"); | |
712 | TH2D* dNdetadphiTrRef = (TH2D*)fOutputList->FindObject("hdNdetadphiSETrRef"); | |
713 | ||
714 | dNdphi->Reset(); | |
715 | dNdetadphi->Reset(); | |
716 | dNdphiTrRef->Reset(); | |
717 | dNdetadphiTrRef->Reset(); | |
718 | ||
719 | // Loads AliFMDFlowCommon and fills histograms and runs analysis. | |
720 | // AOD events also get a TrackRef histogram | |
721 | AliForwardFlowBase* common = new AliForwardFlowBase(fOutputList); | |
722 | ||
723 | if (fAOD) { | |
724 | if (!common->LoopAODmc(fAOD)) return; | |
725 | if (!common->LoopAODtrrefHits(fAOD)) return; | |
9d05ffeb | 726 | // if (!common->LoopMCaddptFlow(fAOD)) return; |
727 | // if (!common->LoopMCaddpdgFlow(fAOD)) return; | |
728 | // if (!common->LoopMCaddetaFlow(fAOD)) return; | |
d2bea14e | 729 | } |
730 | ||
731 | // Run analysis on MC truth | |
732 | for (Int_t n = 1; n <= 4; n++) { | |
733 | if (fv[n]) | |
734 | CumulantsMethod("MC", n); | |
735 | } | |
736 | ||
737 | // Run analysis on TrackRefs | |
738 | for (Int_t n = 1; n <= 4; n++) { | |
739 | if (fv[n]) | |
740 | CumulantsMethod("TrRef", n); | |
741 | } | |
742 | ||
743 | } | |
744 | //_____________________________________________________________________ | |
9d05ffeb | 745 | Double_t AliForwardFlowTaskQC::VarSQ(Double_t wxx2, Double_t x, Double_t wx, |
746 | Double_t wxx, Double_t sqrtwx2) | |
d2bea14e | 747 | { |
748 | // | |
749 | // Small function to compute the variance squared - used by Terminte() | |
750 | // | |
751 | Double_t sx; | |
752 | ||
753 | sx = wxx2 + x*x*wx - (Double_t)2*x*wxx; | |
754 | sx *= (Double_t)1 / wx; | |
755 | sx *= (Double_t)1 / (1 - TMath::Power(sqrtwx2, 2) / (wx*wx)); | |
756 | ||
757 | return sx; | |
758 | } | |
759 | //_____________________________________________________________________ | |
9d05ffeb | 760 | Double_t AliForwardFlowTaskQC::CovXY(Double_t wxwyxy, Double_t wxwy, |
761 | Double_t xy, Double_t wx, Double_t wy) | |
d2bea14e | 762 | { |
763 | // | |
764 | // Small function to compute the covariance between two numbers | |
765 | // - used by Terminate() | |
766 | // | |
767 | Double_t Cov, denominator, numerator; | |
768 | ||
769 | denominator = (wxwyxy / wxwy) - xy; | |
770 | numerator = 1 - (wxwy / (wx * wy)); | |
771 | ||
772 | Cov = denominator / numerator; | |
773 | return Cov; | |
774 | } | |
775 | //_____________________________________________________________________ | |
776 | // | |
777 | // | |
778 | // EOF |