- update track cuts
authorkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 14 Oct 2013 12:18:20 +0000 (12:18 +0000)
committerkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 14 Oct 2013 12:18:20 +0000 (12:18 +0000)
- update event mixing classes
- additional PID histograms
- fix for trigger angle w.r.t. event plane
- add pair cuts

PWGJE/UserTasks/AliAnalysisTaskJetProtonCorr.cxx
PWGJE/UserTasks/AliAnalysisTaskJetProtonCorr.h

index f648868..03d2292 100644 (file)
@@ -5,6 +5,7 @@
 #include "TH2.h"
 #include "TH3.h"
 #include "TFormula.h"
+#include "TRandom.h"
 
 // analysis framework
 #include "AliAnalysisManager.h"
@@ -39,6 +40,9 @@
 
 #include "AliAnalysisTaskJetProtonCorr.h"
 
+#include <iostream>
+#include <cmath>
+
 AliAnalysisTaskJetProtonCorr::AliAnalysisTaskJetProtonCorr(const char *name) :
   AliAnalysisTaskSE(name),
   fMCEvent(0x0),
@@ -51,6 +55,7 @@ AliAnalysisTaskJetProtonCorr::AliAnalysisTaskJetProtonCorr(const char *name) :
   fZvtx(0.),
   fPIDResponse(0x0),
   fEventPlane(5.),
+  fEventPlaneCheck(5.),
   fPrimTrackArray(0x0),
   fJetArray(0x0),
   fPoolMgr(),
@@ -59,13 +64,17 @@ AliAnalysisTaskJetProtonCorr::AliAnalysisTaskJetProtonCorr(const char *name) :
   fOutputList(),
   fHist(),
   fShortTaskId("jet_prot_corr"),
+  fUseStandardCuts(kTRUE),
   fCutsPrim(0x0),
+  fCutsTwoTrackEff(0.02),
   fTrgPartPtMin(6.),
   fTrgPartPtMax(8.),
   fTrgJetPtMin(50.),
   fTrgJetPtMax(80.),
+  fTrgJetLeadTrkPtMin(5.),
   fAssPartPtMin(2.),
-  fAssPartPtMax(4.)
+  fAssPartPtMax(4.),
+  fTrgAngleToEvPlane(TMath::Pi() / 4.)
 {
   // default ctor
 
@@ -82,39 +91,49 @@ AliAnalysisTaskJetProtonCorr::AliAnalysisTaskJetProtonCorr(const char *name) :
   fkEvName[kEvMix]  = "mixed";
 
   // track cuts
-  fCutsPrim = new AliESDtrackCuts();
-
-  // this is taken from PWGJE track cuts
-  TFormula *f1NClustersTPCLinearPtDep = new TFormula("f1NClustersTPCLinearPtDep","70.+30./20.*x");
-  fCutsPrim->SetMinNClustersTPCPtDep(f1NClustersTPCLinearPtDep,20.);
-  fCutsPrim->SetMinNClustersTPC(70);
-  fCutsPrim->SetMaxChi2PerClusterTPC(4);
-  fCutsPrim->SetRequireTPCStandAlone(kTRUE); //cut on NClustersTPC and chi2TPC Iter1
-  fCutsPrim->SetAcceptKinkDaughters(kFALSE);
-  fCutsPrim->SetRequireTPCRefit(kTRUE);
-  fCutsPrim->SetMaxFractionSharedTPCClusters(0.4);
-  // ITS
-  fCutsPrim->SetRequireITSRefit(kTRUE);
-  //accept secondaries
-  fCutsPrim->SetMaxDCAToVertexXY(2.4);
-  fCutsPrim->SetMaxDCAToVertexZ(3.2);
-  fCutsPrim->SetDCAToVertex2D(kTRUE);
-  //reject fakes
-  fCutsPrim->SetMaxChi2PerClusterITS(36);
-  fCutsPrim->SetMaxChi2TPCConstrainedGlobal(36);
-
-  fCutsPrim->SetRequireSigmaToVertex(kFALSE);
-
-  fCutsPrim->SetEtaRange(-0.9,0.9);
-  fCutsPrim->SetPtRange(0.15, 1E+15);
+  if (fUseStandardCuts) {
+    fCutsPrim = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE);
+  } else {
+    fCutsPrim = new AliESDtrackCuts();
+
+    // this is taken from PWGJE track cuts
+    TFormula *f1NClustersTPCLinearPtDep = new TFormula("f1NClustersTPCLinearPtDep","70.+30./20.*x");
+    fCutsPrim->SetMinNClustersTPCPtDep(f1NClustersTPCLinearPtDep,20.);
+    fCutsPrim->SetMinNClustersTPC(70);
+    fCutsPrim->SetMaxChi2PerClusterTPC(4);
+    fCutsPrim->SetRequireTPCStandAlone(kTRUE); //cut on NClustersTPC and chi2TPC Iter1
+    fCutsPrim->SetAcceptKinkDaughters(kFALSE);
+    fCutsPrim->SetRequireTPCRefit(kTRUE);
+    fCutsPrim->SetMaxFractionSharedTPCClusters(0.4);
+    // ITS
+    fCutsPrim->SetRequireITSRefit(kTRUE);
+    //accept secondaries
+    fCutsPrim->SetMaxDCAToVertexXY(2.4);
+    fCutsPrim->SetMaxDCAToVertexZ(3.2);
+    fCutsPrim->SetDCAToVertex2D(kTRUE);
+    //reject fakes
+    fCutsPrim->SetMaxChi2PerClusterITS(36);
+    fCutsPrim->SetMaxChi2TPCConstrainedGlobal(36);
+
+    fCutsPrim->SetRequireSigmaToVertex(kFALSE);
+
+    fCutsPrim->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
+  }
 
-  fCutsPrim->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
+  fCutsPrim->SetEtaRange(-0.9, 0.9);
+  fCutsPrim->SetPtRange(0.15, 1E+15);
 
   // event mixing pool
-  Double_t centralityBins[] = { 0., 10., 40., 80.};
+  Double_t centralityBins[] = {
+    0., 2., 4., 6., 8., 10., // central
+    30., 32., 34., 36., 38., 40., 42., 44., 46., 48., 50., // semi-central
+    90.
+  };
   Int_t nCentralityBins = sizeof(centralityBins)/sizeof(centralityBins[0]);
 
-  Double_t vertexBins[] = { -10., -5., 0., 5., 10.};
+  Double_t vertexBins[] = {
+    -10., -8., -6., -4., -2., 0., 2., 4., 6., 8., 10.
+  };
   Int_t nVertexBins = sizeof(vertexBins)/sizeof(vertexBins[0]);
 
   Double_t psiBins[12];
@@ -127,8 +146,8 @@ AliAnalysisTaskJetProtonCorr::AliAnalysisTaskJetProtonCorr(const char *name) :
       GetPoolMgr((Trg_t) iTrg, (Ass_t) iAss) =
        new AliEventPoolManager(10, 100,
                                nCentralityBins, centralityBins,
-                               nVertexBins, vertexBins,
-                               nPsiBins, psiBins);
+                               nVertexBins, vertexBins);
+                               // nPsiBins, psiBins);
       GetPoolMgr((Trg_t) iTrg, (Ass_t) iAss)->SetTargetValues(100, .1, 1);
     }
   }
@@ -165,6 +184,8 @@ void AliAnalysisTaskJetProtonCorr::UserCreateOutputObjects()
   histStat->GetXaxis()->SetBinLabel(ID(kStatCent));
   histStat->GetXaxis()->SetBinLabel(ID(kStatEvPlane));
   histStat->GetXaxis()->SetBinLabel(ID(kStatPID));
