// ---------------------------------------------------------------------------------
// -----                       PndSttHitProducerReal                           -----
// -----                        Author: S.Costanza                             -----
// ---------------------------------------------------------------------------------

#include <iostream>
#include <fstream>
using namespace std;

//ROOT headers
#include "TClonesArray.h"
#include "TH1F.h"
#include "TMath.h"
#include "TCanvas.h"
#include "TMatrixT.h"
#include "TFolder.h"
#include "TList.h"

//Fair heaeders
#include "FairRootManager.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"

//Stt headers
#include "PndSttRawHit.h"
#include "PndSttHitProducerReal.h"
#include "PndSttHit.h"
#include "PndSttTube.h"
#include "PndSttMapCreator.h"
#include "PndSttGeometryMap.h"

PndSttHitProducerReal::PndSttHitProducerReal()
{
  // init arrays
  fSttParameters = NULL;
  fSttTubeArray  = NULL;

  fPersistence = kTRUE;
  fevtn = 0;
  fVerbose = 0;

}


PndSttHitProducerReal::~PndSttHitProducerReal()
{
}


InitStatus PndSttHitProducerReal::Init()
{
  FairRootManager *ioman = FairRootManager::Instance();
  if ( ! ioman ) {
    cout << "-E- PndSttHitProducerReal::Init: "
	 << "RootManager not instantiated!" << endl;
    return kFATAL;
  }

  // ---> recalling STT map
  PndSttMapCreator *mapper = new PndSttMapCreator(fSttParameters);
  if (!mapper) {
    cout << "-E- PndSttHitProducerReal::Init: "
	 << "mapper not instantiated!" << endl;
    return kFATAL;
  }
  fSttTubeArray = mapper->FillTubeArray();

  fRawData = (TClonesArray*) ioman->GetObject("SttRawHit");
  if (!fRawData) {
    cout << "-E- PndSttHitProducerRealt::Init: "
	 << "No input raw array!" << endl;
    return kERROR;
  }

  // Create and register output array
  fHitArray = new TClonesArray("PndSttHit");
  ioman->Register("STTHit","STT",fHitArray, fPersistence);

  return kSUCCESS;

}


// ----- Retrieve geometry parameters ----------------------------------------
void PndSttHitProducerReal::SetParContainers() {
  FairRuntimeDb* rtdb = FairRunAna::Instance()->GetRuntimeDb();
  fSttParameters = (PndGeoSttPar*) rtdb->getContainer("PndGeoSttPar");
}


void PndSttHitProducerReal::Exec(Option_t *option)
{
  if(fVerbose==0 && fevtn%100==0) 
    cout << "Event Number "<<fevtn<<endl;
  else if(fVerbose >= 3) cout << "Event Number "<<fevtn<<endl;

  fevtn++;
  
  // Reset output array
  if ( ! fHitArray ) Fatal("Exec", "No HitArray");
  
  fHitArray->Clear();

  Int_t detID = 1;    // detectorID
  Int_t layer = 0;
  Int_t tubeID = 0;
  TVector3 pos, dpos; // position and error vectors
  Double_t pulset =-1;
  Double_t radius = 0;
  Double_t driftTime = 0;
  Double_t rawled = 0., rawzct = 0., rawcfd = 0.;
  Double_t depcharge = 0;
  // fast simulation 
  Int_t flag = 2;     // 0) standard curve, from simulation
                      // 1) Juelich exp curve from COSY-TOF (old) (2 bar, 80/20)
                      // 2) Juelich exp curve from COSY-TOF (Feb 2011) (1.25 bar, 80/20)
                      // 3) Juelich exp curve from prototype (Apr 2010)(2 bar, 90/10)

  PndSttRawHit *hit;
  for (Int_t i = 0; i<fRawData->GetEntriesFast(); i++) {
    hit = (PndSttRawHit*) fRawData->At(i);
    
    layer = hit->GetLayer();
    tubeID = hit->GetTubeID();
    PndSttTube* pSttTube = (PndSttTube*) fSttTubeArray->At(tubeID);
    pos = pSttTube->GetPosition();
    dpos.SetXYZ(1./TMath::Sqrt(12), 1./TMath::Sqrt(12), 150./TMath::Sqrt(12));
    if (fVerbose>2)
      cout << "Tube " << tubeID << ", pos = (" << pos.Y() << ", " << pos.Z() <<")" << endl;

    rawled = TMath::Abs(hit->GetRawLed());
    rawzct = TMath::Abs(hit->GetRawZct());
    rawcfd = TMath::Abs(hit->GetRawCfd());

    driftTime = rawled - fTInfo[tubeID][0];
    radius = Calib(fFlag, driftTime); // output in cm
    if (radius>0.5) radius = 0.5;
    if (radius<0) radius = 0.;

    // error calculation according to the curve chosen by flag
    Double_t closestDistanceError = GetError(radius, flag); 
    pulset = hit->GetRawAmpl();
    depcharge = hit->GetRawDe();

    // create hit
    AddHit(detID, tubeID, i, pos, dpos, pulset, radius, closestDistanceError, depcharge);
  }
}


