#include <vector>
using namespace std;

Double_t GetNofZ(Double_t zz, vector<Double_t> &z, vector<Double_t> &nz)
{
  Bool_t found = false;
  UInt_t i;
  for (i=0; i<z.size(); i++) {
    if (z[i]/10.0 >= zz) {
     found = true;
     break;
    }
  }

  if (!found)
    return -1.0;

  if (i == 0)
    return nz[0]*10.0*zz/z[0];

/*
cout << endl << "zz=" << zz << ", i=" << i << ", n(i)=" << nz[i] << ", n(i-1)=" << nz[i-1] << ", n(zz)=" << TMath::Abs(nz[i-1]+(nz[i]-nz[i-1])*(10.0*zz-z[i-1])/(z[i]-z[i-1])) << ", z(i)=" << z[i] << ", z(i-1)=" << z[i-1];
*/
  return TMath::Abs(nz[i-1]+(nz[i]-nz[i-1])*(10.0*zz-z[i-1])/(z[i]-z[i-1]));
}

Int_t FindBkgField(Double_t B, Double_t dB, vector<Double_t> &Bbkg)
{
  UInt_t i;
  Int_t result = -1;
  for (i=0; i<Bbkg.size(); i++) {
    if (TMath::Abs(Bbkg[i]-B) < dB) {
      break;
cout << endl << ">>i=" << i;
    }
  }

  if (i<Bbkg.size())
    result = (Int_t)i;

  return result;
}

void LSCO_test(TString &rge_fln, TString &out_fln, Double_t Bext, Double_t lambda, Double_t d, Double_t w, Double_t widthBkg)
{
  FILE *fp;
  char buffer[128];

  Double_t zz, nzz;
  vector<Double_t> z;
  vector<Double_t> nz;

  // read rge-file
  fp = fopen(rge_fln, "r");
  if (!fp) {
    cout << endl << "Couldn't open rge-file " << rge_fln.Data() << ", will quit.";
    cout << endl;
    return;
  }

  bool header = true;
  while (!feof(fp)) {
    fgets(buffer, sizeof(buffer), fp);
    if (header) {
      memset(buffer, 0, sizeof(buffer));
      header = false;
      continue;
    }
    zz = nzz = 0.0;
    sscanf(buffer, "%lf %lf", &zz, &nzz);
    z.push_back(zz);
    nz.push_back(nzz);
  }
  fclose(fp);

  // calculate pB - Meissner
  Double_t N = TMath::CosH(d/2.0/lambda);
  Double_t Bmin = Bext/N;
  Double_t BB, BBnext;
  Double_t dB = (Bext-Bmin)/1000.0;
  Double_t dz;
  Double_t zm, zp, zNext;
  vector<Double_t> B;
  vector<Double_t> pB;
  Double_t nn;
  for (UInt_t i=0; i<1000; i++) {
    BB = (Double_t)i/1000.0*(Bext-Bmin)+Bmin;
    BBnext = (Double_t)(i+1)/1000.0*(Bext-Bmin)+Bmin;
    B.push_back(BB);
    pB.push_back(0.0);
    // handle zm
    zm = d/2.0-lambda*TMath::ACosH(N*BB/Bext);
    zNext = d/2.0-lambda*TMath::ACosH(N*BBnext/Bext);
    dz = zNext-zm;
    nn = GetNofZ(zm, z, nz);
    if (nn != -1.0)
      pB[i] += nn*TMath::Abs(dz/dB);
/*
    else
      cout << endl << "zm: i=" << i;
*/

    // handle zp
    zp = d/2.0+lambda*TMath::ACosH(N*BB/Bext);
    zNext = d/2.0+lambda*TMath::ACosH(N*BBnext/Bext);
    dz = zNext-zp;
    nn = GetNofZ(zp, z, nz);
    if (nn != -1.0)
      pB[i] += nn*TMath::Abs(dz/dB);
/*
    else
      cout << endl << "zp: i=" << i;
*/
  }

  // normalize pB
  Double_t pBsum = 0.0;
  for (UInt_t i=0; i<B.size(); i++)
    pBsum += pB[i];
  pBsum *= dB;
  for (UInt_t i=0; i<B.size(); i++)
    pB[i] /= pBsum;

  if (w > 0.0) { // only handle bkg if it makes sense
    // calculate pB-Bkg
    vector<Double_t> Bbkg;
    vector<Double_t> pBbkg;
    Int_t range = (Int_t)(5.0*widthBkg/dB);
    for (UInt_t i=0; i<2*range; i++) {
      Bbkg.push_back(Bext+i*dB-5.0*widthBkg);
      pBbkg.push_back(TMath::Exp(-0.5*TMath::Power((Bext-Bbkg[i])/widthBkg, 2.0)));
    }
    // normalize pBbkg
    pBsum = 0.0;
    for (UInt_t i=0; i<pBbkg.size(); i++)
      pBsum += pBbkg[i];
    pBsum *= dB/w;
    for (UInt_t i=0; i<pBbkg.size(); i++)
      pBbkg[i] /= pBsum;
  }

  if (w > 0.0) {
    UInt_t idx;
    TGraph gg(B.size()+Bbkg.size()/2);
    for (UInt_t i=0; i<B.size(); i++) {
      idx = FindBkgField(B[i], dB, Bbkg);
      if (idx == -1) {
        gg.SetPoint(i, B[i], pB[i]);
//cout << endl << "idx == -1: i=" << i << ", B=" << B[i];
      } else {
        gg.SetPoint(i, B[i], pB[i]+pBbkg[idx]);
//cout << endl << "idx != -1: i=" << i << ", B=" << B[i];
      }
    }
    for (UInt_t i=Bbkg.size()/2; i<Bbkg.size(); i++) {
//cout << endl << ">=Bext: i=" << i+B.size()-Bbkg.size()/2 << ", B=" << Bbkg[i];
      gg.SetPoint(i+B.size()-Bbkg.size()/2, Bbkg[i], pBbkg[i]);
    }
  } else {
    TGraph gg(B.size());
    for (UInt_t i=0; i<B.size(); i++) {
      gg.SetPoint(i, B[i], pB[i]);
    }
  }

  gg.DrawClone("a*");

  pBsum = 0.0;
  Double_t valB, valpB;
  for (UInt_t i=0; i<gg.GetN(); i++) {
    gg.GetPoint(i, valB, valpB);
    pBsum += valpB;
  }
  pBsum *= dB;
//  cout << endl << "pBsum = " << pBsum << endl;

  // save pB
  fp = fopen(out_fln, "w");

  fprintf(fp, "# title: Meissner + background\n");
  fprintf(fp, "# parameter: B0 = %lf, lambda = %lf, d = %lf, w = %lf, widthBkg = %lf\n", Bext, lambda, d, w, widthBkg);
  fprintf(fp, "# -----------------------------------------------------\n");
  fprintf(fp, "# B, pB\n");
  for (UInt_t i=0; i<gg.GetN(); i++) {
    gg.GetPoint(i, valB, valpB);
    fprintf(fp, "%lf, %lf\n", valB, valpB/pBsum);
  }
  fprintf(fp, "# end of file\n");
  fclose(fp);

  // clean up
  z.clear();
  nz.clear();
  B.clear();
  pB.clear();
}