Pleasmall bug related to zero inclusive jet PID tasks -
authorChristian Klein-Boesing <Christian.Klein-Boesing@cern.ch>
Thu, 9 Jan 2014 15:16:01 +0000 (16:16 +0100)
committerChristian Klein-Boesing <Christian.Klein-Boesing@cern.ch>
Thu, 9 Jan 2014 15:21:11 +0000 (16:21 +0100)
 Fixed number of generated responses to same (rather high)
 number for all pT and jet pT - Added possibility to set
 (absolute) eta range for PID part - Added TOF dimensions -
 Added debug prints - Fixed end of line - Cleaned up code -
 Moved init of static consts to cxx - Fixed PostData in
 UserCreateOutputObjects

PWGJE/UserTasks/AliAnalysisTaskIDFragmentationFunction.cxx
PWGJE/UserTasks/AliAnalysisTaskPID.cxx
PWGJE/UserTasks/AliAnalysisTaskPID.h

index 3d507fd..0f7d55b 100644 (file)
@@ -2292,7 +2292,7 @@ void AliAnalysisTaskIDFragmentationFunction::UserCreateOutputObjects()
         }
       }
       else {
-        Printf("ERROR: zero jet pid tasks!\n");
+        Printf("WARNING: zero jet pid tasks!\n");
         fUseJetPIDtask = kFALSE;
       }
     }
@@ -2314,8 +2314,8 @@ void AliAnalysisTaskIDFragmentationFunction::UserCreateOutputObjects()
         }
       }
       else {
-        Printf("ERROR: zero inclusive pid tasks!\n");
-        fUseJetPIDtask = kFALSE;
+        Printf("WARNING: zero inclusive pid tasks!\n");
+        fUseInclusivePIDtask = kFALSE;
       }
     }
   }
@@ -2725,8 +2725,10 @@ void AliAnalysisTaskIDFragmentationFunction::UserExec(Option_t *)
           Double_t valuesGenYield[AliAnalysisTaskPID::kGenYieldNumAxes] = { mcID, pT, centPercent, -1, -1, -1, -1 };
 
           for (Int_t i = 0; i < fNumInclusivePIDtasks; i++) {
-            valuesGenYield[fInclusivePIDtask[i]->GetIndexOfChargeAxisGenYield()] = chargeMC;
-            fInclusivePIDtask[i]->FillGeneratedYield(valuesGenYield);
+            if (fInclusivePIDtask[i]->IsInAcceptedEtaRange(TMath::Abs(part->Eta()))) {
+              valuesGenYield[fInclusivePIDtask[i]->GetIndexOfChargeAxisGenYield()] = chargeMC;
+              fInclusivePIDtask[i]->FillGeneratedYield(valuesGenYield);
+            }
           }
           
           Double_t valuesEff[AliAnalysisTaskPID::kEffNumAxes] = { mcID, pT, part->Eta(), chargeMC,
@@ -2839,8 +2841,10 @@ void AliAnalysisTaskIDFragmentationFunction::UserExec(Option_t *)
             for (Int_t i = 0; i < fNumInclusivePIDtasks; i++) {
               if ((!fInclusivePIDtask[i]->GetUseTPCCutMIGeo() && !fInclusivePIDtask[i]->GetUseTPCnclCut()) ||
                   (survivedTPCCutMIGeo && fInclusivePIDtask[i]->GetUseTPCCutMIGeo()) ||
-                  (survivedTPCnclCut && fInclusivePIDtask[i]->GetUseTPCnclCut()))
-                fInclusivePIDtask[i]->ProcessTrack(inclusiveaod, pdg, centPercent, -1); // no jet pT since inclusive spectrum 
+                  (survivedTPCnclCut && fInclusivePIDtask[i]->GetUseTPCnclCut())) {
+                    if (fInclusivePIDtask[i]->IsInAcceptedEtaRange(TMath::Abs(inclusiveaod->Eta())))
+                      fInclusivePIDtask[i]->ProcessTrack(inclusiveaod, pdg, centPercent, -1); // no jet pT since inclusive spectrum 
+              }
             }
             
             if (gentrack) {
@@ -3016,8 +3020,10 @@ void AliAnalysisTaskIDFragmentationFunction::UserExec(Option_t *)
             Double_t valuesGenYield[AliAnalysisTaskPID::kGenYieldNumAxes] = { mcID, trackPt, centPercent, jetPt, z, xi, chargeMC };
             
             for (Int_t i = 0; i < fNumJetPIDtasks; i++) {
-              valuesGenYield[fJetPIDtask[i]->GetIndexOfChargeAxisGenYield()] = chargeMC;
-              fJetPIDtask[i]->FillGeneratedYield(valuesGenYield);
+              if (fJetPIDtask[i]->IsInAcceptedEtaRange(TMath::Abs(part->Eta()))) {
+                valuesGenYield[fJetPIDtask[i]->GetIndexOfChargeAxisGenYield()] = chargeMC;
+                fJetPIDtask[i]->FillGeneratedYield(valuesGenYield);
+              }
             }
             
             
@@ -3213,8 +3219,10 @@ void AliAnalysisTaskIDFragmentationFunction::UserExec(Option_t *)
         for (Int_t i = 0; i < fNumJetPIDtasks; i++) {
           if ((!fJetPIDtask[i]->GetUseTPCCutMIGeo() && !fJetPIDtask[i]->GetUseTPCnclCut()) ||
               (survivedTPCCutMIGeo && fJetPIDtask[i]->GetUseTPCCutMIGeo()) ||
-              (survivedTPCnclCut && fJetPIDtask[i]->GetUseTPCnclCut()))
-            fJetPIDtask[i]->ProcessTrack(aodtrack, pdg, centPercent, jetPt);
+              (survivedTPCnclCut && fJetPIDtask[i]->GetUseTPCnclCut())) {
+                if (fJetPIDtask[i]->IsInAcceptedEtaRange(TMath::Abs(aodtrack->Eta())))
+                  fJetPIDtask[i]->ProcessTrack(aodtrack, pdg, centPercent, jetPt);
+          }
         }
         
         if (fIDFFMode && ((!fJetPIDtask[0]->GetUseTPCCutMIGeo() && !fJetPIDtask[0]->GetUseTPCnclCut()) ||
index 31e4b98..616b413 100644 (file)
@@ -37,9 +37,13 @@ Contact: bhess@cern.ch
 
 ClassImp(AliAnalysisTaskPID)
 
+const Int_t AliAnalysisTaskPID::fgkNumJetAxes = 3; // Number of additional axes for jets
 const Double_t AliAnalysisTaskPID::fgkEpsilon = 1e-8; // Double_t threshold above zero
+const Int_t AliAnalysisTaskPID::fgkMaxNumGenEntries = 500; // Maximum number of generated detector responses per track and delta(Prime) and associated species
+
 const Double_t AliAnalysisTaskPID::fgkOneOverSqrt2 = 0.707106781186547462; // = 1. / TMath::Sqrt2();
-const Double_t AliAnalysisTaskPID::fgkSigmaReferenceForTransitionPars = 0.05; // Reference sigma chosen to calculate transition parameters
+
+const Double_t AliAnalysisTaskPID::fgkSigmaReferenceForTransitionPars = 0.05; // Reference sigma chosen to calculate transition
 
 //________________________________________________________________________
 AliAnalysisTaskPID::AliAnalysisTaskPID()
@@ -63,6 +67,7 @@ AliAnalysisTaskPID::AliAnalysisTaskPID()
   , fkDeltaPrimeLowLimit(0.02)
   , fkDeltaPrimeUpLimit(40.0)
   , fConvolutedGausDeltaPrime(0x0)
+  , fTOFmode(1)
   , fEtaAbsCutLow(0.0)
   , fEtaAbsCutUp(0.9)
   , fDoAnySystematicStudiesOnTheExpectedSignal(kFALSE)
@@ -185,6 +190,7 @@ AliAnalysisTaskPID::AliAnalysisTaskPID(const char *name)
   , fkDeltaPrimeLowLimit(0.02)
   , fkDeltaPrimeUpLimit(40.0)
   , fConvolutedGausDeltaPrime(0x0)
+  , fTOFmode(1)
   , fEtaAbsCutLow(0.0)
   , fEtaAbsCutUp(0.9)
   , fDoAnySystematicStudiesOnTheExpectedSignal(kFALSE)
@@ -554,45 +560,73 @@ void AliAnalysisTaskPID::UserCreateOutputObjects()
   const Double_t xiMin = 0.;
   const Double_t xiMax = 7.;
   
+  const Int_t nTOFpidInfoBins = kNumTOFpidInfoBins;
+  const Double_t tofPIDinfoMin = kNoTOFinfo;
+  const Double_t tofPIDinfoMax = kNoTOFinfo + kNumTOFpidInfoBins;
+  
   // MC PID, SelectSpecies, pT, deltaPrimeSpecies, centrality percentile, jet pT, z = track_pT/jet_pT, xi = log(1/z)
-  Int_t binsNoJets[nBinsNoJets] =    { nMCPIDbins, nSelSpeciesBins,         nPtBins,                 deltaPrimeNBins,  
-                                       nCentBins, nChargeBins};
-  Int_t binsJets[nBinsJets]     =    { nMCPIDbins, nSelSpeciesBins,         nPtBins,                 deltaPrimeNBins,           
-                                       nCentBins,                      nJetPtBins, nZBins, nXiBins, nChargeBins };
+  Int_t binsNoJets[nBinsNoJets] =    { nMCPIDbins,
+                                       nSelSpeciesBins,
+                                       nPtBins,
+                                       deltaPrimeNBins,
+                                       nCentBins,
+                                       nChargeBins,
+                                       nTOFpidInfoBins };
+  
+  Int_t binsJets[nBinsJets]     =    { nMCPIDbins,
+                                       nSelSpeciesBins,
+                                       nPtBins,
+                                       deltaPrimeNBins,           
+                                       nCentBins,
+                                       nJetPtBins,
+                                       nZBins,
+                                       nXiBins,
+                                       nChargeBins,
+                                       nTOFpidInfoBins };
+  
   Int_t *bins = fStoreAdditionalJetInformation ? &binsJets[0] : &binsNoJets[0];
   
-  Double_t xminNoJets[nBinsNoJets] = {   mcPIDmin,   selSpeciesMin,       binsPt[0],               deltaPrimeBins[0],                       
-                                       binsCent[0], binsCharge[0]};
-  Double_t xminJets[nBinsJets] =     {   mcPIDmin,   selSpeciesMin,       binsPt[0],               deltaPrimeBins[0],                       
-                                       binsCent[0],                  binsJetPt[0],   zMin,   xiMin, binsCharge[0] };
+  Double_t xminNoJets[nBinsNoJets] = { mcPIDmin,  
+                                       selSpeciesMin,
+                                       binsPt[0],
+                                       deltaPrimeBins[0],                       
+                                       binsCent[0],
+                                       binsCharge[0],
+                                       tofPIDinfoMin };
+  
+  Double_t xminJets[nBinsJets] =     { mcPIDmin,
+                                       selSpeciesMin,
+                                       binsPt[0],
+                                       deltaPrimeBins[0],                       
+                                       binsCent[0],
+                                       binsJetPt[0],
+                                       zMin,
+                                       xiMin,
+                                       binsCharge[0],
+                                       tofPIDinfoMin };
+  
   Double_t *xmin = fStoreAdditionalJetInformation? &xminJets[0] : &xminNoJets[0];
 
-  Double_t xmaxNoJets[nBinsNoJets] = {   mcPIDmax,   selSpeciesMax, binsPt[nPtBins], deltaPrimeBins[deltaPrimeNBins], 
-                                       binsCent[nCentBins], binsCharge[nChargeBins] };
-  Double_t xmaxJets[nBinsJets] =     {   mcPIDmax,   selSpeciesMax, binsPt[nPtBins], deltaPrimeBins[deltaPrimeNBins], 
-                                       binsCent[nCentBins], binsJetPt[nJetPtBins],   zMax,   xiMax, binsCharge[nChargeBins] };
-  Double_t *xmax = fStoreAdditionalJetInformation? &xmaxJets[0] : &xmaxNoJets[0];
+  Double_t xmaxNoJets[nBinsNoJets] = { mcPIDmax,
+                                       selSpeciesMax,
+                                       binsPt[nPtBins],
+                                       deltaPrimeBins[deltaPrimeNBins], 
+                                       binsCent[nCentBins],
+                                       binsCharge[nChargeBins],
+                                       tofPIDinfoMax };
+  
+  Double_t xmaxJets[nBinsJets] =     { mcPIDmax,
+                                       selSpeciesMax,
+                                       binsPt[nPtBins],
+                                       deltaPrimeBins[deltaPrimeNBins], 
+                                       binsCent[nCentBins],
+                                       binsJetPt[nJetPtBins],
+                                       zMax,
+                                       xiMax,
+                                       binsCharge[nChargeBins],
+                                       tofPIDinfoMax };
   
-  /*OLD with TOF, p_TPC_Inner and p_vertex and deltaSpecies
-  const Int_t nBins = 8;
-  //TODO In case of memory trouble: Remove deltaTOFspecies and p(Vertex) (can be added later, if needed)?
-  //TODO IF everything is working fine: p(TPC_inner) and p(Vertex) can be removed, since everything in the analysis is only a 
-  // function of pT -> Reduces memory consumption significantly!!!
-  
-  // MC PID, SelectSpecies, P(TPC_inner), pT, p(Vertex), deltaSpecies, deltaPrimeSpecies, deltaTOFspecies
-  const Int_t deltaPrimeNBins = 201;
-  
-  const Int_t deltaNBins = 801;
-  const Double_t deltaLowLimit = -200.;
-  const Double_t deltaUpLimit = 200.;
-  
-  Int_t bins[nBins] = 
-    {  5, 4, nPtBins, nPtBins, nPtBins, deltaNBins, deltaPrimeNBins, 501 }; 
-  Double_t xmin[nBins] = 
-    {  0., 0., 0., 0., 0., deltaLowLimit, fkDeltaPrimeLowLimit, -5000.};
-  Double_t xmax[nBins] = 
-    {  5., 4., 50.0, 50.0, 50.0, deltaUpLimit, fkDeltaPrimeUpLimit, 5000.};
-  */
+  Double_t *xmax = fStoreAdditionalJetInformation? &xmaxJets[0] : &xmaxNoJets[0];
   
   fConvolutedGausDeltaPrime->SetNpx(deltaPrimeNBins);
 
@@ -800,7 +834,9 @@ void AliAnalysisTaskPID::UserCreateOutputObjects()
   if(fDebug > 2)
     printf("File: %s, Line: %d: UserCreateOutputObjects -> Posting output data\n", (char*)__FILE__, __LINE__);
   
-  PostOutputData();
+  PostData(1, fOutputContainer);
+  PostData(2, fContainerEff);
+  PostData(3, fPtResolutionContainer);
   
   if(fDebug > 2)
     printf("File: %s, Line: %d: UserCreateOutputObjects -> Done\n", (char*)__FILE__, __LINE__);
@@ -868,7 +904,7 @@ void AliAnalysisTaskPID::UserExec(Option_t *)
         // - Species must be one of those in question (everything else goes to the overflow bin of mcID)
         
         // Geometrie should be the same as on the reconstructed level -> By definition analysis within this eta interval
-        if (TMath::Abs(mcPart->Eta()) < fEtaAbsCutLow || TMath::Abs(mcPart->Eta()) > fEtaAbsCutUp)  continue;
+        if (!IsInAcceptedEtaRange(TMath::Abs(mcPart->Eta()))) continue;
         
         Int_t mcID = PDGtoMCID(mcPart->PdgCode());
         
@@ -954,7 +990,7 @@ void AliAnalysisTaskPID::UserExec(Option_t *)
         // For efficiency: Reconstructed track has survived all cuts on the detector level (excluding acceptance)
         // and has an associated MC track which is a physical primary and was generated inside the acceptance
         if (fMC->IsPhysicalPrimary(TMath::Abs(label)) &&
-            (TMath::Abs(mcTrack->Eta()) >= fEtaAbsCutLow && TMath::Abs(mcTrack->Eta()) <= fEtaAbsCutUp)) {
+            IsInAcceptedEtaRange(TMath::Abs(mcTrack->Eta()))) {
           
           // AliMCParticle->Charge() calls TParticlePDG->Charge(), which returns the charge in units of e0 / 3
           Double_t value[kEffNumAxes] = { mcID, mcTrack->Pt(), mcTrack->Eta(), mcTrack->Charge() / 3., centralityPercentile,
@@ -969,7 +1005,7 @@ void AliAnalysisTaskPID::UserExec(Option_t *)
     }
    
     // Only process tracks inside the desired eta window    
-    if (TMath::Abs(track->Eta()) < fEtaAbsCutLow || TMath::Abs(track->Eta()) > fEtaAbsCutUp)  continue;
+    if (!IsInAcceptedEtaRange(TMath::Abs(track->Eta()))) continue;
    
     if (fDoPID) 
       ProcessTrack(track, pdg, centralityPercentile, -1); // No jet information in this case -> Set jet pT to -1
@@ -1681,6 +1717,57 @@ Double_t AliAnalysisTaskPID::GetMCStrangenessFactorCMS(AliMCEvent* mcEvent, AliM
 
 
 // _________________________________________________________________________________
+AliAnalysisTaskPID::TOFpidInfo AliAnalysisTaskPID::GetTOFType(const AliVTrack* track, Int_t tofMode) const
+{
+  // Get the (locally defined) particle type judged by TOF
+
+  if (!fPIDResponse) {
+    Printf("ERROR: fEvent not available -> Cannot determine TOF type!");
+    return kNoTOFinfo;
+  }
+   
+  // Check kTOFout, kTIME, mismatch
+  const AliPIDResponse::EDetPidStatus tofStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF, track);
+  if (tofStatus != AliPIDResponse::kDetPidOk)
+    return kNoTOFinfo;
+
+  Double_t nsigma[kNumTOFspecies] = { -999., -999., -999. };
+  nsigma[kTOFpion]   = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kPion);
+  nsigma[kTOFkaon]   = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kKaon);
+  nsigma[kTOFproton] = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton);
+
+  Double_t inclusion = -999;
+  Double_t exclusion = -999;
+  
+  if (tofMode == 0) {
+    inclusion = 1.5;
+    exclusion = 3;
+  }
+  else if (tofMode == 1) {
+    inclusion = 2;
+    exclusion = 3;
+  }
+  else if (tofMode == 2) {
+    inclusion = 2.5;
+    exclusion = 3;
+  }
+  else {
+    Printf("ERROR: Bad TOF mode: %d!", tofMode);
+    return kNoTOFinfo;
+  }
+
+  if (TMath::Abs(nsigma[kTOFpion]) < inclusion && TMath::Abs(nsigma[kTOFkaon]) > exclusion && TMath::Abs(nsigma[kTOFproton]) > exclusion)
+    return kTOFpion;
+  if (TMath::Abs(nsigma[kTOFpion]) > exclusion && TMath::Abs(nsigma[kTOFkaon]) < inclusion && TMath::Abs(nsigma[kTOFproton]) > exclusion)
+    return kTOFkaon;
+  if (TMath::Abs(nsigma[kTOFpion]) > exclusion && TMath::Abs(nsigma[kTOFkaon]) > exclusion && TMath::Abs(nsigma[kTOFproton]) < inclusion)
+    return kTOFproton;
+
+  return kNoTOFpid;
+}
+
+
+// _________________________________________________________________________________
 Bool_t AliAnalysisTaskPID::IsSecondaryWithStrangeMotherMC(AliMCEvent* mcEvent, Int_t partLabel)
 {
   // Check whether particle is a secondary with strange mother, i.e. returns kTRUE if a strange mother is found
@@ -1906,6 +1993,7 @@ void AliAnalysisTaskPID::PrintSettings(Bool_t printSystematicsSettings) const
   printf("Use convoluted Gauss: %d\n", GetUseConvolutedGaus());
   printf("Accuracy of non-Gaussian tail: %e\n", GetAccuracyNonGaussianTail());
   printf("Take into account muons: %d\n", GetTakeIntoAccountMuons());
+  printf("TOF mode: %d\n", GetTOFmode());
   printf("\nParams for transition from gauss to asymmetric shape:\n");
   printf("[0]: %e\n", GetConvolutedGaussTransitionPar(0));
   printf("[1]: %e\n", GetConvolutedGaussTransitionPar(1));
@@ -2191,12 +2279,6 @@ Bool_t AliAnalysisTaskPID::ProcessTrack(const AliVTrack* track, Int_t particlePD
                                                               fPIDResponse->UseTPCEtaCorrection(),
                                                               fPIDResponse->UseTPCMultiplicityCorrection()); 
   }
