flyopunta/src/src_reaz/reazione.cpp

/***************************************************************************
                          greazio.cpp  -  description
                             -------------------
    begin                : Sun Nov 26 2000
    copyright            : (C) 2000 by Giuseppe m Pierazzini
    email                : pierazzini@unipi.it
 ***************************************************************************
 *                                                                         *
 *   NA48  simulation program.                                             *
 *                                                                         *
 ***************************************************************************/

// CP violatio       //Ultima revisione 2-7-99
#include "flyoh.h"

#define BEAM    0
#define DECAY2  2
#define DECAY3  3
#define DECAY4  4
#define DECAY5  5
#define URTO2   8
#define URTO3   9
#define NOTDEFINED -999999.

int pattern;

using namespace std;
//==============================================================
//               R E A Z I O N E     Constructor
//==============================================================
Reazione::Reazione ()

{

// Create a new Reaction with links.
  TotRea++;
  id = TotRea;
  if ( last_parz == 0 ) Reaction= this;
  else              last_parz->next = this;
  up = last_parz;
  LastReaction=last_parz = this;
  this->next = 0;
  SelRea=this;
// reset counters for   special Reactions ....


  ss_done=-1; // for thre bodies decay ..partfase ss calculation... -1 == not done; 1 == ok
  Pair_done=Totpair=0;
  Gasint_Done=0;
  Dalitz_done=Totdltz=0;
  Tailsrt_done=Tailsrt=0;
  Tailong_done=Tailong=0;
  Brem_done=Totbrem=0;
  Pipev=Totpart=0;
  avo =  Lastpart=0;
  header="";
  rea_type=rea_pattern=0;
// type =-1 beam, 0=decadimento (default) ,1=urto.  che accade nelle cl derivate
  generati = veri = scrnt=trasmessi=ubeam=0;
  branch=prodrate=rate = prob = 0.0;
  Dnsudp=maxspettro=0.0;
  inc_pinc=vxb=vyb=vzb=inc_teta=0.0;
  slope=wpinc = 0.0;
  rtarg = dztarg = 0.0;
  rbers = dzbers = 0.0;
  ctau_avo_su_mass=0.0;


  if ( up!=0 )
    {
      Ccoll=up->Ccoll;
      rcoll=up->rcoll;
      dzcoll=up->dzcoll;
    }

// controllo dello zero (bias) del sistema del lab. Tutti i rivelatori sono
// spostati di conseguenza..
  if ( Lab_Zero[0]==NOTDEFINED )
    {
      Gout<<"\n \n *** System lab origin NOT defined.. forced to zero!!!";
      Lab_Zero.setvn ( 0.,0.,0. );
    }


// legge le schede di definizione della reazione

  Leggi_schede_rea();

}
//=================   loop input cards ===

