flyopunta/src/src_par/particella.cpp

/***************************************************************************
                             -------------------
    begin                : Sun Nov 26 2000
    copyright            : (C) 2000 by Giuseppe m Pierazzini
    email                : pierazzini@unipi.it
 ***************************************************************************
 *                                                                         *
 *   NA48  simulation program.    1992----2005                                         *
 *                                                                         *
 ***************************************************************************/
#define UNDEF    6
#define ESTERNO  1
#define INTERNO  0
#define INBUCO   -1
#define NEUTRI  1
#define CARICHI 2
#define ALLCHR  3

#include "flyoh.h"
#include "devfluka.h"

using namespace GmpConst;
using namespace std;

//==============================================================
//          P A R T I C E L L A                  Construttor
//=====================================================================

Particella::Particella ()
{

    if ( ( pointer_db=Matter->Get_matter ( Partname, this )) == 0 )
    {
        Gout<<"\n Particella() Partname error....  "<<Partname<<std::endl;
        exit ( 0 );
    }

// Create a new reaction with links.
    SelRea->Totpart++;
    ido = SelRea->Totpart;

    if ( SelRea->Lastpart != 0 )   SelRea->Lastpart->next = this;
    SelRea->Lastpart = this;
    this->next = 0;
    id=0;    //vedi  Set_Id();
    fato = Undef;
    devhit=Ifato =devstory=0;
    if ( charg ==0 )   type = Neut;
    else  type = Char;

    if ( massa==0.0&&width==0.0&&ctau==0.0 ) type=Res;

    nata=0;
    rea_madre=LastReaction;
    Pair_Production=Brems_Production=0;    // 1 attiva la bremsstrhalung, 2 attiva la pair prodution...
    P.Riposo(massa);
    Qpatre=Pburst=Pbar=Gp=P;
    Gx.Reset();

//  Gout<<"\n Part P.mq %12.6lf %12.6lf %12.5lg %12.5lg",massa,P.m,P.mq,Gp.mq);

    pathok=path_max=path_done=last_path=0.0;
    cstar=10.;     // flag esplicitamente voluta...vedi Voldec...

//------------------links----
    up = 0;
    rg = 0;
    dw = 0;
    lf = 0;
    cls = 0;

}

//==============================================================
//                A g g i o r n a _ l i n k
//===============================================================
void Particella::
Aggiorna_link ( Particella * pup, Particella * prg,
                Particella * pdw, Particella * plf )
{
    if ( pup != 0 )     up = pup;
    if ( prg != 0 )     rg = prg;
    if ( pdw != 0 )     dw = pdw;
    if ( plf != 0 )     lf = plf;
}
void Particella::
Aggiusta_link ( Particella * pup, Particella * prg,
                Particella * pdw, Particella * plf )
{
    Aggiorna_link ( pup,prg,pdw,plf );
//restaura i links...
    if ( pup != 0 )  pup->dw=this;
    if ( prg != 0 )  prg->lf=this;
    if ( pdw != 0 )  pdw->up=this;
    if ( plf != 0 )  plf->rg=this;

}
void Particella::
Aggiorna_link_Bullet ( Particella * pbers )
{
    if ( pbers!= 0 )     rg = pbers;
}


Particella * Particella::
Decay_in_2 ( const char *pa_nome,const char *pb_nome )
{
// pattern == DECAY2 :padre=> pa,pb

//  Gout<<"\n Ev Dec1.out "<<Eventi_Fatti<<" "<<nome<<" ==> "
//       <<pa_nome<<" "<<pb_nome<<flush;


    Particella *padre,*pa,*pb;
    padre=this;
    pa = SelRea->NuovaPart ( pa_nome );

    pb = SelRea->NuovaPart ( pb_nome );



//  Aggiorna_link ( up, rg,dw,lf )
    padre->Aggiorna_link ( 0,      0,      pa, 0 );
    pa->Aggiorna_link ( padre,  pb, 0,      0 );
    pb->Aggiorna_link ( padre,  0,      0,      pa );
    pa->Set_Id();
    pb->Set_Id();
//    Gout<<"\n Ev Dec3.out "<<Eventi_Fatti<<" "<<nome<<" ==> "
//       <<pa->nome<<" "<<pb->nome<<" [ "<<pa<<" "<<pb<<" ] "<<flush;

    return pa;
}
//====================
Particella * Particella::
Decay_in_3 (const char *pa_nome, const char *pb_nome,const char *pc_nome  )
{
// three bodies decay....pattern == DECAY3;
//    padre -> pa+res2
//    res2   -> pb + pc;
//  but  the muons...

    Particella *padre,*pa,*pb;
//  const char *word;
    padre=this;

    {
        pa = padre->Decay_in_2 (pa_nome,Res2);
        pb = pa->rg->Decay_in_2 ( pb_nome,pc_nome );
    }
    return pa;
}

