NueSensitivity Class Reference

#include <NueSensitivity.h>

Inheritance diagram for NueSensitivity:
JobCModule

List of all members.

Public Member Functions

 NueSensitivity ()
 ~NueSensitivity ()
void BeginJob ()
JobCResult Ana (const MomNavigator *mom)
void Analysis ()
void EndJob ()
const RegistryDefaultConfig () const
void Config (const Registry &r)
 NueSensitivity ()
 ~NueSensitivity ()
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)

Private Member Functions

void SetPOT ()
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 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

TF1 * nueAppear
TF1 * numuSurvive
TFile * systematicFile
TH2F * systematicHist_allSys
TH2F * systematicHist_oscSys
Double_t systematicHistNorm
Double_t MDCNearToFar
Double_t MDCChallengePOT
Int_t nNuMuFiles
Int_t nNueFiles
Int_t nNuTauFiles
Int_t nNearFiles
Int_t nChallengeNearFiles
Double_t NuMuFilesPOT
Double_t NueFilesPOT
Double_t NuTauFilesPOT
Double_t NearFilesPOT
Double_t ChallengeNearFilesPOT
Int_t nNearUnknownEvents
Int_t nNearNueEvents
Int_t nNearNuMuEvents
Int_t nNearNuTauEvents
Int_t nNearBeamNueEvents
Int_t nNearNCEvents
Int_t nFarUnknownEvents
Int_t nFarNueEvents
Int_t nFarNuMuEvents
Int_t nFarNuTauEvents
Int_t nFarBeamNueEvents
Int_t nFarNCEvents
Int_t nThetaPoints
Double_t * theta13
Int_t nDeltaPoints
Double_t * delta23
Int_t currentRun
float signuenum
TH1D * signuehist
TH1D * reweight
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
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

Detailed Description

Definition at line 19 of file Archive/NueSensitivity.h.


Constructor & Destructor Documentation

NueSensitivity::NueSensitivity (  ) 

Definition at line 41 of file Archive/NueSensitivity.cxx.

References currentRun, delta23, MDCChallengePOT, MDCNearToFar, nChallengeNearFiles, nDeltaPoints, nFarBeamNueEvents, nFarNCEvents, nFarNueEvents, nFarNuMuEvents, nFarNuTauEvents, nFarUnknownEvents, nNearBeamNueEvents, nNearFiles, nNearNCEvents, nNearNueEvents, nNearNuMuEvents, nNearNuTauEvents, nNearUnknownEvents, nNueFiles, nNuMuFiles, nNuTauFiles, nThetaPoints, nueAppear, numuSurvive, systematicFile, systematicHist_allSys, systematicHist_oscSys, systematicHistNorm, and theta13.

00042 {
00043 
00044   nueAppear = NULL;
00045   numuSurvive = NULL;
00046 
00047   systematicHist_oscSys = NULL;
00048   systematicHist_allSys = NULL;
00049   systematicFile = NULL;
00050   systematicHistNorm = 7.4e20; //Trish's normalisation
00051 
00052   MDCChallengePOT = 7.4e20;
00053   MDCNearToFar = 1.1125e-03; //ratio of #near to #far based on MDC
00054 
00055   nNuMuFiles = 0;
00056   nNueFiles = 0;
00057   nNuTauFiles = 0;
00058   nNearFiles = 0;
00059   nChallengeNearFiles = 0;
00060 
00061   nNearUnknownEvents = 0;
00062   nNearNueEvents = 0;
00063   nNearNuMuEvents = 0;
00064   nNearNuTauEvents = 0;
00065   nNearBeamNueEvents = 0;
00066   nNearNCEvents = 0;
00067 
00068   nFarUnknownEvents = 0;
00069   nFarNueEvents = 0;
00070   nFarNuMuEvents = 0;
00071   nFarNuTauEvents = 0;
00072   nFarBeamNueEvents = 0;
00073   nFarNCEvents = 0;
00074 
00075   nDeltaPoints = 60;
00076   nThetaPoints = 40;
00077   theta13 = NULL;
00078   delta23 = NULL;
00079 
00080   currentRun = 0;
00081 
00082 }

NueSensitivity::~NueSensitivity (  ) 

Definition at line 85 of file Archive/NueSensitivity.cxx.

References delta23, nueAppear, numuSurvive, and theta13.

00086 {
00087   delete [] theta13;
00088   delete [] delta23;
00089   delete nueAppear;
00090   delete numuSurvive;
00091 }

NueSensitivity::NueSensitivity (  ) 
NueSensitivity::~NueSensitivity (  ) 

Member Function Documentation

JobCResult NueSensitivity::Ana ( const MomNavigator mom  )  [virtual]

Implement this for read only access to the MomNavigator

Reimplemented from JobCModule.

Definition at line 220 of file Archive/NueSensitivity.cxx.

References ChallengeNearFilesPOT, currentRun, delta23, det, MuELoss::e, HistMan::Fill2d(), MomNavigator::FragmentIter(), VldContext::GetDetector(), SntpHelpers::GetEvent(), SntpHelpers::GetEvent2MCIndex(), RecRecordImp< T >::GetHeader(), SntpHelpers::GetMCTruth(), RecDataHeader::GetRun(), RecHeader::GetVldContext(), NtpSRStripPulseHeight::gev, NtpMCTruth::iaction, NtpMCTruth::inu, NtpMCTruth::inunoosc, isMC, Msg::kDebug, Msg::kError, JobCResult::kFailed, Detector::kFar, Msg::kInfo, Detector::kNear, JobCResult::kPassed, Msg::kVerbose, NtpMCRecord::mc, MDCChallengePOT, MSG, nDeltaPoints, NearFilesPOT, nFarBeamNueEvents, nFarNCEvents, nFarNueEvents, nFarNuMuEvents, nFarNuTauEvents, nFarUnknownEvents, nNearBeamNueEvents, nNearNCEvents, nNearNueEvents, nNearNuMuEvents, nNearNuTauEvents, nNearUnknownEvents, nThetaPoints, nueAppear, NueFilesPOT, NuMuFilesPOT, numuSurvive, NuTauFilesPOT, NtpMCTruth::p4neu, NtpSREvent::ph, Munits::sr, th(), and theta13.

