#include #include #include #include #include #include #include #include #include using namespace std; int N = 1E6; // number of events float Emin = 1.804; // energy threshold float Emax = 8; // max neutrino energy int N_div = 200; // dE double a = 3./100; // energy resolution double p[7]; // parameters double sig[7]; // uncertainties TH1F * datasetNHres; // [0] --> sin^2 2*theta_12 // [1] --> sin^2 2*theta_13 // [2] --> Delta^2 m_21 // [3] --> Delta^2 m_31 // [4] --> Baseline // [5] --> Normalization // [6] --> Resolution // [7] --> Visible energy // x --> Integrated out Energy string fluxExpr = "(0.58 * exp(0.870 - 0.160*x - 0.0910*x^2) +" "0.30 * exp(0.896 - 0.239*x - 0.0981*x^2) +" "0.07 * exp(0.976 - 0.162*x - 0.0790*x^2) +" "0.05 * exp(0.793 - 0.080*x - 0.1085*x^2))"; string crossSecExpr = "(x-1.293) * sqrt( (x-1.293)*(x-1.293) - 0.511*0.511 )"; string PeeExpr = "1 - (pow(1+sqrt(1-[1]),2)/4)*[0]*pow(sin(1.27*[2]*[4]/x),2)" "- (1+sqrt(1-[0])/2)*[1]*pow(sin(1.27*[3]*[4]/x),2)" "- (1-sqrt(1-[0])/2)*[1]*pow(sin(1.27*([3]-[2])*[4]/x),2)"; string gaussExpr = "1./sqrt(2*3.14*[6]*[6]*x) * exp(-0.5*pow( (x-[7])/([6]*sqrt([7])) ,2))"; string intExpr = fluxExpr + "*" + crossSecExpr + "*(" + PeeExpr + ")*" + gaussExpr +"*[5]"; TF1 * integrand = new TF1( "integrand" , intExpr.c_str() , 1.804 , 8 ); void myFCNres( int &npar , double * gin , double &f , double * par , int flag ) { double chi2 = 0; for ( int i = 1 ; i <= datasetNHres->GetSize()-2 ; ++i ) { par[7] = Emin + (Emax-Emin)/(2*N_div) + (i-1)*(Emax-Emin)/N_div; chi2 += pow( ( integrand->Integral( 1.804 , 8 , par , 1E-15) - datasetNHres->GetBinContent(i)) / datasetNHres->GetBinError(i) ,2); } chi2 += pow( (par[0] - p[0]) / sig[0] ,2) + pow( (par[1] - p[1]) / sig[1] ,2) + pow( (par[2] - p[2]) / sig[2] ,2) + pow( (par[3] + p[3]) / sig[3] ,2); f = chi2; } void minuitFitRes_local() { p[0] = 0.857; // sin^2 2*theta_12 // dati dei giappi //0.851 p[1] = 0.089; // sin^2 2*theta_13 //0.0945 p[2] = 7.50E-5; // Delta^2 m_21 //7.54E-5 p[3] = 2.32E-3; // Delta^2 m_31 <---- Segno gerarchia //2.43E-3 p[4] = 52500; // Lenght in m sig[0] = 0.024; // incertezze dai giappi sig[1] = 0.005; sig[2] = 0.2*1E-5; sig[3] = 0.1*1E-3; sig[4] = 0; TStyle * gStyle = new TStyle(); gStyle->cd(); //gStyle->SetOptFit(1111); gStyle->SetOptStat(11); gStyle->SetErrorX(0); // RETRIEVE RESOLUTED PLOTS TDirectory* currentDir = gDirectory; TFile* datasetFile = new TFile( "datasetRes.root" , "READ" ); currentDir->cd(); datasetNHres = dynamic_cast(datasetFile->Get("spectrum_NH_res")->Clone()); //TH1F * datasetIH = dynamic_cast(datasetFile->Get("fluctuating_dataset_IH")->Clone()); delete datasetFile; // FITTING (with NH data) // Correct theory -> spectrumNH TMinuit * gMinuit = new TMinuit(6); gMinuit->SetFCN(myFCNres); gMinuit->SetPrintLevel(-1); double arglist[10]; // vettore con parametri per le istruzioni arglist[0] = 1; // value of Delta_chi2 used // to calculate errors (UP) int ierflg = 0; // flag gMinuit->mnexcm( "SET ERR" , arglist , 1 , ierflg ); // setting of UP // gMinuit->SetErrorDef(1); // same thing // set parameters gMinuit->mnparm( 0 , "sin^2 2*theta_12" , p[0] , 1E-4 , p[0]-sig[0] , p[0]+sig[0] , ierflg ); gMinuit->mnparm( 1 , "sin^2 2*theta_13" , p[1] , 1E-4 , p[1]-sig[1] , p[1]+sig[1] , ierflg ); gMinuit->mnparm( 2 , "Delta^2 m_21 " , p[2] , 1E-8 , p[2]-sig[2] , p[2]+sig[2] , ierflg ); gMinuit->mnparm( 3 , "Delta^2 m_31 " ,-p[3] , 1E-6 ,-p[3]-sig[3] ,-p[3]+sig[3] , ierflg ); gMinuit->mnparm( 4 , "Baseline " , p[4] , 1 , 52000 , 53000 , ierflg ); gMinuit->mnparm( 5 , "Normalization " , 2E4 , 1E3 , 0 , 0 , ierflg ); // fix parameters gMinuit->FixParameter(4); //gMinuit->FixParameter(0); //gMinuit->FixParameter(1); //gMinuit->FixParameter(2); //gMinuit->FixParameter(3); //gMinuit->FixParameter(5); // OR // arglist [0] = 5; // beware of numbering! // gMinuit->mnexcm( "FIX" , arglist , 1 , ierflg ); // SET STRATEGY // 1 standard // 2+ try to improve minimum (slower) arglist [0] = 1; gMinuit->mnexcm( "SET STR" , arglist , 1 , ierflg ); // minimization step arglist[0] = 500; // number of steps arglist[1] = 1; // tolerance (required precision for the value of the minimum) //gMinuit->mnexcm( "SEEK" , arglist , 1 , ierflg ); // MC preliminar search gMinuit->mnexcm( "SIMPLEX" , arglist , 2 , ierflg ); // rough first search //gMinuit->mnexcm( "MIGRAD" , arglist , 2 , ierflg ); arglist[1] = 0; // set parameter numbers for MINOS arglist[2] = 1; arglist[3] = 2; arglist[4] = 3; arglist[5] = 4; arglist[6] = 5; //gMinuit->mnexcm ( "MINOS" , arglist , 7 , ierflg ); // better error estimation // OR // gMinuit->SetMaxIteration(500); // gMinuit->Migrad(); double amin, edm, errdef; // printout results int nvpar, nparx, icstat; gMinuit->mnstat( amin , edm , errdef , nvpar , nparx , icstat ); gMinuit->mnprin( 1 , amin ); double par[8], parErr[6]; // retrieve results gMinuit->GetParameter( 0 , par[0] , parErr[0] ); gMinuit->GetParameter( 1 , par[1] , parErr[1] ); gMinuit->GetParameter( 2 , par[2] , parErr[2] ); gMinuit->GetParameter( 3 , par[3] , parErr[3] ); gMinuit->GetParameter( 4 , par[4] , parErr[4] ); gMinuit->GetParameter( 5 , par[5] , parErr[5] ); par[6] = a; TGraph * plot = new TGraph(); for ( int i = 1 ; i <= 1000 ; ++i ) { par[7]= 1.804 + i*0.0062; plot->SetPoint( i , 1.804 + i*0.0062 , integrand->Integral( 1.804 , 8 , par , 1E-15 )); } TCanvas * c1 = new TCanvas( "c1" , "c1" , 1 ); // print graph c1->SetLogx(); c1->cd(); datasetNHres->Draw("E1"); plot->Draw("SAME"); c1->SaveAs("plotMinuitRes.png"); return; }