flyopunta/src/src_dev62/devlkr.cpp

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

#include "parm.h"
#include "part_db.h"
#include "devlkr.h"
#include "particella.h"
#include "flyorea.h"
#include "evento.h"
#include "triggerbox.h"



#define INTERNO  0
#define EPSR     1.0E-200
using namespace GmpConst;
using namespace std;


PadRaster PAD;

double soglia_pic=0.35; //soglia di .3 gen per selezionare i picchi nelle pad44

//==========================================
// -------  CELLA ------------------

void Cella::Set_index ( double xx ,double yy, int ixx,int iyy )
{
  double Dcell=2.;

  if ( ix>-1&&ix<128&&iy>-1&&iy<128 )
    {
      x=xx;
      y=yy;
      ix=ixx;
      iy=iyy;
      id=128*iy+ix;
      Xcella = ( ix+0.5 ) *Dcell-128;
      Ycella = ( iy+0.5 ) *Dcell-128;
    }
// else pone  id =0;che  corrisponde ad una regione esterna al volume di fiducial del calorimetro
  else
    {
      id=0;
      Xcella=Ycella=0.0;
      ix=0;
      iy=0;
      e=0.;
    }

};


// ============
void Cella::Rilascia_E ( double sig,double e_in )
{
  if ( id<1 ) return;
  double d=1.;  // mezza cella

// if the center  is in the hole... get out....
  if ( ( Xcella*Xcella+Ycella*Ycella ) <64. ) return;


  double u1= ( x+d-Xcella ) / ( sig*SQRT2 );
  double u2= ( x-d-Xcella ) / ( sig*SQRT2 );
  double v1= ( y+d-Ycella ) / ( sig*SQRT2 );
  double v2= ( y-d-Ycella ) / ( sig*SQRT2 );

  double erfx=fabs ( erf ( u2 )-erf ( u1 ) );
  double erfy=fabs ( erf ( v2 )-erf ( v1 ) );
  double erfxy=erfx*erfy;
  if ( erfxy<0.00001 ) erfxy=0.;

  e =0.25*e_in*erfxy;

//    cout<<"\n erfxy  "<<setw(8)<<u1<<" "<<setw(8)<<u2<<" " <<setw(8)<<v1<<" " <<setw(8)<<v2
//    <<" " <<setw(8)<<erfxy<<" " <<setw(8)<<e<<setw(8)<<e_in;

};
void Cella::Cella2Pad()
{
  if ( id<1 ) return;

  int ixp=int ( ix/2 );
  int iyp=int ( iy/8 );
  PAD.padx28[iyp][ixp]+=e;   // rivedre gli indici

  // pads di 4x4
  ixp=int ( ix/4 );
  iyp=int ( iy/4 );
  PAD.padxy44[ixp][iyp]+=e;
//  cout<<" xye "<<ixp<<" " <<iyp<<" "  <<e <<" " << PAD.padxy44[iyp][ixp]  ;

}

//=========================PadRaster==========
void PadRaster::Reset()
{
  epad28=m28x=m28qx=m28y=m28qy=0.0;
  epad44=m44x=m44qx=m44y=m44qy=0.0;
  dist_max=0.;
  padhit=0;

  for ( int i=0;i<1024;i++ )
    {
      pad28[i]=pad44[i]=0.;
      PADflag_l[i]=0;

    }
}

//======

