NueFitter Class Reference

#include <NueFitter.h>

List of all members.

Public Member Functions

 NueFitter ()
 ~NueFitter ()
void Reset ()
void Initialize ()
void AddFile (std::string file)
void AddContour (double cont)
void SetOutputFile (std::string out)
bool PerformFit (double input, std::string outName="dummy.root")
bool PrepareFit ()
bool FitInput (double input, std::string outName)
bool LoadFiles ()
bool BuildChi2Map (double val)
bool BuildContourMaps ()
void WriteToFile ()
void SetFitMethod (int meth)

Private Member Functions

bool InitializeFittingHistograms ()
double CalculateDeltaLog (double nexp, double nobs)
double CalculateDeltaLog (double nexp, double nobs, double nsig)

Private Attributes

std::vector< std::string > fFileList
std::vector< double > fContourLevels
std::map< double, std::map
< double, std::vector< Point > > > 
fDeltaM32MapNormal
std::map< double, std::map
< double, std::vector< Point > > > 
fDeltaM32MapInvert
TH2D * fFitPointChi2N
TH2D * fFitPointChi2I
std::vector< TH2D * > fContourHistsN
std::vector< TH2D * > fContourHistsI
double fMinTotalEvents
double fErrors [10]
std::string fOutFile
int fFitMethod

Detailed Description

Definition at line 27 of file NueFitter.h.


Constructor & Destructor Documentation

NueFitter::NueFitter (  ) 

Definition at line 11 of file NueFitter.cxx.

References fFitMethod, fMinTotalEvents, and fOutFile.

00012 {
00013    fOutFile = "DefaultOut.root";
00014    fMinTotalEvents = -10;
00015    fFitMethod = 0;
00016 }

NueFitter::~NueFitter (  ) 

Member Function Documentation

void NueFitter::AddContour ( double  cont  )  [inline]

Definition at line 37 of file NueFitter.h.

References fContourLevels.

00037 {fContourLevels.push_back(cont);};

void NueFitter::AddFile ( std::string  file  )  [inline]

Definition at line 36 of file NueFitter.h.

References fFileList.

00036 { fFileList.push_back(file); };

bool NueFitter::BuildChi2Map ( double  val  ) 

Definition at line 376 of file NueFitter.cxx.

References CalculateDeltaLog(), fDeltaM32MapInvert, fDeltaM32MapNormal, fFitMethod, fFitPointChi2I, and fFitPointChi2N.

00377 {
00378   double dmDV = fDeltaM32MapNormal.begin()->first;
00379                                                                                 
00380   std::map<double, std::vector<Point> >::iterator delBeg = fDeltaM32MapNormal[dmDV].begin();
00381   std::map<double, std::vector<Point> >::iterator delEnd = fDeltaM32MapNormal[dmDV].end();
00382   fFitPointChi2N->Reset();
00383   fFitPointChi2I->Reset();
00384                                                                                 
00385   while(delBeg != delEnd){
00386     for(unsigned int i = 0; i < delBeg->second.size(); i++){
00387        // at this point
00388        double delta = delBeg->first;
00389        double ss2th = delBeg->second[i].th13;
00390        double nexp  = delBeg->second[i].nExpected;
00391        double nsig =  delBeg->second[i].nSignal;
00392 
00393        double chi2 = 0;
00394        if(fFitMethod == 0) chi2 = CalculateDeltaLog(nexp, val);
00395        if(fFitMethod == 1) chi2 = CalculateDeltaLog(nexp, val, nsig);
00396 
00397        fFitPointChi2N->Fill(ss2th, delta, chi2);
00398     }
00399     delBeg++;
00400   }
00401 
00402   delBeg = fDeltaM32MapInvert[dmDV].begin();
00403   delEnd = fDeltaM32MapInvert[dmDV].end();
00404                                                                                 
00405   while(delBeg != delEnd){
00406     for(unsigned int i = 0; i < delBeg->second.size(); i++){
00407        // at this point
00408        double delta = delBeg->first;
00409        double ss2th = delBeg->second[i].th13;
00410        double nexp  = delBeg->second[i].nExpected;
00411        double nsig =  delBeg->second[i].nSignal;
00412 
00413        double chi2 = 0;
00414        if(fFitMethod == 0) chi2 = CalculateDeltaLog(nexp, val);
00415        if(fFitMethod == 1) chi2 = CalculateDeltaLog(nexp, val, nsig);
00416                                                                                 
00417        fFitPointChi2I->Fill(ss2th, delta, chi2);
00418     }
00419     delBeg++;
00420   }
00421 
00422   cout<<"Completed chi2 map"<<endl;
00423   return true;
00424 }

