MultivariateAnalysis/src/TopTopologicalVariables.cc

/* File TopTopologicalVariables.hpp
 *
 * Created       : Tue Oct  4 16:58:28 CDT 2005
 * Author        : Supriya JAIN, sjain@fnal.gov
 *
 * Purpose      : Class containing generic methods for 
 *                       calculating topological variables 
 *                       (Instantiate using as argument: vector<TMBLorentzVector> 
 *                        of all objects for which you want any variable 
 *                        to be calculated)  
 *
 * Last modified : Amnon Harel, 06 Mar 2006 (const correctness, cache eigenvalues)
 * Comments      : 
 */

#include "TMatrixD.h"
#include "TVectorD.h"
#include "TRandom.h"

#include "LJMet/MultivariateAnalysis/interface/TopTopologicalVariables.h"
#include "LJMet/MultivariateAnalysis/interface/AnglesUtil.h"
#include <iostream>

using namespace std;

namespace top_cafe {
  
  /// Copy constructor
  TopTopologicalVariables::TopTopologicalVariables(const TopTopologicalVariables& that)
    : _myobjects(that._myobjects)
  {
    _pv = 0;
    if (that._pv) _pv = new TVectorD (*(that._pv));
  }

  /// Asignment operator
  TopTopologicalVariables& TopTopologicalVariables::operator= (const TopTopologicalVariables& that)
  {
    if (this != &that)
      {
        _myobjects = that._myobjects;
        if (_pv) delete _pv;
        if (that._pv) _pv = new TVectorD (*(that._pv));
      }
    return *this;
  }

  /// destructor
  TopTopologicalVariables::~TopTopologicalVariables()
  {
    if (_pv) delete _pv;
  }

  /// this method returns the centrality of a group of objects
  double TopTopologicalVariables::Centrality() const 
  {
    return Ht()/H();  
  } // Centrality()

    /// this method returns the aplanarity of a group of objects
    double TopTopologicalVariables::Aplanarity() const 
    {
      ensurePV();
      return 1.5 * (*_pv)[2]; // alternative syntax: ... * _pv->operator[](2)
    } // Aplanarity()
  
  /// this method returns the sphericity of a group of objects
  double TopTopologicalVariables::Sphericity() const
  {
    ensurePV();
    return 1.5 * ( (*_pv)[2] + (*_pv)[1] );
  } // Sphericity()
    
  /// this method returns the C of a group of objects (~momentum elipsiod surface area)
  double TopTopologicalVariables::C() const
  {
    ensurePV();
    return 3 * ( (*_pv)[0] * (*_pv)[1] + (*_pv)[0] * (*_pv)[2] + (*_pv)[1] * (*_pv)[2]);
  } // C()
    
  /// this method returns the D of a group of objects (~momentum elipsiod volume)
  double TopTopologicalVariables::D() const
  {
    ensurePV();
    return 27 * (*_pv)[0] * (*_pv)[1] * (*_pv)[2];
  } // D()
    
    /// this method returns the Momentum Tensor Eigenvalues of a group of objects
    TVectorD TopTopologicalVariables::GetMomentumTensorEigenvalues() const
    {
      ensurePV ();
      return *_pv; // implicit copy and cast
    } // GetMomentumTensorEigenvalues()

  /// internal method to prepare and cache the eigenvectors
  void TopTopologicalVariables::ensurePV () const
  {
    if (_pv) return;

    TMatrixD MomentumTensor(3,3);
    
    Double_t p2_sum=0.0;
    
    // for each _myobjects:
    for ( unsigned int k=0; k<_myobjects.size(); k++ ) {
      
      for ( int i=0; i<3; i++ ) // px, py, pz
        for ( int j=0; j<3; j++ ) // px, py, pz
          MomentumTensor(i,j) += _myobjects[k][i]*_myobjects[k][j];
      
      // add the 3-momentum squared to the sum
      p2_sum += _myobjects[k].Mag32();
    } // Loop over _myobjects
    
      // Divide the sums with the p2 sum
    if ( p2_sum != 0 )
      for ( int i=0; i<3; i++ ) // px, py, pz
        for ( int j=0; j<3; j++ ) // px, py, pz
          MomentumTensor(i,j) /= p2_sum;
    
    _pv = new TVectorD (3);
    MomentumTensor.EigenVectors(*_pv);
  }

