ClassImp(AliPIDResponse);
+Float_t AliPIDResponse::fgTOFmismatchProb = 0.0;
+
AliPIDResponse::AliPIDResponse(Bool_t isMC/*=kFALSE*/) :
TNamed("PIDResponse","PIDResponse"),
fITSResponse(isMC),
fHMPIDPIDParams(NULL),
fEMCALPIDParams(NULL),
fCurrentEvent(NULL),
-fCurrCentrality(0.0)
+fCurrCentrality(0.0),
+fBeamTypeNum(kPP)
{
//
// default ctor
fHMPIDPIDParams(NULL),
fEMCALPIDParams(NULL),
fCurrentEvent(NULL),
-fCurrCentrality(0.0)
+fCurrCentrality(0.0),
+fBeamTypeNum(kPP)
{
//
// copy ctor
fIsMC=other.fIsMC;
fCachePID=other.fCachePID;
fBeamType="PP";
+ fBeamTypeNum=kPP;
fLHCperiod="";
fMCperiodTPC="";
fMCperiodUser=other.fMCperiodUser;
}
// Set up TPC multiplicity for PbPb
- //TODO Will NOT give the desired number for AODs -> Needs new variable/function in future.
- // Fatal, if AOD event and correction enabled
- //printf("DETECTED class: %s (%d)\n\n\n\n", event->IsA()->GetName(), fUseTPCMultiplicityCorrection);//TODO
- if (fUseTPCMultiplicityCorrection && strcmp(event->IsA()->GetName(), "AliESDEvent") != 0) {
- AliFatal("TPC multiplicity correction is enabled, but will NOT work for AOD events, only for ESD => Disabled multiplicity correction!");
- fUseTPCMultiplicityCorrection = kFALSE;
+ if (fUseTPCMultiplicityCorrection) {
+ Int_t numESDtracks = event->GetNumberOfESDTracks();
+ if (numESDtracks < 0) {
+ AliError("Cannot obtain event multiplicity (number of ESD tracks < 0). If you are using AODs, this might be a too old production. Please disable the multiplicity correction to get a reliable PID result!");
+ numESDtracks = 0;
+ }
+ fTPCResponse.SetCurrentEventMultiplicity(numESDtracks);
}
-
- if (fUseTPCMultiplicityCorrection)
- fTPCResponse.SetCurrentEventMultiplicity(event->GetNumberOfTracks());
else
fTPCResponse.SetCurrentEventMultiplicity(0);
fBeamType="";
fBeamType="PP";
+ fBeamTypeNum=kPP;
Bool_t hasProdInfo=(fCurrentFile.BeginsWith("LHC"));
- TPRegexp reg(".*(LHC1[1-3][a-z]+[0-9]+[a-z_]*)/.*");
+ TPRegexp reg(".*(LHC1[1-3][a-z]+[0-9]+[a-z_]*)[/_].*");
if (hasProdInfo) reg=TPRegexp("LHC1[1-2][a-z]+[0-9]+[a-z_]*");
- TPRegexp reg12a17("LHC1[2-3][a-z]");
+ TPRegexp reg12a17("LHC1[2-4][a-z]");
//find the period by run number (UGLY, but not stored in ESD and AOD... )
if (fRun>=114737&&fRun<=117223) { fLHCperiod="LHC10B"; fMCperiodTPC="LHC10D1"; }
// exception for 13d2 and later
if (fCurrentAliRootRev >= 62714) fMCperiodTPC="LHC13D2";
fBeamType="PBPB";
+ fBeamTypeNum=kPBPB;
}
else if (fRun>=139847&&fRun<=146974) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
//TODO: periods 11B (146975-150721), 11C (150722-155837) are not yet treated assume 11d for the moment
fLHCperiod="LHC11H";
fMCperiodTPC="LHC11A10";
fBeamType="PBPB";
+ fBeamTypeNum=kPBPB;
if (reg12a17.MatchB(fCurrentFile)) fMCperiodTPC="LHC12A17";
}
- if (fRun>=170719 && fRun<=177311) { fLHCperiod="LHC12A"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
+ if (fRun>=170719 && fRun<=177311) { fLHCperiod="LHC12A"; fBeamType="PP"; fBeamTypeNum=kPP;/*fMCperiodTPC="";*/ }
// for the moment use LHC12b parameters up to LHC12d
- if (fRun>=177312 /*&& fRun<=179356*/) { fLHCperiod="LHC12B"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
-// if (fRun>=179357 && fRun<=183173) { fLHCperiod="LHC12C"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
-// if (fRun>=183174 && fRun<=186345) { fLHCperiod="LHC12D"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
-// if (fRun>=186346 && fRun<=186635) { fLHCperiod="LHC12E"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
-
-// if (fRun>=186636 && fRun<=188166) { fLHCperiod="LHC12F"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
-// if (fRun >= 188167 && fRun <= 188355 ) { fLHCperiod="LHC12G"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
-// if (fRun >= 188356 && fRun <= 188503 ) { fLHCperiod="LHC12G"; fBeamType="PPB"; /*fMCperiodTPC="";*/ }
+ if (fRun>=177312 /*&& fRun<=179356*/) { fLHCperiod="LHC12B"; fBeamType="PP";fBeamTypeNum=kPP; /*fMCperiodTPC="";*/ }
+// if (fRun>=179357 && fRun<=183173) { fLHCperiod="LHC12C"; fBeamType="PP"; fBeamTypeNum=kPP;/*fMCperiodTPC="";*/ }
+// if (fRun>=183174 && fRun<=186345) { fLHCperiod="LHC12D"; fBeamType="PP"; fBeamTypeNum=kPP;/*fMCperiodTPC="";*/ }
+// if (fRun>=186346 && fRun<=186635) { fLHCperiod="LHC12E"; fBeamType="PP"; fBeamTypeNum=kPP;/*fMCperiodTPC="";*/ }
+
+// if (fRun>=186636 && fRun<=188166) { fLHCperiod="LHC12F"; fBeamType="PP"; fBeamTypeNum=kPP;/*fMCperiodTPC="";*/ }
+// if (fRun >= 188167 && fRun <= 188355 ) { fLHCperiod="LHC12G"; fBeamType="PP"; fBeamTypeNum=kPP;/*fMCperiodTPC="";*/ }
+// if (fRun >= 188356 && fRun <= 188503 ) { fLHCperiod="LHC12G"; fBeamType="PPB"; fBeamTypeNum=kPPB;/*fMCperiodTPC="";*/ }
// for the moment use 12g parametrisation for all full gain runs (LHC12e+)
- if (fRun >= 186346 && fRun < 194480) { fLHCperiod="LHC12G"; fBeamType="PPB"; fMCperiodTPC="LHC12G"; }
+ if (fRun >= 186346 && fRun < 194480) { fLHCperiod="LHC12G"; fBeamType="PPB";fBeamTypeNum=kPPB; fMCperiodTPC="LHC12G"; }
// New parametrisation for 2013 pPb runs
if (fRun >= 194480) {
fLHCperiod="LHC13B";
fBeamType="PPB";
+ fBeamTypeNum=kPPB;
fMCperiodTPC="LHC12G";
if (fCurrentAliRootRev >= 61605)
}
//exception new pp MC productions from 2011 (11a periods have 10f6a splines!)
- if (fBeamType=="PP" && reg.MatchB(fCurrentFile) && !fCurrentFile.Contains("LHC11a")) { fMCperiodTPC="LHC11B2"; fBeamType="PP"; }
+ if (fBeamType=="PP" && reg.MatchB(fCurrentFile) && !fCurrentFile.Contains("LHC11a")) { fMCperiodTPC="LHC11B2"; fBeamType="PP";fBeamTypeNum=kPP; }
// exception for 11f1
if (fCurrentFile.Contains("LHC11f1")) fMCperiodTPC="LHC11F1";
// exception for 12f1a, 12f1b and 12i3
//
- // Setup multiplicity correction
+ // Setup multiplicity correction (only used for non-pp collisions)
//
- if (fUseTPCMultiplicityCorrection && !(fBeamType.CompareTo("PP") == 0)) {
+
+ const Bool_t isPP = (fBeamType.CompareTo("PP") == 0);
+
+ // 2013 pPb data taking at low luminosity
+ const Bool_t isPPb2013LowLuminosity = period.Contains("LHC13B") || period.Contains("LHC13C") || period.Contains("LHC13D");
+ // PbPb 2010, period 10h.pass2
+ //TODO Needs further development const Bool_t is10hpass2 = period.Contains("LHC10H") && recopass == 2;
+
+
+ // In case of MC without(!) tune on data activated for the TPC, don't use the multiplicity correction for the moment
+ Bool_t isMCandNotTPCtuneOnData = fIsMC && !(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC));
+
+ // If correction is available, but disabled (highly NOT recommended!), print warning
+ if (!fUseTPCMultiplicityCorrection && !isPP && !isMCandNotTPCtuneOnData) {
+ //TODO: Needs further development if (is10hpass2 || isPPb2013LowLuminosity) {
+ if (isPPb2013LowLuminosity) {
+ AliWarning("Mulitplicity correction disabled, but correction parameters for this period exist. It is highly recommended to use enable the correction. Otherwise the splines might be off!");
+ }
+ }
+
+ if (fUseTPCMultiplicityCorrection && !isPP && !isMCandNotTPCtuneOnData) {
AliInfo("Multiplicity correction enabled!");
//TODO After testing, load parameters from outside
- /*TODO now correction for MC
+ /*TODO no correction for MC
if (period.Contains("LHC11A10")) {//LHC11A10A
AliInfo("Using multiplicity correction parameters for 11a10!");
fTPCResponse.SetParameterMultiplicityCorrection(0, 6.90133e-06);
fTPCResponse.SetParameterMultiplicitySigmaCorrection(2, -6.36337e-01);
fTPCResponse.SetParameterMultiplicitySigmaCorrection(3, 1.13479e-02);
}
- else*/ if (period.Contains("LHC13B") || period.Contains("LHC13C") || period.Contains("LHC13D")) {// 2013 pPb data taking at low luminosity
- AliInfo("Using multiplicity correction parameters for 13b.pass2!");
+ else*/ if (isPPb2013LowLuminosity) {// 2013 pPb data taking at low luminosity
+ AliInfo("Using multiplicity correction parameters for 13b.pass2 (at least also valid for 13{c,d} and pass 3)!");
fTPCResponse.SetParameterMultiplicityCorrection(0, -5.906e-06);
fTPCResponse.SetParameterMultiplicityCorrection(1, -5.064e-04);
fTPCResponse.SetParameterMultiplicitySigmaCorrection(3, 1.47931e-02);
*/
}
- else if (period.Contains("LHC10H") && recopass == 2) {
+ /*TODO: Needs further development
+ else if (is10hpass2) {
AliInfo("Using multiplicity correction parameters for 10h.pass2!");
fTPCResponse.SetParameterMultiplicityCorrection(0, 3.21636e-07);
fTPCResponse.SetParameterMultiplicityCorrection(1, -6.65876e-04);
fTPCResponse.SetParameterMultiplicitySigmaCorrection(2, -3.20788e-01);
fTPCResponse.SetParameterMultiplicitySigmaCorrection(3, 1.07345e-02);
}
+ */
else {
AliError(Form("Multiplicity correction is enabled, but no multiplicity correction parameters have been found for period %s.pass%d -> Mulitplicity correction DISABLED!", period.Data(), recopass));
fUseTPCMultiplicityCorrection = kFALSE;
// directly and use it for calculations - which should still give valid results, even if
// the multiplicity correction is explicitely enabled in such expert calls.
+ TString reasonForDisabling = "requested by user";
+ if (fUseTPCMultiplicityCorrection) {
+ if (isPP)
+ reasonForDisabling = "pp collisions";
+ else
+ reasonForDisabling = "MC w/o tune on data";
+ }
+
AliInfo(Form("Multiplicity correction %sdisabled (%s)!", fUseTPCMultiplicityCorrection ? "automatically " : "",
- fUseTPCMultiplicityCorrection ? "pp collisions" : "requested by user"));
+ reasonForDisabling.Data()));
fUseTPCMultiplicityCorrection = kFALSE;
fTPCResponse.ResetMultiplicityCorrectionFunctions();
}
- /*
- //TODO NOW start
- for (Int_t i = 0; i <= 4 + 1; i++) {
- printf("parMultCorr: %d, %e\n", i, fTPCResponse.GetMultiplicityCorrectionFunction()->GetParameter(i));
- }
- for (Int_t j = 0; j <= 2 + 1; j++) {
- printf("parMultCorrTanTheta: %d, %e\n", j, fTPCResponse.GetMultiplicityCorrectionFunctionTanTheta()->GetParameter(j));
- }
- for (Int_t j = 0; j <= 3 + 1; j++) {
- printf("parMultSigmaCorr: %d, %e\n", j, fTPCResponse.GetMultiplicitySigmaCorrectionFunction()->GetParameter(j));
+ if (fUseTPCMultiplicityCorrection) {
+ for (Int_t i = 0; i <= 4 + 1; i++) {
+ AliInfo(Form("parMultCorr: %d, %e", i, fTPCResponse.GetMultiplicityCorrectionFunction()->GetParameter(i)));
+ }
+ for (Int_t j = 0; j <= 2 + 1; j++) {
+ AliInfo(Form("parMultCorrTanTheta: %d, %e", j, fTPCResponse.GetMultiplicityCorrectionFunctionTanTheta()->GetParameter(j)));
+ }
+ for (Int_t j = 0; j <= 3 + 1; j++) {
+ AliInfo(Form("parMultSigmaCorr: %d, %e", j, fTPCResponse.GetMultiplicitySigmaCorrectionFunction()->GetParameter(j)));
+ }
}
- //TODO NOW end
- */
-
//
// Setup old resolution parametrisation
//
if (fHMPIDPIDParams) delete fHMPIDPIDParams;
fHMPIDPIDParams=NULL;
- TFile *oadbf = new TFile(Form("%s/COMMON/PID/data/HMPIDPIDParams.root",fOADBPath.Data()));
+ TFile *oadbf;
+ if(!fIsMC) oadbf = new TFile(Form("%s/COMMON/PID/data/HMPIDPIDParams.root",fOADBPath.Data()));
+ else oadbf = new TFile(Form("%s/COMMON/PID/MC/HMPIDPIDParams.root",fOADBPath.Data()));
if (oadbf && oadbf->IsOpen()) {
AliInfo(Form("Loading HMPID Params from %s/COMMON/PID/data/HMPIDPIDParams.root", fOADBPath.Data()));
AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("HMPoadb");
}
//______________________________________________________________________________
-Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel,Double_t centrality,AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const {
+Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack,Double_t efficiencyLevel,Double_t centrality,AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const {
+ // old function for compatibility
+ Int_t ntracklets=0;
+ return IdentifiedAsElectronTRD(vtrack,ntracklets,efficiencyLevel,centrality,PIDmethod);
+}
+
+//______________________________________________________________________________
+Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Int_t &ntracklets,Double_t efficiencyLevel,Double_t centrality,AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const {
//
// Check whether track is identified as electron under a given electron efficiency hypothesis
//
+ // ntracklets is the number of tracklets that has been used to calculate the PID signal
Double_t probs[AliPID::kSPECIES];
- ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs,PIDmethod);
- Int_t ntracklets = vtrack->GetTRDntrackletsPID();
+ ntracklets =CalculateTRDResponse(vtrack,probs,PIDmethod);
+
// Take mean of the TRD momenta in the given tracklets
Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes];
Int_t nmomenta = 0;
//
// Compute PID probabilities for TOF
//
-
+
+ fgTOFmismatchProb = 1E-8;
+
+ // centrality --> fCurrCentrality
+ // Beam type --> fBeamTypeNum
+ // N TOF cluster --> TOF header --> to get the TOF header we need to add a virtual method in AliVTrack extended to ESD and AOD tracks
+ // isMC --> fIsMC
+
+ Int_t nTOFcluster = 0;
+ if(track->GetTOFHeader() && track->GetTOFHeader()->GetTriggerMask()){ // N TOF clusters available
+ nTOFcluster = track->GetTOFHeader()->GetNumberOfTOFclusters();
+ if(fIsMC) nTOFcluster *= 1.5; // +50% in MC
+ }
+ else{
+ switch(fBeamTypeNum){
+ case kPP: // pp 7 TeV
+ nTOFcluster = 50;
+ break;
+ case kPPB: // pPb 5.05 ATeV
+ nTOFcluster = 50 + (100-fCurrCentrality)*50;
+ break;
+ case kPBPB: // PbPb 2.76 ATeV
+ nTOFcluster = 50 + (100-fCurrCentrality)*150;
+ break;
+ }
+ }
+
+ //fTOFResponse.GetMismatchProbability(track->GetTOFsignal(),track->Eta()) * 0.01; // for future implementation of mismatch (i.e. 1% mismatch that should be extended for PbPb, pPb)
+
// set flat distribution (no decision)
for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
const Double_t expTime = fTOFResponse.GetExpectedSignal(track,type);
const Double_t sig = fTOFResponse.GetExpectedSigma(track->P(),expTime,AliPID::ParticleMassZ(type));
- if (TMath::Abs(nsigmas) > (fRange+2)) {
- if(nsigmas < fTOFtail)
- p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
- else
- p[j] = TMath::Exp(-(fRange+2 - fTOFtail*0.5)*fTOFtail)/sig;
- } else{
- if(nsigmas < fTOFtail)
- p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
- else
- p[j] = TMath::Exp(-(nsigmas - fTOFtail*0.5)*fTOFtail)/sig;
- }
+
+ if(nsigmas < fTOFtail)
+ p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
+ else
+ p[j] = TMath::Exp(-(nsigmas - fTOFtail*0.5)*fTOFtail)/sig;
+
+ p[j] += fgTOFmismatchProb;
}
return kDetPidOk;
}
+
+Int_t AliPIDResponse::CalculateTRDResponse(const AliVTrack *track,Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const
+{
+ // new function for backward compatibility
+ // returns number of tracklets PID
+
+ UInt_t TRDslicesForPID[2];
+ SetTRDSlices(TRDslicesForPID,PIDmethod);
+
+ Float_t mom[6]={0.};
+ Double_t dedx[48]={0.}; // Allocate space for the maximum number of TRD slices
+ Int_t nslices = TRDslicesForPID[1] - TRDslicesForPID[0] + 1;
+ AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d", TRDslicesForPID[0], TRDslicesForPID[1], nslices));
+ for(UInt_t ilayer = 0; ilayer < 6; ilayer++){
+ mom[ilayer] = track->GetTRDmomentum(ilayer);
+ for(UInt_t islice = TRDslicesForPID[0]; islice <= TRDslicesForPID[1]; islice++){
+ dedx[ilayer*nslices+islice-TRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice);
+ }
+ }
+
+ return fTRDResponse.GetResponse(nslices, dedx, mom, p,PIDmethod);
+
+}
//______________________________________________________________________________
-AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTRDProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod/*=AliTRDPIDResponse::kLQ1D*/) const
+AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTRDProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const
{
//
// Compute PID probabilities for the TRD
const EDetPidStatus pidStatus=GetTRDPIDStatus(track);
if (pidStatus!=kDetPidOk) return pidStatus;
- UInt_t TRDslicesForPID[2];
- SetTRDSlices(TRDslicesForPID,PIDmethod);
-
- Float_t mom[6]={0.};
- Double_t dedx[48]={0.}; // Allocate space for the maximum number of TRD slices
- Int_t nslices = TRDslicesForPID[1] - TRDslicesForPID[0] + 1;
- AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d", TRDslicesForPID[0], TRDslicesForPID[1], nslices));
- for(UInt_t ilayer = 0; ilayer < 6; ilayer++){
- mom[ilayer] = track->GetTRDmomentum(ilayer);
- for(UInt_t islice = TRDslicesForPID[0]; islice <= TRDslicesForPID[1]; islice++){
- dedx[ilayer*nslices+islice-TRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice);
- }
- }
-
- fTRDResponse.GetResponse(nslices, dedx, mom, p,PIDmethod);
+ CalculateTRDResponse(track,p,PIDmethod);
+
return kDetPidOk;
}
// compute mismatch probability cross-checking at 5 sigmas with TPC
// currently just implemented as a 5 sigma compatibility cut
+ if(!track) return fgTOFmismatchProb;
+
// check pid status
const EDetPidStatus tofStatus=GetTOFPIDStatus(track);
if (tofStatus!=kDetPidOk) return 0.;
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff