//************************************************************************//
//                                                                        //
//  Copyright 2013 Bertram Kopf (bertram@ep1.rub.de)                      //
//               Julian Pychy (julian@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/>.       //
//                                                                        //
//************************************************************************//

// ------------------------------------------------------------------------
// Andreas Pitka 2018
// For Details about the Johnson Function:
// N. L. Johnson - Systems of Frequency Curves Generated by Methods of Translation, Biometrika Vol. 36, No. 1/2 (Jun., 1949), pp. 149-176
// CDF Collaboration: T. Aaltonen, et al, Phys. Rev. D 85, 012009 – Published 27 January 2012 
//  - (DOI: 10.1103/PhysRevD.85.012009)
//  - https://arxiv.org/abs/1111.5023v1
//  (be carefull, the normalization is wrong in this paper, the factor delta/(sigma*sqrt(2pi)) is missing)
// wikipedia: https://en.wikipedia.org/wiki/Johnson%27s_SU-distribution
 

#include "PwaDynamics/Johnson.hh"

#include <cmath>

Johnson::Johnson() {
}

Johnson::~Johnson() {
}

double Johnson::calc(double currentMass, double mass0, double sigma,
		     double delta, double gamma) {
  double result = 0.;
  if ((sigma < 0.) || ( fabs(sigma)<1e-20 ) ) {
    return result;  // Not meant to be for those who want to be thinner than 0
   }
  if ( fabs(sigma) < 1.e-20 ) {
    sigma=1.e-12;
  }

  double z      = (currentMass - mass0)/fabs(sigma);
  double arg    = -0.5 * pow( gamma + fabs(delta) * asinh( z ), 2. ); 

  double norm   = fabs(delta) / ( sqrt(2.*M_PI) * fabs(sigma) * sqrt( 1. + pow(z,2.) ) );

  result = norm * exp( arg );

  return result;
}
