flyopunta/src/src_uti/gvet.h

/***************************************************************************
                        GVET.H
                             -------------------
    revised              : Febbraio settembre 2010    
    copyright            : (C) 2000 by Giuseppe M Pierazzini
    email                : pierazzini@unipi.it
 ***************************************************************************/
#ifndef GVET_H
#define GVET_H
#include <cmath>

//****************************************************
//
//-------------------   Vettori tridimensionali --------------
//----------------------------------------------------
// algebra:
// Operator/Method                          Description
//---- COSTRUCTOR---
//   gvet()         :   empty vector
//   gvet(double a,double b,double c)  : vector= (a,b,c)

//---OPERATOR ------                          -----------
//   operator !    :   to calculate inversion of matrix.
//   operator +=   :   combined subtraction and assignment operator
//   operator -=   :   combined  subctraction  and assignment operator
//   operator -    :   Unary change of sign
//   gvet & operator=(gvet *a) : permette V= *vec
//   void operator=(double *a) : permette V=a;  con a[] = {x,y,z}
//   gvet & operator+(gvet &)  :  V= V1 +   V2
//   gvet & operator+(gvet* a) :  V= V1 +  *V2
//   gvet & operator-(gvet &)  :  V= V1 -  V2
//   gvet & operator-(gvet* a) :  V= V1 -  *V2
//   gvet &  operator&(gvet &) :  Prodotto vettoriale
//   gvet & operator&(gvet *b) :  Prodotto vettoriale
//   gmat & operator&(gmat &)  :  Prodotto vettoriale generalizzato ad una matrice M=VxM1
//   gvet & operator %( gmat &):  Prodotto vettore x matrice in vettore  V= V1.M
//   double operator%(gvet & A):  Prodotto scalare
//   double operator%(gvet * A):  Prodotto scalare
//   gvet & operator*=(double c)  : operatore *= per una costante
//   gvet & operator/=(double c)  : operatore /= per una costante

//----- FUNCTIONS ---
//   gvet &  Verso()   : calcola il versore...
//   gvet & NVerso()   : calcola prima la norma.. just in the case!! +
//   double Norma()    : calcola la Norma
//   double XYNorma()  : calcola la norma solo sul piano xy
//   double YZNorma()  : calcola la norma solo sul piano yz
//   double ZXNorma()  : calcola la norma solo sul piano zx
//   gvet & NxzVerso() : Versore sul piano xz
//   gvet & NyzVerso() : Versore sul piano yz
//   gvet & setvn(double a,double b,double c) : define V e normalizza
//   gvet & setvn(double *a)                  : define V e normalizza
//   gvet & setv(double a,double b,double c)  : define V e NON normalizza
//   gvet & setv(double *a)                   : define V e NON normalizzarriva solo con if(Nome== 
  
//   virtual void Reset()                  : azzera
//   void putv(double &x,double &y,double &z)  : define x,y,z
//   void putv(float &x,float &y,float &z)     : define x,y,z
//   void putvn(double &a,double &b,double &c,double &d)
//   void putvn(float &a,float &b,float &c,float &d)
//   virtual void print(char *)                       : stampa info
//-------------------------------------------------------
class gmatrix;
class gvet
{   
 
protected:
    union
    {
        struct {
            double x, y, z,t;
        };  // t è per combalità con i quadrivettori
        struct {
            double px, py, pz,e;
        };
        double v[4];// componenti del vettore
    };
public:
    union
    {
        struct {
            double norma,normaq;
        };
        struct {
            double mom,momq;
        };
    };
    // norma ... attenzione va aggiornato se cambia v.

    gvet() {
        x=y=z=0.0;
        normaq=norma=0.0;
    };
    gvet ( double a,double b, double c )
    {
        x= a;
        y= b;
        z= c;
        normaq=x*x+y*y+z*z;
        norma=sqrt ( normaq );
    };
// operatori delicati possono generare grane...meglio usare setvn,setv
//   gvet(double * p){x=p[0];y=p[1];z=p[2];normaq=x*x+y*y+z*z;norma=sqrt(normaq);};
//   gvet(double p){x=y=z=p; normaq=3.*p*p; norma=sqrt(normaq);};
//-------------------
    virtual ~gvet() {};