//============
Particella * Particella::
Decay_in_4 ( const char *pa_nome,  const char *pb_nome,const char *pc_nome,
             const char *pd_nome  )
{
    Particella *padre,*pa,*pb;
    padre=this;
//decadimento a 4    padre->  pa+ res3
//decadimento a 3    Res3 -> pb + pc + pd
//
    pa = padre->Decay_in_2 (pa_nome, Res3);
    pb=pa->rg->Decay_in_3(pb_nome,pc_nome,pd_nome);
    return pa;
}


//============
Particella * Particella::
Decay_in_5 ( const char *pa_nome,  const char *pb_nome,const char *pc_nome,
             const char *pd_nome ,const char *pe_nome )
{
    Particella *padre,*pa,*pb;
    padre=this;

//decadimento a 5    padre->  pa+Res4
//decadimento a 4     Res4 ->  pb + pc + pd + pe
    pa = padre->Decay_in_2 (pa_nome, Res4);
    pb=pa->rg->Decay_in_4 (pb_nome,pc_nome,pd_nome,pe_nome);
    return pa;
}

//
void  Particella::Save_Link()
{
    lfsav=lf;
    rgsav=rg;
    upsav=up;
    dwsav=dw;
    if ( next!=0 ) next->Save_Link();
}
void  Particella::Restore_Link()
{
    lf=lfsav;
    rg=rgsav;
    up=upsav;
    dw=dwsav;
    cls=0,nata=0;
    if ( next!=0 ) next->Restore_Link();
}
//==============================================================
//                U p d a t e _ D e v s t o r y
//===============================================================
// salva in formato codificato i devices colpiti dalla particella.
// il bit acceso corrisponde al device colpito
// attenzione memorizza solo i primi 64 device do offset in poi.
void Particella::
Update_Devstory(int idev) {
    static int offset=-1;
    if (offset==-1) {
        offset=0;
        if (collo!=0) offset=collo->idev;
        else if (xch02!=0) offset=xch02->idev;
        else return;
    }
    if (offset<1) return;
    unsigned long long bit_uno=1;
    int shift= idev-offset;
    if (shift>-1&&shift<64)
        devstory = (devstory|(bit_uno<<shift));
    return;
}
//================================================================
//              R e s e t _ stato
//=================================================================

void Particella::
Reset_Stato()
{
// start from the tree top.....


//  Gout<<"\n Reset "<<nome <<" id "<<id ;
    fato = Undef;
    nata=0;
    devhit=devstory=Ifato=0;
    P.Riposo();
    Gp.Riposo();
    X.Reset();
    Gx.Reset();
    if ( Pair_Production>0 )   Pair_Production=1;
    if ( Brems_Production>0 ) Brems_Production=1;
    pathok=path_max=path_done=last_path=0.0;
    cstar=10.;     // flag esplicitamente voluta...vedi Voldec...
    if ( next!=0 ) next->Reset_Stato(); //iterate attenzione meglio di no!!!!

}



//===============================================================
//                R i c e r c a
//==============================================================
Particella *Particella::
Ptr_Particle ( const char *pp )
//   attenzione va bene solo per MakeList...
// cerca la particella di nome pp tra le particelle della lista già create...
{
    Particella *ptr;
    ptr = this;
    while ( ptr != 0 )
    {
        if ( pp == ptr->nome && ptr->dw == 0 )
            return ptr;   // si!
        ptr = ptr->next;
    }
    return ptr;     // no!  return 0

}

//==============================================================
Particella *Particella::
Get_Pointer ( char *pp )
{
//   con procedura iterativa....attenzione va bene solo per MakeList...
// o la dove viene chiamata poche volte per run...
// altrimenti mangia spazio memoria......
    if ( pp==nome ) return this;
    if ( next!=0 ) return Get_Pointer ( pp );
    return 0;
}

