flyopunta/src/src_ana/anak6g.cpp

/***************************************************************************
                          anak6g.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 "anak6g.h"
#include "gmatrix.h"


AnaK6g::AnaK6g()
{
  tipo=2;nfit=6;titol="K--> 3 pi0 -> 6 gam" ;nome="K6gam";
  Gout<<"\n Ridefinito e creato il fit di "<<titol;
  // gmatrix set size;
  Bm64.Set_Size ( 6,4 );Az41.Set_Size ( 4,1 );
  F41.Set_Size ( 4,1 ); R41.Set_Size ( 4,1 );
  De61.Set_Size ( 6,1 );Lmd41.Set_Size ( 4,1 );
  Zm11.Set_Size ( 1,1 );
  Var66.Set_Size ( 6,6 ); InV66.Set_Size ( 6,6 );

  BVB.Set_Size ( 4,4 );
  GB.Set_Size ( 4,4 );
  AGBA.Set_Size ( 1,1 );
  Gx.Set_Size ( 1,1 );
  CHQ.Set_Size ( 1,1 );

}
AnaK6g::~AnaK6g() {}

void AnaK6g::Reset_parm()
{
  AnaK4g::Reset_parm();
  permut_flag=iter_flag=permut_fit=iter_fit=-1;

  dchi2=2000.;
  chisq=101.;
  ss6[0]=ss6[1]=ss6[2]=0.0;
  ws6[0]=ws6[1]=ws6[2]=1.;
}


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

int AnaK6g::Fisica()
{
  
  if ( Debugon==1 ) Write_titol();

  count_call++;
  Reset_parm();
  ngm = mishit=lkry->mhit;
  if ( WrtNt>1 ) Reset_Dummy_Data();  // se si scrivono tutti gli eventi

// also 7 gammas for dalitz....
  if ( mishit<6 ) {getta_gg++;dead6=-1;return -1;}
//if( mishit>6 ) {getta_gg++;dead6=-1;return -1;}   //too many gammas..

  double xi,yi,egx;
  bmax=0.;egamtot=0.;
  int ngood=-1;
  for ( int i=0;i<mishit;i++ )
  {

    xi=lkry->M_Hits[i].Xdev.x;
    yi=lkry->M_Hits[i].Xdev.y;
    egx=lkry->M_Hits[i].e_rivela;


    if ( egx<3.||egx>100. ) {nextra++; continue;}  // giudici 2007 energia test

    // goemetrical cut

    double braccio=sqrt ( xi*xi+yi*yi ); if ( braccio>bmax ) bmax=braccio;
    if ( braccio<15. )               continue;
//    if ( braccio>116. )             continue;   // circolar 116.

    double absx=fabs ( xi ), absy=fabs ( yi );
    if ( absx >116. || absy >116. ) continue;
    if ( absx +absy > 164. )       continue; // rectangular  116.*sqrt(2)
    if ( absx >63.2 && absy >84.7 ) continue;
    if ( absx >52.2 && absy >94.7 ) continue;
    if ( pow ( ( absx-63.2 ),2 ) +pow ( ( absy-94.7 ),2 ) < 121. )  continue;  //octagonal

// save the good gammas in _sav
    ngood++;
    eg_sav[ngood]  =egx;
    xg_sav[ngood]  =xi;
    yg_sav[ngood]  =yi;
    zg_sav[ngood]  =lkry->M_Hits[i].Xlab.z;
    idpr_sav[ngood]=lkry->M_Hits[i].idm;
    egamtot+=egx;
  }


  // cut on gammas  on the base of fit type..here 6 [5 in C++]
  if ( ngood<5 ) {getta_gg++;dead6=-3;return -3;}
  if ( egamtot< 60. ) {getta_en++; dead6=-4;return -4;}


  int Debug_fit=0;  // per test se ==1
  if ( Debug_fit==1 )
  {
    int ngboni=ngood+1;
    Gout<<setprecision ( 3 );
    Gout<<"\n Start   "<<evento_.Gen.Event<<" whit "<<ngboni
    <<" gam  e zv "<<setw ( 8 ) <<evento_.Gen.zv <<" etot "<<setw ( 8 ) <<egamtot
    <<" Lkrz "<<lkry->Get_Cface().z;
    for ( int i=0;i<ngboni;i++ )
      Gout<<"\nG xyze "<<setw ( 8 ) <<xg_sav[i]<<" "
      <<setw ( 8 ) <<yg_sav[i]<<" "
      <<setw ( 8 ) <<eg_sav[i]<<" Idm "<<idpr_sav[i];
    Gout<<endl;
  }

//    FITTING   begin
  nd=0;
  chisq=15.;

// se ci sono sette gamma si itera sulle possibili 6-ple (sette ipotesi)
  int ik =0, ipotesi=1, iok=0;
  if ( ngood==6 ) ipotesi=7;
  while ( ik< ipotesi )
  {
    if ( Get_Data_from_saved ( ik ) <0 ) { ik++;continue;} // select the 6-pla

// first large cuts on neasured values....just to speed up the fitting...gmatrix
    if ( egamtot< 60.||egamtot>180. ) {ik++;continue;}
    if ( cog > 100.0 ) {ik++;continue;}
    ik++;

// debug ======gmatrix
    if ( Debug_fit>=2 )
    {
      for ( int i=0;i<6;i++ )
        Gout<<"\nC xyze "<<setw ( 8 ) <<xgood[i]<<" "
        <<setw ( 8 ) <<ygood[i]<<" "
        <<setw ( 8 ) <<egood[i]<<" Idm "<<idpr[i];
      Gout<<endl;
    }
// end debug =====

//    procedure    F I T T A


    if ( Fitta() ==1 ) iok=1;
  }
  if ( iok==0 ) return -1; // non ci sono fit ok...
  // debug =====

  if ( Debug_fit==1 )
  {
    Gout<<setprecision ( 3 );
    Gout<<"\n Fit ev   "<<evento_.Gen.Event<<" ipot. "<<ik<<" itr "<<iter_fit
    <<" chq "<<setw ( 8 ) <<chisq<<" zeta "<<setw ( 8 ) <<zeta <<" et "<<setw ( 8 ) <<egamt_fit;

    for ( int i=0;i<6;i++ )
      Gout<<"\nF xyze "<<setw ( 8 ) <<xg_fit[i]<<" " <<setw ( 8 ) <<yg_fit[i]
      <<setw ( 8 ) <<eg_fit[i]<<"  Perm "<<permut_fit<<" Idm "<<idpr_fit[i];
    Gout<<endl;

  }

//  ------------fit done ----
  if ( egamt_fit<70. ||egamt_fit>170. ) {getta_en++; dead6=-4;return -4;}

  if ( chisq>15. )  {getta_fit++;dead6=-6;return -6;} // wrong chiq in the fit!
  if ( dist_fit<0.0 ) {getta_gf++;dead6=-7;return -7;} // wrong vx dist_fit
//  if (cog_fit>20.0) {getta_cog++; dead6=-5;return -5;} // wrong cog_fit
  if ( cog_fit>100.0 ) {getta_cog++; dead6=-5;return -5;} // wrong cog_fit  19.07.09
//  time life cut from aks for Kl
//  aks=607.2; // posizione della coincidenza
  double volo=zeta-607.2;

//massa Ks=0.497648;
  cvita= ( ( volo*0.497648 ) /egamt_fit ) /2.6842;    //19.07.09
//  if ( cvita<-2.|| cvita>10. ) {getta_tau++;dead6=-8;return -8;}   //no il 19.07.09
  if ( cvita> ( 0.8+0.06*egamt_fit ) ) {getta_tau++;dead6=-9;return -9;}

// define cinematic variables

  for ( int i=0;i<6;i++ )
  {
    Versg.setvn ( xg_fit[i],yg_fit[i],dist_fit ); Pgam=Versg.NVerso() *eg_fit[i];
    Gq[i].setvn ( Pgam,0.0 );
  }


  Pio01=Gq[0]+Gq[1];mgg[0]=  Pio01.Invarq();
  Pio23=Gq[2]+Gq[3];mgg[1]=  Pio23.Invarq();
  Pio45=Gq[4]+Gq[5];mgg[2]=  Pio45.Invarq();
  
  if(mgg[0]>0.)mgg[0]=sqrt(mgg[0]); else mgg[0]=-sqrt(fabs(mgg[0]));
  if(mgg[1]>0.)mgg[1]=sqrt(mgg[1]); else mgg[1]=-sqrt(fabs(mgg[1]));
  if(mgg[2]>0.)mgg[2]=sqrt(mgg[2]); else mgg[2]=-sqrt(fabs(mgg[2]));
  
  qvet Pm12=Pio23+Pio45;
  qvet Pm02=Pio01+Pio45;
  qvet Pm01=Pio23+Pio01;

  ss6[0]= Pm12.Invarq();
  ss6[1]= Pm02.Invarq();
  ss6[2]= Pm01.Invarq();
  
  vss6[0]= fabs(ss6[1]-ss6[2])*0.577350269;//1./sqrt(3.)
  vss6[1]= fabs(ss6[2]-ss6[0])*0.577350269;
  vss6[2]= fabs(ss6[0]-ss6[1])*0.577350269;

  int Irand=evento_.Gen.Irand_ss;

  if ( Irand==0 )
  {
    ss6_sort=ss6[0];  
    vss6_sort=vss6[0];
  }
  else if ( Irand==1 )
  {
    ss6_sort=ss6[1];  
    vss6_sort=vss6[1];
  }
  else if ( Irand==2 )
  {
    ss6_sort=ss6[2];  
    vss6_sort=vss6[2]; 
  }
  
  
  

//   Pio01.print("p01");  Gout<<setprecision(6)<<" "<<Pio01.m<<"  pi0m "<<pi0mass<<" " <<4.*pi0massq;
//   Pio23.print("p23");  Gout<<setprecision(6)<<" "<<Pio23.m<<" ss6[0] "<<ss6[0]<<" " <<evento_.Gen.ss0;
//   Pio45.print("p45");  Gout<<setprecision(6)<<" "<<Pio45.m<<" dif "<<ss6[0]-evento_.Gen.ss0;


// calcolo il peso per ogni ss
  double wpeso;
  wpeso=Fspace ( 0.49767,Pm12.m, 0.13497 )  *Fspace ( Pm12.m, 0.13497,0.13497 );
  if ( wpeso>0.0 ) ws6[0]=1./wpeso;
  wpeso=Fspace ( 0.49767,Pm02.m, 0.13497 ) *Fspace ( Pm02.m, 0.13497,0.13497 );
  if ( wpeso>0.0 ) ws6[1]=1./wpeso;
  wpeso=Fspace ( 0.49767,Pm01.m, 0.13497 ) *Fspace ( Pm01.m, 0.13497,0.13497 );
  if ( wpeso>0.0 ) ws6[2]=1./wpeso;
  //-------------------------
  count_wnt++; dead6=0;
  return 1;
}


void AnaK6g::Reset_Dummy_Data()
{
  chisq=100.;
  zeta=SelRea->avo->X.z;
  cvita= ( ( zeta-607.5 ) *msAvo/evento_.Gen.pk ) /2.6842;
  egamt_fit=egamtot;
  for ( int i=0;i<mishit;i++ ) eg_fit[i]=lkry->M_Hits[i].e_rivela;
  for ( int i=mishit;i<6;i++ ) eg_fit[i]=-1.0;
  dead6=0;
}






//----------------------------------------------------------
int AnaK6g::Fitta()
{

  double Fsum,chiq=0.0;
  double Fsum_cut=0.00001;
//  double Fsum_cut=0.00001;
  //-----------------  l o o f i t -------
  int iter,ipr;
  int iok=0;
  last_chiq=100.;
  int debugfit=0;
  if ( evento_.Gen.Event>1&&evento_.Gen.Event<20 ) debugfit=1;

  for ( ipr=0;ipr<15;ipr++ )
  {
    Permutazione ( ipr );// get the right permutated parameters (0-1,2-3,4-5)

    InV66=!Var66;
    De61.Reset();

    if ( ( Fsum=Calcolo_Vincoli() ) >1.0 ) {chiq=999.;continue;}

    for ( iter=0;iter<20;iter++ )
    {
      //   Gout<<"\n\n Start  Evento :  "<< evento_.Gen.Event<< "  ipr,iter: "<< ipr<<" "<<iter;
      //  calcolo i vincoli e le derivate (vincoli troppo alti return 10.)

//      if ( ( Fsum=Calcolo_Vincoli() ) >1.0 ) {chiq=999.;break;}


      // algebra di fit.. 9 variabili ed un parametro;
      BVB = ( ~Bm64 ) * ( Var66*Bm64 );

      GB= !BVB;
      if ( GB.okflag<0 )
      {
        Gout<<"\n BVB Okflag error  Ev "<<evento_.Gen.Event;
        break;
      }

      AGBA = ~Az41*GB*Az41;

      //if(evento_.Gen.Event==1392&& ipr==0){AGBA.Det();AGBA.Print("AGBA");}

      Gx=!AGBA;
      if ( Gx.okflag<0 )
      {
        Gout<<"\n AGBA Okflag error  Ev "<<evento_.Gen.Event;
        break;
      }

      Zm11= -Gx*~Az41*GB*R41;
      Lmd41=GB* ( Az41*Zm11+R41 );
      gmatrix daux=Bm64*Lmd41;
      De61=-Var66*daux;

      // updating
      dist_vertex+= Zm11 ( 0,0 );
      distq=dist_vertex*dist_vertex;
      for ( int i =0;i<6;i++ ) ec[i]=De61 ( i,0 ) +em[i];

      //  test of iteration end.
      CHQ= ~De61*InV66*De61;
      chiq=CHQ ( 0,0 );
      if ( iter>0&&chiq>100. ) { break; } // skip very  bad permutation fit
      if ( ( Fsum=Calcolo_Vincoli() ) >1.0 ) {chiq=999.;break;}// skip very  bad permutation fit


// test to save good data...
      if ( Fsum<Fsum_cut &&iter>1 )
      {
        dchi2=abs(last_chiq-chiq);

        if ( dchi2 <20. )
        {
 /*         Gout<<setprecision ( 5 );
          Gout<<"\n Fitta: "<<evento_.Gen.Event<<std::scientific
          <<" chiq " <<setw ( 12 ) <<chiq <<" diffq "<<setw ( 12 ) <<dchi2<<" " <<setw ( 12 )
          <<Fsum<<" ipr " <<ipr<< " itr " << iter;
*/

          return -1;  // Si accettano gli eventi in cui l'accoppiamento non è degenere!!!!
        }
        if ( chiq>chisq ) break;  // if worse than the last one


