PWGJE/UserTasks/AliAnalysisTaskJetProtonCorr.cxx

   1 // ROOT
   2 #include "TFile.h"
   3 #include "TList.h"
   4 #include "TH1.h"
   5 #include "TH2.h"
   6 #include "TH3.h"
   7 #include "TF1.h"
   8 #include "TFormula.h"
   9 #include "TRandom.h"
  10 #include "TSpline.h"
  11
  12 // analysis framework
  13 #include "AliAnalysisManager.h"
  14 #include "AliInputEventHandler.h"
  15 #include "AliVEvent.h"
  16 #include "AliVTrack.h"
  17 #include "AliVTrdTrack.h"
  18 #include "AliVVertex.h"
  19 #include "AliPIDResponse.h"
  20 #include "AliEventPoolManager.h"
  21 #include "AliOADBContainer.h"
  22 #include "AliTOFPIDParams.h"
  23 #include "AliAnalysisTaskVnV0.h"
  24
  25 // MC stuff
  26 #include "AliMCEvent.h"
  27 #include "AliGenPythiaEventHeader.h"
  28
  29 // ESD stuff
  30 #include "AliESDEvent.h"
  31 #include "AliESDInputHandler.h"
  32 #include "AliESDtrack.h"
  33 #include "AliESDtrackCuts.h"
  34 #include "AliESDTrdTrack.h"
  35 #include "AliESDTrdTracklet.h"
  36 #include "AliESDTrdTrigger.h"
  37
  38 // AOD stuff
  39 #include "AliAODEvent.h"
  40 #include "AliAODJet.h"
  41 #include "AliAODTrack.h"
  42
  43 // jet tasks
  44 #include "AliAnalysisTaskJetServices.h"
  45 #include "AliAnalysisHelperJetTasks.h"
  46
  47 #include "AliAnalysisTaskJetProtonCorr.h"
  48
  49 #include <iostream>
  50 #include <cmath>
  51
  52 AliAnalysisTaskJetProtonCorr::AliAnalysisTaskJetProtonCorr(const char *name) :
  53   AliAnalysisTaskSE(name),
  54   fMCEvent(0x0),
  55   fESDEvent(0x0),
  56   fAODEvent(0x0),
  57   fRunNumber(-1),
  58   fOADBContainerTOF(0x0),
  59   fParamsTOF(0x0),
  60   fEventplane(0x0),
  61   fTriggerMask(0),
  62   fClassMask(0),
  63   fCentrality(100.),
  64   fCentralityCheck(100.),
  65   fZvtx(0.),
  66   fPIDResponse(0x0),
  67   fEventPlaneAngle(5.),
  68   fEventPlaneRes(0.),
  69   fEventPlaneAngleCheck(5.),
  70   fEventPlaneAngle3(5.),
  71   fPrimTrackArrayAss(0x0),
  72   fPrimTrackArrayTrg(0x0),
  73   fPrimConstrainedTrackArray(new TClonesArray("AliESDtrack", 100)),
  74   fJetArray(0x0),
  75   fPoolMgr(),
  76   fPool(),
  77   fHistCorr(0x0),
  78   fErrorMsg(10),
  79   fOutputList(),
  80   fHist(),
  81   fShortTaskId("jet_prot_corr"),
  82   fUseStandardCuts(kTRUE),
  83   fUseEvplaneV0(kFALSE),
  84   fCutsPrimTrg(0x0),
  85   fCutsPrimTrgConstrain(new AliESDtrackCuts()),
  86   fCutsPrimAss(0x0),
  87   fCutsTwoTrackEff(0.02),
  88   fAssFilterMask(1 << 10),
  89   fRequirePID(kTRUE),
  90   fTrgJetEtaMax(0.45),
  91   fHadEtaMax(0.8),
  92   fTrgPartPtMin(6.),
  93   fTrgPartPtMax(8.),
  94   fTrgJetPtMin(50.),
  95   fTrgJetPtMax(80.),
  96   fTrgJetLeadTrkPtMin(6.),
  97   fTrgJetLeadTrkPtMax(100.),
  98   fTrgJetAreaMin(0.6 * TMath::Pi() * 0.2*0.2),
  99   fAssPartPtMin(2.),
 100   fAssPartPtMax(4.),
 101   fTrgAngleToEvPlane(TMath::Pi() / 4.),
 102   fToyMeanNoPart(.5),
 103   fToyRadius(.8),
 104   fToySmearPhi(.2),
 105   fTrgJetPhiModCent(new TF1("jetphimodcent", "1 + 2 * [0] * cos(2*x)", 0., 2 * TMath::Pi())),
 106   fTrgJetPhiModSemi(new TF1("jetphimodsemi", "1 + 2 * [0] * cos(2*x)", 0., 2 * TMath::Pi())),
 107   fTrgHadPhiModCent(new TF1("hadphimodcent", "1 + 2 * [0] * cos(2*x)", 0., 2 * TMath::Pi())),
 108   fTrgHadPhiModSemi(new TF1("hadphimodsemi", "1 + 2 * [0] * cos(2*x)", 0., 2 * TMath::Pi())),
 109   fTrgJetV2Cent(0.10),
 110   fTrgJetV2Semi(0.10),
 111   fTrgHadV2Cent(0.04),
 112   fTrgHadV2Semi(0.10),
 113   fSplineEventPlaneRes(0x0)
 114 {
 115   // default ctor
 116
 117   fkCorrTypeName[kCorrHadHad]  = "hh";
 118   fkCorrTypeName[kCorrHadProt] = "hp";
 119   fkCorrTypeName[kCorrJetHad]  = "jh";
 120   fkCorrTypeName[kCorrJetProt] = "jp";
 121   fkCorrTypeName[kCorrRndJetHad]  = "rjh";
 122   fkCorrTypeName[kCorrRndJetProt]  = "rjp";
 123   fkCorrTypeName[kCorrRndHadHad]  = "rhh";
 124   fkCorrTypeName[kCorrRndHadProt]  = "rhp";
 125   fkCorrTypeName[kCorrRndJetExcHad]  = "rjeh";
 126   fkCorrTypeName[kCorrRndJetExcProt]  = "rjep";
 127   fkCorrTypeName[kCorrRndHadExcHad]  = "rheh";
 128   fkCorrTypeName[kCorrRndHadExcProt]  = "rhep";
 129
 130   fkClassName[kClCentral]      = "cent";
 131   fkClassName[kClSemiCentral]  = "semi";
 132   // fkClassName[kClDijet]        = "dijet";
 133
 134   fkEvName[kEvSame] = "same";
 135   fkEvName[kEvMix]  = "mixed";
 136
 137   // track cuts for associates
 138   if (fUseStandardCuts) {
 139     fCutsPrimAss = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE);
 140   } else {
 141     fCutsPrimAss = new AliESDtrackCuts();
 142
 143     // this is taken from PWGJE track cuts
 144     TFormula *f1NClustersTPCLinearPtDep = new TFormula("f1NClustersTPCLinearPtDep","70.+30./20.*x");
 145     fCutsPrimAss->SetMinNClustersTPCPtDep(f1NClustersTPCLinearPtDep,20.);
 146     fCutsPrimAss->SetMinNClustersTPC(70);
 147     fCutsPrimAss->SetMaxChi2PerClusterTPC(4);
 148     fCutsPrimAss->SetRequireTPCStandAlone(kTRUE); //cut on NClustersTPC and chi2TPC Iter1
 149     fCutsPrimAss->SetAcceptKinkDaughters(kFALSE);
 150     fCutsPrimAss->SetRequireTPCRefit(kTRUE);
 151     fCutsPrimAss->SetMaxFractionSharedTPCClusters(0.4);
 152     // ITS
 153     fCutsPrimAss->SetRequireITSRefit(kTRUE);
 154     //accept secondaries
 155     fCutsPrimAss->SetMaxDCAToVertexXY(2.4);
 156     fCutsPrimAss->SetMaxDCAToVertexZ(3.2);
 157     fCutsPrimAss->SetDCAToVertex2D(kTRUE);
 158     //reject fakes
 159     fCutsPrimAss->SetMaxChi2PerClusterITS(36);
 160     fCutsPrimAss->SetMaxChi2TPCConstrainedGlobal(36);
 161
 162     fCutsPrimAss->SetRequireSigmaToVertex(kFALSE);
 163
 164     fCutsPrimAss->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
 165   }
 166
 167   fCutsPrimAss->SetEtaRange(-fHadEtaMax, fHadEtaMax);
 168   fCutsPrimAss->SetPtRange(0.15, 1E+15);
 169
 170   // track cuts for triggers
 171   fCutsPrimTrg = new AliESDtrackCuts(*fCutsPrimAss);
 172
 173   // azimuthal modulation for triggers
 174   fTrgJetPhiModCent->SetParameter(0, fTrgJetV2Cent);
 175   fTrgJetPhiModSemi->SetParameter(0, fTrgJetV2Semi);
 176   fTrgHadPhiModCent->SetParameter(0, fTrgHadV2Cent);
 177   fTrgHadPhiModSemi->SetParameter(0, fTrgHadV2Semi);
 178
 179   // event mixing pool
 180   Double_t centralityBins[] = {
 181     0., 2., 4., 6., 8., 10., // central
 182     30., 32., 34., 36., 38., 40., 42., 44., 46., 48., 50., // semi-central
 183     90.
 184   };
 185   Int_t nCentralityBins = sizeof(centralityBins)/sizeof(centralityBins[0]);
 186
 187   Double_t vertexBins[] = {
 188     -10., -8., -6., -4., -2., 0., 2., 4., 6., 8., 10.
 189   };
 190   Int_t nVertexBins = sizeof(vertexBins)/sizeof(vertexBins[0]);
 191
 192   Double_t psiBins[12];
 193   Int_t nPsiBins = sizeof(psiBins)/sizeof(psiBins[0]);
 194   for (Int_t iBin = 0; iBin < nPsiBins; ++iBin)
 195     psiBins[iBin] = iBin * TMath::Pi()/nPsiBins;
 196
 197   const Int_t poolSize = 10; // unused by current event pool implementation
 198   const Int_t trackDepth = 50000; // number of tracks to maintain in mixing buffer
 199   const Float_t trackDepthFraction = 0.1;
 200   const Int_t targetEvents = 1;
 201   for (Int_t iAss = 0; iAss <= kAssProt; ++iAss) {
 202     GetPoolMgr((Ass_t) iAss) =
 203       new AliEventPoolManager(poolSize, trackDepth,
 204                               nCentralityBins, centralityBins,
 205                               nVertexBins, vertexBins);
 206     GetPoolMgr((Ass_t) iAss)->SetTargetValues(trackDepth, trackDepthFraction, targetEvents);
 207   }
 208
 209   fHistCorr = new AliHistCorr*[kEvLast*kCorrLast*kClLast];
 210
 211   DefineOutput(1, TList::Class());
 212 }
 213
 214 AliAnalysisTaskJetProtonCorr::~AliAnalysisTaskJetProtonCorr()
 215 {
 216   // dtor
 217
 218   // delete [] fHistCorr;
 219 }
 220
 221 void AliAnalysisTaskJetProtonCorr::UserCreateOutputObjects()
 222 {
 223   // create user output objects
 224
 225   // open OADB file and retrieve the OADB container with TOF parameters
 226   TString oadbFileName =
 227     TString::Format("%s/COMMON/PID/data/TOFPIDParams.root", AliAnalysisManager::GetOADBPath());
 228   TFile *oadbFile = TFile::Open(oadbFileName);
 229   if(!oadbFile->IsOpen())
 230     AliFatal(Form("Cannot open OADB file %s", oadbFileName.Data()));
 231   AliOADBContainer *oadbContainer =
 232     (AliOADBContainer*) oadbFile->Get("TOFoadb");
 233   if (!oadbContainer)
 234     AliFatal("Cannot fetch OADB container for VZERO EP selection");
 235   fOADBContainerTOF = new AliOADBContainer(*oadbContainer);
 236   oadbFile->Close();
 237   delete oadbFile;
 238
 239   // setup list
 240   OpenFile(1);
 241   fOutputList = new TList();
 242   fOutputList->SetOwner();
 243
 244   // setup histograms
 245   TH1 *hist;
 246   TH1 *histStat = AddHistogram(ID(kHistStat), "event statistics;;counts",
 247                                kStatLast-1, .5, kStatLast-.5);
 248   histStat->GetXaxis()->SetBinLabel(ID(kStatSeen));
 249   histStat->GetXaxis()->SetBinLabel(ID(kStatTrg));
 250   histStat->GetXaxis()->SetBinLabel(ID(kStatVtx));
 251   histStat->GetXaxis()->SetBinLabel(ID(kStatEvCuts));
 252   histStat->GetXaxis()->SetBinLabel(ID(kStatUsed));
 253   histStat->GetXaxis()->SetBinLabel(ID(kStatCent));
 254   histStat->GetXaxis()->SetBinLabel(ID(kStatEvPlane));
 255   histStat->GetXaxis()->SetBinLabel(ID(kStatPID));
 256   histStat->GetXaxis()->SetBinLabel(ID(kStatCentral));
 257   histStat->GetXaxis()->SetBinLabel(ID(kStatSemiCentral));
 258
 259   AddHistogram(ID(kHistVertexNctb), "number of vertex contributors;N_{ctb};counts",
 260                100, 0., 2000.);
 261   AddHistogram(ID(kHistVertexZ), "z-position of primary vertex;z (cm);counts",
 262                100, -50., 50.);
 263
 264   AddHistogram(ID(kHistCentrality), "centrality;C;counts",
 265                110, -5., 105.);
 266   hist = AddHistogram(ID(kHistCentralityUsed), "centrality used;C;event class",
 267                       110, -5., 105.,
 268                       kClLast, -.5, kClLast-.5);
 269   hist->GetYaxis()->SetBinLabel(LAB(kClCentral));
 270   hist->GetYaxis()->SetBinLabel(LAB(kClSemiCentral));
 271   // hist->GetYaxis()->SetBinLabel(LAB(kClDijet));
 272   AddHistogram(ID(kHistCentralityCheck), "centrality check;C;counts",
 273                110, -5., 105.);
 274   hist = AddHistogram(ID(kHistCentralityCheckUsed), "centrality check used;C;event class",
 275                       110, -5., 105.,
 276                       kClLast, -.5, kClLast-.5);
 277   hist->GetYaxis()->SetBinLabel(LAB(kClCentral));
 278   hist->GetYaxis()->SetBinLabel(LAB(kClSemiCentral));
 279   // hist->GetYaxis()->SetBinLabel(LAB(kClDijet));
 280   AddHistogram(ID(kHistCentralityVsMult), "centrality - multiplicity;centrality percentile (%);N_{prim}",
 281                100, 0., 100.,
 282                100, 0., 2500.);
 283
 284   AddHistogram(ID(kHistSignalTPC), "TPC dE/dx;p (GeV/#it{c});dE/dx (arb. units)",
 285                100, 0., 10., 200, 0., 300.);
 286   AddHistogram(ID(kHistSignalTOF), "TOF time of flight;p_{T} (GeV/#it{c});t (ns)",
 287                100, 0., 10., 200, 0., 50.);
 288   AddHistogram(ID(kHistBetaTOF), "TOF beta;p (GeV/#it{c}); #beta",
 289                100, 0., 10.,
 290                100, 0., 1.);
 291
 292   AddHistogram(ID(kHistDeltaTPC), "TPC dE/dx;p (GeV/#it{c});dE/dx (arb. units)",
 293                100, 0., 10., 200, -100., 100.);
 294   AddHistogram(ID(kHistDeltaTPCSemi), "TPC dE/dx;p (GeV/#it{c});dE/dx (arb. units)",
 295                100, 0., 10., 200, -100., 100.);
 296   AddHistogram(ID(kHistDeltaTOF), "TOF time of flight;p (GeV/#it{c});t (ns)",
 297                100, 0., 10., 200, -2., 2.);
 298   AddHistogram(ID(kHistDeltaTOFSemi), "TOF time of flight;p (GeV/#it{c});t (ns)",
 299                100, 0., 10., 200, -2., 2.);
 300
 301   // Nsigma templates - central
 302   AddHistogram(ID(kHistNsigmaTPCe), "TPC N_{#sigma,p} - e hypothesis;p (GeV/#it{c});N_{#sigma,p}",
 303                100, 0., 10.,
 304                100, -25., 25.);
 305   AddHistogram(ID(kHistNsigmaTPCmu), "TPC N_{#sigma,p} - #mu hypothesis;p (GeV/#it{c});N_{#sigma,p}",
 306                100, 0., 10.,
 307                100, -25., 25.);
 308   AddHistogram(ID(kHistNsigmaTPCpi), "TPC N_{#sigma,p} - #pi hypothesis;p (GeV/#it{c});N_{#sigma,p}",
 309                100, 0., 10.,
 310                100, -25., 25.);
 311   AddHistogram(ID(kHistNsigmaTPCk), "TPC N_{#sigma,p} - K hypothesis;p (GeV/#it{c});N_{#sigma,p}",
 312                100, 0., 10.,
 313                100, -25., 25.);
 314   AddHistogram(ID(kHistNsigmaTPCp), "TPC N_{#sigma,p} - p hypothesis;p (GeV/#it{c});N_{#sigma,p}",
 315                100, 0., 10.,
 316                100, -25., 25.);
 317   AddHistogram(ID(kHistNsigmaTPCd), "TPC N_{#sigma,p} - d hypothesis;p (GeV/#it{c});N_{#sigma,p}",
 318                100, 0., 10.,
 319                100, -25., 25.);
 320   AddHistogram(ID(kHistNsigmaTPCe_e), "TPC N_{#sigma,p} - e hypothesis (id. e);p (GeV/#it{c});N_{#sigma,p}",
 321                100, 0., 10.,
 322                100, -25., 25.);
 323
 324   AddHistogram(ID(kHistNsigmaTOFe), "TOF N_{#sigma,p} - e hypothesis;p (GeV/#it{c});N_{#sigma,p}",
 325                100, 0., 10.,
 326                200, -100., 100.);
 327   AddHistogram(ID(kHistNsigmaTOFmu), "TOF N_{#sigma,p} - #mu hypothesis;p (GeV/#it{c});N_{#sigma,p}",
 328                100, 0., 10.,
 329                200, -100., 100.);
 330   AddHistogram(ID(kHistNsigmaTOFpi), "TOF N_{#sigma,p} - #pi hypothesis;p (GeV/#it{c});N_{#sigma,p}",
 331                100, 0., 10.,
 332                200, -100., 100.);
 333   AddHistogram(ID(kHistNsigmaTOFk), "TOF N_{#sigma,p} - K hypothesis;p (GeV/#it{c});N_{#sigma,p}",
 334                100, 0., 10.,
 335                200, -100., 100.);
 336   AddHistogram(ID(kHistNsigmaTOFp), "TOF N_{#sigma,p} - p hypothesis;p (GeV/#it{c});N_{#sigma,p}",
 337                100, 0., 10.,
 338                200, -100., 100.);
 339   AddHistogram(ID(kHistNsigmaTOFd), "TOF N_{#sigma,p} - d hypothesis;p (GeV/#it{c});N_{#sigma,p}",
 340                100, 0., 10.,
 341                200, -100., 100.);
 342   AddHistogram(ID(kHistNsigmaTOFmismatch), "TOF N_{#sigma,p} - mismatch;p (GeV/#it{c});N_{#sigma,p}",
 343                100, 0., 10.,
 344                200, -100., 100.);
 345
 346   // Nsigma templates - semi-central
 347   AddHistogram(ID(kHistNsigmaTPCeSemi), "TPC N_{#sigma,p} - e hypothesis;p (GeV/#it{c});N_{#sigma,p}",
 348                100, 0., 10.,
 349                100, -25., 25.);
 350   AddHistogram(ID(kHistNsigmaTPCmuSemi), "TPC N_{#sigma,p} - #mu hypothesis;p (GeV/#it{c});N_{#sigma,p}",
 351                100, 0., 10.,
 352                100, -25., 25.);
 353   AddHistogram(ID(kHistNsigmaTPCpiSemi), "TPC N_{#sigma,p} - #pi hypothesis;p (GeV/#it{c});N_{#sigma,p}",
 354                100, 0., 10.,
 355                100, -25., 25.);
 356   AddHistogram(ID(kHistNsigmaTPCkSemi), "TPC N_{#sigma,p} - K hypothesis;p (GeV/#it{c});N_{#sigma,p}",
 357                100, 0., 10.,
 358                100, -25., 25.);
 359   AddHistogram(ID(kHistNsigmaTPCpSemi), "TPC N_{#sigma,p} - p hypothesis;p (GeV/#it{c});N_{#sigma,p}",
 360                100, 0., 10.,
 361                100, -25., 25.);
 362   AddHistogram(ID(kHistNsigmaTPCdSemi), "TPC N_{#sigma,p} - d hypothesis;p (GeV/#it{c});N_{#sigma,p}",
 363                100, 0., 10.,
 364                100, -25., 25.);
 365   AddHistogram(ID(kHistNsigmaTPCe_eSemi), "TPC N_{#sigma,p} - e hypothesis (id. e);p (GeV/#it{c});N_{#sigma,p}",
 366                100, 0., 10.,
 367                100, -25., 25.);
 368
 369   AddHistogram(ID(kHistNsigmaTOFeSemi), "TOF N_{#sigma,p} - e hypothesis;p (GeV/#it{c});N_{#sigma,p}",
 370                100, 0., 10.,
 371                200, -100., 100.);
 372   AddHistogram(ID(kHistNsigmaTOFmuSemi), "TOF N_{#sigma,p} - #mu hypothesis;p (GeV/#it{c});N_{#sigma,p}",
 373                100, 0., 10.,
 374                200, -100., 100.);
 375   AddHistogram(ID(kHistNsigmaTOFpiSemi), "TOF N_{#sigma,p} - #pi hypothesis;p (GeV/#it{c});N_{#sigma,p}",
 376                100, 0., 10.,
 377                200, -100., 100.);
 378   AddHistogram(ID(kHistNsigmaTOFkSemi), "TOF N_{#sigma,p} - K hypothesis;p (GeV/#it{c});N_{#sigma,p}",
 379                100, 0., 10.,
 380                200, -100., 100.);
 381   AddHistogram(ID(kHistNsigmaTOFpSemi), "TOF N_{#sigma,p} - p hypothesis;p (GeV/#it{c});N_{#sigma,p}",
 382                100, 0., 10.,
 383                200, -100., 100.);
 384   AddHistogram(ID(kHistNsigmaTOFdSemi), "TOF N_{#sigma,p} - d hypothesis;p (GeV/#it{c});N_{#sigma,p}",
 385                100, 0., 10.,
 386                200, -100., 100.);
 387   AddHistogram(ID(kHistNsigmaTOFmismatchSemi), "TOF N_{#sigma,p} - mismatch;p (GeV/#it{c});N_{#sigma,p}",
 388                100, 0., 10.,
 389                200, -100., 100.);
 390
 391   // delta templates
 392   AddHistogram(ID(kHistDeltaTOFe), "TOF #Delta;p (GeV/#it{c});t (ns)",
 393                100, 0., 10., 200, -2., 2.);
 394   AddHistogram(ID(kHistDeltaTOFmu), "TOF #Delta;p (GeV/#it{c});t (ns)",
 395                100, 0., 10., 200, -2., 2.);
 396   AddHistogram(ID(kHistDeltaTOFpi), "TOF #Delta;p (GeV/#it{c});t (ns)",
 397                100, 0., 10., 200, -2., 2.);
 398   AddHistogram(ID(kHistDeltaTOFk), "TOF #Delta;p (GeV/#it{c});t (ns)",
 399                100, 0., 10., 200, -2., 2.);
 400   AddHistogram(ID(kHistDeltaTOFp), "TOF #Delta;p (GeV/#it{c});t (ns)",
 401                100, 0., 10., 200, -2., 2.);
 402   AddHistogram(ID(kHistDeltaTOFd), "TOF #Delta;p (GeV/#it{c});t (ns)",
 403                100, 0., 10., 200, -2., 2.);
 404
 405   AddHistogram(ID(kHistDeltaTOFeSemi), "TOF #Delta;p (GeV/#it{c});t (ns)",
 406                100, 0., 10., 200, -2., 2.);
 407   AddHistogram(ID(kHistDeltaTOFmuSemi), "TOF #Delta;p (GeV/#it{c});t (ns)",
 408                100, 0., 10., 200, -2., 2.);
 409   AddHistogram(ID(kHistDeltaTOFpiSemi), "TOF #Delta;p (GeV/#it{c});t (ns)",
 410                100, 0., 10., 200, -2., 2.);
 411   AddHistogram(ID(kHistDeltaTOFkSemi), "TOF #Delta;p (GeV/#it{c});t (ns)",
 412                100, 0., 10., 200, -2., 2.);
 413   AddHistogram(ID(kHistDeltaTOFpSemi), "TOF #Delta;p (GeV/#it{c});t (ns)",
 414                100, 0., 10., 200, -2., 2.);
 415   AddHistogram(ID(kHistDeltaTOFdSemi), "TOF #Delta;p (GeV/#it{c});t (ns)",
 416                100, 0., 10., 200, -2., 2.);
 417
 418   // sigma comparisons
 419   // AddHistogram(ID(kHistExpSigmaTOFe), "TOF time of flight;p (GeV/#it{c});t (ns)",
 420   //           100, 0., 10., 200, 0., .25);
 421   // AddHistogram(ID(kHistExpSigmaTOFmu), "TOF time of flight;p (GeV/#it{c});t (ns)",
 422   //           100, 0., 10., 200, 0., .25);
 423   // AddHistogram(ID(kHistExpSigmaTOFpi), "TOF time of flight;p (GeV/#it{c});t (ns)",
 424   //           100, 0., 10., 200, 0., .25);
 425   // AddHistogram(ID(kHistExpSigmaTOFk), "TOF time of flight;p (GeV/#it{c});t (ns)",
 426   //           100, 0., 10., 200, 0., .25);
 427   // AddHistogram(ID(kHistExpSigmaTOFp), "TOF time of flight;p (GeV/#it{c});t (ns)",
 428   //           100, 0., 10., 200, 0., .25);
 429   // AddHistogram(ID(kHistExpSigmaTOFd), "TOF time of flight;p (GeV/#it{c});t (ns)",
 430   //           100, 0., 10., 200, 0., .25);
 431
 432   // AddHistogram(ID(kHistExpSigmaTOFeSemi), "TOF time of flight;p (GeV/#it{c});t (ns)",
 433   //           100, 0., 10., 200, 0., .25);
 434   // AddHistogram(ID(kHistExpSigmaTOFmuSemi), "TOF time of flight;p (GeV/#it{c});t (ns)",
 435   //           100, 0., 10., 200, 0., .25);
 436   // AddHistogram(ID(kHistExpSigmaTOFpiSemi), "TOF time of flight;p (GeV/#it{c});t (ns)",
 437   //           100, 0., 10., 200, 0., .25);
 438   // AddHistogram(ID(kHistExpSigmaTOFkSemi), "TOF time of flight;p (GeV/#it{c});t (ns)",
 439   //           100, 0., 10., 200, 0., .25);
 440   // AddHistogram(ID(kHistExpSigmaTOFpSemi), "TOF time of flight;p (GeV/#it{c});t (ns)",
 441   //           100, 0., 10., 200, 0., .25);
 442   // AddHistogram(ID(kHistExpSigmaTOFdSemi), "TOF time of flight;p (GeV/#it{c});t (ns)",
 443   //           100, 0., 10., 200, 0., .25);
 444
 445   // AddHistogram(ID(kHistCmpSigmaTOFe), "#sigma comparison;exp #sigma;template #sigma",
 446   //           200, 0., .25, 200, 0., .25);
 447   // AddHistogram(ID(kHistCmpSigmaTOFmu), "#sigma comparison;exp #sigma;template #sigma",
 448   //           200, 0., .25, 200, 0., .25);
 449   // AddHistogram(ID(kHistCmpSigmaTOFpi), "#sigma comparison;exp #sigma;template #sigma",
 450   //           200, 0., .25, 200, 0., .25);
 451   // AddHistogram(ID(kHistCmpSigmaTOFk), "#sigma comparison;exp #sigma;template #sigma",
 452   //           200, 0., .25, 200, 0., .25);
 453   // AddHistogram(ID(kHistCmpSigmaTOFp), "#sigma comparison;exp #sigma;template #sigma",
 454   //           200, 0., .25, 200, 0., .25);
 455   // AddHistogram(ID(kHistCmpSigmaTOFd), "#sigma comparison;exp #sigma;template #sigma",
 456   //           200, 0., .25, 200, 0., .25);
 457
 458   // AddHistogram(ID(kHistCmpSigmaTOFeSemi), "#sigma comparison;exp #sigma;template #sigma",
 459   //           200, 0., .25, 200, 0., .25);
 460   // AddHistogram(ID(kHistCmpSigmaTOFmuSemi), "#sigma comparison;exp #sigma;template #sigma",
 461   //           200, 0., .25, 200, 0., .25);
 462   // AddHistogram(ID(kHistCmpSigmaTOFpiSemi), "#sigma comparison;exp #sigma;template #sigma",
 463   //           200, 0., .25, 200, 0., .25);
 464   // AddHistogram(ID(kHistCmpSigmaTOFkSemi), "#sigma comparison;exp #sigma;template #sigma",
 465   //           200, 0., .25, 200, 0., .25);
 466   // AddHistogram(ID(kHistCmpSigmaTOFpSemi), "#sigma comparison;exp #sigma;template #sigma",
 467   //           200, 0., .25, 200, 0., .25);
 468   // AddHistogram(ID(kHistCmpSigmaTOFdSemi), "#sigma comparison;exp #sigma;template #sigma",
 469   //           200, 0., .25, 200, 0., .25);
 470
 471   // Nsigma distributions
 472   AddHistogram(ID(kHistNsigmaTPCTOF), "N_{#sigma,p} TPC-TOF;p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
 473                100, 0., 10.,
 474                100, -25., 25.,
 475                200, -100., 100.);
 476   AddHistogram(ID(kHistNsigmaTPCTOFPt), "N_{#sigma,p} TPC-TOF;p_{T} (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
 477                100, 0., 10.,
 478                100, -25., 25.,
 479                200, -100., 100.);
 480   AddHistogram(ID(kHistNsigmaTPCTOFUsed), "N_{#sigma,p} TPC-TOF;p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
 481                100, 0., 10.,
 482                100, -25., 25.,
 483                200, -100., 100.);
 484   AddHistogram(ID(kHistNsigmaTPCTOFUsedCentral), "N_{#sigma,p} TPC-TOF (central);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
 485                100, 0., 10.,
 486                100, -25., 25.,
 487                200, -100., 100.);
 488   AddHistogram(ID(kHistNsigmaTPCTOFUsedSemiCentral), "N_{#sigma,p} TPC-TOF (semi-central);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
 489                100, 0., 10.,
 490                100, -25., 25.,
 491                200, -100., 100.);
 492   AddHistogram(ID(kHistNsigmaTPCTOFUsedPt), "N_{#sigma,p} TPC-TOF;p_{T} (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
 493                50, 0., 10.,
 494                100, -25., 25.,
 495                200, -100., 100.);
 496   AddHistogram(ID(kHistNsigmaTPCTOFUsedPtCentral), "N_{#sigma,p} TPC-TOF;p_{T} (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
 497                50, 0., 10.,
 498                100, -25., 25.,
 499                200, -100., 100.);
 500   AddHistogram(ID(kHistNsigmaTPCTOFUsedPtSemiCentral), "N_{#sigma,p} TPC-TOF;p_{T} (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
 501                50, 0., 10.,
 502                100, -25., 25.,
 503                200, -100., 100.);
 504
 505   AddHistogram(ID(kHistNsigmaTPCTOFUsedCentralMCe), "N_{#sigma,p} TPC-TOF (central, MC e);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
 506                100, 0., 10.,
 507                100, -25., 25.,
 508                200, -100., 100.);
 509   AddHistogram(ID(kHistNsigmaTPCTOFUsedCentralMCmu), "N_{#sigma,p} TPC-TOF (central, MC #mu);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
 510                100, 0., 10.,
 511                100, -25., 25.,
 512                200, -100., 100.);
 513   AddHistogram(ID(kHistNsigmaTPCTOFUsedCentralMCpi), "N_{#sigma,p} TPC-TOF (central, MC pi);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
 514                100, 0., 10.,
 515                100, -25., 25.,
 516                200, -100., 100.);
 517   AddHistogram(ID(kHistNsigmaTPCTOFUsedCentralMCk), "N_{#sigma,p} TPC-TOF (central, MC K);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
 518                100, 0., 10.,
 519                100, -25., 25.,
 520                200, -100., 100.);
 521   AddHistogram(ID(kHistNsigmaTPCTOFUsedCentralMCp), "N_{#sigma,p} TPC-TOF (central, MC p);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
 522                100, 0., 10.,
 523                100, -25., 25.,
 524                200, -100., 100.);
 525   AddHistogram(ID(kHistNsigmaTPCTOFUsedCentralMCd), "N_{#sigma,p} TPC-TOF (central, MC d);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
 526                100, 0., 10.,
 527                100, -25., 25.,
 528                200, -100., 100.);
 529
 530   AddHistogram(ID(kHistNsigmaTPCTOFUsedSemiCentralMCe), "N_{#sigma,p} TPC-TOF (semi-central, MC e);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
 531                100, 0., 10.,
 532                100, -25., 25.,
 533                200, -100., 100.);
 534   AddHistogram(ID(kHistNsigmaTPCTOFUsedSemiCentralMCmu), "N_{#sigma,p} TPC-TOF (semi-central, MC #mu);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
 535                100, 0., 10.,
 536                100, -25., 25.,
 537                200, -100., 100.);
 538   AddHistogram(ID(kHistNsigmaTPCTOFUsedSemiCentralMCpi), "N_{#sigma,p} TPC-TOF (semi-central, MC pi);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
 539                100, 0., 10.,
 540                100, -25., 25.,
 541                200, -100., 100.);
 542   AddHistogram(ID(kHistNsigmaTPCTOFUsedSemiCentralMCk), "N_{#sigma,p} TPC-TOF (semi-central, MC K);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
 543                100, 0., 10.,
 544                100, -25., 25.,
 545                200, -100., 100.);
 546   AddHistogram(ID(kHistNsigmaTPCTOFUsedSemiCentralMCp), "N_{#sigma,p} TPC-TOF (semi-central, MC p);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
 547                100, 0., 10.,
 548                100, -25., 25.,
 549                200, -100., 100.);
 550   AddHistogram(ID(kHistNsigmaTPCTOFUsedSemiCentralMCd), "N_{#sigma,p} TPC-TOF (semi-central, MC d);p (GeV/#it{c});N_{#sigma,p}^{TPC};N_{#sigma,p}^{TOF}",
 551                100, 0., 10.,
 552                100, -25., 25.,
 553                200, -100., 100.);
 554
 555   AddHistogram(ID(kHistNevMix), "no. of events for mixing;N_{ev};counts",
 556                20, -.5, 19.5,
 557                2, -.5, 1.5);
 558
 559   AddHistogram(ID(kHistEvPlane), "default event plane;#Psi;counts",
 560                100, -0. * TMath::Pi(), 1. * TMath::Pi());
 561   AddHistogram(ID(kHistEvPlaneRes), "resolution of default event plane;R^{2};centrality (%)",
 562                100, -1., 1.,
 563                50, 0., 100.);
 564   AddHistogram(ID(kHistEvPlaneUsed), "default event plane;#Psi;counts",
 565                100, -0. * TMath::Pi(), 1. * TMath::Pi(),
 566                kClLast, -.5, kClLast-.5);
 567   AddHistogram(ID(kHistEvPlaneCheck), "backup event plane;#Psi;counts",
 568                100, -1. * TMath::Pi(), 1. * TMath::Pi());
 569   AddHistogram(ID(kHistEvPlaneCheckUsed), "backup event plane;#Psi;counts",
 570                100, -0. * TMath::Pi(), 1. * TMath::Pi(),
 571                kClLast, -.5, kClLast-.5);
 572   AddHistogram(ID(kHistEvPlane3), "3rd order event plane;#Psi;counts",
 573                100, -0. * TMath::Pi(), 2./3. * TMath::Pi());
 574   AddHistogram(ID(kHistEvPlaneCorr), "default - backup event plane;#Psi_{def};#Psi_{bak};event class",
 575                100, -0. * TMath::Pi(), 1. * TMath::Pi(),
 576                100, -0. * TMath::Pi(), 1. * TMath::Pi(),
 577                kClLast, -.5, kClLast-.5);
 578   AddHistogram(ID(kHistEvPlaneCross), "default - backup event plane;R_{a} R_{b};centrality",
 579                100, -1., 1.,
 580                50, 0., 100.);
 581   AddHistogram(ID(kHistEvPlaneCorrNoTrgJets), "event plane w/ and w/o trigger jets;#Psi_{full};#Psi_{notrgjets};event class",
 582                100, -0. * TMath::Pi(), 1. * TMath::Pi(),
 583                100, -0. * TMath::Pi(), 1. * TMath::Pi(),
 584                kClLast, -.5, kClLast-.5);
 585   AddHistogram(ID(kHistEvPlaneCorrNoTrgJetsTrgd), "event plane w/ and w/o trigger jets;#Psi_{full};#Psi_{notrgjets};event class",
 586                100, -0. * TMath::Pi(), 1. * TMath::Pi(),
 587                100, -0. * TMath::Pi(), 1. * TMath::Pi(),
 588                kClLast, -.5, kClLast-.5);
 589
 590   AddHistogram(ID(kHistJetPtCentral), "jet spectrum - central;p_{T}^{jet,ch} (GeV/#it{c});counts",
 591                40, 0., 200.);
 592   AddHistogram(ID(kHistJetPtSemi), "jet spectrum - semi-peripheral;p_{T}^{jet,ch} (GeV/#it{c});counts",
 593                40, 0., 200.);
 594
 595   AddHistogram(ID(kHistEtaPhiTrgHad), "trg had;#varphi;#eta",
 596                100, -0. * TMath::Pi(), 2. * TMath::Pi(),
 597                100, -2., 2.);
 598   AddHistogram(ID(kHistEtaPhiTrgJet), "trg jet;#varphi;#eta",
 599                100, -0. * TMath::Pi(), 2. * TMath::Pi(),
 600                100, -2., 2.);
 601   AddHistogram(ID(kHistEtaPhiAssHad), "ass had;#varphi;#eta",
 602                100, -0. * TMath::Pi(), 2. * TMath::Pi(),
 603                100, -2., 2.);
 604   AddHistogram(ID(kHistEtaPhiAssProt), "ass proton;#varphi;#eta",
 605                100, -0. * TMath::Pi(), 2. * TMath::Pi(),
 606                100, -2., 2.);
 607
 608   AddHistogram(ID(kHistPhiTrgJetEvPlane), "trg jet;#varphi - #Psi_{ev};centrality",
 609                100, -0. * TMath::Pi(), 2. * TMath::Pi(),
 610                100, 0., 100.);
 611   AddHistogram(ID(kHistPhiTrgHadEvPlane), "trg had;#varphi - #Psi_{ev};centrality",
 612                100, -0. * TMath::Pi(), 2. * TMath::Pi(),
 613                100, 0., 100.);
 614   AddHistogram(ID(kHistPhiRndTrgJetEvPlane), "rnd trg jet;#varphi - #Psi_{ev};centrality",
 615                100, -0. * TMath::Pi(), 2. * TMath::Pi(),
 616                100, 0., 100.);
 617   AddHistogram(ID(kHistPhiRndTrgHadEvPlane), "rnd trg had;#varphi - #Psi_{ev};centrality",
 618                100, -0. * TMath::Pi(), 2. * TMath::Pi(),
 619                100, 0., 100.);
 620   AddHistogram(ID(kHistPhiAssHadEvPlane), "ass had;#varphi - #Psi_{ev};centrality",
 621                100, -0. * TMath::Pi(), 2. * TMath::Pi(),
 622                100, 0., 100.);
 623   AddHistogram(ID(kHistPhiAssProtEvPlane), "ass prot;#varphi - #Psi_{ev};centrality",
 624                100, -0. * TMath::Pi(), 2. * TMath::Pi(),
 625                100, 0., 100.);
 626   AddHistogram(ID(kHistPhiAssHadVsEvPlane), "ass had;#Psi_{ev};#varphi;centrality",
 627                100, -0. * TMath::Pi(), 1. * TMath::Pi(),
 628                100, -0. * TMath::Pi(), 2. * TMath::Pi(),
 629                100, 0., 100.);
 630
 631   AddHistogram(ID(kHistPhiTrgJetEvPlane3), "trg jet;#varphi - #Psi_{ev};centrality",
 632                100, -0. * TMath::Pi(), 2. * TMath::Pi(),
 633                100, 0., 100.);
 634   AddHistogram(ID(kHistPhiTrgHadEvPlane3), "trg had;#varphi - #Psi_{ev};centrality",
 635                100, -0. * TMath::Pi(), 2. * TMath::Pi(),
 636                100, 0., 100.);
 637   AddHistogram(ID(kHistPhiAssHadEvPlane3), "ass had;#varphi - #Psi_{ev};centrality",
 638                100, -0. * TMath::Pi(), 2. * TMath::Pi(),
 639                100, 0., 100.);
 640   AddHistogram(ID(kHistPhiAssProtEvPlane3), "ass prot;#varphi - #Psi_{ev};centrality",
 641                100, -0. * TMath::Pi(), 2. * TMath::Pi(),
 642                100, 0., 100.);
 643
 644   for (Int_t iCorr = 0; iCorr < kCorrLast; ++iCorr) {
 645     for (Int_t iCl = 0; iCl < kClLast; ++iCl) {
 646       for (Int_t iEv = 0; iEv < kEvLast; ++iEv) {
 647         // we don't need the mixed event histograms for the embedded excess particles
 648         if ((iCorr > kCorrRndJetProt) && (iEv == kEvMix))
 649           continue;
 650
 651         GetHistCorr((CorrType_t) iCorr, (Class_t) iCl, (Ev_t) iEv) =
 652           new AliHistCorr(Form("corr_%s_%s_%s", fkCorrTypeName[iCorr], fkClassName[iCl], fkEvName[iEv]), fOutputList);
 653       }
 654     }
 655   }
 656
 657   PostData(1, fOutputList);
 658 }
 659
 660 Bool_t AliAnalysisTaskJetProtonCorr::Notify()
 661 {
 662   // actions to be taken upon notification about input file change
 663
 664   return AliAnalysisTaskSE::Notify();
 665 }
 666
 667 void AliAnalysisTaskJetProtonCorr::SetParamsTOF()
 668 {
 669   fParamsTOF =
 670     dynamic_cast<AliTOFPIDParams*> (fOADBContainerTOF->GetObject(fRunNumber, "TOFparams"));
 671
 672   if (!fParamsTOF)
 673     AliError(Form("failed to load TOF parameters for run %i", fRunNumber));
 674   else if (fDebug > 2) {
 675     printf("loaded TOF parameters for run %i\n",
 676            fRunNumber);
 677     printf("   intrinsic resolution: %f ps\n",
 678            fParamsTOF->GetTOFresolution());
 679     printf("   tail fraction: %f\n",
 680            fParamsTOF->GetTOFtail());
 681     printf("   start time method: %i\n",
 682            fParamsTOF->GetStartTimeMethod());
 683     printf("   TOF signal parametrization: %f, %f, %f, %f\n",
 684            fParamsTOF->GetSigParams(0), fParamsTOF->GetSigParams(1),
 685            fParamsTOF->GetSigParams(2), fParamsTOF->GetSigParams(3));
 686     printf("   matching loss MC: %f%%\n",
 687            fParamsTOF->GetTOFmatchingLossMC());
 688     printf("   additional mismatch for MC: %f%%\n",
 689            fParamsTOF->GetTOFadditionalMismForMC());
 690     printf("   time offset: %f\n",
 691            fParamsTOF->GetTOFtimeOffset());
 692   }
 693 }
 694
 695 void AliAnalysisTaskJetProtonCorr::UserExec(Option_t * /* option */)
 696 {
 697   // actual work
 698
 699   // setup pointers to input data (null if unavailable)
 700   // mcEvent:  MC input
 701   // esdEvent: ESD input
 702   // outEvent: AOD output
 703   // aodEvent: AOD input if available, otherwise AOD output
 704
 705   fMCEvent   = this->MCEvent();
 706   fESDEvent  = dynamic_cast<AliESDEvent*>(this->InputEvent()); // could also be AOD input
 707   AliAODEvent* outEvent  = this->AODEvent();
 708   fAODEvent  = dynamic_cast<AliAODEvent*> (this->InputEvent());
 709   if (!fAODEvent)
 710     fAODEvent = outEvent;
 711
 712   if ((fDebug > 0) && fESDEvent)
 713     printf("event: %s-%06i\n", CurrentFileName(), fESDEvent->GetEventNumberInFile());
 714
 715   // record number of sampled events and detect trigger contributions
 716   FillH1(kHistStat, kStatSeen);
 717   if (!DetectTriggers()) {
 718     AliError("Failed to detect the triggers");
 719     return;
 720   }
 721
 722   if (!IsTrigger(kTriggerInt))
 723     return;
 724
 725   FillH1(kHistStat, kStatTrg);
 726
 727   // prepare the event
 728   // (make sure it is cleaned up in the end)
 729   if (PrepareEvent()) {
 730     FillH1(kHistStat, kStatUsed);
 731     FillH1(kHistCentrality, fCentrality);
 732     FillH1(kHistCentralityCheck, fCentralityCheck);
 733
 734     FillH1(kHistVertexZ, fZvtx);
 735
 736     // event cuts
 737     if (TMath::Abs(fZvtx) > 10.)
 738       goto stop;
 739     if (GetCentrality() > 90.)
 740       goto stop;
 741
 742     FillH1(kHistStat, kStatEvCuts);
 743
 744     // event category
 745     DetectClasses();
 746     if (IsClass(kClCentral))
 747       FillH1(kHistStat, kStatCentral);
 748     if (IsClass(kClSemiCentral))
 749       FillH1(kHistStat, kStatSemiCentral);
 750
 751     FillH1(kHistEvPlane, fEventPlaneAngle);
 752     FillH2(kHistEvPlaneRes, fEventPlaneRes, fCentrality);
 753     FillH1(kHistEvPlaneCheck, fEventPlaneAngleCheck);
 754     FillH1(kHistEvPlane3, fEventPlaneAngle3);
 755     FillH2(kHistEvPlaneCross, TMath::Cos(2. * (fEventPlaneAngle - fEventPlaneAngleCheck)), GetCentrality());
 756     for (Int_t iClass = 0; iClass < kClLast; ++iClass) {
 757       if (IsClass((Class_t) iClass)) {
 758         FillH2(kHistCentralityUsed, fCentrality, iClass);
 759         FillH2(kHistCentralityCheckUsed, fCentralityCheck, iClass);
 760         FillH2(kHistEvPlaneUsed, fEventPlaneAngle, iClass);
 761         FillH2(kHistEvPlaneCheckUsed, fEventPlaneAngleCheck, iClass);
 762         FillH3(kHistEvPlaneCorr, fEventPlaneAngle, fEventPlaneAngleCheck, iClass);
 763       }
 764     }
 765
 766     Bool_t corrEta  = fPIDResponse->UseTPCEtaCorrection();
 767     Bool_t corrMult = fPIDResponse->UseTPCMultiplicityCorrection();
 768
 769     // select trigger particles and potential associated particles/protons
 770     TObjArray trgArray[kTrgLast];
 771     TObjArray assArray[kAssLast];
 772
 773     // prepare TOF response
 774     TF1 fTOFsignal("fTOFsignal", &AliAnalysisTaskJetProtonCorr::TOFsignal, -2440., 2440., 4);
 775     fTOFsignal.SetParameter(0, 1.);
 776     fTOFsignal.SetParameter(1, 0.);
 777     if (fParamsTOF) {
 778       Float_t res  = fParamsTOF->GetTOFresolution();
 779       Float_t tail = fParamsTOF->GetTOFtail() * res;
 780       fTOFsignal.SetParameter(2, res);
 781       fTOFsignal.SetParameter(3, tail);
 782     }
 783     else {
 784       fTOFsignal.SetParameter(2, 85.);
 785       fTOFsignal.SetParameter(3, 80.);
 786     }
 787
 788     // associate candidates
 789     const Int_t nPrimTracksAss = fPrimTrackArrayAss ? fPrimTrackArrayAss->GetEntries() : 0;
 790     FillH2(kHistCentralityVsMult, fCentrality, nPrimTracksAss);
 791     for (Int_t iTrack = 0; iTrack < nPrimTracksAss; ++iTrack) {
 792       AliVTrack *trk = (AliVTrack*) fPrimTrackArrayAss->At(iTrack);
 793       FillH2(kHistSignalTPC, trk->P(), trk->GetTPCsignal());
 794       // ??? pt or p?
 795       FillH2(kHistSignalTOF, trk->P(), trk->GetTOFsignal() * 1.e-3); // ps -> ns
 796       FillH3(kHistNsigmaTPCTOF,
 797              trk->P(),
 798              fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
 799              fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
 800       FillH3(kHistNsigmaTPCTOFPt,
 801              trk->Pt(),
 802              fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
 803              fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
 804
 805       FillH3(kHistPhiAssHadVsEvPlane, fEventPlaneAngle, trk->Phi(), fCentrality);
 806
 807       if (AcceptAssoc(trk)) {
 808         Float_t phiRel = GetPhiRel2(trk);
 809         Float_t phiRel3 = trk->Phi() - fEventPlaneAngle3;
 810         if (phiRel3 < 0.)
 811           phiRel3 += 2. * TMath::Pi();
 812         assArray[kAssHad].Add(trk);
 813         FillH1(kHistEtaPhiAssHad, trk->Phi(), trk->Eta());
 814         FillH2(kHistPhiAssHadEvPlane, phiRel, fCentrality);
 815         FillH2(kHistPhiAssHadEvPlane3, phiRel3, fCentrality);
 816         FillH3(kHistNsigmaTPCTOFUsed,
 817                trk->P(),
 818                fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
 819                fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
 820         FillH3(kHistNsigmaTPCTOFUsedPt,
 821                trk->Pt(),
 822                fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
 823                fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
 824
 825         // central events
 826         if (IsClass(kClCentral)) {
 827           FillH3(kHistNsigmaTPCTOFUsedCentral,
 828                  trk->P(),
 829                  fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
 830                  fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
 831           FillH3(kHistNsigmaTPCTOFUsedPtCentral,
 832                  trk->Pt(),
 833                  fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
 834                  fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
 835           if (fMCEvent) {
 836             if (AliMCParticle *mcPart = (AliMCParticle*) fMCEvent->GetTrack(TMath::Abs(trk->GetLabel()))) {
 837               for (Int_t iParticle = 0; iParticle <= AliPID::kDeuteron; ++iParticle) {
 838                 if (TMath::Abs(mcPart->Particle()->GetPdgCode()) == AliPID::ParticleCode(iParticle))
 839                   FillH3(kHistNsigmaTPCTOFUsedCentralMCe,
 840                          trk->P(),
 841                          fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
 842                          fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton),
 843                          1., iParticle);
 844               }
 845             }
 846           }
 847         }
 848
 849         // semi-central events
 850         if (IsClass(kClSemiCentral)) {
 851           FillH3(kHistNsigmaTPCTOFUsedSemiCentral,
 852                  trk->P(),
 853                  fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
 854                  fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
 855           FillH3(kHistNsigmaTPCTOFUsedPtSemiCentral,
 856                  trk->Pt(),
 857                  fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
 858                  fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
 859           if (fMCEvent) {
 860             if (AliMCParticle *mcPart = (AliMCParticle*) fMCEvent->GetTrack(TMath::Abs(trk->GetLabel()))) {
 861               for (Int_t iParticle = 0; iParticle <= AliPID::kDeuteron; ++iParticle) {
 862                 if (TMath::Abs(mcPart->Particle()->GetPdgCode()) == AliPID::ParticleCode(iParticle))
 863                   FillH3(kHistNsigmaTPCTOFUsedSemiCentralMCe,
 864                          trk->P(),
 865                          fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
 866                          fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton),
 867                          1., iParticle);
 868               }
 869             }
 870           }
 871         }
 872
 873         // protons
 874         if (IsProton(trk)) {
 875           assArray[kAssProt].Add(trk);
 876           FillH1(kHistEtaPhiAssProt, trk->Phi(), trk->Eta());
 877           FillH2(kHistPhiAssProtEvPlane, phiRel, fCentrality);
 878           FillH2(kHistPhiAssProtEvPlane3, phiRel3, fCentrality);
 879         }
 880
 881         // template generation
 882         if (fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC, trk) == AliPIDResponse::kDetPidOk) {
 883           Double_t deltaTPC;
 884           fPIDResponse->GetSignalDelta(AliPIDResponse::kTPC, trk, AliPID::kProton, deltaTPC);
 885           if (IsClass(kClCentral))
 886             FillH2(kHistDeltaTPC, trk->P(), deltaTPC);
 887           if (IsClass(kClSemiCentral))
 888             FillH2(kHistDeltaTPCSemi, trk->P(), deltaTPC);
 889
 890           Double_t expTPC = fPIDResponse->GetTPCResponse().GetExpectedSignal(trk, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, corrEta, corrMult);
 891           Double_t expSigmaTPC = fPIDResponse->GetTPCResponse().GetExpectedSigma(trk, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, corrEta, corrMult);
 892
 893           // loop over particles
 894           for (Int_t iParticle = 0; iParticle <= AliPID::kDeuteron; ++iParticle) {
 895             Double_t expTPCx = fPIDResponse->GetTPCResponse().GetExpectedSignal(trk, AliPID::EParticleType(iParticle), AliTPCPIDResponse::kdEdxDefault, corrEta, corrMult);
 896             Double_t expSigmaTPCx = fPIDResponse->GetTPCResponse().GetExpectedSigma(trk, AliPID::EParticleType(iParticle), AliTPCPIDResponse::kdEdxDefault, corrEta, corrMult);
 897             Double_t rndTPCx = gRandom->Gaus(expTPCx, expSigmaTPCx);
 898
 899             if(IsClass(kClCentral))
 900               FillH2(kHistNsigmaTPCe, trk->P(), (rndTPCx-expTPC)/expSigmaTPC, 1., iParticle);
 901             if(IsClass(kClSemiCentral))
 902               FillH2(kHistNsigmaTPCeSemi, trk->P(), (rndTPCx-expTPC)/expSigmaTPC, 1., iParticle);
 903           }
 904         }
 905
 906         if (fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF, trk) == AliPIDResponse::kDetPidOk) {
 907           Double_t deltaTOF;
 908           fPIDResponse->GetSignalDelta(AliPIDResponse::kTOF, trk, AliPID::kProton, deltaTOF);
 909
 910           AliTOFPIDResponse &tofResponse = fPIDResponse->GetTOFResponse();
 911
 912           Float_t p = trk->P();
 913           if (const AliExternalTrackParam *param = trk->GetInnerParam())
 914             p = param->GetP();
 915
 916           Double_t expTOF = fPIDResponse->GetTOFResponse().GetExpectedSignal(trk, AliPID::kProton);
 917           Double_t expSigmaTOF = fPIDResponse->GetTOFResponse().GetExpectedSigma(p, expTOF, AliPID::kProton);
 918           Double_t length = trk->GetIntegratedLength() * 1.e-2; // cm -> m
 919           Double_t tof = trk->GetTOFsignal() * 1.e-12; // ps -> s
 920           Double_t beta = length / tof / TMath::C();
 921
 922           FillH2(kHistBetaTOF, p, beta);
 923
 924           // intrinsic TOF smearing
 925           // Double_t signalSigma = TMath::Sqrt(fTOFsignal.Variance(-2000., 2000.));
 926           Double_t signalSmear = fTOFsignal.GetRandom();
 927           // t0 smearing
 928           Float_t  timezeroSigma = tofResponse.GetT0binRes(tofResponse.GetMomBin(p));
 929           Double_t timezeroSmear  = gRandom->Gaus(0., timezeroSigma);
 930           // tracking smearing (default parameters, should be overwritten from OADB)
 931           Double_t fPar[] = { 0.008, 0.008, 0.002, 40.0 };
 932           if (fParamsTOF)
 933             for (Int_t i = 0; i < 4; ++i)
 934               fPar[i] = fParamsTOF->GetSigParams(i);
 935
 936           // loop over particles
 937           for (Int_t iParticle = 0; iParticle <= AliPID::kDeuteron; ++iParticle) {
 938             Double_t expTOFx = fPIDResponse->GetTOFResponse().GetExpectedSignal(trk, AliPID::EParticleType(iParticle));
 939             // Double_t cmpSigmaTOFx = fPIDResponse->GetTOFResponse().GetExpectedSigma(p, expTOFx, AliPID::EParticleType(iParticle));
 940
 941             // tracking smearing
 942             Double_t massx = AliPID::ParticleMassZ(AliPID::EParticleType(iParticle));
 943             Double_t dppx = fPar[0] + fPar[1] * p + fPar[2] * massx / p;
 944             Double_t expSigmaTOFx = dppx * expTOFx / (1.+ p * p / (massx * massx));
 945             Double_t texpSmearx = gRandom->Gaus(0., TMath::Sqrt(expSigmaTOFx * expSigmaTOFx + fPar[3]*fPar[3]/p/p));
 946             // Double_t tmpSigmaTOFx = TMath::Sqrt(signalSigma*signalSigma +
 947             //                                  timezeroSigma*timezeroSigma +
 948             //                                  expSigmaTOFx*expSigmaTOFx + fPar[3]*fPar[3]/p/p);
 949             // printf("sigma comparison %i, %f: %f, %f, %f -> %f vs %f\n",
 950             //     iParticle, expTOFx,
 951             //     signalSigma, // signal
 952             //     timezeroSigma, // timezero
 953             //     expSigmaTOFx, // tracking
 954             //     tmpSigmaTOFx, // total
 955             //     cmpSigmaTOFx); // from PID response
 956
 957             // TOF signal
 958             Double_t rndTOFx = expTOFx + signalSmear + timezeroSmear + texpSmearx;
 959
 960             if (IsClass(kClCentral)) {
 961               FillH2(kHistNsigmaTOFe, trk->P(), (rndTOFx-expTOF)/expSigmaTOF, 1., iParticle);
 962               FillH2(kHistDeltaTOFe, trk->P(), (rndTOFx-expTOF) * 1.e-3, 1., iParticle);
 963               // FillH2(kHistExpSigmaTOFe, p, cmpSigmaTOFx * 1.e-3, 1., iParticle);
 964               // FillH2(kHistCmpSigmaTOFe, cmpSigmaTOFx * 1.e-3, tmpSigmaTOFx * 1.e-3, 1., iParticle);
 965             }
 966             if (IsClass(kClSemiCentral)) {
 967               FillH2(kHistNsigmaTOFeSemi, trk->P(), (rndTOFx-expTOF)/expSigmaTOF, 1., iParticle);
 968               FillH2(kHistDeltaTOFeSemi, trk->P(), (rndTOFx-expTOF) * 1.e-3, 1., iParticle);
 969               // FillH2(kHistCmpSigmaTOFeSemi, cmpSigmaTOFx * 1.e-3, tmpSigmaTOFx * 1.e-3, 1., iParticle);
 970               // FillH2(kHistExpSigmaTOFeSemi, p, cmpSigmaTOFx * 1.e-3, 1., iParticle);
 971             }
 972           }
 973
 974           Double_t rndTOFmismatch = AliTOFPIDResponse::GetMismatchRandomValue(trk->Eta());
 975
 976           if(IsClass(kClCentral)) {
 977             FillH2(kHistNsigmaTOFmismatch, trk->P(), (rndTOFmismatch - expTOF) / expSigmaTOF);
 978             FillH2(kHistDeltaTOF, trk->P(), deltaTOF * 1.e-3); // ps -> ns
 979           }
 980           if(IsClass(kClSemiCentral)) {
 981             FillH2(kHistNsigmaTOFmismatchSemi, trk->P(), (rndTOFmismatch - expTOF) / expSigmaTOF);
 982             FillH2(kHistDeltaTOFSemi, trk->P(), deltaTOF * 1.e-3); // ps -> ns
 983           }
 984         }
 985       }
 986     }
 987
 988     const Int_t nPrimTracksTrg = fPrimTrackArrayTrg ? fPrimTrackArrayTrg->GetEntries() : 0;
 989     for (Int_t iTrack = 0; iTrack < nPrimTracksTrg; ++iTrack) {
 990       AliVTrack *trk = (AliVTrack*) fPrimTrackArrayTrg->At(iTrack);
 991       if (AcceptTrigger(trk)) {
 992         trgArray[kTrgHad].Add(trk);
 993         FillH1(kHistEtaPhiTrgHad, trk->Phi(), trk->Eta());
 994       }
 995     }
 996
 997     // select trigger jet
 998     // and remove them from the Q vector
 999     Int_t nJets = fJetArray ? fJetArray->GetEntries() : 0;
1000     const TVector2 *qVectorOrig = fEventplane->GetQVector();
1001     TVector2 qVector;
1002     if (qVectorOrig)
1003       qVector = *qVectorOrig;
1004
1005     for (Int_t iJet = 0; iJet < nJets; ++iJet) {
1006       AliAODJet *jet = (AliAODJet*) fJetArray->At(iJet);
1007
1008       if (AcceptTrigger(jet)) {
1009         trgArray[kTrgJet].Add(jet);
1010         FillH1(kHistEtaPhiTrgJet, jet->Phi(), jet->Eta());
1011
1012         if (qVectorOrig) {
1013           Int_t nRefTracks = jet->GetRefTracks()->GetEntriesFast();
1014           for (Int_t iTrack = 0; iTrack < nRefTracks; ++iTrack) {
1015             AliVTrack *track = (AliVTrack*) jet->GetRefTracks()->At(iTrack);
1016
1017             if (fEventplane && track &&
1018                 (track->GetID() > -1)) {
1019               TVector2 evplaneContrib(fEventplane->GetQContributionX(track),
1020                                       fEventplane->GetQContributionY(track));
1021               qVector -= evplaneContrib;
1022             }
1023           }
1024         }
1025       }
1026     }
1027     if (qVectorOrig) {
1028       for (Int_t iClass = 0; iClass < kClLast; ++iClass) {
1029         if (IsClass((Class_t) iClass)) {
1030           FillH3(kHistEvPlaneCorrNoTrgJets, qVectorOrig->Phi()/2., qVector.Phi()/2., iClass);
1031           if (trgArray[kTrgJet].GetEntriesFast() > 0)
1032             FillH3(kHistEvPlaneCorrNoTrgJetsTrgd, qVectorOrig->Phi()/2., qVector.Phi()/2., iClass);
1033         }
1034       }
1035     }
1036
1037     // invent a trigger jet/hadron and correlated associates
1038     Float_t pFraction = assArray[kAssHad].GetEntries() > 0 ?
1039       (Float_t) (assArray[kAssProt].GetEntries()) / assArray[kAssHad].GetEntries() :
1040       .5;
1041     GenerateRandom(&trgArray[kTrgJetRnd], &trgArray[kTrgHadRnd],
1042                    &assArray[kAssHadJetExc], &assArray[kAssProtJetExc],
1043                    &assArray[kAssHadHadExc], &assArray[kAssProtHadExc],
1044                    pFraction);
1045
1046     // correlate, both same and mixed event
1047     for (Int_t iClass = 0; iClass < kClLast; ++iClass) {
1048       if (IsClass((Class_t) iClass)) {
1049         for (Int_t iTrg = 0; iTrg < kTrgLast; ++iTrg) {
1050           for (Int_t iAss = 0; iAss <= kAssProt; ++iAss) {
1051             // same event
1052             Correlate((Trg_t) iTrg, (Ass_t) iAss, (Class_t) iClass, kEvSame, &trgArray[iTrg], &assArray[iAss]);
1053
1054             // mixed event
1055             AliEventPool *pool = GetPool((Ass_t) iAss);
1056             if (pool)
1057               FillH2(kHistNevMix, pool->GetCurrentNEvents(), iAss);
1058             if (pool && pool->IsReady()) {
1059               AliDebug(1, Form("----- using pool: %i %i %i -----", iClass, iTrg, iAss));
1060               const Int_t nEvents = pool->GetCurrentNEvents();
1061               for (Int_t iEvent = 0; iEvent < nEvents; ++iEvent) {
1062                 TObjArray *assTracks = pool->GetEvent(iEvent);
1063                 Correlate((Trg_t) iTrg, (Ass_t) iAss, (Class_t) iClass, kEvMix, &trgArray[iTrg], assTracks, 1./nEvents);
1064               }
1065             }
1066           }
1067         }
1068
1069         // fill event pool for mixing
1070         // >= 0: don't require a trigger in the event
1071         // >= 1: require a trigger in the event
1072         // if (trgArray[iTrg].GetEntries() >= 0) {
1073         for (Int_t iAss = 0; iAss <= kAssProt; ++iAss) {
1074           AliEventPool *pool = GetPool((Ass_t) iAss);
1075           if (pool) {
1076             pool->UpdatePool(CloneTracks(&assArray[iAss]));
1077             AliDebug(1, Form("----- updating pool: %i -----", iAss));
1078             if (fDebug > 0)
1079               pool->PrintInfo();
1080           }
1081         }
1082         // }
1083
1084         // correlate artificial triggers and associates
1085         Correlate(kCorrRndJetExcHad,  (Class_t) iClass, kEvSame, &trgArray[kTrgJetRnd], &assArray[kAssHadJetExc]);
1086         Correlate(kCorrRndJetExcProt, (Class_t) iClass, kEvSame, &trgArray[kTrgJetRnd], &assArray[kAssProtJetExc]);
1087         Correlate(kCorrRndHadExcHad,  (Class_t) iClass, kEvSame, &trgArray[kTrgHadRnd], &assArray[kAssHadHadExc]);
1088         Correlate(kCorrRndHadExcProt, (Class_t) iClass, kEvSame, &trgArray[kTrgHadRnd], &assArray[kAssProtHadExc]);
1089       }
1090     }
1091
1092     trgArray[kTrgJetRnd].Delete();
1093     trgArray[kTrgHadRnd].Clear();
1094     assArray[kAssHadJetExc].Delete();
1095     assArray[kAssProtJetExc].Clear();
1096     assArray[kAssHadHadExc].Delete();
1097     assArray[kAssProtHadExc].Clear();
1098   }
1099
1100  stop:
1101   CleanUpEvent();
1102
1103   PostData(1, fOutputList);
1104 }
1105
1106 void AliAnalysisTaskJetProtonCorr::Terminate(const Option_t * /* option */)
1107 {
1108   // actions at task termination
1109
1110 }
1111
1112 void AliAnalysisTaskJetProtonCorr::PrintTask(Option_t *option, Int_t indent) const
1113 {
1114   AliAnalysisTaskSE::PrintTask(option, indent);
1115
1116   std::cout << std::setw(indent) << " " << "using jet branch: " << fJetBranchName << std::endl;
1117 }
1118
1119 Bool_t AliAnalysisTaskJetProtonCorr::DetectTriggers()
1120 {
1121   fTriggerMask = 0;
1122
1123   AliVEvent::EOfflineTriggerTypes physSel =
1124     (AliVEvent::EOfflineTriggerTypes) ((AliInputEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
1125   // TString trgClasses = InputEvent()->GetFiredTriggerClasses();
1126
1127   // physics selection
1128   if (physSel & (AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral))
1129     MarkTrigger(kTriggerInt);
1130
1131   return kTRUE;
1132 }
1133
1134 Bool_t AliAnalysisTaskJetProtonCorr::DetectClasses()
1135 {
1136   fClassMask = 0;
1137
1138   if (IsCentral())
1139     MarkClass(kClCentral);
1140
1141   if (IsSemiCentral())
1142     MarkClass(kClSemiCentral);
1143
1144   return kTRUE;
1145 }
1146
1147 Bool_t AliAnalysisTaskJetProtonCorr::PrepareEvent()
1148 {
1149   Bool_t eventGood = kTRUE;
1150
1151   // check for run change
1152   if (fRunNumber != InputEvent()->GetRunNumber()) {
1153     fRunNumber = InputEvent()->GetRunNumber();
1154     SetParamsTOF();
1155   }
1156
1157   // retrieve z-vertex position
1158   fZvtx = 100.;
1159   const AliVVertex *vtx = InputEvent()->GetPrimaryVertex();
1160   if (vtx) {
1161     FillH1(kHistStat, kStatVtx);
1162     FillH1(kHistVertexNctb, vtx->GetNContributors());
1163     if (vtx->GetNContributors() >= 3.)
1164       fZvtx = vtx->GetZ();
1165   }
1166
1167   // retrieve centrality
1168   AliCentrality *eventCentrality = InputEvent()->GetCentrality();
1169   if (eventCentrality) {
1170     fCentrality = eventCentrality->GetCentralityPercentile("V0M");
1171     fCentralityCheck = eventCentrality->GetCentralityPercentile("TRK");
1172     if (fCentrality >= 0.) {
1173       FillH1(kHistStat, kStatCent);
1174     } else {
1175       // centrality estimation not reliable
1176       eventGood = kFALSE;
1177       fCentrality = 105.;
1178     }
1179   }
1180   else
1181     eventGood = kFALSE;
1182
1183   // retrieve event plane
1184   fEventplane = InputEvent()->GetEventplane();
1185   if (fEventplane) {
1186     fEventPlaneAngle3 = fEventplane->GetEventplane("V0", InputEvent(), 3);
1187
1188     if (fUseEvplaneV0) {
1189       fEventPlaneAngle = fEventplane->GetEventplane("V0", InputEvent());
1190       fEventPlaneAngleCheck = fEventplane->GetEventplane("Q");
1191     }
1192     else {
1193       fEventPlaneAngle = fEventplane->GetEventplane("Q");
1194       fEventPlaneRes   = TMath::Cos(2 * fEventplane->GetQsubRes());
1195       fEventPlaneAngleCheck = fEventplane->GetEventplane("V0", InputEvent());
1196       // use V0 event plane angle from flow task:
1197       // fEventPlaneAngleCheck = AliAnalysisTaskVnV0::GetPsi2V0A();
1198       // printf("V0A evplane = %f\n", fEventPlaneAngleCheck);
1199       // fEventPlaneAngleCheck = AliAnalysisTaskVnV0::GetPsi2V0C();
1200       // printf("V0C evplane = %f\n", fEventPlaneAngleCheck);
1201     }
1202
1203     // ensure angles to be in [0, ...)
1204     if (fEventPlaneAngle3 < 0)
1205       fEventPlaneAngle3 += 2.*TMath::Pi()/3.;
1206     if (fEventPlaneAngle < 0)
1207       fEventPlaneAngle += TMath::Pi();
1208     if (fEventPlaneAngleCheck < 0)
1209       fEventPlaneAngleCheck += TMath::Pi();
1210
1211     FillH1(kHistStat, kStatEvPlane);
1212   }
1213   else
1214     eventGood = kFALSE;
1215
1216   // retrieve PID
1217   fPIDResponse = ((AliInputEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();
1218   if (fPIDResponse)
1219     FillH1(kHistStat, kStatPID);
1220   else
1221     eventGood = kFALSE;
1222
1223   // retrieve primary tracks
1224   if (fESDEvent) {
1225     fPrimTrackArrayAss = fCutsPrimAss->GetAcceptedTracks(fESDEvent);
1226     fPrimTrackArrayTrg = fCutsPrimTrg->GetAcceptedTracks(fESDEvent);
1227     if (fCutsPrimTrgConstrain) {
1228       TIter trkIter(fCutsPrimTrgConstrain->GetAcceptedTracks(fESDEvent));
1229       while (AliESDtrack *trk = (AliESDtrack*) trkIter()) {
1230         if (!fCutsPrimTrg->IsSelected(trk)) {
1231           AliESDtrack *track = (AliESDtrack*) fPrimConstrainedTrackArray->ConstructedAt(fPrimConstrainedTrackArray->GetEntriesFast());
1232           if(trk->GetConstrainedParam()) {
1233             track->Set(trk->GetConstrainedParam()->GetX(),
1234                        trk->GetConstrainedParam()->GetAlpha(),
1235                        trk->GetConstrainedParam()->GetParameter(),
1236                        trk->GetConstrainedParam()->GetCovariance());
1237           }
1238           fPrimTrackArrayTrg->Add(track);
1239         }
1240       }
1241     }
1242   }
1243   else if (fAODEvent) {
1244     // associate candidates
1245     fPrimTrackArrayAss = new TObjArray();
1246     const Int_t nTracksAODAss = fAODEvent->GetNumberOfTracks();
1247     for (Int_t iTrack = 0; iTrack < nTracksAODAss; ++iTrack) {
1248       AliAODTrack *trk = fAODEvent->GetTrack(iTrack);
1249       if (trk->TestFilterBit(fAssFilterMask))
1250         fPrimTrackArrayAss->Add(trk);
1251     }
1252
1253     // trigger candidates
1254     fPrimTrackArrayTrg = new TObjArray();
1255     const Int_t nTracksAODTrg = fAODEvent->GetNumberOfTracks();
1256     for (Int_t iTrack = 0; iTrack < nTracksAODTrg; ++iTrack) {
1257       AliAODTrack *trk = fAODEvent->GetTrack(iTrack);
1258       if (trk->IsHybridGlobalConstrainedGlobal())
1259         fPrimTrackArrayTrg->Add(trk);
1260     }
1261   }
1262   else
1263     eventGood = kFALSE;
1264
1265   // retrieve jet array
1266   if (fAODEvent) {
1267     fJetArray = dynamic_cast<TClonesArray*> (fAODEvent->FindListObject(fJetBranchName));
1268     if (!fJetArray) {
1269       // still try the output event
1270       if (AODEvent())
1271         fJetArray = dynamic_cast<TClonesArray*> (AODEvent()->FindListObject(fJetBranchName));
1272       if (!fJetArray) {
1273         if (fErrorMsg > 0) {
1274           AliError(Form("no jet branch \"%s\" found, in the AODs are:", fJetBranchName));
1275           --fErrorMsg;
1276         }
1277         if (fDebug > 0) {
1278           fAODEvent->GetList()->Print();
1279           if (AODEvent())
1280             AODEvent()->GetList()->Print();
1281         }
1282       }
1283     }
1284   }
1285
1286   // retrieve event pool for the event category
1287   if (eventGood) {
1288     for (Int_t iAss = 0; iAss <= kAssProt; ++iAss) {
1289       AliEventPoolManager *mgr = GetPoolMgr((Ass_t) iAss);
1290       GetPool((Ass_t) iAss) =
1291         // mgr ? mgr->GetEventPool(fCentrality, fZvtx, fEventPlaneAngle) : 0x0;
1292         mgr ? mgr->GetEventPool(fCentrality, fZvtx) : 0x0;
1293     }
1294   }
1295
1296   return eventGood;
1297 }
1298
1299 Float_t AliAnalysisTaskJetProtonCorr::GetPhiRel2(AliVParticle *part) const
1300 {
1301   // calculate the angle to the event plane
1302   // after removing the particle's own contribution to the Q vector
1303
1304   AliVTrack *track = dynamic_cast<AliVTrack*> (part);
1305   AliAODJet *jet   = dynamic_cast<AliAODJet*> (part);
1306
1307   Float_t phiRel = -1.;
1308
1309   if (!fEventplane)
1310     return phiRel;
1311
1312   const TVector2 *qVectorOrig = fEventplane->GetQVector();
1313   if (!qVectorOrig)
1314     return phiRel;
1315
1316   TVector2 qVector = *qVectorOrig;
1317
1318   // ??? protect against negative ID
1319   // but should be handled properly by event plane
1320   if (track && (track->GetID() > -1)) {
1321     TVector2 evplaneContrib(fEventplane->GetQContributionX(track),
1322                             fEventplane->GetQContributionY(track));
1323     qVector -= evplaneContrib;
1324   }
1325   else if (jet) {
1326     const Int_t nRefTracks = jet->GetRefTracks()->GetEntriesFast();
1327     for (Int_t iTrack = 0; iTrack < nRefTracks; ++iTrack) {
1328       AliVTrack *contrib = (AliVTrack*) jet->GetRefTracks()->At(iTrack);
1329
1330       if (contrib && (contrib->GetID() > -1)) {
1331         TVector2 evplaneContrib(fEventplane->GetQContributionX(contrib),
1332                                 fEventplane->GetQContributionY(contrib));
1333         qVector -= evplaneContrib;
1334       }
1335     }
1336   }
1337
1338   // assign phiRel and map to [0, 2 pi)
1339   phiRel = part->Phi() - (qVector.Phi() / 2.);
1340   // if (gRandom->Rndm() > .5)
1341   //    phiRel -= TMath::Pi();
1342   while (phiRel < 0.)
1343     phiRel += TMath::TwoPi();
1344   while (phiRel > TMath::TwoPi())
1345     phiRel -= TMath::TwoPi();
1346
1347   return phiRel;
1348 }
1349
1350 Bool_t AliAnalysisTaskJetProtonCorr::AcceptTrigger(AliVTrack *trg)
1351 {
1352   if (TMath::Abs(trg->Eta()) > fHadEtaMax)
1353     return kFALSE;
1354
1355   // check pt interval
1356   Bool_t acceptPt =
1357     (trg->Pt() >= fTrgPartPtMin) &&
1358     (trg->Pt() <= fTrgPartPtMax);
1359
1360   // for semi-central events check phi w.r.t. event plane
1361   Bool_t acceptOrientation = IsSemiCentral() ?
1362     AcceptAngleToEvPlane(trg->Phi(), GetEventPlaneAngle()) :
1363     kTRUE;
1364
1365   if (acceptPt) {
1366     Float_t phiRel = trg->Phi() - fEventPlaneAngle;
1367     // if (gRandom->Rndm() > .5)
1368     //   phiRel -= TMath::Pi();
1369     if (phiRel < 0.)
1370       phiRel += 2. * TMath::Pi();
1371     Float_t phiRel3 = trg->Phi() - fEventPlaneAngle3;
1372     if (phiRel3 < 0.)
1373       phiRel3 += 2. * TMath::Pi();
1374     FillH2(kHistPhiTrgHadEvPlane, GetPhiRel2(trg), fCentrality);
1375     FillH2(kHistPhiTrgHadEvPlane3, phiRel3, fCentrality);
1376   }
1377
1378   return (acceptPt && acceptOrientation);
1379 }
1380
1381 AliVTrack* AliAnalysisTaskJetProtonCorr::GetLeadingTrack(const AliAODJet *jet) const
1382 {
1383   // return leading track within a jet
1384
1385   // check contributing tracks
1386   Int_t nJetTracks = jet->GetRefTracks()->GetEntriesFast();
1387   Int_t iLeadingTrack = -1;
1388   Float_t ptLeadingTrack = 0.;
1389   for (Int_t iTrack=0; iTrack < nJetTracks; ++iTrack) {
1390     AliAODTrack *track = (AliAODTrack*) jet->GetRefTracks()->At(iTrack);
1391     // find the leading track
1392     if (track->Pt() > ptLeadingTrack) {
1393       ptLeadingTrack = track->Pt();
1394       iLeadingTrack = iTrack;
1395     }
1396   }
1397
1398   // retrieve the leading track
1399   return (AliAODTrack*) jet->GetRefTracks()->At(iLeadingTrack);
1400 }
1401
1402 Bool_t AliAnalysisTaskJetProtonCorr::AcceptTrigger(AliAODJet *trg)
1403 {
1404   // restrict eta
1405   if (TMath::Abs(trg->Eta()) > fTrgJetEtaMax)
1406     return kFALSE;
1407
1408   // reject too small jets
1409   if (trg->EffectiveAreaCharged() < fTrgJetAreaMin)
1410     return kFALSE;
1411
1412   // require hard leading track
1413   // (leading track biased jet sample)
1414   // but reject jets with too high pt constituents
1415   const Float_t ptLeadTrack = GetLeadingTrack(trg)->Pt();
1416   if ((ptLeadTrack < fTrgJetLeadTrkPtMin) ||
1417       (ptLeadTrack > fTrgJetLeadTrkPtMax))
1418     return kFALSE;
1419
1420   // check for jet orientation w.r.t. event plane
1421   // (constrained only for semi-central events)
1422   Bool_t acceptOrientation = IsSemiCentral() ?
1423     AcceptAngleToEvPlane(trg->Phi(), GetEventPlaneAngle()) :
1424     kTRUE;
1425
1426   // check for jet pt
1427   Bool_t acceptPt =
1428     (trg->Pt() >= fTrgJetPtMin) &&
1429     (trg->Pt() <= fTrgJetPtMax);
1430
1431   if (acceptPt) {
1432     // store the azimuthal distribution relative to the event plane
1433     // for the jets in the pt range of interest
1434     Float_t phiRel = trg->Phi() - fEventPlaneAngle;
1435     // if (gRandom->Rndm() > .5)
1436     //   phiRel -= TMath::Pi();
1437     if (phiRel < 0.)
1438       phiRel += 2. * TMath::Pi();
1439     Float_t phiRel3 = trg->Phi() - fEventPlaneAngle3;
1440     if (phiRel3 < 0.)
1441       phiRel3 += 2. * TMath::Pi();
1442     FillH2(kHistPhiTrgJetEvPlane, GetPhiRel2(trg), fCentrality);
1443     FillH2(kHistPhiTrgJetEvPlane3, phiRel3, fCentrality);
1444   }
1445
1446   if (acceptOrientation) {
1447     // spectrum for cross checks
1448     if (IsClass(kClCentral))
1449       FillH1(kHistJetPtCentral, trg->Pt());
1450     if (IsClass(kClSemiCentral))
1451       FillH1(kHistJetPtSemi, trg->Pt());
1452   }
1453
1454   return (acceptPt && acceptOrientation);
1455 }
1456
1457 Bool_t AliAnalysisTaskJetProtonCorr::AcceptAssoc(const AliVTrack *trk) const
1458 {
1459   if ((trk->Pt() > fAssPartPtMin) && (trk->Pt() < fAssPartPtMax) && (TMath::Abs(trk->Eta()) < fHadEtaMax))
1460     if (fRequirePID) {
1461       if ((fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC, trk) == AliPIDResponse::kDetPidOk) &&
1462           (fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF, trk) == AliPIDResponse::kDetPidOk) &&
1463           (trk->GetTPCsignalN() >= 60.) && (trk->GetTPCsignal() >= 10))
1464         return kTRUE;
1465       else
1466         return kFALSE;
1467     }
1468     else
1469       return kTRUE;
1470   else
1471     return kFALSE;
1472 }
1473
1474 Bool_t AliAnalysisTaskJetProtonCorr::IsProton(const AliVTrack *trk) const
1475 {
1476   if ((fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC, trk) != AliPIDResponse::kDetPidOk) ||
1477       (fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF, trk) != AliPIDResponse::kDetPidOk) ||
1478       (trk->GetTPCsignalN() < 60.) || (trk->GetTPCsignal() < 10))
1479     return kFALSE;
1480
1481   Double_t nSigmaProtonTPC = fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton);
1482   Double_t nSigmaProtonTOF = fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton);
1483
1484   if ((TMath::Abs(nSigmaProtonTPC) <= 2.) && (TMath::Abs(nSigmaProtonTOF) <= 2.)) {
1485     return kTRUE;
1486   }
1487
1488   return kFALSE;
1489 }
1490
1491 Bool_t AliAnalysisTaskJetProtonCorr::AcceptAngleToEvPlane(Float_t phi, Float_t psi) const
1492 {
1493   Float_t deltaPhi = phi - psi;
1494   while (deltaPhi < 0.)
1495     deltaPhi += 2 * TMath::Pi();
1496
1497   // map to interval [-pi/2, pi/2)
1498   deltaPhi = std::fmod(deltaPhi + TMath::Pi()/2., TMath::Pi()) - TMath::Pi()/2.;
1499   // printf("delta phi = %f with phi = %f, psi = %f\n", deltaPhi, phi, psi);
1500
1501   if (TMath::Abs(deltaPhi) < fTrgAngleToEvPlane)
1502     return kTRUE;
1503   else
1504     return kFALSE;
1505 }
1506
1507 Float_t AliAnalysisTaskJetProtonCorr::GetDPhiStar(Float_t phi1, Float_t pt1, Float_t charge1,
1508                                                   Float_t phi2, Float_t pt2, Float_t charge2,
1509                                                   Float_t radius, Float_t bSign) const
1510 {
1511   // calculates dphistar
1512   // from AliUEHistograms
1513
1514   Float_t dphistar = phi1 - phi2
1515     - charge1 * bSign * TMath::ASin(0.075 * radius / pt1)
1516     + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
1517
1518   static const Double_t kPi = TMath::Pi();
1519
1520   // circularity
1521   if (dphistar > kPi)
1522     dphistar = kPi * 2 - dphistar;
1523   if (dphistar < -kPi)
1524     dphistar = -kPi * 2 - dphistar;
1525   if (dphistar > kPi) // might look funny but is needed
1526     dphistar = kPi * 2 - dphistar;
1527
1528   return dphistar;
1529 }
1530
1531
1532 Bool_t AliAnalysisTaskJetProtonCorr::AcceptTwoTracks(const AliVParticle *trgPart, const AliVParticle *assPart) const
1533 {
1534   // apply two track pair cut
1535
1536   Float_t phi1 = trgPart->Phi();
1537   Float_t pt1 = trgPart->Pt();
1538   Float_t charge1 = trgPart->Charge();
1539
1540   Float_t phi2 = assPart->Phi();
1541   Float_t pt2 = assPart->Pt();
1542   Float_t charge2 = assPart->Charge();
1543
1544   Float_t deta = trgPart->Eta() - assPart->Eta();
1545
1546   Float_t bSign = (InputEvent()->GetMagneticField() > 0) ? 1 : -1;
1547
1548   // optimization
1549   if (TMath::Abs(deta) < fCutsTwoTrackEff) {
1550     // check first boundaries to see if is worth to loop and find the minimum
1551     Float_t dphistar1 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 0.8, bSign);
1552     Float_t dphistar2 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 2.5, bSign);
1553
1554     const Float_t kLimit = fCutsTwoTrackEff * 3;
1555
1556     Float_t dphistarminabs = 1e5;
1557     // Float_t dphistarmin = 1e5;
1558     if ((TMath::Abs(dphistar1) < kLimit) ||
1559         (TMath::Abs(dphistar2) < kLimit) ||
1560         (dphistar1 * dphistar2 < 0)) {
1561       for (Double_t rad=0.8; rad<2.51; rad+=0.01) {
1562         Float_t dphistarabs = TMath::Abs(GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, rad, bSign));
1563         if (dphistarabs < dphistarminabs)
1564           dphistarminabs = dphistarabs;
1565       }
1566
1567       if (dphistarminabs < fCutsTwoTrackEff)
1568         return kFALSE;
1569     }
1570   }
1571
1572   return kTRUE;
1573 }
1574
1575 Bool_t AliAnalysisTaskJetProtonCorr::Correlate(CorrType_t corr, Class_t cl, Ev_t ev,
1576                                                TCollection *trgArray, TCollection *assArray, Float_t weight)
1577 {
1578   AliHistCorr *histCorr = GetHistCorr(corr, cl, ev);
1579   if (!histCorr) {
1580     AliError(Form("no correlation histograms for corr %i, cl %i, ev %i",
1581                   corr, cl, ev));
1582     return kFALSE;
1583   }
1584
1585   TIter assIter(assArray);
1586   while (TObject *assObj = assIter()) {
1587     if (dynamic_cast<AliVParticle*> (assObj)) {
1588       AliVParticle *assPart = (AliVParticle*) assObj;
1589       histCorr->Ass(assPart->Phi(), assPart->Eta(),
1590                     assPart->Pt() > 0. ? (assPart->Charge() / 3.) / assPart->Pt() : 1.e-4,
1591                     weight);
1592     }
1593     else if (dynamic_cast<TLorentzVector*> (assObj)) {
1594       TLorentzVector *assVec = (TLorentzVector*) assObj;
1595       histCorr->Ass(TVector2::Phi_0_2pi(assVec->Phi()), assVec->Eta(), 0., weight);
1596     }
1597   }
1598
1599   TIter trgIter(trgArray);
1600   while (TObject *trgObj = trgIter()) {
1601     if (AliVParticle *trgPart = dynamic_cast<AliVParticle*> (trgObj)) {
1602       // count the trigger
1603       histCorr->Trigger(trgPart->Phi(), trgPart->Eta(),
1604                         trgPart->Pt() > 0. ? (trgPart->Charge() / 3.) / trgPart->Pt() : 1.e-4,
1605                         weight);
1606
1607       // loop over associates
1608       assIter.Reset();
1609       while (TObject *assObj = assIter()) {
1610         if (AliVParticle *assPart = dynamic_cast<AliVParticle*> (assObj)) {
1611           if (AcceptTwoTracks(trgPart, assPart))
1612             histCorr->Fill(trgPart, assPart, weight);
1613         }
1614         else if (TLorentzVector *assVec = dynamic_cast<TLorentzVector*> (assObj)) {
1615           AliFatal(Form("got %p, but not implemented", assVec));
1616         }
1617       }
1618     }
1619     else if (TLorentzVector *trgVec = dynamic_cast<TLorentzVector*> (trgObj)) {
1620       // count the trigger
1621       histCorr->Trigger(TVector2::Phi_0_2pi(trgVec->Phi()), trgVec->Eta(),
1622                         0.,
1623                         weight);
1624
1625       // loop over associates
1626       assIter.Reset();
1627       while (TObject *assObj = assIter()) {
1628         if (AliVParticle *assPart = dynamic_cast<AliVParticle*> (assObj)) {
1629           histCorr->Fill(trgVec, assPart, weight);
1630         }
1631         else if (TLorentzVector *assVec = dynamic_cast<TLorentzVector*> (assObj)) {
1632           histCorr->Fill(trgVec, assVec, weight);
1633         }
1634       }
1635     }
1636   }
1637
1638   return kTRUE;
1639 }
1640
1641 Bool_t AliAnalysisTaskJetProtonCorr::Correlate(Trg_t trg, Ass_t ass, Class_t cl, Ev_t ev,
1642                                                TCollection *trgArray, TCollection *assArray, Float_t weight)
1643 {
1644   if ((trg < 0) || (trg > kTrgJetRnd))
1645     AliFatal(Form("wrong request for correlation with trigger: %d", trg));
1646
1647   if ((ass < 0) || (ass > kAssProt))
1648     AliFatal(Form("wrong request for correlation with associate: %d", ass));
1649
1650   CorrType_t corr = (CorrType_t) (2 * trg + ass);
1651
1652   return Correlate(corr, cl, ev, trgArray, assArray, weight);
1653 }
1654
1655 Bool_t AliAnalysisTaskJetProtonCorr::GenerateRandom(TCollection *trgJetArray,
1656                                                     TCollection *trgHadArray,
1657                                                     TCollection *assHadJetArray,
1658                                                     TCollection *assProtJetArray,
1659                                                     TCollection *assHadHadArray,
1660                                                     TCollection *assProtHadArray,
1661                                                     Float_t pFraction)
1662 {
1663   const Int_t nPart = gRandom->Poisson(fToyMeanNoPart);
1664
1665   // overcompensate modulation for event plane resolution
1666   if (fSplineEventPlaneRes) {
1667     fTrgJetPhiModCent->SetParameter(0, fTrgJetV2Cent / fSplineEventPlaneRes->Eval(GetCentrality()));
1668     fTrgJetPhiModSemi->SetParameter(0, fTrgJetV2Semi / fSplineEventPlaneRes->Eval(GetCentrality()));
1669     fTrgHadPhiModCent->SetParameter(0, fTrgHadV2Cent / fSplineEventPlaneRes->Eval(GetCentrality()));
1670     fTrgHadPhiModSemi->SetParameter(0, fTrgHadV2Semi / fSplineEventPlaneRes->Eval(GetCentrality()));
1671   }
1672   else {
1673     printf("no event plane resolution spline available\n");
1674     fTrgJetPhiModCent->SetParameter(0, fTrgJetV2Cent);
1675     fTrgJetPhiModSemi->SetParameter(0, fTrgJetV2Semi);
1676     fTrgHadPhiModCent->SetParameter(0, fTrgHadV2Cent);
1677     fTrgHadPhiModSemi->SetParameter(0, fTrgHadV2Semi);
1678   }
1679
1680   // generate random direction
1681   Float_t trgJetEta = gRandom->Uniform(-fTrgJetEtaMax, fTrgJetEtaMax);
1682   Float_t trgJetPhi = 0.;
1683   do {
1684     trgJetPhi = IsClass(kClSemiCentral) ?
1685       fTrgJetPhiModSemi->GetRandom() :
1686       fTrgJetPhiModCent->GetRandom();
1687
1688     trgJetPhi += fEventPlaneAngle;
1689     if (trgJetPhi < 0.)
1690       trgJetPhi += 2. * TMath::Pi();
1691     else if (trgJetPhi > 2. * TMath::Pi())
1692       trgJetPhi -= 2. * TMath::Pi();
1693
1694     Float_t phiRel = trgJetPhi- fEventPlaneAngle;
1695     if (phiRel < 0.)
1696       phiRel += 2. * TMath::Pi();
1697     FillH2(kHistPhiRndTrgJetEvPlane, phiRel, fCentrality);
1698   } while (IsClass(kClSemiCentral) &&
1699            !AcceptAngleToEvPlane(trgJetPhi, GetEventPlaneAngle()));
1700
1701   // generate one trigger particle
1702   TLorentzVector *trgJet = new TLorentzVector();
1703   trgJet->SetPtEtaPhiM(50., trgJetEta, trgJetPhi, .14);
1704   trgJetArray->Add(trgJet);
1705
1706   // generate direction for away side
1707   Float_t trgJetEtaAway = gRandom->Uniform(-fHadEtaMax, fHadEtaMax);
1708   Float_t trgJetPhiAway = trgJetPhi + TMath::Pi() + gRandom->Gaus(0., fToySmearPhi);
1709
1710   // generate associated particles
1711   // with proton/hadron ratio observed in this event
1712   for (Int_t iPart = 0; iPart < nPart; ++iPart) {
1713     Float_t deltaEta, deltaPhi;
1714     const Float_t r = fToyRadius;
1715     do {
1716       deltaEta = gRandom->Uniform(-r, r);
1717       deltaPhi = gRandom->Uniform(-r, r);
1718     } while ((deltaEta * deltaEta + deltaPhi * deltaPhi) > (r * r));
1719
1720     TLorentzVector *assPart = new TLorentzVector();
1721     Float_t eta = trgJetEtaAway + deltaEta;
1722     Float_t phi = trgJetPhiAway + deltaPhi;
1723     if (TMath::Abs(eta) > fHadEtaMax) {
1724       delete assPart;
1725       continue;
1726     }
1727     while (phi < 0.)
1728       phi += 2. * TMath::Pi();
1729     while (phi > 2. * TMath::Pi())
1730       phi -= 2. * TMath::Pi();
1731     assPart->SetPtEtaPhiM(3., eta, phi, .14);
1732     assHadJetArray->Add(assPart);
1733     if (gRandom->Uniform() < pFraction)
1734       assProtJetArray->Add(assPart);
1735   }
1736
1737   // trigger hadron
1738   Float_t trgHadEta = gRandom->Uniform(-fHadEtaMax, fHadEtaMax);
1739   Float_t trgHadPhi = 0.;
1740   do {
1741     trgHadPhi = IsClass(kClSemiCentral) ?
1742       fTrgHadPhiModSemi->GetRandom() :
1743       fTrgHadPhiModCent->GetRandom();
1744
1745     trgHadPhi += fEventPlaneAngle;
1746     if (trgHadPhi < 0.)
1747       trgHadPhi += 2. * TMath::Pi();
1748     else if (trgHadPhi > 2. * TMath::Pi())
1749       trgHadPhi -= 2. * TMath::Pi();
1750
1751     Float_t phiRel = trgHadPhi - fEventPlaneAngle;
1752     if (phiRel < 0.)
1753       phiRel += 2. * TMath::Pi();
1754     FillH2(kHistPhiRndTrgHadEvPlane, phiRel, fCentrality);
1755   } while (IsClass(kClSemiCentral) &&
1756            !AcceptAngleToEvPlane(trgHadPhi, GetEventPlaneAngle()));
1757
1758   // generate one trigger particle
1759   TLorentzVector *trgHad = new TLorentzVector();
1760   trgHad->SetPtEtaPhiM(50., trgHadEta, trgHadPhi, .14);
1761   trgHadArray->Add(trgHad);
1762
1763   // generate direction for away side
1764   Float_t trgHadEtaAway = gRandom->Uniform(-fHadEtaMax, fHadEtaMax);
1765   Float_t trgHadPhiAway = trgHadPhi + TMath::Pi() + gRandom->Gaus(0., fToySmearPhi);
1766
1767   // generate associated particles
1768   // with proton/hadron ratio observed in this event
1769   for (Int_t iPart = 0; iPart < nPart; ++iPart) {
1770     Float_t deltaEta, deltaPhi;
1771     const Float_t r = fToyRadius;
1772     do {
1773       deltaEta = gRandom->Uniform(-r, r);
1774       deltaPhi = gRandom->Uniform(-r, r);
1775     } while ((deltaEta * deltaEta + deltaPhi * deltaPhi) > (r * r));
1776
1777     TLorentzVector *assPart = new TLorentzVector();
1778     Float_t eta = trgHadEtaAway + deltaEta;
1779     Float_t phi = trgHadPhiAway + deltaPhi;
1780     if (TMath::Abs(eta) > fHadEtaMax) {
1781       delete assPart;
1782       continue;
1783     }
1784     while (phi < 0.)
1785       phi += 2. * TMath::Pi();
1786     while (phi > 2. * TMath::Pi())
1787       phi -= 2. * TMath::Pi();
1788     assPart->SetPtEtaPhiM(3., eta, phi, .14);
1789     assHadHadArray->Add(assPart);
1790     if (gRandom->Uniform() < pFraction)
1791       assProtHadArray->Add(assPart);
1792   }
1793
1794   return kTRUE;
1795 }
1796
1797 Bool_t AliAnalysisTaskJetProtonCorr::CleanUpEvent()
1798 {
1799   if (fAODEvent) {
1800     delete fPrimTrackArrayAss;
1801     fPrimTrackArrayAss = 0x0;
1802     delete fPrimTrackArrayTrg;
1803     fPrimTrackArrayTrg = 0x0;
1804   }
1805
1806   fPrimConstrainedTrackArray->Delete();
1807
1808   return kTRUE;
1809 }
1810
1811 TObjArray* AliAnalysisTaskJetProtonCorr::CloneTracks(TObjArray* tracks) const
1812 {
1813   TObjArray* tracksClone = new TObjArray;
1814   tracksClone->SetOwner(kTRUE);
1815
1816   Int_t nTracks = tracks->GetEntriesFast();
1817   for (Int_t i = 0; i < nTracks; i++) {
1818     // WARNING: TObject::Clone() is very!!! expensive, unusable
1819     // tracksClone->Add(particle->Clone());
1820
1821     // if (AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (tracks->At(i)))
1822     //   tracksClone->Add(new AliESDtrack(*esdTrack));
1823     // else if (AliAODTrack* aodTrack = dynamic_cast<AliAODTrack*> (tracks->At(i)))
1824     //   tracksClone->Add(new AliAODTrack(*aodTrack));
1825
1826     if (const AliVTrack *track = dynamic_cast<const AliVTrack*> (tracks->At(i)))
1827       tracksClone->Add(new AliPartCorr(*track));
1828   }
1829
1830   return tracksClone;
1831 }
1832
1833 AliAnalysisTaskJetProtonCorr::AliHistCorr::AliHistCorr(TString name, TList *outputList) :
1834   TNamed(name, name),
1835   fOutputList(outputList),
1836   fHistStat(0x0),
1837   fHistCorrPhi(0x0),
1838   fHistCorrPhi2(0x0),
1839   fHistCorrEtaPhi(0x0),
1840   fHistCorrAvgEtaPhi(0x0),
1841   fHistCorrTrgEtaPhi(0x0),
1842   fHistCorrAssEtaPhi(0x0),
1843   fHistCorrTrgEtaPhiQpt(0x0),
1844   fHistCorrAssEtaPhiQpt(0x0),
1845   fHistDphiLo(-TMath::Pi() / 2.),
1846   fHistDphiNbins(120),
1847   fHistDetaNbins(120)
1848 {
1849   // ctor
1850
1851   fHistStat = new TH1F(Form("%s_stat", name.Data()), "statistics",
1852                        1, .5, 1.5);
1853
1854   fHistCorrPhi = new TH1F(Form("%s_phi", name.Data()), ";#Delta#varphi",
1855                           fHistDphiNbins, fHistDphiLo, 2.*TMath::Pi() + fHistDphiLo);
1856   fHistCorrPhi->Sumw2();
1857   fHistCorrPhi2 = new TH2F(Form("%s_phi2", name.Data()), ";#varphi_{trg};#varphi_{ass}",
1858                           120, 0.*TMath::Pi(), 2.*TMath::Pi(),
1859                           120, 0.*TMath::Pi(), 2.*TMath::Pi());
1860   fHistCorrPhi2->Sumw2();
1861   fHistCorrEtaPhi = new TH2F(Form("%s_etaphi", name.Data()), ";#Delta#varphi;#Delta#eta",
1862                              fHistDphiNbins, fHistDphiLo, 2*TMath::Pi() + fHistDphiLo,
1863                              fHistDetaNbins, -2., 2.);
1864   fHistCorrEtaPhi->Sumw2();
1865   fHistCorrAvgEtaPhi = new TH2F(Form("%s_avgetaphi", name.Data()), ";#Delta#varphi;avg #eta",
1866                                 fHistDphiNbins, fHistDphiLo, 2*TMath::Pi() + fHistDphiLo,
1867                                 fHistDetaNbins, -2., 2.);
1868   fHistCorrAvgEtaPhi->Sumw2();
1869   fHistCorrTrgEtaPhi = new TH2F(Form("%s_trg_etaphi", name.Data()), ";#varphi;#eta",
1870                                 120, 0., 2*TMath::Pi(),
1871                                 120, -1., 1.);
1872   fHistCorrTrgEtaPhi->Sumw2();
1873   fHistCorrAssEtaPhi = new TH2F(Form("%s_ass_etaphi", name.Data()), ";#varphi;#eta",
1874                                 120, 0., 2*TMath::Pi(),
1875                                 120, -1., 1.);
1876   fHistCorrAssEtaPhi->Sumw2();
1877
1878   // fHistCorrTrgEtaPhiQpt = new TH3F(Form("%s_trg_etaphi", name.Data()), ";#varphi;#eta;Q/p_{T}",
1879   //                               120, 0., 2*TMath::Pi(),
1880   //                               120, -1., 1.,
1881   //                               120, -1., 1.);
1882   // fHistCorrTrgEtaPhiQpt->Sumw2();
1883   // fHistCorrAssEtaPhiQpt = new TH3F(Form("%s_ass_etaphi", name.Data()), ";#varphi;#etaQ/p_{T}",
1884   //                               120, 0., 2*TMath::Pi(),
1885   //                               120, -1., 1.,
1886   //                               120, -1., 1.);
1887   // fHistCorrAssEtaPhiQpt->Sumw2();
1888
1889   fOutputList->Add(fHistStat);
1890   fOutputList->Add(fHistCorrPhi);
1891   fOutputList->Add(fHistCorrPhi2);
1892   fOutputList->Add(fHistCorrEtaPhi);
1893   fOutputList->Add(fHistCorrAvgEtaPhi);
1894   fOutputList->Add(fHistCorrTrgEtaPhi);
1895   fOutputList->Add(fHistCorrAssEtaPhi);
1896   // fOutputList->Add(fHistCorrTrgEtaPhiQpt);
1897   // fOutputList->Add(fHistCorrAssEtaPhiQpt);
1898 }
1899
1900 AliAnalysisTaskJetProtonCorr::AliHistCorr::~AliHistCorr()
1901 {
1902   // dtor
1903 }
1904
1905 void AliAnalysisTaskJetProtonCorr::AliHistCorr::Fill(AliVParticle *trgPart, AliVParticle *assPart, Float_t weight)
1906 {
1907   // fill correlation histograms from trigger particle and associate particle
1908
1909   Float_t deltaEta = assPart->Eta() - trgPart->Eta();
1910   Float_t avgEta = (assPart->Eta() + trgPart->Eta()) / 2.;
1911   Float_t deltaPhi = assPart->Phi() - trgPart->Phi();
1912   while (deltaPhi > (2.*TMath::Pi() + fHistDphiLo))
1913     deltaPhi -= 2. * TMath::Pi();
1914   while (deltaPhi < fHistDphiLo)
1915     deltaPhi += 2. * TMath::Pi();
1916
1917   fHistCorrPhi->Fill(deltaPhi, weight);
1918   if (fHistCorrPhi2)
1919     fHistCorrPhi2->Fill(trgPart->Phi(), assPart->Phi(), weight);
1920   fHistCorrEtaPhi->Fill(deltaPhi, deltaEta, weight);
1921   if (fHistCorrAvgEtaPhi)
1922     fHistCorrAvgEtaPhi->Fill(deltaPhi, avgEta, weight);
1923 }
1924
1925 void AliAnalysisTaskJetProtonCorr::AliHistCorr::Fill(TLorentzVector *trgPart, AliVParticle *assPart, Float_t weight)
1926 {
1927   // fill correlation histograms from trigger direction and associate particle
1928
1929   Float_t deltaEta = assPart->Eta() - trgPart->Eta();
1930   Float_t avgEta = (assPart->Eta() + trgPart->Eta()) / 2.;
1931   Float_t deltaPhi = assPart->Phi() - trgPart->Phi();
1932   while (deltaPhi > (2.*TMath::Pi() + fHistDphiLo))
1933     deltaPhi -= 2. * TMath::Pi();
1934   while (deltaPhi < fHistDphiLo)
1935     deltaPhi += 2. * TMath::Pi();
1936
1937   fHistCorrPhi->Fill(deltaPhi, weight);
1938   if (fHistCorrPhi2) {
1939     Float_t trgPhi = trgPart->Phi();
1940     if (trgPhi < 0)
1941       trgPhi += 2. * TMath::Pi();
1942     fHistCorrPhi2->Fill(trgPhi, assPart->Phi(), weight);
1943   }
1944   fHistCorrEtaPhi->Fill(deltaPhi, deltaEta, weight);
1945   if (fHistCorrAvgEtaPhi)
1946     fHistCorrAvgEtaPhi->Fill(deltaPhi, avgEta, weight);
1947 }
1948
1949 void AliAnalysisTaskJetProtonCorr::AliHistCorr::Fill(TLorentzVector *trgPart, TLorentzVector *assPart, Float_t weight)
1950 {
1951   // fill correlation histograms from trigger direction and associate direction
1952
1953   Float_t deltaEta = assPart->Eta() - trgPart->Eta();
1954   Float_t avgEta = (assPart->Eta() + trgPart->Eta()) / 2.;
1955   Float_t deltaPhi = assPart->Phi() - trgPart->Phi();
1956   if (deltaPhi > (2.*TMath::Pi() + fHistDphiLo))
1957     deltaPhi -= 2. * TMath::Pi();
1958   else if (deltaPhi < fHistDphiLo)
1959     deltaPhi += 2. * TMath::Pi();
1960
1961   fHistCorrPhi->Fill(deltaPhi, weight);
1962   if (fHistCorrPhi2) {
1963     Float_t trgPhi = trgPart->Phi();
1964     if (trgPhi < 0)
1965       trgPhi += 2. * TMath::Pi();
1966     Float_t assPhi = assPart->Phi();
1967     if (assPhi < 0)
1968       assPhi += 2. * TMath::Pi();
1969     fHistCorrPhi2->Fill(trgPhi, assPhi, weight);
1970   }
1971   fHistCorrEtaPhi->Fill(deltaPhi, deltaEta, weight);
1972   if (fHistCorrAvgEtaPhi)
1973     fHistCorrAvgEtaPhi->Fill(deltaPhi, avgEta, weight);
1974 }
1975
1976 // ----- histogram management -----
1977 TH1* AliAnalysisTaskJetProtonCorr::AddHistogram(Hist_t hist, const char *hid, TString title,
1978                                                  Int_t xbins, Float_t xmin, Float_t xmax,
1979                                                  Int_t binType)
1980 {
1981   TString hName;
1982   hName.Form("%s_%s", fShortTaskId, hid);
1983   hName.ToLower();
1984   TH1 *h = 0x0;
1985   if (binType == 0)
1986     h = new TH1I(hName.Data(), title,
1987                  xbins, xmin, xmax);
1988   else
1989     h = new TH1F(hName.Data(), title,
1990                  xbins, xmin, xmax);
1991   GetHistogram(hist) = h;
1992   fOutputList->Add(h);
1993   return h;
1994 }
1995
1996 TH2* AliAnalysisTaskJetProtonCorr::AddHistogram(Hist_t hist, const char *hid, TString title,
1997                                                  Int_t xbins, Float_t xmin, Float_t xmax,
1998                                                  Int_t ybins, Float_t ymin, Float_t ymax,
1999                                                  Int_t binType)
2000 {
2001   TString hName;
2002   hName.Form("%s_%s", fShortTaskId, hid);
2003   hName.ToLower();
2004   TH1 *h = 0x0;
2005   if (binType == 0)
2006     h = GetHistogram(hist) = new TH2I(hName.Data(), title,
2007                                      xbins, xmin, xmax,
2008                                      ybins, ymin, ymax);
2009   else
2010     h = GetHistogram(hist) = new TH2F(hName.Data(), title,
2011                                      xbins, xmin, xmax,
2012                                      ybins, ymin, ymax);
2013   fOutputList->Add(h);
2014   return (TH2*) h;
2015 }
2016
2017 TH3* AliAnalysisTaskJetProtonCorr::AddHistogram(Hist_t hist, const char *hid, TString title,
2018                                                  Int_t xbins, Float_t xmin, Float_t xmax,
2019                                                  Int_t ybins, Float_t ymin, Float_t ymax,
2020                                                  Int_t zbins, Float_t zmin, Float_t zmax,
2021                                                  Int_t binType)
2022 {
2023   TString hName;
2024   hName.Form("%s_%s", fShortTaskId, hid);
2025   hName.ToLower();
2026   TH1 *h = 0x0;
2027   if (binType == 0)
2028     h = GetHistogram(hist) = new TH3I(hName.Data(), title,
2029                                      xbins, xmin, xmax,
2030                                      ybins, ymin, ymax,
2031                                      zbins, zmin, zmax);
2032   else
2033     h = GetHistogram(hist) = new TH3F(hName.Data(), title,
2034                                      xbins, xmin, xmax,
2035                                      ybins, ymin, ymax,
2036                                      zbins, zmin, zmax);
2037   fOutputList->Add(h);
2038   return (TH3*) h;
2039 }