//================================================================
//           S t a m p a    q u e s t a   p a r t i c e l l a
//=================================================================
void Particella::
Stampa_questa_particella ()
{
    // Gout<<"\n Ev "<< Eventi_Fatti<<"  Stampa part num: "<<id;
    if ( nome==Dlz && dw == 0 )
    {
        Gout<<"\n Dalitz pair no electrons linking.. stop"<<std::endl;
        exit ( 0 );
    }

//   P.print(nome);
    Gout<<"\n "<<setw ( 4 ) <<ido<<" <"<<setw ( 3 ) <<idm<<"> "<<setw ( 8 ) <<id<<"  "
    <<setw ( 5 ) <<nome<<" ["<<setw ( 6 ) <<fato<<" "<<Ifato<<"] Type: "<<setw ( 6 ) <<type;


    if ( dw != 0 )
    {
        Gout<<" ==>  dw-> "<<setw ( 8 ) <<dw->id<<" <"
        <<setw ( 3 ) <<dw->idm<<"> "
        <<setw ( 5 ) <<dw->nome<<" "<<setw ( 7 ) <<dw->fato<<" ";
        if ( dw->rg != 0 )
            Gout<<" + "<<setw ( 5 ) <<dw->rg->nome<<" "<<setw ( 7 ) <<dw->rg->fato<<" "
            <<setw ( 8 ) <<dw->rg->id<<" <"
            <<setw ( 3 ) <<dw->rg->idm<<">";
    }

    if ( rg!=0 && type==Bullet)
        Gout<<" ==> urta-> "<<setw ( 9 ) <<rg->id<<" <"
        <<setw ( 3 ) <<rg->idm<<"> "
        <<setw ( 4 ) <<rg->nome<<" "<<setw ( 7 ) <<rg->fato;

    if ( cls!=0 )
        Gout<<" ==>  Cls-> "<<setw ( 9 ) <<cls->id<<" <"
        <<setw ( 3 ) <<cls->idm<<"> "
        <<setw ( 4 ) <<cls->nome<<" "<<setw ( 7 ) <<cls->fato;
}

//================================================================
//           S t a m p a  p a r t i c e l l e
//=================================================================

void Particella::
Stampa_particle ( Reazione *rea )
{
    //iterative method  attenzione usare con parsimonia!!!!

    if ( this==rea->avo )
        Gout<<"\n\n    ido  idm     id     nom   stato          stampa Ev "<<evento_.Gen.Event;
    Stampa_questa_particella();
// iterate....
    if ( dw!=0 )
    {
        dw->Stampa_particle ( rea );
        if ( dw->rg!=0 ) dw->rg->Stampa_particle ( rea );
    }
    if ( cls!=0 ) cls->Stampa_particle ( rea );
    if (type==Bullet&&rg!=0) rg->Stampa_particle ( rea );

    return;
}




//================================================================

// redefine the particle id
void Particella::
Set_Id()
{
// The reaction is described as chain of two bodies
// decay.
// The ID is usefull to undertand the particle position.
//
// It is defined as: id = xxx.p.ll:
// ll= the decay order or the generation level.
// p = 0  is the right particle in the deacy,p = 1  to theleft, p = 3  ( cls)  particle belongin to a cluster (see gasint)
// xxx =  ido generation number of the mother pareticle
// so:  ID = ll+100*p+1000*xxx

    if ( up!=0 ) {
        id=(up->id%100)+1+1000*up->ido;

        if (lf!=0 )
        {
            if (type!=Bers)id+=100;
            else id+=400;
        }
        else if (up->cls==this )id+=300;

    }
    else id=0;

}