void PadRaster::Momenti()
{

  // lavora sulle proiezioni..... Padhit deve essere definito come mhit inizialmente..
  if ( padhit==0 ) return;

  for ( int k=0;k<1024;k++ )
    {
      epad28+=pad28[k];
    }
  if ( epad28<0.001 ) return;
  epad44=epad28;

// 28 proiezione x

  for ( int i=0;i<64;i++ )
    {
      double ri= ( i+0.5 ) *4.-128;
      double riq=ri*ri;
      for ( int j=0;j<16;j++ )
        {
          m28x +=padx28[i][j]*ri;
          m28qx+=padx28[i][j]*riq;
        }
    }
// 28 proiezione y

  for ( int i=0;i<16;i++ )
    {
      double ri= ( i+0.5 ) *16.-128;
      double riq=ri*ri;
      for ( int j=0;j<64;j++ )
        {
          m28y +=padx28[j][i]*ri;
          m28qy+=padx28[j][i]*riq;
        }
    }

// proiezioni x e y

  for ( int i=0;i<32;i++ )
    {
      double ri= ( i+0.5 ) *8.-128;
      double riq=ri*ri;
      for ( int j=0;j<32;j++ )
        {
          m44y +=padxy44[j][i]*ri;
          m44qy+=padxy44[j][i]*riq;

          m44x +=padxy44[i][j]*ri;
          m44qx+=padxy44[i][j]*riq;


        }
    }
//mediepadxy44

  m28x/=epad28;
  m28qx/=epad28;
  m28qx-=m28x*m28x;
  m28qx=sqrt ( m28qx );
  m28y/=epad28;
  m28qy/=epad28;
  m28qy-=m28y*m28y;
  m28qy=sqrt ( m28qy );
  m44x/=epad44;
  m44qx/=epad44;
  m44qx-=m44x*m44x;
  m44qx=sqrt ( m44qx );
  m44y/=epad44;
  m44qy/=epad44;
  m44qy-=m44y*m44y;
  m44qy=sqrt ( m44qy );

//quadranti

  quad1=0.0;
  quad2=0.0;
  quad3=0.0;
  quad4=0.0;

  for ( int ix=0;ix<16;ix++ )
    {
      for ( int iy=0;iy<16;iy++ )
        {
          quad1+=padxy44[ix][iy];
          quad2+=padxy44[ix+16][iy];
          quad3+=padxy44[ix+16][iy+16];
          quad4+=padxy44[ix][iy+16];
        }
    }
}
//
void  PadRaster::Get_Peak()
{
  npeak=0;

// main loop to discover peaks: avoid the border channels.
  for ( int ix=1;ix<31;ix++ )
    {
      for ( int iy=1;iy<31;iy++ )
        {
          if ( Peak_mask ( ix,iy ) ==1 )
            {

              Peak[npeak].ix=ix*8-124;
              Peak[npeak].iy=iy*8-124;
              npeak++;
              if ( npeak>9 ) break; //too many peaks

            }
        }
    }

  dist_max=150.;
  if ( npeak>8 ) return;//too many peaks
  dist_max=-1.;
  if ( npeak<2 ) return;  // only one or no peak;

// check max distance between peaks.

  dist_max=0;
  for ( int i=0;i<npeak-1;i++ )
    {
      for ( int j=i+1;j<npeak;j++ )
        {
          int dx =Peak[i].ix-Peak[j].ix;
          int dy =Peak[i].iy-Peak[j].iy;
          double dist=sqrt ( dx*dx+dy*dy ) *4.;
          if ( dist_max<dist ) dist_max=dist;
        }
    }
}
//
int  PadRaster::Peak_mask ( int ix,int iy )
{

// simple peak definition...

// channel allready used...
  if ( PADflag[ix][iy]==1 ) return 0;

  double eseed= padxy44[ix][iy];

  if ( eseed<soglia_pic ) return 0;

  if ( eseed<padxy44[ix-1][iy-1] ||
       eseed<padxy44[ix  ][iy-1] ||
       eseed<padxy44[ix+1][iy-1] ||
       eseed<padxy44[ix-1][iy  ] ||
       eseed<padxy44[ix+1][iy  ] ||
       eseed<padxy44[ix-1][iy+1] ||
       eseed<padxy44[ix  ][iy+1] ||
       eseed<padxy44[ix+1][iy+1] ) return 0;

  double esomm8=  padxy44[ix-1][iy-1]+ padxy44[ix  ][iy-1] + padxy44[ix+1][iy-1];
  esomm8+=        padxy44[ix-1][iy  ]            +           padxy44[ix+1][iy  ];
  esomm8+=        padxy44[ix-1][iy+1]+ padxy44[ix  ][iy+1] + padxy44[ix+1][iy+1];

  if ( eseed<esomm8*.25 ) return 0;

  //  Gout<<"\n ev "<<Eventi_Fatti<<"  xy pad "<<ix<<" "<<iy<< "  eseed "<<eseed<<" etot "<<eseed+esomm8;
  // Trovato !   set flags...
  PADflag[ix-1][iy-1]= PADflag[ix][iy-1]= PADflag[ix+1][iy-1]= 1;
  PADflag[ix-1][iy]=   PADflag[ix][iy]=   PADflag[ix+1][iy]=1;
  PADflag[ix-1][iy+1]= PADflag[ix][iy+1]= PADflag[ix+1][iy+1]=1;

  return 1;
}


