#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),
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);
// 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-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="";*/ }
- // 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>=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";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"; /*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>=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
- if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) { fMCperiodTPC="LHC11B2"; fBeamType="PP"; }
+ //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";
// exception for 12f1a, 12f1b and 12i3
if (fCurrentFile.Contains("LHC12f1") || fCurrentFile.Contains("LHC12i3")) fMCperiodTPC="LHC12F1";
// exception for 12c4
if (fCurrentFile.Contains("LHC12c4")) fMCperiodTPC="LHC12C4";
+ // exception for 12d and 13d pp periods
+ if (fBeamType=="PP" && fCurrentAliRootRev >= 61605) fMCperiodTPC="LHC13D1";
}
//______________________________________________________________________________
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;
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;
(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
//
- // 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") || period.Contains("LHC13E") ||
- period.Contains("LHC13F")) {// 2013 pPb data taking
- 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
//
//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;
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;
//
// 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.;
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;
+}