flyopunta/src/src_dev/device.cpp

/***************************************************************************
                          device.cpp  -  description
                             -------------------
    begin                : Sat May 26 2001
    copyright            : (C) 2001 by Giuseppe Pierazzini
    email                : peppe@unipi.it
  ***************************************************************************
  *                                                                         *
  *   NA48  simulation program.                                             *
  *                                                                         *
  ***************************************************************************/


#define CEPSI 1.0e-6
#include "flyoh.h"

//        Tirrenia   GmP 4.09.
#include <ctype.h>

#define ESTERNO  1
#define INTERNO  0
#define INBUCO   -1
#define EPSR     1.0E-200

using namespace GmpConst;
using namespace std;
static Device *last_device = 0;
static int tot_device = 0;
double sinPhi=0.0, cosPhi=1.0;
void Scrivi_Dati ( int );
// classe Base
Device::Device()
{
    working_buffer=MemoryBuffer;
    Make_working_buffer ( working_buffer );
    Make_measur_buffer ( working_buffer );
    // create e new device
    tot_device++;
    idev = tot_device;
    devtype=TypDevice;  //=0
    if ( last_device== 0 ) Apparato=this;
    else             last_device->next = this;   // linking
    up=last_device;
    last_device = this;
    this->next = 0;   // end linking

    // default orientation and zero setting
    Nx.setvn ( 1.,0.,0. );
    Ny.setvn ( 0.,1.,0. );
    Nz.setvn ( 0.,0.,1. );
    Tot_seen=0;
    for ( int i=0;i<40;i++ ) Part_hit[i]=0;
    nhit =ghcir= gnhit=nvede = -1;
    ghit_last=mhit=0;

    nowpos=prevpos=6;
    wrout=accedi=escludi=rivela=Zint=Zbuc=0;
    registra=1;
    campo=-1;
    reckik=-1. ; // reconstruction Kik non defined
    mattint="";
    mattbuc="";
    maxBrems=maxPair=0;    //if >1 enable bremsstrhalung  and pair production
    Brem_done=Pair_done=0;
    eDev=t_dead_mem =time_now=Ro =Dedx =0.0;
    Roint=Robuc=Dedxint=Dedxbuc=0.0;
    Aint=Abuc=lradi = lradb=lrad= 0.0;
    fun=Nodef;       // set default function at None...

// legge i dati di definizione del dev
    Leggi_schede_dev();
}

