GridGen_Joint Class Reference

#include <GridGen_Joint.h>

Inheritance diagram for GridGen_Joint:
GridGen

List of all members.

Public Member Functions

 GridGen_Joint ()
virtual ~GridGen_Joint ()
void AddExtrapFHC (Extrapolate2D *E)
void AddExtrapRHC (Extrapolate2D *E)
void RunMultiBinOscParErrs (string s="OscParErrDistributions.root")
void RunMultiBin_VaryTheta13 (string s="OscParErrDistributions.root")

Private Attributes

vector< Extrapolate2D * > ExtrapRHC
vector< Extrapolate2D * > ExtrapFHC

Detailed Description

Definition at line 16 of file GridGen_Joint.h.


Constructor & Destructor Documentation

GridGen_Joint::GridGen_Joint (  ) 
GridGen_Joint::~GridGen_Joint (  )  [virtual]

Definition at line 26 of file GridGen_Joint.cxx.

00027 {
00028 }


Member Function Documentation

void GridGen_Joint::AddExtrapFHC ( Extrapolate2D E  ) 

Definition at line 30 of file GridGen_Joint.cxx.

References ExtrapFHC.

00031 {
00032   ExtrapFHC.push_back(E);
00033   return;
00034 }

void GridGen_Joint::AddExtrapRHC ( Extrapolate2D E  ) 

Definition at line 35 of file GridGen_Joint.cxx.

References ExtrapRHC.

00036 {
00037   ExtrapRHC.push_back(E);
00038   return;
00039 }

void GridGen_Joint::RunMultiBin_VaryTheta13 ( string  s = "OscParErrDistributions.root"  )  [virtual]

Reimplemented from GridGen.

Definition at line 491 of file GridGen_Joint.cxx.

References GridGen::AsymGaus(), NueConvention::bnue, GridGen::dDeltaMSq12_dn, GridGen::dDeltaMSq12_up, GridGen::dDeltaMSq23_dn, GridGen::dDeltaMSq23_up, GridGen::DeltaHigh, GridGen::DeltaLow, GridGen::DeltaMSq12, GridGen::DeltaMSq23, GridGen::dTheta12_dn, GridGen::dTheta12_up, GridGen::dTheta13_dn, GridGen::dTheta13_up, GridGen::dTheta23_dn, GridGen::dTheta23_up, ExtrapFHC, ExtrapRHC, Form(), gSystem(), id, Background::kBNueCC, OscPar::kDeltaM12, OscPar::kDeltaM23, Background::kNC, Background::kNueCC, Background::kNuMuCC, Background::kNuTauCC, OscPar::kTh12, OscPar::kTh13, OscPar::kTh23, GridGen::nDeltaSteps, GridGen::NormalHier, NueConvention::nue, GridGen::NumExpts, GridGen::outFileName, GridGen::Theta12, GridGen::Theta13, and GridGen::Theta23.

