#define FadcStuff_cxx
// The class definition in FadcStuff.h has been generated automatically
// by the ROOT utility TTree::MakeSelector(). This class is derived
// from the ROOT class TSelector. For more information on the TSelector
// framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual.

// The following methods are defined in this file:
//    Begin():        called every time a loop on the tree starts,
//                    a convenient place to create your histograms.
//    SlaveBegin():   called after Begin(), when on PROOF called only on the
//                    slave servers.
//    Process():      called for each event, in this function you decide what
//                    to read and fill your histograms.
//    SlaveTerminate: called at the end of the loop on the tree, when on PROOF
//                    called only on the slave servers.
//    Terminate():    called at the end of the loop on the tree,
//                    a convenient place to draw/fit your histograms.
//
// To use this file, try the following session on your Tree T:
//
// Root > T->Process("FadcStuff.C")
// Root > T->Process("FadcStuff.C","some options")
// Root > T->Process("FadcStuff.C+")
//

#include <iostream>

#include "FadcStuff.h"
#include "Filter.h"
#include <TStyle.h>
#include <TCanvas.h>

#ifndef SQR_FUNC
#define SQR_FUNC 1
double sqr(double x) { return x*x; }
#endif

Filter *dFilter=0;

FadcStuff::FadcStuff(TTree * /*tree*/)
{
   // Constructor

   fHist = 0;
   fHist2 = 0;
   fH1runs = 0;
   fH1events = 0;
   fTrace0 = 0;
   fTrace1 = 0;
   fTrace2 = 0;
   fTrace3 = 0;
   fTrace4 = 0;
   fTrace12 = 0;
   fTrace13 = 0;
   fTrace14 = 0;
   fTrace15 = 0;
   fRFclock = 0;
   fThist = 0;
   fTrist = 0;
   fTrig = 0;
   fTfilt = 0;

   // set the pedestal centroids and widths
   mean[0] = 412.5;
   stdev[0] = 34.9;
   mean[1] = 395.5;
   stdev[1] = 40.91;
   mean[2] = 370.7;
   stdev[2] = 44.32;
   mean[3] = 360.7;
   stdev[3] = 45.4;
   mean[4] = 363.9;
   stdev[4] = 30.9;
   mean[12] = 396.1;
   stdev[12] = 44.5;
   mean[13] = 487.6;
   stdev[13] = 34.6;
   mean[14] = 499.3;
   stdev[14] = 37.1;
}

FadcStuff::~FadcStuff()
{
   // Destructor
}

void FadcStuff::Begin(TTree * /*tree*/)
{
   // The Begin() function is called at the start of the query.
   // When running with PROOF Begin() is only called on the client.
   // The tree argument is deprecated (on PROOF 0 is passed).

   TString option = GetOption();
}

