//#include "FullAmp.h"
#include "Intensity.h"
#include "Resonance.h"
#include "Particle.h"
#include <complex>
#include <vector>
#include <TGraph.h>
#include <TAxis.h>
#include <TCanvas.h>
#include "TH2.h"
#include "TH1.h"
#include "TLorentzVector.h"
#include "TLorentzRotation.h"
#include "TGenPhaseSpace.h"
#include <fstream>
#include <iostream> 
#include <string>
#include "TParticlePDG.h"
#include "TDatabasePDG.h"
#include "QNCooking.h"
using namespace TMath;
  //(Momentum, Energy units are GeV/C, GeV)
int main(int __argc,char *__argv[]){ 
  /// Read parameters set by user ----------------------------------
  int sizemass = 10;
  std::string nStr="";

  // decode arguments
  if( __argc>1 && ( strcmp( __argv[1], "-help" ) == 0
		    || strcmp( __argv[1], "--help" ) == 0 ) ){

    std::cout << "This is Dalitz Plot generation with parameters\n"
	      <<"-n Number of events \n"
	      <<"Have fun! \n"
	      << std::endl;
    return 0;
  } 
  while ((optind < (__argc-1) ) && (__argv[optind][0]=='-')) {
    bool found=false;
    std::string sw = __argv[optind];
    if (sw=="-n"){
      optind++;
      nStr = __argv[optind];
      found=true;
    }

    if (!found){
      std::cout<< "Unknown switch: "
	       << __argv[optind] <<std::endl;
      optind++;
    }
  
  while ( (optind < __argc ) && __argv[optind][0]!='-' ) optind++; 
  }

  std::stringstream nSStr(nStr);
  nSStr >> sizemass;
  ///--------------------------------------------------------------

  /// Open\Read file(s) -----------------------------------
  std::ofstream output("DalitzPlotIntencity.txt");
  output<<"-------------------------"<<std::endl;
  output<<"Format: "<<std::endl;
  output<<"p1 (4-vector)"<<std::endl;
  output<<"p2 (4-vector)"<<std::endl;
  output<<"p3 (4-vector)"<<std::endl;
  output<<"intensity (complex)"<<std::endl;
  output<<"-------------------------"<<std::endl;



  //input parameters-------------------------------------
  TDatabasePDG *fdbPDG = TDatabasePDG::Instance();
  int idFP[3];//final particles PDGid
  int idIP;//initial particle PDGid
  double MassDP;//mass of initial particle
  double sDP;//spin of initial particle NB: hasn't used!
  int pDP;//P parity of initial particle NB: hasn't used!
  // double massf1; // mass of 1st finale particle
  // double massf2; // mass of 2nd finale particle
  // double massf3; // mass of 3nd finale particle
  double mFP[3]; // mass of #st finale particle
  double sFP[3]; // spin of #st finale particle
  int pFP[3]; // P parity of #st finale particle
  int idR[24];// resonances IDs
  int idRd[24][2];// IDs resonances decay particles
  double mR[24];// resonances masses
  double wR[24];// resonances widths
  double sR[24];// resonances spin
  int pR[24];// resonances P parity
  for(int ir=0;ir<24;ir++){
    idRd[ir][0]=-1; idRd[ir][1]=-1;
  }
  //-----------------------------------------------------

  // decode arguments from input file
  std::string idStr="";
  std::string pStr="";
  int ip=0;
  ifstream input("DalitzPlotParams.txt");
  std::string sout;
  while(!input.eof()){
      input>>sout;
      if (sout=="Particles"){
	while(sout!="}"){
	  input>>sout;
	  if(sout=="#"){
	    input>>sout;
	    //	    output<<sout<<" decay to:"<<std::endl;
	    output<<sout;
	    input>>sout;
	    output<<"(PDGid="<<sout<<") decay to:"<<std::endl;
	    idStr = sout;
	    std::stringstream idSStr(idStr);
	    idSStr >> idIP;
	    input>>sout;
	    pStr=sout;
	    std::stringstream pSStr(pStr);
	    pSStr>>MassDP;
	    input>>sout;
	    pStr=sout;
	    std::stringstream pSStr1(pStr);
	    pSStr1>>sDP;
	    input>>sout;
	    pStr=sout;
	    std::stringstream pSStr2(pStr);
	    pSStr2>>pDP;
	  }
	  if(sout=="*"){
	    input>>sout;
	    output<<sout<<" ";
	    input>>sout;
	    output<<sout<<std::endl;
	    idStr = sout;
	    std::stringstream idSStr(idStr);
	    idSStr >> idFP[ip];
	    input>>sout;
	    pStr=sout;
	    std::stringstream pSStr(pStr);
	    pSStr>>mFP[ip];
	    input>>sout;
	    pStr=sout;
	    std::stringstream pSStr1(pStr);
	    pSStr1>>sFP[ip];
	    input>>sout;
	    pStr=sout;
	    std::stringstream pSStr2(pStr);
	    pSStr2>>pFP[ip];
	    ip++;
	  }

	}
      }

      if(sout=="Resonances"){
	output<<"Resonance(s):"<<std::endl;
	int ir=0;
	while(sout!="}"){
	  input>>sout;
	  if(sout=="*"){
	    input>>sout;
	    output<<" "<<sout;
	    input>>sout;
	    output<<"(PDGid="<<sout<<") decay to: ";
	    idStr = sout;
	    std::stringstream idSStr(idStr);
	    idSStr >> idR[ir];

	    input>>sout;
	    pStr = sout;
	    std::stringstream pSStr(pStr);
	    pSStr >> mR[ir];

	    input>>sout;
	    pStr = sout;
	    std::stringstream pSStr1(pStr);
	    pSStr1 >> wR[ir];

	    input>>sout;
	    pStr = sout;
	    std::stringstream pSStr2(pStr);
	    pSStr2 >> sR[ir];

	    input>>sout;
	    pStr = sout;
	    std::stringstream pSStr3(pStr);
	    pSStr3 >> pR[ir];

	    input>>sout; 
	    output<<sout<<" and ";
	    input>>sout;
	    idStr = sout;
	    std::stringstream idSStrp1(idStr);
	    idSStrp1 >> idRd[ir][0];
	    input>>sout;
	    output<<sout<<std::endl;
	    input>>sout;
	    idStr = sout;
	    std::stringstream idSStrp2(idStr);
	    idSStrp2 >> idRd[ir][1];
	    ir++;
	  }
	}
      }
  }

  ///--------------------------------------------------------------

  ///DP -> f1 f2 f3
  //initial particle
  //  TParticlePDG *inP = fdbPDG->GetParticle(idIP);
  // inP->Print();
  //final particles
  // TParticlePDG *outP[3] = {fdbPDG->GetParticle(idFP[0]), fdbPDG->GetParticle(idFP[1]), fdbPDG->GetParticle(idFP[2])};
  // outP[0]->Print();  outP[1]->Print();   outP[2]->Print();
  // double MassDP = inP->Mass();//mass of initial particle
  // double massf1 = outP[0]->Mass(); // mass of 1st finale particle
  // double massf2 = outP[1]->Mass(); // mass of 2nd finale particle
  // double massf3 = outP[2]->Mass(); // mass of 3nd finale particle
  Particle pi1(mFP[0],sFP[0],pFP[0]);
  Particle pi2(mFP[1],sFP[1],pFP[1]);
  Particle pi3(mFP[2],sFP[2],pFP[2]);
  pi1.Print();
  pi2.Print();
  pi3.Print();
  /// Setting Resonance(s) and Calculation of Quantum numbers --------------
  ///resonances for pair (1,2) ---------------------------------- //TODO: check coherent\incoheren sum!
  std::vector<Resonance> resonance12;
  for(int ir=0;ir<24;ir++){
    //  std::cout<<"ir#"<<ir<<" decay to: "<<idRd[ir][0]<<" and "<<idRd[ir][1]<<std::endl;
    if(idRd[ir][0]==-1 && idRd[ir][1]==-1) break;
    if((idRd[ir][0]==idFP[0] && idRd[ir][1]==idFP[1]) || 
       (idRd[ir][0]==idFP[1] && idRd[ir][1]==idFP[0])){
      Resonance f0(mR[ir],wR[ir],sR[ir],pR[ir]);
      QNCooking *QNCalc = new QNCooking(f0,pi1,pi2);
      std::vector<double> l;
      QNCalc->AngularMomenta(l);
      f0.setl(l);
      //      int sRes = sR[ir];
      resonance12.push_back(f0);
      //   f0.Print();
      //   std::cout<<"resonance #"<<ir<<" was taken for pair(1,2): "<<idFP[0]<<" and "<<idFP[1]<<std::endl;
    }
  }
  ///-----------------------------------------------------------
  ///resonances for pair (2,3) ---------------------------------- //TODO: check coherent\incoheren sum!
  std::vector<Resonance> resonance23;
  for(int ir=0;ir<24;ir++){
    //   std::cout<<"ir#"<<ir<<" decay to: "<<idRd[ir][1]<<" and "<<idRd[ir][2]<<std::endl;
    if(idRd[ir][0]==-1 && idRd[ir][1]==-1) break;
    if((idRd[ir][0]==idFP[1] && idRd[ir][1]==idFP[2]) || 
       (idRd[ir][0]==idFP[2] && idRd[ir][1]==idFP[1])){
      Resonance f0(mR[ir],wR[ir],sR[ir],pR[ir]);
      QNCooking *QNCalc = new QNCooking(f0,pi2,pi3);
      std::vector<double> l;
      QNCalc->AngularMomenta(l);
      f0.setl(l);
      //      int sRes = sR[ir];
      resonance23.push_back(f0);
      //   f0.Print();
      //    std::cout<<"resonance #"<<ir<<" was taken for pair(2,3): "<<idFP[1]<<" and "<<idFP[2]<<std::endl;
    }
  }
  ///-----------------------------------------------------------
 ///resonances for pair (1,3) ---------------------------------- //TODO: check coherent\incoheren sum!
  std::vector<Resonance> resonance13;
  for(int ir=0;ir<24;ir++){
    //  std::cout<<"ir#"<<ir<<" decay to: "<<idRd[ir][0]<<" and "<<idRd[ir][2]<<std::endl;
    if(idRd[ir][0]==-1 && idRd[ir][1]==-1) break;
    if((idRd[ir][0]==idFP[0] && idRd[ir][1]==idFP[2]) || 
       (idRd[ir][0]==idFP[2] && idRd[ir][1]==idFP[0])){
      Resonance f0(mR[ir],wR[ir],sR[ir],pR[ir]);
      QNCooking *QNCalc = new QNCooking(f0,pi1,pi3);
      std::vector<double> l;
      QNCalc->AngularMomenta(l);
      f0.setl(l);
      //      int sRes = sR[ir];
      resonance13.push_back(f0);
      //  f0.Print();
      // std::cout<<"resonance #"<<ir<<" was taken for pair(1,3): "<<idFP[0]<<" and "<<idFP[2]<<std::endl;
    }
  }
  ///-----------------------------------------------------------

  Intensity *f0_pi_pi_s = new Intensity(resonance12,resonance23, resonance13, pi1, pi2, pi3); 
  //TO DO: check J, Mj of vector<resonance>
 //  f0_pi_pi_s->PrintPar();

  /// Test Histogramm ------------------------------------------
  double minx = TMath::Power((mFP[0]+mFP[1]),2)-0.1;
  double maxx = TMath::Power((MassDP-mFP[2]),2)+0.1;
  double miny = TMath::Power((mFP[1]+mFP[2]),2)-0.1;
  double maxy = TMath::Power((MassDP-mFP[0]),2)+0.1;
  //  std::cout<<minx<<", "<<maxx<<", "<<miny<<", "<<maxy<<std::endl;
  TH2F *h2 = new TH2F("h2","Uniform phase space di stribution;m^{2}_{12}, GeV^{2};m^{2}_{23}, GeV^{2}", 100,minx,maxx,100,miny,maxy); //test
  TH1F *h1m12i = new TH1F("h1m12i",";m_{12};Events",100,TMath::Sqrt(minx+0.1)-0.1,TMath::Sqrt(maxx));
  TH1F *h1m23i = new TH1F("h1m23i",";m_{23};Events",100,TMath::Sqrt(miny+0.1)-0.1,TMath::Sqrt(maxy));
  TH2F *h2r = new TH2F("h2r","Result phase space distribution;m^{2}_{12}, GeV^{2};m^{2}_{23}, GeV^{2}", 100, minx,maxx, 100,miny,maxy); //test
 TH1F *h1m12 = new TH1F("h1m12",";m_{12};Events",100,TMath::Sqrt(minx+0.1)-0.1,TMath::Sqrt(maxx));
 TH1F *h1m23 = new TH1F("h1m23",";m_{23};Events",100,TMath::Sqrt(miny+0.1)-0.1,TMath::Sqrt(maxy));
 minx = TMath::Power((mFP[0]+mFP[2]),2)-0.1;
 maxx = TMath::Power((MassDP-mFP[1]),2)+0.1;
 // std::cout<<minx<<", "<<maxx<<", "<<miny<<", "<<maxy<<std::endl;
 TH2F *h2r2 = new TH2F("h2r2","Result phase space distribution;m^{2}_{13}, GeV^{2};m^{2}_{23}, GeV^{2}", 100, minx,maxx, 100,miny,maxy); //test
   TH1F *h1m13 = new TH1F("h1m13",";m_{13};Events",100,TMath::Sqrt(minx+0.1)-0.1,TMath::Sqrt(maxx));
  TH1F *hinttest = new TH1F("h1r","Intensity;I", 1000,0,10); //test
  TCanvas *c1 = new TCanvas("test");
  c1->Divide(3,3);
  ///--------------------------------------------------------- 
 
 
 
  TLorentzVector W(0.0,0.0,0.0,MassDP);
  //  Double_t masses[3] = {massf1, massf2, massf3} ;
  Double_t *masses = mFP;
  TGenPhaseSpace event;
  event.SetDecay(W, 3, masses);
  // const int sizemass=20000;

  for (Int_t n=0;n<sizemass;n++) {
    /// Generation uniform phase space -------------------------------------
    /// and relative angle calculation in the rest system of resonance -----
    /// (which is opposite system to the one of vector)
    
    double weight = event.Generate();
    TLorentzVector *p1_cur = event.GetDecay(0);
    TLorentzVector *p2_cur = event.GetDecay(1);
    TLorentzVector *p3_cur = event.GetDecay(2);

    TLorentzVector p2_1= *p2_cur;
    TLorentzVector p3_1= *p3_cur;
    TLorentzVector p1_1= *p1_cur;

    TLorentzVector p1_2 = *p1_cur;
    TLorentzVector p3_2 = *p3_cur;
    TLorentzVector p2_2 = *p2_cur;

    TLorentzVector p1_3 = *p1_cur;
    TLorentzVector p2_3 = *p2_cur;
    TLorentzVector p3_3 = *p3_cur;
    // p1_1.Print();
    // p2_1.Print();
    // p3_1.Print();
    p1_1.SetE(p2_1.E()+p3_1.E()+1e-6);
    TLorentzRotation lt1;
    lt1.RotateZ(p1_1.Phi());
    lt1.RotateY(p1_1.Theta());
    p1_1.Transform(lt1);
    lt1.Boost(p1_1.BoostVector());
    p1_1.Boost(-(p1_1.BoostVector()));
    // std::cout<<"p1_1"<<std::endl;
    // p1_1.Print();
    p2_1.Transform(lt1);
    p3_1.Transform(lt1);

    // p2_1.Print();
    // p3_1.Print();
  

  
    p2_2.SetE(p1_2.E()+p3_2.E()+1e-6);
    TLorentzRotation lt2;
    lt2.RotateZ(p2_2.Phi());
    lt2.RotateY(p2_2.Theta());
    p2_2.Transform(lt2);
    lt2.Boost(p2_2.BoostVector());
    p2_2.Boost(-(p2_2.BoostVector()));
    // std::cout<<"p2_2"<<std::endl;
    //  p2_2.Print();
    p1_2.Transform(lt2);
    p3_2.Transform(lt2);

    //    p2_2.Print();
    // p1_2.Print();
    // p3_2.Print();

 
    p3_3.SetE(p1_3.E()+p2_3.E()+1e-6);
    TLorentzRotation lt3;
    lt3.RotateZ(p3_3.Phi());
    lt3.RotateY(p3_3.Theta());
    p3_3.Transform(lt3);
    lt3.Boost(p3_3.BoostVector());
    p3_3.Boost(-(p3_3.BoostVector()));
    // std::cout<<"p3_3"<<std::endl;
    // p3_3.Print();
    p1_3.Transform(lt3);
    p2_3.Transform(lt3);

    // p3_3.Print();
    // p1_3.Print();
    // p2_3.Print();

    TLorentzVector p12 = *p1_cur + *p2_cur;
    TLorentzVector p23 = *p2_cur + *p3_cur;
    TLorentzVector p13 = *p1_cur + *p3_cur;
    double mass12 = p12.M();
    double mass23 = p23.M();
    double mass13 = p13.M();
    h2->Fill(p12.M2(), p23.M2() ,weight);
    h1m12i->Fill(p12.M(),weight);
    //  std::cout<<" h1m12i filled by: "<<p12.M()<<" "<<weight<<std::endl;
    h1m23i->Fill(p23.M(),weight);
    ///---------------------------------------------------------

 

    ///Intensity calculation in helicity formalizm -------------
    ///phi and theta are calculated by boost
     std::complex<double> itensity_sum = f0_pi_pi_s->Calc(p1_3.Phi(),p1_3.Theta(),mass12,p2_1.Phi(),p2_1.Theta(),mass23,
							  p1_2.Phi(),p1_2.Theta(),mass13);
     double intensbin = real(itensity_sum)*weight; 
    // std::complex<double> itensity_sum2 = f0_pi_pi_s->Calc(p2_1.Phi(),p2_1.Theta(),mass23,
    // 							  p1_2.Phi(),p1_2.Theta(),mass13,p1_3.Phi(),p1_3.Theta(),mass12);

    // std::complex<double> itensity_sum3 = f0_pi_pi_s->Calc(p1_2.Phi(),p1_2.Theta(),mass13,p1_3.Phi(),p1_3.Theta(),mass12,p2_1.Phi(),p2_1.Theta(),mass23);
    // double intensbin = real(itensity_sum+itensity_sum2+itensity_sum3)*weight; 

    // for(int idem=0;idem<10;idem++){
    //   if(resonance12[idem].size()!=0 || resonance23[idem].size()!=0 || resonance13[idem].size()!=0){
    // 	Intensity *f0_pi_pi_p = new Intensity(resonance12[idem],resonance23[idem], resonance13[idem], pi1, pi2, pi3);
    // 	itensity_sum+=f0_pi_pi_p->Calc(p1_3.Phi(),p1_3.Theta(),mass12,p2_1.Phi(),p2_1.Theta(),mass23,
    // 				       p1_2.Phi(),p1_2.Theta(),mass13);
    //   }
    // }
    
    //    double intensbin = real(itensity_sum); 
    // std::cout<<"itensity_sum="<<itensity_sum<<" from "<<real(itensity_sum)<<std::endl;
    // std::cout<<"==============================="<<std::endl;
  
    hinttest->Fill(intensbin);
    h2r->Fill(p12.M2(), p23.M2(),intensbin);
    h2r2->Fill(p13.M2(),p23.M2(),intensbin);
    h1m12->Fill(p12.M(),intensbin);
    h1m23->Fill(p23.M(),intensbin);
    h1m13->Fill(p13.M(),intensbin);
    output<<" "<<std::endl;
    output<<(*p1_cur).Px()<<" "<<(*p1_cur).Py()<<" "<<(*p1_cur).Pz()<<" "<<(*p1_cur).E()<<std::endl;
    output<<(*p2_cur).Px()<<" "<<(*p2_cur).Py()<<" "<<(*p2_cur).Pz()<<" "<<(*p2_cur).E()<<std::endl;
    output<<(*p3_cur).Px()<<" "<<(*p3_cur).Py()<<" "<<(*p3_cur).Pz()<<" "<<(*p3_cur).E()<<std::endl;
    output<<itensity_sum<<std::endl;

  ///-----------------------------------------------------------

  }
  c1->cd(1);
  h2->Draw("COLZ");
  c1->cd(2);
  //  h2->ProjectionX()->Draw();
  h1m12i->Draw();
  //h1m12i->DrawCopy();
  c1->cd(3);
  //  h2->ProjectionY()->Draw();
  h1m23i->Draw();
  c1->cd(4);
  h2r->Draw("COLZ");
  c1->cd(5);
  // h2r->ProjectionX()->Draw();
  h1m12->Draw();
  c1->cd(6);
  // h2r->ProjectionY()->Draw();
  h1m23->Draw();
  c1->cd(7);
  h2r2->Draw("COLZ");
  c1->cd(8);
  // h2r2->ProjectionX()->Draw();
  h1m13->Draw();
  c1->cd(9);
  h1m23->Draw();
  // h2r2->ProjectionY()->Draw();
  c1->SaveAs("testDalitzPlot.pdf");
  
  TCanvas *c2 = new TCanvas("DalitzPlot");
  h2r->Draw("COLZ");
  c2->SaveAs("testDalitzPlot_onlyDP.root");

  TCanvas *c3 = new TCanvas("IntencityTest");
  hinttest->Draw();
  c3->SaveAs("testIntenc.root");

}