00221 {
00222   //get all NueRecords from mom 
00223   //may have more than one per go since mom reads in a snarl's worth of data
00224   //so, this is a little more complicated than just asking for a NueRecord
00225   MSG("NueSensitivity",Msg::kVerbose) << "**********IN ANA**********" << endl;
00226 
00227   TObject *obj=0;
00228   TIter objiter = mom->FragmentIter();
00229   NtpSRRecord *sr=0;
00230   NtpMCRecord *mc=0;
00231   NtpTHRecord *th=0;
00232   vector<NuePID *> vpid;
00233   while((obj=objiter.Next())){
00234     const char *cn=obj->ClassName();
00235     if(strcmp(cn,"NtpSRRecord")==0){
00236       sr = dynamic_cast<NtpSRRecord *>(obj);
00237     }
00238     else if(strcmp(cn,"NtpMCRecord")==0){
00239       mc = dynamic_cast<NtpMCRecord *>(obj);
00240     }
00241     else if(strcmp(cn,"NtpTHRecord")==0){
00242       th = dynamic_cast<NtpTHRecord *>(obj);
00243     }
00244     else if(strcmp(cn,"NuePID")==0){
00245       NuePID *npid  = dynamic_cast<NuePID *>(obj);
00246       vpid.push_back(npid);
00247     }
00248     else{
00249       continue;
00250     }
00251   }
00252 
00253   bool isMC = true;
00254   bool isTH = true;
00255   if(!sr){
00256     MSG("NueSensitivity",Msg::kError)<<"No NtpSRRecord in Mom"<<endl;
00257     return JobCResult::kFailed;
00258   }  
00259   if(!mc) isMC = false;
00260   if(!th) isTH = false;
00261   if(vpid.size()==0){
00262     MSG("NueSensitivity",Msg::kError)<<"No NuePID Records in Mom"<<endl;
00263     return JobCResult::kFailed;
00264   }
00265 
00266   if(sr->GetHeader().GetRun()!=currentRun) {
00267     currentRun = sr->GetHeader().GetRun();
00268     MSG("NueSensitivity",Msg::kInfo)<< "Current Run: " 
00269                                     << currentRun << endl;
00270   }
00271 
00272   //use HistMan to plot something for events that pass/fail
00273   HistMan man("sensitivity");
00274   Int_t det = 0;
00275   if(sr->GetHeader().GetVldContext().GetDetector()==Detector::kNear) 
00276     {
00277       nueAppear->SetParameter(0,1.0);
00278       numuSurvive->SetParameter(0,1.0);
00279       det = 1;
00280     }
00281   else if(sr->GetHeader().GetVldContext().GetDetector()==Detector::kFar) 
00282     {
00283       nueAppear->SetParameter(0,735.);
00284       numuSurvive->SetParameter(0,735.);
00285       det = 2;
00286     }
00287 
00288   //so, mom will match up snarls for us,
00289   //but, we have to match up events for ourselves.  
00290   Int_t *nmc = NULL; 
00291   if(isMC){ 
00292     TClonesArray& mcArray = *(mc->mc);
00293     nmc = new Int_t[mcArray.GetEntries()];
00294     for(int i=0;i<mcArray.GetEntries();i++) nmc[i] = 0;
00295   }
00296   for(unsigned int i=0;i<vpid.size();i++){
00297     int evtno = vpid[i]->GetHeader().GetEventNo();
00298     int mcindex=0;
00299     NtpSREvent *evt = NULL;
00300     NtpMCTruth *mcth = NULL;
00301     if(evtno<0){
00302       MSG("NueSensitivity",Msg::kDebug)<< "can't get event "
00303                                        << evtno << endl;
00304       if(isMC && isTH) {
00305         mcth = SntpHelpers::GetMCTruth(mcindex,mc);
00306         if(mcth==0) continue;
00307       }
00308     }
00309     else{
00310       evt = SntpHelpers::GetEvent(evtno,sr);
00311       if(isMC && isTH){
00312         mcindex = SntpHelpers::GetEvent2MCIndex(evtno,th);
00313         if(mcindex>=0) nmc[mcindex] += 1;
00314         mcth = SntpHelpers::GetMCTruth(mcindex,mc);
00315         if(mcth==0){
00316           MSG("NueSensitivity",Msg::kError)<< "can't get mctruth for event "
00317                                            << evtno << endl;
00318           continue;
00319         }
00320       }
00321     }
00322 
00323     //nuIntType:
00324     //0=unknown; 1=nue from oscillated numu; 2=numuCC; 
00325     //3=nutauCC; 4=beam nue; 5=nc      
00326     Int_t nuIntType = 0;
00327     if(isMC) {
00328       if(mcth->iaction==0){
00329         nuIntType=5;
00330         //if(nmc[mcindex]==1) {
00331         if(det==1) nNearNCEvents +=1;
00332         else if(det==2) nFarNCEvents +=1;
00333         //}
00334       }
00335       else if(abs(mcth->inu)==12){
00336         if(abs(mcth->inunoosc)==12){
00337           nuIntType=4;
00338           //if(nmc[mcindex]==1) {
00339           if(det==1) nNearBeamNueEvents +=1;
00340           else if(det==2) nFarBeamNueEvents +=1;
00341           //}
00342         }
00343         else if(abs(mcth->inunoosc)==14){
00344           nuIntType=1;
00345           //if(nmc[mcindex]==1) {
00346           if(det==1) nNearNueEvents +=1;
00347           else if(det==2) nFarNueEvents +=1;
00348           //}
00349         }
00350       }
00351       else if(abs(mcth->inu)==14){
00352         nuIntType=2;
00353         //if(nmc[mcindex]==1) {
00354         if(det==1) nNearNuMuEvents +=1;
00355         else if(det==2) nFarNuMuEvents +=1;
00356         //}
00357       }
00358       else if(abs(mcth->inu)==16){
00359         nuIntType=3;
00360         //if(nmc[mcindex]==1) {
00361         if(det==1) nNearNuTauEvents +=1;
00362         else if(det==2) nFarNuTauEvents +=1;
00363         //}
00364       }
00365     }
00366 
00367     if(nuIntType==0) {
00368       if(!isMC) {
00369         if(det==1) nNearUnknownEvents +=1;
00370         else if(det==2) nFarUnknownEvents +=1;
00371       }
00372     }
00373     
00374     int pass = vpid[i]->IsNue;
00375     if(pass==1){
00376       
00377       if(evt==0) {
00378         MSG("NueSensitivity",Msg::kError)<< "Have PID==pass but evt==0!"
00379                                          << endl;
00380         continue;
00381       }
00382 
00383       //POT norms:
00384       Double_t totPOT = 0;
00385       Double_t NCScale = 1;
00386       Double_t BeamNueScale = 1;
00387       Double_t nuMuScale = 1;
00388       Double_t nuEScale = 1;
00389       Double_t nuTauScale = 1;
00390       
00391       if(det==1){
00392         //for NearDet scale everything to MDC Challenge Set POT:
00393         totPOT       = NearFilesPOT;
00394         NCScale      = MDCChallengePOT/totPOT;
00395         BeamNueScale = MDCChallengePOT/totPOT;
00396         nuMuScale    = MDCChallengePOT/totPOT;
00397         nuEScale     = MDCChallengePOT/totPOT;
00398         nuTauScale   = MDCChallengePOT/totPOT;
00399       }
00400       else if(det==2){
00401         totPOT = NuMuFilesPOT + NueFilesPOT + NuTauFilesPOT;
00402         //normalise NC to MDC Challenge POT:
00403         NCScale = MDCChallengePOT/totPOT;
00404         //normalise beam nue to MDC Challenge POT:
00405         //there are no beam nue's in the nutau files!
00406         BeamNueScale = MDCChallengePOT/(NuMuFilesPOT + NueFilesPOT);
00407         //normalise numu to MDC Challenge POT: 
00408         nuMuScale = MDCChallengePOT/NuMuFilesPOT;
00409         //normalise nue to MDC Challenge POT:
00410         nuEScale = MDCChallengePOT/NueFilesPOT;
00411         //normalise nutau to MDC Challenge POT:
00412         nuTauScale = MDCChallengePOT/NuTauFilesPOT;
00413       }
00414 
00415       if(det==1){
00416         if(nuIntType==1){
00417           man.Fill2d("nearHist",nuIntType,evt->ph.gev,
00418                      nuEScale);
00419         }
00420         else if(nuIntType==2){
00421           man.Fill2d("nearHist",nuIntType,evt->ph.gev,
00422                      nuMuScale);
00423         }
00424         else if(nuIntType==3){
00425           man.Fill2d("nearHist",nuIntType,evt->ph.gev,
00426                      nuTauScale);
00427         }
00428         else if(nuIntType==4){
00429           man.Fill2d("nearHist",nuIntType,evt->ph.gev,
00430                      BeamNueScale);
00431         }
00432         else if(nuIntType==5){
00433           man.Fill2d("nearHist",nuIntType,evt->ph.gev,
00434                      NCScale);
00435         }
00436         else if(nuIntType==0){ //assume that these events are Challenge
00437           man.Fill2d("nearHist",nuIntType,evt->ph.gev,
00438                      MDCChallengePOT/ChallengeNearFilesPOT);  
00439         }
00440       }
00441       else if(det==2){
00442         for(int j=0;j<nThetaPoints;j++){
00443           for(int k=0;k<nDeltaPoints;k++){
00444             //histman names:
00445             char name[256];
00446             sprintf(name,"farHist%.0f_%.0f",1000*theta13[j],10000.*delta23[k]);
00447             
00448             Double_t nueAppearWeight = 1;
00449             Double_t numuSurviveWeight = 1;
00450             if(isMC){
00451               //osc probs:
00452               nueAppear->SetParameter(3,theta13[j]);
00453               numuSurvive->SetParameter(3,theta13[j]);
00454               nueAppear->SetParameter(4,delta23[k]);
00455               numuSurvive->SetParameter(4,delta23[k]);
00456               nueAppearWeight = nueAppear->Eval(mcth->p4neu[3]);
00457               numuSurviveWeight = numuSurvive->Eval(mcth->p4neu[3]);
00458             }
00459             double UE32 = 0.5*(1 - sqrt(1-theta13[j]));
00460             if(((4.*UE32*(1-UE32)) - theta13[j])>1e-7) {
00461               cout << "calc not right: UE32 = " << UE32 
00462                    << " theta13 = " << theta13[j] << endl;
00463             }
00464             
00465 //          cout << "nue appearance prob: " << endl;
00466 //          cout << nueAppearWeight - 
00467 //            Oscillate(12,14,mcth->p4neu[3],
00468 //                      735.,delta23[k],0.65,UE32)
00469 //               << endl;
00470 //          cout << "numu dissappearance prob: " << endl;
00471 //          cout << numuSurviveWeight-
00472 //            Oscillate(14,14,mcth->p4neu[3],
00473 //                      735.,delta23[k],0.65,UE32)
00474 //               << endl;
00475 //          cout << "nutau appearance prob: " << endl;
00476 //          cout << 1.-numuSurviveWeight -
00477 //            Oscillate(16,14,mcth->p4neu[3],
00478 //                            735.,delta23[k],0.65,UE32)
00479 //               << endl;
00480 
00481             if(nuIntType==1){
00482               man.Fill2d(name,nuIntType,evt->ph.gev,
00483                          nueAppearWeight*nuEScale);
00484               //man.Fill2d(name,nuIntType,evt->ph.gev,
00485               //         Oscillate(12,14,mcth->p4neu[3],
00486               //                   735.,delta23[k],0.65,
00487               //                   UE32)*nuEScale);
00488             }
00489             else if(nuIntType==2){
00490               man.Fill2d(name,nuIntType,evt->ph.gev,
00491                  numuSurviveWeight*nuMuScale);
00492               // man.Fill2d(name,nuIntType,evt->ph.gev,
00493               // Oscillate(14,14,mcth->p4neu[3],
00494               //           735.,delta23[k],0.65,
00495               //           UE32)*nuMuScale);
00496             }
00497             else if(nuIntType==3){
00498               man.Fill2d(name,nuIntType,evt->ph.gev,
00499                  (1.-numuSurviveWeight)*nuTauScale);
00500               //man.Fill2d(name,nuIntType,evt->ph.gev,
00501               // Oscillate(16,14,mcth->p4neu[3],
00502               //           735.,delta23[k],0.65,
00503               //           UE32)*nuTauScale);
00504             }
00505             else if(nuIntType==4){
00506               man.Fill2d(name,nuIntType,evt->ph.gev,
00507                  BeamNueScale);
00508               //man.Fill2d(name,nuIntType,evt->ph.gev,
00509               // Oscillate(12,12,mcth->p4neu[3],
00510               //           735.,delta23[k],0.65,
00511               //           UE32)*BeamNueScale);
00512             }
00513             else if(nuIntType==5){
00514               man.Fill2d(name,nuIntType,evt->ph.gev,
00515                          NCScale);
00516             }
00517             else if(nuIntType==0){ //assume that these events are Challenge
00518               man.Fill2d(name,nuIntType,evt->ph.gev);
00519             }
00520           }
00521         }
00522       }
00523     }
00524   }
00525   delete [] nmc;
00526   return JobCResult::kPassed; // kNoDecision, kFailed, etc.
00527 }

void NueSensitivity::Analysis (  ) 

Definition at line 530 of file Archive/NueSensitivity.cxx.

References HistMan::Adopt(), delta23, HistMan::Get(), Msg::kInfo, MDCChallengePOT, MDCNearToFar, MSG, nDeltaPoints, NearFilesPOT, nFarBeamNueEvents, nFarNCEvents, nFarNueEvents, nFarNuMuEvents, nFarNuTauEvents, nFarUnknownEvents, nNearBeamNueEvents, nNearNCEvents, nNearNueEvents, nNearNuMuEvents, nNearNuTauEvents, nNearUnknownEvents, nThetaPoints, NueFilesPOT, NuMuFilesPOT, NuTauFilesPOT, systematicHist_allSys, systematicHist_oscSys, systematicHistNorm, and theta13.

Referenced by EndJob().

