/////////////////////////////////////////////////////////////////////
// 
// Decay Sigma11385_0 -> Lambda e+ e- 
//
// References:
// [L1] Hadronic three-body decays of light vector mesons.
//      S. Leupold, (Frankfurt U.) , M.F.M. Lutz, (Darmstadt, GSI)
//      Jul 2008. 7pp.
//      Published in Eur.Phys.J.A39:205-212,2009.
//      e-Print: arXiv:0807.4686 [hep-ph] 
// 
// Authors:  I. Froehlich, T. Scheib and S. Leupold
/////////////////////////////////////////////////////////////////////


using namespace std;
#include <sstream>
#include <iostream>
#include <iomanip>

#include "PSigmatoLee.h"


ClassImp(PSigmatoLee)


PSigmatoLee::PSigmatoLee() {

    primary = NULL;
    parent = NULL;
    
 
} ;

PSigmatoLee::PSigmatoLee(const Char_t *id, const Char_t *de) :
    PDistribution(id, de) {

    primary = NULL;
    parent = NULL;
    max = -1;
} ;

PDistribution* PSigmatoLee::Clone(const char*delme) const {
    return new PSigmatoLee((const PSigmatoLee &)* this);
};

Bool_t PSigmatoLee::Init(void) {

    
    //looking for primary. This is mandatory
    primary = GetParticle("primary");
    if (!primary) {
	Warning("Init","Primary not found");
	return kFALSE;
    }

    //now get the parent
    for (int i=0; i<position; i++) {
	if (particle_flag[i] == PARTICLE_LIST_PARENT)
	    parent=particle[i];
    }
    if (!parent) {
	Warning("Init","Parent not found");
	return kFALSE;
    }

//    int side = 0;

    for (int i=0; i<position; i++) {
	if ((particle_flag[i] & PARTICLE_LIST_DAUGHTER) 
	    && (particle[i] != primary)) {
	    lambda_ = particle[i];
	}
    }


    return kTRUE;    
};

Bool_t PSigmatoLee::Prepare(void) {
    return kTRUE;
};

Bool_t PSigmatoLee::Finalize(void) {
    return kTRUE;
};

Bool_t PSigmatoLee::CheckAbort(void) {
    return kFALSE;
};

Bool_t PSigmatoLee::IsValid(void) {

    TLorentzVector dl = (*(TLorentzVector *) primary);
    TLorentzVector sigma = (*(TLorentzVector *) parent);
    TLorentzVector lambda = (*(TLorentzVector *) lambda_);

    double m2ll = dl.M2(); 
    double m2sigma = sigma.M2(); 
    double m2lambda = lambda.M2(); 
    double factor = PSigmatoLee::diffgam(m2ll,m2sigma,m2lambda);
    

    if (factor>max) {
	Warning("IsValid","Dalitz factor %f > max %f",factor,max);
	if (max < 0) {
	    Warning("IsValid","Dalitz factor max not set");
	}
	max = factor*1.02;
    } 

    if (factor<0  ) {
	return kFALSE;
    } 

    if (((factor/max))>PUtils::sampleFlat()) {   
	return kTRUE; 
    }

    return kFALSE;
};


double PSigmatoLee::diffgam(double m2ll,double m2sigma,double m2lambda) {
    //Written by S. Leupold
    double pi = TMath::Pi();

    // physical parameters m_V, h_P, h_A, f, pion mass m_pi, omega mass m_omega
    double msigma=sqrt(m2sigma);
    double mlambda=sqrt(m2lambda);
    // phase space, returned variable
    

    double alpha=1./137;
    double me=0.000511;
    double q=sqrt(pow(m2ll-m2lambda+m2ll,2)/(4*m2sigma)-m2ll);
    double G_M=2.7;
 
    double QED_pl=2*alpha/(3*pi*sqrt(m2ll))*sqrt(1-4*me*me/m2ll)*(1+me*me/m2ll);
    double M2_AtoBvg=G_M*G_M*pow((msigma+mlambda),2)*(pow((msigma-mlambda),2)-m2ll)/(4*m2lambda)/(pow((msigma+mlambda),2)-m2ll)*(7.0*m2sigma*m2sigma+14.0*m2sigma*m2ll+3.0*m2ll*m2ll+8.0*m2sigma*msigma*mlambda+2*m2sigma*m2lambda+6*m2ll*m2lambda+3.0*m2lambda*m2lambda);

    double dg = QED_pl*q/(8.0*pi*m2sigma)*M2_AtoBvg;
    
    return dg;
};