// get the input Data-list from file
void Reazione::Leggi_schede_rea()
{
  char line[156];
  string stringline,label;
  double a1,a2,a3,a4,a5,a6;


  Cardcur->getline ( line, 156 );

  while ( !Cardcur->eof() )

    {
      // Decodifico la linea
      stringline=line;
      stringstream scheda;
      scheda<<line;
      scheda>>label>>a1>>a2>>a3>>a4>>a5>>a6;

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

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

// Decodifico la linea

      if ( label.find ( "done" ) ==0||label.find ( "end" ) ==0 )
        {
          type=id*100 + avo->Get_Idm();
          if ( avo->Get_Massa() >0.0 ) ctau_avo_su_mass=avo->Get_Ctau() /avo->Get_Massa(); // reabeam...
          Gout<<"\n Done reazione <"<<id<<"> "<<header<<"   Br.R."<<rate<<endl;
          return;
        }

      if ( label.find ( "rate" ) ==0 ) branch=prodrate=rate = a1;
      else if ( label.find ( "slope" ) ==0 ) slope=a1;
      else if ( label.find ( "pinc" ) ==0||label.find ( "upinc" ) ==0 )
        {
          // beam define by user as gaussian shape (p,wpinc)
          // P is the mean value, wpinc is the sigma.
          if ( label.find ( "upinc" ) ==0 ) ubeam=1; // set flag for user defined beam
          inc_pinc=400.;
          inc_teta=0.;
          Pinc.setvn ( a1,a2,a3 );
          wpinc = a4;
          wpincq=wpinc*wpinc;
          if ( abs ( a5 ) >0.00001 ) inc_teta=a5;
          if ( a6>0.00001 ) inc_pinc=a6;
          Vinc_ppp.setvn ( sin ( inc_teta ), 0., cos ( inc_teta ) );
          Pinc_ppp=Vinc_ppp*inc_pinc;
        }
      else if ( label.find ( "target" ) ==0 )
        {
          Ctarg.setvn ( a1,a2,a3 );
          Ctarg=Ctarg-Lab_Zero;
          Ctarg.Norma();
          rtarg =a4;
          dztarg =a5;
        }
      else if ( label.find ( "coll" ) ==0 )
        {
          Ccoll.setvn ( a1,a2,a3 );
          Ccoll=Ccoll-Lab_Zero;
          Ccoll.Norma();
          rcoll  = a4;
          dzcoll = a5;
        }
      else if ( label.find ( "bersag" ) ==0 )
        {
          Cbers.setvn ( a1,a2,a3 );
          Cbers=Cbers-Lab_Zero;
          Cbers.Norma();
          rbers = a4;
          dzbers= a5;
        }
      else if ( label.find ( "pbers" ) == 0 ) //particella inc da destra
        {
          Pbers.setvn ( a1,a2,a3 );
          wpbers = a4;
          drbers = a5;
        }


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

          header=stringline.substr ( label.size(),stringline.size() );
          Gout<<"\n get header: "<< header;
          Makelist ();

          Particella *pp, *pr=avo;
// and Id setting
          while ( pr!=0 )
            {
              pr->Set_Id();
              pr=pr->next;
            }
// define rea_patern  only beam=0, n= n bodies decay
          if ( avo->dw!=0 )
            {
              int nb=2;
              pp=avo->dw;
              if ( pp->Get_Type() ==Res ) nb++;
              if ( ( pp=avo->dw->rg ) !=0 &&
                   pp->Get_Type() ==Res ) nb++;
              if ( ( pp=avo->dw->lf ) !=0 &&
                   pp->Get_Type() ==Res ) nb++;
              rea_pattern=nb;
            }
// compatibility test for four bodies decay
          if ( rea_pattern==4 )
            {
              if ( avo->dw->Get_Nome() !=avo->dw->rg->Get_Nome() )
                {
                  Gout<<"\n Reaction list error:...<"<<avo->dw->Get_Nome() <<" "<<avo->dw->rg->Get_Nome();
                  Gout<<"\nIn four body decay the two res. must have the same names!"<<std::endl;
                  exit ( 0 );
                }
            }
//    save the particle describing the fix part of the reaction
//    for bremsstrhalung /Paiproduction or see...MkCluster....
          avo->Save_Link();
          LastpartSave=Lastpart;
          TotpartSave=Totpart;
        }


      else
        {
          Gout<<"\n Bad line:  "<<line;
          exit ( 0 );
        }

      Cardcur->getline ( line, 156 );
    }
  Gout<<"\n Warnign ... Reaction  bad ending!!"<<std::endl;

}


