flyopunta/src/src_ana/anak2pi.cpp

/***************************************************************************
                          anak2pi.cpp  -  description
                             -------------------
    begin                : Sat Jan 18 2003
    copyright            : (C) 2003 by Giuseppe Pierazzini
    email                : peppe@unipi.it
  ***************************************************************************
  *                                                                         *
  *   NA48  simulation program.                                             *
  *                                                                         *
  ***************************************************************************/
  // #define MIN(x,y,z) ((x)<(y)?(x)<(z)?(x):(z):(y)<(z)?(y):(z))
  
#include "flyoh.h"
#include "anak2pi.h"
using namespace std;

//-------------------------------------------------------------
//-------------------------------------------------------------

 AnaK2pi::AnaK2pi()
   {tipo=3;nfit=2;titol="K--> pi+ pi-";nome="K2pi";
    nd=np=0;
    Gout<<"\n Attivato il fit di "<<titol;
   }                  
                                                
  AnaK2pi::~AnaK2pi(){};

 //=============== C h a r g e  2 pions    F i t ========

int AnaK2pi::Fisica()
  {

//  Device *dev;
//  gvet Xa,Va,Vb;
  qvet qp[3],qq,qtot;   // crea e azzera
  np=nd=0;
  count_call++;
//  
  np=Le_Tracce->Tot_trk+1;
  if(np!=2) {np=0;return(-1);}
 // attenzione per piu' di tre tracce si dovrebbe iterare sulla meglio coppia
//  vedi anak3pi.cpp 
  charge=Tracc[0]->charge+Tracc[1]->charge;   
  pmis[0]=Tracc[0]->ptrk;
  pmis[1]=Tracc[1]->ptrk;
  qp[0]=  Tracc[0]->Get_P();
  qp[1]=  Tracc[1]->Get_P();
  qtot=qp[0]+qp[1]; 
  qtot.Invar();
  mkmis=qtot.m;
  pkmis=qtot.norma;

   if(pkmis<50.||pkmis>190.){getta_en++;return -2;}
   if(pmis[0]<10.){getta_ch++; return -2;}
   if(pmis[1]<10.){getta_ch++; return -2;}

   nd=2;
   qq=qp[0]+qp[1];
   qq.Invar();qq.Norma();
   mkmis=qq.m;
   pkmis=qq .norma;
   mcut=0.1608e-2+qq[3]*0.8505e-5;
   mcut=fabs(mkmis-0.49367)/mcut;

   if(Vertice()!=1){getta_fit++;return -3;}
   double lvita=vz-607.5;   // - Aks posizione...
   cvita= lvita/(2.676*pkmis/.49367);
  
   double cogx= (pmis[0]*Cam4->M_Hits[0].Xdev.x+pmis[1]*Cam4->M_Hits[1].Xdev.x)/pkmis;
   double cogy= (pmis[0]*Cam4->M_Hits[0].Xdev.y+pmis[1]*Cam4->M_Hits[1].Xdev.y)/pkmis;
 
   cog=sqrt(cogx*cogx+cogy*cogy);
   asp= fabs((pmis[1]-pmis[0])/pkmis);
// calcolo del  mon trsverso rispetto all'asse del fascio.
   gvet Bvers(0.,0.,1.);           //asse KL
   if(SelRea->avo->Get_Idm()!=16)
     {Bvers.setvn(0.,-7.19,12171.4); // asse Ks
      Bvers=Bvers.Verso();
     }
    gvet Pt=Bvers&qq; Pt.Norma();
    cptq= Pt.normaq;
    count_wnt++;
   return 1;

}

//======================================================================================
 int AnaK2pi::Vertice()
   {
// calcola il vertice ed il chiq relativo  per due tracce
     gvet P1,P2,V1,V2;
     gvet Noto,PV;
         
     V1= Tracc[0]->Get_Pc2()-(P1=Tracc[0]->Get_Pc1()); V1.Norma();V1=V1.Verso();
     V2= Tracc[1]->Get_Pc2()-(P2=Tracc[1]->Get_Pc1()); V2.Norma();V2=V2.Verso();
     
     Noto=P1+P2;
     Noto-=(V1*(V1%P1)+V2*(V2%P2));
//     Noto.print("Noto");

// matrice righe MV1,MV2,MV3
     gvet MV1(2.,0.,0.), MV2(0.,2.,0.),MV3(0.,0.,2.);
     MV1-=V1*V1[0] + V2*V2[0];
     MV2-=V1*V1[1] + V2*V2[1];
     MV3-=V1*V1[2] + V2*V2[2];

     gvet VM1,VM2,VM3;    //righe della matrice inversa
     VM1= MV2&MV3;        // nota: prodotto vettoriale
     double det = MV1%VM1; // determinante
     VM1= VM1/det;
     VM2= (MV3&MV1)/det;
     VM3= (MV1&MV2)/det;
     PV.setvn(VM1%Noto,VM2%Noto,VM3%Noto);  // vertice
//calcolo il chiq
     gvet DP1=PV-P1,DP2=PV-P2;
     double t1=V1%DP1, t2=V2%DP2;
     DP1-=V1*t1; DP2-=V2*t2;
     DP1.Norma(); DP2.Norma();
     chisq= DP1.normaq + DP2.normaq;

     if(Debugon>0) {
     PV.print("Vert");
     Gout<<"\n Vertice chiq ="<<chisq<<std::endl;
     }
     /*     PV.setv(cam2->mx[icp]-cam1->mx[icp],
                 cam2->my[icp]-cam1->my[icp],
                 cam2->mz[icp]-cam1->mz[icp]);
                 */

     PV.putv(vx,vy,vz);
//     Gout<<"\n vertex %10.5lf  %10.5lf   %10.5lf", vx,vy,vz);
     if(chisq>200.)return -1;
     return 1;
   }

 //----------------------------------------------------
   void AnaK2pi::print_scale()
    {

     Gout<<"\n -----> "<<titol<<"   Analysis Summary. <-------\n";
     Gout<<"\n Fit entries          "<<count_call;
     Gout<<"\n Track energ to low   "<<getta_ch;
     Gout<<"\n Killed by energy     "<<getta_en;
     Gout<<"\n Tau cut              "<<getta_tau;
     Gout<<"\n No good vertex       "<<getta_fit;
     Gout<<"\n Good events          "<<count_wnt;
     Gout<<"\n\n ======>           d o n e     <============ ";
    }
//-------------------------------------------