00492 {
00493   if((ExtrapRHC.size() + ExtrapFHC.size())==0)
00494   {
00495     cout<<"No Extrapolate2D input.  Quitting..."<<endl;
00496     return;
00497   }
00498   
00499   unsigned int ie;
00500   for(ie=0;ie<ExtrapFHC.size();ie++)
00501   {
00502     ExtrapFHC[ie]->GetPrediction();
00503   }
00504   for(ie=0;ie<ExtrapRHC.size();ie++)
00505   {
00506     ExtrapRHC[ie]->GetPrediction();
00507   }
00508   
00509   for(ie=0;ie<ExtrapFHC.size();ie++)
00510   {
00511     ExtrapFHC[ie]->SetOscPar(OscPar::kTh13,Theta13);
00512     ExtrapFHC[ie]->SetOscPar(OscPar::kTh12,Theta12);
00513     ExtrapFHC[ie]->SetOscPar(OscPar::kTh23,Theta23);
00514     ExtrapFHC[ie]->SetOscPar(OscPar::kDeltaM12,DeltaMSq12);
00515     ExtrapFHC[ie]->SetOscPar(OscPar::kDeltaM23,DeltaMSq23);
00516     if(!NormalHier) ExtrapFHC[ie]->InvertMassHierarchy();
00517     ExtrapFHC[ie]->OscillatePrediction();
00518   }
00519   for(ie=0;ie<ExtrapRHC.size();ie++)
00520   {
00521     ExtrapRHC[ie]->SetOscPar(OscPar::kTh13,Theta13);
00522     ExtrapRHC[ie]->SetOscPar(OscPar::kTh12,Theta12);
00523     ExtrapRHC[ie]->SetOscPar(OscPar::kTh23,Theta23);
00524     ExtrapRHC[ie]->SetOscPar(OscPar::kDeltaM12,DeltaMSq12);
00525     ExtrapRHC[ie]->SetOscPar(OscPar::kDeltaM23,DeltaMSq23);
00526     if(!NormalHier) ExtrapRHC[ie]->InvertMassHierarchy();
00527     ExtrapRHC[ie]->OscillatePrediction();
00528   }
00529   int i,k;
00530   unsigned int j;
00531   
00532   int nbins=0;
00533   int nbinsFHC=0;
00534   int nbinsRHC=0;
00535   if (ExtrapFHC.size()>0) nbinsFHC = ExtrapFHC[0]->Pred_TotalBkgd_VsBinNumber->GetNbinsX();
00536   if (ExtrapRHC.size()>0) nbinsRHC = ExtrapRHC[0]->Pred_TotalBkgd_VsBinNumber->GetNbinsX();
00537   nbins = nbinsRHC + nbinsFHC;
00538   
00539   int nPIDFHC=0;
00540   int nPIDRHC=0;
00541   if (ExtrapFHC.size()>0) nPIDFHC = ExtrapFHC[0]->GetNPID();
00542   if (ExtrapRHC.size()>0) nPIDRHC = ExtrapRHC[0]->GetNPID();
00543   double delta,t13axis;
00544   vector<double> sig, bkgd;
00545   vector< vector<double> > nc,numucc,bnue,tau,nue;
00546   vector<double> oscparerr;
00547   vector<double> oscparerr_offdiag;
00548   int noff;
00549   
00550   t13axis = TMath::Sin(2*Theta13)*TMath::Sin(2*Theta13);
00551   
00552   for(i=0;i<nbins;i++)
00553   {
00554     sig.push_back(0);
00555     bkgd.push_back(0);
00556     oscparerr.push_back(0);
00557     nc.push_back( vector<double>() );
00558     numucc.push_back( vector<double>() );
00559     bnue.push_back( vector<double>() );
00560     tau.push_back( vector<double>() );
00561     nue.push_back( vector<double>() );
00562     for(k=0;k<nbins;k++)
00563     {
00564       if(k>i)
00565       {
00566         oscparerr_offdiag.push_back(0);
00567       }
00568     }
00569     if (i<nbinsFHC){
00570       for(j=0;j<ExtrapFHC.size();j++)
00571       {
00572         nc[i].push_back(0);
00573         numucc[i].push_back(0);
00574         bnue[i].push_back(0);
00575         tau[i].push_back(0);
00576         nue[i].push_back(0);
00577       } 
00578     } else if (i>=nbinsFHC){
00579       for(j=0;j<ExtrapRHC.size();j++)
00580       {
00581         nc[i].push_back(0);
00582         numucc[i].push_back(0);
00583         bnue[i].push_back(0);
00584         tau[i].push_back(0);
00585         nue[i].push_back(0);
00586       }
00587     }
00588   }
00589   
00590   vector<TTree*> ftree;
00591   vector< vector<TTree*> > ftree2;
00592   TTree *ttmp;
00593   noff=0;
00594   for(i=0;i<nbins;i++)
00595   {
00596     ttmp = new TTree(Form("Bin_%i",i),Form("Bin_%i",i));
00597     ttmp->Branch("Delta",&delta,"Delta/D");
00598     ttmp->Branch("Th13Axis",&t13axis,"Th13Axis/D");
00599     ttmp->Branch("Signal",&sig[i],"Signal/D");
00600     ttmp->Branch("Background",&bkgd[i],"Background/D");
00601     ttmp->Branch("DNExp_DOscPars",&oscparerr[i],"DNExp_DOscPars/D");
00602     for(k=0;k<nbins;k++)
00603     {
00604       if(k>i)
00605       {
00606         ttmp->Branch(Form("Bin_%i_Bin_%i",i,k),&oscparerr_offdiag[noff],Form("Bin_%i_Bin_%i/D",i,k));
00607         noff++;
00608       }
00609     }
00610     ftree.push_back(ttmp);
00611     
00612     ttmp->Reset();
00613     
00614     ftree2.push_back( vector<TTree*>() );
00615     
00616     if (i<nbinsFHC){
00617       
00618       for(j=0;j<ExtrapFHC.size();j++)
00619       {
00620         ttmp = new TTree(Form("Bin_%i_Run_%i",i,j),Form("Bin_%i_Run_%i",i,j));
00621         ttmp->Branch("Delta",&delta,"Delta/D");
00622         ttmp->Branch("Th13Axis",&t13axis,"Th13Axis/D");
00623         ttmp->Branch("Signal",&nue[i][j],"Signal/D");
00624         ttmp->Branch("NC",&nc[i][j],"NC/D");
00625         ttmp->Branch("NuMuCC",&numucc[i][j],"NuMuCC/D");
00626         ttmp->Branch("BNueCC",&bnue[i][j],"BNueCC/D");
00627         ttmp->Branch("NuTauCC",&tau[i][j],"NuTauCC/D");
00628         
00629         ftree2[i].push_back(ttmp);
00630         
00631         ttmp->Reset();
00632       }
00633     } else if (i>=nbinsFHC){
00634       
00635       for(j=0;j<ExtrapRHC.size();j++)
00636       {
00637         ttmp = new TTree(Form("Bin_%i_Run_%i",i,j),Form("Bin_%i_Run_%i",i,j));
00638         ttmp->Branch("Delta",&delta,"Delta/D");
00639         ttmp->Branch("Th13Axis",&t13axis,"Th13Axis/D");
00640         ttmp->Branch("Signal",&nue[i][j],"Signal/D");
00641         ttmp->Branch("NC",&nc[i][j],"NC/D");
00642         ttmp->Branch("NuMuCC",&numucc[i][j],"NuMuCC/D");
00643         ttmp->Branch("BNueCC",&bnue[i][j],"BNueCC/D");
00644         ttmp->Branch("NuTauCC",&tau[i][j],"NuTauCC/D");
00645         
00646         ftree2[i].push_back(ttmp);
00647         
00648         ttmp->Reset();
00649       }
00650     }
00651   }
00652   
00653   double delta_increment = 0;
00654   if(nDeltaSteps>0) delta_increment = (DeltaHigh - DeltaLow)/(nDeltaSteps);
00655   
00656   int id,l,u;
00657   int ip,ir;
00658   double theta23,theta12,dm21,dm32,theta13;
00659   vector<double> nexp,nobs;
00660   vector<TH1D*> delnexphist;
00661   vector<TH1D*> delnexphist_offdiag;
00662   for(i=0;i<nbins;i++)
00663   {
00664     delnexphist.push_back(new TH1D(Form("delnexphist_%i",i),"",400,-1,1));
00665     nexp.push_back(0);
00666     nobs.push_back(0);
00667     for(k=0;k<nbins;k++)
00668     {
00669       if(k>i)
00670       {
00671         delnexphist_offdiag.push_back(new TH1D(Form("delnexphist_%i_%i",i,k),"",400,-1,1));
00672       }
00673     }
00674   }
00675   
00676   gRandom->SetSeed(0);
00677   
00678   TFile *f = new TFile(gSystem->ExpandPathName(s.c_str()),"RECREATE");//save the osc par err distributions to this file
00679   
00680   l=0;
00681   for(id=0;id<nDeltaSteps+1;id++)
00682   {
00683     delta = (id*delta_increment + DeltaLow)*TMath::Pi();
00684     
00685     
00686     //get nominal prediction
00687     for(ie=0;ie<ExtrapFHC.size();ie++)
00688     {
00689       ExtrapFHC[ie]->SetOscPar(OscPar::kTh13,Theta13);
00690       ExtrapFHC[ie]->SetOscPar(OscPar::kTh12,Theta12);
00691       ExtrapFHC[ie]->SetOscPar(OscPar::kTh23,Theta23);
00692       ExtrapFHC[ie]->SetOscPar(OscPar::kDeltaM12,DeltaMSq12);
00693       ExtrapFHC[ie]->SetOscPar(OscPar::kDeltaM23,DeltaMSq23);//note that DeltaMSq23 is always positive
00694       if(!NormalHier) ExtrapFHC[ie]->InvertMassHierarchy();
00695       ExtrapFHC[ie]->SetDeltaCP(delta);
00696       ExtrapFHC[ie]->OscillatePrediction();
00697     }
00698     
00699     for(ie=0;ie<ExtrapRHC.size();ie++)
00700     {
00701       ExtrapRHC[ie]->SetOscPar(OscPar::kTh13,Theta13);
00702       ExtrapRHC[ie]->SetOscPar(OscPar::kTh12,Theta12);
00703       ExtrapRHC[ie]->SetOscPar(OscPar::kTh23,Theta23);
00704       ExtrapRHC[ie]->SetOscPar(OscPar::kDeltaM12,DeltaMSq12);
00705       ExtrapRHC[ie]->SetOscPar(OscPar::kDeltaM23,DeltaMSq23);//note that DeltaMSq23 is always positive
00706       if(!NormalHier) ExtrapRHC[ie]->InvertMassHierarchy();
00707       ExtrapRHC[ie]->SetDeltaCP(delta);
00708       ExtrapRHC[ie]->OscillatePrediction();
00709     }
00710     
00711     for(i=0;i<nbins;i++)
00712     {
00713       sig[i]=0;
00714       bkgd[i]=0;
00715       
00716       if (i<nbinsFHC){
00717         
00718         ir = int(i/nPIDFHC);
00719         ip = i%nPIDFHC;
00720         
00721         for(ie=0;ie<ExtrapFHC.size();ie++)
00722         {
00723           bkgd[i] += (ExtrapFHC[ie]->Pred_TotalBkgd_VsBinNumber->GetBinContent(i+1));
00724           sig[i] += (ExtrapFHC[ie]->Pred_Signal_VsBinNumber->GetBinContent(i+1));
00725           
00726           nc[i][ie] = ExtrapFHC[ie]->Pred[Background::kNC]->GetBinContent(ip+1,ir+1);
00727           numucc[i][ie] = ExtrapFHC[ie]->Pred[Background::kNuMuCC]->GetBinContent(ip+1,ir+1);
00728           bnue[i][ie] = ExtrapFHC[ie]->Pred[Background::kBNueCC]->GetBinContent(ip+1,ir+1);
00729           tau[i][ie] = ExtrapFHC[ie]->Pred[Background::kNuTauCC]->GetBinContent(ip+1,ir+1);
00730           nue[i][ie] = ExtrapFHC[ie]->Pred[Background::kNueCC]->GetBinContent(ip+1,ir+1);
00731         }
00732       } else if (i>=nbinsFHC){
00733         
00734         ir = int((i-nbinsFHC)/nPIDRHC);
00735         ip = (i-nbinsFHC)%nPIDRHC;
00736         
00737         for(ie=0;ie<ExtrapRHC.size();ie++)
00738         {
00739           bkgd[i] += (ExtrapRHC[ie]->Pred_TotalBkgd_VsBinNumber->GetBinContent(i+1-nbinsFHC));
00740           sig[i] += (ExtrapRHC[ie]->Pred_Signal_VsBinNumber->GetBinContent(i+1-nbinsFHC));
00741           
00742           nc[i][ie] = ExtrapRHC[ie]->Pred[Background::kNC]->GetBinContent(ip+1,ir+1);
00743           numucc[i][ie] = ExtrapRHC[ie]->Pred[Background::kNuMuCC]->GetBinContent(ip+1,ir+1);
00744           bnue[i][ie] = ExtrapRHC[ie]->Pred[Background::kBNueCC]->GetBinContent(ip+1,ir+1);
00745           tau[i][ie] = ExtrapRHC[ie]->Pred[Background::kNuTauCC]->GetBinContent(ip+1,ir+1);
00746           nue[i][ie] = ExtrapRHC[ie]->Pred[Background::kNueCC]->GetBinContent(ip+1,ir+1);
00747         }
00748       }
00749       
00750       nexp[i] = sig[i]+bkgd[i];
00751      }
00752       
00753       //do pseudo experiments
00754       noff=0;
00755       for(i=0;i<nbins;i++)
00756       {
00757         delnexphist[i]->Reset();
00758         delnexphist[i]->SetName(Form("DeltaNexp_%i_Diag_%i",id,i));
00759         for(k=0;k<nbins;k++)
00760         {
00761           if(k>i)
00762           {
00763             delnexphist_offdiag[noff]->Reset();
00764             delnexphist_offdiag[noff]->SetName(Form("DeltaNexp_%i_OffDiag_%i_%i",id,i,k));
00765             noff++;
00766           }
00767         }
00768       }
00769       
00770       //Do this NumExpts times
00771       for(u=0;u<NumExpts;u++)
00772       {
00773         theta13 = Theta13 + AsymGaus(dTheta13_dn,dTheta13_up);
00774         theta12 = Theta12 + AsymGaus(dTheta12_dn,dTheta12_up);
00775         dm21 = DeltaMSq12 + AsymGaus(dDeltaMSq12_dn,dDeltaMSq12_up);
00776         dm32 = DeltaMSq23 + AsymGaus(dDeltaMSq23_dn,dDeltaMSq23_up);
00777         theta23 = Theta23 + AsymGaus(dTheta23_dn,dTheta23_up);
00778         
00779         for(ie=0;ie<ExtrapFHC.size();ie++)
00780         {
00781           ExtrapFHC[ie]->SetOscPar(OscPar::kTh13,theta13);
00782           ExtrapFHC[ie]->SetOscPar(OscPar::kTh12,theta12);
00783           ExtrapFHC[ie]->SetOscPar(OscPar::kDeltaM12,dm21);
00784           ExtrapFHC[ie]->SetOscPar(OscPar::kDeltaM23,dm32);
00785           if(!NormalHier) ExtrapFHC[ie]->InvertMassHierarchy();
00786           ExtrapFHC[ie]->SetOscPar(OscPar::kTh23,theta23);
00787           ExtrapFHC[ie]->OscillatePrediction();
00788         }
00789         for(ie=0;ie<ExtrapRHC.size();ie++)
00790         {
00791           ExtrapRHC[ie]->SetOscPar(OscPar::kTh13,theta13);
00792           ExtrapRHC[ie]->SetOscPar(OscPar::kTh12,theta12);
00793           ExtrapRHC[ie]->SetOscPar(OscPar::kDeltaM12,dm21);
00794           ExtrapRHC[ie]->SetOscPar(OscPar::kDeltaM23,dm32);
00795           if(!NormalHier) ExtrapRHC[ie]->InvertMassHierarchy();
00796           ExtrapRHC[ie]->SetOscPar(OscPar::kTh23,theta23);
00797           ExtrapRHC[ie]->OscillatePrediction();
00798         }
00799         
00800         noff=0;
00801         for(i=0;i<nbins;i++)
00802         {
00803           nobs[i]=0;
00804           
00805           if (i<nbinsFHC){
00806             for(ie=0;ie<ExtrapFHC.size();ie++)
00807             {
00808               nobs[i] += (ExtrapFHC[ie]->Pred_TotalBkgd_VsBinNumber->GetBinContent(i+1));
00809               nobs[i] += (ExtrapFHC[ie]->Pred_Signal_VsBinNumber->GetBinContent(i+1));
00810             }
00811           } else if (i>=nbinsFHC){
00812             for(ie=0;ie<ExtrapRHC.size();ie++)
00813             {
00814               nobs[i] += (ExtrapRHC[ie]->Pred_TotalBkgd_VsBinNumber->GetBinContent(i+1-nbinsFHC));
00815               nobs[i] += (ExtrapRHC[ie]->Pred_Signal_VsBinNumber->GetBinContent(i+1-nbinsFHC));
00816             }
00817           }
00818           delnexphist[i]->Fill((nobs[i]-nexp[i])/(nexp[i]));
00819         }
00820         
00821         for(i=0;i<nbins;i++)
00822         {
00823           for(k=0;k<nbins;k++)
00824           {
00825             if(k>i)
00826             {
00827               delnexphist_offdiag[noff]->Fill((nobs[i]-nexp[i])*(nobs[k]-nexp[k])/(nexp[i]*nexp[k]));
00828               noff++;
00829             }
00830           }
00831         }
00832       }
00833       
00834       noff=0;
00835       for(i=0;i<nbins;i++)
00836       {
00837         oscparerr[i] = delnexphist[i]->GetRMS();
00838         delnexphist[i]->Write();
00839         for(k=0;k<nbins;k++)
00840         {
00841           if(k>i)
00842           {
00843             oscparerr_offdiag[noff] = delnexphist_offdiag[noff]->GetRMS();
00844             if(delnexphist_offdiag[noff]->GetMean()<0) oscparerr_offdiag[noff] = -1.*oscparerr_offdiag[noff];
00845             delnexphist_offdiag[noff]->Write();
00846             noff++;
00847           }
00848         }
00849       }
00850       f->Close();
00851       
00852       gROOT->cd("/");
00853       
00854       for(i=0;i<nbins;i++)
00855       {
00856         if (i<nbinsFHC){
00857           for(ie=0;ie<ExtrapFHC.size();ie++)
00858           {
00859             ftree2[i][ie]->Fill();
00860           }
00861         } else if (i>=nbinsFHC){
00862           for(ie=0;ie<ExtrapRHC.size();ie++)
00863           {
00864             ftree2[i][ie]->Fill();
00865           }
00866         }
00867         ftree[i]->Fill();
00868       }
00869       
00870       f = new TFile(gSystem->ExpandPathName(s.c_str()),"UPDATE");
00871       
00872       if(l%100==0) cout<<100.*l/(nDeltaSteps+1)<<"% complete"<<endl;
00873       l++;
00874     
00875   }
00876   
00877   f->Close();
00878   
00879   double nPOTNearFHC,nPOTFarFHC;
00880   double nPOTNearRHC,nPOTFarRHC;
00881   
00882   TTree *paramtree = new TTree("paramtree","paramtree");
00883   paramtree->Branch("nearPOTFHC",&nPOTNearFHC,"nearPOTFHC/D");
00884   paramtree->Branch("farPOTFHC",&nPOTFarFHC,"farPOTFHC/D");
00885   paramtree->Branch("nearPOTRHC",&nPOTNearRHC,"nearPOTRHC/D");
00886   paramtree->Branch("farPOTRHC",&nPOTFarRHC,"farPOTRHC/D");
00887   paramtree->Branch("Theta13",&Theta13,"Theta13/D");
00888   paramtree->Branch("Theta12",&Theta12,"Theta12/D");
00889   paramtree->Branch("Theta23",&Theta23,"Theta23/D");
00890   paramtree->Branch("DeltaMSq23",&dm32,"DeltaMSq23/D");
00891   paramtree->Branch("DeltaMSq12",&DeltaMSq12,"DeltaMSq12/D");
00892   
00893   dm32 = DeltaMSq23;
00894   if(!NormalHier) dm32 = -1.*DeltaMSq23;
00895   
00896   for(ie=0;ie<ExtrapFHC.size();ie++)
00897   {
00898     
00899     nPOTNearRHC=0;
00900     nPOTFarRHC=0;
00901     nPOTNearFHC = ExtrapFHC[ie]->GetNearPOT();
00902     nPOTFarFHC = ExtrapFHC[ie]->GetFarPOT();
00903     paramtree->Fill();
00904   }
00905   for(ie=0;ie<ExtrapRHC.size();ie++)
00906   {
00907     nPOTNearFHC=0;
00908     nPOTFarFHC=0;
00909     nPOTNearRHC=0;
00910     nPOTFarRHC=0;
00911     nPOTNearRHC = ExtrapRHC[ie]->GetNearPOT();
00912     nPOTFarRHC = ExtrapRHC[ie]->GetFarPOT();
00913     paramtree->Fill();
00914   }
00915   
00916   TFile *fout = new TFile(gSystem->ExpandPathName(outFileName.c_str()),"RECREATE");
00917   for(i=0;i<nbins;i++)
00918   {
00919     ftree[i]->Write();
00920     if (i<nbinsFHC){
00921       for(ie=0;ie<ExtrapFHC.size();ie++)
00922       {
00923         ftree2[i][ie]->Write();
00924       } 
00925     } else if (i>=nbinsFHC){
00926       for(ie=0;ie<ExtrapRHC.size();ie++)
00927       {
00928         ftree2[i][ie]->Write();
00929       } 
00930     }
00931   }
00932   paramtree->Write();
00933   fout->Close();
00934   
00935   return;
00936 }