//================ooooOOO0OOOooo=   DevLkr    ==ooooOOO0OOOooo= ====================
//================ooooOOO0OOOooo=   DevLkr    ==ooooOOO0OOOooo= ====================
//================ooooOOO0OOOooo=   DevLkr    ==ooooOOO0OOOooo= ====================
void  DevLkr::Pad_Filling( )
{
  //reset parameters also for root
  epad28=m28x=m28qx=m28y=m28qy=0.0;
  epad44=m44x=m44qx=m44y=m44qy=0.0;
  dist_max=0.;
  padhit=0;

  if ( mhit<1 ) return;

  PAD.Reset();
  PAD.padhit=mhit;   //initial setting
  Cella CelSpace;
  double Dcell=2.;


  for ( int ipr=0; ipr<mhit; ipr++ )
    {


      double mx = M_Hits[ipr].Xdev.x;
      double my = M_Hits[ipr].Xdev.y;
      double e_in =M_Hits[ipr].e_rivela;

// larghezza della distribuzione deposito energia.
      double  sig=2.;
      double rap= e_in/M_Hits[ipr].Pdev.e;

      if ( M_Hits[ipr].pnome==Pip||M_Hits[ipr].pnome==Pim )
        {
          sig=1.2;

          if ( rap<0.1 ) sig=0.7;
          if ( rap<0.5 ) sig=1.4;

        }
      else  if ( M_Hits[ipr].pnome==Gam||M_Hits[ipr].pnome==Elec||M_Hits[ipr].pnome==Elep ) sig=1.;
      else if ( M_Hits[ipr].pnome==Mup||M_Hits[ipr].pnome==Mum ) sig=0.5;

// definisco CelSpace come la finestra delle (11x11) celle interessate  dalla particella incidente
// la finestra è di 121 celle;
//  la cella centrale è la (5,5);


      int ix=int ( ( mx+128. ) /Dcell );
      int iy=int ( ( my+128. ) /Dcell );


      int ixshift= ix-5;
      int iyshift= iy-5;

//         cout<<"\n Ev: "<<Eventi_Fatti<< " =====>  ";
//          cout<< setw(8)<<M_Hits[ipr].pnome<<" x "<< setw(8) <<x << " y " << setw(8)<<y <<" e "<< setw(8)<<e_in
//          <<" cel_chs " << setw(8)<<ix << " " << setw(8)<<iy<< flush;


      for ( int i=0;i<11;i++ )
        {
          for ( int j=0;j<11;j++ )
            {
              ix=ixshift+i;
              iy=iyshift+j;
              CelSpace.Set_index ( mx,my,ix,iy );
// deposito dell' energia rilasciata nelle celle dalla pertcella incidente.
              CelSpace.Rilascia_E ( sig,e_in );
              CelSpace.Cella2Pad();
            }
        }
    }


  PAD.Momenti();

// memorizzo per root
  padhit=PAD.padhit;
  epad28=PAD.epad28;
  m28x  =PAD.m28x  ;
  m28qx =PAD.m28qx ;
  m28y  =PAD.m28y  ;
  m28qy =PAD.m28qy ;
  epad44=PAD.epad44;
  m44x  =PAD.m44x  ;
  m44qx =PAD.m44qx ;
  m44y  =PAD.m44y  ;
  m44qy =PAD.m44qy ;
  quad1= PAD.quad1;
  quad2= PAD.quad2;
  quad3= PAD.quad3;
  quad4= PAD.quad4;

// ridefinizione dgli hit visti come peak

  PAD.Get_Peak();
  padhit=PAD.npeak;
  dist_max=PAD.dist_max;


// DEBUG DEBUG.......................if(mhit==1000....) per non eseguire il debug
  if ( mhit==1000&&padhit==2 )
    {


      Gout<<"\n ev "<<Eventi_Fatti<< " hit "<<mhit<<" peak "<<PAD.npeak<<"  "<<PAD.padhit;

      for ( int k=0;k<mhit;k++ )
        {
          double e_in =M_Hits[k].e_rivela;
          double mx = M_Hits[k].Xdev.x;
          double my = M_Hits[k].Xdev.y;

          Gout<<"\n Lkrhit "<<k<<" "<<M_Hits[k].pnome<<"  x,y "<<mx<<" "<<my<< "  e "<<e_in<<" "<<M_Hits[k].Pdev.e;



        }

      for ( int k=0;k<padhit;k++ )
        {

          double mx = PAD.Peak[k].ix;
          double my = PAD.Peak[k].iy;
          Gout<<"\n Padhit "<<k<<"  x,y "<<mx<<" "<<my;
        }

      Gout<<"\n    ";
      for ( int jx=0;jx<32;jx++ )
        {
          Gout<<setw ( 4 ) <<jx*8-124;
        }
      for ( int jy=0;jy<32;jy++ )
        {
          Gout<<"\n "<<setw ( 4 ) <<jy*8-124;
          for ( int jx=0;jx<32;jx++ )
            {

              int epad=int ( PAD.padxy44[jx][jy]*10. );

              if ( epad>0 ) Gout<<" "<<setw ( 3 ) <<epad;
              else    Gout<<"  . ";


            }
        }
    }
// end DEBUG....................

}



