for (unsigned int j=0;j<nGhosts; j++) mtMinusPt.push_back(jmbge->rho(ghosts[j])*ghosts_area[j]);
jmbge->set_jet_density_class(orig_density);
#if FASTJET_VERSION_NUMBER >= 30100
}
#endif
}
else { // a single bge, only rho requested
for (unsigned int j=0;j<nGhosts; j++) mtMinusPt.push_back(1e-200);
#if FASTJET_VERSION_NUMBER >= 30100
// In FJ3.1 and BGE with rho_m support, add a warning, similar to that in Subtractor
double const rho_m_warning_threshold = 1e-5;
if (_bge_rho->has_rho_m() && _bge_rho->rho_m()>rho_m_warning_threshold*_bge_rho->rho()){
_warning_unused_rhom.warn("ConstituentSubtractor:: Background estimator indicates non-zero rho_m, but the constituent subtractor does not use rho_m information; consider calling set_common_bge_for_rho_and_rhom(true) to include the rho_m information");
if (_distance==ConstituentSubtractor::deltaR) max_distance_transformed=pow(_max_distance,2);
if (_distance==ConstituentSubtractor::angle) max_distance_transformed=-cos(_max_distance);
std::vector<std::pair<double,int> > CS_distances; // the first element is the distance, the second element is only the index in the vector used for sorting
if (fabs(_alpha)>1e-5) pt_factor=pow(particles_pt[i],_alpha);
double polarAngle_factor=1.;
if (fabs(_polarAngleExp)>1e-5) polarAngle_factor=pow(particles_pt[i]/sqrt(particles_sorted[i].pt2()+particles_sorted[i].pz()*particles_sorted[i].pz()),_polarAngleExp);
if (backgroundProxies_fraction_of_pt[i]>0) subtracted_pt=backgroundProxies_pt[i]*backgroundProxies_fraction_of_pt[i];
double subtracted_mtMinusPt=0;
if (backgroundProxies_fraction_of_mtMinusPt[i]>0) subtracted_mtMinusPt=(backgroundProxies_mt[i]-backgroundProxies_pt[i])*backgroundProxies_fraction_of_mtMinusPt[i];
if (!_bge_rho) throw Error("ConstituentSubtractor::set_common_bge_for_rho_and_rhom() is not allowed when _bge_rho is not set!");
if (value){
if (_bge_rhom) throw Error("ConstituentSubtractor::set_common_bge_for_rho_and_rhom() is not allowed in the presence of an existing background estimator for rho_m.");
if (_externally_supplied_rho_rhom) throw Error("ConstituentSubtractor::set_common_bge_for_rho_and_rhom() is not allowed when supplying externally the values for rho and rho_m.");
if (!jmbge) throw Error("ConstituentSubtractor::set_common_bge_for_rho_and_rhom() is currently only allowed for background estimators of JetMedianBackgroundEstimator type.");
#if FASTJET_VERSION_NUMBER >= 30100
}
#endif
_common_bge=value;
_do_mass_subtraction=value;
}
// setting this to true will result in rho_m being estimated using bge_rhom->rho_m() instead of bge_rhom->rho()
throw Error("ConnstituentSubtractor::use_rhom_from_bge_rhom() can only be used with FastJet >=3.1.");
#else
if (!_bge_rhom) throw Error("ConstituentSubtractor::use_rhom_from_bge_rhom() requires a background estimator for rho_m.");
if (!(_bge_rhom->has_rho_m())) throw Error("ConstituentSubtractor::use_rhom_from_bge_rhom() requires rho_m support for the background estimator for rho_m.");
+ - added new functions "get_ghosts" and "get_ghosts_area".
+ - now the user can modify the ghosts outside the CS after constructing them with function "construct_ghosts_uniformly" and getting them with "get_ghosts" and "get_ghosts_area"
+ - added new rescaling class useful for heavy ion events. Now the rapidity dependence can be taken from a Root TH1 object.
+
2017/01/20: release of version 1.1.4:
- fixed few bugs concerning the correction of masses
2016/11/14: release of version 1.1.3:
- updated rho rescaling functions
- new function for pileup subtraction: sequential_subtraction
- possibility to change the computation of distance between particles and background proxies
2016/03/11: release of version 1.1.2:
- added two FunctionOfPseudojet<object> classes for rho rescaling: one using root TH1 rapidity parametrization, and the other for heavy ions rapidity-azimuth parametrization
2016/02/09: release of version 1.1.1:
- faster algortihm for subtraction
- added new function for Constituent Subtraction using tracking information
2015/11/05: release of version 1.1.0:
- added FastJet v3.1 rho_m support
- simplified way for Constituent Subtraction of the whole event
- updated the meaning of parameter max_deltaR
- added new parameter to the definition of the distance
- /// with four parameters sigma1, sigma2, a1, and a2.
-
- /// This function is used to rescale the background which is subtracted such that one can correctly account
- /// for the modulation of the UE due to rapidity dependence of the particle production
- /// and/or due to the modulation in the azimuthal angle which is characteristic for heavy ion collisions.
- /// 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).
+ class BackgroundRescalingYFromRootPhi : public FunctionOfPseudoJet<double> {
+ public:
+ /// Construct background rescaling function in rapidity and azimuth using ROOT TH1 histogram bin contents for the rapidity dependence and this parametrization for the azimuth:
+ /// This product of the TH1 histogram and function is used to rescale the background which is subtracted such that one can correctly account
+ /// for the modulation of the UE due to rapidity dependence of the particle production
+ /// and/or due to the modulation in the azimuthal angle which is characteristic for heavy ion collisions.
+ /// 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).
+ 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;
+ _use_rap=false;
+ }
+ else _use_rap=true;
+ }
+
+ /// Turn on or off the rapidity rescaling. Throwing in case true is set and no histogram is provided.
+ void use_rap_term(bool use_rap){
+ _use_rap=use_rap;
+ if (!_hist && _use_rap){
+ std::cout << "ConstituentSubtractor::BackgroundRescalingYFromRootPhi ERROR: Requested rapidity rescaling, but the histogram for rescaling is not defined!!! Throwing." << std::endl;
+ throw;
+ }
+ }
+
+ /// Turn on or off the azimuth rescaling.
+ void use_phi_term(bool use_phi){
+ _use_phi=use_phi;
+ }
+
+ /// Return the rescaling factor associated with this particle
+ /// with four parameters sigma1, sigma2, a1, and a2.
+
+ /// This function is used to rescale the background which is subtracted such that one can correctly account
+ /// for the modulation of the UE due to rapidity dependence of the particle production
+ /// and/or due to the modulation in the azimuthal angle which is characteristic for heavy ion collisions.
+ /// 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).
/// A class to perform subtraction of background, e.g. pileup, from a set of particles at particle-level. The output is a jet or the whole event with corrected constituents.
///
/// This class corrects the input particles for background contamination with the algorithm described in:
/// Peter Berta, Martin Spousta, David W. Miller, Rupert Leitner [arXiv:1403.3108]
///
/// For individual jet background subtraction, see example.cc
/// For whole event background subtraction before jet clustering, see example_whole_event.cc
///
/// The distance used for matching between particle i and ghost k is defined as:
/// Constructor that takes a pointer to a background estimator for rho and optionally a pointer to a background estimator for rho_m. If the latter is not supplied, rho_m is assumed to always be zero (this behaviour can be changed by calling use_common_bge_for_rho_and_rhom).
/// do the constituent subtraction for the input particles using the provided background proxies. The output is a vector with corrected particles - particles with zero corrected pt are removed.
/// do the subtraction of the whole event - more user-friendly approach. The particles with |eta|>max_eta are discarded at the beginning, i.e. they are not used, nor returned. The ghosts are added automatically inside this function up to max_eta.
/// do the subtraction of the whole event using the tracking information for charged particles, i.e. the 4-momenta of charged particles from signal vertex, and 4-momenta of charged particles from background. The user can set the scaling of charged particles from background and signal using parameters charged_background_scale (CBS) and charged_signal_scale (CSS). These scales are useful if one assumes correlation between charged and neutral particles or in case the inputs from calorimeter are miscalibrated wrt tracks. In case CBS=CSS=0, the input charged particles are not used. In case CBS=CSS=1, the input charged particles are not scaled. Recommending to try several combinations for CBS and CSS from range [0.8, 1.5]. It is no more necessary to provide background estimator. The GridMedianBackgroundEstimator is used - probably more flexibility will be added in the future. The rescaling function for background estimator can be also provided - the rescaling function will be used for the event after subtracting charged scaled particles.
/// do the subtraction of the whole event in two steps: in the first step, very small max_distance is used (0.1), and in the second step larger max_distance is used (1) to subtract the remaining background proxies
/// Set the pointer to a background estimator for rho and optionally a pointer to a background estimator for rho_m. If the latter is not supplied, rho_m is assumed to always be zero (this behaviour can be changed by calling use_common_bge_for_rho_and_rhom).
/// when only one background estimator, bge_rho, is specified, calling this function with argument true, causes rho_m to be calculated from the same background estimator as rho, instead of being set to zero. Currently this only works if the estimator is a JetMedianBackgroundEstimator or other estimator which has such function.
/// when two background estimators are used (one for rho, the second for rho_m), setting this to true will result in rho_m being estimated using bge_rhom->rho_m() instead of bge_rhom->rho().
void set_use_bge_rhom_rhom(bool value=true);
/// specify if also the mass term sqrt(pT^2+m^2)-pT should be corrected during the subtraction procedure. This is automatically set when calling functions set_common_bge_for_rho_and_rhom, set_scalar_background_density or set_use_bge_rhom_rhom
void set_do_mass_subtraction(bool value=true);
/// function to change the alpha-parameter figuring in the distance measure deltaR. The larger the alpha, the more are preferred to be corrected the low pt particles. The default value is 0, i.e. by default the standard deltaR definition is used: deltaR=sqrt(deltay^2 + deltaphi^2)
void set_alpha(double alpha);
/// function to change the parameter polarAngleExp
void set_polarAngleExp(double polarAngleExp);
/// function to change the parameter ghost_area
void set_ghost_area(double ghost_area);
/// function to change the distance type
void set_distance_type(Distance distance);
/// function to change the free parameter max_distance. The distance is defined by enum Distance. The particle-ghost pairs with distance>max_distance are not used. When max_distance<=0, the max_distance parameter is not used (no upper limit on distance). The default value is -1, i.e. by default there is no upper limit for possible distance values.
void set_max_distance(double max_distance);
/// function to change the free parameter max_distance_sequential for sequantial subtraction. This is the smaller maximal deltaR distance used in the first subtraction.
/// function to return the maximal deltaR distance used in the subtraction
double get_max_distance();
/// This function returns true if the first argument is smaller than the second argument, otherwise returns false. The comparison is done only on the first element in the two pairs. This function is used to sort in ascending order the deltaR values for each pair particle-ghost while keeping track of particles and ghosts
+2017-07-11 Peter Berta <berta@ipnp.troja.mff.cuni.cz>
+
+ * VERSION:
+ - Set version number to 1.1.6
+
+ * ConstituentSubtractor.cc
+ * ConstituentSubtractor.hh
+ - Implemented suggestion from Marta Verweij:
+ - Added new functions: "get_ghosts" and "get_ghosts_area". They return the vector of ghosts and the corresponding ghost areas used for whole event subtraction.
+ - Changed the function "construct_ghosts_uniformly" to be public
+ - Now the user does not need to provide the BackgroundEstimator to CS, but he/she can do background estimation independently and then assign the estimated background to ghosts. The procedure is the following:
+ - construct the ghosts with function "construct_ghosts_uniformly"
+ - get them with functions "get_ghosts" and "get_ghosts_area"
+ - do background estimation independently and assign the corresponding background estimates to ghosts
+ - use the function "do_subtraction" to subtract the modified ghosts from real particles.
+
+
+ * RescalingClasses.cc
+ * RescalingClasses.hh:
+ - Added new rescaling class template BackgroundRescalingYFromRootPhi - useful mainly for heavy ion events. The rapidity dependance is taken from a Root TH1 object provided by the user. The next parameters are used to parametrize the eliptic flow.
+
+
+ * example_background_rescaling.cc
+ - Added example for usage of the new rescaling class BackgroundRescalingYFromRootPhi (commented at the beginning of the macro).
+
+
+
2017-01-25 Gavin Salam <gavin.salam@cern.ch>
* VERSION:
updated VERSION to 1.1.5
* Makefile:
- resolved issue with make fragile-shared, caused by repeated
SRCS lines in the Makefile. Also removed other duplicates.
2017-01-20 Peter Berta <berta@ipnp.troja.mff.cuni.cz>
* VERSION:
- Set version number to 1.1.4
* ConstituentSubtractor.cc
- fixed rounding issues for zero masses, when doing the correction of term mtMinusPt
- more precise evaluation of corrected 4-momenta, which affected the mass of corrected particles with low masses
- added new function set_do_mass_subtraction
* ConstituentSubtractor.hh
- fixed initialization of member variable "_do_mass_subtraction"
- updated reference file. Added also all corrected particles.
2016-11-14 Peter Berta <berta@ipnp.troja.mff.cuni.cz>
* VERSION:
- Set version number to 1.1.3
* ConstituentSubtractor.cc
* ConstituentSubtractor.hh:
- Fixed bug in "do_subtraction" function for the unsubtracted background proxies.
- Added new function "sequential_subtraction" - for testing purposes so far, example will be added later. The ConstituentSubtraction is performed twice: first with small maximal allowed deltaR distance, then the background proxies are again uniformly redistributed and the ConstituentSubtraction runs second time with higher allowed maximal deltaR
- Added new function "set_distance_type". With this function, the user can change the definition of distance between particles and background proxies. There are two possibilities using enums Distance: deltaR (longitudinal Lorentz invariant distance in rapidity vs. phi) or angle (in Euclidean space). The default distance is deltaR, i.e. the same definition of distance as in previous versions.
* Makefile:
- Removed dependence on Root. Using template in RescalingClasses now.
* RescalingClasses.cc
* RescalingClasses.hh:
- Removed dependence on Root for the rescaling using TH1 histogram. Using template now. Also the name of the class changed. Old name: BackgroundRescalingYTH1. New name: BackgroundRescalingYFromRoot
- Changed the name of the rescaling class for heavy ion events. Old name: BackgroundRescalingYPhiHeavyIon. New name: BackgroundRescalingYPhi. Also the description of the rescaling function is updated.
* example_rescaling_using_TH1.cc
* example_rescaling_using_TH1.ref:
- Removed
* example_background_rescaling.cc
* example_background_rescaling.ref:
- Added new example for usage of Heavy Ion background rescaling function
- The file example_background_rescaling.cc contains also commented example for usage of rapidity rescaling according to root's TH1 distribution
* example_whole_event.cc:
- Updated this example with usage of the "set_distance_type" function.
2016-03-11 Peter Berta <berta@ipnp.troja.mff.cuni.cz>
* VERSION:
set version number to 1.1.2
* Makefile:
- added CXXFLAGS and LDFLAGS in case the "ROOTSYS" variable is defined from ROOT. In case it is not defined, the standard usage still works, just the RescalingClasses are not compiled.
* RescalingClasses.cc:
* RescalingClasses.hh:
- New files: contain classes for rho estimation rescaling. Two classes were added:
- BackgroundRescalingYTH1 - input is a TH1 object (stabdard histogram in ROOT) containg the actual rapidity dependence of backgound pt
- BackgroundRescalingYPhiHeavyIon - input are 8 parameters. The first four parameters parametrize the elliptic flow. The next four parameters parametrize the rapidity dependence using two Gaussian distributions
* example_rescaling_using_TH1.cc:
* example_rescaling_using_TH1.ref:
- New files: example in which the ConstituentSubtractor is used for the whole event subtraction using rapidity rescaling for the rho estimation. The input for rapidity rescaling is a TH1 object
2016-02-09 Peter Berta <berta@ipnp.troja.mff.cuni.cz>
* VERSION:
set version number to 1.1.1
* ConstituentSubtractor.hh:
* ConstituentSubtractor.cc:
- Faster algorithm for subtraction in the "do_subtraction" member function.
- New member function: subtract_event_using_charged_info
- whole event subtraction
- 4 additional inputs: two vectors of PseudoJets for charged background and charged signal, and two doubles, charged background scale (CBS) and charged signal scale (CSS)
- the scales CBS and CSS scale the input charged particles. Useful if one assumes correlation between charged and neutral particles or in case the inputs from calorimeter are miscalibrated wrt tracks. In case CBS=CSS=0, the input charged particles are not used. In case CBS=CSS=1, the input charged particles are not scaled. Recommending to try several combinations for CBS and CSS from range [0.8, 1.5].
- the background estimation is built in this function using GridMedianBackgroundEstimator with possibility to rescale
- see example file "example_whole_event_using_charged_info.cc" for usage
* example_whole_event_using_charged_info.cc:
* example_whole_event_using_charged_info.ref:
- New files: example in which the ConstituentSubtractor is used for the whole event subtraction using tracking information for charged particles
* functions.h:
- Updated function "read_event". Now, obtaining also the vectors of charged particles originating from signal and background.
2015-11-05 Peter Berta <berta@ipnp.troja.mff.cuni.cz>
* VERSION:
set version number to 1.1.0
* Makefile:
- added "-fPIC" CXX flag
* ConstituentSubtractor.hh:
* ConstituentSubtractor.cc:
- updated output for result(PseudoJet): copying the "user_index" from the input jet and its constituents to the output corrected jet and its constituents
- discarded the original max_deltaR parameter, and replaced it by parameter max_standardDeltaR defined as standardDeltaR=sqrt(deltay^2 + deltaphi^2). This new parameter restricts the usage of ghost-particle pairs which have too large distance standardDeltaR, i.e. not the distance deltaR which is the standardDeltaR multiplied by other factors as in the previous version.
- added a new multiplicative factor (sin(theta_i))^polarAngleExp into the definition of distance for a ghost-particle pair. Angle theta_i is the polar angle of partile i. The polarAngleExp is a free parameter with default value polarAngleExp=0, i.e. the default definition of deltaR has not changed wrt previous version. It can be changed using set_polarAngleExp member function
- changed the structure of the package and its member functions
- the standard usage with result(PseudoJet) is the same and it gives the same results as in previous version
- new member functions are available to simplify the whole event subtraction: subtract_event, set_background_estimator, set_scalar_background_density, see example file example_whole_event.cc
- similar changes as in GenericSubtractor due to FJ>=3.1:
- added set_common_bge_for_rho_and_rhom() and removed use_common_bge_for_rho_and_rhom()
- added set_use_bge_rhom_rhom() to tweak how rho_m is obtained when ConstituentSubtractor has been constructed with two estimators
- altered common_bge_for_rho_and_rhom() and added use_bge_rhom_rhom() to retrieve the behaviour
- added a warning when rhom is non-zero and unused
- Updated the terminology:
- ghosts - soft particles distributed with high density
- background proxies - particles used in the "do_subtraction" member function. They can be obtained from ghosts or in other way. Previously, these particles were also called "ghosts"
* example_whole_event.cc:
* example_whole_event.ref:
-new files: example in which the ConstituentSubtractor is used for the whole event
* functions.h:
-new file: contains useful functions for the example files
* example.hh:
-removed. Content moved to file functions.h
2014-04-06 Peter Berta <berta@ipnp.troja.mff.cuni.cz>
+ /// rap_term(y) = bin content of the histogram at position y
+ ///
+ /// You need to set the parameters event-by-event. Example where the TH1 histogram is filled with a random distribution. Do not use this histogram!!! Fill your own histogram with all particles (topoclusters) weighted by their pt using minimum bias events.
+ rescaling.use_rap_term(true); // this is useful to make sure that the histogram with rapidity dependence is defined. Set to false, if you do not want to use the rapidity rescaling
+ */
+
+
+ ///**** rescaling in heavy ion events using parametrized rapidity dependence: ****
// Set the seven parameters for rescaling using function
rescaling.use_rap_term(false); // set to true, if you have derived also the rapidity terms for the rescaling
///**** rescaling using rapidity dependence stored in root TH1 object: ****
// find the rapidity distribution of pileup particles from minimum bias events in a separate run. Fill a root TH1 histogram with this distribution.
// Here as an example where the TH1 histogram is filled with a random distribution. Do not use this!!! Fill it with all particles (topoclusters) weighted by their pt using minimum bias events.
subtractor.set_distance_type(contrib::ConstituentSubtractor::deltaR); // free parameter for the type of distance between particle i and ghost k. There are two options: "deltaR" or "angle" which are defined as deltaR=sqrt((y_i-y_k)^2+(phi_i-phi_k)^2) or Euclidean angle between the momenta
subtractor.set_max_distance(1); // free parameter for the maximal allowed distance between particle i and ghost k
subtractor.set_alpha(2); // free parameter for the distance measure (the exponent of particle pt). The larger the parameter alpha, the more are favoured the lower pt particles in the subtraction process
subtractor.set_ghost_area(0.01); // free parameter for the density of ghosts. The smaller, the better - but also the computation is slower.
// event loop
// read in input particles
vector<PseudoJet> hard_event, full_event;
read_event(hard_event, full_event);
double maxEta=4; // specify the maximal rapidity for the particles used in the subtraction
double maxEta_jet=3; // the maximal rapidity for selected jets
AreaDefinition area_def(active_area_explicit_ghosts,GhostedAreaSpec(maxEta,1)); // the area definiton is used only for the jet backgroud estimator. It is not important for the ConstituentSubtractor when subtracting the whole event - this is not true when subtracting the individual jets
bge_rho.set_jet_density_class(scalarPtDensity); // this changes the computation of pt of patches from vector sum to scalar sum. The scalar sum seems more reasonable.
// setting the TH1 rescaling:
bge_rho.set_rescaling_class(&rescaling);
bge_rho.set_particles(full_event);
subtractor.set_background_estimator(&bge_rho);
// this sets the same background estimator to be used for deltaMass density, rho_m, as for pt density, rho:
subtractor.set_common_bge_for_rho_and_rhom(true); // for massless input particles it does not make any difference (rho_m is always zero)