+  histStat->GetXaxis()->SetBinLabel(ID(kStatCentral));
+  histStat->GetXaxis()->SetBinLabel(ID(kStatSemiCentral));
 
   AddHistogram(ID(kHistCentrality), "centrality;C;counts",
               110, -5., 105.);
@@ -183,18 +204,251 @@ void AliAnalysisTaskJetProtonCorr::UserCreateOutputObjects()
   hist->GetYaxis()->SetBinLabel(LAB(kClSemiCentral));
   hist->GetYaxis()->SetBinLabel(LAB(kClDijet));
 
-  AddHistogram(ID(kHistNsigmaTPCTOF), "N#sigma TPC-TOF;p_{T} (GeV/c);N#sigma_{TPC};N#sigma_{TOF}",
-               100, 0., 10.,
-               100, -5., 5.,
-               100, -5., 5.);
 
-  AddHistogram(ID(kHistEvPlane), "event plane angle;#Psi;counts",
+  AddHistogram(ID(kHistSignalTPC), "TPC dE/dx;p (GeV/c);dE/dx (arb. units)",
+              100, 0., 10., 200, 0., 300.);
+  AddHistogram(ID(kHistSignalTOF), "TOF time of flight;p_{T} (GeV/c);t (ns)",
+              100, 0., 10., 200, 0., 50.);
+  AddHistogram(ID(kHistBetaTOF), "TOF beta;p (GeV/c); #beta",
+              100, 0., 10.,
+              100, 0., 1.);
+  AddHistogram(ID(kHistDeltaTPC), "TPC dE/dx;p (GeV/c);dE/dx (arb. units)",
+              100, 0., 10., 200, -100., 100.);
+  AddHistogram(ID(kHistDeltaTPCSemi), "TPC dE/dx;p (GeV/c);dE/dx (arb. units)",
+              100, 0., 10., 200, -100., 100.);
+  AddHistogram(ID(kHistDeltaTOF), "TOF time of flight;p_{T} (GeV/c);t (ns)",
+              100, 0., 10., 200, -2., 2.);
+  AddHistogram(ID(kHistDeltaTOFSemi), "TOF time of flight;p_{T} (GeV/c);t (ns)",
+              100, 0., 10., 200, -2., 2.);
+
+  // Nsigma templates
+  AddHistogram(ID(kHistNsigmaTPCe), "TPC N#sigma - e hypothesis;p (GeV/c)",
+              100, 0., 10.,
+              100, -25., 25.);
+  AddHistogram(ID(kHistNsigmaTPCmu), "TPC N#sigma - #mu hypothesis;p (GeV/c)",
+              100, 0., 10.,
+              100, -25., 25.);
+  AddHistogram(ID(kHistNsigmaTPCpi), "TPC N#sigma - #pi hypothesis;p (GeV/c)",
+              100, 0., 10.,
+              100, -25., 25.);
+  AddHistogram(ID(kHistNsigmaTPCk), "TPC N#sigma - K hypothesis;p (GeV/c)",
+              100, 0., 10.,
+              100, -25., 25.);
+  AddHistogram(ID(kHistNsigmaTPCp), "TPC N#sigma - p hypothesis;p (GeV/c)",
+              100, 0., 10.,
+              100, -25., 25.);
+  AddHistogram(ID(kHistNsigmaTPCd), "TPC N#sigma - d hypothesis;p (GeV/c)",
+              100, 0., 10.,
+              100, -25., 25.);
+  AddHistogram(ID(kHistNsigmaTPCe_e), "TPC N#sigma - e hypothesis (id. e);p (GeV/c)",
+              100, 0., 10.,
+              100, -25., 25.);
+
+  AddHistogram(ID(kHistNsigmaTOFe), "TOF N#sigma - e hypothesis;p (GeV/c)",
+              100, 0., 10.,
+              200, -50., 50.);
+  AddHistogram(ID(kHistNsigmaTOFmu), "TOF N#sigma - #mu hypothesis;p (GeV/c)",
+              100, 0., 10.,
+              200, -50., 50.);
+  AddHistogram(ID(kHistNsigmaTOFpi), "TOF N#sigma - #pi hypothesis;p (GeV/c)",
+              100, 0., 10.,
+              200, -50., 50.);
+  AddHistogram(ID(kHistNsigmaTOFk), "TOF N#sigma - K hypothesis;p (GeV/c)",
+              100, 0., 10.,
+              200, -50., 50.);
+  AddHistogram(ID(kHistNsigmaTOFp), "TOF N#sigma - p hypothesis;p (GeV/c)",
+              100, 0., 10.,
+              200, -50., 50.);
+  AddHistogram(ID(kHistNsigmaTOFd), "TOF N#sigma - d hypothesis;p (GeV/c)",
+              100, 0., 10.,
+              200, -50., 50.);
+  AddHistogram(ID(kHistNsigmaTOFmismatch), "TOF N#sigma - mismatch;p (GeV/c)",
+              100, 0., 10.,
+              200, -50., 50.);
+  AddHistogram(ID(kHistNsigmaTOFmismatch2), "TOF N#sigma - mismatch;p (GeV/c)",
+              100, 0., 10.,
+              200, -50., 50.);
+
+  // Nsigma templates
+  AddHistogram(ID(kHistNsigmaTPCeSemi), "TPC N#sigma - e hypothesis;p (GeV/c)",
+              100, 0., 10.,
+              100, -25., 25.);
+  AddHistogram(ID(kHistNsigmaTPCmuSemi), "TPC N#sigma - #mu hypothesis;p (GeV/c)",
+              100, 0., 10.,
+              100, -25., 25.);
+  AddHistogram(ID(kHistNsigmaTPCpiSemi), "TPC N#sigma - #pi hypothesis;p (GeV/c)",
+              100, 0., 10.,
+              100, -25., 25.);
+  AddHistogram(ID(kHistNsigmaTPCkSemi), "TPC N#sigma - K hypothesis;p (GeV/c)",
+              100, 0., 10.,
+              100, -25., 25.);
+  AddHistogram(ID(kHistNsigmaTPCpSemi), "TPC N#sigma - p hypothesis;p (GeV/c)",
+              100, 0., 10.,
+              100, -25., 25.);
+  AddHistogram(ID(kHistNsigmaTPCdSemi), "TPC N#sigma - d hypothesis;p (GeV/c)",
+              100, 0., 10.,
+              100, -25., 25.);
+  AddHistogram(ID(kHistNsigmaTPCe_eSemi), "TPC N#sigma - e hypothesis (id. e);p (GeV/c)",
+              100, 0., 10.,
+              100, -25., 25.);
+
+  AddHistogram(ID(kHistNsigmaTOFeSemi), "TOF N#sigma - e hypothesis;p (GeV/c)",
+              100, 0., 10.,
+              200, -50., 50.);
+  AddHistogram(ID(kHistNsigmaTOFmuSemi), "TOF N#sigma - #mu hypothesis;p (GeV/c)",
+              100, 0., 10.,
+              200, -50., 50.);
+  AddHistogram(ID(kHistNsigmaTOFpiSemi), "TOF N#sigma - #pi hypothesis;p (GeV/c)",
+              100, 0., 10.,
+              200, -50., 50.);
+  AddHistogram(ID(kHistNsigmaTOFkSemi), "TOF N#sigma - K hypothesis;p (GeV/c)",
+              100, 0., 10.,
+              200, -50., 50.);
+  AddHistogram(ID(kHistNsigmaTOFpSemi), "TOF N#sigma - p hypothesis;p (GeV/c)",
+              100, 0., 10.,
+              200, -50., 50.);
+  AddHistogram(ID(kHistNsigmaTOFdSemi), "TOF N#sigma - d hypothesis;p (GeV/c)",
+              100, 0., 10.,
+              200, -50., 50.);
+  AddHistogram(ID(kHistNsigmaTOFmismatchSemi), "TOF N#sigma - mismatch;p (GeV/c)",
+              100, 0., 10.,
+              200, -50., 50.);
+  AddHistogram(ID(kHistNsigmaTOFmismatch2Semi), "TOF N#sigma - mismatch;p (GeV/c)",
+              100, 0., 10.,
+              200, -50., 50.);
+
+  // delta templates
+  AddHistogram(ID(kHistDeltaTOFe), "TOF #Delta;p (GeV/c);t (ns)",
+              100, 0., 10., 200, -2., 2.);
+  AddHistogram(ID(kHistDeltaTOFmu), "TOF #Delta;p (GeV/c);t (ns)",
+              100, 0., 10., 200, -2., 2.);
+  AddHistogram(ID(kHistDeltaTOFpi), "TOF #Delta;p (GeV/c);t (ns)",
+              100, 0., 10., 200, -2., 2.);
+  AddHistogram(ID(kHistDeltaTOFk), "TOF #Delta;p (GeV/c);t (ns)",
+              100, 0., 10., 200, -2., 2.);
+  AddHistogram(ID(kHistDeltaTOFp), "TOF #Delta;p (GeV/c);t (ns)",
+              100, 0., 10., 200, -2., 2.);
+  AddHistogram(ID(kHistDeltaTOFd), "TOF #Delta;p (GeV/c);t (ns)",
+              100, 0., 10., 200, -2., 2.);
+
+  AddHistogram(ID(kHistDeltaTOFeSemi), "TOF #Delta;p (GeV/c);t (ns)",
+              100, 0., 10., 200, -2., 2.);
+  AddHistogram(ID(kHistDeltaTOFmuSemi), "TOF #Delta;p (GeV/c);t (ns)",
+              100, 0., 10., 200, -2., 2.);
+  AddHistogram(ID(kHistDeltaTOFpiSemi), "TOF #Delta;p (GeV/c);t (ns)",
+              100, 0., 10., 200, -2., 2.);
+  AddHistogram(ID(kHistDeltaTOFkSemi), "TOF #Delta;p (GeV/c);t (ns)",
+              100, 0., 10., 200, -2., 2.);
+  AddHistogram(ID(kHistDeltaTOFpSemi), "TOF #Delta;p (GeV/c);t (ns)",
+              100, 0., 10., 200, -2., 2.);
+  AddHistogram(ID(kHistDeltaTOFdSemi), "TOF #Delta;p (GeV/c);t (ns)",
+              100, 0., 10., 200, -2., 2.);
+
+  // sigma comparisons
+  AddHistogram(ID(kHistExpSigmaTOFe), "TOF time of flight;p (GeV/c);t (ns)",
+              100, 0., 10., 200, 0., .25);
+  AddHistogram(ID(kHistExpSigmaTOFmu), "TOF time of flight;p (GeV/c);t (ns)",
+              100, 0., 10., 200, 0., .25);
+  AddHistogram(ID(kHistExpSigmaTOFpi), "TOF time of flight;p (GeV/c);t (ns)",
+              100, 0., 10., 200, 0., .25);
+  AddHistogram(ID(kHistExpSigmaTOFk), "TOF time of flight;p (GeV/c);t (ns)",
+              100, 0., 10., 200, 0., .25);
+  AddHistogram(ID(kHistExpSigmaTOFp), "TOF time of flight;p (GeV/c);t (ns)",
+              100, 0., 10., 200, 0., .25);
+  AddHistogram(ID(kHistExpSigmaTOFd), "TOF time of flight;p (GeV/c);t (ns)",
+              100, 0., 10., 200, 0., .25);
+
+  AddHistogram(ID(kHistExpSigmaTOFeSemi), "TOF time of flight;p (GeV/c);t (ns)",
+              100, 0., 10., 200, 0., .25);
+  AddHistogram(ID(kHistExpSigmaTOFmuSemi), "TOF time of flight;p (GeV/c);t (ns)",
+              100, 0., 10., 200, 0., .25);
+  AddHistogram(ID(kHistExpSigmaTOFpiSemi), "TOF time of flight;p (GeV/c);t (ns)",
+              100, 0., 10., 200, 0., .25);
+  AddHistogram(ID(kHistExpSigmaTOFkSemi), "TOF time of flight;p (GeV/c);t (ns)",
+              100, 0., 10., 200, 0., .25);
+  AddHistogram(ID(kHistExpSigmaTOFpSemi), "TOF time of flight;p (GeV/c);t (ns)",
+              100, 0., 10., 200, 0., .25);
+  AddHistogram(ID(kHistExpSigmaTOFdSemi), "TOF time of flight;p (GeV/c);t (ns)",
+              100, 0., 10., 200, 0., .25);
+
+  AddHistogram(ID(kHistCmpSigmaTOFe), "#sigma comparison;exp #sigma;template #sigma",
+              200, 0., .25, 200, 0., .25);
+  AddHistogram(ID(kHistCmpSigmaTOFmu), "#sigma comparison;exp #sigma;template #sigma",
+              200, 0., .25, 200, 0., .25);
+  AddHistogram(ID(kHistCmpSigmaTOFpi), "#sigma comparison;exp #sigma;template #sigma",
+              200, 0., .25, 200, 0., .25);
+  AddHistogram(ID(kHistCmpSigmaTOFk), "#sigma comparison;exp #sigma;template #sigma",
+              200, 0., .25, 200, 0., .25);
+  AddHistogram(ID(kHistCmpSigmaTOFp), "#sigma comparison;exp #sigma;template #sigma",
+              200, 0., .25, 200, 0., .25);
+  AddHistogram(ID(kHistCmpSigmaTOFd), "#sigma comparison;exp #sigma;template #sigma",
+              200, 0., .25, 200, 0., .25);
+
+  AddHistogram(ID(kHistCmpSigmaTOFeSemi), "#sigma comparison;exp #sigma;template #sigma",
+              200, 0., .25, 200, 0., .25);
+  AddHistogram(ID(kHistCmpSigmaTOFmuSemi), "#sigma comparison;exp #sigma;template #sigma",
+              200, 0., .25, 200, 0., .25);
+  AddHistogram(ID(kHistCmpSigmaTOFpiSemi), "#sigma comparison;exp #sigma;template #sigma",
+              200, 0., .25, 200, 0., .25);
+  AddHistogram(ID(kHistCmpSigmaTOFkSemi), "#sigma comparison;exp #sigma;template #sigma",
+              200, 0., .25, 200, 0., .25);
+  AddHistogram(ID(kHistCmpSigmaTOFpSemi), "#sigma comparison;exp #sigma;template #sigma",
+              200, 0., .25, 200, 0., .25);
+  AddHistogram(ID(kHistCmpSigmaTOFdSemi), "#sigma comparison;exp #sigma;template #sigma",
+              200, 0., .25, 200, 0., .25);
+
+  // Nsigma distributions
+  AddHistogram(ID(kHistNsigmaTPCTOF), "N#sigma TPC-TOF;p (GeV/c);N#sigma_{TPC};N#sigma_{TOF}",
+               100, 0., 10.,
+               100, -25., 25.,
+               200, -50., 50.);
+  AddHistogram(ID(kHistNsigmaTPCTOFPt), "N#sigma TPC-TOF;p_{T} (GeV/c);N#sigma_{TPC};N#sigma_{TOF}",
+               100, 0., 10.,
+               100, -25., 25.,
+               200, -50., 50.);
+  AddHistogram(ID(kHistNsigmaTPCTOFUsed), "N#sigma TPC-TOF;p (GeV/c);N#sigma_{TPC};N#sigma_{TOF}",
+               100, 0., 10.,
+               100, -25., 25.,
+               200, -50., 50.);
+  AddHistogram(ID(kHistNsigmaTPCTOFUsedCentral), "N#sigma TPC-TOF (central);p (GeV/c);N#sigma_{TPC};N#sigma_{TOF}",
+               100, 0., 10.,
+               100, -25., 25.,
+               200, -50., 50.);
+  AddHistogram(ID(kHistNsigmaTPCTOFUsedSemiCentral), "N#sigma TPC-TOF (semi-central);p (GeV/c);N#sigma_{TPC};N#sigma_{TOF}",
+               100, 0., 10.,
+               100, -25., 25.,
+               200, -50., 50.);
+  AddHistogram(ID(kHistNsigmaTPCTOFUsedPt), "N#sigma TPC-TOF;p_{T} (GeV/c);N#sigma_{TPC};N#sigma_{TOF}",
+               50, 0., 10.,
+               100, -25., 25.,
+               200, -50., 50.);
+  AddHistogram(ID(kHistNsigmaTPCTOFUsedPtCentral), "N#sigma TPC-TOF;p_{T} (GeV/c);N#sigma_{TPC};N#sigma_{TOF}",
+               50, 0., 10.,
+               100, -25., 25.,
+               200, -50., 50.);
+  AddHistogram(ID(kHistNsigmaTPCTOFUsedPtSemiCentral), "N#sigma TPC-TOF;p_{T} (GeV/c);N#sigma_{TPC};N#sigma_{TOF}",
+               50, 0., 10.,
+               100, -25., 25.,
+               200, -50., 50.);
+
+  AddHistogram(ID(kHistEvPlane), "default event plane;#Psi;counts",
                100, -0. * TMath::Pi(), 1. * TMath::Pi());
-  AddHistogram(ID(kHistEvPlaneUsed), "event plane angle;#Psi;counts",
+  AddHistogram(ID(kHistEvPlaneUsed), "default event plane;#Psi;counts",
+               100, -0. * TMath::Pi(), 1. * TMath::Pi(),
+              kClLast, -.5, kClLast-.5);
+  AddHistogram(ID(kHistEvPlaneCheck), "backup event plane;#Psi;counts",
+               100, -1. * TMath::Pi(), 1. * TMath::Pi());
+  AddHistogram(ID(kHistEvPlaneCheckUsed), "backup event plane;#Psi;counts",
+               100, -0. * TMath::Pi(), 1. * TMath::Pi(),
+              kClLast, -.5, kClLast-.5);
+  AddHistogram(ID(kHistEvPlaneCorr), "default - backup event plane;#Psi_{def};#Psi_{bak};counts",
+               100, -0. * TMath::Pi(), 1. * TMath::Pi(),
                100, -0. * TMath::Pi(), 1. * TMath::Pi(),
               kClLast, -.5, kClLast-.5);
 
-  AddHistogram(ID(kHistJetPt), "jet spectrum",
+  AddHistogram(ID(kHistJetPtCentral), "jet spectrum - central",
+               40, 0., 200.);
+  AddHistogram(ID(kHistJetPtSemi), "jet spectrum - semi-peripheral",
                40, 0., 200.);
 
   AddHistogram(ID(kHistEtaPhiTrgHad), "trg had;#varphi;#eta",
@@ -278,27 +532,48 @@ void AliAnalysisTaskJetProtonCorr::UserExec(Option_t * /* option */)
 
     // event category
     DetectClasses();
+    if (IsClass(kClCentral))
+      FillH1(kHistStat, kStatCentral);
+    if (IsClass(kClSemiCentral))
+      FillH1(kHistStat, kStatSemiCentral);
 
     FillH1(kHistEvPlane, fEventPlane);
+    FillH1(kHistEvPlaneCheck, fEventPlaneCheck);
     for (Int_t iClass = 0; iClass < kClLast; ++iClass) {
       if (IsClass((Class_t) iClass)) {
         FillH2(kHistCentralityUsed, fCentrality, iClass);
         FillH2(kHistCentralityCheckUsed, fCentralityCheck, iClass);
        FillH2(kHistEvPlaneUsed, fEventPlane, iClass);
+       FillH2(kHistEvPlaneCheckUsed, fEventPlaneCheck, iClass);
+       FillH3(kHistEvPlaneCorr, fEventPlane, fEventPlaneCheck, iClass);
       }
     }
 
+    Bool_t corrEta  = fPIDResponse->UseTPCEtaCorrection();
+    Bool_t corrMult = fPIDResponse->UseTPCMultiplicityCorrection();
+
     // select trigger particles and potential associated particles/protons
     TObjArray trgArray[kTrgLast];
     TObjArray assArray[kAssLast];
 
+    TF1 fTOFsignal("fTOFsignal", &AliAnalysisTaskJetProtonCorr::TOFsignal, -2440., 2440., 4);
+    fTOFsignal.SetParameter(0, 1.);
+    fTOFsignal.SetParameter(1, 0.);
+    fTOFsignal.SetParameter(2, 85.);
+    fTOFsignal.SetParameter(3, 80.);
+
     Int_t nPrimTracks = fPrimTrackArray ? fPrimTrackArray->GetEntries() : 0;
     for (Int_t iTrack = 0; iTrack < nPrimTracks; ++iTrack) {
       AliVTrack *trk = (AliVTrack*) fPrimTrackArray->At(iTrack);
       FillH3(kHistNsigmaTPCTOF,
+             trk->P(),
+             fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
+             fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
+      FillH3(kHistNsigmaTPCTOFPt,
              trk->Pt(),
              fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
              fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
+
       if (AcceptTrigger(trk)) {
        trgArray[kTrgHad].Add(trk);
        FillH1(kHistEtaPhiTrgHad, trk->Phi(), trk->Eta());
@@ -306,10 +581,136 @@ void AliAnalysisTaskJetProtonCorr::UserExec(Option_t * /* option */)
       if (AcceptAssoc(trk)) {
        assArray[kAssHad].Add(trk);
        FillH1(kHistEtaPhiAssHad, trk->Phi(), trk->Eta());
+       FillH3(kHistNsigmaTPCTOFUsed,
+              trk->P(),
+              fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
+              fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
+       FillH3(kHistNsigmaTPCTOFUsedPt,
+              trk->Pt(),
+              fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
+              fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
+       if (IsClass(kClCentral)) {
+         FillH3(kHistNsigmaTPCTOFUsedCentral,
+                trk->P(),
+                fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
+                fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
+         FillH3(kHistNsigmaTPCTOFUsedPtCentral,
+                trk->Pt(),
+                fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
+                fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
+       }
+       if (IsClass(kClSemiCentral)) {
+         FillH3(kHistNsigmaTPCTOFUsedSemiCentral,
+                trk->P(),
+                fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
+                fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
+         FillH3(kHistNsigmaTPCTOFUsedPtSemiCentral,
+                trk->Pt(),
+                fPIDResponse->NumberOfSigmasTPC(trk, AliPID::kProton),
+                fPIDResponse->NumberOfSigmasTOF(trk, AliPID::kProton));
+       }
        if (IsProton(trk)) {
          assArray[kAssProt].Add(trk);
          FillH1(kHistEtaPhiAssProt, trk->Phi(), trk->Eta());
        }
+
+       // template generation
+       Double_t deltaTPC;
+       fPIDResponse->GetSignalDelta(AliPIDResponse::kTPC, trk, AliPID::kProton, deltaTPC);
+       FillH2(kHistSignalTPC, trk->P(), trk->GetTPCsignal());
+       FillH2(kHistDeltaTPC, trk->P(), deltaTPC);
+       if (fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC, trk) == AliPIDResponse::kDetPidOk) {
+         Double_t expTPC = fPIDResponse->GetTPCResponse().GetExpectedSignal(trk, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, corrEta, corrMult);
+         Double_t expSigmaTPC = fPIDResponse->GetTPCResponse().GetExpectedSigma(trk, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, corrEta, corrMult);
+
+         // loop over particles
+         for (Int_t iParticle = 0; iParticle <= AliPID::kDeuteron; ++iParticle) {
+           Double_t expTPCx = fPIDResponse->GetTPCResponse().GetExpectedSignal(trk, AliPID::EParticleType(iParticle), AliTPCPIDResponse::kdEdxDefault, corrEta, corrMult);
+           Double_t expSigmaTPCx = fPIDResponse->GetTPCResponse().GetExpectedSigma(trk, AliPID::EParticleType(iParticle), AliTPCPIDResponse::kdEdxDefault, corrEta, corrMult);
+           Double_t rndTPCx = gRandom->Gaus(expTPCx, expSigmaTPCx);
+           if(IsClass(kClCentral))
+             FillH2(kHistNsigmaTPCe, trk->P(), (rndTPCx-expTPC)/expSigmaTPC, 1., iParticle);
+           if(IsClass(kClSemiCentral))
+             FillH2(kHistNsigmaTPCeSemi, trk->P(), (rndTPCx-expTPC)/expSigmaTPC, 1., iParticle);
+         }
+       }
+
+       Double_t deltaTOF;
+       fPIDResponse->GetSignalDelta(AliPIDResponse::kTOF, trk, AliPID::kProton, deltaTOF);
+       FillH2(kHistSignalTOF, trk->Pt(), trk->GetTOFsignal() * 1.e-3); // ps -> ns
+       if (fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF, trk) == AliPIDResponse::kDetPidOk) {
+         AliTOFPIDResponse &tofResponse = fPIDResponse->GetTOFResponse();
+
+         Float_t p = trk->GetInnerParam()->GetP();
+
+         Double_t expTOF = fPIDResponse->GetTOFResponse().GetExpectedSignal(trk, AliPID::kProton);
+         Double_t expSigmaTOF = fPIDResponse->GetTOFResponse().GetExpectedSigma(p, expTOF, AliPID::kProton);
+         Double_t length = trk->GetIntegratedLength() * 1.e-2; // cm -> m
+         Double_t tof = trk->GetTOFsignal() * 1.e-12; // ps -> s
+         Double_t beta = length / tof / TMath::C();
+
+         FillH2(kHistBetaTOF, p, beta);
+
+         // intrinsic TOF smearing
+         Double_t signalSigma = TMath::Sqrt(fTOFsignal.Variance(-2000., 2000.));
+         Double_t signalSmear = fTOFsignal.GetRandom();
+         // t0 smearing
+         Float_t  timezeroSigma = tofResponse.GetT0binRes(tofResponse.GetMomBin(p));
+         Double_t timezeroSmear  = gRandom->Gaus(0., timezeroSigma);
+         // tracking smearing
+         Double_t fPar[] = { 0.008, 0.008, 0.002, 40.0 };
+
+         // loop over particles
+         for (Int_t iParticle = 0; iParticle <= AliPID::kDeuteron; ++iParticle) {
+           Double_t expTOFx = fPIDResponse->GetTOFResponse().GetExpectedSignal(trk, AliPID::EParticleType(iParticle));
+           Double_t cmpSigmaTOFx = fPIDResponse->GetTOFResponse().GetExpectedSigma(p, expTOFx, AliPID::EParticleType(iParticle));
+
+           // tracking smearing
+           Double_t massx = AliPID::ParticleMassZ(AliPID::EParticleType(iParticle));
+           Double_t dppx = fPar[0] + fPar[1] * p + fPar[2] * massx / p;
+           Double_t expSigmaTOFx = dppx * expTOFx / (1.+ p * p / (massx * massx));
+           Double_t texpSmearx = gRandom->Gaus(0., TMath::Sqrt(expSigmaTOFx * expSigmaTOFx + fPar[3]*fPar[3]/p/p));
+           Double_t tmpSigmaTOFx = TMath::Sqrt(signalSigma*signalSigma +
+                                               timezeroSigma*timezeroSigma +
+                                               expSigmaTOFx*expSigmaTOFx + fPar[3]*fPar[3]/p/p);
+           // printf("sigma comparison %i, %f: %f, %f, %f -> %f vs %f\n",
+           //     iParticle, expTOFx,
+           //     signalSigma, // signal
+           //     timezeroSigma, // timezero
+           //     expSigmaTOFx, // tracking
+           //     tmpSigmaTOFx, // total
+           //     cmpSigmaTOFx); // from PID response
+
+           // TOF signal
+           Double_t rndTOFx = expTOFx + signalSmear + timezeroSmear + texpSmearx;
+
+           if (IsClass(kClCentral)) {
+             FillH2(kHistNsigmaTOFe, trk->P(), (rndTOFx-expTOF)/expSigmaTOF, 1., iParticle);
+             FillH2(kHistDeltaTOFe, trk->P(), (rndTOFx-expTOF) * 1.e-3, 1., iParticle);
+             FillH2(kHistExpSigmaTOFe, p, cmpSigmaTOFx * 1.e-3, 1., iParticle);
+             FillH2(kHistCmpSigmaTOFe, cmpSigmaTOFx * 1.e-3, tmpSigmaTOFx * 1.e-3, 1., iParticle);
+           }
+           if (IsClass(kClSemiCentral)) {
+             FillH2(kHistNsigmaTOFeSemi, trk->P(), (rndTOFx-expTOF)/expSigmaTOF, 1., iParticle);
+             FillH2(kHistDeltaTOFeSemi, trk->P(), (rndTOFx-expTOF) * 1.e-3, 1., iParticle);
+             FillH2(kHistCmpSigmaTOFeSemi, cmpSigmaTOFx * 1.e-3, tmpSigmaTOFx * 1.e-3, 1., iParticle);
+             FillH2(kHistExpSigmaTOFeSemi, p, cmpSigmaTOFx * 1.e-3, 1., iParticle);
+           }
+         }
+
+         Double_t rndTOFmismatch = AliTOFPIDResponse::GetMismatchRandomValue(trk->Eta());
+
+         if(IsClass(kClCentral)) {
+           FillH2(kHistNsigmaTOFmismatch, p, (rndTOFmismatch - expTOF) / expSigmaTOF);
+           FillH2(kHistNsigmaTOFmismatch2, p, (rndTOFmismatch - expTOF) / expSigmaTOF);
+           FillH2(kHistDeltaTOF, trk->P(), deltaTOF * 1.e-3); // ps -> ns
+         }
+         if(IsClass(kClSemiCentral)) {
+           FillH2(kHistNsigmaTOFmismatchSemi, p, (rndTOFmismatch - expTOF) / expSigmaTOF);
+           FillH2(kHistNsigmaTOFmismatch2Semi, p, (rndTOFmismatch - expTOF) / expSigmaTOF);
+           FillH2(kHistDeltaTOFSemi, trk->P(), deltaTOF * 1.e-3); // ps -> ns
+         }
+       }
       }
     }
 
@@ -317,7 +718,6 @@ void AliAnalysisTaskJetProtonCorr::UserExec(Option_t * /* option */)
     Int_t nJets = fJetArray ? fJetArray->GetEntries() : 0;
     for (Int_t iJet = 0; iJet < nJets; ++iJet) {
       AliAODJet *jet = (AliAODJet*) fJetArray->At(iJet);
-      FillH1(kHistJetPt, jet->Pt());
       if (AcceptTrigger(jet)) {
        trgArray[kTrgJet].Add(jet);
        FillH1(kHistEtaPhiTrgJet, jet->Phi(), jet->Eta());
@@ -349,7 +749,9 @@ void AliAnalysisTaskJetProtonCorr::UserExec(Option_t * /* option */)
          }
 
          // fill event pool for mixing
-         if (trgArray[iTrg].GetEntries() > 0) {
+         // >= 0: don't require a trigger in the event
+         // >= 1: require a trigger in the event
+         if (trgArray[iTrg].GetEntries() >= 0) {
            for (Int_t iAss = 0; iAss < kAssLast; ++iAss) {
              AliEventPool *pool = GetPool((Class_t) iClass, (Trg_t) iTrg, (Ass_t) iAss);
              if (pool) {
@@ -376,6 +778,13 @@ void AliAnalysisTaskJetProtonCorr::Terminate(const Option_t * /* option */)
 
 }
 
+void AliAnalysisTaskJetProtonCorr::PrintTask(Option_t *option, Int_t indent) const
+{
+  AliAnalysisTaskSE::PrintTask(option, indent);
+
+  std::cout << std::setw(indent) << " " << "using jet branch: " << fJetBranchName << std::endl;
+}
+
 Bool_t AliAnalysisTaskJetProtonCorr::DetectTriggers()
 {
   fTriggerMask = 0;
@@ -427,12 +836,17 @@ Bool_t AliAnalysisTaskJetProtonCorr::PrepareEvent()
       fCentrality = 105.;
     }
   }
+  else
+    eventGood = kFALSE;
 
   // retrieve event plane
   AliEventplane *eventPlane = InputEvent()->GetEventplane();
   if (eventPlane) {
     FillH1(kHistStat, kStatEvPlane);
     fEventPlane = eventPlane->GetEventplane("Q");
+    fEventPlaneCheck = eventPlane->GetEventplane("V0", InputEvent());
+    if (fEventPlaneCheck < 0)
+      fEventPlaneCheck += TMath::Pi();
   }
   else
     eventGood = kFALSE;
@@ -453,7 +867,8 @@ Bool_t AliAnalysisTaskJetProtonCorr::PrepareEvent()
     Int_t nTracksAOD = fAODEvent->GetNumberOfTracks();
     for (Int_t iTrack = 0; iTrack < nTracksAOD; ++iTrack) {
       AliAODTrack *trk = fAODEvent->GetTrack(iTrack);
-      // track cuts esdTrackCutsH ???
+      // 4: track cuts esdTrackCutsH ???
+      // 10: R_AA cuts
       if (trk->TestFilterMask(1 << 4))
         fPrimTrackArray->Add(trk);
     }
@@ -478,7 +893,8 @@ Bool_t AliAnalysisTaskJetProtonCorr::PrepareEvent()
        for (Int_t iAss = 0; iAss < kAssLast; ++iAss) {
          AliEventPoolManager *mgr = GetPoolMgr((Trg_t) iTrg, (Ass_t) iAss);
          GetPool((Class_t) iClass, (Trg_t) iTrg, (Ass_t) iAss) =
-           mgr ? mgr->GetEventPool(fCentrality, fZvtx, fEventPlane) : 0x0;
+           // mgr ? mgr->GetEventPool(fCentrality, fZvtx, fEventPlane) : 0x0;
+           mgr ? mgr->GetEventPool(fCentrality, fZvtx) : 0x0;
        }
       }
     }
@@ -493,44 +909,77 @@ Bool_t AliAnalysisTaskJetProtonCorr::AcceptTrigger(AliVTrack *trg)
     if (IsCentral())
       return kTRUE;
     else if (IsSemiCentral()) {
-      // check calculation ???
-      Float_t deltaPhi = trg->Phi() - GetEventPlane();
-      if (deltaPhi > TMath::Pi())
-        deltaPhi -= 2. * TMath::Pi();
-      else if (deltaPhi < -TMath::Pi())
-        deltaPhi += 2. * TMath::Pi();
-      if (TMath::Abs(TMath::Pi() - TMath::Abs(deltaPhi)) >= (3./8. * TMath::Pi()))
+      if (AcceptAngleToEvPlane(trg->Phi(), GetEventPlane())) {
+       // printf("track accepted with phi = %f, psi = %f\n", trg->Phi(), GetEventPlane());
         return kTRUE;
+      }
     }
   }
 
   return kFALSE;
 }
 
-Bool_t AliAnalysisTaskJetProtonCorr::AcceptTrigger(AliAODJet *trg)
+AliVTrack* AliAnalysisTaskJetProtonCorr::GetLeadingTrack(AliAODJet *jet) const
 {
-  if ((trg->Pt() > fTrgJetPtMin) && (trg->Pt() < fTrgJetPtMax)) {
-    if (IsCentral())
-      return kTRUE;
-    else if (IsSemiCentral()) {
-      // check calculation ???
-      Float_t deltaPhi = trg->Phi() - GetEventPlane();
-      if (deltaPhi > TMath::Pi())
-        deltaPhi -= 2. * TMath::Pi();
-      else if (deltaPhi < -TMath::Pi())
-        deltaPhi += 2. * TMath::Pi();
-      if (TMath::Abs(TMath::Pi() - TMath::Abs(deltaPhi)) >= (3./8. * TMath::Pi()))
-        return kTRUE;
+  // return leading track within a jet
+
+  // check contributing tracks
+  Int_t nJetTracks = jet->GetRefTracks()->GetEntriesFast();
+  Int_t iLeadingTrack = -1;
+  Float_t ptLeadingTrack = 0.;
+  for (Int_t iTrack=0; iTrack < nJetTracks; ++iTrack) {
+    AliAODTrack *track = (AliAODTrack*) jet->GetRefTracks()->At(iTrack);
+    // find the leading track
+    if (track->Pt() > ptLeadingTrack) {
+      ptLeadingTrack = track->Pt();
+      iLeadingTrack = iTrack;
     }
   }
 
-  return kFALSE;
+  // retrieve the leading track
+  return (AliAODTrack*) jet->GetRefTracks()->At(iLeadingTrack);
 }
 
-Bool_t AliAnalysisTaskJetProtonCorr::AcceptAssoc(AliVTrack *trg)
+Bool_t AliAnalysisTaskJetProtonCorr::AcceptTrigger(AliAODJet *trg)
 {
-  if ((trg->Pt() > fAssPartPtMin) && (trg->Pt() < fAssPartPtMax))
-    return kTRUE;
+  // restrict eta
+  if (TMath::Abs(trg->Eta()) > .5)
+    return kFALSE;
+
+  // require hard leading track
+  // (biased jet sample)
+  if (GetLeadingTrack(trg)->Pt() < fTrgJetLeadTrkPtMin)
+    return kFALSE;
+
+  // check for pt and azimuthal orientation
+  if (IsSemiCentral()) {
+    if (!AcceptAngleToEvPlane(trg->Phi(), GetEventPlane()))
+      return kFALSE;
+  }
+
+  // printf("jet accepted with phi = %f, psi = %f\n", trg->Phi(), GetEventPlane());
+
+  // spectrum for cross checks
+  if (IsClass(kClCentral))
+    FillH1(kHistJetPtCentral, trg->Pt());
+  if (IsClass(kClSemiCentral))
+    FillH1(kHistJetPtSemi, trg->Pt());
+
+  if ((trg->Pt() < fTrgJetPtMin) || (trg->Pt() > fTrgJetPtMax))
+    return kFALSE;
+
+  return kTRUE;
+}
+
+Bool_t AliAnalysisTaskJetProtonCorr::AcceptAssoc(AliVTrack *trk)
+{
+  if ((trk->Pt() > fAssPartPtMin) && (trk->Pt() < fAssPartPtMax) &&
+      (fPIDResponse->CheckPIDStatus(AliPIDResponse::kTPC, trk) == AliPIDResponse::kDetPidOk) &&
+      (fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF, trk) == AliPIDResponse::kDetPidOk))
+    if ((trk->GetTPCsignalN() >= 60.) &&
+       (trk->GetTPCsignal() >= 10) &&
+       (TMath::Abs(trk->Eta()) < .9))
+      return kTRUE;
 
   return kFALSE;
 }
@@ -547,31 +996,96 @@ Bool_t AliAnalysisTaskJetProtonCorr::IsProton(AliVTrack *trk)
   return kFALSE;
 }
 
-Bool_t AliAnalysisTaskJetProtonCorr::Correlate(CorrType_t corr, Class_t cl, Ev_t ev,
-                                              TCollection *trgArray, TCollection *assArray, Float_t weight)
+Bool_t AliAnalysisTaskJetProtonCorr::AcceptAngleToEvPlane(Float_t phi, Float_t psi)
 {
-  AliHistCorr *histCorr = GetHistCorr(corr, cl, ev);
+  Float_t deltaPhi = phi - psi;
 
-  TIter trgIter(trgArray);
+  // map to interval [-pi/2, pi/2)
+  deltaPhi = std::fmod(deltaPhi + TMath::Pi()/2., TMath::Pi()) - TMath::Pi()/2.;
+  // printf("delta phi = %f with phi = %f, psi = %f\n", deltaPhi, phi, psi);
 
-  while (AliVParticle *trgPart = (AliVParticle*) trgIter()) {
-    // count the trigger
-    histCorr->Trigger(weight);
+  if (TMath::Abs(deltaPhi) < fTrgAngleToEvPlane)
+    return kTRUE;
+  else
+    return kFALSE;
+}
 
-    // loop over associates
-    TIter assIter(assArray);
-    while (AliVParticle *assPart = (AliVParticle*) assIter()) {
-      histCorr->Fill(trgPart, assPart, weight);
+Float_t AliAnalysisTaskJetProtonCorr::GetDPhiStar(Float_t phi1, Float_t pt1, Float_t charge1,
+                                                 Float_t phi2, Float_t pt2, Float_t charge2,
+                                                 Float_t radius, Float_t bSign)
+{
+  // calculates dphistar
+  // from AliUEHistograms
+
+  Float_t dphistar = phi1 - phi2
+    - charge1 * bSign * TMath::ASin(0.075 * radius / pt1)
+    + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
+
+  static const Double_t kPi = TMath::Pi();
+
+  // circularity
+  if (dphistar > kPi)
+    dphistar = kPi * 2 - dphistar;
+  if (dphistar < -kPi)
+    dphistar = -kPi * 2 - dphistar;
+  if (dphistar > kPi) // might look funny but is needed
+    dphistar = kPi * 2 - dphistar;
+
+  return dphistar;
+}
+
+
+Bool_t AliAnalysisTaskJetProtonCorr::AcceptTwoTracks(AliVParticle *trgPart, AliVParticle *assPart)
+{
+  // apply two track pair cut
+
+  Float_t phi1 = trgPart->Phi();
+  Float_t pt1 = trgPart->Pt();
+  Float_t charge1 = trgPart->Charge();
+
+  Float_t phi2 = assPart->Phi();
+  Float_t pt2 = assPart->Pt();
+  Float_t charge2 = assPart->Charge();
+
+  // check ???
+  Float_t bSign = 1.;
+  Float_t deta = trgPart->Eta() - assPart->Eta();
+
+  // optimization
+  if (TMath::Abs(deta) < fCutsTwoTrackEff * 2.5 * 3) {
+    // check first boundaries to see if is worth to loop and find the minimum
+    Float_t dphistar1 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 0.8, bSign);
+    Float_t dphistar2 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 2.5, bSign);
+
+    const Float_t kLimit = fCutsTwoTrackEff * 3;
+
+    Float_t dphistarminabs = 1e5;
+    // Float_t dphistarmin = 1e5;
+    if ((TMath::Abs(dphistar1) < kLimit) ||
+       (TMath::Abs(dphistar2) < kLimit) ||
+       (dphistar1 * dphistar2 < 0)) {
+      for (Double_t rad=0.8; rad<2.51; rad+=0.01) {
+       Float_t dphistar = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, rad, bSign);
+       Float_t dphistarabs = TMath::Abs(dphistar);
+
+       if (dphistarabs < dphistarminabs) {
+         // dphistarmin = dphistar;
+         dphistarminabs = dphistarabs;
+       }
+      }
+
+      if ((dphistarminabs < fCutsTwoTrackEff) &&
+         (TMath::Abs(deta) < fCutsTwoTrackEff))
+       return kFALSE;
     }
   }
 
   return kTRUE;
 }
 
-Bool_t AliAnalysisTaskJetProtonCorr::Correlate(Trg_t trg, Ass_t ass, Class_t cl, Ev_t ev,
+Bool_t AliAnalysisTaskJetProtonCorr::Correlate(CorrType_t corr, Class_t cl, Ev_t ev,
                                               TCollection *trgArray, TCollection *assArray, Float_t weight)
 {
-  CorrType_t corr = (CorrType_t) (2 * trg + ass);
   AliHistCorr *histCorr = GetHistCorr(corr, cl, ev);
 
   TIter trgIter(trgArray);
@@ -583,13 +1097,22 @@ Bool_t AliAnalysisTaskJetProtonCorr::Correlate(Trg_t trg, Ass_t ass, Class_t cl,
     // loop over associates
     TIter assIter(assArray);
     while (AliVParticle *assPart = (AliVParticle*) assIter()) {
-      histCorr->Fill(trgPart, assPart, weight);
+      if (AcceptTwoTracks(trgPart, assPart))
+       histCorr->Fill(trgPart, assPart, weight);
     }
   }
 
   return kTRUE;
 }
 
+Bool_t AliAnalysisTaskJetProtonCorr::Correlate(Trg_t trg, Ass_t ass, Class_t cl, Ev_t ev,
+                                              TCollection *trgArray, TCollection *assArray, Float_t weight)
+{
+  CorrType_t corr = (CorrType_t) (2 * trg + ass);
+
+  return Correlate(corr, cl, ev, trgArray, assArray, weight);
+}
+
 Bool_t AliAnalysisTaskJetProtonCorr::CleanUpEvent()
 {
   if (fAODEvent) {
index b10889d..02d24c7 100644 (file)
@@ -32,7 +32,7 @@ public:
   virtual void   Terminate(const Option_t *option);
 
   // task configuration
-  void SetJetBranchName(const char* const branchName) { strncpy(fJetBranchName, branchName, fgkStringLength-1); }
+  void SetJetBranchName(const char* const branchName) { strncpy(fJetBranchName, branchName, fgkStringLength); }
   const char* GetJetBranchName() const { return fJetBranchName; }
 
   void SetPtThrPart(Float_t minPt, Float_t maxPt) { fTrgPartPtMin = minPt; fTrgPartPtMax = maxPt; }
@@ -45,6 +45,21 @@ public:
   Float_t GetPtMinAss() const { return fAssPartPtMin; }
   Float_t GetPtMaxAss() const { return fAssPartPtMax; }
 
+  void PrintTask(Option_t *option, Int_t indent) const;
+
+  static Double_t TOFsignal(Double_t *x, Double_t *par)
+  {
+    Double_t norm = par[0];
+    Double_t mean = par[1];
+    Double_t sigma = par[2];
+    Double_t tail = par[3];
+
+    if (x[0] <= (tail + mean))
+      return norm * TMath::Gaus(x[0], mean, sigma);
+    else
+      return norm * TMath::Gaus(tail + mean, mean, sigma) * TMath::Exp(-tail * (x[0] - tail - mean) / (sigma * sigma));
+  }
+
   // histograms
   enum Hist_t {
       kHistStat = 0,
@@ -52,10 +67,94 @@ public:
       kHistCentralityUsed,
       kHistCentralityCheck,
       kHistCentralityCheckUsed,
+      kHistSignalTPC,
+      kHistSignalTOF,
+      kHistBetaTOF,
+      kHistDeltaTPC,
+      kHistDeltaTPCSemi,
+      kHistDeltaTOF,
+      kHistDeltaTOFSemi,
+      kHistExpSigmaTOFe,
+      kHistExpSigmaTOFmu,
+      kHistExpSigmaTOFpi,
+      kHistExpSigmaTOFk,
+      kHistExpSigmaTOFp,
+      kHistExpSigmaTOFd,
+      kHistExpSigmaTOFeSemi,
+      kHistExpSigmaTOFmuSemi,
+      kHistExpSigmaTOFpiSemi,
+      kHistExpSigmaTOFkSemi,
+      kHistExpSigmaTOFpSemi,
+      kHistExpSigmaTOFdSemi,
+      kHistCmpSigmaTOFe,
+      kHistCmpSigmaTOFmu,
+      kHistCmpSigmaTOFpi,
+      kHistCmpSigmaTOFk,
+      kHistCmpSigmaTOFp,
+      kHistCmpSigmaTOFd,
+      kHistCmpSigmaTOFeSemi,
+      kHistCmpSigmaTOFmuSemi,
+      kHistCmpSigmaTOFpiSemi,
+      kHistCmpSigmaTOFkSemi,
+      kHistCmpSigmaTOFpSemi,
+      kHistCmpSigmaTOFdSemi,
+      kHistNsigmaTPCe,
+      kHistNsigmaTPCmu,
+      kHistNsigmaTPCpi,
+      kHistNsigmaTPCk,
+      kHistNsigmaTPCp,
+      kHistNsigmaTPCd,
+      kHistNsigmaTPCe_e,
+      kHistNsigmaTOFe,
+      kHistNsigmaTOFmu,
+      kHistNsigmaTOFpi,
+      kHistNsigmaTOFk,
+      kHistNsigmaTOFp,
+      kHistNsigmaTOFd,
+      kHistNsigmaTOFmismatch,
+      kHistNsigmaTOFmismatch2,
+      kHistDeltaTOFe,
+      kHistDeltaTOFmu,
+      kHistDeltaTOFpi,
+      kHistDeltaTOFk,
+      kHistDeltaTOFp,
+      kHistDeltaTOFd,
+      kHistNsigmaTPCeSemi,
+      kHistNsigmaTPCmuSemi,
+      kHistNsigmaTPCpiSemi,
+      kHistNsigmaTPCkSemi,
+      kHistNsigmaTPCpSemi,
+      kHistNsigmaTPCdSemi,
+      kHistNsigmaTPCe_eSemi,
+      kHistNsigmaTOFeSemi,
+      kHistNsigmaTOFmuSemi,
+      kHistNsigmaTOFpiSemi,
+      kHistNsigmaTOFkSemi,
+      kHistNsigmaTOFpSemi,
+      kHistNsigmaTOFdSemi,
+      kHistNsigmaTOFmismatchSemi,
+      kHistNsigmaTOFmismatch2Semi,
+      kHistDeltaTOFeSemi,
+      kHistDeltaTOFmuSemi,
+      kHistDeltaTOFpiSemi,
+      kHistDeltaTOFkSemi,
+      kHistDeltaTOFpSemi,
+      kHistDeltaTOFdSemi,
       kHistNsigmaTPCTOF,
+      kHistNsigmaTPCTOFPt,
+      kHistNsigmaTPCTOFUsed,
+      kHistNsigmaTPCTOFUsedCentral,
+      kHistNsigmaTPCTOFUsedSemiCentral,
+      kHistNsigmaTPCTOFUsedPt,
+      kHistNsigmaTPCTOFUsedPtCentral,
+      kHistNsigmaTPCTOFUsedPtSemiCentral,
       kHistEvPlane,
       kHistEvPlaneUsed,
-      kHistJetPt,
+      kHistEvPlaneCheck,
+      kHistEvPlaneCheckUsed,
+      kHistEvPlaneCorr,
+      kHistJetPtCentral,
+      kHistJetPtSemi,
       kHistEtaPhiTrgHad,
       kHistEtaPhiTrgJet,
       kHistEtaPhiAssHad,
@@ -67,20 +166,22 @@ public:
   enum Stat_t {
       kStatSeen = 1,
       kStatTrg,
-      kStatUsed,
       kStatCent,
       kStatEvPlane,
       kStatPID,
+      kStatUsed,
       kStatEvCuts,
+      kStatCentral,
+      kStatSemiCentral,
       kStatLast
   };
 
   // trigger conditions
-  enum Trigger_t { 
+  enum Trigger_t {
       kTriggerMB = 0,
       kTriggerInt,
       kTriggerLast
-  };      
+  };
 
   // classification
   enum CorrType_t {
@@ -155,6 +256,7 @@ protected:
   Float_t fZvtx; //!
   AliPIDResponse *fPIDResponse; //!
   Float_t fEventPlane; //!
+  Float_t fEventPlaneCheck; //!
   TObjArray *fPrimTrackArray; //!
   TClonesArray *fJetArray; //!
 
@@ -180,10 +282,18 @@ protected:
   Bool_t IsCentral() { return ((fCentrality >= 0.) && (fCentrality <= 10.)); }
   Bool_t IsSemiCentral() { return ((fCentrality >= 30.) && (fCentrality <= 50.)); }
 
+  AliVTrack* GetLeadingTrack(AliAODJet *jet) const;
+
+  Float_t GetDPhiStar(Float_t phi1, Float_t pt1, Float_t charge1,
+                     Float_t phi2, Float_t pt2, Float_t charge2,
+                     Float_t radius, Float_t bSign);
+
   Bool_t AcceptTrigger(AliVTrack *trg);
   Bool_t AcceptTrigger(AliAODJet *trg);
-  Bool_t AcceptAssoc(AliVTrack *ass);
+  Bool_t AcceptAssoc(AliVTrack *trk);
   Bool_t IsProton(AliVTrack *trk);
+  Bool_t AcceptAngleToEvPlane(Float_t phi, Float_t psi);
+  Bool_t AcceptTwoTracks(AliVParticle *trgPart, AliVParticle *assPart);
 
   TObjArray* CloneTracks(TObjArray *tracks) const;
 
@@ -227,14 +337,19 @@ protected:
   static const Int_t fgkStringLength = 100; // max length for the jet branch name
   char fJetBranchName[fgkStringLength];     // jet branch name
 
+  Bool_t fUseStandardCuts;
+
   AliESDtrackCuts *fCutsPrim;  // track cuts for primary particles
+  Float_t fCutsTwoTrackEff;
 
   Float_t fTrgPartPtMin;
   Float_t fTrgPartPtMax;
   Float_t fTrgJetPtMin;
   Float_t fTrgJetPtMax;
+  Float_t fTrgJetLeadTrkPtMin;
   Float_t fAssPartPtMin;
   Float_t fAssPartPtMax;
+  Float_t fTrgAngleToEvPlane;
 
   // not implemented
   AliAnalysisTaskJetProtonCorr(const AliAnalysisTaskJetProtonCorr &rhs);