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 "TFormula.h"
   8 #include "TRandom.h"
   9
  10 // analysis framework
  11 #include "AliAnalysisManager.h"
  12 #include "AliInputEventHandler.h"
  13 #include "AliVEvent.h"
  14 #include "AliVTrdTrack.h"
  15 #include "AliVVertex.h"
  16 #include "AliPIDResponse.h"
  17 #include "AliEventPoolManager.h"
  18
  19 // MC stuff
  20 #include "AliMCEvent.h"
  21 #include "AliGenPythiaEventHeader.h"
  22
  23 // ESD stuff
  24 #include "AliESDEvent.h"
  25 #include "AliESDInputHandler.h"
  26 #include "AliESDtrack.h"
  27 #include "AliESDtrackCuts.h"
  28 #include "AliESDTrdTrack.h"
  29 #include "AliESDTrdTracklet.h"
  30 #include "AliESDTrdTrigger.h"
  31
  32 // AOD stuff
  33 #include "AliAODEvent.h"
  34 #include "AliAODJet.h"
  35 #include "AliAODTrack.h"
  36
  37 // jet tasks
  38 #include "AliAnalysisTaskJetServices.h"
  39 #include "AliAnalysisHelperJetTasks.h"
  40
  41 #include "AliAnalysisTaskJetProtonCorr.h"
  42
  43 #include <iostream>
  44 #include <cmath>
  45
  46 AliAnalysisTaskJetProtonCorr::AliAnalysisTaskJetProtonCorr(const char *name) :
  47   AliAnalysisTaskSE(name),
  48   fMCEvent(0x0),
  49   fESDEvent(0x0),
  50   fAODEvent(0x0),
  51   fTriggerMask(0),
  52   fClassMask(0),
  53   fCentrality(100.),
  54   fCentralityCheck(100.),
  55   fZvtx(0.),
  56   fPIDResponse(0x0),
  57   fEventPlane(5.),
  58   fEventPlaneCheck(5.),
  59   fPrimTrackArray(0x0),
  60   fJetArray(0x0),
  61   fPoolMgr(),
  62   fPool(),
  63   fHistCorr(0x0),
  64   fOutputList(),
  65   fHist(),
  66   fShortTaskId("jet_prot_corr"),
  67   fUseStandardCuts(kTRUE),
  68   fCutsPrim(0x0),
  69   fCutsTwoTrackEff(0.02),
  70   fTrgPartPtMin(6.),
  71   fTrgPartPtMax(8.),
  72   fTrgJetPtMin(50.),
  73   fTrgJetPtMax(80.),
  74   fTrgJetLeadTrkPtMin(5.),
  75   fAssPartPtMin(2.),
  76   fAssPartPtMax(4.),
  77   fTrgAngleToEvPlane(TMath::Pi() / 4.)
  78 {
  79   // default ctor
  80
  81   fkCorrTypeName[kCorrHadHad]  = "hh";
  82   fkCorrTypeName[kCorrHadProt] = "hp";
  83   fkCorrTypeName[kCorrJetHad]  = "jh";
  84   fkCorrTypeName[kCorrJetProt] = "jp";
  85
  86   fkClassName[kClCentral]      = "cent";
  87   fkClassName[kClSemiCentral]  = "semi";
  88   fkClassName[kClDijet]        = "dijet";
  89
  90   fkEvName[kEvSame] = "same";
  91   fkEvName[kEvMix]  = "mixed";
  92
  93   // track cuts
  94   if (fUseStandardCuts) {
  95     fCutsPrim = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE);
  96   } else {
  97     fCutsPrim = new AliESDtrackCuts();
  98
  99     // this is taken from PWGJE track cuts
 100     TFormula *f1NClustersTPCLinearPtDep = new TFormula("f1NClustersTPCLinearPtDep","70.+30./20.*x");
 101     fCutsPrim->SetMinNClustersTPCPtDep(f1NClustersTPCLinearPtDep,20.);
 102     fCutsPrim->SetMinNClustersTPC(70);
 103     fCutsPrim->SetMaxChi2PerClusterTPC(4);
 104     fCutsPrim->SetRequireTPCStandAlone(kTRUE); //cut on NClustersTPC and chi2TPC Iter1
 105     fCutsPrim->SetAcceptKinkDaughters(kFALSE);
 106     fCutsPrim->SetRequireTPCRefit(kTRUE);
 107     fCutsPrim->SetMaxFractionSharedTPCClusters(0.4);
 108     // ITS
 109     fCutsPrim->SetRequireITSRefit(kTRUE);
 110     //accept secondaries
 111     fCutsPrim->SetMaxDCAToVertexXY(2.4);
 112     fCutsPrim->SetMaxDCAToVertexZ(3.2);
 113     fCutsPrim->SetDCAToVertex2D(kTRUE);
 114     //reject fakes
 115     fCutsPrim->SetMaxChi2PerClusterITS(36);
 116     fCutsPrim->SetMaxChi2TPCConstrainedGlobal(36);
 117
 118     fCutsPrim->SetRequireSigmaToVertex(kFALSE);
 119
 120     fCutsPrim->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
 121   }
 122
 123   fCutsPrim->SetEtaRange(-0.9, 0.9);
 124   fCutsPrim->SetPtRange(0.15, 1E+15);
 125
 126   // event mixing pool
 127   Double_t centralityBins[] = {
 128     0., 2., 4., 6., 8., 10., // central
 129     30., 32., 34., 36., 38., 40., 42., 44., 46., 48., 50., // semi-central
 130     90.
 131   };
 132   Int_t nCentralityBins = sizeof(centralityBins)/sizeof(centralityBins[0]);
 133
 134   Double_t vertexBins[] = {
 135     -10., -8., -6., -4., -2., 0., 2., 4., 6., 8., 10.
 136   };
 137   Int_t nVertexBins = sizeof(vertexBins)/sizeof(vertexBins[0]);
 138
 139   Double_t psiBins[12];
 140   Int_t nPsiBins = sizeof(psiBins)/sizeof(psiBins[0]);
 141   for (Int_t iBin = 0; iBin < nPsiBins; ++iBin)
 142     psiBins[iBin] = iBin * TMath::Pi()/nPsiBins;
 143
 144   for (Int_t iTrg = 0; iTrg < kTrgLast; ++iTrg) {
 145     for (Int_t iAss = 0; iAss < kAssLast; ++iAss) {
 146       GetPoolMgr((Trg_t) iTrg, (Ass_t) iAss) =
 147         new AliEventPoolManager(10, 100,
 148                                 nCentralityBins, centralityBins,
 149                                 nVertexBins, vertexBins);
 150                                 // nPsiBins, psiBins);
 151       GetPoolMgr((Trg_t) iTrg, (Ass_t) iAss)->SetTargetValues(100, .1, 1);
 152     }
 153   }
 154
 155   fHistCorr = new AliHistCorr*[kEvLast*kCorrLast*kClLast];
 156
 157   DefineOutput(1, TList::Class());
 158 }
 159
 160 AliAnalysisTaskJetProtonCorr::~AliAnalysisTaskJetProtonCorr()
 161 {
 162   // dtor
 163
 164   // delete [] fHistCorr;
 165 }
 166
 167 void AliAnalysisTaskJetProtonCorr::UserCreateOutputObjects()
 168 {
 169   // create user output objects
 170
 171   // setup list
 172   OpenFile(1);
 173   fOutputList = new TList();
 174   fOutputList->SetOwner();
 175
 176   // setup histograms
 177   TH1 *hist;
 178   TH1 *histStat = AddHistogram(ID(kHistStat), "event statistics;;counts",
 179                                kStatLast-1, .5, kStatLast-.5);
 180   histStat->GetXaxis()->SetBinLabel(ID(kStatSeen));
 181   histStat->GetXaxis()->SetBinLabel(ID(kStatTrg));
 182   histStat->GetXaxis()->SetBinLabel(ID(kStatEvCuts));
 183   histStat->GetXaxis()->SetBinLabel(ID(kStatUsed));
 184   histStat->GetXaxis()->SetBinLabel(ID(kStatCent));
 185   histStat->GetXaxis()->SetBinLabel(ID(kStatEvPlane));
 186   histStat->GetXaxis()->SetBinLabel(ID(kStatPID));
 187   histStat->GetXaxis()->SetBinLabel(ID(kStatCentral));
 188   histStat->GetXaxis()->SetBinLabel(ID(kStatSemiCentral));
 189
 190   AddHistogram(ID(kHistCentrality), "centrality;C;counts",
 191                110, -5., 105.);
 192   hist = AddHistogram(ID(kHistCentralityUsed), "centrality used;C;event class",
 193                       110, -5., 105.,
 194                       kClLast, -.5, kClLast-.5);
 195   hist->GetYaxis()->SetBinLabel(LAB(kClCentral));
 196   hist->GetYaxis()->SetBinLabel(LAB(kClSemiCentral));
 197   hist->GetYaxis()->SetBinLabel(LAB(kClDijet));
 198   AddHistogram(ID(kHistCentralityCheck), "centrality check;C;counts",
 199                110, -5., 105.);
 200   hist = AddHistogram(ID(kHistCentralityCheckUsed), "centrality check used;C;event class",
 201                       110, -5., 105.,
 202                       kClLast, -.5, kClLast-.5);
 203   hist->GetYaxis()->SetBinLabel(LAB(kClCentral));
 204   hist->GetYaxis()->SetBinLabel(LAB(kClSemiCentral));
 205   hist->GetYaxis()->SetBinLabel(LAB(kClDijet));
 206
 207
 208   AddHistogram(ID(kHistSignalTPC), "TPC dE/dx;p (GeV/c);dE/dx (arb. units)",
 209                100, 0., 10., 200, 0., 300.);
 210   AddHistogram(ID(kHistSignalTOF), "TOF time of flight;p_{T} (GeV/c);t (ns)",
 211                100, 0., 10., 200, 0., 50.);
 212   AddHistogram(ID(kHistBetaTOF), "TOF beta;p (GeV/c); #beta",
 213                100, 0., 10.,
 214                100, 0., 1.);
 215   AddHistogram(ID(kHistDeltaTPC), "TPC dE/dx;p (GeV/c);dE/dx (arb. units)",
 216                100, 0., 10., 200, -100., 100.);
 217   AddHistogram(ID(kHistDeltaTPCSemi), "TPC dE/dx;p (GeV/c);dE/dx (arb. units)",
 218                100, 0., 10., 200, -100., 100.);
 219   AddHistogram(ID(kHistDeltaTOF), "TOF time of flight;p_{T} (GeV/c);t (ns)",
 220                100, 0., 10., 200, -2., 2.);
 221   AddHistogram(ID(kHistDeltaTOFSemi), "TOF time of flight;p_{T} (GeV/c);t (ns)",
 222                100, 0., 10., 200, -2., 2.);
 223
 224   // Nsigma templates
 225   AddHistogram(ID(kHistNsigmaTPCe), "TPC N#sigma - e hypothesis;p (GeV/c)",
 226                100, 0., 10.,
 227                100, -25., 25.);
 228   AddHistogram(ID(kHistNsigmaTPCmu), "TPC N#sigma - #mu hypothesis;p (GeV/c)",
 229                100, 0., 10.,
 230                100, -25., 25.);
 231   AddHistogram(ID(kHistNsigmaTPCpi), "TPC N#sigma - #pi hypothesis;p (GeV/c)",
 232                100, 0., 10.,
 233                100, -25., 25.);
 234   AddHistogram(ID(kHistNsigmaTPCk), "TPC N#sigma - K hypothesis;p (GeV/c)",
 235                100, 0., 10.,
 236                100, -25., 25.);
 237   AddHistogram(ID(kHistNsigmaTPCp), "TPC N#sigma - p hypothesis;p (GeV/c)",
 238                100, 0., 10.,
 239                100, -25., 25.);
 240   AddHistogram(ID(kHistNsigmaTPCd), "TPC N#sigma - d hypothesis;p (GeV/c)",
 241                100, 0., 10.,
 242                100, -25., 25.);
 243   AddHistogram(ID(kHistNsigmaTPCe_e), "TPC N#sigma - e hypothesis (id. e);p (GeV/c)",
 244                100, 0., 10.,
 245                100, -25., 25.);
 246
 247   AddHistogram(ID(kHistNsigmaTOFe), "TOF N#sigma - e hypothesis;p (GeV/c)",
 248                100, 0., 10.,
 249                200, -50., 50.);
 250   AddHistogram(ID(kHistNsigmaTOFmu), "TOF N#sigma - #mu hypothesis;p (GeV/c)",
 251                100, 0., 10.,
 252                200, -50., 50.);
 253   AddHistogram(ID(kHistNsigmaTOFpi), "TOF N#sigma - #pi hypothesis;p (GeV/c)",
 254                100, 0., 10.,
 255                200, -50., 50.);
 256   AddHistogram(ID(kHistNsigmaTOFk), "TOF N#sigma - K hypothesis;p (GeV/c)",
 257                100, 0., 10.,
 258                200, -50., 50.);
 259   AddHistogram(ID(kHistNsigmaTOFp), "TOF N#sigma - p hypothesis;p (GeV/c)",
 260                100, 0., 10.,
 261                200, -50., 50.);
 262   AddHistogram(ID(kHistNsigmaTOFd), "TOF N#sigma - d hypothesis;p (GeV/c)",
 263                100, 0., 10.,
 264                200, -50., 50.);
 265   AddHistogram(ID(kHistNsigmaTOFmismatch), "TOF N#sigma - mismatch;p (GeV/c)",
 266                100, 0., 10.,
 267                200, -50., 50.);
 268   AddHistogram(ID(kHistNsigmaTOFmismatch2), "TOF N#sigma - mismatch;p (GeV/c)",
 269                100, 0., 10.,
 270                200, -50., 50.);
 271
 272   // Nsigma templates
 273   AddHistogram(ID(kHistNsigmaTPCeSemi), "TPC N#sigma - e hypothesis;p (GeV/c)",
 274                100, 0., 10.,
 275                100, -25., 25.);
 276   AddHistogram(ID(kHistNsigmaTPCmuSemi), "TPC N#sigma - #mu hypothesis;p (GeV/c)",
 277                100, 0., 10.,
 278                100, -25., 25.);
 279   AddHistogram(ID(kHistNsigmaTPCpiSemi), "TPC N#sigma - #pi hypothesis;p (GeV/c)",
 280                100, 0., 10.,
 281                100, -25., 25.);
 282   AddHistogram(ID(kHistNsigmaTPCkSemi), "TPC N#sigma - K hypothesis;p (GeV/c)",
 283                100, 0., 10.,
 284                100, -25., 25.);
 285   AddHistogram(ID(kHistNsigmaTPCpSemi), "TPC N#sigma - p hypothesis;p (GeV/c)",
 286                100, 0., 10.,
 287                100, -25., 25.);
 288   AddHistogram(ID(kHistNsigmaTPCdSemi), "TPC N#sigma - d hypothesis;p (GeV/c)",
 289                100, 0., 10.,
 290                100, -25., 25.);
 291   AddHistogram(ID(kHistNsigmaTPCe_eSemi), "TPC N#sigma - e hypothesis (id. e);p (GeV/c)",
 292                100, 0., 10.,
 293                100, -25., 25.);
 294
 295   AddHistogram(ID(kHistNsigmaTOFeSemi), "TOF N#sigma - e hypothesis;p (GeV/c)",
 296                100, 0., 10.,
 297                200, -50., 50.);
 298   AddHistogram(ID(kHistNsigmaTOFmuSemi), "TOF N#sigma - #mu hypothesis;p (GeV/c)",
 299                100, 0., 10.,
 300                200, -50., 50.);
 301   AddHistogram(ID(kHistNsigmaTOFpiSemi), "TOF N#sigma - #pi hypothesis;p (GeV/c)",
 302                100, 0., 10.,
 303                200, -50., 50.);
 304   AddHistogram(ID(kHistNsigmaTOFkSemi), "TOF N#sigma - K hypothesis;p (GeV/c)",
 305                100, 0., 10.,
 306                200, -50., 50.);
 307   AddHistogram(ID(kHistNsigmaTOFpSemi), "TOF N#sigma - p hypothesis;p (GeV/c)",
 308                100, 0., 10.,
 309                200, -50., 50.);
 310   AddHistogram(ID(kHistNsigmaTOFdSemi), "TOF N#sigma - d hypothesis;p (GeV/c)",
 311                100, 0., 10.,
 312                200, -50., 50.);
 313   AddHistogram(ID(kHistNsigmaTOFmismatchSemi), "TOF N#sigma - mismatch;p (GeV/c)",
 314                100, 0., 10.,
 315                200, -50., 50.);
 316   AddHistogram(ID(kHistNsigmaTOFmismatch2Semi), "TOF N#sigma - mismatch;p (GeV/c)",
 317                100, 0., 10.,
 318                200, -50., 50.);
 319
 320   // delta templates
 321   AddHistogram(ID(kHistDeltaTOFe), "TOF #Delta;p (GeV/c);t (ns)",
 322                100, 0., 10., 200, -2., 2.);
 323   AddHistogram(ID(kHistDeltaTOFmu), "TOF #Delta;p (GeV/c);t (ns)",
 324                100, 0., 10., 200, -2., 2.);
 325   AddHistogram(ID(kHistDeltaTOFpi), "TOF #Delta;p (GeV/c);t (ns)",
 326                100, 0., 10., 200, -2., 2.);
 327   AddHistogram(ID(kHistDeltaTOFk), "TOF #Delta;p (GeV/c);t (ns)",
 328                100, 0., 10., 200, -2., 2.);
 329   AddHistogram(ID(kHistDeltaTOFp), "TOF #Delta;p (GeV/c);t (ns)",
 330                100, 0., 10., 200, -2., 2.);
 331   AddHistogram(ID(kHistDeltaTOFd), "TOF #Delta;p (GeV/c);t (ns)",
 332                100, 0., 10., 200, -2., 2.);
 333
 334   AddHistogram(ID(kHistDeltaTOFeSemi), "TOF #Delta;p (GeV/c);t (ns)",
 335                100, 0., 10., 200, -2., 2.);
 336   AddHistogram(ID(kHistDeltaTOFmuSemi), "TOF #Delta;p (GeV/c);t (ns)",
 337                100, 0., 10., 200, -2., 2.);
 338   AddHistogram(ID(kHistDeltaTOFpiSemi), "TOF #Delta;p (GeV/c);t (ns)",
 339                100, 0., 10., 200, -2., 2.);
 340   AddHistogram(ID(kHistDeltaTOFkSemi), "TOF #Delta;p (GeV/c);t (ns)",
 341                100, 0., 10., 200, -2., 2.);
 342   AddHistogram(ID(kHistDeltaTOFpSemi), "TOF #Delta;p (GeV/c);t (ns)",
 343                100, 0., 10., 200, -2., 2.);
 344   AddHistogram(ID(kHistDeltaTOFdSemi), "TOF #Delta;p (GeV/c);t (ns)",
 345                100, 0., 10., 200, -2., 2.);
 346
 347   // sigma comparisons
 348   AddHistogram(ID(kHistExpSigmaTOFe), "TOF time of flight;p (GeV/c);t (ns)",
 349                100, 0., 10., 200, 0., .25);
 350   AddHistogram(ID(kHistExpSigmaTOFmu), "TOF time of flight;p (GeV/c);t (ns)",
 351                100, 0., 10., 200, 0., .25);
 352   AddHistogram(ID(kHistExpSigmaTOFpi), "TOF time of flight;p (GeV/c);t (ns)",
 353                100, 0., 10., 200, 0., .25);
 354   AddHistogram(ID(kHistExpSigmaTOFk), "TOF time of flight;p (GeV/c);t (ns)",
 355                100, 0., 10., 200, 0., .25);
 356   AddHistogram(ID(kHistExpSigmaTOFp), "TOF time of flight;p (GeV/c);t (ns)",
 357                100, 0., 10., 200, 0., .25);
 358   AddHistogram(ID(kHistExpSigmaTOFd), "TOF time of flight;p (GeV/c);t (ns)",
 359                100, 0., 10., 200, 0., .25);
 360
 361   AddHistogram(ID(kHistExpSigmaTOFeSemi), "TOF time of flight;p (GeV/c);t (ns)",
 362                100, 0., 10., 200, 0., .25);
 363   AddHistogram(ID(kHistExpSigmaTOFmuSemi), "TOF time of flight;p (GeV/c);t (ns)",
 364                100, 0., 10., 200, 0., .25);
 365   AddHistogram(ID(kHistExpSigmaTOFpiSemi), "TOF time of flight;p (GeV/c);t (ns)",
 366                100, 0., 10., 200, 0., .25);
 367   AddHistogram(ID(kHistExpSigmaTOFkSemi), "TOF time of flight;p (GeV/c);t (ns)",
 368                100, 0., 10., 200, 0., .25);
 369   AddHistogram(ID(kHistExpSigmaTOFpSemi), "TOF time of flight;p (GeV/c);t (ns)",
 370                100, 0., 10., 200, 0., .25);
 371   AddHistogram(ID(kHistExpSigmaTOFdSemi), "TOF time of flight;p (GeV/c);t (ns)",
 372                100, 0., 10., 200, 0., .25);
 373
 374   AddHistogram(ID(kHistCmpSigmaTOFe), "#sigma comparison;exp #sigma;template #sigma",
 375                200, 0., .25, 200, 0., .25);
 376   AddHistogram(ID(kHistCmpSigmaTOFmu), "#sigma comparison;exp #sigma;template #sigma",
 377                200, 0., .25, 200, 0., .25);
 378   AddHistogram(ID(kHistCmpSigmaTOFpi), "#sigma comparison;exp #sigma;template #sigma",
 379                200, 0., .25, 200, 0., .25);
 380   AddHistogram(ID(kHistCmpSigmaTOFk), "#sigma comparison;exp #sigma;template #sigma",
 381                200, 0., .25, 200, 0., .25);
 382   AddHistogram(ID(kHistCmpSigmaTOFp), "#sigma comparison;exp #sigma;template #sigma",
 383                200, 0., .25, 200, 0., .25);
 384   AddHistogram(ID(kHistCmpSigmaTOFd), "#sigma comparison;exp #sigma;template #sigma",
 385                200, 0., .25, 200, 0., .25);
 386
 387   AddHistogram(ID(kHistCmpSigmaTOFeSemi), "#sigma comparison;exp #sigma;template #sigma",
 388                200, 0., .25, 200, 0., .25);
 389   AddHistogram(ID(kHistCmpSigmaTOFmuSemi), "#sigma comparison;exp #sigma;template #sigma",
 390                200, 0., .25, 200, 0., .25);
 391   AddHistogram(ID(kHistCmpSigmaTOFpiSemi), "#sigma comparison;exp #sigma;template #sigma",
 392                200, 0., .25, 200, 0., .25);
 393   AddHistogram(ID(kHistCmpSigmaTOFkSemi), "#sigma comparison;exp #sigma;template #sigma",
 394                200, 0., .25, 200, 0., .25);
 395   AddHistogram(ID(kHistCmpSigmaTOFpSemi), "#sigma comparison;exp #sigma;template #sigma",
 396                200, 0., .25, 200, 0., .25);
 397   AddHistogram(ID(kHistCmpSigmaTOFdSemi), "#sigma comparison;exp #sigma;template #sigma",
 398                200, 0., .25, 200, 0., .25);
 399
 400   // Nsigma distributions
 401   AddHistogram(ID(kHistNsigmaTPCTOF), "N#sigma TPC-TOF;p (GeV/c);N#sigma_{TPC};N#sigma_{TOF}",
 402                100, 0., 10.,
 403                100, -25., 25.,
 404                200, -50., 50.);
 405   AddHistogram(ID(kHistNsigmaTPCTOFPt), "N#sigma TPC-TOF;p_{T} (GeV/c);N#sigma_{TPC};N#sigma_{TOF}",
 406                100, 0., 10.,
 407                100, -25., 25.,
 408                200, -50., 50.);
 409   AddHistogram(ID(kHistNsigmaTPCTOFUsed), "N#sigma TPC-TOF;p (GeV/c);N#sigma_{TPC};N#sigma_{TOF}",
 410                100, 0., 10.,
 411                100, -25., 25.,
 412                200, -50., 50.);
 413   AddHistogram(ID(kHistNsigmaTPCTOFUsedCentral), "N#sigma TPC-TOF (central);p (GeV/c);N#sigma_{TPC};N#sigma_{TOF}",
 414                100, 0., 10.,
 415                100, -25., 25.,
 416                200, -50., 50.);
 417   AddHistogram(ID(kHistNsigmaTPCTOFUsedSemiCentral), "N#sigma TPC-TOF (semi-central);p (GeV/c);N#sigma_{TPC};N#sigma_{TOF}",
 418                100, 0., 10.,
 419                100, -25., 25.,
 420                200, -50., 50.);
 421   AddHistogram(ID(kHistNsigmaTPCTOFUsedPt), "N#sigma TPC-TOF;p_{T} (GeV/c);N#sigma_{TPC};N#sigma_{TOF}",
 422                50, 0., 10.,
 423                100, -25., 25.,
 424                200, -50., 50.);
 425   AddHistogram(ID(kHistNsigmaTPCTOFUsedPtCentral), "N#sigma TPC-TOF;p_{T} (GeV/c);N#sigma_{TPC};N#sigma_{TOF}",
 426                50, 0., 10.,
 427                100, -25., 25.,
 428                200, -50., 50.);
 429   AddHistogram(ID(kHistNsigmaTPCTOFUsedPtSemiCentral), "N#sigma TPC-TOF;p_{T} (GeV/c);N#sigma_{TPC};N#sigma_{TOF}",
 430                50, 0., 10.,
 431                100, -25., 25.,
 432                200, -50., 50.);
 433
 434   AddHistogram(ID(kHistEvPlane), "default event plane;#Psi;counts",
 435                100, -0. * TMath::Pi(), 1. * TMath::Pi());
 436   AddHistogram(ID(kHistEvPlaneUsed), "default event plane;#Psi;counts",
 437                100, -0. * TMath::Pi(), 1. * TMath::Pi(),
 438                kClLast, -.5, kClLast-.5);
 439   AddHistogram(ID(kHistEvPlaneCheck), "backup event plane;#Psi;counts",
 440                100, -1. * TMath::Pi(), 1. * TMath::Pi());
 441   AddHistogram(ID(kHistEvPlaneCheckUsed), "backup event plane;#Psi;counts",
 442                100, -0. * TMath::Pi(), 1. * TMath::Pi(),
 443                kClLast, -.5, kClLast-.5);
 444   AddHistogram(ID(kHistEvPlaneCorr), "default - backup event plane;#Psi_{def};#Psi_{bak};counts",
 445                100, -0. * TMath::Pi(), 1. * TMath::Pi(),
 446                100, -0. * TMath::Pi(), 1. * TMath::Pi(),
 447                kClLast, -.5, kClLast-.5);
 448
 449   AddHistogram(ID(kHistJetPtCentral), "jet spectrum - central",
 450                40, 0., 200.);
 451   AddHistogram(ID(kHistJetPtSemi), "jet spectrum - semi-peripheral",
 452                40, 0., 200.);
 453
 454   AddHistogram(ID(kHistEtaPhiTrgHad), "trg had;#varphi;#eta",
 455                100, -0. * TMath::Pi(), 2. * TMath::Pi(),
 456                100, -2., 2.);
 457   AddHistogram(ID(kHistEtaPhiTrgJet), "trg jet;#varphi;#eta",
 458                100, -0. * TMath::Pi(), 2. * TMath::Pi(),
 459                100, -2., 2.);
 460   AddHistogram(ID(kHistEtaPhiAssHad), "ass had;#varphi;#eta",
 461                100, -0. * TMath::Pi(), 2. * TMath::Pi(),
 462                100, -2., 2.);
 463   AddHistogram(ID(kHistEtaPhiAssProt), "ass proton;#varphi;#eta",
 464                100, -0. * TMath::Pi(), 2. * TMath::Pi(),
 465                100, -2., 2.);
 466
 467   for (Int_t iCorr = 0; iCorr < kCorrLast; ++iCorr) {
 468     for (Int_t iCl = 0; iCl < kClLast; ++iCl) {
 469       for (Int_t iEv = 0; iEv < kEvLast; ++iEv) {
 470         GetHistCorr((CorrType_t) iCorr, (Class_t) iCl, (Ev_t) iEv) =
 471           new AliHistCorr(Form("corr_%s_%s_%s", fkCorrTypeName[iCorr], fkClassName[iCl], fkEvName[iEv]), fOutputList);
 472       }
 473     }
 474   }
 475
 476   PostData(1, fOutputList);
 477 }
 478
 479 Bool_t AliAnalysisTaskJetProtonCorr::Notify()
 480 {
 481   // actions to be taken upon notification about input file change
 482
 483   return AliAnalysisTaskSE::Notify();
 484 }
 485
 486 void AliAnalysisTaskJetProtonCorr::UserExec(Option_t * /* option */)
 487 {
 488   // actual work
 489
 490   // setup pointers to input data (null if unavailable)
 491   // mcEvent:  MC input
 492   // esdEvent: ESD input
 493   // outEvent: AOD output
 494   // aodEvent: AOD input if available, otherwise AOD output
 495
 496   fMCEvent   = this->MCEvent();
 497   fESDEvent  = dynamic_cast<AliESDEvent*>(this->InputEvent()); // could also be AOD input
 498   AliAODEvent* outEvent  = this->AODEvent();
 499   fAODEvent  = dynamic_cast<AliAODEvent*> (this->InputEvent());
 500   if (!fAODEvent)
 501     fAODEvent = outEvent;
 502
 503   if ((fDebug > 0) && fESDEvent)
 504     printf("event: %s-%06i\n", CurrentFileName(), fESDEvent->GetEventNumberInFile());
 505
 506   // record number of sampled events and detect trigger contributions
 507   FillH1(kHistStat, kStatSeen);
 508   if (!DetectTriggers()) {
 509     AliError("Failed to detect the triggers");
 510     return;
 511   }
 512
 513   if (!IsTrigger(kTriggerInt))
 514     return;
 515
 516   FillH1(kHistStat, kStatTrg);
 517
 518   // prepare the event
 519   // (make sure it is cleaned up in the end)
 520   if (PrepareEvent()) {
 521     FillH1(kHistStat, kStatUsed);
 522     FillH1(kHistCentrality, fCentrality);
 523     FillH1(kHistCentralityCheck, fCentralityCheck);
 524
 525     // event cuts
 526     if (TMath::Abs(fZvtx) > 10.)
 527       goto stop;
 528     if (GetCentrality() > 90.)
 529       goto stop;
 530
 531     FillH1(kHistStat, kStatEvCuts);
 532
 533     // event category
 534     DetectClasses();
 535     if (IsClass(kClCentral))
 536       FillH1(kHistStat, kStatCentral);
 537     if (IsClass(kClSemiCentral))
 538       FillH1(kHistStat, kStatSemiCentral);
 539
 540     FillH1(kHistEvPlane, fEventPlane);
 541     FillH1(kHistEvPlaneCheck, fEventPlaneCheck);
 542     for (Int_t iClass = 0; iClass < kClLast; ++iClass) {
 543       if (IsClass((Class_t) iClass)) {
 544         FillH2(kHistCentralityUsed, fCentrality, iClass);
 545         FillH2(kHistCentralityCheckUsed, fCentralityCheck, iClass);
 546         FillH2(kHistEvPlaneUsed, fEventPlane, iClass);
 547         FillH2(kHistEvPlaneCheckUsed, fEventPlaneCheck, iClass);
 548         FillH3(kHistEvPlaneCorr, fEventPlane, fEventPlaneCheck, iClass);
 549       }
 550     }
 551
 552     Bool_t corrEta  = fPIDResponse->UseTPCEtaCorrection();
 553     Bool_t corrMult = fPIDResponse->UseTPCMultiplicityCorrection();
 554
 555     // select trigger particles and potential associated particles/protons
 556     TObjArray trgArray[kTrgLast];
 557     TObjArray assArray[kAssLast];
 558
 559     TF1 fTOFsignal("fTOFsignal", &AliAnalysisTaskJetProtonCorr::TOFsignal, -2440., 2440., 4);
 560     fTOFsignal.SetParameter(0, 1.);
 561     fTOFsignal.SetParameter(1, 0.);
 562     fTOFsignal.SetParameter(2, 85.);
 563     fTOFsignal.SetParameter(3, 80.);
 564
 565     Int_t nPrimTracks = fPrimTrackArray ? fPrimTrackArray->GetEntries() : 0;
 566     for (Int_t iTrack = 0; iTrack < nPrimTracks; ++iTrack) {
 567       AliVTrack *trk = (AliVTrack*) fPrimTrackArray->At(iTrack);
 568       FillH3(kHistNsigmaTPCTOF,
 569              trk->P(),
 570              fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
 571              fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
 572       FillH3(kHistNsigmaTPCTOFPt,
 573              trk->Pt(),
 574              fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
 575              fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
 576
 577       if (AcceptTrigger(trk)) {
 578         trgArray[kTrgHad].Add(trk);
 579         FillH1(kHistEtaPhiTrgHad, trk->Phi(), trk->Eta());
 580       }
 581       if (AcceptAssoc(trk)) {
 582         assArray[kAssHad].Add(trk);
 583         FillH1(kHistEtaPhiAssHad, trk->Phi(), trk->Eta());
 584         FillH3(kHistNsigmaTPCTOFUsed,
 585                trk->P(),
 586                fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
 587                fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
 588         FillH3(kHistNsigmaTPCTOFUsedPt,
 589                trk->Pt(),
 590                fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
 591                fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
 592         if (IsClass(kClCentral)) {
 593           FillH3(kHistNsigmaTPCTOFUsedCentral,
 594                  trk->P(),
 595                  fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
 596                  fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
 597           FillH3(kHistNsigmaTPCTOFUsedPtCentral,
 598                  trk->Pt(),
 599                  fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
 600                  fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
 601         }
 602         if (IsClass(kClSemiCentral)) {
 603           FillH3(kHistNsigmaTPCTOFUsedSemiCentral,
 604                  trk->P(),
 605                  fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
 606                  fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
 607           FillH3(kHistNsigmaTPCTOFUsedPtSemiCentral,
 608                  trk->Pt(),
 609                  fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
 610                  fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
 611         }
 612         if (IsProton(trk)) {
 613           assArray[kAssProt].Add(trk);
 614           FillH1(kHistEtaPhiAssProt, trk->Phi(), trk->Eta());
 615         }
 616
 617         // template generation
 618         Double_t deltaTPC;
 619         fPIDResponse->GetSignalDelta(AliPIDResponse::kTPC, trk, AliPID::kProton, deltaTPC);
 620         FillH2(kHistSignalTPC, trk->P(), trk->GetTPCsignal());
 621         FillH2(kHistDeltaTPC, trk->P(), deltaTPC);
 622         if (fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC, trk) == AliPIDResponse::kDetPidOk) {
 623           Double_t expTPC = fPIDResponse->GetTPCResponse().GetExpectedSignal(trk, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, corrEta, corrMult);
 624           Double_t expSigmaTPC = fPIDResponse->GetTPCResponse().GetExpectedSigma(trk, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, corrEta, corrMult);
 625
 626           // loop over particles
 627           for (Int_t iParticle = 0; iParticle <= AliPID::kDeuteron; ++iParticle) {
 628             Double_t expTPCx = fPIDResponse->GetTPCResponse().GetExpectedSignal(trk, AliPID::EParticleType(iParticle), AliTPCPIDResponse::kdEdxDefault, corrEta, corrMult);
 629             Double_t expSigmaTPCx = fPIDResponse->GetTPCResponse().GetExpectedSigma(trk, AliPID::EParticleType(iParticle), AliTPCPIDResponse::kdEdxDefault, corrEta, corrMult);
 630             Double_t rndTPCx = gRandom->Gaus(expTPCx, expSigmaTPCx);
 631             if(IsClass(kClCentral))
 632               FillH2(kHistNsigmaTPCe, trk->P(), (rndTPCx-expTPC)/expSigmaTPC, 1., iParticle);
 633             if(IsClass(kClSemiCentral))
 634               FillH2(kHistNsigmaTPCeSemi, trk->P(), (rndTPCx-expTPC)/expSigmaTPC, 1., iParticle);
 635           }
 636         }
 637
 638         Double_t deltaTOF;
 639         fPIDResponse->GetSignalDelta(AliPIDResponse::kTOF, trk, AliPID::kProton, deltaTOF);
 640         FillH2(kHistSignalTOF, trk->Pt(), trk->GetTOFsignal() * 1.e-3); // ps -> ns
 641         if (fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF, trk) == AliPIDResponse::kDetPidOk) {
 642           AliTOFPIDResponse &tofResponse = fPIDResponse->GetTOFResponse();
 643
 644           Float_t p = trk->GetInnerParam()->GetP();
 645
 646           Double_t expTOF = fPIDResponse->GetTOFResponse().GetExpectedSignal(trk, AliPID::kProton);
 647           Double_t expSigmaTOF = fPIDResponse->GetTOFResponse().GetExpectedSigma(p, expTOF, AliPID::kProton);
 648           Double_t length = trk->GetIntegratedLength() * 1.e-2; // cm -> m
 649           Double_t tof = trk->GetTOFsignal() * 1.e-12; // ps -> s
 650           Double_t beta = length / tof / TMath::C();
 651
 652           FillH2(kHistBetaTOF, p, beta);
 653
 654           // intrinsic TOF smearing
 655           Double_t signalSigma = TMath::Sqrt(fTOFsignal.Variance(-2000., 2000.));
 656           Double_t signalSmear = fTOFsignal.GetRandom();
 657           // t0 smearing
 658           Float_t  timezeroSigma = tofResponse.GetT0binRes(tofResponse.GetMomBin(p));
 659           Double_t timezeroSmear  = gRandom->Gaus(0., timezeroSigma);
 660           // tracking smearing
 661           Double_t fPar[] = { 0.008, 0.008, 0.002, 40.0 };
 662
 663           // loop over particles
 664           for (Int_t iParticle = 0; iParticle <= AliPID::kDeuteron; ++iParticle) {
 665             Double_t expTOFx = fPIDResponse->GetTOFResponse().GetExpectedSignal(trk, AliPID::EParticleType(iParticle));
 666             Double_t cmpSigmaTOFx = fPIDResponse->GetTOFResponse().GetExpectedSigma(p, expTOFx, AliPID::EParticleType(iParticle));
 667
 668             // tracking smearing
 669             Double_t massx = AliPID::ParticleMassZ(AliPID::EParticleType(iParticle));
 670             Double_t dppx = fPar[0] + fPar[1] * p + fPar[2] * massx / p;
 671             Double_t expSigmaTOFx = dppx * expTOFx / (1.+ p * p / (massx * massx));
 672             Double_t texpSmearx = gRandom->Gaus(0., TMath::Sqrt(expSigmaTOFx * expSigmaTOFx + fPar[3]*fPar[3]/p/p));
 673             Double_t tmpSigmaTOFx = TMath::Sqrt(signalSigma*signalSigma +
 674                                                 timezeroSigma*timezeroSigma +
 675                                                 expSigmaTOFx*expSigmaTOFx + fPar[3]*fPar[3]/p/p);
 676             // printf("sigma comparison %i, %f: %f, %f, %f -> %f vs %f\n",
 677             //     iParticle, expTOFx,
 678             //     signalSigma, // signal
 679             //     timezeroSigma, // timezero
 680             //     expSigmaTOFx, // tracking
 681             //     tmpSigmaTOFx, // total
 682             //     cmpSigmaTOFx); // from PID response
 683
 684             // TOF signal
 685             Double_t rndTOFx = expTOFx + signalSmear + timezeroSmear + texpSmearx;
 686
 687             if (IsClass(kClCentral)) {
 688               FillH2(kHistNsigmaTOFe, trk->P(), (rndTOFx-expTOF)/expSigmaTOF, 1., iParticle);
 689               FillH2(kHistDeltaTOFe, trk->P(), (rndTOFx-expTOF) * 1.e-3, 1., iParticle);
 690               FillH2(kHistExpSigmaTOFe, p, cmpSigmaTOFx * 1.e-3, 1., iParticle);
 691               FillH2(kHistCmpSigmaTOFe, cmpSigmaTOFx * 1.e-3, tmpSigmaTOFx * 1.e-3, 1., iParticle);
 692             }
 693             if (IsClass(kClSemiCentral)) {
 694               FillH2(kHistNsigmaTOFeSemi, trk->P(), (rndTOFx-expTOF)/expSigmaTOF, 1., iParticle);
 695               FillH2(kHistDeltaTOFeSemi, trk->P(), (rndTOFx-expTOF) * 1.e-3, 1., iParticle);
 696               FillH2(kHistCmpSigmaTOFeSemi, cmpSigmaTOFx * 1.e-3, tmpSigmaTOFx * 1.e-3, 1., iParticle);
 697               FillH2(kHistExpSigmaTOFeSemi, p, cmpSigmaTOFx * 1.e-3, 1., iParticle);
 698             }
 699           }
 700
 701           Double_t rndTOFmismatch = AliTOFPIDResponse::GetMismatchRandomValue(trk->Eta());
 702
 703           if(IsClass(kClCentral)) {
 704             FillH2(kHistNsigmaTOFmismatch, p, (rndTOFmismatch - expTOF) / expSigmaTOF);
 705             FillH2(kHistNsigmaTOFmismatch2, p, (rndTOFmismatch - expTOF) / expSigmaTOF);
 706             FillH2(kHistDeltaTOF, trk->P(), deltaTOF * 1.e-3); // ps -> ns
 707           }
 708           if(IsClass(kClSemiCentral)) {
 709             FillH2(kHistNsigmaTOFmismatchSemi, p, (rndTOFmismatch - expTOF) / expSigmaTOF);
 710             FillH2(kHistNsigmaTOFmismatch2Semi, p, (rndTOFmismatch - expTOF) / expSigmaTOF);
 711             FillH2(kHistDeltaTOFSemi, trk->P(), deltaTOF * 1.e-3); // ps -> ns
 712           }
 713         }
 714       }
 715     }
 716
 717     // select trigger jet
 718     Int_t nJets = fJetArray ? fJetArray->GetEntries() : 0;
 719     for (Int_t iJet = 0; iJet < nJets; ++iJet) {
 720       AliAODJet *jet = (AliAODJet*) fJetArray->At(iJet);
 721       if (AcceptTrigger(jet)) {
 722         trgArray[kTrgJet].Add(jet);
 723         FillH1(kHistEtaPhiTrgJet, jet->Phi(), jet->Eta());
 724       }
 725     }
 726
 727     // correlate, both same and mixed event
 728     for (Int_t iClass = 0; iClass < kClLast; ++iClass) {
 729       if (IsClass((Class_t) iClass)) {
 730         for (Int_t iTrg = 0; iTrg < kTrgLast; ++iTrg) {
 731           for (Int_t iAss = 0; iAss < kAssLast; ++iAss) {
 732             // same event
 733             Correlate((Trg_t) iTrg, (Ass_t) iAss, (Class_t) iClass, kEvSame, &trgArray[iTrg], &assArray[iAss]);
 734
 735             // mixed event
 736             AliEventPool *pool = GetPool((Class_t) iClass, (Trg_t) iTrg, (Ass_t) iAss);
 737             if (pool && pool->IsReady()) {
 738               // printf("----- using pool: %i %i %i -----\n", iClass, iTrg, iAss);
 739               Int_t nEvents = pool->GetCurrentNEvents();
 740               for (Int_t iEvent = 0; iEvent < nEvents; ++iEvent) {
 741                 TObjArray *assTracks = pool->GetEvent(iEvent);
 742                 Correlate((Trg_t) iTrg, (Ass_t) iAss, (Class_t) iClass, kEvMix, &trgArray[iTrg], assTracks, 1./nEvents);
 743               }
 744             }
 745             // if (pool && !pool->IsReady()) {
 746             //   printf("----- pool not ready: %i %i %i -----\n", iClass, iTrg, iAss);
 747             //   pool->PrintInfo();
 748             // }
 749           }
 750
 751           // fill event pool for mixing
 752           // >= 0: don't require a trigger in the event
 753           // >= 1: require a trigger in the event
 754           if (trgArray[iTrg].GetEntries() >= 0) {
 755             for (Int_t iAss = 0; iAss < kAssLast; ++iAss) {
 756               AliEventPool *pool = GetPool((Class_t) iClass, (Trg_t) iTrg, (Ass_t) iAss);
 757               if (pool) {
 758                 pool->UpdatePool(CloneTracks(&assArray[iAss]));
 759                 // printf("----- updating pool: %i %i %i -----\n", iClass, iTrg, iAss);
 760                 // pool->PrintInfo();
 761               }
 762             }
 763           }
 764         }
 765       }
 766     }
 767   }
 768
 769  stop:
 770   CleanUpEvent();
 771
 772   PostData(1, fOutputList);
 773 }
 774
 775 void AliAnalysisTaskJetProtonCorr::Terminate(const Option_t * /* option */)
 776 {
 777   // actions at task termination
 778
 779 }
 780
 781 void AliAnalysisTaskJetProtonCorr::PrintTask(Option_t *option, Int_t indent) const
 782 {
 783   AliAnalysisTaskSE::PrintTask(option, indent);
 784
 785   std::cout << std::setw(indent) << " " << "using jet branch: " << fJetBranchName << std::endl;
 786 }
 787
 788 Bool_t AliAnalysisTaskJetProtonCorr::DetectTriggers()
 789 {
 790   fTriggerMask = 0;
 791
 792   AliVEvent::EOfflineTriggerTypes physSel = (AliVEvent::EOfflineTriggerTypes) ((AliInputEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 793   TString trgClasses = InputEvent()->GetFiredTriggerClasses();
 794
 795   // physics selection
 796   if (physSel & (AliVEvent::kAnyINT | AliVEvent::kCentral | AliVEvent::kSemiCentral))
 797     MarkTrigger(kTriggerInt);
 798
 799   return kTRUE;
 800 }
 801
 802 Bool_t AliAnalysisTaskJetProtonCorr::DetectClasses()
 803 {
 804   fClassMask = 0;
 805
 806   if (IsCentral())
 807     MarkClass(kClCentral);
 808
 809   if (IsSemiCentral())
 810     MarkClass(kClSemiCentral);
 811
 812   return kTRUE;
 813 }
 814
 815 Bool_t AliAnalysisTaskJetProtonCorr::PrepareEvent()
 816 {
 817   Bool_t eventGood = kTRUE;
 818
 819   // retrieve z-vertex position
 820   const AliVVertex *vtx = InputEvent()->GetPrimaryVertex();
 821   if (vtx && (vtx->GetNContributors() >= 3.))
 822     fZvtx = vtx->GetZ();
 823   else
 824     fZvtx = 100.;
 825
 826   // retrieve centrality
 827   AliCentrality *eventCentrality = InputEvent()->GetCentrality();
 828   if (eventCentrality) {
 829     fCentrality = eventCentrality->GetCentralityPercentile("V0M");
 830     fCentralityCheck = eventCentrality->GetCentralityPercentile("TRK");
 831     if (fCentrality >= 0.) {
 832       FillH1(kHistStat, kStatCent);
 833     } else {
 834       // centrality estimation not reliable
 835       eventGood = kFALSE;
 836       fCentrality = 105.;
 837     }
 838   }
 839   else
 840     eventGood = kFALSE;
 841
 842   // retrieve event plane
 843   AliEventplane *eventPlane = InputEvent()->GetEventplane();
 844   if (eventPlane) {
 845     FillH1(kHistStat, kStatEvPlane);
 846     fEventPlane = eventPlane->GetEventplane("Q");
 847     fEventPlaneCheck = eventPlane->GetEventplane("V0", InputEvent());
 848     if (fEventPlaneCheck < 0)
 849       fEventPlaneCheck += TMath::Pi();
 850   }
 851   else
 852     eventGood = kFALSE;
 853
 854   // retrieve PID
 855   fPIDResponse = ((AliInputEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();
 856   if (fPIDResponse)
 857     FillH1(kHistStat, kStatPID);
 858   else
 859     eventGood = kFALSE;
 860
 861   // retrieve primary tracks
 862   if (fESDEvent) {
 863     fPrimTrackArray = fCutsPrim->GetAcceptedTracks(fESDEvent);
 864   }
 865   else if (fAODEvent) {
 866     fPrimTrackArray = new TObjArray();
 867     Int_t nTracksAOD = fAODEvent->GetNumberOfTracks();
 868     for (Int_t iTrack = 0; iTrack < nTracksAOD; ++iTrack) {
 869       AliAODTrack *trk = fAODEvent->GetTrack(iTrack);
 870       // 4: track cuts esdTrackCutsH ???
 871       // 10: R_AA cuts
 872       if (trk->TestFilterMask(1 << 4))
 873         fPrimTrackArray->Add(trk);
 874     }
 875   }
 876   else
 877     eventGood = kFALSE;
 878
 879   // retrieve jet array
 880   if (fAODEvent) {
 881     fJetArray = dynamic_cast<TClonesArray*> (fAODEvent->FindListObject(fJetBranchName));
 882     if (!fJetArray) {
 883       printf("no jet branch \"%s\" found, in the AODs are:\n", fJetBranchName);
 884       if (fDebug > 0)
 885         fAODEvent->GetList()->Print();
 886     }
 887   }
 888
 889   // retrieve event pool for the event category
 890   if (eventGood) {
 891     for (Int_t iClass = 0; iClass < kClLast; ++iClass) {
 892       for (Int_t iTrg = 0; iTrg < kTrgLast; ++iTrg) {
 893         for (Int_t iAss = 0; iAss < kAssLast; ++iAss) {
 894           AliEventPoolManager *mgr = GetPoolMgr((Trg_t) iTrg, (Ass_t) iAss);
 895           GetPool((Class_t) iClass, (Trg_t) iTrg, (Ass_t) iAss) =
 896             // mgr ? mgr->GetEventPool(fCentrality, fZvtx, fEventPlane) : 0x0;
 897             mgr ? mgr->GetEventPool(fCentrality, fZvtx) : 0x0;
 898         }
 899       }
 900     }
 901   }
 902
 903   return eventGood;
 904 }
 905
 906 Bool_t AliAnalysisTaskJetProtonCorr::AcceptTrigger(AliVTrack *trg)
 907 {
 908   if ((trg->Pt() > fTrgPartPtMin) && (trg->Pt() < fTrgPartPtMax)) {
 909     if (IsCentral())
 910       return kTRUE;
 911     else if (IsSemiCentral()) {
 912       if (AcceptAngleToEvPlane(trg->Phi(), GetEventPlane())) {
 913         // printf("track accepted with phi = %f, psi = %f\n", trg->Phi(), GetEventPlane());
 914         return kTRUE;
 915       }
 916     }
 917   }
 918
 919   return kFALSE;
 920 }
 921
 922 AliVTrack* AliAnalysisTaskJetProtonCorr::GetLeadingTrack(AliAODJet *jet) const
 923 {
 924   // return leading track within a jet
 925
 926   // check contributing tracks
 927   Int_t nJetTracks = jet->GetRefTracks()->GetEntriesFast();
 928   Int_t iLeadingTrack = -1;
 929   Float_t ptLeadingTrack = 0.;
 930   for (Int_t iTrack=0; iTrack < nJetTracks; ++iTrack) {
 931     AliAODTrack *track = (AliAODTrack*) jet->GetRefTracks()->At(iTrack);
 932     // find the leading track
 933     if (track->Pt() > ptLeadingTrack) {
 934       ptLeadingTrack = track->Pt();
 935       iLeadingTrack = iTrack;
 936     }
 937   }
 938
 939   // retrieve the leading track
 940   return (AliAODTrack*) jet->GetRefTracks()->At(iLeadingTrack);
 941 }
 942
 943 Bool_t AliAnalysisTaskJetProtonCorr::AcceptTrigger(AliAODJet *trg)
 944 {
 945   // restrict eta
 946   if (TMath::Abs(trg->Eta()) > .5)
 947     return kFALSE;
 948
 949   // require hard leading track
 950   // (biased jet sample)
 951   if (GetLeadingTrack(trg)->Pt() < fTrgJetLeadTrkPtMin)
 952     return kFALSE;
 953
 954   // check for pt and azimuthal orientation
 955   if (IsSemiCentral()) {
 956     if (!AcceptAngleToEvPlane(trg->Phi(), GetEventPlane()))
 957       return kFALSE;
 958   }
 959
 960   // printf("jet accepted with phi = %f, psi = %f\n", trg->Phi(), GetEventPlane());
 961
 962   // spectrum for cross checks
 963   if (IsClass(kClCentral))
 964     FillH1(kHistJetPtCentral, trg->Pt());
 965   if (IsClass(kClSemiCentral))
 966     FillH1(kHistJetPtSemi, trg->Pt());
 967
 968   if ((trg->Pt() < fTrgJetPtMin) || (trg->Pt() > fTrgJetPtMax))
 969     return kFALSE;
 970
 971   return kTRUE;
 972 }
 973
 974 Bool_t AliAnalysisTaskJetProtonCorr::AcceptAssoc(AliVTrack *trk)
 975 {
 976   if ((trk->Pt() > fAssPartPtMin) && (trk->Pt() < fAssPartPtMax) &&
 977       (fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC, trk) == AliPIDResponse::kDetPidOk) &&
 978       (fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF, trk) == AliPIDResponse::kDetPidOk))
 979     if ((trk->GetTPCsignalN() >= 60.) &&
 980         (trk->GetTPCsignal() >= 10) &&
 981         (TMath::Abs(trk->Eta()) < .9))
 982       return kTRUE;
 983
 984   return kFALSE;
 985 }
 986
 987 Bool_t AliAnalysisTaskJetProtonCorr::IsProton(AliVTrack *trk)
 988 {
 989   Double_t nSigmaProtonTPC = fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton);
 990   Double_t nSigmaProtonTOF = fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton);
 991
 992   if ((TMath::Abs(nSigmaProtonTPC) <= 2.) && (TMath::Abs(nSigmaProtonTOF) <= 2.)) {
 993     return kTRUE;
 994   }
 995
 996   return kFALSE;
 997 }
 998
 999 Bool_t AliAnalysisTaskJetProtonCorr::AcceptAngleToEvPlane(Float_t phi, Float_t psi)
1000 {
1001   Float_t deltaPhi = phi - psi;
1002
1003   // map to interval [-pi/2, pi/2)
1004   deltaPhi = std::fmod(deltaPhi + TMath::Pi()/2., TMath::Pi()) - TMath::Pi()/2.;
1005   // printf("delta phi = %f with phi = %f, psi = %f\n", deltaPhi, phi, psi);
1006
1007   if (TMath::Abs(deltaPhi) < fTrgAngleToEvPlane)
1008     return kTRUE;
1009   else
1010     return kFALSE;
1011 }
1012
1013 Float_t AliAnalysisTaskJetProtonCorr::GetDPhiStar(Float_t phi1, Float_t pt1, Float_t charge1,
1014                                                   Float_t phi2, Float_t pt2, Float_t charge2,
1015                                                   Float_t radius, Float_t bSign)
1016 {
1017   // calculates dphistar
1018   // from AliUEHistograms
1019
1020   Float_t dphistar = phi1 - phi2
1021     - charge1 * bSign * TMath::ASin(0.075 * radius / pt1)
1022     + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
1023
1024   static const Double_t kPi = TMath::Pi();
1025
1026   // circularity
1027   if (dphistar > kPi)
1028     dphistar = kPi * 2 - dphistar;
1029   if (dphistar < -kPi)
1030     dphistar = -kPi * 2 - dphistar;
1031   if (dphistar > kPi) // might look funny but is needed
1032     dphistar = kPi * 2 - dphistar;
1033
1034   return dphistar;
1035 }
1036
1037
1038 Bool_t AliAnalysisTaskJetProtonCorr::AcceptTwoTracks(AliVParticle *trgPart, AliVParticle *assPart)
1039 {
1040   // apply two track pair cut
1041
1042   Float_t phi1 = trgPart->Phi();
1043   Float_t pt1 = trgPart->Pt();
1044   Float_t charge1 = trgPart->Charge();
1045
1046   Float_t phi2 = assPart->Phi();
1047   Float_t pt2 = assPart->Pt();
1048   Float_t charge2 = assPart->Charge();
1049
1050   // check ???
1051   Float_t bSign = 1.;
1052   Float_t deta = trgPart->Eta() - assPart->Eta();
1053
1054   // optimization
1055   if (TMath::Abs(deta) < fCutsTwoTrackEff * 2.5 * 3) {
1056     // check first boundaries to see if is worth to loop and find the minimum
1057     Float_t dphistar1 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 0.8, bSign);
1058     Float_t dphistar2 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 2.5, bSign);
1059
1060     const Float_t kLimit = fCutsTwoTrackEff * 3;
1061
1062     Float_t dphistarminabs = 1e5;
1063     // Float_t dphistarmin = 1e5;
1064     if ((TMath::Abs(dphistar1) < kLimit) ||
1065         (TMath::Abs(dphistar2) < kLimit) ||
1066         (dphistar1 * dphistar2 < 0)) {
1067       for (Double_t rad=0.8; rad<2.51; rad+=0.01) {
1068         Float_t dphistar = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, rad, bSign);
1069         Float_t dphistarabs = TMath::Abs(dphistar);
1070
1071         if (dphistarabs < dphistarminabs) {
1072           // dphistarmin = dphistar;
1073           dphistarminabs = dphistarabs;
1074         }
1075       }
1076
1077       if ((dphistarminabs < fCutsTwoTrackEff) &&
1078           (TMath::Abs(deta) < fCutsTwoTrackEff))
1079         return kFALSE;
1080     }
1081   }
1082
1083   return kTRUE;
1084 }
1085
1086 Bool_t AliAnalysisTaskJetProtonCorr::Correlate(CorrType_t corr, Class_t cl, Ev_t ev,
1087                                                TCollection *trgArray, TCollection *assArray, Float_t weight)
1088 {
1089   AliHistCorr *histCorr = GetHistCorr(corr, cl, ev);
1090
1091   TIter trgIter(trgArray);
1092
1093   while (AliVParticle *trgPart = (AliVParticle*) trgIter()) {
1094     // count the trigger
1095     histCorr->Trigger(weight);
1096
1097     // loop over associates
1098     TIter assIter(assArray);
1099     while (AliVParticle *assPart = (AliVParticle*) assIter()) {
1100       if (AcceptTwoTracks(trgPart, assPart))
1101         histCorr->Fill(trgPart, assPart, weight);
1102     }
1103   }
1104
1105   return kTRUE;
1106 }
1107
1108 Bool_t AliAnalysisTaskJetProtonCorr::Correlate(Trg_t trg, Ass_t ass, Class_t cl, Ev_t ev,
1109                                                TCollection *trgArray, TCollection *assArray, Float_t weight)
1110 {
1111   CorrType_t corr = (CorrType_t) (2 * trg + ass);
1112
1113   return Correlate(corr, cl, ev, trgArray, assArray, weight);
1114 }
1115
1116 Bool_t AliAnalysisTaskJetProtonCorr::CleanUpEvent()
1117 {
1118   if (fAODEvent) {
1119     delete fPrimTrackArray;
1120     fPrimTrackArray = 0x0;
1121   }
1122
1123   return kTRUE;
1124 }
1125
1126 TObjArray* AliAnalysisTaskJetProtonCorr::CloneTracks(TObjArray* tracks) const
1127 {
1128   TObjArray* tracksClone = new TObjArray;
1129   tracksClone->SetOwner(kTRUE);
1130
1131   Int_t nTracks = tracks->GetEntriesFast();
1132   for (Int_t i = 0; i < nTracks; i++) {
1133     // tracksClone->Add(new AliDPhiBasicParticle(particle->Eta(), particle->Phi(), particle->Pt(), particle->Charge()));
1134
1135     // WARNING: TObject::Clone() is very!!! expensive, unusable
1136     // tracksClone->Add(particle->Clone());
1137
1138     if (AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (tracks->At(i)))
1139       tracksClone->Add(new AliESDtrack(*esdTrack));
1140     else if (AliAODTrack* aodTrack = dynamic_cast<AliAODTrack*> (tracks->At(i)))
1141       tracksClone->Add(new AliAODTrack(*aodTrack));
1142   }
1143
1144   return tracksClone;
1145 }
1146
1147 AliAnalysisTaskJetProtonCorr::AliHistCorr::AliHistCorr(TString name, TList *outputList) :
1148   TNamed(name, name),
1149   fOutputList(outputList),
1150   fHistStat(0x0),
1151   fHistCorrPhi(0x0),
1152   fHistCorrPhi2(0x0),
1153   fHistCorrEtaPhi(0x0)
1154 {
1155   // ctor
1156
1157   fHistStat = new TH1F(Form("%s_stat", name.Data()), "statistics",
1158                        1, .5, 1.5);
1159
1160   fHistCorrPhi = new TH1F(Form("%s_phi", name.Data()), ";#Delta #phi",
1161                           100, -2.*TMath::Pi(), 2.*TMath::Pi());
1162   fHistCorrPhi2 = new TH2F(Form("%s_phi2", name.Data()), ";#phi_{trg};#phi_{ass}",
1163                           100, 0.*TMath::Pi(), 2.*TMath::Pi(),
1164                           100, 0.*TMath::Pi(), 2.*TMath::Pi());
1165   fHistCorrEtaPhi = new TH2F(Form("%s_etaphi", name.Data()), ";#Delta#phi;#Delta#eta",
1166                              100, -1., 2*TMath::Pi()-1.,
1167                              100, -2., 2.);
1168
1169   fOutputList->Add(fHistStat);
1170   fOutputList->Add(fHistCorrPhi);
1171   fOutputList->Add(fHistCorrPhi2);
1172   fOutputList->Add(fHistCorrEtaPhi);
1173 }
1174
1175 AliAnalysisTaskJetProtonCorr::AliHistCorr::~AliHistCorr()
1176 {
1177   // dtor
1178 }
1179
1180 void AliAnalysisTaskJetProtonCorr::AliHistCorr::Fill(AliVParticle *trgPart, AliVParticle *assPart, Float_t weight)
1181 {
1182   Float_t deltaEta = assPart->Eta() - trgPart->Eta();
1183   Float_t deltaPhi = assPart->Phi() - trgPart->Phi();
1184   if (deltaPhi > (2.*TMath::Pi()-1.))
1185     deltaPhi -= 2. * TMath::Pi();
1186   else if (deltaPhi < -1.)
1187     deltaPhi += 2. * TMath::Pi();
1188   // printf("trg: pt = %5.2f, phi = %5.2f, eta = %5.2f; ass: pt = %5.2f, phi = %5.2f, eta = %5.2f; deltaphi = %5.2f, deltaeta = %5.2f\n",
1189   //     trgPart->Pt(), trgPart->Phi(), trgPart->Eta(), assPart->Pt(), assPart->Phi(), assPart->Eta(), deltaPhi, deltaEta);
1190
1191   fHistCorrPhi->Fill(deltaPhi);
1192   fHistCorrPhi2->Fill(trgPart->Phi(), assPart->Phi(), weight);
1193   fHistCorrEtaPhi->Fill(deltaPhi, deltaEta, weight);
1194 }
1195
1196 // ----- histogram management -----
1197 TH1* AliAnalysisTaskJetProtonCorr::AddHistogram(Hist_t hist, const char *hid, TString title,
1198                                                  Int_t xbins, Float_t xmin, Float_t xmax,
1199                                                  Int_t binType)
1200 {
1201   TString hName;
1202   hName.Form("%s_%s", fShortTaskId, hid);
1203   hName.ToLower();
1204   TH1 *h = 0x0;
1205   if (binType == 0)
1206     h = new TH1I(hName.Data(), title,
1207                  xbins, xmin, xmax);
1208   else
1209     h = new TH1F(hName.Data(), title,
1210                  xbins, xmin, xmax);
1211   GetHistogram(hist) = h;
1212   fOutputList->Add(h);
1213   return h;
1214 }
1215
1216 TH2* AliAnalysisTaskJetProtonCorr::AddHistogram(Hist_t hist, const char *hid, TString title,
1217                                                  Int_t xbins, Float_t xmin, Float_t xmax,
1218                                                  Int_t ybins, Float_t ymin, Float_t ymax,
1219                                                  Int_t binType)
1220 {
1221   TString hName;
1222   hName.Form("%s_%s", fShortTaskId, hid);
1223   hName.ToLower();
1224   TH1 *h = 0x0;
1225   if (binType == 0)
1226     h = GetHistogram(hist) = new TH2I(hName.Data(), title,
1227                                      xbins, xmin, xmax,
1228                                      ybins, ymin, ymax);
1229   else
1230     h = GetHistogram(hist) = new TH2F(hName.Data(), title,
1231                                      xbins, xmin, xmax,
1232                                      ybins, ymin, ymax);
1233   fOutputList->Add(h);
1234   return (TH2*) h;
1235 }
1236
1237 TH3* AliAnalysisTaskJetProtonCorr::AddHistogram(Hist_t hist, const char *hid, TString title,
1238                                                  Int_t xbins, Float_t xmin, Float_t xmax,
1239                                                  Int_t ybins, Float_t ymin, Float_t ymax,
1240                                                  Int_t zbins, Float_t zmin, Float_t zmax,
1241                                                  Int_t binType)
1242 {
1243   TString hName;
1244   hName.Form("%s_%s", fShortTaskId, hid);
1245   hName.ToLower();
1246   TH1 *h = 0x0;
1247   if (binType == 0)
1248     h = GetHistogram(hist) = new TH3I(hName.Data(), title,
1249                                      xbins, xmin, xmax,
1250                                      ybins, ymin, ymax,
1251                                      zbins, zmin, zmax);
1252   else
1253     h = GetHistogram(hist) = new TH3F(hName.Data(), title,
1254                                      xbins, xmin, xmax,
1255                                      ybins, ymin, ymax,
1256                                      zbins, zmin, zmax);
1257   fOutputList->Add(h);
1258   return (TH3*) h;
1259 }