JETAN/fastjet/siscone/momentum.h

   1 // -*- C++ -*-
   2 ///////////////////////////////////////////////////////////////////////////////
   3 // File: momentum.h                                                          //
   4 // Description: header file for 4-momentum class Cmomentum                   //
   5 // This file is part of the SISCone project.                                 //
   6 // For more details, see http://projects.hepforge.org/siscone                //
   7 //                                                                           //
   8 // Copyright (c) 2006 Gavin Salam and Gregory Soyez                          //
   9 //                                                                           //
  10 // This program is free software; you can redistribute it and/or modify      //
  11 // it under the terms of the GNU General Public License as published by      //
  12 // the Free Software Foundation; either version 2 of the License, or         //
  13 // (at your option) any later version.                                       //
  14 //                                                                           //
  15 // This program is distributed in the hope that it will be useful,           //
  16 // but WITHOUT ANY WARRANTY; without even the implied warranty of            //
  17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             //
  18 // GNU General Public License for more details.                              //
  19 //                                                                           //
  20 // You should have received a copy of the GNU General Public License         //
  21 // along with this program; if not, write to the Free Software               //
  22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
  23 //                                                                           //
  24 // $Revision:: 163                                                          $//
  25 // $Date:: 2007-04-26 22:31:02 +0200 (Thu, 26 Apr 2007)                     $//
  26 ///////////////////////////////////////////////////////////////////////////////
  27
  28 #ifndef __VECTOR_H__
  29 #define __VECTOR_H__
  30
  31 #include <vector>
  32 #include <math.h>
  33 #include "reference.h"
  34 #include "geom_2d.h"
  35 #include "defines.h"
  36
  37 namespace siscone{
  38
  39 /**
  40  * \class Cmomentum
  41  * \brief base class for dynamic coordinates management
  42  *
  43  * This class contains the information for particle or group of
  44  * particles management.
  45  * It includes all Lorentz properties as well as tools for summing them.
  46  * Note: 'sums' over phi angles are indeed averages. This allows to
  47  *       deal with periodicity at each step
  48  */
  49 class Cmomentum{
  50  public:
  51   /// default ctor
  52   Cmomentum();
  53
  54   /// ctor with initialisation
  55   Cmomentum(double _px, double _py, double _pz, double _E);
  56
  57   /// ctor with detailed initialisation
  58   Cmomentum(double _eta, double _phi, Creference _ref);
  59
  60   /// default dtor
  61   ~Cmomentum();
  62
  63   /// computes pT
  64   inline double perp() const {return sqrt(perp2());}
  65
  66   /// computes pT^2
  67   inline double perp2() const {return px*px+py*py;}
  68
  69   /// computes m
  70   inline double mass() const {return sqrt(mass2());}
  71
  72   /// computes m^2
  73   inline double mass2() const {return perpmass2()-perp2();}
  74
  75   /// transverse mass, mt = sqrt(pt^2+m^2) = sqrt(E^2 - pz^2)
  76   inline double perpmass() const {return sqrt((E-pz)*(E+pz));}
  77
  78   /// transverse mass squared, mt^2 = pt^2+m^2 = E^2 - pz^2
  79   inline double perpmass2() const {return (E-pz)*(E+pz);}
  80
  81   /// computes transverse energy
  82   inline double Et() const {return E/sqrt(1.0+pz*pz/perp2());}
  83
  84   /// computes transverse energy (squared)
  85   inline double Et2() const {return E*E/(1.0+pz*pz/perp2());}
  86
  87   /// assignment of vectors
  88   Cmomentum& operator = (const Cmomentum &v);
  89
  90   /// addition of vectors
  91   /// !!! WARNING !!! no updating of eta and phi !!!
  92   const Cmomentum operator + (const Cmomentum &v);
  93
  94   /// incrementation of vectors
  95   /// !!! WARNING !!! no updating of eta and phi !!!
  96   Cmomentum& operator += (const Cmomentum &v);
  97
  98   /// decrementation of vectors
  99   /// !!! WARNING !!! no updating of eta and phi !!!
 100   Cmomentum& operator -= (const Cmomentum &v);
 101
 102   /// build eta-phi from 4-momentum info
 103   /// !!!                WARNING                   !!!
 104   /// !!! computing eta and phi is time-consuming  !!!
 105   /// !!! use this whenever you need eta or phi    !!!
 106   /// !!! automatically called for single-particle !!!
 107   void build_etaphi();
 108
 109   double px;        ///< x-momentum
 110   double py;        ///< y-momentum
 111   double pz;        ///< z-momentum
 112   double E;         ///< energy
 113
 114   double eta;       ///< particle pseudo-rapidity
 115   double phi;       ///< particle azimuthal angle
 116   int parent_index; ///< particle number in the parent list
 117   int index;        ///< internal particle number
 118
 119   //////////////////////////////////////////////
 120   // the following part is used for checksums //
 121   //////////////////////////////////////////////
 122   Creference ref;   ///< reference number for the vector
 123 };
 124
 125 /// ordering of two vectors
 126 /// this is by default done w.r.t. their references
 127 bool operator < (const Cmomentum &v1, const Cmomentum &v2);
 128
 129 /// ordering of vectors in eta (e.g. used in collinear tests)
 130 bool momentum_eta_less(const Cmomentum &v1, const Cmomentum &v2);
 131
 132 /// ordering of vectors in pt
 133 bool momentum_pt_less(const Cmomentum &v1, const Cmomentum &v2);
 134
 135
 136 //////////////////////////
 137 // some handy utilities //
 138 //////////////////////////
 139
 140 /// get distance between to eta-phi points
 141 /// \param eta  eta coordinate of first point
 142 /// \param phi  phi coordinate of first point
 143 /// \param v    vector defining the second point
 144 inline double get_distance(double eta, double phi, Cmomentum *v){
 145   double dx, dy;
 146
 147   dx = eta - v->eta;
 148   dy = fabs(phi - v->phi);
 149   if (dy>M_PI)
 150     dy -= twopi;
 151
 152   return dx*dx+dy*dy;
 153 }
 154
 155 }
 156
 157 #endif