flyopunta/src/src_uti/gvet.cpp

/***************************************************************************
                          GPmrg->h  -  description
                             -------------------
    revised              : Febbraio agosto 2009 <===========
    copyright            : (C) 2000 by Giuseppe M Pierazzini
    email                : pierazzini@unipi.it
 ***************************************************************************/
//----------------------   vector -----------------------------
#define GVECTOR
#include "parm.h"
#include "gvet.h"
#include "gmatrix.h"
using namespace std;
//------------------------------------------------------
//--------------------------------------------------------
// note precision bits: double 53 bits, float 24 bits


double normaux;
// costruttori
// vedi gevet.h
//======


//====================================================
// Definisco gli operatori dell'algebra tra  vettori
// attenzione mettere le parentesi per le priorita'
// si usano i registri per permettere operazioni
// complesse sulla stessa linea e con migliore utilizzazione
// della memoria.
//----------------------------------------------------

// operatore unario di cambio segno
gvet  gvet::operator-()
{
  gvet temp;
  temp.x=-x;
  temp.y=-y;
  temp.z=-z;
  temp.norma=norma;
  temp.normaq=normaq;
  return temp;
}

gvet  gvet::operator* ( const double a )
{
  gvet temp;
  temp.x=x*a;
  temp.y=y*a;
  temp.z=z*a;
  temp.norma=norma*fabs ( a );
  temp.normaq=temp.norma*temp.norma;
  return temp;
}

// ritorna il vettore normalizzato sul piano xz
gvet   gvet::NxzVerso()
{
  gvet temp;
  normaux=sqrt ( x*x+z*z );
  temp.x=x/normaux;
  temp.y=0. ;
  temp.z=z/normaux;
  temp.norma=temp.normaq=1.0;
  return temp;
}
// ritorna il vettore normalizzato sul piano yz
gvet  gvet::NyzVerso()
{
  gvet temp;
  normaux=sqrt ( y*y+z*z );
  temp.x=0.;
  temp.y=y/normaux;
  temp.z=z/normaux;
  temp.norma=temp.normaq=1.0;
  return temp;
}
// somma di due vetori V1+V2 or dif
gvet  gvet::operator+ ( const gvet& a )
{
  gvet temp;
  temp.x=x+a.x;
  temp.y=y+a.y;
  temp.z=z+a.z;
  temp.Norma();
  return temp;
}
gvet   gvet::operator- ( const gvet& a )
{
  gvet temp;
  temp.x=x-a.x;
  temp.y=y-a.y;
  temp.z=z-a.z;
  temp.Norma();
  return temp;
}
// prodotto vettoriale  V1 esterno V2
gvet gvet::operator& ( const gvet& b )
{
  gvet temp;
  temp.x=y*b.z-z*b.y;
  temp.y=z*b.x-x*b.z;
  temp.z=x*b.y-y*b.x;
  temp.Norma();
  return temp;
}


gmatrix&  gvet::operator >> (gmatrix & M )
{
//crea il tensore  |this><this|  
  M.Proj(*this);
  return  M;
}

gvet & gvet::operator = (gmatrix & M )

{ 
  // restituisce in this il vettore  colonna 0 della matrice M.
  M.Get_Vcol (*this);
 return *this;
} 


//----------------   output  ---------------
// output per debug...............
void   gvet::print ( const char *t )
{
  Gout.setf ( ios::right,ios::adjustfield );
  Gout<<std::setfill ( ' ' );
  Gout<<std::setprecision ( 3 ) <<fixed;
  Gout<<"\n"<<setw ( 7 ) <<t<<"=["
  <<setw ( 9 ) <<x<<" "
  <<setw ( 9 ) <<y<<" "
  <<setw ( 9 ) <<z<<" ]" <<" |"
  <<setw ( 9 ) <<norma<<"|";
  Gout.setf ( ios::floatfield );

}

void   gvet::print ( long ev,const char *pnome,int pid,const char *devnome )
{
  Gout.setf ( ios::right,ios::adjustfield );
  Gout<<std::setfill ( ' ' );
  Gout<<std::setprecision ( 2 ) <<fixed;
  Gout<<"\n"<<setw ( 7 ) <<ev
  <<"  "<<setw ( 5 ) <<pnome<< " "<<setw ( 3 ) <<pid
  <<"  "<<setw ( 5 ) <<devnome <<"=["
  <<setw ( 9 ) <<x<<" "
  <<setw ( 9 ) <<y<<" "
  <<setw ( 11 ) <<z<<" ]| "
  <<setw ( 9 ) <<norma<<" | ";
  Gout.setf ( ios::floatfield );
}




//------------------------------------------------------------
//============================================================
//------------------------------------------------------------

qvet::qvet ( double massa )
{
  x=0.0;
  y=0.0;
  z=0.0;
  m=massa;
  mq=m*m;
  Energia();
  Norma();
}
// attenzione un po'pericoloso....
qvet::qvet ( double* p )  //da arrey  px,py,pz,massa
{
  x=p[0];
  y=p[1];
  z=p[2];
  m=p[3];
  mq=m*m;
  Energia();
  Norma();
}
qvet::qvet ( gvet P,double massa )  //da vettore  P e massa
{
  P.putv( v );
  m=massa;
  mq=m*m;
  Energia();
  Norma();
}


qvet::qvet ( double px,double py, double pz, double massa )
{
  x=px;
  y=py;
  z=pz;
  m=massa;
  mq=m*m;
  Energia();
  Norma();
}


//====================================================
// Definisco gli operatori dell'algebra tra  vettori.

// double & qvet::operator[](int i) {return v[i];}


// si introducono registri per eseguire senza errori equazioni con  piu
// termini del tipo V*c contemporanei

// operatore unario; cambio solo la parte spaziale....!!!!
qvet  qvet::operator-()
{
  qvet qaux;
  qaux.x=-x;
  qaux.y=-y;
  qaux.z=-z;
  qaux.e=e;
  qaux.norma=norma;
  qaux.normaq=normaq;
  qaux.m =m;
  qaux.mq=mq;
  return qaux;;
}

qvet  qvet::operator* ( const double a )
{
  qvet qaux;
  qaux.x=x*a;
  qaux.y=y*a;
  qaux.z=z*a;
  qaux.e=e*a;
  qaux.norma=norma*a;
  qaux.normaq=normaq*a*a;
  qaux.m =m*a;
  qaux.mq=mq*a*a;
  return qaux;;
}


// ricalcola l'invariante relativistico
double qvet::Invar()
{
  normaq=x*x+y*y+z*z; norma=sqrt ( normaq );
  mq=e*e - normaq;
  if ( mq>0.0 ) m=sqrt ( mq ); else m=-sqrt ( -mq ); //Ritorna la massa negativa!!!!
  return m;
}

double qvet::Invarq()
{
  normaq=x*x+y*y+z*z; norma=sqrt ( normaq );
  mq=e*e - normaq;
  if ( mq>0.0 ) m=sqrt ( mq ); else m=-sqrt ( -mq );//Ritorna la massa negativa!!!!
  return mq;
}



// ricalcola l'energia ...a partire dall'impulso e massa
double qvet::Energia()
{
  mq=m*m;
  normaq=x*x+y*y+z*z; norma=sqrt ( normaq );
  e=sqrt ( mq+normaq );
  return e;
}
//


void qvet::operator= ( double *a )    // attenzione come (px,py,pz,m)
{
  x=a[0];y=a[1];z=a[2];
  m=a[3];
  Energia();
}

qvet & qvet::operator+= ( const qvet& a )              // attenzione non calcola norma..
{
  x+=a.x;y+=a.y;z+=a.z;
  e+=a.e;   return *this;
}

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


// ricalcola l'impulso  ...a partire dall'Energia e massa
// usa come versore di P quello di prima. Attenzione il versore non deve essere nullo.
double qvet::Impulso()
{
  double puls=sqrt ( ( e-m ) * ( e+m ) );
  double rnor=puls/sqrt ( x*x+y*y+z*z );
  x= x*rnor;
  y= y*rnor;
  z= z*rnor;
  norma=puls;
  normaq=puls*puls;
  return norma;
}


//

qvet  qvet::operator+ ( const qvet& a )
{
  qvet qaux;
  qaux.x=x+a.x;
  qaux.y=y+a.y;
  qaux.z=z+a.z;
  qaux.e=e+a.e;
  return qaux;;
}

qvet   qvet::operator- ( const qvet& a )
{
  qvet qaux;
  qaux.x=x-a.x;
  qaux.y=y-a.y;
  qaux.z=z-a.z;
  qaux.e=e-a.e;
  return qaux;;
}


//------------------------------------------------------

qvet& qvet::Lburst ( qvet & Burst )
{
  static gvet Beta;
  static double gam;
  // Burst e' un quadrivettore (P,E) relativo alla particella centro di massa.
  //se m e' nulla si usano i vecchi valori...

  if ( Burst.m>0.0 )
  {
    Beta.setvn ( Burst[0]/Burst.e, Burst[1]/Burst.e, Burst[2]/Burst.e );
    gam=Burst.e/Burst.m;
  }
//    Pnew=(gvet)*this + Beta*(gam*(pb/(1.0+1.0/gam) + e));  // funge perfettamente
//    gvet PL=Beta*(Beta%this)/Beta.normaq;
//    gvet Pnew=(gvet)*this + PL*(gam-1.) +Beta*(gam*e);

  gvet Pnew= ( gvet ) *this + Beta* ( Beta%this ) * ( ( gam-1. ) /Beta.normaq ) +Beta* ( gam*e );
  setvn ( Pnew[0],Pnew[1],Pnew[2],m );
  return *this;
}


// output per debug...............
void   qvet::print ( const char *t )
{

 Gout.setf ( ios::right,ios::adjustfield );
// Gout.setf ( ios::floatfield );
  
  Gout<<setfill ( ' ' ) <<fixed<<setprecision(2)
  <<"\n"<<setw ( 7 ) <<t<<"=["
  <<setw ( 8 ) <<x<<" "
  <<setw ( 8 ) <<y<<" "
  <<setw ( 9 ) <<z<<" "
  <<setw ( 8 ) <<e<<" ] { n="
  <<setw ( 8 ) <<norma<<" "
  <<setw ( 8 ) <<normaq<<" m="
  <<setw ( 8 )<<setprecision(3) <<m<<" "
  <<setw ( 8 ) <<mq<<"}"<<flush;
 
}

void   qvet::print ( long ev,const char *pnome,int pid,const char *devnome )
{
  Gout.setf ( ios::right,ios::adjustfield );
  Gout<<std::setfill ( ' ' );
  Gout<<std::setprecision ( 2 ) <<fixed;
  Gout<<"\n"<<setw ( 7 ) <<ev
  <<"  "<<setw ( 5 ) <<pnome<< " "<<setw ( 3 ) <<pid
  <<"  "<<setw ( 5 ) <<devnome <<"=["
  <<setw ( 8 ) <<x<<" "
  <<setw ( 8 ) <<y<<" "
  <<setw ( 11 ) <<z<<" "
  <<setw ( 8 ) <<e<<" ]"     <<" |"
  <<setw ( 8 ) <<norma<<"| "
  <<setw ( 8 ) <<normaq<<"| "
  <<setw ( 8 ) <<m<<"| "
  <<setw ( 8 ) <<m<<flush;
  Gout.setf ( ios::floatfield );
}