/// Class to build average temperature inside TPC
/// over a particular time period (or run).  Also
/// has facilities to make an animated gif of the
/// TPC temperature over a period of time.
/// Blair Jamieson (C) 2009
#include "ITPCTempAvg.hxx"
#include <algorithm>

/// Initialize TPC temperature averaging object using unix times
ITPCTempAvg::ITPCTempAvg( int maxsamples, long autmin,  long autmax){
  fTimeMin = autmin;
  fTimeMax = autmax;
  fMaxSamples = maxsamples;

  Init();
}
  
/// Initialize TPC temperature averaging object using year,
/// month(0-11), day (1-31), ...
ITPCTempAvg:: ITPCTempAvg(int maxsamples,
				int ayearmin, int amonthmin, int adaymin, int ahourmin, int aminmin, int asecmin,
				int ayearmax, int amonthmax, int adaymax, int ahourmax, int aminmax, int asecmax ){
  fMaxSamples = maxsamples;

  tm mt;  
  // use time structure from time.h
  // to build unix time from year,month, etc.
  mt.tm_year = ayearmin - 1900;
  mt.tm_mon  = amonthmin-1;
  mt.tm_mday = adaymin;
  mt.tm_hour = ahourmin;
  mt.tm_min  = aminmin;
  mt.tm_sec  = asecmin;
  fTimeMin = (long)mktime(&mt);

  mt.tm_year = ayearmax - 1900;
  mt.tm_mon  = amonthmax-1;
  mt.tm_mday = adaymax;
  mt.tm_hour = ahourmax;
  mt.tm_min  = aminmax;
  mt.tm_sec  = asecmax;
  fTimeMax = (long)mktime(&mt);

  std::cout<<"<ITPCTempAvg> reqesting "<<fMaxSamples<<" samples from "
	   <<ayearmin<<"/"<<amonthmin<<"/"<<adaymin<<" "
	   <<ahourmin<<":"<<aminmin<<":"<<asecmin
	   <<" to "
	   <<ayearmax<<"/"<<amonthmax<<"/"<<adaymax<<" "
	   <<ahourmax<<":"<<aminmax<<":"<<asecmax<<std::endl;
  Init();
}