bool NueFitter::BuildContourMaps (  ) 

Definition at line 304 of file NueFitter.cxx.

References fContourHistsI, fContourHistsN, fDeltaM32MapInvert, and fDeltaM32MapNormal.

Referenced by PrepareFit().

00305 {
00306   //This fills the maps for the various chi2 levels
00307 
00308   double dmDV = fDeltaM32MapNormal.begin()->first;
00309 
00310   std::map<double, std::vector<Point> >::iterator delBeg = fDeltaM32MapNormal[dmDV].begin();
00311   std::map<double, std::vector<Point> >::iterator delEnd = fDeltaM32MapNormal[dmDV].end();
00312 
00313   while(delBeg != delEnd){     
00314     double delta = delBeg->first;
00315     for(unsigned int i = 0; i < delBeg->second.size(); i++){
00316        // at this point
00317        double ss2th = delBeg->second[i].th13;
00318   
00319        int nChi2 = delBeg->second[i].cutVals.size();
00320        for(int j = 0; j < nChi2; j++){
00321          double val = delBeg->second[i].cutVals[j].cutval;
00322          fContourHistsN[j]->Fill(ss2th, delta, val);
00323        }
00324     }    
00325     delBeg++;
00326   }
00327 
00328   dmDV = fDeltaM32MapInvert.begin()->first;
00329   delBeg = fDeltaM32MapInvert[dmDV].begin();
00330   delEnd = fDeltaM32MapInvert[dmDV].end();
00331                                                                                 
00332   while(delBeg != delEnd){
00333     for(unsigned int i = 0; i < delBeg->second.size(); i++){
00334        // at this point
00335        double delta = delBeg->first;
00336        double ss2th = delBeg->second[i].th13;
00337                                                                                 
00338        int nChi2 = delBeg->second[i].cutVals.size();
00339        for(int j = 0; j < nChi2; j++){
00340          double val = delBeg->second[i].cutVals[j].cutval;
00341          fContourHistsI[j]->Fill(ss2th, delta, val);
00342        }
00343     }
00344     delBeg++;
00345   }
00346 
00347   cout<<"Completed contour maps"<<endl;
00348   return true;
00349 }

double NueFitter::CalculateDeltaLog ( double  nexp,
double  nobs,
double  nsig 
) [private]

Definition at line 261 of file NueFitter.cxx.

References n, and Munits::rad.

00262 { 
00263     //First calculate the best fit to the particular point mu given these values of b,s,k
00264     // These are the same formulae used in the analytic feldman cousins
00265 
00266     double b = nexp - nsig;
00267     double s = nsig;
00268     double k = 1.0;
00269     double errK = 0.0774118*0.0774118;
00270     double errBg = b*b*0.073739*0.073739;
00271     double n = nobs;
00272 
00273     double sol_A = 1 + errK/errBg*s*s;
00274     double sol_B = errBg + s*k - 2*errK*s*s*b/errBg - b + s*s*errK;
00275     double sol_C = -n*errBg + s*k*errBg - s*s*errK*b - s*k*b + s*s*b*b*errK/errBg;
00276 
00277     double rad = sol_B*sol_B-4*sol_A*sol_C;
00278     if(rad < 0) std::cout<<"negative radical - not sure about this...."<<std::endl;
00279 
00280    double betaHat = (-sol_B + TMath::Sqrt(rad))/(2*sol_A);
00281    double kHat = k - s*errK/errBg*(b-betaHat);
00282 
00283    nexp = betaHat + s*kHat;
00284    double chisq = 2*(nexp - n + n*TMath::Log(n/nexp))
00285                              + (betaHat-b)*(betaHat-b)/errBg + (kHat - k)*(kHat-k)/errK;
00286 
00287    double chiBF = 0.0;
00288 
00289   if(n > b){
00290           chiBF = 0.0;
00291      }else{
00292        double  bBest = 0.5*(b - errBg + TMath::Sqrt( (b-errBg)*(b - errBg) + 4*n*errBg));
00293         // then we have no signal at best fit point
00294 
00295         chiBF = 2*(bBest - n + n*TMath::Log(n/bBest)) + (b-bBest)*(b-bBest)/errBg;
00296     }
00297 
00298 //   if(chiBF > chisq+1e-10)
00299 //       cout<<n<<"  "<<chisq<<"  "<<chiBF<<"  "<<bBest<<"  "<<betaHat<<"  "<<nexp<<"  "<<kHat<<"  "<<k<<"  "<<errBg<<"  "<<errK<<"  "<<s<<endl;
00300 
00301     return (chisq - chiBF);
00302 }

