//
// fana.C - C++ utility classes for the fana pyroot macro
//

#include <fana.h>
#include <cmath>
#include <complex>
#include <TROOT.h>
#include <TTreeFormula.h>
#include <TH1D.h>

unsigned int Fana::SetPoints(TChain *ttagm, const char *varexp,
                                            const char *cond,
                                            int maxentries)
{
   /// load time series vector with values returned by expression varexp
   /// evaluated against the rows of tree ttagm that satisfy conditional
   /// expression cond.

   fVarexp = std::string(varexp);
   TTreeFormula varform("varexp", varexp, ttagm);
   const char defcond[] = "eventtime > 0";
   if (cond == 0 || strlen(cond) == 0)
      cond = defcond;
   TTreeFormula condform("condform", cond, ttagm);

   fTseries.clear();
   int maxentry = ttagm->GetEntries();
   double lastvar = 0;
   TBranch *branch = 0;
   int tree_size = 0;
   int i0 = 0;
   double var;
   for (int i=0; i < maxentry && i < maxentries; ++i) {
      ttagm->LoadTree(i);
      if (i == i0 + tree_size) {
         i0 = i;
#if SIMPLE_LEAF_FORMULA
         ttagm->SetBranchAddress(varexp, &var, &branch);
#else
         condform.UpdateFormulaLeaves();
         varform.UpdateFormulaLeaves();
#endif
         tree_size = ttagm->GetTree()->GetEntries();
      }
#if SIMPLE_LEAF_FORMULA
      branch->GetEntry(i-i0);
      if (true) {
         // assume condition true 
#else
      if (condform.EvalInstance() != 0) {
         double var = varform.EvalInstance();
#endif
         int repeats;
         if (var != lastvar) {
            fTseries.push_back(var);
            repeats = 0;
         }
         else {
            if (repeats > 10000) {
               printf("Error - TChain lost sync, too many repeated values!\n");
               return 0;
            }
            ++repeats;
         }
         lastvar = var;
      }
   }
   return fTseries.size();
}

TH1D *Fana::GetTransform(int n, double *ftmag) {
   /// Compute the Fourier transform of fTseries at regularly-spaced
   /// frequencies passed in the n elements of argument ftmag. Return
   /// the results as a frequency histogram containing the magnitudes
   /// of the Fourier coefficients as the bin contents.

   double f0 = ftmag[0];
   double f1 = ftmag[1] - ftmag[0];
   std::vector<std::complex<double> > phase;
   std::vector<std::complex<double> > efactor;
   const std::complex<double> i2pi(0, 2*M_PI);
   for (unsigned int i=0; i < fTseries.size(); ++i) {
      efactor.push_back(exp(i2pi * fmod(fTseries[i] * f0, 1)));
      phase.push_back(exp(i2pi * fmod(fTseries[i] * f1, 1)));
   }
   std::vector<std::complex<double> > ft;
   for (int p=0; p < n; ++p) {
      std::complex<double> sum;
      for (unsigned int i=0; i < fTseries.size(); ++i) {
         sum += efactor[i];
#if VERBOSE
         printf("p=%d,i=%d,efactor[i]=(%f,%f),sum=(%f,%f)\n",
                p,i,efactor[i].real(),efactor[i].imag(),sum.real(),sum.imag());
#endif
         efactor[i] *= phase[i];
      }
      ft.push_back(sum);
   }

   std::string title("Fourier transform of ");
   title += fVarexp;
   TH1D *fthist = new TH1D("ft", title.c_str(), n, f0, f0 + n * f1);
   for (int p=0; p < n; ++p)
      fthist->SetBinContent(p+1, abs(ft[p]));
   return fthist;
}

TH1I *Fana::Check_t0(TChain *ttagm, const char *cond) {
   TH1I *h = dynamic_cast<TH1I*>(gROOT->FindObject("checkt0"));
   if (h != 0)
      h->Delete();
   h = new TH1I("checkt0", "check of thit - ttrig (ns)", 2000, -500., 500.);
   h->SetBit(TH1::kCanRebin);

   ttagm->LoadTree(1);
   const char defcond[] = "eventtime > 0";
   if (cond == 0 || strlen(cond) == 0)
      cond = defcond;
   TTreeFormula condform("cond", cond, ttagm);

   double F1_clock_period = 32. / 152 * 128 / 232;
   double F1_rollover_period = 3744;
   int maxentry = ttagm->GetEntries();
   
   int row;
   int column;
   int f1trig;
   int rftick;
   int tdctick;
   double rftime;
   double tdctime;
   double trigtime;
   double eventtime;
   ttagm->SetBranchAddress("row", &row);
   ttagm->SetBranchAddress("column", &column);
   ttagm->SetBranchAddress("tick", &tdctick);
   ttagm->SetBranchAddress("rftick", &rftick);
   ttagm->SetBranchAddress("rftime", &rftime);
   ttagm->SetBranchAddress("f1trig", &f1trig);
   ttagm->SetBranchAddress("trigtime", &trigtime);
   ttagm->SetBranchAddress("eventtime", &eventtime);
   ttagm->SetBranchAddress("time", &tdctime);
   for (int i=1; i <= maxentry; ++i) {
      ttagm->GetEvent(i);
      if (condform.EvalInstance() == 0)
         continue;
      double ttrigger = fmod(eventtime, F1_rollover_period);
      double tphase = tdctick * F1_clock_period/2 - ttrigger;
      tphase += (tphase > -1000)? 0 : F1_rollover_period;
      h->Fill(tphase);
   }
   return h;
}