/*****************************************************************************
 * Project: RooFit                                                           *
 * Package: RooFitModels                                                     *
 * @(#)root/roofit:$Id$
 * Authors:                                                                  *
 *   DB, Dieter Best,     UC Irvine,         best@slac.stanford.edu          *
 *   HT, Hirohisa Tanaka  SLAC               tanaka@slac.stanford.edu        *
 *                                                                           *
 *   Updated version with analytical integral                                *
 *   MP, Marko Petric,   J. Stefan Institute,  marko.petric@ijs.si           *
 *                                                                           *
 * Copyright (c) 2000-2013, Regents of the University of California          *
 *                          and Stanford University. All rights reserved.    *
 *                                                                           *
 * Redistribution and use in source and binary forms,                        *
 * with or without modification, are permitted according to the terms        *
 * listed in LICENSE (http://roofit.sourceforge.net/license.txt)             *
 *****************************************************************************/

//////////////////////////////////////////////////////////////////////////////
//
// BEGIN_HTML
//
// RooNovosibirsk implements the Novosibirsk function 
//
// Function taken from H. Ikeda et al. NIM A441 (2000), p. 401 (Belle Collaboration)
// 
// END_HTML
//


#include "RooFit.h"

#include <math.h>
#include "TMath.h"

#include "RooNovosibirsk.h"
#include "RooRealVar.h"

using namespace std;

ClassImp(RooNovosibirsk)



//_____________________________________________________________________________
RooNovosibirsk::RooNovosibirsk(const char *name, const char *title,
			     RooAbsReal& _x,     RooAbsReal& _peak,
			     RooAbsReal& _width, RooAbsReal& _tail) :
  // The two addresses refer to our first dependent variable and
  // parameter, respectively, as declared in the rdl file
  RooAbsPdf(name, title),
  x("x","x",this,_x),
  width("width","width",this,_width),
  peak("peak","peak",this,_peak),
  tail("tail","tail",this,_tail)
{
}



//_____________________________________________________________________________
RooNovosibirsk::RooNovosibirsk(const RooNovosibirsk& other, const char *name):
  RooAbsPdf(other,name),
  x("x",this,other.x),
  width("width",this,other.width),
  peak("peak",this,other.peak),
  tail("tail",this,other.tail)
{
}



//_____________________________________________________________________________
Double_t RooNovosibirsk::evaluate() const
{ 
  //If tail=eta=0 the Belle distribution becomes gaussian
  if (TMath::Abs(tail) < 1.e-7) {
    return TMath::Exp( -0.5 * TMath::Power( ( (x - peak) / width ), 2 ));
  }
  
  Double_t arg = 1.0 - ( x - peak ) * tail / width;
  
  if (arg < 1.e-7) {
    //Argument of logaritem negative. Real continuation -> function equals zero
    return 0.0;
  }
    
  Double_t log = TMath::Log(arg);    
  static const Double_t xi = 2.3548200450309494; // 2 Sqrt( Ln(4) )
    
  Double_t width_zero = ( 2.0 / xi ) * TMath::ASinH( tail * xi * 0.5 );
  Double_t width_zero2 = width_zero * width_zero;
  Double_t exponent = ( -0.5 / (width_zero2) * log * log ) - ( width_zero2 * 0.5 );
    
  return TMath::Exp(exponent) ;
}




//_____________________________________________________________________________
Int_t RooNovosibirsk::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* ) const
{
  if (matchArgs(allVars,analVars,x)) return 1 ;
  if (matchArgs(allVars,analVars,peak)) return 2 ;
    
  //The other two integrals over tali and width are not integrable
    
  return 0 ;
}



//_____________________________________________________________________________
Double_t RooNovosibirsk::analyticalIntegral(Int_t code, const char* rangeName) const
{
  assert(code==1 || code==2) ;
  
  //The range is defined as [A,B]
  
  //Numerical values need for the evaluation of the integral
  static const Double_t sqrt2 = 1.4142135623730950; // Sqrt(2)
  static const Double_t sqlog2 = 0.832554611157697756; //Sqrt( Log(2) )
  static const Double_t sqlog4 = 1.17741002251547469; //Sqrt( Log(4) )
  static const Double_t log4 = 1.38629436111989062; //Log(2)
  static const Double_t rootpiby2 = 1.2533141373155003; // Sqrt(pi/2)
  static const Double_t sqpibylog2 = 2.12893403886245236; //Sqrt( pi/Log(2) )
    
  if (code==1) {
    Double_t A = x.min(rangeName);
    Double_t B = x.max(rangeName);
        
    Double_t result = 0;

        
    //If tail==0 the function becomes gaussian, thus we return a gassian integral
    if (TMath::Abs(tail) < 1.e-7) {

      Double_t xscale = sqrt2*width;
            
      result = rootpiby2*width*(TMath::Erf((B-peak)/xscale)-TMath::Erf((A-peak)/xscale));
            
      return result;
            
    }
        
    // If the range is not defined correctly the function becomes complex
    Double_t log_argument_A = ( (peak - A)*tail + width ) / width ;
    Double_t log_argument_B = ( (peak - B)*tail + width ) / width ;
        
    //lower limit
    if ( log_argument_A < 1.e-7) {
      log_argument_A = 1.e-7;
    }
        
    //upper limit
    if ( log_argument_B < 1.e-7) {
      log_argument_B = 1.e-7;
    }
        
    Double_t term1 =  TMath::ASinH( tail * sqlog4 );
    Double_t term1_2 =  term1 * term1;
        
    //Calculate the error function arguments
    Double_t erf_termA = ( term1_2 - log4 * TMath::Log( log_argument_A ) ) / ( 2 * term1 * sqlog2 );
    Double_t erf_termB = ( term1_2 - log4 * TMath::Log( log_argument_B ) ) / ( 2 * term1 * sqlog2 );
        
    result = 0.5 / tail * width * term1 * ( TMath::Erf(erf_termB) - TMath::Erf(erf_termA)) * sqpibylog2;
        
    return result;
        
  } else if (code==2) {
    Double_t A = x.min(rangeName);
    Double_t B = x.max(rangeName);
        
    Double_t result = 0;
        
        
    //If tail==0 the function becomes gaussian, thus we return a gassian integral
    if (TMath::Abs(tail) < 1.e-7) {
            
      Double_t xscale = sqrt2*width;
            
      result = rootpiby2*width*(TMath::Erf((B-x)/xscale)-TMath::Erf((A-x)/xscale));
            
      return result;
            
    }
        
    // If the range is not defined correctly the function becomes complex
    Double_t log_argument_A = ( (A - x)*tail + width ) / width;
    Double_t log_argument_B = ( (B - x)*tail + width ) / width;
        
    //lower limit
    if ( log_argument_A < 1.e-7) {
      log_argument_A = 1.e-7;
    }
        
    //upper limit
    if ( log_argument_B < 1.e-7) {
      log_argument_B = 1.e-7;
    }
        
    Double_t term1 =  TMath::ASinH( tail * sqlog4 );
    Double_t term1_2 =  term1 * term1;
        
    //Calculate the error function arguments
    Double_t erf_termA = ( term1_2 - log4 * TMath::Log( log_argument_A ) ) / ( 2 * term1 * sqlog2 );
    Double_t erf_termB = ( term1_2 - log4 * TMath::Log( log_argument_B ) ) / ( 2 * term1 * sqlog2 );
        
    result = 0.5 / tail * width * term1 * ( TMath::Erf(erf_termB) - TMath::Erf(erf_termA)) * sqpibylog2;
        
    return result;
        
  }

  // Emit fatal error
  coutF(Eval) << "Error in RooNovosibirsk::analyticalIntegral" << std::endl;  

  // Put dummy return here to avoid compiler warnings
  return 1.0 ;
}
 RooNovosibirsk.cxx:1
 RooNovosibirsk.cxx:2
 RooNovosibirsk.cxx:3
 RooNovosibirsk.cxx:4
 RooNovosibirsk.cxx:5
 RooNovosibirsk.cxx:6
 RooNovosibirsk.cxx:7
 RooNovosibirsk.cxx:8
 RooNovosibirsk.cxx:9
 RooNovosibirsk.cxx:10
 RooNovosibirsk.cxx:11
 RooNovosibirsk.cxx:12
 RooNovosibirsk.cxx:13
 RooNovosibirsk.cxx:14
 RooNovosibirsk.cxx:15
 RooNovosibirsk.cxx:16
 RooNovosibirsk.cxx:17
 RooNovosibirsk.cxx:18
 RooNovosibirsk.cxx:19
 RooNovosibirsk.cxx:20
 RooNovosibirsk.cxx:21
 RooNovosibirsk.cxx:22
 RooNovosibirsk.cxx:23
 RooNovosibirsk.cxx:24
 RooNovosibirsk.cxx:25
 RooNovosibirsk.cxx:26
 RooNovosibirsk.cxx:27
 RooNovosibirsk.cxx:28
 RooNovosibirsk.cxx:29
 RooNovosibirsk.cxx:30
 RooNovosibirsk.cxx:31
 RooNovosibirsk.cxx:32
 RooNovosibirsk.cxx:33
 RooNovosibirsk.cxx:34
 RooNovosibirsk.cxx:35
 RooNovosibirsk.cxx:36
 RooNovosibirsk.cxx:37
 RooNovosibirsk.cxx:38
 RooNovosibirsk.cxx:39
 RooNovosibirsk.cxx:40
 RooNovosibirsk.cxx:41
 RooNovosibirsk.cxx:42
 RooNovosibirsk.cxx:43
 RooNovosibirsk.cxx:44
 RooNovosibirsk.cxx:45
 RooNovosibirsk.cxx:46
 RooNovosibirsk.cxx:47
 RooNovosibirsk.cxx:48
 RooNovosibirsk.cxx:49
 RooNovosibirsk.cxx:50
 RooNovosibirsk.cxx:51
 RooNovosibirsk.cxx:52
 RooNovosibirsk.cxx:53
 RooNovosibirsk.cxx:54
 RooNovosibirsk.cxx:55
 RooNovosibirsk.cxx:56
 RooNovosibirsk.cxx:57
 RooNovosibirsk.cxx:58
 RooNovosibirsk.cxx:59
 RooNovosibirsk.cxx:60
 RooNovosibirsk.cxx:61
 RooNovosibirsk.cxx:62
 RooNovosibirsk.cxx:63
 RooNovosibirsk.cxx:64
 RooNovosibirsk.cxx:65
 RooNovosibirsk.cxx:66
 RooNovosibirsk.cxx:67
 RooNovosibirsk.cxx:68
 RooNovosibirsk.cxx:69
 RooNovosibirsk.cxx:70
 RooNovosibirsk.cxx:71
 RooNovosibirsk.cxx:72
 RooNovosibirsk.cxx:73
 RooNovosibirsk.cxx:74
 RooNovosibirsk.cxx:75
 RooNovosibirsk.cxx:76
 RooNovosibirsk.cxx:77
 RooNovosibirsk.cxx:78
 RooNovosibirsk.cxx:79
 RooNovosibirsk.cxx:80
 RooNovosibirsk.cxx:81
 RooNovosibirsk.cxx:82
 RooNovosibirsk.cxx:83
 RooNovosibirsk.cxx:84
 RooNovosibirsk.cxx:85
 RooNovosibirsk.cxx:86
 RooNovosibirsk.cxx:87
 RooNovosibirsk.cxx:88
 RooNovosibirsk.cxx:89
 RooNovosibirsk.cxx:90
 RooNovosibirsk.cxx:91
 RooNovosibirsk.cxx:92
 RooNovosibirsk.cxx:93
 RooNovosibirsk.cxx:94
 RooNovosibirsk.cxx:95
 RooNovosibirsk.cxx:96
 RooNovosibirsk.cxx:97
 RooNovosibirsk.cxx:98
 RooNovosibirsk.cxx:99
 RooNovosibirsk.cxx:100
 RooNovosibirsk.cxx:101
 RooNovosibirsk.cxx:102
 RooNovosibirsk.cxx:103
 RooNovosibirsk.cxx:104
 RooNovosibirsk.cxx:105
 RooNovosibirsk.cxx:106
 RooNovosibirsk.cxx:107
 RooNovosibirsk.cxx:108
 RooNovosibirsk.cxx:109
 RooNovosibirsk.cxx:110
 RooNovosibirsk.cxx:111
 RooNovosibirsk.cxx:112
 RooNovosibirsk.cxx:113
 RooNovosibirsk.cxx:114
 RooNovosibirsk.cxx:115
 RooNovosibirsk.cxx:116
 RooNovosibirsk.cxx:117
 RooNovosibirsk.cxx:118
 RooNovosibirsk.cxx:119
 RooNovosibirsk.cxx:120
 RooNovosibirsk.cxx:121
 RooNovosibirsk.cxx:122
 RooNovosibirsk.cxx:123
 RooNovosibirsk.cxx:124
 RooNovosibirsk.cxx:125
 RooNovosibirsk.cxx:126
 RooNovosibirsk.cxx:127
 RooNovosibirsk.cxx:128
 RooNovosibirsk.cxx:129
 RooNovosibirsk.cxx:130
 RooNovosibirsk.cxx:131
 RooNovosibirsk.cxx:132
 RooNovosibirsk.cxx:133
 RooNovosibirsk.cxx:134
 RooNovosibirsk.cxx:135
 RooNovosibirsk.cxx:136
 RooNovosibirsk.cxx:137
 RooNovosibirsk.cxx:138
 RooNovosibirsk.cxx:139
 RooNovosibirsk.cxx:140
 RooNovosibirsk.cxx:141
 RooNovosibirsk.cxx:142
 RooNovosibirsk.cxx:143
 RooNovosibirsk.cxx:144
 RooNovosibirsk.cxx:145
 RooNovosibirsk.cxx:146
 RooNovosibirsk.cxx:147
 RooNovosibirsk.cxx:148
 RooNovosibirsk.cxx:149
 RooNovosibirsk.cxx:150
 RooNovosibirsk.cxx:151
 RooNovosibirsk.cxx:152
 RooNovosibirsk.cxx:153
 RooNovosibirsk.cxx:154
 RooNovosibirsk.cxx:155
 RooNovosibirsk.cxx:156
 RooNovosibirsk.cxx:157
 RooNovosibirsk.cxx:158
 RooNovosibirsk.cxx:159
 RooNovosibirsk.cxx:160
 RooNovosibirsk.cxx:161
 RooNovosibirsk.cxx:162
 RooNovosibirsk.cxx:163
 RooNovosibirsk.cxx:164
 RooNovosibirsk.cxx:165
 RooNovosibirsk.cxx:166
 RooNovosibirsk.cxx:167
 RooNovosibirsk.cxx:168
 RooNovosibirsk.cxx:169
 RooNovosibirsk.cxx:170
 RooNovosibirsk.cxx:171
 RooNovosibirsk.cxx:172
 RooNovosibirsk.cxx:173
 RooNovosibirsk.cxx:174
 RooNovosibirsk.cxx:175
 RooNovosibirsk.cxx:176
 RooNovosibirsk.cxx:177
 RooNovosibirsk.cxx:178
 RooNovosibirsk.cxx:179
 RooNovosibirsk.cxx:180
 RooNovosibirsk.cxx:181
 RooNovosibirsk.cxx:182
 RooNovosibirsk.cxx:183
 RooNovosibirsk.cxx:184
 RooNovosibirsk.cxx:185
 RooNovosibirsk.cxx:186
 RooNovosibirsk.cxx:187
 RooNovosibirsk.cxx:188
 RooNovosibirsk.cxx:189
 RooNovosibirsk.cxx:190
 RooNovosibirsk.cxx:191
 RooNovosibirsk.cxx:192
 RooNovosibirsk.cxx:193
 RooNovosibirsk.cxx:194
 RooNovosibirsk.cxx:195
 RooNovosibirsk.cxx:196
 RooNovosibirsk.cxx:197
 RooNovosibirsk.cxx:198
 RooNovosibirsk.cxx:199
 RooNovosibirsk.cxx:200
 RooNovosibirsk.cxx:201
 RooNovosibirsk.cxx:202
 RooNovosibirsk.cxx:203
 RooNovosibirsk.cxx:204
 RooNovosibirsk.cxx:205
 RooNovosibirsk.cxx:206
 RooNovosibirsk.cxx:207
 RooNovosibirsk.cxx:208
 RooNovosibirsk.cxx:209
 RooNovosibirsk.cxx:210
 RooNovosibirsk.cxx:211
 RooNovosibirsk.cxx:212
 RooNovosibirsk.cxx:213
 RooNovosibirsk.cxx:214
 RooNovosibirsk.cxx:215
 RooNovosibirsk.cxx:216
 RooNovosibirsk.cxx:217
 RooNovosibirsk.cxx:218
 RooNovosibirsk.cxx:219
 RooNovosibirsk.cxx:220
 RooNovosibirsk.cxx:221
 RooNovosibirsk.cxx:222
 RooNovosibirsk.cxx:223