flyopunta/src/src_ana/anake2g.cpp

/***************************************************************************
                          anake2g.cpp  - 
   gmp/SGal  4 ottobre 2010 
  ***************************************************************************/
//la successione degli include può essere causa di errore.....

#include "flyoh.h"
#include "anake2g.h"

using namespace std;

gvet  Vert_coor(gvet,gvet,gvet,gvet);

AnaKe2g::AnaKe2g()
{
  tipo=24;
  nfit=1;
  titol="K->e v g" ;
  nome="Ke2g";
  Gout<<"\n Ridefinito e creato il fit di "<<titol<<endl;

  if(Le_Tracce==0)
    {
      Gout<<"\n Error.. "<<titol<<"  analisys ... Tracks non defined!"<<endl;
      exit(0);
    }
  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("elep");
  mase= pm->Get_Massa();

}
AnaKe2g::~AnaKe2g()
{}





//==================================

//utility vectors  for Fisica 
 qvet PK,Pelep,Pke,Pgam;
 gvet Vtxm, Vgam;
 
int AnaKe2g::Fisica()
{
  count_call++;
  
  kp=0.0;
  engm=enlk=miske=miskex=miskey=miskez=0.0;
  ptrske=emisq=0.0;
  vtx=vty=vtz=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
    } 

// test on lkry hits 
      if(lkry->mhit<2)
        {
          getta_en++;
          return -3;// calorimetro con pochi hits
        }
   
        
     

// calcolo la distanza dell'elettrone estrapolato al calorimetro dagli hits visti per trovare
// il gamma e il positrone

// test primo hits 
// controllare che Xdev sia quello giusto con errori!
 int  pr_el=-1 , pr_gm=-1;
 
   engm=enlk=Xtcal= Ytcal=0.0;
   if( ( discal=(Tracc[0]->Xlk-lkry->M_Hits[0].Xdev).XYNorma())<1.)
   {
     pr_el=0;
    pr_gm=1;     // l'altro hit è del gamma
   }
   else if( (discal=(Tracc[0]->Xlk-lkry->M_Hits[1].Xdev).XYNorma())<1.)
   {
   pr_el=1;
    pr_gm=0;
   }  
   else  
   {
      getta_ch++;
      return(-1);// non c'e' matching 
    }       
    Xtcal=Tracc[0]->Xlk[0];
    Ytcal=Tracc[0]->Xlk[1];    

   
//      qvet PK,Pelep,Pgam;
      // K incidente 
      PK.setvn(0.,0.,75.,mask);   // impulso fisso del K (px,py,pz,Massa)
      
      // gamma 
      if(( engm= lkry->M_Hits[pr_gm].e_rivela)<0.5)
      {
          getta_en++;
          return -3;// cut on energia misurata in lkry...
        }
       
       
       
      
       

      
      // traccia elettrone  con soglia
      if((enlk = lkry->M_Hits[pr_el].e_rivela)<0.1)
        {
          getta_en++;
          return -3;// cut on energia misurata in lkry...
        }
      gvet Mome=gvet(Tracc[0]->Get_P()); // impulso della traccia carica calcolato in DevTrk
      Pelep.setvn(Mome, mase);       
      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 vertice in Vert_coor()
        gvet Xkb(0.,0.,8000.),Vkb(0.,0.,1.);
        gvet Xe=Tracc[0]->Get_Pc1();
        gvet Ve= Tracc[0]->Get_Pc2()-Tracc[0]->Get_Pc1(); Ve.Norma();Ve=Ve.Verso();
   
        Vtxm =  Vert_coor(Vkb,Xkb,Ve,Xe);
        Vtxm.putv(vtx,vty,vtz);
        
// calcolo della massa del neutrino in funzione dello zeta del vertice.    
      Pke=PK-Pelep; 
      Pke.Invar();
//      Vgam.setvn(lkry->M_Hits[pr_gm].Xlab.x-vtx,lkry->M_Hits[pr_gm].Xlab.y-vty,lkry->M_Hits[pr_gm].Xlab.z-vtz);
      Vgam.setvn(lkry->M_Hits[pr_gm].Xlab.x,lkry->M_Hits[pr_gm].Xlab.y,lkry->M_Hits[pr_gm].Xlab.z-vtz);
      
      Vgam=Vgam.Verso();
      Pgam.setqn(Vgam*engm,engm);  
      emisq= Pke.mq -2.*(Pke.e*Pgam.e -Pgam%Pke);
 
// define miske....
     miske=Pke.mom;
     Pke.putv(miskex,miskey,miskez);
     ptrske=Pke.XYNorma();
     
         
         
  //*************************************************************************

  count_wnt++;
  return 1 ;
}


// ========gvet  Vert_coor      def ================

// utility matrix for Vert_coor with general scope
gmatrix PJ1(3,3),PJ2(3,3),PJ(3,3),Noto(3,1),Xvt(3,1);
gmatrix Ident(3,1.0);       // oppure gmatrix Ident(3,3);Ident.Unit();
gvet Vtx;

gvet  Vert_coor(gvet v1,gvet x1,gvet v2,gvet x2)
{
// calcolo del vertice a partire dalle de due rette.

gmatrix P1(x1),P2(x2);  //passa da vettore gvet a matrice(3,1):

//  Note
/*
  attenzione "P.Proij( gvet &v)" definisce la matrice P come  un proiettore
 sul vettore: "|v><v|" 
 Il proiettore applicato ad un vettore, per esempio B, estrae la parte 
 B' parallela a v  come  B'=P*B; 
 [ricorda che anche B è una matrice (3X1) ad una colonna ottenuta per esempio 
 generando la matrice a partire dal vettore tridimensionale gvet "b" con
 gmatrix B(b)] 
 
 Alternativamente  in modo più compatto!
 L'operatore ">>"  usato nella forma "b>>P" (dove b è un gvet e P una gmatrix 3X3 vuota) 
 è equivalente alla Proij(). Attenzione a mettere le parentesi per definire la priorità nelle operazioni! 
*/

PJ1= Ident-(v1>>PJ1); PJ2= Ident-(v2>>PJ2);  // attenzione alle parentesi
 
 PJ1.Det(); PJ2.Det();
 PJ = PJ1+PJ2;     Noto= PJ1*P1+PJ2*P2;
 Xvt=!PJ*Noto; 
 // passo da matrice a gvet.
// Xvt.Get_Vcol(Vtx);   //oppure 
 Vtx=Xvt; // equivale a Vtx.x=Xvt(0,0); Vtx.y=Xvt(1,0); Vtx.z=Xvt(2,0);

 return Vtx;
 
}