//This macro is supposed to be used for the Hands on Neutrino Workshop

void FillDetection(TH1D*,int events, double resol,TF1 *spectrum, TF1 *oscil);
void Neutrino_2023()
{
  gStyle->SetHistLineWidth(2); 

  cout<<"Please use this macro to guide you through the neutrino's exercises."<<endl;
  cout<<"Just add numbers and delete commands as you proceed (look for *return* commands)."<<endl<<endl<<endl;

  

  
  /********************** STEP 1 **********************/
  cout<<"Step 1: Here you have an example of a function describing the reactor antineutrino energy spectrum"<<endl;
  cout<<"The reactor antineutrino energy spectrum can be described by a normalized Gaussian having a mean value of 4 MeV and a sigma of 2 MeV"<<endl;
  cout<<"A normalized gaussian function in Root has 2 parameters: f(x) = exp(-0.5*((x-[0])/[1])^2)/([1]*sqrt(2*Pi()))"<<endl<<endl;
  cout<<"Reactor function = TMath::Gaus(x,[0],[1],0)/([1]*sqrt(2.*TMath::Pi())), where [0] = 4.0 and [1] = 2.0"<<endl<<endl<<endl;

  
  // input spectrum as a normalized gaussian - structure of TF1 is name/function/lower x/upper x
  //Once completed the TF1 *reactor function comment the line below
  return;

  TF1 *reactor=new TF1("reactor","---write-function-here---",0.0,10.0);

  reactor->Draw(); 
  return;
  //comment '//' the two lines above and paste STEP 2 to proceed with the exercise

  
}





//This is the FillDetection function -- do NOT modify this function
//The function generates an histogram given Ntotal events, an energy resultion res1, and the two initial functions: the reactor function (orig) and the oscillation function (oscil)
void FillDetection(TH1D *detection,int Ntotal,double res1, TF1 *orig, TF1 *oscil)
{
  double IBDthr = 1.8;//MeV	This is the threhold of your reaction
  double emass = 0.511;//MeV	This is the electron's mass
  TRandom tmp;
  
  for(int iev = 0; iev < Ntotal; iev++)//Loop over the input number of events
    {
      double NeutrinoEnergy=orig->GetRandom(); //obtain the energy from the given spectrum distribution

      //What will we detect?
      if(tmp.Uniform()>oscil->Eval(NeutrinoEnergy)) continue; //oscillated to another neutrino type, do not plot it 
      if(NeutrinoEnergy<IBDthr) continue; //not enough energy to make inverse beta decay, do not plot it
      
      //visible energy
      double Epositron = NeutrinoEnergy - IBDthr + emass;
      double resolutionE = sqrt(Epositron)*res1; //apply the energy resolution to the detected events, from statistical fluctuation on # photo-electrons
      double RecoEnergy=tmp.Gaus(Epositron,resolutionE); //Now smear your detected energy function
      
      //back to neutrino
      double NeutrinoRecoE = RecoEnergy - emass + IBDthr;
      detection->Fill(NeutrinoRecoE,1.); // this is what will be seen for this 1 event
    }
  
  return;   
}