//=======================================
int DevLkr::Traccia ( Particella *pr )
{
// controllo il tracing della particella
// se è un pione
  if ( pr->Get_Nome() ==Pip||pr->Get_Nome() ==Pim )
    {
// ed è interno al lkr ed ha un cammino potenziale maggiore di un cm..
      if ( nowpos==INTERNO&&pr->pathok>1. )
        {
// calcolo in cammino prima di interagire 
          double pathfai = -liint* log ( EPSR + Pran() );
          if ( pathfai<pr->pathok )
            {
// interagisce prima di fare il cammino potenziale
// restante dichiaro la particella come distrutta: Ifato=idead=-3;
              pr->Set_Fato ( Distr,-3 ); 
              pr->pathok=pathfai;  // ridefinisco il cammino potenziale nel dev prima dellinterazione distruttiva...
            }
        }
    }
// se non è un pione traccio seguendo la procedura standard...
// se è un pione traccio la parte di cammino restante prima della distruzione.... 
  int ret=Device::Traccia ( pr );
  return ret;
}

//*********************************************
//*********************************************
DevLkr::~DevLkr()
{}

void DevLkr::Prgeom()
{
  Gout<<"==================  R i v e l a t o r e =================";
  Gout<<"\n --> < "<<nome <<" >  [ "<<fun
  <<" ]  Calorim. ";
  Device::Prgeom();

}
//***********************************************
void DevLkr::DataSmear()
{
// operate on the dev_variables of type (2) M_Hits
//  double x,y,ep,se,sx;


  const char *partnom;
  i=0;
  while ( i<=nhit )
    {
      x=M_Hits[i].Xdev.x;
      y=M_Hits[i].Xdev.y;
      partnom=M_Hits[i].pnome;
      ep=M_Hits[i].e_rivela;     //energia rilasciate nel dev

//

      if ( partnom==Gam|| partnom==Elec|| partnom==Elep )
        {


          //================================================================
          // new model from eta runs fit
          //  double rr, xfr=0.,slope=0.0;
          int Tail_todo=1;
          if ( Tail_todo==1 )
            {

//        if ( Pran() <=0.085 )    // first small  correction..to be checked
              if ( Pran() <=0.090 )    // 19.07.09

                {
                  if ( ep<=30.0 )
                    slope=32.5+1.25*ep;
                  else
                    slope=75.0;
                  rr =Pran();
                  if ( rr>1.e-9 )
//          if ( rr>1.e-9 )      // 19.07.09
                    {
                      ep*= ( 1.0+log ( rr ) /slope );
                      SelRea->Tailsrt_done++;
                    }
                }


              // add long tails from E/P nel 3 per mille
              if ( Pran() <=0.0021 )    // hadronproduction correction..
//        if ( Pran() <=0.003 )    // 19.07.09
                {
                  rr =Pran();
                  if ( rr>=0.000245 )
//         if ( rr>=0.0001 )  // 19.07.09
                    xfr=0.92+log ( rr ) /10.14;
                  else
                    xfr=.1;

                  ep*=xfr;
                  SelRea->Tailong_done++;
                }
// tails done
            }


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


          //  calcolo degli errori per il calorimetro
          //    x=Xdev[i][0]; y=Xdev[i][1];
          //    se= sqrt(0.01 + (0.001024+2.5e-5*ep)*ep);       // 1998
          se= sqrt ( 0.0081 + ( 0.001024+1.764e-5*ep ) *ep );      //1999
          ep=Pgauss ( se,ep );

          sx=sqrt ( 3.6e-3 + 0.1764/ep );   // correzione 2001
//     sx=sqrt ( 2.5e-3 + 0.1806/ep );  // 19.07.09  ... vedere nella ricostruzione....
//

          x=Pgauss ( sx,x );
          y=Pgauss ( sx,y );
        }

      else
        {
          // le altre particelle...pioni o "fuse" ..si usa un errore di 1 cm
          // per i mu si usa  2 mm.

          double sigma=1.; // for handron or "Fuse"
          if ( partnom==Mum||partnom==Mup )
            sigma=0.2;
          x=Pgauss ( sigma,x );
          y=Pgauss ( sigma,y );

          se= sqrt ( 0.0081 + ( 0.001024+1.764e-5*ep ) *ep );      //1999
          ep=Pgauss ( se,ep );

        }
      M_Hits[i].e_rivela=ep;
      M_Hits[i].Xdev.x=x;
      M_Hits[i].Xdev.y=y;
      M_Hits[i].Xlab=Devc2Lab ( M_Hits[i].Xdev );


      i++;
    }

}
//************************************
//=====================Simula il rivelatore =======================
int DevLkr::SimulaDev()
{
  // this procedure performs a filter of the measures
  // to simulate the device response
  // attention for  no gamma particles...

// opera sulla variabili di tipo (2) predifinite   M_Hits


  if ( mhit<1 ) return 0;
  emDev=0.;
  dead=0;
  // attenzione le dimensioni sono legate al massimo numere dellle particelle viste dal dev....
  int flag[30]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
  // hit number of gammas in the associated det. (remember:-1, 0,1,2,3,4,5)


  const char *pnome;
  // ----------
  // set flags and for  pions try to make life a little bit worse into lkr
  i=0;
  while ( i<=nhit )
    {
      flag[i]=1; // set all flags to 1
      pnome=M_Hits[i].pnome;
      if ( pnome==Pip||pnome==Pim )
        {
          M_Hits[i].e_rivela+=E_pio_Relesed( &M_Hits[i] );  //  released pion energy
        }
      i++;
    }


  // cuts


  double Ecutg=0.35;
//    double Ecutg=.35;

  ht=nhit;
  while ( ht>-1 )
    {
      // qui tagli in energia sui gamma  o elettroni
      pnome=M_Hits[ht].pnome;
      if ( ! ( pnome==Gam||pnome==Elec||pnome==Elep ) )
        {
          ht--;
          continue;
        }

      if ( flag[ht]!=0&& ( M_Hits[ht].e_rivela<Ecutg ) )
        {
          flag[ht]=-2;
          dead++;
        } //  geometry
      else if ( ( xyd= M_Hits[ht].Xdev.XYNorma() ) < 10. ) //taglio sul tubo del calorimetro
        {
          flag[ht]=-1;
          dead++;
        } // geometrical  cut:

      else   //if still alive,  then  dead cels test....
        {
          x= M_Hits[ht].Xdev.x;
          y= M_Hits[ht].Xdev.y;

          int  killcel=Get_DeadCel ( x,y );
          if ( killcel==1 )
            {
              flag[ht]=-3;
              dead++;
            }
        }
      ht--;
    }


  //  cout<<"\n Lkr nhit dead "<< nhit<<" "<<dead<<" "<<killcel<<" "<<flag[0];


//    for (int i=0;i<=nhit;i++)
//    Gout<<"\n Lk: "<< evento_.Gen.Event<<setprecision(3)
//    <<" xy "<< setw(8)<<M_Hits[i].Xdev.x<<" "<< setw(8)<<M_Hits[i].Xdev.y
//    <<" e "<< setw(8)<< M_Hits[i].e_rivela<<" pr "<<M_Hits[i].idm<< " flg "<<flag[i];

// Fusion procedure. Default to 1 ... zero for background test....

  if ( Fusione==1 )
    {
      evento_.Gen.Ifus=0;
      double ee1=0.0,ee2=0.0,eet=0.0;
      ht=nhit;
      while ( ht>-1 )
        {
          if ( flag[ht]<1 )
            {
              ht--;
              continue;
            }
          ee1=M_Hits[ht].e_rivela;

          gvet Daux=M_Hits[ht].Xdev;

//      Daux.print("Daux");

          //

          hi=ht-1;
          while ( hi>-1 )
            {
              if ( flag[hi]<1 )
                {
                  hi--;
                  continue;
                } //skip...

              sepm = ( M_Hits[hi].Xdev-Daux ).XYNorma();   // distanza proiettata su xy

//        if ( sepm<10. ) return 0;  // attenzione per Kl->pio0  ????


              if ( sepm<2. )           // allora  fusion sepm<2.
                {
                  ee2=M_Hits[hi].e_rivela;
                  eet=ee1+ee2;
                  if ( eet<0.00001 )
                    {
                      M_Hits[hi].Xdev= ( M_Hits[ht].Xdev+ M_Hits[hi].Xdev ) *0.5;
                      M_Hits[hi].Xlab= ( M_Hits[ht].Xlab+ M_Hits[hi].Xlab ) *0.5;
                      M_Hits[hi].tempo= ( M_Hits[ht].tempo + M_Hits[hi].tempo ) *0.5;
                    }
                  else
                    {
                      M_Hits[hi].Xdev= ( M_Hits[ht].Xdev*ee1+ M_Hits[hi].Xdev*ee2 ) /eet;
                      M_Hits[hi].Xlab= ( M_Hits[ht].Xlab*ee1+ M_Hits[hi].Xlab*ee2 ) /eet;
                      M_Hits[hi].tempo= ( M_Hits[ht].tempo*ee1 + M_Hits[hi].tempo*ee2 ) /eet;


                    }
                  M_Hits[hi].Pdev+=M_Hits[ht].Pdev;
                  M_Hits[hi].Vers=M_Hits[hi].Pdev.Verso();
                  M_Hits[hi].e_rivela=eet;
// to remeber that is an artefact define  idm =-1
                  M_Hits[hi].idm=-1;

                  if ( evento_.Gen.Event<1000 )
                    Gout<<" \n Fuse: "<<setw ( 8 ) <<evento_.Gen.Event<<setprecision ( 3 )
                    <<"  "<<setw ( 8 ) << M_Hits[hi].pnome<<" di "<<setw ( 8 ) <<ee2
                    <<" con " <<M_Hits[ht].pnome
                    <<" di "<<setw ( 8 ) <<ee1<<" sp "<<setw ( 8 ) <<sepm<<" " << M_Hits[hi].idm;

                  M_Hits[hi].pnome=Fuse;

                  flag[ht]=-4;
                  dead++;
                  evento_.Gen.Ifus++;
                  break;
                }
              hi--;
            }
          ht--;
        }
    }
  // end fusione

  // soglia rivelazione nel kripton   400. MeV

  for ( int i=0;i<=nhit;i++ )
    {


      if ( flag[i]>=0 )
        {
          if ( M_Hits[i].e_rivela< 0.35 )
            {
              flag[i]=-5;
              dead++;
            }
          else
            {
              double xcl=M_Hits[i].Xdev.x;
              double ycl=M_Hits[i].Xdev.y;
              if ( ( xcl*xcl+ycl*ycl ) <144. )
                {
                  flag[i]=-6;
                  dead++;
                }
            }
        }
    }


  // compact the mesures....
  int htt=-1;
  for ( int i=0;i<=nhit;i++ )
    {
      if ( flag[i]>=0 )
        {
          htt++;
          emDev+=M_Hits[htt].e_rivela=   M_Hits[i].e_rivela;
          if ( htt==i )  continue;   //avoid to rewrite the same data into itsself...
          M_Hits[htt]=M_Hits[i];
        }
    }

  mhit=htt+1;  //attenzione ridefinizione da 1 a n anziche' da 0 a n-1;/*
  nhit=mhit-1;

//cout<<"\n Lkr mhit "<< nhit<<" "<<mhit;

  // pad filling per trigger zero level

  Pad_Filling( );
//=================================
  if ( Debugon==1 )
    Gout<<"\n CalSimul  "<<evento_.Gen.Event
    <<" mhit "<<setw ( 3 ) <<mhit<<" of "<<nhit+1
    <<"    dead "<<setw ( 3 ) << dead
    <<" fusi "<<setw ( 3 ) << evento_.Gen.Ifus
    <<" flag "<< flag[0]<<" "<<flag[1]<<" "<<flag[2]
    <<" " << flag[3]<<" "<<flag[4]<<" "<<flag[5]<<endl;

  return htt;
}
//******************************************
double DevLkr::E_pio_Relesed( Bufdev *MH)
{
  // nota MH sono i parametri misurati di un hit (vedi BufDev *)
  static double anor=0.0;

  // define the energy left in the hadronic interction if it has flagged idead=-3 
  if ( MH->idead!=-3 ) return 0.0;
  double brnd=Pran();
  if ( brnd<.3 ) return 0.;


  double rp=0.0, fp=0.0;
  //double alfa=10., beta=10.; originali
  double alfa=12., beta=10.;
  if ( anor==0.0 )  // prima volta calcolo l'integrale della fun. di distribuzione.
    {
      while ( rp<1. )
        {
          anor+= ( 1.-exp ( -alfa*rp ) ) *exp ( -beta*rp ) *rp*rp;
          rp+=0.01;
        }
    }
  // normal path(gvet_0)

  double rnd=Pran();
  rp=0.0;
  while ( rp<1. )
    {
      fp+= ( 1.-exp ( -alfa*rp ) ) *exp ( -beta*rp ) *rp*rp/anor;
      if ( fp>rnd ) break;
      rp+=0.01;
    }
  // correzione dell'energia rilasciata in funzione della profondita' di interazione
  double depth= ( 2.*Lout.z-MH->last_path ) /60.;
  double reduci =1.- exp ( -depth*depth );
  return rp* MH->Pdev.e*reduci;

}
//===============================================================

