]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/FORWARD/analysis2/AliForwardFlowTaskQC.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Mega-commit by Alexander - added Event plane to analysis task, and other flow updates...
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliForwardFlowTaskQC.cxx
CommitLineData
d2bea14e 1//
d226802c 2// Calculate flow in the forward and central regions using the Q cumulants method.
d2bea14e 3//
4// Inputs:
5// - AliAODEvent
6//
7// Outputs:
8// - AnalysisResults.root
9//
d2bea14e 10#include <TROOT.h>
11#include <TSystem.h>
12#include <TInterpreter.h>
13#include <TChain.h>
14#include <TFile.h>
15#include <TList.h>
16#include <iostream>
17#include <TMath.h>
58f5fae2 18#include "TH3D.h"
2f9be372 19#include "TProfile2D.h"
d2bea14e 20#include "AliLog.h"
21#include "AliForwardFlowTaskQC.h"
22#include "AliAnalysisManager.h"
23#include "AliAODHandler.h"
24#include "AliAODInputHandler.h"
d2bea14e 25#include "AliAODForwardMult.h"
2b556440 26#include "AliAODCentralMult.h"
58f5fae2 27#include "AliAODEvent.h"
28
d2bea14e 29ClassImp(AliForwardFlowTaskQC)
9d05ffeb 30#if 0
31; // For emacs
32#endif
d2bea14e 33
34AliForwardFlowTaskQC::AliForwardFlowTaskQC()
2b556440 35 : AliAnalysisTaskSE(),
36 fBinsFMD(), // List with FMD flow histos
37 fBinsSPD(), // List with SPD flow histos
38 fSumList(0), // Event sum list
39 fOutputList(0), // Result output list
9d05ffeb 40 fAOD(0), // AOD input event
2b556440 41 fZvertex(1111), // Z vertex coordinate
42 fCent(-1), // Centrality
43 fHistCent(), // Histo for centrality
44 fHistVertexSel(), // Histo for selected vertices
45 fHistVertexAll() // Histo for all vertices
d2bea14e 46{
47 //
48 // Default constructor
49 //
d015ecfe 50 for (Int_t n = 0; n <= 6; n++) fv[n] = kTRUE;
d2bea14e 51}
52//_____________________________________________________________________
2b556440 53AliForwardFlowTaskQC::AliForwardFlowTaskQC(const char* name)
54 : AliAnalysisTaskSE(name),
55 fBinsFMD(), // List with FMD flow histos
56 fBinsSPD(), // List with SPD flow histos
57 fSumList(0), // Event sum list
58 fOutputList(0), // Result output list
59 fAOD(0), // AOD input event
60 fZvertex(1111), // Z vertex coordinate
61 fCent(-1), // Centrality
62 fHistCent(), // Histo for centrality
63 fHistVertexSel(), // Histo for selected vertices
64 fHistVertexAll() // Histo for all vertices
d2bea14e 65{
66 //
67 // Constructor
68 //
69 // Parameters:
70 // name: Name of task
71 //
d226802c 72 for (Int_t n = 0; n <= 6; n++) fv[n] = kTRUE;
2b556440 73
74 fHistCent = new TH1D("hCent", "Centralities", 100, 0, 100);
75 fHistVertexSel = new TH1D("hVertexSel", "Selected vertices", 40, -20, 20);
76 fHistVertexAll = new TH1D("hVertexAll", "All vertices", 40, -20, 20);
77
d2bea14e 78 DefineOutput(1, TList::Class());
2b556440 79 DefineOutput(2, TList::Class());
d2bea14e 80}
81//_____________________________________________________________________
2b556440 82AliForwardFlowTaskQC::AliForwardFlowTaskQC(const AliForwardFlowTaskQC& o)
83 : AliAnalysisTaskSE(o),
84 fBinsFMD(), // List with FMD flow histos
85 fBinsSPD(), // List with SPD flow histos
86 fSumList(o.fSumList), // Event sum list
87 fOutputList(o.fOutputList), // Result output list
88 fAOD(o.fAOD), // AOD input event
89 fZvertex(o.fZvertex), // Z vertex coordinate
90 fCent(o.fCent), // Centrality
91 fHistCent(o.fHistCent), // Histo for centrality
92 fHistVertexSel(o.fHistVertexSel), // Histo for selected vertices
93 fHistVertexAll(o.fHistVertexAll) // Histo for all vertices
d2bea14e 94{
95 //
96 // Copy constructor
97 //
98 // Parameters:
99 // o Object to copy from
100 //
d226802c 101 for (Int_t n = 0; n <= 6; n++) fv[n] = o.fv[n];
2b556440 102}
103//_____________________________________________________________________
104AliForwardFlowTaskQC&
105AliForwardFlowTaskQC::operator=(const AliForwardFlowTaskQC& o)
106{
107 //
108 // Assignment operator
109 //
110 fSumList = o.fSumList;
111 fOutputList = o.fOutputList;
112 fAOD = o.fAOD;
113 fZvertex = o.fZvertex;
114 fCent = o.fCent;
115 fHistCent = o.fHistCent;
116 fHistVertexSel = o.fHistVertexSel;
117 fHistVertexAll = o.fHistVertexAll;
118 fHistVertexAll = o.fHistVertexAll;
119
120 for (Int_t n = 0; n <= 6; n++) fv[n] = o.fv[n];
121 return *this;
d2bea14e 122}
123//_____________________________________________________________________
9453b19e 124void AliForwardFlowTaskQC::UserCreateOutputObjects()
d2bea14e 125{
126 //
127 // Create output objects
2b556440 128 //
129 InitVertexBins();
130 InitHists();
2f9be372 131
2b556440 132 PostData(1, fSumList);
133 PostData(2, fOutputList);
134
135}
136//_____________________________________________________________________
137void AliForwardFlowTaskQC::InitVertexBins()
138{
139 //
140 // Init vertexbin objects for FMD and SPD, and add them to the lists
141 //
142 for(Int_t n = 1; n <= 6; n++) {
143 if (!fv[n]) continue;
144 for (Int_t v = -10; v < 10; v++) {
145 fBinsFMD.Add(new VertexBin(v, v+1, n, "FMD"));
146 fBinsSPD.Add(new VertexBin(v, v+1, n, "SPD"));
d226802c 147 }
148 }
2b556440 149
150}
151//_____________________________________________________________________
152void AliForwardFlowTaskQC::InitHists()
153{
154 //
155 // Init histograms and add vertex bin histograms to the sum list
156 //
157 if (!fSumList)
158 fSumList = new TList();
159 fSumList->SetName("Sums");
160 fSumList->SetOwner();
161
162 TList* dList = new TList();
163 dList->SetName("Diagnostics");
164 dList->Add(fHistCent);
165 dList->Add(fHistVertexSel);
166 dList->Add(fHistVertexAll);
167 fSumList->Add(dList);
168
169 TIter nextFMD(&fBinsFMD);
170 VertexBin* bin = 0;
171 while ((bin = static_cast<VertexBin*>(nextFMD()))) {
172 bin->AddOutput(fSumList);
d226802c 173 }
2b556440 174 TIter nextSPD(&fBinsSPD);
175 while ((bin = static_cast<VertexBin*>(nextSPD()))) {
176 bin->AddOutput(fSumList);
d226802c 177 }
d2bea14e 178}
179//_____________________________________________________________________
180void AliForwardFlowTaskQC::UserExec(Option_t */*option*/)
181{
2b556440 182 //
183 // Calls the analyze function - called every event
d2bea14e 184 //
185 // Parameters:
186 // option: Not used
187 //
188
2b556440 189 Analyze();
190
191 PostData(1, fSumList);
192
193}
194//_____________________________________________________________________
195Bool_t AliForwardFlowTaskQC::Analyze()
196{
197 //
198 // Load FMD and SPD objects from aod tree and call the cumulants
199 // calculation for the correct vertexbin
200 //
201
d226802c 202 // Reset data members
203 fCent = -1;
204 fZvertex = 1111;
205
d2bea14e 206 // Get input event
207 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
2b556440 208 if (!fAOD) return kFALSE;
d2bea14e 209
2b556440 210 AliAODForwardMult* aodfmult = static_cast<AliAODForwardMult*>(fAOD->FindListObject("Forward"));
211 if (!aodfmult) return kFALSE;
212 if (!AODCheck(aodfmult)) return kFALSE;
213 TH2D fmddNdetadphi = aodfmult->GetHistogram();
d2bea14e 214
2b556440 215 AliAODCentralMult* aodcmult = static_cast<AliAODCentralMult*>(fAOD->FindListObject("CentralClusters"));
216 if (!aodcmult) return kFALSE;
217 TH2D spddNdetadphi = aodcmult->GetHistogram();
d226802c 218
2b556440 219 // TODO: remove me!
220// fCent = 0.5;
2f9be372 221
d226802c 222 // Run analysis on FMD and SPD
2b556440 223 TIter nextFMD(&fBinsFMD);
224 VertexBin* bin = 0;
225 while ((bin = static_cast<VertexBin*>(nextFMD()))) {
226 if (bin->CheckVertex(fZvertex)) {
227 if (!bin->FillHists(fmddNdetadphi)) return kFALSE;
228 bin->CumulantsAccumulate(fCent);
d226802c 229 }
2f9be372 230 }
d226802c 231
2b556440 232 TIter nextSPD(&fBinsSPD);
233 while ((bin = static_cast<VertexBin*>(nextSPD()))) {
234 if (bin->CheckVertex(fZvertex)) {
235 if (!bin->FillHists(spddNdetadphi)) return kFALSE;
236 bin->CumulantsAccumulate(fCent);
237 }
d226802c 238 }
d2bea14e 239
2b556440 240 return kTRUE;
d2bea14e 241}
242//_____________________________________________________________________
2b556440 243void AliForwardFlowTaskQC::Terminate(Option_t */*option*/)
d2bea14e 244{
2b556440 245 //
246 // Calls the finalize function, done at the end of the analysis
d2bea14e 247 //
248 // Parameters:
2b556440 249 // option: Not used
d2bea14e 250 //
2b556440 251
252 // Make sure pointers are set to the correct lists
253 fSumList = dynamic_cast<TList*> (GetOutputData(1));
254 if(!fSumList) {
255 AliError("Could not retrieve TList fSumList");
256 return;
257 }
258
259 if (!fOutputList)
260 fOutputList = new TList();
261 fOutputList->SetName("Results");
262 fOutputList->SetOwner();
263
264 // Run finalize on VertexBins
265 Finalize();
266
267 // Collect centralities
268 TProfile2D* hist2D = 0;
269 TList* centList = 0;
270 TH1D* hist1D = 0;
271 TIter nextProfile(fOutputList);
272 while ((hist2D = dynamic_cast<TProfile2D*>(nextProfile()))) {
273 for (Int_t cBin = 1; cBin <= 100; ) {
274 Int_t cMin = cBin - 1;
275 Int_t cMax = (cMin < 80 ? (cMin < 20 ? cMin + 5 : cMin + 10) : cMin + 20);
276 TString name = Form("cent_%d-%d", cMin, cMax);
277 centList = (TList*)fOutputList->FindObject(name.Data());
278 if (!centList) {
279 centList = new TList();
280 centList->SetName(name.Data());
281 fOutputList->Add(centList);
282 }
283 hist1D = hist2D->ProjectionX(Form("%s_%s", hist2D->GetName(), name.Data()),
284 cMin, cMax, "E");
285 hist1D->SetTitle(hist1D->GetName());
286 hist1D->Scale(1./(cMax-cMin));
287 centList->Add(hist1D);
288
289 cBin = cMax+1;
290 }
291 }
292
293 PostData(2, fOutputList);
294
295}
296//_____________________________________________________________________
297void AliForwardFlowTaskQC::Finalize()
298{
299 //
300 // Finalize command, called by Terminate()
301 //
302
303
304 // Reinitiate vertex bins if Terminate is called separately!
305 if (fBinsFMD.GetEntries() == 0) InitVertexBins();
306
307 // Iterate over all vertex bins objects and finalize cumulants
308 // calculations
309 TIter nextFMD(&fBinsFMD);
310 VertexBin* bin = 0;
311 while ((bin = static_cast<VertexBin*>(nextFMD()))) {
312 bin->CumulantsTerminate(fSumList, fOutputList);
313 }
314 TIter nextSPD(&fBinsSPD);
315 while ((bin = static_cast<VertexBin*>(nextSPD()))) {
316 bin->CumulantsTerminate(fSumList, fOutputList);
317 }
318
319}
320// _____________________________________________________________________
321Bool_t AliForwardFlowTaskQC::AODCheck(const AliAODForwardMult* aodfm)
322{
323 //
324 // Function to check that and AOD event meets the cuts
325 //
326 // Parameters:
327 // AliAODForwardMult: forward mult object with trigger and vertex info
328 //
329 // Returns false if there is no trigger or if the centrality or vertex
330 // is out of range. Otherwise true.
331 //
332
333 // First check for trigger
334 if (!aodfm->IsTriggerBits(AliAODForwardMult::kOffline)) return kFALSE;
335
336 // Then check for centrality
337 fCent = (Double_t)aodfm->GetCentrality();
338 if (0. >= fCent || fCent >= 100.) return kFALSE;
339 fHistCent->Fill(fCent);
340
341 // And finally check for vertex
342 fZvertex = aodfm->GetIpZ();
343 fHistVertexAll->Fill(fZvertex);
344 if (TMath::Abs(fZvertex) >= 10.) return kFALSE;
345 fHistVertexSel->Fill(fZvertex);
346
347 return kTRUE;
348
349}
350//_____________________________________________________________________
351AliForwardFlowTaskQC::VertexBin::VertexBin()
352 : TNamed(),
353 fMoment(0), // Flow moment for this vertexbin
354 fVzMin(0), // Vertex z-coordinate min
355 fVzMax(0), // Vertex z-coordinate max
356 fType(), // Data type (FMD/SPD/FMDTR/SPDTR/MC)
357 fCumuRef(), // Histogram for reference flow
358 fCumuDiff(), // Histogram for differential flow
359 fCumuHist(), // Sum histogram for cumulants
360 fHistTwoCorr(), // Diagnostics histogram for <<2>>
361 fHistW2(), // Diagnostics histogram for w_<<2>>
362 fHistFourCorr(), // Diagnostics histogram for <<4>>
363 fHistW4(), // Diagnostics histogram for w_<<4>>
364 fdNdedpAcc() // Diagnostics histogram to make acc. maps
365{
366 //
367 // Default constructor
368 //
369}
370//_____________________________________________________________________
371AliForwardFlowTaskQC::VertexBin::VertexBin(const Int_t vLow, const Int_t vHigh,
372 const Int_t moment, const Char_t* name)
373 : TNamed("", ""),
374 fMoment(moment), // Flow moment for this vertexbin
375 fVzMin(vLow), // Vertex z-coordinate min
376 fVzMax(vHigh), // Vertex z-coordinate max
377 fType(name), // Data type (FMD/SPD/FMDTR/SPDTR/MC)
378 fCumuRef(), // Histogram for reference flow
379 fCumuDiff(), // Histogram for differential flow
380 fCumuHist(), // Sum histogram for cumulants
381 fHistTwoCorr(), // Diagnostics histogram for <<2>>
382 fHistW2(), // Diagnostics histogram for w_<<2>>
383 fHistFourCorr(), // Diagnostics histogram for <<4>>
384 fHistW4(), // Diagnostics histogram for w_<<4>>
385 fdNdedpAcc() // Diagnostics histogram to make acc. maps
386{
387 //
388 // Constructor
389 //
390 // Parameters
391 // vLow: min z-coordinate
392 // vHigh: max z-coordinate
393 // moment: flow moment
394 // name: data type name (FMD/SPD/FMDTR/SPDTR/MC)
395 //
396 SetName(Form("%svertexBin%d_%d_%d", name, moment, vLow, vHigh));
397 SetTitle(Form("%svertexBin%d_%d_%d", name, moment, vLow, vHigh));
398
399}
400//_____________________________________________________________________
401AliForwardFlowTaskQC::VertexBin&
402AliForwardFlowTaskQC::VertexBin::operator=(const AliForwardFlowTaskQC::VertexBin& o)
403{
404 //
405 // Assignment operator
406 //
407 // Parameters
408 // o: AliForwardFlowTaskQC::VertexBin
409 //
410 fType = o.fType;
411 fCumuRef = o.fCumuRef;
412 fCumuDiff = o.fCumuDiff;
413 fCumuHist = o.fCumuHist;
414 fHistTwoCorr = o.fHistTwoCorr;
415 fHistW2 = o.fHistW2;
416 fHistFourCorr = o.fHistFourCorr;
417 fHistW4 = o.fHistW4;
418 fdNdedpAcc = o.fdNdedpAcc;
419
420 return *this;
421}
422//_____________________________________________________________________
423void AliForwardFlowTaskQC::VertexBin::AddOutput(TList* outputlist)
424{
425 //
426 // Add histograms to outputlist
427 //
428 // Parameters
429 // outputlist: list of histograms
430 //
431
432 // First we try to find an outputlist for this vertexbin
433 TList* list = (TList*)outputlist->FindObject(Form("%svertex_%d_%d", fType, fVzMin, fVzMax));
434
435 // If it doesn't exist we make one
436 if (!list) {
437 list = new TList();
438 list->SetName(Form("%svertex_%d_%d", fType, fVzMin, fVzMax));
439 outputlist->Add(list);
440 }
441
442 // We initiate the reference histogram according to an acceptance correction map,
443 // so we don't shift the SPD coverage within a reference bin
444 fCumuRef = new TH2D(Form("%s_v%d_%d_%d_ref", fType, fMoment, fVzMin, fVzMax),
445 Form("%s_v%d_%d_%d_ref", fType, fMoment, fVzMin, fVzMax),
446 24, -6., 6., 5, 0.5, 5.5);
447 TFile acc("$ALICE_ROOT/PWGLF/FORWARD/corrections/FlowCorrections/FlowAccMap.root", "READ");
448 TH1D* accMap = (TH1D*)acc.Get(Form("%saccVertex_%d_%d", fType, fVzMin, fVzMax));
449 if (accMap) {
450 Int_t nBins = accMap->GetNbinsX();
451 Double_t eta[48] = { 0. };
452 Int_t n = 0;
453 Double_t newOcc[48] = { 0. };
454 Double_t prev = -1;
455 for (Int_t i = 0; i < nBins; i++) {
456 Double_t occ = accMap->GetBinContent(i+1);
457 if (prev != occ && (occ > 0.6 || occ == 0)) {
458 eta[n] = i*0.25-6.;
459 newOcc[n] = occ;
460 n++;
461// printf("eta: %f \t occ: %f \t Vertex: %d \n", eta[n-1], occ, fVzMin);
462 }
463 prev = occ;
464 }
465 eta[n] = 6.;
466 fCumuRef->GetXaxis()->Set(n, eta);
467 }
468 acc.Close();
469
470 fCumuRef->Sumw2();
471 list->Add(fCumuRef);
472
473 // We initiate the differential histogram
474 fCumuDiff = new TH2D(Form("%s_v%d_%d_%d_diff", fType, fMoment, fVzMin, fVzMax),
475 Form("%s_v%d_%d_%d_diff", fType, fMoment, fVzMin, fVzMax),
476 48, -6., 6., 5, 0.5, 5.5);
477 fCumuDiff->Sumw2();
478 list->Add(fCumuDiff);
479
480 // Initiate the cumulant sum histogram
481 fCumuHist = new TH3D(Form("%sv%d_vertex_%d_%d", fType, fMoment, fVzMin, fVzMax),
482 Form("%sv%d_vertex_%d_%d", fType, fMoment, fVzMin, fVzMax),
483 48, -6., 6., 100, 0., 100., 26, 0.5, 26.5);
484 fCumuHist->Sumw2();
485
486 list->Add(fCumuHist);
487
488 // We check for diagnostics histograms (only done per type and moment, not vertexbin)
489 // If they are not found we create them.
490 TList* dList = (TList*)outputlist->FindObject("Diagnostics");
491 if (!dList) AliFatal("No diagnostics list found, what kind of game are you running here?!?!");
492
493 // Corr. hists are shared over all vertex bins...
494 fHistTwoCorr = (TH3D*)dList->FindObject(Form("hHistTwoCorr_%s_v%d", fType, fMoment));
495 if (fHistTwoCorr) {
496 fHistW2 = (TH3D*)dList->FindObject(Form("hHistW2_%s_v%d", fType, fMoment));
497 fHistFourCorr = (TH3D*)dList->FindObject(Form("hHistFourCorr_%s_v%d", fType, fMoment));
498 fHistW4 = (TH3D*)dList->FindObject(Form("hHistW4_%s_v%d", fType, fMoment));
499 } else {
500 fHistTwoCorr = new TH3D(Form("hHistTwoCorr_%s_v%d", fType, fMoment),
501 "Two particle correlator: w_{2}<<2>>",
502 48, -6., 6., 100, 0., 150000, 100, 0., 100.);
503 fHistTwoCorr->Sumw2();
504 fHistW2 = new TH3D(Form("hHistW2_%s_v%d", fType, fMoment),
505 "Two particle event weight: w_{2}",
506 48, -6., 6., 100, 0., 2e+7, 100, 0., 100.);
507 fHistW2->Sumw2();
508 fHistFourCorr = new TH3D(Form("hHistFourCorr_%s_v%d", fType, fMoment),
509 "Four particle correlator: w_{4}<<4>>",
510 48, -6., 6., 100, 0., 1e+10, 100, 0., 100.);
511 fHistFourCorr->Sumw2();
512 fHistW4 = new TH3D(Form("hHistW4_%s_v%d", fType, fMoment),
513 "Four particle event weight: w_{4}",
514 48, -6., 6., 100, 0., 4e+14, 100, 0., 100.);
515 fHistW4->Sumw2();
516
517 dList->Add(fHistTwoCorr);
518 dList->Add(fHistW2);
519 dList->Add(fHistFourCorr);
520 dList->Add(fHistW4);
521 }
2f9be372 522
2b556440 523 // Acceptance hists are shared over all moments
524 fdNdedpAcc = (TH2D*)dList->FindObject(Form("h%sdNdedpAcc_%d_%d", fType, fVzMin, fVzMax));
525 if (!fdNdedpAcc) {
526 fdNdedpAcc = new TH2D(Form("h%sdNdedpAcc_%d_%d", fType, fVzMin, fVzMax),
527 Form("%s acceptance map for %d cm < v_{z} < %d cm", fType, fVzMin, fVzMax),
528 48, -6, 6, 20, 0, TMath::TwoPi());
529 fdNdedpAcc->Sumw2();
530 dList->Add(fdNdedpAcc);
531 }
532
533
534}
535//_____________________________________________________________________
536Bool_t AliForwardFlowTaskQC::VertexBin::CheckVertex(Double_t vz)
537{
538 //
539 // We check if this is the correct bin for the current event's vertex
540 //
541 // Parameters:
542 // vZ: Current event vertex
543 //
544 // Returns false if out of range, true otherwise
545 //
546 if ((Double_t)fVzMin > vz) return kFALSE;
547 if ((Double_t)fVzMax <= vz) return kFALSE;
548
549 return kTRUE;
550}
551//_____________________________________________________________________
552Bool_t AliForwardFlowTaskQC::VertexBin::FillHists(TH2D dNdetadphi)
553{
554 //
555 // Fill reference and differential eta-histograms
556 //
557 // Parameters:
558 // dNdetadphi: 2D histogram with input data
559 //
560
561 if (!fCumuRef) AliFatal("You have not called AddOutput() - Terminating!");
562
563 // Fist we reset histograms
564 fCumuRef->Reset();
565 fCumuDiff->Reset();
566
567 // Numbers to cut away bad events and acceptance.
568 Double_t runAvg = 0;
569 Double_t max = 0;
570 Int_t nInAvg = 0;
571 Int_t nBadBins = 0;
572 Int_t nBins = (dNdetadphi.GetNbinsX() * 6) / (fCumuDiff->GetNbinsX() * 5);
573 Int_t nInBin = 0;
574 Int_t nCurBin = 0, nPrevBin = 0;
575
576 // Then we loop over the input and calculate sum cos(k*n*phi)
577 // and fill it in the reference and differential histograms
578 Double_t eta, phi, weight;
579 Double_t dQnRe = 0, dQ2nRe = 0, dQnIm = 0, dQ2nIm = 0;
580 for (Int_t etaBin = 1; etaBin <= dNdetadphi.GetNbinsX(); etaBin++) {
581 eta = dNdetadphi.GetXaxis()->GetBinCenter(etaBin);
582 nCurBin = fCumuDiff->GetXaxis()->FindBin(eta);
583 // If we have moved to a new bin in the flow hist, and less than half the eta
584 // region has been covered by it we cut it away.
585 if (nPrevBin && nCurBin != nPrevBin) {
586 if (nInBin <= nBins/2) {
587 for (Int_t pBin = 1; pBin <= fCumuDiff->GetNbinsY(); pBin++) {
588 fCumuDiff->SetBinContent(nPrevBin, pBin, 0);
589 fCumuDiff->SetBinError(nPrevBin, pBin, 0);
590 }
591 }
592 nInBin = 0;
593 nPrevBin = nCurBin;
594 runAvg = 0;
595 nInAvg = 0;
596 max = 0;
597 }
598 Bool_t data = kFALSE;
599 for (Int_t phiBin = 1; phiBin <= dNdetadphi.GetNbinsY(); phiBin++) {
600 phi = dNdetadphi.GetYaxis()->GetBinCenter(phiBin);
601 weight = dNdetadphi.GetBinContent(etaBin, phiBin);
602 if (!weight) continue;
603 if (!data) data = kTRUE;
604 // We calculate the running average Nch per. bin
605 runAvg *= nInAvg;
606 runAvg += weight;
607 nInAvg++;
608 runAvg /= nInAvg;
609 // if the current bin has more than write the avg we count a bad bin
610 if (weight > max) max = weight;
611
612 dQnRe = weight*TMath::Cos(fMoment*phi);
613 dQnIm = weight*TMath::Sin(fMoment*phi);
614 dQ2nRe = weight*TMath::Cos(2.*fMoment*phi);
615 dQ2nIm = weight*TMath::Sin(2.*fMoment*phi);
616
617 fCumuRef->Fill(eta, kHmult, weight);
618 fCumuRef->Fill(eta, kHQnRe, dQnRe);
619 fCumuRef->Fill(eta, kHQnIm, dQnIm);
620 fCumuRef->Fill(eta, kHQ2nRe, dQ2nRe);
621 fCumuRef->Fill(eta, kHQ2nIm, dQ2nIm);
622
623 fCumuDiff->Fill(eta, kHmult, weight);
624 fCumuDiff->Fill(eta, kHQnRe, dQnRe);
625 fCumuDiff->Fill(eta, kHQnIm, dQnIm);
626 fCumuDiff->Fill(eta, kHQ2nRe, dQ2nRe);
627 fCumuDiff->Fill(eta, kHQ2nIm, dQ2nIm);
628
629 // Fill acc. map
630 fdNdedpAcc->Fill(eta, phi, weight);
631 }
632 if (data) {
633 nInBin++;
634 if (max > 35) nBadBins++;
635 }
636 // If there are too many bad bins we throw the event away!
637 if (nBadBins > 2) return kFALSE;
638 }
639 return kTRUE;
640}
641//_____________________________________________________________________
642void AliForwardFlowTaskQC::VertexBin::CumulantsAccumulate(Double_t cent)
643{
644 //
645 // Calculate the Q cumulant of order fMoment
646 //
647 // Parameters:
648 // cent: Centrality of event
649 //
650
651 if (!fCumuRef) AliFatal("You have not called AddOutput() - Terminating!");
652
d2bea14e 653 // We create the objects needed for the analysis
d226802c 654 Double_t dQnRe = 0, dQ2nRe = 0, dQnIm = 0, dQ2nIm = 0, mult = 0;
58f5fae2 655 Double_t pnRe = 0, p2nRe = 0, qnRe = 0, q2nRe = 0, pnIm = 0, p2nIm = 0, qnIm = 0, q2nIm = 0;
2f9be372 656 Double_t two = 0, four = 0, twoPrime = 0, fourPrime = 0;
d226802c 657 Double_t cosPhi1Phi2 = 0, cosPhi1Phi2Phi3m = 0;
658 Double_t sinPhi1Phi2 = 0, sinPhi1Phi2Phi3m = 0;
659 Double_t cosPsi1Phi2 = 0, cosPsi1Phi2Phi3m = 0, cosPsi1Phi2Phi3p = 0;
660 Double_t sinPsi1Phi2 = 0, sinPsi1Phi2Phi3m = 0, sinPsi1Phi2Phi3p = 0;
2b556440 661 Double_t eta = 0;
d2bea14e 662 Double_t multi = 0, multp = 0, mp = 0, mq = 0;
58f5fae2 663 Double_t w2 = 0, w4 = 0, w2p = 0, w4p = 0;
d226802c 664 Int_t refEtaBin = 0;
2b556440 665 Bool_t eventCount = kFALSE;
666
d2bea14e 667 // We loop over the data 1 time!
2b556440 668 for (Int_t etaBin = 1; etaBin <= fCumuDiff->GetNbinsX(); etaBin++) {
669 eta = fCumuDiff->GetXaxis()->GetBinCenter(etaBin);
670 refEtaBin = fCumuRef->GetXaxis()->FindBin(eta);
58f5fae2 671 // The values for each individual etaBin bins are reset
d2bea14e 672 mp = 0;
673 pnRe = 0;
674 p2nRe = 0;
675 pnIm = 0;
676 p2nIm = 0;
677
d226802c 678 mult = 0;
679 dQnRe = 0;
680 dQnIm = 0;
681 dQ2nRe = 0;
682 dQ2nIm = 0;
d2bea14e 683
2b556440 684 // Reference flow
685 multi = fCumuRef->GetBinContent(refEtaBin, kHmult);
686 dQnRe = fCumuRef->GetBinContent(refEtaBin, kHQnRe);
687 dQnIm = fCumuRef->GetBinContent(refEtaBin, kHQnIm);
688 dQ2nRe = fCumuRef->GetBinContent(refEtaBin, kHQ2nRe);
689 dQ2nIm = fCumuRef->GetBinContent(refEtaBin, kHQ2nIm);
690 mult += multi;
691
692 // For each etaBin bin the necessary values for differential flow
693 // is calculated. Here is the loop over the phi's.
694 multp = fCumuDiff->GetBinContent(etaBin, kHmult);
695 pnRe = fCumuDiff->GetBinContent(etaBin, kHQnRe);
696 pnIm = fCumuDiff->GetBinContent(etaBin, kHQnIm);
697 p2nRe = fCumuDiff->GetBinContent(etaBin, kHQ2nRe);
698 p2nIm = fCumuDiff->GetBinContent(etaBin, kHQ2nIm);
699 mp += multp;
700
701 if (mult <= 3) continue;
d226802c 702
2b556440 703 if (!eventCount) {
d226802c 704 // Count number of events
2b556440 705 fCumuHist->Fill(-7., cent, -0.5, 1.);
706 eventCount = kTRUE;
d226802c 707 }
2b556440 708 if (mp == 0) continue;
d226802c 709 // The reference flow is calculated
d226802c 710
711 // 2-particle
712 w2 = mult * (mult - 1.);
713 two = dQnRe*dQnRe + dQnIm*dQnIm - mult;
2b556440 714
715 fCumuHist->Fill(eta, cent, kW2Two, two);
716 fCumuHist->Fill(eta, cent, kW2, w2);
717
718 fCumuHist->Fill(eta, cent, kQnRe, dQnRe);
719 fCumuHist->Fill(eta, cent, kQnIm, dQnIm);
720 fCumuHist->Fill(eta, cent, kM, mult);
d226802c 721
722 // 4-particle
723 w4 = mult * (mult - 1.) * (mult - 2.) * (mult - 3.);
724
725 four = 2.*mult*(mult-3.) + TMath::Power((TMath::Power(dQnRe,2.)+TMath::Power(dQnIm,2.)),2.)
726 -4.*(mult-2.)*(TMath::Power(dQnRe,2.) + TMath::Power(dQnIm,2.))
727 -2.*(TMath::Power(dQnRe,2.)*dQ2nRe+2.*dQnRe*dQnIm*dQ2nIm-TMath::Power(dQnIm,2.)*dQ2nRe)
728 +(TMath::Power(dQ2nRe,2.)+TMath::Power(dQ2nIm,2.));
2f9be372 729
2b556440 730 fCumuHist->Fill(eta, cent, kW4Four, four);
731 fCumuHist->Fill(eta, cent, kW4, w4);
2f9be372 732
d226802c 733 cosPhi1Phi2 = dQnRe*dQnRe - dQnIm*dQnIm - dQ2nRe;
734 sinPhi1Phi2 = 2.*dQnRe*dQnIm - dQ2nIm;
d2bea14e 735
d226802c 736 cosPhi1Phi2Phi3m = dQnRe*(TMath::Power(dQnRe,2)+TMath::Power(dQnIm,2))-dQnRe*dQ2nRe-dQnIm*dQ2nIm-2.*(mult-1)*dQnRe;
d2bea14e 737
d226802c 738 sinPhi1Phi2Phi3m = -dQnIm*(TMath::Power(dQnRe,2)+TMath::Power(dQnIm,2))+dQnRe*dQ2nIm-dQnIm*dQ2nRe+2.*(mult-1)*dQnIm;
d2bea14e 739
2b556440 740 fCumuHist->Fill(eta, cent, kCosphi1phi2, cosPhi1Phi2);
741 fCumuHist->Fill(eta, cent, kSinphi1phi2, sinPhi1Phi2);
742 fCumuHist->Fill(eta, cent, kCosphi1phi2phi3m, cosPhi1Phi2Phi3m);
743 fCumuHist->Fill(eta, cent, kSinphi1phi2phi3m, sinPhi1Phi2Phi3m);
744 fCumuHist->Fill(eta, cent, kMm1m2, mult*(mult-1.)*(mult-2.));
745
746 // Diagnostics are filled
747 fHistTwoCorr->Fill(eta, two, cent);
748 fHistW2->Fill(eta, w2, cent);
749 fHistFourCorr->Fill(eta, four, cent);
750 fHistW4->Fill(eta, w4, cent);
751
752 // Differential flow calculations for each eta bin bin is done:
d2bea14e 753 mq = mp;
754 qnRe = pnRe;
755 qnIm = pnIm;
58f5fae2 756 q2nRe = p2nRe;
757 q2nIm = p2nIm;
d2bea14e 758
759 // 2-particle differential flow
58f5fae2 760 w2p = mp * mult - mq;
2f9be372 761 twoPrime = pnRe*dQnRe + pnIm*dQnIm - mq;
d2bea14e 762
2b556440 763 fCumuHist->Fill(eta, cent, kw2two, twoPrime);
764 fCumuHist->Fill(eta, cent, kw2, w2p);
765
766 fCumuHist->Fill(eta, cent, kpnRe, pnRe);
767 fCumuHist->Fill(eta, cent, kpnIm, pnIm);
768 fCumuHist->Fill(eta, cent, kmp, mp);
d2bea14e 769
770 // 4-particle differential flow
58f5fae2 771 w4p = (mp * mult - 3.*mq)*(mult - 1.)*(mult - 2.);
d226802c 772
773 fourPrime = (TMath::Power(dQnRe,2.)+TMath::Power(dQnIm,2.))*(pnRe*dQnRe+pnIm*dQnIm)
774 - q2nRe*(TMath::Power(dQnRe,2.)-TMath::Power(dQnIm,2.))
775 - 2.*q2nIm*dQnRe*dQnIm
776 - pnRe*(dQnRe*dQ2nRe+dQnIm*dQ2nIm)
777 + pnIm*(dQnIm*dQ2nRe-dQnRe*dQ2nIm)
778 - 2.*mult*(pnRe*dQnRe+pnIm*dQnIm)
779 - 2.*(TMath::Power(dQnRe,2.)+TMath::Power(dQnIm,2.))*mq
780 + 6.*(qnRe*dQnRe+qnIm*dQnIm)
781 + 1.*(q2nRe*dQ2nRe+q2nIm*dQ2nIm)
782 + 2.*(pnRe*dQnRe+pnIm*dQnIm)
783 + 2.*mq*mult
784 - 6.*mq;
2b556440 785
786 fCumuHist->Fill(eta, cent, kw4four, fourPrime);
787 fCumuHist->Fill(eta, cent, kw4, w4p);
58f5fae2 788
d226802c 789 cosPsi1Phi2 = pnRe*dQnRe - pnIm*dQnIm - q2nRe;
790 sinPsi1Phi2 = pnRe*dQnIm + pnIm*dQnRe - q2nIm;
58f5fae2 791
d226802c 792 cosPsi1Phi2Phi3p = pnRe*(TMath::Power(dQnIm,2.)+TMath::Power(dQnRe,2.)-mult)
793 - 1.*(q2nRe*dQnRe+q2nIm*dQnIm)
794 - mq*dQnRe+2.*qnRe;
795
796 sinPsi1Phi2Phi3p = pnIm*(TMath::Power(dQnIm,2.)+TMath::Power(dQnRe,2.)-mult)
797 - 1.*(q2nIm*dQnRe-q2nRe*dQnIm)
798 - mq*dQnIm+2.*qnIm;
58f5fae2 799
d226802c 800 cosPsi1Phi2Phi3m = pnRe*(TMath::Power(dQnRe,2.)-TMath::Power(dQnIm,2.))+2.*pnIm*dQnRe*dQnIm
801 - 1.*(pnRe*dQ2nRe+pnIm*dQ2nIm)
802 - 2.*mq*dQnRe+2.*qnRe;
803
804 sinPsi1Phi2Phi3m = pnIm*(TMath::Power(dQnRe,2.)-TMath::Power(dQnIm,2.))-2.*pnRe*dQnRe*dQnIm
805 - 1.*(pnIm*dQ2nRe-pnRe*dQ2nIm)
806 + 2.*mq*dQnIm-2.*qnIm;
807
2b556440 808 fCumuHist->Fill(eta, cent, kCospsi1phi2, cosPsi1Phi2);
809 fCumuHist->Fill(eta, cent, kSinpsi1phi2, sinPsi1Phi2);
810 fCumuHist->Fill(eta, cent, kCospsi1phi2phi3m, cosPsi1Phi2Phi3m);
811 fCumuHist->Fill(eta, cent, kSinpsi1phi2phi3m, sinPsi1Phi2Phi3m);
812 fCumuHist->Fill(eta, cent, kmpmq, (mp*mult-2.*mq)*(mult-1.));
813 fCumuHist->Fill(eta, cent, kCospsi1phi2phi3p, cosPsi1Phi2Phi3p);
814 fCumuHist->Fill(eta, cent, kSinpsi1phi2phi3p, sinPsi1Phi2Phi3p);
d2bea14e 815 }
816
817}
818//_____________________________________________________________________
2b556440 819void AliForwardFlowTaskQC::VertexBin::CumulantsTerminate(TList* inlist, TList* outlist)
d2bea14e 820{
821 //
d2bea14e 822 // Finalizes the Q cumulant calculations
823 //
824 // Parameters:
2b556440 825 // inlist: input sumlist
826 // outlist: output result list
d2bea14e 827 //
2f9be372 828
2b556440 829 // Re-find cumulants hist if Terminate is called separately
830 if (!fCumuHist) {
831 TList* list = (TList*)inlist->FindObject(Form("%svertex_%d_%d", fType, fVzMin, fVzMax));
832 fCumuHist = (TH3D*)list->FindObject(Form("%sv%d_vertex_%d_%d", fType, fMoment, fVzMin, fVzMax));
833 }
834
835 // Create result profiles
836 TProfile2D* cumu2 = (TProfile2D*)outlist->FindObject(Form("%sQC2_v%d_unCorr", fType, fMoment));
837 TProfile2D* cumu4 = (TProfile2D*)outlist->FindObject(Form("%sQC4_v%d_unCorr", fType, fMoment));
838 if (!cumu2) {
839 cumu2 = new TProfile2D(Form("%sQC2_v%d_unCorr", fType, fMoment),
840 Form("%sQC2_v%d_unCorr", fType, fMoment),
841 48, -6., 6., 100, 0., 100);
842 outlist->Add(cumu2);
843 }
844 if (!cumu4) {
845 cumu4 = new TProfile2D(Form("%sQC4_v%d_unCorr", fType, fMoment),
846 Form("%sQC4_v%d_unCorr", fType, fMoment),
847 48, -6., 6., 100, 0., 100);
848 outlist->Add(cumu4);
849 }
58f5fae2 850
851 // For flow calculations
2b556440 852 Double_t two = 0, qc2 = 0, vnTwo = 0, four = 0, qc4 = 0, vnFour = 0;
2f9be372 853 Double_t twoPrime = 0, qc2Prime = 0, vnTwoDiff = 0, fourPrime = 0, qc4Prime = 0, vnFourDiff = 0;
854 Double_t w2 = 0, w4 = 0, w2p = 0, w4p = 0;
2f9be372 855 Double_t w2Two = 0, w2pTwoPrime = 0, w4Four = 0, w4pFourPrime = 0;
58f5fae2 856 Double_t cosP1nPhi = 0, sinP1nPhi = 0, mult = 0, cosP1nPhi1P1nPhi2 = 0, sinP1nPhi1P1nPhi2 = 0;
857 Double_t cosP1nPhi1M1nPhi2M1nPhi3 = 0, sinP1nPhi1M1nPhi2M1nPhi3 = 0, multm1m2 = 0;
858 Double_t cosP1nPsi = 0, sinP1nPsi = 0, mp = 0, cosP1nPsi1P1nPhi2 = 0, sinP1nPsi1P1nPhi2 = 0;
859 Double_t cosP1nPsi1M1nPhi2M1nPhi3 = 0, sinP1nPsi1M1nPhi2M1nPhi3 = 0, mpqMult = 0;
860 Double_t cosP1nPsi1P1nPhi2M1nPhi3 = 0, sinP1nPsi1P1nPhi2M1nPhi3 = 0;
2f9be372 861
2b556440 862 // Loop over cumulant histogram for final calculations
863 // Centrality loop
864 for (Int_t c = 1; c <= 80; c++) {
865 Double_t nEv = 0;
866 // Eta loop
867 for (Int_t etaBin = 1; etaBin <= fCumuHist->GetNbinsX(); etaBin++) {
868 Double_t eta = fCumuHist->GetXaxis()->GetBinCenter(etaBin);
869 // 2-particle reference flow
870 w2Two = fCumuHist->GetBinContent(etaBin, c, kW2Two);
871 w2 = fCumuHist->GetBinContent(etaBin, c, kW2);
872 mult = fCumuHist->GetBinContent(etaBin, c, kM);
873 if (!w2 || !mult) continue;
874 cosP1nPhi = fCumuHist->GetBinContent(etaBin, c, kQnRe);
875 sinP1nPhi = fCumuHist->GetBinContent(etaBin, c, kQnIm);
2f9be372 876
2b556440 877 cosP1nPhi /= mult;
878 sinP1nPhi /= mult;
879 two = w2Two / w2;
880 qc2 = two - TMath::Power(cosP1nPhi, 2) - TMath::Power(sinP1nPhi, 2);
881 if (qc2 <= 0) continue;
882 vnTwo = TMath::Sqrt(qc2);
883 // if (!TMath::IsNaN(vnTwo*mult))
884 // cumu2->Fill(eta, vnTwo, fCumuHist->GetBinContent(0,c,0));
885
886 // 4-particle reference flow
887 w4Four = fCumuHist->GetBinContent(etaBin, c, kW4Four);
888 w4 = fCumuHist->GetBinContent(etaBin, c, kW4);
889 multm1m2 = fCumuHist->GetBinContent(etaBin, c, kMm1m2);
890 if (!w4 || !multm1m2) continue;
891 cosP1nPhi1P1nPhi2 = fCumuHist->GetBinContent(etaBin, c, kCosphi1phi2);
892 sinP1nPhi1P1nPhi2 = fCumuHist->GetBinContent(etaBin, c, kSinphi1phi2);
893 cosP1nPhi1M1nPhi2M1nPhi3 = fCumuHist->GetBinContent(etaBin, c, kCosphi1phi2phi3m);
894 sinP1nPhi1M1nPhi2M1nPhi3 = fCumuHist->GetBinContent(etaBin, c, kSinphi1phi2phi3m);
895
896 cosP1nPhi1P1nPhi2 /= w2;
897 sinP1nPhi1P1nPhi2 /= w2;
898 cosP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
899 sinP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
900 four = w4Four / w4;
901 qc4 = four-2.*TMath::Power(two,2.)
902 - 4.*cosP1nPhi*cosP1nPhi1M1nPhi2M1nPhi3
903 + 4.*sinP1nPhi*sinP1nPhi1M1nPhi2M1nPhi3-TMath::Power(cosP1nPhi1P1nPhi2,2.)-TMath::Power(sinP1nPhi1P1nPhi2,2.)
904 + 4.*cosP1nPhi1P1nPhi2*(TMath::Power(cosP1nPhi,2.)-TMath::Power(sinP1nPhi,2.))
905 + 8.*sinP1nPhi1P1nPhi2*sinP1nPhi*cosP1nPhi
906 + 8.*two*(TMath::Power(cosP1nPhi,2.)+TMath::Power(sinP1nPhi,2.))
907 - 6.*TMath::Power((TMath::Power(cosP1nPhi,2.)+TMath::Power(sinP1nPhi,2.)),2.);
908
909 if (qc4 >= 0) continue;
910 vnFour = TMath::Power(-qc4, 0.25);
911 // if (!TMath::IsNaN(vnFour*mult))
912 // cumu4->Fill(eta, vnFour, fCumuHist->GetBinContent(0,c,0));
913
914 // 2-particle differential flow
915 w2pTwoPrime = fCumuHist->GetBinContent(etaBin, c, kw2two);
916 w2p = fCumuHist->GetBinContent(etaBin, c, kw2);
917 mp = fCumuHist->GetBinContent(etaBin, c, kmp);
918 if (!w2p || !mp) continue;
919 cosP1nPsi = fCumuHist->GetBinContent(etaBin, c, kpnRe);
920 sinP1nPsi = fCumuHist->GetBinContent(etaBin, c, kpnIm);
921
922 cosP1nPsi /= mp;
923 sinP1nPsi /= mp;
924 twoPrime = w2pTwoPrime / w2p;
925 qc2Prime = twoPrime - sinP1nPsi*sinP1nPhi - cosP1nPsi*cosP1nPhi;
926
927 vnTwoDiff = qc2Prime / TMath::Sqrt(qc2);
928 if (!TMath::IsNaN(vnTwoDiff*mp)) cumu2->Fill(eta, (Double_t)c-1., vnTwoDiff, fCumuHist->GetBinContent(0,c,0));
929
930 // 4-particle differential flow
931 w4pFourPrime = fCumuHist->GetBinContent(etaBin, c, kw4four);
932 w4p = fCumuHist->GetBinContent(etaBin, c, kw4);
933 mpqMult = fCumuHist->GetBinContent(etaBin, c, kmpmq);
934 if (!w4p || !mpqMult) continue;
935 cosP1nPsi1P1nPhi2 = fCumuHist->GetBinContent(etaBin, c, kCospsi1phi2);
936 sinP1nPsi1P1nPhi2 = fCumuHist->GetBinContent(etaBin, c, kSinpsi1phi2);
937 cosP1nPsi1M1nPhi2M1nPhi3 = fCumuHist->GetBinContent(etaBin, c, kCospsi1phi2phi3m);
938 sinP1nPsi1M1nPhi2M1nPhi3 = fCumuHist->GetBinContent(etaBin, c, kSinpsi1phi2phi3m);
939 cosP1nPsi1P1nPhi2M1nPhi3 = fCumuHist->GetBinContent(etaBin, c, kCospsi1phi2phi3p);
940 sinP1nPsi1P1nPhi2M1nPhi3 = fCumuHist->GetBinContent(etaBin, c, kSinpsi1phi2phi3p);
941
942 cosP1nPsi1P1nPhi2 /= w2p;
943 sinP1nPsi1P1nPhi2 /= w2p;
944 cosP1nPsi1M1nPhi2M1nPhi3 /= mpqMult;
945 sinP1nPsi1M1nPhi2M1nPhi3 /= mpqMult;
946 cosP1nPsi1P1nPhi2M1nPhi3 /= mpqMult;
947 sinP1nPsi1P1nPhi2M1nPhi3 /= mpqMult;
948
949 fourPrime = w4pFourPrime / w4p;
950
951 qc4Prime = fourPrime-2.*twoPrime*two
952 - cosP1nPsi*cosP1nPhi1M1nPhi2M1nPhi3
953 + sinP1nPsi*sinP1nPhi1M1nPhi2M1nPhi3
954 - cosP1nPhi*cosP1nPsi1M1nPhi2M1nPhi3
955 + sinP1nPhi*sinP1nPsi1M1nPhi2M1nPhi3
956 - 2.*cosP1nPhi*cosP1nPsi1P1nPhi2M1nPhi3
957 - 2.*sinP1nPhi*sinP1nPsi1P1nPhi2M1nPhi3
958 - cosP1nPsi1P1nPhi2*cosP1nPhi1P1nPhi2
959 - sinP1nPsi1P1nPhi2*sinP1nPhi1P1nPhi2
960 + 2.*cosP1nPhi1P1nPhi2*(cosP1nPsi*cosP1nPhi-sinP1nPsi*sinP1nPhi)
961 + 2.*sinP1nPhi1P1nPhi2*(cosP1nPsi*sinP1nPhi+sinP1nPsi*cosP1nPhi)
962 + 4.*two*(cosP1nPsi*cosP1nPhi+sinP1nPsi*sinP1nPhi)
963 + 2.*cosP1nPsi1P1nPhi2*(TMath::Power(cosP1nPhi,2.)-TMath::Power(sinP1nPhi,2.))
964 + 4.*sinP1nPsi1P1nPhi2*cosP1nPhi*sinP1nPhi
965 + 4.*twoPrime*(TMath::Power(cosP1nPhi,2.)+TMath::Power(sinP1nPhi,2.))
966 - 6.*(TMath::Power(cosP1nPhi,2.)-TMath::Power(sinP1nPhi,2.))
967 * (cosP1nPsi*cosP1nPhi-sinP1nPsi*sinP1nPhi)
968 - 12.*cosP1nPhi*sinP1nPhi
969 * (sinP1nPsi*cosP1nPhi+cosP1nPsi*sinP1nPhi);
970
971 vnFourDiff = - qc4Prime / TMath::Power(-qc4, 0.75);
972 if (!TMath::IsNaN(vnFourDiff*mp)) cumu4->Fill(eta, (Double_t)c-1., vnFourDiff, fCumuHist->GetBinContent(0,c,0));
973 } // End of eta loop
974 // Number of events:
975 nEv += fCumuHist->GetBinContent(0,c,0);
976 cumu2->Fill(7., (Double_t)c-1., nEv);
977 cumu4->Fill(7., (Double_t)c-1., nEv);
978 } // End of centrality loop
d226802c 979
980 return;
d2bea14e 981}
d2bea14e 982//_____________________________________________________________________
983//
984//
985// EOF