NueFCSensitivity Class Reference

#include <NueFCSensitivity.h>

List of all members.

Public Member Functions

 NueFCSensitivity ()
 ~NueFCSensitivity ()
void Run (std::string input, std::string output, double outPOT=4.0)
void WriteToFile (std::string file)
float Oscillate (double dm23, double t13, double delta, int sign)
void RunFromGrid ()
void RunStandardApproach ()
void SetObserved (int in)
double MinimizationFunction (double *par)
void TestCode ()
void SetDeltaRange (double start, double end)
void SetTheta23 (double in)

Private Member Functions

void Initialize ()
void SetupGridRun ()
void LoadEventsFromFile ()
void GetPoint (int i, float &bg, float &sig)
float OscillateNumber (double Ue32)
float OscillateHist (double dm23, double Ue32, double delta, int sign)
float OscillateFile (double dm23, double Ue32, double delta, int sign)
float CalculateChi2 (int i, float bg, float sig)
float CalculateFitChi2 (int i, float bg, float sig)
int SetupChi2Histograms (NSCDataParam DPst13, NSCDataParam DPdelta, NSCDataParam DPdm23)
void ProduceTree (TTree *&infotree, TTree *&errtree)
double EvaluateOmega (double signal, double background, int i)
double CalculateRank (int n, double s, double b, double errBg, double k, double errK)
bool FindBestFit (int n, double s, double b, double errBg, double k, double errK, double *res)
double OscillateMatter (int nuFlavor, int nonOscNuFlavor, float Energy)
double OscillateMatter (int nuFlavor, int nonOscNuFlavor, float Energy, float dm2, float th13, float delta, int hierarchy)
void SetOscParamBase (float dm2, float ss13, float delta, int hierarchy)

Private Attributes

float signuenum
TH1D * signuehist
TH1D * reweight
double fTheta23
double * fBinCenter
std::vector< mininfomcInfo
double fMeasuredPOT
float currentVal [5]
float * fZeroValBg
float * fZeroValSig
float fBaseLine
int fMethod
int fNumConfig
NueSenseConfignsc
TH3D ** chi2n
TH3D ** chi2i
TH3D ** nsigN
TH3D ** nsigI
TH3D ** nbgN
TH3D ** nbgI
float t13Step
float dStep
float dm23Step
double fDeltaMS12
double fTh12
double fTh23
double fDensity
std::map< double, std::map<
double, std::vector< Point > > > 
fDeltaM32MapNormal
std::map< double, std::map<
double, std::vector< Point > > > 
fDeltaM32MapInvert
OscCalc fOscCalc
int fObserved
bool fObservedIsSet
double fResult
double fFixedFitPar [6]
double fDeltaStart
double fDeltaEnd


Detailed Description

Definition at line 21 of file NueFCSensitivity.h.


Constructor & Destructor Documentation

NueFCSensitivity::NueFCSensitivity (  ) 

Definition at line 47 of file NueFCSensitivity.cxx.

00047                                   {
00048   t13Step = 0.001;
00049   dStep = 0.02;
00050   dm23Step = 0.025;
00051   fBaseLine = 735;
00052   fObservedIsSet = false;
00053   fDeltaStart = -1;
00054   fDeltaEnd = 3;
00055   fTheta23 = -1e4;
00056   for(int i = 0; i < 6; i++)  fFixedFitPar[i] = 0.0;
00057 }

NueFCSensitivity::~NueFCSensitivity (  ) 

Definition at line 60 of file NueFCSensitivity.cxx.

00060                                    {
00061 }


Member Function Documentation

float NueFCSensitivity::CalculateChi2 ( int  i,
float  bg,
float  sig 
) [private]

Definition at line 952 of file NueFCSensitivity.cxx.

References err(), fZeroValBg, fZeroValSig, NueSenseConfig::GetErrorConfig(), and nsc.

00953 {
00954   NSCErrorParam err = nsc->GetErrorConfig(i);
00955 
00956   float nObs = sig+bg;
00957   float nZero = fZeroValBg[i] + fZeroValSig[i];
00958 
00959   double bgerr = err.bg_systematic*fZeroValBg[i];
00960   double sigerr = err.sig_systematic*fZeroValSig[i];
00961 
00962   float syserrsq = bgerr*bgerr + sigerr*sigerr;
00963  
00964   float errScale = nObs/(nObs + syserrsq);
00965                                                                                 
00966   float chi2 = 2*(nZero - nObs + nObs*TMath::Log(nObs/nZero)) * errScale;
00967   return chi2;
00968 }

float NueFCSensitivity::CalculateFitChi2 ( int  i,
float  bg,
float  sig 
) [private]

Definition at line 970 of file NueFCSensitivity.cxx.

References err(), fObserved, fObservedIsSet, fZeroValBg, fZeroValSig, NueSenseConfig::GetErrorConfig(), and nsc.

Referenced by RunStandardApproach().

00971 {
00972   NSCErrorParam err = nsc->GetErrorConfig(i);
00973                                                                                 
00974   float nZero = sig+bg;
00975   float nObs = fObserved;
00976   if(!fObservedIsSet) nObs = fZeroValBg[i] + fZeroValSig[i];
00977            
00978   double bgerr = err.bg_systematic*bg;
00979   double sigerr = err.sig_systematic*sig;
00980 
00981   float syserrsq = bgerr*bgerr + sigerr*sigerr;                                               
00982   float errScale = nZero/(nZero + syserrsq);
00983                                                                                 
00984   float chi2 = 2*(nZero - nObs + nObs*TMath::Log(nObs/nZero)) * errScale;
00985   return chi2;
00986 }

double NueFCSensitivity::CalculateRank ( int  n,
double  s,
double  b,
double  errBg,
double  k,
double  errK 
) [private]

Definition at line 462 of file NueFCSensitivity.cxx.

References FindBestFit().

Referenced by EvaluateOmega().

00463 {
00464    double results[3]; 
00465    //Calculate the numerator
00466    FindBestFit(n, s, b, errBg, k, errK, results);
00467 
00468 //   double numerator = results[2];
00469    
00470    double sBest, kBest, bBest;
00471 
00472    // Now calculate the denominator 
00473    if(n >= b) { // then we have n = sk+b, Beta = b; k = k0// don't have to do anything  
00474         sBest = (n-b)/k; kBest = k; bBest = b;
00475    }
00476    else{   //if(n < b) {
00477         bBest = 0.5*(b - errBg + TMath::Sqrt( (b-errBg)*(b - errBg) + 4*n*errBg)); 
00478         // then we have no signal at best fit point
00479         sBest = 0;  kBest = k;
00480     }
00481 
00482    double betaHat = results[0];
00483    double kHat = results[1];
00484    double logN = n*TMath::Log(s*kHat+betaHat) - (s*kHat+betaHat) - (betaHat-b)*(betaHat-b)/(2*errBg)
00485                       - (kHat-k)*(kHat-k)/(2*errK);
00486    double logD = n*TMath::Log(sBest*kBest+bBest) - (sBest*kBest+bBest) - (bBest-b)*(bBest-b)/(2*errBg)
00487                       - (kBest-k)*(kBest-k)/(2*errK);
00488 
00489     double rank = 1.0;
00490     if(logD != 0) { rank = TMath::Exp(logN-logD); }
00491     else std::cout<<"odd"<<std::endl;
00492 
00493     return rank;
00494 }

double NueFCSensitivity::EvaluateOmega ( double  signal,
double  background,
int  i 
) [private]

Definition at line 297 of file NueFCSensitivity.cxx.

References CalculateRank(), MuELoss::e, err(), FindBestFit(), fObserved, Gaussian(), NueSenseConfig::GetErrorConfig(), nsc, and Poisson().

Referenced by RunFromGrid().

