loon: NuMMParameters.cxx Source File

00001 #include "MessageService/MsgService.h"
00002 #include "NtupleUtils/NuMMParameters.h"
00003 #include "Minuit2/FunctionMinimum.h"
00004 #include "Minuit2/MnUserParameters.h"
00005 
00006 
00007 ClassImp(NuMMParameters)
00008 
00009 CVSID("$Id: NuMMParameters.cxx,v 1.21 2010/01/15 18:56:54 bckhouse Exp $");
00010 
00011 Double_t NuMMParameters::fncBackOneSigma = 0.5;
00012 
00013 //____________________________________________________________________72
00014 NuMMParameters::NuMMParameters()
00015   : fdm2Fixed(false),
00016     fsn2Fixed(false),
00017     fnormFixed(false),
00018     fshwScaleFixed(false),
00019     fncBackFixed(false),
00020     fdm2BarFixed(false),
00021     fsn2BarFixed(false),
00022     ftransitionProbFixed(true),
00023     fDm2Constrained(false),
00024     fDm2BarConstrained(false),
00025     fPhysical(false),
00026     fCPT(false),
00027     fDm2LowLimit(0.0),
00028     fDm2HighLimit(-1.0),
00029     fDm2BarLowLimit(0.0),
00030     fDm2BarHighLimit(-1.0),
00031     fContourPars(0,1),
00032     fminuitPass(false),
00033     fSn2LowLimit(0.0),
00034     fSn2HighLimit(1.0),
00035     fOneSidedMass(false)
00036 {
00037   fParameters.push_back(0);
00038   fParameters.push_back(0);
00039   fParameters.push_back(1.0);
00040   fParameters.push_back(0.0);
00041   fParameters.push_back(1.0);
00042   fParameters.push_back(0);
00043   fParameters.push_back(0);
00044   fParameters.push_back(0);
00045 }
00046 
00047 //____________________________________________________________________72
00048 NuMMParameters::NuMMParameters(const std::vector<double>& pars)
00049   : fdm2Fixed(false),
00050     fsn2Fixed(false),
00051     fnormFixed(false),
00052     fshwScaleFixed(false),
00053     fncBackFixed(false),
00054     fdm2BarFixed(false),
00055     fsn2BarFixed(false),
00056     ftransitionProbFixed(true),
00057     fDm2Constrained(false),
00058     fDm2BarConstrained(false),
00059     fPhysical(false),
00060     fCPT(false),
00061     fDm2LowLimit(0.0),
00062     fDm2HighLimit(-1.0),
00063     fDm2BarLowLimit(0.0),
00064     fDm2BarHighLimit(-1.0),
00065     fContourPars(0,1),
00066     fminuitPass(false),
00067     fSn2LowLimit(0.0),
00068     fSn2HighLimit(1.0),
00069     fOneSidedMass(false)
00070 {
00071     fParameters=pars;
00072 }
00073 
00074 //____________________________________________________________________72
00075 NuMMParameters::NuMMParameters(const ROOT::Minuit2::FunctionMinimum& minPars)
00076   : fdm2Fixed(false),
00077     fsn2Fixed(false),
00078     fnormFixed(false),
00079     fshwScaleFixed(false),
00080     fncBackFixed(false),
00081     fdm2BarFixed(false),
00082     fsn2BarFixed(false),
00083     ftransitionProbFixed(true),
00084     fDm2Constrained(false),
00085     fDm2BarConstrained(false),
00086     fPhysical(false),
00087     fCPT(false),
00088     fDm2LowLimit(0.0),
00089     fDm2HighLimit(-1.0),
00090     fDm2BarLowLimit(0.0),
00091     fDm2BarHighLimit(-1.0),
00092     fContourPars(0,1)
00093 {
00094   for (Int_t i=0; i<8; ++i){
00095     fParameters.push_back(minPars.UserParameters().Parameter(i).Value());
00096   }
00097 }
00098 
00099 //____________________________________________________________________72
00100 NuMMParameters::~NuMMParameters()
00101 {
00102 }
00103 
00104 //____________________________________________________________________72
00105 const Bool_t NuMMParameters::AreAllParametersFixed() const
00106 {
00107   if (fdm2Fixed &&
00108       fsn2Fixed &&
00109       fnormFixed &&
00110       fshwScaleFixed &&
00111       fncBackFixed &&
00112       fdm2BarFixed &&
00113       fsn2BarFixed &&
00114       ftransitionProbFixed){
00115     return true;
00116   }
00117   else {return false;}
00118 }
00119 
00120 //____________________________________________________________________72
00121 const Double_t NuMMParameters::PenaltyTerm() const
00122 {
00123   const Double_t normalisation = this->Normalisation();
00124   const Double_t shwEn = this->ShwEnScale();
00125   const Double_t ncBack = this->NCBackgroundScale();
00126 
00127   Double_t normOneSigma = 0.04;
00128   Double_t shwEnOneSigma = 1.0;
00129   Double_t ncBackOneSigma = fncBackOneSigma;
00130   Double_t penalty = ((normalisation - 1.0)*(normalisation-1.0))/
00131     (normOneSigma*normOneSigma);
00132   penalty += ((shwEn - 0.0)*(shwEn - 0.0))/
00133     (shwEnOneSigma*shwEnOneSigma);
00134   penalty += ((ncBack - 1.0)*(ncBack - 1.0))/
00135     (ncBackOneSigma*ncBackOneSigma);
00136   return penalty;
00137 }
00138 
00139 //____________________________________________________________________72
00140 ROOT::Minuit2::MnUserParameters NuMMParameters::MinuitParameters() const
00141 {
00142   ROOT::Minuit2::MnUserParameters mnPars;
00143   
00144   // Neutrino Parameters
00145   if (fDm2Constrained) mnPars.Add("Dm2",this->Dm2(),0.005e-3,
00146                                   fDm2LowLimit,fDm2HighLimit);
00147   else                 mnPars.Add("Dm2",this->Dm2(),0.005e-3);
00148   
00149   if (fPhysical) mnPars.Add("Sn2",this->Sn2(),0.01,fSn2LowLimit,fSn2HighLimit);
00150   else           mnPars.Add("Sn2",this->Sn2(),0.01);
00151 
00152   // Systematics
00153   mnPars.Add("Norm",this->Normalisation(),0.01);
00154   mnPars.Add("ShwScale",this->ShwEnScale(),1.0);
00155   mnPars.Add("NCBack",this->NCBackgroundScale(),0.05);
00156 
00157   // Antineutrino Parameters
00158   if (fDm2BarConstrained) {
00159     // We want to constrain dm2bar
00160     if (fOneSidedMass) {
00161       // Make the constraint one sided
00162       mnPars.Add("Dm2Bar", this->Dm2Bar(), 0.005e-3);
00163       mnPars.SetLowerLimit("Dm2Bar", fDm2BarLowLimit);
00164     } else {
00165       // A two sided constraint is fine
00166       mnPars.Add("Dm2Bar",this->Dm2Bar(),0.005e-3,
00167                   fDm2BarLowLimit,fDm2BarHighLimit);
00168     }
00169               
00170   } else {
00171       // Do not constrain dm2bar
00172       mnPars.Add("Dm2Bar",this->Dm2Bar(),0.005e-3);
00173   }
00174   
00175   if (fPhysical) mnPars.Add("Sn2Bar",this->Sn2Bar(),0.01,fSn2LowLimit,fSn2HighLimit);
00176   else           mnPars.Add("Sn2Bar",this->Sn2Bar(),0.01);
00177   
00178   if (fPhysical) mnPars.Add("TransProb",this->TransitionProb(), 0.01, 0.0, 1.0);
00179   else           mnPars.Add("TransProb",this->TransitionProb(), 0.01);
00180 
00181 
00182   if (fdm2Fixed){mnPars.Fix("Dm2");}
00183   if (fsn2Fixed){mnPars.Fix("Sn2");}
00184   if (fnormFixed){mnPars.Fix("Norm");}
00185   if (fshwScaleFixed){mnPars.Fix("ShwScale");}
00186   if (fncBackFixed){mnPars.Fix("NCBack");}
00187   if (fdm2BarFixed){mnPars.Fix("Dm2Bar");}
00188   if (fsn2BarFixed){mnPars.Fix("Sn2Bar");}
00189   if (ftransitionProbFixed){mnPars.Fix("TransProb");}
00190   return mnPars;
00191 }
00192 
00193 void NuMMParameters::PrintStatus() const 
00194 {
00195     MSG("NuMMParameters",Msg::kInfo) << endl << "=====NuMMParameters Status====" << endl
00196     << "Dm2    = " << setw(8) << Dm2()               << " (" << (fdm2Fixed?"Fixed":"Free") << ")" << endl
00197     << "Sn2    = " << setw(8) << Sn2()               << " (" << (fsn2Fixed?"Fixed":"Free") << ")" << endl
00198     << "Dm2Bar = " << setw(8) << Dm2Bar()            << " (" << (fdm2BarFixed?"Fixed":"Free") << ")" << endl
00199     << "Sn2Bar = " << setw(8) << Sn2Bar()            << " (" << (fsn2BarFixed?"Fixed":"Free") << ")" << endl
00200     << "Trans  = " << setw(8) << TransitionProb()    << " (" << (ftransitionProbFixed?"Fixed":"Free") << ")" << endl
00201     << "Norm   = " << setw(8) << Normalisation()     << " (" << (fnormFixed?"Fixed":"Free") << ")" << endl
00202     << "ShwEn  = " << setw(8) << ShwEnScale()        << " (" << (fshwScaleFixed?"Fixed":"Free") << ")" << endl
00203     << "NCBkgd = " << setw(8) << NCBackgroundScale() << " (" << (fncBackFixed?"Fixed":"Free") << ")" << endl
00204     << "Constrain to CPT = " << (fCPT ? "Yes" : "No") << endl
00205     << "Constrain to Physical = " << (fPhysical ? "Yes" : "No") << endl
00206     << "Oone-sided mass = " << (fOneSidedMass ? "Yes" : "No") << endl;
00207     if (fDm2Constrained) MSG("NuMMParameters",Msg::kInfo) << "Dm2 Constrained to [" << fDm2LowLimit << ", " << fDm2HighLimit << "]" << endl;
00208     else                 MSG("NuMMParameters",Msg::kInfo) << "Dm2 Unconstrained" << endl;
00209     if (fDm2BarConstrained) MSG("NuMMParameters",Msg::kInfo) << "Dm2Bar Constrained to [" << fDm2BarLowLimit << ", " << fDm2BarHighLimit << "]" << endl;
00210     else                    MSG("NuMMParameters",Msg::kInfo) << "Dm2Bar Unconstrained" << endl;
00211 
00212 }