void GridGen_Joint::RunMultiBinOscParErrs ( string  s = "OscParErrDistributions.root"  )  [virtual]

Reimplemented from GridGen.

Definition at line 42 of file GridGen_Joint.cxx.

References GridGen::AsymGaus(), NueConvention::bnue, GridGen::dDeltaMSq12_dn, GridGen::dDeltaMSq12_up, GridGen::dDeltaMSq23_dn, GridGen::dDeltaMSq23_up, GridGen::DeltaHigh, GridGen::DeltaLow, GridGen::DeltaMSq12, GridGen::DeltaMSq23, GridGen::dTheta12_dn, GridGen::dTheta12_up, GridGen::dTheta23_dn, GridGen::dTheta23_up, ExtrapFHC, ExtrapRHC, Form(), gSystem(), id, Background::kBNueCC, OscPar::kDeltaM12, OscPar::kDeltaM23, Background::kNC, Background::kNueCC, Background::kNuMuCC, Background::kNuTauCC, OscPar::kTh12, OscPar::kTh23, GridGen::nDeltaSteps, GridGen::NormalHier, GridGen::nSinSq2Th13Steps, NueConvention::nue, GridGen::NumExpts, GridGen::outFileName, GridGen::SinSq2Th13High, GridGen::SinSq2Th13Low, GridGen::Theta12, and GridGen::Theta23.