double NueFitter::CalculateDeltaLog ( double  nexp,
double  nobs 
) [private]

Definition at line 244 of file NueFitter.cxx.

References MuELoss::e, and fMinTotalEvents.

Referenced by BuildChi2Map().

00245 {
00246   //I assume that if the observation has a value greater than the minimum number  // of events, than it can be predicted perfectly
00247                                                                                 
00248   double chi2 = 2*(nexp - nobs + nobs*TMath::Log(nobs/nexp));
00249                                                                                 
00250   if(nobs < fMinTotalEvents)
00251     chi2 -= 2*(fMinTotalEvents - nobs + nobs*TMath::Log(nobs/fMinTotalEvents));
00252                                                                                 
00253   if(nexp+1e-4 < fMinTotalEvents){
00254     cout<<"The expected number is less then the minimum number of events - "
00255         <<" something is messed up, please check: min: "<<fMinTotalEvents<<" nexp: "<<nexp<<endl;
00256   }
00257      
00258   return chi2;
00259 }

bool NueFitter::FitInput ( double  input,
std::string  outName 
)
void NueFitter::Initialize (  ) 
bool NueFitter::InitializeFittingHistograms (  )  [private]

Definition at line 152 of file NueFitter.cxx.

References fContourHistsI, fContourHistsN, fContourLevels, fDeltaM32MapNormal, fFitPointChi2I, fFitPointChi2N, and one().

Referenced by PrepareFit().

