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

gStyle->SetHistLineWidth(2); 

cout<<"Please use this macro to guide you through the neutrino exercise."<<endl;
cout<<"It will become useful from point 7 onwards, when you can start editting it."<<endl;
cout<<"Just add numbers and delete commands as you go along (look for return commands :-)."<<endl<<endl<<endl;
 
//delete or comment the following line to go to exercise 7
return;
 
cout<<"Ex: 7: here you have an example of a fixed function for the reactor spectrum,"<<endl;
cout<<"you need to write also the oscillation function, with variable parameters."<<endl;
cout<<"Then you need to enter also the relevant distance for the Far Detector (from 6);"<<endl;
cout<<"and notice functions can be defined has opperations on previously defined functions"<<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.",1.8,8.2); 

cout<<"what is the survival probability function, what does it depend on? TF1 *oscil=??"<<endl<<endl; 
return;








// survival probability function - with parameters distance=[0],mixing=[1],dm2=[2]
TF1 *oscil=new TF1("oscil","1.-[1]*pow(sin(1.27*[2]*[0]*1.e3/x),2.)",1.8,8.2); //- this is the solution!

// fix the oscillation parameters 
oscil->SetParameter(1,0.856); oscil->FixParameter(1,0.856); 
oscil->SetParameter(2,8.e-5); oscil->FixParameter(2,8.e-5);

// set the distance for each detector in km
double dNear=0.5;//km
cout<<" what is the distance (in km) you have chosen? double dFar=?"<<endl<<endl; 
return; 

double dFar=100.; //km -- replace by your solution (you can try testing other values)

cout<<"In each detector we see the spectrum modulated by the survival probability"<<endl; 
oscil->SetParameter(0,dNear); 
TF1 *near=new TF1("near","reactor*oscil",1.8,8.2); near->SetLineColor(kBlue);
oscil->SetParameter(0,dFar); 
TF1 *far=new TF1("far","reactor*oscil",1.8,8.2); far->SetLineColor(kRed);
near->SetTitle("Spectrum; E [MeV]; expected flux");
near->SetMinimum(0.); near->SetMaximum(0.25);
TCanvas *functions=new TCanvas(); near->Draw(); far->Draw("same");

cout<<"And can calculate numbers of events (with Integral)"<<endl;
double nearVisibleFraction=near->Integral(1.8,8.2);
double farVisibleFraction=far->Integral(1.8,8.2);

cout<<endl
    <<nearVisibleFraction<<" ND visible and "<<farVisibleFraction<<" FD visible spectra"<<endl<<endl;

//delete or comment the following line to go to exercise 8
return; 

cout<<"Ex: 8: use FillDetection(DET, N , resolution, spectrum) to fill histogram DET"<<endl;
cout<<"------ with N events from an input spectrum function with a given resolution"<<endl;
cout<<"------ use the spectra functions created above, and try changing N and resolution"<<endl<<endl;

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 value for EventsInNearDetector and resolution"<<endl; 
return;
double EventsInNearDetector=1000; //replace by the value you want to test
double resolution=0.0001; //replace by the value you want to test

//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/100.*farVisibleFraction);

new TCanvas; FillDetection(ND,nearEv,resolution,near); ND->DrawClone("hist,pe"); 
new TCanvas; FillDetection(FD,farEv,resolution,far); FD->DrawClone("hist,pe");

cout<<endl<<ND->Integral()<<" ND events and "<<FD->Integral()<<" FD events"<<endl<<endl<<endl;

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

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,100.,1.); 
division->Draw("pe");

cout<<"What does this division correspond too?"<<endl<<endl;
return;







oscil->SetParameter(0,dFar); oscil->FixParameter(0,dFar); 
oscil->ReleaseParameter(1);oscil->ReleaseParameter(2);

division->Fit("oscil","R","",2.5,7.5); //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<<endl<<endl<<"these are the results for "<<EventsInNearDetector<<" events and "<<resolution*100<<"% resolution"<<endl;

return;
}
void FillDetection(TH1D *detection,int Ntotal,double res1, TF1 *distrib)
{

   TRandom tmp;
   for(int iev=0;iev<Ntotal;iev++)
	{
	  double TrueEnergy=distrib->GetRandom(); //choose energy from the given spectrum distribution
	  double resolutionE=sqrt(TrueEnergy)*res1; //from statistical fluctuation on # photo-electrons
	  double RecoEnergy=tmp.Gaus(TrueEnergy,resolutionE); //gaussian fluctuations on reconstruction
	  detection->Fill(RecoEnergy,1.); // this is what will be seen for this 1 event
	 }
   return;   
}