void MakeHelper(Int_t whichPID,Int_t whichSys) { gROOT->LoadMacro("NueAna/macros/LoadLibs.C"); LoadLibs(); gSystem->Load("libNueAnaExtrapolation.so"); MCReweight *mcr = &MCReweight::Instance(); NeugenWeightCalculator *n = new NeugenWeightCalculator(); mcr->AddWeightCalculator(n); TChain *nearChain = new TChain("ana_nue"); TChain *farChain = new TChain("ana_nue"); nearChain->Add("/unix/minos/cbs/AnaNue/mc/Trimmed/AnaNue-n130*.root"); farChain->Add("/unix/minos/cbs/AnaNue/mc/Trimmed/AnaNue-f21*.root"); Double_t farPOT = 981.*1.02716e20/3.; Double_t nearPOT = 2982*400*2.568e13; /* nearChain->Add("/unix/minos/cbs/AnaNue/mc/near/AnaNue-n13022000_0000_L010185.sntp.R1_18_2.root"); farChain->Add("/unix/minos/cbs/AnaNue/mc/far/AnaNue-f21[0-3]01001_0000_L010185.sntp.R1_18_2.root"); Double_t farPOT = farChain->GetNtrees()*1.02716e20/3.; Double_t nearPOT = nearChain->GetNtrees()*400*2.568e13; */ //configure what changes you want to study: Int_t nNorm = 8; //+/- 4% Double_t Norm[8] = {0.96,0.97,0.98,0.99,1.01,1.02,1.03,1.04}; Int_t nEMEnergy = 4; // +/- 2% on electron energy Double_t EMEnergy[4] = {-0.02,-0.01,0.01,0.02}; Int_t nHadEnergy = 10; //+/- 11% on hadronic system Double_t HadEnergy[10] = {-0.11,-0.09,-0.07,-0.05,-0.03, 0.03,0.05,0.07,0.09,0.11}; Int_t nRelEnergy = 4; //+/- 2% on N->F relative calibration Double_t RelEnergy[4] = {-0.02,-0.01,0.01,0.02}; Int_t nMA_QEL = 10; //+/- 30% Double_t MA_QEL[10] = {0.7,0.76,0.82,0.88,0.94,1.06,1.12,1.18,1.24,1.3}; Int_t nMA_RES = 10; //+/- 30% Double_t MA_RES[10] = {0.7,0.76,0.82,0.88,0.94,1.06,1.12,1.18,1.24,1.3}; Int_t nKNO = 10; //+/- 30% Double_t KNO[10] = {0.7,0.76,0.82,0.88,0.94,1.06,1.12,1.18,1.24,1.3}; Int_t nTrkPlane = 4; //+/- 2 planes Double_t TrkPlane[4] = {-2,-1,1,2}; Int_t nPID = 10; // +/- 5% on pid cut Double_t PID[10] = {-0.05,-0.04,-0.03,-0.02,-0.01,0.01,0.02,0.03,0.04,0.05}; Int_t ndelta23 = 10; // 90% MINOS confidence limits Double_t delta23[10] = {0.0022,0.0023,0.0024,0.0025,0.0026, 0.0028,0.0029,0.0030,0.0031,0.0032}; Int_t ntheta23 = 10; //90% MINOS confidence limits Double_t theta23[10] = {28,32,36,40,44,46,50,54,58,62}; for(int i=0;i<10;i++) theta23[i] *= TMath::Pi()/180.; Int_t nShwDev = 3; //numbers are indexes: 1 = scale to daikon spectrum; 2 = MODBYRS4; 3 = PiZero reweighting Double_t ShwDev[3] = {1,2,3}; Int_t nTauProd = 10; //tau neutrino CC cross-section uncertainty of +/-50% Double_t TauProd[10] = {0.5,0.6,0.7,0.8,0.9,1.1,1.2,1.3,1.4,1.5}; Int_t nPIDSkew = 10; //Skew PID dists to weight high PID events up to +/-50% Double_t PIDSkew[10] = {-0.5,-0.4,-0.3,-0.2,-0.1,0.1,0.2,0.3,0.4,0.5}; Selection::Selection_t sel = whichPID; Systematic::Systematic_t sys = whichSys; //make plot maker helper: NueFNHelper help(200,0.,100.); help.SetChains(nearChain,farChain,nearPOT,farPOT); //add a standard entry: (numu CC osc pars, no nue oscillations; SKZP out to 30GeV) NueSystematic *nueSys = new NueSystematic("Standard"); nueSys->AddSystematic(Systematic::kSKZP,Systematic::GetDefaultValue(Systematic::kSKZP)); nueSys->AddSystematic(Systematic::kOscProb,Systematic::GetDefaultValue(Systematic::kOscProb)); nueSys->SetOscParams(0.554,0.7854,0.,8.2e-5,2.7e-3,0,1); if(sys!=Systematic::kSKZP && sys!=Systematic::kOscProb) nueSys->AddSystematic(sys,Systematic::GetDefaultValue(sys)); help.AddNueSystematic(nueSys); //add systematic tweaks: if(sys==Systematic::kOscProb) { //special case for(int i=0;iAddSystematic(Systematic::kSKZP,1); nueSys->AddSystematic(sys,1); nueSys->SetOscParams(0.554,theta23[i],0.,8.2e-5,2.7e-3,0,1); help.AddNueSystematic(nueSys); } for(int i=0;iAddSystematic(Systematic::kSKZP,1); nueSys->AddSystematic(sys,1); nueSys->SetOscParams(0.554,0.7854,0.,8.2e-5,delta23[i],0,1); help.AddNueSystematic(nueSys); } } else if(sys==Systematic::kSKZP) { //special case //switch off SKZP char sysname[256]; sprintf(sysname,"AlterSKZP_0"); nueSys = new NueSystematic(sysname); nueSys->AddSystematic(Systematic::kSKZP,0); nueSys->AddSystematic(Systematic::kOscProb,1); nueSys->SetOscParams(0.554,0.7854,0.,8.2e-5,2.7e-3,0,1); help.AddNueSystematic(nueSys); //use SKZP_120GeV sprintf(sysname,"AlterSKZP_2"); nueSys = new NueSystematic(sysname); nueSys->AddSystematic(Systematic::kSKZP,2); nueSys->AddSystematic(Systematic::kOscProb,1); nueSys->SetOscParams(0.554,0.7854,0.,8.2e-5,2.7e-3,0,1); help.AddNueSystematic(nueSys); } else { Int_t nmax = 10; for(int i=0;iAddSystematic(Systematic::kSKZP,1); nueSys->AddSystematic(Systematic::kOscProb,1); nueSys->SetOscParams(0.554,0.7854,0.,8.2e-5,2.7e-3,0,1); if(sys==Systematic::kNorm) { nueSys->AddSystematic(sys,Norm[i]); if(i>=nNorm-1) i=nmax; } if(sys==Systematic::kEMCalib) { nueSys->AddSystematic(sys,EMEnergy[i]); if(i>=nEMEnergy-1) i=nmax; } else if(sys==Systematic::kHadCalib) { nueSys->AddSystematic(sys,HadEnergy[i]); if(i>=nHadEnergy-1) i=nmax; } else if(sys==Systematic::kRelCalib) { nueSys->AddSystematic(sys,RelEnergy[i]); if(i>=nRelEnergy-1) i=nmax; } else if(sys==Systematic::kMA_QE) { nueSys->AddSystematic(sys,MA_QEL[i]); if(i>=nMA_QEL-1) i=nmax; } else if(sys==Systematic::kMA_RES) { nueSys->AddSystematic(sys,MA_RES[i]); if(i>=nMA_RES-1) i=nmax; } else if(sys==Systematic::kKNO) { nueSys->AddSystematic(sys,KNO[i]); if(i>=nKNO-1) i=nmax; } else if(sys==Systematic::kTrkPlane) { nueSys->AddSystematic(sys,TrkPlane[i]); if(i>=nTrkPlane-1) i=nmax; } else if(sys==Systematic::kPIDShift) { nueSys->AddSystematic(sys,PID[i]); if(i>=nPID-1) i=nmax; } else if(sys==Systematic::kShwDev) { nueSys->AddSystematic(sys,ShwDev[i]); if(i>=nShwDev-1) i=nmax; } else if(sys==Systematic::kTauProd) { nueSys->AddSystematic(sys,TauProd[i]); if(i>=nTauProd-1) i=nmax; } else if(sys==Systematic::kPIDSkew) { nueSys->AddSystematic(sys,PIDSkew[i]); if(i>=nPIDSkew-1) i=nmax; } help.AddNueSystematic(nueSys); } } //run all the plot making code: help.MakeHelpers(sel); //write out the file: string outfiletag = string(Selection::AsString(sel)) + "_Alter" + string(Systematic::AsString(sys)); help.WriteFile(outfiletag.c_str()); }