//=============================================================
//             D e b u g _ o u t
//==============================================================
// printout for events debug                      Debug_out
void Particella::Debug_out ()
{
    Particella *pr;
    Device *pv;
    pr = this;

    Gout<<"\n\n Ev "
    <<setprecision ( 2 ) <<fixed
    <<setw ( 7 ) << evento_.Gen.Event<<" time="
    <<setw ( 6 ) <<TBurst->Get_Tempo() <<" tl="
    <<setw ( 6 ) <<TBurst->Get_Tlast() <<" "
    <<setw ( 3 ) << pr->id<<" <"
    <<setw ( 3 ) <<pr->idm<<"> "
    <<setw ( 25 ) <<pr->rea_madre->Get_Titolo();

    TBurst->Reset_Tlast();
    while ( pr != 0 )
    {
        if ( pr->fato != Undef )
        {
            Gout<<"\n "
            <<setprecision ( 2 ) <<fixed<<" P=("
            <<setw ( 7 ) << pr->P[0]<<" "
            <<setw ( 7 ) << pr->P[1]<<" "
            <<setw ( 8 ) << pr->P[2]<<") E="
            <<setw ( 6 ) << pr->P[3]<<"   Fato: "
            <<setw ( 6 ) << pr->fato;
            Gout<<"\n Gp=("
            <<setw ( 7 ) << pr->Gp[0]<<" "
            <<setw ( 7 ) << pr->Gp[1]<<" "
            <<setw ( 8 ) << pr->Gp[2]<<") "<<pr->Get_Nome();
            Gout<<"\n Gx=("
            <<setw ( 7 ) << pr->Gx[0]<<" "
            <<setw ( 7 ) << pr->Gx[1]<<" "
            <<setw ( 8 ) << pr->Gx[2]<<")";
            Gout<<"\n x =("
            <<setw ( 7 ) << pr->X[0]<<" "
            <<setw ( 7 ) << pr->X[1]<<" "
            <<setw ( 8 ) << pr->X[2]<<")   Path_done: "
            <<setw ( 8 ) <<pr->path_done<<std::endl;
        }

//      if(pr->nome==Res) resmas=pr->massa;
        pr = pr->next;
    }
    double energia=0.0;
    pv = Apparato;
    while ( pv != 0 )
    {
        energia+=pv->eDev;
//           Gout<<"\n %s %8.3lf", pv->nome,pv->eDev);
        pv=pv->next;
    }
    Gout<<"\n Tot E seen ="<<energia;
    Gout<<"\n ----gmp-----   -------    ------\n"<<std::endl;
    return;
}

//              V o l d e c
//=============================================================
// Decay in flight of pa in pb and pc            Voldec
//-------------------------------------------------------------

int Particella::
Voldec ( Particella * pb, Particella * pc )