-  /*OLD with deltaSpecies
-  Double_t deltaElectron = dEdxTPC - dEdxEl;
-  Double_t deltaKaon = dEdxTPC - dEdxKa;
-  Double_t deltaPion = dEdxTPC - dEdxPi;
-  Double_t deltaProton = dEdxTPC - dEdxPr;
-  */
   
   Double_t deltaPrimeElectron = (dEdxEl > 0) ? dEdxTPC / dEdxEl : -1;
   if (dEdxEl <= 0)  {
@@ -2221,17 +2303,6 @@ Bool_t AliAnalysisTaskPID::ProcessTrack(const AliVTrack* track, Int_t particlePD
     Printf("Error: Expected TPC signal <= 0 for proton hypothesis");
     return kFALSE;
   }
-      
-  /*TODO for TOF 
-  // TOF signal
-  Double_t times[AliPID::kSPECIES];
-  track->GetIntegratedTimes(times);
-  AliTOFPIDResponse tofPIDResponse = fPIDResponse->GetTOFResponse();
-  Float_t electronDeltaTOF = GetDeltaTOF(track, &tofPIDResponse, times, AliPID::kElectron);
-  Float_t pionDeltaTOF = GetDeltaTOF(track, &tofPIDResponse, times, AliPID::kPion);
-  Float_t kaonDeltaTOF = GetDeltaTOF(track, &tofPIDResponse, times, AliPID::kKaon);
-  Float_t protonDeltaTOF = GetDeltaTOF(track, &tofPIDResponse, times, AliPID::kProton);
-  */
   
   if(fDebug > 2)
     printf("File: %s, Line: %d: ProcessTrack -> Compute probabilities\n", (char*)__FILE__, __LINE__);
@@ -2302,6 +2373,7 @@ Bool_t AliAnalysisTaskPID::ProcessTrack(const AliVTrack* track, Int_t particlePD
   pTPC = temp;
   */
   
+  TOFpidInfo tofPIDinfo = GetTOFType(track, fTOFmode);
   
   Double_t entry[fStoreAdditionalJetInformation ? kDataNumAxes : kDataNumAxes - fgkNumJetAxes];
   entry[kDataMCID] = binMC;
@@ -2317,30 +2389,20 @@ Bool_t AliAnalysisTaskPID::ProcessTrack(const AliVTrack* track, Int_t particlePD
   }
   
   entry[GetIndexOfChargeAxisData()] = trackCharge;
-  
-  /*OLD with TOF, p_TPC_Inner and p_vertex and deltaSpecies
-  // MC PID, SelectSpecies, P(TPC_inner), pT, p(Vertex), deltaSpecies, deltaPrimeSpecies, deltaTOFspecies
-  Double_t entry[8] = { binMC, 0, pTPC, pT, p, deltaElectron, deltaPrimeElectron, electronDeltaTOF };
-  */
+  entry[GetIndexOfTOFpidInfoAxisData()] = tofPIDinfo;
   
   fhPIDdataAll->Fill(entry);
   
   entry[kDataSelectSpecies] = 1;
-  //OLD with deltaSpecies entry[5] = deltaKaon;
   entry[kDataDeltaPrimeSpecies] = deltaPrimeKaon;
-  //TODO for TOF entry[7] = kaonDeltaTOF;
   fhPIDdataAll->Fill(entry);
     
   entry[kDataSelectSpecies] = 2;
-  //OLD with deltaSpecies entry[5] = deltaPion;
   entry[kDataDeltaPrimeSpecies] = deltaPrimePion;
-  //TODO for TOF entry[7] = pionDeltaTOF;
   fhPIDdataAll->Fill(entry);
     
   entry[kDataSelectSpecies] = 3;
-  //OLD with deltaSpecies entry[5] = deltaProton;
   entry[kDataDeltaPrimeSpecies] = deltaPrimeProton;
-  //TODO for TOF entry[7] = protonDeltaTOF;
   fhPIDdataAll->Fill(entry);
   
   
@@ -2360,17 +2422,19 @@ Bool_t AliAnalysisTaskPID::ProcessTrack(const AliVTrack* track, Int_t particlePD
   }
   
   genEntry[GetIndexOfChargeAxisGen()] = trackCharge;
+  genEntry[GetIndexOfTOFpidInfoAxisGen()] = tofPIDinfo;
   
-  //OLD with deltaSpecies Double_t genEntry[5] = { binMC, 0, pT, -999, -999 }; // MC PID, SelectSpecies, pT, deltaSpecies, deltaPrimeSpecies
+  // Generate numGenEntries "responses" with fluctuations around the expected values.
+  // fgkMaxNumGenEntries = 500 turned out to give reasonable templates even for highest track pT in 15-20 GeV/c jet pT bin.
+  Int_t numGenEntries = fgkMaxNumGenEntries; // fgkMaxNumGenEntries = 500
   
+  /*OLD: Different number of responses depending on pT and jet pT for fgkMaxNumGenEntries = 1000
+   * => Problem: If threshold to higher number of responses inside a bin (or after rebinning), then the template
+   * is biased to the higher pT.
   // Generate numGenEntries "responses" with fluctuations around the expected values.
   // The higher the (transverse) momentum, the more "responses" will be generated in order not to run out of statistics too fast.
   Int_t numGenEntries = 10;
-  if (pT >= 5) 
-    numGenEntries *= 5;
-  else if (pT >= 2)
-    numGenEntries *= 2;
-  
   // Jets have even worse statistics, therefore, scale numGenEntries further
   if (jetPt >= 40)
     numGenEntries *= 20;
@@ -2380,10 +2444,13 @@ Bool_t AliAnalysisTaskPID::ProcessTrack(const AliVTrack* track, Int_t particlePD
     numGenEntries *= 2;
   
   
+  
   // Do not generate more entries than available in memory!
   if (numGenEntries > fgkMaxNumGenEntries)// fgkMaxNumGenEntries = 1000
     numGenEntries = fgkMaxNumGenEntries;
-      
+  */
+  
+  
   ErrorCode errCode = kNoErrors;
   
   if(fDebug > 2)
@@ -2421,57 +2488,22 @@ Bool_t AliAnalysisTaskPID::ProcessTrack(const AliVTrack* track, Int_t particlePD
   errCode = GenerateDetectorResponse(errCode, dEdxPr / dEdxPi, sigmaPr / dEdxPi, fGenRespPrDeltaPrimePi, numGenEntries);
   errCode = GenerateDetectorResponse(errCode, 1.,              sigmaPr / dEdxPr, fGenRespPrDeltaPrimePr, numGenEntries);
   
-  
-  /*OLD with deltaSpecies 
-  // Delta
-  
-  // Electrons
-  errCode = GenerateDetectorResponse(errCode, 0.,              sigmaEl, fGenRespElDeltaEl, numGenEntries, usePureGausForDelta);
-  errCode = GenerateDetectorResponse(errCode, dEdxEl - dEdxKa, sigmaEl, fGenRespElDeltaKa, numGenEntries, usePureGausForDelta);
-  errCode = GenerateDetectorResponse(errCode, dEdxEl - dEdxPi, sigmaEl, fGenRespElDeltaPi, numGenEntries, usePureGausForDelta);
-  errCode = GenerateDetectorResponse(errCode, dEdxEl - dEdxPr, sigmaEl, fGenRespElDeltaPr, numGenEntries, usePureGausForDelta);
-  
-  // Kaons
-  errCode = GenerateDetectorResponse(errCode, dEdxKa - dEdxEl, sigmaKa, fGenRespKaDeltaEl, numGenEntries, usePureGausForDelta);
-  errCode = GenerateDetectorResponse(errCode, 0.,              sigmaKa, fGenRespKaDeltaKa, numGenEntries, usePureGausForDelta);
-  errCode = GenerateDetectorResponse(errCode, dEdxKa - dEdxPi, sigmaKa, fGenRespKaDeltaPi, numGenEntries, usePureGausForDelta);
-  errCode = GenerateDetectorResponse(errCode, dEdxKa - dEdxPr, sigmaKa, fGenRespKaDeltaPr, numGenEntries, usePureGausForDelta);
-  
-  // Pions
-  errCode = GenerateDetectorResponse(errCode, dEdxPi - dEdxEl, sigmaPi, fGenRespPiDeltaEl, numGenEntries, usePureGausForDelta);
-  errCode = GenerateDetectorResponse(errCode, dEdxPi - dEdxKa, sigmaPi, fGenRespPiDeltaKa, numGenEntries, usePureGausForDelta);
-  errCode = GenerateDetectorResponse(errCode, 0.,              sigmaPi, fGenRespPiDeltaPi, numGenEntries, usePureGausForDelta);
-  errCode = GenerateDetectorResponse(errCode, dEdxPi - dEdxPr, sigmaPi, fGenRespPiDeltaPr, numGenEntries, usePureGausForDelta);
-  
-  // Muons
-  errCode = GenerateDetectorResponse(errCode, dEdxMu - dEdxEl, sigmaMu, fGenRespMuDeltaEl, numGenEntries, usePureGausForDelta);
-  errCode = GenerateDetectorResponse(errCode, dEdxMu - dEdxKa, sigmaMu, fGenRespMuDeltaKa, numGenEntries, usePureGausForDelta);
-  errCode = GenerateDetectorResponse(errCode, dEdxMu - dEdxPi, sigmaMu, fGenRespMuDeltaPi, numGenEntries, usePureGausForDelta);
-  errCode = GenerateDetectorResponse(errCode, dEdxMu - dEdxPr, sigmaMu, fGenRespMuDeltaPr, numGenEntries, usePureGausForDelta);
-  
-  // Protons
-  errCode = GenerateDetectorResponse(errCode, dEdxPr - dEdxEl, sigmaPr, fGenRespPrDeltaEl, numGenEntries, usePureGausForDelta);
-  errCode = GenerateDetectorResponse(errCode, dEdxPr - dEdxKa, sigmaPr, fGenRespPrDeltaKa, numGenEntries, usePureGausForDelta);
-  errCode = GenerateDetectorResponse(errCode, dEdxPr - dEdxPi, sigmaPr, fGenRespPrDeltaPi, numGenEntries, usePureGausForDelta);
-  errCode = GenerateDetectorResponse(errCode, 0.,              sigmaPr, fGenRespPrDeltaPr, numGenEntries, usePureGausForDelta);
-  */
-  
   if (errCode != kNoErrors) {
-    if (errCode == kWarning) {
-      //Printf("Warning: Questionable detector response for track, most likely due to very low number of PID clusters! Debug output (dEdx_expected, sigma_expected):");
+    if (errCode == kWarning && fDebug > 1) {
+      Printf("Warning: Questionable detector response for track, most likely due to very low number of PID clusters! Debug output (dEdx_expected, sigma_expected):");
     }
     else 
       Printf("Error: Failed to generate detector response for track - skipped! Debug output (dEdx_expected, sigma_expected):");
     
-    /*
-    Printf("Pr: %e, %e", dEdxPr, sigmaPr);
-    Printf("Pi: %e, %e", dEdxPi, sigmaPi);
-    Printf("El: %e, %e", dEdxEl, sigmaEl);
-    Printf("Mu: %e, %e", dEdxMu, sigmaMu);
-    Printf("Ka: %e, %e", dEdxKa, sigmaKa);
-    Printf("track: dEdx %f, pTPC %f, eta %f, ncl %d\n", track->GetTPCsignal(), track->GetTPCmomentum(), track->Eta(), 
-           track->GetTPCsignalN());
-    */
+    if (fDebug > 1) {
+      Printf("Pr: %e, %e", dEdxPr, sigmaPr);
+      Printf("Pi: %e, %e", dEdxPi, sigmaPi);
+      Printf("El: %e, %e", dEdxEl, sigmaEl);
+      Printf("Mu: %e, %e", dEdxMu, sigmaMu);
+      Printf("Ka: %e, %e", dEdxKa, sigmaKa);
+      Printf("track: dEdx %f, pTPC %f, eta %f, ncl %d\n", track->GetTPCsignal(), track->GetTPCmomentum(), track->Eta(), 
+            track->GetTPCsignalN());
+    }
     
     if (errCode != kWarning) {
       fhSkippedTracksForSignalGeneration->Fill(track->GetTPCmomentum(), track->GetTPCsignalN());
@@ -2499,108 +2531,88 @@ Bool_t AliAnalysisTaskPID::ProcessTrack(const AliVTrack* track, Int_t particlePD
     
     // Electrons
     genEntry[kGenSelectSpecies] = 0;
-    //OLD with deltaSpecies genEntry[3] = fGenRespElDeltaEl[n];
     genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimeEl[n];
     fhGenEl->Fill(genEntry);
     
     genEntry[kGenSelectSpecies] = 1;
-    //OLD with deltaSpecies genEntry[3] = fGenRespElDeltaKa[n];
     genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimeKa[n];
     fhGenEl->Fill(genEntry);
     
     genEntry[kGenSelectSpecies] = 2;
-    //OLD with deltaSpecies genEntry[3] = fGenRespElDeltaPi[n];
     genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimePi[n];
     fhGenEl->Fill(genEntry);
     
     genEntry[kGenSelectSpecies] = 3;
-    //OLD with deltaSpecies genEntry[3] = fGenRespElDeltaPr[n];
     genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimePr[n];
     fhGenEl->Fill(genEntry);
     
     // Kaons
     genEntry[kGenSelectSpecies] = 0;
-    //OLD with deltaSpecies genEntry[3] = fGenRespKaDeltaEl[n];
     genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimeEl[n];
     fhGenKa->Fill(genEntry);
     
     genEntry[kGenSelectSpecies] = 1;
-    //OLD with deltaSpecies genEntry[3] = fGenRespKaDeltaKa[n];
     genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimeKa[n];
     fhGenKa->Fill(genEntry);
     
     genEntry[kGenSelectSpecies] = 2;
-    //OLD with deltaSpecies genEntry[3] = fGenRespKaDeltaPi[n];
     genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimePi[n];
     fhGenKa->Fill(genEntry);
     
     genEntry[kGenSelectSpecies] = 3;
-    //OLD with deltaSpecies genEntry[3] = fGenRespKaDeltaPr[n];
     genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimePr[n];
     fhGenKa->Fill(genEntry);
     
     // Pions
     genEntry[kGenSelectSpecies] = 0;
-    //OLD with deltaSpecies genEntry[3] = fGenRespPiDeltaEl[n];
     genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimeEl[n];
     fhGenPi->Fill(genEntry);
     
     genEntry[kGenSelectSpecies] = 1;
-    //OLD with deltaSpecies genEntry[3] = fGenRespPiDeltaKa[n];
     genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimeKa[n];
     fhGenPi->Fill(genEntry);
     
     genEntry[kGenSelectSpecies] = 2;
-    //OLD with deltaSpecies genEntry[3] = fGenRespPiDeltaPi[n];
     genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimePi[n];
     fhGenPi->Fill(genEntry);
     
     genEntry[kGenSelectSpecies] = 3;
-    //OLD with deltaSpecies genEntry[3] = fGenRespPiDeltaPr[n];
     genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimePr[n];
     fhGenPi->Fill(genEntry);
     
     if (fTakeIntoAccountMuons) {
       // Muons, if desired
       genEntry[kGenSelectSpecies] = 0;
-      //OLD with deltaSpecies genEntry[3] = fGenRespMuDeltaEl[n];
       genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimeEl[n];
       fhGenMu->Fill(genEntry);
       
       genEntry[kGenSelectSpecies] = 1;
-      //OLD with deltaSpecies genEntry[3] = fGenRespMuDeltaKa[n];
       genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimeKa[n];
       fhGenMu->Fill(genEntry);
       
       genEntry[kGenSelectSpecies] = 2;
-      //OLD with deltaSpecies genEntry[3] = fGenRespMuDeltaPi[n];
       genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimePi[n];
       fhGenMu->Fill(genEntry);
       
       genEntry[kGenSelectSpecies] = 3;
-      //OLD with deltaSpecies genEntry[3] = fGenRespMuDeltaPr[n];
       genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimePr[n];
       fhGenMu->Fill(genEntry);
     }
     
     // Protons
     genEntry[kGenSelectSpecies] = 0;
-    //OLD with deltaSpecies genEntry[3] = fGenRespPrDeltaEl[n];
     genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimeEl[n];
     fhGenPr->Fill(genEntry);
     
     genEntry[kGenSelectSpecies] = 1;
-    //OLD with deltaSpecies genEntry[3] = fGenRespPrDeltaKa[n];
     genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimeKa[n];
     fhGenPr->Fill(genEntry);
     
     genEntry[kGenSelectSpecies] = 2;
-    //OLD with deltaSpecies genEntry[3] = fGenRespPrDeltaPi[n];
     genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimePi[n];
     fhGenPr->Fill(genEntry);
     
     genEntry[kGenSelectSpecies] = 3;
-    //OLD with deltaSpecies genEntry[3] = fGenRespPrDeltaPr[n];
     genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimePr[n];
     fhGenPr->Fill(genEntry);
   }
@@ -2974,6 +2986,14 @@ void AliAnalysisTaskPID::SetUpGenHist(THnSparse* hist, Double_t* binsPt, Double_
   }
   
   hist->GetAxis(GetIndexOfChargeAxisGen())->SetTitle("Charge (e_{0})");
+  
+  hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetTitle("TOF PID Info");
+  // Offset is (TMath::Abs(kNoTOFinfo) + 1), such that bin 1 (= first label) corresponds to kNoTOFinfo (< 0)
+  hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kNoTOFinfo + (TMath::Abs(kNoTOFinfo) + 1), "No TOF Info");
+  hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kNoTOFpid + (TMath::Abs(kNoTOFinfo) + 1), "No TOF PID");
+  hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kTOFpion + (TMath::Abs(kNoTOFinfo) + 1), "#pi");
+  hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kTOFkaon + (TMath::Abs(kNoTOFinfo) + 1), "K");
+  hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kTOFproton + (TMath::Abs(kNoTOFinfo) + 1), "p");
 }
 
 
