/////////////////////////////////////////////////////////////////////////////
//
// EEmcIsoCal
//
// Author: Jason C. Webb
//
// Calibration of an EEMC smd plane using isolated showers.
//
/////////////////////////////////////////////////////////////////////////////

#include "EEmcIsoSmdCal.h"

#include <iostream>

#include <TH1F.h>
#include <TH2F.h>
#include <TProfile.h>

#include "EEezTower.h"
#include "EEezStrip.h"
#include "EEezCluster.h"

#include "EEmcDb.h"
#include "StEEmcDbMaker/EEmcDbItem.h"
#include "StEEmcUtil/EEfeeRaw/EEname2Index.h"

#include "TRandom.h"

ClassImp(EEmcIsoSmdCal);

/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////

 EEmcIsoSmdCal::EEmcIsoSmdCal ( const Char_t *myName, const Char_t *myTitle )
  : TDirectory ( myName, myTitle ) 
{
  // Class constructor

  if ( gRandom == NULL ) gRandom = new TRandom();

}

/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////

 Int_t EEmcIsoSmdCal::Init()
{
  // Books histograms and performs general initialization
  cd();

  mADC[0] = new TH2F("mADC_U","ADC response of highest U-strip beneath isolated cluster/seed tower",1024,0.,1024,288,0.,288.);
  mADC[1] = new TH2F("mADC_V","ADC response of highest V-strip beneath isolated cluster/seed tower",1024,0.,1024,288,0.,288.);

  mMIP[0] = new TH2F("mMIP_U","N_{mip} of highest U-strip beneath isolated cluster/seed tower",400,0.,400.,288,0.,288.);
  mMIP[1] = new TH2F("mMIP_V","N_{mip} of highest V-strip beneath isolated cluster/seed tower",400,0.,400.,288,0.,288.);
  mMIP[2] = new TH2F("mMIP_Uh","N_{mip} of highest U-strip beneath isolated cluster/seed tower, #eta > 1.5",400,0.,400.,288,0.,288.);
  mMIP[3] = new TH2F("mMIP_Vh","N_{mip} of highest V-strip beneath isolated cluster/seed tower, #eta > 1.5",400,0.,400.,288,0.,288.);
  mMIP[4] = new TH2F("mMIP_Ul","N_{mip} of highest U-strip beneath isolated cluster/seed tower, #eta < 1.5",400,0.,400.,288,0.,288.);
  mMIP[5] = new TH2F("mMIP_Vl","N_{mip} of highest V-strip beneath isolated cluster/seed tower, #eta < 1.5",400,0.,400.,288,0.,288.);

  mETow[0] = new TH2F("mETow_U","Cluster energy used for strips",50,0.,10.,288,0.,288.);
  mETow[1] = new TH2F("mETow_V","Cluster energy used for strips",50,0.,10.,288,0.,288.);

  mFreqLarge[0] = new TH1F("mFreqLarge_U","Frequency of strip being classified as largest",288,0.,288.);
  mFreqLarge[1] = new TH1F("mFreqLarge_V","Frequency of strip being classified as largest",288,0.,288.);

  cd("..");

  return 1;
}

