/********************** STEP 5 **********************/
  cout<<"Step 5: Next let's fit the oscillation parameters."<<endl; 

  new TCanvas;
  
  cout<<"Histograms can be divided and fitted with a given function"<<endl;
  cout<<"------ Histogram -> Divide (H1, H2, c1, c2) means Histogram = (c1.H1)/(c2.H2)"<<endl;
  cout<<"------ the free parameters of the function 'oscil' 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./EventRatio); 
  division->Draw("pe");
  
  cout<<"What does this ratio correspond to?"<<endl<<endl;
  
  //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");
  //division->Fit("oscil","R","",2.5,7.5); //you can also 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;

  oscil->Draw("same");

  return;