/// Initialize temperature / pressure averaging object
void ITPCTempAvg::Init(){
  // check if gifsicle program exists
  if ( gSystem->Exec("gifsicle --version") == 0 ){
    fHaveGifSicle = true;
  } else {
    fHaveGifSicle = false;
  }

  // first calculate delta-time for maximum number of samples:
  long adt;
  if ( fMaxSamples >0 ){
    adt = (fTimeMax-fTimeMin)/fMaxSamples;
    if (adt<600) adt = 600; // readings are only every five minutes anyway so set minimum
  } else {
    // default to trying every 600 sec?
    adt = 600; 
  }

  
  std::cout<<"<ITPCTempAvg::Init> Look for temperature readings from time "
	   <<fTimeMin<<" to "<<fTimeMax<<std::endl;
  std::cout<<"<ITPCTempAvg::Init> MaxSamples = "<<fMaxSamples<<" (every "<<adt<<" seconds)"<<std::endl;

  aT = new ITPCTemperature( fTimeMin, fTimeMax );

  // create average/rms histograms
  char ahname[256];
  char ahtitle[256];
  sprintf(ahname,"havgt_%ld_%ld",fTimeMin, fTimeMax );
  sprintf(ahtitle,"TPC T Avg (degC) %ld to %ld", fTimeMin, fTimeMax );
  fhAvgT = (TH2F*)( aT->GetTH2F()->Clone( ahname ) );
  fhAvgT->SetNdivisions(0,"X");
  fhAvgT->SetNdivisions(0,"Y");
  fhAvgT->SetTickLength(0,"X");
  fhAvgT->SetTickLength(0,"Y");
  fhAvgT->SetName( ahname );
  fhAvgT->SetTitle( ahtitle );
  fhAvgT->Clear();

  sprintf(ahname,"hrmst_%ld_%ld",fTimeMin, fTimeMax );
  sprintf(ahtitle,"TPC T RMS (degC) %ld to %ld", fTimeMin, fTimeMax );
  fhRMST = (TH2F*)( aT->GetTH2F()->Clone( ahname ) );
  fhRMST->SetNdivisions(0,"X");
  fhRMST->SetNdivisions(0,"Y");
  fhRMST->SetTickLength(0,"X");
  fhRMST->SetTickLength(0,"Y");
  fhRMST->SetName( ahname );
  fhRMST->SetTitle( ahtitle );
  fhRMST->Clear();

  // setup variables to calculate mean and rms
  int inx = fhAvgT->GetNbinsX();
  int iny = fhAvgT->GetNbinsY();
  for (int ix=0; ix<=inx; ix++){
    asum.push_back ( vector<double>() );
    asum2.push_back ( vector<double>() );
    anum.push_back( vector<int>() );
    for (int iy=0; iy<=iny; iy++){
      asum[ix].push_back(0.);
      asum2[ix].push_back(0.);
      anum[ix].push_back(0);
    }
  }

  // setup counters
  fTMin = 273.;
  fTMax = -273.;
  fRMSTMin = 273.;
  fRMSTMax = -273.;
  fTAvg = 0.;
  fTRMS = 0.;
  fNavg=0;

  // setup for gifsicle    
  TCanvas*tc=new TCanvas("canvas","canvas",1400,700);
  char aname[256];

  long astart_time;
  long aend_time;
  aT->GetTimeInterval(astart_time, aend_time);
  int nquery = 0;
  do {
    if (nquery%1 == 0 ){
      std::cout<<"DB Query"<<nquery
	       <<" time = "<<aT->GetTimeStamp();
    }
    AddUpdate( false );

    if (fHaveGifSicle==true){
      tc->cd();
      tc->Clear();
      sprintf( aname, "ttpctempavg_%06d.gif",nquery);
      TH2F * htmp = aT->GetTH2F();
      htmp->SetNdivisions(0,"X");
      htmp->SetNdivisions(0,"Y");
      htmp->SetTickLength(0,"X");
      htmp->SetTickLength(0,"Y");
      htmp->SetMinimum( 12.0 );// std::max<float>(fTMin,fTAvg-30.0*fTRMS) );
      htmp->SetMaximum( 42.0 ); //std::min<float>(fTMax,fTAvg+30.0*fTRMS) );
      htmp->Draw("colz");
      tc->Print(aname);
    }

    aT->ReadNext( adt );
    nquery++;
    aT->GetTimeInterval(astart_time, aend_time);
  } while ( aend_time < fTimeMax );

  // add final point
  AddUpdate( true );
  fNTemps = nquery+1;

  if (fHaveGifSicle==true){
    tc->cd();
    tc->Clear();
    sprintf( aname, "ttpctempavg_%06d.gif",nquery);
    TH2F * htmp = aT->GetTH2F();
    htmp->SetNdivisions(0,"X");
    htmp->SetNdivisions(0,"Y");
    htmp->SetTickLength(0,"X");
    htmp->SetTickLength(0,"Y");
    htmp->SetMinimum( 12.0 );// std::max<float>(fTMin,fTAvg-30.0*fTRMS) );
    htmp->SetMaximum( 42.0 ); //std::min<float>(fTMax,fTAvg+30.0*fTRMS) );
    htmp->Draw("colz");
    tc->Print(aname);
  }

  return;
}