{
    Particella *pa;
    double eb, pcm,pt,pcmq;
    double rma, rmb, rmc;
    double cost,fi, sinfi, cosfi;
    static qvet Stesso;   //dummy for lorentz buster...
    //============

    pa = this;                // nota giusto per differenziare rispetto a pb e pc
    if ( pa ->fato==Decay );  // particelle che decadono
    else if ( pa->ctau < 0. )
        return -1;        // stable particle. no action is taken*/

    // Test the fato of the father particle, if not normal  */
    // stop the child generation */
// Gout<<"\n voldec "<<nome;


    if ( pa->fato != Nata && pa->fato != Decay )
    {
        pb->fato = Undef;
        pb->Ifato=0;
        pc->fato = Undef;
        pc->Ifato=0;
        pc->X = pb->X = pa->X;
        return -2;  //non va trasportato!
    }
// copia P nel vettore ausiliario Qpatre, può essere alterato da DecParm
    pa->Qpatre=pa->P;


//-----------------------------------------------------------
// produzione coppie nella interazione dei gamma nel materiale...
// produzione bremsstrhalung da leptoni...

    if ( pa->Pair_Production>1 )  pa->DecParm_Pair();
    if ( pa->Brems_Production>1 ) pa->DecParm_Brems();

    pb->fato = Nata;
    pb->Ifato = -1;
    pc->fato = Nata;
    pc->Ifato = -1;

//  get the parameters for particle : FaseSpace resonance == dec isotropici,
//  coppia dalitz, ... munu da K3l...
// does nothing for stable particles ....
    pb->DecParm();
    pc->DecParm();


// define the particles parameters for next use
    rma = Qpatre.m;
    rmb = pb->massa;
    rmc = pc->massa;

    pb->P[3] = eb = ( rma * rma + rmb * rmb - rmc * rmc ) * 0.5 / rma;
    pc->P[3] = rma - eb;
    pcm = sqrt ( ( pcmq= ( eb + rmb ) * ( eb - rmb ) ) );
    pb->P.normaq = pc->P.normaq = pcmq;
    pb->P.norma = pc->P.norma = pcm;

//******************************************************************/
// decdimento isotropo se cstar non e' stato gia' definito!
    if ( pb->cstar>1. )   // decadimento non definito: genero isotropico
    {
        cost = 1.0 -  2.0 * Pran();
        fi = DuePiGreco* Pran();
        cosfi = cos ( fi );
        sinfi = sin ( fi );
    }
    else  // fisica di decadimento definita
    {
        cost = pb->cstar;
        fi   = pb->fistar;
        cosfi = cos ( fi );
        sinfi = sin ( fi );
    }
// Gout<<"\n Voldec ev2 "<<Eventi_Fatti<<" "<<this<<" "<<pb <<" "<<pc<<flush;
//   Gout<<" "<<this->Get_Nome()<<" "<<pb->Get_Nome() <<" "<<pc->Get_Nome()<<flush;
    // Definisco un sistema con l'asse z nella direzione di moto
    // della particella padre nel baricentro dove fu creata.
    // Come creare l'asse x! adesso a random!!

    // Attenzione Pbar va definito cosi' se si tratta di risonanza
    // altrimenti occorre Pbar=pa->P vedere meglio!!!!!!!!!!!
    gvet vaz=pa->Pbar.Verso();          // asse z
// definizione del miglior asse sghembo a vaz.. arbitrario
    gvet vzn ( 1.,0.,0. );
// attenzione rivedere anche se funziona !!!!!!!!!!!!!
    double ax,ay,az;
    vaz.putv( ax,ay,az );
    ax=fabs ( ax );
    ay=fabs ( ay );
    az=fabs ( az );
    if ( ay<=az&&ay<=ax )   vzn.setv ( 0.,1.,0. );
    else if ( az<=ax )      vzn.setv ( 0.,0.,1. );
//-------
    gvet vay=vaz&vzn;
    vay.Norma();
    vay=vay.Verso(); // asse y
    gvet vax=vay&vaz;
    vax.Norma();
    vax=vax.Verso(); // asse x
    pt = pcm * ( sqrt ( 1. - cost * cost ) );// pericolo di negativa sqrt!
    gvet mom ( vax*pt*cosfi + vay*pt*sinfi + vaz*pcm*cost );
    //      mom.Norma();
    //    mom.print("mom");

    pc->P[0]=- ( pb->P[0]=mom[0] );
    pc->P[1]=- ( pb->P[1]=mom[1] );
    pc->P[2]=- ( pb->P[2]=mom[2] );



//ricordo il suo momento nel  baricentro.
    pb->Pbar=pb->P;
    pc->Pbar=pc->P;



    pb->Pburst=Qpatre;
    pc->Pburst=Qpatre;  // ricordo il momento del padre  ...
// vedi qui sotto.. si evitano problemi con le trasformate di lrentz   ..





//Get the lorentz booster



    if ( pa->type==Res )   // decadim in 3 or 4 bodies via Res....
//    pa->nome==Pic2 ,Pi02,P00p,Gg2
    {
// prima un burst nella direzione e velocità che aveva la Pa
// nel baricentro di suo padre
        pb->P=pb->P.Lburst ( pa->Pbar );
        pc->P=pc->P.Lburst ( Stesso );   //sono nel bar del padre della risonanza


// quindi passo al lab con  il burst impresso a Pa  a suo tempo
        pb->P=pb->P.Lburst ( pa->Pburst );
        pc->P=pc->P.Lburst ( Stesso );
    }

// se il padre non e' una risonanza passo direttamente al lab if pa moves    *

    else  if ( Qpatre.norma > 0.0 )
    {
        pb->P=pb->P.Lburst ( Qpatre );
        pc->P=pc->P.Lburst ( Stesso );


    }
//-------------------------------------



    pb->Gp = pb->P;
    pc->Gp = pc->P;

    pb->Vers=pb->P.Verso();
    pc->Vers=pc->P.Verso();
// nascita set to ok...
    pb->nata=pc->nata=1;
    // define the  initial coordinates at the father space point */
    // see the Transport  routine */
    pb->Gx=pb->X = pc->Gx= pc->X = pa->X;

    /*******************   debug   ***/
    if ( Debugon==1 )
//     if (Eventi_Fatti==129)
    {
//   pa->Stampa_questa_particella();

        pa->P.print ("Pa_p");
        pa->X.print ("Pa_X");
        Gout<<" "<<setw(8)<<pa->nome <<" "<<setw(4)<<pa->ido<<" "<<setw(4)<<pa->id  <<"  Ev "<<Eventi_Fatti;
        pb->P.print ("Pb_p");
        pb->X.print ("Pb_X");
        Gout<<" "<<setw(8)<<pb->nome <<" "<<setw(4)<<pb->ido<<" "<<setw(4)<<pb->id;
        pc->P.print ("Pc_p");
        pc->X.print ("Pc_X");
        Gout<<" "<<setw(8)<<pc->nome <<" "<<setw(4)<<pc->ido<<" "<<setw(4)<<pc->id;
        Gout<<endl;
    }

    /************   end debug ********/
    return 1;
}