  /// this method returns the  Pt of a group of objects
  double TopTopologicalVariables::Pt() const {
    // initialize
    TMBLorentzVector Sum(0.0,0.0,0.0,0.0);
    double pt;
    for (unsigned int i=0; i<_myobjects.size(); i++)
      Sum += _myobjects[i];
    pt = Sum.Pt();
    return pt;
  } // Pt

    /// this method returns the Ht: scalar sum of Pt
    double TopTopologicalVariables::Ht() const {
      // initialize
      double ht = 0.0;
      for (unsigned int i=0; i<_myobjects.size(); i++)
        ht += _myobjects[i].Pt();
      return ht;
    } // Ht

  /// this method returns the H: sum of the (scalar) energy, E
  double TopTopologicalVariables::H() const {
    // initialize
    double H = 0.0;
    for (unsigned int i=0; i<_myobjects.size(); i++)
      H+=_myobjects[i].E();
    return H;
  } // H

    /// this method returns the Transverse Mass, Mt(a, b)  
  /// = sqrt( sum{pT^2} - sum{px^2} - sum{px^2});
  double TopTopologicalVariables::TransverseMass() const {
    // Require at least two objects 
    if(_myobjects.size()<2) return -1.0; 
    // initialize
    double Mt = 0.0;
    double sum_pT = 0.0;
    double sum_px = 0.0;
    double sum_py = 0.0;
    // loop over _myobjects
    for (unsigned int i = 0; i < _myobjects.size(); i++) {
      sum_pT += _myobjects[i].Pt();
      sum_px += _myobjects[i].Px();
      sum_py += _myobjects[i].Py();
    } 
    
    Mt = sum_pT*sum_pT - sum_px*sum_px - sum_py*sum_py;
    
    if ( Mt >= 0.0 )
      Mt = sqrt(Mt); 
    else {
      if(fabs(Mt)<1.0e-6) 
        Mt = 0.0;  // Square of Mt is negative, but small
      else {
        std::cout << "In TopTopologicalVariables\n  Error: Square of transverse mass is negative!\n";
        return -1.0;
      } // Square of Mt is negative and large
    } // Square of Mt is negative
    
    return Mt;
    
  } // TransverseMass()
      
    /// this method returns the invariant mass: sqrt(E^2 - ThreeVector{P}^2 )
    double TopTopologicalVariables::M() const {
      // initialize
      TMBLorentzVector Sum(0,0,0,0);
      for (unsigned int i=0; i<_myobjects.size(); i++)
        Sum += _myobjects[i];
    
      if ( Sum.E()*Sum.E() - Sum.Mag32() < 0 ) {
        std::cout << "Error: Square of Invariant_mass is negative!" << std::endl;
        return -1.0;
      }
      return Sum.M();
    } // M()

  double TopTopologicalVariables::Mt() const {      
    double mt = -1.0;

    if(_myobjects.size() == 2 ) {
      double pt1 = _myobjects[0].Pt();
      double pt2 = _myobjects[1].Pt();
      double phi1 = _myobjects[0].Phi();
      double phi2 = _myobjects[1].Phi();

      double delta_phi = kinem::delta_phi(phi1, phi2);

      mt = TMath::Sqrt(2*pt1*pt2*(1-TMath::Cos(delta_phi)));
    }

    return( mt );
  } // Mt()

    /// this method returns the geometric mean of the Pt()s of the objects
    double TopTopologicalVariables::GeometricMeanPt() const {
      if(_myobjects.size() != 2) return -1.0;
      double meansq = 0.0;
      meansq = _myobjects[0].Pt() * _myobjects[1].Pt();
      if (meansq < 0) {
        std::cout << "Error: Square of GeometricMeanPt is negative!" << std::endl;
        return -1.0;
      }
      return sqrt(meansq);
    } // GeometricMeanPt()
  
