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

void FillDetection(TH1D*,int events, double resol,TF1 *spectrum, TF1 *oscil);
void Neutrino_2021()
{

  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;
  
  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
  TF1 *reactor=new TF1("reactor","TMath::Gaus(x,4.,2.,0)/sqrt(2.*TMath::Pi())/2.",0.0,10.0);

  //reactor->Draw(); 
  //return;
  //comment '//' the two lines above to proceed with the exercise.
 
  cout<<"Step 2: You now need to write also the oscillation function, with variable parameters."<<endl;
  cout<<"Go back to your exercise sheet and look for the oscillation probability function you just obtained."<<endl;
  cout<<"Starting from the oscillation probability function, try to write in root the survival probability function as a funcion of the antineutrino energy, x."<<endl;
  cout<<"The 3 parameters that you need are [0] = the distance (in km), [1] = the mixing, and [2] = the Dm^2"<<endl;
  cout<<"Write the function in 'TF1 *oscil' below. Hint: f(x)^n in root is pow(f(x),n)"<<endl<<endl<<endl;

  //Once completed the TF1 *oscil function comment the line below
  //return;

  //Survival probability function - with parameters distance in km=[0],mixing=[1],Dm2=[2]

  TF1 *oscil=new TF1("oscil","1.-[1]*pow(sin(1.27*[2]*[0]*1.e3/x),2.)",0.,10.);
  //TF1 *oscil=new TF1("oscil","write here your function",0.,10.); //- this is the solution for the distance in km!

  //Set and fix the oscillation parameters [1] and [2]
  oscil->SetParameter(1,0.856); oscil->FixParameter(1,0.856);
  oscil->SetParameter(2,8.e-5); oscil->FixParameter(2,8.e-5);
 
  cout<<"Then you need to enter also the relevant distance for the Far Detector."<<endl<<endl<<endl;

  //Far detector
  //Set the distance for the far detector in km
  double dFar= 62; //km -- set here the distance of the far detector
  oscil->SetParameter(0,dFar);

  //oscil->Draw(); 
  //return;
  //comment '//' the two lines above to proceed with the exercise.
  
  cout<<"Step 3: The actual survival probability is defined as the product of the reactor antineutrino energy spectrum (the reactor function) and the survival probability function."<<endl; 
  cout<<"Plot the survival probability function for the near and the far detector."<<endl; 
  cout<<"Do they agree with your expectations? What does happen if you change the near and the far distances? What is the optimal choice?"<<endl<<endl<<endl; 

  //Near detector
  //Set the distance for the near detector in km
  double dNear= 0.1;//km -- set here the distance of the near detector
  oscil->SetParameter(0,dNear);
  
  TF1 *near=new TF1("near","reactor*oscil",0.,10.);
  near->SetLineColor(kRed);
  near->SetTitle("Spectrum; E [MeV]; expected flux");
  near->SetMinimum(0.); near->SetMaximum(0.25);

  oscil->SetParameter(0,dFar);
  TF1 *far=new TF1("far","reactor*oscil",0.,10.);
  far->SetLineColor(kBlue);

  //Draw the functions
  TCanvas *functions=new TCanvas();
  near->Draw();
  far->Draw("same"); 
  reactor->SetLineColor(kBlack); reactor->SetLineStyle(4); reactor->Draw("same");

  //Once done with this step, comment the line below
  //return;

  cout<<"We can now calculate the numbers of events in each spectrum (using the function Integral)"<<endl;
  cout<<"What would be the mininum value of energy, x, to calculate your integral?"<<endl;

  //return;
  //Substitute x with the value of the energy's threshold. Once done comment the line above
  double x = 1.8 ;
  double nearVisibleFraction = near->Integral(x,8.2);
  double farVisibleFraction = far->Integral(x,8.2);

  cout<<"Fraction of ND events "<<nearVisibleFraction <<endl;
  cout<<"Fraction of FD events "<<farVisibleFraction<<endl<<endl<<endl;

  //Comment the following line to proceed with the exercise
  //return;
    
  cout<<"Step 4: Use FillDetection(DET, N , resolution, original spectrum, oscillation function) to fill histogram DET"<<endl;
  cout<<"------ where N is the number of events generated using an input spectrum function with a given resolution"<<endl;
  cout<<"------ use the reactor (original spectrum) and oscilation functions created above, and try changing N and resolution"<<endl<<endl;

  //Let's create two DET hisograms: the near, ND, and the far, FD.
  TH1D *ND=new TH1D("ND","Near Detector; E [MeV]; # events",32,1.8,8.2); ND->SetLineColor(kBlue);
  TH1D *FD=new TH1D("FD","Far Detector; E [MeV]; # events",32,1.8,8.2); FD->SetLineColor(kRed); 
  
  ND->Sumw2(); FD->Sumw2(); //this line is necessary for root to store statistical uncertainties
  
  cout<<"Write in your macro the value for EventsInNearDetector, the resolution and calculate the value for size"<<endl<<endl<<endl;

  //Comment the line below once filled the EventsInNearDetector and resolution below
  //return;
  
  double EventsInNearDetector=10000000; //replace by the value you want to test. These are the events that you want to see in ND
  double resolution= 0.02; //replace by the value you want to test
  
  double size=100; //This should be approximately equal to the ratio (distance_far-distance_near)^2/(mass_far/mass_near)
  double EventsInFarDetector=EventsInNearDetector/size; 
  
  //now get approximate number of events to generate in each detector and how many will be visible
  int nearEv=(int)(EventsInNearDetector/nearVisibleFraction);
  int farEv=(int)(nearEv/size*farVisibleFraction);
  
  
  //Please see at the end of the macro for the FillDetection function
  new TCanvas;
  oscil->SetParameter(0,dNear); 
  FillDetection(ND,EventsInNearDetector,resolution,reactor,oscil);
  ND->DrawClone("hist,pe");
  
  new TCanvas;
  oscil->SetParameter(0,dFar); 
  FillDetection(FD,EventsInFarDetector,resolution,reactor,oscil);
  FD->DrawClone("hist,pe");
  
  cout<<endl<<"Fraction of ND events "<<ND->Integral()<<" and fraction of FD events "<<FD->Integral()<<endl<<endl<<endl;

  //Comment the line below once done.
  //return;

  cout<<"Next let's fit the oscillation parameters."<<endl; 

  new TCanvas;
  
  cout<<"Histograms can be divided and fitted with a given funtion"<<endl;
  cout<<"------ HistFim -> Divide (H1, H2, c1, c2) means HistFim = (c1.H1)/(c2.H2)"<<endl;
  cout<<"------ the free parameters of function are given in order as defined above"<<endl<<endl<<endl;
  
  TH1D *division=(TH1D*)FD->Clone("division");
  division->Reset(); 
  division->SetLineColor(kBlack);
  division->SetTitle("Far Spectrum / Near Spectrum; E [MeV]; FD/ND");
  division->Divide(FD,ND,1.,1./size); 
  division->Draw("pe");
  
  cout<<"What does this ratio correspond too?"<<endl<<endl;
  
  oscil->Draw("same");

  //Comment this line once done
  //return;
  
  
  oscil->SetParameter(0,dFar); oscil->FixParameter(0,dFar); 
  oscil->ReleaseParameter(1); oscil->ReleaseParameter(2);
  oscil->SetParameter(1,0.5); oscil->SetParameter(2,1.e-5); //You can change these parameters here

  //division->Fit("oscil","R","",2.5,7.5); //you can choose the fitting range, here 2.5 MeV to 7.5 MeV
  division->Fit("oscil"); //you can choose the fitting range, here 2.5 MeV to 7.5 MeV
  
  cout<<endl<<endl<<endl;
  cout<<"Oscillation Amplitude = "<<oscil->GetParameter(1)<<" +/- "<<oscil->GetParError(1)<<endl;
  cout<<"Delta Mass Squared = ("<<oscil->GetParameter(2)<<" +/- "<<oscil->GetParError(2)<<") eV^2"<<endl;
  cout<<endl<<"These are the results for "<<EventsInNearDetector<<" events and "<<resolution*100<<"% resolution"<<endl;
  
  return;
}



//This is the FillDetection function -- do NOT modify this function
//The function generate 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 smaer 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;   
}