void Device::Leggi_schede_dev()
{
    // legge  le schede di definizione...
    string stringline,label,label2;
    char line[128];
    double ax,ay,az;
    int centro_dato= 0;

    Cardcur->getline ( line, 128 );
    while ( !Cardcur->eof() )
    {

        stringline=line;
        stringstream scheda;
        scheda<<line;
        scheda>>label;

        if ( Dbginput==1 )
        {
            Gout<<"\n Device label: "<<label<<"  Linea:"<<line;
        }

        if ( stringline.size()      <1  ||
                label.find ( "**" )  ==0  ||
                label.find ( "\\\\" ) ==0  ||
                label.find ( "//" )  ==0 )
        {
            Cardcur->getline ( line, 128 );
            continue;
        } //line empty or comment ...read next line


        if ( label.find ( "device" ) == 0 )
        {
            Gout<<"\n Device card out of sequence !!!\n "<<endl;
            cout<<"\n Device card out of sequence !!!\n "<<endl;
            exit ( 0 );
        }
        else if ( label.find ( "nome" ) == 0 )
        {
            scheda>>nome>>nickname;
        }

        else if ( label.find ( "flags" ) == 0 ||
                  label.find ( "fun" ) == 0 )
        {
            label="";
            scheda>>label;
            if ( label.find ( None ) == 0 )  fun = None;
            else if ( label.find ( Trig ) == 0 )   fun = Trig;
            else if ( label.find ( Veto ) == 0 )  fun = Veto;
            else if ( label.find ( Magn ) == 0 )  fun = Magn;
            else if ( label.find ( Trigh ) == 0 )  fun = Trigh;
            else if ( label.find ( Dead ) == 0 )   fun = Dead;
            else if ( label.find ( Dump ) == 0 )  fun = Dump;
            else if ( label.find ( Gasint ) == 0 ) fun =  Gasint;
            else if ( label.find ( Bersaglio ) == 0 ) fun =  Bersaglio;
            else
            {
                Gout<<"\n Device "<<nome<< "... flags error:  "<<label;
                exit ( 0 );
            }
        }
        else if ( label.find ( "detect" ) == 0 )
        {
            if ( nvede< 20 )
            {
                label="";
                scheda>>label;
                nvede++;
                const char *ptm;
                if ( ( ptm = Matter->Ptbsdat ( label.c_str() ) ) != 0 )
                {
                    vede[nvede] = ptm;
                }
                else if ( label.find ( All ) == 0 )  vede[nvede] = All;
                else if ( label.find ( Char ) == 0 ) vede[nvede] = Char;
                else if ( label.find ( Neut ) == 0 ) vede[nvede] = Neut;
                else
                {
                    Gout<<"\n Device parameter error ==> "<< label<<endl;
                    cout<<"\n Device parameter error ==> "<< label<<endl;
                    exit ( 0 );
                }
                assorb[nvede]=0;
                label="";
                scheda>>label;
                if ( label.size() !=0 )
                {
                    if ( label.find ( "enevd" ) ==0 ) assorb[nvede]=1;
                    else if ( label.find ( "dump" ) ==0 ) assorb[nvede]=2; // dump ma registra
                    else if ( label.find ( "noreg" ) ==0 ) assorb[nvede]=3; // non registra
                    else if ( label.find ( "cdump" ) ==0 ) assorb[nvede]=4; // dump ma non registra
                }
            }
        }
        else if ( label.find ( "esclu" ) == 0 )
        {
            escludi=1;
            accedi=100;
        }
        else if ( label.find ( "writeout" ) == 0 ) wrout=1;

        else if ( label.find ( "buffer" ) == 0 ) // working buffer: standard 30
        {
            scheda>>working_buffer;

        }



        else if ( label.find ( "centro" ) == 0 )
        {
            if ( centro_dato!=0 )
            {
                Gout<<"\n  Error: centro duplicato! in dev "<< idev<<endl;
                exit ( 1 );
            }
            centro_dato=1;
            scheda>>ax>>ay>>az;
            Centro.setvn ( ax,ay,az );
            Centro=Centro-Lab_Zero;
            Centro.Norma();
        }
        else if ( label.find ( "cface" ) == 0 ) // coord.centro superficie di ingresso
        {
            if ( centro_dato!=0 )
            {
                Gout<<"\n  Error: centro duplicato! in dev "<< idev<<endl;
                exit ( 1 );
            }
            centro_dato=2;
            scheda>>ax>>ay>>az;
            Cface.setvn ( ax,ay,az );

        }
        else if ( label.find ( "centrin" ) == 0 )   // centro del buco...coordinate interne
        {
            scheda>>ax>>ay>>az;
            Centrin.setvn ( ax,ay,az );
        }
        else if ( label.find ( "versx" ) == 0 )
        {
            scheda>>ax>>ay>>az;
            Nx.setvn ( ax,ay,az );
        }
        else if ( label.find ( "versy" ) == 0 )
        {
            scheda>>ax>>ay>>az;
            Ny.setvn ( ax,ay,az );
        }
        else if ( label.find ( "versz" ) == 0 )
        {
            scheda>>ax>>ay>>az;
            Nz.setvn ( ax,ay,az );
        }
        else if ( label.find ( "rotap" ) == 0 )  // START RUGGIERO
        {
            double Phi;
            scheda>>Phi;
            Phi*=PiGreco/180.;
            sinPhi = sin ( Phi );
            cosPhi = cos ( Phi );
// la rotazione è attorno all'asse zeta del laboratorio...
// ridefinisco i versori del device rispetto al sistema del lab.
            Nx.setvn ( cosPhi, sinPhi,0.0 );
            Ny.setvn ( -sinPhi, cosPhi,0.0 );
            Nz.setvn ( 0.0,    0.0,   1.0 );
        }                                          // END RUGGIERO

        else if ( label.find ( "sizein" ) == 0 ||
                  label.find ( "rin" ) == 0 )
        {
            scheda>>ax>>ay>>az;
            Lin.setvn ( ax,ay,az );
        }

        else if ( label.find ( "totin" ) == 0 ||
                  label.find ( "dtin" ) == 0 )
        {
            scheda>>ax>>ay>>az;
            Lin.setvn ( ax*.5, ay*.5, az*.5 );
        }

        else if ( label.find ( "sizeout" ) == 0 ||
                  label.find ( "rout" ) == 0 )
        {
            scheda>>ax>>ay>>az;
            Lout.setvn ( ax,ay,az );
            Rqsize=Lout.normaq*1.1;
        }
        else if ( label.find ( "totout" ) == 0 ||
                  label.find ( "dtout" ) == 0 )
        {
            scheda>>ax>>ay>>az;
            Lout.setvn ( ax*.5,ay*.5,az*.5 );
            Rqsize=Lout.normaq*1.1;
        }
        else if ( label.find ( "opens" ) == 0 )
        {
            scheda>> OpenSector;
            OpenSector*= 0.017453293;
        }
        else if ( label.find ( "fortax" ) == 0 )
            // per definire dimensioni semiapertura del foro finale  interno dei tax
        {
            scheda>>ax>>ay>>az;
            Fortax.setvn ( ax,ay,az );
        }

        else if ( label.find ( "field" ) == 0 )
        {
            scheda>>ax>>ay>>az;
            Field.setvn ( ax,ay,az );
            Bref=ax;
            Field_sav=Field;
        }
        else if ( label.find ( "vfiel" ) == 0 )
        {
            //      Bref e posizione valore di riferimento del campo
            //      per il focussing quadrupoli.....
            scheda>>Bref ;
            scheda>>Dbref;
            Field.setvn ( Bref,0.0,0.0 );
            Field_sav=Field;
        }
        else if ( label.find ( "magki" ) == 0 )
            scheda>>reckik;
        else if ( label.find ( "mem" ) == 0 )
        {
            scheda>> t_dead_mem ;
            if ( Tpipe<t_dead_mem ) Tpipe=t_dead_mem;
        }

        else if ( label.find ( "registra" ) == 0 )
        {
            scheda>> label;
            if ( label=="no" ) registra=0;  //default=1;
        }

        else if ( label.find ( "materiali" ) == 0 )
        {

            scheda>>mattint>>mattbuc;
            Get_Materiali();
        }
        else if ( label.find ( "brems" ) == 0 ) scheda>>maxBrems;
        else if ( label.find ( "pair" ) == 0 ) scheda>>maxPair;

        //--------------------------------------------------
        //  the  following procedure is obsolate
        else if ( label.find ( "radi" ) == 0 )
        {
            Gout<<"\n In dev "<<nome
            <<"  \"radi\" too old ... to be subtituted by \"materiali\" card ...thanks " ;
        }
        else if ( label.find ( "dens" ) == 0 )
        {
            Gout<<"\n In dev "<<nome
            <<"  \"dens\" too old ... to be subtituted by \"materiali\" card ...thanks ";
        }
        //----------------------------------------------------


        //exit dev generation
        else if ( label.find ( "end" ) == 0 ||
                  label.find ( "done" ) == 0 )
        {
            if ( centro_dato==2 ) // calcolo il vero centro del dev
            {
                Centro=Cface-Lab_Zero;
                Centro[2]+=Lout[2];
                Centro.Norma();
            }
            else if ( centro_dato==1 )
            {
                Cface=Centro;
                Cface[2]-=Lout[2];
                Cface.Norma();
            }
            else
            {
                Gout<<"\n Error: Missing centro in dev "<<idev<<endl;
                exit ( 1 );
            }

            if ( escludi==1 ) wrout=0;

// attenzione  aggiustamento memoria di lavoro se richieste dimensioni maggiori di 30 !!!!!!!
// device buffers: local, circular, after dev simulation and measured buffer for root using.
// standard dim=30

            Make_working_buffer ( working_buffer );
            Make_measur_buffer ( working_buffer );


            return;
        }
        //  next for debug.....
        else if ( label.find ( "stopr" ) == 0 )
        {
            Gout<<"\n Stop run by operator. "<<endl;
            cout<<"\n Stop run by operator. "<<endl;
            exit ( 0 );
        }
        Cardcur->getline ( line, 120 );
    }




}
//---------------------------------
Device::~Device()
{}
//--------------------------------------------
// making buffer procedures
void  Device::Make_working_buffer ( int dim )
{
    if ( dim>MemoryBuffer )
    {
        delete[] Hits;
        delete[] G_Hits;
        delete[] M_Hits;
        Gout<<"\n In device ** "<<nome <<" **  working_buffer adjusted to ==> "<<dim;
    }

    circular_buffer=dim;
    if ( Pipeline>0 )  circular_buffer=4*dim;
    Hits  = new Bufdev[dim];
    G_Hits= new Bufdev[circular_buffer];
    M_Hits= new Bufdev[dim];

}
void  Device::Make_measur_buffer ( int dim )
{



    if ( dim>MemoryBuffer )
    {
// nel primo dev non di delete.....
//if(idev==2) Gout<<"\n old pointer  "<< index<<" "<< fdx[0]<<" "<<&index<<" "<< &fdx;


        delete[] index;
        delete[] mid;
        delete[] mdead;
        delete[] devstory;
        delete[] mxd;
        delete[] myd;
        delete[] mzd;
        delete[] mdd;
        delete[] me;
        delete[] mxl;
        delete[] myl;
        delete[] mzl;
        delete[] gpx;
        delete[] gpy;
        delete[] gpz;
        delete[] ge;
        delete[] dtemp;
        delete[] structinfo;

        Gout<<"\n In device ** "<<nome <<" **  measur_buffer adjusted to ==> "<<dim;

    }

    index=new int[dim];
    mid =new int[dim];
    mdead=new int[dim];
    devstory= new unsigned long long[dim];
    mxd=new float[dim];
    myd=new float[dim];
    mzd=new float[dim];
    mdd=new float[dim];
    me=new float[dim];
    mxl=new float[dim];
    myl=new float[dim];
    mzl=new float[dim];
    gpx=new float[dim];
    gpy=new float[dim];
    gpz=new float[dim];
    ge=new float[dim];
    dtemp=new double[dim];
    structinfo=new int[dim];


}

