PWGJE/UserTasks/AliAnalysisTaskPID.h

   1 #ifndef ALI_ANALYSIS_TASK_PID_H
   2 #define ALI_ANALYSIS_TASK_PID_H
   3
   4 /*
   5 This task collects PID output from different detectors.
   6 Only tracks fulfilling some standard quality cuts are taken into account.
   7 At the moment, only data from TPC and TOF is collected. But in future,
   8 data from e.g. HMPID is also foreseen.
   9
  10 Class written by Benjamin Hess.
  11 Contact: bhess@cern.ch
  12 */
  13
  14 class TF1;
  15 class TRandom3;
  16 class AliAnalysisFilter;
  17 class AliCFContainer;
  18 class AliESDEvent;
  19 class AliMCEvent;
  20 class AliMCParticle;
  21 class AliPID;
  22 class AliPIDCombined;
  23 class AliPIDResponse;
  24 class AliTOFPIDResponse;
  25 class AliVEvent;
  26 class AliVTrack;
  27
  28 #include "TAxis.h"
  29 #include "TH1D.h"
  30 #include "TH2D.h"
  31 #include "TH3D.h"
  32 #include "THnSparse.h"
  33 #include "TProfile.h"
  34 #include "TString.h"
  35
  36 #include "AliCentrality.h"
  37 #include "AliCFContainer.h"
  38
  39 #include "AliPID.h"
  40 #include "AliAnalysisTaskPIDV0base.h"
  41
  42 class AliAnalysisTaskPID : public AliAnalysisTaskPIDV0base {
  43  public:
  44   AliAnalysisTaskPID();
  45   AliAnalysisTaskPID(const char *name);
  46   virtual ~AliAnalysisTaskPID();
  47
  48   virtual void   UserCreateOutputObjects();
  49   virtual void   UserExec(Option_t *option);
  50   virtual void   Terminate(const Option_t*);
  51
  52   enum ErrorCode { kNoErrors = 1, kWarning = 0, kError = -1};
  53
  54   enum dataAxes { kDataMCID = 0, kDataSelectSpecies = 1, kDataPt = 2, kDataDeltaPrimeSpecies = 3, kDataCentrality = 4,
  55                   kDataJetPt = 5, kDataZ = 6, kDataXi = 7, kDataCharge = 8, kDataTOFpidInfo = 9, kDataNumAxes = 10 };
  56
  57   enum genAxes { kGenMCID = 0, kGenSelectSpecies = 1, kGenPt = 2, kGenDeltaPrimeSpecies = 3, kGenCentrality = 4,
  58                  kGenJetPt = 5, kGenZ = 6, kGenXi = 7, kGenCharge = 8, kGenTOFpidInfo = 9, kGenNumAxes = 10 };
  59
  60   enum genYieldAxes { kGenYieldMCID = 0, kGenYieldPt = 1, kGenYieldCentrality = 2, kGenYieldJetPt = 3, kGenYieldZ = 4, kGenYieldXi = 5,
  61                       kGenYieldCharge = 6, kGenYieldNumAxes = 7 };
  62
  63   enum ptResolutionAxes { kPtResJetPt = 0, kPtResGenPt = 1, kPtResRecPt = 2, kPtResCharge = 3, kPtResCentrality = 4, kPtResNumAxes = 5 };
  64
  65   enum qaSharedClsAxes { kQASharedClsJetPt = 0, kQASharedClsPt = 1, kQASharedClsNumSharedCls = 2, kQASharedClsPadRow = 3,
  66                          kQASharedClsNumAxes = 4 };
  67
  68   enum dEdxCheckAxes { kDeDxCheckPID = 0, kDeDxCheckP = 1, kDeDxCheckJetPt = 2, kDeDxCheckEtaAbs = 3 , kDeDxCheckDeDx = 4,
  69                        kDeDxCheckNumAxes = 5 };
  70
  71   enum binZeroStudyAxes { kBinZeroStudyCentrality = 0, kBinZeroStudyGenPt = 1, kBinZeroStudyGenEta = 2, kBinZeroStudyNumAxes = 3 };
  72
  73   enum efficiencyAxes { kEffMCID = 0, kEffTrackPt = 1, kEffTrackEta = 2, kEffTrackCharge = 3, kEffCentrality = 4, kEffJetPt = 5,
  74                         kEffZ = 6, kEffXi = 7, kEffNumAxes = 8 };
  75
  76   enum EffSteps { kStepGenWithGenCuts = 0, kStepRecWithGenCuts = 1, kStepRecWithGenCutsMeasuredObs = 2,
  77                   kStepRecWithRecCutsMeasuredObs = 3, kStepRecWithRecCutsMeasuredObsPrimaries = 4,
  78                   kStepRecWithRecCutsMeasuredObsStrangenessScaled = 5, kStepRecWithRecCutsPrimaries = 6, kNumSteps = 7};
  79
  80   enum TOFpidInfo { kNoTOFinfo = -2, kNoTOFpid = -1, kTOFpion = 0, kTOFkaon = 1, kTOFproton = 2, kNumTOFspecies = 3,
  81                     kNumTOFpidInfoBins = 5 };
  82
  83   enum EventCounterType { kTriggerSel = 0, kTriggerSelAndVtxCut = 1, kTriggerSelAndVtxCutAndZvtxCutNoPileUpRejection = 2,
  84                           kTriggerSelAndVtxCutAndZvtxCut = 3 };
  85
  86   enum CutHistoType { kMCPtHardCut = 0 };
  87
  88   static Int_t PDGtoMCID(Int_t pdg);
  89
  90   static void GetJetTrackObservables(Double_t trackPt, Double_t jetPt, Double_t& z, Double_t& xi);
  91
  92   static Double_t GetMCStrangenessFactorCMS(Int_t motherPDG, Double_t motherGenPt);
  93   static Double_t GetMCStrangenessFactorCMS(AliMCEvent* mcEvent, AliMCParticle* daughter);
  94
  95   static Bool_t IsSecondaryWithStrangeMotherMC(AliMCEvent* mcEvent, Int_t partLabel);
  96
  97   virtual void ConfigureTaskForCurrentEvent(AliVEvent* event);
  98
  99   Int_t GetIndexOfChargeAxisData() const
 100     { return fStoreAdditionalJetInformation ? kDataCharge : kDataCharge - fgkNumJetAxes; };
 101   Int_t GetIndexOfChargeAxisGen() const
 102     { return fStoreAdditionalJetInformation ? kGenCharge : kGenCharge - fgkNumJetAxes; };
 103   Int_t GetIndexOfChargeAxisGenYield() const
 104     { return fStoreAdditionalJetInformation ? kGenYieldCharge : kGenYieldCharge - fgkNumJetAxes; };
 105
 106   Int_t GetIndexOfTOFpidInfoAxisData() const
 107     { return fStoreAdditionalJetInformation ? kDataTOFpidInfo : kDataTOFpidInfo - fgkNumJetAxes; };
 108   Int_t GetIndexOfTOFpidInfoAxisGen() const
 109     { return fStoreAdditionalJetInformation ? kGenTOFpidInfo : kGenTOFpidInfo - fgkNumJetAxes; };
 110
 111   Bool_t FillXsec(Double_t xsection)
 112     { if (!fh1Xsec) return kFALSE; fh1Xsec->Fill("<#sigma>", xsection); return kTRUE; };
 113   Bool_t FillPythiaTrials(Double_t avgTrials)
 114     { if (!fh1Trials) return kFALSE; fh1Trials->Fill("#sum{ntrials}", avgTrials); return kTRUE; };
 115
 116   Bool_t FillEfficiencyContainer(const Double_t* values, EffSteps step, Double_t weight = 1.0);
 117
 118   Bool_t FillGeneratedYield(const Double_t* values, Double_t weight = 1.0);
 119   Bool_t FillPtResolution(Int_t mcID, const Double_t* values);
 120   Bool_t FillGenJets(Double_t centralityPercentile, Double_t jetPt, Double_t norm = -1.);
 121   Bool_t FillRecJets(Double_t centralityPercentile, Double_t jetPt, Double_t norm = -1.);
 122
 123   Bool_t IncrementEventCounter(Double_t centralityPercentile, EventCounterType type);
 124
 125   Bool_t FillCutHisto(Double_t value, CutHistoType type);
 126
 127   void PostOutputData();
 128
 129   void PrintSettings(Bool_t printSystematicsSettings = kFALSE) const;
 130
 131   void PrintSystematicsSettings() const;
 132
 133   Bool_t ProcessTrack(const AliVTrack* track, Int_t particlePDGcode, Double_t centralityPercentile, Double_t jetPt) ;
 134
 135   ErrorCode GenerateDetectorResponse(ErrorCode errCode, Double_t mean, Double_t sigma, Double_t* responses,
 136                                      Int_t nResponses,
 137                                      Bool_t usePureGaus = kFALSE);
 138   ErrorCode SetParamsForConvolutedGaus(Double_t gausMean, Double_t gausSigma);
 139
 140   const TString GetCentralityEstimator() const { return fCentralityEstimator; };
 141   const TString GetPPCentralityEstimator() const {
 142     TString ppCentEstimator = fCentralityEstimator; ppCentEstimator = ppCentEstimator.ReplaceAll("ppMult", ""); return ppCentEstimator; }
 143   void SetCentralityEstimator(TString estimator) { fCentralityEstimator = estimator; };
 144
 145   Double_t GetCentralityPercentile(AliVEvent* evt) const;
 146
 147   inline Double_t GetConvolutedGaussTransitionPar(Int_t index) const;
 148
 149   Bool_t SetConvolutedGaussLambdaParameter(Double_t lambda);
 150
 151   Bool_t GetInputFromOtherTask() const { return fInputFromOtherTask; };
 152   void SetInputFromOtherTask(Bool_t flag) { fInputFromOtherTask = flag; };
 153
 154   Bool_t GetDoPID() const { return fDoPID; };
 155   void SetDoPID(Bool_t flag) { fDoPID = flag; };
 156
 157   Bool_t GetDoEfficiency() const { return fDoEfficiency; };
 158   void SetDoEfficiency(Bool_t flag) { fDoEfficiency = flag; };
 159
 160   Bool_t GetDoPtResolution() const { return fDoPtResolution; };
 161   void SetDoPtResolution(Bool_t flag) { fDoPtResolution = flag; };
 162
 163   Bool_t GetDoDeDxCheck() const { return fDoDeDxCheck; };
 164   void SetDoDeDxCheck(Bool_t flag) { fDoDeDxCheck = flag; };
 165
 166   Bool_t GetDoBinZeroStudy() const { return fDoBinZeroStudy; };
 167   void SetDoBinZeroStudy(Bool_t flag) { fDoBinZeroStudy = flag; };
 168
 169   Bool_t GetStoreCentralityPercentile() const { return fStoreCentralityPercentile; };
 170   void SetStoreCentralityPercentile(Bool_t flag) { fStoreCentralityPercentile = flag; };
 171
 172   Bool_t GetStoreAdditionalJetInformation() const { return fStoreAdditionalJetInformation; };
 173   void SetStoreAdditionalJetInformation(Bool_t flag) { fStoreAdditionalJetInformation = flag; };
 174
 175   Bool_t GetUseMCidForGeneration() const { return fUseMCidForGeneration; };
 176   void SetUseMCidForGeneration(Bool_t flag) { fUseMCidForGeneration = flag; };
 177
 178   Bool_t GetUseConvolutedGaus() const { return fUseConvolutedGaus; };
 179   void SetUseConvolutedGaus(Bool_t flag) { fUseConvolutedGaus = flag; };
 180
 181   Double_t GetAccuracyNonGaussianTail() const { return fAccuracyNonGaussianTail; };
 182   void SetAccuracyNonGaussianTail(Double_t value) { fAccuracyNonGaussianTail = value; };
 183
 184   Bool_t GetTakeIntoAccountMuons() const { return fTakeIntoAccountMuons; };
 185   void SetTakeIntoAccountMuons(Bool_t flag) { fTakeIntoAccountMuons = flag; };
 186
 187   Int_t GetTOFmode() const { return fTOFmode; };
 188   void SetTOFmode(Int_t tofMode) { fTOFmode = tofMode; };
 189
 190   Bool_t GetUseTPCDefaultPriors() const { return fTPCDefaultPriors; };
 191   void SetUseTPCDefaultPriors(Bool_t flag) { fTPCDefaultPriors = flag; };
 192
 193   Bool_t GetUsePriors() const { return fUsePriors; };
 194   void SetUsePriors(Bool_t flag) { fUsePriors = flag; };
 195
 196   Bool_t GetUseITS() const { return fUseITS; };
 197   void SetUseITS(Bool_t flag) { fUseITS = flag; };
 198
 199   Bool_t GetUseTOF() const { return fUseTOF; };
 200   void SetUseTOF(Bool_t flag) { fUseTOF = flag; };
 201
 202   Double_t GetEtaAbsCutLow() const { return fEtaAbsCutLow; };
 203   Double_t GetEtaAbsCutUp() const { return fEtaAbsCutUp; };
 204   Bool_t SetEtaAbsCutRange(Double_t lowerLimit, Double_t upperLimit);
 205
 206   Bool_t IsInAcceptedEtaRange(Double_t etaAbs) const { return (etaAbs >= fEtaAbsCutLow && etaAbs <= fEtaAbsCutUp); };
 207
 208   AliAnalysisTaskPIDV0base::PileUpRejectionType GetPileUpRejectionType() const { return fPileUpRejectionType; };
 209   void SetPileUpRejectionType(AliAnalysisTaskPIDV0base::PileUpRejectionType newType) { fPileUpRejectionType = newType; };
 210
 211   Double_t GetSystematicScalingSplinesThreshold() const { return fSystematicScalingSplinesThreshold; };
 212   void SetSystematicScalingSplinesThreshold(Double_t threshold) { fSystematicScalingSplinesThreshold = threshold; };
 213
 214   Double_t GetSystematicScalingSplinesBelowThreshold() const { return fSystematicScalingSplinesBelowThreshold; };
 215   void SetSystematicScalingSplinesBelowThreshold(Double_t scaleFactor)
 216     { fSystematicScalingSplinesBelowThreshold = scaleFactor; CheckDoAnyStematicStudiesOnTheExpectedSignal(); };
 217
 218   Double_t GetSystematicScalingSplinesAboveThreshold() const { return fSystematicScalingSplinesAboveThreshold; };
 219   void SetSystematicScalingSplinesAboveThreshold(Double_t scaleFactor)
 220     { fSystematicScalingSplinesAboveThreshold = scaleFactor; CheckDoAnyStematicStudiesOnTheExpectedSignal(); };
 221
 222   Double_t GetSystematicScalingEtaCorrectionMomentumThr() const { return fSystematicScalingEtaCorrectionMomentumThr; };
 223   void SetSystematicScalingEtaCorrectionMomentumThr(Double_t threshold) { fSystematicScalingEtaCorrectionMomentumThr = threshold; };
 224
 225   Double_t GetSystematicScalingEtaCorrectionLowMomenta() const { return fSystematicScalingEtaCorrectionLowMomenta; };
 226   void SetSystematicScalingEtaCorrectionLowMomenta(Double_t scaleFactor)
 227     { fSystematicScalingEtaCorrectionLowMomenta = scaleFactor; CheckDoAnyStematicStudiesOnTheExpectedSignal(); };
 228
 229   Double_t GetSystematicScalingEtaCorrectionHighMomenta() const { return fSystematicScalingEtaCorrectionHighMomenta; };
 230   void SetSystematicScalingEtaCorrectionHighMomenta(Double_t scaleFactor)
 231     { fSystematicScalingEtaCorrectionHighMomenta = scaleFactor; CheckDoAnyStematicStudiesOnTheExpectedSignal(); };
 232
 233   Double_t GetSystematicScalingEtaSigmaParaThreshold() const { return fSystematicScalingEtaSigmaParaThreshold; };
 234   void SetSystematicScalingEtaSigmaParaThreshold(Double_t threshold) { fSystematicScalingEtaSigmaParaThreshold = threshold; };
 235
 236   Double_t GetSystematicScalingEtaSigmaParaBelowThreshold() const { return fSystematicScalingEtaSigmaParaBelowThreshold; };
 237   void SetSystematicScalingEtaSigmaParaBelowThreshold(Double_t scaleFactor)
 238     { fSystematicScalingEtaSigmaParaBelowThreshold = scaleFactor; CheckDoAnyStematicStudiesOnTheExpectedSignal(); };
 239
 240   Double_t GetSystematicScalingEtaSigmaParaAboveThreshold() const { return fSystematicScalingEtaSigmaParaAboveThreshold; };
 241   void SetSystematicScalingEtaSigmaParaAboveThreshold(Double_t scaleFactor)
 242     { fSystematicScalingEtaSigmaParaAboveThreshold = scaleFactor; CheckDoAnyStematicStudiesOnTheExpectedSignal(); };
 243
 244   Double_t GetSystematicScalingMultCorrection() const { return fSystematicScalingMultCorrection; };
 245   void SetSystematicScalingMultCorrection(Double_t scaleFactor)
 246     { fSystematicScalingMultCorrection = scaleFactor; CheckDoAnyStematicStudiesOnTheExpectedSignal(); };
 247
 248   Double_t GetMaxEtaVariation(Double_t dEdxSplines);
 249   Bool_t CalculateMaxEtaVariationMapFromPIDResponse();
 250
 251   void CleanupParticleFractionHistos();
 252   Bool_t GetParticleFraction(Double_t trackPt, Double_t jetPt, Double_t multiplicity,
 253                              AliPID::EParticleType species, Double_t& fraction, Double_t& fractionErrorStat,
 254                              Double_t& fractionErrorSys) const;
 255   Bool_t GetParticleFractions(Double_t trackPt, Double_t jetPt, Double_t centralityPercentile,
 256                               Double_t* prob, Int_t smearSpeciesByError, Int_t takeIntoAccountSpeciesSysError,
 257                               Bool_t uniformSystematicError = kFALSE) const;
 258   const TH3D* GetParticleFractionHisto(Int_t species, Bool_t sysError = kFALSE) const;
 259   Bool_t SetParticleFractionHisto(const TH3D* hist, Int_t species, Bool_t sysError = kFALSE);
 260   Int_t GetParticleFractionHistoNbinsTrackPt() const;
 261   Int_t GetParticleFractionHistoNbinsJetPt() const;
 262   Int_t GetParticleFractionHistoNbinsCentrality() const;
 263   Bool_t SetParticleFractionHistosFromFile(const TString filePathName, Bool_t sysError = kFALSE);
 264   Int_t GetRandomParticleTypeAccordingToParticleFractions(Double_t trackPt, Double_t jetPt,
 265                                                           Double_t centralityPercentile,
 266                                                           Bool_t smearByError,
 267                                                           Bool_t takeIntoAccountSysError = kFALSE) const;
 268
 269   TOFpidInfo GetTOFType(const AliVTrack* track, Int_t tofMode) const;
 270
 271  protected:
 272   void CheckDoAnyStematicStudiesOnTheExpectedSignal();
 273   Double_t ConvolutedGaus(const Double_t* xx, const Double_t* par) const;
 274   inline Double_t FastGaus(Double_t x, Double_t mean, Double_t sigma) const;
 275   inline Double_t FastNormalisedGaus(Double_t x, Double_t mean, Double_t sigma) const;
 276   Int_t FindBinWithinRange(TAxis* axis, Double_t value) const;
 277   Int_t FindFirstBinAboveIn3dSubset(const TH3* hist, Double_t threshold, Int_t yValue, Int_t zValue) const;
 278   Int_t FindLastBinAboveIn3dSubset(const TH3* hist, Double_t threshold, Int_t yValue, Int_t zValue) const;
 279   virtual void SetUpGenHist(THnSparse* hist, Double_t* binsPt, Double_t* binsDeltaPrime, Double_t* binsCent, Double_t* binsJetPt) const;
 280   virtual void SetUpGenYieldHist(THnSparse* hist, Double_t* binsPt, Double_t* binsCent, Double_t* binsJetPt) const;
 281   virtual void SetUpHist(THnSparse* hist, Double_t* binsPt, Double_t* binsDeltaPrime, Double_t* binsCent, Double_t* binsJetPt) const;
 282   virtual void SetUpPtResHist(THnSparse* hist, Double_t* binsPt, Double_t* binsJetPt, Double_t* binsCent) const;
 283   virtual void SetUpSharedClsHist(THnSparse* hist, Double_t* binsPt, Double_t* binsJetPt) const;
 284   virtual void SetUpDeDxCheckHist(THnSparse* hist, const Double_t* binsPt, const Double_t* binsJetPt, const Double_t* binsEtaAbs) const;
 285   virtual void SetUpBinZeroStudyHist(THnSparse* hist, const Double_t* binsCent, const Double_t* binsPt) const;
 286   virtual void SetUpPIDcombined();
 287
 288   static const Int_t fgkNumJetAxes; // Number of additional axes for jets
 289   static const Double_t fgkEpsilon; // Double_t threshold above zero
 290   static const Int_t fgkMaxNumGenEntries; // Maximum number of generated detector responses per track and delta(Prime) and associated species
 291
 292   static const Double_t fgkOneOverSqrt2; // = 1. / TMath::Sqrt2();
 293
 294  private:
 295   Int_t fRun; // Current run number
 296   AliPIDCombined* fPIDcombined; //! PID combined object
 297
 298   Bool_t fInputFromOtherTask; // If set to kTRUE, no events are processed and the input must be fed in from another task. If set to kFALSE, normal event processing
 299
 300   Bool_t fDoPID; // Do PID processing (and post the output), if flag is set to kTRUE
 301   Bool_t fDoEfficiency; // Do efficiency processing (and post the output), if flag is set to kTRUE
 302   Bool_t fDoPtResolution; // Do pT resolution processing (and post the output), if flag is set to kTRUE
 303   Bool_t fDoDeDxCheck; // Check dEdx, if flag set to kTRUE
 304   Bool_t fDoBinZeroStudy; // Do bin zero study, if flag is set to kTRUE
 305
 306   Bool_t fStoreCentralityPercentile; // If set to kTRUE, store centrality percentile for each event. In case of kFALSE (appropriate for pp), centrality percentile will be set to -1 for every event
 307   Bool_t fStoreAdditionalJetInformation; // If set to kTRUE, additional jet information like jetPt, z, xi will be stored in the THnSparses
 308
 309   Bool_t fTakeIntoAccountMuons; // Also take into account muons for the generation of the expected response and the most probable PID
 310   Bool_t fUseITS; // Use ITS for PID combined probabilities
 311   Bool_t fUseTOF; // Use TOF for PID combined probabilities
 312   Bool_t fUsePriors; // Use priors for PID combined probabilities
 313   Bool_t fTPCDefaultPriors; // Use TPC default priors for PID combined probabilities, if priors are enabled
 314
 315   Bool_t fUseMCidForGeneration; // If MC, use MCid instead of PIDcombined to generate the signals
 316
 317   Bool_t fUseConvolutedGaus; // Use convoluted gaus to generate detector response instead of pure gaus
 318   const Int_t fkConvolutedGausNPar; // Number of parameters for convoluted gaus
 319   Double_t fAccuracyNonGaussianTail; // Accuracy of the non-gaussian tail (fraction of the maximum)
 320   const Double_t fkDeltaPrimeLowLimit; // Lower deltaPrime limit
 321   const Double_t fkDeltaPrimeUpLimit; // Upper deltaPrime limit
 322   TF1* fConvolutedGausDeltaPrime; // Gaus convoluted with exponential tail to generate detector response (deltaPrime)
 323
 324   Int_t fTOFmode; // TOF mode used for TOF PID info (affects num sigma inclusion/exclusion)
 325   Double_t fConvolutedGaussTransitionPars[3]; // Parameter for transition from gaussian parameters to asymmetric shape
 326   static const Double_t fgkSigmaReferenceForTransitionPars; // Reference sigma chosen to calculate transition parameters
 327
 328   Double_t fEtaAbsCutLow; // Lower cut value on |eta|
 329   Double_t fEtaAbsCutUp;  // Upper cut value on |eta|
 330
 331   AliAnalysisTaskPIDV0base::PileUpRejectionType fPileUpRejectionType; // Which pile-up rejection is used (if any)
 332
 333   // For systematic studies
 334   Bool_t   fDoAnySystematicStudiesOnTheExpectedSignal; // Internal flag indicating whether any systematic studies are going to be performed
 335   Double_t fSystematicScalingSplinesThreshold;         // beta-gamma threshold for the systematic spline scale factor
 336   Double_t fSystematicScalingSplinesBelowThreshold;        // Systematic scale factor for the splines (1. = no systematics) below threshold
 337   Double_t fSystematicScalingSplinesAboveThreshold;        // Systematic scale factor for the splines (1. = no systematics) above threshold
 338   Double_t fSystematicScalingEtaCorrectionMomentumThr;  // Momentum threshold for the systematic scale factor for the eta correction (separates low-p from high-p
 339   Double_t fSystematicScalingEtaCorrectionLowMomenta;   // Systematic scale factor for the eta correction (1. = no systematics) at low momenta
 340   Double_t fSystematicScalingEtaCorrectionHighMomenta;  // Systematic scale factor for the eta correction (1. = no systematics) at high momenta
 341
 342   Double_t fSystematicScalingEtaSigmaParaThreshold; // dEdx threshold for the systematic scale factor for the parametrisation of the relative signal width
 343   Double_t fSystematicScalingEtaSigmaParaBelowThreshold; // Systematic scale factor for the parametrisation of the relative signal width (1. = no systematics) below threshold
 344   Double_t fSystematicScalingEtaSigmaParaAboveThreshold; // Systematic scale factor for the parametrisation of the relative signal width (1. = no systematics) above threshold
 345   Double_t fSystematicScalingMultCorrection; // Systematic scale factor for the multiplicity correction (1. = no systematics)
 346
 347   TH3D* fFractionHists[AliPID::kSPECIES]; // 3D histos of particle fraction as function  with trackPt, jetPt (-1 for inclusive spectra), centralityPercentile (-1 for pp)
 348   TH3D* fFractionSysErrorHists[AliPID::kSPECIES]; // 3D histos of sys. error of particle fraction as function  with trackPt, jetPt (-1 for inclusive spectra), centralityPercentile (-1 for pp)
 349
 350   TString fCentralityEstimator; // Estimator for the centrality (e.g. V0A, V0M)
 351
 352   THnSparseD* fhPIDdataAll; //! Data histo
 353
 354   // Generated response information
 355   THnSparseD* fhGenEl; //! Generated response for el
 356   THnSparseD* fhGenKa; //! Generated response for ka
 357   THnSparseD* fhGenPi; //! Generated response for pi
 358   THnSparseD* fhGenMu; //! Generated response for mu
 359   THnSparseD* fhGenPr; //! Generated response for pr
 360
 361   // Generated responses for a single track
 362   Double_t* fGenRespElDeltaPrimeEl; //! Generated responses for a single track
 363   Double_t* fGenRespElDeltaPrimeKa; //! Generated responses for a single track
 364   Double_t* fGenRespElDeltaPrimePi; //! Generated responses for a single track
 365   Double_t* fGenRespElDeltaPrimePr; //! Generated responses for a single track
 366   Double_t* fGenRespKaDeltaPrimeEl; //! Generated responses for a single track
 367   Double_t* fGenRespKaDeltaPrimeKa; //! Generated responses for a single track
 368   Double_t* fGenRespKaDeltaPrimePi; //! Generated responses for a single track
 369   Double_t* fGenRespKaDeltaPrimePr; //! Generated responses for a single track
 370   Double_t* fGenRespPiDeltaPrimeEl; //! Generated responses for a single track
 371   Double_t* fGenRespPiDeltaPrimeKa; //! Generated responses for a single track
 372   Double_t* fGenRespPiDeltaPrimePi; //! Generated responses for a single track
 373   Double_t* fGenRespPiDeltaPrimePr; //! Generated responses for a single track
 374   Double_t* fGenRespMuDeltaPrimeEl; //! Generated responses for a single track
 375   Double_t* fGenRespMuDeltaPrimeKa; //! Generated responses for a single track
 376   Double_t* fGenRespMuDeltaPrimePi; //! Generated responses for a single track
 377   Double_t* fGenRespMuDeltaPrimePr; //! Generated responses for a single track
 378   Double_t* fGenRespPrDeltaPrimeEl; //! Generated responses for a single track
 379   Double_t* fGenRespPrDeltaPrimeKa; //! Generated responses for a single track
 380   Double_t* fGenRespPrDeltaPrimePi; //! Generated responses for a single track
 381   Double_t* fGenRespPrDeltaPrimePr; //! Generated responses for a single track
 382   /*
 383   Double_t* fGenRespElDeltaEl; //! Generated responses for a single track
 384   Double_t* fGenRespElDeltaKa; //! Generated responses for a single track
 385   Double_t* fGenRespElDeltaPi; //! Generated responses for a single track
 386   Double_t* fGenRespElDeltaPr; //! Generated responses for a single track
 387   Double_t* fGenRespKaDeltaEl; //! Generated responses for a single track
 388   Double_t* fGenRespKaDeltaKa; //! Generated responses for a single track
 389   Double_t* fGenRespKaDeltaPi; //! Generated responses for a single track
 390   Double_t* fGenRespKaDeltaPr; //! Generated responses for a single track
 391   Double_t* fGenRespPiDeltaEl; //! Generated responses for a single track
 392   Double_t* fGenRespPiDeltaKa; //! Generated responses for a single track
 393   Double_t* fGenRespPiDeltaPi; //! Generated responses for a single track
 394   Double_t* fGenRespPiDeltaPr; //! Generated responses for a single track
 395   Double_t* fGenRespMuDeltaEl; //! Generated responses for a single track
 396   Double_t* fGenRespMuDeltaKa; //! Generated responses for a single track
 397   Double_t* fGenRespMuDeltaPi; //! Generated responses for a single track
 398   Double_t* fGenRespMuDeltaPr; //! Generated responses for a single track
 399   Double_t* fGenRespPrDeltaEl; //! Generated responses for a single track
 400   Double_t* fGenRespPrDeltaKa; //! Generated responses for a single track
 401   Double_t* fGenRespPrDeltaPi; //! Generated responses for a single track
 402   Double_t* fGenRespPrDeltaPr; //! Generated responses for a single track
 403   */
 404
 405   TAxis* fDeltaPrimeAxis; //! Axis holding the deltaPrime binning
 406   TH1D* fhMaxEtaVariation; //! Histo holding the maximum deviation of the eta correction factor from unity vs. 1/dEdx(splines)
 407
 408   TH1D* fhEventsProcessed; //! Histo holding the number of processed events (i.e. passing trigger selection, vtx and zvtx cuts and (if enabled) pile-up rejection)
 409   TH1D* fhEventsTriggerSel; //! Histo holding the number of events passing trigger selection
 410   TH1D* fhEventsTriggerSelVtxCut; //! Histo holding the number of events passing trigger selection and vtx cut
 411   TH1D* fhEventsProcessedNoPileUpRejection; //! Histo holding the number of processed events before pile-up rejection
 412
 413   THnSparseD* fChargedGenPrimariesTriggerSel; //! Histo holding the generated charged primary yields for triggered events
 414   THnSparseD* fChargedGenPrimariesTriggerSelVtxCut; //! Histo holding the generated charged primary yields for triggered events passing vertex cuts
 415   THnSparseD* fChargedGenPrimariesTriggerSelVtxCutZ; //! Histo holding the generated charged primary yields for triggered events passing vertex cuts (including cut on z)
 416   THnSparseD* fChargedGenPrimariesTriggerSelVtxCutZPileUpRej; //! Histo holding the generated charged primary yields for triggered events passing vertex cuts (including cut on z) and pile-up rejection
 417
 418   THnSparseD* fhMCgeneratedYieldsPrimaries; //! Histo holding the generated (no reco, no cuts) primary particle yields in considered eta range
 419
 420   TH2D* fh2FFJetPtRec;            //! Number of reconstructed jets vs. jetPt and centrality
 421   TH2D* fh2FFJetPtGen;            //! Number of generated jets vs. jetPt and centrality
 422
 423   TProfile* fh1Xsec;              //! pythia cross section and trials
 424   TH1D*     fh1Trials;            //! sum of trials
 425   TH1F*     fh1EvtsPtHardCut;     //! Number events before and after the cut on MC pT hard
 426
 427   AliCFContainer* fContainerEff; //! Container for efficiency determination
 428
 429   THnSparseD* fPtResolution[AliPID::kSPECIES]; //! Pt Resolution for the individual species
 430   THnSparseD* fQASharedCls; //! QA for shared clusters
 431
 432   THnSparseD* fDeDxCheck; //! dEdx check
 433
 434   TObjArray* fOutputContainer;  //! output data container
 435
 436   TObjArray* fQAContainer; //! output data container for QA
 437
 438   AliAnalysisTaskPID(const AliAnalysisTaskPID&); // not implemented
 439   AliAnalysisTaskPID& operator=(const AliAnalysisTaskPID&); // not implemented
 440
 441   ClassDef(AliAnalysisTaskPID, 22);
 442 };
 443
 444
 445 //_____________________________________________________________________________
 446 inline Bool_t AliAnalysisTaskPID::FillEfficiencyContainer(const Double_t* values, AliAnalysisTaskPID::EffSteps step,
 447                                                           Double_t weight)
 448 {
 449   // Fill efficiency container at step "step" with the values
 450
 451   if (!fDoEfficiency)
 452     return kFALSE;
 453
 454   if (!fContainerEff) {
 455     AliError("Efficiency container not initialised -> cannot be filled!");
 456     return kFALSE;
 457   }
 458
 459   fContainerEff->Fill(values, step, weight);
 460
 461   return kTRUE;
 462 }
 463
 464
 465 //_____________________________________________________________________________
 466 inline Bool_t AliAnalysisTaskPID::FillGeneratedYield(const Double_t* values, Double_t weight)
 467 {
 468   // Fill histos with generated primary yields with provided values
 469
 470   if (!fDoPID)
 471     return kFALSE;
 472
 473   if (!fhMCgeneratedYieldsPrimaries) {
 474     AliError("Histo for generated primary yield not initialised -> cannot be filled!");
 475     return kFALSE;
 476   }
 477
 478   fhMCgeneratedYieldsPrimaries->Fill(values, weight);
 479
 480   return kTRUE;
 481 }
 482
 483
 484 //_____________________________________________________________________________
 485 inline Bool_t AliAnalysisTaskPID::FillGenJets(Double_t centralityPercentile, Double_t jetPt, Double_t norm)
 486 {
 487   if (!fDoPID && !fDoEfficiency)
 488     return kFALSE;
 489
 490   if (!fh2FFJetPtGen)
 491     return kFALSE;
 492
 493   if (norm > 0.)
 494     fh2FFJetPtGen->Fill(centralityPercentile, jetPt, 1. / norm);
 495   else
 496     fh2FFJetPtGen->Fill(centralityPercentile, jetPt);
 497
 498   return kTRUE;
 499 }
 500
 501
 502 //_____________________________________________________________________________
 503 inline Bool_t AliAnalysisTaskPID::FillRecJets(Double_t centralityPercentile, Double_t jetPt, Double_t norm)
 504 {
 505   if (!fDoPID && !fDoEfficiency)
 506     return kFALSE;
 507
 508   if (!fh2FFJetPtRec)
 509     return kFALSE;
 510
 511   if (norm > 0.)
 512     fh2FFJetPtRec->Fill(centralityPercentile, jetPt, 1. / norm);
 513   else
 514     fh2FFJetPtRec->Fill(centralityPercentile, jetPt);
 515
 516   return kTRUE;
 517 }
 518
 519
 520 //_____________________________________________________________________________
 521 inline Bool_t AliAnalysisTaskPID::FillPtResolution(Int_t mcID, const Double_t* values)
 522 {
 523   // Fill histos with pT resolution with provided values
 524
 525   if (!fDoPtResolution || mcID < 0 || mcID >= AliPID::kSPECIES)
 526     return kFALSE;
 527
 528   if (!fPtResolution[mcID]) {
 529     AliError(Form("Histo for pT resolution (species: %s) not initialised -> cannot be filled!", AliPID::ParticleName(mcID)));
 530     return kFALSE;
 531   }
 532
 533   fPtResolution[mcID]->Fill(values);
 534
 535   return kTRUE;
 536 }
 537
 538
 539 //_____________________________________________________________________________
 540 inline Bool_t AliAnalysisTaskPID::IncrementEventCounter(Double_t centralityPercentile, AliAnalysisTaskPID::EventCounterType type)
 541 {
 542   // Increment the number of events for the given centrality percentile and the corresponding counter type
 543
 544   if (type == kTriggerSel) {
 545     if (!fhEventsTriggerSel) {
 546       AliError("Histogram for number of events (kTriggerSel) not initialised -> cannot be incremented!");
 547       return kFALSE;
 548     }
 549
 550     fhEventsTriggerSel->Fill(centralityPercentile);
 551   }
 552   else if (type == kTriggerSelAndVtxCut) {
 553     if (!fhEventsTriggerSelVtxCut) {
 554       AliError("Histogram for number of events (kTriggerSelAndVtxCut) not initialised -> cannot be incremented!");
 555       return kFALSE;
 556     }
 557
 558     fhEventsTriggerSelVtxCut->Fill(centralityPercentile);
 559   }
 560   else if (type == kTriggerSelAndVtxCutAndZvtxCut) {
 561     if (!fhEventsProcessed) {
 562       AliError("Histogram for number of events (kTriggerSelAndVtxCutAndZvtxCut) not initialised -> cannot be incremented!");
 563       return kFALSE;
 564     }
 565
 566     fhEventsProcessed->Fill(centralityPercentile);
 567   }
 568   else if (type == kTriggerSelAndVtxCutAndZvtxCutNoPileUpRejection) {
 569     if (!fhEventsProcessedNoPileUpRejection) {
 570       AliError("Histogram for number of events (kTriggerSelAndVtxCutAndZvtxCutNoPileUpRejection) not initialised -> cannot be incremented!");
 571       return kFALSE;
 572     }
 573
 574     fhEventsProcessedNoPileUpRejection->Fill(centralityPercentile);
 575   }
 576
 577
 578   return kTRUE;
 579 };
 580
 581
 582 //_____________________________________________________________________________
 583 inline Bool_t AliAnalysisTaskPID::FillCutHisto(Double_t value, AliAnalysisTaskPID::CutHistoType type)
 584 {
 585   // Fill cut histo of corresponding type at given value.
 586   if (type == kMCPtHardCut) {
 587     if (!fh1EvtsPtHardCut) {
 588       AliError("Histogram \"fh1EvtsPtHardCut\" not initialised -> cannot be filled!");
 589       return kFALSE;
 590     }
 591
 592     fh1EvtsPtHardCut->Fill(value);
 593   }
 594
 595   return kTRUE;
 596 }
 597
 598
 599 //_____________________________________________________________________________
 600 inline Bool_t AliAnalysisTaskPID::SetEtaAbsCutRange(Double_t lowerLimit, Double_t upperLimit)
 601 {
 602   if (lowerLimit >= upperLimit) {
 603     AliError(Form("Requested lower |eta| cut limit >= upper |eta| cut limit. Old eta cut range will be used (low %f, high %f).",
 604                   fEtaAbsCutLow, fEtaAbsCutUp));
 605     return kFALSE;
 606   }
 607
 608   fEtaAbsCutLow = lowerLimit;
 609   fEtaAbsCutUp = upperLimit;
 610
 611   return kTRUE;
 612 };
 613
 614
 615 //_____________________________________________________________________________
 616 inline Double_t AliAnalysisTaskPID::GetConvolutedGaussTransitionPar(Int_t index) const
 617 {
 618   if (index < 0 || index >= 3) {
 619     printf("Invalid index %d!\n", index);
 620     return -1;
 621   }
 622   return fConvolutedGaussTransitionPars[index];
 623 }
 624
 625
 626 //_____________________________________________________________________________
 627 inline Int_t AliAnalysisTaskPID::GetParticleFractionHistoNbinsTrackPt() const
 628 {
 629   if (!fFractionHists[AliPID::kPion])
 630     return -1;
 631
 632   return fFractionHists[AliPID::kPion]->GetNbinsX();
 633 }
 634
 635
 636 //_____________________________________________________________________________
 637 inline Int_t AliAnalysisTaskPID::GetParticleFractionHistoNbinsJetPt() const
 638 {
 639   if (!fFractionHists[AliPID::kPion])
 640     return -1;
 641
 642   return fFractionHists[AliPID::kPion]->GetNbinsY();
 643 }
 644
 645
 646 //_____________________________________________________________________________
 647 inline Int_t AliAnalysisTaskPID::GetParticleFractionHistoNbinsCentrality() const
 648 {
 649   if (!fFractionHists[AliPID::kPion])
 650     return -1;
 651
 652   return fFractionHists[AliPID::kPion]->GetNbinsZ();
 653 }
 654
 655
 656 //_____________________________________________________________________________
 657 inline Double_t AliAnalysisTaskPID::GetCentralityPercentile(AliVEvent* evt) const
 658 {
 659   // WARNING: This function may not be used in case of special pp centrality estimators which require different handling
 660   // (and sometimes ESD events)
 661   if (!evt)
 662     return -1;
 663
 664   Double_t centralityPercentile = -1.;
 665   if (fStoreCentralityPercentile)
 666     centralityPercentile = evt->GetCentrality()->GetCentralityPercentile(fCentralityEstimator.Data());
 667
 668   return centralityPercentile;
 669 }
 670
 671
 672 //_____________________________________________________________________________
 673 inline void AliAnalysisTaskPID::PostOutputData()
 674 {
 675   PostData(1, fOutputContainer);
 676
 677   if (fDoEfficiency)
 678     PostData(2, fContainerEff);
 679
 680   if (fDoPtResolution || fDoDeDxCheck)
 681     PostData(3, fQAContainer);
 682 }
 683
 684 #endif