  /// this method returns the weighted RMS of eta of the objects
  double TopTopologicalVariables::WeightedEtaRMS() const {
    double meaneta = 0.0;
    for (unsigned i=0; i<_myobjects.size(); i++) {
      meaneta += _myobjects[i].Pt() * _myobjects[i].Eta();
    }
    meaneta /= Ht();

    double weightsum = 0.0;
    double weightedEtaRMS = 0.0;
    for (unsigned i=0; i<_myobjects.size(); i++) {
      double sigmabg = 1.10 - 1.99e-3*_myobjects[i].Pt()  + 15.8/_myobjects[i].Pt();
      double sigmatt = 0.68 - 1.05e-3*_myobjects[i].Pt()  + 13.6/_myobjects[i].Pt();
      double weight = (sigmatt - sigmabg) / sigmatt;  // Usually this is negative.
      weightsum += weight;
      double deta = _myobjects[i].Eta()-meaneta;
      weightedEtaRMS += weight * deta * deta;
    }
    weightedEtaRMS /= weightsum;

    return weightedEtaRMS;
  } // WeightedEtaRMS()

    /// return minimum pair invariant mass
    double TopTopologicalVariables::MinimumPairMass() const
    {
      double minMass = 100000;
      Iterator end = _myobjects.end();
      for(Iterator itr1 = _myobjects.begin(); itr1 != end; ++itr1)
        {
          Iterator itr2 = itr1 + 1;
          for( ; itr2 != end; ++itr2)
            {
              double tempmass = (*itr1 + *itr2).M();
              if(tempmass < minMass)
                {
                  minMass = tempmass;
                }
            }
        }
      return minMass;
    } // MinimumPairMass()

  /// Note: this is KtMin, **NOT** KtMinPrime
    /// User needs to divide by his/her favorite energy variable to get
    /// KtMinPrime
    double TopTopologicalVariables::KtMin() const
    {
      double ptmin = 0;
      double drmin2 = 100000;
      Iterator end = _myobjects.end();
      for(Iterator itr1 = _myobjects.begin(); itr1 != end; ++itr1)
        {
          Iterator itr2 = itr1 + 1;
          for( ; itr2 != end; ++itr2)
            {
              double dp = kinem::delta_phi(itr1->Phi(), itr2->Phi());
              double de = fabs(itr1->Eta() - itr2->Eta());
              double dr2 = (dp*dp + de*de);
              if (dr2 < drmin2)
                {
                  drmin2 = dr2;
                  ptmin = itr1->Pt() < itr2->Pt() ? itr1->Pt() : itr2->Pt();
                }
            }
        }

      return sqrt(drmin2)*ptmin;
    } // KtMin()
  
  // Note: this is the minimal relative momentum between two objects, 
  //       a clean variation of KtMin. As for Ktmin, the user is
  //       encouraged to divide by an energy and get PtrelMinPrime
  double TopTopologicalVariables::PtrelMin() const
  {
    double ptmin = 999;
    for (unsigned int i=0; i<_myobjects.size(); ++i) {

      for (unsigned int j=0; j<_myobjects.size(); ++j) {

        if (i==j) continue;

        double cur = _myobjects[i].Pt(_myobjects[j]);
        if (cur < ptmin) ptmin = cur;
      }
    }
    return ptmin;
  }
    

  /// Calculate cos(theta*), where theta* is the angle between the
  /// 1st object and the beam axis in the objects' rest frame
  double TopTopologicalVariables::CosThetaStar() const
  {
    if (_myobjects.size() <= 0) return -1;
      
    // initialize
    TMBLorentzVector Sum(0,0,0,0);
    for (unsigned int i=0; i<_myobjects.size(); i++)
      Sum += _myobjects[i];
      
    TVector3 sumBoost (Sum.BoostVector());
      
    TMBLorentzVector lv (_myobjects[0]);
    lv.Boost (-1 * sumBoost); // boost to Sum's rest frame
      
    return cos (lv.Theta());
  } // CosThetaStar()

  /// Calculate cos(theta*), assuming no boost in x and y, where theta*
  /// is the angle between the
  /// 1st object and the beam axis in the objects' rest frame
  double TopTopologicalVariables::CosThetaStarJustZ() const
  {
    if (_myobjects.size() <= 0) return -1;
      
    // initialize
    TMBLorentzVector Sum(0,0,0,0);
    for (unsigned int i=0; i<_myobjects.size(); i++)
      Sum += _myobjects[i];
      
    TVector3 sumBoost (Sum.BoostVector());
      
    TMBLorentzVector lv (_myobjects[0]);
    lv.Boost (0, 0, -sumBoost.Z()); // boost to Sum's rest frame, assuming no x&y boost
      
    return cos (lv.Theta());
  } // CosThetaStar()


