flyopunta/src/src_dev62/devlav.cpp

/***************************************************************************
 *   Copyright (C) 2009 by giuseppe Pierazzini                                      *
 *   giuseppe@pierazzini.it                                                               *
 *                                                                         *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.                  *
 *                                                                         *
 *  Dipartimento di Fisica  E.Fermi / INFN Pisa Italy                      *
 ***************************************************************************/
#include "parm.h"
#include "devlav.h"
#include "part_db.h"

using namespace std;

DevLav::DevLav()
        : DevCrn()
{
    devtype=TypDevLav;
    devclass="Large Angle Veto";         
Gout<<"\n\n < "<<nome <<" > Dev Id "<<idev<<" Device type <"<<devclass<<"> typ "<<devtype; 
}

DevLav::~DevLav()
{
}

int DevLav::SimulaDev()
{
    if(mhit<1) return 0;
    double Thres_gam=.100; 
    double Thres_mu=.300; //threshold in p!
    //double Thres_mu=0.1; //prova per vedere se cambia molto l'efficienza
    double RefIndx=1.3;  // leadglass EFFECTIVE refraction index (cherenkov threshold != real threshold!)
    //double RefIndx=2.; //prova per vedere se cambia molto l'efficienza
    double beta=0.;

// loop sulle particelle viste;
    qvet Eimp;
    int loop=0;
    const char *Part_nome;
    while (loop<mhit)
    {
        Eimp=M_Hits[loop].Pdev;
        Part_nome= M_Hits[loop].pnome;

        if ( Part_nome==Gam|| Part_nome==Elec|| Part_nome==Elep)
        {
            if (M_Hits[loop].Pdev.e <Thres_gam)
                M_Hits[loop].idead=-99999;
            /*uccido la particella*/
        }

        else if ( Part_nome==Mup|| Part_nome==Mum)
        {
            if (M_Hits[loop].Pdev.mom <Thres_mu)
                M_Hits[loop].idead=-99999;
            /*uccido la particella*/
        }

        else if ( Part_nome==Neutrone|| Part_nome==A_Neutrone)
        {
          if ( Eff_Neutron(Part_nome,Eimp)==0) M_Hits[loop].idead=-99999;
            /*uccido la particella*/
        }
        
        
        else if(M_Hits[loop].Pr->Get_Charg()!=0)
        {
          beta =Eimp.mom/Eimp.e;
//        Eimp.print("Pvet"); Gout<<" "<<M_Hits[loop].pnome<<" beta "<< beta ;
          if(beta <1./RefIndx)  M_Hits[loop].idead=-99999;
          /* attenzione  atomi inclusi..*/
        
        }
        
          else if ( M_Hits[loop].Pr->Get_Charg()==0)
        {
          if ( Eff_Neutron(Part_nome,Eimp)==0) M_Hits[loop].idead=-99999;
            /* per il momento trattatano come i neutroni....*/
        }
        
          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;
}
int DevLav::Eff_Neutron( const char *Nome,qvet Eimp)
{
    // attenzione qui si deve arrivare solo con i neutroni... altrimenti con gli atomi....Axy ...
    double eff=.70;
    if (Eimp.mom<0.1) eff=0.05;
    else if (Eimp.mom<1.) eff=.50;
    if (Pran()<eff)return 1;
    return 0;
}
//===============
int DevLav::Get_Structinfo()
{
// ritorna in gradi la posizione nel veto  system
    double x = M_Hits[nhit].Xdev.x;
    double y = M_Hits[nhit].Xdev.y;
    int gradi= atan2(y,x)*57.30;
    return gradi;
}