flyopunta/src/src_dev/devsector.cpp

/***************************************************************************
                          devsector.cpp  -  description
                             -------------------
    begin                : settembre 2008
    copyright            : (C) 2001 by Giuseppe Pierazzini
    email                : pierazzini@unipi.it
  ***************************************************************************
  *                                                                         *
  *   NA48  simulation program.                                             *
  *                                                                         *
  ***************************************************************************/
#include "parm.h"
#include "devsector.h"

//using namespace std;

DevSector::~DevSector()
{
}
//==================  Cilindro=======Erede di RettEdro
DevSector::DevSector()
{
  devtype=TypDevSector;
  devclass="Cili_Sector";
  rin =Lin.y=Lin.x;
  Lin.Norma();
  rinq=rin*rin;
}

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

}

//======================================================
//----------------- -----P o s i z i o n e ---
int DevSector::Posizione()
{

  if ( DevCln::Posizione() >0 ) return 1;  
  //  ritorn 1: il punto e' esterno al cilindro in cui è inscritto il settore
// il centro del cilindro non è detto che sia sull'asse z!

// verifico se il punto è interno al settore cilindrico.
// Il vertice del settore corrisponde al centro coordinate del device.
// L'asse x è mediano al settore. La semiapertura è
// OpenSector= x rad.
// verifico se il punto è interno al settore cilindrico
// che suppongo avere un raggio maggiore rout già verificato
// e un raggio minore rin che verifico ora
  if ( ( X_dev.x*X_dev.x+X_dev.y*X_dev.y ) <rinq ) return 1;
  double fi= atan2 ( X_dev.y,X_dev.x );
  if ( abs ( fi ) > OpenSector ) return 1;
  return 0;    // zero se interno, 1 se esterno.


}
//----------------------------------------------------------
//----------------------C a m E s t e r --------------------
double DevSector::CamEster()
{
  double nv,dxyq,b,c,dltq,fi,t=-1.;
  gvet Vn,Xpos;
  camm=-1.;

// se mi allontano  ritorno -1.
  if ( ( nv=V_dev.z ) >0.0&&Doutd.z>0.0 ) return -1.;

// *****   controllo se colpisce la base
  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;}
// colpisco nei limiti giusti
  Xpos= X_dev+V_dev*camm;
  dxyq=Xpos.x *Xpos.x +Xpos.y*Xpos.y;
  fi= atan2 ( X_dev.y,X_dev.x );
  if ( abs ( fi ) <= OpenSector && dxyq >=rinq && dxyq<=routq ) return camm;
  else camm=-1.0;
  if ( vxyq<=0.0 ) return camm; //per vxyq==0 corre lungo zeta e non colpisce nulla

  if ( ( c=drq-routq ) >0.0 )  // la particella è esterna e può colpire il cilindro lateralmente....
  {
    b=V_dev.x*X_dev.x+V_dev.y*X_dev.y;
    if ( b<0.0 )
    {
      if ( ( dltq=b*b-vxyq*c ) < 0.0 )  return camm; // non colpisce mai  il cilindro... quindi esco con -1.0
      t= - ( b + sqrt ( dltq ) ) /vxyq; // Scelgo la soluzione minore
      if ( t>0.0 )  // ma lo è per forza....verifico i limiti
      {
        Xpos= X_dev+V_dev*t;
        fi= atan2 ( X_dev.y,X_dev.x );
        if ( abs ( fi ) <= OpenSector && fabs ( Xpos.z ) <= Lout.z ) {camm=t; return camm;} // esco contento
      }
    }
  }

// la normale al piano è calcolata uscente dal settore.
// se urta il piano laterale verifico i limiti.
 
  Vn.x=-sin ( OpenSector );
  Vn.y= cos ( OpenSector );
  Vn.z=0.0;
  if ( ( b=Vn%V_dev ) <0. )  t= - ( Vn%X_dev ) /b;
  else
  {
    Vn.x=sin ( OpenSector );
    Vn.y=cos ( OpenSector );
    Vn.z=0.0;
    if ( ( b=Vn%V_dev ) <0. )  t= - ( Vn%X_dev ) /b;
  }
  if ( t>0.0 )  // verifico i limiti....in raggi e zeta 
  {
    Xpos= X_dev+V_dev*t;
    dxyq=Xpos.x *Xpos.x +Xpos.y*Xpos.y;
    if ( dxyq <=routq && dxyq >=rinq &&
            fabs ( Xpos.z ) <= Lout.z ) {camm=t; return camm;} //  ha colpito il lato ed esco contento..
  }

