[ROOT] Relativistic Breit-Wigner Convoluted with a Gaussian

From: M. Muniruzzaman (munir@ucrhi3.ucr.edu)
Date: Mon Nov 11 2002 - 02:21:37 MET 

Hi, 

I have a few line code of a relativistic Breit-Wigner convoluted with a 
Gaussian. But I am not sure if it is doing the correct normalization. The 
code is attached below. I would appreciate if someone could point me of 
anything I am doing wrong or point me if there is any working code.

Best Regards,

-- munir

// -----------------------------------------------------------------------
//  Convoluted Relativistic Breit Wigner function                        |
//  Author : Basanta K. Nandi                                            |
//           M. Muniruzzaman                                             |
//           University of California Riverside                          |
//------------------------------------------------------------------------

// Convolution integral of Relativistic Breit-Wigner and Gaussian by Trapezoidal rule
Double_t rbwconv(Double_t *x, Double_t *par)
{
  Int_t i, Ndiv;
  Double_t fitval;
  Double_t arg = 0;
  Double_t arg1, arg2, arg3, arg4, arg5;
  Double_t delq, delq0;
  Double_t delq02, delq2;
  Double_t GammaM;
  Double_t mass = 0.493677;//PDG 
  Double_t Arg, Arg1, Func1, Func2;
  Double_t ExpMassRes, Xloww, Xhigh, Delta;
  Double_t Sum;
  Double_t invsq2pi = 0.3989422804014; 
  Double_t binwidth = 0.002;

  Double_t F[201];
  
  // ExpMassRes : User has to put the Experimental Mass Resolution
  ExpMassRes  = 0.0012;
  
  Xloww = -4.*ExpMassRes;
  Xhigh =  4.*ExpMassRes;
  Ndiv  = 200;
  Delta = (Xhigh - Xloww)/(Float_t) Ndiv;

  delq2  = x[0]*x[0]/4. - mass*mass;
  delq02 = par[1]*par[1]/4. - mass*mass;

  if (delq2 <=0.0)delq2 =0.0;
  if (delq02<=0.0)delq02=1.0;

  arg1 = sqrt(delq2/delq02);

  arg2 = arg1*arg1*arg1;
  arg3 = (2.*delq02)/(delq2 + delq02);

  GammaM = par[2]*arg2*arg3;

  for (i=0; i<201; i++)
    {
      Arg1 = (Xloww + (Float_t)i*Delta);
      Arg  = (x[0] - Arg1);
      arg4 = (Arg*Arg - par[1]*par[1]);
      arg5 = par[1]*GammaM;
      arg  = arg4*arg4 + arg5*arg5;

      //Func1 : Gaussian Function
      Func1 = TMath::Exp(-0.5*Arg1*Arg1/(ExpMassRes*ExpMassRes));
      //Func2 : Rel. Breit-Wigner
      Func2 = (Arg*arg5/arg)*(arg5/(TMath::Pi()));//not sure of the factor Pi
      F[i] = Func1*Func2;
    }

  Sum = 0.;
  for (i=1; i <Ndiv; i++)
    {
      Sum += F[i];
    }
  
  fitval = par[0] * (invsq2pi/ExpMassRes) * 0.5*Delta*(F[0] + 2.*Sum + F[200]);
  // Do I need to multiply with the bin width somewhere??

  return fitval;

}

---------------------------------------------------------------------------
                              M. Muniruzzaman
---------------------------------------------------------------------------
Graduate Student         |email:                     |Phone:
Department of Physics    |munir@phyun0.ucr.edu       |(909)-787-3084
University of California |                           |
Riverside, CA 92521      |web:                       |
USA                      |www.phenix.bnl.gov/~munir  |
---------------------------------------------------------------------------

This archive was generated by hypermail 2b29 : Sat Jan 04 2003 - 23:51:17 MET