// fitfun.cpp --- Fitfunktionen fuer A und B
// 2004 10 20 Daniela Mund
#include "fitfun.h"

Double_t triggerfun(Double_t *x, Double_t *par)
{
  Double_t q=1.-par[0];
  if (q==0) return 1.;
  Double_t actPhotoElektronen=x[0]*par[1];
  Double_t help=pow(q,actPhotoElektronen);
  Double_t return_value=1.-help*(1.+(actPhotoElektronen*par[1]) /q);
  if (return_value<0) return 0;
  else if (return_value>1) return 1;
  else return return_value;
}

Double_t beta(Double_t EeKin)
{
  // Hierbei ist EeKin die kinetische Energie in keV ohne die Ruhemasse!

  Double_t p = sqrt(EeKin*EeKin + 2*me*EeKin); // Impuls berechnen
  return (p/(EeKin+me));  // relativistische Geschwindigkeit
}

Double_t min_channel(Double_t channel, Double_t sigma)
{
  return channel-3.*sigma;
}

Double_t max_channel(Double_t channel, Double_t sigma)
{
  return channel+3.*sigma;
}

// ---- Verbreiterungsfunktionen

void poisson_broad(TH1F *fromHist, TH1F *toHist,Double_t PhotoElektronen_per_keV)
{
  for (int bin=0; bin<=SPEKLEN; bin++)
  {
    Double_t counts=fromHist->GetBinContent(bin);
    if(counts>0)
    {
      Double_t E_in_keV_bin=fromHist->GetBinCenter(bin); //Energie des theoretischen Spektrums, unabhaengig von meiner Verstaerkung, da THEORIE
      Double_t PhotoElektronenInThisBin=E_in_keV_bin*PhotoElektronen_per_keV;
      //SPEKLEN ist die Achse, Bin SPEKLEN geht von SPEKLEN-1 bis SPEKLEN, BIN SPEKLEN+1 enthaelt Overflows
      for (int PhotoElektronen=0;PhotoElektronen<=PE_SPEKLEN;PhotoElektronen++)
       toHist->Fill(PhotoElektronen,counts*TMath::Poisson(PhotoElektronen,PhotoElektronenInThisBin));
       // 0 PhotoElektronen landen in BIN 1, BIN 0 ist UNDERFLOW
    }
  }
}

void BinomialCorr(TH1F *fromHist)
{
  for (int PhotoElektronenBin=0;PhotoElektronenBin<=PE_SPEKLEN+1;PhotoElektronenBin++)
  {
    // 0 PhotoElektronen landen in BIN 1, BIN 0 ist UNDERFLOW, BIN SPEKLEN+1 enthaelt Overflows
    // => In Bin #Bin sind Bin-1 PhotoElektronen
    Double_t PhotoElektronen=PhotoElektronenBin-1;
    Double_t binomial_corrected=fromHist->GetBinContent(PhotoElektronenBin)*tfun->Eval(PhotoElektronen);
    fromHist->Fill(PhotoElektronen,binomial_corrected); //Fill sucht bin, in dem PhotoElektronen liegt
  }
}


void gaus_broad(TH1F *fromHist, TH1F *toHist,Double_t BinPerPhotoElectron, Double_t SigmaGaus,Double_t sigmaPed, Double_t ChannelOffset)
{
  TH1F *tempHist=(TH1F*)toHist->Clone();
  tempHist->SetNameTitle("temphist","tempHist");
  tempHist->Reset();
  Double_t binWidth=tempHist->GetBinWidth(1);
  for (int PhotoElektronenBin=1;PhotoElektronenBin<=PE_SPEKLEN+1;PhotoElektronenBin++)
  {
    // 0 PhotoElektronen landen in BIN 1, BIN 0 ist UNDERFLOW, BIN SPEKLEN+1 enthaelt Overflows
    // => In Bin #Bin sind Bin-1 PhotoElektronen
    Double_t PhotoElektronen=PhotoElektronenBin-1;
    Double_t sigmaPE=SigmaGaus*sqrt(PhotoElektronen);

    //Umrechnung auf effektives Binning
    Double_t EffectiveChannel=BinPerPhotoElectron*PhotoElektronen;
    Double_t EffectiveMean=EffectiveChannel+ChannelOffset;
    Double_t EffectiveWidth=BinPerPhotoElectron*sigmaPE;
    EffectiveWidth=sqrt(EffectiveWidth*EffectiveWidth+sigmaPed*sigmaPed);

    //Unverbreitert in toHist fuellen
    Double_t intensity=fromHist->GetBinContent(PhotoElektronenBin);
    toHist->Fill(EffectiveMean,intensity);


    for (Double_t channel=min_channel(EffectiveMean,EffectiveWidth); channel<=max_channel(EffectiveMean,EffectiveWidth); channel+=binWidth)
    {
      Double_t myGaus=TMath::Gaus(channel,EffectiveMean,EffectiveWidth,kTRUE);
      tempHist->Fill(channel,myGaus);
    }
  }
  toHist->Multiply(tempHist);
}


// --------- Eichpraeparate
void bi(TH1F *spec,Double_t biscale)
{
    Double_t energy[25] = {     59.13,   67.13,   75.13,   70.44,   78.44,
                               481.65,  540.77,  548.77,  556.77,  552.09,
                               560.09,  554.95,  562.95,  566.54,  569.33,
                               975.64, 1034.76, 1042.76, 1191.64, 1046.08,
                              1054.08, 1048.94, 1056.94, 1060.53, 1063.32};
    Double_t intensity[25] = {0.181, 0.549, 0.415, 0.296, 0.448,
                              1.513, 0.004, 0.013, 0.010, 0.007,
                              0.010, 0.173, 0.261, 0.109, 0.033,
                              8.363, 0.023, 0.070, 0.053, 0.038,
                              0.057, 0.871, 1.318, 0.537, 0.174};


    for (int line=0; line<25; line++) {
 if (line<15) spec->Fill(energy[line],intensity[line]);
	else spec->Fill(energy[line],intensity[line]*biscale);
 }
}

