#include <fstream>
#include <iostream>
#include "TFile.h"
#include "TProfile.h"
#include "TCanvas.h"
#include "TH1.h"
#include "TVector.h"
#include "TObjArray.h"
#include <TLimit.h>
#include <TLimitDataSource.h>
#include <TConfidenceLevel.h>


using namespace std;


double CalculateLambdaAtExpectedCL(, TH1D* BackgroundEnergy, TH1D* SignalEnergy);

void CLc_level_searching(){

  double lumi = 500.0;
  double Pe = 0.0;
  double Pp = 0.0;
  char filename[256];
  char _operator[20];
  double _mass=102;
  double _spin=0.5;
  double _lambda =1580;
  double ExpCL = 0.9; //0.997
  
  
 

    //sprintf(_operator,"Vector");   
    //sprintf(_operator,"Scalar,s-channel");
    sprintf(_operator,"Axial-vector");  
    cout<<_operator<<endl;

    
    sprintf(filename,"Sensitivity_%s_Pe%.1f_Pp%.1f_Lumi%.0f_90percent_CL.root",_operator,Pe,Pp,lumi);
    TFile * f= new TFile(filename,"recreate");
    sprintf(filename,"Sensitivity %s",_operator);
     TH1D* Sensitivity = new TH1D(filename,filename,250,1,250);
     
     
    

for(;_mass<111; _mass++){

  int counter=0; 
    
    //sprintf(filename,"/afs/desy.de/group/flc/pool/achaus/WIMPs/data/Templates_Andrii/Lumi500_%s_full/1D_Data_Signal_Pe%.1f_Pp%.1f_Lumi%.0f.root",_operator,Pe,Pp,lumi);
    sprintf(filename,"%s_Pe%.1f_Pp%.1f_Lumi%.0f.root",_operator,Pe,Pp,lumi);
    TFile * infile_sig = new TFile(filename,"READ");
    infile_sig->cd();
    
    
    sprintf(filename,"M%.0f_S%.1f_%s_Lambda%.0f",_mass,_spin,_operator,_lambda); 
    TH1D* sh=(TH1D*)infile_sig->Get(filename);
    
       
    //sprintf(filename,"/afs/desy.de/group/flc/pool/achaus/WIMPs/data/Templates_Andrii/Data_Background_Pe%.1f_Pp%.1f_Lumi%.0f.root",Pe,Pp,lumi);
    //TFile * infile_bg = new TFile(filename,"READ");
    //infile_bg->cd();

    sprintf(filename,"Background_Pe%.1f_Pp%.1f_Lumi%.0f",Pe,Pp,lumi);
    TH1D* bh=(TH1D*)infile_sig->Get(filename);

   	
	cout<<endl <<" Mass = "<<_mass<<endl;	
	

	sh->Scale(pow(_lambda,4));


  double CLs = 0.0;
  double k = 1000.0;
  double k_last = 1000.0;
  double stepwidth = 100; 

  while(stepwidth > 1) //  normal 0.0001   
    {
      CLs = 0.0;
      while (CLs < (1-ExpCL))
	{ 
	  sh->Scale(1/pow(k,4)); 
	  
	 
	CLs=CalculateLambdaAtExpectedCL( bh, sh);
	 
	  sh->Scale(pow(k,4));
	  
	  
	  if (CLs < (1-ExpCL))
	    k_last = k;
	  
	  cout <<"k1 out "<<k<<endl;
	  k += stepwidth;
 
	    
	    counter++;
	    cout<<"Counter = "<<counter<<endl;
	}
      k = k_last;
      cout <<"k2 out "<<k<<endl;
      cout <<"stepwidth = "<<stepwidth<<endl;
      stepwidth *= 0.5;
    }
    cout <<"k3 out "<<k; 
   
	
    	
 	Sensitivity->Fill(_mass,k); 

	
	infile_sig->Close();
    	//infile_bg->Close();

	
	
	}


  f->cd();
  Sensitivity->Write();
  cout<<"END LOOP"<<endl;
		
  
 f->Close();
 cout<<"SAVE FILE"<<endl; 
}





// #######################################################################
double CalculateLambdaAtExpectedCL(TH1D* BackgroundEnergy, TH1D* SignalEnergy)
{
  double CLs = 0.0;

     	  TVectorD errorb(1);
	  TVectorD errors(1);
	  TObjArray* names = new TObjArray();
	  TObjString name1("bg uncertainty");
	  TObjString name2("sh uncertainty");	  
	  	  names->AddLast(&name1);
		  names->AddLast(&name2);
 	 
	
	  errorb[0]=0.003; // error source 1: 0.3%
	  errors[0]=0.0;   // error source 1: 0%
     	 // errorb[1]=0.0;   // error source 2: 0%
 	 // errors[1]=0.0;   // error source 2: 0%
 	  
	 
	  TLimitDataSource *infile = new TLimitDataSource(SignalEnergy,BackgroundEnergy,BackgroundEnergy,&errors,&errorb,names);
	  TConfidenceLevel *confidence = TLimit::ComputeLimit(infile,10000);


	  CLs = confidence->GetExpectedCLs_b(); 
	  
	  cout<<endl<<"CLs = "<<CLs<<endl;
	  
	  
	  delete confidence;
	  delete infile;
	  

  return CLs;
}