00153 {
00154    double dmDV = fDeltaM32MapNormal.begin()->first;
00155 
00156    std::map<double, std::vector<Point> >::iterator delBeg = fDeltaM32MapNormal[dmDV].begin();
00157    std::map<double, std::vector<Point> >::iterator delEnd = fDeltaM32MapNormal[dmDV].end();
00158 
00159    vector<double> deltaPos;
00160    while(delBeg != delEnd){
00161      deltaPos.push_back(delBeg->first);
00162      delBeg++;
00163    }
00164    
00165    vector<double> sin22th13Pos;
00166    std::vector<Point>::iterator thBeg = fDeltaM32MapNormal[dmDV][0].begin();
00167    std::vector<Point>::iterator thEnd = fDeltaM32MapNormal[dmDV][0].end();
00168                                                                                 
00169    while(thBeg != thEnd){
00170      sin22th13Pos.push_back(thBeg->th13);
00171      thBeg++;
00172    }
00173 
00174    //Loaded in all the positions - now convert to array 
00175    // with offsets
00176    sort(deltaPos.begin(), deltaPos.end());
00177    sort(sin22th13Pos.begin(), sin22th13Pos.end());
00178 
00179    int nDelta = deltaPos.size();
00180    int nTh13 = sin22th13Pos.size();
00181 
00182    if(nDelta <= 1){
00183       cout<<"Only zero/one delta value??"<<std::endl;
00184       return false;
00185    }
00186    if(nTh13 <= 1){
00187       cout<<"Only zero/one theta value??"<<std::endl;
00188       return false;
00189    }
00190 
00191    double* deltaArray = new double[nDelta+1];
00192    double* sinthArray = new double[nTh13+1];   
00193 
00194    double offset = 0.0;
00195    for(int i = 1; i < nDelta; i++){
00196      double prev = deltaPos[i-1];
00197      double curr = deltaPos[i];
00198      offset = (curr - prev)/2.0;
00199      if(i == 1) deltaArray[0] = prev - offset;
00200      deltaArray[i] = curr - offset;
00201    }
00202    deltaArray[nDelta] = deltaArray[nDelta-1] + 2*offset;
00203  
00204    offset = 0.0;
00205    for(int i = 1; i < nTh13; i++){
00206      double prev = sin22th13Pos[i-1];
00207      double curr = sin22th13Pos[i];
00208      cout<<prev<<"  "<<curr<<endl;
00209      offset = (curr - prev)/2.0;
00210      if(i == 1) sinthArray[0] = prev - offset;
00211      sinthArray[i] = curr - offset;
00212    }
00213    sinthArray[nTh13] = sinthArray[nTh13-1] + 2*offset;
00214 
00215    TString n1 = "Chi2HistN";
00216    TString n2 = "Chi2HistI";
00217 
00218    fFitPointChi2N = new TH2D(n1, n1, nTh13, sinthArray, nDelta, deltaArray);
00219    fFitPointChi2I = new TH2D(n2, n2, nTh13, sinthArray, nDelta, deltaArray);
00220 
00221    fFitPointChi2N->SetDirectory(0);
00222    fFitPointChi2I->SetDirectory(0);   
00223    
00224    for(unsigned int i = 0; i < fContourLevels.size(); i++)
00225    {
00226       char num[3];
00227       sprintf(num, "%d", i);
00228       TString one = "Chi2Hist" + string(num) + "N";
00229       TString two = "Chi2Hist" + string(num) + "I";
00230 
00231       TH2D* temp = new TH2D(one, one, nTh13, sinthArray, nDelta, deltaArray);
00232       fContourHistsN.push_back(temp);
00233       temp  = new TH2D(two, two, nTh13, sinthArray, nDelta, deltaArray);
00234       fContourHistsI.push_back(temp);
00235 
00236       fContourHistsN[i]->SetDirectory(0);
00237       fContourHistsI[i]->SetDirectory(0);
00238    }
00239 
00240    cout<<"Finished Initializing Histograms"<<endl;
00241    return true;
00242 }

bool NueFitter::LoadFiles (  ) 

Definition at line 18 of file NueFitter.cxx.

References count, Chi2Cut::cutval, Point::cutVals, fContourLevels, fDeltaM32MapInvert, fDeltaM32MapNormal, fErrors, fFileList, fMinTotalEvents, Munits::m, Point::nExpected, Point::nSignal, Chi2Cut::percent, Point::th13, and total().

Referenced by PrepareFit().

