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