00298 {
00299    int n0 = fObserved;
00300    double b0 = background;
00301    double s0 = signal;
00302    double k0 = 1.0;
00303 
00304    NSCErrorParam err = nsc->GetErrorConfig(i);
00305 
00306    double errK = 1.0;
00307    if(s0 > 0) errK = (err.sig_systematic)*(err.sig_systematic);
00308    if(errK < 1e-5){ errK = 1.0;  std::cout<<"can't do perfect signal"<<std::endl; }
00309    double errBg = err.bg_systematic*b0*err.bg_systematic*b0;
00310  
00311    double rank0 = CalculateRank(n0, s0, b0, errBg, k0, errK);
00312 
00313    double results[3];
00314    FindBestFit(n0, s0, b0, errBg, k0, errK, results);
00315 
00316    double betaHat = results[0];
00317    double kHat = results[1];
00318 
00319    errBg = (betaHat*err.bg_systematic)*(betaHat*err.bg_systematic);
00320    errK = (err.sig_systematic)*(err.sig_systematic)*kHat*kHat;
00321 
00322 //   std::cout<<n0<<"  "<<s0<<"  "<<b0<<"  "<<k0<<"  "<<betaHat<<" "<<kHat<<"  "<<errBg<<" "<<errK<<std::endl;
00323 //   std::cout<<s0*kHat + betaHat<<"  "<<n0<<std::endl;
00324 
00325    double omega = 0;
00326    double omegaBar = Poisson(s0*kHat + betaHat, n0);
00327 
00328    int nBase = int(s0*kHat + betaHat);
00329 
00330    int nHigh = nBase;
00331    int nLow  = nBase - 1;
00332 
00333    double delta = 1.0;
00334    bool filled = false;
00335 
00336    const int NUMB = 30;
00337    const int NUMK = 30;
00338 
00339    static double bVal[NUMB];
00340    static double kVal[NUMK];
00341    static double errBVar[NUMB];
00342    static double errKVar[NUMK];
00343 
00344    static bool first = true;
00345    static double scale = 1.0;
00346 
00347    double bStart = 0.4*b0;   double bStop = 2.2*b0;
00348    double kStart = 0.4;   double kStop = 1.6;
00349 
00350    if(first){
00351       for(int i = 0; i < NUMB; i++) { 
00352           bVal[i] = bStart + i*(bStop - bStart)/float(NUMB); 
00353           errBVar[i] = err.bg_systematic*err.bg_systematic*bVal[i]*bVal[i];
00354       }
00355   
00356       for(int i = 0; i < NUMK; i++) { 
00357           kVal[i] = kStart + i*(kStop - kStart)/float(NUMK);
00358           errKVar[i] = (err.sig_systematic)*(err.sig_systematic)*kVal[i]*kVal[i];
00359        }
00360       first = false; 
00361       scale = (bStop-bStart)*(kStop-kStart)/(NUMK*NUMB);
00362    }
00363    double bkProb[NUMB][NUMK];
00364 
00365    /*   so lets be clear about this, the values for these gaussians are the same
00366     but its the rank that changes for any given value of n (the values themselve differ with mu)
00367     so the first time through i calculate the contribution for each point in the space
00368     then when looping through just look it up in a giant array             */
00369 
00370    // better yet - I can save cycles by building the arrays separately on the first pass 
00371 
00372    double bGauss[NUMB];
00373    double kGauss[NUMK];
00374 
00375    for(int i = 0; i < NUMB; i++) { 
00376       bGauss[i] = Gaussian(betaHat, bVal[i], errBg);
00377    }
00378    for(int i = 0; i < NUMB; i++) { 
00379       kGauss[i] = Gaussian(kHat, kVal[i], errK);
00380    }    
00381 
00382    bool risingH = false, risingL = false;
00383    bool doneH = false;  // some calc savers
00384    bool doneL = false;  // some calc savers
00385 
00386    double ThreshHold = 1e-5;
00387    double pfThresh = ThreshHold*1e-2;
00388 
00389    double slip = 0;   // Error on the numerical integral
00390 
00391    while(delta > ThreshHold || (1 - (omega + omegaBar) > 2*ThreshHold) )
00392    {
00393       if(nHigh == n0) nHigh++;
00394       if(nLow == n0) nLow--;
00395 
00396       double dOmegaH = 0.0, dOmegaBarH = 0.0;
00397       double dOmegaL = 0.0, dOmegaBarL = 0.0;
00398 
00399       double PrefactorH = Poisson(s0*kHat+betaHat, nHigh);
00400       double PrefactorL = Poisson(s0*kHat+betaHat, nLow);
00401 
00402       if(PrefactorH > pfThresh) risingH = true;
00403       if(PrefactorL > pfThresh) risingL = true;
00404       if(PrefactorH < pfThresh && risingH) doneH = true;
00405       if(PrefactorL < pfThresh && risingL) doneL = true;
00406       if(doneH && doneL){
00407           // just a sanity check that the code doesn't think its done too early, this will cause an inf loop
00408           if(1 - (omega + omegaBar) < 1.5*slip) break;
00409           std::cout<<"Thats unexpected:  "<<1 - (omega + omegaBar)<<"  "<<slip<<std::endl;
00410       }
00411 
00412       if(PrefactorH > pfThresh || PrefactorL > pfThresh){  // No need to loop if contribution too small
00413         for(int i = 0; i < NUMB; i++) //  0 to infinity
00414         {
00415           double b = bVal[i];
00416           double eb = errBVar[i];
00417           for(int j = 0 ; j < NUMK; j++) //  -inf to inf
00418           {
00419               double k = kVal[j];
00420               if(!filled) bkProb[i][j] = bGauss[i]*kGauss[j];
00421  
00422               double val = bkProb[i][j];   
00423               double eK = errKVar[j];
00424               if(val < ThreshHold*1e-4) continue; // no point in using these contributions 
00425               
00426               double rank = 1.0;
00427 
00428               if(PrefactorH > pfThresh){
00429                 rank = CalculateRank(nHigh, s0, b, eb, k, eK);
00430                 if(rank > rank0) dOmegaH += val*scale;
00431                 else          dOmegaBarH += val*scale; 
00432               }
00433  
00434               if(PrefactorL > pfThresh){
00435                 if(nLow >= 0){
00436                   rank = CalculateRank(nLow, s0, b, eb, k, eK);
00437                   if(rank > rank0) dOmegaL += val*scale;
00438                   else          dOmegaBarL += val*scale; 
00439                 }
00440               }
00441           }
00442         }
00443         if(!filled) filled = true;          
00444         if(PrefactorH > pfThresh && 1 - (dOmegaH+ dOmegaBarH) > slip) slip = 1 - (dOmegaH+ dOmegaBarH);
00445         if(PrefactorL > pfThresh && 1 - (dOmegaL+ dOmegaBarL) > slip) slip = 1 - (dOmegaL+ dOmegaBarL);
00446       }
00447       delta = TMath::Max(PrefactorH*dOmegaH + PrefactorL*dOmegaL, PrefactorH*dOmegaBarH + PrefactorL*dOmegaBarL);
00448       omega += PrefactorH*dOmegaH + PrefactorL*dOmegaL;
00449       omegaBar += PrefactorH*dOmegaBarH + PrefactorL*dOmegaBarL;
00450       nHigh++; nLow--;
00451    }
00452 
00453    if(slip > ThreshHold){
00454      std::cout<<"Integral error past threshold:  "<<omega<<"  "<<omegaBar<<"  "<<omega+omegaBar<<"  "<<slip<<std::endl;
00455      std::cout<<n0<<"  "<<s0<<"  "<<b0<<"  "<<k0<<"  "<<betaHat<<" "<<kHat<<"  "<<errBg<<" "<<errK<<std::endl;
00456      std::cout<<s0*kHat + betaHat<<"  "<<n0<<std::endl;
00457    }
00458    return omega;
00459 }

bool NueFCSensitivity::FindBestFit ( int  n,
double  s,
double  b,
double  errBg,
double  k,
double  errK,
double *  res 
) [private]

Definition at line 496 of file NueFCSensitivity.cxx.

References Munits::rad.

Referenced by CalculateRank(), EvaluateOmega(), and TestCode().

00497 {
00498 //  At the moment this function is solvable in ~closed form (we ignore the dependance of errBg, errK on khat, betaHat
00499 //  this lets us avoid having a numerical minimization, but the second half of this function is the code that allows for that
00500 
00501    double sol_A = 1 + errK/errBg*s*s;
00502    double sol_B = errBg + s*k - 2*errK*s*s*b/errBg - b + s*s*errK;
00503    double sol_C = -n*errBg + s*k*errBg - s*s*errK*b - s*k*b + s*s*b*b*errK/errBg;
00504 
00505    double rad = sol_B*sol_B-4*sol_A*sol_C;
00506    if(rad < 0) std::cout<<"negative radical - not sure about this...."<<std::endl;
00507 
00508    double betaHat = (-sol_B + TMath::Sqrt(rad))/(2*sol_A);
00509 
00510    //solving for kHat given betaHat
00511    double kHat = k - s*errK/errBg*(b-betaHat);
00512 
00513    res[0] = betaHat;
00514    res[1] = kHat;
00515 
00516 /*
00517    fFixedFitPar[0] = n;
00518    fFixedFitPar[1] = s;
00519    fFixedFitPar[2] = b;
00520    fFixedFitPar[3] = errBg;
00521    fFixedFitPar[4] = k;
00522    fFixedFitPar[5] = errK;
00523 */
00524 //   res[2] = fResult = Poisson(s*kHat+betaHat, n, true)*Gaussian(betaHat, b, errBg)*Gaussian(kHat, k, errK);
00525 /*
00526    static bool first = true;
00527    static TMinuit *min = 0;
00528 
00529    if(first){    
00530       min = new TMinuit(2);
00531       gMinuit = min;
00532       min->SetFCN(WrapperFunction);
00533       min->SetObjectFit(this);
00534       min->DefineParameter(0,"bfit",fObserved,1, 0,60);
00535       min->DefineParameter(1,"kfit",1,0.99, 0,22);
00536 
00537       const double ERRDEF=1.;
00538       min->SetErrorDef(ERRDEF);
00539       min->SetPrintLevel(-1);
00540       first = false;
00541      std::cout<<"Built"<<std::endl;
00542    }
00543    Double_t arglist[10];
00544    Int_t ierflg = 0;
00545 
00546    arglist[0] = 10000;       //max calls
00547    arglist[1] = 0.01;         //tolerance
00548 
00549    min->mnexcm("SIMPLEX",arglist,2,ierflg);
00550 
00551    double errs[2];
00552    for(int i=0;i<2;i++){
00553      min->GetParameter(i,res[i],errs[i]);
00554 //     std::cout<<res[i]<<"  "<<errs[i]<<std::endl;
00555    }
00556 */
00557 //   res[2] = fResult;
00558 
00559    return true;
00560 }