// -----   Private method AddHit   --------------------------------------------
PndSttHit* PndSttHitProducerReal::AddHit(Int_t detID, Int_t tubeID, Int_t iPoint, TVector3& pos, TVector3& dpos, Double_t p, Double_t rsim, Double_t closestDistanceError, Double_t depcharge) {
  // see PndSttHit for hit description
  TClonesArray& clref = *fHitArray;
  Int_t size = clref.GetEntriesFast();
 
  PndSttHit *hitnew = new(clref[size]) PndSttHit(detID, tubeID, iPoint, pos, dpos, p, rsim, closestDistanceError, depcharge);
  return hitnew;

}


// ----- Read t0 offset values from file -------------------------------------
Bool_t PndSttHitProducerReal::ReadToffset(TString filename) 
{
  
  ifstream fIn;
  fIn.open(filename);

  if (fIn.is_open()) {
    if (!fIn.good()) {
      cout << "-E- Input file for t0 offset not found!!!" << endl;
      return kFALSE;
    }

    Int_t t1, t2, t3, t4;
    fTInfo.ResizeTo(fNTubes, 3);
    for (int i=1; i<fNTubes; i++) {
      for (int j=0; j<2; j++)
	fTInfo[i][j] = 0.;
      fIn >> t1 >> t2 >> t3 >> t4;
      fTInfo[i][0] = t2;   // t0
      fTInfo[i][1] = t3;   // deltat
      fTInfo[i][2] = t4;   // noise level
      if (fVerbose>1) 
	cout << "Tube " << i << ", t0 " << fTInfo[i][0] << " ns, deltat " << fTInfo[i][1] << endl;
    }
    fIn.close(); 
  }
  else {
    cout << "-E- Input file for t0 offset not found!!!" << endl;
    return kFALSE;
  }

  return kTRUE;

}


// ----- Read radius-drift time calibration parameters from file ---------------------------------
Bool_t PndSttHitProducerReal::ReadCalib(TString filename) 
{

  ifstream fIn;
  fIn.open(filename);

  if (fIn.is_open()) {
    if (!fIn.good()) {
      cout << "-E- Input file for t0 offset not found!!!" << endl;
      return kFALSE;
    }

   Double_t t0, t1, t2, t3, t4;
    fIn >> t0 >> t1 >> t2 >> t3 >> t4;
    fcalfit[0] = t0;  
    fcalfit[1] = t1;  
    fcalfit[2] = t2;  
    fcalfit[3] = t3;
    fcalfit[4] = t4;
    if (fVerbose>0) 
      cout << "Fit parameters of the r-t curve: " << fcalfit[0] 
	   << "\t" << fcalfit[1] << "\t" << fcalfit[2] << "\t" 
	   << fcalfit[3] << "\t" << fcalfit[4] << endl;

    fIn.close(); 
  }
  else {
    cout << "-E- Input file for t0 offset not found!!!" << endl;
    return kFALSE;
  }

  return kTRUE;

}


