loon: NueSensitivity.cxx Source File

00001 
00002 // $Id: NueSensitivity.cxx,v 1.1 2006/09/18 21:41:22 boehm Exp $
00003 //
00004 // FILL_IN: [Document your code!!]
00005 //
00006 // vahle@hep.ucl.ac.uk
00008 #include "TH1F.h"
00009 #include "TFile.h"
00010 #include "TClonesArray.h"
00011 #include "TF1.h"
00012 #include "TH2F.h"
00013 #include "TGraph.h"
00014 #include "NueAna/NueSensitivity.h"
00015 #include "NueAna/NuePID.h"
00016 #include "MessageService/MsgService.h"
00017 #include "MinosObjectMap/MomNavigator.h"
00018 #include "JobControl/JobCModuleRegistry.h" // For JOBMODULE macro
00019 #include "HistMan/HistMan.h"
00020 #include "MCNtuple/NtpMCRecord.h"
00021 #include "MCNtuple/NtpMCTruth.h"
00022 #include "TruthHelperNtuple/NtpTHRecord.h"
00023 #include "CandNtupleSR/NtpSRRecord.h"
00024 #include "CandNtupleSR/NtpSREvent.h"
00025 #include "CandNtupleSR/NtpSRTrack.h"
00026 #include "CandNtupleSR/NtpSRShower.h"
00027 #include "NueAna/SntpHelpers.h"
00028 #include "NueAna/NueRWHelpers.h"
00029 #include "NueAna/OscProb.h"
00030 
00031 using namespace NueRWHelpers;
00032 
00033 #define POTPERFARFILE 6.5e20
00034 #define POTPERNEARFILE 550*2.4e13
00035 
00036 JOBMODULE(NueSensitivity, "NueSensitivity",
00037           "Reads in ana_nue ntuple and pid ntuple");
00038 CVSID("$Id:");
00039 
00040 //......................................................................
00041 NueSensitivity::NueSensitivity()
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 }
00083 
00084 //......................................................................
00085 NueSensitivity::~NueSensitivity()
00086 {
00087   delete [] theta13;
00088   delete [] delta23;
00089   delete nueAppear;
00090   delete numuSurvive;
00091 }
00092 
00093 //......................................................................
00094 const Registry& NueSensitivity::DefaultConfig() const
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 }
00117 
00118 //......................................................................
00119 void NueSensitivity::Config(const Registry& r)
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 }
00132 
00133 
00134 //......................................................................
00135 void NueSensitivity::SetPOT(){
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 }
00150 
00151 //......................................................................
00152 void NueSensitivity::BeginJob(){
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 }
00218 
00219 //......................................................................
00220 JobCResult NueSensitivity::Ana(const MomNavigator* mom)
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 }
00528 
00529 //......................................................................
00530 void NueSensitivity::Analysis()
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 }
00761 
00762 void NueSensitivity::EndJob()
00763 {
00764   Analysis();
00765   HistMan man("sensitivity");
00766   man.WriteOut("SensitivityFile.root");
00767   systematicFile->Close();
00768 }
00769 
00770
Archive/NueSensitivity.cxx
