//************************************************************************//
//									  //
//  Copyright 2013 Bertram Kopf (bertram@ep1.rub.de)			  //
//  	      	   Julian Pychy (julian@ep1.rub.de)			  //
//          	   - Ruhr-Universität Bochum 				  //
//									  //
//  This file is part of Pawian.					  //
//									  //
//  Pawian is free software: you can redistribute it and/or modify	  //
//  it under the terms of the GNU General Public License as published by  //
//  the Free Software Foundation, either version 3 of the License, or 	  //
//  (at your option) any later version.	 	      	  	   	  //
//									  //
//  Pawian is distributed in the hope that it will be useful,		  //
//  but WITHOUT ANY WARRANTY; without even the implied warranty of	  //
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the	  //
//  GNU General Public License for more details.	      		  //
//									  //
//  You should have received a copy of the GNU General Public License     //
//  along with Pawian.  If not, see <http://www.gnu.org/licenses/>.	  //
//									  //
//************************************************************************//

#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 "PwaDynamics/KPole.hh"
#include "PwaDynamics/KMatrixRel.hh"
#include "PwaDynamics/TMatrixRel.hh"
#include "PwaDynamics/KMatrixNonRel.hh"
#include "PwaDynamics/TMatrixNonRel.hh"
#include "PwaDynamics/AbsPhaseSpace.hh"
#include "PwaDynamics/PhaseSpaceIsobar.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);
  std::shared_ptr<KPole> kPoleLow(new KPole(low_gis, _massLow));

  vector<double> high_gis;
  high_gis.push_back(gHigh);
  // KpoleMatrix kPoleHigh(high_gis, _massHigh, theMassDecPairVec);
  std::shared_ptr<KPole> kPoleHigh(new KPole(high_gis, _massHigh));

  // vector<KpoleMatrix> kPoles;
  // kPoles.push_back(kPoleLow);
  // kPoles.push_back(kPoleHigh);
  std::vector<std::shared_ptr<KPole> > kPoles;
   kPoles.push_back(kPoleLow);
   kPoles.push_back(kPoleHigh);

   std::vector<std::shared_ptr<AbsPhaseSpace> > phpVec;
   std::shared_ptr<AbsPhaseSpace> firstPhp(new PhaseSpaceIsobar(piMass, piMass));
   phpVec.push_back(firstPhp); 

   // TMatrix theTMatrix(kPoles);
   std::shared_ptr<KMatrixRel> theKMatrixRel( new KMatrixRel(kPoles, phpVec));
   TMatrixRel theTMatrixRel(theKMatrixRel);


   std::shared_ptr<KMatrixNonRel> theKMatrixNonRel( new KMatrixNonRel(kPoles, phpVec));
   TMatrixNonRel theTMatrixNonRel(theKMatrixNonRel);  

  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;


    theTMatrixNonRel.evalMatrix(massIt);
    complex<double> currentTmatrix=theTMatrixNonRel(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()));

    theTMatrixRel.evalMatrix(massIt);
    complex<double> currentTmatrixRel=theTMatrixRel(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();
}

