flyopunta/src/src_dev/devotcl.cpp

/***************************************************************************
                          devotcl.cpp  -  description
                             -------------------
    begin                : Sun Dec 7 2003
    copyright            : (C) 2003 by gmp
    email                : peppe@newpeppe
  ***************************************************************************
  *                                                                         *
  *   NA48  simulation program.                                             *
  *                                                                         *
  ***************************************************************************/
#include "parm.h"
#include "devotcl.h"
#include "devrtcl.h"

DevOtCl::~DevOtCl()
{
}


//================== DevOtCl ======Erede di DevOt
DevOtCl::DevOtCl()
{
  devtype=TypDevOtCl;
  devclass="Ottg_Cili";
  rin = Lin[1]=Lin[0];
  Lin.Norma();
  rinq=rin*rin;
}

//----------------- --------
void DevOtCl::Prgeom()
{

  Gout<<"==================  R i v e l a t o r e =================";
  Gout<<"\n --> < "<<nome <<" >  [ "<<fun
  <<" ]  Ottagonale con buco cilindrico  "<<std::endl;
  Device::Prgeom();
  Gout<<std::endl;
}

//======================================================
//----------------- -----P o s i z i o n e ---
int DevOtCl::Posizione()
{
  if ( DevOt::Posizione() >0 ) return 1;
// puo' essere interno o nel buco
  P_Cin=X_dev-Centrin;
  Dind=P_Cin-Lin;
  Dins=- ( Lin+P_Cin );
// se sono interno al rivelatore verifico se siamo nel buco
  if ( ( drinq= P_Cin[0]*P_Cin[0]+P_Cin[1]*P_Cin[1] ) >rinq ) return 0;
//  sono forse nel buco ...controllo meglio i limiti
// in zeta. Ricordo Lin[2] e' la zeta della base destra...
  if ( P_Cin[2]< -Lin[2] || P_Cin[2]> Lin[2] ) return 0; // e' esterno
// ritorno -1  perche' nel buco
  return -1;
}
//----------------------C a m E s t e r --------------------
// usa quella di DevOt  per eredita'
//===========================================================
//----------------------C a m I n t e r --------------------

double DevOtCl::CamInter()
{
//uso quella di DevRt per calcolare il cammino
// interno senza il buco...
  double tratto=DevOt::CamInter();
//poi verifico se colpisce le pareti del buco .....
  double a,b,c,dltq,t,nv,z;
// calcolo il cammino positivo verso le  superfici del cilindro interno.
// e scelgo il cammino maggiore
//  calcolo i coefficienti della equazione intersezione di una retta con cili
  camm=-1.;
// controllo se colpisce la base
  if ( ( nv=V_dev[2] ) >0.0 && Dins[2]>0.0 )  // piani z
    {t=Dins[2]/nv;  if ( t>camm ) camm=t;}
  else if ( nv<0.0 && Dind[2]>0.0 )
    {t=-Dind[2]/nv; if ( t>camm ) camm=t;}
//
  a=V_dev[0]*V_dev[0]+V_dev[1]*V_dev[1];
  if ( a>0.0 )  //revised 31.07.01 gmp
  {
    b=V_dev[0]*P_Cin[0]+V_dev[1]*P_Cin[1];
    if ( b<0. )
    {
      c=drinq-rinq;
      if ( ( dltq=b*b-a*c ) >=0.0 )
      {
        t= - ( b + sqrt ( dltq ) ) /a; // Scelgo la soluzione minore
        if ( t>=0.0 )
        {
          // check limiti in zeta
          z=P_Cin[2]+nv*t;
          if ( -Lin[2]<=z&&z<=Lin[2] ) camm=t;
        }
      }
    }
  }
  // quindi confronto camm con tratto e accetto il cammino piu' breve......
  if ( camm<0.0 ) {camm=tratto; return camm;}
  if ( camm>0.0 && tratto<camm ) camm=tratto;
  return camm;
}
//===============================================================
//--------------- C a m B u c o ---------------------------------
double DevOtCl::CamBuco()
{

  double a,b,c,dltq,t,nv;
  camm=1.0e+10;
// calcolo il cammino verso le  superfici
// quella cilindrica e dalle 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.
  a=V_dev[0]*V_dev[0]+V_dev[1]*V_dev[1];
  if ( a>0.0 )
  {
    b=V_dev[0]*P_Cin[0]+V_dev[1]*P_Cin[1];
    c=drinq-rinq;           // nota sempre negativo internamente al cilindro.
    if ( ( dltq=b*b-a*c ) >=0.0 )
    {
      t= ( -b + sqrt ( dltq ) ) /a; // Scelgo la soluzione positiva
      if ( t>=0.0 ) camm=t;
    }
  }
// controllo anche la zeta....
  if ( ( nv=V_dev[2] ) >0.0 )  // piani z
    {t=-Dind[2]/nv; if ( t<camm ) camm=t; }
  else if ( nv<0.0 )
    {t=Dins[2]/nv;  if ( t<camm ) camm=t; }
  return camm;

}