PWGJE/EMCALJetTasks/AliAnalysisTaskLocalRhoDev.cxx

   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 #include <TArrayI.h>
  28 // aliroot includes
  29 #include <AliAnalysisTask.h>
  30 #include <AliAnalysisManager.h>
  31 #include <AliCentrality.h>
  32 #include <AliVVertex.h>
  33 #include <AliESDEvent.h>
  34 #include <AliAODEvent.h>
  35 #include <AliAODTrack.h>
  36 // emcal jet framework includes
  37 #include <AliPicoTrack.h>
  38 #include <AliEmcalJet.h>
  39 #include <AliRhoParameter.h>
  40 #include <AliLocalRhoParameter.h>
  41 #include <AliAnalysisTaskLocalRhoDev.h>
  42
  43 class AliAnalysisTaskLocalRhoDev;
  44 using namespace std;
  45
  46 ClassImp(AliAnalysisTaskLocalRhoDev)
  47
  48 //_____________________________________________________________________________
  49 AliAnalysisTaskLocalRhoDev::AliAnalysisTaskLocalRhoDev() :
  50   AliAnalysisTaskEmcalJet("AliAnalysisTaskLocalRhoDev", kTRUE),
  51   fDebug(0), fInitialized(0), fAttachToEvent(kTRUE), fFillHistograms(kFALSE), fNoEventWeightsForQC(kTRUE),
  52   fUseScaledRho(0), fCentralityClasses(0), fUserSuppliedV2(0), fUserSuppliedV3(0), fUserSuppliedR2(0),
  53   fUserSuppliedR3(0), fNAcceptedTracks(0), fNAcceptedTracksQCn(0), fInCentralitySelection(-1), fFitModulationType(kCombined),
  54   fFitGoodnessTest(kChi2Poisson), fQCRecovery(kTryFit), fUsePtWeight(kTRUE), fUsePtWeightErrorPropagation(kFALSE), fDetectorType(kTPC),
  55   fFitModulationOptions("WLQI"), fRunModeType(kGrid), fFitModulation(0), fFitControl(0x0), fMinPvalue(0.01), fMaxPvalue(1), fExpectedRuns(0x0),
  56   fExpectedSemiGoodRuns(0x0), fRunNumber(-1), fCachedRho(0x0), fNameSmallRho(""), fLocalJetMinEta(-10), fLocalJetMaxEta(-10),
  57   fLocalJetMinPhi(-10), fLocalJetMaxPhi(-10), fSoftTrackMinPt(0.15),
  58   fSoftTrackMaxPt(5.), fSemiGoodJetMinPhi(0.), fSemiGoodJetMaxPhi(4.), fSemiGoodTrackMinPhi(0.), fSemiGoodTrackMaxPhi(4.),
  59   fHistPvalueCDFROOT(0x0), fHistPvalueCDFROOTCent(0x0), fHistChi2ROOTCent(0x0), fHistPChi2Root(0x0), fHistPvalueCDF(0x0), fHistPvalueCDFCent(0x0),
  60   fHistChi2Cent(0x0), fHistPChi2(0x0), fHistRhoStatusCent(0), fAbsVnHarmonics(kTRUE), fExcludeLeadingJetsFromFit(1.),
  61   fRebinSwapHistoOnTheFly(kTRUE), fPercentageOfFits(10.), fUseV0EventPlaneFromHeader(kTRUE), fOutputList(0),
  62   fOutputListGood(0), fOutputListBad(0), fHistSwap(0), fHistAnalysisSummary(0), fProfV2(0), fProfV2Cumulant(0),
  63   fProfV3(0), fProfV3Cumulant(0)
  64 {
  65   // Default constructor
  66
  67   for(Int_t i(0); i < 10; i++) {
  68     fHistPsi2[i] = 0;
  69     fHistPsi3[i] = 0;
  70   }
  71 }
  72
  73 //_____________________________________________________________________________
  74 AliAnalysisTaskLocalRhoDev::AliAnalysisTaskLocalRhoDev(const char* name, runModeType type) :
  75   AliAnalysisTaskEmcalJet(name, kTRUE),
  76   fDebug(0), fInitialized(0), fAttachToEvent(kTRUE), fFillHistograms(kFALSE), fNoEventWeightsForQC(kTRUE),
  77   fUseScaledRho(0), fCentralityClasses(0), fUserSuppliedV2(0), fUserSuppliedV3(0), fUserSuppliedR2(0),
  78   fUserSuppliedR3(0), fNAcceptedTracks(0), fNAcceptedTracksQCn(0), fInCentralitySelection(-1), fFitModulationType(kCombined),
  79   fFitGoodnessTest(kChi2Poisson), fQCRecovery(kTryFit), fUsePtWeight(kTRUE), fUsePtWeightErrorPropagation(kFALSE), fDetectorType(kTPC),
  80   fFitModulationOptions("WLQI"), fRunModeType(type), fFitModulation(0), fFitControl(0x0), fMinPvalue(0.01), fMaxPvalue(1), fExpectedRuns(0x0),
  81   fExpectedSemiGoodRuns(0x0), fRunNumber(-1), fCachedRho(0x0), fNameSmallRho(""), fLocalJetMinEta(-10), fLocalJetMaxEta(-10),
  82   fLocalJetMinPhi(-10), fLocalJetMaxPhi(-10), fSoftTrackMinPt(0.15),
  83   fSoftTrackMaxPt(5.), fSemiGoodJetMinPhi(0.), fSemiGoodJetMaxPhi(4.), fSemiGoodTrackMinPhi(0.), fSemiGoodTrackMaxPhi(4.),
  84   fHistPvalueCDFROOT(0x0), fHistPvalueCDFROOTCent(0x0), fHistChi2ROOTCent(0x0), fHistPChi2Root(0x0), fHistPvalueCDF(0x0), fHistPvalueCDFCent(0x0),
  85   fHistChi2Cent(0x0), fHistPChi2(0x0), fHistRhoStatusCent(0), fAbsVnHarmonics(kTRUE), fExcludeLeadingJetsFromFit(1.),
  86   fRebinSwapHistoOnTheFly(kTRUE), fPercentageOfFits(10.), fUseV0EventPlaneFromHeader(kTRUE), fOutputList(0),
  87   fOutputListGood(0), fOutputListBad(0), fHistSwap(0), fHistAnalysisSummary(0), fProfV2(0), fProfV2Cumulant(0),
  88   fProfV3(0), fProfV3Cumulant(0)
  89 {
  90   // Constructor
  91   for(Int_t i(0); i < 10; i++) {
  92     fHistPsi2[i] = 0;
  93     fHistPsi3[i] = 0;
  94   }
  95
  96   DefineInput(0, TChain::Class());
  97   DefineOutput(1, TList::Class());
  98   switch (fRunModeType) {
  99   case kLocal : {
 100     gStyle->SetOptFit(1);
 101     DefineOutput(2, TList::Class());
 102     DefineOutput(3, TList::Class());
 103   } break;
 104   default: fDebug = -1;   // suppress debug info explicitely when not running locally
 105   }
 106 }
 107
 108 //_____________________________________________________________________________
 109 AliAnalysisTaskLocalRhoDev::~AliAnalysisTaskLocalRhoDev()
 110 {
 111   // destructor
 112   if(fOutputList)             delete fOutputList;
 113   if(fOutputListGood)         delete fOutputListGood;
 114   if(fOutputListBad)          delete fOutputListBad;
 115   if(fFitModulation)          delete fFitModulation;
 116   if(fHistSwap)               delete fHistSwap;
 117 }
 118
 119 //_____________________________________________________________________________
 120 void AliAnalysisTaskLocalRhoDev::ExecOnce()
 121 {
 122   // Init the analysis
 123   if(fLocalRhoName=="") fLocalRhoName = Form("LocalRhoFrom_%s", GetName());
 124   fLocalRho = new AliLocalRhoParameter(fLocalRhoName.Data(), 0);
 125   // add the local rho to the event if necessary
 126   if(fAttachToEvent) {
 127     if(!(InputEvent()->FindListObject(fLocalRho->GetName()))) {
 128       InputEvent()->AddObject(fLocalRho);
 129     } else {
 130       AliFatal(Form("%s: Container with same name %s already present. Aborting", GetName(), fLocalRho->GetName()));
 131     }
 132   }
 133   AliAnalysisTaskEmcalJet::ExecOnce();        // init the base clas
 134   if(fUseScaledRho) {
 135     // unscaled rho has been retrieved by the parent class, now we retrieve rho scaled
 136     fRho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(Form("%s_Scaled", fRho->GetName())));
 137     if(!fRho) {
 138       AliFatal(Form("%s: Couldn't find container for scaled rho. Aborting !", GetName()));
 139     }
 140   }
 141   if(!GetJetContainer()) AliFatal(Form("%s: Couldn't get jet container. Aborting !", GetName()));
 142 }
 143
 144 //_____________________________________________________________________________
 145 Bool_t AliAnalysisTaskLocalRhoDev::InitializeAnalysis()
 146 {
 147   // Initialize the anaysis
 148
 149   if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 150   if(fLocalJetMinEta > -10 && fLocalJetMaxEta > -10) SetJetEtaLimits(fLocalJetMinEta, fLocalJetMaxEta);
 151   if(fLocalJetMinPhi > -10 && fLocalJetMaxPhi > -10) SetJetPhiLimits(fLocalJetMinPhi, fLocalJetMaxPhi);
 152   switch (fFitModulationType)  {
 153   case kNoFit : { SetModulationFit(new TF1("fit_kNoFit", "[0]", 0, TMath::TwoPi())); } break;
 154   case kV2 : {
 155     SetModulationFit(new TF1("fit_kV2", "[0]*([1]+[2]*[3]*TMath::Cos([2]*(x-[4])))", 0, TMath::TwoPi()));
 156     fFitModulation->SetParameter(0, 0.);        // normalization
 157     fFitModulation->SetParameter(3, 0.2);       // v2
 158     fFitModulation->FixParameter(1, 1.);        // constant
 159     fFitModulation->FixParameter(2, 2.);        // constant
 160   } break;
 161   case kV3: {
 162     SetModulationFit(new TF1("fit_kV3", "[0]*([1]+[2]*[3]*TMath::Cos([2]*(x-[4])))", 0, TMath::TwoPi()));
 163     fFitModulation->SetParameter(0, 0.);        // normalization
 164     fFitModulation->SetParameter(3, 0.2);       // v3
 165     fFitModulation->FixParameter(1, 1.);        // constant
 166     fFitModulation->FixParameter(2, 3.);        // constant
 167   } break;
 168   default : { // for the combined fit, the 'direct fourier series' or the user supplied vn values we use v2 and v3
 169     SetModulationFit(new TF1("fit_kCombined", "[0]*([1]+[2]*([3]*TMath::Cos([2]*(x-[4]))+[7]*TMath::Cos([5]*(x-[6]))))", 0, TMath::TwoPi()));
 170     fFitModulation->SetParameter(0, 0.);       // normalization
 171     fFitModulation->SetParameter(3, 0.2);      // v2
 172     fFitModulation->FixParameter(1, 1.);       // constant
 173     fFitModulation->FixParameter(2, 2.);       // constant
 174     fFitModulation->FixParameter(5, 3.);       // constant
 175     fFitModulation->SetParameter(7, 0.2);      // v3
 176   } break;
 177   }
 178   switch (fRunModeType) {
 179   case kGrid : { fFitModulationOptions += "N0"; } break;
 180   default : break;
 181   }
 182   FillAnalysisSummaryHistogram();
 183   return kTRUE;
 184 }
 185
 186 //_____________________________________________________________________________
 187 void AliAnalysisTaskLocalRhoDev::UserCreateOutputObjects()
 188 {
 189   // create output objects
 190   if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 191   fHistSwap = new TH1F("fHistSwap", "fHistSwap", 20, 0, TMath::TwoPi());
 192   if(!fCentralityClasses) {   // classes must be defined at this point
 193     Double_t c[] = {0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100};
 194     fCentralityClasses = new TArrayD(sizeof(c)/sizeof(c[0]), c);
 195   }
 196   if(!fExpectedRuns) {        // expected runs must be defined at this point
 197     Int_t r[] =  {167813, 167988, 168066, 168068, 168069, 168076, 168104, 168212, 168311, 168322, 168325, 168341, 168361, 168362, 168458, 168460, 168461, 168992, 169091, 169094, 169138, 169143, 169167, 169417, 169835, 169837, 169838, 169846, 169855, 169858, 169859, 169923, 169956, 170027, 170036, 170081, /* up till here original good TPC list */169975, 169981, 170038, 170040, 170083, 170084, 170085, 170088, 170089, 170091, 170152, 170155, 170159, 170163, 170193, 170195, 170203, 170204, 170205, 170228, 170230, 170264, 170268, 170269, 170270, 170306, 170308, 170309, /* original semi-good tpc list */169415, 169411, 169035, 168988, 168984, 168826, 168777, 168512, 168511, 168467, 168464, 168342, 168310, 168115, 168108, 168107, 167987, 167915, 167903, /*new runs, good according to RCT */ 169238, 169160, 169156, 169148, 169145, 169144 /* run swith missing OROC 8 but seem ok in QA */};
 198     fExpectedRuns = new TArrayI(sizeof(r)/sizeof(r[0]), r);
 199   }
 200   if(!fExpectedSemiGoodRuns) {
 201     Int_t r[] = {169975, 169981, 170038, 170040, 170083, 170084, 170085, 170088, 170089, 170091, 170152, 170155, 170159, 170163, 170193, 170195, 170203, 170204, 170205, 170228, 170230, 170264, 170268, 170269, 170270, 170306, 170308, 170309};
 202     fExpectedSemiGoodRuns = new TArrayI(sizeof(r)/sizeof(r[0]), r);
 203   }
 204   fOutputList = new TList();
 205   fOutputList->SetOwner(kTRUE);
 206   // the analysis summary histo which stores all the analysis flags is always written to file
 207   fHistAnalysisSummary = BookTH1F("fHistAnalysisSummary", "flag", 54, -0.5, 54.5);
 208   if(!fFillHistograms) {
 209     PostData(1, fOutputList);
 210     return;
 211   }
 212   for(Int_t i(0); i < fCentralityClasses->GetSize()-1; i++) {
 213     fHistPsi2[i] = BookTH1F("fHistPsi2", "#Psi_{2}", 100, -.5*TMath::Pi(), .5*TMath::Pi(), i);
 214     fHistPsi3[i] = BookTH1F("fHistPsi3", "#Psi_{3}", 100, -1.*TMath::Pi()/3., TMath::Pi()/3., i);
 215   }
 216   // cdf of chisquare distribution
 217   fHistPvalueCDF = BookTH1F("fHistPvalueCDF", "CDF #chi^{2}", 500, 0, 1);
 218   fHistRhoStatusCent = BookTH2F("fHistRhoStatusCent", "centrality", "status [0=ok, 1=failed]", 101, -1, 100, 2, -.5, 1.5);
 219   // vn profiles
 220   Float_t temp[fCentralityClasses->GetSize()];
 221   for(Int_t i(0); i < fCentralityClasses->GetSize(); i++) temp[i] = fCentralityClasses->At(i);
 222   fProfV2 = new TProfile("fProfV2", "fProfV2", fCentralityClasses->GetSize()-1, temp);
 223   fProfV3 = new TProfile("fProfV3", "fProfV3", fCentralityClasses->GetSize()-1, temp);
 224   fOutputList->Add(fProfV2);
 225   fOutputList->Add(fProfV3);
 226   switch (fFitModulationType) {
 227   case kQC2 : {
 228     fProfV2Cumulant = new TProfile("fProfV2Cumulant", "fProfV2Cumulant", fCentralityClasses->GetSize()-1, temp);
 229     fProfV3Cumulant = new TProfile("fProfV3Cumulant", "fProfV3Cumulant", fCentralityClasses->GetSize()-1, temp);
 230     fOutputList->Add(fProfV2Cumulant);
 231     fOutputList->Add(fProfV3Cumulant);
 232   } break;
 233   case kQC4 : {
 234     fProfV2Cumulant = new TProfile("fProfV2Cumulant", "fProfV2Cumulant", fCentralityClasses->GetSize()-1, temp);
 235     fProfV3Cumulant = new TProfile("fProfV3Cumulant", "fProfV3Cumulant", fCentralityClasses->GetSize()-1, temp);
 236     fOutputList->Add(fProfV2Cumulant);
 237     fOutputList->Add(fProfV3Cumulant);
 238   } break;
 239   default : break;
 240   }
 241   if(fUsePtWeight) fHistSwap->Sumw2();
 242   if(fUserSuppliedV2) fOutputList->Add(fUserSuppliedV2);
 243   if(fUserSuppliedV3) fOutputList->Add(fUserSuppliedV3);
 244   if(fUserSuppliedR2) fOutputList->Add(fUserSuppliedR2);
 245   if(fUserSuppliedR3) fOutputList->Add(fUserSuppliedR3);
 246
 247   // increase readability of output list
 248   fOutputList->Sort();
 249
 250   // cdf and pdf of chisquare distribution
 251   fHistPvalueCDF = BookTH1F("fHistPvalueCDF", "CDF #chi^{2}", 50, 0, 1);
 252   fHistPvalueCDFCent = BookTH2F("fHistPvalueCDFCent", "centrality", "p-value", 40, 0, 100, 40, 0, 1);
 253   fHistChi2Cent = BookTH2F("fHistChi2Cent", "centrality", "#tilde{#chi^{2}}", 100, 0, 100, 100, 0, 5);
 254   fHistPChi2 = BookTH2F("fHistPChi2", "p-value", "#tilde{#chi^{2}}", 1000, 0, 1, 100, 0, 5);
 255   fHistPvalueCDFROOT = BookTH1F("fHistPvalueCDFROOT", "CDF #chi^{2} ROOT", 50, 0, 1);
 256   fHistPvalueCDFROOTCent = BookTH2F("fHistPvalueCDFROOTCent", "centrality", "p-value ROOT", 40, 0, 100, 45, 0, 1);
 257   fHistChi2ROOTCent = BookTH2F("fHistChi2ROOTCent", "centrality", "#tilde{#chi^{2}}", 40, 0, 100, 45, 0, 5);
 258   fHistPChi2Root = BookTH2F("fHistPChi2Root", "p-value", "#tilde{#chi^{2}} ROOT", 1000, 0, 1, 100, 0, 5);
 259   fHistRhoStatusCent = BookTH2F("fHistRhoStatusCent", "centrality", "status [-1=lin was better, 0=ok, 1 = failed]", 101, -1, 100, 3, -1.5, 1.5);
 260
 261   PostData(1, fOutputList);
 262   switch (fRunModeType) {
 263   case kLocal : {
 264     fOutputListGood = new TList();
 265     fOutputListGood->SetOwner(kTRUE);
 266     fOutputListBad = new TList();
 267     fOutputListBad->SetOwner(kTRUE);
 268     PostData(2, fOutputListGood);
 269     PostData(3, fOutputListBad);
 270   } break;
 271   default: break;
 272   }
 273 }
 274
 275 //_____________________________________________________________________________
 276 TH1F* AliAnalysisTaskLocalRhoDev::BookTH1F(const char* name, const char* x, Int_t bins, Double_t min, Double_t max, Int_t c, Bool_t append)
 277 {
 278   // Book a TH1F and connect it to the output container
 279
 280   if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 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", (int)fCentralityClasses->At(c), (int)fCentralityClasses->At(1+c));
 286   }
 287   title += Form(";%s;[counts]", x);
 288   TH1F* histogram = new TH1F(name, title.Data(), bins, min, max);
 289   histogram->Sumw2();
 290   if(append) fOutputList->Add(histogram);
 291   return histogram;
 292 }
 293
 294 //_____________________________________________________________________________
 295 TH2F* AliAnalysisTaskLocalRhoDev::BookTH2F(const char* name, const char* x, const char*y, Int_t binsx, Double_t minx, Double_t maxx,
 296                                         Int_t binsy, Double_t miny, Double_t maxy, Int_t c, Bool_t append)
 297 {
 298   // Book a TH2F and connect it to the output container
 299
 300   if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 301   if(!fOutputList) return 0x0;
 302   TString title(name);
 303   if(c!=-1) { // format centrality dependent histograms accordingly
 304     name = Form("%s_%i", name, c);
 305     title += Form("_%i-%i", (int)fCentralityClasses->At(c), (int)fCentralityClasses->At(1+c));
 306   }
 307   title += Form(";%s;%s", x, y);
 308   TH2F* histogram = new TH2F(name, title.Data(), binsx, minx, maxx, binsy, miny, maxy);
 309   histogram->Sumw2();
 310   if(append) fOutputList->Add(histogram);
 311   return histogram;
 312 }
 313
 314 //_____________________________________________________________________________
 315 Bool_t AliAnalysisTaskLocalRhoDev::Run()
 316 {
 317   // Execute once for each event
 318
 319   if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 320   if(!(PassesCuts(InputEvent())||fTracks||fJets||fRho)) return kFALSE;
 321   if(!fInitialized) fInitialized = InitializeAnalysis();
 322   // get the centrality bin (necessary for some control histograms
 323   fInCentralitySelection = -1;
 324   Double_t cent(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"));
 325   for(Int_t i(0); i < fCentralityClasses->GetSize()-1; i++) {
 326     if(cent >= fCentralityClasses->At(i) && cent <= fCentralityClasses->At(1+i)) {
 327       fInCentralitySelection = i;
 328       break; }
 329   }
 330   if(fInCentralitySelection < 0) return kFALSE;
 331   // set the rho value
 332   fLocalRho->SetVal(fRho->GetVal());
 333   // set the correct event plane accordign to the requested reference detector
 334   Double_t psi2(-1), psi3(-1);
 335   switch (fDetectorType) {    // determine the detector type for the rho fit
 336   case kTPC :     {
 337     // [0] psi2         [1] psi3
 338     Double_t tpc[2];
 339     CalculateEventPlaneTPC(tpc);
 340     psi2 = tpc[0];         psi3 = tpc[1];
 341   } break;
 342   case kVZEROA :  {
 343     // [0][0] psi2a     [1,0]   psi2c
 344     // [0][1] psi3a     [1,1]   psi3c
 345     Double_t vzero[2][2];
 346     CalculateEventPlaneVZERO(vzero);
 347     psi2 = vzero[0][0];    psi3 = vzero[0][1];
 348   }   break;
 349   case kVZEROC :  {
 350     // [0][0] psi2a     [1,0]   psi2c
 351     // [0][1] psi3a     [1,1]   psi3c
 352     Double_t vzero[2][2];
 353     CalculateEventPlaneVZERO(vzero);
 354     psi2 = vzero[1][0];    psi3 = vzero[1][1];
 355   }   break;
 356   case kVZEROComb : {
 357     /* for the combined vzero event plane
 358      * [0] psi2         [1] psi3
 359      * not fully implmemented yet, use with caution ! */
 360     Double_t vzeroComb[2];
 361     CalculateEventPlaneCombinedVZERO(vzeroComb);
 362     psi2 = vzeroComb[0]; psi3 = vzeroComb[1];
 363   } break;
 364   default : break;
 365   }
 366   if(fFillHistograms) FillEventPlaneHistograms(psi2, psi3);
 367   switch (fFitModulationType) { // do the fits
 368   case kNoFit : { fFitModulation->FixParameter(0, fLocalRho->GetVal()); } break;
 369   case kV2 : {    // only v2
 370     if(CorrectRho(psi2, psi3)) {
 371       if(fFillHistograms) fProfV2->Fill(fCent, fFitModulation->GetParameter(3));
 372       if(fUserSuppliedR2) {
 373         Double_t r(fUserSuppliedR2->GetBinContent(fUserSuppliedR2->GetXaxis()->FindBin(fCent)));
 374         if(r > 0) fFitModulation->SetParameter(3, fFitModulation->GetParameter(3)/r);
 375       }
 376     }
 377   } break;
 378   case kV3 : {    // only v3
 379     if(CorrectRho(psi2, psi3)) {
 380       if(fUserSuppliedR3) {
 381         Double_t r(fUserSuppliedR3->GetBinContent(fUserSuppliedR3->GetXaxis()->FindBin(fCent)));
 382         if(r > 0) fFitModulation->SetParameter(3, fFitModulation->GetParameter(3)/r);
 383       }
 384       if(fFillHistograms) fProfV3->Fill(fCent, fFitModulation->GetParameter(3));
 385     }
 386   } break;
 387   case kQC2 : {   // qc2 analysis - NOTE: not a wise idea to use this !
 388     if(CorrectRho(psi2, psi3)) {
 389       if(fUserSuppliedR2 && fUserSuppliedR3) {
 390         // note for the qc method, resolution is REVERSED to go back to v2obs
 391         Double_t r2(fUserSuppliedR2->GetBinContent(fUserSuppliedR2->GetXaxis()->FindBin(fCent)));
 392         Double_t r3(fUserSuppliedR3->GetBinContent(fUserSuppliedR3->GetXaxis()->FindBin(fCent)));
 393         if(r2 > 0) fFitModulation->SetParameter(3, fFitModulation->GetParameter(3)*r2);
 394         if(r3 > 0) fFitModulation->SetParameter(7, fFitModulation->GetParameter(7)*r3);
 395       }
 396       if (fUsePtWeight) { // use weighted weights
 397         Double_t dQCnM11 = (fNoEventWeightsForQC) ? 1. : QCnM11();
 398         if(fFillHistograms) {
 399           fProfV2->Fill(fCent, fFitModulation->GetParameter(3), dQCnM11);
 400           fProfV3->Fill(fCent, fFitModulation->GetParameter(7), dQCnM11);
 401         }
 402       } else {
 403         Double_t dQCnM = (fNoEventWeightsForQC) ? 2. : QCnM();
 404         if(fFillHistograms) {
 405           fProfV2->Fill(fCent, fFitModulation->GetParameter(3), dQCnM*(dQCnM-1));
 406           fProfV3->Fill(fCent, fFitModulation->GetParameter(7), dQCnM*(dQCnM-1));
 407         }
 408       }
 409     }
 410   } break;
 411   case kQC4 : {   // NOTE: see comment at kQC2
 412     if(CorrectRho(psi2, psi3)) {
 413       if(fUserSuppliedR2 && fUserSuppliedR3) {
 414         // note for the qc method, resolution is REVERSED to go back to v2obs
 415         Double_t r2(fUserSuppliedR2->GetBinContent(fUserSuppliedR2->GetXaxis()->FindBin(fCent)));
 416         Double_t r3(fUserSuppliedR3->GetBinContent(fUserSuppliedR3->GetXaxis()->FindBin(fCent)));
 417         if(r2 > 0) fFitModulation->SetParameter(3, fFitModulation->GetParameter(3)*r2);
 418         if(r3 > 0) fFitModulation->SetParameter(7, fFitModulation->GetParameter(7)*r3);
 419       }
 420       if (fUsePtWeight) { // use weighted weights
 421         if(fFillHistograms) {
 422           fProfV2->Fill(fCent, TMath::Power(fFitModulation->GetParameter(3),0.5)/*, QCnM1111()*/);
 423           fProfV3->Fill(fCent, TMath::Power(fFitModulation->GetParameter(7),0.5)/*, QCnM1111()*/);
 424         }
 425       } else {
 426         if(fFillHistograms) {
 427           fProfV2->Fill(fCent, TMath::Power(fFitModulation->GetParameter(3),0.5)/*, QCnM()*(QCnM()-1)*(QCnM()-2)*(QCnM()-3)*/);
 428           fProfV3->Fill(fCent, TMath::Power(fFitModulation->GetParameter(7),0.5)/*, QCnM()*(QCnM()-1)*(QCnM()-2)*(QCnM()-3)*/);
 429         }
 430       }
 431     }
 432   } break;
 433   default : {
 434     if(CorrectRho(psi2, psi3)) {
 435       if(fUserSuppliedR2 && fUserSuppliedR3) {
 436         Double_t r2(fUserSuppliedR2->GetBinContent(fUserSuppliedR2->GetXaxis()->FindBin(fCent)));
 437         Double_t r3(fUserSuppliedR3->GetBinContent(fUserSuppliedR3->GetXaxis()->FindBin(fCent)));
 438         if(r2 > 0) fFitModulation->SetParameter(3, fFitModulation->GetParameter(3)/r2);
 439         if(r3 > 0) fFitModulation->SetParameter(7, fFitModulation->GetParameter(7)/r3);
 440       }
 441       if(fFillHistograms) {
 442         fProfV2->Fill(fCent, fFitModulation->GetParameter(3));
 443         fProfV3->Fill(fCent, fFitModulation->GetParameter(7));
 444       }
 445     }
 446   } break;
 447   }
 448   // if all went well, add local rho
 449   fLocalRho->SetLocalRho(fFitModulation);
 450   PostData(1, fOutputList);
 451   return kTRUE;
 452 }
 453
 454 //_____________________________________________________________________________
 455 void AliAnalysisTaskLocalRhoDev::CalculateEventPlaneVZERO(Double_t vzero[2][2]) const
 456 {
 457   // Get the vzero event plane
 458   if(fUseV0EventPlaneFromHeader) {
 459     // use the vzero event plane from the event header
 460     // note: to use the calibrated vzero event plane, run
 461     // $ALICE_ROOT/ANALYSIS/macros/AddTaskVZEROEPSelection.C
 462     // prior to this task (make sure the calibration is available for the dataset
 463     // you want to use)
 464     Double_t a(0), b(0), c(0), d(0), e(0), f(0), g(0), h(0);
 465     vzero[0][0] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 2, a, b);
 466     vzero[1][0] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 2, c, d);
 467     vzero[0][1] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 3, e, f);
 468     vzero[1][1] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 3, g, h);
 469     return;
 470   }
 471   // grab the vzero event plane without recentering
 472   if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 473   Double_t qxa2(0), qya2(0), qxc2(0), qyc2(0);    // for psi2
 474   Double_t qxa3(0), qya3(0), qxc3(0), qyc3(0);    // for psi3
 475   for(Int_t iVZERO(0); iVZERO < 64; iVZERO++) {
 476     Double_t phi(TMath::PiOver4()*(.5+iVZERO%8)), /* eta(0), */ weight(InputEvent()->GetVZEROEqMultiplicity(iVZERO));
 477     //        (iVZERO<32) ? eta = -3.45+.5*(iVZERO/8) : eta = 4.8-.6*((iVZERO/8)-4);
 478     if(iVZERO<32) {
 479       qxa2 += weight*TMath::Cos(2.*phi);
 480       qya2 += weight*TMath::Sin(2.*phi);
 481       qxa3 += weight*TMath::Cos(3.*phi);
 482       qya3 += weight*TMath::Sin(3.*phi);
 483     }
 484     else {
 485       qxc2 += weight*TMath::Cos(2.*phi);
 486       qyc2 += weight*TMath::Sin(2.*phi);
 487       qxc3 += weight*TMath::Cos(3.*phi);
 488       qyc3 += weight*TMath::Sin(3.*phi);
 489     }
 490   }
 491   vzero[0][0] = .5*TMath::ATan2(qya2, qxa2);
 492   vzero[1][0] = .5*TMath::ATan2(qyc2, qxc2);
 493   vzero[0][1] = (1./3.)*TMath::ATan2(qya3, qxa3);
 494   vzero[1][1] = (1./3.)*TMath::ATan2(qyc3, qxc3);
 495 }
 496
 497 //_____________________________________________________________________________
 498 void AliAnalysisTaskLocalRhoDev::CalculateEventPlaneTPC(Double_t* tpc)
 499 {
 500   // Grab the TPC event plane. if parameter fExcludeLeadingJetsFromFit is larger than 0,
 501   // strip in eta of width fExcludeLeadingJetsFromFit * GetJetContainer()->GetJetRadius() around the leading jet (before
 502   // subtraction of rho) will be exluded from the event plane estimate
 503
 504   if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 505   fNAcceptedTracks = 0;                // reset the track counter
 506   Double_t qx2(0), qy2(0);     // for psi2
 507   Double_t qx3(0), qy3(0);     // for psi3
 508   if(fTracks) {
 509     Float_t excludeInEta = -999;
 510     if(fExcludeLeadingJetsFromFit > 0 ) {    // remove the leading jet from ep estimate
 511       AliEmcalJet* leadingJet(GetJetContainer()->GetLeadingJet());
 512       if(leadingJet) excludeInEta = leadingJet->Eta();
 513     }
 514     Int_t iTracks(fTracks->GetEntriesFast());
 515     for(Int_t iTPC(0); iTPC < iTracks; iTPC++) {
 516       AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(iTPC));
 517       if(!PassesCuts(track) || track->Pt() < fSoftTrackMinPt || track->Pt() > fSoftTrackMaxPt) continue;
 518       if(fExcludeLeadingJetsFromFit > 0 &&( (TMath::Abs(track->Eta() - excludeInEta) < GetJetContainer()->GetJetRadius()*fExcludeLeadingJetsFromFit ) || (TMath::Abs(track->Eta()) - GetJetContainer()->GetJetRadius() - GetJetContainer()->GetJetEtaMax() ) > 0 )) continue;
 519       fNAcceptedTracks++;
 520       qx2+= TMath::Cos(2.*track->Phi());
 521       qy2+= TMath::Sin(2.*track->Phi());
 522       qx3+= TMath::Cos(3.*track->Phi());
 523       qy3+= TMath::Sin(3.*track->Phi());
 524     }
 525   }
 526   tpc[0] = .5*TMath::ATan2(qy2, qx2);
 527   tpc[1] = (1./3.)*TMath::ATan2(qy3, qx3);
 528 }
 529
 530 //_____________________________________________________________________________
 531 void AliAnalysisTaskLocalRhoDev::CalculateEventPlaneCombinedVZERO(Double_t* comb) const
 532 {
 533   // Grab the combined vzero event plane
 534
 535   //    if(fUseV0EventPlaneFromHeader) {    // use the vzero from the header
 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   // FIXME the rest of this function isn't impelmented yet (as of 01-07-2013)
 540   // this means a default the combined vzero event plane from the header is used
 541   // to get this value 'by hand', vzeroa and vzeroc event planes have to be combined
 542   // according to their resolution - this will be added ...
 543   //
 544   //    } else {
 545   //        Double_t qx2a(0), qy2a(0), qx2c(0), qy2c(0), qx3a(0), qy3a(0), qx3c(0), qy3c(0);
 546   //        InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 2, qx2a, qy2a);
 547   //        InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 2, qx2c, qy2c);
 548   //        InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 3, qx3a, qy3a);
 549   //        InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 3, qx3c, qy3c);
 550   //        Double_t chi2A(-1), chi2C(-1), chi3A(-1), chi3C(-1);     // get chi from the resolution
 551   //        Double_t qx2(chi2A*chi2A*qx2a+chi2C*chi2C*qx2c);
 552   //        Double_t qy2(chi2A*chi2A*qy2a+chi2C*chi2C*qy2c);
 553   //        Double_t qx3(chi3A*chi3A*qx3a+chi3C*chi3C*qx3c);
 554   //        Double_t qy3(chi3A*chi3A*qy3a+chi3C*chi3C*qy3c);
 555   //        comb[0] = .5*TMath::ATan2(qy2, qx2);
 556   //        comb[1] = (1./3.)*TMath::ATan2(qy3, qx3);
 557   //    }
 558 }
 559
 560 //_____________________________________________________________________________
 561 Double_t AliAnalysisTaskLocalRhoDev::CalculateQC2(Int_t harm)
 562 {
 563   // Get the second order q-cumulant, a -999 return will be caught in the qa routine of CorrectRho
 564
 565   if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 566   Double_t reQ(0), imQ(0), modQ(0), M11(0), M(0);
 567   if(fUsePtWeight) {  // for the weighted 2-nd order q-cumulant
 568     QCnQnk(harm, 1, reQ, imQ);      // get the weighted 2-nd order q-vectors
 569     modQ = reQ*reQ+imQ*imQ;         // get abs Q-squared
 570     M11 = QCnM11();                 // equals S2,1 - S1,2
 571     return (M11 > 0) ? ((modQ - QCnS(1,2))/M11) : -999;
 572   } // else return the non-weighted 2-nd order q-cumulant
 573   QCnQnk(harm, 0, reQ, imQ);          // get the non-weighted 2-nd order q-vectors
 574   modQ = reQ*reQ+imQ*imQ;             // get abs Q-squared
 575   M = QCnM();
 576   return (M > 1) ? (modQ - M)/(M*(M-1)) : -999;
 577 }
 578
 579 //_____________________________________________________________________________
 580 Double_t AliAnalysisTaskLocalRhoDev::CalculateQC4(Int_t harm)
 581 {
 582   // Get the fourth order q-cumulant, a -999 return will be caught in the qa routine of CorrectRho
 583
 584   if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 585   Double_t reQn1(0), imQn1(0), reQ2n2(0), imQ2n2(0), reQn3(0), imQn3(0), M1111(0), M(0);
 586   Double_t a(0), b(0), c(0), d(0), e(0), f(0), g(0);  // terms of the calculation
 587   if(fUsePtWeight) {  // for the weighted 4-th order q-cumulant
 588     QCnQnk(harm, 1, reQn1, imQn1);
 589     QCnQnk(harm*2, 2, reQ2n2, imQ2n2);
 590     QCnQnk(harm, 3, reQn3, imQn3);
 591     // fill in the terms ...
 592     a = (reQn1*reQn1+imQn1*imQn1)*(reQn1*reQn1+imQn1*imQn1);
 593     b = reQ2n2*reQ2n2 + imQ2n2*imQ2n2;
 594     c = -2.*(reQ2n2*reQn1*reQn1-reQ2n2*imQn1*imQn1+2.*imQ2n2*reQn1*imQn1);
 595     d = 8.*(reQn3*reQn1+imQn3*imQn1);
 596     e = -4.*QCnS(1,2)*(reQn1*reQn1+imQn1*imQn1);
 597     f = -6.*QCnS(1,4);
 598     g = 2.*QCnS(2,2);
 599     M1111 = QCnM1111();
 600     return (M1111 > 0) ? (a+b+c+d+e+f+g)/M1111 : -999;
 601   }   // else return the unweighted case
 602   Double_t reQn(0), imQn(0), reQ2n(0), imQ2n(0);
 603   QCnQnk(harm, 0, reQn, imQn);
 604   QCnQnk(harm*2, 0, reQ2n, imQ2n);
 605   // fill in the terms ...
 606   M = QCnM();
 607   if(M < 4) return -999;
 608   a = (reQn*reQn+imQn*imQn)*(reQn*reQn+imQn*imQn);
 609   b = reQ2n*reQ2n + imQ2n*imQ2n;
 610   c = -2.*(reQ2n*reQn*reQn-reQ2n*imQn*imQn+2.*imQ2n*reQn*imQn);
 611   e = -4.*(M-2)*(reQn*reQn+imQn*imQn);
 612   f = 2.*M*(M-3);
 613   return (a+b+c+e+f)/(M*(M-1)*(M-2)*(M-3));
 614 }
 615
 616 //_____________________________________________________________________________
 617 void AliAnalysisTaskLocalRhoDev::QCnQnk(Int_t n, Int_t k, Double_t &reQ, Double_t &imQ)
 618 {
 619   // Get the weighted n-th order q-vector, pass real and imaginary part as reference
 620
 621   if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 622   if(!fTracks) return;
 623   fNAcceptedTracksQCn = 0;
 624   Int_t iTracks(fTracks->GetEntriesFast());
 625   for(Int_t iTPC(0); iTPC < iTracks; iTPC++) {
 626     AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(iTPC));
 627     if(!PassesCuts(track) || track->Pt() < fSoftTrackMinPt || track->Pt() > fSoftTrackMaxPt) continue;
 628     fNAcceptedTracksQCn++;
 629     // for the unweighted case, k equals zero and the weight doesn't contribute to the equation below
 630     reQ += TMath::Power(track->Pt(), k) * TMath::Cos(((double)n)*track->Phi());
 631     imQ += TMath::Power(track->Pt(), k) * TMath::Sin(((double)n)*track->Phi());
 632   }
 633 }
 634
 635 //_____________________________________________________________________________
 636 Double_t AliAnalysisTaskLocalRhoDev::QCnS(Int_t i, Int_t j)
 637 {
 638   // Get the weighted ij-th order autocorrelation correction
 639
 640   if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 641   if(!fTracks || i <= 0 || j <= 0) return -999;
 642   Int_t iTracks(fTracks->GetEntriesFast());
 643   Double_t Sij(0);
 644   for(Int_t iTPC(0); iTPC < iTracks; iTPC++) {
 645     AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(iTPC));
 646     if(!PassesCuts(track) || track->Pt() < fSoftTrackMinPt || track->Pt() > fSoftTrackMaxPt) continue;
 647     Sij+=TMath::Power(track->Pt(), j);
 648   }
 649   return TMath::Power(Sij, i);
 650 }
 651
 652 //_____________________________________________________________________________
 653 Double_t AliAnalysisTaskLocalRhoDev::QCnM()
 654 {
 655   // Get multiplicity for unweighted q-cumulants. function QCnQnk should be called first
 656
 657   if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 658   return (Double_t) fNAcceptedTracksQCn;
 659 }
 660
 661 //_____________________________________________________________________________
 662 Double_t AliAnalysisTaskLocalRhoDev::QCnM11()
 663 {
 664   // Get multiplicity weights for the weighted two particle cumulant
 665
 666   if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 667   return (QCnS(2,1) - QCnS(1,2));
 668 }
 669
 670 //_____________________________________________________________________________
 671 Double_t AliAnalysisTaskLocalRhoDev::QCnM1111()
 672 {
 673   // Get multiplicity weights for the weighted four particle cumulant
 674
 675   if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 676   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));
 677 }
 678
 679 //_____________________________________________________________________________
 680 Bool_t AliAnalysisTaskLocalRhoDev::QCnRecovery(Double_t psi2, Double_t psi3)
 681 {
 682   // Decides how to deal with the situation where c2 or c3 is negative
 683   // Returns kTRUE depending on whether or not a modulated rho is used for the jet background
 684
 685   if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 686   if(TMath::AreEqualAbs(fFitModulation->GetParameter(3), .0, 1e-10) && TMath::AreEqualAbs(fFitModulation->GetParameter(7), .0,1e-10)) {
 687     fFitModulation->SetParameter(7, 0);
 688     fFitModulation->SetParameter(3, 0);
 689     fFitModulation->SetParameter(0, fLocalRho->GetVal());
 690     return kTRUE;   // v2 and v3 have physical null values
 691   }
 692   switch (fQCRecovery) {
 693   case kFixedRho : {      // roll back to the original rho
 694     fFitModulation->SetParameter(7, 0);
 695     fFitModulation->SetParameter(3, 0);
 696     fFitModulation->SetParameter(0, fLocalRho->GetVal());
 697     return kFALSE;       // rho is forced to be fixed
 698   }
 699   case kNegativeVn : {
 700     Double_t c2(fFitModulation->GetParameter(3));
 701     Double_t c3(fFitModulation->GetParameter(7));
 702     if( c2 < 0 ) c2 = -1.*TMath::Sqrt(-1.*c2);
 703     if( c3 < 0 ) c3 = -1.*TMath::Sqrt(-1.*c3);
 704     fFitModulation->SetParameter(3, c2);
 705     fFitModulation->SetParameter(7, c3);
 706     return kTRUE;        // is this a physical quantity ?
 707   }
 708   case kTryFit : {
 709     fitModulationType tempType(fFitModulationType);  // store temporarily
 710     fFitModulationType = kCombined;
 711     fFitModulation->SetParameter(7, 0);
 712     fFitModulation->SetParameter(3, 0);
 713     Bool_t pass(CorrectRho(psi2, psi3));         // do the fit and all quality checks
 714     fFitModulationType = tempType;               // roll back for next event
 715     return pass;
 716   }
 717   default : return kFALSE;
 718   }
 719   return kFALSE;
 720 }
 721
 722 //_____________________________________________________________________________
 723 Bool_t AliAnalysisTaskLocalRhoDev::CorrectRho(Double_t psi2, Double_t psi3)
 724 {
 725   // Get rho' -> rho(phi)
 726   // three routines are available, 1 and 2 can be used with or without pt weights
 727   //  [1] get vn from q-cumulants
 728   //      in case of cumulants, both cumulants and vn values are stored. in both cases, v2 and v3
 729   //      are expected. a check is performed to see if rho has no negative local minimum
 730   //      for full description, see Phys. Rev. C 83, 044913
 731   //      since the cn distribution has negative values, vn = sqrt(cn) can be imaginary sometimes
 732   //      in this case one can either roll back to the 'original' fixed rho, do a fit for vn or take use
 733   //      vn = - sqrt(|cn|) note that because of this, use of q-cumulants is not safe !
 734   //  [2] fitting a fourier expansion to the de/dphi distribution
 735   //      the fit can be done with either v2, v3 or a combination.
 736   //      in all cases, a cut can be made on the p-value of the chi-squared value of the fit
 737   //      and a check can be performed to see if rho has no negative local minimum
 738   //  [3] get v2 and v3 from user supplied histograms
 739   //      in this way, a fixed value of v2 and v3 is subtracted w.r.t. whichever event plane is requested
 740
 741   if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 742   Int_t freeParams(2);                // free parameters of the fit (for NDF)
 743   switch (fFitModulationType) {       // for approaches where no fitting is required
 744     case kQC2 : {
 745       fFitModulation->FixParameter(4, psi2);
 746       fFitModulation->FixParameter(6, psi3);
 747       fFitModulation->FixParameter(3, CalculateQC2(2));   // set here with cn, vn = sqrt(cn)
 748       fFitModulation->FixParameter(7, CalculateQC2(3));
 749       // first fill the histos of the raw cumulant distribution
 750       if (fUsePtWeight) { // use weighted weights
 751         Double_t dQCnM11 = (fNoEventWeightsForQC) ? 1. : QCnM11();
 752         fProfV2Cumulant->Fill(fCent, fFitModulation->GetParameter(3), dQCnM11);
 753         fProfV3Cumulant->Fill(fCent, fFitModulation->GetParameter(7), dQCnM11);
 754       } else {
 755         Double_t dQCnM = (fNoEventWeightsForQC) ? 2. : QCnM();
 756         fProfV2Cumulant->Fill(fCent, fFitModulation->GetParameter(3), dQCnM*(dQCnM-1));
 757         fProfV3Cumulant->Fill(fCent, fFitModulation->GetParameter(7), dQCnM*(dQCnM-1));
 758       }
 759       // then see if one of the cn value is larger than zero and vn is readily available
 760       if(fFitModulation->GetParameter(3) > 0 && fFitModulation->GetParameter(7) > 0) {
 761         fFitModulation->FixParameter(3, TMath::Sqrt(fFitModulation->GetParameter(3)));
 762         fFitModulation->FixParameter(7, TMath::Sqrt(fFitModulation->GetParameter(7)));
 763       } else if (!QCnRecovery(psi2, psi3)) return kFALSE;  // try to recover the cumulant, this will set v2 and v3
 764       if(fAbsVnHarmonics && fFitModulation->GetMinimum(0, TMath::TwoPi()) < 0) {  // general check
 765         fFitModulation->SetParameter(7, 0);
 766         fFitModulation->SetParameter(3, 0);
 767         fFitModulation->SetParameter(0, fLocalRho->GetVal());
 768         return kFALSE;
 769       }
 770       return kTRUE;
 771     } break;
 772     case kQC4 : {
 773       fFitModulation->FixParameter(4, psi2);
 774       fFitModulation->FixParameter(6, psi3);
 775       fFitModulation->FixParameter(3, CalculateQC4(2));   // set here with cn, vn = sqrt(cn)
 776       fFitModulation->FixParameter(7, CalculateQC4(3));
 777       // first fill the histos of the raw cumulant distribution
 778       if (fUsePtWeight) { // use weighted weights
 779         fProfV2Cumulant->Fill(fCent, fFitModulation->GetParameter(3)/*, QCnM1111()*/);
 780         fProfV3Cumulant->Fill(fCent, fFitModulation->GetParameter(7)/*, QCnM1111()*/);
 781       } else {
 782         fProfV2Cumulant->Fill(fCent, fFitModulation->GetParameter(3)/*, QCnM1111()*/);
 783         fProfV3Cumulant->Fill(fCent, fFitModulation->GetParameter(7)/*, QCnM1111()*/);
 784       }
 785       // then see if one of the cn value is larger than zero and vn is readily available
 786       if(fFitModulation->GetParameter(3) > 0 && fFitModulation->GetParameter(7) > 0) {
 787         fFitModulation->FixParameter(3, TMath::Sqrt(fFitModulation->GetParameter(3)));
 788         fFitModulation->FixParameter(7, TMath::Sqrt(fFitModulation->GetParameter(7)));
 789       } else if (!QCnRecovery(psi2, psi3)) return kFALSE;  // try to recover the cumulant, this will set v2 and v3
 790       if(fAbsVnHarmonics && fFitModulation->GetMinimum(0, TMath::TwoPi()) < 0) {  // general check
 791         fFitModulation->SetParameter(7, 0);
 792         fFitModulation->SetParameter(3, 0);
 793         fFitModulation->SetParameter(0, fLocalRho->GetVal());
 794         return kFALSE;
 795       }
 796     } break;
 797     case kIntegratedFlow : {
 798       // use v2 and v3 values from an earlier iteration over the data
 799       fFitModulation->FixParameter(3, fUserSuppliedV2->GetBinContent(fUserSuppliedV2->GetXaxis()->FindBin(fCent)));
 800       fFitModulation->FixParameter(4, psi2);
 801       fFitModulation->FixParameter(6, psi3);
 802       fFitModulation->FixParameter(7, fUserSuppliedV3->GetBinContent(fUserSuppliedV3->GetXaxis()->FindBin(fCent)));
 803       if(fAbsVnHarmonics && fFitModulation->GetMinimum(0, TMath::TwoPi()) < 0) {
 804         fFitModulation->SetParameter(7, 0);
 805         fFitModulation->SetParameter(3, 0);
 806         fFitModulation->SetParameter(0, fLocalRho->GetVal());
 807           return kFALSE;
 808       }
 809       return kTRUE;
 810     }
 811     default : break;
 812   }
 813   TString detector("");
 814   switch (fDetectorType) {
 815     case kTPC : detector+="TPC";
 816       break;
 817     case kVZEROA : detector+="VZEROA";
 818       break;
 819     case kVZEROC : detector+="VZEROC";
 820       break;
 821     case kVZEROComb : detector+="VZEROComb";
 822       break;
 823     default: break;
 824   }
 825   Int_t iTracks(fTracks->GetEntriesFast());
 826   Double_t excludeInEta = -999;
 827   Double_t excludeInPhi = -999;
 828   Double_t excludeInPt  = -999;
 829   if(iTracks <= 0 || fLocalRho->GetVal() <= 0 ) return kFALSE;   // no use fitting an empty event ...
 830   if(fExcludeLeadingJetsFromFit > 0 ) {
 831     AliEmcalJet* leadingJet(GetJetContainer()->GetLeadingJet());
 832     if(PassesCuts(leadingJet)) {
 833       excludeInEta = leadingJet->Eta();
 834       excludeInPhi = leadingJet->Phi();
 835       excludeInPt = leadingJet->Pt();
 836     }
 837   }
 838   // check the acceptance of the track selection that will be used
 839   // if one uses e.g. semi-good tpc tracks, accepance in phi is reduced to 0 < phi < 4
 840   // the defaults (-10 < phi < 10) which accept all, are then overwritten
 841   Double_t lowBound(0.), upBound(TMath::TwoPi());     // bounds for fit
 842   if(GetParticleContainer()->GetParticlePhiMin() > lowBound) lowBound = GetParticleContainer()->GetParticlePhiMin();
 843   if(GetParticleContainer()->GetParticlePhiMax() < upBound) upBound = GetParticleContainer()->GetParticlePhiMax();
 844
 845   fHistSwap->Reset();                 // clear the histogram
 846   TH1F _tempSwap;     // on stack for quick access
 847   TH1F _tempSwapN;    // on stack for quick access, bookkeeping histogram
 848   if(fRebinSwapHistoOnTheFly) {
 849     if(fNAcceptedTracks < 49) fNAcceptedTracks = 49;       // avoid aliasing effects
 850     _tempSwap = TH1F("_tempSwap", "_tempSwap", TMath::CeilNint(TMath::Sqrt(fNAcceptedTracks)), lowBound, upBound);
 851     if(fUsePtWeightErrorPropagation) _tempSwapN = TH1F("_tempSwapN", "_tempSwapN", TMath::CeilNint(TMath::Sqrt(fNAcceptedTracks)), lowBound, upBound);
 852     if(fUsePtWeight) _tempSwap.Sumw2();
 853   }
 854   else _tempSwap = *fHistSwap;         // now _tempSwap holds the desired histo
 855   // non poissonian error when using pt weights
 856   Double_t totalpts(0.), totalptsquares(0.), totalns(0.);
 857   for(Int_t i(0); i < iTracks; i++) {
 858     AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(i));
 859     if(fExcludeLeadingJetsFromFit > 0 &&( (TMath::Abs(track->Eta() - excludeInEta) < GetJetContainer()->GetJetRadius()*fExcludeLeadingJetsFromFit ) || (TMath::Abs(track->Eta()) - GetJetContainer()->GetJetRadius() - GetJetContainer()->GetJetEtaMax() ) > 0 )) continue;
 860     if(!PassesCuts(track) || track->Pt() > fSoftTrackMaxPt || track->Pt() < fSoftTrackMinPt) continue;
 861      if(fUsePtWeight) {
 862        _tempSwap.Fill(track->Phi(), track->Pt());
 863        if(fUsePtWeightErrorPropagation) {
 864          totalpts += track->Pt();
 865          totalptsquares += track->Pt()*track->Pt();
 866          totalns += 1;
 867          _tempSwapN.Fill(track->Phi());
 868        }
 869      }
 870      else _tempSwap.Fill(track->Phi());
 871   }
 872   if(fUsePtWeight && fUsePtWeightErrorPropagation) {
 873     // in the case of pt weights overwrite the poissonian error estimate which is assigned by root by a more sophisticated appraoch
 874     // 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
 875     // 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
 876     // event, for the latter use a poissonian estimate. the two contrubitions are assumed to be uncorrelated
 877     if(totalns < 1) return kFALSE; // not one track passes the cuts
 878     for(Int_t l = 0; l < _tempSwap.GetNbinsX(); l++) {
 879       if(_tempSwapN.GetBinContent(l+1) == 0) {
 880         _tempSwap.SetBinContent(l+1,0);
 881         _tempSwap.SetBinError(l+1,0);
 882       }
 883       else {
 884         Double_t vartimesnsq = totalptsquares*totalns - totalpts*totalpts;
 885         Double_t variance = vartimesnsq/(totalns*(totalns-1.));
 886         Double_t SDOMSq = variance / _tempSwapN.GetBinContent(l+1);
 887         Double_t SDOMSqOverMeanSq = SDOMSq * _tempSwapN.GetBinContent(l+1) * _tempSwapN.GetBinContent(l+1) / (_tempSwapN.GetBinContent(l+1) * _tempSwapN.GetBinContent(l+1));
 888         Double_t poissonfrac = 1./_tempSwapN.GetBinContent(l+1);
 889         Double_t vartotalfrac = SDOMSqOverMeanSq + poissonfrac;
 890         Double_t vartotal = vartotalfrac * _tempSwap.GetBinContent(l+1) * _tempSwap.GetBinContent(l+1);
 891         if(vartotal > 0.0001) _tempSwap.SetBinError(l+1,TMath::Sqrt(vartotal));
 892         else {
 893           _tempSwap.SetBinContent(l+1,0);
 894           _tempSwap.SetBinError(l+1,0);
 895         }
 896       }
 897     }
 898   }
 899
 900   fFitModulation->SetParameter(0, fLocalRho->GetVal());
 901   switch (fFitModulationType) {
 902     case kNoFit : {
 903       fFitModulation->FixParameter(0, fLocalRho->GetVal() );
 904       freeParams = 0;
 905     } break;
 906     case kV2 : {
 907       fFitModulation->FixParameter(4, psi2);
 908       freeParams = 1;
 909     } break;
 910     case kV3 : {
 911       fFitModulation->FixParameter(4, psi3);
 912       freeParams = 1;
 913     } break;
 914     case kCombined : {
 915       fFitModulation->FixParameter(4, psi2);
 916       fFitModulation->FixParameter(6, psi3);
 917       freeParams = 2;
 918     } break;
 919     case kFourierSeries : {
 920       // in this approach, an explicit calculation will be made of vn = sqrt(xn^2+yn^2)
 921       // where x[y] = Integrate[r(phi)cos[sin](n phi)dphi, 0, 2pi]
 922       Double_t cos2(0), sin2(0), cos3(0), sin3(0), sumPt(0);
 923       for(Int_t i(0); i < iTracks; i++) {
 924         AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(i));
 925         if(!PassesCuts(track) || track->Pt() > fSoftTrackMaxPt || track->Pt() < fSoftTrackMinPt) continue;
 926         sumPt += track->Pt();
 927         cos2 += track->Pt()*TMath::Cos(2*PhaseShift(track->Phi()-psi2));
 928         sin2 += track->Pt()*TMath::Sin(2*PhaseShift(track->Phi()-psi2));
 929         cos3 += track->Pt()*TMath::Cos(3*PhaseShift(track->Phi()-psi3));
 930         sin3 += track->Pt()*TMath::Sin(3*PhaseShift(track->Phi()-psi3));
 931       }
 932       fFitModulation->SetParameter(3, TMath::Sqrt(cos2*cos2+sin2*sin2)/fLocalRho->GetVal());
 933       fFitModulation->SetParameter(4, psi2);
 934       fFitModulation->SetParameter(6, psi3);
 935       fFitModulation->SetParameter(7, TMath::Sqrt(cos3*cos3+sin3*sin3)/fLocalRho->GetVal());
 936     } break;
 937     default : break;
 938   }
 939   _tempSwap.Fit(fFitModulation, fFitModulationOptions.Data(), "", lowBound, upBound);
 940   // the quality of the fit is evaluated from 1 - the cdf of the chi square distribution
 941   // three methods are available, all with their drawbacks. all are stored, one is selected to do the cut
 942   Int_t NDF(_tempSwap.GetXaxis()->GetNbins()-freeParams);
 943   if(NDF == 0) return kFALSE;
 944   Double_t CDF(1.-ChiSquareCDF(NDF, ChiSquare(_tempSwap, fFitModulation)));
 945   Double_t CDFROOT(1.-ChiSquareCDF(NDF, fFitModulation->GetChisquare()));
 946   // fill the values and centrality correlation (redundant but easy on the eyes)
 947   fHistPvalueCDF->Fill(CDF);
 948   fHistPvalueCDFCent->Fill(fCent, CDF);
 949   fHistPvalueCDFROOT->Fill(CDFROOT);
 950   fHistPvalueCDFROOTCent->Fill(fCent, CDFROOT);
 951   fHistChi2ROOTCent->Fill(fCent, fFitModulation->GetChisquare()/((float)NDF));
 952   fHistChi2Cent->Fill(fCent, ChiSquare(_tempSwap, fFitModulation)/((float)NDF));
 953   fHistPChi2Root->Fill(CDFROOT, fFitModulation->GetChisquare()/((float)NDF));
 954   fHistPChi2->Fill(CDF, ChiSquare(_tempSwap, fFitModulation)/((float)NDF));
 955
 956   // variable CDF is used for making cuts, so we fill it with the selected p-value
 957   switch (fFitGoodnessTest) {
 958     case kChi2ROOT : {
 959       CDF = CDFROOT;
 960     } break;
 961     case kChi2Poisson : break;      // CDF is already CDF
 962     default: break;
 963   }
 964
 965   if(fFitControl) {
 966     // as an additional quality check, see if fitting a control fit has a higher significance
 967     _tempSwap.Fit(fFitControl, fFitModulationOptions.Data(), "", lowBound, upBound);
 968     Double_t CDFControl(-1.);
 969     switch (fFitGoodnessTest) {
 970       case kChi2ROOT : {
 971         CDFControl = 1.-ChiSquareCDF(fFitControl->GetNDF(), fFitModulation->GetChisquare());
 972       } break;
 973       case kChi2Poisson : {
 974         CDFControl = 1.-ChiSquareCDF(fFitControl->GetNDF(), ChiSquare(_tempSwap, fFitModulation));
 975       } break;
 976       default: break;
 977     }
 978     if(CDFControl > CDF) {
 979       CDF = -1.; // control fit is more significant, so throw out the 'old' fit
 980       fHistRhoStatusCent->Fill(fCent, -1);
 981     }
 982   }
 983   if(CDF >= fMinPvalue && CDF <= fMaxPvalue && ( fAbsVnHarmonics && fFitModulation->GetMinimum(0, TMath::TwoPi()) > 0)) {       // fit quality. not that although with limited acceptance the fit is performed on just
 984     // part of phase space, the requirement that energy desntiy is larger than zero is applied
 985     // to the FULL spectrum
 986     fHistRhoStatusCent->Fill(fCent, 0.);
 987     // for LOCAL didactic purposes, save the  best and the worst fits
 988     // this routine can produce a lot of output histograms (it's not memory 'safe') and will not work on GRID
 989     // since the output will become unmergeable (i.e. different nodes may produce conflicting output)
 990     switch (fRunModeType) {
 991       case kLocal : {
 992         if(gRandom->Uniform(0, 100) > fPercentageOfFits) break;
 993         static Int_t didacticCounterBest(0);
 994         TProfile* didacticProfile = (TProfile*)_tempSwap.Clone(Form("Fit_%i_1-CDF_%.3f_cen_%i_%s", didacticCounterBest, CDF, fInCentralitySelection, detector.Data()));
 995         TF1* didacticFit = (TF1*)fFitModulation->Clone(Form("fit_%i_CDF_%.3f_cen_%i_%s", didacticCounterBest, CDF, fInCentralitySelection, detector.Data()));
 996         switch(fFitModulationType) {
 997           case kCombined : {
 998             // to make a nice picture also plot the separate components (v2 and v3) of the fit
 999             // only done for cobined fit where there are actually components to split ...
1000             TF1* v0(new TF1("dfit_kV2", "[0]", 0, TMath::TwoPi()));
1001             v0->SetParameter(0, didacticFit->GetParameter(0));        // normalization
1002             v0->SetLineColor(kMagenta);
1003             v0->SetLineStyle(7);
1004             didacticProfile->GetListOfFunctions()->Add(v0);
1005             TF1* v2(new TF1("dfit_kV2", "[0]*([1]+[2]*[3]*TMath::Cos([2]*(x-[4])))", 0, TMath::TwoPi()));
1006             v2->SetParameter(0, didacticFit->GetParameter(0));        // normalization
1007             v2->SetParameter(3, didacticFit->GetParameter(3));        // v2
1008             v2->FixParameter(1, 1.);        // constant
1009             v2->FixParameter(2, 2.);        // constant
1010             v2->FixParameter(4, didacticFit->GetParameter(4));        // psi2
1011             v2->SetLineColor(kGreen);
1012             didacticProfile->GetListOfFunctions()->Add(v2);
1013             TF1* v3(new TF1("dfit_kV3", "[0]*([1]+[2]*[3]*TMath::Cos([5]*(x-[4])))", 0, TMath::TwoPi()));
1014             v3->SetParameter(0, didacticFit->GetParameter(0));        // normalization
1015             v3->SetParameter(3, didacticFit->GetParameter(7));        // v3
1016             v3->FixParameter(1, 1.);        // constant
1017             v3->FixParameter(2, 2.);        // constant
1018             v3->FixParameter(4, didacticFit->GetParameter(6));        // psi3
1019             v3->FixParameter(5, 3.);        // constant
1020             v3->SetLineColor(kCyan);
1021             didacticProfile->GetListOfFunctions()->Add(v3);
1022           }
1023           default : break;
1024         }
1025         didacticProfile->GetListOfFunctions()->Add(didacticFit);
1026         didacticProfile->GetYaxis()->SetTitle("#frac{d #sum #it{p}_{T}}{d #varphi} [GeV/#it{c}]");
1027         didacticProfile->GetXaxis()->SetTitle("#varphi");
1028         fOutputListGood->Add(didacticProfile);
1029         didacticCounterBest++;
1030         TH2F* didacticSurface = BookTH2F(Form("surface_%s", didacticProfile->GetName()), "#phi", "#eta", 50, 0, TMath::TwoPi(), 50, -1, 1, -1, kFALSE);
1031         for(Int_t i(0); i < iTracks; i++) {
1032           AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(i));
1033           if(PassesCuts(track)) {
1034             if(fUsePtWeight) didacticSurface->Fill(track->Phi(), track->Eta(), track->Pt());
1035             else didacticSurface->Fill(track->Phi(), track->Eta());
1036           }
1037         }
1038         if(fExcludeLeadingJetsFromFit) {       // visualize the excluded region
1039           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);
1040           f2->SetParameters(excludeInPt/3.,excludeInPhi,.1,excludeInEta,.1);
1041           didacticSurface->GetListOfFunctions()->Add(f2);
1042         }
1043         fOutputListGood->Add(didacticSurface);
1044       } break;
1045       default : break;
1046     }
1047   } else {    // if the fit is of poor quality revert to the original rho estimate
1048       switch (fRunModeType) { // again see if we want to save the fit
1049         case kLocal : {
1050           static Int_t didacticCounterWorst(0);
1051           if(gRandom->Uniform(0, 100) > fPercentageOfFits) break;
1052           TProfile* didacticProfile = (TProfile*)_tempSwap.Clone(Form("Fit_%i_1-CDF_%.3f_cen_%i_%s", didacticCounterWorst, CDF, fInCentralitySelection, detector.Data() ));
1053           TF1* didacticFit = (TF1*)fFitModulation->Clone(Form("fit_%i_p_%.3f_cen_%i_%s", didacticCounterWorst, CDF, fInCentralitySelection, detector.Data()));
1054           didacticProfile->GetListOfFunctions()->Add(didacticFit);
1055           fOutputListBad->Add(didacticProfile);
1056           didacticCounterWorst++;
1057           } break;
1058         default : break;
1059       }
1060       switch (fFitModulationType) {
1061         case kNoFit : break;        // nothing to do
1062         case kCombined : fFitModulation->SetParameter(7, 0);        // no break
1063         case kFourierSeries : fFitModulation->SetParameter(7, 0);   // no break
1064         default : { // needs to be done if there was a poor fit
1065           fFitModulation->SetParameter(3, 0);
1066           fFitModulation->SetParameter(0, fLocalRho->GetVal());
1067         } break;
1068       }
1069     if(CDF > -.5) fHistRhoStatusCent->Fill(fCent, 1.);
1070     return kFALSE;  // return false if the fit is rejected
1071   }
1072   return kTRUE;
1073 }
1074 //_____________________________________________________________________________
1075 void AliAnalysisTaskLocalRhoDev::FillAnalysisSummaryHistogram() const
1076 {
1077   // Fill the analysis summary histrogram, saves all relevant analysis settigns
1078
1079   if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1080   fHistAnalysisSummary->GetXaxis()->SetBinLabel(2, "fJetRadius");
1081   fHistAnalysisSummary->SetBinContent(2, GetJetContainer()->GetJetRadius());
1082   fHistAnalysisSummary->GetXaxis()->SetBinLabel(3, "fJetEtaMin");
1083   fHistAnalysisSummary->SetBinContent(3, GetJetContainer()->GetJetEtaMin());
1084   fHistAnalysisSummary->GetXaxis()->SetBinLabel(4, "fJetEtaMax");
1085   fHistAnalysisSummary->SetBinContent(4, GetJetContainer()->GetJetEtaMax());
1086   fHistAnalysisSummary->GetXaxis()->SetBinLabel(5, "fJetPhiMin");
1087   fHistAnalysisSummary->SetBinContent(5, GetJetContainer()->GetJetPhiMin());
1088   fHistAnalysisSummary->GetXaxis()->SetBinLabel(6, "fJetPhiMax");
1089   fHistAnalysisSummary->SetBinContent(6, GetJetContainer()->GetJetPhiMin());
1090   fHistAnalysisSummary->GetXaxis()->SetBinLabel(34, "fitModulationType");
1091   fHistAnalysisSummary->SetBinContent(34, (int)fFitModulationType);
1092   fHistAnalysisSummary->GetXaxis()->SetBinLabel(35, "runModeType");
1093   fHistAnalysisSummary->SetBinContent(35, (int)fRunModeType);
1094   fHistAnalysisSummary->GetXaxis()->SetBinLabel(37, "iterator");
1095   fHistAnalysisSummary->SetBinContent(37, 1.);
1096   fHistAnalysisSummary->GetXaxis()->SetBinLabel(38, "fMinPvalue");
1097   fHistAnalysisSummary->SetBinContent(38, fMinPvalue);
1098   fHistAnalysisSummary->GetXaxis()->SetBinLabel(39, "fMaxPvalue");
1099   fHistAnalysisSummary->SetBinContent(39, fMaxPvalue);
1100   fHistAnalysisSummary->GetXaxis()->SetBinLabel(40, "fExcludeLeadingJetsFromFit");
1101   fHistAnalysisSummary->SetBinContent(40, fExcludeLeadingJetsFromFit);
1102   fHistAnalysisSummary->GetXaxis()->SetBinLabel(41, "fRebinSwapHistoOnTheFly");
1103   fHistAnalysisSummary->SetBinContent(41, (int)fRebinSwapHistoOnTheFly);
1104   fHistAnalysisSummary->GetXaxis()->SetBinLabel(42, "fUsePtWeight");
1105   fHistAnalysisSummary->SetBinContent(42, (int)fUsePtWeight);
1106   fHistAnalysisSummary->GetXaxis()->SetBinLabel(45, "fLocalJetMinEta");
1107   fHistAnalysisSummary->SetBinContent(45,fLocalJetMinEta );
1108   fHistAnalysisSummary->GetXaxis()->SetBinLabel(46, "fLocalJetMaxEta");
1109   fHistAnalysisSummary->SetBinContent(46, fLocalJetMaxEta);
1110   fHistAnalysisSummary->GetXaxis()->SetBinLabel(47, "fLocalJetMinPhi");
1111   fHistAnalysisSummary->SetBinContent(47, fLocalJetMinPhi);
1112   fHistAnalysisSummary->GetXaxis()->SetBinLabel(48, "fLocalJetMaxPhi");
1113   fHistAnalysisSummary->SetBinContent(48, fLocalJetMaxPhi);
1114   fHistAnalysisSummary->GetXaxis()->SetBinLabel(49, "fSoftTrackMinPt");
1115   fHistAnalysisSummary->SetBinContent(49, fSoftTrackMinPt);
1116   fHistAnalysisSummary->GetXaxis()->SetBinLabel(50, "fSoftTrackMaxPt");
1117   fHistAnalysisSummary->SetBinContent(50, fSoftTrackMaxPt);
1118   fHistAnalysisSummary->GetXaxis()->SetBinLabel(51, "fUseScaledRho");
1119   fHistAnalysisSummary->SetBinContent(51, fUseScaledRho);
1120   fHistAnalysisSummary->GetXaxis()->SetBinLabel(53, "used rho");
1121   fHistAnalysisSummary->GetXaxis()->SetBinLabel(54, "used small rho");
1122 }
1123
1124 //_____________________________________________________________________________
1125 void AliAnalysisTaskLocalRhoDev::FillEventPlaneHistograms(Double_t psi2, Double_t psi3) const
1126 {
1127   // Fill event plane histograms
1128
1129   if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1130   fHistPsi2[fInCentralitySelection]->Fill(psi2);
1131   fHistPsi3[fInCentralitySelection]->Fill(psi3);
1132 }
1133 //_____________________________________________________________________________
1134 Bool_t AliAnalysisTaskLocalRhoDev::PassesCuts(AliVEvent* event)
1135 {
1136     // event cuts
1137     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1138     if(!event) return kFALSE;
1139     // determine the run number to see if the track and jet cuts should be refreshed for semi-good TPC runs
1140     // only done if the runnumber changes, could be moved to a call to AliAnalysisTaskSE::Notify()
1141     if(fRunNumber != InputEvent()->GetRunNumber()) {
1142         fRunNumber = InputEvent()->GetRunNumber();        // set the current run number
1143         if(fDebug > 0) printf("__FUNC__ %s > NEW RUNNUMBER DETECTED \n ", __func__);
1144         // reset the cuts. should be a pointless operation except for the case where the run number changes
1145         // from semi-good back to good on one node, which is not a likely scenario
1146         AliAnalysisTaskEmcal::SetTrackPhiLimits(-10., 10.);
1147         AliAnalysisTaskEmcalJet::SetJetPhiLimits(-10., 10.);
1148         if(fCachedRho) {                // if there's a cached rho, it's the default, so switch back
1149             if(fDebug > 0) printf("__FUNC__ %s > replacing rho with cached rho \n ", __func__);
1150             fRho = fCachedRho;          // reset rho back to cached value. again, should be pointless
1151         }
1152         Bool_t flaggedAsSemiGood(kFALSE);       // not flagged as anything
1153         for(Int_t i(0); i < fExpectedSemiGoodRuns->GetSize(); i++) {
1154             if(fExpectedSemiGoodRuns->At(i) == fRunNumber) { // run is semi-good
1155                if(fDebug > 0) printf("__FUNC__ %s > semi-good tpc run detected, adjusting acceptance \n ", __func__);
1156                 flaggedAsSemiGood = kTRUE;
1157                 AliAnalysisTaskEmcalJet::SetJetPhiLimits(fSemiGoodJetMinPhi, fSemiGoodJetMaxPhi);       // just an acceptance cut, jets are obtained from full azimuth, so no edge effects
1158                 AliAnalysisTaskEmcal::SetTrackPhiLimits(fSemiGoodTrackMinPhi, fSemiGoodTrackMaxPhi);    // only affects vn extraction, NOT jet finding
1159                 // for semi-good runs, also try to get the 'small rho' estimate, if it is available
1160                 AliRhoParameter* tempRho(dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fNameSmallRho.Data())));
1161                 if(tempRho) {
1162                     if(fDebug > 0) printf("__FUNC__ %s > switching to small rho, caching normal rho \n ", __func__);
1163                     fHistAnalysisSummary->SetBinContent(54, 1.);        // bookkeep the fact that small rho is used
1164                     fCachedRho = fRho;          // cache the original rho ...
1165                     fRho = tempRho;             // ... and use the small rho
1166                 }
1167             }
1168         }
1169         if(!flaggedAsSemiGood) fHistAnalysisSummary->SetBinContent(53, 1.);                // bookkeep which rho estimate is used
1170     }
1171     return kTRUE;
1172 }
1173
1174 //_____________________________________________________________________________
1175 void AliAnalysisTaskLocalRhoDev::Terminate(Option_t *)
1176 {
1177   // Terminate
1178 }
1179
1180 //_____________________________________________________________________________
1181 void AliAnalysisTaskLocalRhoDev::SetModulationFit(TF1* fit)
1182 {
1183   // Set function to fit modulation
1184
1185   if (fFitModulation) delete fFitModulation;
1186   fFitModulation = fit;
1187 }