siscone/spherical/momentum.h

// -*- C++ -*-
// File: momentum.h                                                          //
// Description: header file for 4-momentum class Cmomentum                   //
// This file is part of the SISCone project.                                 //
// WARNING: this is not the main SISCone trunk but                           //
//          an adaptation to spherical coordinates                           //
// For more details, see http://projects.hepforge.org/siscone                //
//                                                                           //
// Copyright (c) 2006-2008 Gavin Salam and Gregory Soyez                     //
//                                                                           //
// This program is free software; you can redistribute it and/or modify      //
// it under the terms of the GNU General Public License as published by      //
// the Free Software Foundation; either version 2 of the License, or         //
// (at your option) any later version.                                       //
//                                                                           //
// This program is distributed in the hope that it will be useful,           //
// but WITHOUT ANY WARRANTY; without even the implied warranty of            //
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             //
// GNU General Public License for more details.                              //
//                                                                           //
// You should have received a copy of the GNU General Public License         //
// along with this program; if not, write to the Free Software               //
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
//                                                                           //
// $Revision:: 256                                                          $//
// $Date:: 2008-07-14 13:52:16 +0200 (Mon, 14 Jul 2008)                     $//

#ifndef __SPH_VECTOR_H__
#define __SPH_VECTOR_H__

#include <vector>
#include <math.h>
#include <siscone/reference.h>
#include "geom_2d.h"
#include <siscone/defines.h>

namespace siscone_spherical{

class CSph3vector{
 public:
  CSph3vector();

  CSph3vector(double _px, double _py, double _pz);

  ~CSph3vector();

  CSph3vector& operator = (const CSph3vector &v);

  const CSph3vector operator + (const CSph3vector &v);

  const CSph3vector operator - (const CSph3vector &v);

  const CSph3vector operator / (const double &r);

  CSph3vector& operator += (const CSph3vector &v);

  CSph3vector& operator -= (const CSph3vector &v);

  CSph3vector& operator *= (const double &r);

  CSph3vector& operator /= (const double &r);

  inline double perp() const {return sqrt(perp2());}

  inline double perp2() const {return px*px+py*py;}

  inline double norm() const {return sqrt(px*px+py*py+pz*pz);}

  inline double norm2() const {return px*px+py*py+pz*pz;}

  inline double phi() const {return atan2(py, px);}

  inline double theta() const {return atan2(perp(),pz);}

  void build_norm();

  void build_thetaphi();

  void get_angular_directions(CSph3vector &angular_dir1, CSph3vector &angular_dir2);

  double px;        
  double py;        
  double pz;        

  double _norm;     
  double _theta;    
  double _phi;      

  // the following part is used for checksums //
  siscone::Creference ref;   
};

class CSphmomentum : public CSph3vector{
 public:
  CSphmomentum();

  CSphmomentum(CSph3vector &init, double E=0.0);

  CSphmomentum(double _px, double _py, double _pz, double _E);

  //CSphmomentum(double _eta, double _phi, siscone::Creference _ref);

  ~CSphmomentum();

  inline double mass() const {return sqrt(mass2());}

  inline double mass2() const {return perpmass2()-perp2();}

  inline double perpmass() const {return sqrt((E-pz)*(E+pz));}

  inline double perpmass2() const {return (E-pz)*(E+pz);}

  inline double Et() const {return E/sqrt(1.0+pz*pz/perp2());}

  inline double Et2() const {return E*E/(1.0+pz*pz/perp2());}

  CSphmomentum& operator = (const CSphmomentum &v);

  const CSphmomentum operator + (const CSphmomentum &v);

  CSphmomentum& operator += (const CSphmomentum &v);

  CSphmomentum& operator -= (const CSphmomentum &v);

  double E;         

  int parent_index; 
  int index;        
};

bool operator < (const CSphmomentum &v1, const CSphmomentum &v2);

bool momentum_theta_less(const CSphmomentum &v1, const CSphmomentum &v2);

bool momentum_pt_less(const CSphmomentum &v1, const CSphmomentum &v2);


// some handy utilities //

inline double sqr(double x){return x*x;}

inline double dot_product3(const CSph3vector &v1, const CSph3vector &v2){
  //double tmp = v1.px*v2.px + v1.py*v2.py + v1.pz*v2.pz;
  //if (!isfinite(tmp)){
  //  std::cout << "dot_product inf: " << std::endl;
  //  std::cout << "  angles: " << v1._theta << " " << v1._phi << " and " << v2._theta << " " << v2._phi << std::endl; 
  //  std::cout << "  moms  : " << v1.px << " " << v1.py << " " << v1.pz 
  //          << " and "      << v2.px << " " << v2.py << " " << v2.pz << std::endl;
  //}
  return v1.px*v2.px + v1.py*v2.py + v1.pz*v2.pz;
}

inline CSph3vector cross_product3(const CSph3vector &v1, const CSph3vector &v2){
  //CSph3vector tmp;
  //tmp.px = v1.py*v2.pz-v1.pz*v2.py;
  //tmp.py = v1.pz*v2.px-v1.px*v2.pz;
  //tmp.pz = v1.px*v2.py-v1.py*v2.px;
  //return tmp;
  return CSph3vector(v1.py*v2.pz-v1.pz*v2.py,
                  v1.pz*v2.px-v1.px*v2.pz,
                  v1.px*v2.py-v1.py*v2.px);
}

inline double norm2_cross_product3(const CSph3vector &v1, const CSph3vector &v2){
  return sqr(v1.py*v2.pz-v1.pz*v2.py) + sqr(v1.pz*v2.px-v1.px*v2.pz) + sqr(v1.px*v2.py-v1.py*v2.px);
}

inline double get_tan2_distance(const CSphmomentum &v1, const CSphmomentum &v2){
  return norm2_cross_product3(v1,v2)/sqr(dot_product3(v1,v2));
}

inline double get_distance(const CSph3vector *v1, const CSph3vector *v2){
  return atan2(sqrt(norm2_cross_product3(*v1,*v2)), dot_product3(*v1,*v2));
}

inline bool is_closer(const CSph3vector *v1, const CSph3vector *v2, const double tan2R){
  double dot = dot_product3(*v1,*v2);
  return (dot>=0) && (norm2_cross_product3(*v1,*v2)<=tan2R*dot*dot);
}

inline bool is_closer_safer(const CSph3vector *v1, const CSph3vector *v2, const double cosR){
  return dot_product3(*v1,*v2)>=cosR*sqrt(v1->norm2()*v2->norm2());
  //double dot = dot_product3(*v1,*v2);
  //return (dot>=0) && (norm2_cross_product3(*v1,*v2)<tan2R*dot*dot);
}

inline CSph3vector operator * (const double &r, const CSph3vector &v){
  CSph3vector tmp = v;
  return tmp*=r;
}
}
#endif