flyopunta/src/src_par/partdltz.cpp

/***************************************************************************
                          partdltz.cpp  -  description
                             -------------------
    begin                : Sun May 9 2004
    copyright            : (C) 2004 by gmp
    email                : peppe@newpeppe
  ***************************************************************************
  *                                                                         *
  *   NA48  simulation program.                                             *
  *                                                                         *
  ***************************************************************************/

#define EPSR     1.0E-200

#include "flyoh.h"
using namespace GmpConst;

//             D a l i t z
//--------------------------------------------------------------
// Get the Dalitz pair mass  and decay angles

void PartDltz::
 DecParm()
 {

// Attenzione questa routine vale solo per la
// produzione di dalitz pair nel decadimento del mesone pi0 in gam e e
// si suppone una distribuzione che va come df/dmq=(1-mq/mpiq)**3 *1./mq;
// limiti tra 2 volte la massa dell'elettrone e quella del pi0.
// mass electron= 0.00051 , mass pi0 = 0.13497 (**2=0.0182169009)
// vedi anche  "Selected problems in theoretical physics"

  double wtot=0.0,x=0.0,xq=0.0,aux=0.0;
  int i,j=0;
//==============================================================
// inizializza la forma integrale della distribuzione di massa normalizzata.
  if(dlznorm==0.0)
   {
     xmindlz=0.0075572349;    // = 2 me/mpio
     xmaxdlz=1.0;
     stepdlz=(xmaxdlz-xmindlz)/100000.;
     probdlz[0]=0.0;

     for (i=1;i<100000;i++)
      { x=xmindlz+stepdlz*i;
        xq=x*x;
        aux=1.-xq;
        wtot+= (aux*aux*aux)/x;   // df/dm
        if(i%1000 ==0)
        {j++;
         probdlz[j]=wtot;

        }
      }
// normalizzo....
   probdlz[100]=1.0;
   probdlz[101]=10.0;
   stepdlz=(xmaxdlz-xmindlz)/100.;
   dlznorm=1./wtot;
   for (i=1;i<101;i++) probdlz[i]*= dlznorm;
  }
  
//=============================================================
//------------------------------------------
// qui comincia normalmente... get massa and angles

  double drm;
  double R=Pran();
  i=1;
  while(1)
   {
    if( probdlz[i]>R)break;
    i++;
   }
// interpolazione ..... piu' peso alle basse masse
   drm=(R-probdlz[i-1])/(probdlz[i]-probdlz[i-1]);
   double mm=xmindlz+(i-1.0+drm*drm)*stepdlz;

   P[3] = P.m=massa=massdlz=mm*0.13497;



   SelRea->Dalitz_done++;


// calcolo la direzione di decadimento dei figli.

//   double cost = 1.0 -  2.0 * Pran();
/**********/
   double rdm=2.*Pran()-1.0;    // ==cos(tet)**2
   double cost = exp(log( fabs(rdm)+EPSR)/3.);
   if(rdm<0.0)cost=-cost;
/*************
   double cost = 1.0;
   while (1)
   {
   cost = 1.0 -  2.0 * Pran();
   if(Pran()<(cost*cost)) break;

   }
*********
    double cost = 1.0;
     if(Pran()<0.5)cost=-1;
*******/

   double   fi = DuePiGreco * Pran();
   dw->cstar=cost;     dw->fistar=fi;
   dw->rg->cstar=-cost; dw->rg->fistar=fi+PiGreco;
   return;
 }