flyopunta/src/src_ana/anak3pi.cpp

/***************************************************************************
                          anak3pi.cpp  -  description
                             -------------------
    begin                : Sat Jan 18 2003
    copyright            : (C) 2003 by Giuseppe Pierazzini
    email                : peppe@unipi.it
  ***************************************************************************
  *                                                                         *
  *   NA48  simulation program.                                             *
  *                                                                         *
  ***************************************************************************/
#define EPSR     1.0E-200
#define MPIC   0.13956
#define MPICQ  0.019477
#define MASK   0.49367
#define MASKQ  0.24371
#define SS0    0.100714

#include "flyoh.h"
#include "anak3pi.h"

#define MIN(x,y,z) ((x)<(y)?(x)<(z)?(x):(z):(y)<(z)?(y):(z))

using namespace std;

// if only for three tracks ==> one combination
// but for more tarcks => combinations =  k*(k-1)*(k-2)/3! 

  static  int combine[20][3]=  
              { {0,1,2},{0,1,3},{0,2,3},{1,2,3},
                {0,1,4},{0,2,4},{0,3,4},{1,2,4},
                {1,3,4},{2,3,4},{0,1,5},{0,2,5},
                {0,3,5},{0,4,5},{1,2,5},{1,3,5},
                {1,4,5},{2,3,5},{2,4,5},{3,4,5}};
// for n tracks  use the  combinations
//     3                     1
//     4                     1-4
//     5                     1-10
//     6                     1-20
//-----------------------------------------------
//-----------------------------------------------

   AnaK3pi::AnaK3pi()
   {
    tipo=4;nfit=3;titol="K--> 3 pio_ch";nome="K3pi";
    nd=np=0;
    maspi=.13956;
    Ird0=Ird1=Ird2=0;
    Icb0=Icb1=Icb2=0;
    Gout<<"\n Attivato il fit di "<<titol;
    if(Le_Tracce==0)
     {
       Gout<<"\n Error.. K3pi analisys ... Tracks non defined!"<<std::endl;  
       exit(0); 
     }
   }
   AnaK3pi::~AnaK3pi(){};

//=============== C h a r g e   t r a c k   F i t ========
int AnaK3pi::Fisica()
  {


  np=nd=0;
  np=Le_Tracce->Tot_trk+1; 
//  Gout<<"\n Nhit: "<<np;
  if(np<3||np>6){np=0;getta_ch++;return(-1);}
  nd=3;
  carica=.0;
 // define the coupling number (1,4,10,20)
   int trial=1;
   if(np>3)trial=(np*(np-1)*(np-2))/6;


   chisq=300.;
   for (int i=0;i<trial;i++)
   {
    count_call++;
// get track index 
   Icb0=combine[i][0];
   Icb1=combine[i][1];
   Icb2=combine[i][2];
   
  double etot= Tracc[Icb0]->Get_P()[3]+Tracc[Icb1]->Get_P()[3]+Tracc[Icb2]->Get_P()[3];
  
  if(etot<10. ){getta_en++;continue;}     //return -2;};
  carica=Tracc[Icb0]->charge+Tracc[Icb1]->charge+Tracc[Icb2]->charge;
  // get the index order...Ird0,Ird1 (same charge as above) Ird3 is the odd charge 
   if( carica!=1&&carica!=-1 ) {getta_char++;continue;}

  Ordina(); // ordina le tracce da essere (++-) o (--+)

  P1=Tracc[Ird0]->Get_Pc1();
  P2=Tracc[Ird1]->Get_Pc1();
  P3=Tracc[Ird2]->Get_Pc1();
  if(Vertice()<1){getta_fit++;continue;} //return -3;};
  if(chicomb<chisq)
    {chisq=chicomb;
     charge=carica;
     PVsave=PV;
     Isv0=Ird0; Isv1=Ird1;Isv2=Ird2;
    }
   }
  if(chisq>200.)return -3; // see Vertice()

  if(trial>1)   // get the saved values by the best fit for trial >1
    {Ird0=Isv0;Ird1=Isv1;Ird2=Isv2;
     P1=Tracc[Ird0]->Get_Pc1();
     P2=Tracc[Ird1]->Get_Pc1();
     P3=Tracc[Ird2]->Get_Pc1();
     PV=PVsave;
    }


  
  pmis[0]=Tracc[Ird0]->ptrk;
  pmis[1]=Tracc[Ird1]->ptrk;
  pmis[2]=Tracc[Ird2]->ptrk;

// calcolo il cog at lkr
  Cog3p=(Tracc[Ird0]->Xlk*pmis[0]+Tracc[Ird1]->Xlk*pmis[1] +Tracc[Ird2]->Xlk*pmis[2]);
  Cog3p/=(pmis[0]+pmis[1]+pmis[2]);
  kcog  = Cog3p.XYNorma();
  cogm1x=Cog3p[0];
  cogm1y=Cog3p[1];     
  cogm1 =kcog;
  PV.putv(vx,vy,vz);
            
  V1= P1-PV; V1.Norma();V1=V1.Verso()*pmis[0];
  V2= P2-PV; V2.Norma();V2=V2.Verso()*pmis[1];
  V3= P3-PV; V3.Norma();V3=V3.Verso()*pmis[2];
  
// calcolo della massa invariante
  Q1.setvn(V1,maspi);
  Q2.setvn(V2,maspi);
  Q3.setvn(V3,maspi);
  QT=Q1+Q2+Q3;
  mkmis=QT.Invar();
  pkmis=QT.norma;

  ptksq=(QT&PV).Norma()/PV.norma;        
// Calcolo dei parametri s1,s2,s3,su..
// con rinomalizzazione della massa invariante a quella del k  !!
//
// qvet QC = dQT/dcorr  quadrivettore derivato rispetti alla correzione
 
   QC.setqn(V1+V2+V3,  Q1.normaq/Q1.e+Q2.normaq/Q2.e+Q3.normaq/Q3.e);
//   QC.print("QC");
   double corr=1.- (QT.m-MASK)*QT.m/(QC*QT);  
// Gout<<"\n Correzione = %12.7lf  %12.5lf ",corr,QT*QC);   
   Q1.setvn(V1*corr,maspi);
   Q2.setvn(V2*corr,maspi);
   Q3.setvn(V3*corr,maspi);
   QT=Q1+Q2+Q3;  QT.Invar();
//   QT.print("QTcorr");
  fi3=atan2(Q3[1],Q3[0]);

  ss[0]=(QT-Q1).Invarq();
  ss[1]=(QT-Q2).Invarq();
  ss[2]=(QT-Q3).Invarq();   
  su = (ss[2]-SS0)/MPICQ;  // u
  sv = 0.433013*(ss[1]-ss[0])/MPICQ;  // v  
  count_wnt++;
  return(1);

}