//=====================================================================
//                    M a k e l i s t
//====================================================================
void Reazione::
Makelist ()
{
  char line[156];
  string label,a0,a1,a2,a3,a4,a5;



  // possibili patterns di dati in input:
  //      lb  0  1  2  3  4  5     indice di pattern
  // 0)   pa                     solo beam
  // 2)   pa  =  pb pc           decadimento in due corpi
  // 3)   pa  =  pb pc pd        decadimento in tre corpi
  // 4)   pa  =  pb pc pd pe     decadimento in quattro corpi
  // 5)   pa  =  pb pc pd pe pf   decadimento in cinque  corpi
  // 8)   pa pb   = pc pd        urto in due corpi (solo nella prima linea)
  // 9)   pa pb   = pc pd pe     urto in tre corpi (solo nella prima linea)

  Cardcur->getline ( line, 156 );         // da rivedere....!!!!!!!!!!
  while ( !Cardcur->eof() )
    {
      // Decodifico la linea


      a0=a1=a2=a3=a4=a5="";
      stringstream scheda;
      scheda<<line;
      scheda>>label>>a0>>a1>>a2>>a3>>a4>>a5;
                        
                

      if ( Dbginput==1 ) Gout<<"\n\n Mklist reads:  "<<line;

      if ( label.size()     ==0  ||
           label.find ( "**" )  ==0  ||
           label.find ( "\\\\" ) ==0  ||
           label.find ( "//" )  ==0 )
        {
          Cardcur->getline ( line, 156 );   //line empty or comment ...read next line
          continue;
        }
        
      if ( ( a0.size() !=0 &&a1.size() !=0) && ( a0!= "="&& a1 != "=" ) )


        {  // check the = position....
          cout<<"\n Input card error ...: "<<line<<"  ......exit";
          Gout<<"\n Input card error ...: "<<line<<"  ......exit  "
          << a0<< "  "<<a0.size() << "  " << a1<< "  "<<a1.size() << "  " << a2<< "  "<<a2.size();
          exit ( 0 );
        }

      if ( label.find ( "fine" )  ==0  ||
           label.find ( "end" )  ==0  ||
           label.find ( "done" )  ==0 )
        {
          if ( Dbginput==1 ) Gout<<"\n "<<label;
          return;
        }

// extract particles
      const char *p0=0,*p1=0,*p2=0,*p3=0,*p4=0,*p5=0;

// beam particle
      p0= Matter->Ptbsdat ( label.c_str() );
      if ( p0==0 )
        {
          cout<<"\n Reaction  card error ...:  "<<line<<"  ......exit";
          Gout<<"\n Reaction  card error ...:  "<<line<<"  ......exit";
          exit ( 0 );
        }

      pattern=-1;
// è solo un beam
      if ( a0.size() ==0 ) pattern = BEAM;

      else if ( a0 == "=" ) // è un decadimento
        {
          p1 = Matter->Ptbsdat ( a1.c_str() );
          p2 = Matter->Ptbsdat ( a2.c_str() );
          p3 = Matter->Ptbsdat ( a3.c_str() );
          p4 = Matter->Ptbsdat ( a4.c_str() );
          p5 = Matter->Ptbsdat ( a5.c_str() );
          if ( p1!=0 && p2!=0 ) pattern = DECAY2;
          if ( p1!=0 && p2!=0 && p3 !=0 ) pattern = DECAY3;
          if ( p1!=0 && p2!=0 && p3 !=0 && p4!=0 ) pattern = DECAY4;
          if ( p1!=0 && p2!=0 && p3 !=0 && p4!=0 && p5!=0 ) pattern = DECAY5;


//        Gout<<"\n "<<a0<<" "<<a1<<" "<<a2<<" "<<a3<<" "<<a4<<" "<<a5<<" "<<a5<<pattern<<endl;
        }

      else if ( a1 == "=" ) // è un urto
        {
          p1 = Matter->Ptbsdat ( a0.c_str() );
          p2 = Matter->Ptbsdat ( a2.c_str() );
          p3 = Matter->Ptbsdat ( a3.c_str() );
          p4 = Matter->Ptbsdat ( a4.c_str() );
          if ( p1!=0 && p2!=0 && p3 !=0 ) pattern = URTO2;
          if ( p1!=0 && p2!=0 && p3 !=0 && p4!=0 ) pattern = URTO3;

        }

      if ( Dbginput==1 ) Gout<<"\n Pattern found = "<<pattern;
      Gout.flush();

      if ( pattern == -1 )
        {
          Gout << "\n MakeLista: Pattern errore 2 ... exit";
          exit ( 0 );
        }

      if ( ( pattern == 0 ) && avo != 0 )
        {
          Gout << "\n MakeLista: avo errore 0,4,5... exit";
          exit ( 0 );
        }
// definisco l'avo della lista se non gia' fatto.
// altrimenti definisco il pontatore padre della decadimento...
      Particella *padre, *madre, *ccms;

      if ( avo == 0 )
        { 
          if ( pattern==URTO2||pattern==URTO3 )
            {
                rea_type=1; //urto
              avo = padre = NuovaBullet ( p0 );   //Type == Bullet
            }
          else //decay o beam
            {
              rea_type=0; 
              avo = padre = NuovaPart ( p0 );
              avo->Set_Type ( Beam ); // è la prima particella della lista... è il fascio.
            }

          if ( pattern==BEAM )
            { 
              rea_type=-1;
              Cardcur->getline ( line, 156 );
              continue;
            }
        }
      else padre = avo->Ptr_Particle ( p0 );
      if ( padre == 0 )
        {
          Gout << "\n MakeLista: Wrong name or chaining error < "<<p0<<" > ....exit. "<<endl;
          exit ( 0 );
        }

//======================================
//======================================
      if ( pattern == BEAM );
    
      else if ( pattern == DECAY2 ) padre->Decay_in_2 ( p1, p2 );
      else if ( pattern == DECAY3 ) padre->Decay_in_3 ( p1,p2,p3 );
      else if ( pattern == DECAY4 ) padre->Decay_in_4 ( p1,p2,p3,p4 ); 
      else if ( pattern == DECAY5 ) padre->Decay_in_5 ( p1,p2,p3,p4,p5 );
      
     
      
      //**************************************************************************
//    Sezione urti da rivedere
      else if ( pattern == URTO2 ) // two bodies final status // da controllare....!!!!!
        {
          madre = NuovaPart ( p1 );
          madre->Set_Type ( Bers );
          padre->Aggiorna_link_Bullet ( madre );
// attenzione  Cms deve essere definita  meglio.... Anzi va ridefinita ad ogni interazione...
          ccms = NuovaPart ( Cms );
          ccms->Set_Type ( Bari );
          if ( Dbginput==1 ) Gout<<"\n Define dummy Res: "<<Cms;
          padre->Aggiorna_link ( 0, madre, 0, 0 );
          madre->Aggiorna_link ( 0, 0,ccms, padre );
          ccms->Decay_in_2 ( p2, p3 );
        }
      else if ( pattern == URTO3 ) // three bodies final status

        {
          madre = NuovaPart ( p1 );
          madre->Set_Type ( Bers );
          padre->Aggiorna_link_Bullet ( madre );
          ccms = NuovaPart ( Cms );
          ccms->Set_Type ( Bari );
          if ( Dbginput==1 ) Gout<<"\n Define dummy Res: "<<Cms;
          padre->Aggiorna_link ( 0, madre, 0, 0 );
          madre->Aggiorna_link ( 0, 0,ccms, padre );
          ccms->Decay_in_3 ( p2,p3,p4 );
        }

      Cardcur->getline ( line, 156 );
    }
}