00019 {
00020    for(unsigned int i = 0; i < fFileList.size(); i++)
00021    {
00022       TFile* file = new TFile(fFileList[i].c_str(), "READ");
00023       TTree* Pos;
00024       file->GetObject("PositionTree", Pos);
00025 
00026       double dm32, delta, th13, massH;
00027       char name[100];
00028  
00029       double nevent, onesigma, ninety, threesigma, signal;
00030    
00031       Pos->SetBranchAddress("Th13", &th13);
00032       Pos->SetBranchAddress("DeltaM32", &dm32);
00033       Pos->SetBranchAddress("DeltaCP",&delta);
00034       Pos->SetBranchAddress("Hierarchy",&massH);
00035       Pos->SetBranchAddress("ParamID", &name);
00036       Pos->SetBranchAddress("Total", &nevent);
00037       Pos->SetBranchAddress("SixtyEightCut", &onesigma);
00038       Pos->SetBranchAddress("NinetyCut", &ninety);
00039       Pos->SetBranchAddress("ThreeSigma", &threesigma);
00040       Pos->SetBranchAddress("NueCC", &signal);
00041  
00042       Point temp;
00043        
00044       for(int j = 0; j < Pos->GetEntries(); j++){ 
00045         // For each position in this particular file
00046         Pos->GetEntry(j);
00047         temp.cutVals.clear();
00048         temp.nExpected = nevent;
00049         temp.nSignal = signal;
00050         double thVal = TMath::Sin(2*th13);
00051         thVal *= thVal;        
00052 
00053         temp.th13 = thVal;
00054 
00055         for(unsigned int k = 0; k < fContourLevels.size(); k++){
00056            Chi2Cut tempChi; 
00057            tempChi.percent = fContourLevels[k]; tempChi.cutval = 0;
00058            if(fContourLevels[k] == 68){
00059              tempChi.percent = 68; tempChi.cutval = onesigma; 
00060            }
00061            if(fContourLevels[k] == 90){
00062              tempChi.percent = 90; tempChi.cutval = ninety;
00063            }
00064            if(fContourLevels[k] == 99.73){
00065              tempChi.percent = 99.73; tempChi.cutval = threesigma;
00066            }
00067 
00068            if(tempChi.cutval == 0){
00069             cout<<"Why am I here... "<<fContourLevels[k]<<endl;
00070             //Clever function that takes the deltaLog histogram
00071             // from the file and figures out the cut values 
00072             // for now its the Generator code, but not really 
00073             //efficient if multiple cuts here
00074             TH1D* hist = 0;
00075             
00076             TString hName = string(name) + "/DeltaLogLikely";
00077             file->GetObject(hName, hist);
00078 
00079             double total = hist->GetSum();
00080             double count = 0.0;
00081             for(int m = 0; m < hist->GetNbinsX(); m++){
00082               count += hist->GetBinContent(m);
00083               if(count/total > fContourLevels[k]/100.0){
00084                  tempChi.cutval = hist->GetBinLowEdge(m);
00085                  m = hist->GetNbinsX();
00086               }
00087             }
00088             delete hist;
00089            }
00090 
00091            temp.cutVals.push_back(tempChi);
00092         } //end of loop over contours 
00093 
00094         //store it and move on
00095 
00096         vector<Chi2Cut>(temp.cutVals).swap(temp.cutVals);
00097 
00098         if(massH > 0) 
00099            fDeltaM32MapNormal[dm32][delta].push_back(temp);
00100         else
00101            fDeltaM32MapInvert[dm32][delta].push_back(temp);
00102 
00103         if(temp.th13 < 0 || temp.th13> 1)
00104           cout<<"Pushing "<<dm32<<" "<<delta<<"  "<<temp.th13<<"  "<<massH<<endl;
00105 //        cout<<temp.cutVals.capacity()<<"  "<<fDeltaM32MapNormal[dm32][delta].capacity()<<"  "<<endl;
00106         if(j%10000 == 0) cout<<"Loaded "<<j<<" points..."<<endl;
00107      }
00108      
00109      //Check this file against the others
00110      TTree* fileInfo;
00111      file->GetObject("GenerationTree", fileInfo);
00112                                                          
00113      const int NERR = 6;
00114      double error[NERR];
00115      double eventMin;
00116                                                                                 
00117      fileInfo->SetBranchAddress("MinTotal", &eventMin);
00118      fileInfo->SetBranchAddress("NCErr", &error[0]);
00119      fileInfo->SetBranchAddress("NuMuCCErr",&error[1]);
00120      fileInfo->SetBranchAddress("NueCCErr",&error[2]);
00121      fileInfo->SetBranchAddress("NuTauCCErr", &error[3]); 
00122      fileInfo->SetBranchAddress("BNueCCErr", &error[4]);
00123      fileInfo->SetBranchAddress("TotalBGErr", &error[5]);
00124 
00125      fileInfo->GetEntry(0);
00126      if(fMinTotalEvents < 0){
00127        fMinTotalEvents = eventMin;
00128        for(int m = 0; m < NERR; m++) fErrors[m] = error[m];
00129      }else{
00130        bool good = true;
00131 //       good = good && (fMinTotalEvents == eventMin);
00132        if(fMinTotalEvents > eventMin) fMinTotalEvents = eventMin;
00133        cout<<fMinTotalEvents<<"  "<<eventMin<<endl;
00134        for(int m = 0; m < NERR; m++) {
00135           if(fErrors[m] > 0) cout<<fErrors[m]<<"  "<<error[m]<<endl;
00136           good = good && (fErrors[m] == error[m]);
00137        }
00138 
00139        if(!good)
00140          cout<<"Mismatch between files - check the GenerationTrees!"<<endl;
00141      }
00142   
00143      delete fileInfo;
00144      delete Pos;
00145      delete file;
00146    }
00147    std::cout<<"Finished Loading all files"<<std::endl;
00148 
00149    return true;
00150 }

