PWGLF/FORWARD/analysis2/AliForwardFlowTaskQC.cxx

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