void Device::Reset_Pos()
{
    // start from the tree top
    Device *pdev=Apparato;
    while ( pdev!=0 )
    {
        nowpos=prevpos=6;
        pdev=pdev->next;
    }
}

//***********************  R E S E T ****************************

void Reset_Devices()
{
// attenzione dove si mette questo reset... per la relazine tempoarle
    // start from the tree top...
    Device *pdev=Apparato;
    while ( pdev!=0 )
    {
        pdev->Reset_Dev();
        pdev=pdev->next;
    }
}
void Device::Reset_Dev()
{
// attenzione dove si mette questo reset... si suppone all'inizio evento
// quando si conosce il tempo di generazione del nouvo evento ....

    nowpos=prevpos=6;  //posizione non definita
    accedi=0;
    rivela=0;
    eDev=0.0;
    Brem_done=Pair_done=0;
    nhit=-1;
    mhit = 0;

    if ( Pipeline==0 )
    {
        gnhit = -1;  // hits generati
    }
}

//====================   s c r i v i   c o l p i  ================
void  Device::ScrvColpi ( Particella *pr )
{
//  cout<<"\n Ev Scrv "<<Eventi_Fatti<< " dev "<<nome<< " pr "<<pr->nome<<" colpo "<<colpo<<" visto "<<visto;

    if ( nhit> ( working_buffer-2 ) ) //verifica la size del buffer!
    {
        Gout<<"\n Ev "<< Eventi_Fatti<<"  In dev "<<nome<<" too many particles "<< nhit<< "  get out! "<<endl;
        Scrivi_Dati ( 2 );
        flyoend();
        exit ( 0 );
    }


// aggiorna il contatore di dev per le particelle che lo attraversano
    Part_hit[Idsee]++;
    Tot_seen++;

    pr->devhit=idev;  //memorizzo l'ultimo indice del dev colpito.

    // salva in formato codificato i devices colpiti dalla particell
    // il bit acceso corrisponde al device colpito
    // attenzione memorizza solo i primi 64 device.
    pr->Update_Devstory ( idev );


    if ( rivela>2 ) return;   // non registered
// attenzione si evita anche  di registrare  gli hits per i dev che non sono esplicitamente
// usati nella analisi successiva. registra=(0,1)==(no,si)
    if ( registra==0 ) return;

    nhit++;       // remember nhit starts from -1 that corrisponds to no  hits...

    Bufdev  *bf = &Hits[nhit];
    bf->Pr=             pr;
    bf->tempo=          evento_.Gen.tempo;
    bf->idead=          pr->Ifato;
    bf->devstory=       pr->devstory;
    bf->last_path=      pr->last_path;
    bf->pnome =         pr->nome;
    bf->id=             pr->id;
    bf->idm=            pr->idm;
    bf->Vers=           pr->Vers;
    bf->Xlab=           pr->X;
    bf->Xdev=           Lab2cDev ( pr->X );
    bf->Pdev=           pr->P;
    bf->structinfo= Get_Structinfo();

//  ----------------------
    if ( rivela>0 ) bf->e_rivela=pr->P.e;
    else bf->e_rivela=0.0;   //it is done anyway if multiplescattering if required...

//  ------------- pipeline test -------
    if ( Pipeline==0 ) return;

// pipeline is on! copies the hits in the circular buffer

    gnhit++;       // hits index up to now...
    ghcir=gnhit%circular_buffer;  //attention... Circular buffer
    G_Hits[ghcir]=bf[0];
    G_Hits[ ( ghcir+1 ) %circular_buffer].tempo=-10;  //set to -10 next hit not yet arrived....

//    Gout<<"\n Eve_gen "<<evento_.Gen.Event<<" dev "<<nome<< " hit "<<gnhit<<" "<<ghcir<<" tempo "<< G_Hits[ghcir].tempo
//     <<" " << G_Hits[(ghcir+1)%circular_buffer].tempo;
    ;

};