// ----- Retrieve calib parameters  ----------------------------------------------------------
Double_t PndSttHitProducerReal::Calib(Int_t flag, Double_t time) 
{

  if (flag==1) {

    // radius-drift time calibration from Garfield 
    // ArCO2(90/10) @ 1850 V [calculated by P.Wintz]
    const Double_t calfit[5] = {0.0034801, 0.0534479, -6.32247e-5, -1.11069e-6, 4.31749e-9};
    Double_t radius = calfit[0] + calfit[1]*time + calfit[2]*pow(time,2) 
      + calfit[3]*pow(time,3) + calfit[4]*pow(time,4);
    
    return radius*0.1; // cm
  }
  else if (flag==2) {

    // radius-drift time calibration from PndSttRawAna
    Double_t radius = fcalfit[0] + fcalfit[1]*time + fcalfit[2]*pow(time,2) 
      + fcalfit[3]*pow(time,3) + fcalfit[4]*pow(time,4);
    
    return radius; // cm
  }
}


// ----- Get Error on radius -------------------------------------------------
Double_t PndSttHitProducerReal::GetError(Double_t RadiusCm, Int_t rescurve) {

  Double_t resmic=-1;
  
  // simulation res curve used
  if (rescurve==0) {
    resmic =      +1.06966e+02  
      -4.03073e+03 *RadiusCm 
      +1.60851e+05 *pow(RadiusCm,2)  
      -2.87722e+06 *pow(RadiusCm,3)   
      +2.67581e+07 *pow(RadiusCm,4)   
      -1.43397e+08 *pow(RadiusCm,5) 
      +4.61046e+08 *pow(RadiusCm,6) 
      -8.79170e+08 *pow(RadiusCm,7) 
      +9.17095e+08 *pow(RadiusCm,8)  
      -4.03253e+08 *pow(RadiusCm,9); 
  }

  // Juelich exp curve from COSY-TOF (old) used
  else if (rescurve==1) {
    resmic =     +1.48048e+02
	-3.35951e+02*RadiusCm 
	-1.87575e+03*pow(RadiusCm,2)  
	+1.92910e+04*pow(RadiusCm,3)   
	-6.90036e+04*pow(RadiusCm,4)   
	+1.07960e+05*pow(RadiusCm,5) 
	-5.90064e+04*pow(RadiusCm,6);
  }
  
  // Juelich exp curve from COSY-TOF (Feb 2011) used
  else if (rescurve==2) {

    // data from COSY-TOF (Feb 2011)
    // the parametrization comes from mm vs mm:
    // => RadiusCm must be in mm
    RadiusCm *= 10.; // cm -> mm
    // and resmic will be given in mm ...

    // pol5 parametriz
    resmic = 0.188119 
      + 0.00211993 * RadiusCm
      + 0.00336004 * pow(RadiusCm, 2)
      - 0.0103979  * pow(RadiusCm, 3)
      + 0.0033461  * pow(RadiusCm, 4)
      -0.000315764 * pow(RadiusCm, 5);
    if (fVerbose>2)
      cout << "radius " << RadiusCm << ", resolution (micron) " << resmic << endl;
    // convert resmic to micron and RadiusCm to cm
    resmic *= 1000.;
    RadiusCm *= 0.1;
}
  
  // Juelich exp curve from prototype (Apr 2010) used
  else if (rescurve==3) {
    // RadiusCm needed in mm --> error in mm
    resmic = 4.521331e-01
      -2.087216e-01 *10.*RadiusCm
      +4.911102e-02 *pow(10.*RadiusCm,2)
      -3.934728e-03 *pow(10.*RadiusCm,3);
    resmic = resmic*1000.; // error in micron
  }
  // resmic in micron
  return resmic*0.0001; // --> resmic in cm
}


ClassImp(PndSttHitProducerReal)