void ITPCTempAvg::AddUpdate( bool lastflag){

  // add to sum and sum^2 temperature from histogram bins
  TH2F * htmp = aT->GetTH2F();
  int inx = htmp->GetNbinsX();
  int iny = htmp->GetNbinsY();
  for (int ix=1; ix<=inx; ix++){
    for (int iy=1; iy<=iny; iy++){
      float curt = htmp->GetBinContent(ix,iy);
      if ( curt > 3. && curt < 100. ){ // only avg reasonable temperatures
	(asum[ix][iy]) += curt;
	(asum2[ix][iy]) += curt * curt;
	(anum[ix][iy]) ++;
      }
    }
  }
  // keep track of min/max temperatures
  float tmpmin;
  float tmpmax;
  aT->GetMinMaxT( tmpmin, tmpmax );
  if ( tmpmin < fTMin ) fTMin = tmpmin;
  if ( tmpmax > fTMax ) fTMax = tmpmax;

  // if last call, compute mean, sigma
  if ( lastflag == true ){
    for (int ix=1; ix<=inx; ix++){
      for (int iy=1; iy<=iny; iy++){
	double aTavg =  0.;
	if  ( anum[ix][iy] > 0 ) {
	  aTavg = (asum[ix][iy])/double(anum[ix][iy]);
	}
	fhAvgT->SetBinContent(ix,iy, aTavg );
	double aTrms = 0.;
	if ( anum[ix][iy] > 1 ){
	  aTrms = float( ( (asum2[ix][iy]) - 
			   double(anum[ix][iy])*double(aTavg)*double(aTavg) )/double( (anum[ix][iy])-1) );
	  //std::cout<<"ix="<<ix<<" iy="<<iy<<" Avg="<<aTavg<<" N="<<anum[ix][iy]<<" rms^2="<<aTrms<<std::endl;
	  if (aTrms<0.) aTrms = 0.;
	  else aTrms = sqrt(aTrms);
	} 
	fhRMST->SetBinContent(ix,iy,aTrms);
	if ( aTrms >0. && aTrms < fRMSTMin ) fRMSTMin = aTrms;
	if ( aTrms >0. && aTrms > fRMSTMax ) fRMSTMax = aTrms;
	
	// check if this bin contains a real measurement
	if ( strlen( aT->GetBinName(ix,iy) ) > 3 && aTavg > 5. ){
	  fNavg++;
	  fTAvg += aTavg;
	  fTRMS += aTrms;
	}
      }
    }

    if ( fNavg>0 ){
      fTAvg /= float(fNavg);
      fTRMS /= float(fNavg);
    }
    
    // Set min/max for histograms
    fhAvgT->SetMinimum( fTAvg-7.0 ); // std::max<float>(fTMin,fTAvg-30.0*fTRMS) );
    fhAvgT->SetMaximum( fTAvg+7.0 ); // std::min<float>(fTMax,fTAvg+30.0*fTRMS) );
    fhRMST->SetMinimum( fRMSTMin );
    fhRMST->SetMaximum( fRMSTMax );
  }

  return;
}

/// Method to create gifsicle of temperature readings used
void ITPCTempAvg::MakeGifSicle(){
  if ( fHaveGifSicle==false ){
    std::cout<<"<ITPCTempAvg::MakeGifSicles> WARNING: Cannot make gifsicles; first install gifsicle program"<<std::endl;
    std::cout<<"<ITPCTempAvg::MakeGifSicles> see http://www.lcdf.org/~eddietwo/gifsicle/"<<std::endl;
    return;
  }
  
  char cmdline[256];
  sprintf(cmdline,"gifsicle --delay=10 --loop=5 ttpctempavg_*.gif > ttpctemps_%ld_to_%ld.gif",
	  fTimeMin, fTimeMax );

  gSystem->Exec(cmdline);
  gSystem->Exec("rm -f ttpctempavg_*.gif");

  std::cout<<"<ITPCTempAvg::MakeGifSicle> gifsicle ttpctemps_" 
	   << fTimeMin << "_to_" << fTimeMax << ".gif was created. "<<std::endl;
}


/// Print average and rms temperature results to std::cout
void ITPCTempAvg::Print(){
  std::cout<<"TPC Volume from "
	   <<fTimeMin<<" to "
	   <<fTimeMax << std::endl;
  std::cout<<"Number of Temperature readings: "<<fNTemps<<std::endl;
  std::cout<<"  Avg+-RMS AllTemperatures = "<<fTAvg<<" +- "<<fTRMS<<" degC "<<std::endl;

  int inx = fhAvgT->GetNbinsX();
  int iny = fhAvgT->GetNbinsY();
  for (int ix=1; ix<=inx; ix++){
    for (int iy=1; iy<=iny; iy++){      
      // ignore empty bins
      if ( strlen( aT->GetBinName(ix,iy) ) < 3 ) continue;
      std::cout<<"  "<<ix<<" "<<iy<<"  Avg+-RMS "<<aT->GetBinName(ix,iy)<<" = "
	       <<fhAvgT->GetBinContent(ix,iy)<<" +- "
	       <<fhRMST->GetBinContent(ix,iy)<<" degC"<<std::endl;
    }
  }
  return;
}