fastjet/contrib/RescalingClasses.hh

0001 // ConstituentSubtractor package

0002 // Questions/comments: berta@ipnp.troja.mff.cuni.cz, 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 #ifndef __FASTJET_CONTRIB_CONSTITUENTSUBTRACTOR_RESCALINGCLASSES_HH__
0024 #define __FASTJET_CONTRIB_CONSTITUENTSUBTRACTOR_RESCALINGCLASSES_HH__
0025
0026
0027
0028 #include <fastjet/FunctionOfPseudoJet.hh>
0029 #include <iostream>
0030 #include <vector>
0031
0032 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh

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

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

0042     BackgroundRescalingYFromRoot(): _hist(0) {}
0043     BackgroundRescalingYFromRoot(T* hist=0) {_hist = hist;}
0044
0045     // return the rescaling factor associated with this jet

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

0067     /// construct a background rescaling function using ROOT TH2 histogram bin contents (rapidity vs azimuth binning)

0068     BackgroundRescalingYPhiFromRoot(): _hist(0) {}
0069     BackgroundRescalingYPhiFromRoot(T* hist=0) {_hist = hist;}
0070
0071     // return the rescaling factor associated with this jet

0072     virtual double result(const PseudoJet & particle) const {
0073       if (!_hist){
0074     throw Error("BackgroundRescalingYPhiFromRoot (from ConstituentSubtractor)  The histogram for rescaling not defined! ");
0075       }
0076       double rap = particle.rap();
0077       double phi = particle.phi();
0078       int xbin=1;
0079       if (rap<_hist->GetXaxis()->GetBinLowEdge(1)) xbin=1;
0080       else if (rap>=_hist->GetXaxis()->GetBinUpEdge(_hist->GetNbinsX())) xbin=_hist->GetNbinsX();
0081       else xbin=_hist->GetXaxis()->FindBin(rap);
0082       int ybin=1;
0083       if (phi<_hist->GetYaxis()->GetBinLowEdge(1) || phi>_hist->GetYaxis()->GetBinUpEdge(_hist->GetNbinsY())){
0084     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.");
0085       }
0086       else ybin=_hist->GetYaxis()->FindBin(phi);
0087       return _hist->GetBinContent(xbin,ybin);
0088     }
0089
0090   private:
0091     T* _hist;
0092   };
0093
0094
0095
0096
0097
0098   template<class T>
0099   class BackgroundRescalingYFromRootPhi : public FunctionOfPseudoJet<double> {
0100   public:
0101     ///  Construct background rescaling function in rapidity and azimuth using ROOT TH1 histogram bin contents for the rapidity dependence and this parametrization for the azimuth:

0102
0103     ///  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))

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

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

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

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

0109     ///  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).

0110
0111     BackgroundRescalingYFromRootPhi(): _v2(0), _v3(0), _v4(0), _psi(0), _use_rap(false), _use_phi(false), _hist(0) {}
0112     BackgroundRescalingYFromRootPhi(double v2, double v3, double v4, double psi, T* hist=0){
0113       _v2=v2;
0114       _v3=v3;
0115       _v4=v4;
0116       _psi=psi;
0117       _hist=hist;
0118       _use_phi=true;
0119       if (!_hist){
0120     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;
0121     _use_rap=false;
0122       }
0123       else _use_rap=true;
0124     }
0125
0126     ///

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

0128     void use_rap_term(bool use_rap){
0129       _use_rap=use_rap;
0130       if (!_hist && _use_rap){
0131     throw Error("BackgroundRescalingYFromRootPhi (from ConstituentSubtractor)  Requested rapidity rescaling, but the histogram for rescaling is not defined!");
0132       }
0133     }
0134
0135     ///

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

0137     void use_phi_term(bool use_phi){
0138       _use_phi=use_phi;
0139     }
0140
0141     ///

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

0143     virtual double result(const PseudoJet & particle) const{
0144       double phi_term=1;
0145       if (_use_phi){
0146     double phi=particle.phi();
0147     phi_term=1 + 2*_v2*_v2*cos(2*(phi-_psi)) + 2*_v3*_v3*cos(3*(phi-_psi)) +  2*_v4*_v4*cos(4*(phi-_psi));
0148       }
0149       double rap_term=1;
0150       if (_use_rap){
0151     double y=particle.rap();
0152     int bin=_hist->FindBin(y);
0153     rap_term=_hist->GetBinContent(bin);
0154       }
0155
0156       return phi_term*rap_term;
0157     }
0158
0159   private:
0160     double _v2, _v3, _v4, _psi;
0161     bool _use_rap, _use_phi;
0162     T* _hist;
0163   };
0164
0165
0166
0167
0168
0169
0170   class BackgroundRescalingYPhi : public FunctionOfPseudoJet<double> {
0171   public:
0172     ///

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

0174     ///

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

0176     ///  where

0177     ///

0178     ///  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))

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

0180     ///

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

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

0183     ///

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

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

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

0187     ///  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).

0188
0189     BackgroundRescalingYPhi(): _v2(0), _v3(0), _v4(0), _psi(0), _a1(1), _sigma1(1000), _a2(0), _sigma2(1000), _use_rap(false), _use_phi(false) {}
0190     BackgroundRescalingYPhi(double v2, double v3, double v4, double psi, double a1, double sigma1, double a2, double sigma2);
0191
0192     void use_rap_term(bool use_rap);
0193     void use_phi_term(bool use_phi);
0194
0195     /// return the rescaling factor associated with this jet

0196     virtual double result(const PseudoJet & particle) const;
0197   private:
0198     double _v2, _v3, _v4, _psi, _a1, _sigma1, _a2, _sigma2;
0199     bool _use_rap, _use_phi;
0200   };
0201
0202
0203
0204   class BackgroundRescalingYPhiUsingVectorForY : public FunctionOfPseudoJet<double> {
0205   public:
0206     ///

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

0208     ///

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

0210     ///  where

0211     ///

0212     ///  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))

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

0214     ///

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

0216     ///

0217     /// 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.

0218     ///

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

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

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

0222     ///  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).

0223
0224     BackgroundRescalingYPhiUsingVectorForY(): _v2(0), _v3(0), _v4(0), _psi(0), _values(0), _rap_binning(0), _use_rap(false), _use_phi(false) {}
0225     BackgroundRescalingYPhiUsingVectorForY(double v2, double v3, double v4, double psi, std::vector<double> values, std::vector<double> rap_binning);
0226
0227     void use_rap_term(bool use_rap);
0228     void use_phi_term(bool use_phi);
0229
0230     /// return the rescaling factor associated with this jet

0231     virtual double result(const PseudoJet & particle) const;
0232   private:
0233     double _v2, _v3, _v4, _psi;
0234     std::vector<double> _values;
0235     std::vector<double> _rap_binning;
0236     bool _use_rap, _use_phi;
0237   };
0238
0239
0240
0241   class BackgroundRescalingYPhiUsingVectors : public FunctionOfPseudoJet<double> {
0242   public:
0243     ///

0244     ///  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.

0245     ///

0246
0247     BackgroundRescalingYPhiUsingVectors(): _values(0), _rap_binning(0), _phi_binning(0) {}
0248     BackgroundRescalingYPhiUsingVectors(std::vector<std::vector<double> > values, std::vector<double> rap_binning, std::vector<double> phi_binning);
0249
0250     void use_rap_term(bool use_rap);
0251     void use_phi_term(bool use_phi);
0252
0253     /// return the rescaling factor associate

0254     virtual double result(const PseudoJet & particle) const;
0255   private:
0256     std::vector<std::vector<double> > _values;
0257     std::vector<double> _rap_binning;
0258     std::vector<double> _phi_binning;
0259     bool _use_rap, _use_phi;
0260   };
0261
0262
0263
0264 } // namespace contrib

0265
0266 FASTJET_END_NAMESPACE
0267
0268
0269 #endif   //__FASTJET_CONTRIB_CONSTITUENTSUBTRACTOR_RESCALINGCLASSES_HH__