00531 {
00532   
00533   HistMan man("sensitivity");
00534   
00535   MSG("NueSensitivity",Msg::kInfo) << "==============================" << endl;
00536   MSG("NueSensitivity",Msg::kInfo) 
00537     << "Job Summary: (considering deltam^{2}_{23} = 0.0025 eV^{2})" << endl;
00538   MSG("NueSensitivity",Msg::kInfo) << "------------------------------" << endl;
00539   MSG("NueSensitivity",Msg::kInfo) << "Near Detector: " << endl;
00540   MSG("NueSensitivity",Msg::kInfo) << " Total Events = " 
00541                                    << man.Get<TH2F>("nearHist")->Integral()
00542                                    << endl;
00543   MSG("NueSensitivity",Msg::kInfo) << "------------------------------" << endl;
00544   MSG("NueSensitivity",Msg::kInfo) << "Far Detector: " << endl;
00545   for(int i=0;i<nThetaPoints;i++){
00546     char name[256];
00547     sprintf(name,"farHist%.0f_%.0f",1000.*theta13[i],25.);
00548     MSG("NueSensitivity",Msg::kInfo) << "SinSq(2Theta13) = " << theta13[i] 
00549                                      << " Total Events = " 
00550                                      << man.Get<TH2F>(name)->Integral()
00551                                      << endl;
00552   }
00553 
00554   MSG("NueSensitivity",Msg::kInfo) << "------------------------------" << endl;
00555   MSG("NueSensitivity",Msg::kInfo) << "Event Totals: (Near)  (Far)" << endl;
00556   MSG("NueSensitivity",Msg::kInfo) << "Unknown: " << nNearUnknownEvents << " "
00557                                   << nFarUnknownEvents << endl;
00558   MSG("NueSensitivity",Msg::kInfo) << "Nue: " << nNearNueEvents << " "
00559                                   << nFarNueEvents << endl;
00560   MSG("NueSensitivity",Msg::kInfo) << "NuMu: " << nNearNuMuEvents << " "
00561                                   << nFarNuMuEvents << endl;
00562   MSG("NueSensitivity",Msg::kInfo) << "NuTau: " << nNearNuTauEvents << " "
00563                                   << nFarNuTauEvents << endl;
00564   MSG("NueSensitivity",Msg::kInfo) << "BeamNue: " << nNearBeamNueEvents << " "
00565                                   << nFarBeamNueEvents << endl;
00566   MSG("NueSensitivity",Msg::kInfo) << "NC: " << nNearNCEvents << " "
00567                                   << nFarNCEvents << endl;
00568 
00569 
00570   MSG("NueSensitivity",Msg::kInfo) << "------------------------------" 
00571                                   << endl;
00572   MSG("NueSensitivity",Msg::kInfo) << "POT Summary: (in units of 1e20)" 
00573                                   << endl;
00574   MSG("NueSensitivity",Msg::kInfo) << "Near Total: " 
00575                                   << NearFilesPOT/1e20 << endl;
00576   MSG("NueSensitivity",Msg::kInfo) << "Far Total: " 
00577                                   << (NuMuFilesPOT + 
00578                                       NueFilesPOT + 
00579                                       NuTauFilesPOT)/1e20 
00580                                   << endl;
00581   MSG("NueSensitivity",Msg::kInfo) << "consisting of: " << endl; 
00582   MSG("NueSensitivity",Msg::kInfo) << "numu POT: " 
00583                                   << NuMuFilesPOT/1e20 << endl;
00584   MSG("NueSensitivity",Msg::kInfo) << "nue POT: " 
00585                                   << NueFilesPOT/1e20 << endl;
00586   MSG("NueSensitivity",Msg::kInfo) << "nutau POT: " 
00587                                   << NuTauFilesPOT/1e20 << endl;
00588   MSG("NueSensitivity",Msg::kInfo) << "==============================" 
00589                                   << endl;
00590 
00591 
00592   //Make 2D sensitivity plot:
00593   TH2F *sensitivityHist1 = man.Get<TH2F>("sensitivityHist1");
00594   TH2F *sensitivityHist2 = man.Get<TH2F>("sensitivityHist2");
00595 
00596   for(int i=0;i<nDeltaPoints;i++){
00597     char name[256];
00598     sprintf(name,"farHist%.0f_%.0f",1000.*theta13[0],10000.*delta23[i]);
00599     TH2F *bkgFar = man.Get<TH2F>(name);
00600     TH1D *bkgHist = bkgFar->ProjectionY("bkgHist");
00601     
00602     for(int j=0;j<nThetaPoints;j++){
00603       sprintf(name,"farHist%.0f_%.0f",1000.*theta13[j],10000.*delta23[i]);
00604       TH2F *theFar = man.Get<TH2F>(name);
00605       TH1D *sigHist = theFar->ProjectionY("sigHist");
00606 
00607       Float_t chi2 = 0;
00608       for(int k=1;k<=200;k++){
00609         double sig = sigHist->GetBinContent(k);
00610         double bkg = bkgHist->GetBinContent(k);
00611         if(sig==0||bkg==0) chi2 += 2.*(bkg-sig);
00612         else chi2 += 2.*(bkg-sig)+2.*sig*TMath::Log(sig/bkg);
00613       }
00614 
00615       sensitivityHist1->Fill(theta13[j],delta23[i],chi2);
00616       sensitivityHist2->Fill(theta13[j],delta23[i],
00617                              2.*(bkgFar->Integral()-theFar->Integral()) +
00618                              2.*theFar->Integral() *
00619                              TMath::Log(theFar->Integral() / 
00620                                         bkgFar->Integral()));
00621 
00622       delete sigHist;
00623     }
00624     delete bkgHist;
00625   }
00626 
00627   //for calculating prob given difference in #events:
00628   TF1 *probgaus = new TF1("probgaus","gaus",-10,10);
00629   probgaus->SetParameters(1./TMath::Sqrt(2.*TMath::Pi()),0,1);
00630   
00631   //check that near det is present:
00632   TH2F *nearHist = man.Get<TH2F>("nearHist");
00633   Float_t normNear = nearHist->Integral();  
00634   if(normNear==0) return; //if no near det then don't try to do study
00635   
00636   //for estimating error on FD from ND:
00637   TF1 *gaus = new TF1("gaus","gaus",0,100);
00638 
00639   //calculate simple normalisation for Near to Far:
00640   Float_t normFar = man.Get<TH2F>("farHist0_25")->Integral();
00641   Float_t norm = normFar/normNear;
00642 
00643   //check to see if Near detector files are challenge set
00644   bool isChallenge = false;
00645   TH1D * normNearHistProj = nearHist->ProjectionX("normNearHistProj");
00646   if(normNearHistProj->GetBinContent(1) > 0) isChallenge = true;
00647   delete normNearHistProj;
00648 
00649   //remember some things for different theta_13's:
00650   Double_t *testStat_noSys  = new Double_t[nThetaPoints];
00651   Double_t *testProb_noSys  = new Double_t[nThetaPoints];
00652   Double_t *testStat_oscSys = new Double_t[nThetaPoints];
00653   Double_t *testProb_oscSys = new Double_t[nThetaPoints];
00654   Double_t *testStat_allSys = new Double_t[nThetaPoints];
00655   Double_t *testProb_allSys = new Double_t[nThetaPoints];
00656   for(int i=0;i<nThetaPoints;i++){
00657     testStat_noSys[i]  = 0;    testProb_noSys[i]  = 0;
00658     testStat_oscSys[i] = 0;    testProb_oscSys[i] = 0;
00659     testStat_allSys[i] = 0;    testProb_allSys[i] = 0;
00660   }
00661 
00662   //make TGraphs for a particular delta m^{2} (= 0.0025)
00663   for(int i=0;i<nThetaPoints;i++){
00664     char name[256];
00665     sprintf(name,"farHist%.0f_%.0f",1000.*theta13[i],25.);
00666     
00667     //predict number of far events from near measurement (using simple norm)
00668     Float_t nEventsFromNear = normNear*MDCNearToFar;
00669     
00670     TH2F *farHist = man.Get<TH2F>(name);
00671     Float_t nEvents = farHist->Integral();
00672     //if ND challenge present, use measured event rate to correct far rate
00673     if(isChallenge) nEvents *= MDCNearToFar/norm;
00674     Float_t statError = TMath::Sqrt(nEvents);
00675     
00676     //One-tail hypothesis test no Sys:
00677     testStat_noSys[i] = (nEvents-nEventsFromNear)/statError;
00678     testProb_noSys[i] = probgaus->Integral(testStat_noSys[i],10.);
00679     
00680     if(systematicHist_oscSys){
00681       Int_t bin = 1;
00682       bin = systematicHist_oscSys->GetXaxis()->FindBin(systematicHistNorm * 
00683                                                        nearHist->Integral() / 
00684                                                        MDCChallengePOT);
00685       
00686       TH1D *tempHist = systematicHist_oscSys->ProjectionY("tempHist",bin,bin);
00687       tempHist->Fit("gaus");
00688       Double_t mean=gaus->GetParameter(1)*MDCChallengePOT/systematicHistNorm;;
00689       Double_t sigma=gaus->GetParameter(2)*MDCChallengePOT/systematicHistNorm;;
00690       Double_t rms=tempHist->GetRMS()*MDCChallengePOT/systematicHistNorm;;
00691       //don't trust sigma, if rms is much smaller
00692       //if(sigma>2.*rms)
00693       sigma = rms;
00694       mean = nEventsFromNear;
00695       
00696       //One-tail hypothesis test osc Sys:
00697       testStat_oscSys[i] = (nEvents-mean)/sqrt(statError*statError + 
00698                                                sigma*sigma);
00699       testProb_oscSys[i] = probgaus->Integral(testStat_oscSys[i],10.);
00700       delete tempHist;
00701     }
00702     
00703     if(systematicHist_allSys){
00704       Int_t bin = 1;
00705       bin = systematicHist_allSys->GetXaxis()->FindBin(systematicHistNorm * 
00706                                                        nearHist->Integral() / 
00707                                                        MDCChallengePOT);
00708       
00709       TH1D *tempHist = systematicHist_allSys->ProjectionY("tempHist",bin,bin);
00710       tempHist->Fit("gaus");
00711       Double_t mean=gaus->GetParameter(1)*MDCChallengePOT/systematicHistNorm;
00712       Double_t sigma=gaus->GetParameter(2)*MDCChallengePOT/systematicHistNorm;
00713       Double_t rms=tempHist->GetRMS()*MDCChallengePOT/systematicHistNorm;
00714       //don't trust sigma, if rms is much smaller
00715       if(sigma>2.*rms) sigma = rms;
00716       
00717       //One-tail hypothesis test all Sys:
00718       testStat_allSys[i] = (nEvents-mean)/sqrt(statError*statError + 
00719                                                sigma*sigma);
00720       testProb_allSys[i] = probgaus->Integral(testStat_allSys[i],10.);
00721       delete tempHist;
00722     }
00723 
00724   }
00725 
00726   TGraph *sensitivityGraph1_noSys = new TGraph(40,theta13,testStat_noSys);
00727   sensitivityGraph1_noSys->SetName("sensitivityGraph1_noSys");
00728   sensitivityGraph1_noSys->SetTitle("One-Tailed Hypothesis Test Statistic vs Sin^{2}(2#theta_{13}) - No Systematics");
00729   man.Adopt("",sensitivityGraph1_noSys);
00730 
00731   TGraph *sensitivityGraph2_noSys = new TGraph(40,theta13,testProb_noSys);
00732   sensitivityGraph2_noSys->SetName("sensitivityGraph2_noSys");
00733   sensitivityGraph2_noSys->SetTitle("Probability Data Is Consistent with #theta_{13}=0 vs Sin^{2}(2#theta_{13}) - No Systematics");
00734   man.Adopt("",sensitivityGraph2_noSys);
00735 
00736   TGraph *sensitivityGraph1_oscSys = new TGraph(40,theta13,testStat_oscSys);
00737   sensitivityGraph1_oscSys->SetName("sensitivityGraph1_oscSys");
00738   sensitivityGraph1_oscSys->SetTitle("One-Tailed Hypothesis Test Statistic vs Sin^{2}(2#theta_{13}) - Oscillation Systematics");
00739   man.Adopt("",sensitivityGraph1_oscSys);
00740 
00741   TGraph *sensitivityGraph2_oscSys = new TGraph(40,theta13,testProb_oscSys);
00742   sensitivityGraph2_oscSys->SetName("sensitivityGraph2_oscSys");
00743   sensitivityGraph2_oscSys->SetTitle("Probability Data Is Consistent with #theta_{13}=0 vs Sin^{2}(2#theta_{13}) - Oscillation Systematics");
00744   man.Adopt("",sensitivityGraph2_oscSys);
00745 
00746   TGraph *sensitivityGraph1_allSys = new TGraph(40,theta13,testStat_allSys);
00747   sensitivityGraph1_allSys->SetName("sensitivityGraph1_allSys");
00748   sensitivityGraph1_allSys->SetTitle("One-Tailed Hypothesis Test Statistic vs Sin^{2}(2#theta_{13}) - Oscillation + Neugen Systematics");
00749   man.Adopt("",sensitivityGraph1_allSys);
00750 
00751   TGraph *sensitivityGraph2_allSys = new TGraph(40,theta13,testProb_allSys);
00752   sensitivityGraph2_allSys->SetName("sensitivityGraph2_allSys");
00753   sensitivityGraph2_allSys->SetTitle("Probability Data Is Consistent with #theta_{13}=0 vs Sin^{2}(2#theta_{13}) - Oscillation + Neugen Systematics");
00754   man.Adopt("",sensitivityGraph2_allSys);
00755 
00756   delete [] testStat_noSys;     delete [] testProb_noSys;
00757   delete [] testStat_oscSys;    delete [] testProb_oscSys;
00758   delete [] testStat_allSys;    delete [] testProb_allSys;
00759   
00760 }

