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