flyopunta/src/src_ana/anakl2.cpp

/***************************************************************************
                          anake2.cpp  -  description
                             -------------------
    begin                : Mon Nov 24 2003
    copyright            : (C) 2003 by gmp
    email                : peppe@newpeppe
  ***************************************************************************
  *                                                                         *
  *   NA48  simulation program.                                             *
  *                                                                         *
  ***************************************************************************/
//la successione degli include può essere causa di errore.....

#include "flyoh.h"
#include "anakl2.h"

using namespace std;

AnaKl2::AnaKl2()
{
  tipo=23;
  nfit=1;
  titol="K->l nu" ;
  nome="Kl2";
  Gout<<"\n Ridefinito e creato il fit di "<<titol;

  if(Le_Tracce==0)
    {
      Gout<<"\n Error.. "<<titol<<"  analisys ... Tracks non defined!";
    }
  Part_db *pm;
  pm=Matter->Pntmatter("Kp");
  mask= pm->Get_Massa();
  pm=Matter->Pntmatter("mup");
  masmu= pm->Get_Massa();
  pm=Matter->Pntmatter("pip");
  maspi= pm->Get_Massa();
  pm=Matter->Pntmatter("elec");
  mase= pm->Get_Massa();

}
AnaKl2::~AnaKl2()
{}

int AnaKl2::Fisica()
{
  count_call++;
  //   note:    otagon cal cut not here....

  //  if(lkry->Brem_done<1) return -4;   //solo gli eventi con brems
  kp=kpx=kpy=0.0;
  enlk=pmis=pmisx=pmisy=pmisz=0.0;
  ptrs=emisq=mumisq=pimisq=0.0;


// nota Le_Tracce->First_trk corrispomde a Tracc[0]
  // atenzione forzata ad uno per vedere la statistica......

  nt=Le_Tracce->Tot_trk+1;
  if(nt<1)
    {
      getta_ch++;
      return(-1);
    } // non c'e' la traccia carica

  //   lkry->Pdev[0].print( evento_.Event,lkry->pnome[0],0,"ana");Gout<<"  "<<lkry->mhit;

  if(Errori>0&&Scatter>0)
    {
      if(lkry->mhit==0)
        {
          getta_en++;
          return -3;
        }// calorimetro vuoto
      enlk=0.0;
      if((enlk = lkry->M_Hits[0].e_rivela)<0.1)
        {
          getta_en++;
          return -3;
        }
      ;  // cut on energia misurata in lk...
    }
  else  // to run also if not error simulation.. for debug
    {
      if(lkry->nhit==-1)
        {
          getta_en++;
          return -3;
        }// calorimetro vuoto
      enlk=0.0;
      if((enlk = lkry->M_Hits[0].Pdev.e)<0.1)
        {
          getta_en++;
          return -3;
        }
      ;  // cut on energia misurata in lk...
    }

// calcolo la distanza del pione estrapolato al calorimetro dallo hit visto nel cal;
// si usa solo la prima.... si suppone che gli eventi a più tracce siano da scartare nella analisi..
   discal=(Tracc[0]->Xlk-lkry->M_Hits[0].Xdev).XYNorma();
    Xtcal=Tracc[0]->Xlk[0];
    Ytcal=Tracc[0]->Xlk[1];

   

  // Impulso dell'elettrone * carica  dallo spettrometro ,
  // si potrebbe usare quello dal calor.

  double puls=Tracc[0]->ptrk;


  if(puls>100.)
    puls=100.;

  pq=puls*Tracc[0]->charge;
  rapsel=-1.;
  if(abs(pq)>0.0)
    rapsel=enlk/abs(pq);

  //  Calcolo il quadrimomento del leptone suposto un elettrone
  Mome=gvet(Tracc[0]->Get_P());      // impulso calcolato in DevTrk
  // impulso fisso del K
  Kmom.setvn(0.,0.,75.,.49367);

  // missing mass come elettrone
  Emom.setvn(Mome, mase);    // quadriimpulso come se elettrone
  Qtr=Kmom-Emom;       // quadrimomento mancate o trasferito( neutrino)
  Qtr.Invar();         // Completa il calcolo della massa invariante....
  pmis   = Qtr.norma;  // impulso nel nu  non dipemde dalla massa del leptone....
  pmisx  = Qtr[0];
  pmisy  = Qtr[1];
  pmisz  = Qtr[2];
  ptrs=sqrt(Emom[0]*Emom[0]+Emom[1]*Emom[1]); //impulso  transverso
  emisq  = Qtr.mq;
  // debug
  /*
     Tracc[0]->Get_P().print(evento_.Event,"Emb",0,"Trk");;   
     Mome.print(evento_.Event,"Mome",0,"Anal2");
     Emom.print(evento_.Event,"Emom",0,"Anal2");
     Qtr.print(evento_.Event,"Qtr0",0,"Anal2");
  */


  //  come se muone
  Emom.setvn(Mome, masmu);  // quadriimpulso come se muone
  Qtr=Kmom-Emom;           // quadrimomento mancate o trasferito( neutrino)
  Qtr.Invar();             // Completa il calcolo della massa invariante....
  mumisq  = Qtr.mq;

  //   come se pione
  Emom.setvn(Mome, maspi);    // quadriimpulso come se pione
  Qtr=Kmom-Emom;       // quadrimomento mancate o trasferito( neutrino)
  Qtr.Invar();         // Completa il calcolo della massa invariante....
  pimisq  = Qtr.mq;

  //*************************************************************************
  // using the kabes....***********************************
  if(kabs3==0)
    {
      getta_kabes++;
      return -2;
    }  //non c'e' la kabes 3 si potrebbe uscire..

  //  qui si usano le kabes..solo per la direzione ...da rivedere

  yb1=kabs1->M_Hits[0].Xdev.y;
  yb3=kabs3->M_Hits[0].Xdev.y;

  xb1=kabs1->M_Hits[0].Xdev.x;
  xb3=kabs3->M_Hits[0].Xdev.x;
  //   VK.setvn( xb3-xb1, yb3-yb1, 800.);//stima  fare meglio
  VK.setvn( xb3, yb3, 6378.);

  kp=75.0;
  PK=VK.Verso()*kp;   // vettore impulso del K
  kpx    = PK.x;
  kpy    = PK.y;

  Kmom.setvn(PK,mask);  // quadrimomento del K
  //***

  //  come  se elettrone
  Emom.setvn(Mome, mase);  // quadriimpulso come se muone
  Qtr=Kmom-Emom;           // quadrimomento mancate o trasferito( neutrino)
  Qtr.Invar();             // Completa il calcolo della massa invariante....
  miskbq  = Qtr.mq;
  //    ---------------------------------
  /*
     Kmom.print(evento_.Event,"Pk",0,"ric");
     SelRea->avo->P.print(evento_.Event,"Pk",0,"Gen");;   
   
     Mome.print(evento_.Event,"Mom2",0,"Anal2");
  //   Emom.print(evento_.Event,"Emo2",0,"Anal2");
     Qtr.print(evento_.Event,"Qtr2",0,"Anal2");
    
  */
  count_wnt++;
  return 1 ;
}