void  Device::GetWindowHits ( double trigger_time )
{



    double dtm, adtm;
    int pipe=0;
    nhit=-1;
    mhit=0;


    if ( gnhit<0 ) return;  // gnhit: total hit seen up to now memorized in the circular buffer

    double tmorto=t_dead_mem;  // get dead time of the device
    if ( tmorto>0.0 ) tmorto=t_dead_mem +0.001;
    else tmorto=0.001;   //if zero set to 0.001 ns

    int ghit=ghit_last; // get the last hit
    ghit_last=0;

// Gout<<"\n In Dev " <<nome << " for  trig_time  " <<trigger_time;
// cout<<"\n In Dev " <<nome << " for  trig_time  " <<trigger_time;

    int ih=0;
    while ( ih<circular_buffer )
    {
        ghit=ghit%circular_buffer;
        if ( G_Hits[ghit].tempo<-9. ) break; // test on last hit to break the loop...
        dtm = trigger_time -G_Hits[ghit].tempo;  // difference between trigger_time and hit_time...
        adtm=fabs ( dtm );

        if ( adtm<=tmorto )  //Test if the hit is inside the dead_time....
        {
            pipe=0;
            nhit++;
            if ( adtm>0.001 )
            {
                pipe=1;
                evento_.Gen.Pipev=1;
            }
            M_Hits[nhit]=G_Hits[ghit];   //copy from circular buffer into M_Hits...
            M_Hits[nhit].Pileup=pipe;
            if ( ghit_last==0 ) ghit_last=ghit;

        }
        else if ( dtm>0.0 ) ghit_last=ghit;
        else if ( dtm<0.0 ) break;

//      Gout<<" Lasthit "<<setw(5)<<ghit_last<<"  pipe "<<pipe ;

        if ( nhit> ( working_buffer-2 ) )
        {
            Gout<<"\n Ev "<< evento_.Gen.Event<<"  Too many hits in dev " << nome<< "  see GetWindowHits()... exit" <<endl;
            exit ( 0 );
        }
        ghit++;
        ih++;
    }
    mhit=nhit+1;
    //   if(mhit>0)Gout<<"\n Ev "<< Eventi_Fatti<<"  In dev "<<nome<<"  Collected hits "<<mhit;
}


//***************************************************
//======== fill default dev coordinate type(2) with type(1).====
int Device::Fill_Default()
{
    mhit=nhit+1;
//  Gout<<"\n Default copy Ev "<< evento_.Gen.Event<<" dev "<<nome<< " hits "<< mhit;

    for ( int i=0;i<mhit;i++ )
    {
        M_Hits[i] =Hits[i];
    }
    return mhit;
}

//===============================================================
int Device::Get_DatiOut()
{
    hit_recorded=mhit;

    emDev=0.;
    if ( mhit>0 )
    {
// test on the total hits
        if ( mhit>working_buffer )
        {
            Gout<<"\n Ev "<< Eventi_Fatti<<" <Measur_buffer> too low  in "<<nome<< " Detected hits: "<<mhit
            << " cut to  "<<working_buffer;
            mhit=working_buffer;
        }

        hit_recorded=mhit;
        for ( int hts=0;hts<hit_recorded;hts++ )
        {
            index[hts]=M_Hits[hts].id;
            mid[hts]=M_Hits[hts].idm;
            mdead[hts]=M_Hits[hts].idead;
            structinfo[hts]=M_Hits[hts].structinfo;
            devstory[hts]=M_Hits[hts].devstory;
            M_Hits[hts].Xdev.putv ( mxd[hts],myd[hts],mzd[hts] );
            mdd[hts]= M_Hits[hts].Xdev.XYNorma();
            me[hts]=M_Hits[hts].e_rivela;
            M_Hits[hts].Xlab.putv ( mxl[hts],myl[hts],mzl[hts] );
            M_Hits[hts].Pdev.putq ( gpx[hts],gpy[hts],gpz[hts],ge[hts] );
            dtemp[hts]=M_Hits[hts].tempo-evento_.Gen.tempo;
            emDev+=me[hts];

        }
    }

    return hit_recorded;
}

//*********************************************************

