// 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 ) { fLHCperiod="LHC12G"; fBeamType="PPB"; /*fMCperiodTPC="";*/ }
+ if (fRun >= 186636 ) { fLHCperiod="LHC12G"; fBeamType="PPB"; fMCperiodTPC="LHC12G"; }
//exception new pp MC productions from 2011
- if ( (fBeamType=="PP" || fBeamType=="PPB") && reg.MatchB(fCurrentFile)) { fMCperiodTPC="LHC11B2"; fBeamType="PP"; }
+ if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) { fMCperiodTPC="LHC11B2"; fBeamType="PP"; }
// exception for 11f1
if (fCurrentFile.Contains("LHC11f1/")) fMCperiodTPC="LHC11F1";
+ // exception for 12f1a and 12f1b
+ if (fCurrentFile.Contains("LHC12f1a/") || fCurrentFile.Contains("LHC12f1b/")) fMCperiodTPC="LHC12F1";
}
//______________________________________________________________________________
Double_t upperMapBoundY = h->GetYaxis()->GetBinUpEdge(h->GetYaxis()->GetNbins());
Double_t lowerMapBoundY = h->GetYaxis()->GetBinLowEdge(1);
- Int_t nBinsX = 20;
+ Int_t nBinsX = 30;
// Binning was find to yield good results, if 40 bins are chosen for the range 0.0016 to 0.02. For the new variable range,
// scale the number of bins correspondingly
- Int_t nBinsY = TMath::Nint((upperMapBoundY - lowerMapBoundY) / (0.02 - lowerMapBoundY) * 40);
+ Int_t nBinsY = TMath::Nint((upperMapBoundY - lowerMapBoundY) / (0.02 - 0.0016) * 40);
Int_t nBinsXrefined = nBinsX * refineFactorX;
Int_t nBinsYrefined = nBinsY * refineFactorY;
Double_t centerX = hRefined->GetXaxis()->GetBinCenter(binX);
Double_t centerY = hRefined->GetYaxis()->GetBinCenter(binY);
- /*
+ /*OLD
linExtrapolation->ClearPoints();
// For interpolation: Just take the corresponding bin from the old histo.
}
}
+
+ // Problem: Interpolation does not work before/beyond center of first/last bin (as the name suggests).
+ // Therefore, for each row in dEdx: Take last bin from old map and interpolate values from center and edge.
+ // Assume line through these points and extropolate to last bin of refined map
+ const Double_t firstOldXbinUpEdge = h->GetXaxis()->GetBinUpEdge(1);
+ const Double_t firstOldXbinCenter = h->GetXaxis()->GetBinCenter(1);
+
+ const Double_t oldXbinHalfWidth = firstOldXbinUpEdge - firstOldXbinCenter;
+
+ const Double_t lastOldXbinLowEdge = h->GetXaxis()->GetBinLowEdge(h->GetNbinsX());
+ const Double_t lastOldXbinCenter = h->GetXaxis()->GetBinCenter(h->GetNbinsX());
+
+ for (Int_t binY = 1; binY <= nBinsYrefined; binY++) {
+ Double_t centerY = hRefined->GetYaxis()->GetBinCenter(binY);
+
+ const Double_t interpolatedCenterFirstXbin = h->Interpolate(firstOldXbinCenter, centerY);
+ const Double_t interpolatedUpEdgeFirstXbin = h->Interpolate(firstOldXbinUpEdge, centerY);
+
+ const Double_t extrapolationSlopeFirstXbin = (interpolatedUpEdgeFirstXbin - interpolatedCenterFirstXbin) / oldXbinHalfWidth;
+ const Double_t extrapolationOffsetFirstXbin = interpolatedCenterFirstXbin;
+
+
+ const Double_t interpolatedCenterLastXbin = h->Interpolate(lastOldXbinCenter, centerY);
+ const Double_t interpolatedLowEdgeLastXbin = h->Interpolate(lastOldXbinLowEdge, centerY);
+
+ const Double_t extrapolationSlopeLastXbin = (interpolatedCenterLastXbin - interpolatedLowEdgeLastXbin) / oldXbinHalfWidth;
+ const Double_t extrapolationOffsetLastXbin = interpolatedCenterLastXbin;
+
+ for (Int_t binX = 1; binX <= nBinsXrefined; binX++) {
+ Double_t centerX = hRefined->GetXaxis()->GetBinCenter(binX);
+
+ if (centerX < firstOldXbinCenter) {
+ Double_t extrapolatedValue = extrapolationOffsetFirstXbin + (centerX - firstOldXbinCenter) * extrapolationSlopeFirstXbin;
+ hRefined->SetBinContent(binX, binY, extrapolatedValue);
+ }
+ else if (centerX <= lastOldXbinCenter) {
+ continue;
+ }
+ else {
+ Double_t extrapolatedValue = extrapolationOffsetLastXbin + (centerX - lastOldXbinCenter) * extrapolationSlopeLastXbin;
+ hRefined->SetBinContent(binX, binY, extrapolatedValue);
+ }
+ }
+ }
+
delete linExtrapolation;
return hRefined;
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");
+ 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
+ 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;
//
for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=0.07;fResN2[i]=0.0;}
}
-
+/*TODO remove?
//_________________________________________________________________________
AliTPCPIDResponse::AliTPCPIDResponse(const Double_t *param):
TNamed(),
//
for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=param[1];fResN2[i]=param[2];}
}
-
+*/
//_________________________________________________________________________
AliTPCPIDResponse::~AliTPCPIDResponse()
for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];}
// Copy eta maps
- if (that.fhEtaCorr){
+ if (that.fhEtaCorr) {
fhEtaCorr = new TH2D(*(that.fhEtaCorr));
fhEtaCorr->SetDirectory(0);
}
- if (that.fhEtaSigmaPar1){
+ if (that.fhEtaSigmaPar1) {
fhEtaSigmaPar1 = new TH2D(*(that.fhEtaSigmaPar1));
fhEtaSigmaPar1->SetDirectory(0);
}
delete fhEtaCorr;
fhEtaCorr=0x0;
- if (that.fhEtaCorr){
+ if (that.fhEtaCorr) {
fhEtaCorr = new TH2D(*(that.fhEtaCorr));
fhEtaCorr->SetDirectory(0);
}
delete fhEtaSigmaPar1;
fhEtaSigmaPar1=0x0;
- if (that.fhEtaSigmaPar1){
+ 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,
+ Double_t& dEdx,
+ Int_t& nPoints,
+ ETPCgainScenario& gainScenario,
TSpline3** responseFunction) const
{
// Calculates the right parameters for PID
Float_t& inphi,
Float_t& outphi,
Int_t& in,
- Int_t& out ) const
+ Int_t& out ) const
{
//calculate the sector numbers (equivalent to IROC chamber numbers) a track crosses
//for OROC chamber numbers add 36
}
return kTRUE;
}
-
-
+
+
//_____________________________________________________________________________
Int_t AliTPCPIDResponse::sectorNumber(Double_t phi) const
{
return TMath::Floor(phi/width);
}
+
//_____________________________________________________________________________
void AliTPCPIDResponse::Print(Option_t* /*option*/) const
{
fResponseFunctions.Print();
}
+
//_____________________________________________________________________________
AliTPCPIDResponse::EChamberStatus AliTPCPIDResponse::TrackStatus(const AliVTrack* track, Int_t layer) const
{
/////////////////////////////////////////////////////////////////////////////
//find out where track enters and leaves the layer.
//
- Double_t trackPositionInner[3];
- Double_t trackPositionOuter[3];
-
+ 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 (!sectorNumbersInOut(trackPositionInner,
- trackPositionOuter,
- inphi,
- outphi,
- in,
+ if (!sectorNumbersInOut(trackPositionInner,
+ trackPositionOuter,
+ inphi,
+ outphi,
+ in,
out)) return kChamberInvalid;
/////////////////////////////////////////////////////////////////////////////
- //now we have the location of the track we can check
+ //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;
in=in%18;
Float_t lengthFractionInBadSectors = 0.;
const Float_t sectorWidth = TMath::TwoPi()/18.;
-
+
for (Int_t i=in; i<=out; i++)
{
int j=i;
if (i<0) j+=18; //correct for the negative values
if (!sideA) j+=18; //move to the correct side
-
+
Float_t deltaPhi = 0.;
Float_t phiEdge=sectorWidth*i;
if (inphi>phiEdge) {deltaPhi=phiEdge+sectorWidth-inphi;}
else if ((outphi>=phiEdge) && (outphi<(phiEdge+sectorWidth))) {deltaPhi=outphi-phiEdge;}
else {deltaPhi=sectorWidth;}
-
+
Float_t v = fVoltageMap[(layer==1)?(j):(j+36)];
- if (v<=fBadIROCthreshhold)
- {
- trackLengthInBad+=deltaPhi;
+ if (v<=fBadIROCthreshhold)
+ {
+ trackLengthInBad+=deltaPhi;
lengthFractionInBadSectors=1.;
}
- if (v<=fLowGainIROCthreshold && v>fBadIROCthreshhold)
+ if (v<=fLowGainIROCthreshold && v>fBadIROCthreshhold)
{
trackLengthInLowGain+=deltaPhi;
lengthFractionInBadSectors=1.;
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);
-
- if (lengthFractionInBadSectors>fMaxBadLengthFraction)
+
+ if (lengthFractionInBadSectors>fMaxBadLengthFraction)
{
//printf("%%%%%%%% %s kChamberLowGain\n",(layer==1)?"IROC":"OROC");
return kChamberLowGain;
}
-
+
return kChamberHighGain;
}
+
//_____________________________________________________________________________
Float_t AliTPCPIDResponse::MaxClusterRadius(const AliVTrack* track) const
{
if (!clusterMap) return 0.;
//from AliTPCParam, radius of first IROC row
- const Float_t rfirstIROC = 8.52249984741210938e+01;
+ const Float_t rfirstIROC = 8.52249984741210938e+01;
const Float_t padrowHeightIROC = 0.75;
const Float_t rfirstOROC0 = 1.35100006103515625e+02;
const Float_t padrowHeightOROC0 = 1.0;
return 0.0;
}
+
//_____________________________________________________________________________
Bool_t AliTPCPIDResponse::TrackApex(const AliVTrack* track, Float_t magField, Double_t position[3]) const
{
Double_t p[3];
track->GetPxPyPz(p);
Double_t r = 1./track->OneOverPt()/0.0299792458/magField; //signed - will determine the direction of b
- //printf("b: %.2f, x:%.2f, y:%.2f, pt: %.2f, px:%.2f, py%.2f, r: %.2f\n",magField, x[0],x[1],track->Pt(), p[0],p[1],r);
+ //printf("b: %.2f, x:%.2f, y:%.2f, pt: %.2f, px:%.2f, py%.2f, r: %.2f\n",magField, x[0],x[1],track->Pt(), p[0],p[1],r);
//find orthogonal vector (Gram-Schmidt)
Double_t alpha = (p[0]*x[0] + p[1]*x[1])/(p[0]*p[0] + p[1]*p[1]);
Double_t b[2];
- b[0]=x[0]-alpha*p[0];
- b[1]=x[1]-alpha*p[1];
-
+ b[0]=x[0]-alpha*p[0];
+ b[1]=x[1]-alpha*p[1];
+
Double_t norm = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]);
if (TMath::AreEqualAbs(norm,0.0,1e-10)) return kFALSE;
- b[0]/=norm;
- b[1]/=norm;
- b[0]*=r;
- b[1]*=r;
- b[0]+=x[0];
+ b[0]/=norm;
+ b[1]/=norm;
+ b[0]*=r;
+ b[1]*=r;
+ b[0]+=x[0];
b[1]+=x[1];
//printf("center: x:%.2f, y:%.2f\n",b[0],b[1]);
-
+
norm = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]);
if (TMath::AreEqualAbs(norm,0.0,1e-10)) return kFALSE;
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff