flyopunta/src/src_dev/devcln.cpp

/***************************************************************************
                          devcln.cpp  -  description
                             -------------------
    begin                : Sat May 26 2001
    copyright            : (C) 2001 by Giuseppe Pierazzini
    email                : peppe@unipi.it
  ***************************************************************************
  *                                                                         *
  *   NA48  simulation program.                                             *
  *                                                                         *
  ***************************************************************************/
#include "parm.h"
#include "devcln.h"
#include "evento.h"

//using namespace std;

DevCln::~DevCln()
{
}
//==================  Cilindro=======Erede di RettEdro
DevCln::DevCln()
{
  devtype=TypDevCln;
  devclass="Cilindro";
  rout=Lout.y=Lout.x;  // simmetrizzo le coordinate x,y Lout.x sempre differente da 0 e = rout;
  routq=rout*rout;
  Lout.Norma();
  Rqsize=Lout.normaq;

}

//----------------- --------
void DevCln::Prgeom()
{
  Gout<<"==================  R i v e l a t o r e =================";
  Gout<<"\n --> < "<<nome <<" >  [ "<<fun
  <<" ]  Cilindro di rout "<< rout<<"  senza buco. ";
  Device::Prgeom();

}

//======================================================
//----------------- -----P o s i z i o n e ---
int DevCln::Posizione()
{
  
//  Doutd=X_dev-Lout;     Doutd.Norma();
//  Douts= - ( Lout+X_dev );  Douts.Norma();

  Doutd.z=X_dev.z-Lout.z;   //solo la componente z
  // Distanze dai piani dispari in Douts (per differenza).
  Douts.z=  - ( Lout.z+X_dev.z );   // l'operatore unario "-" funge....

  if ( X_dev.normaq>Rqsize ) return 1;

  drq= X_dev.x*X_dev.x+ X_dev.y*X_dev.y;
  vxyq=V_dev.x*V_dev.x+V_dev.y*V_dev.y;   // per futuro uso
  if ( Douts.z>0.0 || Doutd.z>0.0 || drq>routq ) return 1; //e' fouri
  return 0;    // zero se interno, 1 se esterno.


}
//----------------------------------------------------------
//----------------------C a m E s t e r --------------------
double DevCln::CamEster()
{
  double b,c,dltq,t,nv;
// calcolo il cammino positivo verso la  superfice del cilindro.
// e scelgo il cammino maggiore
// calcolo i coefficienti della equazione intersezione di una retta con cil
  camm=-1.;

// se mi allontano  ritorno -1.


  if ( ( nv=V_dev.z ) >0.0&&Doutd.z>0.0 )return -1.;



// *****   controllo se colpisce una delle  basi

  if ( nv >0.0 && Douts.z>0.0 )  // piani z
    {t=Douts.z/nv;  if ( t>camm ) camm=t;}
  else if ( nv<0.0 && Doutd.z>0.0 )
    {t=-Doutd.z/nv; if ( t>camm ) camm=t;}

// non ha colpito la base ...
//   ora verifico se colpisce il lato nei limiti di zeta...
//  revised  31.07.01 gmp


  if ( vxyq<=0.0 ) return camm; //per vxyq==0 corre lungo zeta /vedi Posizione()
  double z=0.;
  b=V_dev.x*X_dev.x+V_dev.y*X_dev.y;
  if ( b>0. )  return camm;
  c=drq-routq;
  if ( ( dltq=b*b-vxyq*c ) <0.0 ) return camm;
  t= - ( b + sqrt ( dltq ) ) /vxyq; // Scelgo la soluzione minore
  if ( t>=0.0 )
  {
    // check limiti in zeta
    z=X_dev.z+nv*t;
    if ( -Lout.z<=z&&z<=Lout.z ) camm=t;
  }
//***********
//    if ( pr_dev->Get_Idm()==14 ) SpiaPr ( " Ester2 ", 49,49,camm );
  return camm;
}
//===========================================================
//----------------------C a m I n t e r --------------------

double DevCln::CamInter()
{
  double b,c,dltq,t,nv;
  camm=1.0e+10;
// calcolo il cammino verso le  superfici esterne
// quella cilindrica e le basi del cilindro.
// scelgo il cammino positivo minore...
// calcolo il cammino positivo verso la superfice del cilindro interno.
// le soluzioni sono sempre di segno opposto.


  if ( vxyq>0.0 )
  {
    b=V_dev.x*X_dev.x+V_dev.y*X_dev.y;
    c=drq-routq;           // nota sempre negativo internamente al cilindro.
    if ( ( dltq=b*b-vxyq*c ) >=0.0 )
    {
      t= ( -b + sqrt ( dltq ) ) /vxyq; // Scelgo la soluzione positiva
      if ( t>=0.0 ) camm=t;    
    }
  }
// controllo anche la zeta....
  if ( ( nv=V_dev.z ) >0.0 )  // piani z
    {t=-Doutd.z/nv; if ( t<camm ) camm=t; }
  else if ( nv<0.0 )
    {t=Douts.z/nv;  if ( t<camm ) camm=t; }
  return camm;
}