    double& operator[] ( int i ) {
        return v[i];
    };
    double& operator() ( int i ) {
        return v[i];
    };
// operatore unario di cambio segno
    gvet  operator-();

    double Norma()
    {
        normaq=x*x+y*y+z*z;
        norma=sqrt ( normaq );
        return norma;
    };
    double XYNorma() {
        return sqrt ( x*x+y*y );
    };
    double YZNorma() {
        return sqrt ( y*y+z*z );
    };
    double ZXNorma() {
        return sqrt ( z*z+x*x );
    };
    gvet  Verso() {
        return *this/norma;
    };       // calcola il versore...
    gvet NVerso() {
        Norma();    // calcola la norma.. just in the case!!
        return *this/norma;
    }

    gvet  NxzVerso();      // Versore sul piano xz
    gvet  NyzVerso();       // Versore sul piano yz


//     gvet & operator = ( gvet *a )
//     {
//         *this=*a;
//         return *this;
//     }

    gvet & operator = ( gmatrix &m );




    gvet  operator+ ( const gvet & );
    gvet  operator+ ( gvet* a ) {
        return *this+*a;
    }
    gvet  operator- ( const gvet& a );
    gvet  operator- ( gvet* a ) {
        return *this-*a;
    }


    // this+=V; this-=V;
    gvet & operator+= ( const gvet& a )              // attenzione non calcola norma..
    {
        x+=a.x;
        y+=a.y;
        z+=a.z;
        return *this;
    }
    gvet & operator+= ( gvet * a )              // attenzione non calcola norma..
    {
        *this+=*a;
        return *this;
    }

    gvet & operator-= ( const gvet& a ) // attenzione non calcola norma..
    {
        x-=a.x;
        y-=a.y;
        z-=a.z;
        return *this;
    }
    gvet & operator-= ( gvet * a )   // attenzione non calcola norma..
    {
        *this-=*a;
        return *this;
    }


    

//Prodotto vettoriale P=this.x.A
    gvet    operator & ( const gvet & );
    gvet   operator & ( gvet *b ) {
        return *this & *b;
    }
// prodotto esterno per una matrice



// prodotto scalare
    double operator% ( const gvet & A )
    {
        return ( x*A.x+y*A.y+z*A.z );
    };  // prodotto scalare
    double operator% ( gvet * A ) {
        return *this%*A;
    };
    friend double operator% ( gvet * A,gvet & B ) {
        return *A%B;
    } ;

// prodotto per un parametro this*a
    gvet  operator* ( const double );
//   friend  gvet  operator*(double c, gvet & A){ return (A*c);} // friend function
    friend  gvet  operator* ( double c, gvet  A ) {
        return ( A*c );
    }
// divide per un parametro
    gvet  operator/ ( const double a ) {
        double b=1./a;
        return   *this*b;
    }

// operatore *= per una costante
    gvet & operator*= ( double c )
    {
        x*=c;
        y*=c;
        z*=c;
        norma*=fabs ( c );
        normaq=norma*norma;
        return *this;
    }

// operatore /= per una costante
    gvet & operator/= ( double c )
    {
        x/=c;
        y/=c;
        z/=c;
        norma/=fabs ( c );
        normaq=norma*norma;
        return *this;
    }
    
    
    
gmatrix&  operator >> (gmatrix &);
    

// vettore  per una matrice
//   gvet &  operator %( gmat &);
//setting or putting
    virtual void Reset() {
        x=y=z=t=norma=normaq=0.;
    };
    gvet & setvn ( const double a,const double b,const double c ) {
        x=a;
        y=b;
        z=c;
        Norma();
        return *this;
    }
    gvet & setvn ( const float a,const float b,const float c ) {
        x=a;
        y=b;
        z=c;
        Norma();
        return *this;
    }


    gvet & setvn ( double *a ) {
        x=a[0];
        y=a[1];
        z=a[2];
        Norma();
        return *this;
    }
    gvet & setv ( double a,double b,double c ) {
        x=a;
        y=b;
        z=c;
        return *this;
    }
    
      gvet & setvt ( double a,double b,double c,double g) {
        x=a;
        y=b;
        z=c;
                                t=g;   //il tempo !!! 
        return *this;
    }  
    
    
    gvet & setv ( double *a ) {
        x=a[0];
        y=a[1];
        z=a[2];
        return *this;
    }
    
   void putv( double &a,double &b,double &c ) {
        a=x;
        b=y;
        c=z;
    }
    void putv( float &a,float &b,float &c ) {
        a=x;
        b=y;
        c=z;
    }
    void putv( double *a ) {
        a[0]=x;
        a[1]=y;
        a[2]=z;
    }  
    
     void putvn(double &a,double &b,double &c,double &d ) {
        a=x;
        b=y;
        c=z;
        d=norma;
    }
    void putvn( float &a,float &b,float &c,float &d ) {
        a=x;
        b=y;
        c=z;
        d=norma;
    }
    void putvn( double *a) {
        a[0]=x;
        a[1]=y;
        a[2]=z;
        a[3]=norma;
    }  
 
    virtual void print ( const char * );
    virtual void print ( long ev,const char *pnome,int pid,const char *devnome );
};


//****************************************************
//****************************************************
//=====================================================
// The 4 dimensions vector @author Giuseppe Pierazzini
//-----------------------------------------------------
//----------------------------------------------------
// algebra:
// Operator/Method                          Description
//   constructor
// qvet(double );
// qvet(double*);
// qvet(gvet &,double);
// qvet(double,double,double,double);

//      operatori
// qvet & operator=(qvet *a)        ;
// void   operator=(double *a)      ;attenzione come (px,py,pz,m)
// qvet & operator+=(qvet& a)       ;
// qvet & operator+=(qvet * a)      ;
// qvet & operator-=(qvet& a)       ;
// qvet & operator-=(qvet * a)      ;
// qvet & operator-()               ;operatore unario: solo space
// qvet & operator+(qvet &)         ;
// qvet & operator+(qvet* a)        ;
// qvet & operator-(qvet &)         ;
// qvet & operator-(qvet* a)        ;
// double operator*(qvet & A)       ;prodotto scalare tra quadrivettori
// double operator*(qvet * A)       ;
// friend double operator*(qvet * A,qvet & B)
// qvet & operator*(double )        ;prodotto per un parametro
// friend qvet & operator*(double a,qvet & B)
// qvet & operator/(double a)       ;divisione per un parametro

// double Invar();
// double Invarq();
// double Energia();  //Calcola l'energia a partire dalla massa e impulso
// double Impulso();  //Calcola l'impulso  a partire dalla massa e energia

//   qvet & setvn(double a,double b,double c, double d)
//   qvet & setvn(double *a)
//   qvet & setvn(gvet &Vg,double a)
//   qvet & setv(double a,double b,double c,double mass)
//   qvet & setv(double *a)
//   qvet & setv(gvet &Vg,double a)

// set del quadrivettore ; il quarto valore e' l'energia;
//   qvet & setqn(double a,double b,double c, double d)
//   qvet & setqn(double *a)
//   qvet & setqn(gvet &Vg,double a)

//   void Reset()
//   void putq(double &a,double &b,double &c,double &d)

//   void putq(double *a)

//   qvet & Lburst(qvet &);  //lorentz buster....in direzione di Burst...
//   void print(char *);
//-------------------------------------------------------------------------------

class qvet:public gvet
{
    friend class gvet;
protected:
//  double e = componente [3] del vettore  V
public: // norma ... attenzione va aggiornato se cambia v.
    double m,mq;    // invariante relativistico...
    qvet() {
        x=y=z=e=m=mq=normaq=norma=0.0;
    };
    qvet ( double );
    qvet ( double* );
    qvet ( gvet ,double );
    qvet ( double,double,double,double );
    ~qvet() {};
    double Invar();
    double Invarq();
    double Energia();  //Calcola l'energia a partire dalla massa e impulso
    double Impulso();  //Calcola l'impulso  a partire dalla massa e energia

    qvet & operator= ( qvet *a )    {
        *this=*a;
        return *a;
    }
    void   operator= ( double *a );   // attenzione come (px,py,pz,m)
    qvet & operator+= ( const qvet& a );             // attenzione non calcola norma..
    qvet & operator+= ( qvet * a ) {
        *this+=*a;
        return *this;
    }
    qvet & operator-= ( const qvet& a );             // attenzione non calcola norma..
    qvet & operator-= ( qvet * a ) {
        *this-=*a;
        return *this;
    }
    qvet   operator-();                      // operatore unario: solo space