00043 {
00044   if((ExtrapRHC.size() + ExtrapFHC.size())==0)
00045   {
00046     cout<<"No Extrapolate2D input.  Quitting..."<<endl;
00047     return;
00048   }
00049   
00050   unsigned int ie;
00051   for(ie=0;ie<ExtrapFHC.size();ie++)
00052   {
00053     ExtrapFHC[ie]->GetPrediction();
00054   }
00055   for(ie=0;ie<ExtrapRHC.size();ie++)
00056   {
00057     ExtrapRHC[ie]->GetPrediction();
00058   }
00059   
00060   for(ie=0;ie<ExtrapFHC.size();ie++)
00061   {
00062     ExtrapFHC[ie]->SetOscPar(OscPar::kTh12,Theta12);
00063     ExtrapFHC[ie]->SetOscPar(OscPar::kTh23,Theta23);
00064     ExtrapFHC[ie]->SetOscPar(OscPar::kDeltaM12,DeltaMSq12);
00065     ExtrapFHC[ie]->SetOscPar(OscPar::kDeltaM23,DeltaMSq23);
00066     if(!NormalHier) ExtrapFHC[ie]->InvertMassHierarchy();
00067     ExtrapFHC[ie]->OscillatePrediction();
00068   }
00069     for(ie=0;ie<ExtrapRHC.size();ie++)
00070   {
00071     ExtrapRHC[ie]->SetOscPar(OscPar::kTh12,Theta12);
00072     ExtrapRHC[ie]->SetOscPar(OscPar::kTh23,Theta23);
00073     ExtrapRHC[ie]->SetOscPar(OscPar::kDeltaM12,DeltaMSq12);
00074     ExtrapRHC[ie]->SetOscPar(OscPar::kDeltaM23,DeltaMSq23);
00075     if(!NormalHier) ExtrapRHC[ie]->InvertMassHierarchy();
00076     ExtrapRHC[ie]->OscillatePrediction();
00077   }
00078   int i,k;
00079   unsigned int j;
00080 
00081   int nbins=0;
00082   int nbinsFHC=0;
00083   int nbinsRHC=0;
00084   if (ExtrapFHC.size()>0) nbinsFHC = ExtrapFHC[0]->Pred_TotalBkgd_VsBinNumber->GetNbinsX();
00085   if (ExtrapRHC.size()>0) nbinsRHC = ExtrapRHC[0]->Pred_TotalBkgd_VsBinNumber->GetNbinsX();
00086   nbins = nbinsRHC + nbinsFHC;
00087 
00088   int nPIDFHC=0;
00089   int nPIDRHC=0;
00090  if (ExtrapFHC.size()>0) nPIDFHC = ExtrapFHC[0]->GetNPID();
00091  if (ExtrapRHC.size()>0) nPIDRHC = ExtrapRHC[0]->GetNPID();
00092   double delta,t13axis;
00093   vector<double> sig, bkgd;
00094   vector< vector<double> > nc,numucc,bnue,tau,nue;
00095   vector<double> oscparerr;
00096   vector<double> oscparerr_offdiag;
00097   int noff;
00098   
00099   for(i=0;i<nbins;i++)
00100   {
00101     sig.push_back(0);
00102     bkgd.push_back(0);
00103     oscparerr.push_back(0);
00104     nc.push_back( vector<double>() );
00105     numucc.push_back( vector<double>() );
00106     bnue.push_back( vector<double>() );
00107     tau.push_back( vector<double>() );
00108     nue.push_back( vector<double>() );
00109     for(k=0;k<nbins;k++)
00110     {
00111       if(k>i)
00112       {
00113         oscparerr_offdiag.push_back(0);
00114       }
00115     }
00116     if (i<nbinsFHC){
00117       for(j=0;j<ExtrapFHC.size();j++)
00118         {
00119           nc[i].push_back(0);
00120           numucc[i].push_back(0);
00121           bnue[i].push_back(0);
00122           tau[i].push_back(0);
00123           nue[i].push_back(0);
00124         } 
00125     } else if (i>=nbinsFHC){
00126       for(j=0;j<ExtrapRHC.size();j++)
00127         {
00128           nc[i].push_back(0);
00129           numucc[i].push_back(0);
00130           bnue[i].push_back(0);
00131           tau[i].push_back(0);
00132           nue[i].push_back(0);
00133         }
00134     }
00135   }
00136   
00137   vector<TTree*> ftree;
00138   vector< vector<TTree*> > ftree2;
00139   TTree *ttmp;
00140   noff=0;
00141   for(i=0;i<nbins;i++)
00142   {
00143     ttmp = new TTree(Form("Bin_%i",i),Form("Bin_%i",i));
00144     ttmp->Branch("Delta",&delta,"Delta/D");
00145     ttmp->Branch("Th13Axis",&t13axis,"Th13Axis/D");
00146     ttmp->Branch("Signal",&sig[i],"Signal/D");
00147     ttmp->Branch("Background",&bkgd[i],"Background/D");
00148     ttmp->Branch("DNExp_DOscPars",&oscparerr[i],"DNExp_DOscPars/D");
00149     for(k=0;k<nbins;k++)
00150     {
00151       if(k>i)
00152       {
00153         ttmp->Branch(Form("Bin_%i_Bin_%i",i,k),&oscparerr_offdiag[noff],Form("Bin_%i_Bin_%i/D",i,k));
00154         noff++;
00155       }
00156     }
00157     ftree.push_back(ttmp);
00158     
00159     ttmp->Reset();
00160     
00161     ftree2.push_back( vector<TTree*>() );
00162     
00163    if (i<nbinsFHC){
00164  
00165      for(j=0;j<ExtrapFHC.size();j++)
00166        {
00167          ttmp = new TTree(Form("Bin_%i_Run_%i",i,j),Form("Bin_%i_Run_%i",i,j));
00168          ttmp->Branch("Delta",&delta,"Delta/D");
00169          ttmp->Branch("Th13Axis",&t13axis,"Th13Axis/D");
00170          ttmp->Branch("Signal",&nue[i][j],"Signal/D");
00171          ttmp->Branch("NC",&nc[i][j],"NC/D");
00172          ttmp->Branch("NuMuCC",&numucc[i][j],"NuMuCC/D");
00173          ttmp->Branch("BNueCC",&bnue[i][j],"BNueCC/D");
00174          ttmp->Branch("NuTauCC",&tau[i][j],"NuTauCC/D");
00175          
00176          ftree2[i].push_back(ttmp);
00177          
00178          ttmp->Reset();
00179        }
00180    } else if (i>=nbinsFHC){
00181 
00182      for(j=0;j<ExtrapRHC.size();j++)
00183        {
00184          ttmp = new TTree(Form("Bin_%i_Run_%i",i,j),Form("Bin_%i_Run_%i",i,j));
00185          ttmp->Branch("Delta",&delta,"Delta/D");
00186          ttmp->Branch("Th13Axis",&t13axis,"Th13Axis/D");
00187          ttmp->Branch("Signal",&nue[i][j],"Signal/D");
00188          ttmp->Branch("NC",&nc[i][j],"NC/D");
00189          ttmp->Branch("NuMuCC",&numucc[i][j],"NuMuCC/D");
00190          ttmp->Branch("BNueCC",&bnue[i][j],"BNueCC/D");
00191          ttmp->Branch("NuTauCC",&tau[i][j],"NuTauCC/D");
00192          
00193          ftree2[i].push_back(ttmp);
00194          
00195          ttmp->Reset();
00196        }
00197    }
00198   }
00199   
00200   double delta_increment = 0;
00201   if(nDeltaSteps>0) delta_increment = (DeltaHigh - DeltaLow)/(nDeltaSteps);
00202   double ssq2th13_increment = 0;
00203   if(nSinSq2Th13Steps>0) ssq2th13_increment = (SinSq2Th13High - SinSq2Th13Low)/(nSinSq2Th13Steps);
00204   
00205   int id,is,l,u;
00206   int ip,ir;
00207   double theta23,theta12,dm21,dm32,ssq2th13;
00208   vector<double> nexp,nobs;
00209   vector<TH1D*> delnexphist;
00210   vector<TH1D*> delnexphist_offdiag;
00211   for(i=0;i<nbins;i++)
00212   {
00213     delnexphist.push_back(new TH1D(Form("delnexphist_%i",i),"",400,-1,1));
00214     nexp.push_back(0);
00215     nobs.push_back(0);
00216     for(k=0;k<nbins;k++)
00217     {
00218       if(k>i)
00219       {
00220         delnexphist_offdiag.push_back(new TH1D(Form("delnexphist_%i_%i",i,k),"",400,-1,1));
00221       }
00222     }
00223   }
00224   
00225   gRandom->SetSeed(0);
00226   
00227   TFile *f = new TFile(gSystem->ExpandPathName(s.c_str()),"RECREATE");//save the osc par err distributions to this file
00228   
00229   l=0;
00230   for(id=0;id<nDeltaSteps+1;id++)
00231   {
00232     delta = (id*delta_increment + DeltaLow)*TMath::Pi();
00233     
00234     for(is=0;is<nSinSq2Th13Steps+1;is++)
00235     {
00236       t13axis = (is*ssq2th13_increment + SinSq2Th13Low);
00237       ssq2th13 = t13axis/(2.*TMath::Sin(Theta23)*TMath::Sin(Theta23));
00238       
00239       //get nominal prediction
00240       for(ie=0;ie<ExtrapFHC.size();ie++)
00241         {
00242           ExtrapFHC[ie]->SetOscPar(OscPar::kTh12,Theta12);
00243           ExtrapFHC[ie]->SetOscPar(OscPar::kTh23,Theta23);
00244           ExtrapFHC[ie]->SetOscPar(OscPar::kDeltaM12,DeltaMSq12);
00245           ExtrapFHC[ie]->SetOscPar(OscPar::kDeltaM23,DeltaMSq23);//note that DeltaMSq23 is always positive
00246           if(!NormalHier) ExtrapFHC[ie]->InvertMassHierarchy();
00247           ExtrapFHC[ie]->SetDeltaCP(delta);
00248           ExtrapFHC[ie]->SetSinSq2Th13(ssq2th13);
00249           ExtrapFHC[ie]->OscillatePrediction();
00250         }
00251       
00252       for(ie=0;ie<ExtrapRHC.size();ie++)
00253         {
00254           ExtrapRHC[ie]->SetOscPar(OscPar::kTh12,Theta12);
00255           ExtrapRHC[ie]->SetOscPar(OscPar::kTh23,Theta23);
00256           ExtrapRHC[ie]->SetOscPar(OscPar::kDeltaM12,DeltaMSq12);
00257           ExtrapRHC[ie]->SetOscPar(OscPar::kDeltaM23,DeltaMSq23);//note that DeltaMSq23 is always positive
00258           if(!NormalHier) ExtrapRHC[ie]->InvertMassHierarchy();
00259           ExtrapRHC[ie]->SetDeltaCP(delta);
00260           ExtrapRHC[ie]->SetSinSq2Th13(ssq2th13);
00261           ExtrapRHC[ie]->OscillatePrediction();
00262         }
00263 
00264       for(i=0;i<nbins;i++)
00265       {
00266         sig[i]=0;
00267         bkgd[i]=0;
00268 
00269         if (i<nbinsFHC){
00270 
00271           ir = int(i/nPIDFHC);
00272           ip = i%nPIDFHC;
00273 
00274           for(ie=0;ie<ExtrapFHC.size();ie++)
00275             {
00276               bkgd[i] += (ExtrapFHC[ie]->Pred_TotalBkgd_VsBinNumber->GetBinContent(i+1));
00277               sig[i] += (ExtrapFHC[ie]->Pred_Signal_VsBinNumber->GetBinContent(i+1));
00278               
00279               nc[i][ie] = ExtrapFHC[ie]->Pred[Background::kNC]->GetBinContent(ip+1,ir+1);
00280               numucc[i][ie] = ExtrapFHC[ie]->Pred[Background::kNuMuCC]->GetBinContent(ip+1,ir+1);
00281               bnue[i][ie] = ExtrapFHC[ie]->Pred[Background::kBNueCC]->GetBinContent(ip+1,ir+1);
00282               tau[i][ie] = ExtrapFHC[ie]->Pred[Background::kNuTauCC]->GetBinContent(ip+1,ir+1);
00283               nue[i][ie] = ExtrapFHC[ie]->Pred[Background::kNueCC]->GetBinContent(ip+1,ir+1);
00284             }
00285         } else if (i>=nbinsFHC){
00286 
00287           ir = int((i-nbinsFHC)/nPIDRHC);
00288           ip = (i-nbinsFHC)%nPIDRHC;
00289 
00290           for(ie=0;ie<ExtrapRHC.size();ie++)
00291             {
00292               bkgd[i] += (ExtrapRHC[ie]->Pred_TotalBkgd_VsBinNumber->GetBinContent(i+1-nbinsFHC));
00293               sig[i] += (ExtrapRHC[ie]->Pred_Signal_VsBinNumber->GetBinContent(i+1-nbinsFHC));
00294               
00295               nc[i][ie] = ExtrapRHC[ie]->Pred[Background::kNC]->GetBinContent(ip+1,ir+1);
00296               numucc[i][ie] = ExtrapRHC[ie]->Pred[Background::kNuMuCC]->GetBinContent(ip+1,ir+1);
00297               bnue[i][ie] = ExtrapRHC[ie]->Pred[Background::kBNueCC]->GetBinContent(ip+1,ir+1);
00298               tau[i][ie] = ExtrapRHC[ie]->Pred[Background::kNuTauCC]->GetBinContent(ip+1,ir+1);
00299               nue[i][ie] = ExtrapRHC[ie]->Pred[Background::kNueCC]->GetBinContent(ip+1,ir+1);
00300             }
00301         }
00302 
00303         nexp[i] = sig[i]+bkgd[i];
00304       }
00305       
00306       //do pseudo experiments
00307       noff=0;
00308       for(i=0;i<nbins;i++)
00309       {
00310         delnexphist[i]->Reset();
00311         delnexphist[i]->SetName(Form("DeltaNexp_%i_%i_Diag_%i",id,is,i));
00312         for(k=0;k<nbins;k++)
00313         {
00314           if(k>i)
00315           {
00316             delnexphist_offdiag[noff]->Reset();
00317             delnexphist_offdiag[noff]->SetName(Form("DeltaNexp_%i_%i_OffDiag_%i_%i",id,is,i,k));
00318             noff++;
00319           }
00320         }
00321       }
00322 
00323       //Do this NumExpts times
00324       for(u=0;u<NumExpts;u++)
00325       {
00326         theta12 = Theta12 + AsymGaus(dTheta12_dn,dTheta12_up);
00327         dm21 = DeltaMSq12 + AsymGaus(dDeltaMSq12_dn,dDeltaMSq12_up);
00328         dm32 = DeltaMSq23 + AsymGaus(dDeltaMSq23_dn,dDeltaMSq23_up);
00329         //theta23 = DrawTheta23(dTheta23_up);
00330         theta23 = Theta23 + AsymGaus(dTheta23_dn,dTheta23_up);
00331         
00332         ssq2th13 = t13axis/(2.*TMath::Sin(theta23)*TMath::Sin(theta23));
00333 
00334         for(ie=0;ie<ExtrapFHC.size();ie++)
00335         {
00336           ExtrapFHC[ie]->SetOscPar(OscPar::kTh12,theta12);
00337           ExtrapFHC[ie]->SetOscPar(OscPar::kDeltaM12,dm21);
00338           ExtrapFHC[ie]->SetOscPar(OscPar::kDeltaM23,dm32);
00339           if(!NormalHier) ExtrapFHC[ie]->InvertMassHierarchy();
00340           ExtrapFHC[ie]->SetOscPar(OscPar::kTh23,theta23);
00341           ExtrapFHC[ie]->SetSinSq2Th13(ssq2th13);
00342           ExtrapFHC[ie]->OscillatePrediction();
00343         }
00344         for(ie=0;ie<ExtrapRHC.size();ie++)
00345         {
00346           ExtrapRHC[ie]->SetOscPar(OscPar::kTh12,theta12);
00347           ExtrapRHC[ie]->SetOscPar(OscPar::kDeltaM12,dm21);
00348           ExtrapRHC[ie]->SetOscPar(OscPar::kDeltaM23,dm32);
00349           if(!NormalHier) ExtrapRHC[ie]->InvertMassHierarchy();
00350           ExtrapRHC[ie]->SetOscPar(OscPar::kTh23,theta23);
00351           ExtrapRHC[ie]->SetSinSq2Th13(ssq2th13);
00352           ExtrapRHC[ie]->OscillatePrediction();
00353         }
00354         
00355         noff=0;
00356         for(i=0;i<nbins;i++)
00357         {
00358           nobs[i]=0;
00359           
00360           if (i<nbinsFHC){
00361             for(ie=0;ie<ExtrapFHC.size();ie++)
00362               {
00363                 nobs[i] += (ExtrapFHC[ie]->Pred_TotalBkgd_VsBinNumber->GetBinContent(i+1));
00364                 nobs[i] += (ExtrapFHC[ie]->Pred_Signal_VsBinNumber->GetBinContent(i+1));
00365               }
00366           } else if (i>=nbinsFHC){
00367             for(ie=0;ie<ExtrapRHC.size();ie++)
00368               {
00369                 nobs[i] += (ExtrapRHC[ie]->Pred_TotalBkgd_VsBinNumber->GetBinContent(i+1-nbinsFHC));
00370                 nobs[i] += (ExtrapRHC[ie]->Pred_Signal_VsBinNumber->GetBinContent(i+1-nbinsFHC));
00371               }
00372           }
00373           delnexphist[i]->Fill((nobs[i]-nexp[i])/(nexp[i]));
00374         }
00375 
00376         for(i=0;i<nbins;i++)
00377         {
00378           for(k=0;k<nbins;k++)
00379           {
00380             if(k>i)
00381             {
00382               delnexphist_offdiag[noff]->Fill((nobs[i]-nexp[i])*(nobs[k]-nexp[k])/(nexp[i]*nexp[k]));
00383               noff++;
00384             }
00385           }
00386         }
00387       }
00388       
00389       noff=0;
00390       for(i=0;i<nbins;i++)
00391       {
00392         oscparerr[i] = delnexphist[i]->GetRMS();
00393         delnexphist[i]->Write();
00394         for(k=0;k<nbins;k++)
00395         {
00396           if(k>i)
00397           {
00398             oscparerr_offdiag[noff] = delnexphist_offdiag[noff]->GetRMS();
00399             if(delnexphist_offdiag[noff]->GetMean()<0) oscparerr_offdiag[noff] = -1.*oscparerr_offdiag[noff];
00400             delnexphist_offdiag[noff]->Write();
00401             noff++;
00402           }
00403         }
00404       }
00405       f->Close();
00406       
00407       gROOT->cd("/");
00408       
00409       for(i=0;i<nbins;i++)
00410       {
00411         if (i<nbinsFHC){
00412           for(ie=0;ie<ExtrapFHC.size();ie++)
00413             {
00414               ftree2[i][ie]->Fill();
00415             }
00416         } else if (i>=nbinsFHC){
00417           for(ie=0;ie<ExtrapRHC.size();ie++)
00418             {
00419               ftree2[i][ie]->Fill();
00420             }
00421         }
00422         ftree[i]->Fill();
00423       }
00424       
00425       f = new TFile(gSystem->ExpandPathName(s.c_str()),"UPDATE");
00426       
00427       if(l%100==0) cout<<100.*l/((nDeltaSteps+1)*(nSinSq2Th13Steps+1))<<"% complete"<<endl;
00428       l++;
00429     }
00430   }
00431   
00432   f->Close();
00433   
00434   double nPOTNearFHC,nPOTFarFHC;
00435   double nPOTNearRHC,nPOTFarRHC;
00436   
00437   TTree *paramtree = new TTree("paramtree","paramtree");
00438   paramtree->Branch("nearPOTFHC",&nPOTNearFHC,"nearPOTFHC/D");
00439   paramtree->Branch("farPOTFHC",&nPOTFarFHC,"farPOTFHC/D");
00440   paramtree->Branch("nearPOTRHC",&nPOTNearRHC,"nearPOTRHC/D");
00441   paramtree->Branch("farPOTRHC",&nPOTFarRHC,"farPOTRHC/D");
00442   paramtree->Branch("Theta12",&Theta12,"Theta12/D");
00443   paramtree->Branch("Theta23",&Theta23,"Theta23/D");
00444   paramtree->Branch("DeltaMSq23",&dm32,"DeltaMSq23/D");
00445   paramtree->Branch("DeltaMSq12",&DeltaMSq12,"DeltaMSq12/D");
00446   
00447   dm32 = DeltaMSq23;
00448   if(!NormalHier) dm32 = -1.*DeltaMSq23;
00449   
00450   for(ie=0;ie<ExtrapFHC.size();ie++)
00451   {
00452 
00453     nPOTNearRHC=0;
00454     nPOTFarRHC=0;
00455     nPOTNearFHC = ExtrapFHC[ie]->GetNearPOT();
00456     nPOTFarFHC = ExtrapFHC[ie]->GetFarPOT();
00457     paramtree->Fill();
00458   }
00459   for(ie=0;ie<ExtrapRHC.size();ie++)
00460   {
00461     nPOTNearFHC=0;
00462     nPOTFarFHC=0;
00463     nPOTNearRHC=0;
00464     nPOTFarRHC=0;
00465     nPOTNearRHC = ExtrapRHC[ie]->GetNearPOT();
00466     nPOTFarRHC = ExtrapRHC[ie]->GetFarPOT();
00467     paramtree->Fill();
00468   }
00469   
00470   TFile *fout = new TFile(gSystem->ExpandPathName(outFileName.c_str()),"RECREATE");
00471   for(i=0;i<nbins;i++)
00472   {
00473     ftree[i]->Write();
00474     if (i<nbinsFHC){
00475       for(ie=0;ie<ExtrapFHC.size();ie++)
00476         {
00477           ftree2[i][ie]->Write();
00478         } 
00479     } else if (i>=nbinsFHC){
00480       for(ie=0;ie<ExtrapRHC.size();ie++)
00481         {
00482           ftree2[i][ie]->Write();
00483         } 
00484     }
00485   }
00486   paramtree->Write();
00487   fout->Close();
00488   
00489   return;
00490 }


Member Data Documentation

Definition at line 29 of file GridGen_Joint.h.

Referenced by AddExtrapFHC(), RunMultiBin_VaryTheta13(), and RunMultiBinOscParErrs().

Definition at line 28 of file GridGen_Joint.h.

Referenced by AddExtrapRHC(), RunMultiBin_VaryTheta13(), and RunMultiBinOscParErrs().


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

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1