//======================================================================================
 int AnaK3pi::Vertice()
   {
// calcola il vertice ed il chiq relativo
//     gvet P1,P2,P3,V1,V2,V3;
//     gvet Noto,PV;
     
     V1= Tracc[Ird0]->Get_Pc2()-P1; V1.Norma();V1=V1.Verso();
     V2= Tracc[Ird1]->Get_Pc2()-P2; V2.Norma();V2=V2.Verso();
     V3= Tracc[Ird2]->Get_Pc2()-P3; V3.Norma();V3=V3.Verso();
     Noto=P1+P2+P3;
     Noto-=(V1*(V1%P1)+V2*(V2%P2)+V3*(V3%P3));
//     Noto.print("Noto");

// matrice righe MV1,MV2,MV3
     gvet MV1(3.,0.,0.), MV2(0.,3.,0.),MV3(0.,0.,3.);
     MV1-=V1*V1[0] + V2*V2[0] + V3*V3[0];
     MV2-=V1*V1[1] + V2*V2[1] + V3*V3[1];
     MV3-=V1*V1[2] + V2*V2[2] + V3*V3[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
// nota:  definizione del chiq
//   chiq = somma distanze quadrate del vertice dalle rette.
//   per normalizzare , ciascuna distamza va divisa per il suo
//   errore. qui non si normalizza..
//   Ricordo comunque che l'errore deriva dagli errori nelle camere 1 e 2
//   quindi sul punto nella camera 1 e nel calcolo del versore di ciascuna
//   retta che dipende dai punti nelle due camere. Prendo gli errori x e y nelle
//   camere tutti uguali, segue:
//    Err=  Sgma*(sqrt( (1+dz/dc)**2 + (dz/dc)**2)   con dz distanza del
//   vertice dalla camera 1(~90 mt) e dc distanza tra le camere 1 e 2 (~11 mt ).
//    approssimativamente   err=sigm*8.2*sqrt(2)
//    con un errore nelle camere di 0.012 cm  ==> err= 0.14 cm

//     gvet DP1,DP2,DP3;
     DP1=PV-P1,DP2=PV-P2,DP3=PV-P3;
     double t1=V1%DP1, t2=V2%DP2, t3=V3%DP3;
     DP1-=V1*t1; DP2-=V2*t2; DP3-=V3*t3;
     DP1.Norma(); DP2.Norma(); DP3.Norma();
     chicomb= DP1.normaq + DP2.normaq+ DP3.normaq;

     if(Debugon==1) {
     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]);
                 */
   
     if(chicomb>200.)return -3;
     return 1;
   }

 //----------------------------------------------------
   void AnaK3pi::print_scale()
    {
     Gout<<"\n -----> "<<titol<<" Analysis Summary <-------\n";
     Gout<<"\n Fit entries          "<<count_call;
     Gout<<"\n Tracks <3 or >6      "<<getta_ch;
     Gout<<"\n Wrong charge         "<<getta_char;
     Gout<<"\n Killed by energy     "<<getta_en;
     Gout<<"\n No good vertex       "<<getta_fit;
     Gout<<"\n Good events          "<<count_wnt;
     Gout<<"\n\n ======>           d o n e     <============ "<<std::endl;
    }
//-------------------------------------------
  void AnaK3pi::Ordina()
  {
   Ird0=Icb0; Ird1=Icb1; Ird2=Icb2;
   if(Tracc[Icb0]->charge!=charge) {Ird0=Icb2; Ird1=Icb1; Ird2=Icb0;}
   else if(Tracc[Icb1]->charge!=charge) {Ird0=Icb0; Ird1=Icb2; Ird2=Icb1;}
   
  }