void NueFCSensitivity::GetPoint ( int  i,
float &  bg,
float &  sig 
) [private]

Definition at line 732 of file NueFCSensitivity.cxx.

References NueConvention::bnue, currentVal, err(), NueSenseConfig::GetErrorConfig(), NueConvention::NC, nsc, NueConvention::nue, NueConvention::numu, and NueConvention::nutau.

Referenced by RunStandardApproach().

00733 {
00734    NSCErrorParam err = nsc->GetErrorConfig(i);  
00735 
00736    float bgt =  currentVal[ClassType::NC]  * (1.0 + err.scale[ClassType::NC])
00737              + currentVal[ClassType::numu]* (1.0 + err.scale[ClassType::numu])
00738              + currentVal[ClassType::bnue]* (1.0 + err.scale[ClassType::bnue])
00739              + currentVal[ClassType::nutau]*(1.0 + err.scale[ClassType::nutau]);
00740 
00741    bg = bgt;
00742    sig = currentVal[ClassType::nue]*(1.0 + err.scale[ClassType::nue]);
00743 }

void NueFCSensitivity::Initialize (  )  [private]

Definition at line 745 of file NueFCSensitivity.cxx.

References NueConvention::bnue, currentVal, MuELoss::e, fBinCenter, fDeltaMS12, fDensity, fMeasuredPOT, fMethod, fTh12, fTh23, NueSenseConfig::GetDeltaMS12(), NueSenseConfig::GetDensity(), NueSenseConfig::GetNueHistFile(), NueSenseConfig::GetNueHistName(), NueSenseConfig::GetNumber(), NueSenseConfig::GetOldDeltaMSquare(), NueSenseConfig::GetOldUe3Square(), NueSenseConfig::GetPOT(), NueSenseConfig::GetSinS2Th12(), NueSenseConfig::GetSinS2Th23(), OscPar::invKmToeV, NueSenseConfig::kAllNumbers, NueSenseConfig::kAnaNueFiles, NueSenseConfig::kNueHist, LoadEventsFromFile(), NueConvention::NC, nsc, NueConvention::nue, NueConvention::numu, NueConvention::nutau, reweight, NueSenseConfig::ShouldDeOsc(), signuehist, and signuenum.

Referenced by Run().

00746 {  
00747    float inPOT = nsc->GetPOT();
00748    float scale = fMeasuredPOT/inPOT;
00749 
00750    //Initilize the background
00751    if(fMethod == NueSenseConfig::kAllNumbers || fMethod == NueSenseConfig::kNueHist){
00752      currentVal[ClassType::NC] = nsc->GetNumber(ClassType::NC)*scale;
00753      currentVal[ClassType::numu] = nsc->GetNumber(ClassType::numu)*scale;
00754      currentVal[ClassType::bnue] = nsc->GetNumber(ClassType::bnue)*scale;
00755      currentVal[ClassType::nutau] = nsc->GetNumber(ClassType::nutau)*scale;
00756    }
00757 
00758    if(fMethod == NueSenseConfig::kAllNumbers){
00759      float nuenumber = nsc->GetNumber(ClassType::nue);
00760      if(nsc->ShouldDeOsc()){
00761         //Assumed originally Oscillated to:
00762         //  pow(TMath::Sin(theta23),2)*4.*UE32*(1-UE32)*pow(oscterm,2);
00763         float UE32 = nsc->GetOldUe3Square();
00764         float deoscweight = UE32*(1-UE32);
00765         signuenum = nuenumber/deoscweight*scale;
00766      }else{
00767         signuenum = nuenumber*scale;
00768      }
00769    }
00770 
00771    if(fMethod == NueSenseConfig::kNueHist){
00772       TFile inf(nsc->GetNueHistFile().c_str());
00773       TH1D *h = dynamic_cast<TH1D*>(inf.Get(nsc->GetNueHistName().c_str()));
00774       if(h == NULL){
00775         TH1D *h2 = dynamic_cast<TH1D*>(inf.Get(nsc->GetNueHistName().c_str()));
00776         if(h2 == NULL)
00777           std::cout<<"Failure to load signal histogram"<<std::endl;
00778         //This isn't really setup yet, but ideally I would convert it over to the right type
00779       }
00780 
00781       TH1D* hist = dynamic_cast<TH1D*>(h->Clone());
00782       hist->SetDirectory(0);
00783 
00784       if(nsc->ShouldDeOsc()){
00785         hist->Reset();
00786 
00787         float UE32 = nsc->GetOldUe3Square();
00788         float dms23 = nsc->GetOldDeltaMSquare();                                                
00789         for(int i = 0; i < h->GetNbinsX()+1; i++){
00790           float E = h->GetBinCenter(i);
00791           float num = h->GetBinContent(i);
00792 
00793           double invKmToeV = 1.97e-10; //convert 1/km to eV
00794           double Delta13Old = dms23*1e-3*735/(4.*E*1e9*invKmToeV);        
00795           double unweight = UE32*(1-UE32)*4*(1.0/2.0)
00796                              *TMath::Power(TMath::Sin(Delta13Old),2);
00797 
00798           hist->Fill(E, num/unweight*scale);
00799         }
00800      }else{
00801        hist->Scale(scale);
00802      }
00803 
00804      signuehist = dynamic_cast<TH1D*>(hist->Clone());
00805      signuehist->SetDirectory(0);
00806      reweight = dynamic_cast<TH1D*>(signuehist->Clone());
00807      reweight->SetDirectory(0);
00808 
00809      fBinCenter = new double[reweight->GetNbinsX()+1];
00810      for(int i = 0; i < reweight->GetNbinsX()+1; i++){
00811          fBinCenter[i] = reweight->GetBinCenter(i);
00812      }
00813      
00814    }
00815 
00816    if(fMethod == NueSenseConfig::kAnaNueFiles)  
00817       LoadEventsFromFile();
00818 
00819 
00820    fDeltaMS12 = nsc->GetDeltaMS12();
00821    float temp = nsc->GetSinS2Th12();
00822    fTh12 = TMath::ASin(TMath::Sqrt(temp))/2.;
00823    temp = nsc->GetSinS2Th23();
00824    fTh23 = TMath::ASin(TMath::Sqrt(temp))/2.;
00825    
00826    fDensity = nsc->GetDensity();
00827 
00828 }

void NueFCSensitivity::LoadEventsFromFile (  )  [private]

Definition at line 830 of file NueFCSensitivity.cxx.

References mininfo::energy, fMeasuredPOT, NueSenseConfig::GetFile(), NueSenseConfig::GetNumFiles(), mcInfo, nsc, mininfo::nuClass, mininfo::NuFlavor, mininfo::NuFlavorBeforeOsc, NuePOT::pot, and mininfo::weight.

Referenced by Initialize().

00831 {
00832   for(int i = 0; i < nsc->GetNumFiles(); i++){
00833    
00834      TChain *selected = new TChain("ana_nue");
00835      TChain *pots = new TChain("pottree");
00836     
00837      selected->Add(nsc->GetFile(i).c_str());
00838      pots->Add(nsc->GetFile(i).c_str());
00839  
00840      int nonOscFlavor, nuFlavor, type;
00841      float trueE;
00842      double skzpweight;        
00843 
00844      selected->SetMakeClass(1);
00845      selected->SetBranchStatus("*", 0);
00846      selected->SetBranchStatus("mctrue*", 1);
00847      selected->SetBranchStatus("fluxweight.totbeamweight", 1);
00848      selected->SetBranchAddress("mctrue.nonOscNuFlavor",&nonOscFlavor);
00849      selected->SetBranchAddress("mctrue.nuFlavor",&nuFlavor);
00850      selected->SetBranchAddress("mctrue.fNueClass",&type);
00851      selected->SetBranchAddress("mctrue.nuEnergy", &trueE);
00852      selected->SetBranchAddress("fluxweight.totbeamweight",&skzpweight);
00853         
00854      NuePOT *np;                                                             
00855      pots->SetBranchAddress("NuePOT", &np);
00856      double filePOT, scale;
00857      filePOT = 0.0;                                                                     
00858 
00859      for(int j=0; j < pots->GetEntries(); j++) {
00860          pots->GetEntry(j);
00861          filePOT += np->pot;
00862      }
00863      scale = filePOT/fMeasuredPOT;
00864         
00865      for(int j = 0; j < selected->GetEntries(); j++){
00866         selected->GetEntry(j);
00867         mininfo newInfo;
00868         newInfo.energy = trueE;
00869         newInfo.NuFlavor = nuFlavor;
00870         newInfo.NuFlavorBeforeOsc = nonOscFlavor;
00871         newInfo.nuClass = type;
00872         if(skzpweight < 0) skzpweight = 1.0;
00873         newInfo.weight = skzpweight*scale;
00874         mcInfo.push_back(newInfo);
00875      }  //End of this file
00876 
00877      delete selected;  
00878      delete pots; 
00879   }  //Done loading all files and POT Normalization is taken care of
00880 }