    qvet   operator+ ( const  qvet & );
    qvet   operator+ ( qvet* a ) {
        return *this+*a;
    };
    qvet   operator- ( const qvet & );
    qvet   operator- ( qvet* a ) {
        return *this-*a;
    }

// prodotto vettoriale... eredita!
// prodotto scalare   quadrivettoriale

    double operator* ( const qvet & A ) {
        return ( e*A.e- ( x*A.x+y*A.y+z*A.z ) );
    }
    double operator* ( qvet * A ) {
        return *this**A;    // idem
    }
    friend double operator* ( qvet * A,qvet & B ) {
        return *A*B;
    }


    qvet   operator* ( double );    // prodotto per un parametro
//   friend qvet  operator*(double a,qvet & B){return B*a;};
    friend qvet  operator* ( double a,qvet  B ) {
        return B*a;
    };

    qvet  operator/ ( double a ) {
        a=1./a;
        return   *this*a;
    }


// set del vettore parte spaziale; il quarto valore e' la massa
// calcola la norma  e l'energia

    qvet & setvn ( const double a,const double b,const double c, const double massa )
    {
        x=a;
        y=b;
        z=c;
        m=massa;
        mq=m*m;
        normaq=x*x+y*y+z*z;
        norma=sqrt(normaq);
        e=sqrt(normaq+mq);
        return *this;
    }
    qvet & setvn ( double *a )
    {
        x=a[0];
        y=a[1];
        z=a[2];
        m=a[3];
        mq=m*m;
        normaq=x*x+y*y+z*z;
        norma=sqrt(normaq);
        e=sqrt(normaq+mq);     
        return *this;
    }
    qvet & setvn ( gvet Vg,const double massa )
    {
        x=Vg[0];
        y=Vg[1];
        z=Vg[2];
        m=massa;
        mq=m*m;
        normaq=x*x+y*y+z*z;
        norma=sqrt(normaq);
        e=sqrt(normaq+mq);  
        return *this;
    }

// set del vettore parte spaziale; il quarto valore e' la massa; non normalizza

    qvet & setv ( double a,double b,double c,double massa )
    {
        x=a;
        y=b;
        z=c;
        m=massa;
        return *this;
    }
    qvet & setv ( double *a )
    {
        x=a[0];
        y=a[1];
        z=a[2];
        m=a[3];
        return *this;
    }
    qvet & setv ( gvet &Vg,double massa )
    {
        x=Vg[0];
        y=Vg[1];
        z=Vg[2];
        m=massa;
        return *this;
    }

// set del quadrivettore ; il quarto valore e' l'energia;
// calcolo della norma e dell'invariante
    qvet & setqn ( double a,double b,double c, double eg )
    {
        x=a;
        y=b;
        z=c;
        e=eg;
        Invar();
        return *this;
    }
    qvet & setqn ( double *a )
    {
        x=a[0];
        y=a[1];
        z=a[2];
        e=a[3];
        Invar();
        return *this;
    }
    qvet & setqn ( gvet Vg,double eg )
    {
        x=Vg[0];
        y=Vg[1];
        z=Vg[2];
        e=eg;
        Invar();
        return *this;
    }

    void Reset() {
        x=y=z=e=norma=normaq=m=mq=0.;
    };
    
    
    void putq( double &a,double &b,double &c,double &d ) {
        a=x;
        b=y;
        c=z;
        d=e;
    }
    void putq( float &a,float &b,float &c,float &d ) {
        a=x;
        b=y;
        c=z;
        d=e;
    }

    void putq( double *a ) {
        a[0]=x;
        a[1]=y;
        a[2]=z;
        a[3]=e;
    }
    
    
    
    
    void Riposo() {   // con massa già definita
        x=y=z=norma=normaq=0.;
        e=m;
        mq=m*m;
    };
    void Riposo(double massa) {// con massa in input 
        x=y=z=norma=normaq=0.;
        e=m=massa;
        mq=m*m;
    };
    qvet & Lburst ( qvet & Burst );  //lorentz buster....in direzione del Burst...
    void print ( const char * );
    void print ( long ev,const char *pnome,int pid,const char *devnome );
};

#endif