//=============================================================
//              T R A S P O R T
//=============================================================
int Particella::
Trasport ()
{

    int botta = 0, loop = 1000;



// do nothing if .....
    if ( P.norma == 0.0 || ctau == 0.0 || ( fato != Nata ) ) return botta;


//     Gout<<"\n Entro in Transport "<<evento_.Gen.Event;
//     Gout<<"\n Trsa: "<<setprecision(2)
//     <<" Ev "<<evento_.Gen.Event<<" id "<< id<<" "<<fato<<" P "<<P.norma<<" ctau "<<ctau<<" see "<<see;



    if ( see == 0 ) {
        fato = Novede;
        Ifato=-10;
        return botta;
    }

    double lontano=Zona - X[2];

//   Gout<<"\nEv "<<Eventi_Fatti<<" Trans. " <<nome<<" "<< ido<< " " <<type;
//   if(dw!=0) Gout<< " "<<dw->nome;

    if ( type != Beam && type!=Bullet)
    {
        if ( ctau > 0.0 )   // particella che decade seguita da .....anche nulla..
        {
            double volo = ctau * P.norma / massa;
            path_max = -volo * log ( EPSI + Pran() );
            if ( path_max> lontano ) path_max = lontano;

//          if(Eventi_Fatti>71&Eventi_Fatti<91)
//           Gout<<"\nEv "<<setw(8)<<Eventi_Fatti<<" Trans. "<<setw(6)<<nome<<"  "<<setw(5)<< ido<<setprecision(2)
//           <<"  "<<setw(9)<< path_max<<" pz "<<setw(9)<<P.z<<" Rand "<<setw(9)<<Pran()<<flush;
        }
        else path_max = lontano;  // for stable particles...or with decay not defined...


    }

//    Gout<<"\nEv "<<Eventi_Fatti<<" Trans. "<<nome<<"  " << path_max<<" pz "<<P.z;
// set the devices accessible to the flighing particle..

    Set_Rivela_in_Dev ( this );
// Apparato->Set_Rivela(this);
//=====================================================
//=====================================================
// start  the loop on the apparatus

    // Vers = flight direction defined in particle object....

    path_done = pathok = 0.0;
    Device *pv,*Dvok;
    Devloop = Apparato;
    Apparato->Reset_Pos();

    while ( --loop >= 0 )
    {
        pv=Devloop;  // starting Dev... can be changed by Get_camm.
//      if(idm>4&&idm<9)
//        Gout<<"\n Ev "<<setw(5)<<Eventi_Fatti<< " "<<setw(8)<<nome << "   "<<setw(8)<<pv->nome;

        Dvok=pv->Get_camm ( this );

        if ( Dvok == 0 ) break; //end transport: passati tutti i dev...
        pathok=Dvok->camm;

//      if(Dvok->idev==2&&ido==1)
//       Gout<<"\n Ev "<<setw(5)<<Eventi_Fatti<< " "<<setw(8)<<nome << "  devloop "<<setw(8)<<pv->nome
//           << "  devok  "<<setw(8)<<Dvok->nome<<" " <<setw(8)<<Dvok->camm<< "   "<<Dvok->nowpos;


        /*
                if(Dvok->idev==2&&ido==1)
                  cout<<"\n Ev "<<setw(5)<<Eventi_Fatti<< " "<<setw(8)<<nome
                  << "   "<<setw(8)<<Dvok->nome<<" pth "<<setw(8)<<pathok
                 << " pos "<<setw(2)<<Dvok->nowpos<<" prev "<<setw(8)<<Dvok->prevpos;*/


        // potential path redefined if bigger then dacay path.
        if ( ctau>0.0 && pathok > ( path_max-path_done ) )
        {
            pathok= ( path_max-path_done+0.000001 );
            fato=Decay;
            Ifato=-2;


            /*
                                                if ( this==SelRea->avo && Dvok->fun==Gasint )
                      {
                        MkCluster(); fato=Gasint; Ifato=-20;   //Qui si crea un cluster.....
                        Gout<<"\n Ev0  "<<evento_.Gen.Event<<" "<<nome
                        <<" in "<<Dvok->nome
                        <<" "<<fato
                        <<" path "<<setw ( 8 ) << pathok
                        << " Z "<<setw ( 8 ) << X.z;

                      }
            */

        }





//==============================================================================
//==============================================================================

        if ( Dvok->nowpos == ESTERNO );  // do nothing


//==============================================================================
//==============================================================================


        else if ( Dvok->nowpos == INTERNO )

        {
        
            botta++;

//      if ( idm==14 ) Spia ( "Interno",pv,Dvok,47,52 );

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

// destino secondo la funzione

            if ( Dvok->fun == Trig    ||
                    Dvok->fun == None ||
                    Dvok->fun == Nodef ); // nessuna decisione .. si continua a tracciare

// altrimenti le cause per terminare il tracing...si cambia il fato

            else if ( Dvok->fun == Dead ||
                      Dvok->fun == Trigh ) {
                fato = Dvok->nome;
                Ifato=Dvok->idev;
            }

            if ( Dvok->rivela==2||   // assorbe tutta l'energia e muore
                    Dvok->rivela==4 )    {
                fato = Dvok->nome;
                Ifato=Dvok->idev;
            }

            else if ( Dvok->fun == Magn || Dvok->fun == Dump )
            {
                // nota il frame del magnete (che corrisponde all' INTERNO)  fa da dump. Il campo è sempre  nel "buco"
                if ( nome!=Mum && nome!=Mup ) //except the muons..
                {
                    fato = Dvok->nome;
                    Ifato=Dvok->idev;
                }
            }

            else if ( Dvok->fun == Veto )
            {
                // nel caso del "veto" si uccide addirittura l'evento con il flag botta =-1000
                // invece si termina  solo il tracing della particella se si vuole scrivere tutto
                if ( WrtNt!=WRTALL )
                {
                    fato = Dvok->nome;
                    Ifato=Dvok->idev;
                    botta=-1000;
                }
                else   {
                    fato = Dvok->nome;
                    Ifato=Dvok->idev;
                }
            }
            // debug  quale decisione?
            /*
            if ( Dvok->idev==2 )
            {
                X.print ( "X1" );
                Gout<<setprecision(3)
                    <<" ev "<< Eventi_Fatti<<" "<< nome
                    <<" pos " <<Dvok->nowpos<<" "<<Dvok->prevpos
                    <<" dh "<<devhit
                    <<"    ft "<<fato
                    << "    pathok " <<setw(9)<<pathok;
            }
             */
//----------------------------------------------------
//----------------------------------------------------
// scrivo i colpi (cioè memorizza gli hits nei devices)  quando

            if ( Dvok->prevpos==ESTERNO    ||   // è appena entrata dall'esterno
                    Dvok->prevpos==INBUCO  ||   // è entrata dal "buco"
                    Dvok->prevpos>=UNDEF   ||   // appena nata o devhit=0 (vedi ScrvColpi() )
                    devhit==0 )    Dvok->ScrvColpi ( this );                 // prima volta vista


        }
// ---------- fine interno -----------
//==============================================================================
//==============================================================================
//------------ E' nel buco
        else if ( Dvok->nowpos == INBUCO )
        {
            if ( Dvok->fun==Magn  &&          // entra nel "buco" del mag o nasce nel mag!
                    ( devhit==0           ||
                      Dvok->prevpos==UNDEF||
                      Dvok->prevpos==ESTERNO ) )  {
                Dvok->ScrvColpi ( this );
            }

            //pro memoria ...dall'INTERNO va  nel buco..per ora si lascia tracciare
            //if ( Dvok->prevpos== INTERNO );

        }
//--------Fine in buco------------
//==============================================================================
//==============================================================================
//.... check the fato

        if ( fato!=Nata&&fato!=Decay )
        {
            return botta;  //stop tracking....se il fato è diverso da normale o da "decay"
        }
// qui solo per particelle con fato normale o in decadimento.
//***************************************************************************
//----------------------------- real trasport  -------------------------------
//  now move the particle in the  space doing  also
//  pair production  or bremsstrhalung .or clusters..or move in magnetic field..and multiplescattering.
//


        if (Dvok->nowpos == ESTERNO)
        {
            // free flight

            Move ( pathok );
            last_path=pathok;
        }
        else
        {
            int ift=0;
            if ( ( ift=Dvok->Traccia ( this )) > 0 )  // muore nel tracing ...
            {
                fato=Dvok->nome;
                Ifato=Dvok->idev;
                return botta;
            }
            else if (ift < 0 )  // muore nel tracing ...
            {
                fato=Dvok->nome;
                Ifato=ift;
                return botta;
            }



        }

//==============================================================================
//==============================================================================
        if ( fato==Decay ) //stop tracking....anche per gam->pair prod and brems ...
        {
            if ( X.Norma() >Zona )  {
                fato=Decfar;
                Ifato=-8;
            }


// Si verifica se la particella beam  deve fare una interazione con il gas
// attenzione  i dev sono annidati....

            if ( this==SelRea->avo &&  Dvok->nowpos == INBUCO )  //Sono nel buco del dev
            {
// quindi guardo se sono anche nel buco del dev Dectubo ....
                if ( dectubo!=0 )
                {
                    if ( dectubo->Posizione() == INBUCO &&
                            dectubo->fun==Gasint )  // attenzione correzione per analisi fondi particolare...
                    {
// allora genero un cluster da   interazione in gas
                        MkCluster();
                        fato=Gasint;
                        Ifato=-20;   //Qui si crea un cluster.....
                    }
                }
            }
            return botta;
        }
        if ( P.norma<0.0001 )   // stop tracking  per fine energia
        {
            fato = Dvok->nome;
            Ifato=Dvok->idev;
            if ( Debugon==1 )
                Gout<<"\nTout2: Ev  "<<evento_.Gen.Event<<" "<<nome<<" "<<fato<<" "<<X.z;
            return botta;

        }
    }
//==============================================================================
//==============================================================================
//................  e n d    t r a c i n g   lo o p ........
//  here if fato == norm .....  rivedere qui....

    if ( ctau < 0.0 )  // particella stabile
    {
        if ( X.z<.90*Zona )
        {
            // X.Norma();X.print(nome);
            pathok=path_max-path_done;
            X+=Vers*pathok;  //transport
            // X.Norma();X.print(nome);
        }
        fato = Lost;
        Ifato=-9;
    }

// particella che decade......
    else if ( ctau > 0.0 && P.norma>0.001 )
    {
        pathok=path_max-path_done;
        X+=Vers*pathok;  //transport
        fato= Decay;
        Ifato=-2;

        if ( X.Norma() >Zona )
        {
            fato=Decfar;
            Ifato=-8;
        }
    }
    return ( botta );
}
//====================================
void Particella::Move ( double passo )
{
// move free along  remaining tratto
    X+=Vers*passo;
    X.Norma();
    path_done+=passo;
}

//================================================================
//                                   S p i a
//=================================================================

void Particella::Spia ( char * tit, Device *stdev, Device * loopdev,int start,int stop )
{
    if ( loopdev->idev<start|| loopdev->idev>stop ) return;
    Gout.setf ( ios::fixed,ios::floatfield );
    Gout<<setprecision ( 2 );
    Gout<<"\n Trans "<<tit<<" "<<evento_.Gen.Event
    <<" StartDev: "<<setw ( 10 ) <<stdev->nome
    <<" at z"<<setw ( 8 ) << stdev->Cface.z
    <<" "<<nome<<" at Z "<< X.z
    << " Path "<<setw ( 8 ) <<path_max
    << "  fatto "<<setw ( 8 ) <<path_done
    <<" in "<<setw ( 8 ) <<loopdev->nome
    <<" at z"<<setw ( 8 ) << loopdev->Cface.z
    <<" prev " <<setw ( 8 ) <<loopdev->prevpos
    <<" nowpos " <<setw ( 8 ) <<loopdev->nowpos
    << "   pathok "<<setw ( 8 ) <<pathok;

}