<commit_before>/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/


#include "AliDecayerExodus.h"
#include <Riostream.h>
#include <TMath.h>
#include <AliLog.h>
#include <TH1.h>
#include <TRandom.h>
#include <TParticle.h>
#include <TDatabasePDG.h>
#include <TPDGCode.h>
#include <TLorentzVector.h>
#include <TClonesArray.h>
#include <TF1.h>


ClassImp(AliDecayerExodus)

//---------------------------------------------------------------------------------------------------
//                                 
// Generate electron-pair mass distributions for Dalitz decays according
// to the Kroll-Wada parametrization: <NAME>, <NAME>: Phys. Rev 98(1955)1355
// and generate electron-pair mass distributions for resonances according
// to the Gounaris-Sakurai parametrization: G.J. Gounaris, <NAME>: Phys.Rev.Lett. 21(1968)244 
//
// For the electromagnetic form factor the parameterization from
// Lepton-G is used: L.G. Landsberg et al.: Phys. Rep. 128(1985)301
//
// <NAME> (R.<EMAIL>) 
// <NAME>  (<EMAIL>)
//
//---------------------------------------------------------------------------------------------------


AliDecayerExodus::AliDecayerExodus():
    AliDecayer(),
    fEPMassPion(0),
    fEPMassEta(0),
    fEPMassEtaPrime(0),
    fEPMassRho(0),
    fEPMassOmega(0),
    fEPMassOmegaDalitz(0),
    fEPMassPhi(0),
    fEPMassPhiDalitz(0),
    fEPMassJPsi(0),
    fPol(new TF1("dsigdcostheta","1.+[0]*x*x",-1.,1.)), /* Polarization Function for resonances */
    fInit(0)

{
// Constructor
}