/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////

 void EEmcIsoSmdCal::Fill ( EEezCluster *cluster ) 
{
  // Given a cluster, fill histograms

  //-- Get the seed tower
  EEezTower *seed = cluster -> getSeedTower();

  Int_t nu = seed -> getNUStrips();
  Int_t nv = seed -> getNVStrips();
  Int_t etabin = seed -> getEtabin();

  //--
  //-- Must have at least 3 strips in each plane
  //--
  if ( nu < 3 || nv < 3 ) return;
  //-- ------------------------------------------------------------- CUT --/

  //--
  //-- Find the strip with the largest mip response
  //--
  EEezStrip *umax = seed -> getUStrip(0);
  EEezStrip *vmax = seed -> getVStrip(0);
  
  for ( Int_t i = 1; i < nu; i++ ) {

    EEezStrip *strip = seed -> getUStrip(i);
    if ( strip -> getNumMips() > umax -> getNumMips() ) umax = strip;

  }

  for ( Int_t i = 1; i < nv; i++ ) {

    EEezStrip *strip = seed -> getVStrip(i);
    if ( strip -> getNumMips() > vmax -> getNumMips() ) vmax = strip;

  }

  //--
  //-- A max strip was not found (not possible for edge strip to
  //-- be maximum since shape-limits apply.)
  //--
  if ( umax == seed -> getUStrip(0) || 
       vmax == seed -> getVStrip(0) ) return;


  //  umax -> print();
  //  vmax -> print();

  mADC[0] -> Fill ( (Float_t)umax -> getRawAdc(),  (Float_t)umax -> getIndex()+0.5 );
  mMIP[0] -> Fill ( (Float_t)umax -> getNumMips(), (Float_t)umax -> getIndex()+0.5 );
  mETow[0] -> Fill ( cluster -> getEnergy(), (Float_t)umax -> getIndex()+0.5 );
  mFreqLarge[0] -> Fill ( umax -> getIndex()+0.5 );
  if ( etabin > 5 ) {
    mMIP[2] -> Fill ( (Float_t)umax -> getNumMips(), (Float_t)umax -> getIndex()+0.5 );
  }
  else {
    mMIP[4] -> Fill ( (Float_t)umax -> getNumMips(), (Float_t)umax -> getIndex()+0.5 );
  }
  
  mADC[1] -> Fill ( (Float_t)vmax -> getRawAdc(),  (Float_t)vmax -> getIndex()+0.5 );
  mMIP[1] -> Fill ( (Float_t)vmax -> getNumMips(), (Float_t)vmax -> getIndex()+0.5 );
  mETow[1] -> Fill ( cluster -> getEnergy(), (Float_t)vmax -> getIndex()+0.5 );
  mFreqLarge[1] -> Fill ( vmax -> getIndex()+0.5 );
  if ( etabin > 5 ) {
    mMIP[3] -> Fill ( (Float_t)vmax -> getNumMips(), (Float_t)vmax -> getIndex()+0.5 );
  }
  else {
    mMIP[5] -> Fill ( (Float_t)vmax -> getNumMips(), (Float_t)vmax -> getIndex()+0.5 );
  }

}