void cer(TH1F *spec)
{
    Double_t energy[14] = { 027.19,      32., 126.926, 154.118, 181.308,
                            186.118, 158.926, 190.928, 159.622, 186.812,
                            191.622, 164.777, 191.967, 196.777};
    Double_t intensity[30] = {  2.92,   1.5935,   14.696,   1.4624, 0.0352,
                              0.4176,   1.1146,   0.0105,   1.7268, 0.09686,
                              0.0522, 0.663, 0.0258, 0.0066};
    for (int line=0; line<14; line++) spec->Fill(energy[line],intensity[line]);
}


void sn(TH1F *spec)
{
  Double_t energy[14] = {19.97, 22.87, 25.77, 22.81, 25.71,
                        363.75,383.72,386.62,389.52,386.56,
                        389.46,387.72,390.62,391.07};
  Double_t intensity[14] = { 0.036,1.035,7.446,0.273,3.922,
                            24.628,0.012,0.354,2.546,0.093,
                             1.341,0.356,5.117,1.323};

  for (int line=0; line<14; line++)
    spec->Fill(energy[line],intensity[line]);
}



// fitfun erstellen

TH1F *fitspec( Double_t *par)
{
  Double_t BinPerPhotoElectron=par[0];
  Double_t PhotoElektronen_per_keV=par[1];
  Double_t SigmaGaus=par[2];
  Double_t sigmaPed=par[3];
  Double_t ChannelOffset=par[4];
  Double_t biscale=par[5];
  Double_t norm=par[6];

  TH1F *theor_spec=new TH1F ("theor_spec","Theoretisches Spektrum",SPEKLEN,0,MAX_ENERGY);

  TH1F *poisson_spec=new TH1F ("poisson_spec","Theoretisches Spektrum mit Poisson-Verbreiterung", PE_SPEKLEN, 0, PE_SPEKLEN);
  TH1F *gaus_spec=new TH1F ("gaus_spec","Theoretisches Spektrum mit allen Verbreiterungen, x skaliert", SPEKLEN, 0, SPEKLEN);
  bool again=false;
  do
  {
    switch (FitFlag)
    {
      case 1:  bi(theor_spec,biscale);again=false;break;
      case 2:  cer(theor_spec);again=false;break;
      case 3:  sn(theor_spec);again=false;break;
      default: cout<<"Falsche Eingabe fuer Fitspektrum, waehle neues:";
              cin>>FitFlag;
              again=true;
    }
  } while(again);

  poisson_broad(theor_spec,poisson_spec,PhotoElektronen_per_keV);
  BinomialCorr(poisson_spec);
  gaus_broad(poisson_spec, gaus_spec,BinPerPhotoElectron,SigmaGaus,sigmaPed,ChannelOffset);
  gaus_spec->Scale(norm);
  delete poisson_spec;
  delete theor_spec;

  return gaus_spec;
}


// Klasse  MyFit zum Erstellen einer Fitfunktion


MyFit::MyFit(Double_t *parameter, int FunFlag)
{
  for (int par=0;par<PAR_NUM;par++)
  {
    old_parameter[par]=parameter[par];
    fix_parameter[par]=true;
  }

  theor_spec=new TH1F ("theor_spec","Theoretisches Spektrum",SPEKLEN,0,MAX_ENERGY);
  pe_spec=new TH1F ("pe_spec","Theoretisches Photoelektornrn-Spektrum", PE_SPEKLEN, 0, PE_SPEKLEN);
  effective_spec=new TH1F ("effective_spec","Theoretisches Spektrum mit allen Verbreiterungen, x skaliert", SPEKLEN, 0, SPEKLEN);
  SelectTheory(FunFlag,old_parameter);
  ChangeHistos(old_parameter);
}


Double_t MyFit::GetHistoBin(Double_t *x, Double_t *parameter)
{
  bool change=true;
  for (int par=0;par<PAR_NUM;par++)
  {
    if(old_parameter[par]==parameter[par])
      change=false;
    if(change)
    {
      ChangeHistos(parameter);
      break;
    }

  }
  if (change)
    for (int par=0;par<PAR_NUM;par++)
      old_parameter[par]=parameter[par];
  Double_t val=effective_spec->GetBinContent(x[0]);
  return val;
}


void MyFit::ChangeHistos(Double_t *parameter)
{
  theor_spec->Reset();
  pe_spec->Reset();
  effective_spec->Reset();
  SelectTheory(flag,parameter);
  poisson_broad(theor_spec,pe_spec,parameter[1]);
  BinomialCorr(pe_spec);
  gaus_broad(pe_spec,effective_spec,parameter[0],parameter[2],parameter[3],parameter[4]);
  effective_spec->Scale(parameter[6]);
}


void MyFit::SelectTheory(int MyFlag, Double_t *newParameter)
{
  flag=MyFlag;
  bool again=false;
  do
  {
    switch (flag)
    {
      case 1:  bi(theor_spec,newParameter[5]);again=false;break;
      case 2:  cer(theor_spec);again=false;break;
      case 3:  sn(theor_spec);again=false;break;
      default: cout<<"Falsche Eingabe fuer Fitspektrum, waehle neues:";
               cin>>flag;
               again=true;
    }
  } while(again);
}

TF1 *MyFit::GetFunction(void)
{
  TF1 *fitfun=new TF1("fitfun",GetHistoBin,0.,(double)SPEKLEN,PAR_NUM);
  return fitfun;
}