double NueFCSensitivity::MinimizationFunction ( double *  par  ) 

Definition at line 63 of file NueFCSensitivity.cxx.

References fFixedFitPar, and fResult.

Referenced by WrapperFunction().

00064 {
00065    int n = (int) fFixedFitPar[0];
00066    double s = fFixedFitPar[1];
00067    double b = fFixedFitPar[2];
00068    double sigBg = fFixedFitPar[3];
00069    double k = fFixedFitPar[4];
00070    double sigK = fFixedFitPar[5];
00071 
00072 //   std::cout<<n<<"  "<<b<<"  "<<s<<std::endl;
00073 
00074    double bFit = par[0];
00075    double kFit = par[1];
00076 
00077    double tot = s*kFit +bFit;
00078    if(tot < 0) std::cout<<s<<" "<<kFit<<"  "<<bFit<<" serious error "<<std::endl;
00079 
00080    double lnF = -n*TMath::Log(s*kFit + bFit) + (s*kFit + bFit) + 0.5*1/sigBg*(bFit-b)*(bFit-b);
00081           lnF +=  0.5*1/sigK*(kFit-k)*(kFit-k);
00082 
00083 //   std::cout<<"res: "<<kFit<<"  "<<bFit<<"  "<<lnF<<std::endl;
00084    fResult = lnF;
00085 
00086    return  lnF;
00087 }

float NueFCSensitivity::Oscillate ( double  dm23,
double  t13,
double  delta,
int  sign 
)

Definition at line 882 of file NueFCSensitivity.cxx.

References fMethod, NueSenseConfig::kAllNumbers, NueSenseConfig::kAnaNueFiles, NueSenseConfig::kNueHist, OscillateFile(), OscillateHist(), and OscillateNumber().

Referenced by RunStandardApproach().

00883 {
00884   if(fMethod == NueSenseConfig::kAllNumbers)
00885      return OscillateNumber(Ue32);
00886 
00887   if(fMethod == NueSenseConfig::kNueHist)
00888      return OscillateHist(dm23, Ue32, delta, sign);
00889    
00890   if(fMethod == NueSenseConfig::kAnaNueFiles)
00891      return OscillateFile(dm23, Ue32, delta, sign);
00892 
00893   return -1;
00894 }

float NueFCSensitivity::OscillateFile ( double  dm23,
double  Ue32,
double  delta,
int  sign 
) [private]

Definition at line 926 of file NueFCSensitivity.cxx.

References currentVal, mcInfo, NueConvention::nue, OscillateMatter(), and total().

Referenced by Oscillate().

00927 {
00928   float total[5];
00929 
00930   for(int i = 0; i < 5; i++) total[i] = 0;
00931 
00932   for(unsigned int i = 0; i < mcInfo.size(); i++){
00933     float E = mcInfo[i].energy;
00934     int inu = mcInfo[i].NuFlavor;
00935     int inunoosc = mcInfo[i].NuFlavorBeforeOsc;
00936     int nuClass = mcInfo[i].nuClass;
00937     double weight = mcInfo[i].weight;
00938 
00939     double oscweight = 
00940       OscillateMatter(inu,inunoosc,E,dm23,Ue32,delta,sign);
00941 
00942     total[nuClass] += weight*oscweight;
00943   }
00944 
00945   for(int i = 0; i < 5; i++){
00946     currentVal[i] = total[i];
00947   }
00948 
00949   return currentVal[ClassType::nue];
00950 }

float NueFCSensitivity::OscillateHist ( double  dm23,
double  Ue32,
double  delta,
int  sign 
) [private]

Definition at line 907 of file NueFCSensitivity.cxx.

References currentVal, fBinCenter, NueConvention::nue, OscillateMatter(), reweight, and signuehist.

Referenced by Oscillate().

00908 {
00909   reweight->Reset("ICE");
00910 
00911   for(int i = 0; i < signuehist->GetNbinsX()+1; i++){
00912      float E = fBinCenter[i];
00913      if(E > 10) break;
00914      float num = signuehist->GetBinContent(i);
00915      float weight = 1.0;
00916      if(num > 0)    weight =OscillateMatter(12,14,E); //,dm23,th13,delta,sign);
00917                                                                                 
00918      reweight->Fill(E, num*weight);
00919   }
00920  
00921   currentVal[ClassType::nue] = reweight->GetSum();
00922 //  std::cout<<dm23<<"  "<<Ue32<<"  "<<delta<<"  "<<sign<<"  "<<reweight->GetSum()<<std::endl;
00923   return reweight->GetSum();
00924 }                                                                                

double NueFCSensitivity::OscillateMatter ( int  nuFlavor,
int  nonOscNuFlavor,
float  Energy,
float  dm2,
float  th13,
float  delta,
int  hierarchy 
) [private]

Definition at line 1080 of file NueFCSensitivity.cxx.

References MuELoss::e, fBaseLine, fDeltaMS12, fDensity, fOscCalc, fTh12, fTh23, OscCalc::Oscillate(), and OscCalc::SetOscParam().

01083 {
01084  Double_t x[1] = {};
01085   x[0] = Energy;
01086   Double_t dm2_12 = fDeltaMS12*1e-5; //best fit SNO
01087   Double_t dm2_23 = dm2;
01088                                                                                 
01089                                                                                 
01090   Double_t par[9] = {0};
01091   par[0] = fBaseLine;
01092   par[1] = fTh23;
01093   par[2] = fTh12;
01094   par[3] = th13;
01095   par[4] = hierarchy*dm2_23;
01096   par[5] = dm2_12;
01097   par[6] = fDensity; //standard rock density
01098   par[7] = delta;
01099   par[8] = 1;
01100   if(nonOscNuFlavor < 0) par[8] = -1;
01101 
01102   //std::cout<<"About to call "<<dm2<<"  "<<ss13<<"  "<<delta<<std::endl;
01103   fOscCalc.SetOscParam(par);
01104   return fOscCalc.Oscillate(nuFlavor, nonOscNuFlavor, Energy);
01105 }

double NueFCSensitivity::OscillateMatter ( int  nuFlavor,
int  nonOscNuFlavor,
float  Energy 
) [private]

Definition at line 1107 of file NueFCSensitivity.cxx.

References fOscCalc, OscPar::kNuAntiNu, OscCalc::Oscillate(), and OscCalc::SetOscParam().

Referenced by OscillateFile(), and OscillateHist().

01109 {
01110 
01111   if(nonOscNuFlavor > 0) fOscCalc.SetOscParam(OscPar::kNuAntiNu, 1);
01112   if(nonOscNuFlavor < 0) fOscCalc.SetOscParam(OscPar::kNuAntiNu, -1);
01113 
01114   return fOscCalc.Oscillate(nuFlavor, nonOscNuFlavor, Energy);
01115 }

float NueFCSensitivity::OscillateNumber ( double  Ue32  )  [private]

Definition at line 896 of file NueFCSensitivity.cxx.

References currentVal, NueConvention::nue, and signuenum.

Referenced by Oscillate().

00897 {
00898 //  double Ue3 = TMath::Sin(t13);
00899 
00900   //Assumed originally Oscillated to:
00901   //  pow(TMath::Sin(theta23),2)*4.*UE32*(1-UE32)*pow(oscterm,2);
00902   currentVal[ClassType::nue] = signuenum*Ue32*(1-Ue32);
00903                                                                                            
00904   return currentVal[ClassType::nue];
00905 }

void NueFCSensitivity::ProduceTree ( TTree *&  infotree,
TTree *&  errtree 
) [private]

Definition at line 1138 of file NueFCSensitivity.cxx.

References NueConvention::bnue, MuELoss::e, NSCDataParam::end, err(), fBaseLine, fDensity, NueSenseConfig::GetDelta(), NueSenseConfig::GetDeltaMS12(), NueSenseConfig::GetDeltaMS23(), NueSenseConfig::GetErrorConfig(), NueSenseConfig::GetNumberConfig(), NueSenseConfig::GetSinS2Th12(), NueSenseConfig::GetSinS2Th13(), NueSenseConfig::GetSinS2Th23(), NueConvention::NC, nsc, NueConvention::nue, NueConvention::numu, NueConvention::nutau, and NSCDataParam::start.

Referenced by WriteToFile().