bool NueFitter::PerformFit ( double  input,
std::string  outName = "dummy.root" 
)
bool NueFitter::PrepareFit (  ) 

Definition at line 357 of file NueFitter.cxx.

References BuildContourMaps(), InitializeFittingHistograms(), and LoadFiles().

00358 {
00359    LoadFiles();
00360    if(!InitializeFittingHistograms()) return false;
00361    BuildContourMaps();
00362 
00363    return true;
00364 }

void NueFitter::Reset (  ) 
void NueFitter::SetFitMethod ( int  meth  )  [inline]

Definition at line 48 of file NueFitter.h.

References fFitMethod.

00048 {fFitMethod = meth;};

void NueFitter::SetOutputFile ( std::string  out  )  [inline]

Definition at line 38 of file NueFitter.h.

References fOutFile.

00038 {fOutFile = out;};

void NueFitter::WriteToFile (  ) 

Definition at line 433 of file NueFitter.cxx.

References fContourHistsI, fContourHistsN, fContourLevels, fFitPointChi2I, fFitPointChi2N, and fOutFile.

00434 {
00435   TFile* file = 0;  
00436 
00437   file = new TFile(fOutFile.c_str(),"RECREATE");
00438   file->cd();
00439 
00440   fFitPointChi2N->Write(); 
00441   fFitPointChi2I->Write();
00442                                                                                 
00443   for(unsigned int j = 0; j < fContourLevels.size(); j++){
00444     fContourHistsN[j]->Write();
00445     fContourHistsI[j]->Write();
00446   }
00447   file->Close();
00448 }


Member Data Documentation

std::vector<TH2D*> NueFitter::fContourHistsI [private]

Definition at line 67 of file NueFitter.h.

Referenced by BuildContourMaps(), InitializeFittingHistograms(), and WriteToFile().

std::vector<TH2D*> NueFitter::fContourHistsN [private]

Definition at line 66 of file NueFitter.h.

Referenced by BuildContourMaps(), InitializeFittingHistograms(), and WriteToFile().

std::vector<double> NueFitter::fContourLevels [private]

Definition at line 56 of file NueFitter.h.

Referenced by AddContour(), InitializeFittingHistograms(), LoadFiles(), and WriteToFile().

std::map<double, std::map<double, std::vector<Point> > > NueFitter::fDeltaM32MapInvert [private]

Definition at line 59 of file NueFitter.h.

Referenced by BuildChi2Map(), BuildContourMaps(), and LoadFiles().

std::map<double, std::map<double, std::vector<Point> > > NueFitter::fDeltaM32MapNormal [private]

Definition at line 58 of file NueFitter.h.

Referenced by BuildChi2Map(), BuildContourMaps(), InitializeFittingHistograms(), and LoadFiles().

double NueFitter::fErrors[10] [private]

Definition at line 70 of file NueFitter.h.

Referenced by LoadFiles().

std::vector<std::string> NueFitter::fFileList [private]

Definition at line 55 of file NueFitter.h.

Referenced by AddFile(), and LoadFiles().

int NueFitter::fFitMethod [private]

Definition at line 74 of file NueFitter.h.

Referenced by BuildChi2Map(), NueFitter(), and SetFitMethod().

TH2D* NueFitter::fFitPointChi2I [private]

Definition at line 64 of file NueFitter.h.

Referenced by BuildChi2Map(), InitializeFittingHistograms(), and WriteToFile().

TH2D* NueFitter::fFitPointChi2N [private]

Definition at line 63 of file NueFitter.h.

Referenced by BuildChi2Map(), InitializeFittingHistograms(), and WriteToFile().

double NueFitter::fMinTotalEvents [private]

Definition at line 69 of file NueFitter.h.

Referenced by CalculateDeltaLog(), LoadFiles(), and NueFitter().

std::string NueFitter::fOutFile [private]

Definition at line 72 of file NueFitter.h.

Referenced by NueFitter(), SetOutputFile(), and WriteToFile().


The documentation for this class was generated from the following files:

Generated on 22 Nov 2017 for loon by  doxygen 1.6.1