NuTransitionFitterMinuit.cxx

Go to the documentation of this file.

// $Id: NuTransitionFitterMinuit.cxx,v 1.4 2008/03/18 21:14:22 ahimmel Exp $
//
// Performs Matrix Method analyses using Minuit2.
//
// Alexander Himmel
// ahimmel@caltech.edu
//
// $Log: NuTransitionFitterMinuit.cxx,v $
// Revision 1.4  2008/03/18 21:14:22  ahimmel
//
//
// Some small changes to make the fitter happy.
//
// Revision 1.3  2008/03/07 23:51:07  ahimmel
//
//
// Minor changes to keep Minuit happy.
//
// Revision 1.2  2008/03/03 22:05:00  ahimmel
//
//
// Some bug fixes.
//
// Revision 1.1  2008/03/01 00:29:33  ahimmel
//
//
// A minuit fitter for transitions.  Should speed up the fitting stage.  Also, more careful memory management should increase the maximum number of trial experiments.
//

#include "TFile.h"
#include "TFitterMinuit.h"
#include "TH1D.h"
#include "TH1F.h"
#include "TMath.h"
#include "TROOT.h"

#include "Conventions/Munits.h"
#include "MessageService/MsgService.h"
#include "NtupleUtils/NuTransitionFitterMinuit.h"
#include "NtupleUtils/NuMatrixInput.h"
#include "NtupleUtils/NuMatrixOutput.h"
#include "NtupleUtils/NuFluctuator.h"
#include "NtupleUtils/NuXMLConfig.h"

ClassImp(NuTransitionFitterMinuit)

CVSID("$Id: NuTransitionFitterMinuit.cxx,v 1.4 2008/03/18 21:14:22 ahimmel Exp $");


NuTransitionFitterMinuit::NuTransitionFitterMinuit()
{
    mmInput = 0;
    fdData = 0;
    fdFitData = 0;
    fxml = 0;
    
    outputFileName = "";
    
    fFitter = new TFitterMinuit();
    fFitter->CreateMinimizer();
    fFitter->SetMinuitFCN(this);
}


NuTransitionFitterMinuit::NuTransitionFitterMinuit(const TString mmFile, 
                                                   const TString fdFile, 
                                                   const TString outFile)
{   
    // mmInput, fxml
    TFile *infile = new TFile(mmFile);
    mmInput = new NuMatrixInput(infile);
    TString objName = "NuMatrixOutput";
    fxml = (NuXMLConfig*)gROOT->FindObject("NuXMLConfig");
    infile->Close();   
    
    // fdData
    TFile *fdDataFile = new TFile(fdFile,"READ");
    fdData = (TH1D*) fdDataFile->Get("RecoEnergyPQ_FD");
    fdData->SetDirectory(0);
    fdData->Sumw2();
    fdDataFile->Close();
    
    outputFileName = outFile;
    
    fdFitData = 0;
    
    fFitter = new TFitterMinuit();
    fFitter->CreateMinimizer();
    fFitter->SetMinuitFCN(this);
}


NuTransitionFitterMinuit::~NuTransitionFitterMinuit()
{
    if (mmInput) delete mmInput;
    if (fdData) delete fdData;
    if (fdFitData) delete fdFitData;
    if (fxml) delete fxml;    
    
    if (fFitter) delete fFitter;
}


void NuTransitionFitterMinuit::SetupPars()
{
    // Fitted Physics Parameters
    fFitter->SetParameter(0, "prob", 
                          0.1, 0.01, 
                          1.0, 0.0);
    
    // Fixed Physics Parameters
    fFitter->SetParameter(1, "dm2", 
                          fxml->DM2Nu(), 0.1,
                          1.0, 0.0);
    if (!fFitter->IsFixed(1)) fFitter->FixParameter(1);  
    
    fFitter->SetParameter(2, "sn2", 
                          fxml->SN2Nu(), 0.01,
                          1.0,0.0);
    if (!fFitter->IsFixed(2)) fFitter->FixParameter(2);
    
    if (fxml->DM2Nu() !=  fxml->DM2Bar() || fxml->SN2Nu() !=  fxml->SN2Bar()) {
        MAXMSG("NuTransitionFitterMinuit",Msg::kWarning,1) 
        << "Nu and NuBar oscillation parameters don't match.  This is not accounted for." << endl;
    }
    
    /*
    // Systematic Parameters (fixed for now)
    fFitter->SetParameter(3, "Normalisation", 
                          1.0, 0.01,
                          1.0, 0.0);
    if (!fFitter->IsFixed(3)) fFitter->FixParameter(3);
    */
}


const Double_t NuTransitionFitterMinuit::CalculateChi2(const TH1D* fdPred, const TH1D* fdData) const
{
    Double_t chi2 = 0;
    Double_t mnu, dnu;
    for (Int_t i=1; i<=fdPred->GetNbinsX(); ++i){
        mnu = fdPred->GetBinContent(i);
        dnu = fdData->GetBinContent(i);
        if (mnu<0.0001) mnu = 0.0001;
        chi2 += (dnu-mnu)*(dnu-mnu)/mnu;
    }
    return chi2;
}


const Double_t NuTransitionFitterMinuit::CalculateLikelihood(const TH1D* fdPred, const TH1D* fdData) const
{
    Double_t like = 0;   
    Double_t mnu, dnu;
    for (Int_t i=6; i<=fdPred->GetNbinsX(); ++i){
        mnu = fdPred->GetBinContent(i);
        dnu = fdData->GetBinContent(i);
        if (mnu<0.0001) mnu = 0.0001;
        like += 2*(mnu - dnu);
        if (dnu) like += 2*dnu*log(dnu/mnu);
    }
    return like;
}


/*
const Double_t NuTransitionFitterMinuit::CalculateUnbinnedLikelihood(const TH1D* fdPred, vector<double> fdData) const
{
    // Just a placeholder right now
}
*/


void NuTransitionFitterMinuit::TransitionFit(const Int_t experiments)
{
  // Initialize
    TString name = "";
    TString objName = "NuMatrixOutput";
    
    MSG("NuTransitionFitterMinuit",Msg::kInfo) << "Recreating file " << outputFileName << endl;
    TFile *outputFile = new TFile(outputFileName,"RECREATE");
    
    NuMatrixOutput *mmOutput = 0;
    NuFluctuator *fl = new NuFluctuator();
    
    TH1F *bestFits = new TH1F("bestFits","Best Fit Points",175,-0.5,1.25);
    
    
    // Perform an unfluctuated fit
    mmOutput = DoTransitionFit(*fdData);
    
    if (fxml) mmOutput->XMLConfig(fxml);    
    name = objName + "_";
    name += "nofluct";
    outputFile->WriteObject(mmOutput,name);
    
    delete mmOutput;
    
    // Perform fluctuated fits
    for (int i = 0; i < experiments; i++) {
        mmOutput = DoTransitionFit(fl->Fluctuate(fdData));

        bestFits->Fill(mmOutput->bestTransitionProb);
        
        if (fxml) mmOutput->XMLConfig(fxml);
        name = objName + "_";
        name += i;
        outputFile->WriteObject(mmOutput,name);
        
        delete mmOutput;
    }
    
    // Write Best Fits and close the file
    outputFile->WriteObject(bestFits,bestFits->GetName());
    outputFile->Close();
    MSG("NuTransitionFitterMinuit",Msg::kInfo) 
    << "Don't mind the segfault.  Everything is done and it appears to be harmless." << endl;
}


NuMatrixOutput* NuTransitionFitterMinuit::DoTransitionFit(const TH1D &hDataToFit)
{
    // Setup
    this->SetupPars();
    fdFitData = new TH1D(hDataToFit);
    
    for (int i = 1; i < 3; i++) {
        if (!fFitter->IsFixed(i)) {
            cout << "Refixing parameter " << i << endl;
            fFitter->FixParameter(i);
        }
    }

    // Perform Fit
    fFitter->Minimize();
    
    // Get Results
    //fFitter->PrintResults(0,0);
    
    Double_t prob = fFitter->GetParameter(0);
    Double_t dm2 = fFitter->GetParameter(1);
    Double_t sn2 = fFitter->GetParameter(2);

    // Create Output Object
    NuMatrixOutput *mmOutput = new NuMatrixOutput();
    mmOutput->NuMuBarFDData(fdFitData);    
    mmOutput->NuMuBarBestFitFDPrediction(*mmInput->PredictFDSpectrumTransition(dm2, sn2, prob));
    mmOutput->NuMuBarNoTransPrediction(*mmInput->PredictFDSpectrumTransition(dm2, sn2, 0.));
    mmOutput->NuMuBarNoOscPrediction(*mmInput->PredictFDSpectrumTransition(0., 0., 0.));
    mmOutput->BestFitPoint(sn2, dm2, prob);
    
    // Clean up and return
    delete fdFitData; fdFitData = 0;
    return mmOutput;
}


double NuTransitionFitterMinuit::operator () (const vector<double>& pars) const
{
    // Track Fitter Calls
    static Int_t fitcounter = 0;
    if (!(fitcounter%50)){
        MSG("NuTransitionFitterMinuit",Msg::kDebug) << "Minuit TF call " << fitcounter << endl;
    }
    ++fitcounter;

    // Get Current Parameters
    Double_t prob = pars[0];
    Double_t dm2 = pars[1];
    Double_t sn2 = pars[2];
    //Double_t normalisation = pars[3];
    
    // Predict the FD Spectrum for these parameters
    TH1D *fdFitPred = mmInput->PredictFDSpectrumTransition(dm2, sn2, prob);
    
    // Apply systematics
    // fdFitPred->Scale(normalisation);

    // Calculate our goodness-of-fit test, T
    double T = this->CalculateLikelihood(fdFitPred,fdFitData);
    
    // Apply penalty terms
    // Example of normalisation penalty term for center = 1, error = 0.04
    // T += (normalisation - 1)*(normalisation - 1) / (0.04 * 0.04); 
    
    // Clean Up
    delete fdFitPred; fdFitPred = 0;

    return T;
}

---