// save the best data


        chisq=chiq;
        permut_flag = ipr;
        iter_flag=iter;
        Save_good_fit();
        last_chiq=chisq;
        iok=1;



// debug....====

        if ( debugfit==1 )
        {

          double Smas=0.0;
          for ( int i=0;i<5;i++ )
          {
            for ( int j=i+1;j<6;j++ )
              { Smas+= ec[i]*ec[j]*dq[i][j];}
          }
          Gout<<setprecision ( 5 );
          Gout<<"\n Fitta out  Ev "<<evento_.Gen.Event<< " chq "<<chisq <<" Fsum: "<< Fsum
          <<" ipr "<<ipr<< " itr "<<iter
          <<"  mas_k  "<< sqrt ( Smas/distq ) <<endl;
//          for ( int i=0;i<6;i++ ) Gout<<" [ "<<em[i]<<" "<<ec[i]<<" ]";Gout<<endl;
        }
        break;

      }
    }  //iter end....



  } // end permutazioni
  return iok;
}


double AnaK6g::Calcolo_Vincoli ()
{
  //atenzione si usano le matrici.....
  // dei vincoli e delle sue derivate.
  static double  deltz= 0.01,deltg=0.0001;
  double Fsum=0.0;
  int i4,ifu;
  for ( i4=0;i4<4;i4++ )
  {
    F41 ( i4,0 ) = F_Vinc ( i4,ec,dist_vertex,-1 ,0.0 );
    Fsum+=abs ( F41 ( i4,0 ) );
  }
  if ( Fsum>1. ) return 10.;  // vincoli non soddisfatti
  // derivate
  for ( i4=0;i4<4;i4++ )
  {
    Az41 ( i4,0 ) =0.5* ( F_Vinc ( i4,ec,dist_vertex,10,deltz )-F_Vinc ( i4,ec,dist_vertex,10,-deltz ) ) /deltz;
    for ( int i=0; i<6;i++ )
    {
      ifu=i;
      Bm64 ( i,i4 ) =0.5* ( F_Vinc ( i4,ec,dist_vertex,ifu,deltg )-F_Vinc ( i4,ec,dist_vertex,ifu,-deltg ) ) /deltg;
    }
  }
  R41= ( F41- ( ~Bm64 ) *De61 );// residui ...algebra matriciale....in gmatrix.h
  return Fsum;
}
//------------------------------------
double AnaK6g::F_Vinc ( int i4,double *ec ,double dist,int id,double delta )
{
  // procedura per calcolare i valori dei vincoli
  double zd,ecd[6];
  double  zdq,fsum=0.0;

  // ----------------------------
  // ridefisce i parametri per il cacolo anche delle derivate!
  for ( int i=0;i<6;i++ )  // energia e zeta
  {
    ecd[i]=ec[i];
    if ( id==i ) ecd[i]+=delta;
  }
  //la zeta del vertice
  zd=dist;
  if ( id==10 ) zd+=delta;
  zdq=zd*zd;
  //--------------------------------
  // calcolo dei vincoli normalizzati alla massa.
  switch ( i4 )
  {
    case 0:
      fsum= ecd[0]*ecd[1]*dq[0][1];
      return ( fsum/zdq-pi0massq );
//      return abs ( fsum/zdq-pi0massq );  ma perchè!!!!!!!!!!!!!!!!!!!buttare.....

      break;
    case 1:
      fsum= ecd[2]*ecd[3]*dq[2][3];
      return ( fsum/zdq-pi0massq );
//      return abs ( fsum/zdq-pi0massq );

      break;
    case 2:
      fsum= ecd[4]*ecd[5]*dq[4][5];
      return ( fsum/zdq-pi0massq );
//      return abs ( fsum/zdq-pi0massq );
      break;
    case 3:
      fsum=0.0;
      for ( int i=0;i<5;i++ )
      {
        for ( int j=i+1;j<6;j++ )
          { fsum+= ecd[i]*ecd[j]*dq[i][j];}
      }
//    if(evento_.Gen.Event==1392)
//      Gout<<"\n F_Vinc "<<setprecision(5)<<" s "<<fsum <<" z "<<zd<<" mkq "<< msqAvo <<"
//   "<<fsum/(zdq)<<" " <<id;
      return ( fsum/zdq -msqAvo );
//      return abs ( fsum/zdq -msqAvo );
      break;
    default:
      printf ( "errore in F_Vinc\n" ); exit ( -1 );
      break;
  }
  return 0.0;
}
//*****
void AnaK6g::Permutazione ( int ipr )
{
  // perumtazion ipr (0-14);
  // ricava le informazioni di base.
  //  Ditanze tra i gamma e previosione di dist_vertex
  int jc;
  int * Perm=Get_Permuta ( ipr );
  double corr_errorq= 1.21;
//    se= sqrt(0.01 + (0.001024+2.5e-5*ep)*ep);       // 1998
//    se= sqrt ( 0.0081 + ( 0.001024+1.764e-5*ep ) *ep );      //1999

  for ( int i =0;i<6;i++ )
  {
    jc=Perm[i];
    em[i]=ec[i]=egood[jc];
    Var66 ( i,i ) = ( 0.0081 + ( 0.001024+1.764e-5*egood[jc] ) *egood[jc] ) *corr_errorq;

    //  Var66 ( i,i ) = ( 0.01 + ( 0.001024+2.5e-5*egood[jc] ) *egood[jc] );


    xg[i]=xgood[jc];
    yg[i]=ygood[jc];

  }

  //Distanze tra gli hits dei gamma.
  double dx,dy;
  for ( int i =0;i<6;i++ )
  {
    for ( int j =i+1;j<6;j++ )
    {
      dx=xg[i]-xg[j];
      dy=yg[i]-yg[j];
      dq[i][j]=dq[j][i]=dx*dx+dy*dy;
    }
  }

  // inizializzo la distanza del vertice dal cal.
  double fsum=0.0;
  for ( int i=0;i<6-1;i++ )
  {
    for ( int j=i+1;j<6;j++ )
      { fsum+= ec[i]*ec[j]*dq[i][j];}
  }
  dist_vertex=sqrt ( fsum/msqAvo );


}
//---------------------------------
int * AnaK6g::Get_Permuta ( int ipr )
{
  // permutation indexes...
  static int  iperm[15][6] =
  {
    {0,1,2,3,4,5},{0,1,2,4,3,5},{0,1,2,5,4,3},
    {0,2,1,3,4,5},{0,2,1,4,3,5},{0,2,1,5,4,3},
    {0,3,2,1,4,5},{0,3,2,4,1,5},{0,3,2,5,4,1},
    {0,4,2,3,1,5},{0,4,2,1,3,5},{0,4,2,5,1,3},
    {0,5,2,3,4,1},{0,5,2,4,3,1},{0,5,2,1,4,3}
  };

  return iperm[ipr];
}
//-------------------------
int AnaK6g::Get_Data_from_saved ( int k )
{
// sceglie i gamma in una lista ordinata di 7 gamma
// parte dal gamma  k e copia i primi sei i successione tornando a zero dopo il settimo..
  int kora=k;
  double xi,yi,xcog=0.0,ycog=0.0,egx;
  egamtot=0.;
  for ( int i=0;i<6;i++ )
  {
    kora=kora%7;
    egood[i]=egx=eg_sav[kora];
    xgood[i]=xi=xg_sav[kora];
    ygood[i]=yi=yg_sav[kora];
    idpr[i] =idpr_sav[kora];
    xcog+=egx*xi;
    ycog+=egx*yi;
    egamtot+=egx;
    kora++;

  }
// check the separation between gammas (importantissimo non dimenticare....)
  double x,y;
  int ii;
  for ( int i=0;i<6;i++ )
  {
    x=xgood[i];
    y=ygood[i];
    ii=i+1;
    for ( int j=ii;j<6;j++ )
    {
      if ( ( pow ( ( x-xgood[j] ),2 ) +pow ( ( y-ygood[j] ),2 ) ) <100. ) return -1;
    }
  }
  xcog/=egamtot; ycog/=egamtot;
  cog=sqrt ( xcog*xcog+ycog*ycog );

  return 1;
}

// ----------------------
void  AnaK6g::Save_good_fit()
{
  double xi,yi,xcog=0.0,ycog=0.0,egx;
//  dist_fit=dist_vertex+15.;
  dist_fit=dist_vertex;
  zeta=lkry->Get_Cface().z-dist_fit;
  permut_fit=permut_flag;
  iter_fit=iter_flag;
  egamt_fit=0.0;
  for ( int i=0;i<6;i++ )
  {
    egamt_fit+=eg_fit[i]=egx=ec[i];
    xg_fit[i]=xi=xg[i];
    yg_fit[i]=yi=yg[i];
    idpr_fit[i]=idpr[i];
    xcog+=egx*xi;
    ycog+=egx*yi;
  }
  xcog/=egamt_fit; ycog/=egamt_fit;
  cog_fit=sqrt ( xcog*xcog+ycog*ycog );
}