void AliDecayerExodus::Init()
{
 
// Initialisation
//
   Int_t ibin, nbins;
   Double_t min, maxpion, maxeta, maxomega, maxetaprime, maxphi, binwidth_pion, binwidth_eta, binwidth_omega, binwidth_etaprime, binwidth_phi;
   Double_t pionmass, etamass, omegamass, etaprimemass, phimass, emass, proton_mass, omasspion, omasseta, omassgamma;
   Double_t epsilon_pion, epsilon_eta, epsilon_omega, epsilon_etaprime, epsilon_phi;
   Double_t delta_pion, delta_eta, delta_omega, delta_etaprime, delta_phi;
   Double_t mLL_pion, mLL_eta, mLL_omega, mLL_etaprime, mLL_phi;
   Double_t q_pion, q_eta, q_omega, q_etaprime, q_phi;
   Double_t kwHelp_pion, kwHelp_eta, kwHelp_omega, kwHelp_etaprime, kwHelp_phi;
   Double_t krollWada_pion, krollWada_eta, krollWada_omega, krollWada_etaprime, krollWada_phi;
   Double_t formFactor_pion, formFactor_eta, formFactor_omega, formFactor_etaprime, formFactor_phi;
   Double_t weight_pion, weight_eta, weight_omega_dalitz, weight_etaprime, weight_phi_dalitz;

   Float_t  binwidth;
   Float_t  mass_bin, mass_min, mass_max;
   Double_t vmass_rho, vmass_omega, vmass_phi, vmass_jpsi, vwidth_rho, vwidth_omega, vwidth_phi, vwidth_jpsi;
   Double_t weight_rho, weight_omega, weight_phi, weight_jpsi;

//================================================================================//
//          Create electron pair mass histograms from dalitz decays               //
//================================================================================//

    // Get the particle masses
    // parent
    nbins = 1000;
    mass_min = 0.;
    mass_max = 0.;
    pionmass     = (TDatabasePDG::Instance()->GetParticle(111))->Mass();
    etamass      = (TDatabasePDG::Instance()->GetParticle(221))->Mass();  
    omegamass    = (TDatabasePDG::Instance()->GetParticle(223))->Mass();  
    etaprimemass = (TDatabasePDG::Instance()->GetParticle(331))->Mass();  
    phimass      = (TDatabasePDG::Instance()->GetParticle(333))->Mass();
    // child - electron
    emass = (TDatabasePDG::Instance()->GetParticle(11))->Mass();
    // child - other : third childs from Dalitz decays   
    omasspion  = pionmass;
    omasseta   = etamass;
    omassgamma = 0.;
       
    min         = 2.0 * emass;
    maxpion     = pionmass - omassgamma;
    maxeta      = etamass - omassgamma;
    maxomega    = omegamass - pionmass;
    maxetaprime = etaprimemass - omassgamma;
    maxphi      = phimass - omasseta; 

    binwidth_pion     = (maxpion - min) / (Double_t)nbins;
    binwidth_eta      = (maxeta - min) / (Double_t)nbins;
    binwidth_omega    = (maxomega - min) / (Double_t)nbins;
    binwidth_etaprime = (maxetaprime - min) / (Double_t)nbins;
    binwidth_phi      = (maxphi - min) / (Double_t)nbins;


    epsilon_pion     = (emass / pionmass) * (emass / pionmass);
    epsilon_eta      = (emass / etamass) * (emass / etamass);
    epsilon_omega    = (emass / omegamass) * (emass / omegamass);
    epsilon_etaprime = (emass / etaprimemass) * (emass / etaprimemass);
    epsilon_phi      = (emass / phimass) * (emass / phimass);    


    delta_pion       = (omassgamma / pionmass) * (omassgamma / pionmass);
    delta_eta        = (omassgamma / etamass) * (omassgamma / etamass);
    delta_omega      = (omasspion / omegamass) * (omasspion / omegamass);
    delta_etaprime   = (omassgamma / etaprimemass) * (omassgamma / etaprimemass);
    delta_phi        = (omasseta / phimass) * (omasseta / phimass);    

    // create pair mass histograms for Dalitz decays of pi0, eta, omega, eta' and phi
    if (!fEPMassPion)          {delete fEPMassPion;        fEPMassPion          = new TH1F("fEPMassPion", "Dalitz electron pair from pion", nbins, min, maxpion); }
    if (!fEPMassEta)           {delete fEPMassEta;         fEPMassEta           = new TH1F("fEPMassEta", "Dalitz electron pair from eta", nbins, min, maxeta);}
    if (!fEPMassOmegaDalitz)   {delete fEPMassOmegaDalitz; fEPMassOmegaDalitz   = new TH1F("fEPMassOmegaDalitz", "Dalitz electron pair from omega ", nbins, min, maxomega);}
    if (!fEPMassEtaPrime)      {delete fEPMassEtaPrime;    fEPMassEtaPrime      = new TH1F("fEPMassEtaPrime", "Dalitz electron pair from etaprime", nbins, min, maxetaprime);}
    if (!fEPMassPhiDalitz)     {delete fEPMassPhiDalitz;   fEPMassPhiDalitz     = new TH1F("fEPMassPhiDalitz", "Dalitz electron pair from phi ", nbins, min, maxphi);}


    mLL_pion =  mLL_eta = mLL_omega = mLL_etaprime = mLL_phi = 0.;

    for (ibin = 1; ibin <= nbins; ibin++ )
        {
         mLL_pion     = min + (Double_t)(ibin - 1) * binwidth_pion + binwidth_pion / 2.0;
         mLL_eta      = min + (Double_t)(ibin - 1) * binwidth_eta + binwidth_eta / 2.0;
         mLL_omega    = min + (Double_t)(ibin - 1) * binwidth_omega + binwidth_omega / 2.0;
         mLL_etaprime = min + (Double_t)(ibin - 1) * binwidth_etaprime + binwidth_etaprime / 2.0;
         mLL_phi      = min + (Double_t)(ibin - 1) * binwidth_phi + binwidth_phi / 2.0;


         q_pion        = (mLL_pion / pionmass) * (mLL_pion / pionmass);
         q_eta         = (mLL_eta / etamass) * (mLL_eta / etamass);
         q_omega       = (mLL_omega / omegamass)*(mLL_omega / omegamass);
         q_etaprime    = (mLL_etaprime / etaprimemass) * (mLL_etaprime / etaprimemass);
         q_phi         = (mLL_phi / phimass) * (mLL_phi / phimass);

    if ( q_pion <= 4.0 * epsilon_pion || q_eta <= 4.0 * epsilon_eta || q_omega <= 4.0 * epsilon_omega || q_etaprime <= 4.0 * epsilon_etaprime || q_phi <= 4.0 * epsilon_phi )
       {
         AliFatal("Error in calculating Dalitz mass histogram binning!");
       }
  

    kwHelp_pion     = (1.0 + q_pion /  (1.0 - delta_pion)) * (1.0 + q_pion / (1.0 - delta_pion))
                                     - 4.0 * q_pion / ((1.0 - delta_pion) * (1.0 - delta_pion));

    kwHelp_eta      = (1.0 + q_eta /  (1.0 - delta_eta)) * (1.0 + q_eta / (1.0 - delta_eta))
                                    - 4.0 * q_eta / ((1.0 - delta_eta) * (1.0 - delta_eta));

    kwHelp_omega    = (1.0 + q_omega /  (1.0 - delta_omega)) * (1.0 + q_omega / (1.0 - delta_omega))
                                      - 4.0 * q_omega / ((1.0 - delta_omega) * (1.0 - delta_omega));

    kwHelp_etaprime = (1.0 + q_etaprime /  (1.0 - delta_etaprime)) * (1.0 + q_etaprime / (1.0 - delta_etaprime))
                                         - 4.0 * q_etaprime / ((1.0 - delta_etaprime) * (1.0 - delta_etaprime));

    kwHelp_phi      = (1.0 + q_phi /  (1.0 - delta_phi)) * (1.0 + q_phi / (1.0 - delta_phi))
                                    - 4.0 * q_phi / ((1.0 - delta_phi) * (1.0 - delta_phi));




    if ( kwHelp_pion <= 0.0 || kwHelp_eta <= 0.0 || kwHelp_omega <= 0.0 || kwHelp_etaprime <= 0.0 || kwHelp_phi <= 0.0 )
       {
         AliFatal("Error in calculating Dalitz mass histogram binning!");
    
       }


 // Invariant mass distributions of electron pairs from Dalitz decays
 // using Kroll-Wada function   

      krollWada_pion     = (2.0 / mLL_pion) * TMath::Exp(1.5 * TMath::Log(kwHelp_pion))
                                            * TMath::Sqrt(1.0 - 4.0 * epsilon_pion / q_pion)
                                            * (1.0 + 2.0 * epsilon_pion / q_pion);
    
     
      krollWada_eta      = (2.0 / mLL_eta) * TMath::Exp(1.5 * TMath::Log(kwHelp_eta))
                                           * TMath::Sqrt(1.0 - 4.0 * epsilon_eta / q_eta)
                                           * (1.0 + 2.0 * epsilon_eta / q_eta);
   
   
      krollWada_omega    = (2.0 / mLL_omega) * TMath::Exp(1.5 * TMath::Log(kwHelp_omega))
                                             * TMath::Sqrt(1.0 - 4.0 * epsilon_omega / q_omega)
                                             * (1.0 + 2.0 * epsilon_omega / q_omega);
   
   
      krollWada_etaprime = (2.0 / mLL_etaprime) * TMath::Exp(1.5 * TMath::Log(kwHelp_etaprime))
                                                * TMath::Sqrt(1.0 - 4.0 * epsilon_etaprime / q_etaprime)
                                                * (1.0 + 2.0 * epsilon_etaprime / q_etaprime);
   
      krollWada_phi      = (2.0 / mLL_phi) * TMath::Exp(1.5 * TMath::Log(kwHelp_phi))
                                           * TMath::Sqrt(1.0 - 4.0 * epsilon_phi / q_phi)
                                           * (1.0 + 2.0 * epsilon_phi / q_phi);   



    // Form factors from Lepton-G  
    formFactor_pion     = 1.0/(1.0-5.5*mLL_pion*mLL_pion);
    formFactor_eta      = 1.0/(1.0-1.9*mLL_eta*mLL_eta);
    formFactor_omega    = (TMath::Power(TMath::Power(0.6519,2),2))
                          /(TMath::Power(TMath::Power(0.6519,2)-TMath::Power(mLL_omega, 2), 2)
                          + TMath::Power(0.04198, 2)*TMath::Power(0.6519, 2));
    formFactor_etaprime = (TMath::Power(TMath::Power(0.764,2),2))
                          /(TMath::Power(TMath::Power(0.764,2)-TMath::Power(mLL_etaprime, 2), 2)
                          + TMath::Power(0.1020, 2)*TMath::Power(0.764, 2));
    formFactor_phi      = 1.0; 




    weight_pion         = krollWada_pion * formFactor_pion * formFactor_pion;
    weight_eta          = krollWada_eta * formFactor_eta * formFactor_eta;
    weight_omega_dalitz = krollWada_omega * formFactor_omega;
    weight_etaprime     = krollWada_etaprime * formFactor_etaprime;
    weight_phi_dalitz   = krollWada_phi * formFactor_phi * formFactor_phi;


    // Fill histograms of electron pair masses from dalitz decays
    fEPMassPion       ->AddBinContent(ibin, weight_pion);
    fEPMassEta        ->AddBinContent(ibin, weight_eta);
    fEPMassOmegaDalitz->AddBinContent(ibin, weight_omega_dalitz);
    fEPMassEtaPrime   ->AddBinContent(ibin, weight_etaprime);
    fEPMassPhiDalitz  ->AddBinContent(ibin, weight_phi_dalitz);
    }


   

//===================================================================================//
//         Create electron pair mass histograms from resonance decays                //
//===================================================================================//

   Double_t pimass = 0.13956995;

   // Get the particle masses
   // parent
   vmass_rho   = (TDatabasePDG::Instance()->GetParticle(113))->Mass();  
   vmass_omega = (TDatabasePDG::Instance()->GetParticle(223))->Mass();  
   vmass_phi   = (TDatabasePDG::Instance()->GetParticle(333))->Mass();  
   vmass_jpsi  = (TDatabasePDG::Instance()->GetParticle(443))->Mass();
   // Get the particle widths
   // parent
   vwidth_rho   = (TDatabasePDG::Instance()->GetParticle(113))->Width();   
   vwidth_omega = (TDatabasePDG::Instance()->GetParticle(223))->Width();  
   vwidth_phi   = (TDatabasePDG::Instance()->GetParticle(333))->Width();  
   vwidth_jpsi  = (TDatabasePDG::Instance()->GetParticle(443))->Width();


       if ( mass_min == 0. && mass_max == 0. )
          {
           mass_min  = 2.*pimass;
           mass_max  = 5.;
          }

     binwidth  = (mass_max-mass_min)/(Double_t)nbins;

     // create pair mass histograms for resonances of rho, omega, phi and jpsi
     if (!fEPMassRho)   {delete fEPMassRho;   fEPMassRho    = new TH1F("fEPMassRho","mass rho",nbins,mass_min,mass_max);}
     if (!fEPMassOmega) {delete fEPMassOmega; fEPMassOmega  = new TH1F("fEPMassOmega","mass omega",nbins,mass_min,mass_max);}
     if (!fEPMassPhi)   {delete fEPMassPhi;   fEPMassPhi    = new TH1F("fEPMassPhi","mass phi",nbins,mass_min,mass_max);}
     if (!fEPMassJPsi)  {delete fEPMassJPsi;  fEPMassJPsi   = new TH1F("fEPMassJPsi","mass jpsi",nbins,mass_min,mass_max);}

     for (ibin=1; ibin<=nbins; ibin++ )
     {
     mass_bin = mass_min+(Double_t)(ibin-1)*binwidth+binwidth/2.0;

     weight_rho     = (Float_t)GounarisSakurai(mass_bin,vmass_rho,vwidth_rho,emass);
     weight_omega   = (Float_t)GounarisSakurai(mass_bin,vmass_omega,vwidth_omega,emass);
     weight_phi     = (Float_t)GounarisSakurai(mass_bin,vmass_phi,vwidth_phi,emass); 
     weight_jpsi    = (Float_t)Lorentz(mass_bin,vmass_jpsi,vwidth_jpsi);

     // Fill histograms of electron pair masses from resonance decays
     fEPMassRho  ->AddBinContent(ibin,weight_rho);
     fEPMassOmega->AddBinContent(ibin,weight_omega);
     fEPMassPhi  ->AddBinContent(ibin,weight_phi);
     fEPMassJPsi ->AddBinContent(ibin,weight_jpsi);
    }  

}

