flyopunta/src/src_par/partp00p.cpp

/***************************************************************************
                          partp00p.cpp  -  description
                             -------------------
    begin                : Tue Apr 30 2002
    copyright            : (C) 2002 by Giuseppe Pierazzini
    email                : peppe@unipi.it
  ***************************************************************************
  *                                                                         *
  *   NA48  simulation program.                                             *
  *                                                                         *
  ***************************************************************************/


#include "flyoh.h"

using namespace std;
// decadimento del K+ e K- in pic+ p00p 
// p00p -> pi0 pi0
// particella virtuale  che attiva questa procedura   p00p

// constructor
Partp00p::Partp00p(){
  mpadre=0.;
}


//---------------
 void  Partp00p::Get_MatrixParm()
{
  // decay  K->pic p00p ; p00p-> pi0 pi0
  // ovvero    mpadre -> mfrate mchd1 mchd2 
  par_a0 =  0.220;        // +/- 0.005 (articolo Cabibbo)
        par_a2 = -0.0444;       // +/- 0.0010   ibid.
  madiff=0.66667*(par_a0 -par_a2);
  par_A0v= 1.;      
  par_Apv= 2.*par_A0v;      
        par_gzero = 0.652;
  par_gplus = -0.2154;
  Ethres   = 4.*mfrateq;
  mfrateq2 = 2.*mfrateq;
  Mtxthr=par_Apv*(1. + par_gplus*(mpadreq-9.*mfrateq)/(12.*mfrateq));
//--------
  normprb = 0;
// prob max finding to speed up the  generation...
   wmtr_max=0.0;
   for (mresq=mminq+0.000001;mresq<mminq+dmq;mresq+=dmq*0.01)
   {Get_MatrixValue();
    if(prb>wmtr_max) wmtr_max=prb;

//    mres=sqrt(mresq);    //per debug
//    double prbf = Fspace(mres, mchd1, mchd2) * Fspace(mpadre, mfrate, mres);
//    double prbb =prbf*prb;
//    double prb0 =Mtx0*Mtx0*prbf;

   }
  Gout<<"\n p00p Wmtr_max  "<< wmtr_max<<" "<<SelRea->Get_Titolo();
  normprb=1;    // switch on = norm prob n done
 }
//---------------


 int  Partp00p::Get_MatrixValue()
{ // see  cabibbo matrx
 double s3=mresq;    // mresq and s0 from  PartFase::Get_PhaseSpace() and Get_PhaseParm()
 Mtx0= par_A0v*(1. + par_gzero*(s3-s0)/mfrateq2);

 if((sdum= (Ethres -s3)/s3 )>0.0)      // select on threshold
    {
      par_j= sqrt(sdum);
      Mtx1 =-madiff*Mtxthr*par_j;    // nota qui  sempre negativo.
      prb = Mtx0*Mtx0 + Mtx1*Mtx1 + 2.*Mtx0*Mtx1 ;  // S3< 4.mpic**2
    }
 else
    { par_j= sqrt(-sdum);     // in verita'  par_j= (-img)*PIGREC*sqrt(-sdum);
      Mtx1 = madiff*Mtxthr*par_j;    // = (+img)*madiff*Mtxthr*par_j   immag positivo  ==>
                                     // prb= Mtx0*Mtx0 + [(+img)*Mtx1]**2 ;
      prb = Mtx0*Mtx0 + Mtx1*Mtx1 ;                // S3> 4.mpic**2
    }                                              

 if(prb>wmtr_max&& normprb==1)
    {
    Gout<<"\n Wmtr off: Ev "
        <<setprecision(4)<<fixed
        <<setw(6)<<evento_.Gen.Event<<" iter="
        <<setw(4)<<itermtr<<" pb="
        <<setw(7)<<prb<<" pbx="
        <<setw(7)<<wmtr_max <<" e0="
        <<setprecision(2)<<fixed
        <<setw(6)<<e0<<" e1="
        <<setw(6)<<e1<<" "
        <<setw(25)<<SelRea->Get_Titolo()<<std::endl;  
        
        
      wmtr_max=prb*1.01; // aggiorno ma informo.... dovrebbe  accadere poche volte!
    }
    
 if((wmtr_max*Pran())<prb) return 1;
 else return 0;
}