//                    N u o v a p a r t
Particella *Reazione::NuovaPart ( const  char *pnome )
{

  Partname=pnome;
  Particella *punta=0;

  if ( pnome==Res )         punta = new PartFase();     // due corpi solo spazio fasi
  else if ( pnome==Res2 )   punta = new PartFase();    // per tre corpi solo spazio fasi
  else if ( pnome==Res3 )   punta = new Partres4;   // per il decadimento a 4
  else if ( pnome==Res4 )   punta = new Partres5;   // per il decadimento a 5 corpi secondo lo spazio delle fasi...
  else if ( pnome==Pic2 )   punta = new Part3pic ( slope );   // decadimento K in 3 pic
  else if ( pnome==Pi02 )   punta = new Part3pi0 ( slope );  // dec in tre pi0
  else if ( pnome==P00p )   punta = new Partp00p;      //decadimento in pi0 pi0 pic
  else if ( pnome==Pvv )    punta = new PartPvv;  //Risonanza nu nu per decadimento Kl->pi0 nu nu
  else if ( pnome==Dlz )    punta = new PartDltz;   // dalitz
  else if ( pnome==Munu )   punta = new PartMunu;   // munu in K3l
  else if ( Abilita_Brems>0  &&
            ( pnome==Mum || pnome==Mup||
              pnome==Elec|| pnome==Elep ) )  punta = new PartLepton(); // tiene conto del bremsstrhalung
  else if ( Abilita_Pair>0 && pnome==Gam )  punta = new PartGam(); //with pair production..

  else    punta=  new Particella();

  if ( Dbginput==1 ) Gout<<"\n Creato " << pnome;


  return punta;


}
Particella *Reazione::NuovaBullet ( const  char *pnome )
{
  Partname=pnome;
  Particella *punta=0;
  punta=  new PartBullet();
  if ( Dbginput==1 ) Gout<<"\n Creato bullet " << pnome;
  return punta;
}



