/*****************************************
histogram fit function.
 ****************************************/

//returns value of last paramter....as i can't set 
//it as the fit parameter,
//after the fit has taken place. YUK.

Double_t lumiFit( TH1 *hist , TString fitName,Int_t colNum){
  TF1 *fit= new TF1(fitName.Data(),parFit1,0.9,0.999,3);
  fit->SetParameters(1.3e-3,1,-1);
  //fit->FixParameter(4,0);
  fit->SetParNames("a","b","c");
  fit->SetLineColor(colNum);
  hist->Fit(fitName.Data(),"R");

  //calculate d
  cout <<"\n last bin: "<< hist->GetBinContent(102)<< endl;
  cout <<"no of bins: " << hist->GetNbinsX()<< endl;
  cout <<"integral, full: " <<hist->Integral()<<" to last bin: "<<hist->Integral(0,hist->GetNbinsX()-1) << endl;

  // as integral is normalised, d is 1-integral of all but last
  Double_t d1 = 1-(hist->Integral(0, hist->GetNbinsX()-1 ));

  fit->SetParameter(4,d1);
  cout <<" d is " <<d1 <<endl;
  return d1;
}


//this one uses extrapolation
Double_t lumiFit_old( TH1 *hist , TString fitName,Int_t colNum){
  TF1 *fit= new TF1(fitName.Data(),parFit,0.9,0.999,4);
  fit->SetParameters(1.3e-3,1,-1,0);
  //fit->FixParameter(4,0);
  fit->SetParNames("a","b","c","d");
  fit->SetLineColor(colNum);
  hist->Fit(fitName.Data(),"R");

  //calculate d
  Double_t x1[]={0.9985};
  Double_t x2[]={0.9995};
  Double_t mypars[]={fit->GetParameter(0),
           fit->GetParameter(1),
           fit->GetParameter(2),
           0};

  Double_t extrap=parFit(x2,mypars)+(parFit(x2,mypars)-parFit(x1,mypars))/(x2-x1);

  cout <<"\n last bin: "<< hist->GetBinContent(101)<< endl;
  cout <<" extrap: "<< extrap << endl;

  Double_t d1 = hist->GetBinContent(101)- extrap;

  fit->SetParameter(4,d1);
  cout <<"d = " <<fit->GetParameter(3) <<"d is " <<d1 <<endl;
  return d1;
}


/******************************************
User-defined functions for fitting.
 *****************************************/
Double_t expFit(Double_t *x, Double_t par[4]){
  Double_t retval= par[0]*TMath::Exp(par[1]*x[0])+par[2]*TMath::Exp(par[3]*x[0])+par[4];
  if(x[0]<1) return retval;
  else return 0;
}

Double_t expFit1(Double_t *x, Double_t par[1]){
  Double_t retval= par[0]*TMath::Exp(par[1]*x[0])+par[1];
  return retval;
}

Double_t parFit1(Double_t *x, Double_t par[2]){
  Double_t val=par[0]*pow(x[0],par[1])*pow(1-x[0],par[2]);
  return val;
}

Double_t parFit(Double_t *x, Double_t par[3]){
  Double_t epsilon=0;
  Double_t val=par[0]*pow(x[0],par[1])*pow(1+epsilon-x[0],par[2])+par[3]*diracDelta(1-x[0]);
  return val;
}

Double_t diracDelta(Double_t y){
  if(y==0)return 1;
  else return 0;
}

/*****************************************
User-defined stats moments function
 *****************************************/
Double_t histMoment(TH1 *h,Int_t m){
  const Int_t n=h.GetNbinsX();
  Double_t mom=0.0;
  Double_t s=0.0;
  const Double_t xmean=h.GetMean();
  //get standard dev
  for(Int_t j=0; j< n;j++){
      s+=pow(h.GetBin(j)-xmean,2);
  }
  s=sqrt(s/n);

  //now can get any moment
  if (m==1){
    return xmean;
  }else{
    for(Int_t j=0; j< n;j++){
      mom+=pow(h.GetBin(j)-xmean,m)/s;
    }
    return mom/n;
  }
}


/****************************************
Some filtering functions:
photFilter functions take last values in
bhwide output.
 ***************************************/

//to get events where only one has had 
//beamsstrahlung
bool beamsstrahlFilter(Double_t a, Double_t b){
  if(a>0.99 && b==1.0)return kFALSE;
  else if (b>0.99 && a==1.0)return kFALSE;
  else return kTRUE;
}

//to get events where no photon 
//is emitted
bool photFilter0(Double_t a, Double_t b){
  if( (a== 0) && (b == 0) ) return kFALSE;
  else return kTRUE;
}

//to get events where only one photon is 
//emitted in one direction
bool photFilter1(Double_t a, Double_t b){
  if( (a==0) && (b <= 1)) return kFALSE;
  else if ( (b==0) && (a <= 1)) return kFALSE;
  else return kTRUE;
}

//to get events where one photon is 
//emitted in each direction  
bool photFilter2(Double_t a, Double_t b){
  if( (a<=1) && (b<=1) ) return kFALSE;
  else return kTRUE;
}

//get rid of any events where theta is 
//outside range
bool angleRange(Double_t a, Double_t b,Double_t low,Double_t hi){
  if ((a<low) || (a>hi)) return kTRUE;
  else if((b<low) || (b>hi)) return kTRUE;
  else return  kFALSE;
}
bool phiRange(Double_t ph1, Double_t ph2,Double_t hi){
  if (abs(abs(ph1-ph2)-TMath::Pi()) <= hi) return kFALSE;
  else return kTRUE;
}

bool angleRange1(Double_t a){
  if ((a<0.099) || (a>0.301)) return kTRUE;
  else return  kFALSE;
}

//format numbers using 
//scientific notation, 
//3 decimal places
TString numFormat(Double_t a){
  Char_t buffer[1024];
  sprintf(buffer,"%.4E", a );
  return TString(buffer);
  //return buffer1;
}