void FadcStuff::SlaveBegin(TTree * /*tree*/)
{
   // The SlaveBegin() function is called after the Begin() function.
   // When running with PROOF SlaveBegin() is called on each slave server.
   // The tree argument is deprecated (on PROOF 0 is passed).

   TString option = GetOption();
   TObjArray *options = option.Tokenize(",");
   TObjString *ostring;
   if (options->GetEntries() > 0) {
     ostring = (TObjString *)(*options)[0];
     xchan = ostring->GetString().Atoi();
   }
   else {
     xchan = 0;
   }
   if (options->GetEntries() > 1) {
     ostring = (TObjString *)(*options)[1];
     ychan = ostring->GetString().Atoi();
   }
   else {
     ychan = xchan;
   }
   delete options;

   if (xchan == ychan) {
     fHist = new TH1F("hist","pulse amplitude",512,0,2048);
     TString title("ADC channel ");
     fHist->GetXaxis()->SetTitle(title+=xchan);
     fHist->GetYaxis()->SetTitle("hits");
     fOutput->Add(fHist);
   }
   else {
     fHist2 = new TH2F("hist2","pulse correlations",400,0,1600,400,0,1600);
     TString xtitle("ADC channel ");
     fHist2->GetXaxis()->SetTitle(xtitle+=xchan);
     TString ytitle("ADC channel ");
     fHist2->GetYaxis()->SetTitle(ytitle+=ychan);
     fHist2->GetYaxis()->SetTitleOffset(1.8);
     fOutput->Add(fHist2);
   }
#if USE_EXTENDED_FADCSTUFF_TREE
   fH1runs = new TH1F("hruns","run number",200,0,200);
   fOutput->Add(fH1runs);
   fH1events = new TH1F("hevents","event number",500,0,5000000);
   fOutput->Add(fH1events);
   fRFclock = new TH1F("rf clock","rf clock",16000,-8000*0.06,8000*0.06);
   fOutput->Add(fRFclock);
   fThist = new TH1F("thist","hit-trig time",16000,-8000*0.06,8000*0.06);
   fOutput->Add(fThist);
   fTrist = new TH1F("trist","hit-rf time",16000,-8000*0.06,8000*0.06);
   fOutput->Add(fTrist);
   fTrig = new TH1F("trig","trig-rf time",16000,-8000*0.06,8000*0.06);
   fOutput->Add(fTrig);
   fTfilt = new TH1F("tfilt","hit-trig filtered time",8000,-100,100);
   fOutput->Add(fTfilt);
#endif
}