//===============================================
void Reazione:: AddGen()
{
  generati++;
//   remember the last event type...
  evento_.Gen.lastype=evento_.Gen.type;
  evento_.Gen.type =type;
}

//===============================================================
//                 S t a m p a _ r e a z i o n i
//===============================================================
void Reazione::
Stampa ()
{
//  using namespace std;
  Reazione *pts;
  pts = Reaction;
  while ( pts != 0 )
    {
      Gout<<"\n\n[Stampa_reazioni] \n  Reaction n. "<< setw ( 5 ) <<pts->id<<" [ "<<pts->header<<" ]";
      Gout<<"\n  Reaction pattern type:         "<<pts->rea_pattern;
      Gout<<"\n  Interaction parameter (slope): "<<pts->slope;
      Gout<<"\n  Tipo di reazione: "<< pts->rea_type
      <<"\n  Beam:   "<< ( pts->ubeam?"Gaussian profile by user":"Standard" );
      Gout<<"\n  Tot.rate = "
      <<setprecision ( 4 ) <<fixed
      <<setw ( 8 ) << pts->rate;
      Gout<<"\n  Prod target "
      <<setprecision ( 2 ) <<fixed
      <<setw ( 8 ) <<pts->Ctarg[0]<<" "
      <<setw ( 8 ) <<pts->Ctarg[1]<<" "
      <<setw ( 8 ) <<pts->Ctarg[2]<<" r="
      <<setw ( 8 ) << pts->rtarg<<" dz="
      <<setw ( 8 ) << pts->dztarg;

      pts->avo->Stampa_particle ( this );
      pts = pts->next;
    }

  Gout<<endl;
}
//=================================================================
void Reazione::
PrintPinc ()
{
  Gout<<" P =("
  <<setprecision ( 3 ) <<fixed
  <<setw ( 8 ) <<Pinc[0]<<" "
  <<setw ( 8 ) <<Pinc[1]<<" "
  <<setw ( 8 ) <<Pinc[2]<<") dp="
  <<setw ( 8 ) <<wpinc<<endl;
  Gout.setf ( ios::floatfield );
}