void NueSensitivity::BeginJob ( void   )  [virtual]

Implement for notification of begin of job

Reimplemented from JobCModule.

Definition at line 152 of file Archive/NueSensitivity.cxx.

References HistMan::Book(), delta23, MuELoss::e, nDeltaPoints, nThetaPoints, nueAppear, numuSurvive, SetPOT(), systematicFile, systematicHist_allSys, systematicHist_oscSys, and theta13.

00152                              {
00153 
00154   Double_t baseline = 735.0;
00155   nueAppear = new TF1("nueAppear",ElecAppear,0.05,100,10);
00156   nueAppear->SetParameter(0,baseline); //baseline (km)
00157   nueAppear->SetParameter(1,0.6); //sinsq_2th23
00158   nueAppear->SetParameter(2,0.554); //sinsq_2th12
00159   nueAppear->SetParameter(3,0.0); //sinsq_2th13
00160   nueAppear->SetParameter(4,0.002); //dmsq23
00161   nueAppear->SetParameter(5,8.2e-5); //dmsq12
00162   nueAppear->SetParameter(6,0); //density
00163   nueAppear->SetParameter(7,0); //cp phase
00164   nueAppear->SetParameter(8,1); //anti-nu
00165 
00166   numuSurvive = new TF1("numuSurvive",MuSurvive,0.05,100,10);
00167   numuSurvive->SetParameter(0,baseline); //baseline (km)
00168   numuSurvive->SetParameter(1,0.6); //sinsq_2th23
00169   numuSurvive->SetParameter(2,0.554); //sinsq_2th12
00170   numuSurvive->SetParameter(3,0.0); //sinsq_2th13
00171   numuSurvive->SetParameter(4,0.002); //dmsq23
00172   numuSurvive->SetParameter(5,8.2e-5); //dmsq12
00173   numuSurvive->SetParameter(6,0); //density
00174   numuSurvive->SetParameter(7,0); //cp phase
00175   numuSurvive->SetParameter(8,1); //anti-nu
00176 
00177   systematicFile = new TFile("sysHys/moneyplot_newest.root","READ");
00178   if(systematicFile->IsOpen() && !systematicFile->IsZombie()){
00179     systematicHist_oscSys = (TH2F*) systematicFile->Get("hfdvnd");
00180     systematicHist_allSys = (TH2F*) systematicFile->Get("hfdvnd");
00181   }
00182 
00183   theta13 = new Double_t[nThetaPoints];
00184   delta23 = new Double_t[nDeltaPoints];
00185 
00186   for(int i=0;i<nThetaPoints;i++){
00187     theta13[i] = 0.000 + 0.005*Double_t(i);
00188   }
00189 
00190   for(int i=0;i<nDeltaPoints;i++){
00191     delta23[i] = 0.001 + 0.0001*Double_t(i);
00192   }
00193 
00194   HistMan man("sensitivity");
00195   man.Book<TH2F>("nearHist","nearHist",6,0,6,200,0,50);
00196   for(int i=0;i<nThetaPoints;i++){
00197     for(int j=0;j<nDeltaPoints;j++){
00198       char name[256];
00199       sprintf(name,"farHist%.0f_%.0f",1000.*theta13[i],10000.*delta23[j]);
00200       man.Book<TH2F>(name,name,6,0,6,200,0,50);
00201     }
00202   }
00203 
00204   man.Book<TH2F>("sensitivityHist1","sensitivityHist1",
00205                  nThetaPoints,theta13[0]-0.0025,
00206                  theta13[nThetaPoints-1]+0.0025,
00207                  nDeltaPoints,delta23[0]-0.00005,
00208                  delta23[nDeltaPoints-1]+0.00005);
00209   man.Book<TH2F>("sensitivityHist2","sensitivityHist2",
00210                  nThetaPoints,theta13[0]-0.0025,
00211                  theta13[nThetaPoints-1]+0.0025,
00212                  nDeltaPoints,delta23[0]-0.00005,
00213                  delta23[nDeltaPoints-1]+0.00005);
00214 
00215 
00216   SetPOT();
00217 }

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

Definition at line 512 of file Sensitivity/NueSensitivity.cxx.

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

00513 {
00514   NSCErrorParam err = nsc->GetErrorConfig(i);
00515 
00516   float nObs = sig+bg;
00517   float nZero = fZeroValBg[i] + fZeroValSig[i];
00518 
00519   double bgerr = err.bg_systematic*fZeroValBg[i];
00520   double sigerr = err.sig_systematic*fZeroValSig[i];
00521 
00522   float syserrsq = bgerr*bgerr + sigerr*sigerr;
00523  
00524   float errScale = nObs/(nObs + syserrsq);
00525                                                                                 
00526   float chi2 = 2*(nZero - nObs + nObs*TMath::Log(nObs/nZero)) * errScale;
00527   return chi2;
00528 }

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

Definition at line 530 of file Sensitivity/NueSensitivity.cxx.

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

Referenced by RunFromGrid(), and RunStandardApproach().

00531 {
00532   NSCErrorParam err = nsc->GetErrorConfig(i);
00533                                                                                 
00534   float nZero = sig+bg;
00535   float nObs = fObserved;
00536   if(!fObservedIsSet) nObs = fZeroValBg[i] + fZeroValSig[i];
00537 
00538   double bgerr = err.bg_systematic*bg;
00539   double sigerr = err.sig_systematic*sig;
00540 
00541   float syserrsq = bgerr*bgerr + sigerr*sigerr;                                               
00542   float errScale = nZero/(nZero + syserrsq);
00543                                                                                 
00544   float chi2 = 2*(nZero - nObs + nObs*TMath::Log(nObs/nZero)) * errScale;
00545   return chi2;
00546 }

void NueSensitivity::Config ( const Registry r  )  [virtual]

Return the actual configuration. If your module directly pulls its configuration from the fConfig Registry, you don't need to override this. Override if you have local config variables.

Reimplemented from JobCModule.

Definition at line 119 of file Archive/NueSensitivity.cxx.

References Registry::Get(), nChallengeNearFiles, nNearFiles, nNueFiles, nNuMuFiles, and nNuTauFiles.

00120 {
00121 //======================================================================
00122 // Configure the module given the Registry r
00123 //======================================================================
00124   
00125   int tmpi;
00126   if(r.Get("nNuMuFiles",tmpi))           { nNuMuFiles  = tmpi;}
00127   if(r.Get("nNueFiles",tmpi))            { nNueFiles   = tmpi;}
00128   if(r.Get("nNuTauFiles",tmpi))          { nNuTauFiles = tmpi;}
00129   if(r.Get("nNearFiles",tmpi))           { nNearFiles  = tmpi;}
00130   if(r.Get("nChallengeNearFiles",tmpi))  { nChallengeNearFiles  = tmpi;}
00131 }

const Registry & NueSensitivity::DefaultConfig ( void   )  const [virtual]

Get the default configuration registry. This should normally be overridden. One useful idiom is to implement it like:

const Registry& MyModule::DefaultConfig() const { static Registry cfg; // never is destroyed if (cfg.Size()) return cfg; // already filled it // set defaults: cfg.Set("TheAnswer",42); cfg.Set("Units","unknown"); return cfg; }

Reimplemented from JobCModule.

Definition at line 94 of file Archive/NueSensitivity.cxx.

References JobCModule::GetName(), Registry::LockValues(), Registry::Set(), and Registry::UnLockValues().

00095 {
00096 //======================================================================
00097 // Supply the default configuration for the module
00098 //======================================================================
00099   static Registry r; // Default configuration for module
00100  
00101   // Set name of config
00102   std::string name = this->GetName();
00103   name += ".config.default";
00104   r.SetName(name.c_str());
00105  
00106   //Set values in configuration
00107   r.UnLockValues();
00108   r.Set("nNuMuFiles",0);
00109   r.Set("nNueFiles",0);
00110   r.Set("nNuTauFiles",0);
00111   r.Set("nNearFiles",0);
00112   r.Set("nChallengeNearFiles",0);
00113   r.LockValues();
00114  
00115   return r;
00116 }

void NueSensitivity::EndJob (  )  [virtual]

Implement for notification of end of job

Reimplemented from JobCModule.

Definition at line 762 of file Archive/NueSensitivity.cxx.

References Analysis(), systematicFile, and HistMan::WriteOut().

00763 {
00764   Analysis();
00765   HistMan man("sensitivity");
00766   man.WriteOut("SensitivityFile.root");
00767   systematicFile->Close();
00768 }

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

Definition at line 305 of file Sensitivity/NueSensitivity.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().

00306 {  
00307    float inPOT = nsc->GetPOT();
00308    float scale = fMeasuredPOT/inPOT;
00309 
00310    //Initilize the background
00311    if(fMethod == NueSenseConfig::kAllNumbers || fMethod == NueSenseConfig::kNueHist){
00312      currentVal[ClassType::NC] = nsc->GetNumber(ClassType::NC)*scale;
00313      currentVal[ClassType::numu] = nsc->GetNumber(ClassType::numu)*scale;
00314      currentVal[ClassType::bnue] = nsc->GetNumber(ClassType::bnue)*scale;
00315      currentVal[ClassType::nutau] = nsc->GetNumber(ClassType::nutau)*scale;
00316    }
00317 
00318    if(fMethod == NueSenseConfig::kAllNumbers){
00319      float nuenumber = nsc->GetNumber(ClassType::nue);
00320      if(nsc->ShouldDeOsc()){
00321         //Assumed originally Oscillated to:
00322         //  pow(TMath::Sin(theta23),2)*4.*UE32*(1-UE32)*pow(oscterm,2);
00323         float UE32 = nsc->GetOldUe3Square();
00324         float deoscweight = UE32*(1-UE32);
00325         signuenum = nuenumber/deoscweight*scale;
00326      }else{
00327         signuenum = nuenumber*scale;
00328      }
00329    }
00330 
00331    if(fMethod == NueSenseConfig::kNueHist){
00332       TFile inf(nsc->GetNueHistFile().c_str());
00333       TH1D *h = dynamic_cast<TH1D*>(inf.Get(nsc->GetNueHistName().c_str()));
00334       if(h == NULL){
00335         TH1D *h2 = dynamic_cast<TH1D*>(inf.Get(nsc->GetNueHistName().c_str()));
00336         if(h2 == NULL)
00337           std::cout<<"Failure to load signal histogram"<<std::endl;
00338         //This isn't really setup yet, but ideally I would convert it over to the right type
00339       }
00340 
00341       TH1D* hist = dynamic_cast<TH1D*>(h->Clone());
00342       hist->SetDirectory(0);
00343 
00344       if(nsc->ShouldDeOsc()){
00345         hist->Reset();
00346 
00347         float UE32 = nsc->GetOldUe3Square();
00348         float dms23 = nsc->GetOldDeltaMSquare();                                                
00349         for(int i = 0; i < h->GetNbinsX()+1; i++){
00350           float E = h->GetBinCenter(i);
00351           float num = h->GetBinContent(i);
00352 
00353           double invKmToeV = 1.97e-10; //convert 1/km to eV
00354           double Delta13Old = dms23*1e-3*735/(4.*E*1e9*invKmToeV);        
00355           double unweight = UE32*(1-UE32)*4*(1.0/2.0)
00356                              *TMath::Power(TMath::Sin(Delta13Old),2);
00357 
00358           hist->Fill(E, num/unweight*scale);
00359         }
00360      }else{
00361        hist->Scale(scale);
00362      }
00363 
00364      signuehist = dynamic_cast<TH1D*>(hist->Clone());
00365      signuehist->SetDirectory(0);
00366      reweight = dynamic_cast<TH1D*>(signuehist->Clone());
00367      reweight->SetDirectory(0);
00368 
00369      fBinCenter = new double[reweight->GetNbinsX()+1];
00370      for(int i = 0; i < reweight->GetNbinsX()+1; i++){
00371          fBinCenter[i] = reweight->GetBinCenter(i);
00372      }
00373      
00374    }
00375 
00376    if(fMethod == NueSenseConfig::kAnaNueFiles)  
00377       LoadEventsFromFile();
00378 
00379 
00380    fDeltaMS12 = nsc->GetDeltaMS12();
00381    float temp = nsc->GetSinS2Th12();
00382    fTh12 = TMath::ASin(TMath::Sqrt(temp))/2.;
00383    temp = nsc->GetSinS2Th23();
00384    fTh23 = TMath::ASin(TMath::Sqrt(temp))/2.;
00385    
00386    fDensity = nsc->GetDensity();
00387 
00388 }

void NueSensitivity::LoadEventsFromFile (  )  [private]

Definition at line 390 of file Sensitivity/NueSensitivity.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().

00391 {
00392   for(int i = 0; i < nsc->GetNumFiles(); i++){
00393    
00394      TChain *selected = new TChain("ana_nue");
00395      TChain *pots = new TChain("pottree");
00396     
00397      selected->Add(nsc->GetFile(i).c_str());
00398      pots->Add(nsc->GetFile(i).c_str());
00399  
00400      int nonOscFlavor, nuFlavor, type;
00401      float trueE;
00402      double skzpweight;        
00403 
00404      selected->SetMakeClass(1);
00405      selected->SetBranchStatus("*", 0);
00406      selected->SetBranchStatus("mctrue*", 1);
00407      selected->SetBranchStatus("fluxweight.totbeamweight", 1);
00408      selected->SetBranchAddress("mctrue.nonOscNuFlavor",&nonOscFlavor);
00409      selected->SetBranchAddress("mctrue.nuFlavor",&nuFlavor);
00410      selected->SetBranchAddress("mctrue.fNueClass",&type);
00411      selected->SetBranchAddress("mctrue.nuEnergy", &trueE);
00412      selected->SetBranchAddress("fluxweight.totbeamweight",&skzpweight);
00413         
00414      NuePOT *np;                                                             
00415      pots->SetBranchAddress("NuePOT", &np);
00416      double filePOT, scale;
00417      filePOT = 0.0;                                                                     
00418 
00419      for(int j=0; j < pots->GetEntries(); j++) {
00420          pots->GetEntry(j);
00421          filePOT += np->pot;
00422      }
00423      scale = filePOT/fMeasuredPOT;
00424         
00425      for(int j = 0; j < selected->GetEntries(); j++){
00426         selected->GetEntry(j);
00427         mininfo newInfo;
00428         newInfo.energy = trueE;
00429         newInfo.NuFlavor = nuFlavor;
00430         newInfo.NuFlavorBeforeOsc = nonOscFlavor;
00431         newInfo.nuClass = type;
00432         if(skzpweight < 0) skzpweight = 1.0;
00433         newInfo.weight = skzpweight*scale;
00434         mcInfo.push_back(newInfo);
00435      }  //End of this file
00436 
00437      delete selected;  
00438      delete pots; 
00439   }  //Done loading all files and POT Normalization is taken care of
00440 }

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

Definition at line 442 of file Sensitivity/NueSensitivity.cxx.

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

Referenced by RunStandardApproach().

00443 {
00444   if(fMethod == NueSenseConfig::kAllNumbers)
00445      return OscillateNumber(Ue32);
00446 
00447   if(fMethod == NueSenseConfig::kNueHist)
00448      return OscillateHist(dm23, Ue32, delta, sign);
00449    
00450   if(fMethod == NueSenseConfig::kAnaNueFiles)
00451      return OscillateFile(dm23, Ue32, delta, sign);
00452 
00453   return -1;
00454 }

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

Definition at line 486 of file Sensitivity/NueSensitivity.cxx.

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

Referenced by Oscillate().

00487 {
00488   float total[5];
00489 
00490   for(int i = 0; i < 5; i++) total[i] = 0;
00491 
00492   for(unsigned int i = 0; i < mcInfo.size(); i++){
00493     float E = mcInfo[i].energy;
00494     int inu = mcInfo[i].NuFlavor;
00495     int inunoosc = mcInfo[i].NuFlavorBeforeOsc;
00496     int nuClass = mcInfo[i].nuClass;
00497     double weight = mcInfo[i].weight;
00498 
00499     double oscweight = 
00500       OscillateMatter(inu,inunoosc,E);
00501 
00502     total[nuClass] += weight*oscweight;
00503   }
00504 
00505   for(int i = 0; i < 5; i++){
00506     currentVal[i] = total[i];
00507   }
00508 
00509   return currentVal[ClassType::nue];
00510 }

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

Definition at line 467 of file Sensitivity/NueSensitivity.cxx.

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

Referenced by Oscillate().

00468 {
00469   reweight->Reset("ICE");
00470 
00471   for(int i = 0; i < signuehist->GetNbinsX()+1; i++){
00472      float E = fBinCenter[i];
00473      if(E > 10) break;
00474      float num = signuehist->GetBinContent(i);
00475      float weight = 1.0;
00476      if(num > 0)    weight =OscillateMatter(12,14,E); //,dm23,th13,delta,sign);
00477                                                                                 
00478      reweight->Fill(E, num*weight);
00479   }
00480  
00481   currentVal[ClassType::nue] = reweight->GetSum();
00482 //  std::cout<<dm23<<"  "<<Ue32<<"  "<<delta<<"  "<<sign<<"  "<<reweight->GetSum()<<std::endl;
00483   return reweight->GetSum();
00484 }                                                                                

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

Definition at line 612 of file Sensitivity/NueSensitivity.cxx.

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

00615 {
00616  Double_t x[1] = {};
00617   x[0] = Energy;
00618   Double_t dm2_12 = fDeltaMS12*1e-5; //best fit SNO
00619   Double_t dm2_23 = dm2;
00620                                                                                 
00621                                                                                 
00622   Double_t par[9] = {0};
00623   par[0] = fBaseLine;
00624   par[1] = fTh23;
00625   par[2] = fTh12;
00626   par[3] = th13;
00627   par[4] = hierarchy*dm2_23;
00628   par[5] = dm2_12;
00629   par[6] = fDensity; //standard rock density
00630   par[7] = delta;
00631   par[8] = 1;
00632   if(nonOscNuFlavor < 0) par[8] = -1;
00633 
00634   //std::cout<<"About to call "<<dm2<<"  "<<ss13<<"  "<<delta<<std::endl;
00635   fOscCalc.SetOscParam(par);
00636   return fOscCalc.Oscillate(nuFlavor, nonOscNuFlavor, Energy);
00637 }

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

Definition at line 639 of file Sensitivity/NueSensitivity.cxx.

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

Referenced by OscillateFile(), and OscillateHist().

00641 {
00642 
00643   if(nonOscNuFlavor > 0) fOscCalc.SetOscParam(OscPar::kNuAntiNu, 1);
00644   if(nonOscNuFlavor < 0) fOscCalc.SetOscParam(OscPar::kNuAntiNu, -1);
00645 
00646   return fOscCalc.Oscillate(nuFlavor, nonOscNuFlavor, Energy);
00647 }

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

Definition at line 456 of file Sensitivity/NueSensitivity.cxx.

References currentVal, NueConvention::nue, and signuenum.

Referenced by Oscillate().

00457 {
00458 //  double Ue3 = TMath::Sin(t13);
00459 
00460   //Assumed originally Oscillated to:
00461   //  pow(TMath::Sin(theta23),2)*4.*UE32*(1-UE32)*pow(oscterm,2);
00462   currentVal[ClassType::nue] = signuenum*Ue32*(1-Ue32);
00463                                                                                            
00464   return currentVal[ClassType::nue];
00465 }

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

Definition at line 670 of file Sensitivity/NueSensitivity.cxx.

References NSCErrorParam::bg_systematic, 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, NSCErrorParam::scale, NSCErrorParam::sig_systematic, and NSCDataParam::start.

Referenced by WriteToFile().

00671 {
00672   infotree = new TTree("OscPar","OscPar");
00673 
00674   double Th23, Th12, dm2_12;
00675   dm2_12 = nsc->GetDeltaMS12()*1e-5;
00676   Th12 = nsc->GetSinS2Th12();
00677   Th23 = nsc->GetSinS2Th23();
00678 
00679   infotree->Branch("L",&fBaseLine, "L/F");
00680   infotree->Branch("Sin2_2Theta23",&Th23, "Sin2_2Theta23/D");
00681   infotree->Branch("Sin2_2Theta12",&Th12, "Sin2_2Theta12/D");
00682   infotree->Branch("DeltaMS12",&dm2_12, "DeltaMS12/D");
00683   infotree->Branch("Density",&fDensity, "Density/D");
00684 
00685   NSCDataParam DPst13 = nsc->GetSinS2Th13();
00686   NSCDataParam DPdelta = nsc->GetDelta();
00687   NSCDataParam DPdm23 = nsc->GetDeltaMS23();
00688 
00689   infotree->Branch("DeltaMS23_start",&DPst13.start, "DeltaMS23_start/F");
00690   infotree->Branch("DeltaMS23_end",&DPst13.end, "DeltaMS23_end/F");
00691   infotree->Branch("Delta_start",&DPdelta.start, "Delta_start/F");
00692   infotree->Branch("Delta_end",&DPdelta.end, "Delta_end/F");
00693   infotree->Branch("Sin2_2Theta13_start",&DPdm23.start, "Sin2_2Theta13_start/F");
00694   infotree->Branch("Sin2_2Theta13_end",&DPdm23.end, "Sin2_2Theta13_end/F");
00695 
00696   infotree->Fill();
00697 
00698   errtree = new TTree("ErrPar","ErrPar");
00699   
00700   double bgsys, sigsys, nc, numu, nue, nutau, bnue;
00701   errtree->Branch("BG_systematic",&bgsys,"BG_systematic/D");
00702   errtree->Branch("SIG_systematic",&sigsys, "SIG_systematic/D");
00703   errtree->Branch("nc_scale",&nc, "nc_scale/D");
00704   errtree->Branch("numu_scale",&numu, "numu_scale/D");
00705   errtree->Branch("nue_scale",&nue, "nue_scale/D");
00706   errtree->Branch("nutau_scale",&nutau, "nutau_scale/D");
00707   errtree->Branch("bnue_scale",&bnue, "bnue_scale/D");
00708 
00709   for(int i = 0; i < nsc->GetNumberConfig(); i++){
00710       NSCErrorParam err = nsc->GetErrorConfig(i);
00711       bgsys = err.bg_systematic;
00712       sigsys = err.sig_systematic;
00713       nc = err.scale[ClassType::NC];
00714       numu = err.scale[ClassType::numu];
00715       bnue =err.scale[ClassType::bnue];
00716       nutau = err.scale[ClassType::nutau];
00717       nue = err.scale[ClassType::nue];
00718       
00719       errtree->Fill();
00720   }
00721 }

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

Definition at line 18 of file Sensitivity/NueSensitivity.cxx.

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

00019 {
00020    //Loading in the data
00021    nsc = new NueSenseConfig(input);
00022    if(!nsc->CheckConfig()) return;
00023    std::cout<<"Configuration file read and confirmed"<<std::endl;
00024 
00025    fMeasuredPOT = outPOT;
00026    // if just numbers roll forward, if taking form a chain pull out
00027    //  just enough information to do the oscillations, if hist grab it
00028 
00029    fMethod = nsc->GetDataMethod();
00030    fNumConfig = nsc->GetNumberConfig();
00031 
00032    Initialize();    //Load in any values and prepare listings
00033    std::cout<<"Initialization complete"<<std::endl;
00034 
00035    if(fMethod != 4)  RunStandardApproach();
00036    else              RunFromGrid();
00037  //Now we have all the data points
00038    // now we write it out to file
00039    WriteToFile(output);
00040 }

void NueSensitivity::RunFromGrid (  ) 

Definition at line 114 of file Sensitivity/NueSensitivity.cxx.

References CalculateFitChi2(), chi2i, chi2n, fDeltaM32MapInvert, fDeltaM32MapNormal, fNumConfig, and SetupGridRun().

Referenced by Run().

00115 { 
00116    SetupGridRun();
00117    std::cout<<"Setup is Complete"<<std::endl;
00118 
00119    double dmDV = fDeltaM32MapNormal.begin()->first;
00120 
00121    std::map<double, std::vector<Point> >::iterator delBeg = fDeltaM32MapNormal[dmDV].begin();
00122    std::map<double, std::vector<Point> >::iterator delEnd = fDeltaM32MapNormal[dmDV].end();
00123 
00124    while(delBeg != delEnd){
00125      double delta = delBeg->first;
00126      for(unsigned int i = 0; i < delBeg->second.size(); i++){
00127         // at this point
00128         double ss2th = delBeg->second[i].th13;
00129 
00130         double signal = delBeg->second[i].nsignal;
00131         double bg = delBeg->second[i].nbg;
00132 
00133         for(int j = 0; j < fNumConfig; j++){
00134           double val = CalculateFitChi2(j, bg, signal);
00135           chi2n[j]->Fill(ss2th, delta, dmDV, val);
00136         }
00137     }
00138     delBeg++;
00139   }
00140   dmDV = fDeltaM32MapInvert.begin()->first;
00141   delBeg = fDeltaM32MapInvert[dmDV].begin();
00142   delEnd = fDeltaM32MapInvert[dmDV].end();
00143   
00144   while(delBeg != delEnd){
00145     double delta = delBeg->first;
00146     for(unsigned int i = 0; i < delBeg->second.size(); i++){
00147        // at this point
00148        double ss2th = delBeg->second[i].th13;
00149        double signal = delBeg->second[i].nsignal;
00150        double bg = delBeg->second[i].nbg;
00151 
00152        for(int j = 0; j < fNumConfig; j++){
00153           double val = CalculateFitChi2(j, bg, signal);
00154           chi2i[j]->Fill(ss2th, delta, dmDV, val);
00155        }
00156     }
00157     delBeg++;
00158   }
00159 }

void NueSensitivity::RunStandardApproach (  ) 

Definition at line 42 of file Sensitivity/NueSensitivity.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().

00043 {
00044    NSCDataParam DPst13 = nsc->GetSinS2Th13();
00045    NSCDataParam DPdelta = nsc->GetDelta();
00046    NSCDataParam DPdm23 = nsc->GetDeltaMS23();
00047 
00048    fZeroValBg = new float[fNumConfig];
00049    fZeroValSig   = new float[fNumConfig];
00050 
00051    //int total = 
00052    SetupChi2Histograms(DPst13, DPdelta, DPdm23);
00053 
00054    SetOscParamBase(0.0024, 0.0, 0, 1);
00055 
00056    //Determine the number of background at Ue3 = 0
00057    Oscillate(0.0024, 0.0, 0, 1);
00058    for(int i = 0; i < fNumConfig; i++){
00059       fZeroValBg[i] = fZeroValSig[i] = 0.0;
00060       GetPoint(i, fZeroValBg[i], fZeroValSig[i]);
00061       std::cout<<"For configuration "<<i<<" starting with "
00062           <<fZeroValBg[i]<<", "<<fZeroValSig[i]<<" events."<<std::endl;
00063    } 
00064 
00065    //And now the big run over all the numbers
00066    bool d2,d3;  //some bools to keep the loops right
00067    float bg, sig;
00068    double pi = 3.1415926;
00069    int count = 0;
00070 
00071    for(double sins2t13 = DPst13.start; sins2t13 <= DPst13.end; sins2t13 += t13Step)
00072    {
00073       double Th13 = TMath::ASin(TMath::Sqrt(sins2t13))/2.;
00074       fOscCalc.SetOscParam(OscPar::kTh13, Th13);
00075 
00076       d2 = DPdelta.isfixed;
00077       for(double delta = DPdelta.start; delta <= DPdelta.end || d2 ; delta += dStep)
00078       {
00079          d3 = DPdm23.isfixed;
00080 
00081          fOscCalc.SetOscParam(OscPar::kDelta, delta*pi);
00082          for(double dm23 = DPdm23.start; dm23 <= DPdm23.end || d3; dm23+= dm23Step)
00083          {
00084             int sign = 1;
00085             fOscCalc.SetOscParam(OscPar::kDeltaM23, dm23*1e-3);
00086 
00087             Oscillate(dm23*1e-3, Th13, delta*pi, sign);   //Change the values
00088             for(int i = 0; i < fNumConfig; i++){
00089               GetPoint(i, bg, sig);
00090               float chi2 = CalculateFitChi2(i, bg, sig);
00091               chi2n[i]->Fill(sins2t13, delta, dm23, chi2);
00092             }
00093 
00094             if(fMethod != NueSenseConfig::kAllNumbers){
00095               sign = -1;
00096               fOscCalc.SetOscParam(OscPar::kDeltaM23, sign*dm23*1e-3);
00097               Oscillate(dm23*1e-3, Th13, delta*pi, sign);   //Change the values
00098 
00099               for(int i = 0; i < fNumConfig; i++){
00100                 GetPoint(i, bg, sig);
00101                 float chi2 = CalculateFitChi2(i, bg, sig);
00102                 chi2i[i]->Fill(sins2t13, delta, dm23, chi2);
00103               }
00104             }
00105             if(DPdm23.isfixed){ d3 = false;  dm23 = DPdm23.end + 1; }
00106             count++;
00107             if(count%100000 == 0) std::cout<<"On iteration "<<count<<std::endl;
00108          }//end of dm23
00109          if(DPdelta.isfixed){ d2 = false;  delta = DPdelta.end + 1; }
00110       }//end of delta
00111    }//end of sins2theta13
00112 }

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

Definition at line 43 of file Sensitivity/NueSensitivity.h.

References fObserved, and fObservedIsSet.

00043 {fObserved = in; fObservedIsSet = true;}

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

Definition at line 649 of file Sensitivity/NueSensitivity.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().

00650                                                                {
00651 
00652   Double_t dm2_12 = fDeltaMS12*1e-5; //best fit SNO
00653   Double_t dm2_23 = dm2;
00654                                                                                 
00655   Double_t par[9] = {0};
00656   par[OscPar::kL] = fBaseLine;
00657   par[OscPar::kTh23] = fTh23;
00658   par[OscPar::kTh12] = fTh12;
00659   par[OscPar::kTh13] = ss13; // TMath::ASin(TMath::Sqrt(ss2th13))/2.;
00660   par[OscPar::kDeltaM23] = hierarchy*dm2_23;
00661   par[OscPar::kDeltaM12] = dm2_12;
00662   par[OscPar::kDensity] = fDensity; //standard rock density
00663   par[OscPar::kDelta] = delta;
00664   par[OscPar::kNuAntiNu] = 1;
00665                                                                                 
00666 //  std::cout<<"About to call "<<dm2<<"  "<<ss13<<"  "<<delta<<std::endl;
00667   fOscCalc.SetOscParam(par);
00668 }

void NueSensitivity::SetPOT (  )  [private]

Definition at line 135 of file Archive/NueSensitivity.cxx.

References ChallengeNearFilesPOT, Msg::kInfo, MSG, nChallengeNearFiles, NearFilesPOT, nNearFiles, nNueFiles, nNuMuFiles, nNuTauFiles, NueFilesPOT, NuMuFilesPOT, NuTauFilesPOT, POTPERFARFILE, and POTPERNEARFILE.

Referenced by BeginJob().

00135                            {
00136 
00137   NuMuFilesPOT          = float(nNuMuFiles)*POTPERFARFILE;
00138   NueFilesPOT           = float(nNueFiles)*POTPERFARFILE;
00139   NuTauFilesPOT         = float(nNuTauFiles)*POTPERFARFILE;
00140   NearFilesPOT          = float(nNearFiles)*POTPERNEARFILE;  
00141   ChallengeNearFilesPOT = float(nChallengeNearFiles)*POTPERNEARFILE;
00142 
00143   MSG("NueSensitivity",Msg::kInfo) << "POT: " << NuMuFilesPOT << " " 
00144                                    << NueFilesPOT << " "
00145                                    << NuTauFilesPOT << " " 
00146                                    << NearFilesPOT << " " 
00147                                    << ChallengeNearFilesPOT << endl;
00148 
00149 }

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

Definition at line 549 of file Sensitivity/NueSensitivity.cxx.

References chi2i, chi2n, dm23Step, dStep, NSCDataParam::end, fNumConfig, NSCDataParam::isfixed, NSCDataParam::start, and t13Step.

Referenced by RunStandardApproach().

00550 {
00551    chi2n = new TH3D*[fNumConfig];
00552    chi2i = new TH3D*[fNumConfig];
00553                                                                                            
00554    int t13Bins = int((DPst13.end - DPst13.start)/t13Step) + 1;
00555    int delBins = int((DPdelta.end - DPdelta.start)/dStep) + 1;
00556    int dm23Bins = int((DPdm23.end - DPdm23.start)/dm23Step) + 1;
00557    
00558    double xstart, xend, ystart, yend, zstart, zend;
00559    xstart = xend = ystart = yend = zstart = zend = 0;
00560                                                                                         
00561    if(DPst13.isfixed){ // i have no idea why you are running the code
00562    }else{
00563       xstart = DPst13.start - t13Step/2; xend = xstart + t13Bins*t13Step;
00564    }
00565    if(DPdelta.isfixed){
00566       ystart = DPdelta.start - 3*dStep/2;
00567       delBins = 5;
00568    }else{
00569       ystart = DPdelta.start - dStep/2;
00570    }
00571    yend = ystart + delBins*dStep;
00572 
00573    if(DPdm23.isfixed){
00574       zstart = DPdm23.start - 3*dm23Step/2;
00575       dm23Bins = 5;
00576    }else{
00577       zstart = DPdm23.start - dm23Step/2;
00578    }
00579    zend = zstart + dm23Bins*dm23Step;
00580 
00581    for(int i = 0; i < fNumConfig; i++){
00582       char temp[20];
00583       sprintf(temp, "chi2normal%d", i);
00584       chi2n[i] = new TH3D(temp, temp, t13Bins, xstart, xend, delBins,
00585                              ystart, yend, dm23Bins, zstart, zend);
00586       sprintf(temp, "chi2inverted%d", i);
00587       chi2i[i] = new TH3D(temp, temp, t13Bins, xstart, xend, delBins,
00588                              ystart, yend, dm23Bins, zstart, zend);
00589    }
00590    std::cout<<"Chi2 Histograms created ("<<t13Bins<<", "<<delBins<<", "<<dm23Bins<<"): " 
00591             <<t13Bins*delBins*dm23Bins<<" iterations to perform"<<std::endl;
00592    return t13Bins*delBins*dm23Bins;
00593 }

void NueSensitivity::SetupGridRun (  )  [private]

Definition at line 161 of file Sensitivity/NueSensitivity.cxx.

References chi2i, chi2n, fDeltaM32MapInvert, fDeltaM32MapNormal, fMeasuredPOT, fNumConfig, NueSenseConfig::GetDataFiles(), NueSenseConfig::GetPOT(), Point::nbg, nsc, Point::nsignal, and Point::th13.

Referenced by RunFromGrid().

00162 { 
00163   std::vector<std::string> datafiles = nsc->GetDataFiles();
00164   float inPOT = nsc->GetPOT();
00165   float scale = fMeasuredPOT/inPOT;
00166  
00167   double dmVal = 0.0; 
00168   for(unsigned int i = 0; i < datafiles.size(); i++){
00169      double th13, delta, mass, ni, bg[5], sig, dummy;
00170      std::ifstream ins(datafiles[i].c_str());
00171 
00172      ins>>th13>>delta>>mass>>ni>>bg[0]>>bg[1]>>bg[2]>>bg[3]>>sig>>dummy;
00173      if(i == 0) dmVal = mass;
00174      while(!ins.eof()){
00175        Point temp;
00176        double temp2 = TMath::Sin(2*th13);
00177 
00178        temp.th13 = temp2*temp2;
00179        temp.nsignal = sig*scale;
00180        temp.nbg = (dummy-sig)*scale;
00181 
00182        if(ni > 0)  fDeltaM32MapNormal[mass][delta].push_back(temp);
00183        else        fDeltaM32MapInvert[mass][delta].push_back(temp);
00184 
00185        ins>>th13>>delta>>mass>>ni>>bg[0]>>bg[1]>>bg[2]>>bg[3]>>sig>>dummy;
00186      }
00187      ins.clear();
00188    }
00189 
00190    std::map<double, std::vector<Point> >::iterator delBeg = fDeltaM32MapNormal[dmVal].begin();
00191    std::map<double, std::vector<Point> >::iterator delEnd = fDeltaM32MapNormal[dmVal].end();
00192 
00193    std::vector<double> deltaM2Pos;
00194    std::map<double, std::map<double, std::vector<Point> > >::iterator dmBeg = fDeltaM32MapNormal.begin();
00195    std::map<double, std::map<double, std::vector<Point> > >::iterator dmEnd = fDeltaM32MapNormal.end();
00196 
00197    while(dmBeg != dmEnd){
00198      deltaM2Pos.push_back(dmBeg->first);
00199      dmBeg++;
00200    }
00201 
00202    std::vector<double> deltaPos;
00203    while(delBeg != delEnd){
00204      deltaPos.push_back(delBeg->first);
00205      delBeg++;
00206    }
00207 
00208    std::vector<double> sin22th13Pos;
00209    std::vector<Point>::iterator thBeg = fDeltaM32MapNormal[dmVal][0].begin();
00210    std::vector<Point>::iterator thEnd = fDeltaM32MapNormal[dmVal][0].end();
00211 
00212    while(thBeg != thEnd){
00213      sin22th13Pos.push_back(thBeg->th13);
00214      thBeg++;
00215    }
00216 
00217    sort(deltaM2Pos.begin(), deltaM2Pos.end());
00218    sort(deltaPos.begin(), deltaPos.end());
00219    sort(sin22th13Pos.begin(), sin22th13Pos.end());
00220 
00221    int nDM2 = deltaM2Pos.size();
00222    int nDelta = deltaPos.size();
00223    int nTh13 = sin22th13Pos.size();
00224 
00225    if(nDelta <= 1){
00226       std::cout<<"Only zero/one delta value??"<<std::endl;
00227       return;
00228    }
00229    if(nTh13 <= 1){
00230       std::cout<<"Only zero/one theta value??"<<std::endl;
00231       return;
00232    }
00233 
00234    double* deltaM2Array = new double[nDM2+1];
00235    double* deltaArray = new double[nDelta+1];
00236    double* sinthArray = new double[nTh13+1];
00237 
00238    double offset = 0.0;
00239    for(int i = 1; i < nDelta; i++){
00240      double prev = deltaPos[i-1];
00241      double curr = deltaPos[i];
00242      offset = (curr - prev)/2.0;
00243      if(i == 1) deltaArray[0] = prev - offset;
00244      deltaArray[i] = curr - offset;
00245    }
00246    deltaArray[nDelta] = deltaArray[nDelta-1] + 2*offset;
00247 
00248    offset = 0.0;
00249    for(int i = 1; i < nTh13; i++){
00250      double prev = sin22th13Pos[i-1];
00251      double curr = sin22th13Pos[i];
00252      offset = (curr - prev)/2.0;
00253      if(i == 1) sinthArray[0] = prev - offset;
00254      sinthArray[i] = curr - offset;
00255    }
00256    sinthArray[nTh13] = sinthArray[nTh13-1] + 2*offset;
00257 
00258    offset = 0.00002;
00259    deltaM2Array[0] = deltaM2Pos[0];
00260    for(int i = 1; i < nDM2; i++){
00261      double prev = deltaM2Pos[i-1];
00262      double curr = deltaM2Pos[i];
00263      offset = (curr - prev)/2.0;
00264      if(i == 1) deltaM2Array[0] = prev - offset;
00265      deltaM2Array[i] = curr - offset;
00266    }
00267    deltaM2Array[nDM2] = deltaM2Array[0] + 2*offset;
00268 
00269    TString n1 = "Chi2HistN";
00270    TString n2 = "Chi2HistI";
00271 
00272    chi2n = new TH3D*[fNumConfig];
00273    chi2i = new TH3D*[fNumConfig];
00274 
00275    for(int i = 0; i < fNumConfig; i++){
00276      char temp[20];
00277      sprintf(temp, "chi2normal%d", i);
00278      chi2n[i] = new TH3D(temp, temp, nTh13, sinthArray, nDelta, deltaArray, nDM2, deltaM2Array);
00279      sprintf(temp, "chi2inverted%d", i);
00280      chi2i[i] = new TH3D(temp, temp, nTh13, sinthArray, nDelta, deltaArray, nDM2, deltaM2Array);
00281 
00282      chi2n[i]->SetDirectory(0);
00283      chi2i[i]->SetDirectory(0);
00284    }
00285    std::cout<<"Chi2 Histograms created ("<<nTh13<<", "<<nDelta<<", "<<nDM2<<"): Running Grid Method"<<std::endl;
00286 
00287    //So now I have built all the histograms as appropriate 
00288    // - now I just loop over the data
00289 
00290 }

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

Definition at line 595 of file Sensitivity/NueSensitivity.cxx.

References chi2i, chi2n, err(), fNumConfig, info, and ProduceTree().

Referenced by Run().

00596 {
00597    TFile out(file.c_str(), "RECREATE");
00598    out.cd();
00599 
00600    for(int i = 0; i < fNumConfig; i++){
00601       chi2n[i]->Write();
00602       chi2i[i]->Write();
00603    }
00604    TTree *info, *err;
00605 
00606    ProduceTree(info, err);
00607    info->Write();
00608    err->Write();
00609    out.Close();
00610 }


Member Data Documentation

Definition at line 63 of file Archive/NueSensitivity.h.

Referenced by Ana(), and SetPOT().

TH3D** NueSensitivity::chi2i [private]
TH3D** NueSensitivity::chi2n [private]
Int_t NueSensitivity::currentRun [private]

Definition at line 85 of file Archive/NueSensitivity.h.

Referenced by Ana(), and NueSensitivity().

float NueSensitivity::currentVal[5] [private]
Double_t* NueSensitivity::delta23 [private]

Definition at line 83 of file Archive/NueSensitivity.h.

Referenced by Ana(), Analysis(), BeginJob(), NueSensitivity(), and ~NueSensitivity().

float NueSensitivity::dm23Step [private]

Definition at line 97 of file Sensitivity/NueSensitivity.h.

Referenced by RunStandardApproach(), and SetupChi2Histograms().

float NueSensitivity::dStep [private]

Definition at line 96 of file Sensitivity/NueSensitivity.h.

Referenced by RunStandardApproach(), and SetupChi2Histograms().

float NueSensitivity::fBaseLine [private]

Definition at line 86 of file Sensitivity/NueSensitivity.h.

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

double* NueSensitivity::fBinCenter [private]

Definition at line 72 of file Sensitivity/NueSensitivity.h.

Referenced by Initialize(), and OscillateHist().

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

Definition at line 106 of file Sensitivity/NueSensitivity.h.

Referenced by RunFromGrid(), and SetupGridRun().

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

Definition at line 105 of file Sensitivity/NueSensitivity.h.

Referenced by RunFromGrid(), and SetupGridRun().

double NueSensitivity::fDeltaMS12 [private]

Definition at line 100 of file Sensitivity/NueSensitivity.h.

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

double NueSensitivity::fDensity [private]
double NueSensitivity::fMeasuredPOT [private]

Definition at line 75 of file Sensitivity/NueSensitivity.h.

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

int NueSensitivity::fMethod [private]

Definition at line 88 of file Sensitivity/NueSensitivity.h.

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

Definition at line 110 of file Sensitivity/NueSensitivity.h.

Referenced by CalculateFitChi2(), and SetObserved().

Definition at line 111 of file Sensitivity/NueSensitivity.h.

Referenced by CalculateFitChi2(), and SetObserved().

double NueSensitivity::fTh12 [private]

Definition at line 101 of file Sensitivity/NueSensitivity.h.

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

double NueSensitivity::fTh23 [private]

Definition at line 102 of file Sensitivity/NueSensitivity.h.

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

float* NueSensitivity::fZeroValBg [private]
float* NueSensitivity::fZeroValSig [private]
std::vector<mininfo> NueSensitivity::mcInfo [private]

Definition at line 74 of file Sensitivity/NueSensitivity.h.

Referenced by LoadEventsFromFile(), and OscillateFile().

Double_t NueSensitivity::MDCChallengePOT [private]

Definition at line 51 of file Archive/NueSensitivity.h.

Referenced by Ana(), Analysis(), and NueSensitivity().

Double_t NueSensitivity::MDCNearToFar [private]

Definition at line 48 of file Archive/NueSensitivity.h.

Referenced by Analysis(), and NueSensitivity().

Definition at line 57 of file Archive/NueSensitivity.h.

Referenced by Config(), NueSensitivity(), and SetPOT().

Definition at line 82 of file Archive/NueSensitivity.h.

Referenced by Ana(), Analysis(), BeginJob(), and NueSensitivity().

Double_t NueSensitivity::NearFilesPOT [private]

Definition at line 62 of file Archive/NueSensitivity.h.

Referenced by Ana(), Analysis(), and SetPOT().

Definition at line 76 of file Archive/NueSensitivity.h.

Referenced by Ana(), Analysis(), and NueSensitivity().

Definition at line 77 of file Archive/NueSensitivity.h.

Referenced by Ana(), Analysis(), and NueSensitivity().

Definition at line 73 of file Archive/NueSensitivity.h.

Referenced by Ana(), Analysis(), and NueSensitivity().

Definition at line 74 of file Archive/NueSensitivity.h.

Referenced by Ana(), Analysis(), and NueSensitivity().

Definition at line 75 of file Archive/NueSensitivity.h.

Referenced by Ana(), Analysis(), and NueSensitivity().

Definition at line 72 of file Archive/NueSensitivity.h.

Referenced by Ana(), Analysis(), and NueSensitivity().

Definition at line 69 of file Archive/NueSensitivity.h.

Referenced by Ana(), Analysis(), and NueSensitivity().

Int_t NueSensitivity::nNearFiles [private]

Definition at line 56 of file Archive/NueSensitivity.h.

Referenced by Config(), NueSensitivity(), and SetPOT().

Definition at line 70 of file Archive/NueSensitivity.h.

Referenced by Ana(), Analysis(), and NueSensitivity().

Definition at line 66 of file Archive/NueSensitivity.h.

Referenced by Ana(), Analysis(), and NueSensitivity().

Definition at line 67 of file Archive/NueSensitivity.h.

Referenced by Ana(), Analysis(), and NueSensitivity().

Definition at line 68 of file Archive/NueSensitivity.h.

Referenced by Ana(), Analysis(), and NueSensitivity().

Definition at line 65 of file Archive/NueSensitivity.h.

Referenced by Ana(), Analysis(), and NueSensitivity().

Int_t NueSensitivity::nNueFiles [private]

Definition at line 54 of file Archive/NueSensitivity.h.

Referenced by Config(), NueSensitivity(), and SetPOT().

Int_t NueSensitivity::nNuMuFiles [private]

Definition at line 53 of file Archive/NueSensitivity.h.

Referenced by Config(), NueSensitivity(), and SetPOT().

Int_t NueSensitivity::nNuTauFiles [private]

Definition at line 55 of file Archive/NueSensitivity.h.

Referenced by Config(), NueSensitivity(), and SetPOT().

Definition at line 79 of file Archive/NueSensitivity.h.

Referenced by Ana(), Analysis(), BeginJob(), and NueSensitivity().

TF1* NueSensitivity::nueAppear [private]

Definition at line 40 of file Archive/NueSensitivity.h.

Referenced by Ana(), BeginJob(), NueSensitivity(), and ~NueSensitivity().

Double_t NueSensitivity::NueFilesPOT [private]

Definition at line 60 of file Archive/NueSensitivity.h.

Referenced by Ana(), Analysis(), and SetPOT().

Double_t NueSensitivity::NuMuFilesPOT [private]

Definition at line 59 of file Archive/NueSensitivity.h.

Referenced by Ana(), Analysis(), and SetPOT().

Definition at line 41 of file Archive/NueSensitivity.h.

Referenced by Ana(), BeginJob(), NueSensitivity(), and ~NueSensitivity().

Double_t NueSensitivity::NuTauFilesPOT [private]

Definition at line 61 of file Archive/NueSensitivity.h.

Referenced by Ana(), Analysis(), and SetPOT().

TH1D* NueSensitivity::reweight [private]

Definition at line 70 of file Sensitivity/NueSensitivity.h.

Referenced by Initialize(), and OscillateHist().

TH1D* NueSensitivity::signuehist [private]

Definition at line 69 of file Sensitivity/NueSensitivity.h.

Referenced by Initialize(), and OscillateHist().

float NueSensitivity::signuenum [private]

Definition at line 68 of file Sensitivity/NueSensitivity.h.

Referenced by Initialize(), and OscillateNumber().

Definition at line 43 of file Archive/NueSensitivity.h.

Referenced by BeginJob(), EndJob(), and NueSensitivity().

Definition at line 44 of file Archive/NueSensitivity.h.

Referenced by Analysis(), BeginJob(), and NueSensitivity().

Definition at line 45 of file Archive/NueSensitivity.h.

Referenced by Analysis(), BeginJob(), and NueSensitivity().

Definition at line 46 of file Archive/NueSensitivity.h.

Referenced by Analysis(), and NueSensitivity().

float NueSensitivity::t13Step [private]

Definition at line 95 of file Sensitivity/NueSensitivity.h.

Referenced by RunStandardApproach(), and SetupChi2Histograms().

Double_t* NueSensitivity::theta13 [private]

Definition at line 80 of file Archive/NueSensitivity.h.

Referenced by Ana(), Analysis(), BeginJob(), NueSensitivity(), and ~NueSensitivity().


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

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1