AliPhysics  ff1d528 (ff1d528)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskLocalRho.cxx
Go to the documentation of this file.
1 // $Id$
2 //
3 // Analysis task to estimate an event's local energy density
4 //
5 // This task is part of the emcal jet framework and should be run in the emcaljet train
6 // The following extensions to an accepted AliVEvent are expected:
7 // - (anti-kt) jets -> necessary if one wants to exclude leading jet contribution to the event plane
8 // - background estimate of rho -> this task estimates modulation, not rho itself
9 // - pico tracks -> a uniform track selection is necessary to estimate the contribution of v_n harmonics
10 // aod's and esd's are handled transparently
11 // The task will estimates a phi-dependent background density rho
12 // which is added to the event as a AliLocalRhoParamter object
13 //
14 // Author: Redmer Alexander Bertens, Utrecht Univeristy, Utrecht, Netherlands
15 // (rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl)
16 
17 // root includes
18 #include <TStyle.h>
19 #include <TRandom3.h>
20 #include <TChain.h>
21 #include <TMath.h>
22 #include <TF1.h>
23 #include <TF2.h>
24 #include <TH1F.h>
25 #include <TH2F.h>
26 #include <TProfile.h>
27 // aliroot includes
28 #include <AliAnalysisTask.h>
29 #include <AliAnalysisManager.h>
30 #include <AliCentrality.h>
31 #include <AliVVertex.h>
32 #include <AliESDEvent.h>
33 #include <AliAODEvent.h>
34 #include <AliAODTrack.h>
35 // emcal jet framework includes
36 #include <AliPicoTrack.h>
37 #include <AliEmcalJet.h>
38 #include <AliRhoParameter.h>
39 #include <AliLocalRhoParameter.h>
41 #include "AliMultSelection.h"
42 
44 using namespace std;
45 
47 
48 //_____________________________________________________________________________
60  fProfV3(0), fProfV3Cumulant(0)
61 {
62  // Default constructor
63 
64  for(Int_t i(0); i < 10; i++) {
65  fHistPsi2[i] = 0;
66  fHistPsi3[i] = 0;
67  }
68 }
69 
70 //_____________________________________________________________________________
72  AliAnalysisTaskEmcalJet(name, kTRUE),
82  fProfV3(0), fProfV3Cumulant(0)
83 {
84  // Constructor
85  for(Int_t i(0); i < 10; i++) {
86  fHistPsi2[i] = 0;
87  fHistPsi3[i] = 0;
88  }
89 
90  DefineInput(0, TChain::Class());
91  DefineOutput(1, TList::Class());
92  switch (fRunModeType) {
93  case kLocal : {
94  gStyle->SetOptFit(1);
95  DefineOutput(2, TList::Class());
96  DefineOutput(3, TList::Class());
97  } break;
98  default: break; // suppress debug info explicitely when not running locally
99  }
100 }
101 
102 //_____________________________________________________________________________
104 {
105  // destructor
106  if(fOutputList) delete fOutputList;
107  if(fOutputListGood) delete fOutputListGood;
108  if(fOutputListBad) delete fOutputListBad;
109  if(fFitModulation) delete fFitModulation;
110  if(fHistSwap) delete fHistSwap;
111 }
112 
113 //_____________________________________________________________________________
115 {
116  #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
117  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
118  #endif
119  // Init the analysis
120  if(fLocalRhoName=="") fLocalRhoName = Form("LocalRhoFrom_%s", GetName());
121  fLocalRho = new AliLocalRhoParameter(fLocalRhoName.Data(), 0);
122  // add the local rho to the event if necessary
123  if(fAttachToEvent) {
124  if(!(InputEvent()->FindListObject(fLocalRho->GetName()))) {
125  InputEvent()->AddObject(fLocalRho);
126  } else {
127  AliFatal(Form("%s: Container with same name %s already present. Aborting", GetName(), fLocalRho->GetName()));
128  }
129  }
130  AliAnalysisTaskEmcalJet::ExecOnce(); // init the base clas
131  if(fUseScaledRho) {
132  // unscaled rho has been retrieved by the parent class, now we retrieve rho scaled
133  fRho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(Form("%s_Scaled", fRho->GetName())));
134  if(!fRho) {
135  AliFatal(Form("%s: Couldn't find container for scaled rho. Aborting !", GetName()));
136  }
137  }
138  if(!GetJetContainer()) AliFatal(Form("%s: Couldn't get jet container. Aborting !", GetName()));
139 }
140 
141 //_____________________________________________________________________________
143 {
144  // Initialize the anaysis
145  #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
146  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
147  #endif
150  switch (fFitModulationType) {
151  case kNoFit : { SetModulationFit(new TF1("fit_kNoFit", "[0]", 0, TMath::TwoPi())); } break;
152  case kV2 : {
153  SetModulationFit(new TF1("fit_kV2", "[0]*([1]+[2]*[3]*TMath::Cos([2]*(x-[4])))", 0, TMath::TwoPi()));
154  fFitModulation->SetParameter(0, 0.); // normalization
155  fFitModulation->SetParameter(3, 0.2); // v2
156  fFitModulation->FixParameter(1, 1.); // constant
157  fFitModulation->FixParameter(2, 2.); // constant
158  } break;
159  case kV3: {
160  SetModulationFit(new TF1("fit_kV3", "[0]*([1]+[2]*[3]*TMath::Cos([2]*(x-[4])))", 0, TMath::TwoPi()));
161  fFitModulation->SetParameter(0, 0.); // normalization
162  fFitModulation->SetParameter(3, 0.2); // v3
163  fFitModulation->FixParameter(1, 1.); // constant
164  fFitModulation->FixParameter(2, 3.); // constant
165  } break;
166  default : { // for the combined fit, the 'direct fourier series' or the user supplied vn values we use v2 and v3
167  SetModulationFit(new TF1("fit_kCombined", "[0]*([1]+[2]*([3]*TMath::Cos([2]*(x-[4]))+[7]*TMath::Cos([5]*(x-[6]))))", 0, TMath::TwoPi()));
168  fFitModulation->SetParameter(0, 0.); // normalization
169  fFitModulation->SetParameter(3, 0.2); // v2
170  fFitModulation->FixParameter(1, 1.); // constant
171  fFitModulation->FixParameter(2, 2.); // constant
172  fFitModulation->FixParameter(5, 3.); // constant
173  fFitModulation->SetParameter(7, 0.2); // v3
174  } break;
175  }
176  switch (fRunModeType) {
177  case kGrid : { fFitModulationOptions += "N0"; } break;
178  default : break;
179  }
181  return kTRUE;
182 }
183 
184 //_____________________________________________________________________________
186 {
187  // create output objects
188  #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
189  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
190  #endif
191  fHistSwap = new TH1F("fHistSwap", "fHistSwap", 20, 0, TMath::TwoPi());
192  if(!fCentralityClasses) { // classes must be defined at this point
193  Int_t c[] = {0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100};
194  fCentralityClasses = new TArrayI(sizeof(c)/sizeof(c[0]), c);
195  }
196  fOutputList = new TList();
197  fOutputList->SetOwner(kTRUE);
198  // the analysis summary histo which stores all the analysis flags is always written to file
199  fHistAnalysisSummary = BookTH1F("fHistAnalysisSummary", "flag", 51, -0.5, 51.5);
200  if(!fFillHistograms) {
201  PostData(1, fOutputList);
202  return;
203  }
204  for(Int_t i(0); i < fCentralityClasses->GetSize()-1; i++) {
205  fHistPsi2[i] = BookTH1F("fHistPsi2", "#Psi_{2}", 100, -.5*TMath::Pi(), .5*TMath::Pi(), i);
206  fHistPsi3[i] = BookTH1F("fHistPsi3", "#Psi_{3}", 100, -1.*TMath::Pi()/3., TMath::Pi()/3., i);
207  }
208  // cdf of chisquare distribution
209  fHistPvalueCDF = BookTH1F("fHistPvalueCDF", "CDF #chi^{2}", 500, 0, 1);
210  fHistRhoStatusCent = BookTH2F("fHistRhoStatusCent", "centrality", "status [0=ok, 1=failed]", 101, -1, 100, 2, -.5, 1.5);
211  // vn profiles
212  Float_t temp[fCentralityClasses->GetSize()];
213  for(Int_t i(0); i < fCentralityClasses->GetSize(); i++) temp[i] = fCentralityClasses->At(i);
214  fProfV2 = new TProfile("fProfV2", "fProfV2", fCentralityClasses->GetSize()-1, temp);
215  fProfV3 = new TProfile("fProfV3", "fProfV3", fCentralityClasses->GetSize()-1, temp);
216  fOutputList->Add(fProfV2);
217  fOutputList->Add(fProfV3);
218  switch (fFitModulationType) {
219  case kQC2 : {
220  fProfV2Cumulant = new TProfile("fProfV2Cumulant", "fProfV2Cumulant", fCentralityClasses->GetSize()-1, temp);
221  fProfV3Cumulant = new TProfile("fProfV3Cumulant", "fProfV3Cumulant", fCentralityClasses->GetSize()-1, temp);
223  fOutputList->Add(fProfV3Cumulant);
224  } break;
225  case kQC4 : {
226  fProfV2Cumulant = new TProfile("fProfV2Cumulant", "fProfV2Cumulant", fCentralityClasses->GetSize()-1, temp);
227  fProfV3Cumulant = new TProfile("fProfV3Cumulant", "fProfV3Cumulant", fCentralityClasses->GetSize()-1, temp);
229  fOutputList->Add(fProfV3Cumulant);
230  } break;
231  default : break;
232  }
233  if(fUsePtWeight) fHistSwap->Sumw2();
238  // increase readability of output list
239  fOutputList->Sort();
240  PostData(1, fOutputList);
241  switch (fRunModeType) {
242  case kLocal : {
243  fOutputListGood = new TList();
244  fOutputListGood->SetOwner(kTRUE);
245  fOutputListBad = new TList();
246  fOutputListBad->SetOwner(kTRUE);
247  PostData(2, fOutputListGood);
248  PostData(3, fOutputListBad);
249  } break;
250  default: break;
251  }
252 }
253 
254 //_____________________________________________________________________________
255 TH1F* AliAnalysisTaskLocalRho::BookTH1F(const char* name, const char* x, Int_t bins, Double_t min, Double_t max, Int_t c, Bool_t append)
256 {
257  // Book a TH1F and connect it to the output container
258  #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
259  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
260  #endif
261  if(!fOutputList) return 0x0;
262  TString title(name);
263  if(c!=-1) { // format centrality dependent histograms accordingly
264  name = Form("%s_%i", name, c);
265  title += Form("_%i-%i", fCentralityClasses->At(c), fCentralityClasses->At(1+c));
266  }
267  title += Form(";%s;[counts]", x);
268  TH1F* histogram = new TH1F(name, title.Data(), bins, min, max);
269  histogram->Sumw2();
270  if(append) fOutputList->Add(histogram);
271  return histogram;
272 }
273 
274 //_____________________________________________________________________________
275 TH2F* AliAnalysisTaskLocalRho::BookTH2F(const char* name, const char* x, const char*y, Int_t binsx, Double_t minx, Double_t maxx,
276  Int_t binsy, Double_t miny, Double_t maxy, Int_t c, Bool_t append)
277 {
278  // Book a TH2F and connect it to the output container
279  #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
280  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
281  #endif
282  if(!fOutputList) return 0x0;
283  TString title(name);
284  if(c!=-1) { // format centrality dependent histograms accordingly
285  name = Form("%s_%i", name, c);
286  title += Form("_%i-%i", fCentralityClasses->At(c), fCentralityClasses->At(1+c));
287  }
288  title += Form(";%s;%s", x, y);
289  TH2F* histogram = new TH2F(name, title.Data(), binsx, minx, maxx, binsy, miny, maxy);
290  histogram->Sumw2();
291  if(append) fOutputList->Add(histogram);
292  return histogram;
293 }
294 
295 //_____________________________________________________________________________
297 {
298  // Execute once for each event
299  #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
300  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
301  #endif
302  if(!(InputEvent()||fTracks||fJets||fRho)) return kFALSE;
304  // get the centrality bin (necessary for some control histograms
305  if(InputEvent()->GetRunNumber()>200000)
306  {
307  Float_t lPercentile = 300;
308  AliMultSelection *MultSelection = 0x0;
309  MultSelection = (AliMultSelection*)InputEvent()->FindListObject("MultSelection");
310  if(!MultSelection) {
311  //If you get this warning (and lPercentiles 300) please check that the AliMultSelectionTask actually ran (before your task)
312  AliWarning("AliMultSelection object not found!");
313  }else{
314  lPercentile = MultSelection->GetMultiplicityPercentile("V0M");
315  }
316  fCent = lPercentile;
317  }
319 
320  for(Int_t i(0); i < fCentralityClasses->GetSize()-1; i++) {
321  if(fCent >= fCentralityClasses->At(i) && fCent <= fCentralityClasses->At(1+i)) {
323  break; }
324  }
325  if(fInCentralitySelection < 0) return kFALSE;
326  // set the rho value
327  fLocalRho->SetVal(fRho->GetVal());
328  // set the correct event plane accordign to the requested reference detector
329  Double_t psi2(-1), psi3(-1);
330  switch (fDetectorType) { // determine the detector type for the rho fit
331  case kTPC : {
332  // [0] psi2 [1] psi3
333  Double_t tpc[2];
335  psi2 = tpc[0]; psi3 = tpc[1];
336  } break;
337  case kVZEROA : {
338  // [0][0] psi2a [1,0] psi2c
339  // [0][1] psi3a [1,1] psi3c
340  Double_t vzero[2][2];
342  psi2 = vzero[0][0]; psi3 = vzero[0][1];
343  } break;
344  case kVZEROC : {
345  // [0][0] psi2a [1,0] psi2c
346  // [0][1] psi3a [1,1] psi3c
347  Double_t vzero[2][2];
349  psi2 = vzero[1][0]; psi3 = vzero[1][1];
350  } break;
351  case kVZEROComb : {
352  /* for the combined vzero event plane
353  * [0] psi2 [1] psi3
354  * not fully implmemented yet, use with caution ! */
355  Double_t vzeroComb[2];
357  psi2 = vzeroComb[0]; psi3 = vzeroComb[1];
358  } break;
359  default : break;
360  }
362  switch (fFitModulationType) { // do the fits
363  case kNoFit : { fFitModulation->FixParameter(0, fLocalRho->GetVal()); } break;
364  case kV2 : { // only v2
365  if(CorrectRho(psi2, psi3)) {
366  if(fFillHistograms) fProfV2->Fill(fCent, fFitModulation->GetParameter(3));
367  if(fUserSuppliedR2) {
368  Double_t r(fUserSuppliedR2->GetBinContent(fUserSuppliedR2->GetXaxis()->FindBin(fCent)));
369  if(r > 0) fFitModulation->SetParameter(3, fFitModulation->GetParameter(3)/r);
370  }
371  }
372  } break;
373  case kV3 : { // only v3
374  if(CorrectRho(psi2, psi3)) {
375  if(fUserSuppliedR3) {
376  Double_t r(fUserSuppliedR3->GetBinContent(fUserSuppliedR3->GetXaxis()->FindBin(fCent)));
377  if(r > 0) fFitModulation->SetParameter(3, fFitModulation->GetParameter(3)/r);
378  }
379  if(fFillHistograms) fProfV3->Fill(fCent, fFitModulation->GetParameter(3));
380  }
381  } break;
382  case kQC2 : { // qc2 analysis - NOTE: not a wise idea to use this !
383  if(CorrectRho(psi2, psi3)) {
385  // note for the qc method, resolution is REVERSED to go back to v2obs
386  Double_t r2(fUserSuppliedR2->GetBinContent(fUserSuppliedR2->GetXaxis()->FindBin(fCent)));
387  Double_t r3(fUserSuppliedR3->GetBinContent(fUserSuppliedR3->GetXaxis()->FindBin(fCent)));
388  if(r2 > 0) fFitModulation->SetParameter(3, fFitModulation->GetParameter(3)*r2);
389  if(r3 > 0) fFitModulation->SetParameter(7, fFitModulation->GetParameter(7)*r3);
390  }
391  if (fUsePtWeight) { // use weighted weights
392  Double_t dQCnM11 = (fNoEventWeightsForQC) ? 1. : QCnM11();
393  if(fFillHistograms) {
394  fProfV2->Fill(fCent, fFitModulation->GetParameter(3), dQCnM11);
395  fProfV3->Fill(fCent, fFitModulation->GetParameter(7), dQCnM11);
396  }
397  } else {
398  Double_t dQCnM = (fNoEventWeightsForQC) ? 2. : QCnM();
399  if(fFillHistograms) {
400  fProfV2->Fill(fCent, fFitModulation->GetParameter(3), dQCnM*(dQCnM-1));
401  fProfV3->Fill(fCent, fFitModulation->GetParameter(7), dQCnM*(dQCnM-1));
402  }
403  }
404  }
405  } break;
406  case kQC4 : { // NOTE: see comment at kQC2
407  if(CorrectRho(psi2, psi3)) {
409  // note for the qc method, resolution is REVERSED to go back to v2obs
410  Double_t r2(fUserSuppliedR2->GetBinContent(fUserSuppliedR2->GetXaxis()->FindBin(fCent)));
411  Double_t r3(fUserSuppliedR3->GetBinContent(fUserSuppliedR3->GetXaxis()->FindBin(fCent)));
412  if(r2 > 0) fFitModulation->SetParameter(3, fFitModulation->GetParameter(3)*r2);
413  if(r3 > 0) fFitModulation->SetParameter(7, fFitModulation->GetParameter(7)*r3);
414  }
415  if (fUsePtWeight) { // use weighted weights
416  if(fFillHistograms) {
417  fProfV2->Fill(fCent, TMath::Power(fFitModulation->GetParameter(3),0.5)/*, QCnM1111()*/);
418  fProfV3->Fill(fCent, TMath::Power(fFitModulation->GetParameter(7),0.5)/*, QCnM1111()*/);
419  }
420  } else {
421  if(fFillHistograms) {
422  fProfV2->Fill(fCent, TMath::Power(fFitModulation->GetParameter(3),0.5)/*, QCnM()*(QCnM()-1)*(QCnM()-2)*(QCnM()-3)*/);
423  fProfV3->Fill(fCent, TMath::Power(fFitModulation->GetParameter(7),0.5)/*, QCnM()*(QCnM()-1)*(QCnM()-2)*(QCnM()-3)*/);
424  }
425  }
426  }
427  } break;
428  default : {
429  if(CorrectRho(psi2, psi3)) {
431  Double_t r2(fUserSuppliedR2->GetBinContent(fUserSuppliedR2->GetXaxis()->FindBin(fCent)));
432  Double_t r3(fUserSuppliedR3->GetBinContent(fUserSuppliedR3->GetXaxis()->FindBin(fCent)));
433  if(r2 > 0) fFitModulation->SetParameter(3, fFitModulation->GetParameter(3)/r2);
434  if(r3 > 0) fFitModulation->SetParameter(7, fFitModulation->GetParameter(7)/r3);
435  }
436  if(fFillHistograms) {
437  fProfV2->Fill(fCent, fFitModulation->GetParameter(3));
438  fProfV3->Fill(fCent, fFitModulation->GetParameter(7));
439  }
440  }
441  } break;
442  }
443  // if all went well, add local rho
445  PostData(1, fOutputList);
446  return kTRUE;
447 }
448 
449 //_____________________________________________________________________________
451 {
452  #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
453  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
454  #endif
455  // Get the vzero event plane
457  // use the vzero event plane from the event header
458  // note: to use the calibrated vzero event plane, run
459  // $ALICE_ROOT/ANALYSIS/macros/AddTaskVZEROEPSelection.C
460  // prior to this task (make sure the calibration is available for the dataset
461  // you want to use)
462  Double_t a(0), b(0), c(0), d(0), e(0), f(0), g(0), h(0);
463  vzero[0][0] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 2, a, b);
464  vzero[1][0] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 2, c, d);
465  vzero[0][1] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 3, e, f);
466  vzero[1][1] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 3, g, h);
467  return;
468  }
469  // grab the vzero event plane without recentering
470  if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
471  Double_t qxa2(0), qya2(0), qxc2(0), qyc2(0); // for psi2
472  Double_t qxa3(0), qya3(0), qxc3(0), qyc3(0); // for psi3
473  for(Int_t iVZERO(0); iVZERO < 64; iVZERO++) {
474  Double_t phi(TMath::PiOver4()*(.5+iVZERO%8)), /* eta(0), */ weight(InputEvent()->GetVZEROEqMultiplicity(iVZERO));
475  // (iVZERO<32) ? eta = -3.45+.5*(iVZERO/8) : eta = 4.8-.6*((iVZERO/8)-4);
476  if(iVZERO<32) {
477  qxa2 += weight*TMath::Cos(2.*phi);
478  qya2 += weight*TMath::Sin(2.*phi);
479  qxa3 += weight*TMath::Cos(3.*phi);
480  qya3 += weight*TMath::Sin(3.*phi);
481  }
482  else {
483  qxc2 += weight*TMath::Cos(2.*phi);
484  qyc2 += weight*TMath::Sin(2.*phi);
485  qxc3 += weight*TMath::Cos(3.*phi);
486  qyc3 += weight*TMath::Sin(3.*phi);
487  }
488  }
489  vzero[0][0] = .5*TMath::ATan2(qya2, qxa2);
490  vzero[1][0] = .5*TMath::ATan2(qyc2, qxc2);
491  vzero[0][1] = (1./3.)*TMath::ATan2(qya3, qxa3);
492  vzero[1][1] = (1./3.)*TMath::ATan2(qyc3, qxc3);
493 }
494 
495 //_____________________________________________________________________________
497 {
498  // Grab the TPC event plane. if parameter fExcludeLeadingJetsFromFit is larger than 0,
499  // strip in eta of width fExcludeLeadingJetsFromFit * GetJetContainer()->GetJetRadius() around the leading jet (before
500  // subtraction of rho) will be exluded from the event plane estimate
501  #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
502  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
503  #endif
504  fNAcceptedTracks = 0; // reset the track counter
505  Double_t qx2(0), qy2(0); // for psi2
506  Double_t qx3(0), qy3(0); // for psi3
507  if(fTracks) {
508  Float_t excludeInEta = -999;
509  if(fExcludeLeadingJetsFromFit > 0 ) { // remove the leading jet from ep estimate
510  AliEmcalJet* leadingJet(GetJetContainer()->GetLeadingJet());
511  if(leadingJet) leadingJet->Eta();
512  }
513  Int_t iTracks(fTracks->GetEntriesFast());
514  for(Int_t iTPC(0); iTPC < iTracks; iTPC++) {
515  AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(iTPC));
516  if(!PassesCuts(track) || track->Pt() < fSoftTrackMinPt || track->Pt() > fSoftTrackMaxPt) continue;
517  if(fExcludeLeadingJetsFromFit > 0 &&( (TMath::Abs(track->Eta() - excludeInEta) < GetJetContainer()->GetJetRadius()*fExcludeLeadingJetsFromFit ) || (TMath::Abs(track->Eta()) - GetJetContainer()->GetJetRadius() - GetJetContainer()->GetJetEtaMax() ) > 0 )) continue;
519  qx2+= TMath::Cos(2.*track->Phi());
520  qy2+= TMath::Sin(2.*track->Phi());
521  qx3+= TMath::Cos(3.*track->Phi());
522  qy3+= TMath::Sin(3.*track->Phi());
523  }
524  }
525  tpc[0] = .5*TMath::ATan2(qy2, qx2);
526  tpc[1] = (1./3.)*TMath::ATan2(qy3, qx3);
527 }
528 
529 //_____________________________________________________________________________
531 {
532  // Grab the combined vzero event plane
533  #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
534  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
535  #endif
536  Double_t a(0), b(0), c(0), d(0);
537  comb[0] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 10, 2, a, b);
538  comb[1] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 10, 3, c, d);
539 }
540 
541 //_____________________________________________________________________________
543 {
544  // Get the second order q-cumulant, a -999 return will be caught in the qa routine of CorrectRho
545  #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
546  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
547  #endif
548  Double_t reQ(0), imQ(0), modQ(0), M11(0), M(0);
549  if(fUsePtWeight) { // for the weighted 2-nd order q-cumulant
550  QCnQnk(harm, 1, reQ, imQ); // get the weighted 2-nd order q-vectors
551  modQ = reQ*reQ+imQ*imQ; // get abs Q-squared
552  M11 = QCnM11(); // equals S2,1 - S1,2
553  return (M11 > 0) ? ((modQ - QCnS(1,2))/M11) : -999;
554  } // else return the non-weighted 2-nd order q-cumulant
555  QCnQnk(harm, 0, reQ, imQ); // get the non-weighted 2-nd order q-vectors
556  modQ = reQ*reQ+imQ*imQ; // get abs Q-squared
557  M = QCnM();
558  return (M > 1) ? (modQ - M)/(M*(M-1)) : -999;
559 }
560 
561 //_____________________________________________________________________________
563 {
564  // Get the fourth order q-cumulant, a -999 return will be caught in the qa routine of CorrectRho
565  #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
566  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
567  #endif
568  Double_t reQn1(0), imQn1(0), reQ2n2(0), imQ2n2(0), reQn3(0), imQn3(0), M1111(0), M(0);
569  Double_t a(0), b(0), c(0), d(0), e(0), f(0), g(0); // terms of the calculation
570  if(fUsePtWeight) { // for the weighted 4-th order q-cumulant
571  QCnQnk(harm, 1, reQn1, imQn1);
572  QCnQnk(harm*2, 2, reQ2n2, imQ2n2);
573  QCnQnk(harm, 3, reQn3, imQn3);
574  // fill in the terms ...
575  a = (reQn1*reQn1+imQn1*imQn1)*(reQn1*reQn1+imQn1*imQn1);
576  b = reQ2n2*reQ2n2 + imQ2n2*imQ2n2;
577  c = -2.*(reQ2n2*reQn1*reQn1-reQ2n2*imQn1*imQn1+2.*imQ2n2*reQn1*imQn1);
578  d = 8.*(reQn3*reQn1+imQn3*imQn1);
579  e = -4.*QCnS(1,2)*(reQn1*reQn1+imQn1*imQn1);
580  f = -6.*QCnS(1,4);
581  g = 2.*QCnS(2,2);
582  M1111 = QCnM1111();
583  return (M1111 > 0) ? (a+b+c+d+e+f+g)/M1111 : -999;
584  } // else return the unweighted case
585  Double_t reQn(0), imQn(0), reQ2n(0), imQ2n(0);
586  QCnQnk(harm, 0, reQn, imQn);
587  QCnQnk(harm*2, 0, reQ2n, imQ2n);
588  // fill in the terms ...
589  M = QCnM();
590  if(M < 4) return -999;
591  a = (reQn*reQn+imQn*imQn)*(reQn*reQn+imQn*imQn);
592  b = reQ2n*reQ2n + imQ2n*imQ2n;
593  c = -2.*(reQ2n*reQn*reQn-reQ2n*imQn*imQn+2.*imQ2n*reQn*imQn);
594  e = -4.*(M-2)*(reQn*reQn+imQn*imQn);
595  f = 2.*M*(M-3);
596  return (a+b+c+e+f)/(M*(M-1)*(M-2)*(M-3));
597 }
598 
599 //_____________________________________________________________________________
601 {
602  // Get the weighted n-th order q-vector, pass real and imaginary part as reference
603  #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
604  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
605  #endif
606  if(!fTracks) return;
608  Int_t iTracks(fTracks->GetEntriesFast());
609  for(Int_t iTPC(0); iTPC < iTracks; iTPC++) {
610  AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(iTPC));
611  if(!PassesCuts(track) || track->Pt() < fSoftTrackMinPt || track->Pt() > fSoftTrackMaxPt) continue;
613  // for the unweighted case, k equals zero and the weight doesn't contribute to the equation below
614  reQ += TMath::Power(track->Pt(), k) * TMath::Cos(((double)n)*track->Phi());
615  imQ += TMath::Power(track->Pt(), k) * TMath::Sin(((double)n)*track->Phi());
616  }
617 }
618 
619 //_____________________________________________________________________________
621 {
622  // Get the weighted ij-th order autocorrelation correction
623  #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
624  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
625  #endif
626  if(!fTracks || i <= 0 || j <= 0) return -999;
627  Int_t iTracks(fTracks->GetEntriesFast());
628  Double_t Sij(0);
629  for(Int_t iTPC(0); iTPC < iTracks; iTPC++) {
630  AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(iTPC));
631  if(!PassesCuts(track) || track->Pt() < fSoftTrackMinPt || track->Pt() > fSoftTrackMaxPt) continue;
632  Sij+=TMath::Power(track->Pt(), j);
633  }
634  return TMath::Power(Sij, i);
635 }
636 
637 //_____________________________________________________________________________
639 {
640  // Get multiplicity for unweighted q-cumulants. function QCnQnk should be called first
641  #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
642  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
643  #endif
644  return (Double_t) fNAcceptedTracksQCn;
645 }
646 
647 //_____________________________________________________________________________
649 {
650  // Get multiplicity weights for the weighted two particle cumulant
651  #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
652  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
653  #endif
654  return (QCnS(2,1) - QCnS(1,2));
655 }
656 
657 //_____________________________________________________________________________
659 {
660  // Get multiplicity weights for the weighted four particle cumulant
661  #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
662  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
663  #endif
664  return (QCnS(4,1)-6*QCnS(1,2)*QCnS(2,1)+8*QCnS(1,3)*QCnS(1,1)+3*QCnS(2,2)-6*QCnS(1,4));
665 }
666 
667 //_____________________________________________________________________________
669 {
670  // Decides how to deal with the situation where c2 or c3 is negative
671  // Returns kTRUE depending on whether or not a modulated rho is used for the jet background
672  #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
673  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
674  #endif
675  if(TMath::AreEqualAbs(fFitModulation->GetParameter(3), .0, 1e-10) && TMath::AreEqualAbs(fFitModulation->GetParameter(7), .0,1e-10)) {
676  fFitModulation->SetParameter(7, 0);
677  fFitModulation->SetParameter(3, 0);
678  fFitModulation->SetParameter(0, fLocalRho->GetVal());
679  return kTRUE; // v2 and v3 have physical null values
680  }
681  switch (fQCRecovery) {
682  case kFixedRho : { // roll back to the original rho
683  fFitModulation->SetParameter(7, 0);
684  fFitModulation->SetParameter(3, 0);
685  fFitModulation->SetParameter(0, fLocalRho->GetVal());
686  return kFALSE; // rho is forced to be fixed
687  }
688  case kNegativeVn : {
689  Double_t c2(fFitModulation->GetParameter(3));
690  Double_t c3(fFitModulation->GetParameter(7));
691  if( c2 < 0 ) c2 = -1.*TMath::Sqrt(-1.*c2);
692  if( c3 < 0 ) c3 = -1.*TMath::Sqrt(-1.*c3);
693  fFitModulation->SetParameter(3, c2);
694  fFitModulation->SetParameter(7, c3);
695  return kTRUE; // is this a physical quantity ?
696  }
697  case kTryFit : {
698  fitModulationType tempType(fFitModulationType); // store temporarily
700  fFitModulation->SetParameter(7, 0);
701  fFitModulation->SetParameter(3, 0);
702  Bool_t pass(CorrectRho(psi2, psi3)); // do the fit and all quality checks
703  fFitModulationType = tempType; // roll back for next event
704  return pass;
705  }
706  default : return kFALSE;
707  }
708  return kFALSE;
709 }
710 
711 //_____________________________________________________________________________
713 {
714  // Get rho' -> rho(phi)
715  // three routines are available, 1 and 2 can be used with or without pt weights
716  // [1] get vn from q-cumulants
717  // in case of cumulants, both cumulants and vn values are stored. in both cases, v2 and v3
718  // are expected. a check is performed to see if rho has no negative local minimum
719  // for full description, see Phys. Rev. C 83, 044913
720  // since the cn distribution has negative values, vn = sqrt(cn) can be imaginary sometimes
721  // in this case one can either roll back to the 'original' fixed rho, do a fit for vn or take use
722  // vn = - sqrt(|cn|) note that because of this, use of q-cumulants is not safe !
723  // [2] fitting a fourier expansion to the de/dphi distribution
724  // the fit can be done with either v2, v3 or a combination.
725  // in all cases, a cut can be made on the p-value of the chi-squared value of the fit
726  // and a check can be performed to see if rho has no negative local minimum
727  // [3] get v2 and v3 from user supplied histograms
728  // in this way, a fixed value of v2 and v3 is subtracted w.r.t. whichever event plane is requested
729  #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
730  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
731  #endif
732  switch (fFitModulationType) { // for approaches where no fitting is required
733  case kQC2 : {
734  fFitModulation->FixParameter(4, psi2);
735  fFitModulation->FixParameter(6, psi3);
736  fFitModulation->FixParameter(3, CalculateQC2(2)); // set here with cn, vn = sqrt(cn)
737  fFitModulation->FixParameter(7, CalculateQC2(3));
738  // first fill the histos of the raw cumulant distribution
739  if (fUsePtWeight) { // use weighted weights
740  Double_t dQCnM11 = (fNoEventWeightsForQC) ? 1. : QCnM11();
741  if(fFillHistograms) {
742  fProfV2Cumulant->Fill(fCent, fFitModulation->GetParameter(3), dQCnM11);
743  fProfV3Cumulant->Fill(fCent, fFitModulation->GetParameter(7), dQCnM11);
744  }
745  } else {
746  Double_t dQCnM = (fNoEventWeightsForQC) ? 2. : QCnM();
747  if(fFillHistograms) {
748  fProfV2Cumulant->Fill(fCent, fFitModulation->GetParameter(3), dQCnM*(dQCnM-1));
749  fProfV3Cumulant->Fill(fCent, fFitModulation->GetParameter(7), dQCnM*(dQCnM-1));
750  }
751  }
752  // then see if one of the cn value is larger than zero and vn is readily available
753  if(fFitModulation->GetParameter(3) > 0 && fFitModulation->GetParameter(7) > 0) {
754  fFitModulation->FixParameter(3, TMath::Sqrt(fFitModulation->GetParameter(3)));
755  fFitModulation->FixParameter(7, TMath::Sqrt(fFitModulation->GetParameter(7)));
756  } else if (!QCnRecovery(psi2, psi3)) return kFALSE; // try to recover the cumulant, this will set v2 and v3
757  if(fAbsVnHarmonics && fFitModulation->GetMinimum(0, TMath::TwoPi()) < 0) { // general check
758  fFitModulation->SetParameter(7, 0);
759  fFitModulation->SetParameter(3, 0);
760  fFitModulation->SetParameter(0, fLocalRho->GetVal());
761  return kFALSE;
762  }
763  return kTRUE;
764  } break;
765  case kQC4 : {
766  fFitModulation->FixParameter(4, psi2);
767  fFitModulation->FixParameter(6, psi3);
768  fFitModulation->FixParameter(3, CalculateQC4(2)); // set here with cn, vn = sqrt(cn)
769  fFitModulation->FixParameter(7, CalculateQC4(3));
770  // first fill the histos of the raw cumulant distribution
771  if (fUsePtWeight) { // use weighted weights
772  if(fFillHistograms) {
773  fProfV2Cumulant->Fill(fCent, fFitModulation->GetParameter(3)/*, QCnM1111()*/);
774  fProfV3Cumulant->Fill(fCent, fFitModulation->GetParameter(7)/*, QCnM1111()*/);
775  }
776  } else {
777  if(fFillHistograms) {
778  fProfV2Cumulant->Fill(fCent, fFitModulation->GetParameter(3)/*, QCnM1111()*/);
779  fProfV3Cumulant->Fill(fCent, fFitModulation->GetParameter(7)/*, QCnM1111()*/);
780  }
781  }
782  // then see if one of the cn value is larger than zero and vn is readily available
783  if(fFitModulation->GetParameter(3) > 0 && fFitModulation->GetParameter(7) > 0) {
784  fFitModulation->FixParameter(3, TMath::Sqrt(fFitModulation->GetParameter(3)));
785  fFitModulation->FixParameter(7, TMath::Sqrt(fFitModulation->GetParameter(7)));
786  } else if (!QCnRecovery(psi2, psi3)) return kFALSE; // try to recover the cumulant, this will set v2 and v3
787  if(fAbsVnHarmonics && fFitModulation->GetMinimum(0, TMath::TwoPi()) < 0) { // general check
788  fFitModulation->SetParameter(7, 0);
789  fFitModulation->SetParameter(3, 0);
790  fFitModulation->SetParameter(0, fLocalRho->GetVal());
791  return kFALSE;
792  }
793  } break;
794  case kIntegratedFlow : {
795  // use v2 and v3 values from an earlier iteration over the data
796  fFitModulation->FixParameter(3, fUserSuppliedV2->GetBinContent(fUserSuppliedV2->GetXaxis()->FindBin(fCent)));
797  fFitModulation->FixParameter(4, psi2);
798  fFitModulation->FixParameter(6, psi3);
799  fFitModulation->FixParameter(7, fUserSuppliedV3->GetBinContent(fUserSuppliedV3->GetXaxis()->FindBin(fCent)));
800  if(fAbsVnHarmonics && fFitModulation->GetMinimum(0, TMath::TwoPi()) < 0) {
801  fFitModulation->SetParameter(7, 0);
802  fFitModulation->SetParameter(3, 0);
803  fFitModulation->SetParameter(0, fLocalRho->GetVal());
804  return kFALSE;
805  }
806  return kTRUE;
807  }
808  default : break;
809  }
810  TString detector("");
811  switch (fDetectorType) {
812  case kTPC : detector+="TPC";
813  break;
814  case kVZEROA : detector+="VZEROA";
815  break;
816  case kVZEROC : detector+="VZEROC";
817  break;
818  case kVZEROComb : detector+="VZEROComb";
819  break;
820  default: break;
821  }
822  Int_t iTracks(fTracks->GetEntriesFast());
823  Double_t excludeInEta = -999;
824  Double_t excludeInPhi = -999;
825  Double_t excludeInPt = -999;
826  if(iTracks <= 0 || fLocalRho->GetVal() <= 0 ) return kFALSE; // no use fitting an empty event ...
827  if(fExcludeLeadingJetsFromFit > 0 ) {
828  AliEmcalJet* leadingJet = GetJetContainer()->GetLeadingJet();
829  if(PassesCuts(leadingJet)) {
830  excludeInEta = leadingJet->Eta();
831  excludeInPhi = leadingJet->Phi();
832  excludeInPt = leadingJet->Pt();
833  }
834  }
835  fHistSwap->Reset(); // clear the histogram
836  TH1F _tempSwap; // on stack for quick access
837  TH1F _tempSwapN; // on stack for quick access, bookkeeping histogram
839  if(fNAcceptedTracks < 49) fNAcceptedTracks = 49; // avoid aliasing effects
840  _tempSwap = TH1F("_tempSwap", "_tempSwap", TMath::CeilNint(TMath::Sqrt(fNAcceptedTracks)), 0, TMath::TwoPi());
841  if(fUsePtWeightErrorPropagation) _tempSwapN = TH1F("_tempSwapN", "_tempSwapN", TMath::CeilNint(TMath::Sqrt(fNAcceptedTracks)), 0, TMath::TwoPi());
842  if(fUsePtWeight) _tempSwap.Sumw2();
843  }
844  else _tempSwap = *fHistSwap; // now _tempSwap holds the desired histo
845  // non poissonian error when using pt weights
846  Double_t totalpts(0.), totalptsquares(0.), totalns(0.);
847  for(Int_t i(0); i < iTracks; i++) {
848  AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(i));
849  if(fExcludeLeadingJetsFromFit > 0 &&( (TMath::Abs(track->Eta() - excludeInEta) < GetJetContainer()->GetJetRadius()*fExcludeLeadingJetsFromFit ) || (TMath::Abs(track->Eta()) - GetJetContainer()->GetJetRadius() - GetJetContainer()->GetJetEtaMax() ) > 0 )) continue;
850  if(!PassesCuts(track) || track->Pt() > fSoftTrackMaxPt || track->Pt() < fSoftTrackMinPt) continue;
851  if(fUsePtWeight) {
852  _tempSwap.Fill(track->Phi(), track->Pt());
854  totalpts += track->Pt();
855  totalptsquares += track->Pt()*track->Pt();
856  totalns += 1;
857  _tempSwapN.Fill(track->Phi());
858  }
859  }
860  else _tempSwap.Fill(track->Phi());
861  }
863  // in the case of pt weights overwrite the poissonian error estimate which is assigned by root by a more sophisticated appraoch
864  // the assumption here is that the bin error will be dominated by the uncertainty in the mean pt in a bin and in the uncertainty
865  // of the number of tracks in a bin, the first of which will be estimated from the sample standard deviation of all tracks in the
866  // event, for the latter use a poissonian estimate. the two contrubitions are assumed to be uncorrelated
867  if(totalns < 1) return kFALSE; // not one track passes the cuts
868  for(Int_t l = 0; l < _tempSwap.GetNbinsX(); l++) {
869  if(_tempSwapN.GetBinContent(l+1) == 0) {
870  _tempSwap.SetBinContent(l+1,0);
871  _tempSwap.SetBinError(l+1,0);
872  }
873  else {
874  Double_t vartimesnsq = totalptsquares*totalns - totalpts*totalpts;
875  Double_t variance = vartimesnsq/(totalns*(totalns-1.));
876  Double_t SDOMSq = variance / _tempSwapN.GetBinContent(l+1);
877  Double_t SDOMSqOverMeanSq = SDOMSq * _tempSwapN.GetBinContent(l+1) * _tempSwapN.GetBinContent(l+1) / (_tempSwapN.GetBinContent(l+1) * _tempSwapN.GetBinContent(l+1));
878  Double_t poissonfrac = 1./_tempSwapN.GetBinContent(l+1);
879  Double_t vartotalfrac = SDOMSqOverMeanSq + poissonfrac;
880  Double_t vartotal = vartotalfrac * _tempSwap.GetBinContent(l+1) * _tempSwap.GetBinContent(l+1);
881  if(vartotal > 0.0001) _tempSwap.SetBinError(l+1,TMath::Sqrt(vartotal));
882  else {
883  _tempSwap.SetBinContent(l+1,0);
884  _tempSwap.SetBinError(l+1,0);
885  }
886  }
887  }
888  }
889 
890  fFitModulation->SetParameter(0, fLocalRho->GetVal());
891  switch (fFitModulationType) {
892  case kNoFit : { fFitModulation->FixParameter(0, fLocalRho->GetVal() );
893  } break;
894  case kV2 : {
895  fFitModulation->FixParameter(4, psi2);
896  } break;
897  case kV3 : {
898  fFitModulation->FixParameter(4, psi3);
899  } break;
900  case kCombined : {
901  fFitModulation->FixParameter(4, psi2);
902  fFitModulation->FixParameter(6, psi3);
903  } break;
904  case kFourierSeries : {
905  // in this approach, an explicit calculation will be made of vn = sqrt(xn^2+yn^2)
906  // where x[y] = Integrate[r(phi)cos[sin](n phi)dphi, 0, 2pi]
907  Double_t cos2(0), sin2(0), cos3(0), sin3(0), sumPt(0);
908  for(Int_t i(0); i < iTracks; i++) {
909  AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(i));
910  if(!PassesCuts(track) || track->Pt() > fSoftTrackMaxPt || track->Pt() < fSoftTrackMinPt) continue;
911  sumPt += track->Pt();
912  cos2 += track->Pt()*TMath::Cos(2*PhaseShift(track->Phi()-psi2));
913  sin2 += track->Pt()*TMath::Sin(2*PhaseShift(track->Phi()-psi2));
914  cos3 += track->Pt()*TMath::Cos(3*PhaseShift(track->Phi()-psi3));
915  sin3 += track->Pt()*TMath::Sin(3*PhaseShift(track->Phi()-psi3));
916  }
917  fFitModulation->SetParameter(3, TMath::Sqrt(cos2*cos2+sin2*sin2)/fLocalRho->GetVal());
918  fFitModulation->SetParameter(4, psi2);
919  fFitModulation->SetParameter(6, psi3);
920  fFitModulation->SetParameter(7, TMath::Sqrt(cos3*cos3+sin3*sin3)/fLocalRho->GetVal());
921  } break;
922  default : break;
923  }
924  _tempSwap.Fit(fFitModulation, fFitModulationOptions.Data(), "", 0, TMath::TwoPi());
925  // the quality of the fit is evaluated from 1 - the cdf of the chi square distribution
926  Double_t CDF(1.-ChiSquareCDF(fFitModulation->GetNDF(), fFitModulation->GetChisquare()));
927  if(fFillHistograms) fHistPvalueCDF->Fill(CDF);
928  if(CDF > fMinPvalue && CDF < fMaxPvalue && ( fAbsVnHarmonics && fFitModulation->GetMinimum(0, TMath::TwoPi()) > 0)) { // fit quality
930  // for LOCAL didactic purposes, save the best and the worst fits
931  // this routine can produce a lot of output histograms (it's not memory 'safe') and will not work on GRID
932  // since the output will become unmergeable (i.e. different nodes may produce conflicting output)
933  switch (fRunModeType) {
934  case kLocal : {
935  if(gRandom->Uniform(0, 100) > fPercentageOfFits) break;
936  static Int_t didacticCounterBest(0);
937  TProfile* didacticProfile = (TProfile*)_tempSwap.Clone(Form("Fit_%i_1-CDF_%.3f_cen_%i_%s", didacticCounterBest, CDF, fInCentralitySelection, detector.Data()));
938  TF1* didactifFit = (TF1*)fFitModulation->Clone(Form("fit_%i_CDF_%.3f_cen_%i_%s", didacticCounterBest, CDF, fInCentralitySelection, detector.Data()));
939  didacticProfile->GetListOfFunctions()->Add(didactifFit);
940  fOutputListGood->Add(didacticProfile);
941  didacticCounterBest++;
942  TH2F* didacticSurface = BookTH2F(Form("surface_%s", didacticProfile->GetName()), "#phi", "#eta", 50, 0, TMath::TwoPi(), 50, -1, 1, -1, kFALSE);
943  for(Int_t i(0); i < iTracks; i++) {
944  AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(i));
945  if(PassesCuts(track)) {
946  if(fUsePtWeight) didacticSurface->Fill(track->Phi(), track->Eta(), track->Pt());
947  else didacticSurface->Fill(track->Phi(), track->Eta());
948  }
949  }
950  if(fExcludeLeadingJetsFromFit) { // visualize the excluded region
951  TF2 *f2 = new TF2(Form("%s_LJ", didacticSurface->GetName()),"[0]*TMath::Gaus(x,[1],[2])*TMath::Gaus(y,[3],[4])", 0, TMath::TwoPi(), -1, 1);
952  f2->SetParameters(excludeInPt/3.,excludeInPhi,.1,excludeInEta,.1);
953  didacticSurface->GetListOfFunctions()->Add(f2);
954  }
955  fOutputListGood->Add(didacticSurface);
956  } break;
957  default : break;
958  }
959  } else { // if the fit is of poor quality revert to the original rho estimate
961  switch (fRunModeType) { // again see if we want to save the fit
962  case kLocal : {
963  static Int_t didacticCounterWorst(0);
964  if(gRandom->Uniform(0, 100) > fPercentageOfFits) break;
965  TProfile* didacticProfile = (TProfile*)_tempSwap.Clone(Form("Fit_%i_1-CDF_%.3f_cen_%i_%s", didacticCounterWorst, CDF, fInCentralitySelection, detector.Data() ));
966  TF1* didactifFit = (TF1*)fFitModulation->Clone(Form("fit_%i_p_%.3f_cen_%i_%s", didacticCounterWorst, CDF, fInCentralitySelection, detector.Data()));
967  didacticProfile->GetListOfFunctions()->Add(didactifFit);
968  fOutputListBad->Add(didacticProfile);
969  didacticCounterWorst++;
970  } break;
971  default : break;
972  }
973  switch (fFitModulationType) {
974  case kNoFit : break; // nothing to do
975  case kCombined : fFitModulation->SetParameter(7, 0); // no break
976  case kFourierSeries : fFitModulation->SetParameter(7, 0); // no break
977  default : { // needs to be done if there was a poor fit
978  fFitModulation->SetParameter(3, 0);
979  fFitModulation->SetParameter(0, fLocalRho->GetVal());
980  } break;
981  }
982  return kFALSE; // return false if the fit is rejected
983  }
984  return kTRUE;
985 }
986 
987 //_____________________________________________________________________________
989 {
990  // Fill the analysis summary histrogram, saves all relevant analysis settigns
991  #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
992  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
993  #endif
994  fHistAnalysisSummary->GetXaxis()->SetBinLabel(2, "fJetRadius");
995  fHistAnalysisSummary->SetBinContent(2, GetJetContainer()->GetJetRadius());
996  fHistAnalysisSummary->GetXaxis()->SetBinLabel(3, "fJetEtaMin");
997  fHistAnalysisSummary->SetBinContent(3, GetJetContainer()->GetJetEtaMin());
998  fHistAnalysisSummary->GetXaxis()->SetBinLabel(4, "fJetEtaMax");
999  fHistAnalysisSummary->SetBinContent(4, GetJetContainer()->GetJetEtaMax());
1000  fHistAnalysisSummary->GetXaxis()->SetBinLabel(5, "fJetPhiMin");
1001  fHistAnalysisSummary->SetBinContent(5, GetJetContainer()->GetJetPhiMin());
1002  fHistAnalysisSummary->GetXaxis()->SetBinLabel(6, "fJetPhiMax");
1003  fHistAnalysisSummary->SetBinContent(6, GetJetContainer()->GetJetPhiMin());
1004  fHistAnalysisSummary->GetXaxis()->SetBinLabel(34, "fitModulationType");
1005  fHistAnalysisSummary->SetBinContent(34, (int)fFitModulationType);
1006  fHistAnalysisSummary->GetXaxis()->SetBinLabel(35, "runModeType");
1007  fHistAnalysisSummary->SetBinContent(35, (int)fRunModeType);
1008  fHistAnalysisSummary->GetXaxis()->SetBinLabel(37, "iterator");
1009  fHistAnalysisSummary->SetBinContent(37, 1.);
1010  fHistAnalysisSummary->GetXaxis()->SetBinLabel(38, "fMinPvalue");
1011  fHistAnalysisSummary->SetBinContent(38, fMinPvalue);
1012  fHistAnalysisSummary->GetXaxis()->SetBinLabel(39, "fMaxPvalue");
1013  fHistAnalysisSummary->SetBinContent(39, fMaxPvalue);
1014  fHistAnalysisSummary->GetXaxis()->SetBinLabel(40, "fExcludeLeadingJetsFromFit");
1016  fHistAnalysisSummary->GetXaxis()->SetBinLabel(41, "fRebinSwapHistoOnTheFly");
1017  fHistAnalysisSummary->SetBinContent(41, (int)fRebinSwapHistoOnTheFly);
1018  fHistAnalysisSummary->GetXaxis()->SetBinLabel(42, "fUsePtWeight");
1019  fHistAnalysisSummary->SetBinContent(42, (int)fUsePtWeight);
1020  fHistAnalysisSummary->GetXaxis()->SetBinLabel(45, "fLocalJetMinEta");
1021  fHistAnalysisSummary->SetBinContent(45,fLocalJetMinEta );
1022  fHistAnalysisSummary->GetXaxis()->SetBinLabel(46, "fLocalJetMaxEta");
1023  fHistAnalysisSummary->SetBinContent(46, fLocalJetMaxEta);
1024  fHistAnalysisSummary->GetXaxis()->SetBinLabel(47, "fLocalJetMinPhi");
1025  fHistAnalysisSummary->SetBinContent(47, fLocalJetMinPhi);
1026  fHistAnalysisSummary->GetXaxis()->SetBinLabel(48, "fLocalJetMaxPhi");
1027  fHistAnalysisSummary->SetBinContent(48, fLocalJetMaxPhi);
1028  fHistAnalysisSummary->GetXaxis()->SetBinLabel(49, "fSoftTrackMinPt");
1029  fHistAnalysisSummary->SetBinContent(49, fSoftTrackMinPt);
1030  fHistAnalysisSummary->GetXaxis()->SetBinLabel(50, "fSoftTrackMaxPt");
1031  fHistAnalysisSummary->SetBinContent(50, fSoftTrackMaxPt);
1032  fHistAnalysisSummary->GetXaxis()->SetBinLabel(51, "fUseScaledRho");
1033  fHistAnalysisSummary->SetBinContent(51, fUseScaledRho);
1034 }
1035 
1036 //_____________________________________________________________________________
1038 {
1039  // Fill event plane histograms
1040  #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
1041  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1042  #endif
1043  fHistPsi2[fInCentralitySelection]->Fill(psi2);
1044  fHistPsi3[fInCentralitySelection]->Fill(psi3);
1045 }
1046 
1047 //_____________________________________________________________________________
1049 {
1050  #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
1051  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1052  #endif
1053  // Terminate
1054 }
1055 
1056 //_____________________________________________________________________________
1058 {
1059  // Set function to fit modulation
1060  #ifdef ALIANALYSISTASKLOCALRHO_DEBUG_FLAG_0
1061  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1062  #endif
1063  if (fFitModulation) delete fFitModulation;
1064  fFitModulation = fit;
1065 }
Double_t QCnS(Int_t i, Int_t j)
void CalculateEventPlaneCombinedVZERO(Double_t *comb) const
virtual Bool_t Run()
Run function. This is the core function of the analysis and contains the user code. Therefore users have to implement this function.
void CalculateEventPlaneTPC(Double_t *tpc)
virtual void Terminate(Option_t *option)
double Double_t
Definition: External.C:58
Definition: External.C:236
const char * title
Definition: MakeQAPdf.C:26
Bool_t PassesCuts(AliVTrack *track) const
AliJetContainer * GetJetContainer(Int_t i=0) const
Bool_t QCnRecovery(Double_t psi2, Double_t psi3)
Bool_t fAbsVnHarmonics
status of rho vs centrality
Double_t Eta() const
Definition: AliEmcalJet.h:114
Double_t Phi() const
Definition: AliEmcalJet.h:110
Double_t ChiSquareCDF(Int_t ndf, Double_t x) const
Double_t GetJetEtaMax() const
TH1F * fHistAnalysisSummary
swap histogram
TCanvas * c
Definition: TestFitELoss.C:172
TH1F * fHistPsi2[10]
v3 cumulant
TRandom * gRandom
ClassImp(AliAnalysisTaskLocalRho) AliAnalysisTaskLocalRho
Double_t GetJetRadius(Int_t i=0) const
void CalculateEventPlaneVZERO(Double_t vzero[2][2]) const
TString fLocalRhoName
name for local rho
Bool_t CorrectRho(Double_t psi2, Double_t psi3)
Int_t fNAcceptedTracksQCn
number of accepted tracks
AliEmcalJet * GetLeadingJet(const char *opt="")
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
void SetJetEtaLimits(Float_t min, Float_t max, Int_t c=0)
AliRhoParameter * fRho
! event rho
TList * fOutputListBad
output list for local analysis
void QCnQnk(Int_t n, Int_t k, Double_t &reQ, Double_t &imQ)
void ExecOnce()
Perform steps needed to initialize the analysis.
Double_t fCent
!event centrality
AliLocalRhoParameter * fLocalRho
! local event rho
fitModulationType fFitModulationType
centrality bin, only for QA plots
TClonesArray * fJets
! jets
void SetJetPhiLimits(Float_t min, Float_t max, Int_t c=0)
Double_t Pt() const
Definition: AliEmcalJet.h:102
TList * fOutputListGood
output list
TH2F * fHistRhoStatusCent
cdf value of chisquare p
TProfile * fProfV3
v2 cumulant
Float_t GetJetRadius() const
void FillEventPlaneHistograms(Double_t psi2, Double_t psi3) const
TH2F * BookTH2F(const char *name, const char *x, const char *y, Int_t binsx, Double_t minx, Double_t maxx, Int_t binsy, Double_t miny, Double_t maxy, Int_t c=-1, Bool_t append=kTRUE)
TClonesArray * fTracks
!tracks
TProfile * fProfV3Cumulant
extracted v3
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
TH1F * fHistSwap
output list for local analysis
Int_t fInCentralitySelection
accepted tracks for QCn
Int_t GetRunNumber(TString)
Definition: PlotMuonQA.C:2235
const char Option_t
Definition: External.C:48
TProfile * fProfV2Cumulant
extracted v2
bool Bool_t
Definition: External.C:53
Double_t PhaseShift(Double_t x) const
TH1F * BookTH1F(const char *name, const char *x, Int_t bins, Double_t min, Double_t max, Int_t c=-1, Bool_t append=kTRUE)
Bool_t fAttachToEvent
is the analysis initialized?