#include <TArrayI.h>
#include <TArrayF.h>
#include <TLinearFitter.h>
+#include <TSystem.h>
+#include <TMD5.h>
#include <AliVEvent.h>
#include <AliVTrack.h>
ClassImp(AliPIDResponse);
+Float_t AliPIDResponse::fgTOFmismatchProb = 0.0;
+
AliPIDResponse::AliPIDResponse(Bool_t isMC/*=kFALSE*/) :
TNamed("PIDResponse","PIDResponse"),
fITSResponse(isMC),
fMCperiodTPC(),
fMCperiodUser(),
fCurrentFile(),
+fCurrentAliRootRev(-1),
fRecoPass(0),
fRecoPassUser(-1),
fRun(-1),
fArrPidResponseMaster(NULL),
fResolutionCorrection(NULL),
fOADBvoltageMaps(NULL),
-fUseTPCEtaCorrection(kFALSE),//TODO: In future, default kTRUE
+fUseTPCEtaCorrection(kFALSE),
+fUseTPCMultiplicityCorrection(kFALSE),
fTRDPIDResponseObject(NULL),
-fTOFtail(1.1),
+fTOFtail(0.9),
fTOFPIDParams(NULL),
fHMPIDPIDParams(NULL),
fEMCALPIDParams(NULL),
fCurrentEvent(NULL),
-fCurrCentrality(0.0)
+fCurrCentrality(0.0),
+fBeamTypeNum(kPP)
{
//
// default ctor
fMCperiodTPC(),
fMCperiodUser(other.fMCperiodUser),
fCurrentFile(),
+fCurrentAliRootRev(other.fCurrentAliRootRev),
fRecoPass(0),
fRecoPassUser(other.fRecoPassUser),
fRun(-1),
fResolutionCorrection(NULL),
fOADBvoltageMaps(NULL),
fUseTPCEtaCorrection(other.fUseTPCEtaCorrection),
+fUseTPCMultiplicityCorrection(other.fUseTPCMultiplicityCorrection),
fTRDPIDResponseObject(NULL),
-fTOFtail(1.1),
+fTOFtail(0.9),
fTOFPIDParams(NULL),
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;
fCurrentFile="";
+ fCurrentAliRootRev=other.fCurrentAliRootRev;
fRecoPass=0;
fRecoPassUser=other.fRecoPassUser;
fRun=-1;
fResolutionCorrection=NULL;
fOADBvoltageMaps=NULL;
fUseTPCEtaCorrection=other.fUseTPCEtaCorrection;
+ fUseTPCMultiplicityCorrection=other.fUseTPCMultiplicityCorrection;
fTRDPIDResponseObject=NULL;
fEMCALPIDParams=NULL;
- fTOFtail=1.1;
+ fTOFtail=0.9;
fTOFPIDParams=NULL;
fHMPIDPIDParams=NULL;
fCurrentEvent=other.fCurrentEvent;
+
}
return *this;
}
//get number of sigmas according the selected TPC gain configuration scenario
const AliVTrack *track=static_cast<const AliVTrack*>(vtrack);
-// return 0.;
- Float_t nSigma=fTPCResponse.GetNumberOfSigmas(track, type, dedxSource, fUseTPCEtaCorrection);
+ Float_t nSigma=fTPCResponse.GetNumberOfSigmas(track, type, dedxSource, fUseTPCEtaCorrection, fUseTPCMultiplicityCorrection);
return nSigma;
}
fTPCResponse.SetSigma(3.79301e-03*corrSigma, 2.21280e+04);
}
+ // Set up TPC multiplicity for PbPb
+ 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);
+ }
+ else
+ fTPCResponse.SetCurrentEventMultiplicity(0);
+
//TOF resolution
SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFPIDParams->GetStartTimeMethod());
else{
fCurrCentrality = -1;
}
+
+ // Set centrality percentile for EMCAL
+ fEMCALResponse.SetCentrality(fCurrCentrality);
+
+ // switch off some TOF channel according to OADB to match data TOF matching eff
+ if (fTuneMConData && ((fTuneMConDataMask & kDetTOF) == kDetTOF) && fTOFPIDParams->GetTOFmatchingLossMC() > 0.01){
+ Int_t ntrk = event->GetNumberOfTracks();
+ for(Int_t i=0;i < ntrk;i++){
+ AliVParticle *trk = event->GetTrack(i);
+ Int_t channel = GetTOFResponse().GetTOFchannel(trk);
+ Int_t swoffEachOfThem = Int_t(100./fTOFPIDParams->GetTOFmatchingLossMC() + 0.5);
+ if(!(channel%swoffEachOfThem)) ((AliVTrack *) trk)->ResetStatus(AliVTrack::kTOFout);
+ }
+ }
+
}
//______________________________________________________________________________
fBeamType="";
fBeamType="PP";
+ fBeamTypeNum=kPP;
+
+ Bool_t hasProdInfo=(fCurrentFile.BeginsWith("LHC"));
- TPRegexp reg(".*(LHC1[1-2][a-z]+[0-9]+[a-z_]*)/.*");
- TPRegexp reg12a17("LHC1[2-3][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-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"; }
fLHCperiod="LHC10H";
fMCperiodTPC="LHC10H8";
if (reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11A10";
+ // 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="";*/ }
- // for the moment use LHC12b parameters up to LHC12e
- 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="";*/ }
-// for the moment use 12g parametrisation for all full gain runs (LHC12f+)
- if (fRun >= 186636 && fRun < 194480) { fLHCperiod="LHC12G"; fBeamType="PPB"; fMCperiodTPC="LHC12G"; }
- if (fRun >= 194480) { fLHCperiod="LHC13B"; fBeamType="PPB"; fMCperiodTPC="LHC12G"; }
-
- //exception new pp MC productions from 2011
- if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) { fMCperiodTPC="LHC11B2"; fBeamType="PP"; }
+ if (fRun>=170719 && fRun<=177311) {
+ fLHCperiod="LHC12A";
+ fBeamType="PP";
+ fBeamTypeNum=kPP;
+ fMCperiodTPC="LHC10F6A";
+ if (fCurrentAliRootRev >= 62714)
+ fMCperiodTPC="LHC13B2_FIXn1";
+ }
+ // for the moment use LHC12b parameters up to LHC12d
+ if (fRun>=177312 /*&& fRun<=179356*/) {
+ fLHCperiod="LHC12B";
+ fBeamType="PP";
+ fBeamTypeNum=kPP;
+ fMCperiodTPC="LHC10F6A";
+ if (fCurrentAliRootRev >= 62714)
+ fMCperiodTPC="LHC13B2_FIXn1";
+ }
+// 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";fBeamTypeNum=kPPB; fMCperiodTPC="LHC12G"; }
+
+ // New parametrisation for 2013 pPb runs
+ if (fRun >= 194480) {
+ fLHCperiod="LHC13B";
+ fBeamType="PPB";
+ fBeamTypeNum=kPPB;
+ fMCperiodTPC="LHC12G";
+
+ if (fCurrentAliRootRev >= 61605)
+ fMCperiodTPC="LHC13B2_FIX";
+ if (fCurrentAliRootRev >= 62714)
+ fMCperiodTPC="LHC13B2_FIXn1";
+
+ // High luminosity pPb runs require different parametrisations
+ if (fRun >= 195875 && fRun <= 197411) {
+ fLHCperiod="LHC13F";
+ }
+ }
+
+ //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";fBeamTypeNum=kPP; }
// exception for 11f1
- if (fCurrentFile.Contains("LHC11f1/")) fMCperiodTPC="LHC11F1";
+ if (fCurrentFile.Contains("LHC11f1")) fMCperiodTPC="LHC11F1";
// exception for 12f1a, 12f1b and 12i3
- if (fCurrentFile.Contains("LHC12f1a/") || fCurrentFile.Contains("LHC12f1b/")
- || fCurrentFile.Contains("LHC12i3/")) fMCperiodTPC="LHC12F1";
+ if (fCurrentFile.Contains("LHC12f1") || fCurrentFile.Contains("LHC12i3")) fMCperiodTPC="LHC12F1";
// exception for 12c4
- if (fCurrentFile.Contains("LHC12c4/")) fMCperiodTPC="LHC12C4";
+ if (fCurrentFile.Contains("LHC12c4")) fMCperiodTPC="LHC12C4";
+ // exception for 12d and 13d pp periods
+ if (fBeamType=="PP" && fCurrentAliRootRev >= 61605) fMCperiodTPC="LHC13D1";
}
//______________________________________________________________________________
TString period = fLHCperiod.IsNull() ? "No period information" : fLHCperiod;
if (fIsMC) {
- if (!fTuneMConData) {
+ if (!(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))) {
period=fMCperiodTPC;
dataType="MC";
}
fRecoPass = 1;
- if (!fTuneMConData && fMCperiodTPC.IsNull()) {
+ if (!(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC)) && fMCperiodTPC.IsNull()) {
AliFatal("MC detected, but no MC period set -> Not changing eta maps!");
return;
}
Form("TPCetaMaps_%s_pass%d", dataType.Data(), recopass));
if (statusCont) {
AliError("Failed initializing TPC eta correction maps from OADB -> Disabled eta correction");
+ fUseTPCEtaCorrection = kFALSE;
}
else {
AliInfo(Form("Loading TPC eta correction map from %s/COMMON/PID/data/TPCetaMaps.root", fOADBPath.Data()));
TH2D* etaMap = 0x0;
- if (fIsMC && !fTuneMConData) {
+ if (fIsMC && !(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))) {
TString searchMap = Form("TPCetaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), recopass);
etaMap = dynamic_cast<TH2D *>(etaMapsCont.GetDefaultObject(searchMap.Data()));
if (!etaMap) {
if (!etaMap) {
AliError(Form("TPC eta correction map not found for run %d and also no default map found -> Disabled eta correction!!!", fRun));
+ fUseTPCEtaCorrection = kFALSE;
}
else {
TH2D* etaMapRefined = RefineHistoViaLinearInterpolation(etaMap, refineFactorMapX, refineFactorMapY);
if (!fTPCResponse.SetEtaCorrMap(etaMapRefined)) {
AliError(Form("Failed to set TPC eta correction map for run %d -> Disabled eta correction!!!", fRun));
fTPCResponse.SetEtaCorrMap(0x0);
+ fUseTPCEtaCorrection = kFALSE;
}
else {
- AliInfo(Form("Loaded TPC eta correction map (refine factors %.2f/%.2f) from %s/COMMON/PID/data/TPCetaMaps.root: %s",
- refineFactorMapX, refineFactorMapY, fOADBPath.Data(), fTPCResponse.GetEtaCorrMap()->GetTitle()));
+ AliInfo(Form("Loaded TPC eta correction map (refine factors %.2f/%.2f) from %s/COMMON/PID/data/TPCetaMaps.root: %s (MD5(map) = %s)",
+ refineFactorMapX, refineFactorMapY, fOADBPath.Data(), fTPCResponse.GetEtaCorrMap()->GetTitle(),
+ GetChecksum(fTPCResponse.GetEtaCorrMap()).Data()));
}
delete etaMapRefined;
}
else {
AliError(Form("Failed to set TPC eta correction map for run %d (map was loaded, but couldn't be refined) -> Disabled eta correction!!!", fRun));
+ fUseTPCEtaCorrection = kFALSE;
}
}
}
+ // If there was some problem loading the eta maps, it makes no sense to load the sigma maps (that require eta corrected data)
+ if (fUseTPCEtaCorrection == kFALSE) {
+ AliError("Failed to load TPC eta correction map required by sigma maps -> Using old parametrisation for sigma");
+ return;
+ }
+
// Load the sigma parametrisation (1/dEdx vs tanTheta_local (~eta))
AliOADBContainer etaSigmaMapsCont(Form("TPCetaSigmaMaps_%s_pass%d", dataType.Data(), recopass));
TObjArray* etaSigmaPars = 0x0;
- if (fIsMC && !fTuneMConData) {
+ if (fIsMC && !(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))) {
TString searchMap = Form("TPCetaSigmaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), recopass);
etaSigmaPars = dynamic_cast<TObjArray *>(etaSigmaMapsCont.GetDefaultObject(searchMap.Data()));
if (!etaSigmaPars) {
fTPCResponse.SetSigmaParams(0x0, 0);
}
else {
- AliInfo(Form("Loaded TPC sigma correction map (refine factors %.2f/%.2f) from %s/COMMON/PID/data/TPCetaMaps.root: %s",
- refineFactorSigmaMapX, refineFactorSigmaMapY, fOADBPath.Data(), fTPCResponse.GetSigmaPar1Map()->GetTitle()));
+ AliInfo(Form("Loaded TPC sigma correction map (refine factors %.2f/%.2f) from %s/COMMON/PID/data/TPCetaMaps.root: %s (MD5(map) = %s, sigmaPar0 = %f)",
+ refineFactorSigmaMapX, refineFactorSigmaMapY, fOADBPath.Data(), fTPCResponse.GetSigmaPar1Map()->GetTitle(),
+ GetChecksum(fTPCResponse.GetSigmaPar1Map()).Data(), sigmaPar0));
}
delete etaSigmaPar1MapRefined;
TString datatype="DATA";
//in case of mc fRecoPass is per default 1
if (fIsMC) {
- if(!fTuneMConData) datatype="MC";
+ if(!(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))) datatype="MC";
fRecoPass=1;
}
// period
TString period=fLHCperiod;
- if (fIsMC && !fTuneMConData) period=fMCperiodTPC;
+ if (fIsMC && !(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))) period=fMCperiodTPC;
Int_t recopass = fRecoPass;
- if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC) ) recopass = fRecoPassUser;
+ if(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC)) recopass = fRecoPassUser;
AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
Bool_t found=kFALSE;
(AliPID::EParticleType)ispec,
(AliTPCPIDResponse::ETPCgainScenario)igainScenario );
fTPCResponse.SetUseDatabase(kTRUE);
- AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,responseFunction->GetName()));
+ AliInfo(Form("Adding graph: %d %d - %s (MD5(spline) = %s)",ispec,igainScenario,responseFunction->GetName(),
+ GetChecksum((TSpline3*)responseFunction).Data()));
found=kTRUE;
break;
}
(AliPID::EParticleType)ispec,
(AliTPCPIDResponse::ETPCgainScenario)igainScenario );
fTPCResponse.SetUseDatabase(kTRUE);
- AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,responseFunctionPion->GetName()));
+ AliInfo(Form("Adding graph: %d %d - %s (MD5(spline) = %s)",ispec,igainScenario,responseFunctionPion->GetName(),
+ GetChecksum((TSpline3*)responseFunctionPion).Data()));
found=kTRUE;
}
else if (grAll) {
(AliPID::EParticleType)ispec,
(AliTPCPIDResponse::ETPCgainScenario)igainScenario );
fTPCResponse.SetUseDatabase(kTRUE);
- AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,grAll->GetName()));
+ AliInfo(Form("Adding graph: %d %d - %s (MD5(spline) = %s)",ispec,igainScenario,grAll->GetName(),
+ GetChecksum((TSpline3*)grAll).Data()));
found=kTRUE;
}
//else
(AliPID::EParticleType)ispec,
(AliTPCPIDResponse::ETPCgainScenario)igainScenario );
fTPCResponse.SetUseDatabase(kTRUE);
- AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,responseFunctionProton->GetName()));
+ AliInfo(Form("Adding graph: %d %d - %s (MD5(spline) = %s)",ispec,igainScenario,responseFunctionProton->GetName(),
+ GetChecksum((TSpline3*)responseFunctionProton).Data()));
found=kTRUE;
}
else if (grAll) {
(AliPID::EParticleType)ispec,
(AliTPCPIDResponse::ETPCgainScenario)igainScenario );
fTPCResponse.SetUseDatabase(kTRUE);
- AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,grAll->GetName()));
+ AliInfo(Form("Adding graph: %d %d - %s (MD5(spline) = %s)",ispec,igainScenario,grAll->GetName(),
+ GetChecksum((TSpline3*)grAll).Data()));
found=kTRUE;
}
//else
AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
}
+
//
- // Setup resolution parametrisation
+ // Setup multiplicity correction (only used for non-pp collisions)
+ //
+
+ 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 no correction for MC
+ if (period.Contains("LHC11A10")) {//LHC11A10A
+ AliInfo("Using multiplicity correction parameters for 11a10!");
+ fTPCResponse.SetParameterMultiplicityCorrection(0, 6.90133e-06);
+ fTPCResponse.SetParameterMultiplicityCorrection(1, -1.22123e-03);
+ fTPCResponse.SetParameterMultiplicityCorrection(2, 1.80220e-02);
+ fTPCResponse.SetParameterMultiplicityCorrection(3, 0.1);
+ fTPCResponse.SetParameterMultiplicityCorrection(4, 6.45306e-03);
+
+ fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(0, -2.85505e-07);
+ fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(1, -1.31911e-06);
+ fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(2, -0.5);
+
+ fTPCResponse.SetParameterMultiplicitySigmaCorrection(0, -4.29665e-05);
+ fTPCResponse.SetParameterMultiplicitySigmaCorrection(1, 1.37023e-02);
+ fTPCResponse.SetParameterMultiplicitySigmaCorrection(2, -6.36337e-01);
+ fTPCResponse.SetParameterMultiplicitySigmaCorrection(3, 1.13479e-02);
+ }
+ 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.SetParameterMultiplicityCorrection(2, -3.521e-02);
+ fTPCResponse.SetParameterMultiplicityCorrection(3, 2.469e-02);
+ fTPCResponse.SetParameterMultiplicityCorrection(4, 0);
+
+ fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(0, -5.32e-06);
+ fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(1, 1.177e-05);
+ fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(2, -0.5);
+
+ fTPCResponse.SetParameterMultiplicitySigmaCorrection(0, 0.);
+ fTPCResponse.SetParameterMultiplicitySigmaCorrection(1, 0.);
+ fTPCResponse.SetParameterMultiplicitySigmaCorrection(2, 0.);
+ fTPCResponse.SetParameterMultiplicitySigmaCorrection(3, 0.);
+
+ /* Not too bad, but far from perfect in the details
+ fTPCResponse.SetParameterMultiplicityCorrection(0, -6.27187e-06);
+ fTPCResponse.SetParameterMultiplicityCorrection(1, -4.60649e-04);
+ fTPCResponse.SetParameterMultiplicityCorrection(2, -4.26450e-02);
+ fTPCResponse.SetParameterMultiplicityCorrection(3, 2.40590e-02);
+ fTPCResponse.SetParameterMultiplicityCorrection(4, 0);
+
+ fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(0, -5.338e-06);
+ fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(1, 1.220e-05);
+ fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(2, -0.5);
+
+ fTPCResponse.SetParameterMultiplicitySigmaCorrection(0, 7.89237e-05);
+ fTPCResponse.SetParameterMultiplicitySigmaCorrection(1, -1.30662e-02);
+ fTPCResponse.SetParameterMultiplicitySigmaCorrection(2, 8.91548e-01);
+ fTPCResponse.SetParameterMultiplicitySigmaCorrection(3, 1.47931e-02);
+ */
+ }
+ /*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.SetParameterMultiplicityCorrection(2, 1.28786e-03);
+ fTPCResponse.SetParameterMultiplicityCorrection(3, 1.47677e-02);
+ fTPCResponse.SetParameterMultiplicityCorrection(4, 0);
+
+ fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(0, 7.23591e-08);
+ fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(1, 2.7469e-06);
+ fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(2, -0.5);
+
+ fTPCResponse.SetParameterMultiplicitySigmaCorrection(0, -1.22590e-05);
+ fTPCResponse.SetParameterMultiplicitySigmaCorrection(1, 6.88888e-03);
+ 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;
+ fTPCResponse.ResetMultiplicityCorrectionFunctions();
+ }
+ }
+ else {
+ // Just set parameters such that overall correction factor is 1, i.e. no correction.
+ // This is just a reasonable choice for the parameters for safety reasons. Disabling
+ // the multiplicity correction will anyhow skip the calculation of the corresponding
+ // correction factor inside THIS class. Nevertheless, experts can access the TPCPIDResponse
+ // 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 " : "",
+ reasonForDisabling.Data()));
+
+ fUseTPCMultiplicityCorrection = kFALSE;
+ fTPCResponse.ResetMultiplicityCorrectionFunctions();
+ }
+
+ 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)));
+ }
+ }
+
+ //
+ // Setup old resolution parametrisation
//
//default
fTPCResponse.SetSigma(3.79301e-03, 2.21280e+04);
- if (fRun>=122195){
+ if (fRun>=122195){ //LHC10d
fTPCResponse.SetSigma(2.30176e-02, 5.60422e+02);
}
-
- if (fRun>=186636){
-// if (fRun>=188356){
+
+ if (fRun>=170719){ // LHC12a
+ fTPCResponse.SetSigma(2.95714e-03, 1.01953e+05);
+ }
+
+ if (fRun>=177312){ // LHC12b
+ fTPCResponse.SetSigma(3.74633e-03, 7.11829e+04 );
+ }
+
+ if (fRun>=186346){ // LHC12e
fTPCResponse.SetSigma(8.62022e-04, 9.08156e+05);
}
if (fArrPidResponseMaster)
fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
- if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s",fResolutionCorrection->GetName()));
+ if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s (MD5(corr function) = %s)",
+ fResolutionCorrection->GetName(), GetChecksum(fResolutionCorrection).Data()));
//read in the voltage map
TVectorF* gsm = 0x0;
//______________________________________________________________________________
void AliPIDResponse::SetTRDSlices(UInt_t TRDslicesForPID[2],AliTRDPIDResponse::ETRDPIDMethod method) const{
- if(fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
+ if(fLHCperiod.Contains("LHC10D") || fLHCperiod.Contains("LHC10E")){
// backward compatibility for setting with 8 slices
TRDslicesForPID[0] = 0;
TRDslicesForPID[1] = 7;
AliInfo(Form(" StartTime method %d",fTOFPIDParams->GetStartTimeMethod()));
AliInfo(Form(" TOF res. mom. params: %5.2f %5.2f %5.2f %5.2f",
fTOFPIDParams->GetSigParams(0),fTOFPIDParams->GetSigParams(1),fTOFPIDParams->GetSigParams(2),fTOFPIDParams->GetSigParams(3)));
-
+ AliInfo(Form(" Fraction of tracks within gaussian behaviour: %6.4f",fTOFPIDParams->GetTOFtail()));
+ AliInfo(Form(" MC: Fraction of tracks (percentage) to cut to fit matching in data: %6.2f%%",fTOFPIDParams->GetTOFmatchingLossMC()));
+ AliInfo(Form(" MC: Fraction of random hits (percentage) to add to fit mismatch in data: %6.2f%%",fTOFPIDParams->GetTOFadditionalMismForMC()));
+ AliInfo(Form(" Start Time Offset %6.2f ps",fTOFPIDParams->GetTOFtimeOffset()));
+
for (Int_t i=0;i<4;i++) {
fTOFResponse.SetTrackParameter(i,fTOFPIDParams->GetSigParams(i));
}
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;
// Set TOF response function
// Input option for event_time used
//
-
+
Float_t t0spread = 0.; //vevent->GetEventTimeSpread();
if(t0spread < 10) t0spread = 80;
- // T0 from TOF algorithm
+ // T0-FILL and T0-TO offset (because of TOF misallignment
+ Float_t starttimeoffset = 0;
+ if(fTOFPIDParams && !(fIsMC)) starttimeoffset=fTOFPIDParams->GetTOFtimeOffset();
+ if(fTOFPIDParams){
+ fTOFtail = fTOFPIDParams->GetTOFtail();
+ GetTOFResponse().SetTOFtail(fTOFtail);
+ }
+ // T0 from TOF algorithm
Bool_t flagT0TOF=kFALSE;
Bool_t flagT0T0=kFALSE;
Float_t *startTime = new Float_t[fTOFResponse.GetNmomBins()];
startTime[i]=tofHeader->GetDefaultEventTimeVal();
startTimeRes[i]=tofHeader->GetDefaultEventTimeRes();
if(startTimeRes[i] < 1.e-5) startTimeRes[i] = t0spread;
+
+ if(startTimeRes[i] > t0spread - 10 && TMath::Abs(startTime[i]) < 0.001) startTime[i] = -starttimeoffset; // apply offset for T0-fill
}
TArrayI *ibin=(TArrayI*)tofHeader->GetNvalues();
startTime[icurrent]=t0Bin->GetAt(j);
startTimeRes[icurrent]=t0ResBin->GetAt(j);
if(startTimeRes[icurrent] < 1.e-5) startTimeRes[icurrent] = t0spread;
+ if(startTimeRes[icurrent] > t0spread - 10 && TMath::Abs(startTime[icurrent]) < 0.001) startTime[icurrent] = -starttimeoffset; // apply offset for T0-fill
}
}
if(option == kFILL_T0){ // T0-FILL is used
for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
- estimatedT0event[i]=0.0;
+ estimatedT0event[i]=0.0-starttimeoffset;
estimatedT0resolution[i]=t0spread;
}
fTOFResponse.SetT0event(estimatedT0event);
}
else{
for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
- estimatedT0event[i]=0.0;
+ estimatedT0event[i]=0.0-starttimeoffset;
estimatedT0resolution[i]=t0spread;
fTOFResponse.SetT0binMask(i,startTimeMask[i]);
}
Float_t t0A=-10000;
Float_t t0C=-10000;
if(flagT0T0){
- t0A= vevent->GetT0TOF()[1];
- t0C= vevent->GetT0TOF()[2];
+ t0A= vevent->GetT0TOF()[1] - starttimeoffset;
+ t0C= vevent->GetT0TOF()[2] - starttimeoffset;
// t0AC= vevent->GetT0TOF()[0];
- t0AC= t0A/resT0A/resT0A + t0C/resT0C/resT0C;
- resT0AC= TMath::Sqrt(1./resT0A/resT0A + 1./resT0C/resT0C);
- t0AC /= resT0AC*resT0AC;
+ t0AC= t0A/resT0A/resT0A + t0C/resT0C/resT0C;
+ resT0AC= TMath::Sqrt(1./resT0A/resT0A + 1./resT0C/resT0C);
+ t0AC /= resT0AC*resT0AC;
}
Float_t t0t0Best = 0;
estimatedT0resolution[i]=t0t0BestRes;
}
else{
- estimatedT0event[i]=0.0;
+ estimatedT0event[i]=0.0-starttimeoffset;
estimatedT0resolution[i]=t0spread;
}
}
Float_t t0A=-10000;
Float_t t0C=-10000;
if(flagT0T0){
- t0A= vevent->GetT0TOF()[1];
- t0C= vevent->GetT0TOF()[2];
+ t0A= vevent->GetT0TOF()[1] - starttimeoffset;
+ t0C= vevent->GetT0TOF()[2] - starttimeoffset;
// t0AC= vevent->GetT0TOF()[0];
- t0AC= t0A/resT0A/resT0A + t0C/resT0C/resT0C;
- resT0AC= TMath::Sqrt(1./resT0A/resT0A + 1./resT0C/resT0C);
- t0AC /= resT0AC*resT0AC;
+ t0AC= t0A/resT0A/resT0A + t0C/resT0C/resT0C;
+ resT0AC= TMath::Sqrt(1./resT0A/resT0A + 1./resT0C/resT0C);
+ t0AC /= resT0AC*resT0AC;
}
if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
}
else{
for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
- estimatedT0event[i]=0.0;
+ estimatedT0event[i]= 0.0 - starttimeoffset;
estimatedT0resolution[i]=t0spread;
fTOFResponse.SetT0binMask(i,0);
}
fTOFResponse.SetT0event(estimatedT0event);
fTOFResponse.SetT0resolution(estimatedT0resolution);
}
+
delete [] startTime;
delete [] startTimeRes;
delete [] startTimeMask;
// the following call is needed in order to fill the transient data member
// fTPCsignalTuned which is used in the TPCPIDResponse to judge
// if using tuned on data
- if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC) ) this->GetTPCsignalTunedOnData(track);
+ if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))
+ this->GetTPCsignalTunedOnData(track);
- return fTPCResponse.GetNumberOfSigmas(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
+ return fTPCResponse.GetNumberOfSigmas(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection, fUseTPCMultiplicityCorrection);
}
//______________________________________________________________________________
// the following call is needed in order to fill the transient data member
// fTPCsignalTuned which is used in the TPCPIDResponse to judge
// if using tuned on data
- if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC) )
+ if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))
this->GetTPCsignalTunedOnData(track);
- val=fTPCResponse.GetSignalDelta(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection, ratio);
+ val=fTPCResponse.GetSignalDelta(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection, fUseTPCMultiplicityCorrection, ratio);
return GetTPCPIDStatus(track);
}
//
AliVTrack *track=(AliVTrack*)vtrack;
val=GetSignalDeltaTOFold(track, type, ratio);
+
return GetTOFPIDStatus(track);
}
Double_t dedx=track->GetTPCsignal();
Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
- if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC) ) dedx = this->GetTPCsignalTunedOnData(track);
+ if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC)) dedx = this->GetTPCsignalTunedOnData(track);
Double_t bethe = 0.;
Double_t sigma = 0.;
for (Int_t j=0; j<nSpecies; j++) {
AliPID::EParticleType type=AliPID::EParticleType(j);
- bethe=fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
- sigma=fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
+ bethe=fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection, fUseTPCMultiplicityCorrection);
+ sigma=fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection, fUseTPCMultiplicityCorrection);
if (TMath::Abs(dedx-bethe) > fRange*sigma) {
p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
//
// 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 EDetPidStatus pidStatus=GetTOFPIDStatus(track);
if (pidStatus!=kDetPidOk) return pidStatus;
- const Double_t meanCorrFactor = 0.11/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1)
+ const Double_t meanCorrFactor = 0.07/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1)
for (Int_t j=0; j<nSpecies; j++) {
AliPID::EParticleType type=AliPID::EParticleType(j);
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.;
return kDetNoSignal;
}
+
+//______________________________________________________________________________
+TString AliPIDResponse::GetChecksum(const TObject* obj) const
+{
+ // Return the checksum for an object obj (tested to work properly at least for histograms and TSplines).
+
+ TString fileName = Form("tempChecksum.C"); // File name must be fixed for data type "TSpline3", since the file name will end up in the file content!
+
+ // For parallel processing, a unique file pathname is required. Uniqueness can be guaranteed by using a unique directory name
+ UInt_t index = 0;
+ TString uniquePathName = Form("tempChecksum_%u", index);
+
+ // To get a unique path name, increase the index until no directory
+ // of such a name exists.
+ // NOTE: gSystem->AccessPathName(...) returns kTRUE, if the access FAILED!
+ while (!gSystem->AccessPathName(uniquePathName.Data()))
+ uniquePathName = Form("tempChecksum_%u", ++index);
+
+ if (gSystem->mkdir(uniquePathName.Data()) < 0) {
+ AliError("Could not create temporary directory to store temp file for checksum determination!");
+ return "ERROR";
+ }
+
+ TString option = "";
+
+ // Save object as a macro, which will be deleted immediately after the checksum has been computed
+ // (does not work for desired data types if saved as *.root for some reason) - one only wants to compare the content, not
+ // the modification time etc. ...
+ if (dynamic_cast<const TH1*>(obj))
+ option = "colz"; // Histos need this option, since w/o this option, a counter is added to the filename
+
+
+ // SaveAs must be called with the fixed fileName only, since the first argument goes into the file content
+ // for some object types. Thus, change the directory, save the file and then go back
+ TString oldDir = gSystem->pwd();
+ gSystem->cd(uniquePathName.Data());
+ obj->SaveAs(fileName.Data(), option.Data());
+ gSystem->cd(oldDir.Data());
+
+ // Use the file to calculate the MD5 checksum
+ TMD5* md5 = TMD5::FileChecksum(Form("%s/%s", uniquePathName.Data(), fileName.Data()));
+ TString checksum = md5->AsString();
+
+ // Clean up
+ delete md5;
+ gSystem->Exec(Form("rm -rf %s", uniquePathName.Data()));
+
+ return checksum;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff