/*****************************************************************************
 * Project: RooFit                                                           *
 * Package: RooFitModels                                                     *
 * @(#)root/roofit:$Id$
 * Authors:                                                                  *
 *   WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu       *
 *   DK, David Kirkby,    UC Irvine,         dkirkby@uci.edu                 *
 *                                                                           *
 * Copyright (c) 2000-2005, 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
// P.d.f implementing the Crystall Ball line shape
// END_HTML
//

#include "RooFit.h"

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

#include "RooCBShape.h"
#include "RooAbsReal.h"
#include "RooRealVar.h"
#include "RooMath.h"
#include "TMath.h"

#include "TError.h"

using namespace std;

ClassImp(RooCBShape)
;

  

//_____________________________________________________________________________
Double_t RooCBShape::ApproxErf(Double_t arg) const 
{
  static const double erflim = 5.0;
  if( arg > erflim )
    return 1.0;
  if( arg < -erflim )
    return -1.0;
  
  return RooMath::erf(arg);
}



//_____________________________________________________________________________
RooCBShape::RooCBShape(const char *name, const char *title,
		       RooAbsReal& _m, RooAbsReal& _m0, RooAbsReal& _sigma,
		       RooAbsReal& _alpha, RooAbsReal& _n) :
  RooAbsPdf(name, title),
  m("m", "Dependent", this, _m),
  m0("m0", "M0", this, _m0),
  sigma("sigma", "Sigma", this, _sigma),
  alpha("alpha", "Alpha", this, _alpha),
  n("n", "Order", this, _n)
{
}


//_____________________________________________________________________________
RooCBShape::RooCBShape(const RooCBShape& other, const char* name) :
  RooAbsPdf(other, name), m("m", this, other.m), m0("m0", this, other.m0),
  sigma("sigma", this, other.sigma), alpha("alpha", this, other.alpha),
  n("n", this, other.n)
{
}


//_____________________________________________________________________________
Double_t RooCBShape::evaluate() const {

  Double_t t = (m-m0)/sigma;
  if (alpha < 0) t = -t;

  Double_t absAlpha = fabs((Double_t)alpha);

  if (t >= -absAlpha) {
    return exp(-0.5*t*t);
  }
  else {
    Double_t a =  TMath::Power(n/absAlpha,n)*exp(-0.5*absAlpha*absAlpha);
    Double_t b= n/absAlpha - absAlpha; 

    return a/TMath::Power(b - t, n);
  }
}


//_____________________________________________________________________________
Int_t RooCBShape::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
{
  if( matchArgs(allVars,analVars,m) )
    return 1 ;
  
  return 0;
}



//_____________________________________________________________________________
Double_t RooCBShape::analyticalIntegral(Int_t code, const char* rangeName) const
{
  static const double sqrtPiOver2 = 1.2533141373;
  static const double sqrt2 = 1.4142135624;

  R__ASSERT(code==1);
  double result = 0.0;
  bool useLog = false;
  
  if( fabs(n-1.0) < 1.0e-05 )
    useLog = true;
  
  double sig = fabs((Double_t)sigma);
  
  double tmin = (m.min(rangeName)-m0)/sig;
  double tmax = (m.max(rangeName)-m0)/sig;
  
  if(alpha < 0) {
    double tmp = tmin;
    tmin = -tmax;
    tmax = -tmp;
  }

  double absAlpha = fabs((Double_t)alpha);
  
  if( tmin >= -absAlpha ) {
    result += sig*sqrtPiOver2*(   ApproxErf(tmax/sqrt2)
                                - ApproxErf(tmin/sqrt2) );
  }
  else if( tmax <= -absAlpha ) {
    double a = TMath::Power(n/absAlpha,n)*exp(-0.5*absAlpha*absAlpha);
    double b = n/absAlpha - absAlpha;
    
    if(useLog) {
      result += a*sig*( log(b-tmin) - log(b-tmax) );
    }
    else {
      result += a*sig/(1.0-n)*(   1.0/(TMath::Power(b-tmin,n-1.0))
                                - 1.0/(TMath::Power(b-tmax,n-1.0)) );
    }
  }
  else {
    double a = TMath::Power(n/absAlpha,n)*exp(-0.5*absAlpha*absAlpha);
    double b = n/absAlpha - absAlpha;
    
    double term1 = 0.0;
    if(useLog) {
      term1 = a*sig*(  log(b-tmin) - log(n/absAlpha));
    }
    else {
      term1 = a*sig/(1.0-n)*(   1.0/(TMath::Power(b-tmin,n-1.0))
                              - 1.0/(TMath::Power(n/absAlpha,n-1.0)) );
    }
    
    double term2 = sig*sqrtPiOver2*(   ApproxErf(tmax/sqrt2)
                                     - ApproxErf(-absAlpha/sqrt2) );
    
    
    result += term1 + term2;
  }
  
  return result;
}



//_____________________________________________________________________________
Int_t RooCBShape::getMaxVal(const RooArgSet& vars) const 
{
  // Advertise that we know the maximum of self for given (m0,alpha,n,sigma)
  RooArgSet dummy ;

  if (matchArgs(vars,dummy,m)) {
    return 1 ;  
  }
  return 0 ;  
}



//_____________________________________________________________________________
Double_t RooCBShape::maxVal(Int_t code) const
{
  R__ASSERT(code==1) ;

  // The maximum value for given (m0,alpha,n,sigma)
  return 1.0 ;
}


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