void Reazione::
PrintPbers ()
{
  Gout<<" P =("
  <<setprecision ( 3 ) <<fixed
  <<setw ( 8 ) <<Pbers[0]<<" "
  <<setw ( 8 ) <<Pbers[1]<<" "
  <<setw ( 8 ) <<Pbers[2]<<") dp="
  <<setw ( 8 ) <<wpbers<<" drad="
  <<setw ( 8 ) <<drbers<<endl;
  Gout.setf ( ios::floatfield );
}
//===============================================================
//                 F i n e _ r u n
//===============================================================
void Reazione::
Fine_run ()
{
  int gen=0,totnt=0,Vtot=0,VTtot=0,Tpair=0,Tlsrt=0,Tlong=0, Totpipev=0 ;
  int tgasdone=0;
  int Tdlz=0,Tbrem=0;
  Reazione *pts;
// calcolo dei protoni equivalenti su targhetta.
// vedi nota 01/GmP  Pisa 19.6.99
  double totime=TBurst->Get_TTempo();
  Gout<<"\n\n ----------   F i n e _ R u n "<<evento_.Gen.Run<<"  -----------\n";
  Gout.setf ( ios::floatfield , ios::scientific );
  Gout<<setprecision ( 5 )
  <<"\n\nTot PPP "<<Protoni<<"   in "<< totime <<" ns"<<endl;



  Gout.setf ( ios::left, ios::adjustfield );
  pts = Reaction;
  Gout<<"\nBeam  type Id   Gen       Gas      Dlz   Pair  Brem  Hsr     Hlg   Piped   Writ       Titolo";

  while ( pts != 0 )
    {
      Gout<<"\n <"
      <<setw ( 4 ) <<pts->avo->Get_Nome() <<" "
      <<setw ( 4 ) <<pts->type<<" "
      <<setw ( 3 ) <<pts->id<<"> "
      <<setw ( 8 ) <<pts->generati<<"   "
      <<setw ( 8 ) <<pts->Gasint_Done<<" "
      <<setw ( 5 ) <<pts->Totdltz<<" "
      <<setw ( 5 ) <<pts->Totpair<<" "
      <<setw ( 5 ) <<pts->Totbrem<<" "
      <<setw ( 5 ) <<pts->Tailsrt<<" "
      <<setw ( 5 ) <<pts->Tailong<<" "
      <<setw ( 7 ) <<pts->Pipev<<" "
      <<setw ( 7 ) <<pts->scrnt<<" "
      <<setw ( 25 ) <<pts->header;

      gen+=   pts->generati;
      tgasdone+=pts->Gasint_Done;
      totnt+= pts->scrnt;
      Vtot+=  pts->veri;
      VTtot+= pts->trasmessi;
      Tlsrt+= pts->Tailsrt;
      Tlong+= pts->Tailong;
      Tpair+= pts->Totpair;
      Tdlz+=  pts->Totdltz;
      Tbrem+= pts->Totbrem;
      Totpipev+= pts->Pipev;
      pts =   pts->next;
    }

  Gout<<"\n  T o t a l i :  "
  <<setw ( 8 ) <<gen<<"   "
  <<setw ( 8 ) <<tgasdone<<" "
  <<setw ( 5 ) <<Tdlz<<" "
  <<setw ( 5 ) <<Tpair<<" "
  <<setw ( 5 ) <<Tbrem<<" "
  <<setw ( 5 ) <<Tlsrt<<" "
  <<setw ( 5 ) <<Tlong<<" "
  <<setw ( 7 ) <<Totpipev<<" "
  <<setw ( 7 ) <<totnt<<" "
  <<" All Reactions";
  Gout<<"\n\nTau per ppp produced   "<<totime/Protoni<<" Tau per evento written "<< totime/totnt <<endl;
  Gout<<"\n\nLast event ="<<setw ( 8 ) <<evento_.Gen.Event<<" Normal end run....."<<std::endl;
  return;
}

