#include <getopt.h>
#include <fstream>
#include <sstream>
#include <string>
#include "Examples/Tutorial/LineShapes/TwoPoles.hh"
#include "qft++/matrix/KpoleMatrix.hh"
#include "qft++/matrix/TMatrix.hh"
#include "TFile.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TMath.h"


#include "ErrLogger/ErrLogger.hh"

TwoPoles::TwoPoles(double MassRes1, double Width1, double MassRes2, double Width2, double deltaMass) :
  _theTFile(0)
{

  if (MassRes1 < MassRes2){
    _massLow = MassRes1;
    _widthLow = Width1;
    _massHigh = MassRes2;
    _widthHigh = Width2;
    
  }
  else{
    _massLow = MassRes2;
    _widthLow = Width2;
    _massHigh = MassRes1;
    _widthHigh = Width1;
    
  }

  _theTFile=new TFile("TwoPoles.root","recreate");
 
  _massShapeHist= new TH1F("_massShapeHist", "BreitWigner",301, _massLow-deltaMass, _massHigh+deltaMass);
  _massShapeKmatrHist= new TH1F("_massShapeKmatrHist", "mass shape Kmatrix",301, _massLow-deltaMass, _massHigh+deltaMass);
  _massShapeKmatrHistRel= new TH1F("_massShapeKmatrHistRel", "mass shape rel. Kmatrix",301, _massLow-deltaMass, _massHigh+deltaMass);

  _argandHist=new TH2F("_argandHist","Argand plot",301, -1., 1., 301, 0., 1.3);
  _argandHist->SetXTitle("Re(Bw)");
  _argandHist->SetYTitle("Im(Bw)");
  _argandHist->SetMarkerStyle(6);
  
  _argandKmatrHist=new TH2F("_argandKmatrHist","Argand plot K matrix",301, -1., 1., 301, 0., 1.3);
  _argandKmatrHist->SetXTitle("Re(Bw)");
  _argandKmatrHist->SetYTitle("Im(Bw)");
  _argandKmatrHist->SetMarkerStyle(6);

  _argandKmatrHistRel=new TH2F("_argandKmatrHistRel","Argand plot rel. K matrix",301, -1., 1., 301, 0., 1.3);
  _argandKmatrHistRel->SetXTitle("Re(Bw)");
  _argandKmatrHistRel->SetYTitle("Im(Bw)");
  _argandKmatrHistRel->SetMarkerStyle(6);

  
  //initialize here the phase shift histogramms _phaseHistMap
  
  _phaseHist=new TH2F("_phaseHist", "phase shift",301, _massLow-deltaMass, _massHigh+deltaMass, 301, 0., 3.1415);
  _phaseHist->SetXTitle("mass");
  _phaseHist->SetYTitle("phase #delta(#phi)");
  _phaseHist->SetMarkerStyle(6);
  
  _phaseKmatrHist=new TH2F("_phaseKmatrHist", "phase shift Kmatrix",301, _massLow-deltaMass, _massHigh+deltaMass, 301, 0., 3.1415);
  _phaseKmatrHist->SetXTitle("mass");
  _phaseKmatrHist->SetYTitle("phase #delta(#phi)");
  _phaseKmatrHist->SetMarkerStyle(6);
  
  _phaseKmatrHistRel=new TH2F("_phaseKmatrHistRel", "phase shift rel Kmatrix",301, _massLow-deltaMass, _massHigh+deltaMass, 301, 0., 3.1415);
  _phaseKmatrHistRel->SetXTitle("mass");
  _phaseKmatrHistRel->SetYTitle("phase #delta(#phi)");
  _phaseKmatrHistRel->SetMarkerStyle(6);
  // break up momentum q0

  double piMass=0.135;
//   double rhoLow=sqrt( ( _massLow*_massLow-(piMass+piMass)*(piMass+piMass) ) 
// 		    * ( _massLow*_massLow - (piMass-piMass) * (piMass-piMass) ) ) / (_massLow);
  
//   double rhoHigh=sqrt( ( _massHigh*_massHigh-(piMass+piMass)*(piMass+piMass) ) 
// 		     * ( _massHigh*_massHigh - (piMass-piMass) * (piMass-piMass) ) ) / (_massHigh);

  double rhoLow=sqrt( (1-(2.*piMass/_massLow)*(2.*piMass/_massLow)) ); 
  
  double rhoHigh=sqrt( (1-(2.*piMass/_massHigh)*(2.*piMass/_massHigh)) ); 

  double gLow=sqrt(_massLow*_widthLow/rhoLow);
  double gHigh=sqrt(_massHigh*_widthHigh/rhoHigh);  

  pair <double, double> theMassDecPair(piMass, piMass);
  vector <pair <double, double> > theMassDecPairVec;
  theMassDecPairVec.push_back(theMassDecPair);

  vector<double> low_gis;
  low_gis.push_back(gLow);
  KpoleMatrix kPoleLow(low_gis, _massLow, theMassDecPairVec);

  vector<double> high_gis;
  high_gis.push_back(gHigh);
  KpoleMatrix kPoleHigh(high_gis, _massHigh, theMassDecPairVec);

  vector<KpoleMatrix> kPoles;
  kPoles.push_back(kPoleLow);
  kPoles.push_back(kPoleHigh);

  TMatrix theTMatrix(kPoles);
  
  double stepSize=( (_massHigh+deltaMass)-(_massLow-deltaMass) )/300.;  




  for (double massIt=_massLow-deltaMass; massIt<_massHigh+deltaMass; massIt+=stepSize){
    
    Vector4<double>  res4V(massIt, 0., 0., 0.);
    
    complex<double> currentBWs=BreitWigner(res4V, _massLow, _widthLow)+BreitWigner(res4V, _massHigh, _widthHigh);
    double weight=norm(currentBWs);
    _massShapeHist->Fill(massIt,weight);
    
    _argandHist->Fill(currentBWs.real(),currentBWs.imag());
    
    //fill here the phase shift histogramms: hint "atan2(imag,real)"
    
    _phaseHist->Fill(massIt, atan2(currentBWs.imag(), currentBWs.real()));

    //Kmatrix
  // break up momentum q0
 //    double piMass=0.135;
//     double qlow=sqrt( ( _massLow*_massLow-(piMass+piMass)*(piMass+piMass) ) 
// 		      * ( _massLow*_massLow - (piMass-piMass) * (piMass-piMass) ) ) / (_massLow);

//     double qhigh=sqrt( ( _massHigh*_massHigh-(piMass+piMass)*(piMass+piMass) ) 
// 		       * ( _massHigh*_massHigh - (piMass-piMass) * (piMass-piMass) ) ) / (_massHigh);

    // break up momentum qab
 
    // phase space factor rhoLow
//     double  rhoLow=2*qlow/_massLow;

//    double  rhoHigh=2*qhigh/_massHigh;
    // phase space factor rhoMab
//     double  rhoMab=2*qab/massIt;


    theTMatrix.updateMatrix(massIt);
    complex<double> currentTmatrix=theTMatrix(0,0);
    cout << "TMatrix:\t" <<  currentTmatrix << endl;
    weight=norm(currentTmatrix); 
    _massShapeKmatrHist->Fill(massIt,weight);
    _argandKmatrHist->Fill(currentTmatrix.real(),currentTmatrix.imag()); 
    _phaseKmatrHist->Fill(massIt, atan2(currentTmatrix.imag(), currentTmatrix.real()));

    theTMatrix.updateMatrixRel(massIt);
    complex<double> currentTmatrixRel=theTMatrix(0,0);
    cout << "TMatrixRel:\t" <<  currentTmatrixRel << endl;
    weight=norm(currentTmatrixRel); 
    _massShapeKmatrHistRel->Fill(massIt,weight);
    _argandKmatrHistRel->Fill(currentTmatrixRel.real(),currentTmatrixRel.imag()); 
    _phaseKmatrHistRel->Fill(massIt, atan2(currentTmatrixRel.imag(), currentTmatrixRel.real()));
    
  }


}

TwoPoles::~TwoPoles()
{
   _theTFile->Write();
   _theTFile->Close();
}