//------------------------------------------
void  Device::Get_Materiali()
{
    MattPar *Elem;
    // nel interno del dev
    //   Gout<<"\n Elemeto interno "<< mattint;
    if ( mattint.size() !=0 )
    {
        Elem= Elemento->PntElem ( mattint );
        Zint   =Elem->Get_Z();   // carica atomica
        Aint   =Elem->Get_A();   // numero atomico
        liint  =Elem->Get_Intl();// lungh. interazione cm
        Roint  =Elem->Get_Ro();  // densita' g/cm3
        Dedxint=Elem->Get_Dedxl();// perdita energia Gev/cm
        lradi  =Elem->Get_Radl();//  rad in cm
    }

    if ( mattbuc.size() !=0 )
    {
        Elem= Elemento->PntElem ( mattbuc );
        Zbuc   =Elem->Get_Z();   // carica atomica
        Abuc   =Elem->Get_A();   // numero atomico
        libuc  =Elem->Get_Intl();// lungh. interazione cm
        Robuc  =Elem->Get_Ro();  // densita' g/cm3
        Dedxbuc=Elem->Get_Dedxl();// perdita energia Gev/cm
        lradb  =Elem->Get_Radl();//  rad in cm
    }
}
//------------------------------------------
void  Device::Get_Pointers()
{


    // to be updated according to the device names....
    // get dev pointers
    Gout<<"\n Device::Get_pointers() : get the device pointers to speed up the simulation..."
    <<"\n The user must adjust the device labels occording to names as defined in  .app file \n";


    collo=   Apparato->Get_Pointer ( "Collo" );
    xch02=   Apparato->Get_Pointer ( "XCH02" );
    spione0 = Apparato->Get_Pointer ( "spione0" );
    spione1 = Apparato->Get_Pointer ( "spione1" );
    vtappo = Apparato->Get_Pointer ( "Vtappo" );

    dectubo= Apparato->Get_Pointer ( "dectubo" );
    kabs1= (Devkabes*)  Apparato->Get_Pointer ( "Spibes1" );
    kabs2= (Devkabes*)  Apparato->Get_Pointer ( "Spibes2" );
    kabs3= (Devkabes*)  Apparato->Get_Pointer ( "Spibes3" );

    giga1=  (DevGtk*) Apparato->Get_Pointer ( "Gtk1" );
    giga2=  (DevGtk*) Apparato->Get_Pointer ( "Gtk2" );
    giga3=  (DevGtk*) Apparato->Get_Pointer ( "Gtk3" );
         cedar=  (DevCedar*) Apparato->Get_Pointer ( "cedar" );
    mgb1 =   Apparato->Get_Pointer ( "MCB80" );   // magnet pointers of acromat 2
    mgb2 =   Apparato->Get_Pointer ( "MCB84" );
    mgb3 =   Apparato->Get_Pointer ( "MCB95" );
    mgb4 =   Apparato->Get_Pointer ( "MCB97" );
    mgb5 =   Apparato->Get_Pointer ( "MDX99" );

    Chanti =   (DevChanti*) Apparato->Get_Pointer ( "chanti" );

    // get anticounter pointers...

    anti[0]=  (DevLav*) Apparato->Get_Pointer ( "Ant1" );
    anti[1]=  (DevLav*) Apparato->Get_Pointer ( "Ant2" );
    anti[2]=  (DevLav*) Apparato->Get_Pointer ( "Ant3" );
    anti[3]=  (DevLav*) Apparato->Get_Pointer ( "Ant4" );
    anti[4]=  (DevLav*) Apparato->Get_Pointer ( "Ant5" );
    anti[5]=  (DevLav*) Apparato->Get_Pointer ( "Ant6" );
    anti[6]=  (DevLav*) Apparato->Get_Pointer ( "Ant7" );
    anti[7]=  (DevLav*) Apparato->Get_Pointer ( "Ant8" );
    anti[8]=  (DevLav*) Apparato->Get_Pointer ( "Ant9" );
    anti[9]=  (DevLav*) Apparato->Get_Pointer ( "Ant10" );
    anti[10]= (DevLav*) Apparato->Get_Pointer ( "Ant11" );
    anti[11]= (DevLav*) Apparato->Get_Pointer ( "Ant12" );


    Cam1= (DevMwc*)  Apparato->Get_Pointer ( "Cam1" );
    Cam2= (DevMwc*)  Apparato->Get_Pointer ( "Cam2" );
    Cam3= (DevMwc*)  Apparato->Get_Pointer ( "Cam3" );
    Cam4= (DevMwc*)  Apparato->Get_Pointer ( "Cam4" );

    Rich= (DevRich*)   Apparato->Get_Pointer ( "Rich" );

    St1x= (DevStraw*)  Apparato->Get_Pointer ( "St1x" );  // la coordinata è parallela ai fili attenzione
    St1y= (DevStraw*)  Apparato->Get_Pointer ( "St1y" );
    St1u= (DevStraw*)  Apparato->Get_Pointer ( "St1u" );
    St1v= (DevStraw*)  Apparato->Get_Pointer ( "St1v" );
    St2x= (DevStraw*)  Apparato->Get_Pointer ( "St2x" );
    St2y= (DevStraw*)  Apparato->Get_Pointer ( "St2y" );
    St2u= (DevStraw*)  Apparato->Get_Pointer ( "St2u" );
    St2v= (DevStraw*)  Apparato->Get_Pointer ( "St2v" );
    St3x= (DevStraw*)  Apparato->Get_Pointer ( "St3x" );
    St3y= (DevStraw*)  Apparato->Get_Pointer ( "St3y" );
    St3u= (DevStraw*)  Apparato->Get_Pointer ( "St3u" );
    St3v= (DevStraw*)  Apparato->Get_Pointer ( "St3v" );
    St4x= (DevStraw*)  Apparato->Get_Pointer ( "St4x" );
    St4y= (DevStraw*)  Apparato->Get_Pointer ( "St4y" );
    St4u= (DevStraw*)  Apparato->Get_Pointer ( "St4u" );
    St4v= (DevStraw*)  Apparato->Get_Pointer ( "St4v" );

    mnp1 = (DevMagn*) Apparato->Get_Pointer ( "Mnp1" );
    mnp2 = (DevMagn*)  Apparato->Get_Pointer ( "Mnp2" );
    Mag48= (DevMagn*)  Apparato->Get_Pointer ( "Magn" );
    Mnp33= (DevMagn*)  Apparato->Get_Pointer ( "Mnp33" );
    hodo =  (DevChod*) Apparato->Get_Pointer ( "Hodo" );
    hodox = (DevChod*) Apparato->Get_Pointer ( "Hodx" );
    hodoy = (DevChod*) Apparato->Get_Pointer ( "Hody" );

    hadc =  (DevHac*)  Apparato->Get_Pointer ( "HACM" );
    iron  =  Apparato->Get_Pointer ( "Ironwall" );
    muvet1= (DevMuv*) Apparato->Get_Pointer ( "MUVET1" );
    muvet2= (DevMuv*) Apparato->Get_Pointer ( "MUVET2" );
    muvet3= (DevMuv*) Apparato->Get_Pointer ( "MUVET3" );
    monit = Apparato->Get_Pointer ( "Moni" );
    mamud=   Apparato->Get_Pointer ( "Mamu" );

    sacdev= (DevSac*)Apparato->Get_Pointer ( "Sac" );
    ircdev=(DevIrc*)Apparato->Get_Pointer ( "IRC" );
    lkry=   (DevLkr*) Apparato->Get_Pointer ( "Lkry" );
    if ( lkry==0 ) lkry= (DevLkr*)Apparato->Get_Pointer ( "Lkr" );


}
//------static Device *Walls,*lkry;
//  static Device *Coll, *Updet,*Veto,*Buco,*Walls,*lkry;

