//Program to generate a PLUTO root file with the option of selecting the decay of titled particle
//Author Michael C. Kunkel
//Branching is the branching ratio found in pdg. I.e. for eta' 1 is eta' decay to pi+ pi- eta while 2 is rho gamma decay of eta'
//The decimal is used for the decay of eta' daughter particle, i.e. eta' -> w gamma where w -> pi+ pi- pi0 would be 4.1 because 4 is the w gamma decay of eta' and .1 is the first decay of w.
//If channels want to be added, please contact mkunkel@jlab.org

#include "TGraphAsymmErrors.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TChain.h"
#include "TCanvas.h"
#include "TF1.h"
#include "TGraphAsymmErrors.h"
#include "TFile.h"
#include "TSystem.h"

#include "PParticle.h"
#include "PInclusiveModel.h"
#include "PReaction.h"
#include "PBeamSmearing.h"
#include "PDalitzDecay.h"
#include "PAnyDistribution.h"
#include "plugins/scatter_mod/PScatterCrossSection.h"
//#include "plugins/sigmatoLee/PSigmatoLee.h"


void tmp_YtoLambdaee() //int input, int events
{
  
  int events = 10000;
  int input = 10;
  double ebeam = 11.0;
  bool xsection = 1;
  gROOT->Reset();
  PUtils::SetSeed(input*10);
  

makeStaticData()->AddParticle(170,"Sigma13850", 1.3837);	
makeStaticData()->SetParticleTotalWidth("Sigma13850",0.0360);
makeStaticData()->SetParticleBaryon("Sigma13850",1);
makeStaticData()->SetParticleLMass("Sigma13850",1.115);
makeStaticData()->AddDecay("Sigma13850 --> Lambda + dilepton", "Sigma13850", "Lambda, dilepton", 1.0);
//makeStaticData()->AddDecay("dilepton --> e+ + e-", "dilepton", "e+, e-", 1.0);


makeStaticData()->AddParticle(171,"Delta1910+", 1.910);	
makeStaticData()->SetParticleTotalWidth("Delta1910+",0.340);
makeStaticData()->SetParticleBaryon("Delta1910+",1);
makeStaticData()->SetParticleLMass("Delta1910+",1.89);
makeStaticData()->AddDecay("Delta1910+ --> Sigma13850 + K+", "Delta1910+", "Sigma13850, K+", 1.0);


//~~~~//makeStaticData()->AddDecay("Lambda --> p + pi- + dilepton", "Lambda", "p, pi-, dilepton",1.0);
//~~~~  
//~~~~  

//  PAnyDistribution* decay =
//  new PAnyDistribution("v_slope", "A function to add a new t-slope");
//  decay->Add("q,     parent");
//  decay->Add("e-,     daughter");
//  decay->Add("Delta1910+,     daughter");
//  TH1F * cache  = new TH1F ("cache","v cache",400,0,11.0);
//  decay->AddEquation(cache,"beam = _parent->GetBeam(); beam->Boost(_parent) ; t1 = (beam - [e-,1])->E(); _x = t1;");
//  decay->AddEquation("_f = exp(-1000*t1);");
//  makeDistributionManager()->Add(decay);





  PScatterCrossSection *model = new PScatterCrossSection("mymodel","My cross section");
  model->Add("e-,grandparent,beam");
  model->Add("p,grandparent,target");
  model->Add("q,parent");
  model->Add("e-,daughter,spectator");
  model->Add("Delta1910+,daughter,primary");
  
  double M_P = 0.938272;   //Proton mass
  double M_el = 0.000510999;   //Electron mass
  
  //Define the range of the c.m. sampling
  double En_beam = TMath::Sqrt(ebeam*ebeam + M_el*M_el);
  double q= TMath::Sqrt((2.*En_beam*M_P)+(M_P*M_P) + (M_el*M_el));
  double q_min = q - q/10000.;
  double q_max = q + q/10000.;
  //Now let set up the model
  
  model->SetRange(q_min,q_max);
 
   TFile *fhist = new TFile("Hists_Williams.root");
                    
   TH1D *example42 = (TH1D*)fhist->Get("hist42");
   model->AddHistogram(example42,"if (_y > 1.49 ) value = Eval(_x); _f =value*pow((1.49*1.49),7)/pow((_y*_y),7)");

  model->SetNpy(250);
  model->SetNpx(250);  //Use this after SetNpy analysis
  makeDistributionManager()->Add(model);

  //Set up reaction as usual
  char lund_file[75];
  
  cout<<"#########################################################################"<<endl;
  cout<<"################## Cross-Section Genertation is occuring ################"<<endl;
  cout<<"#########################################################################"<<endl;
  sprintf(lund_file,"%s","etaPXSection.lund");
/////////disribution of Sigma1385_0 to Lambda dileption
//PSigmatoLee  *sigmadecaytoLee = new PSigmatoLee("sigma13850_to_e+_e-_lambda_matrix","sigma13850 to e+e-_Lambda Distribution");
//sigmadecaytoLee->Add("Sigma13850,parent");
//sigmadecaytoLee->Add("Lambda,daughter");
//sigmadecaytoLee->Add("dilepton,primary");
//sigmadecaytoLee->SetMax(0.2);
//makeDistributionManager()->Add(sigmadecaytoLee);

PDalitzDecay *pmodel = new PDalitzDecay("sigma13850dalitz_decay@Sigma13850_to_Lambda_dilepton","Sigma13850 Dalitz decay",-1);
makeDistributionManager()->Add(pmodel);

  
  
  
  //Define reaction
  
 // PReaction *my_reaction = new PReaction("_P1 = 11.0","e-","p","e- Delta1910+ [K+ Sigma13850 [dilepton [e+ e-] Lambda [dilepton [e+ e-] p pi-]]]","YtosigmaK",0,0,0,0);
  //PReaction *my_reaction = new PReaction("_P1 = 11.0","e-","p","e- Delta1910+ [K+ Sigma13850 [dilepton [e+ e-] Lambda]]","YtosigmaK",0,0,0,0);
  PReaction *my_reaction = new PReaction("_P1 = 11.0","e-","p","e- Delta1910+ [Sigma13850 [dilepton [e+ e-] Lambda] K+]","YtosigmaK",0,0,0,0);
//  PReaction *my_reaction = new PReaction("_P1 = 11.0","e-","p","e- Delta1910+ [Sigma13850 K+]","YtosigmaK",0,0,0,0);
//  PReaction *my_reaction = new PReaction("_P1 = 11.0","e-","p","e- Delta1910+","YtosigmaK",0,0,0,0);
//~~~~
//~~~~  
  my_reaction->Preheating(500);
  my_reaction->Print();
  my_reaction->Loop(events);
//  
  
  
  
}