///////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////Bukin Class: BukinFn                                                                        
///////                                                                                                            
///////This is Bukin's complete function for the shape of  the omega and  
////// phi peak.                                                                                              
///////                                                                                                            
///////It has 8 parameters: N, mu, c, d, Mx, Gammax, sigma and phi.      
///////It also has a background function with another 5 parameters.          
///////                                                                                                                                                
////////////////////////////////////////////////////////////////////////////////////////////////////////////// 

#include "TMath.h" 
#include "TComplex.h" 
#include <iostream>

//Define class Bukin 
class Bukin { 

  //Constructor 


  Double_t oMASS; 
  Double_t pMASS; 
  Double_t oWIDTH; 
  Double_t pWIDTH; 
  
  Double_t rho1; 
  Double_t rho2; 
  Double_t rho; 
  
  Double_t bukin; 
  Double_t bg; 

  Double_t xMASS; 
  Double_t xWIDTH; 
  
  Double_t N; 
  Double_t mu; 
  Double_t c; 
  Double_t d; 
  Double_t sigma; 
  Double_t phi; 
  
  Double_t b0; 
  Double_t b1; 
  Double_t b2; 
  Double_t b3; 
  Double_t b4; 
  Double_t b5; 

  Double_t rhosinsq; 
  Double_t aRe; 
  Double_t bRe; 

  Double_t x1; 
  Double_t x2; 
  Double_t x3; 
  Double_t x4; 
  Double_t x5; 
  
//I don't know why but it crashes if 'all' is defined here. 
//Instead it's defined in Load(). 
  TComplex all;  
  TComplex a; 
  TComplex b; 
  TComplex po; 
  TComplex px; 
  TComplex op; 
  TComplex ox; 
  TComplex xp; 
  TComplex xo; 
  TComplex ptemp; 
  TComplex otemp; 
  TComplex xtemp; 
  TComplex cd; 
  TComplex cminusd; 
  TComplex i; 
  TComplex sqbracket1; 
  TComplex sqbracket2; 
  TComplex sqbracket3; 


public: 


   Bukin() { 
      std::cout << "constructor done" << std::endl;  
   }

  //Function Load(): Does anything that doesn't depend on the parameters. 
  //Load() is performed once before the fitting begins. 
  void Load(){ 
    
    oMASS = 0.782; 
    pMASS = 1.020; 
    oWIDTH = 0.008*0.5; 
    pWIDTH = 0.004*0.5; 

    rho1 = (4.0*TMath::Power(pMASS-oMASS,2.0))+TMath::Power(pWIDTH+oWIDTH,2.0); 
    rho2 = (2.0*(pWIDTH+oWIDTH)*TMath::Sqrt(pWIDTH*oWIDTH)); 
    rho = (rho1+rho2)/(rho1-rho2); 
  
    bukin = 0.0; 
    all = TComplex(); 
    bg = 0.0; 
    
  } 

  //Function BukinFn(): returns the y-value corresponding to a given x. 
  //The function describes the shape of the omega and phi peaks. 
  Double_t BukinFn(Double_t *x, Double_t *par) { 
    
    xMASS = par[4]; 
    xWIDTH = par[5]*0.5; 
    
    N = par[0]; 
    mu = par[1]; 
    c = par[2]; 
    d = par[3]; 
    sigma = par[6]; 
    phi = par[7]; 

    b0 = par[8]; 
    b1 = par[9]; 
    b2 = par[10]; 
    b3 = par[11]; 
    b4 = par[12]; 
    b5 = par[13]; 

    rhosinsq = rho*TMath::Sin(mu)*TMath::Sin(mu); 
    
    aRe = TMath::Sqrt(pWIDTH/8.0*TMath::Pi())*(TMath::Sqrt(1.0+rhosinsq)+TMath::Cos(mu)); 
    bRe = TMath::Sqrt(oWIDTH/8.0*TMath::Pi())*(TMath::Sqrt(1.0+rhosinsq)-TMath::Cos(mu)); 
    
    b = TComplex(bRe, 0); 
    a = TComplex(aRe*TMath::Cos(phi), aRe*TMath::Sin(phi)); 
    
    po = TComplex(pMASS-oMASS, pWIDTH+oWIDTH); 
    px = TComplex(pMASS-xMASS, pWIDTH+xWIDTH); 
    op = TComplex(oMASS-pMASS, oWIDTH+pWIDTH); 
    ox = TComplex(oMASS-xMASS, oWIDTH+xWIDTH); 
    xp = TComplex(xMASS-pMASS, xWIDTH+pWIDTH); 
    xo = TComplex(xMASS-oMASS, xWIDTH+oWIDTH); 
    
    ptemp = TComplex((1.17741*sigma)+pWIDTH, x[0]-pMASS); 
    otemp = TComplex((1.17741*sigma)+oWIDTH, x[0]-oMASS); 
    xtemp = TComplex((1.17741*sigma)+xWIDTH, x[0]-xMASS); 
    
    cd = TComplex(c, d); 
    cminusd = TComplex(c, d*(-1)); 

    i = TComplex(0, 1); 

    sqbracket1 = TComplex(); 
    sqbracket2 = TComplex(); 
    sqbracket3 = TComplex();  
        
    sqbracket1 = (a/(2*pWIDTH)) - ((i*b)/po) + ((i*cd)/px); 
    sqbracket2 = ((i*a)/op) - (b/(2*oWIDTH)) + ((i*cd)/ox); 
    sqbracket3 = ((i*a)/xp) - ((i*b)/xo) + (cd/(2*xWIDTH)); 
        
    all = ((a/ptemp)*sqbracket1) - ((b/otemp)*sqbracket2) + ((cminusd/xtemp)/sqbracket3); 
        
    bukin = N*2*all.Re(); 
          
    //Background function: 6 parameters; b0, b1, b2, b3, b4, b5. 
    
    x1 = x[0]-1; 
    x2 = x1*x1; 
    x3 = x2*x1; 
    x4 = x3*x1; 
    x5 = x4*x1; 
    
    bg = TMath::Exp(b0 + (b1*x1) + (b2*x2) + (b3*x3) + (b4*x4) + (b5*x5)); 
        
    return bukin + bg; 
    
  } 
};