Float_t dedx = track->GetTPCsignalTunedOnData();
if(dedx > 0) return dedx;
- Double_t mom = t->GetTPCmomentum();
-
dedx = t->GetTPCsignal();
track->SetTPCsignalTunedOnData(dedx);
kGood = kFALSE;
if(kGood){
- Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
- Double_t sigma=fTPCResponse.GetExpectedSigma(mom,t->GetTPCsignalN(),type);
- dedx = gRandom->Gaus(bethe,sigma);
-
+ //TODO maybe introduce different dEdxSources?
+ Double_t bethe = fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, this->UseTPCEtaCorrection());
+ Double_t sigma = fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, this->UseTPCEtaCorrection());
+ dedx = gRandom->Gaus(bethe,sigma);
+
if(iS == AliPID::ParticleCode(AliPID::kHe3) || iS == AliPID::ParticleCode(AliPID::kAlpha)) dedx *= 5;
}
if(dedx > 0) return dedx;
- Double_t mom = t->GetTPCmomentum();
dedx = t->GetTPCsignal();
track->SetTPCsignalTunedOnData(dedx);
if(dedx < 20) return dedx;
kGood = kFALSE;
if(kGood){
- Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
- Double_t sigma=fTPCResponse.GetExpectedSigma(mom,t->GetTPCsignalN(),type);
+ //TODO maybe introduce different dEdxSources?
+ Double_t bethe = fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, this->UseTPCEtaCorrection());
+ Double_t sigma = fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, this->UseTPCEtaCorrection());
dedx = gRandom->Gaus(bethe,sigma);
if(iS == AliPID::ParticleCode(AliPID::kHe3) || iS == AliPID::ParticleCode(AliPID::kAlpha)) dedx *= 5;
}
if (fRun >= 186636 ) { fLHCperiod="LHC12G"; fBeamType="PPB"; /*fMCperiodTPC="";*/ }
//exception new pp MC productions from 2011
- if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2";
+ if ( (fBeamType=="PP" || fBeamType=="PPB") && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2";
// exception for 11f1
if (fCurrentFile.Contains("LHC11f1/")) fMCperiodTPC="LHC11F1";
}
Double_t centerX = hRefined->GetXaxis()->GetBinCenter(binX);
Double_t centerY = hRefined->GetYaxis()->GetBinCenter(binY);
- /*TODO NOW NOW
+ /*
linExtrapolation->ClearPoints();
// For interpolation: Just take the corresponding bin from the old histo.
Double_t interpolatedValue = linExtrapolation->GetParameter(0) + linExtrapolation->GetParameter(1) * centerX
+ linExtrapolation->GetParameter(2) * centerY;
*/
- Double_t interpolatedValue = h->Interpolate(centerX, centerY) ;//TODO NOW NOW NOW
+ Double_t interpolatedValue = h->Interpolate(centerX, centerY) ;
hRefined->SetBinContent(binX, binY, interpolatedValue);
}
}
// Load the TPC eta correction maps from the OADB
//
+ if (fUseTPCEtaCorrection == kFALSE) {
+ // Disable eta correction via setting no maps
+ if (!fTPCResponse.SetEtaCorrMap(0x0))
+ AliInfo("Request to disable TPC eta correction -> Eta correction has been disabled");
+ else
+ AliError("Request to disable TPC eta correction -> Some error occured when unloading the correction maps");
+
+ if (!fTPCResponse.SetSigmaParams(0x0, 0))
+ AliInfo("Request to disable TPC eta correction -> Using old parametrisation for sigma");
+ else
+ AliError("Request to disable TPC eta correction -> Some error occured when unloading the sigma maps");
+
+ return;
+ }
+
TString dataType = "DATA";
TString period = fLHCperiod.IsNull() ? "No period information" : fLHCperiod;
if (fIsMC) {
- dataType = "MC";
+ if (!fTuneMConData) {
+ period=fMCperiodTPC;
+ dataType="MC";
+ }
fRecoPass = 1;
- if (fMCperiodTPC.IsNull()) {
+ if (!fTuneMConData && fMCperiodTPC.IsNull()) {
AliFatal("MC detected, but no MC period set -> Not changing eta maps!");
return;
}
-
- period = fMCperiodTPC;
}
+
+ Int_t recopass = fRecoPass;
+ if (fTuneMConData)
+ recopass = fRecoPassUser;
- TString defaultObj = Form("Default_%s_pass%d", dataType.Data(), fRecoPass);
+ TString defaultObj = Form("Default_%s_pass%d", dataType.Data(), recopass);
- AliInfo(Form("Current period and reco pass: %s.pass%d", period.Data(), fRecoPass));
- if (fUseTPCEtaCorrection == kFALSE) {
- // Disable eta correction via setting no maps
- if (!fTPCResponse.SetEtaCorrMap(0x0))
- AliInfo("Request to disable TPC eta correction -> Eta correction has been disabled");
- else
- AliError("Request to disable TPC eta correction -> Some error occured when unloading the correction maps");
-
- if (!fTPCResponse.SetSigmaParams(0x0, 0))
- AliInfo("Request to disable TPC eta correction -> Using old parametrisation for sigma");
- else
- AliError("Request to disable TPC eta correction -> Some error occured when unloading the sigma maps");
-
- return;
- }
+ AliInfo(Form("Current period and reco pass: %s.pass%d", period.Data(), recopass));
// Invalidate old maps
fTPCResponse.SetEtaCorrMap(0x0);
fTPCResponse.SetSigmaParams(0x0, 0);
// Load the eta correction maps
- AliOADBContainer etaMapsCont(Form("TPCetaMaps_%s_pass%d", dataType.Data(), fRecoPass));
+ AliOADBContainer etaMapsCont(Form("TPCetaMaps_%s_pass%d", dataType.Data(), recopass));
Int_t statusCont = etaMapsCont.InitFromFile(Form("%s/COMMON/PID/data/TPCetaMaps.root", fOADBPath.Data()),
- Form("TPCetaMaps_%s_pass%d", dataType.Data(), fRecoPass));
+ Form("TPCetaMaps_%s_pass%d", dataType.Data(), recopass));
if (statusCont) {
AliError("Failed initializing TPC eta correction maps from OADB -> Disabled eta correction");
}
TH2D* etaMap = 0x0;
- if (fIsMC) {
- TString searchMap = Form("TPCetaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), fRecoPass);
+ if (fIsMC && !fTuneMConData) {
+ TString searchMap = Form("TPCetaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), recopass);
etaMap = dynamic_cast<TH2D *>(etaMapsCont.GetDefaultObject(searchMap.Data()));
if (!etaMap) {
// Try default object
}
// Load the sigma parametrisation (1/dEdx vs tanTheta_local (~eta))
- AliOADBContainer etaSigmaMapsCont(Form("TPCetaSigmaMaps_%s_pass%d", dataType.Data(), fRecoPass));
+ AliOADBContainer etaSigmaMapsCont(Form("TPCetaSigmaMaps_%s_pass%d", dataType.Data(), recopass));
statusCont = etaSigmaMapsCont.InitFromFile(Form("%s/COMMON/PID/data/TPCetaMaps.root", fOADBPath.Data()),
- Form("TPCetaSigmaMaps_%s_pass%d", dataType.Data(), fRecoPass));
+ Form("TPCetaSigmaMaps_%s_pass%d", dataType.Data(), recopass));
if (statusCont) {
AliError("Failed initializing TPC eta sigma maps from OADB -> Using old sigma parametrisation");
}
TObjArray* etaSigmaPars = 0x0;
- if (fIsMC) {
- TString searchMap = Form("TPCetaSigmaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), fRecoPass);
+ if (fIsMC && !fTuneMConData) {
+ TString searchMap = Form("TPCetaSigmaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), recopass);
etaSigmaPars = dynamic_cast<TObjArray *>(etaSigmaMapsCont.GetDefaultObject(searchMap.Data()));
if (!etaSigmaPars) {
// Try default object
TString period=fLHCperiod;
if (fIsMC && !fTuneMConData) period=fMCperiodTPC;
- AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
+ Int_t recopass = fRecoPass;
+ if(fTuneMConData) recopass = fRecoPassUser;
+
+ AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
Bool_t found=kFALSE;
//
//set the new PID splines
//
if (fArrPidResponseMaster){
- Int_t recopass = fRecoPass;
- if(fTuneMConData) recopass = fRecoPassUser;
//for MC don't use period information
//if (fIsMC) period="[A-Z0-9]*";
//for MC use MC period information
//pattern for the default entry (valid for all particles)
TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN(_*)([A-Z1-9]*)",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
- //find particle id ang gain scenario
+ //find particle id and gain scenario
for (Int_t igainScenario=0; igainScenario<AliTPCPIDResponse::fgkNumberOfGainScenarios; igainScenario++)
{
TObject *grAll=NULL;
else AliInfo("no fArrPidResponseMaster");
if (!found){
- AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
+ AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
}
//
}
if (fArrPidResponseMaster)
- fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
+ 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()));
Double_t nSigma = -999.;
- //TODO: TPCsignalTunedOnData not taken into account for the new TPCPIDresponse functions, so take the old functions
- if (fTuneMConData) {
- Double_t mom = track->GetTPCmomentum();
- Double_t sig = this->GetTPCsignalTunedOnData(track);
- UInt_t sigN = track->GetTPCsignalN();
+ if (fTuneMConData)
+ this->GetTPCsignalTunedOnData(track);
- if (sigN>0)
- nSigma = fTPCResponse.GetNumberOfSigmas(mom, sig, sigN, type);
- }
- else {
- nSigma = fTPCResponse.GetNumberOfSigmas(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
- }
+ nSigma = fTPCResponse.GetNumberOfSigmas(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
return nSigma;
}
for (Int_t j=0; j<AliPID::kSPECIESC; j++) {
AliPID::EParticleType type=AliPID::EParticleType(j);
- //TODO: TPCsignalTunedOnData not taken into account for the new TPCPIDresponse functions, so take the old functions
- if (fTuneMConData) {
- Double_t mom = track->GetTPCmomentum();
- bethe=fTPCResponse.GetExpectedSignal(mom,type);
- sigma=fTPCResponse.GetExpectedSigma(mom,track->GetTPCsignalN(),type);
- }
- else {
- bethe=fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
- sigma=fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
- }
+ bethe=fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
+ sigma=fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
+
if (TMath::Abs(dedx-bethe) > fRange*sigma) {
p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
} else {
for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];}
delete fhEtaCorr;
+ fhEtaCorr=0x0;
if (that.fhEtaCorr){
fhEtaCorr = new TH2D(*(that.fhEtaCorr));
fhEtaCorr->SetDirectory(0);
}
delete fhEtaSigmaPar1;
+ fhEtaSigmaPar1=0x0;
if (that.fhEtaSigmaPar1){
fhEtaSigmaPar1 = new TH2D(*(that.fhEtaSigmaPar1));
fhEtaSigmaPar1->SetDirectory(0);
Bool_t AliTPCPIDResponse::ResponseFunctiondEdxN( const AliVTrack* track,
AliPID::EParticleType species,
ETPCdEdxSource dedxSource,
- Double_t& dEdx, Int_t& nPoints, ETPCgainScenario& gainScenario, TSpline3** responseFunction) const
+ Double_t& dEdx,
+ Int_t& nPoints,
+ ETPCgainScenario& gainScenario,
+ TSpline3** responseFunction) const
{
// Calculates the right parameters for PID
// dEdx parametrization for the proper gain scenario, dEdx
if (dedxSource == kdEdxDefault) {
// Fast handling for default case. In addition: Keep it simple (don't call additional functions) to
// avoid possible bugs
- dEdx = track->GetTPCsignal();
+
+ // GetTPCsignalTunedOnData will be non-positive, if it has not been set (i.e. in case of MC NOT tuned to data).
+ // If this is the case, just take the normal signal
+ dEdx = track->GetTPCsignalTunedOnData();
+ if (dEdx <= 0) {
+ dEdx = track->GetTPCsignal();
+ }
+
nPoints = track->GetTPCsignalN();
gainScenario = kDefault;
return kTRUE;
}
+ //TODO Proper handle of tuneMConData for other dEdx sources
Double32_t signal[4]; //0: IROC, 1: OROC medium, 2:OROC long, 3: OROC all (def. truncation used)
Char_t ncl[3]; //same
if (dEdxInfo) dEdxInfo->GetTPCSignalRegionInfo(signal,ncl,nrows);
//check if we cross a bad OROC in which case we reject
- EChamberStatus trackStatus = TrackStatus(track,2);
- if (trackStatus==kChamberOff || trackStatus==kChamberLowGain)
+ EChamberStatus trackOROCStatus = TrackStatus(track,2);
+ if (trackOROCStatus==kChamberOff || trackOROCStatus==kChamberLowGain)
{
return kFALSE;
}
{
case kdEdxOROC:
{
+ if (trackOROCStatus==kChamberInvalid) return kFALSE; //never reached OROC
dEdx = signal[3];
nPoints = ncl[2]+ncl[1];
gainScenario = kOROChigh;
//_________________________________________________________________________
-Bool_t AliTPCPIDResponse::sectorNumbersInOut(const AliVTrack* track,
- Double_t innerRadius,
- Double_t outerRadius,
- Float_t& inphi,
+Bool_t AliTPCPIDResponse::sectorNumbersInOut(Double_t* trackPositionInner,
+ Double_t* trackPositionOuter,
+ Float_t& inphi,
Float_t& outphi,
- Int_t& in,
- Int_t& out ) const
+ Int_t& in,
+ Int_t& out ) const
{
//calculate the sector numbers (equivalent to IROC chamber numbers) a track crosses
//for OROC chamber numbers add 36
//returned angles are between (0,2pi)
-
- Double_t trackPositionInner[3];
- Double_t trackPositionOuter[3];
-
- Bool_t trackAtInner=kTRUE;
- Bool_t trackAtOuter=kTRUE;
- const AliExternalTrackParam* ip = track->GetInnerParam();
-
- //if there is no inner param this could mean we're using the AOD track,
- //we still can extrapolate from the vertex - so use those params.
- if (ip) track=ip;
-
- if (!track->GetXYZAt(innerRadius, fMagField, trackPositionInner))
- trackAtInner=kFALSE;
-
- if (!track->GetXYZAt(outerRadius, fMagField, trackPositionOuter))
- trackAtOuter=kFALSE;
-
- if (!trackAtInner)
- {
- //if we dont even enter inner radius we do nothing
- inphi=0.0;
- outphi=0.0;
- in=0;
- out=0;
- return kFALSE;
- }
-
- if (!trackAtOuter)
- {
- //if we don't reach the outer radius check that the apex is indeed within the outer radius and use apex position
- Bool_t haveApex = TrackApex(track, fMagField, trackPositionOuter);
- Float_t apexRadius = TMath::Sqrt(trackPositionOuter[0]*trackPositionOuter[0]+trackPositionOuter[1]*trackPositionOuter[1]);
- if ( haveApex && apexRadius<=outerRadius && apexRadius>innerRadius)
- {
- //printf("pt: %.2f, apexRadius: %.2f(%s), x: %.2f, y: %.2f\n",track->Pt(),apexRadius,(haveApex)?"OK":"BAD",trackPositionOuter[0],trackPositionOuter[1]);
- }
- else
- {
- inphi=0.0;
- outphi=0.0;
- in=0;
- out=0;
- return kFALSE;
- }
- }
inphi = TMath::ATan2(trackPositionInner[1],trackPositionInner[0]);
outphi = TMath::ATan2(trackPositionOuter[1], trackPositionOuter[0]);
}
return kTRUE;
}
-
+
+
//_____________________________________________________________________________
Int_t AliTPCPIDResponse::sectorNumber(Double_t phi) const
{
Float_t outphi=0.;
Float_t innerRadius = (layer==1)?83.0:133.7;
Float_t outerRadius = (layer==1)?133.5:247.7;
- if (!sectorNumbersInOut(track,
- innerRadius,
- outerRadius,
+
+ /////////////////////////////////////////////////////////////////////////////
+ //find out where track enters and leaves the layer.
+ //
+ Double_t trackPositionInner[3];
+ Double_t trackPositionOuter[3];
+
+ //if there is no inner param this could mean we're using the AOD track,
+ //we still can extrapolate from the vertex - so use those params.
+ const AliExternalTrackParam* ip = track->GetInnerParam();
+ if (ip) track=ip;
+
+ Bool_t trackAtInner = track->GetXYZAt(innerRadius, fMagField, trackPositionInner);
+ Bool_t trackAtOuter = track->GetXYZAt(outerRadius, fMagField, trackPositionOuter);
+
+ if (!trackAtInner)
+ {
+ //if we dont even enter inner radius we do nothing and return invalid
+ inphi=0.0;
+ outphi=0.0;
+ in=0;
+ out=0;
+ return kChamberInvalid;
+ }
+
+ if (!trackAtOuter)
+ {
+ //if we don't reach the outer radius check that the apex is indeed within the outer radius and use apex position
+ Bool_t haveApex = TrackApex(track, fMagField, trackPositionOuter);
+ Float_t apexRadius = TMath::Sqrt(trackPositionOuter[0]*trackPositionOuter[0]+trackPositionOuter[1]*trackPositionOuter[1]);
+ if ( haveApex && apexRadius<=outerRadius && apexRadius>innerRadius)
+ {
+ //printf("pt: %.2f, apexRadius: %.2f(%s), x: %.2f, y: %.2f\n",track->Pt(),apexRadius,(haveApex)?"OK":"BAD",trackPositionOuter[0],trackPositionOuter[1]);
+ }
+ else
+ {
+ inphi=0.0;
+ outphi=0.0;
+ in=0;
+ out=0;
+ return kChamberInvalid;
+ }
+ }
+
+
+ if (!sectorNumbersInOut(trackPositionInner,
+ trackPositionOuter,
inphi,
outphi,
in,
- out)) return kChamberOff;
+ out)) return kChamberInvalid;
+ /////////////////////////////////////////////////////////////////////////////
+ //now we have the location of the track we can check
+ //if it is in a good/bad chamber
+ //
Bool_t sideA = kTRUE;
if (((in/18)%2==1) && ((out/18)%2==1)) sideA=kFALSE;
Float_t trackLengthInBad = 0.;
Float_t trackLengthInLowGain = 0.;
Float_t trackLengthTotal = TMath::Abs(outphi-inphi);
+ Float_t lengthFractionInBadSectors = 0.;
const Float_t sectorWidth = TMath::TwoPi()/18.;
else {deltaPhi=sectorWidth;}
Float_t v = fVoltageMap[(layer==1)?(j):(j+36)];
- if (v<=fBadIROCthreshhold) trackLengthInBad+=deltaPhi;
- if (v<=fLowGainIROCthreshold && v>fBadIROCthreshhold) trackLengthInLowGain+=deltaPhi;
+ if (v<=fBadIROCthreshhold)
+ {
+ trackLengthInBad+=deltaPhi;
+ lengthFractionInBadSectors=1.;
+ }
+ if (v<=fLowGainIROCthreshold && v>fBadIROCthreshhold)
+ {
+ trackLengthInLowGain+=deltaPhi;
+ lengthFractionInBadSectors=1.;
+ }
}
//for now low gain and bad (off) chambers are treated equally
- Float_t lengthFractionInBadSectors = (trackLengthInLowGain+trackLengthInBad)/trackLengthTotal;
+ if (trackLengthTotal>0)
+ lengthFractionInBadSectors = (trackLengthInLowGain+trackLengthInBad)/trackLengthTotal;
//printf("### side: %s, pt: %.2f, pz: %.2f, in: %i, out: %i, phiIN: %.2f, phiOUT: %.2f, rIN: %.2f, rOUT: %.2f\n",(sideA)?"A":"C",track->Pt(),track->Pz(),in,out,inphi,outphi,innerRadius,outerRadius);
position[1]=b[1]+b[1]*TMath::Abs(r)/norm;
position[2]=0.;
return kTRUE;
-}
\ No newline at end of file
+}
AliPID::EParticleType species,
ETPCdEdxSource dedxSource,
Double_t& dEdx, Int_t& nPoints, ETPCgainScenario& gainScenario, TSpline3** responseFunction) const;
- Bool_t sectorNumbersInOut(const AliVTrack* track,
- Double_t innerRadius, Double_t outerRadius,
+ Bool_t sectorNumbersInOut(Double_t* trackPositionInner,
+ Double_t* trackPositionOuter,
Float_t& phiIn, Float_t& phiOut,
Int_t& in, Int_t& out ) const;
AliTPCPIDResponse::EChamberStatus TrackStatus(const AliVTrack* track, Int_t layer) const;
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff