00001 #include "TwoFlavourFitter.h"
00002
00003 #include "TemplateReader.h"
00004 #include "TemplateWriter.h"
00005 #include "TemplateCalculator.h"
00006 #include "TemplateGrid.h"
00007
00008 #include "ParametersUtil.h"
00009 #include "Configuration.h"
00010 #include "Likelihood.h"
00011
00012 #include "TString.h"
00013 #include "TRandom.h"
00014 #include "TMarker.h"
00015 #include "TGraph.h"
00016 #include "TH2.h"
00017 #include "TH1.h"
00018
00019 #include "TFile.h"
00020 #include "TTree.h"
00021 #include "TDirectory.h"
00022
00023 #include <iostream>
00024 #include <cmath>
00025 #include <ctime>
00026 #include <cassert>
00027
00028 using namespace OscFit;
00029
00030 static TwoFlavourFitter* fgTwoFlavourFitter = 0;
00031
00032 static void TwoFlavourLikelihoodFunction(Int_t&, Double_t*, Double_t& f, Double_t* par, Int_t)
00033 {
00034 Bool_t printDebug = 0;
00035
00036 Double_t osc_dmsq = fgTwoFlavourFitter->fOscDmsq;
00037 Double_t osc_sinsq = fgTwoFlavourFitter->fOscSinsq;
00038 Double_t beam_norm = fgTwoFlavourFitter->fBeamNorm;
00039 Double_t beam_ncbkg = fgTwoFlavourFitter->fBeamNCBkg;
00040 Double_t beam_shwen = fgTwoFlavourFitter->fBeamShwEn;
00041 Double_t beam_trken = fgTwoFlavourFitter->fBeamTrkEn;
00042 Double_t beam_trken_exit = fgTwoFlavourFitter->fBeamTrkEnExit;
00043 Double_t atmos_norm_cv = fgTwoFlavourFitter->fAtmosNormCV;
00044 Double_t atmos_norm_rock = fgTwoFlavourFitter->fAtmosNormRock;
00045 Double_t atmos_chg_cv = fgTwoFlavourFitter->fAtmosChgCV;
00046 Double_t atmos_chg_rock = fgTwoFlavourFitter->fAtmosChgRock;
00047 Double_t atmos_nue = fgTwoFlavourFitter->fAtmosNue;
00048 Double_t atmos_ncbkg = fgTwoFlavourFitter->fAtmosNCBkg;
00049 Double_t atmos_spec_cv_nu = fgTwoFlavourFitter->fAtmosSpecNuCV;
00050 Double_t atmos_spec_cv_nubar = fgTwoFlavourFitter->fAtmosSpecNuBarCV;
00051 Double_t atmos_spec_rock_nu = fgTwoFlavourFitter->fAtmosSpecNuRock;
00052 Double_t atmos_spec_rock_nubar = fgTwoFlavourFitter->fAtmosSpecNuBarRock;
00053 Double_t atmos_zenith = fgTwoFlavourFitter->fAtmosZenith;
00054 Double_t atmos_trken = fgTwoFlavourFitter->fAtmosTrkEn;
00055 Double_t atmos_trken_exit = fgTwoFlavourFitter->fAtmosTrkEnExit;
00056 Double_t atmos_shwen = fgTwoFlavourFitter->fAtmosShwEn;
00057
00058 if( fgTwoFlavourFitter->fFittingDmsq ) osc_dmsq = par[0];
00059 if( fgTwoFlavourFitter->fFittingSinsq ) osc_sinsq = par[1];
00060 if( fgTwoFlavourFitter->fFittingBeamNorm ) beam_norm = par[2];
00061 if( fgTwoFlavourFitter->fFittingBeamNCBkg ) beam_ncbkg = par[3];
00062 if( fgTwoFlavourFitter->fFittingBeamShwEn ) beam_shwen = par[4];
00063 if( fgTwoFlavourFitter->fFittingBeamTrkEn ) beam_trken = par[5];
00064 if( fgTwoFlavourFitter->fFittingBeamTrkEnExit ) beam_trken_exit = par[6];
00065 if( fgTwoFlavourFitter->fFittingAtmosNormCV ) atmos_norm_cv = par[7];
00066 if( fgTwoFlavourFitter->fFittingAtmosNormRock ) atmos_norm_rock = par[8];
00067 if( fgTwoFlavourFitter->fFittingAtmosChgCV ) atmos_chg_cv = par[9];
00068 if( fgTwoFlavourFitter->fFittingAtmosChgRock ) atmos_chg_rock = par[10];
00069 if( fgTwoFlavourFitter->fFittingAtmosNue ) atmos_nue = par[11];
00070 if( fgTwoFlavourFitter->fFittingAtmosNCBkg ) atmos_ncbkg = par[12];
00071 if( fgTwoFlavourFitter->fFittingAtmosSpecNuCV ) atmos_spec_cv_nu = par[13];
00072 if( fgTwoFlavourFitter->fFittingAtmosSpecNuBarCV ) atmos_spec_cv_nubar = par[14];
00073 if( fgTwoFlavourFitter->fFittingAtmosSpecNuRock ) atmos_spec_rock_nu = par[15];
00074 if( fgTwoFlavourFitter->fFittingAtmosSpecNuBarRock ) atmos_spec_rock_nubar = par[16];
00075 if( fgTwoFlavourFitter->fFittingAtmosZenith ) atmos_zenith = par[17];
00076 if( fgTwoFlavourFitter->fFittingAtmosTrkEn ) atmos_trken = par[18];
00077 if( fgTwoFlavourFitter->fFittingAtmosTrkEnExit ) atmos_trken_exit = par[19];
00078 if( fgTwoFlavourFitter->fFittingAtmosShwEn ) atmos_shwen = par[20];
00079
00080 fgTwoFlavourFitter->GetLnL( osc_dmsq, osc_sinsq,
00081 beam_norm, beam_ncbkg,
00082 beam_shwen, beam_trken,
00083 beam_trken_exit,
00084 atmos_norm_cv, atmos_norm_rock,
00085 atmos_chg_cv, atmos_chg_rock,
00086 atmos_nue, atmos_ncbkg,
00087 atmos_spec_cv_nu, atmos_spec_cv_nubar,
00088 atmos_spec_rock_nu, atmos_spec_rock_nubar,
00089 atmos_zenith,
00090 atmos_trken, atmos_trken_exit,
00091 atmos_shwen,
00092 f );
00093
00094 if( printDebug ){
00095 std::cout << " --- [" << fgTwoFlavourFitter->get_counter() << "]: " << std::endl;
00096 std::cout << " osc_dmsq=" << osc_dmsq << ", osc_sinsq=" << osc_sinsq << std::endl;
00097 std::cout << " beam_norm=" << beam_norm << ", beam_ncbkg=" << beam_ncbkg << ", beam_shwen=" << beam_shwen << ", beam_trken=" << beam_trken << ", beam_trken_exit=" << beam_trken_exit << std::endl;
00098 std::cout << " atmos_norm_cv=" << atmos_norm_cv << ", atmos_norm_rock=" << atmos_norm_rock << ", atmos_chg_cv=" << atmos_chg_cv << ", atmos_chg_rock=" <<atmos_chg_rock << ", atmos_nue=" << atmos_nue << ", atmos_ncbkg=" << atmos_ncbkg << std::endl;
00099 std::cout << " atmos_spec_cv_nu=" << atmos_spec_cv_nu << " atmos_spec_cv_nubar=" << atmos_spec_cv_nubar << ", atmos_spec_rock_nu=" << atmos_spec_rock_nu << ", atmos_spec_rock_nubar=" << atmos_spec_rock_nubar << ", atmos_zenith=" << atmos_zenith << std::endl;
00100 std::cout << " atmos_trken=" << atmos_trken << ", atmos_trken_exit=" << atmos_trken_exit << ", atmos_shwen=" << atmos_shwen << std::endl;
00101 std::cout << " f=" << f << std::endl;
00102 }
00103
00104 fgTwoFlavourFitter->itr();
00105
00106 return;
00107 }
00108
00109 TwoFlavourFitter* TwoFlavourFitter::Instance()
00110 {
00111 if( !fgTwoFlavourFitter ){
00112 fgTwoFlavourFitter = new TwoFlavourFitter();
00113 }
00114
00115 return fgTwoFlavourFitter;
00116 }
00117
00118 TwoFlavourFitter::TwoFlavourFitter()
00119 {
00120 std::cout << " *** Starting Two Flavour Fitter *** " << std::endl;
00121
00122 fLnL68 = 1.15;
00123 fLnL90 = 2.30;
00124 fLnL99 = 4.61;
00125
00126 fLnL90_1D = 1.35;
00127
00128 fSinsqBins = 0;
00129 fSinsqMin = 0.0;
00130 fSinsqMax = 0.0;
00131
00132 fDmsqBins = 0;
00133 fDmsqMin = 0.0;
00134 fDmsqMax = 0.0;
00135
00136 fInputDmsq = 2.40e-3;
00137 fInputSinsq = 0.96;
00138 fInputBeamNorm = 0.0;
00139 fInputBeamNCBkg = 0.0;
00140 fInputBeamShwEn = 0.0;
00141 fInputBeamTrkEn = 0.0;
00142 fInputBeamTrkEnExit = 0.0;
00143 fInputAtmosNormCV = 0.0;
00144 fInputAtmosNormRock = 0.0;
00145 fInputAtmosChgCV = 0.0;
00146 fInputAtmosChgRock = 0.0;
00147 fInputAtmosNue = 0.0;
00148 fInputAtmosNCBkg = 0.0;
00149 fInputAtmosSpecNuCV = 0.0;
00150 fInputAtmosSpecNuBarCV = 0.0;
00151 fInputAtmosSpecNuRock = 0.0;
00152 fInputAtmosSpecNuBarRock = 0.0;
00153 fInputAtmosZenith = 0.0;
00154 fInputAtmosTrkEn = 0.0;
00155 fInputAtmosTrkEnExit = 0.0;
00156 fInputAtmosShwEn = 0.0;
00157
00158 fSeedDmsq = 2.32e-3;
00159 fSeedSinsq = 1.0;
00160 fSeedBeamNorm = 0.0;
00161 fSeedBeamNCBkg = 0.0;
00162 fSeedBeamShwEn = 0.0;
00163 fSeedBeamTrkEn = 0.0;
00164 fSeedBeamTrkEnExit = 0.0;
00165 fSeedAtmosNormCV = 0.0;
00166 fSeedAtmosNormRock = 0.0;
00167 fSeedAtmosChgCV = 0.0;
00168 fSeedAtmosChgRock = 0.0;
00169 fSeedAtmosNue = 0.0;
00170 fSeedAtmosNCBkg = 0.0;
00171 fSeedAtmosSpecNuCV = 0.0;
00172 fSeedAtmosSpecNuBarCV = 0.0;
00173 fSeedAtmosSpecNuRock = 0.0;
00174 fSeedAtmosSpecNuBarRock = 0.0;
00175 fSeedAtmosZenith = 0.0;
00176 fSeedAtmosTrkEn = 0.0;
00177 fSeedAtmosTrkEnExit = 0.0;
00178 fSeedAtmosShwEn = 0.0;
00179
00180 fFitDmsq = 0;
00181 fFitSinsq = 0;
00182 fFitBeamNorm = 0;
00183 fFitBeamNCBkg = 0;
00184 fFitBeamShwEn = 0;
00185 fFitBeamTrkEn = 0;
00186 fFitBeamTrkEnExit = 0;
00187 fFitAtmosNormCV = 0;
00188 fFitAtmosNormRock = 0;
00189 fFitAtmosChgCV = 0;
00190 fFitAtmosChgRock = 0;
00191 fFitAtmosNue = 0;
00192 fFitAtmosNCBkg = 0;
00193 fFitAtmosSpecNuCV = 0;
00194 fFitAtmosSpecNuBarCV = 0;
00195 fFitAtmosSpecNuRock = 0;
00196 fFitAtmosSpecNuBarRock = 0;
00197 fFitAtmosZenith = 0;
00198 fFitAtmosTrkEn = 0;
00199 fFitAtmosTrkEnExit = 0;
00200 fFitAtmosShwEn = 0;
00201
00202 fMaskDmsq = 1;
00203 fMaskSinsq = 1;
00204 fMaskBeamNorm = 1;
00205 fMaskBeamNCBkg = 1;
00206 fMaskBeamShwEn = 1;
00207 fMaskBeamTrkEn = 1;
00208 fMaskBeamTrkEnExit = 1;
00209 fMaskAtmosNormCV = 1;
00210 fMaskAtmosNormRock = 1;
00211 fMaskAtmosChgCV = 1;
00212 fMaskAtmosChgRock = 1;
00213 fMaskAtmosNue = 1;
00214 fMaskAtmosNCBkg = 1;
00215 fMaskAtmosSpecNuCV = 1;
00216 fMaskAtmosSpecNuBarCV = 1;
00217 fMaskAtmosSpecNuRock = 1;
00218 fMaskAtmosSpecNuBarRock = 1;
00219 fMaskAtmosZenith = 1;
00220 fMaskAtmosTrkEn = 1;
00221 fMaskAtmosTrkEnExit = 1;
00222 fMaskAtmosShwEn = 1;
00223
00224 fFittingDmsq = 0;
00225 fFittingSinsq = 0;
00226 fFittingBeamNorm = 0;
00227 fFittingBeamNCBkg = 0;
00228 fFittingBeamShwEn = 0;
00229 fFittingBeamTrkEn = 0;
00230 fFittingBeamTrkEnExit = 0;
00231 fFittingAtmosNormCV = 0;
00232 fFittingAtmosNormRock = 0;
00233 fFittingAtmosChgCV = 0;
00234 fFittingAtmosChgRock = 0;
00235 fFittingAtmosNue = 0;
00236 fFittingAtmosNCBkg = 0;
00237 fFittingAtmosSpecNuCV = 0;
00238 fFittingAtmosSpecNuBarCV = 0;
00239 fFittingAtmosSpecNuRock = 0;
00240 fFittingAtmosSpecNuBarRock = 0;
00241 fFittingAtmosZenith = 0;
00242 fFittingAtmosTrkEn = 0;
00243 fFittingAtmosTrkEnExit = 0;
00244 fFittingAtmosShwEn = 0;
00245
00246 fMockData = 1;
00247 fMockDataWithFluctuations = 0;
00248
00249 fCorrelatedSystematics = 0;
00250
00251 fUseSimplex = 0;
00252
00253 fPseudoThreeFlavours = 0;
00254
00255 fPenaltyTerms = 1;
00256 fPhysicalBoundary = 1;
00257 fDmsqBoundary = 1;
00258
00259 fTemplateMC = new Template();
00260 fTemplateExpt = new Template();
00261
00262 fWriteTemplates = 0;
00263 fWriteContours = 0;
00264
00265 fNumPoints = 150;
00266
00267 fOutFile = "oscfit.root";
00268
00269 gRandom->SetSeed(0);
00270
00271 fMinuit = new TMinuit();
00272 fMinuit->SetPrintLevel(-1);
00273 fMinuit->SetMaxIterations(5000);
00274
00275 fMinuitStatus = 999;
00276
00277 fBigNumber = 1.0e10;
00278
00279 fDebug = 0;
00280 fCtr = 0;
00281 }
00282
00283 TwoFlavourFitter::~TwoFlavourFitter()
00284 {
00285 delete fMinuit;
00286
00287 delete fTemplateMC;
00288 delete fTemplateExpt;
00289 }
00290
00291 void TwoFlavourFitter::Reset()
00292 {
00293 std::cout << " *** TwoFlavourFitter::Reset() *** " << std::endl;
00294
00295 ResetExpt();
00296
00297 return;
00298 }
00299
00300 void TwoFlavourFitter::PrintSettings()
00301 {
00302 Bool_t usingBeamData = UsingBeamData();
00303 Bool_t usingAtmosData = UsingAtmosData();
00304
00305 std::cout << " *** TwoFlavourFitter::PrintSettings() *** " << std::endl;
00306
00307 std::cout << " Data Sets: " << std::endl;
00308 std::cout << " UseBeamData = " << UsingBeamData() << std::endl;
00309 std::cout << " UseAtmosData = " << UsingAtmosData() << std::endl;
00310 std::cout << " Configuration: " << std::endl;
00311 std::cout << " MockData = " << fMockData << std::endl;
00312 std::cout << " (Fluctuations = " << fMockDataWithFluctuations << ") " << std::endl;
00313 std::cout << " CorrelatedSystematics = " << fCorrelatedSystematics << std::endl;
00314 std::cout << " PenaltyTerms = " << fPenaltyTerms << std::endl;
00315 std::cout << " PhysicalBoundary = " << fPhysicalBoundary << std::endl;
00316 std::cout << " PseudoThreeFlavours = " << fPseudoThreeFlavours << std::endl;
00317
00318 std::cout << " Inputs for Building Expt: " << std::endl;
00319 std::cout << " osc_dmsq = " << fInputDmsq << std::endl;
00320 std::cout << " osc_sinsq = " << fInputSinsq << std::endl;
00321 if( usingBeamData ){
00322 std::cout << " beam_norm = " << fInputBeamNorm << std::endl;
00323 std::cout << " beam_ncbkg = " << fInputBeamNCBkg << std::endl;
00324 std::cout << " beam_shwen = " << fInputBeamShwEn << std::endl;
00325 std::cout << " beam_trken = " << fInputBeamTrkEn << std::endl;
00326 std::cout << " beam_trken_exit = " << fInputBeamTrkEnExit << std::endl;
00327 }
00328 if( usingAtmosData ){
00329 std::cout << " atmos_norm_cv = " << fInputAtmosNormCV << std::endl;
00330 std::cout << " atmos_norm_rock = " << fInputAtmosNormRock << std::endl;
00331 std::cout << " atmos_chg_cv = " << fInputAtmosChgCV << std::endl;
00332 std::cout << " atmos_chg_rock = " << fInputAtmosChgRock << std::endl;
00333 std::cout << " atmos_nue = " << fInputAtmosNue << std::endl;
00334 std::cout << " atmos_ncbkg = " << fInputAtmosNCBkg << std::endl;
00335 std::cout << " atmos_spec_cv_nu = " << fInputAtmosSpecNuCV << std::endl;
00336 std::cout << " atmos_spec_cv_nubar = " << fInputAtmosSpecNuBarCV << std::endl;
00337 std::cout << " atmos_spec_rock_nu = " << fInputAtmosSpecNuRock << std::endl;
00338 std::cout << " atmos_spec_rock_nubar = " << fInputAtmosSpecNuBarRock << std::endl;
00339 std::cout << " atmos_zenith = " << fInputAtmosZenith << std::endl;
00340 std::cout << " atmos_trken = " << fInputAtmosTrkEn << std::endl;
00341 std::cout << " atmos_trken_exit = " << fInputAtmosTrkEnExit << std::endl;
00342 std::cout << " atmos_shwen = " << fInputAtmosShwEn << std::endl;
00343 }
00344 std::cout << " Inputs for Seeding Fit: " << std::endl;
00345 std::cout << " osc_dmsq = " << fSeedDmsq << std::endl;
00346 std::cout << " osc_sinsq = " << fSeedSinsq << std::endl;
00347 if( usingBeamData ){
00348 std::cout << " beam_norm = " << fSeedBeamNorm << std::endl;
00349 std::cout << " beam_ncbkg = " << fSeedBeamNCBkg << std::endl;
00350 std::cout << " beam_shwen = " << fSeedBeamShwEn << std::endl;
00351 std::cout << " beam_trken = " << fSeedBeamTrkEn << std::endl;
00352 std::cout << " beam_trken_exit = " << fSeedBeamTrkEnExit << std::endl;
00353 }
00354 if( usingAtmosData ){
00355 std::cout << " atmos_norm_cv = " << fSeedAtmosNormCV << std::endl;
00356 std::cout << " atmos_norm_rock = " << fSeedAtmosNormRock << std::endl;
00357 std::cout << " atmos_chg_cv = " << fSeedAtmosChgCV << std::endl;
00358 std::cout << " atmos_chg_rock = " << fSeedAtmosChgRock << std::endl;
00359 std::cout << " atmos_nue = " << fSeedAtmosNue << std::endl;
00360 std::cout << " atmos_ncbkg = " << fSeedAtmosNCBkg << std::endl;
00361 std::cout << " atmos_spec_cv_nu = " << fSeedAtmosSpecNuCV << std::endl;
00362 std::cout << " atmos_spec_cv_nubar = " << fSeedAtmosSpecNuBarCV << std::endl;
00363 std::cout << " atmos_spec_rock_nu = " << fSeedAtmosSpecNuRock << std::endl;
00364 std::cout << " atmos_spec_rock_nubar = " << fSeedAtmosSpecNuBarRock << std::endl;
00365 std::cout << " atmos_zenith = " << fSeedAtmosZenith << std::endl;
00366 std::cout << " atmos_trken = " << fSeedAtmosTrkEn << std::endl;
00367 std::cout << " atmos_trken_exit = " << fSeedAtmosTrkEnExit << std::endl;
00368 std::cout << " atmos_shwen = " << fSeedAtmosShwEn << std::endl;
00369 }
00370
00371 std::cout << " Fitting Parameters: " << std::endl;
00372 std::cout << " osc_dmsq = " << fFitDmsq << std::endl;
00373 std::cout << " osc_sinsq = " << fFitSinsq << std::endl;
00374 if( usingBeamData ){
00375 std::cout << " beam_norm = " << fFitBeamNorm << std::endl;
00376 std::cout << " beam_ncbkg = " << fFitBeamNCBkg << std::endl;
00377 std::cout << " beam_shwen = " << fFitBeamShwEn << std::endl;
00378 std::cout << " beam_trken = " << fFitBeamTrkEn << std::endl;
00379 std::cout << " beam_trken_exit = " << fFitBeamTrkEnExit << std::endl;
00380 }
00381 if( usingAtmosData ){
00382 std::cout << " atmos_norm_cv = " << fFitAtmosNormCV << std::endl;
00383 std::cout << " atmos_norm_rock = " << fFitAtmosNormRock << std::endl;
00384 std::cout << " atmos_chg_cv = " << fFitAtmosChgCV << std::endl;
00385 std::cout << " atmos_chg_rock = " << fFitAtmosChgRock << std::endl;
00386 std::cout << " atmos_nue = " << fFitAtmosNue << std::endl;
00387 std::cout << " atmos_ncbkg = " << fFitAtmosNCBkg << std::endl;
00388 std::cout << " atmos_spec_cv_nu = " << fFitAtmosSpecNuCV << std::endl;
00389 std::cout << " atmos_spec_cv_nubar = " << fFitAtmosSpecNuBarCV << std::endl;
00390 std::cout << " atmos_spec_rock_nu = " << fFitAtmosSpecNuRock << std::endl;
00391 std::cout << " atmos_spec_rock_nubar = " << fFitAtmosSpecNuBarRock << std::endl;
00392 std::cout << " atmos_zenith = " << fFitAtmosZenith << std::endl;
00393 std::cout << " atmos_trken = " << fFitAtmosTrkEn << std::endl;
00394 std::cout << " atmos_trken_exit = " << fFitAtmosTrkEnExit << std::endl;
00395 std::cout << " atmos_shwen = " << fFitAtmosShwEn << std::endl;
00396 }
00397
00398 return;
00399 }
00400
00401 Bool_t TwoFlavourFitter::UsingBeamData()
00402 {
00403 return ( Configuration::Instance()->UsingBeamData()
00404 && TemplateReader::Instance()->TouchBeamData() );
00405 }
00406
00407 Bool_t TwoFlavourFitter::UsingAtmosData()
00408 {
00409 return ( Configuration::Instance()->UsingAtmosData()
00410 && TemplateReader::Instance()->TouchAtmosData() );
00411 }
00412
00413 Double_t TwoFlavourFitter::GetDmsq( Int_t n )
00414 {
00415 if( fDmsqBins<=1 ){
00416 return fDmsqMin;
00417 }
00418 else if( n>=0 && n<fDmsqBins ){
00419 return fDmsqMin + (fDmsqMax-fDmsqMin)*((double)n)/(fDmsqBins-1.0);
00420 }
00421 else{
00422 return 0.0;
00423 }
00424 }
00425
00426 Double_t TwoFlavourFitter::GetSinsq( Int_t n )
00427 {
00428 if( fSinsqBins<=1 ){
00429 return fSinsqMin;
00430 }
00431 else if( n>=0 && n<fSinsqBins ){
00432 return fSinsqMin + (fSinsqMax-fSinsqMin)*((double)n)/(fSinsqBins-1.0);
00433 }
00434 else{
00435 return 0.0;
00436 }
00437 }
00438
00439 void TwoFlavourFitter::BuildGridDmsq()
00440 {
00441 Int_t n = TemplateGrid::Instance()->GetDmsqBins()-1;
00442 Double_t min = TemplateGrid::Instance()->GetDmsq(1);
00443 Double_t max = TemplateGrid::Instance()->GetDmsq(n);
00444
00445 BuildGridDmsq( n, min, max );
00446 }
00447
00448 void TwoFlavourFitter::BuildGridDmsq( Int_t n, Double_t min, Double_t max )
00449 {
00450 std::cout << " *** TwoFlavourFitter::BuildGridDmsq(...) *** " << std::endl;
00451 std::cout << " bins=" << n << " (min=" << min << ", max=" << max << ") " << std::endl;
00452
00453 fDmsqBins = n;
00454 fDmsqMin = min;
00455 fDmsqMax = max;
00456
00457 return;
00458 }
00459
00460 void TwoFlavourFitter::BuildGridSinsq( Int_t n, Double_t min, Double_t max )
00461 {
00462 std::cout << " *** TwoFlavourFitter::BuildGridSinsq(...) *** " << std::endl;
00463 std::cout << " bins=" << n << " (min=" << min << ", max=" << max << ") " << std::endl;
00464
00465 fSinsqBins = n;
00466 fSinsqMin = min;
00467 fSinsqMax = max;
00468
00469 return;
00470 }
00471
00472 void TwoFlavourFitter::BuildGrid( Int_t n, Double_t min, Double_t max )
00473 {
00474 BuildGridDmsq();
00475 BuildGridSinsq( n, min, max );
00476 }
00477
00478 void TwoFlavourFitter::PrintGrid()
00479 {
00480 std::cout << " *** TwoFlavourFitter::PrintGrid() *** " << std::endl;
00481
00482 std::cout << " Dmsq Grid: " << std::endl;
00483 for( Int_t i=0; i<GetDmsqBins(); i++ ){
00484 std::cout << " [" << i << "] " << GetDmsq(i) << std::endl;
00485 }
00486
00487 std::cout << " Sinsq Grid: " << std::endl;
00488 for( Int_t i=0; i<GetSinsqBins(); i++ ){
00489 std::cout << " [" << i << "] " << GetSinsq(i) << std::endl;
00490 }
00491
00492 return;
00493 }
00494
00495 Int_t TwoFlavourFitter::GenSeed()
00496 {
00497 const time_t theTime = time(0);
00498 tm* theTimeStruc = gmtime(&theTime);
00499 int seconds = 60*theTimeStruc->tm_min
00500 + theTimeStruc->tm_sec;
00501 int seed = seconds;
00502 return seed;
00503 }
00504
00505 void TwoFlavourFitter::SetSeed(int x)
00506 {
00507 Int_t seed = x;
00508
00509 if( seed<0 ) seed = this->GenSeed();
00510
00511 gRandom->SetSeed(seed);
00512 }
00513
00514 void TwoFlavourFitter::BuildExpt()
00515 {
00516 if( fMockData ) BuildExptFromMC();
00517 else BuildExptFromData();
00518 }
00519
00520 void TwoFlavourFitter::BuildExptFromData()
00521 {
00522
00523 if( fMockData==1 ) return;
00524
00525 std::cout << " *** TwoFlavourFitter::BuildExptFromData() *** " << std::endl;
00526
00527
00528 this->ResetExpt();
00529
00530
00531 if( TemplateReader::Instance()->TouchData()==0 ) {
00532 std::cout << " <warning> The real data is not loaded " << std::endl;
00533 assert(0);
00534 }
00535
00536
00537 TemplateCalculator::Instance()->GetData( fTemplateExpt );
00538
00539
00540 this->PrintExpt();
00541
00542
00543 assert( TouchExpt() );
00544
00545 return;
00546 }
00547
00548 void TwoFlavourFitter::BuildExptFromMC()
00549 {
00550
00551 if( fMockData==0 ) return;
00552
00553 std::cout << " *** TwoFlavourFitter::BuildExptFromMC() *** " << std::endl;
00554
00555
00556
00557 this->ResetExpt();
00558
00559
00560 if( TemplateReader::Instance()->TouchMC()==0 ) {
00561 std::cout << " <warning> The Monte Carlo is not loaded " << std::endl;
00562 assert(0);
00563 }
00564
00565
00566
00567 Double_t inputDmsq = fInputDmsq;
00568 Double_t inputSinsq = fInputSinsq;
00569 Double_t inputBeamNorm = 0.0;
00570 Double_t inputBeamNCBkg = 0.0;
00571 Double_t inputBeamShwEn = 0.0;
00572 Double_t inputBeamTrkEn = 0.0;
00573 Double_t inputBeamTrkEnExit = 0.0;
00574 Double_t inputAtmosNormCV = 0.0;
00575 Double_t inputAtmosNormRock = 0.0;
00576 Double_t inputAtmosChgCV = 0.0;
00577 Double_t inputAtmosChgRock = 0.0;
00578 Double_t inputAtmosNue = 0.0;
00579 Double_t inputAtmosNCBkg = 0.0;
00580 Double_t inputAtmosSpecNuCV = 0.0;
00581 Double_t inputAtmosSpecNuBarCV = 0.0;
00582 Double_t inputAtmosSpecNuRock = 0.0;
00583 Double_t inputAtmosSpecNuBarRock = 0.0;
00584 Double_t inputAtmosZenith = 0.0;
00585 Double_t inputAtmosTrkEn = 0.0;
00586 Double_t inputAtmosTrkEnExit = 0.0;
00587 Double_t inputAtmosShwEn = 0.0;
00588
00589
00590
00591
00592 if( UsingBeamData() ){
00593
00594 inputBeamNorm = fInputBeamNorm;
00595 inputBeamNCBkg = fInputBeamNCBkg;
00596 inputBeamShwEn = fInputBeamShwEn;
00597 inputBeamTrkEn = fInputBeamTrkEn;
00598 inputBeamTrkEnExit = fInputBeamTrkEnExit;
00599
00600
00601 if( fMockDataWithFluctuations ){
00602
00603 Bool_t flag = 0;
00604 Double_t deltaInputBeamNorm = 0.0;
00605 Double_t deltaInputBeamNCBkg = 0.0;
00606 Double_t deltaInputBeamShwEn = 0.0;
00607 Double_t deltaInputBeamTrkEn = 0.0;
00608 Double_t deltaInputBeamTrkEnExit = 0.0;
00609
00610
00611 flag = 0;
00612 while( !(deltaInputBeamNorm>-2.0 && deltaInputBeamNorm<+2.0 && flag) ){
00613 deltaInputBeamNorm = gRandom->Gaus(0.0,1.0); flag = 1;
00614 }
00615
00616
00617 flag = 0;
00618 while( !(deltaInputBeamNCBkg>-2.0 && deltaInputBeamNCBkg<+2.0 && flag) ){
00619 deltaInputBeamNCBkg = gRandom->Gaus(0.0,1.0); flag = 1;
00620 }
00621
00622
00623 flag = 0;
00624 while( !(deltaInputBeamShwEn>-2.0 && deltaInputBeamShwEn<+2.0 && flag) ){
00625 deltaInputBeamShwEn = gRandom->Gaus(0.0,1.0); flag = 1;
00626 }
00627
00628
00629 flag = 0;
00630 while( !(deltaInputBeamTrkEn>-2.0 && deltaInputBeamTrkEn<+2.0 && flag) ){
00631 deltaInputBeamTrkEn = gRandom->Gaus(0.0,1.0); flag = 1;
00632 }
00633
00634
00635 flag = 0;
00636 while( !(deltaInputBeamTrkEnExit>-2.0 && deltaInputBeamTrkEnExit<+2.0 && flag) ){
00637 deltaInputBeamTrkEnExit = gRandom->Gaus(0.0,1.0); flag = 1;
00638 }
00639
00640
00641 inputBeamNorm += deltaInputBeamNorm;
00642 inputBeamNCBkg += deltaInputBeamNCBkg;
00643 inputBeamShwEn += deltaInputBeamShwEn;
00644 inputBeamTrkEn += deltaInputBeamTrkEn;
00645 inputBeamTrkEnExit += deltaInputBeamTrkEnExit;
00646
00647 }
00648 }
00649
00650
00651
00652 if( UsingAtmosData() ){
00653
00654 inputAtmosNormCV = fInputAtmosNormCV;
00655 inputAtmosNormRock = fInputAtmosNormRock;
00656 inputAtmosChgCV = fInputAtmosChgCV;
00657 inputAtmosChgRock = fInputAtmosChgRock;
00658 inputAtmosNue = fInputAtmosNue;
00659 inputAtmosNCBkg = fInputAtmosNCBkg;
00660 inputAtmosSpecNuCV = fInputAtmosSpecNuCV;
00661 inputAtmosSpecNuBarCV = fInputAtmosSpecNuBarCV;
00662 inputAtmosSpecNuRock = fInputAtmosSpecNuRock;
00663 inputAtmosSpecNuBarRock = fInputAtmosSpecNuBarRock;
00664 inputAtmosZenith = fInputAtmosZenith;
00665 inputAtmosTrkEn = fInputAtmosTrkEn;
00666 inputAtmosTrkEnExit = fInputAtmosTrkEnExit;
00667 inputAtmosShwEn = fInputAtmosShwEn;
00668
00669
00670 if( fMockDataWithFluctuations ){
00671
00672 Bool_t flag = 0;
00673 Double_t deltaInputAtmosNormCV = 0.0;
00674 Double_t deltaInputAtmosNormRock = 0.0;
00675 Double_t deltaInputAtmosChgCV = 0.0;
00676 Double_t deltaInputAtmosChgRock = 0.0;
00677 Double_t deltaInputAtmosNue = 0.0;
00678 Double_t deltaInputAtmosNCBkg = 0.0;
00679 Double_t deltaInputAtmosSpecNuCV = 0.0;
00680 Double_t deltaInputAtmosSpecNuBarCV = 0.0;
00681 Double_t deltaInputAtmosSpecNuRock = 0.0;
00682 Double_t deltaInputAtmosSpecNuBarRock = 0.0;
00683 Double_t deltaInputAtmosZenith = 0.0;
00684 Double_t deltaInputAtmosTrkEn = 0.0;
00685 Double_t deltaInputAtmosTrkEnExit = 0.0;
00686 Double_t deltaInputAtmosShwEn = 0.0;
00687
00688
00689 flag = 0;
00690 while( !(deltaInputAtmosNormCV>-2.0 && deltaInputAtmosNormCV<+2.0 && flag) ){
00691 deltaInputAtmosNormCV = gRandom->Gaus(0.0,1.0); flag = 1;
00692 }
00693
00694
00695 flag = 0;
00696 while( !(deltaInputAtmosNormRock>-2.0 && deltaInputAtmosNormRock<+2.0 && flag) ){
00697 deltaInputAtmosNormRock = gRandom->Gaus(0.0,1.0); flag = 1;
00698 }
00699
00700
00701 flag = 0;
00702 while( !(deltaInputAtmosChgCV>-2.0 && deltaInputAtmosChgCV<+2.0 && flag) ){
00703 deltaInputAtmosChgCV = gRandom->Gaus(0.0,1.0); flag = 1;
00704 }
00705
00706
00707 flag = 0;
00708 while( !(deltaInputAtmosChgRock>-2.0 && deltaInputAtmosChgRock<+2.0 && flag) ){
00709 deltaInputAtmosChgRock = gRandom->Gaus(0.0,1.0); flag = 1;
00710 }
00711
00712
00713 flag = 0;
00714 while( !(deltaInputAtmosNue>-2.0 && deltaInputAtmosNue<+2.0 && flag) ){
00715 deltaInputAtmosNue = gRandom->Gaus(0.0,1.0); flag = 1;
00716 }
00717
00718
00719 flag = 0;
00720 while( !(deltaInputAtmosNCBkg>-2.0 && deltaInputAtmosNCBkg<+2.0 && flag) ){
00721 deltaInputAtmosNCBkg = gRandom->Gaus(0.0,1.0); flag = 1;
00722 }
00723
00724
00725 flag = 0;
00726 while( !(deltaInputAtmosSpecNuCV>-2.0 && deltaInputAtmosSpecNuCV<+2.0 && flag) ){
00727 deltaInputAtmosSpecNuCV = gRandom->Gaus(0.0,1.0); flag = 1;
00728 }
00729
00730
00731 flag = 0;
00732 while( !(deltaInputAtmosSpecNuBarCV>-2.0 && deltaInputAtmosSpecNuBarCV<+2.0 && flag) ){
00733 deltaInputAtmosSpecNuBarCV = gRandom->Gaus(0.0,1.0); flag = 1;
00734 }
00735
00736
00737 flag = 0;
00738 while( !(deltaInputAtmosSpecNuRock>-2.0 && deltaInputAtmosSpecNuRock<+2.0 && flag) ){
00739 deltaInputAtmosSpecNuRock = gRandom->Gaus(0.0,1.0); flag = 1;
00740 }
00741
00742
00743 flag = 0;
00744 while( !(deltaInputAtmosSpecNuBarRock>-2.0 && deltaInputAtmosSpecNuBarRock<+2.0 && flag) ){
00745 deltaInputAtmosSpecNuBarRock = gRandom->Gaus(0.0,1.0); flag = 1;
00746 }
00747
00748
00749 flag = 0;
00750 while( !(deltaInputAtmosZenith>-2.0 && deltaInputAtmosZenith<+2.0 && flag) ){
00751 deltaInputAtmosZenith = gRandom->Gaus(0.0,1.0); flag = 1;
00752 }
00753
00754
00755 flag = 0;
00756 while( !(deltaInputAtmosTrkEn>-2.0 && deltaInputAtmosTrkEn<+2.0 && flag) ){
00757 deltaInputAtmosTrkEn = gRandom->Gaus(0.0,1.0); flag = 1;
00758 }
00759
00760
00761 flag = 0;
00762 while( !(deltaInputAtmosTrkEnExit>-2.0 && deltaInputAtmosTrkEnExit<+2.0 && flag) ){
00763 deltaInputAtmosTrkEnExit = gRandom->Gaus(0.0,1.0); flag = 1;
00764 }
00765
00766
00767 flag = 0;
00768 while( !(deltaInputAtmosShwEn>-2.0 && deltaInputAtmosShwEn<+2.0 && flag) ){
00769 deltaInputAtmosShwEn = gRandom->Gaus(0.0,1.0); flag = 1;
00770 }
00771
00772
00773 inputAtmosNormCV += deltaInputAtmosNormCV;
00774 inputAtmosNormRock += deltaInputAtmosNormRock;
00775 inputAtmosChgCV += deltaInputAtmosChgCV;
00776 inputAtmosChgRock += deltaInputAtmosChgRock;
00777 inputAtmosNue += deltaInputAtmosNue;
00778 inputAtmosNCBkg += deltaInputAtmosNCBkg;
00779 inputAtmosSpecNuCV += deltaInputAtmosSpecNuCV;
00780 inputAtmosSpecNuBarCV += deltaInputAtmosSpecNuBarCV;
00781 inputAtmosSpecNuRock += deltaInputAtmosSpecNuRock;
00782 inputAtmosSpecNuBarRock += deltaInputAtmosSpecNuBarRock;
00783 inputAtmosZenith += deltaInputAtmosZenith;
00784 inputAtmosTrkEn += deltaInputAtmosTrkEn;
00785 inputAtmosTrkEnExit += deltaInputAtmosTrkEnExit;
00786 inputAtmosShwEn += deltaInputAtmosShwEn;
00787 }
00788 }
00789
00790
00791
00792
00793 if( fCorrelatedSystematics ){
00794
00795 if( Configuration::Instance()->GetExperiment()==OscFit::kMINOS ){
00796 if( UsingBeamData() ){
00797 inputBeamTrkEnExit = inputBeamTrkEn;
00798 if( UsingAtmosData() ){
00799 inputAtmosTrkEn = inputBeamTrkEn;
00800 inputAtmosShwEn = 0.7071*(inputAtmosShwEn+inputBeamShwEn);
00801 }
00802 }
00803 }
00804
00805 else if( Configuration::Instance()->GetExperiment()==OscFit::kLBNE ){
00806 if( UsingBeamData() ){
00807 if( UsingAtmosData() ){
00808 inputAtmosTrkEn = inputBeamTrkEn;
00809 inputAtmosTrkEnExit = inputBeamTrkEnExit;
00810 inputAtmosShwEn = inputBeamShwEn;
00811 }
00812 }
00813 }
00814 }
00815
00816
00817
00818
00819 std::cout << " Input Parameters for Expt: " << std::endl;
00820 std::cout << " osc_dmsq = " << inputDmsq << std::endl;
00821 std::cout << " osc_sinsq = " << inputSinsq << std::endl;
00822
00823 if( UsingBeamData() ){
00824 std::cout << " beam_norm = " << inputBeamNorm << std::endl;
00825 std::cout << " beam_ncbkg = " << inputBeamNCBkg << std::endl;
00826 std::cout << " beam_shwen = " << inputBeamShwEn << std::endl;
00827 std::cout << " beam_trken = " << inputBeamTrkEn << std::endl;
00828 std::cout << " beam_trken_exit = " << inputBeamTrkEnExit << std::endl;
00829 }
00830
00831 if( UsingAtmosData() ){
00832 std::cout << " atmos_norm_cv = " << inputAtmosNormCV << std::endl;
00833 std::cout << " atmos_norm_rock = " << inputAtmosNormRock << std::endl;
00834 std::cout << " atmos_chg_cv = " << inputAtmosChgCV << std::endl;
00835 std::cout << " atmos_chg_rock = " << inputAtmosChgRock << std::endl;
00836 std::cout << " atmos_nue = " << inputAtmosNue << std::endl;
00837 std::cout << " atmos_ncbkg = " << inputAtmosNCBkg << std::endl;
00838 std::cout << " atmos_spec_cv_nu = " << inputAtmosSpecNuCV << std::endl;
00839 std::cout << " atmos_spec_cv_nubar = " << inputAtmosSpecNuBarCV << std::endl;
00840 std::cout << " atmos_spec_rock_nu = " << inputAtmosSpecNuRock << std::endl;
00841 std::cout << " atmos_spec_rock_nubar = " << inputAtmosSpecNuBarRock << std::endl;
00842 std::cout << " atmos_zenith = " << inputAtmosZenith << std::endl;
00843 std::cout << " atmos_trken = " << inputAtmosTrkEn << std::endl;
00844 std::cout << " atmos_trken_exit = " << inputAtmosTrkEnExit << std::endl;
00845 std::cout << " atmos_shwen = " << inputAtmosShwEn << std::endl;
00846 }
00847
00848
00849
00850 if( fMockDataWithFluctuations ){
00851 TemplateCalculator::Instance()->GetExperiment( inputDmsq, inputSinsq,
00852 inputBeamNorm, inputBeamNCBkg,
00853 inputBeamShwEn, inputBeamTrkEn,
00854 inputBeamTrkEnExit,
00855 inputAtmosNormCV, inputAtmosNormRock,
00856 inputAtmosChgCV, inputAtmosChgRock,
00857 inputAtmosNue, inputAtmosNCBkg,
00858 inputAtmosSpecNuCV, inputAtmosSpecNuBarCV,
00859 inputAtmosSpecNuRock, inputAtmosSpecNuBarRock,
00860 inputAtmosZenith,
00861 inputAtmosTrkEn, inputAtmosTrkEnExit,
00862 inputAtmosShwEn,
00863 fTemplateExpt );
00864 }
00865
00866
00867 else{
00868 TemplateCalculator::Instance()->GetExpectation( inputDmsq, inputSinsq,
00869 inputBeamNorm, inputBeamNCBkg,
00870 inputBeamShwEn, inputBeamTrkEn,
00871 inputBeamTrkEnExit,
00872 inputAtmosNormCV, inputAtmosNormRock,
00873 inputAtmosChgCV, inputAtmosChgRock,
00874 inputAtmosNue, inputAtmosNCBkg,
00875 inputAtmosSpecNuCV, inputAtmosSpecNuBarCV,
00876 inputAtmosSpecNuRock, inputAtmosSpecNuBarRock,
00877 inputAtmosZenith,
00878 inputAtmosTrkEn, inputAtmosTrkEnExit,
00879 inputAtmosShwEn,
00880 fTemplateExpt);
00881 }
00882
00883
00884 this->PrintExpt();
00885
00886
00887 assert( TouchExpt() );
00888
00889 return;
00890 }
00891
00892 void TwoFlavourFitter::WriteExpt()
00893 {
00894 std::cout << " *** TwoFlavourFitter::WriteExpt() *** " << std::endl;
00895
00896 TString filename = fOutFile;
00897 TString tag = ".expt";
00898
00899 if( filename.EndsWith(".root") ){
00900 filename.Insert(filename.Length()-5,tag);
00901 }
00902 else{
00903 filename.Append(tag);
00904 }
00905
00906 return TemplateWriter::Instance()->Write( fTemplateExpt,
00907 filename.Data() );
00908 }
00909
00910 void TwoFlavourFitter::PrintExpt()
00911 {
00912 std::cout << " *** TwoFlavourFitter::PrintExpt() *** " << std::endl;
00913
00914 std::cout << " Bin Contents: " << std::endl;
00915
00916 for( Int_t irun=0; irun<OscFit::GetNumRuns(); irun++ ){
00917 for( Int_t icont=0; icont<OscFit::GetNumContainments(); icont++ ){
00918 for( Int_t iflav=0; iflav<OscFit::GetNumFlavours(); iflav++ ){
00919 for( Int_t icharge=0; icharge<OscFit::GetNumCharges(); icharge++ ){
00920
00921 Run_t run = OscFit::GetRun(irun);
00922 Containment_t containment = OscFit::GetContainment(icont);
00923 Flavour_t flavour = OscFit::GetFlavour(iflav);
00924 Charge_t charge = OscFit::GetCharge(icharge);
00925
00926 Double_t content = fTemplateExpt->GetBinContent( run, containment,
00927 flavour, charge );
00928
00929 if( content>0.0 ) std::cout << " [" << OscFit::AsString(run) << "][" << OscFit::AsString(containment) << "][" << OscFit::AsString(flavour) << "][" << OscFit::AsString(charge) << "] Total=" << content << std::endl;
00930
00931 }
00932 }
00933 }
00934 }
00935
00936
00937 std::cout << " Total = " << fTemplateExpt->GetTotalContent() << std::endl;
00938 std::cout << " Printing Histograms: " << std::endl;
00939
00940 return fTemplateExpt->PrintHistograms();
00941 }
00942
00943 void TwoFlavourFitter::ResetExpt()
00944 {
00945 return fTemplateExpt->Reset();
00946 }
00947
00948 Bool_t TwoFlavourFitter::TouchExpt()
00949 {
00950 if( fTemplateExpt->GetTotalContent()>0 ) return 1; else return 0;
00951 }
00952
00953 void TwoFlavourFitter::Run( Int_t ifirst, Int_t ilast )
00954 {
00955 this->RunOscFit( ifirst, ilast );
00956 }
00957
00958 void TwoFlavourFitter::RunTemplates()
00959 {
00960 std::cout << " *** TwoFlavourFitter::RunTemplates() *** " << std::endl;
00961
00962
00963
00964 if( this->TouchExpt()==0 ){
00965 this->BuildExpt();
00966 }
00967
00968
00969
00970 TString filename = fOutFile;
00971 TString tag = "";
00972
00973
00974 filename = fOutFile;
00975 tag = ".expt";
00976
00977 if( filename.EndsWith(".root") )
00978 filename.Insert(filename.Length()-5,tag);
00979 else filename.Append(tag);
00980
00981 TemplateWriter::Instance()->Write( fTemplateExpt,
00982 filename.Data() );
00983
00984
00985 filename = fOutFile;
00986 tag = ".noosc";
00987
00988 if( filename.EndsWith(".root") )
00989 filename.Insert(filename.Length()-5,tag);
00990 else filename.Append(tag);
00991
00992 TemplateWriter::Instance()->WriteExpectation( filename.Data() );
00993
00994
00995 filename = fOutFile;
00996 tag = ".ncbkg";
00997
00998 if( filename.EndsWith(".root") )
00999 filename.Insert(filename.Length()-5,tag);
01000 else filename.Append(tag);
01001
01002 TemplateWriter::Instance()->WriteNCBackground( filename.Data() );
01003
01004
01005 filename = fOutFile;
01006 tag = ".cosmic";
01007
01008 if( filename.EndsWith(".root") )
01009 filename.Insert(filename.Length()-5,tag);
01010 else filename.Append(tag);
01011
01012 TemplateWriter::Instance()->WriteCosmicBackground( filename.Data() );
01013
01014 return;
01015 }
01016
01017 void TwoFlavourFitter::RunOscFit( Int_t ifirst, Int_t ilast )
01018 {
01019 std::cout << " *** TwoFlavourFitter::StartOscFit() *** " << std::endl;
01020
01021
01022
01023 this->PrintGrid();
01024
01025
01026
01027 this->PrintSettings();
01028
01029
01030
01031 if( this->TouchExpt()==0 ){
01032 this->BuildExpt();
01033 }
01034
01035
01036
01037 if( fWriteTemplates ){
01038 this->RunTemplates();
01039 }
01040
01041 std::cout << " *** TwoFlavourFitter::RunOscFit() *** " << std::endl;
01042
01043
01044
01045
01046
01047
01048 Int_t nx = GetSinsqBins();
01049 Double_t xlow = GetSinsqMin();
01050 Double_t xhigh = GetSinsqMax();
01051 Double_t dx = (xhigh-xlow)/((double)nx-1.0);
01052 Double_t xmin = xlow - 0.50*dx;
01053 Double_t xmax = xhigh + 0.50*dx;
01054
01055 Int_t ny = GetDmsqBins();
01056 Double_t ylow = GetDmsqMin();
01057 Double_t yhigh = GetDmsqMax();
01058 Double_t dy = (yhigh-ylow)/((double)ny-1.0);
01059 Double_t ymin = ylow - 0.50*dy;
01060 Double_t ymax = yhigh + 0.50*dy;
01061
01062 TH2D* hDmsq = new TH2D("hDmsq", "", nx,xmin,xmax, ny,ymin,ymax);
01063 TH2D* hSinsq = new TH2D("hSinsq", "", nx,xmin,xmax, ny,ymin,ymax);
01064
01065 TH2D* hBeamNorm = new TH2D("hBeamNorm", "", nx,xmin,xmax, ny,ymin,ymax);
01066 hBeamNorm->GetZaxis()->SetTitle("Shift / #sigma");
01067 hBeamNorm->GetYaxis()->SetTitle("#Delta m^{2} / eV^{2}");
01068 hBeamNorm->GetXaxis()->SetTitle("sin^{2}2#theta");
01069 hBeamNorm->GetZaxis()->CenterTitle();
01070 hBeamNorm->GetYaxis()->CenterTitle();
01071 hBeamNorm->GetXaxis()->CenterTitle();
01072
01073 TH2D* hBeamNCBkg = new TH2D("hBeamNCBkg","", nx,xmin,xmax, ny,ymin,ymax);
01074 hBeamNCBkg->GetZaxis()->SetTitle("Shift / #sigma");
01075 hBeamNCBkg->GetYaxis()->SetTitle("#Delta m^{2} / eV^{2}");
01076 hBeamNCBkg->GetXaxis()->SetTitle("sin^{2}2#theta");
01077 hBeamNCBkg->GetZaxis()->CenterTitle();
01078 hBeamNCBkg->GetYaxis()->CenterTitle();
01079 hBeamNCBkg->GetXaxis()->CenterTitle();
01080
01081 TH2D* hBeamShwEn = new TH2D("hBeamShwEn","", nx,xmin,xmax, ny,ymin,ymax);
01082 hBeamShwEn->GetZaxis()->SetTitle("Shift / #sigma");
01083 hBeamShwEn->GetYaxis()->SetTitle("#Delta m^{2} / eV^{2}");
01084 hBeamShwEn->GetXaxis()->SetTitle("sin^{2}2#theta");
01085 hBeamShwEn->GetZaxis()->CenterTitle();
01086 hBeamShwEn->GetYaxis()->CenterTitle();
01087 hBeamShwEn->GetXaxis()->CenterTitle();
01088
01089 TH2D* hBeamTrkEn = new TH2D("hBeamTrkEn","", nx,xmin,xmax, ny,ymin,ymax);
01090 hBeamTrkEn->GetZaxis()->SetTitle("Shift / #sigma");
01091 hBeamTrkEn->GetYaxis()->SetTitle("#Delta m^{2} / eV^{2}");
01092 hBeamTrkEn->GetXaxis()->SetTitle("sin^{2}2#theta");
01093 hBeamTrkEn->GetZaxis()->CenterTitle();
01094 hBeamTrkEn->GetYaxis()->CenterTitle();
01095 hBeamTrkEn->GetXaxis()->CenterTitle();
01096
01097 TH2D* hBeamTrkEnExit = new TH2D("hBeamTrkEnExit","", nx,xmin,xmax, ny,ymin,ymax);
01098 hBeamTrkEnExit->GetZaxis()->SetTitle("Shift / #sigma");
01099 hBeamTrkEnExit->GetYaxis()->SetTitle("#Delta m^{2} / eV^{2}");
01100 hBeamTrkEnExit->GetXaxis()->SetTitle("sin^{2}2#theta");
01101 hBeamTrkEnExit->GetZaxis()->CenterTitle();
01102 hBeamTrkEnExit->GetYaxis()->CenterTitle();
01103 hBeamTrkEnExit->GetXaxis()->CenterTitle();
01104
01105 TH2D* hAtmosNormCV = new TH2D("hAtmosNormCV","", nx,xmin,xmax, ny,ymin,ymax);
01106 hAtmosNormCV->GetZaxis()->SetTitle("Shift / #sigma");
01107 hAtmosNormCV->GetYaxis()->SetTitle("#Delta m^{2} / eV^{2}");
01108 hAtmosNormCV->GetXaxis()->SetTitle("sin^{2}2#theta");
01109 hAtmosNormCV->GetZaxis()->CenterTitle();
01110 hAtmosNormCV->GetYaxis()->CenterTitle();
01111 hAtmosNormCV->GetXaxis()->CenterTitle();
01112
01113 TH2D* hAtmosNormRock = new TH2D("hAtmosNormRock","", nx,xmin,xmax, ny,ymin,ymax);
01114 hAtmosNormRock->GetZaxis()->SetTitle("Shift / #sigma");
01115 hAtmosNormRock->GetYaxis()->SetTitle("#Delta m^{2} / eV^{2}");
01116 hAtmosNormRock->GetXaxis()->SetTitle("sin^{2}2#theta");
01117 hAtmosNormRock->GetZaxis()->CenterTitle();
01118 hAtmosNormRock->GetYaxis()->CenterTitle();
01119 hAtmosNormRock->GetXaxis()->CenterTitle();
01120
01121 TH2D* hAtmosChgCV = new TH2D("hAtmosChgCV","", nx,xmin,xmax, ny,ymin,ymax);
01122 hAtmosChgCV->GetZaxis()->SetTitle("Shift / #sigma");
01123 hAtmosChgCV->GetYaxis()->SetTitle("#Delta m^{2} / eV^{2}");
01124 hAtmosChgCV->GetXaxis()->SetTitle("sin^{2}2#theta");
01125 hAtmosChgCV->GetZaxis()->CenterTitle();
01126 hAtmosChgCV->GetYaxis()->CenterTitle();
01127 hAtmosChgCV->GetXaxis()->CenterTitle();
01128
01129 TH2D* hAtmosChgRock = new TH2D("hAtmosChgRock","", nx,xmin,xmax, ny,ymin,ymax);
01130 hAtmosChgRock->GetZaxis()->SetTitle("Shift / #sigma");
01131 hAtmosChgRock->GetYaxis()->SetTitle("#Delta m^{2} / eV^{2}");
01132 hAtmosChgRock->GetXaxis()->SetTitle("sin^{2}2#theta");
01133 hAtmosChgRock->GetZaxis()->CenterTitle();
01134 hAtmosChgRock->GetYaxis()->CenterTitle();
01135 hAtmosChgRock->GetXaxis()->CenterTitle();
01136
01137 TH2D* hAtmosNue = new TH2D("hAtmosNue","", nx,xmin,xmax, ny,ymin,ymax);
01138 hAtmosNue->GetZaxis()->SetTitle("Shift / #sigma");
01139 hAtmosNue->GetYaxis()->SetTitle("#Delta m^{2} / eV^{2}");
01140 hAtmosNue->GetXaxis()->SetTitle("sin^{2}2#theta");
01141 hAtmosNue->GetZaxis()->CenterTitle();
01142 hAtmosNue->GetYaxis()->CenterTitle();
01143 hAtmosNue->GetXaxis()->CenterTitle();
01144
01145 TH2D* hAtmosNCBkg = new TH2D("hAtmosNCBkg","", nx,xmin,xmax, ny,ymin,ymax);
01146 hAtmosNCBkg->GetZaxis()->SetTitle("Shift / #sigma");
01147 hAtmosNCBkg->GetYaxis()->SetTitle("#Delta m^{2} / eV^{2}");
01148 hAtmosNCBkg->GetXaxis()->SetTitle("sin^{2}2#theta");
01149 hAtmosNCBkg->GetZaxis()->CenterTitle();
01150 hAtmosNCBkg->GetYaxis()->CenterTitle();
01151 hAtmosNCBkg->GetXaxis()->CenterTitle();
01152
01153 TH2D* hAtmosSpecNuCV = new TH2D("hAtmosSpecNuCV","", nx,xmin,xmax, ny,ymin,ymax);
01154 hAtmosSpecNuCV->GetZaxis()->SetTitle("Shift / #sigma");
01155 hAtmosSpecNuCV->GetYaxis()->SetTitle("#Delta m^{2} / eV^{2}");
01156 hAtmosSpecNuCV->GetXaxis()->SetTitle("sin^{2}2#theta");
01157 hAtmosSpecNuCV->GetZaxis()->CenterTitle();
01158 hAtmosSpecNuCV->GetYaxis()->CenterTitle();
01159 hAtmosSpecNuCV->GetXaxis()->CenterTitle();
01160
01161 TH2D* hAtmosSpecNuBarCV = new TH2D("hAtmosSpecNuBarCV","", nx,xmin,xmax, ny,ymin,ymax);
01162 hAtmosSpecNuBarCV->GetZaxis()->SetTitle("Shift / #sigma");
01163 hAtmosSpecNuBarCV->GetYaxis()->SetTitle("#Delta m^{2} / eV^{2}");
01164 hAtmosSpecNuBarCV->GetXaxis()->SetTitle("sin^{2}2#theta");
01165 hAtmosSpecNuBarCV->GetZaxis()->CenterTitle();
01166 hAtmosSpecNuBarCV->GetYaxis()->CenterTitle();
01167 hAtmosSpecNuBarCV->GetXaxis()->CenterTitle();
01168
01169 TH2D* hAtmosSpecNuRock = new TH2D("hAtmosSpecNuRock","", nx,xmin,xmax, ny,ymin,ymax);
01170 hAtmosSpecNuRock->GetZaxis()->SetTitle("Shift / #sigma");
01171 hAtmosSpecNuRock->GetYaxis()->SetTitle("#Delta m^{2} / eV^{2}");
01172 hAtmosSpecNuRock->GetXaxis()->SetTitle("sin^{2}2#theta");
01173 hAtmosSpecNuRock->GetZaxis()->CenterTitle();
01174 hAtmosSpecNuRock->GetYaxis()->CenterTitle();
01175 hAtmosSpecNuRock->GetXaxis()->CenterTitle();
01176
01177 TH2D* hAtmosSpecNuBarRock = new TH2D("hAtmosSpecNuBarRock","", nx,xmin,xmax, ny,ymin,ymax);
01178 hAtmosSpecNuBarRock->GetZaxis()->SetTitle("Shift / #sigma");
01179 hAtmosSpecNuBarRock->GetYaxis()->SetTitle("#Delta m^{2} / eV^{2}");
01180 hAtmosSpecNuBarRock->GetXaxis()->SetTitle("sin^{2}2#theta");
01181 hAtmosSpecNuBarRock->GetZaxis()->CenterTitle();
01182 hAtmosSpecNuBarRock->GetYaxis()->CenterTitle();
01183 hAtmosSpecNuBarRock->GetXaxis()->CenterTitle();
01184
01185 TH2D* hAtmosZenith = new TH2D("hAtmosZenith","", nx,xmin,xmax, ny,ymin,ymax);
01186 hAtmosZenith->GetZaxis()->SetTitle("Shift / #sigma");
01187 hAtmosZenith->GetYaxis()->SetTitle("#Delta m^{2} / eV^{2}");
01188 hAtmosZenith->GetXaxis()->SetTitle("sin^{2}2#theta");
01189 hAtmosZenith->GetZaxis()->CenterTitle();
01190 hAtmosZenith->GetYaxis()->CenterTitle();
01191 hAtmosZenith->GetXaxis()->CenterTitle();
01192
01193 TH2D* hAtmosTrkEn = new TH2D("hAtmosTrkEn","", nx,xmin,xmax, ny,ymin,ymax);
01194 hAtmosTrkEn->GetZaxis()->SetTitle("Shift / #sigma");
01195 hAtmosTrkEn->GetYaxis()->SetTitle("#Delta m^{2} / eV^{2}");
01196 hAtmosTrkEn->GetXaxis()->SetTitle("sin^{2}2#theta");
01197 hAtmosTrkEn->GetZaxis()->CenterTitle();
01198 hAtmosTrkEn->GetYaxis()->CenterTitle();
01199 hAtmosTrkEn->GetXaxis()->CenterTitle();
01200
01201 TH2D* hAtmosTrkEnExit = new TH2D("hAtmosTrkEnExit","", nx,xmin,xmax, ny,ymin,ymax);
01202 hAtmosTrkEnExit->GetZaxis()->SetTitle("Shift / #sigma");
01203 hAtmosTrkEnExit->GetYaxis()->SetTitle("#Delta m^{2} / eV^{2}");
01204 hAtmosTrkEnExit->GetXaxis()->SetTitle("sin^{2}2#theta");
01205 hAtmosTrkEnExit->GetZaxis()->CenterTitle();
01206 hAtmosTrkEnExit->GetYaxis()->CenterTitle();
01207 hAtmosTrkEnExit->GetXaxis()->CenterTitle();
01208
01209 TH2D* hAtmosShwEn = new TH2D("hAtmosShwEn","", nx,xmin,xmax, ny,ymin,ymax);
01210 hAtmosShwEn->GetZaxis()->SetTitle("Shift / #sigma");
01211 hAtmosShwEn->GetYaxis()->SetTitle("#Delta m^{2} / eV^{2}");
01212 hAtmosShwEn->GetXaxis()->SetTitle("sin^{2}2#theta");
01213 hAtmosShwEn->GetZaxis()->CenterTitle();
01214 hAtmosShwEn->GetYaxis()->CenterTitle();
01215 hAtmosShwEn->GetXaxis()->CenterTitle();
01216
01217 TH2D* hLnLTotal = new TH2D("hLnLTotal", "", nx,xmin,xmax, ny,ymin,ymax);
01218 hLnLTotal->GetYaxis()->SetTitle("#Delta m^{2} / eV^{2}");
01219 hLnLTotal->GetXaxis()->SetTitle("sin^{2}2#theta");
01220 hLnLTotal->GetYaxis()->CenterTitle();
01221 hLnLTotal->GetXaxis()->CenterTitle();
01222
01223 TH2D* hLnL = new TH2D("hLnL", "", nx,xmin,xmax, ny,ymin,ymax);
01224 hLnL->GetYaxis()->SetTitle("#Delta m^{2} / eV^{2}");
01225 hLnL->GetXaxis()->SetTitle("sin^{2}2#theta");
01226 hLnL->GetYaxis()->CenterTitle();
01227 hLnL->GetXaxis()->CenterTitle();
01228
01229
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246 Double_t input_dmsq = 0.0;
01247 Double_t input_sinsq = 0.0;
01248
01249 Double_t osc_dmsq = 0.0;
01250 Double_t osc_sinsq = 0.0;
01251 Double_t beam_norm = 0.0;
01252 Double_t beam_ncbkg = 0.0;
01253 Double_t beam_shwen = 0.0;
01254 Double_t beam_trken = 0.0;
01255 Double_t beam_trken_exit = 0.0;
01256 Double_t atmos_norm_cv = 0.0;
01257 Double_t atmos_norm_rock = 0.0;
01258 Double_t atmos_chg_cv = 0.0;
01259 Double_t atmos_chg_rock = 0.0;
01260 Double_t atmos_nue = 0.0;
01261 Double_t atmos_ncbkg = 0.0;
01262 Double_t atmos_spec_cv_nu = 0.0;
01263 Double_t atmos_spec_cv_nubar = 0.0;
01264 Double_t atmos_spec_rock_nu = 0.0;
01265 Double_t atmos_spec_rock_nubar = 0.0;
01266 Double_t atmos_zenith = 0.0;
01267 Double_t atmos_trken = 0.0;
01268 Double_t atmos_trken_exit = 0.0;
01269 Double_t atmos_shwen = 0.0;
01270
01271 Double_t LnL = 0.0;
01272 Double_t minLnL = fBigNumber;
01273
01274 Int_t ibest = -1;
01275 Int_t jbest = -1;
01276
01277 Int_t first = ifirst;
01278 Int_t last = ilast;
01279
01280 if( first<0 ) first = 0;
01281 if( last<0 ) last = GetDmsqBins();
01282 if( last>GetDmsqBins() ) last = GetDmsqBins();
01283
01284
01285 for( Int_t i=first; i<last; i++ ){
01286 std::cout << " ---- RunFit: [" << i << "] [" << first << "|" << last << "] --- " << std::endl;
01287
01288
01289 for( Int_t j=0; j<GetSinsqBins(); j++ ){
01290 std::cout << " [" << i << "] [" << j << "]" << std::endl;
01291
01292
01293 input_dmsq = GetDmsq( i );
01294 input_sinsq = GetSinsq( j );
01295
01296
01297 if( fPseudoThreeFlavours ){
01298 input_sinsq = 4.0*input_sinsq*(1.0-input_sinsq);
01299 }
01300
01301
01302 if( fFitDmsq ) input_dmsq = fSeedDmsq;
01303 if( fFitSinsq ) input_sinsq = fSeedSinsq;
01304
01305
01306 this->RunFit( input_dmsq, input_sinsq,
01307 osc_dmsq, osc_sinsq,
01308 beam_norm, beam_ncbkg,
01309 beam_shwen, beam_trken,
01310 beam_trken_exit,
01311 atmos_norm_cv, atmos_norm_rock,
01312 atmos_chg_cv, atmos_chg_rock,
01313 atmos_nue, atmos_ncbkg,
01314 atmos_spec_cv_nu, atmos_spec_cv_nubar,
01315 atmos_spec_rock_nu, atmos_spec_rock_nubar,
01316 atmos_zenith,
01317 atmos_trken, atmos_trken_exit,
01318 atmos_shwen,
01319 LnL );
01320
01321
01322 hDmsq->SetBinContent( j+1, i+1, osc_dmsq );
01323 hSinsq->SetBinContent( j+1, i+1, osc_sinsq );
01324
01325 hBeamNorm->SetBinContent( j+1, i+1, beam_norm );
01326 hBeamNCBkg->SetBinContent( j+1, i+1, beam_ncbkg );
01327 hBeamShwEn->SetBinContent( j+1, i+1, beam_shwen );
01328 hBeamTrkEn->SetBinContent( j+1, i+1, beam_trken );
01329 hBeamTrkEnExit->SetBinContent( j+1, i+1, beam_trken_exit );
01330
01331 hAtmosNormCV->SetBinContent( j+1, i+1, atmos_norm_cv );
01332 hAtmosNormRock->SetBinContent( j+1, i+1, atmos_norm_rock );
01333 hAtmosChgCV->SetBinContent( j+1, i+1,atmos_chg_cv );
01334 hAtmosChgRock->SetBinContent( j+1, i+1,atmos_chg_rock );
01335 hAtmosNue->SetBinContent( j+1, i+1, atmos_nue );
01336 hAtmosNCBkg->SetBinContent( j+1, i+1, atmos_ncbkg );
01337 hAtmosSpecNuCV->SetBinContent( j+1, i+1,atmos_spec_cv_nu );
01338 hAtmosSpecNuBarCV->SetBinContent( j+1, i+1, atmos_spec_cv_nubar );
01339 hAtmosSpecNuRock->SetBinContent( j+1, i+1,atmos_spec_rock_nu );
01340 hAtmosSpecNuBarRock->SetBinContent( j+1, i+1, atmos_spec_rock_nubar );
01341 hAtmosZenith->SetBinContent( j+1, i+1, atmos_zenith );
01342 hAtmosTrkEn->SetBinContent( j+1, i+1, atmos_trken );
01343 hAtmosTrkEnExit->SetBinContent( j+1, i+1, atmos_trken_exit );
01344 hAtmosShwEn->SetBinContent( j+1, i+1,atmos_shwen );
01345
01346 hLnL->SetBinContent( j+1, i+1, LnL );
01347 hLnLTotal->SetBinContent( j+1, i+1, LnL );
01348
01349
01350 if( LnL<minLnL ){
01351 ibest = i;
01352 jbest = j;
01353 minLnL = LnL;
01354 }
01355
01356 }
01357 }
01358
01359
01360
01361
01362 if( ibest>=0 && jbest>=0 ){
01363
01364
01365
01366
01367
01368
01369
01370
01371
01372
01373 for( Int_t i=0; i<hLnL->GetYaxis()->GetNbins(); i++ ){
01374 for( Int_t j=0; j<hLnL->GetXaxis()->GetNbins(); j++ ){
01375 Double_t myLnL = hLnL->GetBinContent( j+1, i+1 );
01376 hLnL->SetBinContent( j+1, i+1, myLnL - minLnL );
01377 }
01378 }
01379
01380 osc_dmsq = hDmsq->GetBinContent( jbest+1, ibest+1 );
01381 osc_sinsq = hSinsq->GetBinContent( jbest+1, ibest+1 );
01382
01383 beam_norm = hBeamNorm->GetBinContent( jbest+1, ibest+1 );
01384 beam_ncbkg = hBeamNCBkg->GetBinContent( jbest+1, ibest+1 );
01385 beam_shwen = hBeamShwEn->GetBinContent( jbest+1, ibest+1 );
01386 beam_trken = hBeamTrkEn->GetBinContent( jbest+1, ibest+1 );
01387 beam_trken_exit = hBeamTrkEnExit->GetBinContent( jbest+1, ibest+1 );
01388
01389 atmos_norm_cv = hAtmosNormCV->GetBinContent( jbest+1, ibest+1 );
01390 atmos_norm_rock = hAtmosNormRock->GetBinContent( jbest+1, ibest+1 );
01391 atmos_chg_cv = hAtmosChgCV->GetBinContent( jbest+1, ibest+1 );
01392 atmos_chg_rock = hAtmosChgRock->GetBinContent( jbest+1, ibest+1 );
01393 atmos_nue = hAtmosNue->GetBinContent( jbest+1, ibest+1 );
01394 atmos_ncbkg = hAtmosNCBkg->GetBinContent( jbest+1, ibest+1 );
01395 atmos_spec_cv_nu = hAtmosSpecNuCV->GetBinContent( jbest+1, ibest+1 );
01396 atmos_spec_cv_nubar = hAtmosSpecNuBarCV->GetBinContent( jbest+1, ibest+1 );
01397 atmos_spec_rock_nu = hAtmosSpecNuRock->GetBinContent( jbest+1, ibest+1 );
01398 atmos_spec_rock_nubar = hAtmosSpecNuBarRock->GetBinContent( jbest+1, ibest+1 );
01399 atmos_zenith = hAtmosZenith->GetBinContent( jbest+1, ibest+1 );
01400 atmos_trken = hAtmosTrkEn->GetBinContent( jbest+1, ibest+1 );
01401 atmos_trken_exit = hAtmosTrkEnExit->GetBinContent( jbest+1, ibest+1 );
01402 atmos_shwen = hAtmosShwEn->GetBinContent( jbest+1, ibest+1 );
01403
01404 std::cout << " *** TwoFlavourFitter::PrintResults() *** " << std::endl;
01405 std::cout << " OSC FIT RESULTS: " << std::endl;
01406
01407 std::cout << " seed parameters: " << std::endl;
01408 std::cout << " osc_dmsq = " << fSeedDmsq << std::endl;
01409 std::cout << " osc_sinsq = " << fSeedSinsq << std::endl;
01410 if( UsingBeamData() ){
01411 std::cout << " beam_norm = " << fSeedBeamNorm << std::endl;
01412 std::cout << " beam_ncbkg = " << fSeedBeamNCBkg << std::endl;
01413 std::cout << " beam_shwen = " << fSeedBeamShwEn << std::endl;
01414 std::cout << " beam_trken = " << fSeedBeamTrkEn << std::endl;
01415 std::cout << " beam_trken_exit = " << fSeedBeamTrkEnExit << std::endl;
01416 }
01417 if( UsingAtmosData() ){
01418 std::cout << " atmos_norm_cv = " << fSeedAtmosNormCV << std::endl;
01419 std::cout << " atmos_norm_rock = " << fSeedAtmosNormRock << std::endl;
01420 std::cout << " atmos_chg_cv = " << fSeedAtmosChgCV << std::endl;
01421 std::cout << " atmos_chg_rock = " << fSeedAtmosChgRock << std::endl;
01422 std::cout << " atmos_nue = " << fSeedAtmosNue << std::endl;
01423 std::cout << " atmos_ncbkg = " << fSeedAtmosNCBkg << std::endl;
01424 std::cout << " atmos_spec_cv_nu = " << fSeedAtmosSpecNuCV << std::endl;
01425 std::cout << " atmos_spec_cv_nubar = " << fSeedAtmosSpecNuBarCV << std::endl;
01426 std::cout << " atmos_spec_rock_nu = " << fSeedAtmosSpecNuRock << std::endl;
01427 std::cout << " atmos_spec_rock_nubar = " << fSeedAtmosSpecNuBarRock << std::endl;
01428 std::cout << " atmos_zenith = " << fSeedAtmosZenith << std::endl;
01429 std::cout << " atmos_trken = " << fSeedAtmosTrkEn << std::endl;
01430 std::cout << " atmos_trken_exit = " << fSeedAtmosTrkEnExit << std::endl;
01431 std::cout << " atmos_shwen = " << fSeedAtmosShwEn << std::endl;
01432 }
01433
01434 std::cout << " best fit parameters: " << std::endl;
01435 std::cout << " osc_dmsq = " << osc_dmsq << std::endl;
01436 std::cout << " osc_sinsq = " << osc_sinsq << std::endl;
01437 if( UsingBeamData() ){
01438 std::cout << " beam_norm = " << beam_norm << std::endl;
01439 std::cout << " beam_ncbkg = " << beam_ncbkg << std::endl;
01440 std::cout << " beam_shwen = " << beam_shwen << std::endl;
01441 std::cout << " beam_trken = " << beam_trken << std::endl;
01442 std::cout << " beam_trken_exit = " << beam_trken_exit << std::endl;
01443 }
01444 if( UsingAtmosData() ){
01445 std::cout << " atmos_norm_cv = " << atmos_norm_cv << std::endl;
01446 std::cout << " atmos_norm_rock = " << atmos_norm_rock << std::endl;
01447 std::cout << " atmos_chg_cv = " << atmos_chg_cv << std::endl;
01448 std::cout << " atmos_chg_rock = " << atmos_chg_rock << std::endl;
01449 std::cout << " atmos_nue = " << atmos_nue << std::endl;
01450 std::cout << " atmos_ncbkg = " << atmos_ncbkg << std::endl;
01451 std::cout << " atmos_spec_cv_nu = " << atmos_spec_cv_nu << std::endl;
01452 std::cout << " atmos_spec_cv_nubar = " << atmos_spec_cv_nubar << std::endl;
01453 std::cout << " atmos_spec_rock_nu = " << atmos_spec_rock_nu << std::endl;
01454 std::cout << " atmos_spec_rock_nubar = " << atmos_spec_rock_nubar << std::endl;
01455 std::cout << " atmos_zenith = " << atmos_zenith << std::endl;
01456 std::cout << " atmos_trken = " << atmos_trken << std::endl;
01457 std::cout << " atmos_trken_exit = " << atmos_trken_exit << std::endl;
01458 std::cout << " atmos_shwen = " << atmos_shwen << std::endl;
01459 }
01460
01461 std::cout << " goodness of fit: " << std::endl;
01462 std::cout << " LnLbest=" << minLnL << std::endl;
01463
01464
01465 if( fWriteTemplates ){
01466
01467 TString filename = fOutFile;
01468 TString tag = ".bestfit";
01469
01470 if( filename.EndsWith(".root") )
01471 filename.Insert(filename.Length()-5,tag);
01472 else filename.Append(tag);
01473
01474 TemplateWriter::Instance()->WriteExpectation( osc_dmsq, osc_sinsq,
01475 beam_norm, beam_ncbkg,
01476 beam_shwen, beam_trken,
01477 beam_trken_exit,
01478 atmos_norm_cv, atmos_norm_rock,
01479 atmos_chg_cv, atmos_chg_rock,
01480 atmos_nue, atmos_ncbkg,
01481 atmos_spec_cv_nu, atmos_spec_cv_nubar,
01482 atmos_spec_rock_nu, atmos_spec_rock_nubar,
01483 atmos_zenith,
01484 atmos_trken, atmos_trken_exit,
01485 atmos_shwen,
01486 filename.Data() );
01487 }
01488 }
01489
01490
01491
01492
01493 TH2D* hLnL68 = (TH2D*)(hLnL->Clone("hLnL68"));
01494 hLnL68->SetContour(1);
01495 hLnL68->SetContourLevel(0,fLnL68);
01496
01497 TH2D* hLnL90 = (TH2D*)(hLnL->Clone("hLnL90"));
01498 hLnL90->SetContour(1);
01499 hLnL90->SetContourLevel(0,fLnL90);
01500
01501 TH2D* hLnL99 = (TH2D*)(hLnL->Clone("hLnL99"));
01502 hLnL99->SetContour(1);
01503 hLnL99->SetContourLevel(0,fLnL99);
01504
01505
01506
01507
01508 if( 1
01509 || fWriteContours ){
01510 TString filename = fOutFile;
01511 TString tag = ".output";
01512
01513 if( filename.EndsWith(".root") )
01514 filename.Insert(filename.Length()-5,tag);
01515 else filename.Append(tag);
01516
01517 std::cout << " *** TwoFlavourFitter::WriteContours() *** " << std::endl;
01518 std::cout << " Writing Contours to File: " << filename.Data() << std::endl;
01519
01520 TDirectory* tmpd = gDirectory;
01521 TFile* file = new TFile(filename.Data(),"recreate");
01522
01523 hDmsq->Write();
01524 hSinsq->Write();
01525 hBeamNorm->Write();
01526 hBeamNCBkg->Write();
01527 hBeamShwEn->Write();
01528 hBeamTrkEn->Write();
01529 hBeamTrkEnExit->Write();
01530 hAtmosNormCV->Write();
01531 hAtmosNormRock->Write();
01532 hAtmosChgCV->Write();
01533 hAtmosChgRock->Write();
01534 hAtmosNue->Write();
01535 hAtmosNCBkg->Write();
01536 hAtmosSpecNuCV->Write();
01537 hAtmosSpecNuBarCV->Write();
01538 hAtmosSpecNuRock->Write();
01539 hAtmosSpecNuBarRock->Write();
01540 hAtmosZenith->Write();
01541 hAtmosTrkEn->Write();
01542 hAtmosTrkEnExit->Write();
01543 hAtmosShwEn->Write();
01544
01545 hLnL->Write();
01546 hLnLTotal->Write();
01547 hLnL68->Write();
01548 hLnL90->Write();
01549 hLnL99->Write();
01550
01551 file->Close();
01552 tmpd->cd();
01553 }
01554
01555 return;
01556 }
01557
01558 void TwoFlavourFitter::RunFit()
01559 {
01560 std::cout << " *** TwoFlavourFitter::StartFit() *** " << std::endl;
01561
01562
01563
01564 this->PrintSettings();
01565
01566
01567
01568 if( this->TouchExpt()==0 ){
01569 this->BuildExpt();
01570 }
01571
01572
01573
01574 if( fWriteTemplates ){
01575 this->RunTemplates();
01576 }
01577
01578
01579
01580 Double_t input_dmsq = fSeedDmsq;
01581 Double_t input_sinsq = fSeedSinsq;
01582
01583 Double_t osc_dmsq = 0.0;
01584 Double_t osc_sinsq = 0.0;
01585 Double_t beam_norm = 0.0;
01586 Double_t beam_ncbkg = 0.0;
01587 Double_t beam_shwen = 0.0;
01588 Double_t beam_trken = 0.0;
01589 Double_t beam_trken_exit = 0.0;
01590 Double_t atmos_norm_cv = 0.0;
01591 Double_t atmos_norm_rock = 0.0;
01592 Double_t atmos_chg_cv = 0.0;
01593 Double_t atmos_chg_rock = 0.0;
01594 Double_t atmos_nue = 0.0;
01595 Double_t atmos_ncbkg = 0.0;
01596 Double_t atmos_spec_cv_nu = 0.0;
01597 Double_t atmos_spec_cv_nubar = 0.0;
01598 Double_t atmos_spec_rock_nu = 0.0;
01599 Double_t atmos_spec_rock_nubar = 0.0;
01600 Double_t atmos_zenith = 0.0;
01601 Double_t atmos_trken = 0.0;
01602 Double_t atmos_trken_exit = 0.0;
01603 Double_t atmos_shwen = 0.0;
01604
01605 Double_t LnL = 0.0;
01606
01607 this->RunFit( input_dmsq, input_sinsq,
01608 osc_dmsq, osc_sinsq,
01609 beam_norm, beam_ncbkg,
01610 beam_shwen, beam_trken,
01611 beam_trken_exit,
01612 atmos_norm_cv, atmos_norm_rock,
01613 atmos_chg_cv, atmos_chg_rock,
01614 atmos_nue, atmos_ncbkg,
01615 atmos_spec_cv_nu, atmos_spec_cv_nubar,
01616 atmos_spec_rock_nu, atmos_spec_rock_nubar,
01617 atmos_zenith,
01618 atmos_trken, atmos_trken_exit,
01619 atmos_shwen,
01620 LnL );
01621
01622 std::cout << " *** TwoFlavourFitter::PrintResults() *** " << std::endl;
01623 std::cout << " OSC FIT RESULTS: " << std::endl;
01624
01625 std::cout << " seed parameters: " << std::endl;
01626 std::cout << " osc_dmsq = " << fSeedDmsq << std::endl;
01627 std::cout << " osc_sinsq = " << fSeedSinsq << std::endl;
01628 if( UsingBeamData() ){
01629 std::cout << " beam_norm = " << fSeedBeamNorm << std::endl;
01630 std::cout << " beam_ncbkg = " << fSeedBeamNCBkg << std::endl;
01631 std::cout << " beam_shwen = " << fSeedBeamShwEn << std::endl;
01632 std::cout << " beam_trken = " << fSeedBeamTrkEn << std::endl;
01633 std::cout << " beam_trken_exit = " << fSeedBeamTrkEnExit << std::endl;
01634 }
01635 if( UsingAtmosData() ){
01636 std::cout << " atmos_norm_cv = " << fSeedAtmosNormCV << std::endl;
01637 std::cout << " atmos_norm_rock = " << fSeedAtmosNormRock << std::endl;
01638 std::cout << " atmos_chg_cv = " << fSeedAtmosChgCV << std::endl;
01639 std::cout << " atmos_chg_rock = " << fSeedAtmosChgRock << std::endl;
01640 std::cout << " atmos_nue = " << fSeedAtmosNue << std::endl;
01641 std::cout << " atmos_ncbkg = " << fSeedAtmosNCBkg << std::endl;
01642 std::cout << " atmos_spec_cv_nu = " << fSeedAtmosSpecNuCV << std::endl;
01643 std::cout << " atmos_spec_cv_nubar = " << fSeedAtmosSpecNuBarCV << std::endl;
01644 std::cout << " atmos_spec_rock_nu = " << fSeedAtmosSpecNuRock << std::endl;
01645 std::cout << " atmos_spec_rock_nubar = " << fSeedAtmosSpecNuBarRock << std::endl;
01646 std::cout << " atmos_zenith = " << fSeedAtmosZenith << std::endl;
01647 std::cout << " atmos_trken = " << fSeedAtmosTrkEn << std::endl;
01648 std::cout << " atmos_trken_exit = " << fSeedAtmosTrkEnExit << std::endl;
01649 std::cout << " atmos_shwen = " << fSeedAtmosShwEn << std::endl;
01650 }
01651
01652 std::cout << " best fit parameters: " << std::endl;
01653 std::cout << " osc_dmsq = " << osc_dmsq << std::endl;
01654 std::cout << " osc_sinsq = " << osc_sinsq << std::endl;
01655 if( UsingBeamData() ){
01656 std::cout << " beam_norm = " << beam_norm << std::endl;
01657 std::cout << " beam_ncbkg = " << beam_ncbkg << std::endl;
01658 std::cout << " beam_shwen = " << beam_shwen << std::endl;
01659 std::cout << " beam_trken = " << beam_trken << std::endl;
01660 std::cout << " beam_trken_exit = " << beam_trken_exit << std::endl;
01661 }
01662 if( UsingAtmosData() ){
01663 std::cout << " atmos_norm_cv = " << atmos_norm_cv << std::endl;
01664 std::cout << " atmos_norm_rock = " << atmos_norm_rock << std::endl;
01665 std::cout << " atmos_chg_cv = " << atmos_chg_cv << std::endl;
01666 std::cout << " atmos_chg_rock = " << atmos_chg_rock << std::endl;
01667 std::cout << " atmos_nue = " << atmos_nue << std::endl;
01668 std::cout << " atmos_ncbkg = " << atmos_ncbkg << std::endl;
01669 std::cout << " atmos_spec_cv_nu = " << atmos_spec_cv_nu << std::endl;
01670 std::cout << " atmos_spec_cv_nubar = " << atmos_spec_cv_nubar << std::endl;
01671 std::cout << " atmos_spec_rock_nu = " << atmos_spec_rock_nu << std::endl;
01672 std::cout << " atmos_spec_rock_nubar = " << atmos_spec_rock_nubar << std::endl;
01673 std::cout << " atmos_zenith = " << atmos_zenith << std::endl;
01674 std::cout << " atmos_trken = " << atmos_trken << std::endl;
01675 std::cout << " atmos_trken_exit = " << atmos_trken_exit << std::endl;
01676 std::cout << " atmos_shwen = " << atmos_shwen << std::endl;
01677 }
01678
01679 std::cout << " goodness of fit: " << std::endl;
01680 std::cout << " LnLbest=" << LnL << std::endl;
01681
01682
01683
01684 if( fWriteTemplates ){
01685 TString filename = fOutFile;
01686 TString tag = ".bestfit";
01687
01688 if( filename.EndsWith(".root") )
01689 filename.Insert(filename.Length()-5,tag);
01690 else filename.Append(tag);
01691
01692 TemplateWriter::Instance()->WriteExpectation( osc_dmsq, osc_sinsq,
01693 beam_norm, beam_ncbkg,
01694 beam_shwen, beam_trken,
01695 beam_trken_exit,
01696 atmos_norm_cv, atmos_norm_rock,
01697 atmos_chg_cv, atmos_chg_rock,
01698 atmos_nue, atmos_ncbkg,
01699 atmos_spec_cv_nu, atmos_spec_cv_nubar,
01700 atmos_spec_rock_nu, atmos_spec_rock_nubar,
01701 atmos_zenith,
01702 atmos_trken, atmos_trken_exit,
01703 atmos_shwen,
01704 filename.Data() );
01705 }
01706
01707
01708
01709 if( fWriteContours ) this->RunContours();
01710
01711 return;
01712 }
01713
01714 void TwoFlavourFitter::RunPseudoExpts( Int_t numExpts )
01715 {
01716 std::cout << " *** TwoFlavourFitter::RunPseudoExpts() *** " << std::endl;
01717
01718
01719
01720 this->PrintSettings();
01721
01722
01723
01724
01725 Double_t input_dmsq = fSeedDmsq;
01726 Double_t input_sinsq = fSeedSinsq;
01727 Double_t osc_dmsq = 0.0;
01728 Double_t osc_sinsq = 0.0;
01729 Double_t beam_norm = 0.0;
01730 Double_t beam_ncbkg = 0.0;
01731 Double_t beam_shwen = 0.0;
01732 Double_t beam_trken = 0.0;
01733 Double_t beam_trken_exit = 0.0;
01734 Double_t atmos_norm_cv = 0.0;
01735 Double_t atmos_norm_rock = 0.0;
01736 Double_t atmos_chg_cv = 0.0;
01737 Double_t atmos_chg_rock = 0.0;
01738 Double_t atmos_nue = 0.0;
01739 Double_t atmos_ncbkg = 0.0;
01740 Double_t atmos_spec_cv_nu = 0.0;
01741 Double_t atmos_spec_cv_nubar = 0.0;
01742 Double_t atmos_spec_rock_nu = 0.0;
01743 Double_t atmos_spec_rock_nubar = 0.0;
01744 Double_t atmos_zenith = 0.0;
01745 Double_t atmos_trken = 0.0;
01746 Double_t atmos_trken_exit = 0.0;
01747 Double_t atmos_shwen = 0.0;
01748 Double_t LnL = 0.0;
01749
01750
01751
01752
01753 TString filename = fOutFile;
01754 TDirectory* tmpd = 0;
01755
01756 tmpd = gDirectory;
01757 TFile* file = new TFile(filename.Data(),"recreate");
01758 TTree* tree = new TTree("PseudoExpts","fits to pseudo-experiments");
01759 tree->Branch("input_dmsq",&input_dmsq,"input_dmsq/D");
01760 tree->Branch("input_sinsq",&input_sinsq,"input_sinsq/D");
01761 tree->Branch("osc_dmsq",&osc_dmsq,"osc_dmsq/D");
01762 tree->Branch("osc_sinsq",&osc_sinsq,"osc_sinsq/D");
01763 tree->Branch("beam_norm",&beam_norm,"beam_norm/D");
01764 tree->Branch("beam_ncbkg",&beam_ncbkg,"beam_ncbkg/D");
01765 tree->Branch("beam_shwen",&beam_shwen,"beam_shwen/D");
01766 tree->Branch("beam_trken",&beam_trken,"beam_trken/D");
01767 tree->Branch("beam_trken_exit",&beam_trken_exit,"beam_trken_exit/D");
01768 tree->Branch("atmos_norm_cv",&atmos_norm_cv,"atmos_norm_cv/D");
01769 tree->Branch("atmos_norm_rock",&atmos_norm_rock,"atmos_norm_rock/D");
01770 tree->Branch("atmos_chg_cv",&atmos_chg_cv,"atmos_chg_cv/D");
01771 tree->Branch("atmos_chg_rock",&atmos_chg_rock,"atmos_chg_rock/D");
01772 tree->Branch("atmos_nue",&atmos_nue,"atmos_nue/D");
01773 tree->Branch("atmos_ncbkg",&atmos_ncbkg,"atmos_ncbkg/D");
01774 tree->Branch("atmos_spec_cv_nu",&atmos_spec_cv_nu,"atmos_spec_cv_nu/D");
01775 tree->Branch("atmos_spec_cv_nubar",&atmos_spec_cv_nubar,"atmos_spec_cv_nubar/D");
01776 tree->Branch("atmos_spec_rock_nu",&atmos_spec_rock_nu,"atmos_spec_rock_nu/D");
01777 tree->Branch("atmos_spec_rock_nubar",&atmos_spec_rock_nubar,"atmos_spec_rock_nubar/D");
01778 tree->Branch("atmos_zenith",&atmos_zenith,"atmos_zenith/D");
01779 tree->Branch("atmos_trken",&atmos_trken,"atmos_trken/D");
01780 tree->Branch("atmos_trken_exit",&atmos_trken_exit,"atmos_trken_exit/D");
01781 tree->Branch("atmos_shwen",&atmos_shwen,"atmos_shwen/D");
01782 tree->Branch("LnL",&LnL,"LnL/D");
01783 gDirectory = tmpd;
01784
01785
01786
01787
01788 for( Int_t n=0; n<numExpts; n++ ){
01789 std::cout << " --- NEW EXPERIMENT [" << n << "] --- " << std::endl;
01790
01791 this->BuildExpt();
01792
01793 this->RunFit( input_dmsq, input_sinsq,
01794 osc_dmsq, osc_sinsq,
01795 beam_norm, beam_ncbkg,
01796 beam_shwen, beam_trken,
01797 beam_trken_exit,
01798 atmos_norm_cv, atmos_norm_rock,
01799 atmos_chg_cv, atmos_chg_rock,
01800 atmos_nue, atmos_ncbkg,
01801 atmos_spec_cv_nu, atmos_spec_cv_nubar,
01802 atmos_spec_rock_nu, atmos_spec_rock_nubar,
01803 atmos_zenith,
01804 atmos_trken, atmos_trken_exit,
01805 atmos_shwen,
01806 LnL );
01807
01808 std::cout << " --- RESULTS: dmsq=" << osc_dmsq << " sinsq=" << osc_sinsq << " LnL=" << LnL << std::endl;
01809
01810 tmpd = gDirectory;
01811 file->cd();
01812 tree->Fill();
01813 gDirectory = tmpd;
01814 }
01815
01816
01817
01818
01819 tmpd = gDirectory;
01820 file->cd();
01821 tree->Write();
01822 file->Close();
01823 gDirectory = tmpd;
01824
01825 return;
01826 }
01827
01828 void TwoFlavourFitter::RunFit( Double_t input_dmsq, Double_t input_sinsq, Double_t& osc_dmsq, Double_t& osc_sinsq, Double_t& beam_norm, Double_t& beam_ncbkg, Double_t& beam_shwen, Double_t& beam_trken, Double_t& beam_trken_exit, Double_t& atmos_norm_cv, Double_t& atmos_norm_rock, Double_t& atmos_chg_cv, Double_t& atmos_chg_rock, Double_t& atmos_nue, Double_t& atmos_ncbkg, Double_t& atmos_spec_cv_nu, Double_t& atmos_spec_cv_nubar, Double_t& atmos_spec_rock_nu, Double_t& atmos_spec_rock_nubar, Double_t& atmos_zenith, Double_t& atmos_trken, Double_t& atmos_trken_exit, Double_t& atmos_shwen, Double_t& LnL )
01829 {
01830
01831
01832
01833 fOscDmsq = input_dmsq;
01834 fOscSinsq = input_sinsq;
01835
01836 fBeamNorm = fSeedBeamNorm;
01837 fBeamNCBkg = fSeedBeamNCBkg;
01838 fBeamShwEn = fSeedBeamShwEn;
01839 fBeamTrkEn = fSeedBeamTrkEn;
01840 fBeamTrkEnExit = fSeedBeamTrkEnExit;
01841 fAtmosNormCV = fSeedAtmosNormCV;
01842 fAtmosNormRock = fSeedAtmosNormRock;
01843 fAtmosChgCV = fSeedAtmosChgCV;
01844 fAtmosChgRock = fSeedAtmosChgRock;
01845 fAtmosNue = fSeedAtmosNue;
01846 fAtmosNCBkg = fSeedAtmosNCBkg;
01847 fAtmosSpecNuCV = fSeedAtmosSpecNuCV;
01848 fAtmosSpecNuBarCV = fSeedAtmosSpecNuBarCV;
01849 fAtmosSpecNuRock = fSeedAtmosSpecNuRock;
01850 fAtmosSpecNuBarRock = fSeedAtmosSpecNuBarRock;
01851 fAtmosZenith = fSeedAtmosZenith;
01852 fAtmosTrkEn = fSeedAtmosTrkEn;
01853 fAtmosTrkEnExit = fSeedAtmosTrkEnExit;
01854 fAtmosShwEn = fSeedAtmosShwEn;
01855
01856 RunFitWithMinuitInSteps();
01857
01858 osc_dmsq = fOscDmsq;
01859 osc_sinsq = fOscSinsq;
01860 beam_norm = fBeamNorm;
01861 beam_ncbkg = fBeamNCBkg;
01862 beam_shwen = fBeamShwEn;
01863 beam_trken = fBeamTrkEn;
01864 beam_trken_exit = fBeamTrkEnExit;
01865 atmos_norm_cv = fAtmosNormCV;
01866 atmos_norm_rock = fAtmosNormRock;
01867 atmos_chg_cv = fAtmosChgCV;
01868 atmos_chg_rock = fAtmosChgRock;
01869 atmos_nue = fAtmosNue;
01870 atmos_ncbkg = fAtmosNCBkg;
01871 atmos_spec_cv_nu = fAtmosSpecNuCV;
01872 atmos_spec_cv_nubar = fAtmosSpecNuBarCV;
01873 atmos_spec_rock_nu = fAtmosSpecNuRock;
01874 atmos_spec_rock_nubar = fAtmosSpecNuBarRock;
01875 atmos_zenith = fAtmosZenith;
01876 atmos_trken = fAtmosTrkEn;
01877 atmos_trken_exit = fAtmosTrkEnExit;
01878 atmos_shwen = fAtmosShwEn;
01879
01880 LnL = fOutputLnL;
01881
01882 return;
01883 }
01884
01885 void TwoFlavourFitter::RunFitWithMinuitInSteps()
01886 {
01887 if( fDebug ){ std::cout << " *** TwoFlavourFitter::RunFitWithMinuitInSteps() *** " << std::endl; }
01888
01889
01890 fMinuitStatus = 999;
01891
01892
01893 fMaskDmsq = 1;
01894 fMaskSinsq = 1;
01895 fMaskBeamNorm = 1;
01896 fMaskBeamNCBkg = 0;
01897 fMaskBeamShwEn = 0;
01898 fMaskBeamTrkEn = 0;
01899 fMaskBeamTrkEnExit = 0;
01900 fMaskAtmosNormCV = 1;
01901 fMaskAtmosNormRock = 1;
01902 fMaskAtmosChgCV = 0;
01903 fMaskAtmosChgRock = 0;
01904 fMaskAtmosNue = 0;
01905 fMaskAtmosNCBkg = 0;
01906 fMaskAtmosSpecNuCV = 0;
01907 fMaskAtmosSpecNuBarCV = 0;
01908 fMaskAtmosSpecNuRock = 0;
01909 fMaskAtmosSpecNuBarRock = 0;
01910 fMaskAtmosZenith = 0;
01911 fMaskAtmosTrkEn = 0;
01912 fMaskAtmosTrkEnExit = 0;
01913 fMaskAtmosShwEn = 0;
01914
01915 if( fDebug ){ std::cout << " -*- oscillations (simplex) -*- " << std::endl; }
01916 fUseSimplex = 1; RunFitWithMinuit();
01917
01918 if( fDebug ){ std::cout << " -*- oscillations (migrad) -*- " << std::endl; }
01919 fUseSimplex = 0; RunFitWithMinuit();
01920
01921
01922
01923 fMaskDmsq = 0;
01924 fMaskSinsq = 0;
01925 fMaskBeamNorm = 1;
01926 fMaskBeamNCBkg = 1;
01927 fMaskBeamShwEn = 1;
01928 fMaskBeamTrkEn = 1;
01929 fMaskBeamTrkEnExit = 1;
01930 fMaskAtmosNormCV = 1;
01931 fMaskAtmosNormRock = 1;
01932 fMaskAtmosChgCV = 1;
01933 fMaskAtmosChgRock = 1;
01934 fMaskAtmosNue = 0;
01935 fMaskAtmosNCBkg = 0;
01936 fMaskAtmosSpecNuCV = 1;
01937 fMaskAtmosSpecNuBarCV = 1;
01938 fMaskAtmosSpecNuRock = 1;
01939 fMaskAtmosSpecNuBarRock = 1;
01940 fMaskAtmosZenith = 0;
01941 fMaskAtmosTrkEn = 0;
01942 fMaskAtmosTrkEnExit = 0;
01943 fMaskAtmosShwEn = 0;
01944
01945 if( fDebug ){ std::cout << " -*- major systematics (simplex) -*- " << std::endl; }
01946 fUseSimplex = 1; RunFitWithMinuit();
01947
01948 if( fDebug ){ std::cout << " -*- major systematics (migrad) -*- " << std::endl; }
01949 fUseSimplex = 0; RunFitWithMinuit();
01950
01951
01952
01953 fMaskDmsq = 0;
01954 fMaskSinsq = 0;
01955 fMaskBeamNorm = 0;
01956 fMaskBeamNCBkg = 0;
01957 fMaskBeamShwEn = 0;
01958 fMaskBeamTrkEn = 0;
01959 fMaskBeamTrkEnExit = 0;
01960 fMaskAtmosNormCV = 0;
01961 fMaskAtmosNormRock = 0;
01962 fMaskAtmosChgCV = 0;
01963 fMaskAtmosChgRock = 0;
01964 fMaskAtmosNue = 1;
01965 fMaskAtmosNCBkg = 1;
01966 fMaskAtmosSpecNuCV = 0;
01967 fMaskAtmosSpecNuBarCV = 0;
01968 fMaskAtmosSpecNuRock = 0;
01969 fMaskAtmosSpecNuBarRock = 0;
01970 fMaskAtmosZenith = 1;
01971 fMaskAtmosTrkEn = 1;
01972 fMaskAtmosTrkEnExit = 1;
01973 fMaskAtmosShwEn = 1;
01974
01975 if( fDebug ){ std::cout << " -*- minor systematics (migrad) -*- " << std::endl; }
01976 fUseSimplex = 0; RunFitWithMinuit();
01977
01978
01979
01980 fMaskDmsq = 1;
01981 fMaskSinsq = 1;
01982 fMaskBeamNorm = 1;
01983 fMaskBeamNCBkg = 1;
01984 fMaskBeamShwEn = 1;
01985 fMaskBeamTrkEn = 1;
01986 fMaskBeamTrkEnExit = 1;
01987 fMaskAtmosNormCV = 1;
01988 fMaskAtmosNormRock = 1;
01989 fMaskAtmosChgCV = 1;
01990 fMaskAtmosChgRock = 1;
01991 fMaskAtmosNue = 1;
01992 fMaskAtmosNCBkg = 1;
01993 fMaskAtmosSpecNuCV = 1;
01994 fMaskAtmosSpecNuBarCV = 1;
01995 fMaskAtmosSpecNuRock = 1;
01996 fMaskAtmosSpecNuBarRock = 1;
01997 fMaskAtmosZenith = 1;
01998 fMaskAtmosTrkEn = 1;
01999 fMaskAtmosTrkEnExit = 1;
02000 fMaskAtmosShwEn = 1;
02001
02002 if( fDebug ){ std::cout << " -*- all parameters (migrad) -*- " << std::endl; }
02003 fUseSimplex = 0; RunFitWithMinuit();
02004
02005 return;
02006 }
02007
02008
02009 void TwoFlavourFitter::RunFitWithMinuit()
02010 {
02011
02012 fFittingDmsq = 0;
02013 fFittingSinsq = 0;
02014 fFittingBeamNorm = 0;
02015 fFittingBeamNCBkg = 0;
02016 fFittingBeamShwEn = 0;
02017 fFittingBeamTrkEn = 0;
02018 fFittingBeamTrkEnExit = 0;
02019 fFittingAtmosNormCV = 0;
02020 fFittingAtmosNormRock = 0;
02021 fFittingAtmosChgCV = 0;
02022 fFittingAtmosChgRock = 0;
02023 fFittingAtmosNue = 0;
02024 fFittingAtmosNCBkg = 0;
02025 fFittingAtmosSpecNuCV = 0;
02026 fFittingAtmosSpecNuBarCV = 0;
02027 fFittingAtmosSpecNuRock = 0;
02028 fFittingAtmosSpecNuBarRock = 0;
02029 fFittingAtmosZenith = 0;
02030 fFittingAtmosTrkEn = 0;
02031 fFittingAtmosTrkEnExit = 0;
02032 fFittingAtmosShwEn = 0;
02033
02034
02035 fFittingDmsq = ( fMaskDmsq && fFitDmsq );
02036 fFittingSinsq = ( fMaskSinsq && fFitSinsq );
02037
02038 if( UsingBeamData() ){
02039 fFittingBeamNorm = ( fMaskBeamNorm && fFitBeamNorm );
02040 fFittingBeamNCBkg = ( fMaskBeamNCBkg && fFitBeamNCBkg
02041 && TemplateGrid::Instance()->TouchSystematic( OscFit::kBeamData, OscFit::kNCbkg ) );
02042 fFittingBeamShwEn = ( fMaskBeamShwEn && fFitBeamShwEn
02043 && TemplateGrid::Instance()->TouchSystematic( OscFit::kBeamData, OscFit::kShwEn ) );
02044 fFittingBeamTrkEn = ( fMaskBeamTrkEn && fFitBeamTrkEn
02045 && TemplateGrid::Instance()->TouchSystematic( OscFit::kBeamData, OscFit::kTrkEn ) );
02046 fFittingBeamTrkEnExit = ( fMaskBeamTrkEnExit && fFitBeamTrkEnExit
02047 && TemplateGrid::Instance()->TouchSystematic( OscFit::kBeamData, OscFit::kTrkEnExit ) );
02048 }
02049
02050 if( UsingAtmosData() ){
02051 fFittingAtmosNormCV = ( fMaskAtmosNormCV && fFitAtmosNormCV );
02052 fFittingAtmosChgCV = ( fMaskAtmosChgCV && fFitAtmosChgCV );
02053 fFittingAtmosNormRock = ( fMaskAtmosNormRock && fFitAtmosNormRock
02054 && Configuration::Instance()->UsingContainedVertexOnly()==false
02055 && TemplateGrid::Instance()->TouchAtmosRockMC()==true );
02056 fFittingAtmosChgRock = ( fMaskAtmosChgRock && fFitAtmosChgRock
02057 && Configuration::Instance()->UsingContainedVertexOnly()==false
02058 && TemplateGrid::Instance()->TouchAtmosRockMC()==true );
02059 fFittingAtmosNue = ( fMaskAtmosNue && fFitAtmosNue
02060 && Configuration::Instance()->UsingAtmosMuonNeutrinosOnly()==false );
02061 fFittingAtmosNCBkg = ( fMaskAtmosNCBkg && fFitAtmosNCBkg );
02062 fFittingAtmosSpecNuCV = ( fMaskAtmosSpecNuCV && fFitAtmosSpecNuCV
02063 && fMaskAtmosSpecNuBarCV && fFitAtmosSpecNuBarCV
02064 && TemplateGrid::Instance()->TouchSystematic( OscFit::kAtmosData, OscFit::kCV, OscFit::kSpec ) );
02065 fFittingAtmosSpecNuBarCV = ( fMaskAtmosSpecNuCV && fFitAtmosSpecNuCV
02066 && fMaskAtmosSpecNuBarCV && fFitAtmosSpecNuBarCV
02067 && TemplateGrid::Instance()->TouchSystematic( OscFit::kAtmosData, OscFit::kCV, OscFit::kSpec ) );
02068 fFittingAtmosSpecNuRock = ( fMaskAtmosSpecNuRock && fFitAtmosSpecNuRock
02069 && fMaskAtmosSpecNuBarRock && fFitAtmosSpecNuBarRock
02070 && Configuration::Instance()->UsingContainedVertexOnly()==false
02071 && Configuration::Instance()->DoingRockSpectrumSystematic()==true
02072 && TemplateGrid::Instance()->TouchAtmosRockMC()==true
02073 && TemplateGrid::Instance()->TouchSystematic( OscFit::kAtmosData, OscFit::kRock, OscFit::kSpec ) );
02074 fFittingAtmosSpecNuBarRock = ( fMaskAtmosSpecNuRock && fFitAtmosSpecNuRock
02075 && fMaskAtmosSpecNuBarRock && fFitAtmosSpecNuBarRock
02076 && Configuration::Instance()->UsingContainedVertexOnly()==false
02077 && Configuration::Instance()->DoingRockSpectrumSystematic()==true
02078 && TemplateGrid::Instance()->TouchAtmosRockMC()==true
02079 && TemplateGrid::Instance()->TouchSystematic( OscFit::kAtmosData, OscFit::kRock, OscFit::kSpec ) );
02080 fFittingAtmosZenith = ( fMaskAtmosZenith && fFitAtmosZenith
02081 && TemplateGrid::Instance()->TouchSystematic( OscFit::kAtmosData, OscFit::kZenith ) );
02082 fFittingAtmosTrkEn = ( fMaskAtmosTrkEn && fFitAtmosTrkEn
02083 && TemplateGrid::Instance()->TouchSystematic( OscFit::kAtmosData, OscFit::kTrkEn ) );
02084 fFittingAtmosTrkEnExit = ( fMaskAtmosTrkEnExit && fFitAtmosTrkEnExit
02085 && TemplateGrid::Instance()->TouchSystematic( OscFit::kAtmosData, OscFit::kTrkEnExit ) );
02086 fFittingAtmosShwEn = ( fMaskAtmosShwEn && fFitAtmosShwEn
02087 && TemplateGrid::Instance()->TouchSystematic( OscFit::kAtmosData, OscFit::kShwEn ) );
02088 }
02089
02090
02091 if( fCorrelatedSystematics ){
02092
02093 if( Configuration::Instance()->GetExperiment()==OscFit::kMINOS ){
02094 if( UsingBeamData() ){
02095 fFittingBeamTrkEnExit = 0;
02096 if( UsingAtmosData() ){
02097 fFittingAtmosTrkEn = 0;
02098 }
02099 }
02100 }
02101
02102 else if( Configuration::Instance()->GetExperiment()==OscFit::kLBNE ){
02103 if( UsingBeamData() ){
02104 if( UsingAtmosData() ){
02105 fFittingAtmosTrkEn = 0;
02106 fFittingAtmosTrkEnExit = 0;
02107 fFittingAtmosShwEn = 0;
02108 }
02109 }
02110 }
02111 }
02112
02113
02114 GetLnL( fOscDmsq, fOscSinsq,
02115 fBeamNorm, fBeamNCBkg,
02116 fBeamShwEn, fBeamTrkEn,
02117 fBeamTrkEnExit,
02118 fAtmosNormCV, fAtmosNormRock,
02119 fAtmosChgCV, fAtmosChgRock,
02120 fAtmosNue, fAtmosNCBkg,
02121 fAtmosSpecNuCV, fAtmosSpecNuBarCV,
02122 fAtmosSpecNuRock, fAtmosSpecNuBarRock,
02123 fAtmosZenith,
02124 fAtmosTrkEn, fAtmosTrkEnExit,
02125 fAtmosShwEn,
02126 fOutputLnL );
02127
02128
02129
02130 Bool_t runFit = 0;
02131
02132 if( fFittingDmsq ) runFit = 1;
02133 if( fFittingSinsq ) runFit = 1;
02134 if( fFittingBeamNorm ) runFit = 1;
02135 if( fFittingBeamNCBkg ) runFit = 1;
02136 if( fFittingBeamShwEn ) runFit = 1;
02137 if( fFittingBeamTrkEn ) runFit = 1;
02138 if( fFittingBeamTrkEnExit ) runFit = 1;
02139 if( fFittingAtmosNormCV ) runFit = 1;
02140 if( fFittingAtmosNormRock ) runFit = 1;
02141 if( fFittingAtmosChgCV ) runFit = 1;
02142 if( fFittingAtmosChgRock ) runFit = 1;
02143 if( fFittingAtmosNue ) runFit = 1;
02144 if( fFittingAtmosNCBkg ) runFit = 1;
02145 if( fFittingAtmosSpecNuCV ) runFit = 1;
02146 if( fFittingAtmosSpecNuBarCV ) runFit = 1;
02147 if( fFittingAtmosSpecNuRock ) runFit = 1;
02148 if( fFittingAtmosSpecNuBarRock ) runFit = 1;
02149 if( fFittingAtmosZenith ) runFit = 1;
02150 if( fFittingAtmosTrkEn ) runFit = 1;
02151 if( fFittingAtmosTrkEnExit ) runFit = 1;
02152 if( fFittingAtmosShwEn ) runFit = 1;
02153
02154 if( runFit==0 ) return;
02155
02156
02157
02158 Double_t osc_dmsq = fOscDmsq;
02159 Double_t osc_sinsq = fOscSinsq;
02160 Double_t beam_norm = fBeamNorm;
02161 Double_t beam_ncbkg = fBeamNCBkg;
02162 Double_t beam_shwen = fBeamShwEn;
02163 Double_t beam_trken = fBeamTrkEn;
02164 Double_t beam_trken_exit = fBeamTrkEnExit;
02165 Double_t atmos_norm_cv = fAtmosNormCV;
02166 Double_t atmos_norm_rock = fAtmosNormRock;
02167 Double_t atmos_chg_cv = fAtmosChgCV;
02168 Double_t atmos_chg_rock = fAtmosChgRock;
02169 Double_t atmos_nue = fAtmosNue;
02170 Double_t atmos_ncbkg = fAtmosNCBkg;
02171 Double_t atmos_spec_cv_nu = fAtmosSpecNuCV;
02172 Double_t atmos_spec_cv_nubar = fAtmosSpecNuBarCV;
02173 Double_t atmos_spec_rock_nu = fAtmosSpecNuRock;
02174 Double_t atmos_spec_rock_nubar = fAtmosSpecNuBarRock;
02175 Double_t atmos_zenith = fAtmosZenith;
02176 Double_t atmos_trken = fAtmosTrkEn;
02177 Double_t atmos_trken_exit = fAtmosTrkEnExit;
02178 Double_t atmos_shwen = fAtmosShwEn;
02179
02180 Double_t osc_dmsq_err = 0.0;
02181 Double_t osc_sinsq_err = 0.0;
02182 Double_t beam_norm_err = 0.0;
02183 Double_t beam_ncbkg_err = 0.0;
02184 Double_t beam_shwen_err = 0.0;
02185 Double_t beam_trken_err = 0.0;
02186 Double_t beam_trken_exit_err = 0.0;
02187 Double_t atmos_norm_cv_err = 0.0;
02188 Double_t atmos_norm_rock_err = 0.0;
02189 Double_t atmos_chg_cv_err = 0.0;
02190 Double_t atmos_chg_rock_err = 0.0;
02191 Double_t atmos_nue_err = 0.0;
02192 Double_t atmos_ncbkg_err = 0.0;
02193 Double_t atmos_spec_cv_nu_err = 0.0;
02194 Double_t atmos_spec_cv_nubar_err = 0.0;
02195 Double_t atmos_spec_rock_nu_err = 0.0;
02196 Double_t atmos_spec_rock_nubar_err = 0.0;
02197 Double_t atmos_zenith_err = 0.0;
02198 Double_t atmos_trken_err = 0.0;
02199 Double_t atmos_trken_exit_err = 0.0;
02200 Double_t atmos_shwen_err = 0.0;
02201
02202
02203 Int_t err = 0;
02204 Int_t flag = 0;
02205 Int_t npari = 0;
02206 Int_t nparx = 0;
02207 Int_t istat = 0;
02208 Double_t lnl = 0.0;
02209 Double_t fedm = 0.0;
02210 Double_t errdef = 0.0;
02211
02212 Double_t* arglist = new Double_t[10];
02213
02214
02215
02216 this->reset_counter();
02217
02218
02219 fMinuit->mncler();
02220
02221
02222 fMinuit->SetFCN( TwoFlavourLikelihoodFunction );
02223
02224
02225 arglist[0] = 0.5;
02226
02227 fMinuit->mnexcm("SET STR",arglist,1,err);
02228
02229
02230 arglist[0] = 1;
02231
02232 fMinuit->mnexcm("SET STR",arglist,1,err);
02233
02234
02235
02236 if( fFittingDmsq ) fMinuit->mnparm(0, "dmsq", fOscDmsq, 0.0002, 0.0, +0.1, err);
02237 if( fFittingSinsq ){
02238 if( !fPhysicalBoundary ) fMinuit->mnparm(1, "sinsq", fOscSinsq, 0.1, 0.0, +1.5, err);
02239 else fMinuit->mnparm(1, "sinsq", fOscSinsq, 0.1, 0.0, +1.0, err);
02240 }
02241
02242 if( UsingBeamData() ){
02243 if( fFittingBeamNorm ) fMinuit->mnparm(2, "beam_norm", fBeamNorm, 1.0, -3.0, +3.0, err);
02244 if( fFittingBeamNCBkg ) fMinuit->mnparm(3, "beam_ncbkg", fBeamNCBkg, 1.0, -2.5, +2.5, err);
02245 if( fFittingBeamShwEn ) fMinuit->mnparm(4, "beam_shwen", fBeamShwEn, 1.0, -2.5, +2.5, err);
02246 if( fFittingBeamTrkEn ) fMinuit->mnparm(5, "beam_trken", fBeamTrkEn, 1.0, -2.5, +2.5, err);
02247 if( fFittingBeamTrkEnExit ) fMinuit->mnparm(6, "beam_trken_exit", fBeamTrkEnExit, 1.0, -2.5, +2.5, err);
02248 }
02249
02250 if( UsingAtmosData() ){
02251 if( fFittingAtmosNormCV ) fMinuit->mnparm(7, "atmos_norm_cv", fAtmosNormCV, 1.0, -3.0, +3.0, err);
02252 if( fFittingAtmosNormRock ) fMinuit->mnparm(8, "atmos_norm_rock", fAtmosNormRock, 1.0, -3.0, +3.0, err );
02253 if( fFittingAtmosChgCV ) fMinuit->mnparm(9, "atmos_chg_cv", fAtmosChgCV, 1.0, -2.5, +2.5, err );
02254 if( fFittingAtmosChgRock ) fMinuit->mnparm(10, "atmos_chg_rock", fAtmosChgRock, 1.0, -2.5, +2.5, err);
02255 if( fFittingAtmosNue ) fMinuit->mnparm(11, "atmos_nue", fAtmosNue, 1.0, -2.5, +2.5, err);
02256 if( fFittingAtmosNCBkg ) fMinuit->mnparm(12, "atmos_ncbkg", fAtmosNCBkg, 1.0, -2.5, +2.5, err);
02257 if( fFittingAtmosSpecNuCV ) fMinuit->mnparm(13, "atmos_spec_cv_nu", fAtmosSpecNuCV, 1.0, -2.5, +2.5, err);
02258 if( fFittingAtmosSpecNuBarCV ) fMinuit->mnparm(14, "atmos_spec_cv_nubar", fAtmosSpecNuBarCV, 1.0, -2.5, +2.5, err );
02259 if( fFittingAtmosSpecNuRock ) fMinuit->mnparm(15, "atmos_spec_rock_nu", fAtmosSpecNuRock, 1.0, -2.5, +2.5, err);
02260 if( fFittingAtmosSpecNuBarRock ) fMinuit->mnparm(16, "atmos_spec_rock_nubar", fAtmosSpecNuBarRock, 1.0, -2.5, +2.5, err );
02261 if( fFittingAtmosZenith ) fMinuit->mnparm(17, "atmos_zenith", fAtmosZenith, 1.0, -2.5, +2.5, err);
02262 if( fFittingAtmosTrkEn ) fMinuit->mnparm(18, "atmos_trken", fAtmosTrkEn, 1.0, -2.5, +2.5, err);
02263 if( fFittingAtmosTrkEnExit ) fMinuit->mnparm(19, "atmos_trken_exit", fAtmosTrkEnExit, 1.0, -2.5, +2.5, err);
02264 if( fFittingAtmosShwEn ) fMinuit->mnparm(20, "atmos_shwen", fAtmosShwEn, 1.0, -2.5, +2.5, err);
02265 }
02266
02267
02268 if( fUseSimplex ){
02269 if( fDebug ) std::cout << " Fitting With Minuit (Simplex) " << std::endl;
02270 fMinuit->mnsimp(); flag = 0;
02271 }
02272 else{
02273 if( fDebug ) std::cout << " Fitting With Minuit (Migrad) " << std::endl;
02274 flag = fMinuit->Migrad();
02275 }
02276
02277 if( fMinuit->GetNumPars()>0 ){
02278 fMinuitStatus = flag;
02279 }
02280
02281
02282 if( fFittingDmsq ) fMinuit->GetParameter(0, osc_dmsq, osc_dmsq_err);
02283 if( fFittingSinsq ) fMinuit->GetParameter(1, osc_sinsq, osc_sinsq_err);
02284 if( fFittingBeamNorm ) fMinuit->GetParameter(2, beam_norm, beam_norm_err);
02285 if( fFittingBeamNCBkg ) fMinuit->GetParameter(3, beam_ncbkg, beam_ncbkg_err);
02286 if( fFittingBeamShwEn ) fMinuit->GetParameter(4, beam_shwen, beam_shwen_err);
02287 if( fFittingBeamTrkEn ) fMinuit->GetParameter(5, beam_trken, beam_trken_err);
02288 if( fFittingBeamTrkEnExit ) fMinuit->GetParameter(6, beam_trken_exit, beam_trken_exit_err);
02289 if( fFittingAtmosNormCV ) fMinuit->GetParameter(7, atmos_norm_cv , atmos_norm_cv_err);
02290 if( fFittingAtmosNormRock ) fMinuit->GetParameter(8, atmos_norm_rock, atmos_norm_rock_err);
02291 if( fFittingAtmosChgCV ) fMinuit->GetParameter(9, atmos_chg_cv, atmos_chg_cv_err);
02292 if( fFittingAtmosChgRock ) fMinuit->GetParameter(10, atmos_chg_rock, atmos_chg_rock_err);
02293 if( fFittingAtmosNue ) fMinuit->GetParameter(11, atmos_nue, atmos_nue_err);
02294 if( fFittingAtmosNCBkg ) fMinuit->GetParameter(12, atmos_ncbkg, atmos_ncbkg_err);
02295 if( fFittingAtmosSpecNuCV ) fMinuit->GetParameter(13, atmos_spec_cv_nu, atmos_spec_cv_nu_err);
02296 if( fFittingAtmosSpecNuBarCV ) fMinuit->GetParameter(14, atmos_spec_cv_nubar, atmos_spec_cv_nubar_err);
02297 if( fFittingAtmosSpecNuRock ) fMinuit->GetParameter(15, atmos_spec_rock_nu, atmos_spec_rock_nu_err);
02298 if( fFittingAtmosSpecNuBarRock ) fMinuit->GetParameter(16, atmos_spec_rock_nubar, atmos_spec_rock_nubar_err);
02299 if( fFittingAtmosZenith ) fMinuit->GetParameter(17, atmos_zenith, atmos_zenith_err);
02300 if( fFittingAtmosTrkEn ) fMinuit->GetParameter(18, atmos_trken, atmos_trken_err);
02301 if( fFittingAtmosTrkEnExit ) fMinuit->GetParameter(19, atmos_trken_exit, atmos_trken_exit_err);
02302 if( fFittingAtmosShwEn ) fMinuit->GetParameter(20, atmos_shwen, atmos_shwen_err);
02303
02304
02305
02306 fMinuit->mnstat(lnl,fedm,errdef,npari,nparx,istat);
02307
02308
02309 fOscDmsq = osc_dmsq;
02310 fOscSinsq = osc_sinsq;
02311 fBeamNorm = beam_norm;
02312 fBeamNCBkg = beam_ncbkg;
02313 fBeamShwEn = beam_shwen;
02314 fBeamTrkEn = beam_trken;
02315 fBeamTrkEnExit = beam_trken_exit;
02316 fAtmosNormCV = atmos_norm_cv;
02317 fAtmosNormRock = atmos_norm_rock;
02318 fAtmosChgCV = atmos_chg_cv;
02319 fAtmosChgRock = atmos_chg_rock;
02320 fAtmosNue = atmos_nue;
02321 fAtmosNCBkg = atmos_ncbkg;
02322 fAtmosSpecNuCV = atmos_spec_cv_nu;
02323 fAtmosSpecNuBarCV = atmos_spec_cv_nubar;
02324 fAtmosSpecNuRock = atmos_spec_rock_nu;
02325 fAtmosSpecNuBarRock = atmos_spec_rock_nubar;
02326 fAtmosZenith = atmos_zenith;
02327 fAtmosTrkEn = atmos_trken;
02328 fAtmosTrkEnExit = atmos_trken_exit;
02329 fAtmosShwEn = atmos_shwen;
02330
02331 fOutputLnL = lnl;
02332
02333 delete [] arglist;
02334
02335
02336 if( fDebug ){
02337 std::cout << " Iterations: " << std::endl;
02338 std::cout << " Ctr=" << this->get_counter() << std::endl;
02339 std::cout << " Oscillations: " << std::endl;
02340 std::cout << " osc_dmsq: " << osc_dmsq << "\t+/-\t" << osc_dmsq_err << std::endl;
02341 std::cout << " osc_sinsq: " << osc_sinsq << "\t+/-\t" << osc_sinsq_err << std::endl;
02342 if( UsingBeamData() ){
02343 std::cout << " Beam Systematics: " << std::endl;
02344 std::cout << " beam_norm: " << beam_norm << "\t+/-\t" << beam_norm_err << std::endl;
02345 std::cout << " beam_ncbkg: " << beam_ncbkg << "\t+/-\t" << beam_ncbkg_err << std::endl;
02346 std::cout << " beam_shwen: " << beam_shwen << "\t+/-\t" << beam_shwen_err << std::endl;
02347 std::cout << " beam_trken: " << beam_trken << "\t+/-\t" << beam_trken_err << std::endl;
02348 std::cout << " beam_trken_exit: " << beam_trken_exit << "\t+/-\t" << beam_trken_exit_err << std::endl;
02349 }
02350 if( UsingAtmosData() ){
02351 std::cout << " Atmos Systematics: " << std::endl;
02352 std::cout << " atmos_norm_cv: " << atmos_norm_cv << "\t+/-\t" << atmos_norm_cv_err << std::endl;
02353 std::cout << " atmos_norm_rock: " << atmos_norm_rock << "\t+/-\t" << atmos_norm_rock_err << std::endl;
02354 std::cout << " atmos_chg_cv: " << atmos_chg_cv << "\t+/-\t" << atmos_chg_cv_err << std::endl;
02355 std::cout << " atmos_chg_rock: " << atmos_chg_rock << "\t+/-\t" << atmos_chg_rock_err << std::endl;
02356 std::cout << " atmos_nue: " << atmos_nue << "\t+/-\t" << atmos_nue_err << std::endl;
02357 std::cout << " atmos_ncbkg: " << atmos_ncbkg << "\t+/-\t" << atmos_ncbkg_err << std::endl;
02358 std::cout << " atmos_spec_cv_nu: " << atmos_spec_cv_nu << "\t+/-\t" << atmos_spec_cv_nu_err << std::endl;
02359 std::cout << " atmos_spec_cv_nubar: " << atmos_spec_cv_nubar << "\t+/-\t" << atmos_spec_cv_nubar_err << std::endl;
02360 std::cout << " atmos_spec_rock_nu: " << atmos_spec_rock_nu << "\t+/-\t" << atmos_spec_rock_nu_err << std::endl;
02361 std::cout << " atmos_spec_rock_nubar: " << atmos_spec_rock_nubar << "\t+/-\t" << atmos_spec_rock_nubar_err << std::endl;
02362 std::cout << " atmos_zenith: " << atmos_zenith << "\t+/-\t" << atmos_zenith_err << std::endl;
02363 std::cout << " atmos_trken: " << atmos_trken << "\t+/-\t" << atmos_trken_err << std::endl;
02364 std::cout << " atmos_trken_exit: " << atmos_trken_exit << "\t+/-\t" << atmos_trken_exit_err << std::endl;
02365 std::cout << " atmos_shwen: " << atmos_shwen << "\t+/-\t" << atmos_shwen_err << std::endl;
02366 }
02367 std::cout << " Likelihood: " << std::endl;
02368 std::cout << " LnL=" << lnl << std::endl;
02369 }
02370
02371 return;
02372 }
02373
02374 void TwoFlavourFitter::RunContours()
02375 {
02376 std::cout << " *** TwoFlavourFitter::RunContours() *** " << std::endl;
02377
02378
02379
02380 if( fMinuitStatus==999 ){
02381 std::cout << " <warning> Cannot generate contours: fit failed to converge [return] " << std::endl;
02382 return;
02383 }
02384
02385
02386
02387
02388
02389
02390
02391
02392 Int_t Npts = fNumPoints;
02393
02394 Double_t epsilon = 0.0005;
02395
02396 std::cout << " Generating Contours using " << Npts << " points " << std::endl;
02397
02398 TGraph* myContourNuLnL68 = 0;
02399 TGraph* myContourNuLnL90 = 0;
02400 TGraph* myContourNuLnL90_1D = 0;
02401
02402
02403
02404 if( fFittingDmsq && fFittingSinsq ){
02405
02406
02407 std::cout << " Generating Contour... [Dmsq, Sinsq, 68%]" << std::endl;
02408 fMinuit->SetErrorDef( fLnL68 );
02409 TGraph* graphNuLnL68 = (TGraph*)(fMinuit->Contour(Npts-1,1,0));
02410 if( graphNuLnL68 ){
02411 Int_t npts = graphNuLnL68->GetN();
02412 Double_t* x = graphNuLnL68->GetX();
02413 Double_t* y = graphNuLnL68->GetY();
02414 for( Int_t n=0; n<npts; n++ ){
02415 if( fPhysicalBoundary ){
02416 if( x[n]<0.0+epsilon ) x[n] = -0.001;
02417 if( x[n]>1.0-epsilon ) x[n] = +1.001;
02418 }
02419 }
02420 myContourNuLnL68 = new TGraph(npts,x,y);
02421 }
02422
02423
02424 std::cout << " Generating Contour... [Dmsq, Sinsq, 90%]" << std::endl;
02425 fMinuit->SetErrorDef( fLnL90 );
02426 TGraph* graphNuLnL90 = (TGraph*)(fMinuit->Contour(Npts-1,1,0));
02427 if( graphNuLnL90 ){
02428 Int_t npts = graphNuLnL90->GetN();
02429 Double_t* x = graphNuLnL90->GetX();
02430 Double_t* y = graphNuLnL90->GetY();
02431 for( Int_t n=0; n<npts; n++ ){
02432 if( fPhysicalBoundary ){
02433 if( x[n]<0.0+epsilon ) x[n] = -0.001;
02434 if( x[n]>1.0-epsilon ) x[n] = +1.001;
02435 }
02436 }
02437 myContourNuLnL90 = new TGraph(npts,x,y);
02438 }
02439
02440
02441 std::cout << " Generating Contour... [Dmsq, Sinsq, 90% 1D]" << std::endl;
02442 fMinuit->SetErrorDef( fLnL90_1D );
02443 TGraph* graphNuLnL90_1D = (TGraph*)(fMinuit->Contour(Npts-1,1,0));
02444 if( graphNuLnL90_1D ){
02445 Int_t npts = graphNuLnL90_1D->GetN();
02446 Double_t* x = graphNuLnL90_1D->GetX();
02447 Double_t* y = graphNuLnL90_1D->GetY();
02448 for( Int_t n=0; n<npts; n++ ){
02449 if( fPhysicalBoundary ){
02450 if( x[n]<0.0+epsilon ) x[n] = -0.001;
02451 if( x[n]>1.0-epsilon ) x[n] = +1.001;
02452 }
02453 }
02454 myContourNuLnL90_1D = new TGraph(npts,x,y);
02455 }
02456
02457 }
02458
02459
02460
02461 TString filename = fOutFile;
02462 TString tag = ".contours";
02463
02464 if( filename.EndsWith(".root") )
02465 filename.Insert(filename.Length()-5,tag);
02466 else filename.Append(tag);
02467
02468 std::cout << " *** TwoFlavourFitter::WriteContours() *** " << std::endl;
02469 std::cout << " Writing Contours to File: " << filename.Data() << std::endl;
02470
02471 TDirectory* tmpd = gDirectory;
02472 TFile* fileNew = new TFile(filename.Data(),"recreate");
02473
02474 if( myContourNuLnL68 ) myContourNuLnL68->Write("gContourNuLnL68");
02475 if( myContourNuLnL90 ) myContourNuLnL90->Write("gContourNuLnL90");
02476 if( myContourNuLnL90_1D ) myContourNuLnL90_1D->Write("gContourNuLnL90_1D");
02477
02478 fileNew->Close();
02479 tmpd->cd();
02480
02481 return;
02482 }
02483
02484 Double_t TwoFlavourFitter::GetLikelihood(Double_t input_dmsq, Double_t input_sinsq)
02485 {
02486 Double_t osc_dmsq = 0.0;
02487 Double_t osc_sinsq = 0.0;
02488 Double_t beam_norm = 0.0;
02489 Double_t beam_ncbkg = 0.0;
02490 Double_t beam_shwen = 0.0;
02491 Double_t beam_trken = 0.0;
02492 Double_t beam_trken_exit = 0.0;
02493 Double_t atmos_norm_cv = 0.0;
02494 Double_t atmos_norm_rock = 0.0;
02495 Double_t atmos_chg_cv = 0.0;
02496 Double_t atmos_chg_rock = 0.0;
02497 Double_t atmos_nue = 0.0;
02498 Double_t atmos_ncbkg = 0.0;
02499 Double_t atmos_spec_cv_nu = 0.0;
02500 Double_t atmos_spec_cv_nubar = 0.0;
02501 Double_t atmos_spec_rock_nu = 0.0;
02502 Double_t atmos_spec_rock_nubar = 0.0;
02503 Double_t atmos_zenith = 0.0;
02504 Double_t atmos_trken = 0.0;
02505 Double_t atmos_trken_exit = 0.0;
02506 Double_t atmos_shwen = 0.0;
02507
02508 Double_t LnL = 0.0;
02509
02510 this->GetLnL( input_dmsq, input_sinsq,
02511 osc_dmsq, osc_sinsq,
02512 beam_norm, beam_ncbkg,
02513 beam_shwen, beam_trken,
02514 beam_trken_exit,
02515 atmos_norm_cv, atmos_norm_rock,
02516 atmos_chg_cv, atmos_chg_rock,
02517 atmos_nue, atmos_ncbkg,
02518 atmos_spec_cv_nu, atmos_spec_cv_nubar,
02519 atmos_spec_rock_nu, atmos_spec_rock_nubar,
02520 atmos_zenith,
02521 atmos_trken, atmos_trken_exit,
02522 atmos_shwen,
02523 LnL );
02524 return LnL;
02525 }
02526
02527 void TwoFlavourFitter::GetLnL( Double_t input_dmsq, Double_t input_sinsq, Double_t& osc_dmsq, Double_t& osc_sinsq, Double_t& beam_norm, Double_t& beam_ncbkg, Double_t& beam_shwen, Double_t& beam_trken, Double_t& beam_trken_exit, Double_t& atmos_norm_cv, Double_t& atmos_norm_rock, Double_t& atmos_chg_cv, Double_t& atmos_chg_rock, Double_t& atmos_nue, Double_t& atmos_ncbkg, Double_t& atmos_spec_cv_nu, Double_t& atmos_spec_cv_nubar, Double_t& atmos_spec_rock_nu, Double_t& atmos_spec_rock_nubar, Double_t& atmos_zenith, Double_t& atmos_trken, Double_t& atmos_trken_exit, Double_t& atmos_shwen, Double_t& LnL )
02528 {
02529 return RunFit( input_dmsq, input_sinsq,
02530 osc_dmsq, osc_sinsq,
02531 beam_norm, beam_ncbkg,
02532 beam_shwen, beam_trken,
02533 beam_trken_exit,
02534 atmos_norm_cv, atmos_norm_rock,
02535 atmos_chg_cv, atmos_chg_rock,
02536 atmos_nue, atmos_ncbkg,
02537 atmos_spec_cv_nu, atmos_spec_cv_nubar,
02538 atmos_spec_rock_nu, atmos_spec_rock_nubar,
02539 atmos_zenith,
02540 atmos_trken, atmos_trken_exit,
02541 atmos_shwen,
02542 LnL );
02543 }
02544
02545 void TwoFlavourFitter::GetLnL(Double_t osc_dmsq, Double_t osc_sinsq, Double_t beam_norm, Double_t beam_ncbkg, Double_t beam_shwen, Double_t beam_trken, Double_t beam_trken_exit, Double_t atmos_norm_cv, Double_t atmos_norm_rock, Double_t atmos_chg_cv, Double_t atmos_chg_rock, Double_t atmos_nue, Double_t atmos_ncbkg, Double_t atmos_spec_cv_nu, Double_t atmos_spec_cv_nubar, Double_t atmos_spec_rock_nu, Double_t atmos_spec_rock_nubar, Double_t atmos_zenith, Double_t atmos_trken, Double_t atmos_trken_exit, Double_t atmos_shwen, Double_t& LnL)
02546 {
02547
02548
02549 LnL = 0.0;
02550
02551
02552
02553
02554 Double_t input_dmsq = osc_dmsq;
02555 Double_t input_sinsq = osc_sinsq;
02556
02557
02558 Bool_t checkOsc = 0;
02559 if( fabs(input_dmsq)>0.0 ) checkOsc = 1;
02560
02561
02562
02563
02564 Double_t input_beam_shwen = beam_shwen;
02565 Double_t input_beam_trken = beam_trken;
02566 Double_t input_beam_trken_exit = beam_trken_exit;
02567 Double_t input_atmos_trken = atmos_trken;
02568 Double_t input_atmos_trken_exit = atmos_trken_exit;
02569 Double_t input_atmos_shwen = atmos_shwen;
02570
02571
02572 if( fCorrelatedSystematics ){
02573
02574 if( Configuration::Instance()->GetExperiment()==OscFit::kMINOS ){
02575 if( UsingBeamData() ){
02576 input_beam_trken_exit = input_beam_trken;
02577 if( UsingAtmosData() ){
02578 input_atmos_trken = input_beam_trken;
02579 input_atmos_shwen = 0.7071*(input_atmos_shwen+input_beam_shwen);
02580 }
02581 }
02582 }
02583
02584 else if( Configuration::Instance()->GetExperiment()==OscFit::kLBNE ){
02585 if( UsingBeamData() ){
02586 if( UsingAtmosData() ){
02587 input_atmos_trken = input_beam_trken;
02588 input_atmos_trken_exit = input_beam_trken_exit;
02589 input_atmos_shwen = input_beam_shwen;
02590 }
02591 }
02592 }
02593 }
02594
02595
02596
02597
02598 if( checkOsc ){
02599
02600
02601 assert( TemplateGrid::Instance()->TouchGrid() );
02602
02603
02604 if( fDmsqBoundary ){
02605 Double_t delta_dmsq = 0.0;
02606 Double_t min_dmsq = TemplateGrid::Instance()->GetDmsqMin();
02607 Double_t max_dmsq = TemplateGrid::Instance()->GetDmsqMax();
02608 if( input_dmsq>max_dmsq ){
02609 delta_dmsq = max_dmsq - input_dmsq;
02610 input_dmsq = max_dmsq;
02611 }
02612 if( input_dmsq<min_dmsq ){
02613 delta_dmsq = input_dmsq - min_dmsq;
02614 input_dmsq = min_dmsq;
02615 }
02616 LnL += 1.0e7*(delta_dmsq*delta_dmsq);
02617 }
02618
02619
02620 if( fPhysicalBoundary ){
02621 Double_t delta_sinsq = 0.0;
02622 if( input_sinsq>1.0 ){
02623 delta_sinsq = input_sinsq - 1.0;
02624 }
02625 if( input_sinsq<0.0 ){
02626 delta_sinsq = 0.0 - input_sinsq;
02627 }
02628 LnL += 1.0e4*(delta_sinsq*delta_sinsq);
02629 }
02630
02631 }
02632
02633
02634
02635
02636
02637
02638
02639
02640
02641
02642
02643 TemplateCalculator::Instance()->GetExpectation( input_dmsq, input_sinsq,
02644 beam_norm, beam_ncbkg,
02645 input_beam_shwen, input_beam_trken,
02646 input_beam_trken_exit,
02647 atmos_norm_cv, atmos_norm_rock,
02648 atmos_chg_cv, atmos_chg_rock,
02649 atmos_nue, atmos_ncbkg,
02650 atmos_spec_cv_nu, atmos_spec_cv_nubar,
02651 atmos_spec_rock_nu, atmos_spec_rock_nubar,
02652 atmos_zenith,
02653 input_atmos_trken, input_atmos_trken_exit,
02654 input_atmos_shwen,
02655 fTemplateMC );
02656
02657 LnL += OscFit::CalcLikelihood( fTemplateMC, fTemplateExpt );
02658
02659
02660
02661
02662 if( fPenaltyTerms ){
02663
02664
02665 LnL += 0.5*( beam_norm*beam_norm
02666 + beam_ncbkg*beam_ncbkg
02667 + beam_shwen*beam_shwen
02668 + beam_trken*beam_trken
02669 + beam_trken_exit*beam_trken_exit
02670 + atmos_norm_cv*atmos_norm_cv
02671 + atmos_norm_rock*atmos_norm_rock
02672 + atmos_chg_cv*atmos_chg_cv
02673 + atmos_chg_rock*atmos_chg_rock
02674 + atmos_nue*atmos_nue
02675 + atmos_ncbkg*atmos_ncbkg
02676 + atmos_spec_cv_nu*atmos_spec_cv_nu
02677 + atmos_spec_cv_nubar*atmos_spec_cv_nubar
02678 + atmos_spec_rock_nu*atmos_spec_rock_nu
02679 + atmos_spec_rock_nubar*atmos_spec_rock_nubar
02680 + atmos_zenith*atmos_zenith
02681 + atmos_trken*atmos_trken
02682 + atmos_trken_exit*atmos_trken_exit
02683 + atmos_shwen*atmos_shwen );
02684
02685
02686 if( fCorrelatedSystematics ){
02687
02688 if( Configuration::Instance()->GetExperiment()==OscFit::kMINOS ){
02689 if( UsingBeamData() ){
02690 LnL -= 0.5*beam_trken_exit*beam_trken_exit;
02691 if( UsingAtmosData() ){
02692 LnL -= 0.5*( atmos_shwen*atmos_shwen
02693 + atmos_trken*atmos_trken );
02694 LnL += 0.5*input_atmos_shwen*input_atmos_shwen;
02695 }
02696 }
02697 }
02698
02699 else if( Configuration::Instance()->GetExperiment()==OscFit::kLBNE ){
02700 if( UsingBeamData() ){
02701 if( UsingAtmosData() ){
02702 LnL -= 0.5*( atmos_shwen*atmos_shwen
02703 + atmos_trken_exit*atmos_trken_exit
02704 + atmos_trken*atmos_trken );
02705 }
02706 }
02707 }
02708 }
02709 }
02710
02711
02712
02713 assert( LnL<fBigNumber );
02714
02715 return;
02716 }
02717