void DevLkr::Reset_RootData()
{
  hit_recorded=padhit=0;
  emDev=0.0;
  epad28=m28x=m28qx=m28y=m28qy=0.0;
  epad44=m44x=m44qx=m44y=m44qy=0.0;
  dist_max=0.;

}


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


const float DevLkr::Deadcl[90][2] =
{
  {  -124.3518     ,   33.55526     },{  -104.6134     ,   21.71223     },
  {  -102.6396     ,  4.0000000E-07 },{  -96.71809     ,  -23.68606     },
  {  -88.82274     ,  -78.95355     },{  -82.90123     ,   61.18901     },
  {  -69.08435     ,   75.00587     },{  -67.11052     ,   45.39830     },
  {  -67.11052     ,   43.42446     },{  -67.11052     ,  -33.55526     },
  {  -67.11052     ,  -65.13667     },{  -65.13667     ,  -3.947674     },
  {  -59.21516     ,   114.4827     },{  -57.24132     ,   67.11052     },
  {  -53.29364     ,  -25.65990     },{  -51.31981     ,  -122.3780     },
  {  -49.34597     ,   112.5088     },{  -47.37213     ,  -1.973835     },
  {  -45.39829     ,   124.3519     },{  -45.39829     ,   122.3780     },
  {  -39.47677     ,  -3.947674     },{  -33.55526     ,   51.31981     },
  {  -25.65990     ,  -61.18900     },{  -23.68606     ,   23.68607     },
  {  -21.71222     ,  -100.6658     },{  -19.73838     ,   61.18901     },
  {  -17.76455     ,   106.5873     },{  -15.79071     ,   92.77043     },
  {  -5.921513     ,   65.13669     },{  -1.973835     ,   63.16285     },
  {   1.973843     ,   80.92740     },{   1.973843     ,  -33.55526     },
  {   3.947681     ,   19.73839     },{   9.869198     ,  -100.6658     },
  {   19.73839     ,  -96.71809     },{   21.71223     ,  -35.52909     },
  {   23.68607     ,   82.90123     },{   23.68607     ,  -33.55526     },
  {   25.65991     ,   122.3780     },{   25.65991     ,   5.921520     },
  {   27.63375     ,  -96.71809     },{   31.58142     ,   27.63375     },
  {   31.58142     ,   27.63375     },{   33.55526     ,   124.3519     },
  {   33.55526     ,   122.3780     },{   33.55526     ,  -7.895351     },
  {   33.55526     ,  -7.895351     },{   33.55526     ,  -23.68606     },
  {   35.52910     ,   102.6396     },{   35.52910     ,   96.71811     },
  {   35.52910     ,   59.21517     },{   35.52910     ,   53.29365     },
  {   35.52910     ,   47.37214     },{   35.52910     ,  -19.73838     },
  {   37.50294     ,   118.4303     },{   37.50294     ,   104.6135     },
  {   37.50294     ,   59.21517     },{   37.50294     ,  -65.13667     },
  {   39.47678     ,   53.29365     },{   41.45062     ,   39.47678     },
  {   41.45062     ,  -7.895351     },{   43.42446     ,   102.6396     },
  {   43.42446     ,   69.08437     },
  {   43.42446     ,   63.16285     },
  {   45.39830     ,   122.3780     },{   45.39830     ,   118.4303     },
  {   45.39830     ,   71.05820     },{   45.39830     ,   69.08437     },
  {   49.34597     ,   124.3519     },{   49.34597     ,   122.3780     },
  {   51.31981     ,   37.50294     },{   51.31981     ,   27.63375     },
  {   51.31981     ,  -55.26748     },{   55.26749     ,   73.03204     },
  {   55.26749     ,  -90.79659     },{   61.18901     ,   86.84891     },
  {   63.16285     ,  -19.73838     },{   65.13669     ,  -61.18900     },
  {   73.03204     ,  -90.79659     },{   76.97972     ,  -43.42445     },
  {   80.92740     ,  -35.52909     },{   90.79659     ,  -51.31981     },
  {   90.79659     ,  -51.31981     },
  {   92.77043     ,   80.92740     },
  {   92.77043     ,  -61.18900     },{   98.69194     ,   75.00587     },
  {   110.5350     ,  -1.973835     },{   116.4565     ,  -51.31981     },
  {   118.4303     ,  -17.76455     },{   118.4303     ,  -47.37213     }

};

int DevLkr::Get_DeadCel ( double x,double y )
{
  // Lkr dead cell finding
  int  killcel=0;
  for ( int l=0;l<90;l++ )
    {
      if ( fabs ( x-Deadcl[l][0] ) <2. && fabs ( y-Deadcl[l][1] ) <2. )
        {
          killcel=1;
          break;
        }
    }
  // killcel =0 (x,y)not in a dead cell
  // killcel=1 (x,y)  is in a dead cell
  return killcel;
};