01139 {
01140   infotree = new TTree("OscPar","OscPar");
01141 
01142   double Th23, Th12, dm2_12;
01143   dm2_12 = nsc->GetDeltaMS12()*1e-5;
01144   Th12 = nsc->GetSinS2Th12();
01145   Th23 = nsc->GetSinS2Th23();
01146 
01147   infotree->Branch("L",&fBaseLine, "L/F");
01148   infotree->Branch("Sin2_2Theta23",&Th23, "Sin2_2Theta23/D");
01149   infotree->Branch("Sin2_2Theta12",&Th12, "Sin2_2Theta12/D");
01150   infotree->Branch("DeltaMS12",&dm2_12, "DeltaMS12/D");
01151   infotree->Branch("Density",&fDensity, "Density/D");
01152 
01153   NSCDataParam DPst13 = nsc->GetSinS2Th13();
01154   NSCDataParam DPdelta = nsc->GetDelta();
01155   NSCDataParam DPdm23 = nsc->GetDeltaMS23();
01156 
01157   infotree->Branch("DeltaMS23_start",&DPst13.start, "DeltaMS23_start/F");
01158   infotree->Branch("DeltaMS23_end",&DPst13.end, "DeltaMS23_end/F");
01159   infotree->Branch("Delta_start",&DPdelta.start, "Delta_start/F");
01160   infotree->Branch("Delta_end",&DPdelta.end, "Delta_end/F");
01161   infotree->Branch("Sin2_2Theta13_start",&DPdm23.start, "Sin2_2Theta13_start/F");
01162   infotree->Branch("Sin2_2Theta13_end",&DPdm23.end, "Sin2_2Theta13_end/F");
01163 
01164   infotree->Fill();
01165 
01166   errtree = new TTree("ErrPar","ErrPar");
01167   
01168   double bgsys, sigsys, nc, numu, nue, nutau, bnue;
01169   errtree->Branch("BG_systematic",&bgsys,"BG_systematic/D");
01170   errtree->Branch("SIG_systematic",&sigsys, "SIG_systematic/D");
01171   errtree->Branch("nc_scale",&nc, "nc_scale/D");
01172   errtree->Branch("numu_scale",&numu, "numu_scale/D");
01173   errtree->Branch("nue_scale",&nue, "nue_scale/D");
01174   errtree->Branch("nutau_scale",&nutau, "nutau_scale/D");
01175   errtree->Branch("bnue_scale",&bnue, "bnue_scale/D");
01176 
01177   for(int i = 0; i < nsc->GetNumberConfig(); i++){
01178       NSCErrorParam err = nsc->GetErrorConfig(i);
01179       bgsys = err.bg_systematic;
01180       sigsys = err.sig_systematic;
01181       nc = err.scale[ClassType::NC];
01182       numu = err.scale[ClassType::numu];
01183       bnue =err.scale[ClassType::bnue];
01184       nutau = err.scale[ClassType::nutau];
01185       nue = err.scale[ClassType::nue];
01186       
01187       errtree->Fill();
01188   }
01189 }

void NueFCSensitivity::Run ( std::string  input,
std::string  output,
double  outPOT = 4.0 
)

Definition at line 124 of file NueFCSensitivity.cxx.

References NueSenseConfig::CheckConfig(), fMeasuredPOT, fMethod, fNumConfig, NueSenseConfig::GetDataMethod(), NueSenseConfig::GetNumberConfig(), Initialize(), nsc, RunFromGrid(), RunStandardApproach(), and WriteToFile().

00125 {
00126    //Loading in the data
00127    nsc = new NueSenseConfig(input);
00128    if(!nsc->CheckConfig()) return;
00129    std::cout<<"Configuration file read and confirmed"<<std::endl;
00130 
00131    fMeasuredPOT = outPOT;
00132    // if just numbers roll forward, if taking form a chain pull out
00133    //  just enough information to do the oscillations, if hist grab it
00134 
00135    fMethod = nsc->GetDataMethod();
00136    fNumConfig = nsc->GetNumberConfig();
00137 
00138    Initialize();    //Load in any values and prepare listings
00139    std::cout<<"Initialization complete"<<std::endl;
00140 
00141    if(fMethod != 4)  RunStandardApproach();
00142    else              RunFromGrid();
00143  //Now we have all the data points
00144    // now we write it out to file
00145    WriteToFile(output);
00146 }

void NueFCSensitivity::RunFromGrid (  ) 

Definition at line 220 of file NueFCSensitivity.cxx.

References chi2i, chi2n, MuELoss::e, EvaluateOmega(), fDeltaEnd, fDeltaM32MapInvert, fDeltaM32MapNormal, fDeltaStart, fNumConfig, fTheta23, nbgI, nbgN, nsigI, nsigN, and SetupGridRun().

Referenced by Run().

00221 {
00222    SetupGridRun();
00223    std::cout<<"Setup is Complete"<<std::endl;  
00224  
00225    double dmDV = fDeltaM32MapNormal.begin()->first;
00226 
00227    std::map<double, std::vector<Point> >::iterator delBeg = fDeltaM32MapNormal[dmDV].begin();
00228    std::map<double, std::vector<Point> >::iterator delEnd = fDeltaM32MapNormal[dmDV].end();
00229 
00230    int ndCount = 0;
00231    double th23scale = 1.0;
00232 
00233    if(fTheta23 > -100) th23scale = 2*TMath::Sin(fTheta23)*TMath::Sin(fTheta23);
00234 
00235    while(delBeg != delEnd){
00236      double delta = delBeg->first;
00237      if(delta - fDeltaStart > -1e-4 && delta - fDeltaEnd < -1e-4){ 
00238        if(ndCount%10 == 0) std::cout<<"Running: "<<delta<<"  "<<fDeltaStart<<"  "<<fDeltaEnd<<std::endl;
00239 
00240        for(unsigned int i = 0; i < delBeg->second.size(); i++){
00241           // at this point
00242           double ss2th = delBeg->second[i].th13 * th23scale;
00243            
00244  
00245           double signal = delBeg->second[i].nsignal; 
00246           double bg = delBeg->second[i].nbg;
00247 
00248 //          if(ss2th < 0.426) continue;
00249 //          std::cout<<ss2th<<"  "<<signal<<"  "<<std::endl;
00250         // For each 
00251           for(int j = 0; j < fNumConfig; j++){
00252             double omega = EvaluateOmega(signal, bg, j);
00253 //          double val = CalculateFitChi2(j, bg, signal);
00254             chi2n[j]->Fill(ss2th, delta, dmDV, omega);
00255             nsigN[j]->Fill(ss2th, delta, dmDV, signal);
00256             nbgN[j]->Fill(ss2th, delta, dmDV, bg);
00257           }
00258        }
00259        ndCount++;
00260     }
00261     delBeg++;
00262   }
00263 
00264    std::cout<<"Finished Normal"<<std::endl;
00265 
00266   dmDV = fDeltaM32MapInvert.begin()->first;
00267   delBeg = fDeltaM32MapInvert[dmDV].begin();
00268   delEnd = fDeltaM32MapInvert[dmDV].end();
00269 
00270   ndCount = 0;
00271 
00272   while(delBeg != delEnd){
00273     double delta = delBeg->first;
00274     if(delta - fDeltaStart > -1e-4 && delta - fDeltaEnd < -1e-4){
00275       for(unsigned int i = 0; i < delBeg->second.size(); i++){
00276          // at this point
00277          double ss2th = delBeg->second[i].th13 * th23scale;
00278          double signal = delBeg->second[i].nsignal;
00279          double bg = delBeg->second[i].nbg;
00280 
00281          for(int j = 0; j < fNumConfig; j++){
00282             double omega =  EvaluateOmega(signal, bg, j);
00283             // CalculateFitChi2(j, bg, signal);
00284             chi2i[j]->Fill(ss2th, delta, dmDV, omega);
00285             nsigI[j]->Fill(ss2th, delta, dmDV, signal);
00286             nbgI[j]->Fill(ss2th, delta, dmDV, bg);
00287          }  
00288       }
00289       ndCount++;
00290       if(ndCount%10 == 0) std::cout<<"Finished inv delta = "<<delBeg->first<<std::endl;
00291     }
00292     delBeg++;
00293     
00294   }
00295 }

void NueFCSensitivity::RunStandardApproach (  ) 

Definition at line 148 of file NueFCSensitivity.cxx.

References CalculateFitChi2(), chi2i, chi2n, count, dm23Step, dStep, MuELoss::e, NSCDataParam::end, fMethod, fNumConfig, fOscCalc, fZeroValBg, fZeroValSig, NueSenseConfig::GetDelta(), NueSenseConfig::GetDeltaMS23(), GetPoint(), NueSenseConfig::GetSinS2Th13(), NSCDataParam::isfixed, NueSenseConfig::kAllNumbers, OscPar::kDelta, OscPar::kDeltaM23, OscPar::kTh13, nsc, Oscillate(), Mphysical::pi, OscCalc::SetOscParam(), SetOscParamBase(), SetupChi2Histograms(), NSCDataParam::start, and t13Step.

Referenced by Run().