Bool_t FadcStuff::Process(Long64_t entry)
{
   // The Process() function is called for each entry in the tree (or possibly
   // keyed object in the case of PROOF) to be processed. The entry argument
   // specifies which entry in the currently loaded tree is to be processed.
   // It can be passed to either FadcStuff::GetEntry() or TBranch::GetEntry()
   // to read either all or the required parts of the data. When processing
   // keyed objects with PROOF, the object is already loaded and is available
   // via the fObject pointer.
   //
   // This function should contain the "body" of the analysis. It can contain
   // simple or elaborate selection criteria, run algorithms on the data
   // of the event and typically fill histograms.
   //
   // The processing can be stopped by calling Abort().
   //
   // Use fStatus to set the return value of TTree::Process().
   //
   // The return value is currently not used.

   b_amax->GetEntry(entry);

   Int_t tcut=0;
   Int_t aminCut=0;
   Int_t amaxCut=9999;
   Filter *aFilter=0;
   if (fInput != 0) {
      TIter iter(fInput);
      for (iter=iter.Begin(); iter!=iter.End(); ++iter) {
         if (strcmp((*iter)->GetName(),"tcut") == 0) {
            tcut = 1;
         }
         else if (strcmp((*iter)->GetName(),"amin") == 0) {
            aminCut = atoi((*iter)->GetTitle());
         }
         else if (strcmp((*iter)->GetName(),"amax") == 0) {
            amaxCut = atoi((*iter)->GetTitle());
         }
         else if (strcmp((*iter)->GetName(),"afilter") == 0) {
            aFilter = (Filter*)*iter;
            if (dFilter == 0) {
               dFilter = new Filter();
               dFilter->SetEdgeFilter();
            }
         }
         else if ((*iter)->InheritsFrom("FadcStuff::Cut2d")) {
            Cut2d *cut = (Cut2d *)(*iter);
            Double_t dist = cut->Dist(amax[cut->chan[0]],amax[cut->chan[1]]);
            if ((cut->reject && (dist < 1)) || (!cut->reject && (dist > 1))) {
               return kFALSE;
            }
         }
      }
   }

   if (amax[xchan] < aminCut || amax[xchan] > amaxCut) {
      return kFALSE;
   }

#if USE_EXTENDED_FADCSTUFF_TREE
   b_tcount->GetEntry(entry);
   b_tchan->GetEntry(entry);
   b_tchip->GetEntry(entry);
   b_ttime->GetEntry(entry);
   Double_t rftlast = -1;
   Double_t triglast = -1;
   Double_t thitlast = -1;
   for (int it=0; it<tcount; ++it) {
      if ((tchip[it] == 6) && ((tchan[it]== 4) || (tchan[it] == 5))) {
         Double_t tdcval = (unsigned short int)ttime[it];
         if (rftlast > -1) {
            Double_t tdiff = tdcval-rftlast;
            fRFclock->Fill(tdiff*0.060+0.03,1);
         }
         rftlast = tdcval;
      }
   }
   for (int it=0; it<tcount; ++it) {
      if ((tchip[it] == 6) && ((tchan[it]== 2) || (tchan[it] == 3))) {
         Double_t tdcval = (unsigned short int)ttime[it];
         if (rftlast > -1) {
            Double_t tdiff = tdcval-rftlast;
            fTrig->Fill(tdiff*0.060+0.03,1);
         }
         triglast = tdcval;
      }
   }
   for (int it=0; it<tcount; ++it) {
      switch (xchan) {
       case 0:
         if ((tchip[it] == 5) && ((tchan[it]== 0) || (tchan[it] == 1))) {
            Double_t tdcval = (unsigned short int)ttime[it];
            if (triglast > 0) {
               Double_t tdiff = tdcval-triglast;
               fThist->Fill(tdiff*0.060+0.03,1);
            }
            if (rftlast > 0) {
               Double_t tdiff = tdcval-rftlast;
               fTrist->Fill(tdiff*0.060+0.03,1);
            }
            thitlast = tdcval;
         }
         break;
       case 1:
         if ((tchip[it] == 5) && ((tchan[it]== 2) || (tchan[it] == 3))) {
            Double_t tdcval = (unsigned short int)ttime[it];
            if (triglast > 0) {
               Double_t tdiff = tdcval-triglast;
               fThist->Fill(tdiff*0.060+0.03,1);
            }
            if (rftlast > 0) {
               Double_t tdiff = tdcval-rftlast;
               fTrist->Fill(tdiff*0.060+0.03,1);
            }
            thitlast = tdcval;
         }
         break;
       case 2:
         if ((tchip[it] == 5) && ((tchan[it]== 4) || (tchan[it] == 5))) {
            Double_t tdcval = (unsigned short int)ttime[it];
            if (triglast > 0) {
               Double_t tdiff = tdcval-triglast;
               fThist->Fill(tdiff*0.060+0.03,1);
            }
            if (rftlast > 0) {
               Double_t tdiff = tdcval-rftlast;
               fTrist->Fill(tdiff*0.060+0.03,1);
            }
            thitlast = tdcval;
         }
         break;
       case 3:
         if ((tchip[it] == 5) && ((tchan[it]== 6) || (tchan[it] == 7))) {
            Double_t tdcval = (unsigned short int)ttime[it];
            if (triglast > 0) {
               Double_t tdiff = tdcval-triglast;
               fThist->Fill(tdiff*0.060+0.03,1);
            }
            if (rftlast > 0) {
               Double_t tdiff = tdcval-rftlast;
               fTrist->Fill(tdiff*0.060+0.03,1);
            }
            thitlast = tdcval;
         }
         break;
       case 4:
         if ((tchip[it] == 6) && ((tchan[it]== 0) || (tchan[it] == 1))) {
            Double_t tdcval = (unsigned short int)ttime[it];
            if (triglast > 0) {
               Double_t tdiff = tdcval-triglast;
               fThist->Fill(tdiff*0.060+0.03,1);
            }
            if (rftlast > 0) {
               Double_t tdiff = tdcval-rftlast;
               fTrist->Fill(tdiff*0.060+0.03,1);
            }
            thitlast = tdcval;
         }
         break;
       case 5:
         if ((tchip[it] == 5) && ((tchan[it]== 0) || (tchan[it] == 1))) {
            Double_t tdcval = (unsigned short int)ttime[it];
            if (triglast > 0) {
               Double_t tdiff = tdcval-triglast;
               fThist->Fill(tdiff*0.060+0.03,1);
            }
            if (rftlast > 0) {
               Double_t tdiff = tdcval-rftlast;
               fTrist->Fill(tdiff*0.060+0.03,1);
            }
            thitlast = tdcval;
         }
         break;
       case 12:
         if ((tchip[it] == 4) && ((tchan[it]== 2) || (tchan[it] == 3))) {
            Double_t tdcval = (unsigned short int)ttime[it];
            if (triglast > 0) {
               Double_t tdiff = tdcval-triglast;
               fThist->Fill(tdiff*0.060+0.03,1);
            }
            if (rftlast > 0) {
               Double_t tdiff = tdcval-rftlast;
               fTrist->Fill(tdiff*0.060+0.03,1);
            }
            thitlast = tdcval;
         }
         break;
       case 13:
         if ((tchip[it] == 4) && ((tchan[it]== 4) || (tchan[it] == 5))) {
            Double_t tdcval = (unsigned short int)ttime[it];
            if (triglast > 0) {
               Double_t tdiff = tdcval-triglast;
               fThist->Fill(tdiff*0.060+0.03,1);
            }
            if (rftlast > 0) {
               Double_t tdiff = tdcval-rftlast;
               fTrist->Fill(tdiff*0.060+0.03,1);
            }
            thitlast = tdcval;
         }
         break;
       case 14:
         if ((tchip[it] == 4) && ((tchan[it]== 6) || (tchan[it] == 7))) {
            Double_t tdcval = (unsigned short int)ttime[it];
            if (triglast > 0) {
               Double_t tdiff = tdcval-triglast;
               fThist->Fill(tdiff*0.060+0.03,1);
            }
            if (rftlast > 0) {
               Double_t tdiff = tdcval-rftlast;
               fTrist->Fill(tdiff*0.060+0.03,1);
            }
            thitlast = tdcval;
         }
         break;
      }
   }

   if (tcut && thitlast == -1) {
      return kFALSE;
   }

   b_runNo->GetEntry(entry);
   b_eventNo->GetEntry(entry);
   fH1runs->Fill(runNo,1);
   fH1events->Fill(eventNo,1);
#endif

   Double_t tzero[16];
   for (Int_t i=0; i<16; ++i) {
      tzero[i] = -9999;
   }
   if (aFilter) {
      b_trace0->GetEntry(entry);
      aFilter->Transform(trace0,125);
      Filter::PeakIterator it0 = aFilter->GetPeakIterator(amax[0]/2);
      amax[0] = (Short_t)(*(++it0)).mag/4;
      tzero[0] = (*it0).zero;
      b_trace1->GetEntry(entry);
      aFilter->Transform(trace1,125);
      Filter::PeakIterator it1 = aFilter->GetPeakIterator(amax[1]/2);
      amax[1] = (Short_t)(*(++it1)).mag/4;
      tzero[1] = (*it1).zero;
      b_trace2->GetEntry(entry);
      aFilter->Transform(trace2,125);
      Filter::PeakIterator it2 = aFilter->GetPeakIterator(amax[2]/2);
      amax[2] = (Short_t)(*(++it2)).mag/4;
      tzero[2] = (*it2).zero;
      b_trace3->GetEntry(entry);
      aFilter->Transform(trace3,125);
      Filter::PeakIterator it3 = aFilter->GetPeakIterator(amax[3]/2);
      amax[3] = (Short_t)(*(++it3)).mag/4;
      tzero[3] = (*it3).zero;
      b_trace4->GetEntry(entry);
      aFilter->Transform(trace4,125);
      Filter::PeakIterator it4 = aFilter->GetPeakIterator(amax[4]/2);
      amax[4] = (Short_t)(*(++it4)).mag/4;
      tzero[4] = (*it4).zero;
      b_trace12->GetEntry(entry);
      aFilter->Transform(trace12,125);
      Filter::PeakIterator it12 = aFilter->GetPeakIterator(amax[12]/2);
      amax[12] = (Short_t)(*(++it12)).mag/4;
      tzero[12] = (*it12).zero;
      b_trace13->GetEntry(entry);
      aFilter->Transform(trace13,125);
      Filter::PeakIterator it13 = aFilter->GetPeakIterator(amax[13]/2);
      amax[13] = (Short_t)(*(++it13)).mag/4;
      tzero[13] = (*it13).zero;
      b_trace14->GetEntry(entry);
      aFilter->Transform(trace14,125);
      Filter::PeakIterator it14 = aFilter->GetPeakIterator(amax[14]/2);
      amax[14] = (Short_t)(*(++it14)).mag/4;
      tzero[14] = (*it14).zero;
      b_trace15->GetEntry(entry);
      dFilter->Transform(trace15,125);
      Filter::PeakIterator it15 = dFilter->GetPeakIterator(amax[15]/2);
      amax[15] = (Short_t)(*(++it15)).mag/4;
      tzero[15] = (*it15).zero;

      // The following fix is introduced to iron out bias non-uniformities
      // in the value of zero returned by the default edge filter.  It will
      // only work for the default edge filter of length 12.
      Int_t tbin = int(tzero[15]+0.300);
      Double_t dt = tzero[15]-tbin;
      Double_t dtcor = (((((-6.251)*dt 
                            +1.634)*dt 
                            +2.808)*dt
                            +0.4308)*dt
                            +0.4728)*dt
                            -0.1469;
      tzero[15] = tbin+dtcor;
   }

   // The following fix is introduced to iron out bias non-uniformities
   // in the value of zero returned by the default pulse filter.  It will
   // only work for the default pulse filter of length 12.
   Int_t tbin = int(tzero[xchan]);
   Double_t dt = tzero[xchan]-tbin;
   Double_t dtcor = (((-0.2108)*dt 
                       +0.4711)*dt
                       +0.7397)*dt;
   tzero[xchan] = tbin+dtcor;

   if (xchan == ychan) {
      fHist->Fill(amax[xchan]+0.5,1);
      fTfilt->Fill(tzero[xchan],1);
   }
   else if ((sqr((amax[xchan]-mean[xchan])/stdev[xchan])+sqr((amax[ychan]-mean[ychan])/stdev[ychan])) > 9) {
      fHist2->Fill(amax[xchan]+0.5,amax[ychan]+0.5,1);
      fTfilt->Fill(tzero[xchan]-tzero[ychan],1);
   }

#if USE_EXTENDED_FADCSTUFF_TREE
   if (fTrace0 == 0) {
      fTrace0 = new TH1F(Form("trace0_%d_%d",runNo,eventNo),
                         Form("run %d event %d channel 0 trace",runNo,eventNo),
                         125,0,500.);
      fTrace1 = new TH1F(Form("trace1_%d_%d",runNo,eventNo),
                         Form("run %d event %d channel 1 trace",runNo,eventNo),
                         125,0,500.);
      fTrace2 = new TH1F(Form("trace2_%d_%d",runNo,eventNo),
                         Form("run %d event %d channel 2 trace",runNo,eventNo),
                         125,0,500.);
      fTrace3 = new TH1F(Form("trace3_%d_%d",runNo,eventNo),
                         Form("run %d event %d channel 3 trace",runNo,eventNo),
                         125,0,500.);
      fTrace4 = new TH1F(Form("trace4_%d_%d",runNo,eventNo),
                         Form("run %d event %d channel 4 trace",runNo,eventNo),
                         125,0,500.);
      fTrace12 = new TH1F(Form("trace12_%d_%d",runNo,eventNo),
                          Form("run %d event %d channel 12 trace",runNo,eventNo),
                          125,0,500.);
      fTrace13 = new TH1F(Form("trace13_%d_%d",runNo,eventNo),
                          Form("run %d event %d channel 13 trace",runNo,eventNo),
                          125,0,500.);
      fTrace14 = new TH1F(Form("trace14_%d_%d",runNo,eventNo),
                          Form("run %d event %d channel 14 trace",runNo,eventNo),
                          125,0,500.);
      fTrace15 = new TH1F(Form("trace15_%d_%d",runNo,eventNo),
                          Form("run %d event %d channel 15 trace",runNo,eventNo),
                          125,0,500.);
      b_trace0->GetEntry(entry);
      b_trace1->GetEntry(entry);
      b_trace2->GetEntry(entry);
      b_trace3->GetEntry(entry);
      b_trace4->GetEntry(entry);
      b_trace12->GetEntry(entry);
      b_trace13->GetEntry(entry);
      b_trace14->GetEntry(entry);
      b_trace15->GetEntry(entry);
      for (int i=0; i<125; ++i) {
         Double_t t=i*4+1;
         fTrace0->Fill(t,trace0[i]);
         fTrace1->Fill(t,trace1[i]);
         fTrace2->Fill(t,trace2[i]);
         fTrace3->Fill(t,trace3[i]);
         fTrace4->Fill(t,trace4[i]);
         fTrace12->Fill(t,trace12[i]);
         fTrace13->Fill(t,trace13[i]);
         fTrace14->Fill(t,trace14[i]);
         fTrace15->Fill(t,trace15[i]);
      }
      fOutput->Add(fTrace0);
      fOutput->Add(fTrace1);
      fOutput->Add(fTrace2);
      fOutput->Add(fTrace3);
      fOutput->Add(fTrace4);
      fOutput->Add(fTrace12);
      fOutput->Add(fTrace13);
      fOutput->Add(fTrace14);
      fOutput->Add(fTrace15);
   }
#endif
   return kTRUE;
}

