flyopunta/src/src_ana/anak4g.cpp

/***************************************************************************
                          anak4g.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

#include "flyoh.h"
#include "anak4g.h"

using namespace std;

AnaK4g::AnaK4g()
{
  tipo=1;nfit=4;titol="K--> 2 pi0 -> 4 gam" ;nome="K4gam";
  nd=0;
  pi0mass=ptr_Pi0->Get_Massa();
  twopi0mass=2.*pi0mass,pi0massq=pi0mass*pi0mass;
  Utot=msqAvo+3.*pi0massq;
  Gout<<"\n Definito il fit "<<titol<<" con K="<<msAvo<<" piomas="<<pi0mass;
}

//============
int AnaK4g::kappa(int i,int j)
//---------------nota on index ...
//define k from 0 to 5 (01,02,03,12,13,23)
//   k= i*3+j    with i<j    no i==j!!!
//   if(i>0)k++;
//   if(i=4)k=15;
//--------------------------------

{ if(i>j)   //swap
  { int jj=j; j=i; i=jj;}
  int k=i*2+j-1;
  if(i==2)k=5;
  //   Gout<<"\n i j k %2d %2d %2d",i,j,k);
  return k;
}



//=========  N e u t r a l    n pi0   F i t =====

//------------------------------------
double  AnaK4g::Distvtx()
{
  //  Decay  K -> pi0 pi0
  // procedure for the  calculation of distance
  // of decay vertice from calorimetr + 
  // the gg masses
  double qij,xi,xj,yi,yj,dx,dy,ei,ej,cxij,cyij;
  double dq=0.0;
  int i=0,j,k=-1;
  Kmax=((nfit-1)*nfit)/2;
  // number indipendente pairings
  
  while(i<nfit)
  {
    ei=egood[i];
    xi=xgood[i];
    yi=ygood[i];

//    Gout<<"\n Dstvtx"<< nfit<<"   "<<idpr[i]<<"  "<<lkry->Pdev[i].e<<" "<<ei<<" "<<xi<<" "<<yi;


    j=i+1;
    while(j<nfit)
    {
      ej=egood[j];
      xj=xgood[j];
      yj=ygood[j];

      //  k index to speed up pairing ..
      //  Kmax= (mishit**2-mishit)/2= elementi off diagonal/2.
      //  for 4 gamma:
      //  define k from 0 to 5 (01,02,03,12,13,23)
      //  for 6 gamma:
      //  define k from 0 to 15 (01,02,03,04,05,12,13,14,15,23,24,25,34,35,45)
      k++;

      dxij[k]=dx=xi-xj;
      dyij[k]=dy=yi-yj;

      if ((dx*dx+dy*dy) < 100.) return -2.;   // attenzione???????? distanza cluster

      qij=ei*ej;
      cxij=exij[k]=qij*dx*dx;
      cyij=eyij[k]=qij*dy*dy;
      dq+=cxij+cyij;
      j++;
    }
    i++;
  }
  distq=dq/msqAvo;
  dist=sqrt(distq);

  Particella *pr;
  pr=SelRea->avo; 

   
  // define the massees
  for(k=0;k<Kmax;k++)
  {
    egij[k]=exij[k]+eyij[k];
    mij[k]=sqrt(egij[k]/distq);

    if(Debugon==1)
    {
      Gout<<"\nEv "<<evento_.Gen.Event<<" Distpiu Masse k "<<k<<" ms="<<mij[k]<<std::endl;
    }
  }
  return dist;
}

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

double AnaK4g::Degij(int k,int i)
{
  // calcola la derivata della prodotto  egij(k)=el*ej*dqlj rispetto
  //alle variabili v(i)= x1,..x4,y1,..y4,e1,..e4
  // matrice per la definizione della metrica...
  static double cij[4][6] =    {{2., 2., 2., 0., 0., 0.},
                                {-2., 0., 0., 2., 2., 0.},
                                { 0.,-2., 0.,-2., 0., 2.},
                                { 0., 0.,-2., 0.,-2.,-2.}} ;

  int ii=i%4;
  double derv= cij[ii][k];
  if(derv==0.0) return 0.0;
  // x1..x4  y1...y4   e1...e4

  if(i<4)      derv*=    exij[k]/dxij[k];
  else if(i<8) derv*=    eyij[k]/dyij[k];
  else if(i>7) derv =    egij[k]/egood[ii];
  return derv;
}
//--------------------------------------------------

void AnaK4g::Ddistq()
{

  // errors   on the (x1,..x4(x6),y1,..y4(y6), e1,..e4(e6))
  int mis2=2*nfit;
  for(int i=0;i<nfit;i++)
  {
    double ep=egood[i];
    //    se= sqrt(0.01 + (0.001024+2.5e-5*ep)*ep);       // 1998
    sv[i+mis2]=sqrt(0.0081 + (0.001024+1.764e-5*ep)*ep);  // 1999
    //    sv[i] =sv[i+nfit] =sqrt(3.6e-3 + 0.001764/ep);
    sv[i] =sv[i+nfit] =sqrt(3.6e-3 + 0.1764/ep);      // correzione 2001

  }

  // derivate di distq rispetto alle variabili
  // nvar allready defined nvar=3*nfit
  //var(i)i (x1...y1...e1...e6) in ddistq(nvar)
  double der=0.0;
  int i=0,k;
  while(i<nvar)
  {
    der=0.0;
    k=0;
    while(k<Kmax)
    {
      der+=Degij(k,i);
      k++;
    }
    dersq[i]=der/msqAvo;   // attenzione per il pionio....

    i++;
  }
}

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

void AnaK4g::Simq(int k1, int k2)
{
  // calcolo sigma sulle variabili del chiq
  // v_i (x1...y1...e1...e4) in ddistq(12)
  double de1,de2,des, def;
  int i=0;                        //??????
  double m1=mij[k1], m2=mij[k2];
  dsumq=0.;
  diffq=0.;
  while(i<12)
  {

    de1 =0.9*sv[i]*(Degij(k1,i)- m1*m1*dersq[i])/(2.*distq*m1);
    de2 =0.9*sv[i]*(Degij(k2,i)- m2*m2*dersq[i])/(2.*distq*m2);

    des=de1+de2;
    def=de1-de2;

    dsumq+=des*des;
    diffq+=def*def;

    //     if(Debugon==1)
    //     Gout<<"\n err %3ld %1d %1d %11.7lf %11.7lf %11.7lf  %11.7lf %11.7lf  %11.7lf",
    //     Event,k1,k2,m1,m2,de1,de2,dv1,dv2);

    i++;
  }

  if(Debugon==1)
    Gout<<"\n S m: "<<m1<<" dsmq="<< dsumq<<" ddfq= "<< diffq<<diffq/dsumq;

}

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

  Gout<<"\n -----> "<<titol<<" Analysis Summary <-------\n";

  Gout<<"\n Fit entries          "<<count_call;
  Gout<<"\n Too few gams         "<<getta_gg;
  Gout<<"\n Too few charge       "<<getta_ch;
  Gout<<"\n Killed by energy     "<<getta_en;
  Gout<<"\n Killed by cog        "<<getta_cog;
  Gout<<"\n Vertex               "<<getta_gf;
  Gout<<"\n Tau cut              "<<getta_tau;
  Gout<<"\n No good fit          "<<getta_fit;
  Gout<<"\n Good events          "<<count_wnt;
  Gout<<"\n\n ======>           d o n e     <============ "<<std::endl;
}
//-------------------------------------------
void AnaK4g::Reset_parm()
{
  nd=0;
  mgg[0]=mgg[1]=mgg[2]=-1.;
  egamtot=cog=cvita=peso=0.0; chisq=3000.;
  zeta=-12020.;
  nextra=0;
  for(int i=0;i<8;i++){e[i]=egood[i]=0.0;}
  nvar=3*nfit;
}
//=========================

int AnaK4g::Fisica()
{
  // neutral procedure  for gammas in lkr calorimeter
  // cuts: acceptance, distance between ganmmas,
  //       energy,
  // calculate ekappa,egamtot,cog, cvita,masse,chisq
  // make cuts and evaluate the best chisq fit
  //----------------------------------------------------
  //----------------- reset----------------------------
  double egx;
  count_call++;
  if(Debugon==1)Write_titol();
  Reset_parm();

   if(tipo==2&&lkry->mhit!=6){getta_gg++;return -2;}
   else if(tipo==1&& lkry->mhit!=4){getta_gg++;return -2;}
//  altri test per altre reazioni

  mishit=lkry->mhit;

  // get coordinates  and xcog
  double xi,yi,xcog=0.0,  ycog=0.0;

  int ngood=-1;
  for(int i=0;i<mishit;i++)
  {
    egx=lkry->M_Hits[i].e_rivela;
    if (egx<=1.5)  nextra++;       // fiorini dixit


    if(egx<3.||egx>100.)   continue;  // giudici 2007 energia test
    xi=lkry->M_Hits[i].Xdev.x;
    yi=lkry->M_Hits[i].Xdev.y;


//      if (xi*xi+yi*yi<225.)           continue     // old ottagono cut
    if ((xi*xi+yi*yi)<215.)               continue;    // ottagono cut
    if (fabs(xi)>116. || fabs(yi)>116.) continue; 
    if (fabs(xi)+fabs(yi) > 164.)       continue; 
    if (fabs(xi)>63.2 && fabs(yi)>84.7) continue; 
    if (fabs(xi)>52.2 && fabs(yi)>84.7) continue; 


    ngood++;
    egood[ngood]=egx;
    xgood[ngood]=xi;
    ygood[ngood]=yi;
    xcog+=egx*xi;
    ycog+=egx*yi;
    egamtot+=egx;
    idpr[ngood]=lkry->M_Hits[i].idm%100;

  }
  // cut on gammas on the base of fit type...
  if( ngood!=nfit-1)      {getta_gg++;return -2;}
  // energy cut
  if(egamtot< 60.||egamtot>180.){getta_en++; return -5;}  // giudici cut 2007

  xcog/=egamtot; ycog/=egamtot;
  cog=sqrt(xcog*xcog+ycog*ycog);
// cog=(fabs(xcog)+fabs(ycog))/2.;
  // cut on cog...for all tipo.
  if(tipo>5 && tipo!=21)
  {if(cog > 100.){ getta_cog++; return -7;}}
  else   if(cog > 10.0 ){ getta_cog++; return -6;}



  //Vertex.  Distance from the calorimeter and dij,mij
  if((dist=Distvtx())<0.0) {getta_gf++;return -3;}

  zeta=lkry->Get_Cface().z-dist;
  if(Debugon==1)
  {
    Gout<<"\n Ana_1: Ev "<<evento_.Gen.Event<<" mis="<<mishit<<" Ek="
    <<egamtot<<" cog="<<cog<<" zv="<<zeta<<std::endl;
  }

  //  time life cut from aks for Kl
  if(tipo==1||tipo==2 )
  {
    static double aks=607.4; // posizione della coincidenza
    double volo=zeta-aks;
   
    
   
    double eta =sqrt(egamtot*egamtot-msqAvo)/msAvo;
    cvita=volo/(eta*2.676); // in ks time_life
    //attenzione taglio sbagliato   Corretto 4.5
    if(cvita<-2.|| cvita>10.) {getta_tau++;return -7;}
    if(cvita>(0.8+0.06*egamtot) ){getta_tau++;return -7;}

    peso= exp( -cvita +volo/(eta*1551.0));
  }



  Ddistq();  // define errors and all the derivatives


  pair[0]= pair[1]= pair[2]=0;

  // choice of the best mass fit
  int stato=MassFit();

  if(stato==1) count_wnt++;

  else getta_fit++;

  return stato;
}
//----------------------------------------------------------
int AnaK4g::MassFit()
{
  int  err=-10;  //set fit error level to -10...
//  chisq=200.01;         // nota ?????????
  chisq=100.01;         // nota ?????????
  //    if(mishit!=4) return err;
  g0=0;g1=1;g2=2;g3=3;
  if(MassChiq()==1) err=1;
  g0=0;g1=2;g2=1;g3=3;
  if(MassChiq()==1) err=1;
  g0=0;g1=3;g2=1;g3=2;
  if(MassChiq()==1) err=1;
  return err;
}


//------------------------------------------------
int AnaK4g::MassChiq()
{
  int err=-1;
  //---------------------------------------
  // masses and chisq  for each permutation...
  //---------------nota on index ...
  //   see procedure  kappa()
  //--------------------------------


  //  masse g0,g1 - g2,g3 - g4,g5


  double ma,mb,chiq;
  double dm,sm;
  int k01,k23;
  k01=kappa(g0,g1);
  k23=kappa(g2,g3);
  ma=mij[k01]; mb=mij[k23];
  Simq(k01,k23);
  dm=(ma-mb); sm=(ma+mb-twopi0mass);

  chiq=(dm*dm)/diffq + (sm*sm)/dsumq;

  // unal chiq

  // Gout<<"\n Ev %10ld Chi %10.3lf pair %4d %4d  %10.3lf %10.3lf ";
  //  if(Debugon==1)Gout<<"\n Ev %10ld Chi %10.3lf pair %4d %4d  %10.3lf %10.3lf ",
  //  Event,chiq,idpr[g0]+100*idpr[g1],idpr[g2]+100*idpr[g3],
  //           ma,mb);

  if(chiq<chisq)
  {
    mgg[0]=ma; mgg[1]=mb;chisq=chiq;
    pair[0]=idpr[g0]+100*idpr[g1];
    pair[1]=idpr[g2]+100*idpr[g3];
    dif1=dm/sqrt(diffq); dif2=sm/sqrt(dsumq);
    if(Debugon==1&&fabs(ma-pi0mass)<0.05)
      Gout<<"\nEV "<<evento_.Gen.Event<<" m:"<< sm<<" "<< dm<<" df1="
      << dif1<<" df2="<< dif2<<" df1/df2"<<dif1/dif2;
    nd=2;
    err=1;
  }
  return err;
}