  /// Returns the lowest pT of all the objects
  double TopTopologicalVariables::SoftestPt() const
  {
    double minPt = 9999;
    for (unsigned int i=0; i<_myobjects.size(); ++i) {
      if (_myobjects[i].Pt() < minPt) minPt = _myobjects[i].Pt();
    }
    return minPt;
  }

  /// Find the lowest delta R between two objects
  double TopTopologicalVariables::MinDR() const
  {
    double minDr = 9999;
    for (unsigned int i=0; i<_myobjects.size()-1; ++i) {
      for (unsigned int j=i+1; j<_myobjects.size(); ++j) {
        double curDr = _myobjects[i].DeltaR (_myobjects[j]);
        if (curDr < minDr) minDr = curDr;
      }
    }
    return minDr;
  }

  /// Find the largest delta R between two objects
  double TopTopologicalVariables::MaxDR() const
  {
    double maxDr = -1;
    for (unsigned int i=0; i<_myobjects.size()-1; ++i) {
      for (unsigned int j=i+1; j<_myobjects.size(); ++j) {
        double curDr = _myobjects[i].DeltaR (_myobjects[j]);
        if (curDr > maxDr) maxDr = curDr;
      }
    }
    return maxDr;
  }

} // namespace top_cafe 


Revision:	1.1
Committed:	Tue Nov 11 23:01:22 2008 UTC (16 years, 5 months ago) by kukartse
Content type:	text/plain
Branch:	MAIN
CVS Tags:	V00-03-01, ZMorph_BASE_20100408, gak040610_morphing, V00-02-02, gak011410, gak010310, ejterm2010_25nov2009, V00-02-01, V00-02-00, gak112409, CMSSW_22X_branch_base, segala101609, V00-01-15, V00-01-14, V00-01-13, V00-01-12, V00-01-11, V00-01-10, gak031009, gak030509, gak022309, gak021209, gak040209, gak012809, V00-01-09, V00-01-08, V00-01-07, V00-01-06, V00-01-05, V00-01-04, V00-00-07, V00-00-06, V00-00-05, V00-00-04, V00-01-03, V00-00-02, V00-00-01, HEAD
Branch point for:	ZMorph-V00-03-01, CMSSW_22X_branch
Log Message:	initial creation of a package for LJMET multivariate analysis
#	User	Rev	Content
1	kukartse	1.1	/* File TopTopologicalVariables.hpp
2			*
3			* Created : Tue Oct 4 16:58:28 CDT 2005
4			* Author : Supriya JAIN, sjain@fnal.gov
5			*
6			* Purpose : Class containing generic methods for
7			* calculating topological variables
8			* (Instantiate using as argument: vector<TMBLorentzVector>
9			* of all objects for which you want any variable
10			* to be calculated)
11			*
12			* Last modified : Amnon Harel, 06 Mar 2006 (const correctness, cache eigenvalues)
13			* Comments :
14			*/
15
16			#include "TMatrixD.h"
17			#include "TVectorD.h"
18			#include "TRandom.h"
19
20			#include "LJMet/MultivariateAnalysis/interface/TopTopologicalVariables.h"
21			#include "LJMet/MultivariateAnalysis/interface/AnglesUtil.h"
22			#include <iostream>
23
24			using namespace std;
25
26			namespace top_cafe {
27
28			/// Copy constructor
29			TopTopologicalVariables::TopTopologicalVariables(const TopTopologicalVariables& that)
30			: _myobjects(that._myobjects)
31			{
32			_pv = 0;
33			if (that._pv) _pv = new TVectorD (*(that._pv));
34			}
35
36			/// Asignment operator
37			TopTopologicalVariables& TopTopologicalVariables::operator= (const TopTopologicalVariables& that)
38			{
39			if (this != &that)
40			{
41			_myobjects = that._myobjects;
42			if (_pv) delete _pv;
43			if (that._pv) _pv = new TVectorD (*(that._pv));
44			}
45			return *this;
46			}
47
48			/// destructor
49			TopTopologicalVariables::~TopTopologicalVariables()
50			{
51			if (_pv) delete _pv;
52			}
53
54			/// this method returns the centrality of a group of objects
55			double TopTopologicalVariables::Centrality() const
56			{
57			return Ht()/H();
58			} // Centrality()
59
60			/// this method returns the aplanarity of a group of objects
61			double TopTopologicalVariables::Aplanarity() const
62			{
63			ensurePV();
64			return 1.5 * (_pv)[2]; // alternative syntax: ... _pv->operator[](2)
65			} // Aplanarity()
66
67			/// this method returns the sphericity of a group of objects
68			double TopTopologicalVariables::Sphericity() const
69			{
70			ensurePV();
71			return 1.5 * ( (_pv)[2] + (_pv)[1] );
72			} // Sphericity()
73
74			/// this method returns the C of a group of objects (~momentum elipsiod surface area)
75			double TopTopologicalVariables::C() const
76			{
77			ensurePV();
78			return 3 * ( (_pv)[0] (_pv)[1] + (_pv)[0] * (_pv)[2] + (_pv)[1] * (*_pv)[2]);
79			} // C()
80
81			/// this method returns the D of a group of objects (~momentum elipsiod volume)
82			double TopTopologicalVariables::D() const
83			{
84			ensurePV();
85			return 27 * (_pv)[0] (_pv)[1] (*_pv)[2];
86			} // D()
87
88			/// this method returns the Momentum Tensor Eigenvalues of a group of objects
89			TVectorD TopTopologicalVariables::GetMomentumTensorEigenvalues() const
90			{
91			ensurePV ();
92			return *_pv; // implicit copy and cast
93			} // GetMomentumTensorEigenvalues()
94
95			/// internal method to prepare and cache the eigenvectors
96			void TopTopologicalVariables::ensurePV () const
97			{
98			if (_pv) return;
99
100			TMatrixD MomentumTensor(3,3);
101
102			Double_t p2_sum=0.0;
103
104			// for each _myobjects:
105			for ( unsigned int k=0; k<_myobjects.size(); k++ ) {
106
107			for ( int i=0; i<3; i++ ) // px, py, pz
108			for ( int j=0; j<3; j++ ) // px, py, pz
109			MomentumTensor(i,j) += _myobjects[k][i]*_myobjects[k][j];
110
111			// add the 3-momentum squared to the sum
112			p2_sum += _myobjects[k].Mag32();
113			} // Loop over _myobjects
114
115			// Divide the sums with the p2 sum
116			if ( p2_sum != 0 )
117			for ( int i=0; i<3; i++ ) // px, py, pz
118			for ( int j=0; j<3; j++ ) // px, py, pz
119			MomentumTensor(i,j) /= p2_sum;
120
121			_pv = new TVectorD (3);
122			MomentumTensor.EigenVectors(*_pv);
123			}
124
125			/// this method returns the Pt of a group of objects
126			double TopTopologicalVariables::Pt() const {
127			// initialize
128			TMBLorentzVector Sum(0.0,0.0,0.0,0.0);
129			double pt;
130			for (unsigned int i=0; i<_myobjects.size(); i++)
131			Sum += _myobjects[i];
132			pt = Sum.Pt();
133			return pt;
134			} // Pt
135
136			/// this method returns the Ht: scalar sum of Pt
137			double TopTopologicalVariables::Ht() const {
138			// initialize
139			double ht = 0.0;
140			for (unsigned int i=0; i<_myobjects.size(); i++)
141			ht += _myobjects[i].Pt();
142			return ht;
143			} // Ht
144
145			/// this method returns the H: sum of the (scalar) energy, E
146			double TopTopologicalVariables::H() const {
147			// initialize
148			double H = 0.0;
149			for (unsigned int i=0; i<_myobjects.size(); i++)
150			H+=_myobjects[i].E();
151			return H;
152			} // H
153
154			/// this method returns the Transverse Mass, Mt(a, b)
155			/// = sqrt( sum{pT^2} - sum{px^2} - sum{px^2});
156			double TopTopologicalVariables::TransverseMass() const {
157			// Require at least two objects
158			if(_myobjects.size()<2) return -1.0;
159			// initialize
160			double Mt = 0.0;
161			double sum_pT = 0.0;
162			double sum_px = 0.0;
163			double sum_py = 0.0;
164			// loop over _myobjects
165			for (unsigned int i = 0; i < _myobjects.size(); i++) {
166			sum_pT += _myobjects[i].Pt();
167			sum_px += _myobjects[i].Px();
168			sum_py += _myobjects[i].Py();
169			}
170
171			Mt = sum_pTsum_pT - sum_pxsum_px - sum_py*sum_py;
172
173			if ( Mt >= 0.0 )
174			Mt = sqrt(Mt);
175			else {
176			if(fabs(Mt)<1.0e-6)
177			Mt = 0.0; // Square of Mt is negative, but small
178			else {
179			std::cout << "In TopTopologicalVariables\n Error: Square of transverse mass is negative!\n";
180			return -1.0;
181			} // Square of Mt is negative and large
182			} // Square of Mt is negative
183
184			return Mt;
185
186			} // TransverseMass()
187
188			/// this method returns the invariant mass: sqrt(E^2 - ThreeVector{P}^2 )
189			double TopTopologicalVariables::M() const {
190			// initialize
191			TMBLorentzVector Sum(0,0,0,0);
192			for (unsigned int i=0; i<_myobjects.size(); i++)
193			Sum += _myobjects[i];
194
195			if ( Sum.E()*Sum.E() - Sum.Mag32() < 0 ) {
196			std::cout << "Error: Square of Invariant_mass is negative!" << std::endl;
197			return -1.0;
198			}
199			return Sum.M();
200			} // M()
201
202			double TopTopologicalVariables::Mt() const {
203			double mt = -1.0;
204
205			if(_myobjects.size() == 2 ) {
206			double pt1 = _myobjects[0].Pt();
207			double pt2 = _myobjects[1].Pt();
208			double phi1 = _myobjects[0].Phi();
209			double phi2 = _myobjects[1].Phi();
210
211			double delta_phi = kinem::delta_phi(phi1, phi2);
212
213			mt = TMath::Sqrt(2pt1pt2*(1-TMath::Cos(delta_phi)));
214			}
215
216			return( mt );
217			} // Mt()
218
219			/// this method returns the geometric mean of the Pt()s of the objects
220			double TopTopologicalVariables::GeometricMeanPt() const {
221			if(_myobjects.size() != 2) return -1.0;
222			double meansq = 0.0;
223			meansq = _myobjects[0].Pt() * _myobjects[1].Pt();
224			if (meansq < 0) {
225			std::cout << "Error: Square of GeometricMeanPt is negative!" << std::endl;
226			return -1.0;
227			}
228			return sqrt(meansq);
229			} // GeometricMeanPt()
230
231			/// this method returns the weighted RMS of eta of the objects
232			double TopTopologicalVariables::WeightedEtaRMS() const {
233			double meaneta = 0.0;
234			for (unsigned i=0; i<_myobjects.size(); i++) {
235			meaneta += _myobjects[i].Pt() * _myobjects[i].Eta();
236			}
237			meaneta /= Ht();
238
239			double weightsum = 0.0;
240			double weightedEtaRMS = 0.0;
241			for (unsigned i=0; i<_myobjects.size(); i++) {
242			double sigmabg = 1.10 - 1.99e-3*_myobjects[i].Pt() + 15.8/_myobjects[i].Pt();
243			double sigmatt = 0.68 - 1.05e-3*_myobjects[i].Pt() + 13.6/_myobjects[i].Pt();
244			double weight = (sigmatt - sigmabg) / sigmatt; // Usually this is negative.
245			weightsum += weight;
246			double deta = _myobjects[i].Eta()-meaneta;
247			weightedEtaRMS += weight * deta * deta;
248			}
249			weightedEtaRMS /= weightsum;
250
251			return weightedEtaRMS;
252			} // WeightedEtaRMS()
253
254			/// return minimum pair invariant mass
255			double TopTopologicalVariables::MinimumPairMass() const
256			{
257			double minMass = 100000;
258			Iterator end = _myobjects.end();
259			for(Iterator itr1 = _myobjects.begin(); itr1 != end; ++itr1)
260			{
261			Iterator itr2 = itr1 + 1;
262			for( ; itr2 != end; ++itr2)
263			{
264			double tempmass = (itr1 + itr2).M();
265			if(tempmass < minMass)
266			{
267			minMass = tempmass;
268			}
269			}
270			}
271			return minMass;
272			} // MinimumPairMass()
273
274			/// Note: this is KtMin, NOT KtMinPrime
275			/// User needs to divide by his/her favorite energy variable to get
276			/// KtMinPrime
277			double TopTopologicalVariables::KtMin() const
278			{
279			double ptmin = 0;
280			double drmin2 = 100000;
281			Iterator end = _myobjects.end();
282			for(Iterator itr1 = _myobjects.begin(); itr1 != end; ++itr1)
283			{
284			Iterator itr2 = itr1 + 1;
285			for( ; itr2 != end; ++itr2)
286			{
287			double dp = kinem::delta_phi(itr1->Phi(), itr2->Phi());
288			double de = fabs(itr1->Eta() - itr2->Eta());
289			double dr2 = (dpdp + dede);
290			if (dr2 < drmin2)
291			{
292			drmin2 = dr2;
293			ptmin = itr1->Pt() < itr2->Pt() ? itr1->Pt() : itr2->Pt();
294			}
295			}
296			}
297
298			return sqrt(drmin2)*ptmin;
299			} // KtMin()
300
301			// Note: this is the minimal relative momentum between two objects,
302			// a clean variation of KtMin. As for Ktmin, the user is
303			// encouraged to divide by an energy and get PtrelMinPrime
304			double TopTopologicalVariables::PtrelMin() const
305			{
306			double ptmin = 999;
307			for (unsigned int i=0; i<_myobjects.size(); ++i) {
308
309			for (unsigned int j=0; j<_myobjects.size(); ++j) {
310
311			if (i==j) continue;
312
313			double cur = _myobjects[i].Pt(_myobjects[j]);
314			if (cur < ptmin) ptmin = cur;
315			}
316			}
317			return ptmin;
318			}
319
320
321			/// Calculate cos(theta), where theta is the angle between the
322			/// 1st object and the beam axis in the objects' rest frame
323			double TopTopologicalVariables::CosThetaStar() const
324			{
325			if (_myobjects.size() <= 0) return -1;
326
327			// initialize
328			TMBLorentzVector Sum(0,0,0,0);
329			for (unsigned int i=0; i<_myobjects.size(); i++)
330			Sum += _myobjects[i];
331
332			TVector3 sumBoost (Sum.BoostVector());
333
334			TMBLorentzVector lv (_myobjects[0]);
335			lv.Boost (-1 * sumBoost); // boost to Sum's rest frame
336
337			return cos (lv.Theta());
338			} // CosThetaStar()
339
340			/// Calculate cos(theta), assuming no boost in x and y, where theta
341			/// is the angle between the
342			/// 1st object and the beam axis in the objects' rest frame
343			double TopTopologicalVariables::CosThetaStarJustZ() const
344			{
345			if (_myobjects.size() <= 0) return -1;
346
347			// initialize
348			TMBLorentzVector Sum(0,0,0,0);
349			for (unsigned int i=0; i<_myobjects.size(); i++)
350			Sum += _myobjects[i];
351
352			TVector3 sumBoost (Sum.BoostVector());
353
354			TMBLorentzVector lv (_myobjects[0]);
355			lv.Boost (0, 0, -sumBoost.Z()); // boost to Sum's rest frame, assuming no x&y boost
356
357			return cos (lv.Theta());
358			} // CosThetaStar()
359
360
361			/// Returns the lowest pT of all the objects
362			double TopTopologicalVariables::SoftestPt() const
363			{
364			double minPt = 9999;
365			for (unsigned int i=0; i<_myobjects.size(); ++i) {
366			if (_myobjects[i].Pt() < minPt) minPt = _myobjects[i].Pt();
367			}
368			return minPt;
369			}
370
371			/// Find the lowest delta R between two objects
372			double TopTopologicalVariables::MinDR() const
373			{
374			double minDr = 9999;
375			for (unsigned int i=0; i<_myobjects.size()-1; ++i) {
376			for (unsigned int j=i+1; j<_myobjects.size(); ++j) {
377			double curDr = _myobjects[i].DeltaR (_myobjects[j]);
378			if (curDr < minDr) minDr = curDr;
379			}
380			}
381			return minDr;
382			}
383
384			/// Find the largest delta R between two objects
385			double TopTopologicalVariables::MaxDR() const
386			{
387			double maxDr = -1;
388			for (unsigned int i=0; i<_myobjects.size()-1; ++i) {
389			for (unsigned int j=i+1; j<_myobjects.size(); ++j) {
390			double curDr = _myobjects[i].DeltaR (_myobjects[j]);
391			if (curDr > maxDr) maxDr = curDr;
392			}
393			}
394			return maxDr;
395			}
396
397			} // namespace top_cafe
398
399
400