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