flyopunta/src/src_ana/anakmunu.cpp

/***************************************************************************
                          anakmunu.cpp  -  description
                             -------------------
    begin                : Sun Jan 18 2004
    copyright            : (C) 2004 by gmp
    email                : peppe@newpeppe
  ***************************************************************************
  *                                                                         *
  *   NA48  simulation program.                                             *
  *                                                                         *
  ***************************************************************************/

#include <flyoh.h>
#include <anakmunu.h>

using namespace std;


AnaKmunu::AnaKmunu()
{
    tipo=29;
    nfit=1;
    titol="K--> mu nu " ;
    nome="Kmunu";
    mas_Kp= ptr_Kapp->Get_Massa();
    mas_mu=ptr_Mup->Get_Massa();
                
          Pk_ric =new GigaBeam();  
  
                
    Gout<<"\n Attivato il fit di "<<titol;

}
AnaKmunu::~AnaKmunu()
{
}
//GigaBeam *Pk_ric=0;

int AnaKmunu::Fisica()
{
        
                
    count_call++;

    np=0;
    pxkp=pykp=pzkp=momk=0.;
    pxmu=pymu=pzmu=mommu=0.;
    emu =masnu= mqnu=0.;
    pxnu=pynu=pznu=momnu=0.;
    xnu=ynu=znu =0.;
 
    //taglio sul vertice generato (taglia solo gli eventi che finiscono sul collimatore)da togliere!!!
    if (evento_.Gen.zv < 10000.) {
        counter[0]++;
        return -1;
    }

    count_decay++;
    np=0;
    int npt = 0 ;
    if (Trkstraw!=0)  npt = Trkstraw->ntrk_done;
    if (npt<1) {
        counter[1]++;
        return -1;
    }
    
    // ricavo l'impulso del K dal gigatrak
    Pk_mom=Pk_ric->Beam_Rec();   // impulso del beam dal gigatrack
    Pkp.setvn ( Pk_mom,mas_Kp ); // definizione del quadrimomento
    Pkp.putvn( pxkp,pykp,pzkp,momk );  // salvo per root


// Dato che si selezionano  eventi con solo una traccia, uso  l'indirizzo della traccia [0]

    Track  Trk=Trkstraw->Traccia[0];
    charge=Trk.charge;
    chiq_trk=Trk.chiq_fit;

    if (Trk.chiq_fit>2.) {
        counter[6]++;
        return -1;
    }


    if ( charge !=1. ) { //even for wrong sign study, positive cut later!
        counter[2]++;
        return -1;
    }

// impulso del mu+ dallo spettromentro
    Pmu.setvn(Trk.P_in,mas_mu);
    Pmu.putvn( pxmu,pymu,pzmu,mommu );
// taglio sull'impulso del mup
    if (mommu<15. || mommu>75.) {
        counter[3]++;
        return -1;
    }

    // 4 camere coinvolte
    if (Trk.camere!=4) {
        counter[4]++;
        return -1;
    }
// taglio radiale sulle strawchambers; d < 7.5 || d > 110
    gvet XCamera,XCamera1;
    double dq=0.;
    XCamera= Trk.Xmg + Trk.P_in*((St1x->Get_Cface().z-Trk.Xmg.z)/Trk.P_in.z);
    XCamera1=XCamera;
    XCamera=St1x->Lab2cDev(XCamera) ;
    dq=XCamera.XYNorma();
    if (dq<7.5||dq>110.) {
        counter[5]++;
        return -1;
    }
    XCamera= Trk.Xmg + Trk.P_in*((St2x->Get_Cface().z-Trk.Xmg.z)/Trk.P_in.z);
    XCamera=St2x->Lab2cDev(XCamera) ;
    dq=XCamera.XYNorma();
    if (dq<7.5||dq>110.) {
        counter[5]++;
        return -1;
    }
    XCamera= Trk.Xmg + Trk.P_out*((St3x->Get_Cface().z-Trk.Xmg.z)/Trk.P_out.z);
    XCamera=St3x->Lab2cDev(XCamera) ;
    dq=XCamera.XYNorma();
    if (dq<7.5||dq>110.) {
        counter[5]++;
        return -1;
    }
    XCamera= Trk.Xmg + Trk.P_out*((St4x->Get_Cface().z-Trk.Xmg.z)/Trk.P_out.z);
    XCamera=St4x->Lab2cDev(XCamera) ;
    dq=XCamera.XYNorma();
    if (dq<7.5||dq>110.) {
        counter[5]++;
        return -1;
    }



    emu=Pmu*Pkp/mas_Kp;            // energia del pi nel sistema Kappa
    Pnu=Pkp-Pmu;
    mqnu=Pnu.Invarq();
    masnu=Pnu.m;

    Pnu.putvn( pxnu,pynu,pznu,momnu );
    enu=Pnu*Pkp/mas_Kp;

// calcolo dell'impulso trasverso del mu
    gvet ptrasv;
    ptrasv= Pmu&Pkp;
    ptrasv.Norma();
    tmom=ptrasv.mom/Pkp.mom;



     
    
    Vtxstruct Xdecay(2);
    Xdecay.P[0].setvn(0.,0.,10240.); // beam at G3 (0.,0.,10240.); Fare meglio
    Xdecay.V[0]=Pkp.NVerso();       //  versore del beam...
    Xdecay.P[1]=Trkstraw->Traccia[0].Xmg; // coordinate della traccia 1 al centro del magnete (definizione)
    Xdecay.V[1]=Trkstraw->Traccia[0].P_in.NVerso(); // versore della traccia 1
    Xdecay.Vertice();
    chiq_vert=Xdecay.chiq;
    Xdecay.Vtx.putv( xnu,ynu,znu );
 
        
                
//     //taglio sul vertice ricostruito (livello 0)
    if (znu<10250.||znu>17000.) {
        counter[7]++;
        return -1;
    }

//taglio su xnu e ynu
//per il segnale taglia circa 2 eventi su 10^4 che hanno passato la selezione, per il gas il 10%!
//collimatore a 10000 con raggio 3.5 (fascio generato in 0)
//=> r_min = 3.5, z_min = 10000 => r_max = 5.95 se z_max = 17000

    gvet beam0 (0.,0.,10240.) ;   // beam point at gigatrack 3
    gvet beamc (-10.7,0.,18366.) ;// beam point at st1x chamber
    gvet bvers;
    bvers = beamc-beam0 ;
    bvers = bvers.NVerso() ;
    double dbeam = (Xdecay.Vtx - beam0)%bvers;
    gvet tbeam;
    tbeam =(Xdecay.Vtx - beam0)-bvers*dbeam;
    double dtrv = tbeam.Norma();
    if (dtrv > dbeam*0.00035) {
        counter[8]++;
        return -1;
    }

    // taglio sui cluster del calorimetro
    if ( lkry->hit_recorded>1) {
        counter[9]++;
        return -1;
    }

    // calcolo le coordinate estrapolate della traccia al piano del calorimeter

    gvet X_at_lkr;
    X_at_lkr= Trk.Xmg + Trk.P_out*((lkry->Get_Cface().z-Trk.Xmg.z)/Trk.P_out.z);
    X_at_lkr.putv(xlac,ylac,zlac);

    double rqcl= X_at_lkr.XYNorma();
    rqcl = rqcl*rqcl;
    // accettanza ottagonale
    if (sqrt(rqcl)<=15. || fabs( X_at_lkr.x) >= 113. ||fabs( X_at_lkr.y) >= 113. ||
            (fabs(X_at_lkr.x)+fabs(X_at_lkr.y)) >= 113.*sqrt(2.)) {
        counter[10]++;
        return -1;
    }

    //distanza traccia cluster < 10 cm
    if (lkry->hit_recorded>0) {
        double dx = lkry->mxl[0] - X_at_lkr.x;
        double dy = lkry->myl[0] - X_at_lkr.y;
        if (dx*dx + dy*dy > 100) {
            counter[11]++;
            return -1;
        }
    }

// accettanza dei muv1,2,3

    gvet X_at_muv;
    X_at_muv= Trk.Xmg + Trk.P_out*((muvet1->Get_Cface().z-Trk.Xmg.z)/Trk.P_out.z);
    if (fabs(X_at_muv.x)<10.&&fabs(X_at_muv.y)<10.) { //taglio interno
        counter[12]++;
        return -1;
    }
    if (fabs(X_at_muv.x)>120.||fabs(X_at_muv.y)>120) { //taglio esterno
        counter[13]++;
        return -1;
    }

    // muvet2
    X_at_muv= Trk.Xmg + Trk.P_out*((muvet2->Get_Cface().z-Trk.Xmg.z)/Trk.P_out.z);
    if (fabs(X_at_muv.x)<10.&&fabs(X_at_muv.y)<10.) { //taglio interno
        counter[14]++;
        return -1;
    }
    if (fabs(X_at_muv.x)>120.||fabs(X_at_muv.y)>120) { //taglio esterno
        counter[15]++;
        return -1;
    }

    // muvet3
    X_at_muv= Trk.Xmg + Trk.P_out*((muvet3->Get_Cface().z-Trk.Xmg.z)/Trk.P_out.z);
    if (fabs(X_at_muv.x)<10.&&fabs(X_at_muv.y)<10.) { //taglio interno
        counter[16]++;
        return -1;
    }
    if (fabs(X_at_muv.x)>120.||fabs(X_at_muv.y)>120) { //taglio esterno
        counter[17]++;
        return -1;
    }


//*************************************************
    count_wnt++;
    np=1;
                
                Set_tempo_parms();
                timedif=ana_dtim;
    return 1;

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


//----------------------------------------------------
void AnaKmunu::print_scale()
{
    Gout<<"\n -----> "<<titol<<" Analysis Summary <-------\n";
    Gout<<"\n Fit entries          "<<count_call;
    Gout<<"\n Good decays          "<<count_decay << " (after zv cut)";

    for (int i=0;i<30;i++)
    {
        if (counter[i]>0)  Gout<<"\n Rejection "<<setw(3)<<i<<"         "<<counter[i];

    }

    Gout<<"\n Good events          "<<count_wnt;


    Gout<<"\n\n ======>           d o n e     <============ "<<std::endl;
}