00149 {
00150    NSCDataParam DPst13 = nsc->GetSinS2Th13();
00151    NSCDataParam DPdelta = nsc->GetDelta();
00152    NSCDataParam DPdm23 = nsc->GetDeltaMS23();
00153 
00154    fZeroValBg = new float[fNumConfig];
00155    fZeroValSig   = new float[fNumConfig];
00156 
00157    //int total = 
00158    SetupChi2Histograms(DPst13, DPdelta, DPdm23);
00159 
00160    SetOscParamBase(0.0024, 0.0, 0, 1);
00161 
00162    //Determine the number of background at Ue3 = 0
00163    Oscillate(0.0024, 0.0, 0, 1);
00164    for(int i = 0; i < fNumConfig; i++){
00165       fZeroValBg[i] = fZeroValSig[i] = 0.0;
00166       GetPoint(i, fZeroValBg[i], fZeroValSig[i]);
00167       std::cout<<"For configuration "<<i<<" starting with "
00168           <<fZeroValBg[i]<<", "<<fZeroValSig[i]<<" events."<<std::endl;
00169    } 
00170 
00171    //And now the big run over all the numbers
00172    bool d2,d3;  //some bools to keep the loops right
00173    float bg, sig;
00174    double pi = 3.1415926;
00175    int count = 0;
00176 
00177    for(double sins2t13 = DPst13.start; sins2t13 <= DPst13.end; sins2t13 += t13Step)
00178    {
00179       double Th13 = TMath::ASin(TMath::Sqrt(sins2t13))/2.;
00180       fOscCalc.SetOscParam(OscPar::kTh13, Th13);
00181 
00182       d2 = DPdelta.isfixed;
00183       for(double delta = DPdelta.start; delta <= DPdelta.end || d2 ; delta += dStep)
00184       {
00185          d3 = DPdm23.isfixed;
00186 
00187          fOscCalc.SetOscParam(OscPar::kDelta, delta*pi);
00188          for(double dm23 = DPdm23.start; dm23 <= DPdm23.end || d3; dm23+= dm23Step)
00189          {
00190             int sign = 1;
00191             fOscCalc.SetOscParam(OscPar::kDeltaM23, dm23*1e-3);
00192 
00193             Oscillate(dm23*1e-3, Th13, delta*pi, sign);   //Change the values
00194             for(int i = 0; i < fNumConfig; i++){
00195               GetPoint(i, bg, sig);
00196               float chi2 = CalculateFitChi2(i, bg, sig);
00197               chi2n[i]->Fill(sins2t13, delta, dm23, chi2);
00198             }
00199 
00200             if(fMethod != NueSenseConfig::kAllNumbers){
00201               sign = -1;
00202               fOscCalc.SetOscParam(OscPar::kDeltaM23, sign*dm23*1e-3);
00203               Oscillate(dm23*1e-3, Th13, delta*pi, sign);   //Change the values
00204 
00205               for(int i = 0; i < fNumConfig; i++){
00206                 GetPoint(i, bg, sig);
00207                 float chi2 = CalculateFitChi2(i, bg, sig);
00208                 chi2i[i]->Fill(sins2t13, delta, dm23, chi2);
00209               }
00210             }
00211             if(DPdm23.isfixed){ d3 = false;  dm23 = DPdm23.end + 1; }
00212             count++;
00213             if(count%100000 == 0) std::cout<<"On iteration "<<count<<std::endl;
00214          }//end of dm23
00215          if(DPdelta.isfixed){ d2 = false;  delta = DPdelta.end + 1; }
00216       }//end of delta
00217    }//end of sins2theta13
00218 }

void NueFCSensitivity::SetDeltaRange ( double  start,
double  end 
)

Definition at line 118 of file NueFCSensitivity.cxx.

References MuELoss::e, fDeltaEnd, and fDeltaStart.

00119 {
00120    if(start > -1e-3) fDeltaStart = start;
00121    if(end > -1e-3)   fDeltaEnd = end;
00122 }

void NueFCSensitivity::SetObserved ( int  in  )  [inline]

Definition at line 34 of file NueFCSensitivity.h.

References fObserved, and fObservedIsSet.

00034 {fObserved = in; fObservedIsSet = true;}

void NueFCSensitivity::SetOscParamBase ( float  dm2,
float  ss13,
float  delta,
int  hierarchy 
) [private]

Definition at line 1117 of file NueFCSensitivity.cxx.

References MuELoss::e, fBaseLine, fDeltaMS12, fDensity, fOscCalc, fTh12, fTh23, OscPar::kDelta, OscPar::kDeltaM12, OscPar::kDeltaM23, OscPar::kDensity, OscPar::kL, OscPar::kNuAntiNu, OscPar::kTh12, OscPar::kTh13, OscPar::kTh23, and OscCalc::SetOscParam().

Referenced by RunStandardApproach().

01118                                                                {
01119 
01120   Double_t dm2_12 = fDeltaMS12*1e-5; //best fit SNO
01121   Double_t dm2_23 = dm2;
01122                                                                                 
01123   Double_t par[9] = {0};
01124   par[OscPar::kL] = fBaseLine;
01125   par[OscPar::kTh23] = fTh23;
01126   par[OscPar::kTh12] = fTh12;
01127   par[OscPar::kTh13] = ss13; // TMath::ASin(TMath::Sqrt(ss2th13))/2.;
01128   par[OscPar::kDeltaM23] = hierarchy*dm2_23;
01129   par[OscPar::kDeltaM12] = dm2_12;
01130   par[OscPar::kDensity] = fDensity; //standard rock density
01131   par[OscPar::kDelta] = delta;
01132   par[OscPar::kNuAntiNu] = 1;
01133                                                                                 
01134 //  std::cout<<"About to call "<<dm2<<"  "<<ss13<<"  "<<delta<<std::endl;
01135   fOscCalc.SetOscParam(par);
01136 }

void NueFCSensitivity::SetTheta23 ( double  in  )  [inline]

Definition at line 40 of file NueFCSensitivity.h.

References fTheta23.

00040 {fTheta23 = in; }

int NueFCSensitivity::SetupChi2Histograms ( NSCDataParam  DPst13,
NSCDataParam  DPdelta,
NSCDataParam  DPdm23 
) [private]

Definition at line 989 of file NueFCSensitivity.cxx.

References chi2i, chi2n, dm23Step, dStep, NSCDataParam::end, fNumConfig, NSCDataParam::isfixed, nbgI, nbgN, nsigI, nsigN, NSCDataParam::start, and t13Step.

Referenced by RunStandardApproach().

00990 {
00991    chi2n = new TH3D*[fNumConfig];
00992    chi2i = new TH3D*[fNumConfig];
00993 
00994    nsigN = new TH3D*[fNumConfig];
00995    nsigI = new TH3D*[fNumConfig];
00996    nbgN = new TH3D*[fNumConfig];
00997    nbgI = new TH3D*[fNumConfig];
00998                                                                                            
00999    int t13Bins = int((DPst13.end - DPst13.start)/t13Step) + 1;
01000    int delBins = int((DPdelta.end - DPdelta.start)/dStep) + 1;
01001    int dm23Bins = int((DPdm23.end - DPdm23.start)/dm23Step) + 1;
01002    
01003    double xstart, xend, ystart, yend, zstart, zend;
01004    xstart = xend = ystart = yend = zstart = zend = 0;
01005                                                                                         
01006    if(DPst13.isfixed){ // i have no idea why you are running the code
01007    }else{
01008       xstart = DPst13.start - t13Step/2; xend = xstart + t13Bins*t13Step;
01009    }
01010    if(DPdelta.isfixed){
01011       ystart = DPdelta.start - 3*dStep/2;
01012       delBins = 5;
01013    }else{
01014       ystart = DPdelta.start - dStep/2;
01015    }
01016    yend = ystart + delBins*dStep;
01017 
01018    if(DPdm23.isfixed){
01019       zstart = DPdm23.start - 3*dm23Step/2;
01020       dm23Bins = 5;
01021    }else{
01022       zstart = DPdm23.start - dm23Step/2;
01023    }
01024    zend = zstart + dm23Bins*dm23Step;
01025 
01026    for(int i = 0; i < fNumConfig; i++){
01027       char temp[20];
01028       sprintf(temp, "chi2normal%d", i);
01029       chi2n[i] = new TH3D(temp, temp, t13Bins, xstart, xend, delBins,
01030                              ystart, yend, dm23Bins, zstart, zend);
01031       sprintf(temp, "chi2inverted%d", i);
01032       chi2i[i] = new TH3D(temp, temp, t13Bins, xstart, xend, delBins,
01033                              ystart, yend, dm23Bins, zstart, zend);
01034 
01035      chi2n[i]->SetDirectory(0);
01036      chi2i[i]->SetDirectory(0);
01037 
01038      sprintf(temp, "signormal%d", i);
01039      nsigN[i] = new TH3D(temp, temp, t13Bins, xstart, xend, delBins,
01040                              ystart, yend, dm23Bins, zstart, zend);
01041 
01042      sprintf(temp, "siginverted%d", i);
01043      nsigI[i] = new TH3D(temp, temp, t13Bins, xstart, xend, delBins,
01044                              ystart, yend, dm23Bins, zstart, zend);
01045 
01046      sprintf(temp, "bgnormal%d", i);
01047      nbgN[i] = new TH3D(temp, temp, t13Bins, xstart, xend, delBins,
01048                              ystart, yend, dm23Bins, zstart, zend);
01049  
01050      sprintf(temp, "bginverted%d", i); 
01051      nbgI[i] = new TH3D(temp, temp, t13Bins, xstart, xend, delBins,
01052                              ystart, yend, dm23Bins, zstart, zend);
01053    }
01054    std::cout<<"Chi2 Histograms created ("<<t13Bins<<", "<<delBins<<", "<<dm23Bins<<"): " 
01055             <<t13Bins*delBins*dm23Bins<<" iterations to perform"<<std::endl;
01056    return t13Bins*delBins*dm23Bins;
01057 }