/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////

 void EEmcIsoSmdCal::RemakeGains ( Int_t sector, EEmcDb *db ) 
{
  // For the specified sector (1..12) overwrite the SMD strip 
  // gains with the relative gains determined by this analysis.

  std::cout << "RemakeGains() called" << std::endl << std::flush;

  cd();
  
  mMIP_pfy[0] = mMIP[0] -> ProfileY();
  mMIP_pfy[1] = mMIP[1] -> ProfileY();

  mMIP[0] -> GetXaxis() -> SetRangeUser(1.,400.);
  mMIP[1] -> GetXaxis() -> SetRangeUser(1.,400.);

  Float_t umean = mMIP[0] -> GetMean();
  Float_t vmean = mMIP[1] -> GetMean();

  mMIP[0] -> GetXaxis() -> SetRange(0,0);
  mMIP[1] -> GetXaxis() -> SetRange(0,0);

  TString uname = (sector<10)?"0":""; uname += sector; uname += "U";
  TString vname = (sector<10)?"0":""; vname += sector; vname += "V";


  std::cout << "Booking some new histograms for remaking gains" << std::endl << std::flush;

  TH1F *mURelGain = new TH1F("mURelGain","Relative gains calculated",100,0.,4.);
  TH1F *mVRelGain = new TH1F("mVRelGain","Relative gains calculated",100,0.,4.);
  TH1F *mURelGainApplied = new TH1F("mURelGainApplied","Relative gains applied",100,0.,4.);
  TH1F *mVRelGainApplied = new TH1F("mVRelGainApplied","Relative gains applied",100,0.,4.);

  TH1F *mURelGainByTube = new TH1F("mURelGainByTube","Relative gains sorted by tube",1,0.,1.);
  TH1F *mVRelGainByTube = new TH1F("mVRelGainByTube","Relative gains sorted by tube",1,0.,1.);
  mURelGainByTube -> SetBit( TH1::kCanRebin );
  mVRelGainByTube -> SetBit( TH1::kCanRebin );


  //-- U plane
  std::cout << "Processing U plane" << std::endl << std::flush;

  for ( Int_t i = 0; i < 288; i++ ) {

    Int_t stripID = i + 1;
    TString strip = "";
    if ( stripID < 100 ) strip += "0";
    if ( stripID < 10  ) strip += "0";
    strip += stripID;    

    std::cout << uname+strip << std::endl << std::flush;
    Int_t uindex = EEname2Index( uname+strip );

    Float_t ugain = (db -> byIndex+uindex) -> gain;
 
    Float_t urel = mMIP_pfy[0] -> GetBinContent(stripID) / umean;

    Float_t urelerr = mMIP_pfy[0] -> GetBinError(stripID) / umean;

    Int_t nclust = (Int_t)mFreqLarge[0] -> GetBinContent(stripID);

    mURelGain -> Fill ( urel );

    //-- Sanity checks
    if ( urel == 0. ) continue;                // duh
    if ( urelerr / urel > 0.2   ) continue;    // 20% uncertainty
    if ( nclust < 25 ) continue;               // Sanity check on above
    
    mURelGainApplied -> Fill ( urel );
    mURelGainByTube -> Fill ( (db->byIndex+uindex)->tube, urel );

    ugain *= urel;
    (db -> byIndex+uindex) -> gain = ugain;
   
  }



  //-- V plane
  std::cout << "Processing V plane" << std::endl;

  for ( Int_t i = 0; i < 288; i++ ) {

    Int_t stripID = i + 1;
    TString strip = "";
    if ( stripID < 100 ) strip += "0";
    if ( stripID < 10  ) strip += "0";
    strip += stripID;

    std::cout << vname+strip << std::endl << std::flush;
    Int_t vindex = EEname2Index( vname+strip );

    Float_t vgain = (db -> byIndex+vindex) -> gain;

    Float_t vrel = mMIP_pfy[1] -> GetBinContent(stripID) / vmean;

    Float_t vrelerr = mMIP_pfy[1] -> GetBinError(stripID) / vmean;

    Int_t nclust = (Int_t)mFreqLarge[1] -> GetBinContent(stripID);

    mVRelGain -> Fill ( vrel );       

    //-- Sanity checks
    if ( vrel == 0. ) continue;                // duh
    if ( vrelerr / vrel > 0.2   ) continue;    // 20% uncertainty
    if ( nclust < 25 ) continue;               // Sanity check on above
 
    mVRelGainApplied -> Fill ( vrel );
    mVRelGainByTube -> Fill ( (db->byIndex+vindex)->tube, vrel );

    vgain *= vrel;
    (db -> byIndex+vindex) -> gain = vgain;

  }

  
  cd("..");

}

/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////

 void EEmcIsoSmdCal::RandomizeGains( Int_t sector, Float_t sigma, EEmcDb *db ) 
{
  // For the specified sector (0..11), introduce a gaussian
  // smearing with the specified width.  A floor of 0.2 is 
  // hardcoded.
  
  gRandom -> SetSeed( gRandom -> GetSeed() + sector );

  TString mySector;
  if ( sector+1 < 10 ) mySector += "0";
  mySector += sector+1;

  //-- Loop over all strips in both planes
  Char_t uv[] = {'U','V'};

  for ( Int_t ip = 0; ip < 2; ip++ ) {

    TString myPlane = mySector + uv[ip];

    for ( Int_t is = 1; is <= 288; is++ ) {

      TString myStrip = myPlane;
      if ( is < 100 ) myStrip += "0";
      if ( is <  10 ) myStrip += "0";
      myStrip += is;

      Int_t index = EEname2Index( myStrip );

      Float_t gain = db -> getByIndex( index ) -> gain;

    RAND:
      gain = gRandom -> Gaus ( gain, sigma );
      if ( gain < 0.2 ) goto RAND;

    }

  }

}

/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////