//===============================================================
//                 S e l e c t _ r e a t i o n
//===============================================================
// Select a Reaction according to the rate.
Reazione *Reazione::
Select_reation ()
{
  Reazione *pts;
  double prb = Pran();

  pts = Reaction;
  while ( pts != 0 )
    {
      if ( prb <= pts->prob )
        {
          return ( pts );
        }
      pts = pts->next;
    }

  return ( pts );
}
//==========================================================================
//                 N o r m a _ R a t e
//===============================================================
// Normalizza i rate di decadimento.
void Reazione::
Norma_Rate()
{
  Reazione *pts;

  double tprob=0.0,sigmat=0.0;
  Gout<<"\n\n  Rates normalization....";
  pts = Reaction;
  while ( pts != 0 )
    {
      sigmat+=pts->prodrate;
      pts = pts->next;
    }
  pts=Reaction;
        
        // rinormnalizzazione della sezione d'urto totale.....
         // sigmat*=10.;
        
  while ( pts != 0 )
    {

      // calcolo del rate di produzione normalizzato
      pts->rate=pts->prodrate/sigmat;
      tprob+=pts->rate;
      pts->prob=tprob;  // probabilita' integrale... per la generazione.





      Gout.setf ( ios::floatfield , ios::scientific );

      Gout<<"\n Rea: "
      <<setprecision ( 3 ) <<setw ( 2 ) <<pts->id
      <<"   Ev/PPP " <<setw ( 11 ) <<pts->prodrate
      <<"   Normalizzato "<<setw ( 8 ) <<pts->rate
      <<"  Integrale " <<setw ( 8 ) <<pts->prob
      <<"   " <<setw ( 8 ) <<pts->header;
      pts = pts->next;

    }
  //  fattore di normalizzazione della sezione d'urto.....
  double Aggiusta_Protoni=0.05;
  Protoni=Aggiusta_Protoni*evento_.Mxev/sigmat;
  double tgen=TBurst->Get_Tburst() * ( Protoni/ppp );
  double tmedio=tgen/evento_.Mxev;
  Gout<<"\n\n For all reactions: "
  <<setprecision ( 5 ) <<setw ( 2 )
  <<" Ev/PPP " <<setw ( 11 ) <<sigmat
  <<"\n Needed Protons  "<<setw ( 8 ) <<Protoni
  <<"   for  "<<setw ( 8 ) <<evento_.Mxev<<" Events "
  <<"\n Production time   "<<setw ( 6 ) << tgen<<" ns  with Tau "
  <<setw ( 6 ) << tmedio<<" ns between two events.   "
        <<setw ( 6 ) << 1000000./tmedio<<" kfreq  Eventi secondo\n-----------------------";

//separazione temporale in ns tra due protoni no sterili.....

  TBurst->Set_Tmed ( tmedio );

}
//================================================================
//              R e s e t _ R e a z i  o n e
//=================================================================
void Reazione::
Reset_Rea()
{

  Particella *pr= LastpartSave->next,*questa=0;
  Lastpart->next=0;
  Lastpart= LastpartSave;

  LastpartSave->next=0;
  Totpart=TotpartSave;
  avo->Restore_Link();
  avo->Reset_Stato();

// free the memory.... or destroy the partucle....
  while ( pr!=0 )
    {
      questa=pr;
      pr=pr->next;
      delete questa;  // here delete the particle
    }

  ss_done=-1; //ss calculation in partfase for three bodies...
  Pair_done=0;
  Brem_done=0;
  Dalitz_done=0;
  Tailsrt_done=0;
  Tailong_done=0;
  mbrem=mdltz=mspair=eeb=eleb=edltz1=edltz2=0.0;
  GasIdo=0;
}
//--------------------------------------------------------------*/
//               Particle simple list                           */
//--------------------------------------------------------------*/

void Reazione::
Simple_List()
{
  // itaration method ...do not use this procedura in the loop....
  Particella *loopx;


  Gout<<"\n \nEV  "<<Eventi_Fatti<<"  Reaction particle list: "<< header;
  loopx = avo;
  Gout
  <<"\n     n        id   idm  nome  type    ch       mas    width    t_life          <  G.Z       Gp   >    <   Dead at: X           >      path      fato    Ifato";

  while ( loopx != 0 )
    {
      Gout<<"\n"
      <<setprecision ( 2 ) <<fixed
      <<setw ( 8 ) << loopx->Get_Ido() <<" "
      <<setw ( 8 ) << loopx->Get_Id() <<" "
      <<setw ( 3 ) << loopx->Get_Idm() <<" "
      <<setw ( 6 ) << loopx->Get_Nome() <<" "
      <<setw ( 6 ) << loopx->Get_Type() <<" "
      <<setw ( 6 ) << loopx->Get_Charg() <<" "
      <<setw ( 8 ) << loopx->Get_Massa() <<" "
      <<setw ( 8 ) << loopx->Get_Width() <<" "
      <<setw ( 8 ) << loopx->Get_Ctau() <<"      "
      <<setw ( 9 ) << loopx->Get_Gx().z <<" "
      <<setw ( 9 ) << loopx->Get_Gp().norma<<"     "
      <<setw ( 9 ) << loopx->X.x <<" "
      <<setw ( 9 ) << loopx->X.y <<" "
      <<setw ( 9 ) << loopx->X.z <<"   "
      <<setw ( 9 ) << loopx->path_done<<"   "
      <<setw ( 8 ) << loopx->Get_Fato() <<" "
      <<setw ( 6 ) << loopx->Get_Ifato();
      loopx = loopx->next;
    }
  Gout<<endl;
  if ( next!=0 ) next->Simple_List(); //iterate... attention
}