#include "AliESDEvent.h"
#include "AliMCEvent.h"
-#include "AliESDInputHandler.h"
+#include "AliESDtrackCuts.h"
+#include "AliAnalysisFilter.h"
#include "AliInputEventHandler.h"
#include "AliVVertex.h"
-#include "AliAnalysisFilter.h"
#include "AliPID.h"
#include "AliPIDCombined.h"
#include "AliPIDResponse.h"
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
+
//________________________________________________________________________
AliAnalysisTaskPID::AliAnalysisTaskPID()
: AliAnalysisTaskPIDV0base()
+ , fRun(-1)
, fPIDcombined(new AliPIDCombined())
, fInputFromOtherTask(kFALSE)
, fDoPID(kTRUE)
, fDoEfficiency(kTRUE)
- , fDoPtResolution(kTRUE)
+ , fDoPtResolution(kFALSE)
+ , fDoDeDxCheck(kFALSE)
+ , fDoBinZeroStudy(kFALSE)
, fStoreCentralityPercentile(kFALSE)
, fStoreAdditionalJetInformation(kFALSE)
, fTakeIntoAccountMuons(kFALSE)
, fkDeltaPrimeLowLimit(0.02)
, fkDeltaPrimeUpLimit(40.0)
, fConvolutedGausDeltaPrime(0x0)
+ , fTOFmode(1)
, fEtaAbsCutLow(0.0)
, fEtaAbsCutUp(0.9)
+ , fPileUpRejectionType(AliAnalysisTaskPIDV0base::kPileUpRejectionOff)
, fDoAnySystematicStudiesOnTheExpectedSignal(kFALSE)
- , fSystematicScalingSplines(1.0)
+ , fSystematicScalingSplinesThreshold(50.)
+ , fSystematicScalingSplinesBelowThreshold(1.0)
+ , fSystematicScalingSplinesAboveThreshold(1.0)
, fSystematicScalingEtaCorrectionMomentumThr(0.35)
, fSystematicScalingEtaCorrectionLowMomenta(1.0)
, fSystematicScalingEtaCorrectionHighMomenta(1.0)
- , fSystematicScalingEtaSigmaPara(1.0)
+ , fSystematicScalingEtaSigmaParaThreshold(250.)
+ , fSystematicScalingEtaSigmaParaBelowThreshold(1.0)
+ , fSystematicScalingEtaSigmaParaAboveThreshold(1.0)
, fSystematicScalingMultCorrection(1.0)
, fCentralityEstimator("V0M")
, fhPIDdataAll(0x0)
, fGenRespPrDeltaPi(new Double_t[fgkMaxNumGenEntries])
, fGenRespPrDeltaPr(new Double_t[fgkMaxNumGenEntries])
*/
+ , fDeltaPrimeAxis(0x0)
+ , fhMaxEtaVariation(0x0)
, fhEventsProcessed(0x0)
- , fhSkippedTracksForSignalGeneration(0x0)
+ , fhEventsTriggerSel(0x0)
+ , fhEventsTriggerSelVtxCut(0x0)
+ , fhEventsProcessedNoPileUpRejection(0x0)
+ , fChargedGenPrimariesTriggerSel(0x0)
+ , fChargedGenPrimariesTriggerSelVtxCut(0x0)
+ , fChargedGenPrimariesTriggerSelVtxCutZ(0x0)
+ , fChargedGenPrimariesTriggerSelVtxCutZPileUpRej(0x0)
, fhMCgeneratedYieldsPrimaries(0x0)
, fh2FFJetPtRec(0x0)
, fh2FFJetPtGen(0x0)
, fh1Xsec(0x0)
, fh1Trials(0x0)
+ , fh1EvtsPtHardCut(0x0)
, fContainerEff(0x0)
+ , fQASharedCls(0x0)
+ , fDeDxCheck(0x0)
, fOutputContainer(0x0)
- , fPtResolutionContainer(0x0)
+ , fQAContainer(0x0)
{
// default Constructor
fkDeltaPrimeLowLimit, fkDeltaPrimeUpLimit,
fkConvolutedGausNPar, "AliAnalysisTaskPID", "ConvolutedGaus");
+ // Set some arbitrary parameteres, such that the function call will not crash
+ // (although it should not be called with these parameters...)
+ fConvolutedGausDeltaPrime->SetParameter(0, 0);
+ fConvolutedGausDeltaPrime->SetParameter(1, 1);
+ fConvolutedGausDeltaPrime->SetParameter(2, 2);
+
+
// Initialisation of translation parameters is time consuming.
// Therefore, default values will only be initialised if they are really needed.
// Otherwise, it is left to the user to set the parameter properly.
//________________________________________________________________________
AliAnalysisTaskPID::AliAnalysisTaskPID(const char *name)
: AliAnalysisTaskPIDV0base(name)
+ , fRun(-1)
, fPIDcombined(new AliPIDCombined())
, fInputFromOtherTask(kFALSE)
, fDoPID(kTRUE)
, fDoEfficiency(kTRUE)
- , fDoPtResolution(kTRUE)
+ , fDoPtResolution(kFALSE)
+ , fDoDeDxCheck(kFALSE)
+ , fDoBinZeroStudy(kFALSE)
, fStoreCentralityPercentile(kFALSE)
, fStoreAdditionalJetInformation(kFALSE)
, fTakeIntoAccountMuons(kFALSE)
, fkDeltaPrimeLowLimit(0.02)
, fkDeltaPrimeUpLimit(40.0)
, fConvolutedGausDeltaPrime(0x0)
+ , fTOFmode(1)
, fEtaAbsCutLow(0.0)
, fEtaAbsCutUp(0.9)
+ , fPileUpRejectionType(AliAnalysisTaskPIDV0base::kPileUpRejectionOff)
, fDoAnySystematicStudiesOnTheExpectedSignal(kFALSE)
- , fSystematicScalingSplines(1.0)
+ , fSystematicScalingSplinesThreshold(50.)
+ , fSystematicScalingSplinesBelowThreshold(1.0)
+ , fSystematicScalingSplinesAboveThreshold(1.0)
, fSystematicScalingEtaCorrectionMomentumThr(0.35)
, fSystematicScalingEtaCorrectionLowMomenta(1.0)
, fSystematicScalingEtaCorrectionHighMomenta(1.0)
- , fSystematicScalingEtaSigmaPara(1.0)
+ , fSystematicScalingEtaSigmaParaThreshold(250.)
+ , fSystematicScalingEtaSigmaParaBelowThreshold(1.0)
+ , fSystematicScalingEtaSigmaParaAboveThreshold(1.0)
, fSystematicScalingMultCorrection(1.0)
, fCentralityEstimator("V0M")
, fhPIDdataAll(0x0)
, fGenRespPrDeltaPi(new Double_t[fgkMaxNumGenEntries])
, fGenRespPrDeltaPr(new Double_t[fgkMaxNumGenEntries])
*/
+ , fDeltaPrimeAxis(0x0)
+ , fhMaxEtaVariation(0x0)
, fhEventsProcessed(0x0)
- , fhSkippedTracksForSignalGeneration(0x0)
+ , fhEventsTriggerSel(0x0)
+ , fhEventsTriggerSelVtxCut(0x0)
+ , fhEventsProcessedNoPileUpRejection(0x0)
+ , fChargedGenPrimariesTriggerSel(0x0)
+ , fChargedGenPrimariesTriggerSelVtxCut(0x0)
+ , fChargedGenPrimariesTriggerSelVtxCutZ(0x0)
+ , fChargedGenPrimariesTriggerSelVtxCutZPileUpRej(0x0)
, fhMCgeneratedYieldsPrimaries(0x0)
, fh2FFJetPtRec(0x0)
, fh2FFJetPtGen(0x0)
, fh1Xsec(0x0)
, fh1Trials(0x0)
+ , fh1EvtsPtHardCut(0x0)
, fContainerEff(0x0)
+ , fQASharedCls(0x0)
+ , fDeDxCheck(0x0)
, fOutputContainer(0x0)
- , fPtResolutionContainer(0x0)
+ , fQAContainer(0x0)
{
// Constructor
fkDeltaPrimeLowLimit, fkDeltaPrimeUpLimit,
fkConvolutedGausNPar, "AliAnalysisTaskPID", "ConvolutedGaus");
+ // Set some arbitrary parameteres, such that the function call will not crash
+ // (although it should not be called with these parameters...)
+ fConvolutedGausDeltaPrime->SetParameter(0, 0);
+ fConvolutedGausDeltaPrime->SetParameter(1, 1);
+ fConvolutedGausDeltaPrime->SetParameter(2, 2);
+
+
// Initialisation of translation parameters is time consuming.
// Therefore, default values will only be initialised if they are really needed.
// Otherwise, it is left to the user to set the parameter properly.
delete fOutputContainer;
fOutputContainer = 0x0;
- delete fPtResolutionContainer;
- fPtResolutionContainer = 0x0;
+ delete fQAContainer;
+ fQAContainer = 0x0;
delete fConvolutedGausDeltaPrime;
fConvolutedGausDeltaPrime = 0x0;
+ delete fDeltaPrimeAxis;
+ fDeltaPrimeAxis = 0x0;
+
delete [] fGenRespElDeltaPrimeEl;
delete [] fGenRespElDeltaPrimeKa;
delete [] fGenRespElDeltaPrimePi;
fGenRespPrDeltaPrimePi = 0x0;
fGenRespPrDeltaPrimePr = 0x0;
+ delete fhMaxEtaVariation;
+ fhMaxEtaVariation = 0x0;
+
/*OLD with deltaSpecies
delete [] fGenRespElDeltaEl;
delete [] fGenRespElDeltaKa;
{
// Initialise the PIDcombined object
- if (!fDoPID)
+ if (!fDoPID && !fDoDeDxCheck)
return;
if(fDebug > 1)
}
+//________________________________________________________________________
+Bool_t AliAnalysisTaskPID::CalculateMaxEtaVariationMapFromPIDResponse()
+{
+ // Calculate the maximum deviation from unity of the eta correction factors for each row in 1/dEdx(splines)
+ // from the eta correction map of the TPCPIDResponse. The result is stored in fhMaxEtaVariation.
+
+ if (!fPIDResponse) {
+ AliError("No PID response!");
+ return kFALSE;
+ }
+
+ delete fhMaxEtaVariation;
+
+ const TH2D* hEta = fPIDResponse->GetTPCResponse().GetEtaCorrMap();
+ if (!hEta) {
+ AliError("No eta correction map!");
+ return kFALSE;
+ }
+
+ // Take binning from hEta in Y for fhMaxEtaVariation
+ fhMaxEtaVariation = hEta->ProjectionY("hMaxEtaVariation");
+ fhMaxEtaVariation->SetDirectory(0);
+ fhMaxEtaVariation->Reset();
+
+ // For each bin in 1/dEdx, loop of all tanTheta bins and find the maximum deviation from unity.
+ // Store the result in fhMaxEtaVariation
+
+ for (Int_t binY = 1; binY <= fhMaxEtaVariation->GetNbinsX(); binY++) {
+ Double_t maxAbs = -1;
+ for (Int_t binX = 1; binX <= hEta->GetNbinsX(); binX++) {
+ Double_t curr = TMath::Abs(hEta->GetBinContent(binX, binY) - 1.);
+ if (curr > maxAbs)
+ maxAbs = curr;
+ }
+
+ if (maxAbs < 1e-12) {
+ AliError(Form("Maximum deviation from unity is zero for 1/dEdx = %f (bin %d)", hEta->GetYaxis()->GetBinCenter(binY), binY));
+ delete fhMaxEtaVariation;
+ return kFALSE;
+ }
+
+ fhMaxEtaVariation->SetBinContent(binY, maxAbs);
+ }
+
+ printf("AliAnalysisTaskPID: Calculated max eta variation from map \"%s\".\n", hEta->GetTitle());
+
+ return kTRUE;
+}
+
+
//________________________________________________________________________
void AliAnalysisTaskPID::UserCreateOutputObjects()
{
if(fDebug > 1)
printf("File: %s, Line: %d: UserCreateOutputObjects\n", (char*)__FILE__, __LINE__);
- SetUpPIDcombined();
-
- // Input handler
- AliAnalysisManager* man = AliAnalysisManager::GetAnalysisManager();
- AliInputEventHandler* inputHandler = dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler());
+ // Setup basic things, like PIDResponse
+ AliAnalysisTaskPIDV0base::UserCreateOutputObjects();
- if (!inputHandler)
- AliFatal("Input handler needed");
- else {
- // PID response object
- fPIDResponse = inputHandler->GetPIDResponse();
- if (!fPIDResponse)
- AliFatal("PIDResponse object was not created");
- }
+ if (!fPIDResponse)
+ AliFatal("PIDResponse object was not created");
- if(fDebug > 2)
- printf("File: %s, Line: %d: UserCreateOutputObjects -> Retrieved PIDresponse object\n", (char*)__FILE__, __LINE__);
+ SetUpPIDcombined();
OpenFile(1);
11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0, 20.0, 22.0, 24.0,
26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 40.0, 45.0, 50.0 };
- const Int_t nCentBins = 12;
- //-1 for pp; 90-100 has huge electro-magnetic impurities
- Double_t binsCent[nCentBins+1] = {-1, 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 };
+ const Bool_t useITSTPCtrackletsCentEstimatorWithNewBinning = fCentralityEstimator.CompareTo("ITSTPCtracklets", TString::kIgnoreCase) == 0
+ && fStoreCentralityPercentile;
+
+ const Int_t nCentBinsGeneral = 12;
+ const Int_t nCentBinsNewITSTPCtracklets = 17;
+
+ const Int_t nCentBins = useITSTPCtrackletsCentEstimatorWithNewBinning ? nCentBinsNewITSTPCtracklets : nCentBinsGeneral;
+
+ Double_t binsCent[nCentBins+1];
+ for (Int_t i = 0; i < nCentBins + 1; i++)
+ binsCent[i] = -1;
+
+ //-1 for pp (unless explicitely requested); 90-100 has huge electro-magnetic impurities
+ Double_t binsCentV0[nCentBinsGeneral+1] = {-1, 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 };
+
+ // These centrality estimators deal with integers! This implies that the ranges are always [lowlim, uplim - 1]
+ Double_t binsCentITSTPCTracklets[nCentBinsNewITSTPCtracklets+1] = { -9999, 0, 1, 4, 7, 10, 15, 20, 25, 30, 40, 50, 60, 70, 80, 90, 100, 9999 };
+ Double_t binsCentITSTPCTrackletsOldPreliminary[nCentBinsGeneral+1] = { 0, 7, 13, 20, 29, 40, 50, 60, 72, 83, 95, 105, 115 };
+
+ // Special centrality binning for pp
+ Double_t binsCentpp[nCentBinsGeneral+1] = { 0, 0.01, 0.1, 1, 5, 10, 15, 20, 30, 40, 50, 70, 100};
+
+ if (fCentralityEstimator.CompareTo("ITSTPCtrackletsOldPreliminaryBinning", TString::kIgnoreCase) == 0 && fStoreCentralityPercentile) {
+ // Special binning for this centrality estimator; but keep number of bins!
+ for (Int_t i = 0; i < nCentBinsGeneral+1; i++)
+ binsCent[i] = binsCentITSTPCTrackletsOldPreliminary[i];
+ }
+ else if (fCentralityEstimator.CompareTo("ITSTPCtracklets", TString::kIgnoreCase) == 0 && fStoreCentralityPercentile) {
+ // Special binning for this centrality estimator and different number of bins!
+ for (Int_t i = 0; i < nCentBinsNewITSTPCtracklets+1; i++)
+ binsCent[i] = binsCentITSTPCTracklets[i];
+ }
+ else if (fCentralityEstimator.Contains("ppMult", TString::kIgnoreCase) && fStoreCentralityPercentile) {
+ // Special binning for this pp centrality estimator; but keep number of bins!
+ for (Int_t i = 0; i < nCentBinsGeneral+1; i++)
+ binsCent[i] = binsCentpp[i];
+ }
+ else {
+ // Take default binning for VZERO
+ for (Int_t i = 0; i < nCentBinsGeneral+1; i++)
+ binsCent[i] = binsCentV0[i];
+ }
const Int_t nJetPtBins = 11;
Double_t binsJetPt[nJetPtBins+1] = {0, 2, 5, 10, 15, 20, 30, 40, 60, 80, 120, 200};
deltaPrimeBins[i] = factor * deltaPrimeBins[i - 1];
}
+ fDeltaPrimeAxis = new TAxis(deltaPrimeNBins, deltaPrimeBins);
+
const Int_t nMCPIDbins = 5;
const Double_t mcPIDmin = 0.;
const Double_t mcPIDmax = 5.;
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);
fhGenPr = new THnSparseD("hGenPr", "", nGenBins, genBins, genXmin, genXmax);
SetUpGenHist(fhGenPr, binsPt, deltaPrimeBins, binsCent, binsJetPt);
fOutputContainer->Add(fhGenPr);
-
-
- fhEventsProcessed = new TH1D("fhEventsProcessed","Number of processed events;Centrality percentile", nCentBins,
- binsCent);
- fhEventsProcessed->Sumw2();
- fOutputContainer->Add(fhEventsProcessed);
-
- fhSkippedTracksForSignalGeneration = new TH2D("fhSkippedTracksForSignalGeneration",
- "Number of tracks skipped for the signal generation;P_{T}^{gen} (GeV/c);TPC signal N",
- nPtBins, binsPt, 161, -0.5, 160.5);
- fhSkippedTracksForSignalGeneration->Sumw2();
- fOutputContainer->Add(fhSkippedTracksForSignalGeneration);
}
+ fhEventsProcessed = new TH1D("fhEventsProcessed",
+ "Number of events passing trigger selection, vtx and zvtx cuts and pile-up rejection;Centrality Percentile",
+ nCentBins, binsCent);
+ fhEventsProcessed->Sumw2();
+ fOutputContainer->Add(fhEventsProcessed);
+
+ fhEventsTriggerSelVtxCut = new TH1D("fhEventsTriggerSelVtxCut",
+ "Number of events passing trigger selection and vtx cut;Centrality Percentile",
+ nCentBins, binsCent);
+ fhEventsTriggerSelVtxCut->Sumw2();
+ fOutputContainer->Add(fhEventsTriggerSelVtxCut);
+
+ fhEventsTriggerSel = new TH1D("fhEventsTriggerSel",
+ "Number of events passing trigger selection;Centrality Percentile",
+ nCentBins, binsCent);
+ fOutputContainer->Add(fhEventsTriggerSel);
+ fhEventsTriggerSel->Sumw2();
+
+
+ fhEventsProcessedNoPileUpRejection = new TH1D("fhEventsProcessedNoPileUpRejection",
+ "Number of events passing trigger selection, vtx and zvtx cuts;Centrality Percentile",
+ nCentBins, binsCent);
+ fOutputContainer->Add(fhEventsProcessedNoPileUpRejection);
+ fhEventsProcessedNoPileUpRejection->Sumw2();
+
+
// Generated yields within acceptance
const Int_t nBinsGenYields = fStoreAdditionalJetInformation ? kGenYieldNumAxes : kGenYieldNumAxes - 3;
Int_t genYieldsBins[kGenYieldNumAxes] = { nMCPIDbins, nPtBins, nCentBins, nJetPtBins, nZBins, nXiBins,
const Int_t nEffDims = fStoreAdditionalJetInformation ? kEffNumAxes : kEffNumAxes - 3; // Number of dimensions for the efficiency
const Int_t nMCIDbins = AliPID::kSPECIES;
- Double_t binsMCID[nMCIDbins];
+ Double_t binsMCID[nMCIDbins + 1];
- for(Int_t i = 0; i < nMCIDbins; i++) {
+ for(Int_t i = 0; i <= nMCIDbins; i++) {
binsMCID[i]= i;
}
}
fContainerEff->SetVarTitle(kEffMCID,"MC ID");
- fContainerEff->SetVarTitle(kEffTrackPt,"P_{T} (GeV/c)");
+ fContainerEff->SetVarTitle(kEffTrackPt,"p_{T} (GeV/c)");
fContainerEff->SetVarTitle(kEffTrackEta,"#eta");
fContainerEff->SetVarTitle(kEffTrackCharge,"Charge (e_{0})");
fContainerEff->SetVarTitle(kEffCentrality, "Centrality Percentile");
if (fStoreAdditionalJetInformation) {
- fContainerEff->SetVarTitle(kEffJetPt, "P_{T}^{jet} (GeV/c)");
- fContainerEff->SetVarTitle(kEffZ, "z = P_{T}^{track} / P_{T}^{jet}");
- fContainerEff->SetVarTitle(kEffXi, "#xi = ln(P_{T}^{jet} / P_{T}^{track})");
+ fContainerEff->SetVarTitle(kEffJetPt, "p_{T}^{jet} (GeV/c)");
+ fContainerEff->SetVarTitle(kEffZ, "z = p_{T}^{track} / p_{T}^{jet}");
+ fContainerEff->SetVarTitle(kEffXi, "#xi = ln(p_{T}^{jet} / p_{T}^{track})");
}
// Define clean MC sample
if (fDoPID || fDoEfficiency) {
// Generated jets
- fh2FFJetPtRec = new TH2D("fh2FFJetPtRec", "Number of reconstructed jets;Centrality Percentile;P_{T}^{jet} (GeV/c)",
+ fh2FFJetPtRec = new TH2D("fh2FFJetPtRec", "Number of reconstructed jets;Centrality Percentile;p_{T}^{jet} (GeV/c)",
nCentBins, binsCent, nJetPtBins, binsJetPt);
fh2FFJetPtRec->Sumw2();
fOutputContainer->Add(fh2FFJetPtRec);
- fh2FFJetPtGen = new TH2D("fh2FFJetPtGen", "Number of generated jets;Centrality Percentile;P_{T}^{jet} (GeV/c)",
+ fh2FFJetPtGen = new TH2D("fh2FFJetPtGen", "Number of generated jets;Centrality Percentile;p_{T}^{jet} (GeV/c)",
nCentBins, binsCent, nJetPtBins, binsJetPt);
fh2FFJetPtGen->Sumw2();
fOutputContainer->Add(fh2FFJetPtGen);
fh1Trials->Sumw2();
fh1Trials->GetXaxis()->SetBinLabel(1, "#sum{ntrials}");
+ fh1EvtsPtHardCut = new TH1F("fh1EvtsPtHardCut", "#events before and after MC #it{p}_{T,hard} cut;;Events",2,0,2);
+ fh1EvtsPtHardCut->Sumw2();
+ fh1EvtsPtHardCut->GetXaxis()->SetBinLabel(1, "All");
+ fh1EvtsPtHardCut->GetXaxis()->SetBinLabel(2, "#it{p}_{T,hard}");
+
fOutputContainer->Add(fh1Xsec);
fOutputContainer->Add(fh1Trials);
+ fOutputContainer->Add(fh1EvtsPtHardCut);
- if (fDoPtResolution) {
+ if (fDoDeDxCheck || fDoPtResolution) {
OpenFile(3);
+ fQAContainer = new TObjArray(1);
+ fQAContainer->SetName(Form("%s_QA", GetName()));
+ fQAContainer->SetOwner(kTRUE);
if(fDebug > 2)
- printf("File: %s, Line: %d: UserCreateOutputObjects -> OpenFile(3) successful\n", (char*)__FILE__, __LINE__);
-
- fPtResolutionContainer = new TObjArray(1);
- fPtResolutionContainer->SetName(Form("%s_PtResolution", GetName()));
- fPtResolutionContainer->SetOwner(kTRUE);
+ printf("File: %s, Line: %d: UserCreateOutputObjects -> OpenFile(3) successful\n", (char*)__FILE__, __LINE__);
+ }
+
+ if (fDoPtResolution) {
+ const Int_t nPtBinsRes = 100;
+ Double_t pTbinsRes[nPtBinsRes + 1];
+
+ const Double_t fromLowPtRes = 0.15;
+ const Double_t toHighPtRes = 50.;
+ const Double_t factorPtRes = TMath::Power(toHighPtRes/fromLowPtRes, 1./nPtBinsRes);
+ // Log binning for whole pT range
+ pTbinsRes[0] = fromLowPtRes;
+ for (Int_t i = 0 + 1; i <= nPtBinsRes; i++) {
+ pTbinsRes[i] = factorPtRes * pTbinsRes[i - 1];
+ }
const Int_t nBinsPtResolution = kPtResNumAxes;
- Int_t ptResolutionBins[kPtResNumAxes] = { nJetPtBins, nPtBins, nPtBins };
- Double_t ptResolutionXmin[kPtResNumAxes] = { binsJetPt[0], binsPt[0], binsPt[0] };
- Double_t ptResolutionXmax[kPtResNumAxes] = { binsJetPt[nJetPtBins], binsPt[nPtBins], binsPt[nPtBins] };
+ Int_t ptResolutionBins[kPtResNumAxes] = { nJetPtBins, nPtBinsRes, nPtBinsRes,
+ nChargeBins, nCentBins };
+ Double_t ptResolutionXmin[kPtResNumAxes] = { binsJetPt[0], pTbinsRes[0], pTbinsRes[0],
+ binsCharge[0], binsCent[0] };
+ Double_t ptResolutionXmax[kPtResNumAxes] = { binsJetPt[nJetPtBins], pTbinsRes[nPtBinsRes], pTbinsRes[nPtBinsRes],
+ binsCharge[nChargeBins], binsCent[nCentBins] };
for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
fPtResolution[i] = new THnSparseD(Form("fPtResolution_%s", AliPID::ParticleShortName(i)),
Form("Pt resolution for primaries, %s", AliPID::ParticleLatexName(i)),
nBinsPtResolution, ptResolutionBins, ptResolutionXmin, ptResolutionXmax);
- SetUpPtResHist(fPtResolution[i], binsPt, binsJetPt);
- fPtResolutionContainer->Add(fPtResolution[i]);
+ SetUpPtResHist(fPtResolution[i], pTbinsRes, binsJetPt, binsCent);
+ fQAContainer->Add(fPtResolution[i]);
}
+
+
+ // Besides the pT resolution, also perform check on shared clusters
+ const Int_t nBinsQASharedCls = kQASharedClsNumAxes;
+ Int_t qaSharedClsBins[kQASharedClsNumAxes] = { nJetPtBins, nPtBinsRes, 160, 160 };
+ Double_t qaSharedClsXmin[kQASharedClsNumAxes] = { binsJetPt[0], pTbinsRes[0], 0, -1 };
+ Double_t qaSharedClsXmax[kQASharedClsNumAxes] = { binsJetPt[nJetPtBins], pTbinsRes[nPtBinsRes], 160, 159 };
+
+ fQASharedCls = new THnSparseD("fQASharedCls", "QA shared clusters", nBinsQASharedCls, qaSharedClsBins, qaSharedClsXmin, qaSharedClsXmax);
+
+ SetUpSharedClsHist(fQASharedCls, pTbinsRes, binsJetPt);
+ fQAContainer->Add(fQASharedCls);
+ }
+
+
+
+ if (fDoDeDxCheck) {
+ const Int_t nEtaAbsBins = 9;
+ const Double_t binsEtaAbs[nEtaAbsBins+1] = { 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 };
+
+ const Double_t dEdxMin = 20;
+ const Double_t dEdxMax = 110;
+ const Int_t nDeDxBins = (Int_t) ((dEdxMax - dEdxMin) / 0.02);
+ const Int_t nBinsDeDxCheck= kDeDxCheckNumAxes;
+ Int_t dEdxCheckBins[kDeDxCheckNumAxes] = { nSelSpeciesBins, nPtBins, nJetPtBins, nEtaAbsBins, nDeDxBins };
+ Double_t dEdxCheckXmin[kDeDxCheckNumAxes] = { selSpeciesMin, binsPt[0], binsJetPt[0], binsEtaAbs[0], dEdxMin };
+ Double_t dEdxCheckXmax[kDeDxCheckNumAxes] = { selSpeciesMax, binsPt[nPtBins], binsJetPt[nJetPtBins], binsEtaAbs[nEtaAbsBins], dEdxMax };
+
+ fDeDxCheck = new THnSparseD("fDeDxCheck", "dEdx check", nBinsDeDxCheck, dEdxCheckBins, dEdxCheckXmin, dEdxCheckXmax);
+ SetUpDeDxCheckHist(fDeDxCheck, binsPt, binsJetPt, binsEtaAbs);
+ fQAContainer->Add(fDeDxCheck);
+ }
+
+ if (fDoBinZeroStudy) {
+ const Double_t etaLow = -0.9;
+ const Double_t etaUp = 0.9;
+ const Int_t nEtaBins = 18;
+
+ const Int_t nBinsBinZeroStudy = kBinZeroStudyNumAxes;
+ Int_t binZeroStudyBins[nBinsBinZeroStudy] = { nCentBins, nPtBins, nEtaBins };
+ Double_t binZeroStudyXmin[nBinsBinZeroStudy] = { binsCent[0], binsPt[0], etaLow };
+ Double_t binZeroStudyXmax[nBinsBinZeroStudy] = { binsCent[nCentBins], binsPt[nPtBins], etaUp };
+
+ fChargedGenPrimariesTriggerSel = new THnSparseD("fChargedGenPrimariesTriggerSel", "Trigger sel.", nBinsBinZeroStudy, binZeroStudyBins,
+ binZeroStudyXmin, binZeroStudyXmax);
+ SetUpBinZeroStudyHist(fChargedGenPrimariesTriggerSel, binsCent, binsPt);
+ fOutputContainer->Add(fChargedGenPrimariesTriggerSel);
+
+ fChargedGenPrimariesTriggerSelVtxCut = new THnSparseD("fChargedGenPrimariesTriggerSelVtxCut", "Vertex cut", nBinsBinZeroStudy,
+ binZeroStudyBins, binZeroStudyXmin, binZeroStudyXmax);
+ SetUpBinZeroStudyHist(fChargedGenPrimariesTriggerSelVtxCut, binsCent, binsPt);
+ fOutputContainer->Add(fChargedGenPrimariesTriggerSelVtxCut);
+
+ fChargedGenPrimariesTriggerSelVtxCutZ = new THnSparseD("fChargedGenPrimariesTriggerSelVtxCutZ", "Vertex #it{z} cut", nBinsBinZeroStudy,
+ binZeroStudyBins, binZeroStudyXmin, binZeroStudyXmax);
+ SetUpBinZeroStudyHist(fChargedGenPrimariesTriggerSelVtxCutZ, binsCent, binsPt);
+ fOutputContainer->Add(fChargedGenPrimariesTriggerSelVtxCutZ);
+
+ fChargedGenPrimariesTriggerSelVtxCutZPileUpRej = new THnSparseD("fChargedGenPrimariesTriggerSelVtxCutZPileUpRej", "Vertex #it{z} cut",
+ nBinsBinZeroStudy, binZeroStudyBins, binZeroStudyXmin, binZeroStudyXmax);
+ SetUpBinZeroStudyHist(fChargedGenPrimariesTriggerSelVtxCutZPileUpRej, binsCent, binsPt);
+ fOutputContainer->Add(fChargedGenPrimariesTriggerSelVtxCutZPileUpRej);
}
if(fDebug > 2)
printf("File: %s, Line: %d: UserCreateOutputObjects -> Posting output data\n", (char*)__FILE__, __LINE__);
- PostOutputData();
+ PostData(1, fOutputContainer);
+ if (fDoEfficiency)
+ PostData(2, fContainerEff);
+ if (fDoDeDxCheck || fDoPtResolution)
+ PostData(3, fQAContainer);
if(fDebug > 2)
printf("File: %s, Line: %d: UserCreateOutputObjects -> Done\n", (char*)__FILE__, __LINE__);
{
// Main loop
// Called for each event
-
+
if(fDebug > 1)
printf("File: %s, Line: %d: UserExec\n", (char*)__FILE__, __LINE__);
return;
}
+ ConfigureTaskForCurrentEvent(fEvent);
+
fMC = dynamic_cast<AliMCEvent*>(MCEvent());
if (!fPIDResponse || !fPIDcombined)
return;
- if (!GetVertexIsOk(fEvent))
- return;
+ Double_t centralityPercentile = -1;
+ if (fStoreCentralityPercentile) {
+ if (fCentralityEstimator.Contains("ITSTPCtracklets", TString::kIgnoreCase)) {
+ // Special pp centrality estimator
+ AliESDEvent* esdEvent = dynamic_cast<AliESDEvent*>(fEvent);
+ if (!esdEvent) {
+ AliError("Not esd event -> Cannot use tracklet multiplicity estimator!");
+ centralityPercentile = -1;
+ }
+ else
+ centralityPercentile = AliESDtrackCuts::GetReferenceMultiplicity(esdEvent, AliESDtrackCuts::kTrackletsITSTPC, fEtaAbsCutUp);
+ }
+ else if (fCentralityEstimator.Contains("ppMult", TString::kIgnoreCase)) {
+ // Another special pp centrality estimator
+ centralityPercentile = fAnaUtils->GetMultiplicityPercentile(fEvent, GetPPCentralityEstimator().Data());
+ }
+ else {
+ // Ordinary centrality estimator
+ centralityPercentile = fEvent->GetCentrality()->GetCentralityPercentile(fCentralityEstimator.Data());
+ }
+ }
+
+ // Check if vertex is ok, but don't apply cut on z position
+ const Bool_t passedVertexSelection = GetVertexIsOk(fEvent, kFALSE);
+ // Now check again, but also require z position to be in desired range
+ const Bool_t passedVertexZSelection = GetVertexIsOk(fEvent, kTRUE);
+ // Check pile-up
+ const Bool_t isPileUp = GetIsPileUp(fEvent, fPileUpRejectionType);
+
- fESD = dynamic_cast<AliESDEvent*>(fEvent);
- const AliVVertex* primaryVertex = fESD ? fESD->GetPrimaryVertexTracks() : fEvent->GetPrimaryVertex();
- if (!primaryVertex)
+
+ if (fDoBinZeroStudy && fMC) {
+ for (Int_t iPart = 0; iPart < fMC->GetNumberOfTracks(); iPart++) {
+ AliMCParticle *mcPart = dynamic_cast<AliMCParticle*>(fMC->GetTrack(iPart));
+
+ if (!mcPart)
+ continue;
+
+ if (!fMC->IsPhysicalPrimary(iPart))
+ continue;
+
+ const Double_t etaGen = mcPart->Eta();
+ const Double_t ptGen = mcPart->Pt();
+
+ Double_t values[kBinZeroStudyNumAxes] = { 0. };
+ values[kBinZeroStudyCentrality] = centralityPercentile;
+ values[kBinZeroStudyGenPt] = ptGen;
+ values[kBinZeroStudyGenEta] = etaGen;
+
+ fChargedGenPrimariesTriggerSel->Fill(values);
+ if (passedVertexSelection) {
+ fChargedGenPrimariesTriggerSelVtxCut->Fill(values);
+ if (passedVertexZSelection) {
+ fChargedGenPrimariesTriggerSelVtxCutZ->Fill(values);
+ if (!isPileUp) {
+ fChargedGenPrimariesTriggerSelVtxCutZPileUpRej->Fill(values);
+ }
+ }
+ }
+ }
+ }
+
+
+
+ // Event counters for trigger selection, vertex cuts and pile-up rejection
+ IncrementEventCounter(centralityPercentile, kTriggerSel);
+
+ if (!passedVertexSelection)
return;
- if(primaryVertex->GetNContributors() <= 0)
+ IncrementEventCounter(centralityPercentile, kTriggerSelAndVtxCut);
+
+ if (!passedVertexZSelection)
return;
- Double_t magField = fEvent->GetMagneticField();
+ IncrementEventCounter(centralityPercentile, kTriggerSelAndVtxCutAndZvtxCutNoPileUpRejection);
+ // ATTENTION: Is this the right place for the pile-up rejection? Important to have still the proper bin-0 correction,
+ // which is done solely with sel and selVtx, since the zvtx selection does ~not change the spectra. The question is whether the pile-up
+ // rejection changes the spectra. If not, then it is perfectly fine to put it here and keep the usual histo for the normalisation to number
+ // of events. But if it does change the spectra, this must somehow be corrected for.
+ // NOTE: multiplicity >= 0 usually implies a properly reconstructed vertex. Hence, the bin-0 correction cannot be done in multiplicity bins.
+ // Furthermore, there seems to be no MC simulation with pile-up rejection, so the bin-0 correction cannot be extracted with it. Pile-up
+ // rejection has only a minor impact, so maybe there is no need to dig further.
+ if (isPileUp)
+ return;
- //OLD with DeltaSpecies const Bool_t usePureGausForDelta = kTRUE;
+ IncrementEventCounter(centralityPercentile, kTriggerSelAndVtxCutAndZvtxCut);
-
- Double_t centralityPercentile = -1;
- if (fStoreCentralityPercentile)
- centralityPercentile = fEvent->GetCentrality()->GetCentralityPercentile(fCentralityEstimator.Data());
+ Double_t magField = fEvent->GetMagneticField();
if (fMC) {
if (fDoPID || fDoEfficiency) {
// - 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());
continue;
if (fDoPID) {
- Double_t valuesGenYield[kGenYieldNumAxes] = { mcID, mcPart->Pt(), centralityPercentile, -1, -1, -1, -1 };
+ Double_t valuesGenYield[kGenYieldNumAxes] = { static_cast<Double_t>(mcID), mcPart->Pt(), centralityPercentile, -1, -1, -1, -1 };
valuesGenYield[GetIndexOfChargeAxisGenYield()] = chargeMC;
fhMCgeneratedYieldsPrimaries->Fill(valuesGenYield);
if (fDoEfficiency) {
- Double_t valueEff[kEffNumAxes] = { mcID, mcPart->Pt(), mcPart->Eta(), chargeMC, centralityPercentile,
+ Double_t valueEff[kEffNumAxes] = { static_cast<Double_t>(mcID), mcPart->Pt(), mcPart->Eta(), chargeMC, centralityPercentile,
-1, -1, -1 };
fContainerEff->Fill(valueEff, kStepGenWithGenCuts);
// Apply detector level track cuts
- //if (track->GetTPCsignalN() < 60)
- // continue;
-
- Double_t dEdxTPC = fPIDResponse->IsTunedOnData() ? fPIDResponse->GetTPCsignalTunedOnData(track) : track->GetTPCsignal();
+ const Bool_t tuneOnDataTPC = fPIDResponse->IsTunedOnData() &&
+ ((fPIDResponse->GetTunedOnDataMask() & AliPIDResponse::kDetTPC) == AliPIDResponse::kDetTPC);
+ Double_t dEdxTPC = tuneOnDataTPC ? fPIDResponse->GetTPCsignalTunedOnData(track) : track->GetTPCsignal();
if (dEdxTPC <= 0)
continue;
if(fTrackFilter && !fTrackFilter->IsSelected(track))
continue;
- if (fUseTPCCutMIGeo) {
+ if (GetUseTPCCutMIGeo()) {
if (!TPCCutMIGeo(track, fEvent))
continue;
}
+ else if (GetUseTPCnclCut()) {
+ if (!TPCnclCut(track))
+ continue;
+ }
if(fUsePhiCut) {
if (!PhiPrimeCut(track, magField))
// 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,
+ Double_t value[kEffNumAxes] = { static_cast<Double_t>(mcID), mcTrack->Pt(), mcTrack->Eta(), mcTrack->Charge() / 3., centralityPercentile,
-1, -1, -1 };
fContainerEff->Fill(value, kStepRecWithGenCuts);
- Double_t valueMeas[kEffNumAxes] = { mcID, track->Pt(), track->Eta(), track->Charge(), centralityPercentile,
+ Double_t valueMeas[kEffNumAxes] = { static_cast<Double_t>(mcID), track->Pt(), track->Eta(), static_cast<Double_t>(track->Charge()), centralityPercentile,
-1, -1, -1 };
fContainerEff->Fill(valueMeas, kStepRecWithGenCutsMeasuredObs);
}
}
// 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)
+ if (fDoPID || fDoDeDxCheck || fDoPtResolution)
ProcessTrack(track, pdg, centralityPercentile, -1); // No jet information in this case -> Set jet pT to -1
if (fDoPtResolution) {
if (mcTrack && fMC->IsPhysicalPrimary(TMath::Abs(label))) {
- Double_t valuePtRes[kPtResNumAxes] = { -1, mcTrack->Pt(), track->Pt() };
- fPtResolution[mcID]->Fill(valuePtRes);
+ // AliMCParticle->Charge() calls TParticlePDG->Charge(), which returns the charge in units of e0 / 3
+ Double_t valuePtRes[kPtResNumAxes] = { -1, mcTrack->Pt(), track->Pt(), mcTrack->Charge() / 3., centralityPercentile };
+ FillPtResolution(mcID, valuePtRes);
}
}
if (fDoEfficiency) {
if (mcTrack) {
- Double_t valueRecAllCuts[kEffNumAxes] = { mcID, track->Pt(), track->Eta(), track->Charge(), centralityPercentile,
+ Double_t valueRecAllCuts[kEffNumAxes] = { static_cast<Double_t>(mcID), track->Pt(), track->Eta(), static_cast<Double_t>(track->Charge()), centralityPercentile,
-1, -1, -1 };
fContainerEff->Fill(valueRecAllCuts, kStepRecWithRecCutsMeasuredObs);
fContainerEff->Fill(valueRecAllCuts, kStepRecWithRecCutsMeasuredObsStrangenessScaled, weight);
// AliMCParticle->Charge() calls TParticlePDG->Charge(), which returns the charge in units of e0 / 3
- Double_t valueGenAllCuts[kEffNumAxes] = { mcID, mcTrack->Pt(), mcTrack->Eta(), mcTrack->Charge() / 3.,
+ Double_t valueGenAllCuts[kEffNumAxes] = { static_cast<Double_t>(mcID), mcTrack->Pt(), mcTrack->Eta(), mcTrack->Charge() / 3.,
centralityPercentile, -1, -1, -1 };
if (fMC->IsPhysicalPrimary(TMath::Abs(label))) {
fContainerEff->Fill(valueRecAllCuts, kStepRecWithRecCutsMeasuredObsPrimaries);
}
} //track loop
- IncrementEventsProcessed(centralityPercentile);
-
if(fDebug > 2)
printf("File: %s, Line: %d: UserExec -> Processing done\n", (char*)__FILE__, __LINE__);
fDoAnySystematicStudiesOnTheExpectedSignal = kFALSE;
- if (TMath::Abs(fSystematicScalingSplines - 1.0) > fgkEpsilon) {
+
+ if ((TMath::Abs(fSystematicScalingSplinesBelowThreshold - 1.0) > fgkEpsilon) ||
+ (TMath::Abs(fSystematicScalingSplinesAboveThreshold - 1.0) > fgkEpsilon)) {
fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
return;
}
return;
}
- if (TMath::Abs(fSystematicScalingEtaSigmaPara - 1.0) > fgkEpsilon) {
+ if ((TMath::Abs(fSystematicScalingEtaSigmaParaBelowThreshold - 1.0) > fgkEpsilon) ||
+ (TMath::Abs(fSystematicScalingEtaSigmaParaAboveThreshold - 1.0) > fgkEpsilon)) {
fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
return;
}
}
+//_____________________________________________________________________________
+void AliAnalysisTaskPID::ConfigureTaskForCurrentEvent(AliVEvent* event)
+{
+ // Configure the task for the current event. In particular, this is needed if the run number changes
+
+ if (!event) {
+ AliError("Could not set up task: no event!");
+ return;
+ }
+
+ Int_t run = event->GetRunNumber();
+
+ if (run != fRun){
+ // If systematics on eta is investigated, need to calculate the maxEtaVariationMap
+ if ((TMath::Abs(fSystematicScalingEtaCorrectionLowMomenta - 1.0) > fgkEpsilon) ||
+ (TMath::Abs(fSystematicScalingEtaCorrectionHighMomenta - 1.0) > fgkEpsilon)) {
+ if (!CalculateMaxEtaVariationMapFromPIDResponse())
+ AliFatal("Systematics on eta correction requested, but failed to calculate max eta varation map!");
+ }
+ }
+
+ fRun = run;
+}
+
+
//_____________________________________________________________________________
Int_t AliAnalysisTaskPID::PDGtoMCID(Int_t pdg)
{
//_____________________________________________________________________________
-void AliAnalysisTaskPID::GetJetTrackObservables(const Double_t trackPt, const Double_t jetPt, Double_t& z, Double_t& xi)
+void AliAnalysisTaskPID::GetJetTrackObservables(Double_t trackPt, Double_t jetPt, Double_t& z, Double_t& xi)
{
// Uses trackPt and jetPt to obtain z and xi.
//_____________________________________________________________________________
-inline Double_t AliAnalysisTaskPID::FastGaus(const Double_t x, const Double_t mean, const Double_t sigma) const
+inline Double_t AliAnalysisTaskPID::FastGaus(Double_t x, Double_t mean, Double_t sigma) const
{
// Calculate an unnormalised gaussian function with mean and sigma.
//_____________________________________________________________________________
-inline Double_t AliAnalysisTaskPID::FastNormalisedGaus(const Double_t x, const Double_t mean, const Double_t sigma) const
+inline Double_t AliAnalysisTaskPID::FastNormalisedGaus(Double_t x, Double_t mean, Double_t sigma) const
{
// Calculate a normalised (divided by sqrt(2*Pi)*sigma) gaussian function with mean and sigma.
//_____________________________________________________________________________
-Int_t AliAnalysisTaskPID::FindFirstBinAboveIn3dSubset(const TH3* hist, const Double_t threshold, const Int_t yBin,
- const Int_t zBin) const
+Int_t AliAnalysisTaskPID::FindFirstBinAboveIn3dSubset(const TH3* hist, Double_t threshold, Int_t yBin,
+ Int_t zBin) const
{
// Kind of projects a TH3 to 1 bin combination in y and z
// and looks for the first x bin above a threshold for this projection.
//_____________________________________________________________________________
-Int_t AliAnalysisTaskPID::FindLastBinAboveIn3dSubset(const TH3* hist, const Double_t threshold, const Int_t yBin,
- const Int_t zBin) const
+Int_t AliAnalysisTaskPID::FindLastBinAboveIn3dSubset(const TH3* hist, Double_t threshold, Int_t yBin,
+ Int_t zBin) const
{
// Kind of projects a TH3 to 1 bin combination in y and z
// and looks for the last x bin above a threshold for this projection.
//_____________________________________________________________________________
-Bool_t AliAnalysisTaskPID::GetParticleFraction(const Double_t trackPt, const Double_t jetPt, const Double_t centralityPercentile,
- const AliPID::EParticleType species,
+Bool_t AliAnalysisTaskPID::GetParticleFraction(Double_t trackPt, Double_t jetPt, Double_t centralityPercentile,
+ AliPID::EParticleType species,
Double_t& fraction, Double_t& fractionErrorStat, Double_t& fractionErrorSys) const
{
// Computes the particle fraction for the corresponding species for the given trackPt, jetPt and centrality.
else {
Double_t x0 = 0., x1 = 0., y0 = 0., y1 = 0.;
Double_t y0errStat = 0., y1errStat = 0., y0errSys = 0., y1errSys = 0.;
- Int_t trackPtBin = xAxis->FindBin(trackPt);
+ Int_t trackPtBin = xAxis->FindFixBin(trackPt);
// Linear interpolation between nearest neighbours in trackPt
if (trackPt <= xAxis->GetBinCenter(trackPtBin)) {
//_____________________________________________________________________________
-Bool_t AliAnalysisTaskPID::GetParticleFractions(const Double_t trackPt, const Double_t jetPt, const Double_t centralityPercentile,
- Double_t* prob, const Int_t smearSpeciesByError,
- const Int_t takeIntoAccountSpeciesSysError, const Bool_t uniformSystematicError) const
+Bool_t AliAnalysisTaskPID::GetParticleFractions(Double_t trackPt, Double_t jetPt, Double_t centralityPercentile,
+ Double_t* prob, Int_t smearSpeciesByError,
+ Int_t takeIntoAccountSpeciesSysError, Bool_t uniformSystematicError) const
{
// Fills the particle fractions for the given trackPt, jetPt and centrality into "prob".
// Use jetPt = -1 for inclusive spectra and centralityPercentile = -1 for pp.
//_____________________________________________________________________________
-const TH3D* AliAnalysisTaskPID::GetParticleFractionHisto(const Int_t species, const Bool_t sysError) const
+const TH3D* AliAnalysisTaskPID::GetParticleFractionHisto(Int_t species, Bool_t sysError) const
{
if (species < AliPID::kElectron || species > AliPID::kProton)
return 0x0;
}
if (absMotherPDG == 3122) { // Lambda
+ //if (absMotherPDG == 3122 || absMotherPDG == 3112 || absMotherPDG == 3222) { // Lambda / Sigma- / Sigma+
if (0.00 <= motherGenPt && motherGenPt < 0.20) fac = 0.645162;
else if(0.20 <= motherGenPt && motherGenPt < 0.40) fac = 0.627431;
else if(0.40 <= motherGenPt && motherGenPt < 0.60) fac = 0.457136;
}
+// _________________________________________________________________________________
+AliAnalysisTaskPID::TOFpidInfo AliAnalysisTaskPID::GetTOFType(const AliVTrack* track, Int_t tofMode) const
+{
+ // Get the (locally defined) particle type judged by TOF
+
+ if (!fPIDResponse) {
+ Printf("ERROR: fPIDResponse not available -> Cannot determine TOF type!");
+ return kNoTOFinfo;
+ }
+
+ /*TODO still needs some further thinking how to implement it....
+ // TOF PID status (kTOFout, kTIME) automatically checked by return value of ComputeTPCProbability;
+ // also, probability array will be set there (no need to initialise here)
+ Double_t p[AliPID::kSPECIES];
+ const AliPIDResponse::EDetPidStatus tofStatus = fPIDResponse->ComputeTPCProbability(track, AliPID::kSPECIES, p);
+ if (tofStatus != AliPIDResponse::kDetPidOk)
+ return kNoTOFinfo;
+
+ // Do not consider muons
+ p[AliPID::kMuon] = 0.;
+
+ // Probabilities are not normalised yet
+ Double_t sum = 0.;
+ for (Int_t i = 0; i < AliPID::kSPECIES; i++)
+ sum += p[i];
+
+ if (sum <= 0.)
+ return kNoTOFinfo;
+
+ for (Int_t i = 0; i < AliPID::kSPECIES; i++)
+ p[i] /= sum;
+
+ Double_t probThreshold = -999.;
+
+ // If there is only one distribution, the threshold corresponds to...
+ if (tofMode == 0) {
+ probThreshold = ;
+ }
+ else if (tofMode == 1) { // default
+ probThreshold = 0.9973; // a 3-sigma inclusion cut
+ }
+ else if (tofMode == 2) {
+ inclusion = 3.;
+ exclusion = 3.5;
+ }
+ else {
+ Printf("ERROR: Bad TOF mode: %d!", tofMode);
+ return kNoTOFinfo;
+ }
+
+ */
+
+ ///* OLD: cut with n-sigmas (ATTENTION: Does not take into account properly the TOF tail!)
+ // Check kTOFout, kTIME, mismatch
+ const AliPIDResponse::EDetPidStatus tofStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF, track);
+ if (tofStatus != AliPIDResponse::kDetPidOk)
+ return kNoTOFinfo;
+
+ Double_t nsigma[kNumTOFspecies + 1] = { -999., -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 = 3.;
+ exclusion = 2.5;
+ }
+ else if (tofMode == 1) { // default
+ inclusion = 3.;
+ exclusion = 3.;
+ }
+ else if (tofMode == 2) {
+ inclusion = 3.;
+ exclusion = 3.5;
+ }
+ else {
+ Printf("ERROR: Bad TOF mode: %d!", tofMode);
+ return kNoTOFinfo;
+ }
+
+ // Do not cut on nSigma electron because this would also remove pions in a large pT range.
+ // The electron contamination of the pion sample can be treated by TPC dEdx fits afterwards.
+ 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;
+ //*/
+
+ // There are no TOF electrons selected because the purity is rather bad, even for small momenta
+ // (also a small mismatch probability significantly affects electrons because their fraction
+ // is small). There is no need for TOF electrons anyway, since the dEdx distribution of kaons and
+ // protons in a given pT bin (also at the dEdx crossings) is very different from that of electrons.
+ // This is due to the steeply falling dEdx of p and K with momentum, whereas the electron dEdx stays
+ // rather constant.
+ // As a result, the TPC fit yields a more accurate electron fraction than the TOF selection can do.
+
+ return kNoTOFpid;
+}
+
+
// _________________________________________________________________________________
Bool_t AliAnalysisTaskPID::IsSecondaryWithStrangeMotherMC(AliMCEvent* mcEvent, Int_t partLabel)
{
//_____________________________________________________________________________
-Bool_t AliAnalysisTaskPID::SetParticleFractionHisto(const TH3D* hist, const Int_t species, const Bool_t sysError)
+Bool_t AliAnalysisTaskPID::SetParticleFractionHisto(const TH3D* hist, Int_t species, Bool_t sysError)
{
// Store a clone of hist (containing the particle fractions of the corresponding species with statistical error (sysError = kFALSE)
// or systematic error (sysError = kTRUE), respectively), internally
//_____________________________________________________________________________
-Bool_t AliAnalysisTaskPID::SetParticleFractionHistosFromFile(const TString filePathName, const Bool_t sysError)
+Bool_t AliAnalysisTaskPID::SetParticleFractionHistosFromFile(const TString filePathName, Bool_t sysError)
{
// Loads particle fractions for all species from the desired file and returns kTRUE on success.
// The maps are assumed to be of Type TH3D, to sit in the main directory and to have names
//_____________________________________________________________________________
-Int_t AliAnalysisTaskPID::GetRandomParticleTypeAccordingToParticleFractions(const Double_t trackPt, const Double_t jetPt,
- const Double_t centralityPercentile,
- const Bool_t smearByError,
- const Bool_t takeIntoAccountSysError) const
+Double_t AliAnalysisTaskPID::GetMaxEtaVariation(Double_t dEdxSplines)
+{
+ // Returns the maximum eta variation (i.e. deviation of eta correction factor from unity) for the
+ // given (spline) dEdx
+
+ if (dEdxSplines < 1. || !fhMaxEtaVariation) {
+ Printf("Error GetMaxEtaVariation: No map or invalid dEdxSplines (%f)!", dEdxSplines);
+ return 999.;
+ }
+
+ Int_t bin = fhMaxEtaVariation->GetXaxis()->FindFixBin(1. / dEdxSplines);
+
+ if (bin == 0)
+ bin = 1;
+ if (bin > fhMaxEtaVariation->GetXaxis()->GetNbins())
+ bin = fhMaxEtaVariation->GetXaxis()->GetNbins();
+
+ return fhMaxEtaVariation->GetBinContent(bin);
+}
+
+
+//_____________________________________________________________________________
+Int_t AliAnalysisTaskPID::GetRandomParticleTypeAccordingToParticleFractions(Double_t trackPt, Double_t jetPt,
+ Double_t centralityPercentile,
+ Bool_t smearByError,
+ Bool_t takeIntoAccountSysError) const
{
// Uses the stored histograms with the particle fractions to generate a random particle type according to these fractions.
// In case of problems (e.g. histo missing), AliPID::kUnknown is returned.
//_____________________________________________________________________________
-AliAnalysisTaskPID::ErrorCode AliAnalysisTaskPID::GenerateDetectorResponse(const AliAnalysisTaskPID::ErrorCode errCode,
- const Double_t mean, const Double_t sigma,
- Double_t* responses, const Int_t nResponses,
- const Bool_t usePureGaus)
+AliAnalysisTaskPID::ErrorCode AliAnalysisTaskPID::GenerateDetectorResponse(AliAnalysisTaskPID::ErrorCode errCode,
+ Double_t mean, Double_t sigma,
+ Double_t* responses, Int_t nResponses,
+ Bool_t usePureGaus)
{
// Generate detector response. If a previous generation was not successful or there is something wrong with this signal generation,
// the function will return kFALSE
printf("Eta cut: %.2f <= |eta| <= %.2f\n", GetEtaAbsCutLow(), GetEtaAbsCutUp());
printf("Phi' cut: %d\n", GetUsePhiCut());
printf("TPCCutMIGeo: %d\n", GetUseTPCCutMIGeo());
+ if (GetUseTPCCutMIGeo()) {
+ printf("GetCutGeo: %f\n", GetCutGeo());
+ printf("GetCutNcr: %f\n", GetCutNcr());
+ printf("GetCutNcl: %f\n", GetCutNcl());
+ }
+ printf("TPCnclCut: %d\n", GetUseTPCnclCut());
+ if (GetUseTPCnclCut()) {
+ printf("GetCutPureNcl: %d\n", GetCutPureNcl());
+ }
printf("\n");
printf("\n");
+ printf("Pile up rejection type: %d\n", (Int_t)fPileUpRejectionType);
+
+ printf("\n");
+
printf("Use MC-ID for signal generation: %d\n", GetUseMCidForGeneration());
printf("Use ITS: %d\n", GetUseITS());
printf("Use TOF: %d\n", GetUseTOF());
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));
printf("Do PID: %d\n", fDoPID);
printf("Do Efficiency: %d\n", fDoEfficiency);
printf("Do PtResolution: %d\n", fDoPtResolution);
+ printf("Do dEdxCheck: %d\n", fDoDeDxCheck);
+ printf("Do binZeroStudy: %d\n", fDoBinZeroStudy);
printf("\n");
// Print current settings for systematic studies.
printf("\n\nSettings for systematics for task %s:\n", GetName());
- printf("Splines:\t%f\n", GetSystematicScalingSplines());
+ printf("SplinesThr:\t%f\n", GetSystematicScalingSplinesThreshold());
+ printf("SplinesBelowThr:\t%f\n", GetSystematicScalingSplinesBelowThreshold());
+ printf("SplinesAboveThr:\t%f\n", GetSystematicScalingSplinesAboveThreshold());
printf("EtaCorrMomThr:\t%f\n", GetSystematicScalingEtaCorrectionMomentumThr());
printf("EtaCorrLowP:\t%f\n", GetSystematicScalingEtaCorrectionLowMomenta());
printf("EtaCorrHighP:\t%f\n", GetSystematicScalingEtaCorrectionHighMomenta());
- printf("SigmaPara:\t%f\n", GetSystematicScalingEtaSigmaPara());
+ printf("SigmaParaThr:\t%f\n", GetSystematicScalingEtaSigmaParaThreshold());
+ printf("SigmaParaBelowThr:\t%f\n", GetSystematicScalingEtaSigmaParaBelowThreshold());
+ printf("SigmaParaAboveThr:\t%f\n", GetSystematicScalingEtaSigmaParaAboveThreshold());
printf("MultCorr:\t%f\n", GetSystematicScalingMultCorrection());
+ printf("TOF mode: %d\n", GetTOFmode());
printf("\n\n");
}
if(fDebug > 1)
printf("File: %s, Line: %d: ProcessTrack\n", (char*)__FILE__, __LINE__);
- if (!fDoPID)
+ if (!fDoPID && !fDoDeDxCheck && !fDoPtResolution)
return kFALSE;
if(fDebug > 2)
// Momenta
//Double_t p = track->GetP();
- //Double_t pTPC = track->GetTPCmomentum();
+ const Double_t pTPC = track->GetTPCmomentum();
Double_t pT = track->Pt();
Double_t z = -1, xi = -1;
Double_t trackCharge = track->Charge();
// TPC signal
- Double_t dEdxTPC = fPIDResponse->IsTunedOnData() ? fPIDResponse->GetTPCsignalTunedOnData(track) : track->GetTPCsignal();
+ const Bool_t tuneOnDataTPC = fPIDResponse->IsTunedOnData() &&
+ ((fPIDResponse->GetTunedOnDataMask() & AliPIDResponse::kDetTPC) == AliPIDResponse::kDetTPC);
+ Double_t dEdxTPC = tuneOnDataTPC ? fPIDResponse->GetTPCsignalTunedOnData(track) : track->GetTPCsignal();
if (dEdxTPC <= 0) {
- Printf("Skipping track with strange dEdx value: dEdx %f, pTPC %f, eta %f, ncl %d\n", track->GetTPCsignal(), track->GetTPCmomentum(),
- track->Eta(), track->GetTPCsignalN());
+ if (fDebug > 1)
+ Printf("Skipping track with strange dEdx value: dEdx %f, pTPC %f, eta %f, ncl %d\n", track->GetTPCsignal(), pTPC,
+ track->Eta(), track->GetTPCsignalN());
return kFALSE;
}
+ // Completely remove tracks from light nuclei defined by el and pr <dEdx> as function of pT.
+ // A very loose cut to be sure not to throw away electrons and/or protons
+ Double_t nSigmaPr = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton);
+ Double_t nSigmaEl = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
+ if ((pTPC >= 0.3 && (nSigmaPr > 10 && nSigmaEl > 10)) ||
+ (pTPC < 0.3 && (nSigmaPr > 15 && nSigmaEl > 15))) {
+ if (fDebug > 1)
+ Printf("Skipping track which seems to be a light nucleus: dEdx %f, pTPC %f, pT %f, eta %f, ncl %d, nSigmaPr %f, nSigmaEl %f\n",
+ track->GetTPCsignal(), pTPC, pT, track->Eta(), track->GetTPCsignalN(), nSigmaPr, nSigmaEl);
+ return kFALSE;
+ }
Double_t dEdxEl, dEdxKa, dEdxPi, dEdxMu, dEdxPr;
dEdxPr = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
// Scale splines, if desired
- if (TMath::Abs(fSystematicScalingSplines - 1.0) > fgkEpsilon) {
- dEdxEl *= fSystematicScalingSplines;
- dEdxKa *= fSystematicScalingSplines;
- dEdxPi *= fSystematicScalingSplines;
- dEdxMu *= fTakeIntoAccountMuons ? fSystematicScalingSplines : 1.;
- dEdxPr *= fSystematicScalingSplines;
+ if ((TMath::Abs(fSystematicScalingSplinesBelowThreshold - 1.0) > fgkEpsilon) ||
+ (TMath::Abs(fSystematicScalingSplinesAboveThreshold - 1.0) > fgkEpsilon)) {
+
+ // Tune turn-on of correction for pions (not so relevant for the other species, since at very large momenta)
+ Double_t bg = 0;
+ Double_t scaleFactor = 1.;
+ Double_t usedSystematicScalingSplines = 1.;
+
+ bg = pTPC / AliPID::ParticleMass(AliPID::kElectron);
+ scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
+ usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
+ + fSystematicScalingSplinesAboveThreshold * scaleFactor;
+ dEdxEl *= usedSystematicScalingSplines;
+
+ bg = pTPC / AliPID::ParticleMass(AliPID::kKaon);
+ scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
+ usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
+ + fSystematicScalingSplinesAboveThreshold * scaleFactor;
+ dEdxKa *= usedSystematicScalingSplines;
+
+ bg = pTPC / AliPID::ParticleMass(AliPID::kPion);
+ scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
+ usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
+ + fSystematicScalingSplinesAboveThreshold * scaleFactor;
+ dEdxPi *= usedSystematicScalingSplines;
+
+ if (fTakeIntoAccountMuons) {
+ bg = pTPC / AliPID::ParticleMass(AliPID::kMuon);
+ scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
+ usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
+ + fSystematicScalingSplinesAboveThreshold * scaleFactor;
+ dEdxMu *= usedSystematicScalingSplines;
+ }
+
+
+ bg = pTPC / AliPID::ParticleMass(AliPID::kProton);
+ scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
+ usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
+ + fSystematicScalingSplinesAboveThreshold * scaleFactor;
+ dEdxPr *= usedSystematicScalingSplines;
}
// Get the eta correction factors for the (modified) expected dEdx
Double_t etaCorrEl = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxEl) : 1.;
Double_t etaCorrKa = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxKa) : 1.;
Double_t etaCorrPi = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxPi) : 1.;
- Double_t etaCorrMu = fTakeIntoAccountMuons && !fPIDResponse->UseTPCEtaCorrection() ?
+ Double_t etaCorrMu = fTakeIntoAccountMuons && fPIDResponse->UseTPCEtaCorrection() ?
fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxMu) : 1.;
Double_t etaCorrPr = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxPr) : 1.;
// Due to additional azimuthal effects, there is an additional eta dependence for low momenta which is not corrected successfully so far.
// One can assign a different (higher) systematic scale factor for this low-p region and a threshold which separates low- and high-p.
- // An ERF will be used to get (as a function of P_TPC) from one correction factor to the other within roughly 0.2 GeV/c
+ // An ERF will be used to get (as a function of p_TPC) from one correction factor to the other within roughly 0.2 GeV/c
Double_t usedSystematicScalingEtaCorrection = fSystematicScalingEtaCorrectionHighMomenta;
if (TMath::Abs(fSystematicScalingEtaCorrectionHighMomenta - fSystematicScalingEtaCorrectionLowMomenta) > fgkEpsilon) {
- const Double_t pTPC = track->GetTPCmomentum();
const Double_t fractionHighMomentumScaleFactor = 0.5 * (1. + TMath::Erf((pTPC - fSystematicScalingEtaCorrectionMomentumThr) / 0.1));
usedSystematicScalingEtaCorrection = fSystematicScalingEtaCorrectionLowMomenta * (1 - fractionHighMomentumScaleFactor)
+ fSystematicScalingEtaCorrectionHighMomenta * fractionHighMomentumScaleFactor;
}
+ Double_t maxEtaVariationEl = GetMaxEtaVariation(dEdxEl);
+ etaCorrEl = etaCorrEl * (1.0 + (usedSystematicScalingEtaCorrection - 1.) * (etaCorrEl - 1.0) / maxEtaVariationEl);
+
+ Double_t maxEtaVariationKa = GetMaxEtaVariation(dEdxKa);
+ etaCorrKa = etaCorrKa * (1.0 + (usedSystematicScalingEtaCorrection - 1.) * (etaCorrKa - 1.0) / maxEtaVariationKa);
+
+ Double_t maxEtaVariationPi = GetMaxEtaVariation(dEdxPi);
+ etaCorrPi = etaCorrPi * (1.0 + (usedSystematicScalingEtaCorrection - 1.) * (etaCorrPi - 1.0) / maxEtaVariationPi);
+
+ if (fTakeIntoAccountMuons) {
+ Double_t maxEtaVariationMu = GetMaxEtaVariation(dEdxMu);
+ etaCorrMu = etaCorrMu * (1.0 + (usedSystematicScalingEtaCorrection - 1.) * (etaCorrMu - 1.0) / maxEtaVariationMu);
+ }
+ else
+ etaCorrMu = 1.0;
+
+ Double_t maxEtaVariationPr = GetMaxEtaVariation(dEdxPr);
+ etaCorrPr = etaCorrPr * (1.0 + (usedSystematicScalingEtaCorrection - 1.) * (etaCorrPr - 1.0) / maxEtaVariationPr);
+
+
+ /*OLD
etaCorrEl = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrEl - 1.0);
etaCorrKa = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrKa - 1.0);
etaCorrPi = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrPi - 1.0);
etaCorrMu = fTakeIntoAccountMuons ? (1.0 + usedSystematicScalingEtaCorrection * (etaCorrMu - 1.0)) : 1.0;
etaCorrPr = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrPr - 1.0);
+ */
}
// Get the multiplicity correction factors for the (modified) expected dEdx
// This means there is no extra parameter for the multiplicitySigmaCorrFactor, but only for the sigma map itself.
// This is valid, since it appears only as a product. One has to assume a larger systematic shift in case of additional
// multiplicity dependence....
- Double_t sigmaRelEl = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
- / fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
- * fSystematicScalingEtaSigmaPara * multiplicityCorrSigmaEl;
- Double_t sigmaRelKa = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
- / fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
- * fSystematicScalingEtaSigmaPara * multiplicityCorrSigmaKa;
+ Bool_t doSigmaSystematics = (TMath::Abs(fSystematicScalingEtaSigmaParaBelowThreshold - 1.0) > fgkEpsilon) ||
+ (TMath::Abs(fSystematicScalingEtaSigmaParaAboveThreshold - 1.0) > fgkEpsilon);
- Double_t sigmaRelPi = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
- / fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
- * fSystematicScalingEtaSigmaPara * multiplicityCorrSigmaPi;
- Double_t sigmaRelMu = fTakeIntoAccountMuons ?
- fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
- / fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
- * fSystematicScalingEtaSigmaPara * multiplicityCorrSigmaMu
- : 999.;
+ Double_t dEdxElExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
+ kTRUE, kFALSE);
+ Double_t systematicScalingEtaSigmaParaEl = 1.;
+ if (doSigmaSystematics) {
+ Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxElExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
+ systematicScalingEtaSigmaParaEl = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
+ + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
+ }
+ Double_t sigmaRelEl = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
+ kTRUE, kFALSE)
+ / dEdxElExpected * systematicScalingEtaSigmaParaEl * multiplicityCorrSigmaEl;
+
+
+ Double_t dEdxKaExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
+ kTRUE, kFALSE);
+ Double_t systematicScalingEtaSigmaParaKa = 1.;
+ if (doSigmaSystematics) {
+ Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxKaExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
+ systematicScalingEtaSigmaParaKa = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
+ + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
+ }
+ Double_t sigmaRelKa = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
+ kTRUE, kFALSE)
+ / dEdxKaExpected * systematicScalingEtaSigmaParaKa * multiplicityCorrSigmaKa;
+
+
+ Double_t dEdxPiExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
+ kTRUE, kFALSE);
+ Double_t systematicScalingEtaSigmaParaPi = 1.;
+ if (doSigmaSystematics) {
+ Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxPiExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
+ systematicScalingEtaSigmaParaPi = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
+ + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
+ }
+ Double_t sigmaRelPi = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
+ kTRUE, kFALSE)
+ / dEdxPiExpected * systematicScalingEtaSigmaParaPi * multiplicityCorrSigmaPi;
+
+
+ Double_t sigmaRelMu = 999.;
+ if (fTakeIntoAccountMuons) {
+ Double_t dEdxMuExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault,
+ kTRUE, kFALSE);
+ Double_t systematicScalingEtaSigmaParaMu = 1.;
+ if (doSigmaSystematics) {
+ Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxMuExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
+ systematicScalingEtaSigmaParaMu = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
+ + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
+ }
+ sigmaRelMu = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
+ / dEdxMuExpected * systematicScalingEtaSigmaParaMu * multiplicityCorrSigmaMu;
+ }
- Double_t sigmaRelPr = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
- / fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
- * fSystematicScalingEtaSigmaPara * multiplicityCorrSigmaPr;
+
+ Double_t dEdxPrExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
+ kTRUE, kFALSE);
+ Double_t systematicScalingEtaSigmaParaPr = 1.;
+ if (doSigmaSystematics) {
+ Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxPrExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
+ systematicScalingEtaSigmaParaPr = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
+ + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
+ }
+ Double_t sigmaRelPr = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
+ kTRUE, kFALSE)
+ / dEdxPrExpected * systematicScalingEtaSigmaParaPr * multiplicityCorrSigmaPr;
// Now scale the (possibly modified) spline values with the (possibly modified) correction factors
dEdxEl *= etaCorrEl * multiplicityCorrEl;
}
else {
- // No systematic studies on expected signal - just take it as it comve from the TPCPIDResponse
+ // No systematic studies on expected signal - just take it as it comes from the TPCPIDResponse
dEdxEl = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
fPIDResponse->UseTPCEtaCorrection(),
fPIDResponse->UseTPCMultiplicityCorrection());
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) {
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__);
+
// Use probabilities to weigh the response generation for the different species.
// Also determine the most probable particle type.
Double_t prob[AliPID::kSPECIESC];
fPIDcombined->ComputeProbabilities(track, fPIDResponse, prob);
+ // If muons are not to be taken into account, just set their probability to zero and normalise the remaining probabilities later
+ if (!fTakeIntoAccountMuons)
+ prob[AliPID::kMuon] = 0;
+
+ // Priors might cause problems in case of mismatches, i.e. mismatch particles will be identified as pions.
+ // Can be easily caught for low p_TPC by just setting pion (and muon) probs to zero, if dEdx is larger than
+ // expected dEdx of electrons and kaons
+ if (pTPC < 1. && dEdxTPC > dEdxEl && dEdxTPC > dEdxKa) {
+ prob[AliPID::kMuon] = 0;
+ prob[AliPID::kPion] = 0;
+ }
+
+
// Bug: One can set the number of species for PIDcombined, but PIDcombined will call PIDresponse, which writes without testing
// the probs for kSPECIESC (including light nuclei) into the array.
// In this case, when only kSPECIES are considered, the probabilities have to be rescaled!
- for (Int_t i = AliPID::kSPECIES; i < AliPID::kSPECIESC; i++)
- prob[i] = 0;
- // If muons are not to be taken into account, just set their probability to zero and normalise the remaining probabilities
- if (!fTakeIntoAccountMuons)
- prob[AliPID::kMuon] = 0;
+ // Exceptions:
+ // 1. If the most probable PID is a light nucleus and above expected dEdx + 5 sigma of el to be sure not to mix up things in the
+ // high-pT region (also contribution there completely negligable!)
+ // 2. If dEdx is larger than the expected dEdx + 5 sigma of el and pr, it is most likely a light nucleus
+ // (if dEdx is more than 5 sigma away from expectation, flat TPC probs are set... and the light nuclei splines are not
+ // too accurate)
+ // In these cases:
+ // The highest expected dEdx is either for pr or for el. Assign 100% probability to the species with highest expected dEdx then.
+ // In all other cases: Throw away light nuclei probs and rescale other probs to 100%.
+
+ // Find most probable species for the ID
+ Int_t maxIndex = TMath::LocMax(AliPID::kSPECIESC, prob);
+
+ if ((prob[maxIndex] > 0 && maxIndex >= AliPID::kSPECIES && dEdxTPC > dEdxEl + 5. * sigmaEl) ||
+ (dEdxTPC > dEdxEl + 5. * sigmaEl && dEdxTPC > dEdxPr + 5. * sigmaPr)) {
+ for (Int_t i = 0; i < AliPID::kSPECIESC; i++)
+ prob[i] = 0;
+
+ if (dEdxEl > dEdxPr) {
+ prob[AliPID::kElectron] = 1.;
+ maxIndex = AliPID::kElectron;
+ }
+ else {
+ prob[AliPID::kProton] = 1.;
+ maxIndex = AliPID::kProton;
+ }
+ }
+ else {
+ for (Int_t i = AliPID::kSPECIES; i < AliPID::kSPECIESC; i++)
+ prob[i] = 0;
+
+ Double_t probSum = 0.;
+ for (Int_t i = 0; i < AliPID::kSPECIES; i++)
+ probSum += prob[i];
+
+ if (probSum > 0) {
+ for (Int_t i = 0; i < AliPID::kSPECIES; i++)
+ prob[i] /= probSum;
+ }
+
+ // If all probabilities are equal (should only happen if no priors are used), set most probable PID to pion
+ // because LocMax returns just the first maximal value (i.e. AliPID::kElectron)
+ if (TMath::Abs(prob[maxIndex] - 1. / AliPID::kSPECIES) < 1e-5)
+ maxIndex = AliPID::kPion;
+
+ // In all other cases, the maximum remains untouched from the scaling (and is < AliPID::kSPECIES by construction)
+ }
- Double_t probSum = 0.;
- for (Int_t i = 0; i < AliPID::kSPECIES; i++)
- probSum += prob[i];
+ if (fDoDeDxCheck) {
+ // Generate the expected responses in DeltaPrime and translate to real dEdx. Only take responses for exactly that species
+ // (i.e. with pre-PID)
+
+ Int_t numGenEntries = fgkMaxNumGenEntries; // fgkMaxNumGenEntries = 500
+
+ ErrorCode errCode = kNoErrors;
- if (probSum > 0) {
- for (Int_t i = 0; i < AliPID::kSPECIES; i++)
- prob[i] /= probSum;
+ if(fDebug > 2)
+ printf("File: %s, Line: %d: ProcessTrack -> Generate Responses for dEdx check\n", (char*)__FILE__, __LINE__);
+
+ // Electrons
+ errCode = GenerateDetectorResponse(errCode, 1., sigmaEl / dEdxEl, fGenRespElDeltaPrimeEl, numGenEntries);
+
+ // Kaons
+ errCode = GenerateDetectorResponse(errCode, 1., sigmaKa / dEdxKa, fGenRespKaDeltaPrimeKa, numGenEntries);
+
+ // Pions
+ errCode = GenerateDetectorResponse(errCode, 1., sigmaPi / dEdxPi, fGenRespPiDeltaPrimePi, numGenEntries);
+
+ // Protons
+ errCode = GenerateDetectorResponse(errCode, 1., sigmaPr / dEdxPr, fGenRespPrDeltaPrimePr, numGenEntries);
+
+ if (errCode == kNoErrors) {
+ Double_t genEntry[kDeDxCheckNumAxes];
+ genEntry[kDeDxCheckJetPt] = jetPt;
+ genEntry[kDeDxCheckEtaAbs] = TMath::Abs(track->Eta());
+ genEntry[kDeDxCheckP] = pTPC;
+
+
+ for (Int_t n = 0; n < numGenEntries; n++) {
+ // If no MC info is available or shall not be used, use weighting with priors to generate entries for the different species
+ Double_t rnd = fRandom->Rndm(); // Produce uniformly distributed floating point in ]0, 1]
+
+ // Consider generated response as originating from...
+ if (rnd <= prob[AliPID::kElectron]) {
+ genEntry[kDeDxCheckPID] = 0; // ... an electron
+ genEntry[kDeDxCheckDeDx] = fGenRespElDeltaPrimeEl[n] * dEdxEl;
+ }
+ else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon]) {
+ genEntry[kDeDxCheckPID] = 1; // ... a kaon
+ genEntry[kDeDxCheckDeDx] = fGenRespKaDeltaPrimeKa[n] * dEdxKa;
+ }
+ else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon] + prob[AliPID::kMuon] + prob[AliPID::kPion]) {
+ genEntry[kDeDxCheckPID] = 2; // ... a pion (or a muon)
+ genEntry[kDeDxCheckDeDx] = fGenRespPiDeltaPrimePi[n] * dEdxPi;
+ }
+ else {
+ genEntry[kDeDxCheckPID] = 3; // ... a proton
+ genEntry[kDeDxCheckDeDx] = fGenRespPrDeltaPrimePr[n] * dEdxPr;
+ }
+
+ fDeDxCheck->Fill(genEntry);
+ }
+ }
+
+ if(fDebug > 2)
+ printf("File: %s, Line: %d: ProcessTrack -> Generate Responses for dEdx check done\n", (char*)__FILE__, __LINE__);
}
- if (!isMC) {
- // If there is no MC information, take the most probable species for the ID
- Float_t max = 0.;
- Int_t maxIndex = -1;
- for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
- if (prob[i] > max) {
- max = prob[i];
- maxIndex = i;
+ if (fDoPtResolution) {
+ // Check shared clusters, which is done together with the pT resolution
+ Double_t qaEntry[kQASharedClsNumAxes];
+ qaEntry[kQASharedClsJetPt] = jetPt;
+ qaEntry[kQASharedClsPt] = pT;
+ qaEntry[kDeDxCheckP] = pTPC;
+ qaEntry[kQASharedClsNumSharedCls] = track->GetTPCSharedMapPtr()->CountBits();
+
+ Int_t iRowInd = -1;
+ // iRowInd == -1 for "all rows w/o multiple counting"
+ qaEntry[kQASharedClsPadRow] = iRowInd;
+ fQASharedCls->Fill(qaEntry);
+
+ // Fill hist for every pad row with shared cluster
+ for (iRowInd = 0; iRowInd < 159; iRowInd++) {
+ if (track->GetTPCSharedMapPtr()->TestBitNumber(iRowInd)) {
+ qaEntry[kQASharedClsPadRow] = iRowInd;
+ fQASharedCls->Fill(qaEntry);
}
}
-
+ }
+
+ if (!fDoPID)
+ return kTRUE;
+
+ if (!isMC) {
+ // Clean up the most probable PID for low momenta (not an issue for the probabilities, since fractions small,
+ // but to get clean (and nice looking) templates (if most probable PID is used instead of generated responses), it should be done).
+ // Problem: If more than 5 sigma away (and since priors also included), most probable PID for dEdx around protons could be pions.
+ // Idea: Only clean selection for K and p possible. So, just take for the most probable PID the probs from TPC only calculated
+ // by hand without restriction to 5 sigma and without priors. This will not affect the template generation, but afterwards one
+ // can replace the K and p templates with the most probable PID distributions, so all other species templates remain the same.
+ // I.e. it doesn't hurt if the most probable PID is distributed incorrectly between el, pi, mu.
+ Bool_t maxIndexSetManually = kFALSE;
+
+ if (pTPC < 1.) {
+ Double_t probManualTPC[AliPID::kSPECIES];
+ for (Int_t i = 0; i < AliPID::kSPECIES; i++)
+ probManualTPC[i] = 0;
+
+ probManualTPC[AliPID::kElectron] = TMath::Exp(-0.5*(dEdxTPC-dEdxEl)*(dEdxTPC-dEdxEl)/(sigmaEl*sigmaEl))/sigmaEl;
+ // Muons are not important here, just ignore and save processing time
+ probManualTPC[AliPID::kMuon] = 0.;//TMath::Exp(-0.5*(dEdxTPC-dEdxMu)*(dEdxTPC-dEdxMu)/(sigmaMu*sigmaMu))/sigmaMu;
+ probManualTPC[AliPID::kPion] = TMath::Exp(-0.5*(dEdxTPC-dEdxPi)*(dEdxTPC-dEdxPi)/(sigmaPi*sigmaPi))/sigmaPi;
+ probManualTPC[AliPID::kKaon] = TMath::Exp(-0.5*(dEdxTPC-dEdxKa)*(dEdxTPC-dEdxKa)/(sigmaKa*sigmaKa))/sigmaKa;
+ probManualTPC[AliPID::kProton] = TMath::Exp(-0.5*(dEdxTPC-dEdxPr)*(dEdxTPC-dEdxPr)/(sigmaPr*sigmaPr))/sigmaPr;
+
+ const Int_t maxIndexManualTPC = TMath::LocMax(AliPID::kSPECIES, probManualTPC);
+ // Should always be > 0, but in case the dEdx is far off such that the machine accuracy sets every prob to zero,
+ // better take the "old" result
+ if (probManualTPC[maxIndexManualTPC] > 0.)
+ maxIndex = maxIndexManualTPC;
+
+ maxIndexSetManually = kTRUE;
+ }
+
+
// Translate from AliPID numbering to numbering of this class
- if (max > 0) {
+ if (prob[maxIndex] > 0 || maxIndexSetManually) {
if (maxIndex == AliPID::kElectron)
binMC = 0;
else if (maxIndex == AliPID::kKaon)
pTPC = temp;
*/
+ TOFpidInfo tofPIDinfo = GetTOFType(track, fTOFmode);
Double_t entry[fStoreAdditionalJetInformation ? kDataNumAxes : kDataNumAxes - fgkNumJetAxes];
entry[kDataMCID] = binMC;
}
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);
}
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;
numGenEntries *= 2;
+
// Do not generate more entries than available in memory!
if (numGenEntries > fgkMaxNumGenEntries)// fgkMaxNumGenEntries = 1000
numGenEntries = fgkMaxNumGenEntries;
-
+ */
+
+
ErrorCode errCode = kNoErrors;
if(fDebug > 2)
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());
return kFALSE;// Skip generated response in case of error
}
}
// 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);
}
//_____________________________________________________________________________
-AliAnalysisTaskPID::ErrorCode AliAnalysisTaskPID::SetParamsForConvolutedGaus(const Double_t gausMean, const Double_t gausSigma)
+AliAnalysisTaskPID::ErrorCode AliAnalysisTaskPID::SetParamsForConvolutedGaus(Double_t gausMean, Double_t gausSigma)
{
// Set parameters for convoluted gauss using parameters for a pure gaussian.
// If SetConvolutedGaussLambdaParameter has not been called before to initialise the translation parameters,
rangeEnd = fConvolutedGausDeltaPrime->GetX(maximumFraction, upBoundSearchBoundLow, upBoundSearchBoundUp, (maxX < 0.4) ? 1e-5 : 0.001);
fConvolutedGausDeltaPrime->SetRange(rangeStart,rangeEnd);
- fConvolutedGausDeltaPrime->SetNpx(fhPIDdataAll->GetAxis(kDataDeltaPrimeSpecies)->FindBin(rangeEnd)
- - fhPIDdataAll->GetAxis(kDataDeltaPrimeSpecies)->FindBin(rangeStart) + 1);
+ fConvolutedGausDeltaPrime->SetNpx(fDeltaPrimeAxis->FindFixBin(rangeEnd) - fDeltaPrimeAxis->FindFixBin(rangeStart) + 1);
+
+ fConvolutedGausDeltaPrime->SetNpx(fDeltaPrimeAxis->FindFixBin(rangeEnd) - fDeltaPrimeAxis->FindFixBin(rangeStart) + 1);
+ //fConvolutedGausDeltaPrime->SetNpx(fhPIDdataAll->GetAxis(kDataDeltaPrimeSpecies)->FindFixBin(rangeEnd)
+ // - fhPIDdataAll->GetAxis(kDataDeltaPrimeSpecies)->FindFixBin(rangeStart) + 1);
//fConvolutedGausDeltaPrime->SetNpx((rangeEnd - rangeStart) / fDeltaPrimeBinWidth + 1);
}
hist->GetAxis(kGenSelectSpecies)->SetBinLabel(3, "#pi");
hist->GetAxis(kGenSelectSpecies)->SetBinLabel(4, "p");
- hist->GetAxis(kGenPt)->SetTitle("P_{T} (GeV/c)");
+ hist->GetAxis(kGenPt)->SetTitle("p_{T} (GeV/c)");
hist->GetAxis(kGenDeltaPrimeSpecies)->SetTitle("TPC #Delta'_{species} (arb. unit)");
- hist->GetAxis(kGenCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
+ hist->GetAxis(kGenCentrality)->SetTitle(Form("Centrality Percentile (%s)", GetPPCentralityEstimator().Data()));
if (fStoreAdditionalJetInformation) {
- hist->GetAxis(kGenJetPt)->SetTitle("P_{T}^{jet} (GeV/c)");
+ hist->GetAxis(kGenJetPt)->SetTitle("p_{T}^{jet} (GeV/c)");
- hist->GetAxis(kGenZ)->SetTitle("z = P_{T}^{track} / P_{T}^{jet}");
+ hist->GetAxis(kGenZ)->SetTitle("z = p_{T}^{track} / p_{T}^{jet}");
- hist->GetAxis(kGenXi)->SetTitle("#xi = ln(P_{T}^{jet} / P_{T}^{track})");
+ hist->GetAxis(kGenXi)->SetTitle("#xi = ln(p_{T}^{jet} / p_{T}^{track})");
}
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");
}
// Set axes titles
hist->GetAxis(kGenYieldMCID)->SetTitle("MC PID");
- hist->GetAxis(kGenYieldPt)->SetTitle("P_{T}^{gen} (GeV/c)");
- hist->GetAxis(kGenYieldCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
+ hist->GetAxis(kGenYieldPt)->SetTitle("p_{T}^{gen} (GeV/c)");
+ hist->GetAxis(kGenYieldCentrality)->SetTitle(Form("Centrality Percentile (%s)", GetPPCentralityEstimator().Data()));
if (fStoreAdditionalJetInformation) {
- hist->GetAxis(kGenYieldJetPt)->SetTitle("P_{T}^{jet, gen} (GeV/c)");
+ hist->GetAxis(kGenYieldJetPt)->SetTitle("p_{T}^{jet, gen} (GeV/c)");
- hist->GetAxis(kGenYieldZ)->SetTitle("z = P_{T}^{track} / P_{T}^{jet}");
+ hist->GetAxis(kGenYieldZ)->SetTitle("z = p_{T}^{track} / p_{T}^{jet}");
- hist->GetAxis(kGenYieldXi)->SetTitle("#xi = ln(P_{T}^{jet} / P_{T}^{track})");
+ hist->GetAxis(kGenYieldXi)->SetTitle("#xi = ln(p_{T}^{jet} / p_{T}^{track})");
}
hist->GetAxis(GetIndexOfChargeAxisGenYield())->SetTitle("Charge (e_{0})");
hist->GetAxis(kDataSelectSpecies)->SetBinLabel(3, "#pi");
hist->GetAxis(kDataSelectSpecies)->SetBinLabel(4, "p");
- hist->GetAxis(kDataPt)->SetTitle("P_{T} (GeV/c)");
+ hist->GetAxis(kDataPt)->SetTitle("p_{T} (GeV/c)");
hist->GetAxis(kDataDeltaPrimeSpecies)->SetTitle("TPC #Delta'_{species} (arb. unit)");
- hist->GetAxis(kDataCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
+ hist->GetAxis(kDataCentrality)->SetTitle(Form("Centrality Percentile (%s)", GetPPCentralityEstimator().Data()));
if (fStoreAdditionalJetInformation) {
- hist->GetAxis(kDataJetPt)->SetTitle("P_{T}^{jet} (GeV/c)");
+ hist->GetAxis(kDataJetPt)->SetTitle("p_{T}^{jet} (GeV/c)");
- hist->GetAxis(kDataZ)->SetTitle("z = P_{T}^{track} / P_{T}^{jet}");
+ hist->GetAxis(kDataZ)->SetTitle("z = p_{T}^{track} / p_{T}^{jet}");
- hist->GetAxis(kDataXi)->SetTitle("#xi = ln(P_{T}^{jet} / P_{T}^{track})");
+ hist->GetAxis(kDataXi)->SetTitle("#xi = ln(p_{T}^{jet} / p_{T}^{track})");
}
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(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");
+}
+
+
+//________________________________________________________________________
+void AliAnalysisTaskPID::SetUpPtResHist(THnSparse* hist, Double_t* binsPt, Double_t* binsJetPt, Double_t* binsCent) const
+{
+ // Sets bin limits for axes which are not standard binned and the axes titles.
+
+ hist->SetBinEdges(kPtResJetPt, binsJetPt);
+ hist->SetBinEdges(kPtResGenPt, binsPt);
+ hist->SetBinEdges(kPtResRecPt, binsPt);
+ hist->SetBinEdges(kPtResCentrality, binsCent);
- 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");
+ // Set axes titles
+ hist->GetAxis(kPtResJetPt)->SetTitle("p_{T}^{jet, rec} (GeV/c)");
+ hist->GetAxis(kPtResGenPt)->SetTitle("p_{T}^{gen} (GeV/c)");
+ hist->GetAxis(kPtResRecPt)->SetTitle("p_{T}^{rec} (GeV/c)");
- 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(kPtResCharge)->SetTitle("Charge (e_{0})");
+ hist->GetAxis(kPtResCentrality)->SetTitle(Form("Centrality Percentile (%s)", GetPPCentralityEstimator().Data()));
+}
+
+
+//________________________________________________________________________
+void AliAnalysisTaskPID::SetUpSharedClsHist(THnSparse* hist, Double_t* binsPt, Double_t* binsJetPt) const
+{
+ // Sets bin limits for axes which are not standard binned and the axes titles.
- hist->GetAxis(5)->SetTitle("TPC #Delta_{species} (arb. unit)");
+ hist->SetBinEdges(kQASharedClsJetPt, binsJetPt);
+ hist->SetBinEdges(kQASharedClsPt, binsPt);
- hist->GetAxis(6)->SetTitle("TPC #Delta'_{species} (arb. unit)");
+ // Set axes titles
+ hist->GetAxis(kQASharedClsJetPt)->SetTitle("#it{p}_{T}^{jet} (GeV/#it{c})");
+ hist->GetAxis(kQASharedClsPt)->SetTitle("#it{p}_{T} (GeV/#it{c})");
+ hist->GetAxis(kQASharedClsNumSharedCls)->SetTitle("#it{N}_{shared}^{cls}");
+ hist->GetAxis(kQASharedClsPadRow)->SetTitle("Pad row");
- hist->GetAxis(7)->SetTitle("#Delta TOF_{species} (ps)");
- */
}
//________________________________________________________________________
-void AliAnalysisTaskPID::SetUpPtResHist(THnSparse* hist, Double_t* binsPt, Double_t* binsJetPt) const
+void AliAnalysisTaskPID::SetUpDeDxCheckHist(THnSparse* hist, const Double_t* binsPt, const Double_t* binsJetPt, const Double_t* binsEtaAbs) const
{
// Sets bin limits for axes which are not standard binned and the axes titles.
+ hist->SetBinEdges(kDeDxCheckP, binsPt);
+ hist->SetBinEdges(kDeDxCheckJetPt, binsJetPt);
+ hist->SetBinEdges(kDeDxCheckEtaAbs, binsEtaAbs);
- hist->SetBinEdges(kPtResJetPt, binsJetPt);
- hist->SetBinEdges(kPtResGenPt, binsPt);
- hist->SetBinEdges(kPtResRecPt, binsPt);
// Set axes titles
- hist->GetAxis(kPtResJetPt)->SetTitle("P_{T}^{jet, rec} (GeV/c)");
- hist->GetAxis(kPtResGenPt)->SetTitle("P_{T}^{gen} (GeV/c)");
- hist->GetAxis(kPtResRecPt)->SetTitle("P_{T}^{rec} (GeV/c)");
+ hist->GetAxis(kDeDxCheckPID)->SetTitle("PID");
+ hist->GetAxis(kDeDxCheckPID)->SetBinLabel(1, "e");
+ hist->GetAxis(kDeDxCheckPID)->SetBinLabel(2, "K");
+ hist->GetAxis(kDeDxCheckPID)->SetBinLabel(3, "#pi");
+ hist->GetAxis(kDeDxCheckPID)->SetBinLabel(4, "p");
+
+
+ hist->GetAxis(kDeDxCheckJetPt)->SetTitle("p_{T}^{jet} (GeV/c)");
+ hist->GetAxis(kDeDxCheckEtaAbs)->SetTitle("|#eta|");
+ hist->GetAxis(kDeDxCheckP)->SetTitle("p_{TPC} (GeV/c)");
+ hist->GetAxis(kDeDxCheckDeDx)->SetTitle("TPC dE/dx (arb. unit)");
}
+
+
+//________________________________________________________________________
+void AliAnalysisTaskPID::SetUpBinZeroStudyHist(THnSparse* hist, const Double_t* binsCent, const Double_t* binsPt) const
+{
+ // Sets bin limits for axes which are not standard binned and the axes titles.
+ hist->SetBinEdges(kBinZeroStudyCentrality, binsCent);
+ hist->SetBinEdges(kBinZeroStudyGenPt, binsPt);
+
+ // Set axes titles
+ hist->GetAxis(kBinZeroStudyCentrality)->SetTitle(Form("Centrality Percentile (%s)", GetPPCentralityEstimator().Data()));
+ hist->GetAxis(kBinZeroStudyGenEta)->SetTitle("#it{#eta}^{gen}");
+ hist->GetAxis(kBinZeroStudyGenPt)->SetTitle("#it{p}_{T}^{gen} (GeV/#it{c})");
+}
\ No newline at end of file