//************************************************************************//
//									  //
//  Copyright 2021 Bertram Kopf (bertram@ep1.rub.de)			  //
//          	   - Ruhr-Universität Bochum 				  //
//									  //
//  This file is part of Pawian.					  //
//									  //
//  Pawian is free software: you can redistribute it and/or modify	  //
//  it under the terms of the GNU General Public License as published by  //
//  the Free Software Foundation, either version 3 of the License, or 	  //
//  (at your option) any later version.	 	      	  	   	  //
//									  //
//  Pawian is distributed in the hope that it will be useful,		  //
//  but WITHOUT ANY WARRANTY; without even the implied warranty of	  //
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the	  //
//  GNU General Public License for more details.	      		  //
//									  //
//  You should have received a copy of the GNU General Public License     //
//  along with Pawian.  If not, see <http://www.gnu.org/licenses/>.	  //
//									  //
//************************************************************************//

// VoigtBlattWRelDynamics class definition file. -*- C++ -*-
// Copyright 20213Bertram Kopf

#include <getopt.h>
#include <fstream>
#include <string>
#include <mutex>

#include "PwaUtils/VoigtBlattWRelDynamics.hh"
#include "ErrLogger/ErrLogger.hh"
#include "Particle/Particle.hh"
#include "PwaDynamics/BreitWignerFunction.hh"
#include "Utils/IdStringMapRegistry.hh"
#include "Utils/MathUtils.hh"
#include "Utils/Faddeeva.hh"
#include "qft++Extension/PawianUtils.hh"

VoigtBlattWRelDynamics::VoigtBlattWRelDynamics(std::string& name, std::vector<Particle*>& fsParticles, Particle* mother, std::vector<Particle*>& fsParticlesDaughter1, std::vector<Particle*>& fsParticlesDaughter2, unsigned int prodOrbMom, unsigned int decOrbMom,double qR) :
  BreitWignerBlattWRelDynamics(name, fsParticles, mother, fsParticlesDaughter1, fsParticlesDaughter2, qR)
  ,_massSigmaKey("defaultMassSigmaKey")
  ,_currentSigma(1.)
  ,_prodOrbMom(prodOrbMom)
  ,_decOrbMom(decOrbMom)
{
  _isLdependent=true;
}

VoigtBlattWRelDynamics::~VoigtBlattWRelDynamics()
{
}

complex<double> VoigtBlattWRelDynamics::eval(EvtData* theData, AbsXdecAmp* grandmaAmp, Spin OrbMom){

    complex<double> result(0.,0.);
    //int nConv = 100;

    complex<double> currentProdVal=_grandMaDyn->eval(theData, grandmaAmp, _prodOrbMom); 
    
    double evtMass=theData->DoubleMassId.at(_dynId);
    double evtMassDaughter1=theData->DoubleMassId.at(_dynMassIdDaughter1);
    double evtMassDaughter2=theData->DoubleMassId.at(_dynMassIdDaughter2);
    double xMin = evtMass-3.*_currentWidth;
    double xMax = evtMass+3.*_currentWidth;
    
    int nConv = std::ceil(6.*_currentWidth/0.002);
    
    double step = (xMax-xMin)/((double)nConv);
    double mean = 0.; 

    double y1 = std::norm(currentProdVal*BreitWignerFunction::BlattWRel(_decOrbMom, xMin, _currentMass, _currentWidth, evtMassDaughter1, evtMassDaughter2, _qR)) * MathUtils::Gauss(evtMass-xMin, mean, _currentSigma);
    double y2=0.;
    double xx=0.;
    
   for(int i=0; i<nConv; i++){

     xx = xMin+(i+1)*step;

     y2 = std::norm(currentProdVal*BreitWignerFunction::BlattWRel(_decOrbMom, xx, _currentMass, _currentWidth, evtMassDaughter1, evtMassDaughter2, _qR)) * MathUtils::Gauss(evtMass-xx, mean, _currentSigma);
     result+= complex<double>(0.5*(y2+y1)*step,0.);
	y1=y2;
    }

  return std::sqrt(result);
}

void  VoigtBlattWRelDynamics::fillDefaultParams(std::shared_ptr<AbsPawianParameters> fitPar){
  BreitWignerBlattWRelDynamics::fillDefaultParams(fitPar);
    //fill sigma width
  fitPar->Add(_massSigmaKey, 0.01, 0.4*0.01);
  fitPar->SetLimits(_massSigmaKey, 0., 0.06);
}

void VoigtBlattWRelDynamics::fillParamNameList(){
  _paramNameList.clear();
  BreitWignerBlattWRelDynamics::fillParamNameList();

  //fill sigma width
  _paramNameList.push_back(_massSigmaKey);
}

void VoigtBlattWRelDynamics::updateFitParams(std::shared_ptr<AbsPawianParameters> fitPar){
  BreitWignerBlattWRelDynamics::updateFitParams(fitPar);
  _currentSigma=fitPar->Value(_massSigmaKey);
}

void VoigtBlattWRelDynamics::setMassKey(std::string& theMassKey){
  BreitWignerBlattWRelDynamics::setMassKey(theMassKey);
  _massSigmaKey=theMassKey+"Sigma";
}