void NueFCSensitivity::SetupGridRun (  )  [private]

Definition at line 562 of file NueFCSensitivity.cxx.

References chi2i, chi2n, dummy, err(), fDeltaM32MapInvert, fDeltaM32MapNormal, fMeasuredPOT, fNumConfig, fTheta23, NueSenseConfig::GetDataFiles(), NueSenseConfig::GetErrorConfig(), NueSenseConfig::GetPOT(), nbgI, nbgN, nsc, nsigI, and nsigN.

Referenced by RunFromGrid().

00563 {
00564   std::vector<std::string> datafiles = nsc->GetDataFiles();
00565   float inPOT = nsc->GetPOT();
00566   float scale = fMeasuredPOT/inPOT;
00567 
00568   NSCErrorParam err = nsc->GetErrorConfig(0);
00569   bool first = true;
00570 
00571   double hold = 2.32e-3;
00572 
00573   for(unsigned int i = 0; i < datafiles.size(); i++){
00574      std::cout<<"Openning file: "<<datafiles[i]<<std::endl;
00575      double th13, delta, mass, ni, bg[5], sig, dummy;
00576      std::ifstream ins(datafiles[i].c_str());
00577 
00578      ins>>th13>>delta>>mass>>ni>>bg[0]>>bg[1]>>bg[2]>>bg[3]>>sig>>dummy;
00579      hold = mass;
00580      while(!ins.eof()){
00581        Point temp;
00582        double temp2 = TMath::Sin(2*th13);
00583 
00584 //       bg[0] *=  1.0 + err.scale[ClassType::NC];
00585 //       bg[1] *=  1.0 + err.scale[ClassType::numu];
00586 //       bg[2] *=  1.0 + err.scale[ClassType::bnue];
00587 //       bg[3] *=  1.0 + err.scale[ClassType::nutau];
00588 //       sig *=  1.0 + err.scale[ClassType::nue];
00589        dummy = bg[0]+bg[1]+bg[2]+bg[3]+sig;
00590        temp.th13 = temp2*temp2;
00591        temp.nsignal = sig*scale;
00592        temp.nbg = (dummy-sig)*scale;
00593  
00594        if(i == 0 && first){ std::cout<< temp.nbg<<std::endl; first = false; }
00595 
00596        if(ni > 0)  fDeltaM32MapNormal[mass][delta].push_back(temp);
00597        else        fDeltaM32MapInvert[mass][delta].push_back(temp);
00598 
00599        ins>>th13>>delta>>mass>>ni>>bg[0]>>bg[1]>>bg[2]>>bg[3]>>sig>>dummy;
00600      }
00601      ins.clear();
00602    }
00603    
00604 
00605    std::map<double, std::vector<Point> >::iterator delBeg = fDeltaM32MapNormal[hold].begin();
00606    std::map<double, std::vector<Point> >::iterator delEnd = fDeltaM32MapNormal[hold].end();
00607 
00608    std::vector<double> deltaM2Pos;
00609    std::map<double, std::map<double, std::vector<Point> > >::iterator dmBeg = fDeltaM32MapNormal.begin();
00610    std::map<double, std::map<double, std::vector<Point> > >::iterator dmEnd = fDeltaM32MapNormal.end();
00611 
00612    while(dmBeg != dmEnd){
00613      deltaM2Pos.push_back(dmBeg->first);
00614      dmBeg++;
00615    }
00616 
00617    std::vector<double> deltaPos;
00618    while(delBeg != delEnd){
00619      deltaPos.push_back(delBeg->first);
00620      delBeg++;
00621    }
00622 
00623    std::vector<double> sin22th13Pos;
00624    std::vector<Point>::iterator thBeg = fDeltaM32MapNormal[hold][0].begin();
00625    std::vector<Point>::iterator thEnd = fDeltaM32MapNormal[hold][0].end();
00626 
00627    while(thBeg != thEnd){
00628      sin22th13Pos.push_back(thBeg->th13);
00629      thBeg++;
00630    }
00631 
00632    sort(deltaM2Pos.begin(), deltaM2Pos.end());
00633    sort(deltaPos.begin(), deltaPos.end());
00634    sort(sin22th13Pos.begin(), sin22th13Pos.end());
00635 
00636    int nDM2 = deltaM2Pos.size();
00637    int nDelta = deltaPos.size();
00638    int nTh13 = sin22th13Pos.size();
00639 
00640    if(nDelta <= 1){
00641       std::cout<<"Only zero/one delta value??"<<std::endl;
00642       return;
00643    }
00644    if(nTh13 <= 1){
00645       std::cout<<"Only zero/one theta value??"<<std::endl;
00646       return;
00647    }
00648 
00649    double* deltaM2Array = new double[nDM2+1];
00650    double* deltaArray = new double[nDelta+1];
00651    double* sinthArray = new double[nTh13+1];
00652 
00653    double offset = 0.0;
00654    for(int i = 1; i < nDelta; i++){
00655      double prev = deltaPos[i-1];
00656      double curr = deltaPos[i];
00657      offset = (curr - prev)/2.0;
00658      if(i == 1) deltaArray[0] = prev - offset;
00659      deltaArray[i] = curr - offset;
00660    }
00661    deltaArray[nDelta] = deltaArray[nDelta-1] + 2*offset;
00662 
00663    offset = 0.0;
00664    double th23scale = 1.0;
00665    if(fTheta23 > -100) th23scale = 2*TMath::Sin(fTheta23)*TMath::Sin(fTheta23);
00666 
00667    for(int i = 1; i < nTh13; i++){
00668      double prev = sin22th13Pos[i-1]*th23scale;
00669      double curr = sin22th13Pos[i]*th23scale;
00670      offset = (curr - prev)/2.0;
00671      if(i == 1) sinthArray[0] = prev - offset;
00672      sinthArray[i] = curr - offset;
00673    }
00674    sinthArray[nTh13] = sinthArray[nTh13-1] + 2*offset;
00675 
00676    offset = 0.00002;
00677    deltaM2Array[0] = deltaM2Pos[0] - 0.00001;
00678    deltaM2Array[1] = deltaM2Pos[0] + 0.00001; 
00679 /*   for(int i = 1; i < nDM2; i++){
00680      double prev = deltaM2Pos[i-1];
00681      double curr = deltaM2Pos[i];
00682      offset = (curr - prev)/2.0;
00683      if(i == 1) deltaM2Array[0] = prev - offset;
00684      deltaM2Array[i] = curr - offset;
00685    }*/ 
00686 //   deltaM2Array[nDM2] = deltaM2Array[0] + 2*offset;
00687 
00688    TString n1 = "Chi2HistN";
00689    TString n2 = "Chi2HistI";
00690 
00691    chi2n = new TH3D*[fNumConfig];
00692    chi2i = new TH3D*[fNumConfig];
00693    nsigN = new TH3D*[fNumConfig];
00694    nsigI = new TH3D*[fNumConfig];
00695    nbgN = new TH3D*[fNumConfig];
00696    nbgI = new TH3D*[fNumConfig];  
00697 
00698    for(int i = 0; i < fNumConfig; i++){
00699      char temp[20];
00700      sprintf(temp, "chi2normal%d", i);
00701      chi2n[i] = new TH3D(temp, temp, nTh13, sinthArray, nDelta, deltaArray, nDM2, deltaM2Array);
00702      sprintf(temp, "chi2inverted%d", i);
00703      chi2i[i] = new TH3D(temp, temp, nTh13, sinthArray, nDelta, deltaArray, nDM2, deltaM2Array);
00704  
00705      chi2n[i]->SetDirectory(0);
00706      chi2i[i]->SetDirectory(0);
00707      sprintf(temp, "signormal%d", i);
00708      nsigN[i] = new TH3D(temp, temp, nTh13, sinthArray, nDelta, deltaArray, nDM2, deltaM2Array);
00709 
00710      sprintf(temp, "siginverted%d", i);
00711      nsigI[i] = new TH3D(temp, temp, nTh13, sinthArray, nDelta, deltaArray, nDM2, deltaM2Array);
00712 
00713      sprintf(temp, "bgnormal%d", i);
00714      nbgN[i] = new TH3D(temp, temp, nTh13, sinthArray, nDelta, deltaArray, nDM2, deltaM2Array);
00715 
00716      sprintf(temp, "bginverted%d", i);
00717      nbgI[i] = new TH3D(temp, temp, nTh13, sinthArray, nDelta, deltaArray, nDM2, deltaM2Array);
00718 
00719      nsigN[i]->SetDirectory(0);
00720      nsigI[i]->SetDirectory(0);
00721      nbgN[i]->SetDirectory(0);
00722      nbgI[i]->SetDirectory(0);
00723    }
00724    std::cout<<"Chi2 Histograms created ("<<nTh13<<", "<<nDelta<<", "<<nDM2<<"): Running Grid Method"<<std::endl;
00725 
00726    //So now I have built all the histograms as appropriate 
00727    // - now I just loop over the data
00728 
00729 }

void NueFCSensitivity::TestCode (  ) 

Definition at line 89 of file NueFCSensitivity.cxx.

