// example of use of TGenPhaseSpace for pi->enuee (E Frlez, 06/18/2014)

using namespace std;

#define matrix_pi_enuee matrix_pi_enuee_
extern "C" float matrix_pi_enuee(float P[], float &FIN, float &AIN, float &RIN);

#include "TGenPhaseSpace.h"
#include "TMath.h"
#include "TStyle.h"
#include "TColor.h"
#include "TLatex.h"
#include <TCanvas.h>
#include <TMultiLayerPerceptron.h>
#include <TTree.h>
#include <TF1.h>
#include <TH1.h>
#include <TH2.h>
#include <TBranch.h>
#include <TFile.h>
#include <TCanvas.h>
#include <TGraph2D.h>
#include <TSystem.h>
#include <fstream>
#include <iostream>
#include <TRandom.h>
#include <TRandom2.h>
#include <TPad.h>
#include <TSystem.h>

#include <stdio.h>
#include <stdlib.h>

#include <fcntl.h>

//void gen_pienuee(float P, float FIN=0, float AIN=0, float RIN=0) {
  void gen_pienuee( ) {

  //if (!gROOT->GetClass("TGenPhaseSpace"));

   gSystem->Load("libPhysics");
   gSystem->Load("libMATRIX_PI_ENUEE");

   gStyle->SetPalette(1); // rainbow pallete or below:

   TLorentzVector target(0.0, 0.0, 0.0, 0.13957018);
   TLorentzVector beam(0.0, 0.0, 0.0, 0.0);
   TLorentzVector W = beam + target;

   //(Momentum, Energy units are Gev/C, GeV)
   Double_t masses[4] = { 0.000510998910, 0.0, 0.000510998910, 0.000510998910} ;

   TGenPhaseSpace event;
   event.SetDecay(W, 4, masses);

   TFile *outfile  = new TFile(Form("gen_pienuee.root"),"RECREATE");
   TTree* TreePI3E   =  new TTree("PI3E","PI3E Tree");

   TH1F *h1_ee1 = new TH1F("h1_ee1","ee1", 100,0,140);
   TH1F *h1_ee2 = new TH1F("h1_ee2","ee2", 100,0,140);
   TH1F *h1_ee3 = new TH1F("h1_ee3","ee3", 100,0,140);
   TH1F *h1_SQME = new TH1F("h1_SQME","SQME", 1000,0,1.0);

   Float_t SQME;
   Double_t weight=1; 
   TreePI3E->Branch("weight",   &weight, "weight/F");

   float FIN=1,AIN=1,RIN=1;
   float P[5][18];

   for (Int_t n=0;n<10000;n++) {

      Int_t acc=0;

      weight = event.Generate();

      TLorentzVector *pe1 = event.GetDecay(0);
      TLorentzVector *pnu = event.GetDecay(1);
      TLorentzVector *pe2 = event.GetDecay(2);
      TLorentzVector *pe3 = event.GetDecay(3);

      //P[0][0]=(float) pe1(0); P[1][0]=(float) pe1(1); P[2][0]=(float) pe1(2); P[3][0]=(float) pe1(3);
      //P[0][2]=(float) pnu(0); P[1][1]=(float) pnu(1); P[2][1]=(float) pnu(2); P[3][1]=(float) pnu(3);
      //P[0][2]=(float) pe2(0); P[1][2]=(float) pe2(1); P[2][2]=(float) pe2(2); P[3][2]=(float) pe2(3);
      //P[0][3]=(float) pe3(0); P[1][3]=(float) pe3(1); P[2][3]=(float) pe3(2); P[3][3]=(float) pe3(3);
      P[0][0]=pe1->Px(); P[1][0]=pe1->Py(); P[2][0]=pe1->Pz(); P[3][0]=pe1->Pz();
      P[0][1]=pnu->Px(); P[1][1]=pnu->Py(); P[2][1]=pnu->Pz(); P[3][1]=pnu->Pz();
      P[0][2]=pe2->Px(); P[1][2]=pe2->Py(); P[2][2]=pe2->Pz(); P[3][2]=pe2->Pz();
      P[0][3]=pe3->Px(); P[1][3]=pe3->Py(); P[2][3]=pe3->Pz(); P[3][3]=pe3->Pz();

      float SQME = matrix_pi_enuee(*P, FIN, AIN, RIN);

      if (n%1000==0) printf(" \n\n");

      Float_t tht_e1=pe1->Theta();
      Float_t phi_e1=pe1->Phi();
      Float_t tht_e2=pe2->Theta();
      Float_t phi_e2=pe2->Phi();
      Float_t tht_e3=pe3->Theta();
      Float_t phi_e3=pe3->Phi();
      Float_t ctht12=sin(tht_e1)*sin(tht_e2)+cos(phi_e1-phi_e2)+cos(tht_e1)*cos(tht_e2);
      Float_t tht12=acos(ctht12);
      Float_t ctht13=sin(tht_e1)*sin(tht_e3)+cos(phi_e1-phi_e3)+cos(tht_e1)*cos(tht_e3);
      Float_t tht13=acos(ctht13);
      Float_t ctht23=sin(tht_e2)*sin(tht_e3)+cos(phi_e2-phi_e3)+cos(tht_e2)*cos(tht_e3);
      Float_t tht23=acos(ctht23);

      h1_ee1->Fill(1000*pe1->E(),weight);
      h1_ee2->Fill(1000*pe2->E(),weight);
      h1_ee3->Fill(1000*pe3->E(),weight);

      if (tht_e1>30 && tht_e2>30 && tht_e3>30 ) acc=1;

      TreePI3E->Fill();

      if (n%1000==0) printf(" event: %d\n",n); 
  }

   //gStyle->SetMarkerSize(0.1);
   //h2_SDp->SetStats(0); h2_SDp->SetTitle(0); h2_SDp->SetLineStyle(1);  
   //h2_SDp->SetLineColor(kRed); h2_SDp->SetXTitle(Form("E_{#gamma} (m_{#pi}/2)")); 
   //h2_SDp->SetYTitle("E_{e+} (m_{#pi}/2)");

   TLatex l;
   l.SetTextFont(72);
   l.SetTextSize(0.06);
   l.SetTextAlign(12);

   gStyle->SetPadLeftMargin(0.14);
   gStyle->SetPadRightMargin(0.14);

   /*
   TCanvas* c1 = new TCanvas("c1","RPD",10,10,700,700);
   h2_SDp->Draw("colz");
   c1->Print("SDp.pdf");
   l.DrawLatex (0.2,0.2,Form("SD^{+}"));
   c1->Update();
   h2_SDm->Draw("colz");
   l.DrawLatex (0.2,0.2,Form("SD^{-}"));
   c1->Update();
   c1->Print("SDm.pdf");
   c1->Update();
   gPad->SetLogz(1);
   h2_IB->Draw("colz");
   l.DrawLatex (0.2,0.2,Form("IB"));
   c1->Update();   
   c1->Print("IB.pdf");
   gPad->SetLogz(1);
   h2_INT->Draw("colz");
   l.DrawLatex (0.2,0.2,Form("INT"));
   c1->Update();
   c1->Print("INT.pdf");
   */

   outfile->Write();
   return;
}