//
// PICASSO Universite de Montreal Analysis class : PVar and FVar that went into the 2009 publication; no localisation (yet!)
//
// Sujeewa Kumaratunga
//


#include <iomanip>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <fstream>
#include "math.h"
#include <sstream>

#include <TCanvas.h>
#include <TGraph.h>
#include <TNamed.h>
#include <TMinuit.h>
#include <TTree.h>
#include <TFile.h>
#include <TObject.h>
#include <TMath.h>

#include "PUdeMAnalysis.hh"

using namespace std;

ClassImp(PUdeMAnalysis);

PUdeMAnalysis *global_excl_cheaT3;

PUdeMAnalysis::PUdeMAnalysis(){
  
  numchaN = 0;

}

PUdeMAnalysis::PUdeMAnalysis(Pio *sssdata, Int_t VmE, Int_t RuN, unsigned int MetA){

  sssdata = new Pio(VmE, RuN, MetA); 
  //TTree *sdtree = sdata->GetDataTree();
  //sdtree->SetCacheSize(0);
  numchaN = 0;
  Int_t temparraysize = Int_t(sssdata->GetData(0)->GetELen());  
  buf_raw = new Double_t[temparraysize];
  buf_hp = new Double_t[temparraysize];
  
}

PUdeMAnalysis::~PUdeMAnalysis(){
  delete [] buf_raw;
  delete [] buf_hp;
}