void FadcStuff::SlaveTerminate()
{
   // The SlaveTerminate() function is called after all entries or objects
   // have been processed. When running with PROOF SlaveTerminate() is called
   // on each slave server.
}

void FadcStuff::Terminate()
{
   // The Terminate() function is the last function to be called during
   // a query. It always runs on the client, it can be used to present
   // the results graphically or save the results to file.

   if ((fHist = dynamic_cast<TH1F *>(fOutput->FindObject(Form("hist"))))) {
      TCanvas *c1 = new TCanvas("c1","FadcStuff canvas");
      fHist->Draw("h");
      c1->SetLogy(true);
      c1->cd();
      c1->Update();
   }
   else if ((fHist2 = dynamic_cast<TH2F *>(fOutput->FindObject(Form("hist2"))))) {
      TCanvas *c1 = new TCanvas("c1","FadcStuff canvas",80,10,500,500);
      c1->SetLeftMargin(0.15);
      fHist2->Draw("col");
      c1->SetLogz();
      c1->cd();
      c1->Update();
   }
   else {
      Warning("Terminate", "no hist object returned");
   }

   // Make a persistent copy of all output objects.

   TIter next(fOutput);
   TObject *obj;
   while ((obj = next())) {
      obj->Clone();
   }
}

void FadcStuff::Streamer(TBuffer &b)
{
   if (b.IsReading()) {
      b >> fHist >> fHist2 >> fH1runs >> fH1events;
      for (Int_t chan=0; chan < 16; ++chan) {
         b >> mean[chan] >> stdev[chan];
      }
      b >> xchan >> ychan;
   }
   else if (b.IsWriting()) {
      b << fHist << fHist2 << fH1runs << fH1events;
      for (Int_t chan=0; chan < 16; ++chan) {
         b << mean[chan] << stdev[chan];
      }
      b << xchan << ychan;
   }
}

FadcStuff::Cut2d::Cut2d()
{
   chan[0] = -1;
   chan[1] = -1;
   reject = 0;
}

FadcStuff::Cut2d::Cut2d(Int_t c0, Int_t c1,
                        Double_t x0, Double_t y0,
                        Double_t amajor, Double_t aminor,
                        Double_t theta, Int_t veto) 
 : EllipticalCut(x0,y0,amajor,aminor,theta)
{
   chan[0] = c0;
   chan[1] = c1;
   reject = veto;
}

void FadcStuff::Cut2d::Streamer(TBuffer &b)
{
   if (b.IsReading()) {
      EllipticalCut::Streamer(b);
      b >> chan[0] >> chan[1];
      b >> reject;
   }
   else if (b.IsWriting()) {
      EllipticalCut::Streamer(b);
      b << chan[0] << chan[1];
      b << reject;
   }

}