@@ -3049,36 +3069,13 @@ void AliAnalysisTaskPID::SetUpHist(THnSparse* hist, Double_t* binsPt, Double_t*
   
   hist->GetAxis(GetIndexOfChargeAxisData())->SetTitle("Charge (e_{0})");
   
-  /*OLD with TOF, p_TPC_Inner and p_vertex
-  // MC PID, SelectSpecies, P(TPC_inner), pT, p(Vertex), deltaSpecies, deltaPrimeSpecies, deltaTOFspecies
-  hist->SetBinEdges(2, binsPt);
-  hist->SetBinEdges(3, binsPt);
-  hist->SetBinEdges(4, binsPt);
-                          
-  // Set axes titles
-  hist->GetAxis(0)->SetTitle("MC PID");
-  hist->GetAxis(0)->SetBinLabel(1, "e");
-  hist->GetAxis(0)->SetBinLabel(2, "K");
-  hist->GetAxis(0)->SetBinLabel(3, "#mu");
-  hist->GetAxis(0)->SetBinLabel(4, "#pi");
-  hist->GetAxis(0)->SetBinLabel(5, "p");
-  
-  hist->GetAxis(1)->SetTitle("Select Species");
-  hist->GetAxis(1)->SetBinLabel(1, "e");
-  hist->GetAxis(1)->SetBinLabel(2, "K");
-  hist->GetAxis(1)->SetBinLabel(3, "#pi");
-  hist->GetAxis(1)->SetBinLabel(4, "p");
-  
-  hist->GetAxis(2)->SetTitle("P_{TPC_inner} (GeV/c)");
-  hist->GetAxis(3)->SetTitle("P_{T} (GeV/c)");
-  hist->GetAxis(4)->SetTitle("P_{vertex} (GeV/c)");
-  
-  hist->GetAxis(5)->SetTitle("TPC #Delta_{species} (arb. unit)");
-  
-  hist->GetAxis(6)->SetTitle("TPC #Delta'_{species} (arb. unit)");
-  
-  hist->GetAxis(7)->SetTitle("#Delta TOF_{species} (ps)");
-  */
+  hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetTitle("TOF PID Info");
+  // Offset is (TMath::Abs(kNoTOFinfo) + 1), such that bin 1 (= first label) corresponds to kNoTOFinfo (< 0)
+  hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kNoTOFinfo + (TMath::Abs(kNoTOFinfo) + 1), "No TOF Info");
+  hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kNoTOFpid + (TMath::Abs(kNoTOFinfo) + 1), "No TOF PID");
+  hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kTOFpion + (TMath::Abs(kNoTOFinfo) + 1), "#pi");
+  hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kTOFkaon + (TMath::Abs(kNoTOFinfo) + 1), "K");
+  hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kTOFproton + (TMath::Abs(kNoTOFinfo) + 1), "p");
 }
 
 