Int_t PUdeMAnalysis::udemdiscr(Int_t vmE, Int_t ruN, unsigned int metA,
			       Int_t starTcnt, Int_t enDcnt,
			       Int_t numDE,
			       Int_t deT1,Int_t deT2,Int_t deT3,Int_t deT4,
			       Int_t deT5,Int_t deT6,Int_t deT7)
{

  /*
    This is the main function in the UdeM Analysis class
  */

  //   initialize some variables
  t_abs=0.0; 
  run_id = 0; vme = 0; meta_id = 0;
  Tpcs = 0;
  toTpzo = 0;
  idDet = 0;
  activeMass = 0.;
  activeMasssigma = 0.;
  
  //temperarory variable declaration
  Double_t max_amplitude = 0.; 
  Double_t temPesig; 
  Double_t temPedelta; 
  Double_t preSSuSigma;
  Int_t olddetno = -99;
  Double_t progress = 0.;
  Int_t histcount = 0;
  Int_t alldetevts = 0, treefillcount =0;
  Double_t *r1ymean;
  Double_t *r2ymean;
  Double_t *r1r2yratio;
  Int_t det_no = 0;
  Int_t detind = 0;
  activeMass = 0.;
  activeMasssigma = 0.;


  char filname1 [] = "MontrealLoc";
  char filname3 [] = ".root";
  char filname [100];
  sprintf(filname,"%s%d%d%d%s",filname1,vmE,ruN,metA,filname3);
    
  TFile *sIGnal;
  sIGnal = new TFile (filname,"RECREATE");

  TTree *tree;
  tree = new TTree("MtreeVar","Mtuple tree Variables, Physics");

  tree->SetDirectory(sIGnal);

  MontrealLoc *stkphy = new MontrealLoc;
  //tree->SetCacheSize(0);

  tree->Branch("PhysVariables",&stkphy);

  Int_t mac =  vmE; 

  //read raw data file
  data = new Pio(vmE, ruN, metA);
  //TTree *sdtree = data->GetDataTree();
  //sdtree->SetCacheSize(0);


  //the PAnalysis class let's me override the regular .psf file with my own; the only change I do here is to use a different high pass filter of 18-300kHz
  //This class is used to store all the analysis settings used in PBubble and Pfft. All options that are available in the analysis of Picasso data should be controllable with the contents of this class. The default way of initializing the class is with the file 'settings/std_ana.psf'.
  PAnalysisSettings *pic_set = new PAnalysisSettings("std_ana_montreal.psf");

  //This class (The Picasso Analysis Base Class) is used to initialize the settings for Pfft and PBubble. Both of these classes need a valid pointer to this class in order to work correctly. 
  PAB *pab = new PAB(data->GetPicasso(), pic_set);

  //PBubble is really Queens's Bubble class, BUT, unless I call the IsBubble memeber function from the PBubble class, I will not infact get a Queens bubble. Within this class they do a high pass filter, so I stole this class from them, just so I can modify my .psf file settings and apply my own high pass filter.
  PBubble *bb = new PBubble(pab);

  //If you have a Picasso .root file generated with Prootify you can access the data with the PData class
  PData *dd = 0;
  //dd->SetCacheSize(0);

  Long64_t e_cnt = data->GetEntries();
  cout << " There are " << e_cnt << " events in Fab file." << endl;

  sIGnal->cd();

  run_id=ruN;
  vme=vmE;
  meta_id = Int_t(metA);

  //get some run info
  //run start time, length, etc...
  Double_t ruNStart = data->GetRun()->GetStartDateUnix();
  Double_t ruNLeN = data->GetRun()->GetRunLength();
  
  //optional start and end event indices: allows you control over a small event range testing
  Int_t StarTcnt = 0; Int_t EnDcnt = e_cnt; 
  if(starTcnt != -1 )
    {
      StarTcnt = starTcnt;
    }
  if(enDcnt != -1)
    {
      EnDcnt = enDcnt;
    } 

		
  for (Long64_t ient = StarTcnt ; ient < EnDcnt ; ient++) {
    if ((Double_t)100 * ient/(Double_t)e_cnt - progress >= 10.){
      progress += 10.;
      cout << " Progress: " << setprecision(3) << progress << "%" << endl;
    }
    
    //declare some variables for each event
    Double_t *max_amplitude_ch, *max_amplitude_hp;
    Double_t *PowerDiffInt;
    Double_t ratE;

    bb->Analyse(data->GetEvent(ient));

    dd = data->GetEvent(ient); 
    if (data->GetEvent(ient) == 0) continue;

    //optional detector numbers: you can pick which detectors you want to analyze and you can pick up to 7 such detectors, or else you can just analyze all detectors with the first argument of 32.
    if( numDE == 32 || 
	dd->GetProductionNumber() == deT1 || dd->GetProductionNumber() == deT2 ||
	dd->GetProductionNumber() == deT3 || dd->GetProductionNumber() == deT4 ||
	dd->GetProductionNumber() == deT5 || dd->GetProductionNumber() == deT6 ||
	dd->GetProductionNumber() == deT7 
	)
      {
	alldetevts++;
	 
	//define some variables
	Double_t temPe; Double_t temPesig; Double_t temPedelta; Double_t preSSu; Double_t preSSuSigma;


	max_amplitude_ch= new Double_t[9];
	max_amplitude_hp= new Double_t[9];
	PowerDiffInt = new Double_t[9];
	 
	for(Int_t ijj=0; ijj<9; ijj++) 
	  {
	    PowerDiffInt[ijj]=-9.9;
	    max_amplitude_ch[ijj]=-9999.9;
	    max_amplitude_hp[ijj]=-9999.9;
	    buf_raw[ijj]=-10000.;
	    buf_hp[ijj]=-10000.;
	  }
	 
	max_amplitude=0.0;
	 
	//now get some event specific information
	det_no = dd->GetProductionNumber();
	detind = data->GetPAB()->GetIndexByProdNo(det_no);
	ratE = data->GetPAB()->GetSamplingFrequency()*1000.;// should be 400000.;
	t_abs = Double_t(dd->GetTrigTime(0) - (data->GetRun()->GetEvOffset())/(data->GetPAB()->GetSamplingFrequency()*1000.)); //trigger time is constant for all pzo, so just call pzo 0. 

	//assign values if we are in a new detector
	if (det_no != olddetno)
	  {
	    //T,P etc
	    Tpcs = data->GetPAB()->GetTPCSByProdNo(det_no);
	    data->GetDetAvgTemp(Tpcs, temPe, temPesig, temPedelta, preSSu, preSSuSigma);
	    //number of piezos
	    toTpzo = dd->GetChannels();

	    //get birth certificates for active mass etc
	    PBirthCertificateList *bc = new PBirthCertificateList;
	    idDet = bc->GetIndexByNumber(det_no);
	    
	    activeMass = bc->GetActiveMass(idDet);
	    activeMasssigma = bc->GetActiveMassError(idDet);

	    delete bc; //bc = NULL;

	  }
	olddetno = det_no;
	 
	//define a few variables for FVar
	r1ymean = new Double_t[9];
	r2ymean = new Double_t[9];
	r1r2yratio = new Double_t[9];
	 
	for(Int_t ijj=0; ijj<9; ijj++)
	  {
	    r1ymean[ijj] = -9.9;
	    r2ymean[ijj] = -9.9;
	    r1r2yratio[ijj] = -9.9;
	  }
	 
	for (Int_t k=0;k<dd->GetChannels();k++) { //opened piezos loop
	  if (dd->GetData(k) == 0) continue;
	  histcount++;

	  numchaN = dd->GetData(k)->GetNbinsX(); //number of channels; ie those little 4096 dots that make up the signal, should be 4096 (old daq) or 8192 (new daq)



	  max_amplitude = 0.;max_amplitude_ch[k]=0.0; max_amplitude_hp[k]=0.0;

	  PowerDiffInt[k] = 0.;
	  /*
	    we make a loop over all the samples (the amp dots that make the signal) . 
	  */
	  for (Int_t chanj=0; chanj<numchaN; chanj++)
	    {
	      buf_raw[chanj] = dd->GetData(k)->GetBinContent(chanj+1);
	      //buffer[chanj] = buf_raw[chanj];
	
	      if (fabs(buf_raw[chanj])>max_amplitude_ch[k]) { 
		max_amplitude_ch[k]=fabs(buf_raw[chanj]);
	      }

	      //now let's get the high-pass-filtered signal:
	      buf_hp[chanj] = bb->GetPSignal(k)->GetBinContent(chanj+1);
	      if (fabs(buf_hp[chanj])>max_amplitude_hp[k]) { 
		max_amplitude_hp[k]=fabs(buf_hp[chanj]);
	      }

	    } //closed the chanj loop 
	  
	  if (max_amplitude_ch[k]>max_amplitude) {
	    max_amplitude=max_amplitude_ch[k];
	  }	
	 	
	  //call PUdeMPVar class
	  PUdeMPVar *pudempvar = 0;
	  pudempvar = new PUdeMPVar(buf_hp, numchaN);
	  Int_t evofftopvar = Int_t((data->GetRun()->GetEvOffset()));
	  Double_t PowerDifffromPVar = 0.;
	  pudempvar->udempvar(evofftopvar, PowerDifffromPVar);
	  PowerDiffInt[k] = PowerDifffromPVar;
	  
	  //delete PUdeMPVar to free up memory
	  pudempvar->Clear();
	  delete pudempvar;  pudempvar = NULL;


	  //FVar calculation begin
	  //CAUTION!!! this is not what went into our publication.
	  //the FVar here is FFT after high pass, the FVar in the 2009 publication was FFT w/o high pass
	  //For this version, I am going to leave it simple and leave the high passed FFT, to avoid recasting PBubble
	  //Also I leave it here because we have not studied the difference between high pass and no high pass FFT.
	  TH1D *hUdeMfft = ((TH1D*) (bb->GetFFT(k)));

	  //call PUdeMFVar class
	  PUdeMFVar *pudemfvar = 0;
	  pudemfvar = new PUdeMFVar(hUdeMfft);
	  Double_t r1ymeanfromFV=0., r2ymeanfromFV=0., r1r2yratiofromFV=0.;
	  pudemfvar->udemfvar(r1ymeanfromFV, r2ymeanfromFV,  r1r2yratiofromFV);

	  r1ymean[k] = r1ymeanfromFV;
	  r2ymean[k] = r2ymeanfromFV;
	  r1r2yratio[k] = r1r2yratiofromFV;

	  //delete PUdeMFVar to free up memory
	  pudemfvar->Clear();
	  delete pudemfvar;  pudemfvar = NULL;


	  hUdeMfft->Reset();
	  //FVar Calculation end

	}//closed piezo loop, ie k loop

	//fill the tree; for all variables relating to localisation, I am going to fill a dummy of -100.0
	//I am keeping these leaves in the tree for posterity.

	stkphy->runindex = run_id; 
	stkphy->evindex = ient;
	stkphy->totevents = alldetevts;

	stkphy->IsMLocal = 100;
	;
	stkphy->udemposition[0] = -100.;
	stkphy->udemposition[1] = -100.;
	stkphy->udemposition[2] = -100.;
	      
	stkphy->udemsigmaplus[0] = -100.; 
	stkphy->udemsigmaplus[1] = -100.; 
	stkphy->udemsigmaplus[2] = -100.; 
	stkphy->udemsigmaminus[0] = -100.; 
	stkphy->udemsigmaminus[1] = -100.; 
	stkphy->udemsigmaminus[2] = -100.;
 
	for (int ijkk=0; ijkk<9; ijkk++) {
	  stkphy->udemtzero[ijkk] = -100.;
	  stkphy->udemsigmatzero[ijkk] = -100.;

	  stkphy->MaxAmplitudeRaw[ijkk] = max_amplitude_ch[ijkk];
	  stkphy->MaxAmplitudeHPFilt[ijkk] = max_amplitude_hp[ijkk];
	  
	  stkphy->PowerCutG[ijkk] = PowerDiffInt[ijkk];

	  for(Int_t xyzz=0; xyzz<5; xyzz++)
	    {
	      stkphy->WhyIsLocalFail[5*ijkk+xyzz] = -100.;
	    }
	}
	

	stkphy->velocity = -100.; 
	stkphy->sigmavelocity = -100.; 
	stkphy->chisquared = -100.; 
	stkphy->nuctime = -100.; 
	stkphy->numpzos = 100; 

	//does the event pass Montreal Bubble cut?
	Int_t isBuB = 0;
	if(bb->IsBubbleEvent()){ isBuB = 1;}
	stkphy->MontBubble = isBuB;

	for (int ijkk=0; ijkk<9; ijkk++) {
	  stkphy->goodpzos[ijkk] = -1;
	}
		
	stkphy->detname = dd->GetDetName();
	stkphy->prodnum = dd->GetProductionNumber();

	stkphy->runtype = data->GetRun()->GetRunType();
	stkphy->tabsol = t_abs;
	
	stkphy->ActiveMasS = activeMass; 
	stkphy->ActiveMasSSigma = activeMasssigma;
	stkphy->RunStarTTime = ruNStart; 
	stkphy->RunLengtH = ruNLeN; 
	stkphy->TotPzos = toTpzo; 
	stkphy->TpcsNum = Tpcs; 
	stkphy->TemperaturE = temPe; 
	stkphy->TemperaturESigma = temPesig; 
	stkphy->PressurE = preSSu;
	stkphy->PressurESigma = preSSuSigma; 


	for (int ijkk=0; ijkk<9; ijkk++) {
	  stkphy->R1Ymean[ijkk] = r1ymean[ijkk];
	  stkphy->R2Ymean[ijkk] = r2ymean[ijkk];
	  stkphy->R1YOverR2YRatio[ijkk] = r1r2yratio[ijkk];
	}

	treefillcount++;
   
	tree->Fill();

	delete [] r1ymean; r1ymean =  NULL;
	delete [] r2ymean; r2ymean = NULL;
	delete [] r1r2yratio; r1r2yratio = NULL;	

      } //close det loop
     
     
    delete [] max_amplitude_ch; max_amplitude_ch = NULL;
    delete [] max_amplitude_hp; max_amplitude_hp = NULL;
    delete [] PowerDiffInt; PowerDiffInt = NULL;
     
    dd->Clear(); 
    bb->Clear(); 
    stkphy->Clear(); 
		
    //} //close detector number loop, i.e.GetProductionNumber.
	
  } //close ient loop for events
  
  cout<<"there were "<<alldetevts<<" events and "<<"\t"<<treefillcount<<" treefill counts"<<endl;

  sIGnal->cd();
  tree->SetDirectory(sIGnal);
  sIGnal->Write();
  sIGnal->Close();

  return 0;
}