flyopunta/src/src_dev/devsveto.cpp

/***************************************************************************
 *   Copyright (C) 2009 by giuseppe Pierazzini                                      *
 *   giuseppe@pierazzini.it                                                               *
 *                                                                         *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.                  *
 *                                                                         *
 *  Dipartimento di Fisica  E.Fermi / INFN Pisa Italy                      *
 ***************************************************************************/
#include "parm.h"
#include "evento.h"
#include "devsveto.h"
using namespace GmpConst;
using namespace std;


DevSveto::DevSveto()
        : DevCrn()
{

//     double latox,latoy,latoz, rmedia;
//   gvet *Ubl,*Use; //Versori degli assi dei blocchi e dei settori.
//   gvet *C_bl,*C_sec; // centri dei blocchi e dei settori
//   int Blocchi, iblocco,isector;
//   double Fasebase;

    devtype=TypDevSveto;
    devclass="Sveto_Crn";
    rin =Lin.y=Lin.x;
    Lin.Norma();
    rinq=rin*rin;

    if (nickname[0]=='A' && nickname[1]=='n') { //serve per controllo, ma è superfluo
        if (nickname[2]!='9' && nickname[3]=='\0') { //anti 1-8 ('\0' è il null character!)
            Strati=5; //numero di strati per veto
            if (nickname[2]=='6'||nickname[2]=='7'||nickname[2]=='8') //anti 6-8
                Blocchi=48;   // numero di blocchi per strato
            else //anti 1-5
                Blocchi=32;
        }
        else { //anti 9-12
            Strati=4; //numero di strati per veto
            if (nickname[3]=='2') // anti 12
                Blocchi=64;   // numero di blocchi per strato
            else // blocchi 9-11
                Blocchi=60;   // numero di blocchi per strato
        }
    }
                
    iblocco=-1;
    isector=-1;
    istrato=-1;

    //  Lout.z;  // = semispessore della corona in z
    latoy=latoz= Lout.z/Strati; // semispessori del leadglass ... (attenzione a Spectrometer09.app!)
    latox= rout-rin;   // lunghezza del leadglass=spessore della corona= rout-rin

//calcolo i versori degli assi  e i centri  dei blocchi
//il primo blocco ha l'asse orientato secondo la FaseBase.
//il primo settore è spostato di Dfi5 in senso antiorario se guardato da downstream

    int TBlocchi=Strati*Blocchi;

    Dfi = 2.* PiGreco/Blocchi;
    U_bl=new gvet[TBlocchi];
    U_sc=new gvet[TBlocchi];
    C_bl=new gvet[TBlocchi];
    C_sc=new gvet[TBlocchi];
    double rmed=(rout+rin)*.5;
    double Dfi5=Dfi*0.5;
    for (int i=0;i<TBlocchi;i++)
    {
        istrato=(i/Blocchi);
        FaseBase[istrato]=istrato*(Dfi/Strati);
        double fi = FaseBase[istrato]+i*Dfi;
        // se Fasebase=0. il primo blocco è lungo l'asse x
        double Centro_z_strato=2.*latoz*istrato+latoz -Lout.z;

//      Gout<<"\n C: "<< Centro_z_strato<<" " << i<< " "<<strato;
        U_bl[i].setvn(cos(fi),sin(fi),0);
        C_bl[i]= U_bl[i]*rmed;
        C_bl[i].z=Centro_z_strato;
        fi+=Dfi5;
        U_sc[i].setvn(cos(fi),sin(fi),0);
        C_sc[i]= U_sc[i]*rmed;
        C_sc[i].z=Centro_z_strato;
    }
    
    
    
    
}

DevSveto::~DevSveto()
{
}
//======================================================
//----------------- --------
void DevSveto::Prgeom()
{
    Gout<<"==================  R i v e l a t o r e =================";
    Gout<<"\n --> < "<<nome <<" >  [ "<<fun
    <<" ]  Contenitore a forma di corona cilindrica . ";
    Device::Prgeom();

}
int DevSveto::Get_Structinfo()
{
    int code=-1;
    if ( istrato ==-1) return -1;
    if ( iblocco>-1) code=iblocco;
    if ( isector>-1) code=isector+1000;
    return code;
}