// e verifico i limiti in zeta e angolo  (note vxyq>0.)
// la posizione della particella può essere dentro o fuori il cerchio di raggio rin

  b=V_dev.x*X_dev.x+V_dev.y*X_dev.y;
  if( (c=drq-rinq) >0.0) return camm;  
// return camm=-1.0:  punto di partenza più distante di rin dall'asse ... non posso colpire l'interno 

  if ( ( dltq=b*b-vxyq*c ) <0.0 ) return camm;  // se il punto è interno è sempre positivo....
  t= - ( b - sqrt ( dltq ) ) /vxyq; // devo scegliere comunque la soluzione positiva 
  if ( t>0.0 )
  {
    Xpos= X_dev+V_dev*t;
    fi= atan2 ( X_dev.y,X_dev.x );
    if ( abs ( fi ) > OpenSector ) return camm;
    if ( fabs ( Xpos.z ) > Lout.z ) return camm;
    camm=t;
  }

  return camm;
}
//===========================================================
//----------------------C a m I n t e r --------------------

double DevSector::CamInter()
{
  double b,c,dltq,t=1.0e+11;

  gvet Vn,Xpos;
  camm=DevCln::CamInter();


  Vn.x=-sin ( OpenSector );
  Vn.y= cos ( OpenSector );
  Vn.z=0.0;

  if ( ( b=Vn%V_dev ) >0. ) t= - ( Vn%X_dev ) /b;

  else
  {
    Vn.x=sin ( OpenSector );
    Vn.y=cos ( OpenSector );
    if ( ( b=Vn%V_dev ) >0. ) t= - ( Vn%X_dev ) /b;

  }
  if ( t<camm )
  {
    Xpos= X_dev+V_dev*t;
    double dxyq=Xpos.x *Xpos.x +Xpos.y*Xpos.y;
    if ( dxyq <=routq && dxyq >=rinq &&
            fabs ( Xpos.z ) <= Lout.z ) camm=t; //  ha colpito un piano laterale
  }

  if ( vxyq>0.0 )
  {
    b=V_dev.x*X_dev.x+V_dev.y*X_dev.y;
    if ( b<0. )   // revised 31.07.01 gmp
    {
      c=drq-rinq;
      if ( ( dltq=b*b-vxyq*c ) >=0.0 )
      {
        t= - ( b + sqrt ( dltq ) ) /vxyq; // Scelgo la soluzione minore
        if ( t<camm )
        {
          Xpos= X_dev+V_dev*t;
          double fi= atan2 ( X_dev.y,X_dev.x );
          if ( abs ( fi ) > OpenSector ) return camm;
          if ( fabs ( Xpos.z ) > Lin.z ) return camm;
          camm=t;
        }
      }
    }
  }


  return camm;
}
//--------------------------------------------
gvet& DevSector::Lab2cDev ( gvet &Xin )
{
  // trasformazione di vettori al dev con rototraslazione...
// la procedura qui è invertita rispetto al metodo di base...


static gvet _Xdc,_Xdev;
  _Xdc.setvn (Xin%Nx,Xin%Ny,Xin%Nz );
 _Xdev=_Xdc-Centro;
  return _Xdev;

}
//--------------------------------------------
gvet& DevSector::Devc2Lab ( gvet &Xdev )
{
  // trasformazione di vettori dal Dev al lab con traslazione...
  // non calcola la norma...
/*
  static gvet _Xlab;
  _Xlab=Nx*Xdev[0] + Ny*Xdev[1] + Nz*Xdev[2];
  _Xlab=_Xlab+Centro;
  return  _Xlab;
*/

static gvet _Xlab, _Xaux;
_Xaux=Xdev+Centro;
  _Xlab=Nx*_Xaux[0] + Ny*_Xaux[1] + Nz*_Xaux[2];
  return  _Xlab;


}