//###########################################
//
//
// Datei: double_decay.C
//
// Funktion: Es wird die Messreihe eines Doppelzerfalls analysiert
//
//
// Autor: Oliver Schitthelm
//
//
//###########################################

#include <iostream>
#include <fstream>      // for file handling
#include <string>       // for string handling

const int ndata_max     = 410; // maximum number of (x,y) pairs
                                //    allowed
int ndata         = 0;
double t_data[ndata_max + 1];
double y_data[ndata_max + 1];


// bekannt von früher: Einlesen einer Messreihe
int read_input_file(string filename)
{
  double t_file;    // x value found on file
  double y_file;    // y value found on file

  ndata = 0;

// define input stream "in", open input file

  ifstream in(filename.c_str()); // this can be understood only later...
  if (!in)
  {
    return 1;
  }

// read all data and store numbers in arrays

  while(1)
  {
    in  >> t_file >> y_file;
    if (in.eof()) break;
    ndata++;
    if(ndata <= ndata_max)
    {
      t_data[ndata] = t_file;
      y_data[ndata] = y_file;
    }
  }

  if (ndata>ndata_max)
  {
    return 2;
  }

// close input stream

   in.close();

   return 0;
}





void double_decay() {

  int result = read_input_file("data.txt");

  if(1 == result)
  {
    cout << " ERROR: Inputfile not available "
         << input_file << endl;
    return 1;
  }


  // die Werte sind in der Form (t,A) abgespeichert.
  // in t stehen die Zeiten, in y die Aktivitäten

  // visualisieren sie die Messdaten. Dies kann mittels
  // TGraph oder einem Histogramm geschehen
  TH1I* decay = new TH1I("decay", "Aktivitaet", ndata, 0, ndata );

  for ( int i = 0 ; i < ndata ; i++ ) {
    decay->Fill( t_data[i], y_data[i] );
  }

  decay->Draw();


  // betrachten sie anschließend das Programm decay_fit1.C aus
  // der Vorlesung. Versuchen sie an die Messdaten mittels
  // Chi^2 Methode eine Bestimmung der Zerfallszeit Tau durchzuführen

  double t_max = 400;
  double t_min = 0;
  double delta_t = (t_max-t_min)/ndata;

  double integral = decay->Integral();


  //c1->Update();  // Updates contents of (default) canvas


  //
  // calculate and display chi2 as a function of tau
  // -----------------------------------------------
  //

  int i_max = 10000;     // number of bins in chi2 histogram

  double tau_min = 0.;    // range of values to be
  double tau_max = 200.;  //   considered

  double delta_tau = (tau_max - tau_min)/i_max;

  TH1F* histogram1 = new TH1F("chi2","chi2",i_max,tau_min,tau_max);

  double chi2;

  double t, y, d, f;
  double tau;
  double factor;

  for (int i=1; i<=i_max; i++)      // loop over chi2 histogram
  {
    tau = tau_min + (i-0.5)*delta_tau;
    // would need to check for denominator = 0 !
    factor = integral / tau * delta_t
                      / (exp(-t_min/tau)-exp(-t_max/tau));
    chi2 = 0.;
    for (int n=1; n<=ndata; n++)    // loop over decay histogram
    {
      t = t_min+(n-0.5)*delta_t;
      y = decay->GetBinContent(n);
      d = 0;
      f = factor * exp(-t/tau);
      // need to check if d = 0 !
      chi2 += pow( (y-f)/d, 2 );
    }
    histogram1->Fill(tau,chi2);
  }

  //c2->cd();
  histogram1->Draw();
}