Here the comment:
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 10 Dec 2012 15:28:27 +0000 (15:28 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 10 Dec 2012 15:28:27 +0000 (15:28 +0000)
- Added extrapolation methods to get better eta maps for |eta| > ~0.85
- Improved splines and maps for existing periods + added more periods:
Data: 11a.pass2, 11a.pass4, 11a.pass1 (7TeV), 10{b,c,d,e}.pass2
MC: 12f1a, 12f1b, 10d1, 10f6a, 11b2 + periods with same maps/splines as e.g. 11b10a, 12gX
Jens

OADB/COMMON/PID/data/TPCPIDResponse.root
OADB/COMMON/PID/data/TPCetaMaps.root
STEER/STEERBase/AliPIDResponse.cxx
STEER/STEERBase/AliTPCPIDResponse.cxx
STEER/STEERBase/AliTPCPIDResponse.h

index 7a642ce0dedebcf5e8bbd8b1bafc0d43141198c0..3f6ce03a5179b129d7c498aa5f6f953794dff8d3 100644 (file)
Binary files a/OADB/COMMON/PID/data/TPCPIDResponse.root and b/OADB/COMMON/PID/data/TPCPIDResponse.root differ
index 91871000ca8da12591da6fefa011f164c90ab753..f59e6dc08b20f7415badb8bfc0449d22f5b19d10 100644 (file)
Binary files a/OADB/COMMON/PID/data/TPCetaMaps.root and b/OADB/COMMON/PID/data/TPCetaMaps.root differ
index 3beb0302ddb8954099c671e3f9864ecdd0577aee..80e55ea29e25af214a5bf49668c00091e321a674 100644 (file)
@@ -715,12 +715,14 @@ void AliPIDResponse::SetRecoInfo()
 //   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";
 }
 
 //______________________________________________________________________________
@@ -761,10 +763,10 @@ TH2D* AliPIDResponse::RefineHistoViaLinearInterpolation(TH2D* h, Double_t refine
   
   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; 
   
@@ -780,7 +782,7 @@ TH2D* AliPIDResponse::RefineHistoViaLinearInterpolation(TH2D* h, Double_t refine
       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.
@@ -852,6 +854,51 @@ TH2D* AliPIDResponse::RefineHistoViaLinearInterpolation(TH2D* h, Double_t refine
     }
   } 
   
+  
+  // 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;
@@ -868,18 +915,18 @@ void AliPIDResponse::SetTPCEtaMaps(Double_t refineFactorMapX, Double_t refineFac
   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;
   
index 9b8c515c736706e7ef19e007f007d8787b7c19a2..3faa7e377b4bd6eae8b1bf328d02182bce61b6f6 100644 (file)
@@ -78,7 +78,7 @@ AliTPCPIDResponse::AliTPCPIDResponse():
   //
   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(),
@@ -108,7 +108,7 @@ AliTPCPIDResponse::AliTPCPIDResponse(const Double_t *param):
   //
   for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=param[1];fResN2[i]=param[2];}
 }
-
+*/
 
 //_________________________________________________________________________
 AliTPCPIDResponse::~AliTPCPIDResponse()
@@ -153,12 +153,12 @@ AliTPCPIDResponse::AliTPCPIDResponse(const AliTPCPIDResponse& that):
   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);
   }
@@ -189,14 +189,14 @@ AliTPCPIDResponse& AliTPCPIDResponse::operator=(const AliTPCPIDResponse& that)
 
   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);
   }
@@ -536,9 +536,9 @@ Float_t AliTPCPIDResponse::GetNumberOfSigmas(const AliVTrack* track,
 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
@@ -871,7 +871,7 @@ Bool_t AliTPCPIDResponse::sectorNumbersInOut(Double_t* trackPositionInner,
                                              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
@@ -894,8 +894,8 @@ Bool_t AliTPCPIDResponse::sectorNumbersInOut(Double_t* trackPositionInner,
   }
   return kTRUE;
 }
-    
-    
+
+
 //_____________________________________________________________________________
 Int_t AliTPCPIDResponse::sectorNumber(Double_t phi) const
 {
@@ -904,6 +904,7 @@ Int_t AliTPCPIDResponse::sectorNumber(Double_t phi) const
   return TMath::Floor(phi/width);
 }
 
+
 //_____________________________________________________________________________
 void AliTPCPIDResponse::Print(Option_t* /*option*/) const
 {
@@ -911,6 +912,7 @@ void AliTPCPIDResponse::Print(Option_t* /*option*/) const
   fResponseFunctions.Print();
 }
 
+
 //_____________________________________________________________________________
 AliTPCPIDResponse::EChamberStatus AliTPCPIDResponse::TrackStatus(const AliVTrack* track, Int_t layer) const
 {
@@ -927,9 +929,9 @@ AliTPCPIDResponse::EChamberStatus AliTPCPIDResponse::TrackStatus(const AliVTrack
   /////////////////////////////////////////////////////////////////////////////
   //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();
@@ -968,19 +970,19 @@ AliTPCPIDResponse::EChamberStatus AliTPCPIDResponse::TrackStatus(const AliVTrack
   }
 
 
-  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;
@@ -1019,26 +1021,26 @@ AliTPCPIDResponse::EChamberStatus AliTPCPIDResponse::TrackStatus(const AliVTrack
   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.;
@@ -1050,16 +1052,17 @@ AliTPCPIDResponse::EChamberStatus AliTPCPIDResponse::TrackStatus(const AliVTrack
     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
 {
@@ -1069,7 +1072,7 @@ 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;
@@ -1084,6 +1087,7 @@ Float_t AliTPCPIDResponse::MaxClusterRadius(const AliVTrack* track) const
   return 0.0;
 }
 
+
 //_____________________________________________________________________________
 Bool_t AliTPCPIDResponse::TrackApex(const AliVTrack* track, Float_t magField, Double_t position[3]) const
 {
@@ -1093,23 +1097,23 @@ Bool_t AliTPCPIDResponse::TrackApex(const AliVTrack* track, Float_t magField, Do
   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;
  
index a39ff6e54344efa10be929b006099adb9e982ee6..1816551be898719147baf6e2d9a3a7a37dd4af5d 100644 (file)
@@ -30,7 +30,7 @@ class TSpline3;
 class AliTPCPIDResponse: public TNamed {
 public:
   AliTPCPIDResponse();
-  AliTPCPIDResponse(const Double_t *param);
+  //TODO Remove? AliTPCPIDResponse(const Double_t *param);
   AliTPCPIDResponse(const AliTPCPIDResponse&);
   AliTPCPIDResponse& operator=(const AliTPCPIDResponse&);
   virtual ~AliTPCPIDResponse();
@@ -67,6 +67,7 @@ public:
                                Double_t kp4,
                                Double_t kp5
                                );
+  //Better prevent user from setting fMIP != 50. because fMIP set fix to 50 for much other code:
   void SetMip(Float_t mip) { fMIP = mip; } // Set overall normalisation; mean dE/dx for MIP
   Double_t Bethe(Double_t bg) const;
   void SetUseDatabase(Bool_t useDatabase) { fUseDatabase = useDatabase;}