PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskRhoVnModulation.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /*
  17  * analysis task for jet flow preparation
  18  *
  19  * this task is part of the emcal jet framework and should be run in the emcaljet train
  20  * the following extensions to an accepted AliVEvent are expected:
  21  *      - (anti-kt) jets
  22  *      - background estimate rho
  23  *      - pico tracks
  24  *      aod's and esd's are handled transparently
  25  * the task will attempt to estimate a phi-dependent background density rho
  26  * by fitting vn harmonics
  27  *
  28  * author: Redmer Alexander Bertens, Utrecht Univeristy, Utrecht, Netherlands
  29  * rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl
  30  */
  31
  32 #include <TStyle.h>
  33 #include <TRandom3.h>
  34 #include <TChain.h>
  35 #include <TMath.h>
  36 #include <TF1.h>
  37 #include <TH1F.h>
  38 #include <TH2F.h>
  39 #include <TProfile.h>
  40
  41 #include <AliAnalysisTask.h>
  42 #include <AliAnalysisManager.h>
  43 #include <AliCentrality.h>
  44 #include <AliVVertex.h>
  45 #include <AliESDEvent.h>
  46 #include <AliAODEvent.h>
  47
  48 #include <AliPicoTrack.h>
  49 #include <AliEmcalJet.h>
  50 #include <AliRhoParameter.h>
  51
  52 #include "AliAnalysisTaskRhoVnModulation.h"
  53
  54
  55 class AliAnalysisTaskRhoVnModulation;
  56 using namespace std;
  57
  58 ClassImp(AliAnalysisTaskRhoVnModulation)
  59
  60 AliAnalysisTaskRhoVnModulation::AliAnalysisTaskRhoVnModulation() : AliAnalysisTaskEmcalJet("AliAnalysisTaskRhoVnModulation", kTRUE),
  61     fDebug(0), fInitialized(0), fFillQAHistograms(kTRUE), fFitModulationType(kNoFit), fDetectorType(kTPC), fFitModulationOptions("Q"), fRunModeType(kGrid), fDataType(kESD), fRandom(0), fMappedRunNumber(0), fInCentralitySelection(-1), fFitModulation(0), fNameJetClones(0), fNamePicoTrackClones(0), fNameRho(0), fAbsVertexZ(10), fHistCentrality(0), fHistVertexz(0), fHistRunnumbersPhi(0), fHistRunnumbersEta(0), fMinDisanceRCtoLJ(0), fRandomConeRadius(0.4), fOutputList(0), fOutputListGood(0), fOutputListBad(0), fHistAnalysisSummary(0), fHistSwap(0), fProfVn(0), fHistPsi2(0), fHistPsi2Spread(0), fHistPsiVZEROA(0), fHistPsiVZEROC(0), fHistPsiTPC(0),
  62    fHistRhoVsMult(0), fHistRhoVsCent(0), fHistRhoAVsMult(0), fHistRhoAVsCent(0) {
  63     for(Int_t i(0); i < 10; i++) {
  64         fHistPicoTrackPt[i] = 0;
  65         fHistPicoTrackPhi[i] = 0;
  66         fHistPicoTrackEta[i] = 0;
  67         fHistPicoCat1[i] = 0;
  68         fHistPicoCat2[i] = 0;
  69         fHistPicoCat3[i] = 0;
  70         /* fHistClusterPt[i] = 0; */
  71         /* fHistClusterPhi[i] = 0; */
  72         /* fHistClusterEta[i] = 0; */
  73         /* fHistClusterCorrPt[i] = 0; */
  74         /* fHistClusterCorrPhi[i] = 0; */
  75         /* fHistClusterCorrEta[i] = 0; */
  76         fHistRhoPackage[i] = 0;
  77         fHistRho[i] = 0;
  78         fHistRCPhiEta[i] = 0;
  79         fHistRhoVsRCPt[i] = 0;
  80         fHistRCPt[i] = 0;
  81         fHistDeltaPtRC[i] = 0;
  82         fHistRCPhiEtaExLJ[i] = 0;
  83         fHistRhoVsRCPtExLJ[i] = 0;
  84         fHistRCPtExLJ[i] = 0;
  85         fHistDeltaPtRCExLJ[i] = 0;
  86         fHistRCPhiEtaRand[i] = 0;
  87         fHistRhoVsRCPtRand[i] = 0;
  88         fHistRCPtRand[i] = 0;
  89         fHistDeltaPtRCRand[i] = 0;
  90         fHistJetPtRaw[i] = 0;
  91         fHistJetPt[i] = 0;
  92         fHistJetEtaPhi[i] = 0;
  93         fHistJetPtArea[i] = 0;
  94         fHistJetPtConstituents[i] = 0;
  95         fHistJetPtInPlaneVZEROA[i] = 0;
  96         fHistJetPtOutPlaneVZEROA[i] = 0;
  97         fHistJetPtMidPlaneVZEROA[i] = 0;
  98         fHistJetPtInPlaneVZEROC[i] = 0;
  99         fHistJetPtOutPlaneVZEROC[i] = 0;
 100         fHistJetPtMidPlaneVZEROC[i] = 0;
 101         fHistJetPtInPlaneTPC[i] = 0;
 102         fHistJetPtOutPlaneTPC[i] = 0;
 103         fHistJetPtMidPlaneTPC[i] = 0;
 104         fHistJetPsiTPCPt[i] = 0;
 105         fHistJetPsiVZEROAPt[i] = 0;
 106         fHistJetPsiVZEROCPt[i] = 0;
 107         fHistDeltaPhiVZEROA[i] = 0;
 108         fHistDeltaPhiVZEROC[i] = 0;
 109         fHistDeltaPhiTPC[i] = 0;
 110     }
 111     // default constructor
 112 }
 113 //_____________________________________________________________________________
 114 AliAnalysisTaskRhoVnModulation::AliAnalysisTaskRhoVnModulation(const char* name, runModeType type) : AliAnalysisTaskEmcalJet(name, kTRUE),
 115   fDebug(0), fInitialized(0), fFillQAHistograms(kTRUE), fFitModulationType(kNoFit), fDetectorType(kTPC), fFitModulationOptions("Q"), fRunModeType(type), fDataType(kESD), fRandom(0), fMappedRunNumber(0), fInCentralitySelection(-1), fFitModulation(0), fNameJetClones(0), fNamePicoTrackClones(0), fNameRho(0), fAbsVertexZ(10), fHistCentrality(0), fHistVertexz(0), fHistRunnumbersPhi(0), fHistRunnumbersEta(0), fMinDisanceRCtoLJ(0), fRandomConeRadius(0.4), fOutputList(0), fOutputListGood(0), fOutputListBad(0), fHistAnalysisSummary(0), fHistSwap(0), fProfVn(0), fHistPsi2(0), fHistPsi2Spread(0), fHistPsiVZEROA(0), fHistPsiVZEROC(0), fHistPsiTPC(0),
 116    fHistRhoVsMult(0), fHistRhoVsCent(0), fHistRhoAVsMult(0), fHistRhoAVsCent(0) {
 117     for(Int_t i(0); i < 10; i++) {
 118         fHistPicoTrackPt[i] = 0;
 119         fHistPicoTrackPhi[i] = 0;
 120         fHistPicoTrackEta[i] = 0;
 121         fHistPicoCat1[i] = 0;
 122         fHistPicoCat2[i] = 0;
 123         fHistPicoCat3[i] = 0;
 124         /* fHistClusterPt[i] = 0; */
 125         /* fHistClusterPhi[i] = 0; */
 126         /* fHistClusterEta[i] = 0; */
 127         /* fHistClusterCorrPt[i] = 0; */
 128         /* fHistClusterCorrPhi[i] = 0; */
 129         /* fHistClusterCorrEta[i] = 0; */
 130         fHistRhoPackage[i] = 0;
 131         fHistRho[i] = 0;
 132         fHistRCPhiEta[i] = 0;
 133         fHistRhoVsRCPt[i] = 0;
 134         fHistRCPt[i] = 0;
 135         fHistDeltaPtRC[i] = 0;
 136         fHistRCPhiEtaExLJ[i] = 0;
 137         fHistRhoVsRCPtExLJ[i] = 0;
 138         fHistRCPtExLJ[i] = 0;
 139         fHistDeltaPtRCExLJ[i] = 0;
 140         fHistRCPhiEtaRand[i] = 0;
 141         fHistRhoVsRCPtRand[i] = 0;
 142         fHistRCPtRand[i] = 0;
 143         fHistDeltaPtRCRand[i] = 0;
 144         fHistJetPtRaw[i] = 0;
 145         fHistJetPt[i] = 0;
 146         fHistJetEtaPhi[i] = 0;
 147         fHistJetPtArea[i] = 0;
 148         fHistJetPtConstituents[i] = 0;
 149         fHistJetPtInPlaneVZEROA[i] = 0;
 150         fHistJetPtOutPlaneVZEROA[i] = 0;
 151         fHistJetPtMidPlaneVZEROA[i] = 0;
 152         fHistJetPtInPlaneVZEROC[i] = 0;
 153         fHistJetPtOutPlaneVZEROC[i] = 0;
 154         fHistJetPtMidPlaneVZEROC[i] = 0;
 155         fHistJetPtInPlaneTPC[i] = 0;
 156         fHistJetPtOutPlaneTPC[i] = 0;
 157         fHistJetPtMidPlaneTPC[i] = 0;
 158         fHistJetPsiTPCPt[i] = 0;
 159         fHistJetPsiVZEROAPt[i] = 0;
 160         fHistJetPsiVZEROCPt[i] = 0;
 161         fHistDeltaPhiVZEROA[i] = 0;
 162         fHistDeltaPhiVZEROC[i] = 0;
 163         fHistDeltaPhiTPC[i] = 0;
 164     }
 165     // constructor
 166     DefineInput(0, TChain::Class());
 167     DefineOutput(1, TList::Class());
 168     switch (fRunModeType) {
 169         case kLocal : {
 170             gStyle->SetOptFit(1);
 171             DefineOutput(2, TList::Class());
 172             DefineOutput(3, TList::Class());
 173         } break;
 174         default: fDebug = -1;   // suppress debug info explicitely when not running locally
 175     }
 176 }
 177 //_____________________________________________________________________________
 178 AliAnalysisTaskRhoVnModulation::~AliAnalysisTaskRhoVnModulation()
 179 {
 180     // destructor
 181     if(fOutputList)     delete fOutputList;
 182     if(fOutputListGood) delete fOutputListGood;
 183     if(fOutputListBad)  delete fOutputListBad;
 184     if(fFitModulation)  delete fFitModulation;
 185     if(fHistSwap)       delete fHistSwap;
 186 }
 187 //_____________________________________________________________________________
 188 Bool_t AliAnalysisTaskRhoVnModulation::InitializeAnalysis()
 189 {
 190     // initialize the anaysis
 191     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 192     if(fMinDisanceRCtoLJ==0) fMinDisanceRCtoLJ = .5*fJetRadius;
 193     if(dynamic_cast<AliAODEvent*>(InputEvent())) fDataType = kAOD; // determine the datatype
 194     else if(dynamic_cast<AliESDEvent*>(InputEvent())) fDataType = kESD;
 195     if(!fRandom) fRandom = new TRandom3(0);  // get a randomized if one hasn't been user-supplied
 196     switch (fFitModulationType)  {
 197         case kNoFit : { SetModulationFit(new TF1("fix_kNoFit", "[0]", 0, TMath::TwoPi())); } break;
 198         case kV2 : {
 199             SetModulationFit(new TF1("fit_kV2", "[0]*([1]+[2]*[3]*TMath::Cos([2]*(x-[4])))", 0, TMath::TwoPi()));
 200             fFitModulation->SetParameter(0, 0.);        // normalization
 201             fFitModulation->SetParameter(3, 0.2);       // v2
 202             fFitModulation->FixParameter(1, 1.);        // constant
 203             fFitModulation->FixParameter(2, 2.);        // constant
 204         } break;
 205         case kV3: {
 206             SetModulationFit(new TF1("fit_kV3", "[0]*([1]+[2]*[3]*TMath::Cos([2]*(x-[4])))", 0, TMath::TwoPi()));
 207             fFitModulation->SetParameter(0, 0.);        // normalization
 208             fFitModulation->SetParameter(3, 0.2);       // v3
 209             fFitModulation->FixParameter(1, 1.);        // constant
 210             fFitModulation->FixParameter(2, 3.);        // constant
 211         } break;
 212         case kCombined : {
 213              SetModulationFit(new TF1("fit_kCombined", "[0]*([1]+[2]*([3]*TMath::Cos([2]*(x-[4]))+[7]*TMath::Cos([5]*(x-[6]))))"));
 214              fFitModulation->SetParameter(0, 0.);       // normalization
 215              fFitModulation->SetParameter(3, 0.2);      // v2
 216              fFitModulation->FixParameter(1, 1.);       // constant
 217              fFitModulation->FixParameter(2, 2.);       // constant
 218              fFitModulation->FixParameter(5, 3.);       // constant
 219              fFitModulation->SetParameter(7, 0.2);      // v3
 220         } break;
 221         default: break;
 222     }
 223     switch (fRunModeType) {
 224         case kGrid : { fFitModulationOptions += "N0"; } break;
 225         default : break;
 226     }
 227     return kTRUE;
 228 }
 229 //_____________________________________________________________________________
 230 TH1F* AliAnalysisTaskRhoVnModulation::BookTH1F(const char* name, const char* x, Int_t bins, Double_t min, Double_t max, Int_t c)
 231 {
 232     // book a TH1F and connect it to the output container
 233     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 234     if(!fOutputList) return 0x0;
 235     if(c!=-1) name = Form("%s_%i", name, c);
 236     TH1F* histogram = new TH1F(name, name, bins, min, max);
 237     histogram->GetXaxis()->SetTitle(x);
 238     histogram->GetYaxis()->SetTitle("[counts]");
 239     histogram->Sumw2();
 240     fOutputList->Add(histogram);
 241     return histogram;
 242 }
 243 //_____________________________________________________________________________
 244 TH2F* AliAnalysisTaskRhoVnModulation::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)
 245 {
 246     // book a TH2F and connect it to the output container
 247     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 248     if(!fOutputList) return 0x0;
 249     if(c!=-1) name = Form("%s_%i", name, c);
 250     TH2F* histogram = new TH2F(name, name, binsx, minx, maxx, binsy, miny, maxy);
 251     histogram->GetXaxis()->SetTitle(x);
 252     histogram->GetYaxis()->SetTitle(y);
 253     histogram->Sumw2();
 254     fOutputList->Add(histogram);
 255     return histogram;
 256 }
 257 //_____________________________________________________________________________
 258 void AliAnalysisTaskRhoVnModulation::UserCreateOutputObjects()
 259 {
 260     // create output objects
 261     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 262     fOutputList = new TList();
 263     fOutputList->SetOwner(kTRUE);
 264     // global QA
 265     fHistCentrality =           BookTH1F("fHistCentrality", "centrality \%", 102, -2, 100);
 266     fHistVertexz =              BookTH1F("fHistVertexz", "vertex z (cm)", 100, -12, 12);
 267
 268     // pico track kinematics
 269     for(Int_t i(0); i < 10; i++) {
 270         fHistPicoTrackPt[i] =          BookTH1F("fHistPicoTrackPt", "p_{t} [GeV/c]", 100, 0, 50, i);
 271         fHistPicoTrackPhi[i] =         BookTH1F("fHistPicoTrackPhi", "#phi", 100, 0, TMath::TwoPi(), i);
 272         fHistPicoTrackEta[i] =         BookTH1F("fHistPicoTrackEta", "#eta", 100, -1, 1, i);
 273         if(fFillQAHistograms) {
 274             fHistPicoCat1[i] =             BookTH2F("fHistPicoCat1", "#eta", "#phi", 100, -1, 1, 100, 0, TMath::TwoPi(), i);
 275             fHistPicoCat2[i] =             BookTH2F("fHistPicoCat2", "#eta", "#phi", 100, -1, 1, 100, 0, TMath::TwoPi(), i);
 276             fHistPicoCat3[i] =             BookTH2F("fHistPicoCat3", "#eta", "#phi", 100, -1, 1, 100, 0, TMath::TwoPi(), i);
 277         }
 278         // emcal kinematics
 279         /* fHistClusterPt[i] =            BookTH1F("fHistClusterPt", "p_{t} [GeV/c]", 100, 0, 100, i); */
 280         /* fHistClusterPhi[i] =           BookTH1F("fHistClusterPhi", "#phi", 100, 0, TMath::TwoPi(), i); */
 281         /* fHistClusterEta[i] =           BookTH1F("fHistClusterEta", "#eta", 100, -5, 5); */
 282
 283         // emcal kinematics after hadronic correction
 284         /* fHistClusterCorrPt[i] =        BookTH1F("fHistClusterCorrPt", "p_{t} [GeV/c]", 100, 0, 100, i); */
 285         /* fHistClusterCorrPhi[i] =       BookTH1F("fHistClusterCorrPhi", "#phi", 100, 0, TMath::TwoPi(), i); */
 286         /* fHistClusterCorrEta[i] =       BookTH1F("fHistClusterCorrEta", "#eta", 100, -5, 5, i); */
 287     }
 288
 289     // event plane estimates and quality
 290     fHistPsi2 =                 new TProfile("fHistPsi2", "fHistPsi2", 3, 0, 3);
 291     fHistPsi2->Sumw2();
 292     fHistPsi2Spread =           new TProfile("fHistPsi2Spread", "fHistPsi2Spread", 3, 0, 3);
 293     fHistPsi2Spread->Sumw2();
 294     fHistPsi2->GetXaxis()->SetBinLabel(1, "<#Psi_{2, VZEROA}>");
 295     fHistPsi2->GetXaxis()->SetBinLabel(2, "<#Psi_{2, VZEROC}>");
 296     fHistPsi2->GetXaxis()->SetBinLabel(3, "<#Psi_{2, TPC}>");
 297     fHistPsi2Spread->GetXaxis()->SetBinLabel(1, "<#Psi_{2, VZEROA} - #Psi_{2, VZEROC}>");
 298     fHistPsi2Spread->GetXaxis()->SetBinLabel(2, "<#Psi_{2, VZEROC} - #Psi_{2, TPC}>");
 299     fHistPsi2Spread->GetXaxis()->SetBinLabel(3, "<#Psi_{2, VZEROC} - #Psi_{2, TPC}>");
 300     fOutputList->Add(fHistPsi2);
 301     fOutputList->Add(fHistPsi2Spread);
 302     fHistPsiVZEROA =            BookTH1F("fHistPsiVZEROA", "#Psi_{VZEROA}", 100, -.5*TMath::Pi(), .5*TMath::Pi());
 303     fHistPsiVZEROC =            BookTH1F("fHistPsiVZEROC", "#Psi_{VZEROC}", 100, -.5*TMath::Pi(), .5*TMath::Pi());
 304     fHistPsiTPC =               BookTH1F("fHistPsiTPC", "#Psi_{TPC}", 100, -.5*TMath::Pi(), .5*TMath::Pi());
 305
 306     // background
 307     for(Int_t i(0); i < 10; i ++) {
 308         fHistRhoPackage[i] =           BookTH1F("fHistRhoPackage",  "#rho [GeV/c]", 100, 0, 150, i);
 309         fHistRho[i] =                  BookTH1F("fHistRho", "#rho [GeV/c]", 100, 0, 150, i);
 310     }
 311     fHistRhoVsMult =            BookTH2F("fHistRhoVsMult", "multiplicity", "#rho [GeV/c]", 100, 0, 4000, 100, 0, 250);
 312     fHistRhoVsCent =            BookTH2F("fHistRhoVsCent", "centrality", "#rho [GeV/c]", 100, 0, 100, 100, 0, 250);
 313     fHistRhoAVsMult =           BookTH2F("fHistRhoAVsMult", "multiplicity", "#rho * A (jet) [GeV/c]", 100, 0, 4000, 100, 0, 50);
 314     fHistRhoAVsCent =           BookTH2F("fHistRhoAVsCent", "centrality", "#rho * A (jet) [GeV/c]", 100, 0, 100, 100, 0, 50);
 315
 316     // delta pt distributions
 317     for(Int_t i(0); i < 10; i ++) {
 318         fHistRCPhiEta[i] =             BookTH2F("fHistRCPhiEta", "#phi (RC)", "#eta (RC)", 100, 0, TMath::TwoPi(), 100, -1, 1, i);
 319         fHistRhoVsRCPt[i] =            BookTH2F("fHistRhoVsRCPt", "p_{t} (RC) [GeV/c]", "#rho * A (RC) [GeV/c]", 100, 0, 300, 100, 0, 250, i);
 320         fHistRCPt[i] =                 BookTH1F("fHistRCPt", "p_{t} (RC) [GeV/c]", 130, -20, 150, i);
 321         fHistDeltaPtRC[i] =            BookTH1F("fHistDeltaPtRC", "#delta p_{t} [GeV/c]", 180, -50, 150, i);
 322         fHistRCPhiEtaExLJ[i] =         BookTH2F("fHistRCPhiEtaExLJ", "#phi (RC)", "#eta (RC)", 100, 0, TMath::TwoPi(), 100, -1, 1, i);
 323         fHistRhoVsRCPtExLJ[i] =        BookTH2F("fHistRhoVsRCPtExLJ", "p_{t} (RC) [GeV/c]", "#rho * A (RC) [GeV/c]", 100, 0, 300, 100, 0, 250, i);
 324         fHistRCPtExLJ[i] =             BookTH1F("fHistRCPtExLJ", "p_{t} (RC) [GeV/c]", 130, -20, 150, i);
 325         fHistDeltaPtRCExLJ[i] =        BookTH1F("fHistDeltaPtRCExLJ", "#delta p_{t} [GeV/c]", 180, -50, 150, i);
 326         fHistRCPhiEtaRand[i] =         BookTH2F("fHistRCPhiEtaRand", "#phi (RC)", "#eta (RC)", 100, 0, TMath::TwoPi(), 100, -1, 1, i);
 327         fHistRhoVsRCPtRand[i] =        BookTH2F("fHistRhoVsRCPtRand", "p_{t} (RC) [GeV/c]", "#rho * A (RC) [GeV/c]", 100, 0, 300, 100, 0, 250, i);
 328         fHistRCPtRand[i] =             BookTH1F("fHistRCPtRand", "p_{t} (RC) [GeV/c]", 130, -20, 150, i);
 329         fHistDeltaPtRCRand[i] =        BookTH1F("fHistDeltaPtRCRand", "#delta p_{t} [GeV/c]", 180, -50, 150, i);
 330
 331         // jet histograms (after kinematic cuts)
 332         fHistJetPtRaw[i] =             BookTH1F("fHistJetPtRaw", "p_{t} RAW [GeV/c]", 200, -50, 150, i);
 333         fHistJetPt[i] =                BookTH1F("fHistJetPt", "p_{t} [GeV/c]", 200, -50, 150, i);
 334         fHistJetEtaPhi[i] =            BookTH2F("fHistJetEtaPhi", "#eta", "#phi", 100, -1, 1, 100, 0, TMath::TwoPi(), i);
 335         fHistJetPtArea[i] =            BookTH2F("fHistJetPtArea", "p_{t} [GeV/c]", "Area", 200, -50, 150, 60, 0, 0.3, i);
 336         fHistJetPtConstituents[i] =    BookTH2F("fHistJetPtConstituents", "p_{t} [GeV/c]", "Area", 200, -50, 150, 60, 0, 150, i);
 337
 338         // in plane and out of plane spectra
 339         fHistJetPtInPlaneVZEROA[i] =   BookTH1F("fHistJetPtInPlaneVZEROA", "p_{t} [GeV/c]", 200, -50, 150, i);
 340         fHistJetPtOutPlaneVZEROA[i] =  BookTH1F("fHistJetPtOutPlaneVZEROA", "p_{t} [GeV/c]", 200, -50, 150, i);
 341         fHistJetPtMidPlaneVZEROA[i] =  BookTH1F("fHistJetPtMidPlaneVZEROA", "p_{t} [GeV/c]", 200, -50, 150, i);
 342         fHistJetPtInPlaneVZEROC[i] =   BookTH1F("fHistJetPtInPlaneVZEROC", "p_{t} [GeV/c]", 200, -50, 150, i);
 343         fHistJetPtOutPlaneVZEROC[i] =  BookTH1F("fHistJetPtOutPlaneVZEROC", "p_{t} [GeV/c]", 200, -50, 150, i);
 344         fHistJetPtMidPlaneVZEROC[i] =  BookTH1F("fHistJetPtMidPlaneVZEROC", "p_{t} [GeV/c]", 200, -50, 150, i);
 345         fHistJetPtInPlaneTPC[i] =      BookTH1F("fHistJetPtInPlaneTPC", "p_{t} [GeV/c]", 200, -50, 150, i);
 346         fHistJetPtOutPlaneTPC[i] =     BookTH1F("fHistJetPtOutPlaneTPC", "p_{t} [GeV/c]", 200, -50, 150, i);
 347         fHistJetPtMidPlaneTPC[i] =     BookTH1F("fHistJetPtMidPlaneTPC", "p_{t} [GeV/c]", 200, -50, 150, i);
 348         fHistJetPsiTPCPt[i] =          BookTH2F("fHistJetPsiTPCPt", "#phi_{jet} - #Psi_{TPC}", "p_{t} [GeV/c]", 100, 0., TMath::TwoPi(), 100, -50, 100, i);
 349         fHistJetPsiVZEROAPt[i] =       BookTH2F("fHistJetPsiVZEROAPt", "#phi_{jet} - #Psi_{VZEROA}", "p_{t} [GeV/c]", 100, 0., TMath::TwoPi(), 100, -50, 100, i);
 350         fHistJetPsiVZEROCPt[i] =       BookTH2F("fHistJetPsiVZEROCPt", "#phi_{jet} - #Psi_{VZEROC}", "p_{t} [GeV/c]", 100, 0., TMath::TwoPi(), 100, -50, 100, i);
 351
 352         // phi minus psi
 353         fHistDeltaPhiVZEROA[i] =       BookTH1F("fHistDeltaPhiVZEROA", "#phi_{jet} - #Psi_{VZEROA}", 100, 0, TMath::TwoPi(), i);
 354         fHistDeltaPhiVZEROC[i] =       BookTH1F("fHistDeltaPhiVZEROC", "#phi_{jet} - #Psi_{VZEROC}", 100, 0, TMath::TwoPi(), i);
 355         fHistDeltaPhiTPC[i] =          BookTH1F("fHistDeltaPhiTPC", "#phi_{jet} - #Psi_{TPC}", 100, 0, TMath::TwoPi(), i);
 356     }
 357
 358     // analysis summary histrogram, saves all relevant analysis settigns
 359     fHistAnalysisSummary = BookTH1F("fHistAnalysisSummary", "flag", 37, -0.5, 37.5);
 360     fHistAnalysisSummary->GetXaxis()->SetBinLabel(1, "fjetRadius");
 361     fHistAnalysisSummary->SetBinContent(1, fJetRadius);
 362     fHistAnalysisSummary->GetXaxis()->SetBinLabel(2, "fPtBiasJetTrack");
 363     fHistAnalysisSummary->SetBinContent(2, fPtBiasJetTrack);
 364     fHistAnalysisSummary->GetXaxis()->SetBinLabel(3, "fPtBiasJetClus");
 365     fHistAnalysisSummary->SetBinContent(3, fPtBiasJetClus);
 366     fHistAnalysisSummary->GetXaxis()->SetBinLabel(4, "fJetPtCut");
 367     fHistAnalysisSummary->SetBinContent(4, fJetPtCut);
 368     fHistAnalysisSummary->GetXaxis()->SetBinLabel(5, "fJetAreaCut");
 369     fHistAnalysisSummary->SetBinContent(5, fJetAreaCut);
 370     fHistAnalysisSummary->GetXaxis()->SetBinLabel(6, "fPercAreaCut");
 371     fHistAnalysisSummary->SetBinContent(6, fPercAreaCut);
 372     fHistAnalysisSummary->GetXaxis()->SetBinLabel(7, "fAreaEmcCut");
 373     fHistAnalysisSummary->SetBinContent(7, fAreaEmcCut);
 374     fHistAnalysisSummary->GetXaxis()->SetBinLabel(8, "fJetMinEta");
 375     fHistAnalysisSummary->SetBinContent(8, fJetMinEta);
 376     fHistAnalysisSummary->GetXaxis()->SetBinLabel(9, "fJetMaxEta");
 377     fHistAnalysisSummary->SetBinContent(9, fJetMaxEta);
 378     fHistAnalysisSummary->GetXaxis()->SetBinLabel(10, "fJetMinPhi");
 379     fHistAnalysisSummary->SetBinContent(10, fJetMinPhi);
 380     fHistAnalysisSummary->GetXaxis()->SetBinLabel(11, "fJetMaxPhi");
 381     fHistAnalysisSummary->SetBinContent(11, fJetMaxPhi);
 382     fHistAnalysisSummary->GetXaxis()->SetBinLabel(12, "fMaxClusterPt");
 383     fHistAnalysisSummary->SetBinContent(12, fMaxClusterPt);
 384     fHistAnalysisSummary->GetXaxis()->SetBinLabel(13, "fMaxTrackPt");
 385     fHistAnalysisSummary->SetBinContent(13, fMaxTrackPt);
 386     fHistAnalysisSummary->GetXaxis()->SetBinLabel(14, "fLeadingHadronType");
 387     fHistAnalysisSummary->SetBinContent(14, fLeadingHadronType);
 388     fHistAnalysisSummary->GetXaxis()->SetBinLabel(15, "fAnaType");
 389     fHistAnalysisSummary->SetBinContent(15, fAnaType);
 390     fHistAnalysisSummary->GetXaxis()->SetBinLabel(16, "fForceBeamType");
 391     fHistAnalysisSummary->SetBinContent(16, fForceBeamType);
 392     fHistAnalysisSummary->GetXaxis()->SetBinLabel(17, "fMinCent");
 393     fHistAnalysisSummary->SetBinContent(17, fMinCent);
 394     fHistAnalysisSummary->GetXaxis()->SetBinLabel(18, "fMaxCent");
 395     fHistAnalysisSummary->SetBinContent(18, fMaxCent);
 396     fHistAnalysisSummary->GetXaxis()->SetBinLabel(19, "fMinVz");
 397     fHistAnalysisSummary->SetBinContent(19, fMinVz);
 398     fHistAnalysisSummary->GetXaxis()->SetBinLabel(20, "fMaxVz");
 399     fHistAnalysisSummary->SetBinContent(20, fMaxVz);
 400     fHistAnalysisSummary->GetXaxis()->SetBinLabel(21, "fOffTrigger");
 401     fHistAnalysisSummary->SetBinContent(21, fOffTrigger);
 402     fHistAnalysisSummary->GetXaxis()->SetBinLabel(22, "fClusPtCut");
 403     fHistAnalysisSummary->SetBinContent(22, fClusPtCut);
 404     fHistAnalysisSummary->GetXaxis()->SetBinLabel(23, "fTrackPtCut");
 405     fHistAnalysisSummary->SetBinContent(23, fTrackPtCut);
 406     fHistAnalysisSummary->GetXaxis()->SetBinLabel(24, "fTrackMinEta");
 407     fHistAnalysisSummary->SetBinContent(24, fTrackMinEta);
 408     fHistAnalysisSummary->GetXaxis()->SetBinLabel(25, "fTrackMaxEta");
 409     fHistAnalysisSummary->SetBinContent(25, fTrackMaxEta);
 410     fHistAnalysisSummary->GetXaxis()->SetBinLabel(26, "fTrackMinPhi");
 411     fHistAnalysisSummary->SetBinContent(26, fTrackMinPhi);
 412     fHistAnalysisSummary->GetXaxis()->SetBinLabel(27, "fTrackMaxPhi");
 413     fHistAnalysisSummary->SetBinContent(27, fTrackMaxPhi);
 414     fHistAnalysisSummary->GetXaxis()->SetBinLabel(28, "fClusTimeCutLow");
 415     fHistAnalysisSummary->SetBinContent(28, fClusTimeCutLow);
 416     fHistAnalysisSummary->GetXaxis()->SetBinLabel(29, "fClusTimeCutUp");
 417     fHistAnalysisSummary->SetBinContent(29, fClusTimeCutUp);
 418     fHistAnalysisSummary->GetXaxis()->SetBinLabel(30, "fMinPtTrackInEmcal");
 419     fHistAnalysisSummary->SetBinContent(30, fMinPtTrackInEmcal);
 420     fHistAnalysisSummary->GetXaxis()->SetBinLabel(31, "fEventPlaneVsEmcal");
 421     fHistAnalysisSummary->SetBinContent(31, fEventPlaneVsEmcal);
 422     fHistAnalysisSummary->GetXaxis()->SetBinLabel(32, "fMinEventPlane");
 423     fHistAnalysisSummary->SetBinContent(32, fMaxEventPlane);
 424     fHistAnalysisSummary->GetXaxis()->SetBinLabel(33, "fRandomConeRadius");
 425     fHistAnalysisSummary->SetBinContent(33, fRandomConeRadius);
 426     fHistAnalysisSummary->GetXaxis()->SetBinLabel(34, "fitModulationType");
 427     fHistAnalysisSummary->SetBinContent(34, (int)fFitModulationType);
 428     fHistAnalysisSummary->GetXaxis()->SetBinLabel(35, "runModeType");
 429     fHistAnalysisSummary->SetBinContent(35, (int)fRunModeType);
 430     fHistAnalysisSummary->GetXaxis()->SetBinLabel(36, "data type");
 431     fHistAnalysisSummary->SetBinContent(36, (int)fDataType);
 432     fHistAnalysisSummary->GetXaxis()->SetBinLabel(37, "iterator");
 433     fHistAnalysisSummary->SetBinContent(37, 1.);
 434
 435     if(fFillQAHistograms) {
 436         fHistRunnumbersEta = new TH2F("fHistRunnumbersEta", "fHistRunnumbersEta", 100, -.5, 99.5, 100, -1.1, 1.1);
 437         fHistRunnumbersEta->Sumw2();
 438         fOutputList->Add(fHistRunnumbersEta);
 439         fHistRunnumbersPhi = new TH2F("fHistRunnumbersPhi", "fHistRunnumbersPhi", 100, -.5, 99.5, 100, -0.2, TMath::TwoPi()+0.2);
 440         fHistRunnumbersPhi->Sumw2();
 441         fOutputList->Add(fHistRunnumbersPhi);
 442     }
 443
 444     fHistSwap = new TH1F("fHistSwap", "fHistSwap", 20, 0, TMath::TwoPi());
 445     fHistSwap->Sumw2();
 446     fProfVn = new TProfile("fProfVn", "fProfVn", 2, -0.5, 1.5);
 447     fProfVn->GetXaxis()->SetBinLabel(1, "v_{2}(EBYE)");
 448     fProfVn->GetXaxis()->SetBinLabel(2, "v_{2}(EBYE)");
 449
 450     fOutputList->Add(fProfVn);
 451     PostData(1, fOutputList);
 452
 453     switch (fRunModeType) {
 454         case kLocal : {
 455             fOutputListGood = new TList();
 456             fOutputListGood->SetOwner(kTRUE);
 457             fOutputListBad = new TList();
 458             fOutputListBad->SetOwner(kTRUE);
 459             PostData(2, fOutputListGood);
 460             PostData(3, fOutputListBad);
 461         } break;
 462         default: break;
 463     }
 464
 465
 466 }
 467 //_____________________________________________________________________________
 468 Bool_t AliAnalysisTaskRhoVnModulation::Run()
 469 {
 470     // user exec: execute once for each event
 471     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 472     if(!fInitialized) fInitialized = InitializeAnalysis();
 473     // reject the event if expected data is missing
 474     if(!PassesCuts(InputEvent())) return kFALSE;
 475     if(!(fTracks||fJets||fRho)) return kFALSE;
 476     if(!fCaloClusters && fDebug > 0) printf(" > warning: couldn't retreive calo clusters! < \n");
 477
 478     // [0][0] psi2a     [1,0]   psi2c
 479     // [0][1] psi3a     [1,1]   psi3c
 480     Double_t vzero[2][2];
 481     CalculateEventPlaneVZERO(vzero);
 482     // [0] psi2         [1] psi3
 483     Double_t tpc[2];
 484     CalculateEventPlaneTPC(tpc);
 485
 486     // arrays which will hold the fit parameters
 487     Double_t fitParameters[] = {0,0,0,0,0,0,0,0,0};
 488     Double_t psi2(-1), psi3(-1);
 489     switch (fDetectorType) {    // determine the detector type for the rho fit
 490         case kTPC :     { psi2 = tpc[0];        psi3 = tpc[1]; }
 491         case kVZEROA :  { psi2 = vzero[0][0];   psi3 = vzero[0][1]; }
 492         case kVZEROC :  { psi2 = vzero[1][0];   psi3 = vzero[1][1]; }
 493         default : break;
 494     }
 495
 496     switch (fFitModulationType) { // do the fits
 497         case kNoFit : { fFitModulation->FixParameter(0, RhoVal()); } break;
 498         case kV2 : {
 499             CorrectRho(fitParameters, psi2, psi3);
 500             fProfVn->Fill((double)0, fFitModulation->GetParameter(3));
 501         } break;
 502         case kV3 : {
 503             CorrectRho(fitParameters, psi2, psi3);
 504             fProfVn->Fill((double)1, fFitModulation->GetParameter(3));
 505         } break;
 506         case kCombined : {
 507             CorrectRho(fitParameters, psi2, psi3);
 508             fProfVn->Fill((double)0, fFitModulation->GetParameter(3));
 509             fProfVn->Fill((double)1, fFitModulation->GetParameter(7));
 510         } break;
 511         default : break;
 512     }
 513     // fill a number of histograms
 514     FillHistogramsAfterSubtraction(vzero, tpc);
 515
 516     // send the output to the connected output container
 517     PostData(1, fOutputList);
 518     switch (fRunModeType) {
 519         case kLocal : {
 520             PostData(2, fOutputListGood);
 521             PostData(3, fOutputListBad);
 522         } break;
 523         default: break;
 524     }
 525     return kTRUE;
 526 }
 527 //_____________________________________________________________________________
 528 void AliAnalysisTaskRhoVnModulation::CalculateEventPlaneVZERO(Double_t vzero[2][2]) const
 529 {
 530     // grab the UNCALIBRATED vzero event plane
 531     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 532     Double_t qxa2(0), qya2(0), qxc2(0), qyc2(0);    // for psi2
 533     Double_t qxa3(0), qya3(0), qxc3(0), qyc3(0);    // for psi3
 534     for(Int_t iVZERO(0); iVZERO < 64; iVZERO++) {
 535         Double_t phi(TMath::PiOver4()*(.5+iVZERO%8)), /* eta(0), */ weight(InputEvent()->GetVZEROEqMultiplicity(iVZERO));
 536 //        (iVZERO<32) ? eta = -3.45+.5*(iVZERO/8) : eta = 4.8-.6*((iVZERO/8)-4);
 537         if(iVZERO<32) {
 538             qxa2 += weight*TMath::Cos(2.*phi);
 539             qya2 += weight*TMath::Sin(2.*phi);
 540             qxa3 += weight*TMath::Cos(3.*phi);
 541             qya3 += weight*TMath::Sin(3.*phi);
 542         }
 543         else {
 544             qxc2 += weight*TMath::Cos(2.*phi);
 545             qyc2 += weight*TMath::Sin(2.*phi);
 546             qxc3 += weight*TMath::Cos(3.*phi);
 547             qyc3 += weight*TMath::Sin(3.*phi);
 548        }
 549     }
 550     vzero[0][0] = .5*TMath::ATan2(qya2, qxa2);
 551     vzero[1][0] = .5*TMath::ATan2(qyc2, qxc2);
 552     vzero[0][1] = (1./3.)*TMath::ATan2(qya3, qxa3);
 553     vzero[1][1] = (1./3.)*TMath::ATan2(qyc3, qxc3);
 554 }
 555 //_____________________________________________________________________________
 556 void AliAnalysisTaskRhoVnModulation::CalculateEventPlaneTPC(Double_t* tpc) const
 557 {
 558    // grab the TPC event plane
 559    if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 560    Double_t qx2(0), qy2(0);     // for psi2
 561    Double_t qx3(0), qy3(0);     // for psi3
 562    if(fTracks) {
 563        Int_t iTracks(fTracks->GetEntriesFast());
 564        for(Int_t iTPC(0); iTPC < iTracks; iTPC++) {
 565            AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(iTPC));
 566            if(!PassesCuts(track)) continue;
 567            qx2+= TMath::Cos(2.*track->Phi());
 568            qy2+= TMath::Sin(2.*track->Phi());
 569            qx3+= TMath::Cos(3.*track->Phi());
 570            qy3+= TMath::Sin(3.*track->Phi());
 571        }
 572    }
 573    tpc[0] = .5*TMath::ATan2(qy2, qx2);
 574    tpc[1] = (1./3.)*TMath::ATan2(qy3, qx3);
 575 }
 576 //_____________________________________________________________________________
 577 void AliAnalysisTaskRhoVnModulation::CalculateRandomCone(Float_t &pt, Float_t &eta, Float_t &phi,
 578         AliEmcalJet* jet, Bool_t randomize) const
 579 {
 580     // get a random cone
 581     if(fDebug > 1) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 582     pt = 0; eta = 0; phi = 0;
 583     Float_t etaJet(999), phiJet(999), dJet(999);        // no jet: same as jet very far away
 584     if(jet) { // if a leading jet is given, use its kinematic properties
 585         etaJet = jet->Eta();
 586         phiJet = jet->Phi();
 587     }
 588     // force the random cones to at least be within detector acceptance
 589     Float_t minPhi(fJetMinPhi), maxPhi(fJetMaxPhi);
 590     if(maxPhi > TMath::TwoPi()) maxPhi = TMath::TwoPi();
 591     if(minPhi < 0 ) minPhi = 0;
 592     Float_t diffRcRJR(TMath::Abs(fRandomConeRadius-fJetRadius));
 593     // construct a random cone and see if it's far away enough from the leading jet
 594     Int_t attempts(1000);
 595     while(kTRUE) {
 596         attempts--;
 597         eta = gRandom->Uniform(fJetMinEta+diffRcRJR, fJetMaxEta-diffRcRJR);
 598         phi = gRandom->Uniform(minPhi, maxPhi);
 599
 600         dJet = TMath::Sqrt((etaJet-eta)*(etaJet-eta)+(phiJet-phi)*(phiJet-phi));
 601         if(dJet > fMinDisanceRCtoLJ) break;
 602         else if (attempts == 0) {
 603             printf(" > No random cone after 1000 tries, giving up ... !\n");
 604             return;
 605         }
 606     }
 607     if(fTracks) {
 608         Int_t iTracks(fTracks->GetEntriesFast());
 609         for(Int_t i(0); i < iTracks; i++) {
 610             AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(i));
 611             if(!PassesCuts(track)) continue;
 612             Float_t etaTrack(track->Eta()), phiTrack(track->Phi()), ptTrack(track->Pt());
 613             // if requested, randomize eta and phi to destroy any correlated fluctuations
 614             if(randomize) {
 615                 etaTrack = gRandom->Uniform(fTrackMinEta, fTrackMaxEta);
 616                 phiTrack = gRandom->Uniform(minPhi, maxPhi);
 617             }
 618             // get distance from cone
 619             if(TMath::Abs(phiTrack-phi) > TMath::Abs(phiTrack - phi + TMath::TwoPi())) phiTrack+=TMath::TwoPi();
 620             if(TMath::Abs(phiTrack-phi) > TMath::Abs(phiTrack - phi - TMath::TwoPi())) phiTrack-=TMath::TwoPi();
 621             if(TMath::Sqrt(TMath::Abs((etaTrack-eta)*(etaTrack-eta)+(phiTrack-phi)*(phiTrack-phi))) <= fRandomConeRadius) pt+=ptTrack;
 622         }
 623     }
 624 }
 625 //_____________________________________________________________________________
 626 void AliAnalysisTaskRhoVnModulation::CorrectRho(Double_t* params, Double_t psi2, Double_t psi3) const
 627 {
 628     // get rho' -> rho(phi)
 629     // the fit is constrained based on the switch of fFitModulationType which is called
 630     // in the initialization of this class.
 631     // after fitting, an array of fit parameters is set which can be re-created at
 632     // any point from the TF1 pointer fFitModulation
 633     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 634     TString detector("");
 635     switch (fDetectorType) {
 636         case kTPC : detector+="TPC";
 637             break;
 638         case kVZEROA : detector+="VZEROA";
 639             break;
 640         case kVZEROC : detector+="VZEROC";
 641             break;
 642         default: break;
 643     }
 644     Int_t iTracks(fTracks->GetEntriesFast());
 645     fHistSwap->Reset();     // clear the histogram
 646     for(Int_t i(0); i < iTracks; i++) {
 647             AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(i));
 648             if(!PassesCuts(track) || track->Pt() > 5 || track->Pt() < 0.15) continue;
 649             fHistSwap->Fill(track->Phi(), track->Pt());
 650     }
 651     fFitModulation->SetParameter(0, RhoVal());
 652     switch (fFitModulationType) {
 653         case kNoFit : { fFitModulation->FixParameter(0, RhoVal() ); }
 654         case kV2 : { fFitModulation->FixParameter(4, psi2); } break;
 655         case kV3 : { fFitModulation->FixParameter(4, psi3); } break;
 656         case kCombined : {
 657             fFitModulation->FixParameter(4, psi2);
 658             fFitModulation->FixParameter(6, psi3);
 659         } break;
 660         default : break;
 661     }
 662     fHistSwap->Fit(fFitModulation, fFitModulationOptions.Data(), "", 0, TMath::TwoPi());
 663     for(Int_t i(0); i < fFitModulation->GetNpar(); i++) params[i] = fFitModulation->GetParameter(i);
 664     // for LOCAL didactic purposes, save the  best and the worst fits
 665     // this routine can produce a lot of output histograms and will not work on GRID
 666     // since the output will become unmergeable
 667     switch (fRunModeType) {
 668         case kGrid : break;
 669         case kLocal : {
 670             static Int_t didacticCounterBest(0);
 671             static Int_t didacticCounterWorst(0);
 672             static Double_t bestFitP(.05);      // threshold for significance
 673             static Double_t worstFitP(.05);
 674             Double_t p(ChiSquare(fFitModulation->GetNDF(), fFitModulation->GetChisquare()));
 675             if(p > bestFitP || p > 0.12) {
 676                 TProfile* didacticProfile = (TProfile*)fHistSwap->Clone(Form("Fit_%i_p_%.3f_cen_%i_%s", didacticCounterBest, p, fInCentralitySelection, detector.Data()));
 677                 fOutputListGood->Add(didacticProfile);
 678                 didacticCounterBest++;
 679                 bestFitP = p;
 680              }
 681              else if(p < worstFitP) {
 682                 TProfile* didacticProfile = (TProfile*)fHistSwap->Clone(Form("Fit_%i_p_%.3f_cen_%i_%s", didacticCounterWorst, p, fInCentralitySelection, detector.Data() ));
 683                 fOutputListBad->Add(didacticProfile);
 684                 didacticCounterWorst++;
 685                 worstFitP = p;
 686              }
 687          } break;
 688          default : break;
 689     }
 690 }
 691 //_____________________________________________________________________________
 692 Bool_t AliAnalysisTaskRhoVnModulation::PassesCuts(AliVEvent* event)
 693 {
 694     // event cuts
 695     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 696     if(!event) return kFALSE;
 697     if(TMath::Abs(InputEvent()->GetPrimaryVertex()->GetZ()) > 10.) return kFALSE;
 698     // aod and esd specific checks
 699     switch (fDataType) {
 700        case kESD: {
 701             AliESDEvent* esdEvent = static_cast<AliESDEvent*>(InputEvent());
 702             if( (!esdEvent) || (TMath::Abs(esdEvent->GetPrimaryVertexSPD()->GetZ() - esdEvent->GetPrimaryVertex()->GetZ()) > .5) ) return kFALSE;
 703        } break;
 704        case kAOD: {
 705             AliAODEvent* aodEvent = static_cast<AliAODEvent*>(InputEvent());
 706             if( (!aodEvent) || (TMath::Abs(aodEvent->GetPrimaryVertexSPD()->GetZ() - aodEvent->GetPrimaryVertex()->GetZ()) > .5) ) return kFALSE;
 707        } break;
 708        default: break;
 709     }
 710     fCent = InputEvent()->GetCentrality()->GetCentralityPercentile("V0M");
 711     if(fCent <= 0 || fCent >= 100 || TMath::Abs(fCent-InputEvent()->GetCentrality()->GetCentralityPercentile("TRK")) > 5.) return kFALSE;
 712     fInCentralitySelection = TMath::FloorNint(fCent/10.);
 713     if(fFillQAHistograms) FillQAHistograms(event);
 714     return kTRUE;
 715 }
 716 //_____________________________________________________________________________
 717 Bool_t AliAnalysisTaskRhoVnModulation::PassesCuts(const AliVCluster* cluster) const
 718 {
 719     // cluster cuts
 720     if(fDebug > 1) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 721     if(!cluster) return kFALSE;
 722     return kTRUE;
 723 }
 724 //_____________________________________________________________________________
 725 Bool_t AliAnalysisTaskRhoVnModulation::PassesCuts(const AliEmcalJet* jet) const
 726 {
 727     // jet cuts
 728     if(fDebug > 1) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 729     if(!jet) return kFALSE;
 730     if(TMath::Abs(jet->Eta()) > .5) return kFALSE;
 731     return kTRUE;
 732 }
 733 //_____________________________________________________________________________
 734 void AliAnalysisTaskRhoVnModulation::FillHistogramsAfterSubtraction(Double_t vzero[2][2], Double_t* tpc) const
 735 {
 736     // fill histograms
 737     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 738     FillTrackHistograms();
 739     /* FillClusterHistograms(); */
 740     FillJetHistograms(vzero, tpc);
 741     /* FillCorrectedClusterHistograms(); */
 742     FillEventPlaneHistograms(vzero, tpc);
 743     FillRhoHistograms();
 744     FillDeltaPtHistograms();
 745     FillDeltaPhiHistograms(vzero, tpc);
 746 }
 747 //_____________________________________________________________________________
 748 void AliAnalysisTaskRhoVnModulation::FillTrackHistograms() const
 749 {
 750     // fill track histograms
 751     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 752     Int_t iTracks(fTracks->GetEntriesFast());
 753     for(Int_t i(0); i < iTracks; i++) {
 754         AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(i));
 755         if(!PassesCuts(track)) continue;
 756         fHistPicoTrackPt[fInCentralitySelection]->Fill(track->Pt());
 757         fHistPicoTrackEta[fInCentralitySelection]->Fill(track->Eta());
 758         fHistPicoTrackPhi[fInCentralitySelection]->Fill(track->Phi());
 759         if(fFillQAHistograms) FillQAHistograms(track);
 760     }
 761     return;
 762 }
 763 //_____________________________________________________________________________
 764 void AliAnalysisTaskRhoVnModulation::FillClusterHistograms() const
 765 {
 766     // fill cluster histograms
 767     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 768     /* Int_t iClusters(fCaloClusters->GetEntriesFast());
 769     for(Int_t i(0); i < iClusters; i++) {
 770         AliVCluster* cluster = static_cast<AliVCluster*>(fCaloClusters->At(iClusters));
 771         if (!PassesCuts(cluster)) continue;
 772         TLorentzVector clusterLorentzVector;
 773         cluster->GetMomentum(clusterLorentzVector, const_cast<Double_t*>(fVertex));
 774         fHistClusterPt[fInCentralitySelection]->Fill(clusterLorentzVector.Pt());
 775         fHistClusterEta[fInCentralitySelection]->Fill(clusterLorentzVector.Eta());
 776         fHistClusterPhi[fInCentralitySelection]->Fill(clusterLorentzVector.Phi());
 777     }
 778     return; */
 779 }
 780 //_____________________________________________________________________________
 781 void AliAnalysisTaskRhoVnModulation::FillCorrectedClusterHistograms() const
 782 {
 783     // fill clusters after hadronic correction FIXME implement
 784     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 785     return;
 786 }
 787 //_____________________________________________________________________________
 788 void AliAnalysisTaskRhoVnModulation::FillEventPlaneHistograms(Double_t vzero[2][2], Double_t* tpc) const
 789 {
 790     // fill event plane histograms
 791     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 792     fHistPsi2->Fill(0.5, vzero[0][0]);
 793     fHistPsi2->Fill(1.5, vzero[1][0]);
 794     fHistPsi2->Fill(2.5, tpc[0]);
 795     fHistPsiVZEROA->Fill(vzero[0][0]);
 796     fHistPsiVZEROC->Fill(vzero[1][0]);
 797     fHistPsiTPC->Fill(tpc[0]);
 798     fHistPsi2Spread->Fill(0.5, vzero[0][0]-vzero[1][0]);
 799     fHistPsi2Spread->Fill(1.5, vzero[0][0]-tpc[0]);
 800     fHistPsi2Spread->Fill(2.5, vzero[1][0]-tpc[0]);
 801     return;
 802 }
 803 //_____________________________________________________________________________
 804 void AliAnalysisTaskRhoVnModulation::FillRhoHistograms() const
 805 {
 806     // fill rho histograms
 807     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 808     fHistRhoPackage[fInCentralitySelection]->Fill(RhoVal());    // save the rho estimate from the emcal jet package
 809     // get multiplicity FIXME inefficient
 810     Int_t iTracks(fTracks->GetEntriesFast()), mult(0), iJets(fJets->GetEntriesFast());
 811     for(Int_t i(0); i < iTracks; i ++) { if(PassesCuts(static_cast<AliVTrack*>(fTracks->At(i)))) mult++; }
 812     Double_t rho(RhoVal(TMath::Pi(), TMath::Pi(), fRho->GetVal()));
 813     fHistRho[fInCentralitySelection]->Fill(rho);
 814     fHistRhoVsMult->Fill(mult, rho);
 815     fHistRhoVsCent->Fill(fCent, rho);
 816     for(Int_t i(0); i < iJets; i++) {
 817         AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(i));
 818         if(!PassesCuts(jet)) continue;
 819         fHistRhoAVsMult->Fill(mult, rho * jet->Area());
 820         fHistRhoAVsCent->Fill(fCent, rho * jet->Area());
 821     }
 822     return;
 823 }
 824 //_____________________________________________________________________________
 825 void AliAnalysisTaskRhoVnModulation::FillDeltaPtHistograms() const
 826 {
 827     // fill delta pt histograms
 828     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 829     static Int_t sJets[9999] = {-1};
 830     GetSortedArray(sJets, fJets);
 831 //    if(sJets[0] <= 0) return;
 832     Int_t i(0), maxCones(20);
 833     AliEmcalJet* leadingJet(0x0);
 834     do { // get the leading jet
 835         leadingJet = static_cast<AliEmcalJet*>(fJets->At(sJets[i]));
 836         i++;
 837     }
 838     while (!PassesCuts(leadingJet)&&i<fJets->GetEntriesFast());
 839     if(!leadingJet && fDebug > 0) printf(" > failed to retrieve leading jet ! < \n");
 840     const Float_t areaRC = fRandomConeRadius*fRandomConeRadius*TMath::Pi();
 841     // we're retrieved the leading jet, now get a random cone
 842     for(i = 0; i < maxCones; i++) {
 843        Float_t pt(0), eta(0), phi(0);
 844        // get a random cone without constraints on leading jet position
 845        CalculateRandomCone(pt, eta, phi, 0x0);
 846        if(pt > 0) {
 847            fHistRCPhiEta[fInCentralitySelection]->Fill(phi, eta);
 848            fHistRhoVsRCPt[fInCentralitySelection]->Fill(pt, RhoVal(phi, fJetRadius, fRho->GetVal())*areaRC);
 849            fHistRCPt[fInCentralitySelection]->Fill(pt);
 850            fHistDeltaPtRC[fInCentralitySelection]->Fill(pt - areaRC*RhoVal(phi, fJetRadius, fRho->GetVal()));
 851        }
 852        // get a random cone excluding leading jet area
 853        CalculateRandomCone(pt, eta, phi, leadingJet);
 854        if(pt > 0) {
 855            fHistRCPhiEtaExLJ[fInCentralitySelection]->Fill(phi, eta);
 856            fHistRhoVsRCPtExLJ[fInCentralitySelection]->Fill(pt, RhoVal(phi, fJetRadius, fRho->GetVal())*areaRC);
 857            fHistRCPtExLJ[fInCentralitySelection]->Fill(pt);
 858            fHistDeltaPtRCExLJ[fInCentralitySelection]->Fill(pt - areaRC*RhoVal(phi, fJetRadius, fRho->GetVal()));
 859        }
 860        // get a random cone in an event with randomized phi and eta
 861        CalculateRandomCone(pt, eta, phi, 0x0, kTRUE);
 862        if( pt > 0) {
 863            fHistRCPhiEtaRand[fInCentralitySelection]->Fill(phi, eta);
 864            fHistRhoVsRCPtRand[fInCentralitySelection]->Fill(pt, RhoVal(phi, fJetRadius, fRho->GetVal())*areaRC);
 865            fHistRCPtRand[fInCentralitySelection]->Fill(pt);
 866            fHistDeltaPtRCRand[fInCentralitySelection]->Fill(pt - areaRC*RhoVal(phi, fJetRadius, fRho->GetVal()));
 867        }
 868     }
 869 }
 870 //_____________________________________________________________________________
 871 void AliAnalysisTaskRhoVnModulation::FillJetHistograms(Double_t vzero[2][2], Double_t* tpc) const
 872 {
 873     // fill jet histograms
 874     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 875     Int_t iJets(fJets->GetEntriesFast());
 876     for(Int_t i(0); i < iJets; i++) {
 877         AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(i));
 878         if(!PassesCuts(jet)) continue;
 879         Double_t pt(jet->Pt()), area(jet->Area()), eta(jet->Eta()), phi(jet->Phi());
 880         Double_t rho(RhoVal(phi, fJetRadius, fRho->GetVal()));
 881         fHistJetPtRaw[fInCentralitySelection]->Fill(pt);
 882         fHistJetPt[fInCentralitySelection]->Fill(pt-area*rho);
 883         fHistJetEtaPhi[fInCentralitySelection]->Fill(eta, phi);
 884         fHistJetPtArea[fInCentralitySelection]->Fill(pt-area*rho, area);
 885         Double_t dPhiA = PhaseShift(phi-vzero[0][0]);
 886         Double_t dPhiC = PhaseShift(phi-vzero[1][0]);
 887         Double_t dPhiTPC = PhaseShift(phi-tpc[0]);
 888         Double_t PiE = TMath::PiOver4()/2.;
 889         fHistJetPsiTPCPt[fInCentralitySelection]->Fill(dPhiTPC, pt-area*rho);
 890         fHistJetPsiVZEROAPt[fInCentralitySelection]->Fill(dPhiA, pt-area*rho);
 891         fHistJetPsiVZEROCPt[fInCentralitySelection]->Fill(dPhiC, pt-area*rho);
 892         if(dPhiA > 15.*PiE && dPhiA < PiE) fHistJetPtInPlaneVZEROA[fInCentralitySelection]->Fill(pt-area*rho);
 893         else if(dPhiA > 7.*PiE && dPhiA < 9.*PiE) fHistJetPtInPlaneVZEROA[fInCentralitySelection]->Fill(pt-area*rho);
 894         else if(dPhiA > 3.*PiE && dPhiA < 5.*PiE) fHistJetPtOutPlaneVZEROA[fInCentralitySelection]->Fill(pt-area*rho);
 895         else if(dPhiA > 11.*PiE && dPhiA < 13.*PiE) fHistJetPtOutPlaneVZEROA[fInCentralitySelection]->Fill(pt-area*rho);
 896         else fHistJetPtMidPlaneVZEROA[fInCentralitySelection]->Fill(pt-area*rho);
 897         if(dPhiC > 15.*PiE && dPhiC < PiE) fHistJetPtInPlaneVZEROC[fInCentralitySelection]->Fill(pt-area*rho);
 898         else if(dPhiC > 7.*PiE && dPhiC < 9.*PiE) fHistJetPtInPlaneVZEROC[fInCentralitySelection]->Fill(pt-area*rho);
 899         else if(dPhiC > 3.*PiE && dPhiC < 5.*PiE) fHistJetPtOutPlaneVZEROC[fInCentralitySelection]->Fill(pt-area*rho);
 900         else if(dPhiC > 11.*PiE && dPhiC < 13.*PiE) fHistJetPtOutPlaneVZEROC[fInCentralitySelection]->Fill(pt-area*rho);
 901         else fHistJetPtMidPlaneVZEROC[fInCentralitySelection]->Fill(pt-area*rho);
 902         if(dPhiTPC > 15.*PiE && dPhiTPC < PiE) fHistJetPtInPlaneTPC[fInCentralitySelection]->Fill(pt-area*rho);
 903         else if(dPhiTPC > 7.*PiE && dPhiTPC < 9.*PiE) fHistJetPtInPlaneTPC[fInCentralitySelection]->Fill(pt-area*rho);
 904         else if(dPhiTPC > 3.*PiE && dPhiTPC < 5.*PiE) fHistJetPtOutPlaneTPC[fInCentralitySelection]->Fill(pt-area*rho);
 905         else if(dPhiTPC > 11.*PiE && dPhiTPC < 13.*PiE) fHistJetPtOutPlaneTPC[fInCentralitySelection]->Fill(pt-area*rho);
 906         else fHistJetPtMidPlaneTPC[fInCentralitySelection]->Fill(pt-area*rho);
 907         // last but not least, set the subtracted pt
 908         jet->SetPtSub(jet->PtSub(rho));
 909         fHistJetPtConstituents[fInCentralitySelection]->Fill(jet->PtSub(), jet->Nch());
 910     }
 911 }
 912 //_____________________________________________________________________________
 913 void AliAnalysisTaskRhoVnModulation::FillDeltaPhiHistograms(Double_t vzero[2][2], Double_t* tpc) const
 914 {
 915    // fill phi minus psi histograms
 916    if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 917    if(fTracks) {
 918        Int_t iTracks(fTracks->GetEntriesFast());
 919        for(Int_t iTPC(0); iTPC < iTracks; iTPC++) {
 920            AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(iTPC));
 921            if(!PassesCuts(track)) continue;
 922            fHistDeltaPhiVZEROA[fInCentralitySelection]->Fill(PhaseShift(track->Phi()-vzero[0][0]));
 923            fHistDeltaPhiVZEROC[fInCentralitySelection]->Fill(PhaseShift(track->Phi()-vzero[1][0]));
 924            fHistDeltaPhiTPC[fInCentralitySelection]->Fill(PhaseShift(track->Phi()-tpc[0]));
 925        }
 926    }
 927 }
 928 //_____________________________________________________________________________
 929 void AliAnalysisTaskRhoVnModulation::FillQAHistograms(AliVTrack* vtrack) const
 930 {
 931     // fill qa histograms for pico tracks
 932     if(!vtrack) return;
 933     AliPicoTrack* track = static_cast<AliPicoTrack*>(vtrack);
 934     fHistRunnumbersPhi->Fill(fMappedRunNumber, track->Phi());
 935     fHistRunnumbersEta->Fill(fMappedRunNumber, track->Eta());
 936     if(!track) return;
 937     Int_t type((int)(track->GetTrackType()));
 938     switch (type) {
 939         case 0:
 940            fHistPicoCat1[fInCentralitySelection]->Fill(track->Eta(), track->Phi());
 941            break;
 942         case 1:
 943            fHistPicoCat2[fInCentralitySelection]->Fill(track->Eta(), track->Phi());
 944            break;
 945         case 2:
 946            fHistPicoCat3[fInCentralitySelection]->Fill(track->Eta(), track->Phi());
 947            break;
 948         default: break;
 949     }
 950 }
 951 //_____________________________________________________________________________
 952 void AliAnalysisTaskRhoVnModulation::FillQAHistograms(AliVEvent* vevent)
 953 {
 954     // fill qa histograms for events
 955     if(!vevent) return;
 956     fHistVertexz->Fill(vevent->GetPrimaryVertex()->GetZ());
 957     fHistCentrality->Fill(fCent);
 958     Int_t runNumber(InputEvent()->GetRunNumber());
 959     Int_t runs[] = {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, 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};
 960     for(fMappedRunNumber = 0; fMappedRunNumber < 65; fMappedRunNumber++) {
 961         if(runs[fMappedRunNumber]==runNumber) break;
 962     }
 963 }
 964 //_____________________________________________________________________________
 965 void AliAnalysisTaskRhoVnModulation::Terminate(Option_t *)
 966 {
 967     // terminate
 968     switch (fRunModeType) {
 969         case kLocal : {
 970         printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 971         if(fFillQAHistograms) {
 972             Int_t runs[] = {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, 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};
 973             for(Int_t i(0); i < 64; i++) {
 974                 fHistRunnumbersPhi->GetXaxis()->SetBinLabel(i+1, Form("%i", runs[i]));
 975                 fHistRunnumbersEta->GetXaxis()->SetBinLabel(i+1, Form("%i", runs[i]));
 976             }
 977             fHistRunnumbersPhi->GetXaxis()->SetBinLabel(65, "undetermined");
 978             fHistRunnumbersEta->GetXaxis()->SetBinLabel(65, "undetermined");
 979         }
 980         AliAnalysisTaskRhoVnModulation::Dump();
 981         for(Int_t i(0); i < fHistAnalysisSummary->GetXaxis()->GetNbins(); i++) printf( " > flag: %s \t content %.2f \n", fHistAnalysisSummary->GetXaxis()->GetBinLabel(1+i), fHistAnalysisSummary->GetBinContent(1+i));
 982         } break;
 983         default : break;
 984     }
 985 }
 986 //_____________________________________________________________________________