Device * Device::Get_Pointer ( const  char *word )
{
    Device *dv,*dvok=0;
    dv=Apparato;
    while ( dv!=0 )
    {
        if ( dv->escludi==0 &&
                strcmp ( dv->nome,word ) ==0 ) dvok= dv;
        dv=dv->next;
    }
    Gout<<"\n Device "<< setw ( 9 ) <<word<<" pointer "<< ( dvok?" found":" === not ===  found or excluded" );
    return dvok;
}



//------------------------------------------
//  accedi=escludi=0;     trasformazioni lab -> dev -> lab
//--------------------------------------------
gvet& Device::Lab2Dev ( gvet &Vlab )
{
    // trasformazioni di vettori direzione lab-> dev.
    static gvet _Vdev;
    _Vdev.setvn ( Vlab%Nx,Vlab%Ny,Vlab%Nz );
    return _Vdev;
}
//--------------------------------------------
gvet& Device::Lab2cDev ( gvet &Xin )
{
    // trasformazione di vettori al dev con traslazione...

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

//--------------------------------------------
gvet& Device::Dev2Lab ( gvet &Vdev )
{
    // trasformazioni di vettori direzione dev->lab.
    // calcola la norma.
    static gvet _Vlab;
    _Vlab=Nx*Vdev[0] + Ny*Vdev[1] + Nz*Vdev[2];
    _Vlab.normaq= _Vlab[0]*_Vlab[0]+ _Vlab[1]*_Vlab[1]+_Vlab[2]*_Vlab[2];
    _Vlab.norma= sqrt ( _Vlab.normaq );
    return _Vlab;
}
//--------------------------------------------
gvet& Device::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;
}

//================   T R A C C I A ====================================
//================   T R A C C I A ====================================

int Device::Traccia ( Particella *pr )
{
    int iter=0;
    // impulso ==0 no tracing ... return;
    if ( pr->P.norma == 0.0 ) return 1;

    lint=nowpos?libuc:liint; //get the interaction length


    if ( fun==Bersaglio&& pr->rg!=0&&lint>0.0&&nowpos<1 )
    {   // se qui è un caso particolare... particelle si distrugge  nel dev 
        // Il dev fa da bersaglio 
        // ed è interna o nel buco ( è qui c'è di default)
        // e la lunghezza di interazione è definita non nulla!
        // si ridefisce il cammino se interagisce nel bersaglio
        pr->Get_Bersaglio ( this,lint );
        if ( pr->Ifato==-3 )
        { // la particella è stata distrutta!!!
            pr->Move ( pr->pathok );  // nuovo path (ma vedere meglio)
            pr->last_path=pr->pathok;
            pr->rg->dw->X=pr->rg->X=pr->X;    // definisco la posizione del C.M.
            pr->rg->dw->Gx=pr->rg->Gx=pr->X;
            return -3;   // la bullet muore nel dev
        }
    }

    double tratok=pr->pathok;
    lrad=nowpos?lradb:lradi; //get the radiation  length
    if ( lrad==0. )     // no radiation...
    {
        // free flight
        pr->Move ( tratok );
        pr->last_path=tratok;
        return 0;
    }

// qui lrad != 0

    if ( pr->charg ==0.0 )
    {

// se è un gamma fa Pair_Production
        if ( pr->Pair_Production>0&& Pair_done<maxPair ) //        Pair production activated?
        {
            if ( pr->P.e<0.002 ) // Energy cutoff .... fare meglio
            {
                Hits[nhit].e_rivela+=pr->P.e;  // aggiorno l'energia persa...
                Hits[nhit].idead=-18;
                pr->Ifato=-18;
                pr->P.setqn ( 0.,0.,0.,0. );
                pr->pathok=0.0;
                return -18;
            }
            pr->Get_Pair ( this,tratok,lrad );
            Hits[nhit].idead=pr->Ifato;
            pr->pathok=tratok; //redefine (if nedeed) free path before pair prod..
        }
        pr->Move ( tratok );
        pr->last_path=tratok;
        return 0;
    }
//  qui lrad != 0   ed è carica 
    int crepa;
    double passdone=0.0;
    double passo= lrad;

// fa bremsstrahlung

    while ( passdone<tratok )
    {
        double diff=tratok-passdone;
        if ( diff <passo ) passo =diff+0.001;

        if ( pr->Brems_Production>0 && Brem_done<maxBrems )
        {
            if ( pr->Get_Brems ( this,passo,lrad ) !=0 )
            {
                pr->Move ( passo );  // divido per 2  ??? media
                Mscatter ( pr,passo );
                pr->last_path=passdone+passo;
                Hits[nhit].idead=pr->Ifato;
                return 0;
            }
        }
        // qui se non c'è stato brems...
        pr->Move ( passo );   //forse passo/2  ???!|!!
        if ( ( crepa=Mscatter ( pr,passo ) ) !=0 )
        {
            pr->last_path=passdone+passo;

            Hits[nhit].idead=crepa;
            return crepa;
        }
        passdone+=passo;
        iter++;
    }
    Hits[nhit].idead=pr->Ifato;
    Hits[nhit].last_path =pr->last_path=tratok;
    return 0;

}

//***************************************************
//           Set_Rivela
//--------------------------------------------
void Set_Rivela_in_Dev ( Particella *pr )
{
    // start from the tree top...
    Device *pdev=Apparato;
    while ( pdev!=0 )
    {
        pdev->Set_Rivela ( pr );
        pdev=pdev->next;
    }
}

