flyopunta/src/src_ana/anakpivv.cpp

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

#include <flyoh.h>
#include <anakpivv.h>

using namespace std;


AnaKpivv::AnaKpivv()
{
    tipo=25;
    nfit=1;
    titol="K--> pip nu nu " ;
    nome="Kpivv";
    mas_Kp= ptr_Kapp->Get_Massa();
    mas_pip=ptr_Pip->Get_Massa();
    Gout<<"\n Attivato il fit di "<<titol;
    Pk_ric =new GigaBeam();  
}
AnaKpivv::~AnaKpivv()
{
}


int AnaKpivv::Fisica()
{

    count_call++;

    np=0;
    pxkp=pykp=pzkp=momk=0.;
    pxpi=pypi=pzpi=mompi=0.;
    epi =masvv= mqvv=0.;
    pxvv=pyvv=pzvv=momvv=0.;
    xvv=yvv=zvv =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 pione+ dallo spettromentro
    Ppi.setvn(Trk.P_in,mas_pip);
    Ppi.putvn( pxpi,pypi,pzpi,mompi );
// taglio sull'impulso del pip
    if (mompi<15. || mompi>35.) {
        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;
    }



    epi=Ppi*Pkp/mas_Kp;            // energia del pi nel sistema Kappa
    Pvv=Pkp-Ppi;
    mqvv=Pvv.Invarq();
    masvv=Pvv.m;

    Pvv.putvn( pxvv,pyvv,pzvv,momvv );
    evv=Pvv*Pkp/mas_Kp;

// calcolo dell'impulso trasverso del pi
    gvet ptrasv;
    ptrasv= Ppi&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( xvv,yvv,zvv );
 
     
//     //taglio sul vertice ricostruito (livello 0)
    if (zvv<10250.||zvv>17000.) {
        counter[7]++;
        return -1;
    }

//taglio su xvv e yvv
//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;
    return 1;

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


//----------------------------------------------------
void AnaKpivv::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;
}