{

//tutto relativo
//nbins dispari a partire nbins/sigma
//
//
//
// 
// 
//
//

gROOT->Reset();

Int_t mode=1;
Int_t niter=500; //n of iterations

Int_t i,c;
Int_t nbins=20;
Double_t integral;
const Double_t PI=3.1416;
Double_t nevents=10000;
Double_t sigma_nom=1;//mm
Double_t sigma=sigma_nom; //mm
Double_t *x2;
Double_t naux=0.6;
Double_t nbins_s;
Double_t size;
Double_t dx,ss,mm,ss_fit,mm_fit;
TF1 *profile, *profile1, *profile2;
TGraphErrors *gprofile, *gprofile2;
TMultiGraph *mg;
TH1D *hsigma;
Double_t x1,y1,ex,ey;
Char_t title[5];
TList *hlist, *glist;
TF1 *fg;
Int_t j;
Char_t s[40];
TFile *froot=new TFile();
//cout << "Input NEVENTS: " << endl;
//cin >> nevents;
//cout << "Input MODE: " << endl;
//cin >> mode;

sprintf(title,"simulate_mode%d_nev%d.out",mode,nevents);
ofstream fallout(title, ios::out);
naux=-0.5;
sprintf(title,"simulate_mode%d_nev%d.root",mode,nevents);
froot->Open(title,"recreate");
for(Int_t nb=0;nb<4;nb++) {
	 nbins_s=(pow(10,naux));
	 nbins=(Int_t)(nbins_s*16);
	 if((nbins % 2)==0) nbins++;
	 size=15; //window size mm
	 dx=size/nbins; //mm per bin
	 naux=naux+0.05;
	 cout << "NBINS= " << nbins << " NBINS/S= " << nbins_s << endl;
	 
	 x2=new  Double_t[nbins+1]; //vector of bin centers 
	 x2[0]=-(Double_t)(nbins)/2;
	 
	 for(i=1;i<nbins+1;i++) 
		  x2[i]=x2[i-1]+1;
	 profile=new TF1("profile","gaus",-nbins/2,nbins/2);//ideal profile
	 profile2=new TF1("profile2","gaus",-nbins/2,nbins/2);//ideal profile to be shifted
	 gprofile=new TGraphErrors();
	 gprofile2=new TGraphErrors();
	 profile ->SetParameters(nevents/(sqrt(2*PI)*nbins_s),0,nbins_s); 
	 profile2->SetParameters(nevents/(sqrt(2*PI)*nbins_s),0+1/2,nbins_s);	 
	 
	 for(i=0;i<(int)nbins;i++) {
		  integral=profile->Integral(x2[i],x2[i+1]);
		  gprofile->SetPoint(i,x2[i]+0.5,integral);
		  gprofile->SetPointError(i,0,sqrt(integral));
	 }
	 for(i=0;i<(int)nbins;i++) {
		  integral=profile2->Integral(x2[i],x2[i+1]);
		  gprofile2->SetPoint(i,x2[i],integral);
		  gprofile2->SetPointError(i,0,sqrt(integral));
	 }
	 //in gprofile there is the gaussian with the mean values to 
	 //be given to each point...
  
	 hlist=new TList();
 
	 for(i=0;i<(int)nbins;i++) {
		  if(mode==1) {
			   gprofile->GetPoint(i,x1,y1);
			   ey=gprofile->GetEY();
		  }
		  else {
			   gprofile2->GetPoint(i,x1,y1);
			   ey=gprofile2->GetEY();
		  }
		  fg=new TF1("fg","gaus",y1-5*sqrt(y1),y1+5*sqrt(y1));
		  fg->SetParameters(1,y1,sqrt(y1));
		  sprintf(title,"h%d",i);
		  hlist->Add(new TH1D(title,title,niter,y1-5*sqrt(y1),y1+5*sqrt(y1)));
		  ((TH1D *)hlist->Last())->FillRandom("fg",niter);
 
	 }
 
	 glist=new TList();
	 hsigma=new TH1D("hsigma","hsigma",100,nbins_s*(1-0.4),nbins_s*(1+0.4));
	 mg=new TMultiGraph();
	 c=1;
	 for(i=0;i<niter;i++) {
		  glist->Add(new TGraph());  
		  for(j=0;j<(int)nbins;j++) {
			   ((TGraph *)glist->Last())->SetPoint(j,x2[j],((TH1D *)hlist->At(j))->GetRandom());
		  }
		  ((TGraph *)glist->Last())->SetMarkerColor(c);
		  ((TGraph *)glist->Last())->SetLineColor(c++);
		  if(c==10) c=1;
		  profile1=new TF1("profile1","gaus",-size/2,size/2);//ideal profile used for fit
		  ((TGraph *)glist->Last())->Fit("profile1","q");
    
		  hsigma->Fill(profile1->GetParameter(2));
	
		  ((TGraph *)glist->Last())->GetFunction("profile1")->Delete();
		  ((TGraph *)glist->Last())->SetMarkerSize(.1);
		  mg->Add(((TGraph *)glist->Last()));
	
	 }
 
	 sprintf(s,"gprofile, N events %d",nevents);
	 gprofile->SetNameTitle(s,s);
	 gprofile->Write();
  
	 sprintf(s,"MG, N events %d",nevents);
	 mg->SetNameTitle(s,s);
	 mg->Write();
  
	 ss=hsigma->GetRMS(1);
	 mm=hsigma->GetMean();

	 hsigma->Fit("gaus","q");
	 ss_fit=hsigma->GetFunction("gaus")->GetParameter(2);
	 mm_fit=hsigma->GetFunction("gaus")->GetParameter(1);
	 sprintf(s,"hsigma_nbins_%d",nb); 
	 hsigma->SetTitle(s);
	 hsigma->Write();

	 fallout.setf(ios::left);
	 fallout.precision(4);
	 fallout << 
		  nevents     << '\t' <<
		  niter       << '\t' <<
		  size        << '\t' <<
		  nbins       << '\t' <<
		  nbins_s     << '\t' <<
		  sigma_nom   << '\t' <<
		  mm       << '\t' <<
		  ss       << '\t' <<
		  mm_fit   << '\t' <<
		  ss_fit <<
		  endl;
	    
	 delete [] x2;
	 for(i=0;i<hlist->GetSize();i++) ((TH1D*)hlist->At(i))->Delete();
	 delete hlist;
	 delete hsigma;
	 delete profile; delete fg; delete gprofile; delete glist;
	 delete profile1; delete profile2; delete gprofile2;
	 delete mg;
  
	 cout << nb << ": Bins/sig = " << nbins_s << " Mean = " << mm << " Spread = " << ss << endl;
	 cout << nb << ": Bins/sig = " << nbins_s << " Mean = " << mm_fit << " Spread = " << ss_fit << endl;
}
fallout.close();
froot->Close();
}