//--------------------------------------------
void Device::Set_Rivela ( Particella *pr )
{
    int isee;
    const char *label;



    if ( escludi!=1 )
    {
        // dev visibile
        accedi=0;
        rivela=-1;
        // test if the dev sees the particle ...first by name
        isee = nvede;
        while ( isee >= 0 )
        {
            label = vede[isee];
            if ( label == pr->nome )
            {
                rivela=assorb[isee];
                Idsee=isee;
                break;
            }
            isee--;
        }
        // ... and now  by type
        if ( rivela==-1 )

        {
            isee = nvede;
            while ( isee >= 0 )
            {
                label = vede[isee];
                if ( label == All ||label == pr->type ||
                        ( label == Char && pr->charg!=0. ) )
                {
                    rivela=assorb[isee];
                    Idsee=isee;
                    break;
                }
                isee--;
            }
        }
    }

// if(idev>1&&idev<4)
// Gout<<"\n Riv "<< nome <<" "<<pr->nome<<" " <<accedi << " " <<rivela;


}

//--------------------------------------------
void Device::Dev_summary()
{
    if ( Tot_seen!=0 && escludi!=1 )
    {
        int ns;
        for ( int i=0;i<40;i++ ) if ( ( ns=Part_hit[i] ) !=0 )
            {
                Gout.setf ( ios::fixed, ios::floatfield );
                Gout<<"\n Dev "
                <<setw ( 10 ) <<nome
                <<" "<<setw ( 3 ) <<idev
                <<" sees "<<setw ( 9 ) <<ns
                <<"  "<<setw ( 6 ) <<vede[i]
                <<"   Tot: "<<setw ( 9 ) <<Tot_seen<<setprecision ( 2 )
                <<"  Z_lim.["<<setw ( 8 ) <<Centro.z-Lout.z
                <<" <---> "<<setw ( 8 ) <<Centro.z+Lout.z<<"]"<<flush;
            }
    }
    if ( next!=0 ) next->Dev_summary(); //iterate !

}

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

void Device::PrgeomId()
{

    if ( this==Apparato )
    {
        Gout<<"\n\n==================  Detectors block =====================\n";
        Gout<<"\n    Id     nome   niknome  vs.    Z_in        Z_out      Wrk_buff          type \n";
    }
    Gout.setf ( ios::fixed, ios::floatfield );
    Gout<<setprecision ( 2 ) <<"\n "<<setw ( 3 ) <<idev
    <<setw ( 10 ) <<nome<<"   ["<<setw ( 6 ) <<nickname<<"] "
    <<" "<<escludi
    << "    "<<setw ( 8 ) << Cface.z
    << "    "<<setw ( 8 ) << Cface.z+2.*Lout.z
    << "    Wrb "<<setw ( 4 ) <<working_buffer
    << "   "<<setw ( 4 ) <<devtype<<" <=> "<<devclass;
    if ( next!=0 ) next->PrgeomId();
    else Gout<<"\n===========  End apparatus list   ============\n \n \n";
}


void Device::Prgeom()
{


    Gout<<"\n Dev   Id ="<<setw ( 4 ) <<idev
    <<"\n Dev Type ="<<setw ( 4 ) <<devtype
    <<"\n Nome = "<<nome<< "    Code: "<<nickname;
    if ( escludi==1 ) Gout<<"\n\n Rivelatore  E S C L U S O  dall'utente";

    Gout<<endl;
    if ( Cface.norma>0. ) Cface.print ( "S.Lab" );
    Centro.print ( "C_lab" );
    gvet Cfine=Cface;
    Cfine[2]+=Lout[2]*2.;
    Cfine.print ( "E.Lab" );
    Lout.print ( "Lout" );
    if ( Lin.norma>0. )
    {
        if ( Centrin.norma>0. ) Centrin.print ( "C_in" );
        Lin.print ( "Lin" );
    }
    if ( Fortax.norma>0. )
    {
        Fortax.print ( "Linend" );
    }
    Gout<<endl;
    Nx.print ( "Nx" );
    Ny.print ( "Ny" );
    Nz.print ( "Nz" );
    //campo.....

    if ( campo==0 ) Field.print ( "Field" );
    if ( campo==1 ) Gout<<"\n Quad Bref                "<<setprecision ( 3 ) <<fixed<<setw ( 11 ) <<Bref;
    if ( campo==2 ) Gout<<"\n Campo magnetico variabile!";
    if ( campo==0 ) Gout<<"\n Reconstruction magic kik "<< setw ( 7 ) <<reckik;
    Gout<<"\n Memoria in ns             "<<setw ( 6 ) <<t_dead_mem;
    Gout<<"\n Working buffers           "<<setw ( 6 ) <<working_buffer;
    if ( wrout==1 ) Gout<<"\n Writeout ok";
    else Gout<<"\n Writeout off...";

    Gout<<endl;
    Gout.setf ( ios::floatfield );
    if ( lradi != 0.0 || lradb != 0.0 )
    {
        Gout.setf ( ios::right,ios::adjustfield );
        Gout<<setfill ( ' ' ) <<setprecision ( 5 )
        <<"\n Dev  stuff:           "<<setw ( 10 ) <<mattint<<"  in hole: "<<setw ( 10 ) <<mattbuc<<"  "
        <<"\n DensInt, DensBuco     "<<setw ( 10 ) <<Roint<<"  "<<setw ( 10 ) <<Robuc<<"  g/cm3"
        <<"\n Zint, Zbuc            "<<setw ( 10 ) <<Zint <<"  "<<setw ( 10 ) <<Zbuc<<"  Atom. charges"
        <<"\n Aint, Abuc            "<<setw ( 10 ) <<Aint <<"  "<<setw ( 10 ) <<Abuc<<"  Atom. Numbers"
        <<"\n Radin,Radb            "<<setw ( 10 ) <<lradi<<"  "<<setw ( 10 ) <<lradb<<"  "
        << ( Scatter==0?"M.S. off":"M.S. on " )
        <<"\n InterInt,InterBuco    "<<setw ( 10 ) <<liint<<"  "<<setw ( 10 ) <<libuc
        <<"\n Dedxint,Dedxbuc       "<<setw ( 10 ) <<Dedxint<<"  "<<setw ( 10 ) <<Dedxbuc<<"  Gev/cm  "
        << ( Perdita==0?"Loss off":"Loss on" )
        <<"\n MaxBrems, MaxPair     "<<setw ( 10 ) << maxBrems<<"  "<<setw ( 10 ) <<maxPair;




    }
    Gout.setf ( ios::floatfield );

    Gout<<"\n";
    int vst;
    int idt=0;;
    if ( nvede>-1 ) Gout<<"\n Lista delle particelle rivelate.\n";
    while ( idt<=nvede )
    {
        vst=assorb[idt];
        Gout<<"\n"<<idt<<" ->   Detect: < "<<setw ( 7 ) <<vede[idt]<<" >  "<<setw ( 7 ) <<
        ( vst==0?"Seen": vst==1?"Dump":vst==2?"Seen but not reg.":"Dump but not reg." );
        idt++;
    }
    Gout<<"\n---------------\n";
}