index ad53ea2..bb80b1e 100644 (file)
-#ifndef ALI_ANALYSIS_TASK_PID_H\r
-#define ALI_ANALYSIS_TASK_PID_H\r
-\r
-/*\r
-This task collects PID output from different detectors.\r
-Only tracks fulfilling some standard quality cuts are taken into account.\r
-At the moment, only data from TPC and TOF is collected. But in future,\r
-data from e.g. HMPID is also foreseen.\r
-\r
-Class written by Benjamin Hess.\r
-Contact: bhess@cern.ch\r
-*/\r
-\r
-class TF1;\r
-class TRandom3;\r
-class AliAnalysisFilter;\r
-class AliCFContainer;\r
-class AliESDEvent;\r
-class AliMCEvent;\r
-class AliMCParticle;\r
-class AliPID;\r
-class AliPIDCombined;\r
-class AliPIDResponse;\r
-class AliTOFPIDResponse;\r
-class AliVEvent;\r
-class AliVTrack;\r
-\r
-#include "TH1D.h"\r
-#include "TH2D.h"\r
-#include "TH3D.h"\r
-#include "THnSparse.h"\r
-#include "TProfile.h"\r
-#include "TString.h"\r
-\r
-#include "AliCentrality.h"\r
-#include "AliCFContainer.h"\r
-\r
-#include "AliPID.h"\r
-#include "AliAnalysisTaskPIDV0base.h"\r
-\r
-class AliAnalysisTaskPID : public AliAnalysisTaskPIDV0base {\r
- public:\r
-  AliAnalysisTaskPID();\r
-  AliAnalysisTaskPID(const char *name);\r
-  virtual ~AliAnalysisTaskPID();\r
-  \r
-  virtual void   UserCreateOutputObjects();\r
-  virtual void   UserExec(Option_t *option);\r
-  virtual void   Terminate(const Option_t*);\r
-  \r
-  enum ErrorCode { kNoErrors = 1, kWarning = 0, kError = -1};\r
-  \r
-  enum dataAxes { kDataMCID = 0, kDataSelectSpecies = 1, kDataPt = 2, kDataDeltaPrimeSpecies = 3, kDataCentrality = 4,\r
-                  kDataJetPt = 5, kDataZ = 6, kDataXi = 7, kDataCharge = 8, kDataNumAxes = 9 };\r
-\r
-  enum genAxes { kGenMCID = 0, kGenSelectSpecies = 1, kGenPt = 2, kGenDeltaPrimeSpecies = 3, kGenCentrality = 4,\r
-                 kGenJetPt = 5, kGenZ = 6, kGenXi = 7, kGenCharge = 8, kGenNumAxes = 9 };\r
-  \r
-  enum genYieldAxes { kGenYieldMCID = 0, kGenYieldPt = 1, kGenYieldCentrality = 2, kGenYieldJetPt = 3, kGenYieldZ = 4, kGenYieldXi = 5,\r
-                      kGenYieldCharge = 6, kGenYieldNumAxes = 7 };\r
-  \r
-  enum ptResolutionAxes { kPtResJetPt = 0, kPtResGenPt = 1, kPtResRecPt = 2, kPtResCharge = 3, kPtResCentrality = 4, kPtResNumAxes = 5 };\r
-  \r
-  enum efficiencyAxes { kEffMCID = 0, kEffTrackPt = 1, kEffTrackEta = 2, kEffTrackCharge = 3, kEffCentrality = 4, kEffJetPt = 5,\r
-                        kEffZ = 6, kEffXi = 7, kEffNumAxes = 8 };\r
-  \r
-  enum EffSteps { kStepGenWithGenCuts = 0, kStepRecWithGenCuts = 1, kStepRecWithGenCutsMeasuredObs = 2,\r
-                  kStepRecWithRecCutsMeasuredObs = 3, kStepRecWithRecCutsMeasuredObsPrimaries = 4,\r
-                  kStepRecWithRecCutsMeasuredObsStrangenessScaled = 5, kStepRecWithRecCutsPrimaries = 6, kNumSteps = 7};\r
-  \r
-  static Int_t PDGtoMCID(Int_t pdg);\r
-  \r
-  static void GetJetTrackObservables(Double_t trackPt, Double_t jetPt, Double_t& z, Double_t& xi);\r
-  \r
-  static Double_t GetMCStrangenessFactorCMS(Int_t motherPDG, Double_t motherGenPt);\r
-  static Double_t GetMCStrangenessFactorCMS(AliMCEvent* mcEvent, AliMCParticle* daughter);\r
-  \r
-  static Bool_t IsSecondaryWithStrangeMotherMC(AliMCEvent* mcEvent, Int_t partLabel);\r
-  \r
-  Int_t GetIndexOfChargeAxisData() const\r
-    { return fStoreAdditionalJetInformation ? kDataCharge : kDataCharge - fgkNumJetAxes; };\r
-  Int_t GetIndexOfChargeAxisGen() const\r
-    { return fStoreAdditionalJetInformation ? kGenCharge : kGenCharge - fgkNumJetAxes; };\r
-  Int_t GetIndexOfChargeAxisGenYield() const\r
-    { return fStoreAdditionalJetInformation ? kGenYieldCharge : kGenYieldCharge - fgkNumJetAxes; };\r
-  \r
-  Bool_t FillXsec(Double_t xsection)\r
-    { if (!fh1Xsec) return kFALSE; fh1Xsec->Fill("<#sigma>", xsection); return kTRUE; };\r
-  Bool_t FillPythiaTrials(Double_t avgTrials)\r
-    { if (!fh1Trials) return kFALSE; fh1Trials->Fill("#sum{ntrials}", avgTrials); return kTRUE; };\r
-    \r
-  Bool_t FillEfficiencyContainer(const Double_t* values, EffSteps step, Double_t weight = 1.0);\r
-  \r
-  Bool_t FillGeneratedYield(const Double_t* values, Double_t weight = 1.0);\r
-  Bool_t FillPtResolution(Int_t mcID, const Double_t* values);\r
-  Bool_t FillGenJets(Double_t centralityPercentile, Double_t jetPt, Double_t norm = -1.);\r
-  Bool_t FillRecJets(Double_t centralityPercentile, Double_t jetPt, Double_t norm = -1.);\r
-  \r
-  Bool_t IncrementEventsProcessed(Double_t centralityPercentile);\r
-  \r
-  void PostOutputData();\r
-  \r
-  void PrintSettings(Bool_t printSystematicsSettings = kFALSE) const;\r
-  \r
-  void PrintSystematicsSettings() const;\r
-  \r
-  Bool_t ProcessTrack(const AliVTrack* track, Int_t particlePDGcode, Double_t centralityPercentile, Double_t jetPt) ;\r
-  \r
-  ErrorCode GenerateDetectorResponse(ErrorCode errCode, Double_t mean, Double_t sigma, Double_t* responses,\r
-                                     Int_t nResponses,\r
-                                     Bool_t usePureGaus = kFALSE);\r
-  ErrorCode SetParamsForConvolutedGaus(Double_t gausMean, Double_t gausSigma);\r
-  \r
-  const TString GetCentralityEstimator() const { return fCentralityEstimator; };\r
-  void SetCentralityEstimator(TString estimator) { fCentralityEstimator = estimator; };\r
-  \r
-  Double_t GetCentralityPercentile(AliVEvent* evt) const;\r
-  \r
-  inline Double_t GetConvolutedGaussTransitionPar(Int_t index) const;\r
-  \r
-  Bool_t SetConvolutedGaussLambdaParameter(Double_t lambda);\r
-\r
-  Bool_t GetInputFromOtherTask() const { return fInputFromOtherTask; };\r
-  void SetInputFromOtherTask(Bool_t flag) { fInputFromOtherTask = flag; };\r
-  \r
-  Bool_t GetDoPID() const { return fDoPID; };\r
-  void SetDoPID(Bool_t flag) { fDoPID = flag; };\r
-  \r
-  Bool_t GetDoEfficiency() const { return fDoEfficiency; };\r
-  void SetDoEfficiency(Bool_t flag) { fDoEfficiency = flag; };\r
-  \r
-  Bool_t GetDoPtResolution() const { return fDoPtResolution; };\r
-  void SetDoPtResolution(Bool_t flag) { fDoPtResolution = flag; };\r
-  \r
-  Bool_t GetStoreCentralityPercentile() const { return fStoreCentralityPercentile; };\r
-  void SetStoreCentralityPercentile(Bool_t flag) { fStoreCentralityPercentile = flag; };\r
-  \r
-  Bool_t GetStoreAdditionalJetInformation() const { return fStoreAdditionalJetInformation; };\r
-  void SetStoreAdditionalJetInformation(Bool_t flag) { fStoreAdditionalJetInformation = flag; };\r
-  \r
-  Bool_t GetUseMCidForGeneration() const { return fUseMCidForGeneration; };\r
-  void SetUseMCidForGeneration(Bool_t flag) { fUseMCidForGeneration = flag; };\r
-  \r
-  Bool_t GetUseConvolutedGaus() const { return fUseConvolutedGaus; };\r
-  void SetUseConvolutedGaus(Bool_t flag) { fUseConvolutedGaus = flag; };\r
-  \r
-  Double_t GetAccuracyNonGaussianTail() const { return fAccuracyNonGaussianTail; };\r
-  void SetAccuracyNonGaussianTail(Double_t value) { fAccuracyNonGaussianTail = value; };\r
-  \r
-  Bool_t GetTakeIntoAccountMuons() const { return fTakeIntoAccountMuons; };\r
-  void SetTakeIntoAccountMuons(Bool_t flag) { fTakeIntoAccountMuons = flag; };\r
-  \r
-  Bool_t GetUseTPCDefaultPriors() const { return fTPCDefaultPriors; };\r
-  void SetUseTPCDefaultPriors(Bool_t flag) { fTPCDefaultPriors = flag; };\r
-  \r
-  Bool_t GetUsePriors() const { return fUsePriors; };\r
-  void SetUsePriors(Bool_t flag) { fUsePriors = flag; };\r
-  \r
-  Bool_t GetUseITS() const { return fUseITS; };\r
-  void SetUseITS(Bool_t flag) { fUseITS = flag; };\r
-  \r
-  Bool_t GetUseTOF() const { return fUseTOF; };\r
-  void SetUseTOF(Bool_t flag) { fUseTOF = flag; };\r
-  \r
-  Double_t GetEtaAbsCutLow() const { return fEtaAbsCutLow; };\r
-  Double_t GetEtaAbsCutUp() const { return fEtaAbsCutUp; };\r
-  Bool_t SetEtaAbsCutRange(Double_t lowerLimit, Double_t upperLimit);\r
-  \r
-  Double_t GetSystematicScalingSplines() const { return fSystematicScalingSplines; };\r
-  void SetSystematicScalingSplines(Double_t scaleFactor) \r
-    { fSystematicScalingSplines = scaleFactor; CheckDoAnyStematicStudiesOnTheExpectedSignal(); };\r
-  \r
-  Double_t GetSystematicScalingEtaCorrectionMomentumThr() const { return fSystematicScalingEtaCorrectionMomentumThr; };\r
-  void SetSystematicScalingEtaCorrectionMomentumThr(Double_t threshold) { fSystematicScalingEtaCorrectionMomentumThr = threshold; };\r
-  \r
-  Double_t GetSystematicScalingEtaCorrectionLowMomenta() const { return fSystematicScalingEtaCorrectionLowMomenta; };\r
-  void SetSystematicScalingEtaCorrectionLowMomenta(Double_t scaleFactor) \r
-    { fSystematicScalingEtaCorrectionLowMomenta = scaleFactor; CheckDoAnyStematicStudiesOnTheExpectedSignal(); };\r
-  \r
-  Double_t GetSystematicScalingEtaCorrectionHighMomenta() const { return fSystematicScalingEtaCorrectionHighMomenta; };\r
-  void SetSystematicScalingEtaCorrectionHighMomenta(Double_t scaleFactor) \r
-    { fSystematicScalingEtaCorrectionHighMomenta = scaleFactor; CheckDoAnyStematicStudiesOnTheExpectedSignal(); };\r
-  \r
-  Double_t GetSystematicScalingEtaSigmaPara() const { return fSystematicScalingEtaSigmaPara; };\r
-  void SetSystematicScalingEtaSigmaPara(Double_t scaleFactor)\r
-    { fSystematicScalingEtaSigmaPara = scaleFactor; CheckDoAnyStematicStudiesOnTheExpectedSignal(); };\r
-  \r
-  Double_t GetSystematicScalingMultCorrection() const { return fSystematicScalingMultCorrection; };\r
-  void SetSystematicScalingMultCorrection(Double_t scaleFactor) \r
-    { fSystematicScalingMultCorrection = scaleFactor; CheckDoAnyStematicStudiesOnTheExpectedSignal(); };\r
-  \r
-  void CleanupParticleFractionHistos();\r
-  Bool_t GetParticleFraction(Double_t trackPt, Double_t jetPt, Double_t multiplicity,\r
-                             AliPID::EParticleType species, Double_t& fraction, Double_t& fractionErrorStat,\r
-                             Double_t& fractionErrorSys) const;\r
-  Bool_t GetParticleFractions(Double_t trackPt, Double_t jetPt, Double_t centralityPercentile,\r
-                              Double_t* prob, Int_t smearSpeciesByError, Int_t takeIntoAccountSpeciesSysError,\r
-                              Bool_t uniformSystematicError = kFALSE) const;\r
-  const TH3D* GetParticleFractionHisto(Int_t species, Bool_t sysError = kFALSE) const;\r
-  Bool_t SetParticleFractionHisto(const TH3D* hist, Int_t species, Bool_t sysError = kFALSE);\r
-  Int_t GetParticleFractionHistoNbinsTrackPt() const;\r
-  Int_t GetParticleFractionHistoNbinsJetPt() const;\r
-  Int_t GetParticleFractionHistoNbinsCentrality() const;\r
-  Bool_t SetParticleFractionHistosFromFile(const TString filePathName, Bool_t sysError = kFALSE);\r
-  Int_t GetRandomParticleTypeAccordingToParticleFractions(Double_t trackPt, Double_t jetPt, \r
-                                                          Double_t centralityPercentile,\r
-                                                          Bool_t smearByError,\r
-                                                          Bool_t takeIntoAccountSysError = kFALSE) const;\r
-  \r
-  \r
- protected:\r
-  void CheckDoAnyStematicStudiesOnTheExpectedSignal();\r
-  Double_t ConvolutedGaus(const Double_t* xx, const Double_t* par) const;\r
-  inline Double_t FastGaus(Double_t x, Double_t mean, Double_t sigma) const;\r
-  inline Double_t FastNormalisedGaus(Double_t x, Double_t mean, Double_t sigma) const;\r
-  Int_t FindBinWithinRange(TAxis* axis, Double_t value) const;\r
-  Int_t FindFirstBinAboveIn3dSubset(const TH3* hist, Double_t threshold, Int_t yValue, Int_t zValue) const;\r
-  Int_t FindLastBinAboveIn3dSubset(const TH3* hist, Double_t threshold, Int_t yValue, Int_t zValue) const;\r
-  virtual void SetUpGenHist(THnSparse* hist, Double_t* binsPt, Double_t* binsDeltaPrime, Double_t* binsCent, Double_t* binsJetPt) const;\r
-  virtual void SetUpGenYieldHist(THnSparse* hist, Double_t* binsPt, Double_t* binsCent, Double_t* binsJetPt) const;\r
-  virtual void SetUpHist(THnSparse* hist, Double_t* binsPt, Double_t* binsDeltaPrime, Double_t* binsCent, Double_t* binsJetPt) const;\r
-  virtual void SetUpPtResHist(THnSparse* hist, Double_t* binsPt, Double_t* binsJetPt, Double_t* binsCent) const;\r
-  virtual void SetUpPIDcombined();\r
-  \r
-  static const Int_t fgkNumJetAxes = 3; // Number of additional axes for jets\r
-  static const Double_t fgkEpsilon; // Double_t threshold above zero\r
-  static const Int_t fgkMaxNumGenEntries = 1000; // Maximum number of generated detector responses per track and delta(Prime) and associated species\r
-\r
- private:\r
-  static const Double_t fgkOneOverSqrt2; // = 1. / TMath::Sqrt2();\r
-  \r
-  AliPIDCombined* fPIDcombined; //! PID combined object\r
-  \r
-  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\r
-  \r
-  Bool_t fDoPID; // Only do PID processing (and post the output), if flag is set to kTRUE\r
-  Bool_t fDoEfficiency; // Only do efficiency processing (and post the output), if flag is set to kTRUE\r
-  Bool_t fDoPtResolution; // Only do pT resolution processing (and post the output), if flag is set to kTRUE\r
-  \r
-  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\r
-  Bool_t fStoreAdditionalJetInformation; // If set to kTRUE, additional jet information like jetPt, z, xi will be stored in the THnSparses\r
-\r
-  Bool_t fTakeIntoAccountMuons; // Also take into account muons for the generation of the expected response and the most probable PID\r
-  Bool_t fUseITS; // Use ITS for PID combined probabilities\r
-  Bool_t fUseTOF; // Use TOF for PID combined probabilities\r
-  Bool_t fUsePriors; // Use priors for PID combined probabilities\r
-  Bool_t fTPCDefaultPriors; // Use TPC default priors for PID combined probabilities, if priors are enabled\r
-    \r
-  Bool_t fUseMCidForGeneration; // If MC, use MCid instead of PIDcombined to generate the signals\r
-  \r
-  Bool_t fUseConvolutedGaus; // Use convoluted gaus to generate detector response instead of pure gaus  \r
-  const Int_t fkConvolutedGausNPar; // Number of parameters for convoluted gaus\r
-  Double_t fAccuracyNonGaussianTail; // Accuracy of the non-gaussian tail (fraction of the maximum)\r
-  const Double_t fkDeltaPrimeLowLimit; // Lower deltaPrime limit\r
-  const Double_t fkDeltaPrimeUpLimit; // Upper deltaPrime limit\r
-  TF1* fConvolutedGausDeltaPrime; // Gaus convoluted with exponential tail to generate detector response (deltaPrime)\r
-  \r
-  Double_t fConvolutedGaussTransitionPars[3]; // Parameter for transition from gaussian parameters to asymmetric shape\r
-  static const Double_t fgkSigmaReferenceForTransitionPars; // Reference sigma chosen to calculate transition parameters\r
-  \r
-  Double_t fEtaAbsCutLow; // Lower cut value on |eta|\r
-  Double_t fEtaAbsCutUp;  // Upper cut value on |eta|\r
-  \r
-  // For systematic studies\r
-  Bool_t   fDoAnySystematicStudiesOnTheExpectedSignal; // Internal flag indicating whether any systematic studies are going to be performed\r
-  Double_t fSystematicScalingSplines;        // Systematic scale factor for the splines (1. = no systematics) \r
-  Double_t fSystematicScalingEtaCorrectionMomentumThr;  // Momentum threshold for the systematic scale factor for the eta correction (separates low-p from high-p\r
-  Double_t fSystematicScalingEtaCorrectionLowMomenta;   // Systematic scale factor for the eta correction (1. = no systematics) at low momenta\r
-  Double_t fSystematicScalingEtaCorrectionHighMomenta;  // Systematic scale factor for the eta correction (1. = no systematics) at high momenta\r
-  Double_t fSystematicScalingEtaSigmaPara;   // Systematic scale factor for the parametrisation of the relative signal width (1. = no systematics) \r
-  Double_t fSystematicScalingMultCorrection; // Systematic scale factor for the multiplicity correction (1. = no systematics) \r
-  \r
-  TH3D* fFractionHists[AliPID::kSPECIES]; // 3D histos of particle fraction as function  with trackPt, jetPt (-1 for inclusive spectra), centralityPercentile (-1 for pp)\r
-  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)\r
-  \r
-  TString fCentralityEstimator; // Estimator for the centrality (e.g. V0A, V0M)\r
-  \r
-  THnSparseD* fhPIDdataAll; //! Data histo\r
-  \r
-  // Generated response information\r
-  THnSparseD* fhGenEl; //! Generated response for el\r
-  THnSparseD* fhGenKa; //! Generated response for ka\r
-  THnSparseD* fhGenPi; //! Generated response for pi\r
-  THnSparseD* fhGenMu; //! Generated response for mu\r
-  THnSparseD* fhGenPr; //! Generated response for pr\r
-  \r
-  // Generated responses for a single track\r
-  Double_t* fGenRespElDeltaPrimeEl; //! Generated responses for a single track\r
-  Double_t* fGenRespElDeltaPrimeKa; //! Generated responses for a single track\r
-  Double_t* fGenRespElDeltaPrimePi; //! Generated responses for a single track\r
-  Double_t* fGenRespElDeltaPrimePr; //! Generated responses for a single track\r
-  Double_t* fGenRespKaDeltaPrimeEl; //! Generated responses for a single track\r
-  Double_t* fGenRespKaDeltaPrimeKa; //! Generated responses for a single track\r
-  Double_t* fGenRespKaDeltaPrimePi; //! Generated responses for a single track\r
-  Double_t* fGenRespKaDeltaPrimePr; //! Generated responses for a single track\r
-  Double_t* fGenRespPiDeltaPrimeEl; //! Generated responses for a single track\r
-  Double_t* fGenRespPiDeltaPrimeKa; //! Generated responses for a single track\r
-  Double_t* fGenRespPiDeltaPrimePi; //! Generated responses for a single track\r
-  Double_t* fGenRespPiDeltaPrimePr; //! Generated responses for a single track\r
-  Double_t* fGenRespMuDeltaPrimeEl; //! Generated responses for a single track\r
-  Double_t* fGenRespMuDeltaPrimeKa; //! Generated responses for a single track\r
-  Double_t* fGenRespMuDeltaPrimePi; //! Generated responses for a single track\r
-  Double_t* fGenRespMuDeltaPrimePr; //! Generated responses for a single track\r
-  Double_t* fGenRespPrDeltaPrimeEl; //! Generated responses for a single track\r
-  Double_t* fGenRespPrDeltaPrimeKa; //! Generated responses for a single track\r
-  Double_t* fGenRespPrDeltaPrimePi; //! Generated responses for a single track\r
-  Double_t* fGenRespPrDeltaPrimePr; //! Generated responses for a single track\r
-  /*\r
-  Double_t* fGenRespElDeltaEl; //! Generated responses for a single track\r
-  Double_t* fGenRespElDeltaKa; //! Generated responses for a single track\r
-  Double_t* fGenRespElDeltaPi; //! Generated responses for a single track\r
-  Double_t* fGenRespElDeltaPr; //! Generated responses for a single track\r
-  Double_t* fGenRespKaDeltaEl; //! Generated responses for a single track\r
-  Double_t* fGenRespKaDeltaKa; //! Generated responses for a single track\r
-  Double_t* fGenRespKaDeltaPi; //! Generated responses for a single track\r
-  Double_t* fGenRespKaDeltaPr; //! Generated responses for a single track\r
-  Double_t* fGenRespPiDeltaEl; //! Generated responses for a single track\r
-  Double_t* fGenRespPiDeltaKa; //! Generated responses for a single track\r
-  Double_t* fGenRespPiDeltaPi; //! Generated responses for a single track\r
-  Double_t* fGenRespPiDeltaPr; //! Generated responses for a single track\r
-  Double_t* fGenRespMuDeltaEl; //! Generated responses for a single track\r
-  Double_t* fGenRespMuDeltaKa; //! Generated responses for a single track\r
-  Double_t* fGenRespMuDeltaPi; //! Generated responses for a single track\r
-  Double_t* fGenRespMuDeltaPr; //! Generated responses for a single track\r
-  Double_t* fGenRespPrDeltaEl; //! Generated responses for a single track\r
-  Double_t* fGenRespPrDeltaKa; //! Generated responses for a single track\r
-  Double_t* fGenRespPrDeltaPi; //! Generated responses for a single track\r
-  Double_t* fGenRespPrDeltaPr; //! Generated responses for a single track\r
-  */\r
-  \r
-  TH1D* fhEventsProcessed; //! Histo holding the number of processed events\r
-  TH2D* fhSkippedTracksForSignalGeneration; //! Number of tracks that have been skipped for the signal generation\r
-  THnSparseD* fhMCgeneratedYieldsPrimaries; //! Histo holding the generated (no reco, no cuts) primary particle yields in considered eta range\r
-  \r
-  TH2D* fh2FFJetPtRec;            //! Number of reconstructed jets vs. jetPt and centrality\r
-  TH2D* fh2FFJetPtGen;            //! Number of generated jets vs. jetPt and centrality\r
-  \r
-  TProfile* fh1Xsec;              //! pythia cross section and trials\r
-  TH1D*     fh1Trials;            //! sum of trials\r
-  \r
-  AliCFContainer* fContainerEff; //! Container for efficiency determination\r
-  \r
-  THnSparseD* fPtResolution[AliPID::kSPECIES]; //! Pt Resolution for the individual species\r
-  \r
-  TObjArray* fOutputContainer;  //! output data container\r
-  \r
-  TObjArray* fPtResolutionContainer;  //! output data container for pt resolution\r
-  \r
-  AliAnalysisTaskPID(const AliAnalysisTaskPID&); // not implemented\r
-  AliAnalysisTaskPID& operator=(const AliAnalysisTaskPID&); // not implemented\r
-  \r
-  ClassDef(AliAnalysisTaskPID, 13);\r
-};\r
-\r
-\r
-//_____________________________________________________________________________\r
-inline Bool_t AliAnalysisTaskPID::FillEfficiencyContainer(const Double_t* values, AliAnalysisTaskPID::EffSteps step,\r
-                                                          Double_t weight) \r
-{\r
-  // Fill efficiency container at step "step" with the values\r
-  \r
-  if (!fDoEfficiency)\r
-    return kFALSE;\r
-  \r
-  if (!fContainerEff) {\r
-    AliError("Efficiency container not initialised -> cannot be filled!");\r
-    return kFALSE;\r
-  }\r
-  \r
-  fContainerEff->Fill(values, step, weight);    \r
-  \r
-  return kTRUE;\r
-}\r
-\r
-\r
-//_____________________________________________________________________________\r
-inline Bool_t AliAnalysisTaskPID::FillGeneratedYield(const Double_t* values, Double_t weight)\r
-{\r
-  // Fill histos with generated primary yields with provided values\r
-  \r
-  if (!fDoPID)\r
-    return kFALSE;\r
-  \r
-  if (!fhMCgeneratedYieldsPrimaries) {\r
-    AliError("Histo for generated primary yield not initialised -> cannot be filled!");\r
-    return kFALSE;\r
-  }\r
-  \r
-  fhMCgeneratedYieldsPrimaries->Fill(values, weight);\r
-    \r
-  return kTRUE;\r
-}\r
-\r
-\r
-//_____________________________________________________________________________\r
-inline Bool_t AliAnalysisTaskPID::FillGenJets(Double_t centralityPercentile, Double_t jetPt, Double_t norm)\r
-{\r
-  if (!fDoPID && !fDoEfficiency)\r
-    return kFALSE;\r
-  \r
-  if (!fh2FFJetPtGen)\r
-    return kFALSE;\r
-  \r
-  if (norm > 0.)\r
-    fh2FFJetPtGen->Fill(centralityPercentile, jetPt, 1. / norm);\r
-  else\r
-    fh2FFJetPtGen->Fill(centralityPercentile, jetPt);\r
-  \r
-  return kTRUE;\r
-}\r
-\r
-\r
-//_____________________________________________________________________________\r
-inline Bool_t AliAnalysisTaskPID::FillRecJets(Double_t centralityPercentile, Double_t jetPt, Double_t norm)\r
-{\r
-  if (!fDoPID && !fDoEfficiency)\r
-    return kFALSE;\r
-  \r
-  if (!fh2FFJetPtRec)\r
-    return kFALSE;\r
-  \r
-  if (norm > 0.)\r
-    fh2FFJetPtRec->Fill(centralityPercentile, jetPt, 1. / norm);\r
-  else\r
-    fh2FFJetPtRec->Fill(centralityPercentile, jetPt);\r
-  \r
-  return kTRUE;\r
-}\r
-\r
-\r
-//_____________________________________________________________________________\r
-inline Bool_t AliAnalysisTaskPID::FillPtResolution(Int_t mcID, const Double_t* values)\r
-{\r
-  // Fill histos with pT resolution with provided values\r
-  \r
-  if (!fDoPtResolution || mcID < 0 || mcID >= AliPID::kSPECIES)\r
-    return kFALSE;\r
-  \r
-  if (!fPtResolution[mcID]) {\r
-    AliError(Form("Histo for pT resolution (species: %s) not initialised -> cannot be filled!", AliPID::ParticleName(mcID)));\r
-    return kFALSE;\r
-  }\r
-  \r
-  fPtResolution[mcID]->Fill(values);\r
-    \r
-  return kTRUE;\r
-}\r
\r
-\r
-//_____________________________________________________________________________\r
-inline Bool_t AliAnalysisTaskPID::IncrementEventsProcessed(Double_t centralityPercentile)\r
-{\r
-  // Increment the number of processed events for the given centrality percentile\r
-  \r
-  if (!fDoPID)\r
-    return kFALSE;\r
-  \r
-  if (!fhEventsProcessed) {\r
-    AliError("Histogram for number of events not initialised -> cannot be incremented!");\r
-    return kFALSE;\r
-  }\r
-  \r
-  fhEventsProcessed->Fill(centralityPercentile);\r
-  return kTRUE;\r
-};\r
-\r
-\r
-//_____________________________________________________________________________\r
-inline Bool_t AliAnalysisTaskPID::SetEtaAbsCutRange(Double_t lowerLimit, Double_t upperLimit)\r
-{\r
-  if (lowerLimit >= upperLimit) {\r
-    AliError(Form("Requested lower |eta| cut limit >= upper |eta| cut limit. Old eta cut range will be used (low %f, high %f).",\r
-                  fEtaAbsCutLow, fEtaAbsCutUp));\r
-    return kFALSE;\r
-  }\r
-  \r
-  fEtaAbsCutLow = lowerLimit;\r
-  fEtaAbsCutUp = upperLimit;\r
-  \r
-  return kTRUE;\r
-};\r
-\r
-\r
-//_____________________________________________________________________________\r
-inline Double_t AliAnalysisTaskPID::GetConvolutedGaussTransitionPar(Int_t index) const\r
-{\r
-  if (index < 0 || index >= 3) {\r
-    printf("Invalid index %d!\n", index);\r
-    return -1;\r
-  }\r
-  return fConvolutedGaussTransitionPars[index];\r
-}\r
-\r
-\r
-//_____________________________________________________________________________\r
-inline Int_t AliAnalysisTaskPID::GetParticleFractionHistoNbinsTrackPt() const\r
-{\r
-  if (!fFractionHists[AliPID::kPion])\r
-    return -1;\r
-  \r
-  return fFractionHists[AliPID::kPion]->GetNbinsX();\r
-}\r
-\r
-\r
-//_____________________________________________________________________________\r
-inline Int_t AliAnalysisTaskPID::GetParticleFractionHistoNbinsJetPt() const\r
-{\r
-  if (!fFractionHists[AliPID::kPion])\r
-    return -1;\r
-  \r
-  return fFractionHists[AliPID::kPion]->GetNbinsY();\r
-}\r
-\r
-\r
-//_____________________________________________________________________________\r
-inline Int_t AliAnalysisTaskPID::GetParticleFractionHistoNbinsCentrality() const\r
-{\r
-  if (!fFractionHists[AliPID::kPion])\r
-    return -1;\r
-  \r
-  return fFractionHists[AliPID::kPion]->GetNbinsZ();\r
-}\r
-\r
-\r
-//_____________________________________________________________________________\r
-inline Double_t AliAnalysisTaskPID::GetCentralityPercentile(AliVEvent* evt) const\r
-{\r
-  if (!evt)\r
-    return -1;\r
-  \r
-  Double_t centralityPercentile = -1.;\r
-  if (fStoreCentralityPercentile)\r
-    centralityPercentile = evt->GetCentrality()->GetCentralityPercentile(fCentralityEstimator.Data());\r
-  \r
-  return centralityPercentile;\r
-}\r
-\r
-\r
-//_____________________________________________________________________________\r
-inline void AliAnalysisTaskPID::PostOutputData()\r
-{\r
-  PostData(1, fOutputContainer);\r
-  \r
-  if (fDoEfficiency)\r
-    PostData(2, fContainerEff);\r
-  \r
-  if (fDoPtResolution)\r
-    PostData(3, fPtResolutionContainer);\r
-}\r
-\r
-#endif\r
+#ifndef ALI_ANALYSIS_TASK_PID_H
+#define ALI_ANALYSIS_TASK_PID_H
+
+/*
+This task collects PID output from different detectors.
+Only tracks fulfilling some standard quality cuts are taken into account.
+At the moment, only data from TPC and TOF is collected. But in future,
+data from e.g. HMPID is also foreseen.
+
+Class written by Benjamin Hess.
+Contact: bhess@cern.ch
+*/
+
+class TF1;
+class TRandom3;
+class AliAnalysisFilter;
+class AliCFContainer;
+class AliESDEvent;
+class AliMCEvent;
+class AliMCParticle;
+class AliPID;
+class AliPIDCombined;
+class AliPIDResponse;
+class AliTOFPIDResponse;
+class AliVEvent;
+class AliVTrack;
+
+#include "TH1D.h"
+#include "TH2D.h"
+#include "TH3D.h"
+#include "THnSparse.h"
+#include "TProfile.h"
+#include "TString.h"
+
+#include "AliCentrality.h"
+#include "AliCFContainer.h"
+
+#include "AliPID.h"
+#include "AliAnalysisTaskPIDV0base.h"
+
+class AliAnalysisTaskPID : public AliAnalysisTaskPIDV0base {
+ public:
+  AliAnalysisTaskPID();
+  AliAnalysisTaskPID(const char *name);
+  virtual ~AliAnalysisTaskPID();
+  
+  virtual void   UserCreateOutputObjects();
+  virtual void   UserExec(Option_t *option);
+  virtual void   Terminate(const Option_t*);
+  
+  enum ErrorCode { kNoErrors = 1, kWarning = 0, kError = -1};
+  
+  enum dataAxes { kDataMCID = 0, kDataSelectSpecies = 1, kDataPt = 2, kDataDeltaPrimeSpecies = 3, kDataCentrality = 4,
+                  kDataJetPt = 5, kDataZ = 6, kDataXi = 7, kDataCharge = 8, kDataTOFpidInfo = 9, kDataNumAxes = 10 };
+
+  enum genAxes { kGenMCID = 0, kGenSelectSpecies = 1, kGenPt = 2, kGenDeltaPrimeSpecies = 3, kGenCentrality = 4,
+                 kGenJetPt = 5, kGenZ = 6, kGenXi = 7, kGenCharge = 8, kGenTOFpidInfo = 9, kGenNumAxes = 10 };
+  
+  enum genYieldAxes { kGenYieldMCID = 0, kGenYieldPt = 1, kGenYieldCentrality = 2, kGenYieldJetPt = 3, kGenYieldZ = 4, kGenYieldXi = 5,
+                      kGenYieldCharge = 6, kGenYieldNumAxes = 7 };
+  
+  enum ptResolutionAxes { kPtResJetPt = 0, kPtResGenPt = 1, kPtResRecPt = 2, kPtResCharge = 3, kPtResCentrality = 4, kPtResNumAxes = 5 };
+  
+  enum efficiencyAxes { kEffMCID = 0, kEffTrackPt = 1, kEffTrackEta = 2, kEffTrackCharge = 3, kEffCentrality = 4, kEffJetPt = 5,
+                        kEffZ = 6, kEffXi = 7, kEffNumAxes = 8 };
+  
+  enum EffSteps { kStepGenWithGenCuts = 0, kStepRecWithGenCuts = 1, kStepRecWithGenCutsMeasuredObs = 2,
+                  kStepRecWithRecCutsMeasuredObs = 3, kStepRecWithRecCutsMeasuredObsPrimaries = 4,
+                  kStepRecWithRecCutsMeasuredObsStrangenessScaled = 5, kStepRecWithRecCutsPrimaries = 6, kNumSteps = 7};
+  
+  enum TOFpidInfo { kNoTOFinfo = -2, kNoTOFpid = -1, kTOFpion = 0, kTOFkaon = 1, kTOFproton = 2, kNumTOFspecies = 3,
+                    kNumTOFpidInfoBins = 5 };
+  
+  static Int_t PDGtoMCID(Int_t pdg);
+  
+  static void GetJetTrackObservables(Double_t trackPt, Double_t jetPt, Double_t& z, Double_t& xi);
+  
+  static Double_t GetMCStrangenessFactorCMS(Int_t motherPDG, Double_t motherGenPt);
+  static Double_t GetMCStrangenessFactorCMS(AliMCEvent* mcEvent, AliMCParticle* daughter);
+  
+  static Bool_t IsSecondaryWithStrangeMotherMC(AliMCEvent* mcEvent, Int_t partLabel);
+  
+  Int_t GetIndexOfChargeAxisData() const
+    { return fStoreAdditionalJetInformation ? kDataCharge : kDataCharge - fgkNumJetAxes; };
+  Int_t GetIndexOfChargeAxisGen() const
+    { return fStoreAdditionalJetInformation ? kGenCharge : kGenCharge - fgkNumJetAxes; };
+  Int_t GetIndexOfChargeAxisGenYield() const
+    { return fStoreAdditionalJetInformation ? kGenYieldCharge : kGenYieldCharge - fgkNumJetAxes; };
+  
+  Int_t GetIndexOfTOFpidInfoAxisData() const
+    { return fStoreAdditionalJetInformation ? kDataTOFpidInfo : kDataTOFpidInfo - fgkNumJetAxes; };
+  Int_t GetIndexOfTOFpidInfoAxisGen() const
+    { return fStoreAdditionalJetInformation ? kGenTOFpidInfo : kGenTOFpidInfo - fgkNumJetAxes; };
+  
+  Bool_t FillXsec(Double_t xsection)
+    { if (!fh1Xsec) return kFALSE; fh1Xsec->Fill("<#sigma>", xsection); return kTRUE; };
+  Bool_t FillPythiaTrials(Double_t avgTrials)
+    { if (!fh1Trials) return kFALSE; fh1Trials->Fill("#sum{ntrials}", avgTrials); return kTRUE; };
+    
+  Bool_t FillEfficiencyContainer(const Double_t* values, EffSteps step, Double_t weight = 1.0);
+  
+  Bool_t FillGeneratedYield(const Double_t* values, Double_t weight = 1.0);
+  Bool_t FillPtResolution(Int_t mcID, const Double_t* values);
+  Bool_t FillGenJets(Double_t centralityPercentile, Double_t jetPt, Double_t norm = -1.);
+  Bool_t FillRecJets(Double_t centralityPercentile, Double_t jetPt, Double_t norm = -1.);
+  
+  Bool_t IncrementEventsProcessed(Double_t centralityPercentile);
+  
+  void PostOutputData();
+  
+  void PrintSettings(Bool_t printSystematicsSettings = kFALSE) const;
+  
+  void PrintSystematicsSettings() const;
+  
+  Bool_t ProcessTrack(const AliVTrack* track, Int_t particlePDGcode, Double_t centralityPercentile, Double_t jetPt) ;
+  
+  ErrorCode GenerateDetectorResponse(ErrorCode errCode, Double_t mean, Double_t sigma, Double_t* responses,
+                                     Int_t nResponses,
+                                     Bool_t usePureGaus = kFALSE);
+  ErrorCode SetParamsForConvolutedGaus(Double_t gausMean, Double_t gausSigma);
+  
+  const TString GetCentralityEstimator() const { return fCentralityEstimator; };
+  void SetCentralityEstimator(TString estimator) { fCentralityEstimator = estimator; };
+  
+  Double_t GetCentralityPercentile(AliVEvent* evt) const;
+  
+  inline Double_t GetConvolutedGaussTransitionPar(Int_t index) const;
+  
+  Bool_t SetConvolutedGaussLambdaParameter(Double_t lambda);
+
+  Bool_t GetInputFromOtherTask() const { return fInputFromOtherTask; };
+  void SetInputFromOtherTask(Bool_t flag) { fInputFromOtherTask = flag; };
+  
+  Bool_t GetDoPID() const { return fDoPID; };
+  void SetDoPID(Bool_t flag) { fDoPID = flag; };
+  
+  Bool_t GetDoEfficiency() const { return fDoEfficiency; };
+  void SetDoEfficiency(Bool_t flag) { fDoEfficiency = flag; };
+  
+  Bool_t GetDoPtResolution() const { return fDoPtResolution; };
+  void SetDoPtResolution(Bool_t flag) { fDoPtResolution = flag; };
+  
+  Bool_t GetStoreCentralityPercentile() const { return fStoreCentralityPercentile; };
+  void SetStoreCentralityPercentile(Bool_t flag) { fStoreCentralityPercentile = flag; };
+  
+  Bool_t GetStoreAdditionalJetInformation() const { return fStoreAdditionalJetInformation; };
+  void SetStoreAdditionalJetInformation(Bool_t flag) { fStoreAdditionalJetInformation = flag; };
+  
+  Bool_t GetUseMCidForGeneration() const { return fUseMCidForGeneration; };
+  void SetUseMCidForGeneration(Bool_t flag) { fUseMCidForGeneration = flag; };
+  
+  Bool_t GetUseConvolutedGaus() const { return fUseConvolutedGaus; };
+  void SetUseConvolutedGaus(Bool_t flag) { fUseConvolutedGaus = flag; };
+  
+  Double_t GetAccuracyNonGaussianTail() const { return fAccuracyNonGaussianTail; };
+  void SetAccuracyNonGaussianTail(Double_t value) { fAccuracyNonGaussianTail = value; };
+  
+  Bool_t GetTakeIntoAccountMuons() const { return fTakeIntoAccountMuons; };
+  void SetTakeIntoAccountMuons(Bool_t flag) { fTakeIntoAccountMuons = flag; };
+  
+  Int_t GetTOFmode() const { return fTOFmode; };
+  void SetTOFmode(Int_t tofMode) { fTOFmode = tofMode; };
+  
+  Bool_t GetUseTPCDefaultPriors() const { return fTPCDefaultPriors; };
+  void SetUseTPCDefaultPriors(Bool_t flag) { fTPCDefaultPriors = flag; };
+  
+  Bool_t GetUsePriors() const { return fUsePriors; };
+  void SetUsePriors(Bool_t flag) { fUsePriors = flag; };
+  
+  Bool_t GetUseITS() const { return fUseITS; };
+  void SetUseITS(Bool_t flag) { fUseITS = flag; };
+  
+  Bool_t GetUseTOF() const { return fUseTOF; };
+  void SetUseTOF(Bool_t flag) { fUseTOF = flag; };
+  
+  Double_t GetEtaAbsCutLow() const { return fEtaAbsCutLow; };
+  Double_t GetEtaAbsCutUp() const { return fEtaAbsCutUp; };
+  Bool_t SetEtaAbsCutRange(Double_t lowerLimit, Double_t upperLimit);
+  
+  Bool_t IsInAcceptedEtaRange(Double_t etaAbs) const { return (etaAbs >= fEtaAbsCutLow && etaAbs <= fEtaAbsCutUp); };
+  
+  Double_t GetSystematicScalingSplines() const { return fSystematicScalingSplines; };
+  void SetSystematicScalingSplines(Double_t scaleFactor) 
+    { fSystematicScalingSplines = scaleFactor; CheckDoAnyStematicStudiesOnTheExpectedSignal(); };
+  
+  Double_t GetSystematicScalingEtaCorrectionMomentumThr() const { return fSystematicScalingEtaCorrectionMomentumThr; };
+  void SetSystematicScalingEtaCorrectionMomentumThr(Double_t threshold) { fSystematicScalingEtaCorrectionMomentumThr = threshold; };
+  
+  Double_t GetSystematicScalingEtaCorrectionLowMomenta() const { return fSystematicScalingEtaCorrectionLowMomenta; };
+  void SetSystematicScalingEtaCorrectionLowMomenta(Double_t scaleFactor) 
+    { fSystematicScalingEtaCorrectionLowMomenta = scaleFactor; CheckDoAnyStematicStudiesOnTheExpectedSignal(); };
+  
+  Double_t GetSystematicScalingEtaCorrectionHighMomenta() const { return fSystematicScalingEtaCorrectionHighMomenta; };
+  void SetSystematicScalingEtaCorrectionHighMomenta(Double_t scaleFactor) 
+    { fSystematicScalingEtaCorrectionHighMomenta = scaleFactor; CheckDoAnyStematicStudiesOnTheExpectedSignal(); };
+  
+  Double_t GetSystematicScalingEtaSigmaPara() const { return fSystematicScalingEtaSigmaPara; };
+  void SetSystematicScalingEtaSigmaPara(Double_t scaleFactor)
+    { fSystematicScalingEtaSigmaPara = scaleFactor; CheckDoAnyStematicStudiesOnTheExpectedSignal(); };
+  
+  Double_t GetSystematicScalingMultCorrection() const { return fSystematicScalingMultCorrection; };
+  void SetSystematicScalingMultCorrection(Double_t scaleFactor) 
+    { fSystematicScalingMultCorrection = scaleFactor; CheckDoAnyStematicStudiesOnTheExpectedSignal(); };
+  
+  void CleanupParticleFractionHistos();
+  Bool_t GetParticleFraction(Double_t trackPt, Double_t jetPt, Double_t multiplicity,
+                             AliPID::EParticleType species, Double_t& fraction, Double_t& fractionErrorStat,
+                             Double_t& fractionErrorSys) const;
+  Bool_t GetParticleFractions(Double_t trackPt, Double_t jetPt, Double_t centralityPercentile,
+                              Double_t* prob, Int_t smearSpeciesByError, Int_t takeIntoAccountSpeciesSysError,
+                              Bool_t uniformSystematicError = kFALSE) const;
+  const TH3D* GetParticleFractionHisto(Int_t species, Bool_t sysError = kFALSE) const;
+  Bool_t SetParticleFractionHisto(const TH3D* hist, Int_t species, Bool_t sysError = kFALSE);
+  Int_t GetParticleFractionHistoNbinsTrackPt() const;
+  Int_t GetParticleFractionHistoNbinsJetPt() const;
+  Int_t GetParticleFractionHistoNbinsCentrality() const;
+  Bool_t SetParticleFractionHistosFromFile(const TString filePathName, Bool_t sysError = kFALSE);
+  Int_t GetRandomParticleTypeAccordingToParticleFractions(Double_t trackPt, Double_t jetPt, 
+                                                          Double_t centralityPercentile,
+                                                          Bool_t smearByError,
+                                                          Bool_t takeIntoAccountSysError = kFALSE) const;
+  
+  TOFpidInfo GetTOFType(const AliVTrack* track, Int_t tofMode) const;
+  
+ protected:
+  void CheckDoAnyStematicStudiesOnTheExpectedSignal();
+  Double_t ConvolutedGaus(const Double_t* xx, const Double_t* par) const;
+  inline Double_t FastGaus(Double_t x, Double_t mean, Double_t sigma) const;
+  inline Double_t FastNormalisedGaus(Double_t x, Double_t mean, Double_t sigma) const;
+  Int_t FindBinWithinRange(TAxis* axis, Double_t value) const;
+  Int_t FindFirstBinAboveIn3dSubset(const TH3* hist, Double_t threshold, Int_t yValue, Int_t zValue) const;
+  Int_t FindLastBinAboveIn3dSubset(const TH3* hist, Double_t threshold, Int_t yValue, Int_t zValue) const;
+  virtual void SetUpGenHist(THnSparse* hist, Double_t* binsPt, Double_t* binsDeltaPrime, Double_t* binsCent, Double_t* binsJetPt) const;
+  virtual void SetUpGenYieldHist(THnSparse* hist, Double_t* binsPt, Double_t* binsCent, Double_t* binsJetPt) const;
+  virtual void SetUpHist(THnSparse* hist, Double_t* binsPt, Double_t* binsDeltaPrime, Double_t* binsCent, Double_t* binsJetPt) const;
+  virtual void SetUpPtResHist(THnSparse* hist, Double_t* binsPt, Double_t* binsJetPt, Double_t* binsCent) const;
+  virtual void SetUpPIDcombined();
+  
+  static const Int_t fgkNumJetAxes; // Number of additional axes for jets
+  static const Double_t fgkEpsilon; // Double_t threshold above zero
+  static const Int_t fgkMaxNumGenEntries; // Maximum number of generated detector responses per track and delta(Prime) and associated species
+
+ private:
+  static const Double_t fgkOneOverSqrt2; // = 1. / TMath::Sqrt2();
+  
+  AliPIDCombined* fPIDcombined; //! PID combined object
+  
+  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
+  
+  Bool_t fDoPID; // Only do PID processing (and post the output), if flag is set to kTRUE
+  Bool_t fDoEfficiency; // Only do efficiency processing (and post the output), if flag is set to kTRUE
+  Bool_t fDoPtResolution; // Only do pT resolution processing (and post the output), if flag is set to kTRUE
+  
+  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
+  Bool_t fStoreAdditionalJetInformation; // If set to kTRUE, additional jet information like jetPt, z, xi will be stored in the THnSparses
+
+  Bool_t fTakeIntoAccountMuons; // Also take into account muons for the generation of the expected response and the most probable PID
+  Bool_t fUseITS; // Use ITS for PID combined probabilities
+  Bool_t fUseTOF; // Use TOF for PID combined probabilities
+  Bool_t fUsePriors; // Use priors for PID combined probabilities
+  Bool_t fTPCDefaultPriors; // Use TPC default priors for PID combined probabilities, if priors are enabled
+    
+  Bool_t fUseMCidForGeneration; // If MC, use MCid instead of PIDcombined to generate the signals
+  
+  Bool_t fUseConvolutedGaus; // Use convoluted gaus to generate detector response instead of pure gaus  
+  const Int_t fkConvolutedGausNPar; // Number of parameters for convoluted gaus
+  Double_t fAccuracyNonGaussianTail; // Accuracy of the non-gaussian tail (fraction of the maximum)
+  const Double_t fkDeltaPrimeLowLimit; // Lower deltaPrime limit
+  const Double_t fkDeltaPrimeUpLimit; // Upper deltaPrime limit
+  TF1* fConvolutedGausDeltaPrime; // Gaus convoluted with exponential tail to generate detector response (deltaPrime)
+  
+  Int_t fTOFmode; // TOF mode used for TOF PID info (affects num sigma inclusion/exclusion)
+  Double_t fConvolutedGaussTransitionPars[3]; // Parameter for transition from gaussian parameters to asymmetric shape
+  static const Double_t fgkSigmaReferenceForTransitionPars; // Reference sigma chosen to calculate transition parameters
+  
+  Double_t fEtaAbsCutLow; // Lower cut value on |eta|
+  Double_t fEtaAbsCutUp;  // Upper cut value on |eta|
+  
+  // For systematic studies
+  Bool_t   fDoAnySystematicStudiesOnTheExpectedSignal; // Internal flag indicating whether any systematic studies are going to be performed
+  Double_t fSystematicScalingSplines;        // Systematic scale factor for the splines (1. = no systematics) 
+  Double_t fSystematicScalingEtaCorrectionMomentumThr;  // Momentum threshold for the systematic scale factor for the eta correction (separates low-p from high-p
+  Double_t fSystematicScalingEtaCorrectionLowMomenta;   // Systematic scale factor for the eta correction (1. = no systematics) at low momenta
+  Double_t fSystematicScalingEtaCorrectionHighMomenta;  // Systematic scale factor for the eta correction (1. = no systematics) at high momenta
+  Double_t fSystematicScalingEtaSigmaPara;   // Systematic scale factor for the parametrisation of the relative signal width (1. = no systematics) 
+  Double_t fSystematicScalingMultCorrection; // Systematic scale factor for the multiplicity correction (1. = no systematics) 
+  
+  TH3D* fFractionHists[AliPID::kSPECIES]; // 3D histos of particle fraction as function  with trackPt, jetPt (-1 for inclusive spectra), centralityPercentile (-1 for pp)
+  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)
+  
+  TString fCentralityEstimator; // Estimator for the centrality (e.g. V0A, V0M)
+  
+  THnSparseD* fhPIDdataAll; //! Data histo
+  
+  // Generated response information
+  THnSparseD* fhGenEl; //! Generated response for el
+  THnSparseD* fhGenKa; //! Generated response for ka
+  THnSparseD* fhGenPi; //! Generated response for pi
+  THnSparseD* fhGenMu; //! Generated response for mu
+  THnSparseD* fhGenPr; //! Generated response for pr
+  
+  // Generated responses for a single track
+  Double_t* fGenRespElDeltaPrimeEl; //! Generated responses for a single track
+  Double_t* fGenRespElDeltaPrimeKa; //! Generated responses for a single track
+  Double_t* fGenRespElDeltaPrimePi; //! Generated responses for a single track
+  Double_t* fGenRespElDeltaPrimePr; //! Generated responses for a single track
+  Double_t* fGenRespKaDeltaPrimeEl; //! Generated responses for a single track
+  Double_t* fGenRespKaDeltaPrimeKa; //! Generated responses for a single track
+  Double_t* fGenRespKaDeltaPrimePi; //! Generated responses for a single track
+  Double_t* fGenRespKaDeltaPrimePr; //! Generated responses for a single track
+  Double_t* fGenRespPiDeltaPrimeEl; //! Generated responses for a single track
+  Double_t* fGenRespPiDeltaPrimeKa; //! Generated responses for a single track
+  Double_t* fGenRespPiDeltaPrimePi; //! Generated responses for a single track
+  Double_t* fGenRespPiDeltaPrimePr; //! Generated responses for a single track
+  Double_t* fGenRespMuDeltaPrimeEl; //! Generated responses for a single track
+  Double_t* fGenRespMuDeltaPrimeKa; //! Generated responses for a single track
+  Double_t* fGenRespMuDeltaPrimePi; //! Generated responses for a single track
+  Double_t* fGenRespMuDeltaPrimePr; //! Generated responses for a single track
+  Double_t* fGenRespPrDeltaPrimeEl; //! Generated responses for a single track
+  Double_t* fGenRespPrDeltaPrimeKa; //! Generated responses for a single track
+  Double_t* fGenRespPrDeltaPrimePi; //! Generated responses for a single track
+  Double_t* fGenRespPrDeltaPrimePr; //! Generated responses for a single track
+  /*
+  Double_t* fGenRespElDeltaEl; //! Generated responses for a single track
+  Double_t* fGenRespElDeltaKa; //! Generated responses for a single track
+  Double_t* fGenRespElDeltaPi; //! Generated responses for a single track
+  Double_t* fGenRespElDeltaPr; //! Generated responses for a single track
+  Double_t* fGenRespKaDeltaEl; //! Generated responses for a single track
+  Double_t* fGenRespKaDeltaKa; //! Generated responses for a single track
+  Double_t* fGenRespKaDeltaPi; //! Generated responses for a single track
+  Double_t* fGenRespKaDeltaPr; //! Generated responses for a single track
+  Double_t* fGenRespPiDeltaEl; //! Generated responses for a single track
+  Double_t* fGenRespPiDeltaKa; //! Generated responses for a single track
+  Double_t* fGenRespPiDeltaPi; //! Generated responses for a single track
+  Double_t* fGenRespPiDeltaPr; //! Generated responses for a single track
+  Double_t* fGenRespMuDeltaEl; //! Generated responses for a single track
+  Double_t* fGenRespMuDeltaKa; //! Generated responses for a single track
+  Double_t* fGenRespMuDeltaPi; //! Generated responses for a single track
+  Double_t* fGenRespMuDeltaPr; //! Generated responses for a single track
+  Double_t* fGenRespPrDeltaEl; //! Generated responses for a single track
+  Double_t* fGenRespPrDeltaKa; //! Generated responses for a single track
+  Double_t* fGenRespPrDeltaPi; //! Generated responses for a single track
+  Double_t* fGenRespPrDeltaPr; //! Generated responses for a single track
+  */
+  
+  TH1D* fhEventsProcessed; //! Histo holding the number of processed events
+  TH2D* fhSkippedTracksForSignalGeneration; //! Number of tracks that have been skipped for the signal generation
+  THnSparseD* fhMCgeneratedYieldsPrimaries; //! Histo holding the generated (no reco, no cuts) primary particle yields in considered eta range
+  
+  TH2D* fh2FFJetPtRec;            //! Number of reconstructed jets vs. jetPt and centrality
+  TH2D* fh2FFJetPtGen;            //! Number of generated jets vs. jetPt and centrality
+  
+  TProfile* fh1Xsec;              //! pythia cross section and trials
+  TH1D*     fh1Trials;            //! sum of trials
+  
+  AliCFContainer* fContainerEff; //! Container for efficiency determination
+  
+  THnSparseD* fPtResolution[AliPID::kSPECIES]; //! Pt Resolution for the individual species
+  
+  TObjArray* fOutputContainer;  //! output data container
+  
+  TObjArray* fPtResolutionContainer;  //! output data container for pt resolution
+  
+  AliAnalysisTaskPID(const AliAnalysisTaskPID&); // not implemented
+  AliAnalysisTaskPID& operator=(const AliAnalysisTaskPID&); // not implemented
+  
+  ClassDef(AliAnalysisTaskPID, 14);
+};
+
+
+//_____________________________________________________________________________
+inline Bool_t AliAnalysisTaskPID::FillEfficiencyContainer(const Double_t* values, AliAnalysisTaskPID::EffSteps step,
+                                                          Double_t weight) 
+{
+  // Fill efficiency container at step "step" with the values
+  
+  if (!fDoEfficiency)
+    return kFALSE;
+  
+  if (!fContainerEff) {
+    AliError("Efficiency container not initialised -> cannot be filled!");
+    return kFALSE;
+  }
+  
+  fContainerEff->Fill(values, step, weight);    
+  
+  return kTRUE;
+}
+
+
+//_____________________________________________________________________________
+inline Bool_t AliAnalysisTaskPID::FillGeneratedYield(const Double_t* values, Double_t weight)
+{
+  // Fill histos with generated primary yields with provided values
+  
+  if (!fDoPID)
+    return kFALSE;
+  
+  if (!fhMCgeneratedYieldsPrimaries) {
+    AliError("Histo for generated primary yield not initialised -> cannot be filled!");
+    return kFALSE;
+  }
+  
+  fhMCgeneratedYieldsPrimaries->Fill(values, weight);
+    
+  return kTRUE;
+}
+
+
+//_____________________________________________________________________________
+inline Bool_t AliAnalysisTaskPID::FillGenJets(Double_t centralityPercentile, Double_t jetPt, Double_t norm)
+{
+  if (!fDoPID && !fDoEfficiency)
+    return kFALSE;
+  
+  if (!fh2FFJetPtGen)
+    return kFALSE;
+  
+  if (norm > 0.)
+    fh2FFJetPtGen->Fill(centralityPercentile, jetPt, 1. / norm);
+  else
+    fh2FFJetPtGen->Fill(centralityPercentile, jetPt);
+  
+  return kTRUE;
+}
+
+
+//_____________________________________________________________________________
+inline Bool_t AliAnalysisTaskPID::FillRecJets(Double_t centralityPercentile, Double_t jetPt, Double_t norm)
+{
+  if (!fDoPID && !fDoEfficiency)
+    return kFALSE;
+  
+  if (!fh2FFJetPtRec)
+    return kFALSE;
+  
+  if (norm > 0.)
+    fh2FFJetPtRec->Fill(centralityPercentile, jetPt, 1. / norm);
+  else
+    fh2FFJetPtRec->Fill(centralityPercentile, jetPt);
+  
+  return kTRUE;
+}
+
+
+//_____________________________________________________________________________
+inline Bool_t AliAnalysisTaskPID::FillPtResolution(Int_t mcID, const Double_t* values)
+{
+  // Fill histos with pT resolution with provided values
+  
+  if (!fDoPtResolution || mcID < 0 || mcID >= AliPID::kSPECIES)
+    return kFALSE;
+  
+  if (!fPtResolution[mcID]) {
+    AliError(Form("Histo for pT resolution (species: %s) not initialised -> cannot be filled!", AliPID::ParticleName(mcID)));
+    return kFALSE;
+  }
+  
+  fPtResolution[mcID]->Fill(values);
+    
+  return kTRUE;
+}
+
+//_____________________________________________________________________________
+inline Bool_t AliAnalysisTaskPID::IncrementEventsProcessed(Double_t centralityPercentile)
+{
+  // Increment the number of processed events for the given centrality percentile
+  
+  if (!fDoPID)
+    return kFALSE;
+  
+  if (!fhEventsProcessed) {
+    AliError("Histogram for number of events not initialised -> cannot be incremented!");
+    return kFALSE;
+  }
+  
+  fhEventsProcessed->Fill(centralityPercentile);
+  return kTRUE;
+};
+
+
+//_____________________________________________________________________________
+inline Bool_t AliAnalysisTaskPID::SetEtaAbsCutRange(Double_t lowerLimit, Double_t upperLimit)
+{
+  if (lowerLimit >= upperLimit) {
+    AliError(Form("Requested lower |eta| cut limit >= upper |eta| cut limit. Old eta cut range will be used (low %f, high %f).",
+                  fEtaAbsCutLow, fEtaAbsCutUp));
+    return kFALSE;
+  }
+  
+  fEtaAbsCutLow = lowerLimit;
+  fEtaAbsCutUp = upperLimit;
+  
+  return kTRUE;
+};
+
+
+//_____________________________________________________________________________
+inline Double_t AliAnalysisTaskPID::GetConvolutedGaussTransitionPar(Int_t index) const
+{
+  if (index < 0 || index >= 3) {
+    printf("Invalid index %d!\n", index);
+    return -1;
+  }
+  return fConvolutedGaussTransitionPars[index];
+}
+
+
+//_____________________________________________________________________________
+inline Int_t AliAnalysisTaskPID::GetParticleFractionHistoNbinsTrackPt() const
+{
+  if (!fFractionHists[AliPID::kPion])
+    return -1;
+  
+  return fFractionHists[AliPID::kPion]->GetNbinsX();
+}
+
+
+//_____________________________________________________________________________
+inline Int_t AliAnalysisTaskPID::GetParticleFractionHistoNbinsJetPt() const
+{
+  if (!fFractionHists[AliPID::kPion])
+    return -1;
+  
+  return fFractionHists[AliPID::kPion]->GetNbinsY();
+}
+
+
+//_____________________________________________________________________________
+inline Int_t AliAnalysisTaskPID::GetParticleFractionHistoNbinsCentrality() const
+{
+  if (!fFractionHists[AliPID::kPion])
+    return -1;
+  
+  return fFractionHists[AliPID::kPion]->GetNbinsZ();
+}
+
+
+//_____________________________________________________________________________
+inline Double_t AliAnalysisTaskPID::GetCentralityPercentile(AliVEvent* evt) const
+{
+  if (!evt)
+    return -1;
+  
+  Double_t centralityPercentile = -1.;
+  if (fStoreCentralityPercentile)
+    centralityPercentile = evt->GetCentrality()->GetCentralityPercentile(fCentralityEstimator.Data());
+  
+  return centralityPercentile;
+}
+
+
+//_____________________________________________________________________________
+inline void AliAnalysisTaskPID::PostOutputData()
+{
+  PostData(1, fOutputContainer);
+  
+  if (fDoEfficiency)
+    PostData(2, fContainerEff);
+  
+  if (fDoPtResolution)
+    PostData(3, fPtResolutionContainer);
+}
+
+#endif