fastjet/contrib/RescalingClasses.hh

0001 // ConstituentSubtractor package

0002 // Questions/comments: peter.berta@cern.ch, Martin.Spousta@cern.ch, David.W.Miller@uchicago.edu, Rupert.Leitner@mff.cuni.cz

0003 //

0004 // Copyright (c) 2014-, Peter Berta, Martin Spousta, David W. Miller, Rupert Leitner

0005 //

0006 //----------------------------------------------------------------------

0007 // This file is part of FastJet contrib.

0008 //

0009 // It is free software; you can redistribute it and/or modify it under

0010 // the terms of the GNU General Public License as published by the

0011 // Free Software Foundation; either version 2 of the License, or (at

0012 // your option) any later version.

0013 //

0014 // It is distributed in the hope that it will be useful, but WITHOUT

0015 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY

0016 // or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public

0017 // License for more details.

0018 //

0019 // You should have received a copy of the GNU General Public License

0020 // along with this code. If not, see <http://www.gnu.org/licenses/>.

0021 //----------------------------------------------------------------------

0022
0023
0024 #ifndef __FASTJET_CONTRIB_CONSTITUENTSUBTRACTOR_RESCALINGCLASSES_HH__
0025 #define __FASTJET_CONTRIB_CONSTITUENTSUBTRACTOR_RESCALINGCLASSES_HH__
0026
0027
0028
0029 #include <fastjet/FunctionOfPseudoJet.hh>
0030 #include <iostream>
0031 #include <vector>
0032
0033 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh

0034
0035 namespace contrib{
0036
0037
0038   template<class T>
0039   class BackgroundRescalingYFromRoot : public FunctionOfPseudoJet<double> {
0040   public:
0041     ///

0042     /// construct a background rescaling function using ROOT TH1 histogram bin contents (rapidity binning)

0043     BackgroundRescalingYFromRoot() {}
0044     BackgroundRescalingYFromRoot(T* hist=0, bool interpolate=false) {_hist = hist; _interpolate=interpolate;}
0045
0046     // return the rescaling factor associated with this jet

0047     virtual double result(const PseudoJet & particle) const {
0048       if (!_hist){
0049     throw Error("BackgroundRescalingYFromRoot (from ConstituentSubtractor)  The histogram for rescaling not defined! ");
0050       }
0051       double rap = particle.rap();
0052       if (_interpolate) return _hist->Interpolate(rap);
0053       else{
0054         if (rap<_hist->GetXaxis()->GetBinLowEdge(1)) return _hist->GetBinContent(1);
0055         if (rap>=_hist->GetXaxis()->GetBinUpEdge(_hist->GetNbinsX())) return _hist->GetBinContent(_hist->GetNbinsX());
0056         int bin=_hist->FindBin(rap);
0057         return _hist->GetBinContent(bin);
0058       }
0059     }
0060
0061   private:
0062     T* _hist=0;
0063     bool _interpolate=false;
0064   };
0065
0066
0067
0068   template<class T>
0069   class BackgroundRescalingYPhiFromRoot : public FunctionOfPseudoJet<double> {
0070   public:
0071     ///

0072     /// construct a background rescaling function using ROOT TH2 histogram bin contents (rapidity on the x-axis and azimuth on the y-axis)

0073     BackgroundRescalingYPhiFromRoot() {}
0074     BackgroundRescalingYPhiFromRoot(T* hist=0, bool interpolate=false) {_hist = hist; _interpolate=interpolate;}
0075
0076     // return the rescaling factor associated with this jet

0077     virtual double result(const PseudoJet & particle) const {
0078       if (!_hist){
0079     throw Error("BackgroundRescalingYPhiFromRoot (from ConstituentSubtractor)  The histogram for rescaling not defined! ");
0080       }
0081       double rap = particle.rap();
0082       double phi = particle.phi();
0083       if (_interpolate) return _hist->Interpolate(rap,phi);
0084       else{
0085         int xbin=1;
0086         if (rap<_hist->GetXaxis()->GetBinLowEdge(1)) xbin=1;
0087         else if (rap>=_hist->GetXaxis()->GetBinUpEdge(_hist->GetNbinsX())) xbin=_hist->GetNbinsX();
0088         else xbin=_hist->GetXaxis()->FindBin(rap);
0089         int ybin=1;
0090         if (phi<_hist->GetYaxis()->GetBinLowEdge(1) || phi>_hist->GetYaxis()->GetBinUpEdge(_hist->GetNbinsY())){
0091           throw Error("BackgroundRescalingYPhiFromRoot (from ConstituentSubtractor)  The phi range of the histogram does not correspond to the phi range of the particles! Change the phi range of the histogram.");
0092         }
0093         else ybin=_hist->GetYaxis()->FindBin(phi);
0094         return _hist->GetBinContent(xbin,ybin);
0095       }
0096     }
0097
0098   private:
0099     T* _hist=0;
0100     bool _interpolate=false;
0101   };
0102
0103
0104
0105
0106
0107   template<class T>
0108   class BackgroundRescalingYFromRootPhi : public FunctionOfPseudoJet<double> {
0109   public:
0110     ///  Construct background rescaling function in rapidity and azimuth using ROOT TH1 histogram bin contents for the rapidity dependence and this parametrization for the azimuth:

0111
0112     ///  phi_term(phi) = 1 + 2 * v2^2 * cos(2*(phi-psi)) + 2 * v3^2 * cos(3*(phi-psi)) +  2 * v4^2 * cos(4*(phi-psi))

0113     ///  with four parameters v2, v3, v4, and psi.

0114
0115     ///  This product of the TH1 histogram and function is used to rescale the background which is subtracted such that one can correctly account

0116     ///  for the modulation of the UE due to rapidity dependence of the particle production

0117     ///  and/or due to the modulation in the azimuthal angle which is characteristic for heavy ion collisions.

0118     ///  The overall normalization of the rescaling function is arbitrary since it divides out in the calculation of position dependent rho (background is first demodulated to obtain unbiased position independent rho, and then it is modulated to obtain position dependent rho, see fastjet classes GridMedianBackgroundEstimator and JetMedianBackgroundEstimator for detailed calculation).

0119
0120     BackgroundRescalingYFromRootPhi() {}
0121     BackgroundRescalingYFromRootPhi(double v2, double v3, double v4, double psi, T* hist=0, bool interpolate=false){
0122       _v2=v2;
0123       _v3=v3;
0124       _v4=v4;
0125       _psi=psi;
0126       _hist=hist;
0127       _use_phi=true;
0128       _interpolate=interpolate;
0129       if (!_hist){
0130     std::cout << std::endl << std::endl << "ConstituentSubtractor::BackgroundRescalingYFromRootPhi WARNING: The histogram for rapidity rescaling is not defined!!! Not performing rapidity rescaling." << std::endl << std::endl << std::endl;
0131     _use_rap=false;
0132       }
0133       else _use_rap=true;
0134     }
0135
0136     ///

0137     /// Turn on or off the rapidity rescaling. Throwing in case true is set and no histogram is provided.

0138     void use_rap_term(bool use_rap){
0139       _use_rap=use_rap;
0140       if (!_hist && _use_rap){
0141     throw Error("BackgroundRescalingYFromRootPhi (from ConstituentSubtractor)  Requested rapidity rescaling, but the histogram for rescaling is not defined!");
0142       }
0143     }
0144
0145     ///

0146     /// Turn on or off the azimuth rescaling.

0147     void use_phi_term(bool use_phi){
0148       _use_phi=use_phi;
0149     }
0150
0151     ///

0152     /// Return the rescaling factor associated with this particle

0153     virtual double result(const PseudoJet & particle) const{
0154       double phi_term=1;
0155       if (_use_phi){
0156     double phi=particle.phi();
0157     phi_term=1 + 2*_v2*_v2*cos(2*(phi-_psi)) + 2*_v3*_v3*cos(3*(phi-_psi)) +  2*_v4*_v4*cos(4*(phi-_psi));
0158       }
0159       double rap_term=1;
0160       if (_use_rap){
0161     double y=particle.rap();
0162         if (_interpolate) rap_term=_hist->Interpolate(y);
0163         else{
0164           int bin=_hist->FindBin(y);
0165           rap_term=_hist->GetBinContent(bin);
0166         }
0167       }
0168
0169       return phi_term*rap_term;
0170     }
0171
0172   private:
0173     double _v2=0, _v3=0, _v4=0, _psi=0;
0174     bool _use_rap=false, _use_phi=false;
0175     T* _hist=0;
0176     bool _interpolate=false;
0177   };
0178
0179
0180
0181
0182
0183
0184   class BackgroundRescalingYPhi : public FunctionOfPseudoJet<double> {
0185   public:
0186     ///

0187     ///  Construct background rescaling function in rapidity and azimuth using this parameterization:

0188     ///

0189     ///  f(y,phi) = phi_term(phi) * rap_term(y)

0190     ///  where

0191     ///

0192     ///  phi_term(phi) = 1 + 2 * v2^2 * cos(2*(phi-psi)) + 2 * v3^2 * cos(3*(phi-psi)) +  2 * v4^2 * cos(4*(phi-psi))

0193     ///  with four parameters v2, v3, v4, and psi.

0194     ///

0195     ///  rap_term(y) = a1*exp(-pow(y,2)/(2*sigma1^2)) + a2*exp(-pow(y,2)/(2*sigma2^2))

0196     ///  with four parameters sigma1, sigma2, a1, and a2.

0197     ///

0198     ///  This function is used to rescale the background which is subtracted such that one can correctly account

0199     ///  for the modulation of the UE due to rapidity dependence of the particle production

0200     ///  and/or due to the modulation in the azimuthal angle which is characteristic for heavy ion collisions.

0201     ///  The overall normalization of function f is arbitrary since it divides out in the calculation of position dependent rho (background is first demodulated to obtain unbiased position independent rho, and then it is modulated to obtain position dependent rho, see fastjet classes GridMedianBackgroundEstimator and JetMedianBackgroundEstimator for detailed calculation).

0202
0203     BackgroundRescalingYPhi(): _v2(0), _v3(0), _v4(0), _psi(0), _a1(1), _sigma1(1000), _a2(0), _sigma2(1000), _use_rap(false), _use_phi(false) {}
0204     BackgroundRescalingYPhi(double v2, double v3, double v4, double psi, double a1, double sigma1, double a2, double sigma2);
0205
0206     void use_rap_term(bool use_rap);
0207     void use_phi_term(bool use_phi);
0208
0209     /// return the rescaling factor associated with this jet

0210     virtual double result(const PseudoJet & particle) const;
0211   private:
0212     double _v2, _v3, _v4, _psi, _a1, _sigma1, _a2, _sigma2;
0213     bool _use_rap, _use_phi;
0214   };
0215
0216
0217
0218   class BackgroundRescalingYPhiUsingVectorForY : public FunctionOfPseudoJet<double> {
0219   public:
0220     ///

0221     ///  Construct background rescaling function in rapidity and azimuth using this parameterization:

0222     ///

0223     ///  f(y,phi) = phi_term(phi) * rap_term(y)

0224     ///  where

0225     ///

0226     ///  phi_term(phi) = 1 + 2 * v2^2 * cos(2*(phi-psi)) + 2 * v3^2 * cos(3*(phi-psi)) +  2 * v4^2 * cos(4*(phi-psi))

0227     ///  with four parameters v2, v3, v4, and psi.

0228     ///

0229     ///  rap_term(y) = provided in a vector.

0230     ///

0231     /// The size of the input vector "values" for rapidity dependence is N bins and the corresponding binning should be specified in a separate input vector "rap_binning" of size (N+1). The bin boundaries must be in increasing order.

0232     ///

0233     ///  This function is used to rescale the background which is subtracted such that one can correctly account

0234     ///  for the modulation of the UE due to rapidity dependence of the particle production

0235     ///  and/or due to the modulation in the azimuthal angle which is characteristic for heavy ion collisions.

0236     ///  The overall normalization of function f is arbitrary since it divides out in the calculation of position dependent rho (background is first demodulated to obtain unbiased position independent rho, and then it is modulated to obtain position dependent rho, see fastjet classes GridMedianBackgroundEstimator and JetMedianBackgroundEstimator for detailed calculation).

0237
0238     BackgroundRescalingYPhiUsingVectorForY() {}
0239     BackgroundRescalingYPhiUsingVectorForY(double v2, double v3, double v4, double psi, std::vector<double> values, std::vector<double> rap_binning, bool interpolate=true);
0240
0241     void use_rap_term(bool use_rap);
0242     void use_phi_term(bool use_phi);
0243
0244     /// return the rescaling factor associated with this jet

0245     virtual double result(const PseudoJet & particle) const;
0246   private:
0247     double _v2=0, _v3=0, _v4=0, _psi=0;
0248     std::vector<double> _values;
0249     std::vector<double> _rap_binning;
0250     bool _use_rap=false, _use_phi=false;
0251     bool _interpolate=true;
0252   };
0253
0254
0255
0256   class BackgroundRescalingYPhiUsingVectors : public FunctionOfPseudoJet<double> {
0257   public:
0258     ///

0259     ///  Construct background rescaling function in rapidity and azimuth using dependence recorded in input object vector<vector<double> > values. Its size is N bins for rapidity and M bins for azimuth. The binning of the rapidity should be specified in a separate vector "rap_binning" of size (N+1), and similarly the binning of the azimuth  should be specified in a separate vector "phi_binning" of size (M+1). The bin boundaries must be in increasing order.

0260     ///

0261
0262     BackgroundRescalingYPhiUsingVectors(): _values(0), _rap_binning(0), _phi_binning(0) {}
0263     BackgroundRescalingYPhiUsingVectors(std::vector<std::vector<double> > values, std::vector<double> rap_binning, std::vector<double> phi_binning);
0264
0265     void use_rap_term(bool use_rap);
0266     void use_phi_term(bool use_phi);
0267
0268     /// return the rescaling factor associate

0269     virtual double result(const PseudoJet & particle) const;
0270   private:
0271     std::vector<std::vector<double> > _values;
0272     std::vector<double> _rap_binning;
0273     std::vector<double> _phi_binning;
0274     bool _use_rap, _use_phi;
0275   };
0276
0277 } // namespace contrib

0278
0279 FASTJET_END_NAMESPACE
0280
0281
0282 #endif   //__FASTJET_CONTRIB_CONSTITUENTSUBTRACTOR_RESCALINGCLASSES_HH__