Double_t AliDecayerExodus::GounarisSakurai(Float_t mass, Double_t vmass, Double_t vwidth, Double_t emass)
{
// Invariant mass distributions of electron pairs from resonance decays
// of rho, omega and phi
// using Gounaris-Sakurai function

  Double_t corr = 0.;
  Double_t epsilon = 0.;
  Double_t weight = 0.;

  Double_t pimass = 0.13956995;
 
  if(mass>pimass){
  corr = vwidth*(vmass/mass)*exp(1.5*log((mass*mass/4.0-pimass*pimass)
         /(vmass*vmass/4.0-pimass*pimass)));
  }

  epsilon = (emass/mass)*(emass/mass);
       
  if ( 1.0-4.0*epsilon>=0.0 )
  {
   weight = sqrt(1.0-4.0*epsilon)*(1.0+2.0*epsilon)/
                ((vmass*vmass-mass*mass)*(vmass*vmass-mass*mass)+
                (vmass*corr)*(vmass*corr));
  }
  return weight;  
}


Double_t AliDecayerExodus::Lorentz(Float_t mass, Double_t vmass, Double_t vwidth)
{
// Invariant mass distributions of electron pairs from resonance decay
// of jpsi (and it can also be used for other particles except rho, omega and phi) 
// using Lorentz function

  Double_t weight;
  
  weight = (vwidth*vwidth/4.0)/(vwidth*vwidth/4.0+(vmass-mass)*(vmass-mass));

  return weight;

}

