//Global variables [units: MeV, cm]

//beam energy [MeV]
double E0=10E3;

//resonance energy and width [MeV]
double ER=2E3;
double W=1;

//energy loss, MeV/cm
double alpha=10;

//int. probability per unit length at resonance peak [1/cm]
double lambdaPeak=0.1;
double Range=E0/alpha;

double stepMax=1;

//dummy process for transportation only, int. probability per unit length
double sigmaTransport(){
  return 1;
}

double sigma(double E){
  double num=lambdaPeak*W*W;
  double den=((E-ER)*(E-ER)+W*W);
  return num/den;
}

double sigmaMean(double Ei,double Ef){

  double dE=Ei-Ef;
  
  double a=(lambdaPeak*W)/dE;

  double b=TMath::ATan((Ei-ER)/W);
  double c=TMath::ATan((Ef-ER)/W);

  return a*(b-c);
}

double sigmaMax(double E){
  if (E<ER) return sigma(E);
  else{
    double Emax=stepMax*alpha;
    double xi=ER/(ER+Emax);
    if (E<Emax/xi){
      return lambdaPeak;
    }
    else{
      double Estar=(ER > xi*E ? ER : xi*E);
      return sigma(Estar);
    }
  }
}

double EatX(double x){
  return E0-alpha*x;
}

TRandom3 rrand(0);


double funTeo(double *x,double *par){
  double xx=x[0];
  double f1=exp(-lambdaPeak*(W/alpha)*TMath::ATan((E0-ER)/W));
  double f2=exp(lambdaPeak*(W/alpha)*TMath::ATan((E0-ER-alpha*xx)/W));
  double f3=(lambdaPeak)/(1+(E0-ER-alpha*xx)/W*(E0-ER-alpha*xx)/W);

  return f1*f2*f3;


}

void test(){
  const int Ntry=50000;


  TH1D *hSigmaEf=new TH1D("hSigmaEf","hSigmaEf",50000,0,Range);
  TH1D *hSigmaEmean=new TH1D("hSigmaEmean","hSigmaEmean",50000,0,Range);

  double x=0;
  double maxV=0;
  double sigmaF;
  double Ei,Ef;
  double prob,randTransport,randInteraction;

  cout<<"E0 is: "<<E0<<endl;
  cout<<"Range is: "<<Range<<endl;
  cout<<"Lambda at peak is: "<<lambdaPeak<<endl;
  for (int ii=0;ii<Ntry;ii++){
    x=0;
    while (x<Range){
      Ei=EatX(x);
      maxV=sigmaMax(Ei);
      randInteraction=rrand.Exp(1./maxV);
      //randTransport=rrand.Exp(1./sigmaTransport());
      // cout<<ii<<" : "<<x<<" "<<Ei<<" "<<maxV<<" "<<randInteraction<<" "<<randTransport<<endl;

      randTransport=stepMax;
      if (randTransport<randInteraction){
	x=x+randTransport;
	continue;
      }
      else{
	x=x+randInteraction;
	Ef=EatX(x);
	sigmaF=sigma(Ef);
	prob=sigmaF/maxV;
	//cout<<ii<<" :: "<<maxV<<" "<<sigmaF<<" "<<prob<<endl;
	if (rrand.Uniform()<prob){ //interaction
	  hSigmaEf->Fill(x);
	  break;
	}
      }
    }
    if (ii%1000==0)cout<<"Particle "<<ii<<" done: "<<x<<endl;
  }


  for (int ii=0;ii<Ntry;ii++){
    x=0;
    while (x<Range){
      Ei=EatX(x);
      maxV=sigmaMax(Ei);
      randInteraction=rrand.Exp(1./maxV);
      //randTransport=rrand.Exp(1./sigmaTransport());
      // cout<<ii<<" : "<<x<<" "<<Ei<<" "<<maxV<<" "<<randInteraction<<" "<<randTransport<<endl;

      randTransport=stepMax;
      if (randTransport<randInteraction){
	x=x+randTransport;
	continue;
      }
      else{
	x=x+randInteraction;
	Ef=EatX(x);
	sigmaF=sigmaMean(Ei,Ef);
	prob=sigmaF/maxV;
	//cout<<ii<<" :: "<<maxV<<" "<<sigmaF<<" "<<prob<<endl;
	if (rrand.Uniform()<prob){ //interaction
	  hSigmaEmean->Fill(x);
	  break;
	}
      }
    }
    if (ii%1000==0)cout<<"Particle "<<ii<<" done: "<<x<<endl;
  }

  TF1 *fTeo=new TF1("fTeo",funTeo,0,ER-W,ER+W);

  
  hSigmaEf->Scale(1./Ntry,"width");
  hSigmaEmean->Scale(1./Ntry,"width");

  TCanvas *c=new TCanvas("c","c");
  c->SetLogy();
  hSigmaEf->Draw();
  hSigmaEmean->SetLineColor(2);
  hSigmaEmean->Draw("SAMES");

  fTeo->SetNpx(1000000);
  fTeo->SetLineColor(3);
  fTeo->SetLineWidth(2);
  fTeo->Draw("SAME");
}