flyopunta/src/src_dev62/devsac.cpp

/***************************************************************************
                          devsac.cpp  -  description  2010
                             -------------------

  ***************************************************************************
  *                                                                         *
  *   NA62  simulation program.                                             *
  *                                                                         *
  ***************************************************************************/
#include "parm.h"
#include "evento.h"
#include "devsac.h"
#include "particella.h"
#include "part_db.h"

DevSac::DevSac()
{
    devtype=TypDevSac;
    devclass="Small Angle Calorimeter";
Gout<<"\n\n < "<<nome <<" > Dev Id "<<idev<<" Device type <"<<devclass<<"> typ "<<devtype; 
         
}
DevSac::~DevSac()
{}
//----------------- --------
void DevSac::Prgeom()
{

    Gout<<"==================  R i v e l a t o r e =================";
    Gout<<"\n --> < "<<nome <<" >  [ "<<fun
    <<" ]  Small Angle  Calorimeter:  Rettangolare"<<std::endl;
    Device::Prgeom();
}
//=======-----------------------------------
int DevSac::SimulaDev()
{
          if(mhit<1)return 0;
    double Thres_gam=.100;
    double Thres_neu=.100; //threshold su p!

// loop sulle particelle viste;
    qvet Eimp;
    int loop=0;
    const char *Part_nome;
    //   double Probint= 1.-exp(-2.*Lout.z/liint);
    //double Probint= .99752;
double Probint= .715418; //10.5 cm Pb + 10.5 cm Scintillator (70 strips da 1..5 +1.5 mm l'una)
    while (loop<mhit)
    {
        Eimp=M_Hits[loop].Pdev;
        Part_nome= M_Hits[loop].pnome;

        if ( Part_nome==Gam)
        {

// nota i gamma sono dampati.... quindi  l'energia è automaticamente memorizzata in me
            if ( Eimp.e<Thres_gam)  M_Hits[loop].idead=-99999;
            /*uccido la particella sotto soglia*/
        }
        else if ( Part_nome==Neutrone|| Part_nome==A_Neutrone)
        {
            if (Pran()>Probint|| Eimp.mom<Thres_neu)  M_Hits[loop].idead=-99999;
            /*uccido la particella se non interagisce o sotto soglia*/
        }

        else if (M_Hits[loop].Pr->Get_Charg()==0) //tutte le altre particelle neutre
        {
            //  per ora come per il neutrone
            if (Pran()>Probint|| Eimp.mom<Thres_neu)  M_Hits[loop].idead=-99999;
            /*uccido la particella*/
        }

        loop++;

    }
//==============  La parte che segue è fissa!
// elimino le paricelle morte, compattizzo i dati e  ridefinisco nhit e mhit;
//  Gout<<"\n Prima  del coppattamento..."<<nome<<" "<< nhit<<" " <<mhit<< "  "<<M_Hits[nhit].pnome;

    loop =0;
    int newhit=mhit;
    while (loop<newhit)
    {
        if (M_Hits[loop].idead!=-99999)loop++;
        else if (newhit==1)newhit=0;
        else {
            newhit--;
            for (int i=loop;i<newhit;i++)
            {
                M_Hits[i]=M_Hits[i+1];
            }
        }
    }
    mhit=newhit;
    nhit=mhit-1;
//     Gout<<"\n Dopo il coppattamento..."<<nome<<" "<<nhit<<" " <<mhit;
//     if (nhit>-1)  Gout<<"  "<<M_Hits[nhit].pnome;
    return 0;





    return 1.;
}