void AliDecayerExodus::Decay(Int_t idpart, TLorentzVector* pparent)
{
 
    if (!fInit) {
        Init();
        fInit=1;  
    }

   //local variables for dalitz/2-body decay:
   Double_t pmass, epmass, realp_mass, e1, p1, e3, p3;
   Double_t wp_res, mp_res, md_res, epmass_res, Ed_res, pd_res;
   Double_t PolPar;
   TLorentzVector fProducts_res[2], fProducts_dalitz[3];
   Int_t idRho=113;
   Int_t idOmega=223;<commit_msg>Adding protection in case width used for Lorentz function is 0 -> set it to 1e-6, i.e. 1 keV<commit_after>/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/


#include "AliDecayerExodus.h"
#include <Riostream.h>
#include <TMath.h>
#include <AliLog.h>
#include <TH1.h>
#include <TRandom.h>
#include <TParticle.h>
#include <TDatabasePDG.h>
#include <TPDGCode.h>
#include <TLorentzVector.h>
#include <TClonesArray.h>
#include <TF1.h>


ClassImp(AliDecayerExodus)

//---------------------------------------------------------------------------------------------------
//                                 
// Generate electron-pair mass distributions for Dalitz decays according
// to the Kroll-Wada parametrization: <NAME>, <NAME>: Phys. Rev 98(1955)1355
// and generate electron-pair mass distributions for resonances according
// to the Gounaris-Sakurai parametrization: <NAME>, <NAME>: Phys.Rev.Lett. 21(1968)244 
//
// For the electromagnetic form factor the parameterization from
// Lepton-G is used: L.G. Landsberg et al.: Phys. Rep. 128(1985)301
//
// <NAME> (<EMAIL>) 
// <NAME>  (<EMAIL>)
//
//---------------------------------------------------------------------------------------------------


AliDecayerExodus::AliDecayerExodus():
    AliDecayer(),
    fEPMassPion(0),
    fEPMassEta(0),
    fEPMassEtaPrime(0),
    fEPMassRho(0),
    fEPMassOmega(0),
    fEPMassOmegaDalitz(0),
    fEPMassPhi(0),
    fEPMassPhiDalitz(0),
    fEPMassJPsi(0),
    fPol(new TF1("dsigdcostheta","1.+[0]*x*x",-1.,1.)), /* Polarization Function for resonances */
    fInit(0)

{
// Constructor
}


void AliDecayerExodus::Init()
{
 
// Initialisation
//
   Int_t ibin, nbins;
   Double_t min, maxpion, maxeta, maxomega, maxetaprime, maxphi, binwidth_pion, binwidth_eta, binwidth_omega, binwidth_etaprime, binwidth_phi;
   Double_t pionmass, etamass, omegamass, etaprimemass, phimass, emass, proton_mass, omasspion, omasseta, omassgamma;
   Double_t epsilon_pion, epsilon_eta, epsilon_omega, epsilon_etaprime, epsilon_phi;
   Double_t delta_pion, delta_eta, delta_omega, delta_etaprime, delta_phi;
   Double_t mLL_pion, mLL_eta, mLL_omega, mLL_etaprime, mLL_phi;
   Double_t q_pion, q_eta, q_omega, q_etaprime, q_phi;
   Double_t kwHelp_pion, kwHelp_eta, kwHelp_omega, kwHelp_etaprime, kwHelp_phi;
   Double_t krollWada_pion, krollWada_eta, krollWada_omega, krollWada_etaprime, krollWada_phi;
   Double_t formFactor_pion, formFactor_eta, formFactor_omega, formFactor_etaprime, formFactor_phi;
   Double_t weight_pion, weight_eta, weight_omega_dalitz, weight_etaprime, weight_phi_dalitz;

   Float_t  binwidth;
   Float_t  mass_bin, mass_min, mass_max;
   Double_t vmass_rho, vmass_omega, vmass_phi, vmass_jpsi, vwidth_rho, vwidth_omega, vwidth_phi, vwidth_jpsi;
   Double_t weight_rho, weight_omega, weight_phi, weight_jpsi;

//================================================================================//
//          Create electron pair mass histograms from dalitz decays               //
//================================================================================//

    // Get the particle masses
    // parent
    nbins = 1000;
    mass_min = 0.;
    mass_max = 0.;
    pionmass     = (TDatabasePDG::Instance()->GetParticle(111))->Mass();
    etamass      = (TDatabasePDG::Instance()->GetParticle(221))->Mass();  
    omegamass    = (TDatabasePDG::Instance()->GetParticle(223))->Mass();  
    etaprimemass = (TDatabasePDG::Instance()->GetParticle(331))->Mass();  
    phimass      = (TDatabasePDG::Instance()->GetParticle(333))->Mass();
    // child - electron
    emass = (TDatabasePDG::Instance()->GetParticle(11))->Mass();
    // child - other : third childs from Dalitz decays   
    omasspion  = pionmass;
    omasseta   = etamass;
    omassgamma = 0.;
       
    min         = 2.0 * emass;
    maxpion     = pionmass - omassgamma;
    maxeta      = etamass - omassgamma;
    maxomega    = omegamass - pionmass;
    maxetaprime = etaprimemass - omassgamma;
    maxphi      = phimass - omasseta; 

    binwidth_pion     = (maxpion - min) / (Double_t)nbins;
    binwidth_eta      = (maxeta - min) / (Double_t)nbins;
    binwidth_omega    = (maxomega - min) / (Double_t)nbins;
    binwidth_etaprime = (maxetaprime - min) / (Double_t)nbins;
    binwidth_phi      = (maxphi - min) / (Double_t)nbins;


    epsilon_pion     = (emass / pionmass) * (emass / pionmass);
    epsilon_eta      = (emass / etamass) * (emass / etamass);
    epsilon_omega    = (emass / omegamass) * (emass / omegamass);
    epsilon_etaprime = (emass / etaprimemass) * (emass / etaprimemass);
    epsilon_phi      = (emass / phimass) * (emass / phimass);    


    delta_pion       = (omassgamma / pionmass) * (omassgamma / pionmass);
    delta_eta        = (omassgamma / etamass) * (omassgamma / etamass);
    delta_omega      = (omasspion / omegamass) * (omasspion / omegamass);
    delta_etaprime   = (omassgamma / etaprimemass) * (omassgamma / etaprimemass);
    delta_phi        = (omasseta / phimass) * (omasseta / phimass);    

    // create pair mass histograms for Dalitz decays of pi0, eta, omega, eta' and phi
    if (!fEPMassPion)          {delete fEPMassPion;        fEPMassPion          = new TH1F("fEPMassPion", "Dalitz electron pair from pion", nbins, min, maxpion); }
    if (!fEPMassEta)           {delete fEPMassEta;         fEPMassEta           = new TH1F("fEPMassEta", "Dalitz electron pair from eta", nbins, min, maxeta);}
    if (!fEPMassOmegaDalitz)   {delete fEPMassOmegaDalitz; fEPMassOmegaDalitz   = new TH1F("fEPMassOmegaDalitz", "Dalitz electron pair from omega ", nbins, min, maxomega);}
    if (!fEPMassEtaPrime)      {delete fEPMassEtaPrime;    fEPMassEtaPrime      = new TH1F("fEPMassEtaPrime", "Dalitz electron pair from etaprime", nbins, min, maxetaprime);}
    if (!fEPMassPhiDalitz)     {delete fEPMassPhiDalitz;   fEPMassPhiDalitz     = new TH1F("fEPMassPhiDalitz", "Dalitz electron pair from phi ", nbins, min, maxphi);}


    mLL_pion =  mLL_eta = mLL_omega = mLL_etaprime = mLL_phi = 0.;

    for (ibin = 1; ibin <= nbins; ibin++ )
        {
         mLL_pion     = min + (Double_t)(ibin - 1) * binwidth_pion + binwidth_pion / 2.0;
         mLL_eta      = min + (Double_t)(ibin - 1) * binwidth_eta + binwidth_eta / 2.0;
         mLL_omega    = min + (Double_t)(ibin - 1) * binwidth_omega + binwidth_omega / 2.0;
         mLL_etaprime = min + (Double_t)(ibin - 1) * binwidth_etaprime + binwidth_etaprime / 2.0;
         mLL_phi      = min + (Double_t)(ibin - 1) * binwidth_phi + binwidth_phi / 2.0;


         q_pion        = (mLL_pion / pionmass) * (mLL_pion / pionmass);
         q_eta         = (mLL_eta / etamass) * (mLL_eta / etamass);
         q_omega       = (mLL_omega / omegamass)*(mLL_omega / omegamass);
         q_etaprime    = (mLL_etaprime / etaprimemass) * (mLL_etaprime / etaprimemass);
         q_phi         = (mLL_phi / phimass) * (mLL_phi / phimass);

    if ( q_pion <= 4.0 * epsilon_pion || q_eta <= 4.0 * epsilon_eta || q_omega <= 4.0 * epsilon_omega || q_etaprime <= 4.0 * epsilon_etaprime || q_phi <= 4.0 * epsilon_phi )
       {
         AliFatal("Error in calculating Dalitz mass histogram binning!");
       }
  

    kwHelp_pion     = (1.0 + q_pion /  (1.0 - delta_pion)) * (1.0 + q_pion / (1.0 - delta_pion))
                                     - 4.0 * q_pion / ((1.0 - delta_pion) * (1.0 - delta_pion));

    kwHelp_eta      = (1.0 + q_eta /  (1.0 - delta_eta)) * (1.0 + q_eta / (1.0 - delta_eta))
                                    - 4.0 * q_eta / ((1.0 - delta_eta) * (1.0 - delta_eta));

    kwHelp_omega    = (1.0 + q_omega /  (1.0 - delta_omega)) * (1.0 + q_omega / (1.0 - delta_omega))
                                      - 4.0 * q_omega / ((1.0 - delta_omega) * (1.0 - delta_omega));

    kwHelp_etaprime = (1.0 + q_etaprime /  (1.0 - delta_etaprime)) * (1.0 + q_etaprime / (1.0 - delta_etaprime))
                                         - 4.0 * q_etaprime / ((1.0 - delta_etaprime) * (1.0 - delta_etaprime));

    kwHelp_phi      = (1.0 + q_phi /  (1.0 - delta_phi)) * (1.0 + q_phi / (1.0 - delta_phi))
                                    - 4.0 * q_phi / ((1.0 - delta_phi) * (1.0 - delta_phi));




    if ( kwHelp_pion <= 0.0 || kwHelp_eta <= 0.0 || kwHelp_omega <= 0.0 || kwHelp_etaprime <= 0.0 || kwHelp_phi <= 0.0 )
       {
         AliFatal("Error in calculating Dalitz mass histogram binning!");
    
       }


 // Invariant mass distributions of electron pairs from Dalitz decays
 // using Kroll-Wada function   

      krollWada_pion     = (2.0 / mLL_pion) * TMath::Exp(1.5 * TMath::Log(kwHelp_pion))
                                            * TMath::Sqrt(1.0 - 4.0 * epsilon_pion / q_pion)
                                            * (1.0 + 2.0 * epsilon_pion / q_pion);
    
     
      krollWada_eta      = (2.0 / mLL_eta) * TMath::Exp(1.5 * TMath::Log(kwHelp_eta))
                                           * TMath::Sqrt(1.0 - 4.0 * epsilon_eta / q_eta)
                                           * (1.0 + 2.0 * epsilon_eta / q_eta);
   
   
      krollWada_omega    = (2.0 / mLL_omega) * TMath::Exp(1.5 * TMath::Log(kwHelp_omega))
                                             * TMath::Sqrt(1.0 - 4.0 * epsilon_omega / q_omega)
                                             * (1.0 + 2.0 * epsilon_omega / q_omega);
   
   
      krollWada_etaprime = (2.0 / mLL_etaprime) * TMath::Exp(1.5 * TMath::Log(kwHelp_etaprime))
                                                * TMath::Sqrt(1.0 - 4.0 * epsilon_etaprime / q_etaprime)
                                                * (1.0 + 2.0 * epsilon_etaprime / q_etaprime);
   
      krollWada_phi      = (2.0 / mLL_phi) * TMath::Exp(1.5 * TMath::Log(kwHelp_phi))
                                           * TMath::Sqrt(1.0 - 4.0 * epsilon_phi / q_phi)
                                           * (1.0 + 2.0 * epsilon_phi / q_phi);   



    // Form factors from Lepton-G  
    formFactor_pion     = 1.0/(1.0-5.5*mLL_pion*mLL_pion);
    formFactor_eta      = 1.0/(1.0-1.9*mLL_eta*mLL_eta);
    formFactor_omega    = (TMath::Power(TMath::Power(0.6519,2),2))
                          /(TMath::Power(TMath::Power(0.6519,2)-TMath::Power(mLL_omega, 2), 2)
                          + TMath::Power(0.04198, 2)*TMath::Power(0.6519, 2));
    formFactor_etaprime = (TMath::Power(TMath::Power(0.764,2),2))
                          /(TMath::Power(TMath::Power(0.764,2)-TMath::Power(mLL_etaprime, 2), 2)
                          + TMath::Power(0.1020, 2)*TMath::Power(0.764, 2));
    formFactor_phi      = 1.0; 




    weight_pion         = krollWada_pion * formFactor_pion * formFactor_pion;
    weight_eta          = krollWada_eta * formFactor_eta * formFactor_eta;
    weight_omega_dalitz = krollWada_omega * formFactor_omega;
    weight_etaprime     = krollWada_etaprime * formFactor_etaprime;
    weight_phi_dalitz   = krollWada_phi * formFactor_phi * formFactor_phi;


    // Fill histograms of electron pair masses from dalitz decays
    fEPMassPion       ->AddBinContent(ibin, weight_pion);
    fEPMassEta        ->AddBinContent(ibin, weight_eta);
    fEPMassOmegaDalitz->AddBinContent(ibin, weight_omega_dalitz);
    fEPMassEtaPrime   ->AddBinContent(ibin, weight_etaprime);
    fEPMassPhiDalitz  ->AddBinContent(ibin, weight_phi_dalitz);
    }


   

//===================================================================================//
//         Create electron pair mass histograms from resonance decays                //
//===================================================================================//

   Double_t pimass = 0.13956995;

   // Get the particle masses
   // parent
   vmass_rho   = (TDatabasePDG::Instance()->GetParticle(113))->Mass();  
   vmass_omega = (TDatabasePDG::Instance()->GetParticle(223))->Mass();  
   vmass_phi   = (TDatabasePDG::Instance()->GetParticle(333))->Mass();  
   vmass_jpsi  = (TDatabasePDG::Instance()->GetParticle(443))->Mass();
   // Get the particle widths
   // parent
   vwidth_rho   = (TDatabasePDG::Instance()->GetParticle(113))->Width();   
   vwidth_omega = (TDatabasePDG::Instance()->GetParticle(223))->Width();  
   vwidth_phi   = (TDatabasePDG::Instance()->GetParticle(333))->Width();  
   vwidth_jpsi  = (TDatabasePDG::Instance()->GetParticle(443))->Width();


       if ( mass_min == 0. && mass_max == 0. )
          {
           mass_min  = 2.*pimass;
           mass_max  = 5.;
          }

     binwidth  = (mass_max-mass_min)/(Double_t)nbins;

     // create pair mass histograms for resonances of rho, omega, phi and jpsi
     if (!fEPMassRho)   {delete fEPMassRho;   fEPMassRho    = new TH1F("fEPMassRho","mass rho",nbins,mass_min,mass_max);}
     if (!fEPMassOmega) {delete fEPMassOmega; fEPMassOmega  = new TH1F("fEPMassOmega","mass omega",nbins,mass_min,mass_max);}
     if (!fEPMassPhi)   {delete fEPMassPhi;   fEPMassPhi    = new TH1F("fEPMassPhi","mass phi",nbins,mass_min,mass_max);}
     if (!fEPMassJPsi)  {delete fEPMassJPsi;  fEPMassJPsi   = new TH1F("fEPMassJPsi","mass jpsi",nbins,mass_min,mass_max);}

     for (ibin=1; ibin<=nbins; ibin++ )
     {
     mass_bin = mass_min+(Double_t)(ibin-1)*binwidth+binwidth/2.0;

     weight_rho     = (Float_t)GounarisSakurai(mass_bin,vmass_rho,vwidth_rho,emass);
     weight_omega   = (Float_t)GounarisSakurai(mass_bin,vmass_omega,vwidth_omega,emass);
     weight_phi     = (Float_t)GounarisSakurai(mass_bin,vmass_phi,vwidth_phi,emass); 
     weight_jpsi    = (Float_t)Lorentz(mass_bin,vmass_jpsi,vwidth_jpsi);

     // Fill histograms of electron pair masses from resonance decays
     fEPMassRho  ->AddBinContent(ibin,weight_rho);
     fEPMassOmega->AddBinContent(ibin,weight_omega);
     fEPMassPhi  ->AddBinContent(ibin,weight_phi);
     fEPMassJPsi ->AddBinContent(ibin,weight_jpsi);
    }  

}

Double_t AliDecayerExodus::GounarisSakurai(Float_t mass, Double_t vmass, Double_t vwidth, Double_t emass)
{
// Invariant mass distributions of electron pairs from resonance decays
// of rho, omega and phi
// using Gounaris-Sakurai function

  Double_t corr = 0.;
  Double_t epsilon = 0.;
  Double_t weight = 0.;

  Double_t pimass = 0.13956995;
 
  if(mass>pimass){
  corr = vwidth*(vmass/mass)*exp(1.5*log((mass*mass/4.0-pimass*pimass)
         /(vmass*vmass/4.0-pimass*pimass)));
  }

  epsilon = (emass/mass)*(emass/mass);
       
  if ( 1.0-4.0*epsilon>=0.0 )
  {
   weight = sqrt(1.0-4.0*epsilon)*(1.0+2.0*epsilon)/
                ((vmass*vmass-mass*mass)*(vmass*vmass-mass*mass)+
                (vmass*corr)*(vmass*corr));
  }
  return weight;  
}


Double_t AliDecayerExodus::Lorentz(Float_t mass, Double_t vmass, Double_t vwidth)
{
// Invariant mass distributions of electron pairs from resonance decay
// of jpsi (and it can also be used for other particles except rho, omega and phi) 
// using Lorentz function

  Double_t weight;

  if( vwidth < 1e-6 ){
    AliWarning(Form("Width = %f (Mass = %f) < 1e-6, set to 1e-6",vwidth,vmass));
    vwidth = 1e-6;
  }
  
  weight = (vwidth*vwidth/4.0)/(vwidth*vwidth/4.0+(vmass-mass)*(vmass-mass));

  return weight;

}

void AliDecayerExodus::Decay(Int_t idpart, TLorentzVector* pparent)
{
 
    if (!fInit) {
        Init();
        fInit=1;  
    }

   //local variables for dalitz/2-body decay:
   Double_t pmass, epmass, realp_mass, e1, p1, e3, p3;
   Double_t wp_res, mp_res, md_res, epmass_res, Ed_res, pd_res;
   Double_t PolPar;
   TLorentzVector fProducts_res[2], fProducts_dalitz[3];
   Int_t idRho=113;
   Int_t idOmega=223;
{'fraction_non_alphanumeric': 0.16330460297084173, 'fraction_numerical': 0.025673940949935817, 'mean_word_length': 2.468935538592027, 'pattern_counts': {'":': 0, '<': 66, '<?xml version=': 0, '>': 119, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 66, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3364324', 'n_tokens_mistral': 12303, 'n_tokens_neox': 11560, 'n_words': 2222}