//----------------- -----P o s i z i o n e ---
int DevSveto::Posizione()
{
    int pos=0;
    iblocco=-1;
    isector=-1;
    istrato=-1;
    if ( ( pos=DevCrn::Posizione())!=0 ) return pos; // esce se esterno o nel "buco" della corona.
//
// qui ... la particella  è interno alla corona
// in quale strato siamo
    istrato = int(0.5* (X_dev.z+Lout.z)/latoz); // 0<=istrato<5

// verifico se il punto particella  è interno anche ad un blocco leadglass
// Ricordo:
// L'asse centrale di ogni blocco, pensato come un  parallelepipedo,  punta al centro della corona.
// Il suo versore è individuato da Ubl[i] (calcolato prima) tipo Ubl[i]=(vx,vy,0).
// Se il punto è nella corona la sua posizione in zeta e in alto-basso ( secondo il versore Ubl) è ok.
// Dobbiamo calcolare la sua distanza dall'asse centrale sul piano  x,y
//  e verificare che è minore della senilarghezza, laty,  dal blocco.
    double fdf,fi;
    gvet Xpunto=X_dev;
    Xpunto.z=0;
    fi=fdf =atan2 ( Xpunto.y,Xpunto.x )-FaseBase[istrato];  // attenzione deve essere sempre >0!!
    if (fi<0.)fi+=2.*PiGreco;
    iblocco=int(fi/Dfi+0.5) %Blocchi;// è il blocco più vicino al punto... aggiungo lo strato
    iblocco+=istrato*Blocchi;

//  if(Eventi_Fatti<100) Gout<<"\n Pos1 "<<setw(8)<<pr_dev->Get_Ido()<<" "<<setw(8)<<X_dev.z<<" "
//       <<setw(8)<<fi<<" " <<setw(3)<<iblocco<<" "<<setw(3) <<isector<< " " <<setw(3)<<istrato ;

    double dist= (Xpunto-U_bl[iblocco]*(U_bl[iblocco]%Xpunto)).Norma();
    if ( dist< latoy) return 0;  // è dentro il blocco "iblocco" con isector=-1;


// non è nel blocco,ma è nel settore esterno al blocco
    iblocco=-1;
    fi=fdf-Dfi*0.5;
    if (fi<0.)fi+=2.*PiGreco;

    isector= int(fi/Dfi+0.5) %Blocchi + istrato*Blocchi;

//    if(Eventi_Fatti<100) Gout<<"\n Pos2 "<<setw(8)<<pr_dev->Get_Ido()<<" "<<setw(8)<<X_dev.z<<" "
//       <<setw(8)<<fi<<" " <<setw(3)<<iblocco<<" "<<setw(3) <<isector<< " " <<setw(3)<<istrato ;

    return -1;  // La particella è nel "buco", o meglio, nel settore  isector con iblocco=-1!

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

double DevSveto::CamInter()
{
    /*
     siamo in un blocco leadglass  iblocco.
    Procedura:
    1)Calcolo il cammino che dovrebbe fare la particella per uscire dallla corona e lo memorizzo in camm..
    2) Verifico il cammino per uscire lateralmente dal blocco.. se minore di camm ridefinisco camm.
      Per questo
     Calcolo un sistema interno (destrorso)  al blocco con:
    - l'asse z parallelo allo Z generale del dev
    - l'asse x nella direzione  del'asse del blocco
    - l'asse y perpendicolare ai due di prima.
    */

    DevCrn::CamInter();  // definisco il cammino "camm" nella corona

//    Controllo nel blocco
    gvet Bx = U_bl[iblocco];
    gvet Bz(0.,0.,1.);
    gvet By = Bz&Bx;
    By=By.NVerso();

    gvet X_p = X_dev-C_bl[iblocco];  // vettore distanza del punto dal centro blocco


// piani laterali paralleli in y  e lati in zeta
    double t=0.;
    double  bv=By%V_dev;
    if (fabs(bv)>0.)
    {
        if ( (t=(latoy-By%X_p)/bv) >0.  && t<camm) camm=t;
        else  if ( (t=-(latoy+By%X_p)/bv) >0. && t<camm) camm=t;
    }
    // e lati zeta
    bv=Bz%V_dev;
    if (fabs(bv)>0.)
    {
        if ( (t=(latoz-Bz%X_p)/bv) >0.  && t<camm) camm=t;
        else  if ( (t=-(latoz+Bz%X_p)/bv) >0. && t<camm) camm=t;
    }
    return camm;
}


double DevSveto::CamBuco()
{

    if (isector==-1) // è nel buco della corona
        return DevCrn::CamBuco();

// siamo dentro un settore tra due blocchi e all'interno della corona.

    DevCrn::CamInter(); // calcolo camm all'interno della corona

// verifico se colpisco le superfici del settore adiacenti ai leadglass.
// assi di riferimento e centro del settore
    gvet Sx = U_sc[isector];
    gvet Sz(0.,0.,1.);
    gvet Sy = Sz&Sx;
    Sy=Sy.NVerso();

// calcolo il cammino verso il piano di destra e confronto
    gvet X_p = X_dev-C_sc[isector];
    int ib=isector;
    gvet Bx = U_bl[ib];
    gvet By = Sz&Bx;
    By=By.NVerso();     // normale al piano di destra
    gvet Xpd= C_bl[ib]+ By*latoy-C_sc[isector];  // un punto sul piano di destra riferito al centro settore.

    double t=0.;
    double sv=Sy%V_dev;
    if (fabs(sv)>0. && (t=(Sy%Xpd-Sy%X_p)/sv) >0. && t<camm) camm=t;

// calcolo il cammino verso il piano di sinistra e confronto
    ib=(isector+1)%Blocchi;
    Bx = U_bl[ib];
    By = Sz&Bx;
    By=By.NVerso();     // normale al piano di destra
    Xpd= C_bl[ib]- By*latoy-C_sc[isector];  // un punto sul piano di destra riferito al centro settore.

    sv=Sy%V_dev;
    if (fabs(sv)>0. && (t=(Sy%Xpd-Sy%X_p)/sv) >0. && t<camm) camm=t;

// e lati z dello strato
    sv=Sz%V_dev;
    if (fabs(sv)>0.)
    {
        if ( (t=(latoz-Sz%X_p)/sv) >0.  && t<camm) camm=t;
        else  if ( (t=-(latoz+Sz%X_p)/sv) >0. && t<camm) camm=t;
    }

    return camm;


}