//==================Get_camm ==================

Device *Device::Get_camm ( Particella *pr )
{
    //  double tratto,path;  ...in device.h
    //  int posok,pos;
    //  double Last_distq=1.e+12;
    Device *Dvk,*Dev;
    //
    path=1.e+10;
    Dvk=0;
    Dev=this;    // device di partenza per la loop
    posok=3;


    while ( Dev!=0 )
    {


        if ( Dev->accedi>0|| Dev->escludi>0||
                Dev->rivela<0 )      // dev non visibile o non vede
        {

            if ( Dev == Devloop ) Devloop = Devloop->next; //  aggiorno Devloop
            Dev=Dev->next;
            continue;
        }



        //--------------------=================----------------
        Dev->pr_dev=pr;
        Dev->X_dev=Dev->Lab2cDev ( pr->X ); // trasla al centro dev e ruota
        Dev->V_dev=Dev->Lab2Dev ( pr->Vers ); // routa il versore.... att. aggiornato!

        Dev->X_dev.Norma();
        if ( Dev->X_dev.normaq>Dev->Rqsize )  // test grossolano ma veloce  per dev esterno
        {
            dcs=Dev->X_dev%Dev->V_dev;   // dist. dal centro dev proiet. sul cammino
            if ( dcs>=0.0 )  // dev superato
            {
                if ( Dev == Devloop ) Devloop = Devloop->next; //  aggiorno Devloop
                Dev->accedi=13;

                Dev=Dev->next;
                continue;
            }
            else if ( Dev->X_dev.normaq-dcs*dcs>Dev->Rqsize )
            {
                //  Dev->accedi=12;   //no! se lo scattering o un campo cambia la direzione di volo!

                Dev=Dev->next;
                continue;
            }
        }
        // check position carefullly.....


        pos = Dev->Posizione();

//    if ( pr->idm==14 ) Dev->SpiaPr ( " Posiz ", 49,49 ,-99.);





        if ( pos==1 )              //Esterno
        {
            // calcolo il cammino per colpire il dev.....
            Dev->tratto=Dev->CamEster();
            //    if ( pr->idm==14 ) Dev->SpiaPr ( " Esterno ", 49,49 , Dev->tratto);



            if ( Dev->tratto>Dev->X_dev.norma ) Dev->tratto=Dev->X_dev.norma;
//      if ( pr->idm==14 ) Dev->SpiaPr ( " Ester ", 48,48 );

            if ( Dev->tratto<0.0001 )          // Dev not accessible or just crossed
            {
                Dev->accedi = 14;
                if ( Dev == Devloop ) Devloop = Devloop->next; //  aggiorno Devloop
                Dev=Dev->next;
                continue;
            }
            else if ( Dev->tratto<path )    // viaggia verso il dev...
            {
                path=Dev->tratto+0.01;
                if ( posok>0 )
                {
                    Dvk=Dev;
                    posok=pos;
                }
            }
        }

        else if ( pos==0 )               //Interno ..
        {
            // calcolo la dist dalle parete di uscita. o da un dev  interno
            Dev->tratto=Dev->CamInter();
            if ( Dev->tratto<path ) path=Dev->tratto+0.01;  // corretto
            Dvk=Dev;
            posok=pos;
//      if ( pr->idm==14 ) Dev->SpiaPr ( " Inter ", 48,48 );

            //qui cambiato??????????????????????????????????????
            //       break;  se non ci sono device interni si potrebbe uscire subito

        }
        else if ( pos==-1 )   //nel buco...  anxhe qui attenzione ai dev interni ......
        {
            Dev->tratto=Dev->CamBuco();
            if ( Dev->nowpos>-1 )
            {
                Dvk=Dev;
                posok=pos;
            }
            if ( Dev->tratto<path )  // corretto
            {
                path=Dev->tratto+0.01;
                if ( posok>-1 )
                {
                    Dvk=Dev;
                    posok=pos;
                }
            }

        }


        Dev=Dev->next;
    }       // end loop

    if ( Dvk!=0 )
    {


        Dvk->pr_dev=pr;   // ricordo l'indirizzo della particella in tracing..
        Dvk->camm= path+CEPSI;
        Dvk->prevpos=Dvk->nowpos;
        Dvk->nowpos=posok;
//  Dvk->SpiaPr ( " Scelt ", 2,4, Dvk->camm );
    }
    return Dvk;
}

void Device::SpiaPr ( char *Tit, int start,int stop, double corri )
{
    if ( idev<start|| idev>stop ) return;
    Gout.setf ( ios::fixed,ios::floatfield );
    Gout<<setprecision ( 2 );
    Gout<<"\n  Ev "<<setw ( 8 ) <<evento_.Gen.Event<<setw ( 10 ) <<Tit
    <<" " <<setw ( 10 ) << nome
    <<" at z "<<setw ( 10 ) <<Cface.z
    <<" pr "<<setw ( 6 ) <<pr_dev->nome
    <<"   at_Z "<<setw ( 10 ) <<X_dev.z
    <<" acc "<<setw ( 4 ) <<accedi
    <<" es "<<setw ( 4 ) <<escludi
    <<" Rv "<<setw ( 4 ) <<rivela
    <<" prev "<<setw ( 4 ) <<prevpos
    <<" pos "<<setw ( 4 ) <<nowpos
    <<" tratto "<<setw ( 10 ) <<corri;

}