References FindBestFit(), fObserved, and min.

00090 {
00091 
00092    double results[3];
00093    fObserved = 35;
00094    FindBestFit(35, 1, 26.5, 4, 1, 0.07, results);
00095 
00096    std::cout<<results[0]<<"  "<<results[1]<<"  "<<results[2]<<std::endl;
00097   
00098    TMinuit* min = gMinuit;
00099    min->DefineParameter(0,"bfit",35,0.01, 0,60);
00100    min->DefineParameter(1,"kfit",1,0.01, 0,22);
00101 
00102    FindBestFit(35, 1, 26.5, 4, 1, 0.07, results);
00103 
00104    std::cout<<results[0]<<"  "<<results[1]<<"  "<<results[2]<<std::endl;
00105 
00106    min->DefineParameter(0,"bfit",35,3, 0,60);
00107    min->DefineParameter(1,"kfit",1,3, 0,22);
00108  
00109    FindBestFit(35, 1, 26.5, 4, 1, 0.07, results);
00110  
00111    std::cout<<results[0]<<"  "<<results[1]<<"  "<<results[2]<<std::endl;
00112 
00113 
00114 
00115 }

void NueFCSensitivity::WriteToFile ( std::string  file  ) 

Definition at line 1059 of file NueFCSensitivity.cxx.

References chi2i, chi2n, err(), fNumConfig, nbgI, nbgN, nsigI, nsigN, and ProduceTree().

Referenced by Run().

01060 {
01061    TFile out(file.c_str(), "RECREATE");
01062    out.cd();
01063 
01064    for(int i = 0; i < fNumConfig; i++){
01065       chi2n[i]->Write();
01066       chi2i[i]->Write();
01067       nsigN[i]->Write();
01068       nbgN[i]->Write();
01069       nsigI[i]->Write();
01070       nbgI[i]->Write();
01071    }
01072    TTree *info, *err;
01073 
01074    ProduceTree(info, err);
01075    info->Write();
01076    err->Write();
01077    out.Close();
01078 }


Member Data Documentation

TH3D** NueFCSensitivity::chi2i [private]

Definition at line 95 of file NueFCSensitivity.h.

Referenced by RunFromGrid(), RunStandardApproach(), SetupChi2Histograms(), SetupGridRun(), and WriteToFile().

TH3D** NueFCSensitivity::chi2n [private]

Definition at line 94 of file NueFCSensitivity.h.

Referenced by RunFromGrid(), RunStandardApproach(), SetupChi2Histograms(), SetupGridRun(), and WriteToFile().

float NueFCSensitivity::currentVal[5] [private]

Definition at line 79 of file NueFCSensitivity.h.

Referenced by GetPoint(), Initialize(), OscillateFile(), OscillateHist(), and OscillateNumber().

float NueFCSensitivity::dm23Step [private]

Definition at line 104 of file NueFCSensitivity.h.

Referenced by RunStandardApproach(), and SetupChi2Histograms().

float NueFCSensitivity::dStep [private]

Definition at line 103 of file NueFCSensitivity.h.

Referenced by RunStandardApproach(), and SetupChi2Histograms().

float NueFCSensitivity::fBaseLine [private]

Definition at line 88 of file NueFCSensitivity.h.

Referenced by OscillateMatter(), ProduceTree(), and SetOscParamBase().

double* NueFCSensitivity::fBinCenter [private]

Definition at line 74 of file NueFCSensitivity.h.

Referenced by Initialize(), and OscillateHist().

double NueFCSensitivity::fDeltaEnd [private]

Definition at line 124 of file NueFCSensitivity.h.

Referenced by RunFromGrid(), and SetDeltaRange().

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

Definition at line 113 of file NueFCSensitivity.h.

Referenced by RunFromGrid(), and SetupGridRun().

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

Definition at line 112 of file NueFCSensitivity.h.

Referenced by RunFromGrid(), and SetupGridRun().

double NueFCSensitivity::fDeltaMS12 [private]

Definition at line 107 of file NueFCSensitivity.h.

Referenced by Initialize(), OscillateMatter(), and SetOscParamBase().

double NueFCSensitivity::fDeltaStart [private]

Definition at line 123 of file NueFCSensitivity.h.

Referenced by RunFromGrid(), and SetDeltaRange().

double NueFCSensitivity::fDensity [private]

Definition at line 110 of file NueFCSensitivity.h.

Referenced by Initialize(), OscillateMatter(), ProduceTree(), and SetOscParamBase().

double NueFCSensitivity::fFixedFitPar[6] [private]

Definition at line 121 of file NueFCSensitivity.h.

Referenced by MinimizationFunction().

double NueFCSensitivity::fMeasuredPOT [private]

Definition at line 77 of file NueFCSensitivity.h.

Referenced by Initialize(), LoadEventsFromFile(), Run(), and SetupGridRun().

int NueFCSensitivity::fMethod [private]

Definition at line 90 of file NueFCSensitivity.h.

Referenced by Initialize(), Oscillate(), Run(), and RunStandardApproach().

int NueFCSensitivity::fNumConfig [private]

Definition at line 91 of file NueFCSensitivity.h.

Referenced by Run(), RunFromGrid(), RunStandardApproach(), SetupChi2Histograms(), SetupGridRun(), and WriteToFile().

int NueFCSensitivity::fObserved [private]

Definition at line 117 of file NueFCSensitivity.h.

Referenced by CalculateFitChi2(), EvaluateOmega(), SetObserved(), and TestCode().

bool NueFCSensitivity::fObservedIsSet [private]

Definition at line 118 of file NueFCSensitivity.h.

Referenced by CalculateFitChi2(), and SetObserved().

OscCalc NueFCSensitivity::fOscCalc [private]

Definition at line 115 of file NueFCSensitivity.h.

Referenced by OscillateMatter(), RunStandardApproach(), and SetOscParamBase().

double NueFCSensitivity::fResult [private]

Definition at line 120 of file NueFCSensitivity.h.

Referenced by MinimizationFunction().

double NueFCSensitivity::fTh12 [private]

Definition at line 108 of file NueFCSensitivity.h.

Referenced by Initialize(), OscillateMatter(), and SetOscParamBase().

double NueFCSensitivity::fTh23 [private]

Definition at line 109 of file NueFCSensitivity.h.

Referenced by Initialize(), OscillateMatter(), and SetOscParamBase().

double NueFCSensitivity::fTheta23 [private]

Definition at line 73 of file NueFCSensitivity.h.

Referenced by RunFromGrid(), SetTheta23(), and SetupGridRun().

float* NueFCSensitivity::fZeroValBg [private]

Definition at line 82 of file NueFCSensitivity.h.

Referenced by CalculateChi2(), CalculateFitChi2(), and RunStandardApproach().

float* NueFCSensitivity::fZeroValSig [private]

Definition at line 83 of file NueFCSensitivity.h.

Referenced by CalculateChi2(), CalculateFitChi2(), and RunStandardApproach().

std::vector<mininfo> NueFCSensitivity::mcInfo [private]

Definition at line 76 of file NueFCSensitivity.h.

Referenced by LoadEventsFromFile(), and OscillateFile().

TH3D** NueFCSensitivity::nbgI [private]

Definition at line 100 of file NueFCSensitivity.h.

Referenced by RunFromGrid(), SetupChi2Histograms(), SetupGridRun(), and WriteToFile().

TH3D** NueFCSensitivity::nbgN [private]

Definition at line 99 of file NueFCSensitivity.h.

Referenced by RunFromGrid(), SetupChi2Histograms(), SetupGridRun(), and WriteToFile().

NueSenseConfig* NueFCSensitivity::nsc [private]

Definition at line 92 of file NueFCSensitivity.h.

Referenced by CalculateChi2(), CalculateFitChi2(), EvaluateOmega(), GetPoint(), Initialize(), LoadEventsFromFile(), ProduceTree(), Run(), RunStandardApproach(), and SetupGridRun().

TH3D** NueFCSensitivity::nsigI [private]

Definition at line 98 of file NueFCSensitivity.h.

Referenced by RunFromGrid(), SetupChi2Histograms(), SetupGridRun(), and WriteToFile().

TH3D** NueFCSensitivity::nsigN [private]

Definition at line 97 of file NueFCSensitivity.h.

Referenced by RunFromGrid(), SetupChi2Histograms(), SetupGridRun(), and WriteToFile().

TH1D* NueFCSensitivity::reweight [private]

Definition at line 71 of file NueFCSensitivity.h.

Referenced by Initialize(), and OscillateHist().

TH1D* NueFCSensitivity::signuehist [private]

Definition at line 70 of file NueFCSensitivity.h.

Referenced by Initialize(), and OscillateHist().

float NueFCSensitivity::signuenum [private]

Definition at line 69 of file NueFCSensitivity.h.

Referenced by Initialize(), and OscillateNumber().

float NueFCSensitivity::t13Step [private]

Definition at line 102 of file NueFCSensitivity.h.

Referenced by RunStandardApproach(), and SetupChi2Histograms().


The documentation for this class was generated from the following files:
Generated on Mon Sep 1 00:52:09 2014 for loon by  doxygen 1.4.7