#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 "PAnyDistribution.h"
#include "plugins/scatter_mod/PScatterCrossSection.h"
#include "plugins/sigmatoLee/PSigmatoLee.h"


void tmp_YtoLambdaee() //int input, int events
{

  int events = 100;
  int input = 10;
  double ebeam = 11.0;
  bool xsection = 1;
  gROOT->Reset();
  PUtils::SetSeed(input*1000);
makeStaticData()->AddParticle(70,"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(71,"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);
  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);

  PReaction *my_reaction = new PReaction("_P1 = 11.0","e-","p","e- Delta1910+ [Sigma13850 [dilepton Lambda] K+]","YtosigmaK",0,0,0,0);

  my_reaction->Preheating(500);
  my_reaction->Print();
  my_reaction->Loop(events);
//  

}