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