loon: NueMatrixHelper.cxx Source File

00001 #include "NueMatrixHelper.h"
00002 #include "TCut.h"
00003 #include "TGraph.h"
00004 #include "TRandom.h"
00005 #include <vector>
00006 #include <fstream>
00007 #include <iostream>
00008 #include "StandardNtuple/NtpStRecord.h"
00009 #include "MCNtuple/NtpMCTruth.h"
00010 #include "MCNtuple/NtpMCStdHep.h"
00011 #include "MCReweight/NuParent.h"
00012 #include "MCReweight/Zfluk.h"
00013 #include "MCReweight/Zbeam.h"
00014 #include "Registry/Registry.h"
00015 #include "MCReweight/MCEventInfo.h"
00016 #include "MCReweight/MCReweight.h"
00017 #include "MCReweight/NeugenWeightCalculator.h"
00018 //#include "reweight.C"
00019 //#include "mcinfo.C"
00020 
00021 using namespace std;
00022 
00023 NueMatrixHelper::NueMatrixHelper(Int_t nbins,Double_t lowx,Double_t uppx) :
00024   NueExtrapHelper(nbins,lowx,uppx)
00025 {
00026 }
00027 
00028 NueMatrixHelper::NueMatrixHelper(Int_t nbins,Double_t *bins) :
00029   NueExtrapHelper(nbins,bins)
00030 {
00031 }
00032 
00033 void NueMatrixHelper::AddNueSystematic(NueSystematic *nueSys)
00034 {
00035   fMatrixHists[nueSys];
00036   Int_t max_bg_index = 0;
00037   while(strcmp(Background::
00038                AsString(Background::EBackground(max_bg_index)),
00039                "?Unknown?")!=0) {
00040     gDirectory->cd("/");
00041     string mh_name = string(Background::
00042                             AsString(Background::EBackground(max_bg_index)));
00043     mh_name += "_" + string(nueSys->GetName());
00044     MatrixHists *mh = new MatrixHists(mh_name.c_str());
00045     mh->fDirectory->cd();
00046     
00047     //helper histos:
00048     mh->fRecoVsTrueEnergy_ND = 
00049       new TH2D("RecoVsTrueEnergy_ND",
00050                "Reco Vs True Energy (NearDet)",
00051                fNXBins,fXBins,fNXBins,fXBins);
00052     mh->fRecoVsTrueEnergy_ND->Sumw2();
00053 
00054     mh->fRecoVsTrueEnergy_FD = 
00055       new TH2D("RecoVsTrueEnergy_FD",
00056                "Reco Vs True Energy (FarDet)",
00057                fNXBins,fXBins,fNXBins,fXBins);
00058     mh->fRecoVsTrueEnergy_FD->Sumw2();
00059 
00060     mh->fEfficiency_ND = 
00061       new TH1D("Efficiency_ND",
00062                "NuMu CC Selection Efficiency with True Energy (NearDet)",
00063                fNXBins,fXBins);
00064     mh->fEfficiency_ND->Sumw2();
00065 
00066     mh->fEfficiency_FD = 
00067       new TH1D("Efficiency_FD",
00068                "NuMu CC Selection Efficiency with True Energy (FarDet)",
00069                fNXBins,fXBins);
00070     mh->fEfficiency_FD->Sumw2();
00071 
00072     mh->fPurity_ND = 
00073       new TH1D("Purity_ND",
00074                "NuMu CC Selection Purity with True Energy (NearDet)",
00075                fNXBins,fXBins);
00076     mh->fPurity_ND->Sumw2();
00077 
00078     mh->fPurity_FD = 
00079       new TH1D("Purity_FD",
00080                "NuMu CC Selection Purity with True Energy (FarDet)",
00081                fNXBins,fXBins);
00082     mh->fPurity_FD->Sumw2();
00083   
00084     mh->fFDVsNDMatrix = 
00085       new TH2D("FDVsNDMatrix",
00086                "Number of FD Vs ND Events with True Energy",
00087                fNXBins,fXBins,fNXBins,fXBins);
00088     mh->fFDVsNDMatrix->Sumw2();
00089 
00090     mh->fFDVsNDMatrixRW = 
00091       new TH2D("FDVsNDMatrixRW",
00092                "Number of FD Vs ND Events with True Energy (with Near Reweight)",
00093                fNXBins,fXBins,fNXBins,fXBins);
00094     mh->fFDVsNDMatrixRW->Sumw2();
00095 
00096     mh->fFDVsNDMatrixXSec = 
00097       new TH2D("FDVsNDMatrixXSec",
00098                "Number of FD Vs ND Events with True Energy (with XSec)",
00099                fNXBins,fXBins,fNXBins,fXBins);
00100     mh->fFDVsNDMatrixXSec->Sumw2();
00101     
00102     mh->fFDVsNDMatrixXSecRW = 
00103       new TH2D("FDVsNDMatrixXSecRW",
00104                "Number of FD Vs ND Events with True Energy (with XSec + Near Reweight)",
00105                fNXBins,fXBins,fNXBins,fXBins);
00106     mh->fFDVsNDMatrixXSecRW->Sumw2();
00107     
00108     mh->fXSec_CC = 
00109       new TH1D("XSec_CC","NuMu CC XSection with True Energy",
00110                fNXBins,fXBins);
00111     mh->fXSec_CC->Sumw2();
00112     
00113     mh->fFracErrOnPred = 
00114       new TH1D("FracErrOnPred",
00115                "Fractional Error on Energy Spectrum with Reco Energy",
00116                fNXBins,fXBins);
00117     mh->fFracErrOnPred->Sumw2();
00118     
00119     //Check list histos:
00120     mh->fRecoEnergyAllEvents_ND = 
00121       new TH1D("RecoEnergyAllEvents_ND",
00122                "Selected Events with Reco Energy (NearDet)",
00123                fNXBins,fXBins);
00124     mh->fRecoEnergyAllEvents_ND->Sumw2();
00125     
00126     mh->fRecoEnergyCCOnlyEvents_ND = 
00127       new TH1D("RecoEnergyCCOnlyEvents_ND",
00128                "NuMu CC Selected Events with Reco Energy (NearDet)",
00129                fNXBins,fXBins);
00130     mh->fRecoEnergyCCOnlyEvents_ND->Sumw2();
00131     
00132     mh->fTrueEnergyCCOnlyEvents_ND = 
00133       new TH1D("TrueEnergyCCOnlyEvents_ND",
00134                "NuMu CC Selected Events with True Energy (NearDet)",
00135                fNXBins,fXBins);
00136     mh->fTrueEnergyCCOnlyEvents_ND->Sumw2();
00137     
00138     mh->fTrueEnergyTrueCCFidEvents_ND = 
00139       new TH1D("TrueEnergyTrueCCFidEvents_ND",
00140                "True Fid NuMu CC Events with True Energy (NearDet)",
00141                fNXBins,fXBins);
00142     mh->fTrueEnergyTrueCCFidEvents_ND->Sumw2();
00143     
00144     mh->fTrueEnergyNuFlux_ND = 
00145       new TH1D("TrueEnergyNuFlux_ND",
00146                "Neutrino Flux with True Energy (NearDet)",
00147                fNXBins,fXBins);
00148     mh->fTrueEnergyNuFlux_ND->Sumw2();
00149     
00150     mh->fTrueEnergyNuFluxRW_ND = 
00151       new TH1D("TrueEnergyNuFluxRW_ND",
00152                "Neutrino Flux with True Energy (NearDet with Reweighting)",
00153                fNXBins,fXBins);
00154     mh->fTrueEnergyNuFluxRW_ND->Sumw2();
00155     
00156     mh->fTrueEnergyNuFlux_FD = 
00157       new TH1D("TrueEnergyNuFlux_FD",
00158                "Neutrino Flux with True Energy (FarDet)",
00159                fNXBins,fXBins);
00160     mh->fTrueEnergyNuFlux_FD->Sumw2();
00161     
00162     mh->fTrueEnergyNuFluxRW_FD = 
00163       new TH1D("TrueEnergyNuFluxRW_FD",
00164                "Neutrino Flux with True Energy (FarDet with Reweighting)",
00165                fNXBins,fXBins);
00166     mh->fTrueEnergyNuFluxRW_FD->Sumw2();
00167     
00168     mh->fTrueEnergyCCFlux_ND = 
00169       new TH1D("TrueEnergyCCFlux_ND",
00170                "NuMu CC Flux with True Energy (NearDet)",
00171                fNXBins,fXBins);
00172     mh->fTrueEnergyCCFlux_ND->Sumw2();
00173     
00174     mh->fTrueEnergyCCFluxRW_ND = 
00175       new TH1D("TrueEnergyCCFluxRW_ND",
00176                "NuMu CC Flux with True Energy with (NearDet with Reweighting)",
00177                fNXBins,fXBins);
00178     mh->fTrueEnergyCCFluxRW_ND->Sumw2();
00179     
00180     mh->fTrueEnergyCCFlux_FD = 
00181       new TH1D("TrueEnergyCCFlux_FD",
00182                "NuMu CC Flux with True Energy (FarDet)",
00183                fNXBins,fXBins);
00184     mh->fTrueEnergyCCFlux_FD->Sumw2();
00185     
00186     mh->fTrueEnergyCCFluxRW_FD = 
00187       new TH1D("TrueEnergyCCFluxRW_FD",
00188                "NuMu CC Flux with True Energy (FarDet with Reweighing)",
00189                fNXBins,fXBins);
00190     mh->fTrueEnergyCCFluxRW_FD->Sumw2();
00191     
00192     mh->fTrueEnergyTrueCCFidEvents_FD = 
00193       new TH1D("TrueEnergyTrueCCFidEvents_FD",
00194                "True Fid NuMu CC Events with True Energy (FarDet)",
00195                fNXBins,fXBins);
00196     mh->fTrueEnergyTrueCCFidEvents_FD->Sumw2();
00197     
00198     mh->fTrueEnergyCCOnlyEvents_FD = 
00199       new TH1D("TrueEnergyCCOnlyEvents_FD",
00200                "NuMu CC Selected Events with True Energy (FarDet)",
00201                fNXBins,fXBins);
00202     mh->fTrueEnergyCCOnlyEvents_FD->Sumw2();
00203     
00204     mh->fRecoEnergyCCOnlyEvents_FD = 
00205       new TH1D("RecoEnergyCCOnlyEvents_FD",
00206                "NuMu CC Selected Events with Reco Energy (FarDet)",
00207                fNXBins,fXBins);
00208     mh->fRecoEnergyCCOnlyEvents_FD->Sumw2();
00209   
00210     mh->fRecoEnergyAllEvents_FD = 
00211       new TH1D("RecoEnergyAllEvents_FD",
00212              "Selected Events with Reco Energy (FarDet)",
00213                fNXBins,fXBins);
00214     mh->fRecoEnergyAllEvents_FD->Sumw2();
00215     (fMatrixHists[nueSys])[Background::EBackground(max_bg_index)] = mh;
00216     max_bg_index++;
00217   }
00218   gDirectory->cd("/");
00219 }
00220 
00221 NueMatrixHelper::~NueMatrixHelper()
00222 {
00223 }
00224 
00225 void NueMatrixHelper::MakeANANUEPlots(Selection::Selection_t sel)
00226 {
00227   if(!fFarChain || !fNearChain) return;  
00228   fCurSel = sel;
00229 
00230   std::map<NueSystematic*,
00231     std::map<Background::Background_t,MatrixHists*> >::iterator mapBeg = 
00232     fMatrixHists.begin();
00233   std::map<NueSystematic*,
00234     std::map<Background::Background_t,MatrixHists*> >::iterator mapEnd = 
00235     fMatrixHists.end();
00236   if(mapBeg==mapEnd) return;
00237 
00238   //ND:
00239   Int_t ndEntries = fNearChain->GetEntries();
00240   for(int i=0;i<ndEntries;i++){
00241     if(i%100000==0) cout << "Processing event " << i << " / " << ndEntries << endl;
00242     mapBeg = fMatrixHists.begin();
00243     while(mapBeg!=mapEnd) {
00244       fNearChain->GetEvent(i);
00245       Double_t totWeight = 1;
00246       if((mapBeg->first)) totWeight = (mapBeg->first)->UpdateRecord(fRecord,sel);
00247       if(this->PassCuts(sel)) {
00248         Double_t recoEnergy = this->GetNueEnergy(sel);
00249         //fill histos:
00250         std::map<Background::Background_t,MatrixHists*>::iterator mhBeg = mapBeg->second.begin();
00251         std::map<Background::Background_t,MatrixHists*>::iterator mhEnd = mapBeg->second.end();
00252         while(mhBeg!=mhEnd){
00253           mhBeg->second->fRecoEnergyAllEvents_ND->Fill(recoEnergy,totWeight);
00254           mhBeg++;
00255         }
00256         Background::Background_t bg = 
00257           Background::TranslateFromMC(fRecord->mctrue.interactionType,
00258                                       fRecord->mctrue.nuFlavor,
00259                                       fRecord->mctrue.nonOscNuFlavor,
00260                                       fRecord->fluxinfo.tptype);                
00261         (mapBeg->second)[bg]->fRecoEnergyCCOnlyEvents_ND->Fill(recoEnergy,totWeight);
00262         (mapBeg->second)[bg]->fTrueEnergyCCOnlyEvents_ND->Fill(fRecord->mctrue.nuEnergy,totWeight);
00263         (mapBeg->second)[bg]->fPurity_ND->Fill(recoEnergy,totWeight);
00264         (mapBeg->second)[bg]->fEfficiency_ND->Fill(fRecord->mctrue.nuEnergy,totWeight);
00265         (mapBeg->second)[bg]->fRecoVsTrueEnergy_ND->Fill(fRecord->mctrue.nuEnergy,
00266                                                          recoEnergy,totWeight); 
00267         Background::Background_t bg2 = 
00268           Background::TranslateFromMC(fRecord->mctrue.interactionType,
00269                                       fRecord->mctrue.nuFlavor,
00270                                       fRecord->mctrue.nonOscNuFlavor);
00271         if(bg2!=bg) {
00272           (mapBeg->second)[bg2]->fRecoEnergyCCOnlyEvents_ND->Fill(recoEnergy,totWeight);
00273           (mapBeg->second)[bg2]->fTrueEnergyCCOnlyEvents_ND->Fill(fRecord->mctrue.nuEnergy,totWeight);
00274           (mapBeg->second)[bg2]->fPurity_ND->Fill(recoEnergy,totWeight);
00275           (mapBeg->second)[bg2]->fEfficiency_ND->Fill(fRecord->mctrue.nuEnergy,totWeight);
00276           (mapBeg->second)[bg2]->fRecoVsTrueEnergy_ND->Fill(fRecord->mctrue.nuEnergy,
00277                                                             recoEnergy,totWeight);
00278         }
00279       }
00280     }
00281   }
00282   
00283   //FD:
00284   Int_t fdEntries = fFarChain->GetEntries();
00285   for(int i=0;i<fdEntries;i++){
00286     if(i%100000==0) cout << "Processing event " << i << " / " << ndEntries << endl;
00287     mapBeg = fMatrixHists.begin();
00288     while(mapBeg!=mapEnd) {
00289       fFarChain->GetEvent(i);
00290       Double_t totWeight = 1;
00291       if((mapBeg->first)) totWeight = (mapBeg->first)->UpdateRecord(fRecord,sel);
00292       if(this->PassCuts(sel)) {
00293         Double_t recoEnergy = this->GetNueEnergy(sel);
00294         //fill histos:
00295         std::map<Background::Background_t,MatrixHists*>::iterator mhBeg = mapBeg->second.begin();
00296         std::map<Background::Background_t,MatrixHists*>::iterator mhEnd = mapBeg->second.end();
00297         while(mhBeg!=mhEnd){
00298           mhBeg->second->fRecoEnergyAllEvents_FD->Fill(recoEnergy,totWeight);
00299           mhBeg++;
00300         }
00301         Background::Background_t bg = 
00302           Background::TranslateFromMC(fRecord->mctrue.interactionType,
00303                                       fRecord->mctrue.nuFlavor,
00304                                       fRecord->mctrue.nonOscNuFlavor,
00305                                       fRecord->fluxinfo.tptype);                
00306         (mapBeg->second)[bg]->fRecoEnergyCCOnlyEvents_FD->Fill(recoEnergy,totWeight);
00307         (mapBeg->second)[bg]->fTrueEnergyCCOnlyEvents_FD->Fill(fRecord->mctrue.nuEnergy,totWeight);
00308         (mapBeg->second)[bg]->fPurity_FD->Fill(recoEnergy,totWeight);
00309         (mapBeg->second)[bg]->fEfficiency_FD->Fill(fRecord->mctrue.nuEnergy,totWeight);
00310         (mapBeg->second)[bg]->fRecoVsTrueEnergy_FD->Fill(fRecord->mctrue.nuEnergy,
00311                                                          recoEnergy,totWeight); 
00312         Background::Background_t bg2 = 
00313           Background::TranslateFromMC(fRecord->mctrue.interactionType,
00314                                       fRecord->mctrue.nuFlavor,
00315                                       fRecord->mctrue.nonOscNuFlavor);
00316         if(bg2!=bg) {
00317           (mapBeg->second)[bg2]->fRecoEnergyCCOnlyEvents_FD->Fill(recoEnergy,totWeight);
00318           (mapBeg->second)[bg2]->fTrueEnergyCCOnlyEvents_FD->Fill(fRecord->mctrue.nuEnergy,totWeight);
00319           (mapBeg->second)[bg2]->fPurity_FD->Fill(recoEnergy,totWeight);
00320           (mapBeg->second)[bg2]->fEfficiency_FD->Fill(fRecord->mctrue.nuEnergy,totWeight);
00321           (mapBeg->second)[bg2]->fRecoVsTrueEnergy_FD->Fill(fRecord->mctrue.nuEnergy,
00322                                                             recoEnergy,totWeight);
00323         }
00324       }
00325     }
00326   }
00327   
00328   mapBeg = fMatrixHists.begin();
00329   while(mapBeg!=mapEnd) {
00330     std::map<Background::Background_t,MatrixHists*>::iterator mhBeg = mapBeg->second.begin();
00331     std::map<Background::Background_t,MatrixHists*>::iterator mhEnd = mapBeg->second.end();
00332     while(mhBeg!=mhEnd) {
00333       for(int i=1;i<=fNXBins;i++){ //loop over true bins
00334         for(int j=1;j<=fNXBins+1;j++){ //loop over reco bins
00335           if(mhBeg->second->fRecoEnergyCCOnlyEvents_ND->GetBinContent(j)>0 && 
00336              mhBeg->second->fRecoVsTrueEnergy_ND->GetBinContent(i,j)>0) {
00337             Float_t error = ( mhBeg->second->fRecoVsTrueEnergy_ND->GetBinError(i,j) / 
00338                               mhBeg->second->fRecoVsTrueEnergy_ND->GetBinContent(i,j) );
00339             mhBeg->second->fRecoVsTrueEnergy_ND->
00340               SetBinContent(i,j,mhBeg->second->fRecoVsTrueEnergy_ND->GetBinContent(i,j)/
00341                             mhBeg->second->fRecoEnergyCCOnlyEvents_ND->GetBinContent(j));
00342             mhBeg->second->fRecoVsTrueEnergy_ND->
00343               SetBinError(i,j,error*mhBeg->second->fRecoVsTrueEnergy_ND->GetBinContent(i,j));
00344           }
00345           else {
00346             mhBeg->second->fRecoVsTrueEnergy_ND->SetBinContent(i,j,0);
00347             mhBeg->second->fRecoVsTrueEnergy_ND->SetBinError(i,j,0);
00348           }
00349           if(mhBeg->second->fRecoEnergyCCOnlyEvents_FD->GetBinContent(j)>0 && 
00350              mhBeg->second->fRecoVsTrueEnergy_FD->GetBinContent(i,j)>0) {
00351             Float_t error = ( mhBeg->second->fRecoVsTrueEnergy_FD->GetBinError(i,j) / 
00352                               mhBeg->second->fRecoVsTrueEnergy_FD->GetBinContent(i,j) );
00353             mhBeg->second->fRecoVsTrueEnergy_FD->
00354               SetBinContent(i,j,mhBeg->second->fRecoVsTrueEnergy_FD->GetBinContent(i,j)/
00355                             mhBeg->second->fRecoEnergyCCOnlyEvents_FD->GetBinContent(j));
00356             mhBeg->second->fRecoVsTrueEnergy_FD->
00357               SetBinError(i,j,error*mhBeg->second->fRecoVsTrueEnergy_FD->GetBinContent(i,j));
00358           }
00359           else {
00360             mhBeg->second->fRecoVsTrueEnergy_FD->SetBinContent(i,j,0);
00361             mhBeg->second->fRecoVsTrueEnergy_FD->SetBinError(i,j,0);
00362           }
00363         }
00364       }
00365     }
00366   }  
00367 }
00368 
00369 /*
00370 Bool_t NueMatrixHelper::MakeSNTPPlots(TChain *ndChain,TChain *fdChain)
00371 {
00372   cout << "In MakeSNTPPlots" << endl;
00373 
00374   Zbeam beamWeighter;
00375   //skzp1:
00376   //Double_t beamWeighterPars[6] = {-0.441888,0.0303515,-0.0494623,
00377   //                              -2.63948e-07,-2.20597,0.0};
00378   //skzp2:
00379   //  Double_t beamWeighterPars[6] = {0.74087,0.23129,-0.0084731,
00380   //-2.014e-07,-2.9149,0.0};
00381 
00382   Zfluk reWeighter;
00383   reWeighter.UseParameterization("SKZP");
00384   //skzp1:
00385   //Double_t reWeightPars[5] = {0.535219,0.485084,1.00619,1.58987,0.870674};
00386   //skzp2:
00387   //  Double_t reWeightPars[7] = {0.55227,-4.3194,1.0035,0.49904,
00388   //0.92193,1.8081,0.86612};
00389   
00390   MCReweight *mcr = &MCReweight::Instance();
00391   NeugenWeightCalculator *n = new NeugenWeightCalculator();
00392   mcr->AddWeightCalculator(n);
00393   Registry *rwtconfig = new Registry();
00394   rwtconfig->UnLockValues();
00395   rwtconfig->UnLockKeys();
00396   rwtconfig->Set("neugen:config_name","MODBYRS");
00397   rwtconfig->Set("neugen:config_no",3);
00398   rwtconfig->LockValues();
00399   rwtconfig->LockKeys();
00400   mcr->ComputeWeight(0,rwtconfig);
00401   NuParent *np=0;
00402 
00403   if(ndChain && fdChain){
00404     NtpStRecord *strecord = NULL;
00405     ndChain->SetBranchAddress("NtpStRecord",&strecord);
00406     Int_t ndEntries = ndChain->GetEntries();
00407     cout << "Starting ND" << endl;
00408     for(int i=0;i<ndEntries;i++) {
00409       if(i%100000==0) cout << "Processing event " << i << " / " << ndEntries << endl;
00410       ndChain->GetEvent(i);
00411       TClonesArray& mcArray = *(strecord->mc);      
00412       Int_t mcEntries = mcArray.GetEntriesFast();
00413       for(int j=0;j<mcEntries;j++) {
00414         NtpMCTruth *ntpTruth = dynamic_cast<NtpMCTruth *>(mcArray[j]);  
00415         if(ntpTruth->vtxz>=1.0 && ntpTruth->vtxz<=5.0 && 
00416            sqrt(pow(ntpTruth->vtxx-1.4885,2) + 
00417                 pow(ntpTruth->vtxy-0.1397,2))<=1.0) {
00418           if(ntpTruth->iaction==1 && ntpTruth->inu==14) {
00419             Double_t flux_weight = reWeighter.GetWeight(ntpTruth->flux.tptype,
00420                                                         ntpTruth->flux.tpx,
00421                                                         ntpTruth->flux.tpy,
00422                                                         ntpTruth->flux.tpz,reWeightPars);
00423             for(int j=0;j<6;j++)
00424               flux_weight *= beamWeighter.GetWeight(1,2,j+1,beamWeighterPars[j],
00425                                                     ntpTruth->p4neu[3]);
00426             if(!doSKZP) flux_weight = 1.;
00427             MCEventInfo ei;
00428             ei.UseStoredXSec(false); 
00429             ei.nuE=ntpTruth->p4neu[3]; ei.nuPx=ntpTruth->p4neu[0]; 
00430             ei.nuPy=ntpTruth->p4neu[1]; ei.nuPz=ntpTruth->p4neu[2];
00431             ei.tarE=ntpTruth->p4tgt[3]; ei.tarPx=ntpTruth->p4tgt[0]; 
00432             ei.tarPy=ntpTruth->p4tgt[1]; ei.tarPz=ntpTruth->p4tgt[2];
00433             ei.y=ntpTruth->y; ei.x=ntpTruth->x; ei.q2=ntpTruth->q2; ei.w2=ntpTruth->w2;
00434             ei.iaction=ntpTruth->iaction; ei.inu=ntpTruth->inu;
00435             ei.iresonance=ntpTruth->iresonance; 
00436             ei.initial_state=GetInitialState(strecord,ntpTruth);
00437             ei.nucleus=GetNucleus(int(ntpTruth->z),int(ntpTruth->a));
00438             ei.had_fs=GetHadronicFinalState(strecord,ntpTruth);
00439             Double_t xsec_weight = 1.;
00440             if(doMODBYRS3) xsec_weight = mcr->ComputeWeight(&ei,np);
00441 
00442             fTrueEnergyTrueCCFidEvents_ND->Fill(ntpTruth->p4neu[3],
00443                                                 flux_weight*xsec_weight);
00444           }
00445         }
00446       }
00447     }
00448  
00449     fdChain->SetBranchAddress("NtpStRecord",&strecord);
00450     Int_t fdEntries = fdChain->GetEntries();
00451     cout << "Starting FD" << endl;
00452     for(int i=0;i<fdEntries;i++) {
00453       if(i%100000==0) cout << "Processing event " << i << " / " << fdEntries << endl;
00454       fdChain->GetEvent(i);
00455       TClonesArray& mcArray = *(strecord->mc);      
00456       Int_t mcEntries = mcArray.GetEntriesFast();
00457       for(int j=0;j<mcEntries;j++) {
00458         NtpMCTruth *ntpTruth = dynamic_cast<NtpMCTruth *>(mcArray[j]);
00459         if(ntpTruth->vtxz>=0.5 && ntpTruth->vtxz<=28.0 && 
00460            (ntpTruth->vtxz>=16.2 || ntpTruth->vtxz<=14.3) && 
00461            (ntpTruth->vtxx*ntpTruth->vtxx + 
00462             ntpTruth->vtxy*ntpTruth->vtxy)<=14.0) {
00463           if(ntpTruth->iaction==1 && ntpTruth->inu==14) {
00464             Double_t flux_weight = reWeighter.GetWeight(ntpTruth->flux.tptype,
00465                                                         ntpTruth->flux.tpx,
00466                                                         ntpTruth->flux.tpy,
00467                                                         ntpTruth->flux.tpz,reWeightPars);
00468             for(int j=0;j<6;j++)
00469               flux_weight *= beamWeighter.GetWeight(2,2,j+1,beamWeighterPars[j],
00470                                                     ntpTruth->p4neu[3]);
00471             if(!doSKZP) flux_weight = 1.;
00472             MCEventInfo ei;
00473             ei.UseStoredXSec(false); 
00474             ei.nuE=ntpTruth->p4neu[3]; ei.nuPx=ntpTruth->p4neu[0]; 
00475             ei.nuPy=ntpTruth->p4neu[1]; ei.nuPz=ntpTruth->p4neu[2];
00476             ei.tarE=ntpTruth->p4tgt[3]; ei.tarPx=ntpTruth->p4tgt[0]; 
00477             ei.tarPy=ntpTruth->p4tgt[1]; ei.tarPz=ntpTruth->p4tgt[2];
00478             ei.y=ntpTruth->y; ei.x=ntpTruth->x; ei.q2=ntpTruth->q2; ei.w2=ntpTruth->w2;
00479             ei.iaction=ntpTruth->iaction; ei.inu=ntpTruth->inu;
00480             ei.iresonance=ntpTruth->iresonance; 
00481             ei.initial_state=GetInitialState(strecord,ntpTruth);
00482             ei.nucleus=GetNucleus(int(ntpTruth->z),int(ntpTruth->a)); 
00483             ei.had_fs=GetHadronicFinalState(strecord,ntpTruth);
00484             Double_t xsec_weight = 1.;
00485             if(doMODBYRS3) xsec_weight = mcr->ComputeWeight(&ei,np);
00486 
00487             fTrueEnergyTrueCCFidEvents_FD->Fill(ntpTruth->p4neu[3],
00488                                                 flux_weight*xsec_weight);
00489           }
00490         }
00491       }
00492     }
00493     return true;
00494   }
00495   return false;
00496 }
00497 
00498 Bool_t NueMatrixHelper::MakeFLUXPlots(TChain *fluxChain,TGraph *xsec_Graph,
00499                                    TH1F *mikehist)
00500 {
00501   cout << "In MakeFLUXPlots" << endl;
00502 
00503   Zbeam beamWeighter;
00504   //skzp1:
00505   //Double_t beamWeighterPars[6] = {-0.441888,0.0303515,-0.0494623,
00506   //                              -2.63948e-07,-2.20597,0.0};
00507   //skzp2:
00508   //  Double_t beamWeighterPars[6] = {0.74087,0.23129,-0.0084731,
00509   //-2.014e-07,-2.9149,0.0};
00510 
00511   Zfluk reWeighter;
00512   reWeighter.UseParameterization("SKZP");
00513   //skzp1:
00514   //Double_t reWeightPars[5] = {0.535219,0.485084,1.00619,1.58987,0.870674};
00515   //skzp2:
00516   //  Double_t reWeightPars[7] = {0.55227,-4.3194,1.0035,0.49904,
00517   //0.92193,1.8081,0.86612};
00518 
00519   
00520   if(fluxChain) {
00521     Int_t          run;
00522     Int_t          evtno;
00523     Float_t        Ndxdz;
00524     Float_t        Ndydz;
00525     Float_t        Npz;
00526     Float_t        Nenergy;
00527     Float_t        Ndxdznea;
00528     Float_t        Ndydznea;
00529     Float_t        Nenergyn;
00530     Float_t        Nwtnear;
00531     Float_t        Ndxdzfar;
00532     Float_t        Ndydzfar;
00533     Float_t        Nenergyf;
00534     Float_t        Nwtfar;
00535     Int_t          Norig;
00536     Int_t          Ndecay;
00537     Int_t          Ntype;
00538     Float_t        Vx;
00539     Float_t        Vy;
00540     Float_t        Vz;
00541     Float_t        pdpx;
00542     Float_t        pdpy;
00543     Float_t        pdpz;
00544     Float_t        ppdxdz;
00545     Float_t        ppdydz;
00546     Float_t        pppz;
00547     Float_t        ppenergy;
00548     Int_t          ppmedium;
00549     Int_t          ptype;
00550     Float_t        ppvx;
00551     Float_t        ppvy;
00552     Float_t        ppvz;
00553     Float_t        muparpx;
00554     Float_t        muparpy;
00555     Float_t        muparpz;
00556     Float_t        mupare;
00557     Float_t        Necm;
00558     Float_t        Nimpwt;
00559     Float_t        xpoint;
00560     Float_t        ypoint;
00561     Float_t        zpoint;
00562     Float_t        tvx;
00563     Float_t        tvy;
00564     Float_t        tvz;
00565     Float_t        tpx;
00566     Float_t        tpy;
00567     Float_t        tpz;
00568     Int_t          tptype;
00569     Int_t          tgen;
00570     Int_t          tgptype;
00571     Float_t        tgppx;
00572     Float_t        tgppy;
00573     Float_t        tgppz;
00574     Float_t        tprivx;
00575     Float_t        tprivy;
00576     Float_t        tprivz;
00577     Float_t        beamx;
00578     Float_t        beamy;
00579     Float_t        beamz;
00580     Float_t        beampx;
00581     Float_t        beampy;
00582     Float_t        beampz;
00583     fluxChain->SetBranchAddress("run",&run);
00584     fluxChain->SetBranchAddress("evtno",&evtno);
00585     fluxChain->SetBranchAddress("Ndxdz",&Ndxdz);
00586     fluxChain->SetBranchAddress("Ndydz",&Ndydz);
00587     fluxChain->SetBranchAddress("Npz",&Npz);
00588     fluxChain->SetBranchAddress("Nenergy",&Nenergy);
00589     fluxChain->SetBranchAddress("Ndxdznea",&Ndxdznea);
00590     fluxChain->SetBranchAddress("Ndydznea",&Ndydznea);
00591     fluxChain->SetBranchAddress("Nenergyn",&Nenergyn);
00592     fluxChain->SetBranchAddress("Nwtnear",&Nwtnear);
00593     fluxChain->SetBranchAddress("Ndxdzfar",&Ndxdzfar);
00594     fluxChain->SetBranchAddress("Ndydzfar",&Ndydzfar);
00595     fluxChain->SetBranchAddress("Nenergyf",&Nenergyf);
00596     fluxChain->SetBranchAddress("Nwtfar",&Nwtfar);
00597     fluxChain->SetBranchAddress("Norig",&Norig);
00598     fluxChain->SetBranchAddress("Ndecay",&Ndecay);
00599     fluxChain->SetBranchAddress("Ntype",&Ntype);
00600     fluxChain->SetBranchAddress("Vx",&Vx);
00601     fluxChain->SetBranchAddress("Vy",&Vy);
00602     fluxChain->SetBranchAddress("Vz",&Vz);
00603     fluxChain->SetBranchAddress("pdpx",&pdpx);
00604     fluxChain->SetBranchAddress("pdpy",&pdpy);
00605     fluxChain->SetBranchAddress("pdpz",&pdpz);
00606     fluxChain->SetBranchAddress("ppdxdz",&ppdxdz);
00607     fluxChain->SetBranchAddress("ppdydz",&ppdydz);
00608     fluxChain->SetBranchAddress("pppz",&pppz);
00609     fluxChain->SetBranchAddress("ppenergy",&ppenergy);
00610     fluxChain->SetBranchAddress("ppmedium",&ppmedium);
00611     fluxChain->SetBranchAddress("ptype",&ptype);
00612     fluxChain->SetBranchAddress("ppvx",&ppvx);
00613     fluxChain->SetBranchAddress("ppvy",&ppvy);
00614     fluxChain->SetBranchAddress("ppvz",&ppvz);
00615     fluxChain->SetBranchAddress("muparpx",&muparpx);
00616     fluxChain->SetBranchAddress("muparpy",&muparpy);
00617     fluxChain->SetBranchAddress("muparpz",&muparpz);
00618     fluxChain->SetBranchAddress("mupare",&mupare);
00619     fluxChain->SetBranchAddress("Necm",&Necm);
00620     fluxChain->SetBranchAddress("Nimpwt",&Nimpwt);
00621     fluxChain->SetBranchAddress("xpoint",&xpoint);
00622     fluxChain->SetBranchAddress("ypoint",&ypoint);
00623     fluxChain->SetBranchAddress("zpoint",&zpoint);
00624     fluxChain->SetBranchAddress("tvx",&tvx);
00625     fluxChain->SetBranchAddress("tvy",&tvy);
00626     fluxChain->SetBranchAddress("tvz",&tvz);
00627     fluxChain->SetBranchAddress("tpx",&tpx);
00628     fluxChain->SetBranchAddress("tpy",&tpy);
00629     fluxChain->SetBranchAddress("tpz",&tpz);
00630     fluxChain->SetBranchAddress("tptype",&tptype);
00631     fluxChain->SetBranchAddress("tgen",&tgen);
00632     fluxChain->SetBranchAddress("tgptype",&tgptype);
00633     fluxChain->SetBranchAddress("tgppx",&tgppx);
00634     fluxChain->SetBranchAddress("tgppy",&tgppy);
00635     fluxChain->SetBranchAddress("tgppz",&tgppz);
00636     fluxChain->SetBranchAddress("tprivx",&tprivx);
00637     fluxChain->SetBranchAddress("tprivy",&tprivy);
00638     fluxChain->SetBranchAddress("tprivz",&tprivz);
00639     fluxChain->SetBranchAddress("beamx",&beamx);
00640     fluxChain->SetBranchAddress("beamy",&beamy);
00641     fluxChain->SetBranchAddress("beamz",&beamz);
00642     fluxChain->SetBranchAddress("beampx",&beampx);
00643     fluxChain->SetBranchAddress("beampy",&beampy);
00644     fluxChain->SetBranchAddress("beampz",&beampz);
00645     Int_t entries = fluxChain->GetEntries();
00646     for(int i=0;i<entries;i++){
00647       if(i%100000==0) cout << "Processing event " << i 
00648                            << " / " << entries << endl;
00649       fluxChain->GetEvent(i);
00650       if(Ntype!=56) continue;
00651       Double_t flux_weight_near = reWeighter.GetWeight(tptype,tpx,tpy,tpz,
00652                                                        reWeightPars);
00653       Double_t flux_weight_far = reWeighter.GetWeight(tptype,tpx,tpy,tpz,
00654                                                       reWeightPars);
00655 
00656       for(int j=0;j<6;j++) {
00657         flux_weight_near *= beamWeighter.GetWeight(1,2,j+1,beamWeighterPars[j],
00658                                                   Nenergyn);
00659         flux_weight_far *= beamWeighter.GetWeight(2,2,j+1,beamWeighterPars[j],
00660                                                   Nenergyf);
00661       }
00662       if(!doSKZP) {
00663         flux_weight_near = 1.;
00664         flux_weight_far = 1.;
00665       }
00666 
00667       Double_t xsec_nd = double(xsec_Graph->Eval(Nenergyn,0,""));
00668       Double_t xsec_fd = double(xsec_Graph->Eval(Nenergyf,0,""));
00669       if(xsec_nd<=0.0) xsec_nd = 0.0;
00670       if(xsec_fd<=0.0) xsec_fd = 0.0;
00671 
00672       fTrueEnergyNuFlux_ND->Fill(Nenergyn,Nimpwt*Nwtnear*flux_weight_near);
00673       fTrueEnergyNuFlux_FD->Fill(Nenergyf,Nimpwt*Nwtfar*flux_weight_far);      
00674       fFDVsNDMatrix->Fill(Nenergyn,Nenergyf,Nimpwt*Nwtfar*flux_weight_far);
00675 
00676       fTrueEnergyCCFlux_ND->Fill(Nenergyn,Nimpwt*Nwtnear*xsec_nd*flux_weight_near);
00677       fTrueEnergyCCFlux_FD->Fill(Nenergyf,Nimpwt*Nwtfar*xsec_fd*flux_weight_far);
00678       fFDVsNDMatrixXSec->Fill(Nenergyn,Nenergyf,Nimpwt*Nwtfar*xsec_fd*flux_weight_far);
00679 
00680       //do reweighting for front face of ND:
00681       BeamVar beamvar;
00682       beamvar.brun = run;
00683       beamvar.bevtno = evtno;
00684       beamvar.bNdxdz = Ndxdz;
00685       beamvar.bNdydz = Ndydz;
00686       beamvar.bNpz = Npz;
00687       beamvar.bNenergy = Nenergy;
00688       beamvar.bNdxdznea = Ndxdznea;
00689       beamvar.bNdydznea = Ndydznea;
00690       beamvar.bNenergyn = Nenergyn;
00691       beamvar.bNwtnear = Nwtnear;
00692       beamvar.bNdxdzfar = Ndxdzfar;
00693       beamvar.bNdydzfar = Ndydzfar;
00694       beamvar.bNenergyf = Nenergyf;
00695       beamvar.bNwtfar = Nwtfar;
00696       beamvar.bNorig = Norig;
00697       beamvar.bNdecay = Ndecay;
00698       beamvar.bNtype = Ntype;
00699       beamvar.bVx = Vx;
00700       beamvar.bVy = Vy;
00701       beamvar.bVz = Vz;
00702       beamvar.bpdpx = pdpx;
00703       beamvar.bpdpy = pdpy;
00704       beamvar.bpdpz = pdpz;
00705       beamvar.bppdxdz = ppdxdz;
00706       beamvar.bppdydz = ppdydz;
00707       beamvar.bpppz = pppz;
00708       beamvar.bppenergy = ppenergy;
00709       beamvar.bppmedium = ppmedium;
00710       beamvar.bptype = ptype;
00711       beamvar.bppvx = ppvx;
00712       beamvar.bppvy = ppvy;
00713       beamvar.bppvz = ppvz;
00714       beamvar.bmuparpx = muparpx;
00715       beamvar.bmuparpy = muparpy;
00716       beamvar.bmuparpz = muparpz;
00717       beamvar.bmupare = mupare;
00718       beamvar.bNecm = Necm;
00719       beamvar.bNimpwt = Nimpwt;
00720       beamvar.bxpoint = xpoint;
00721       beamvar.bypoint = ypoint;
00722       beamvar.bzpoint = zpoint;
00723       beamvar.btvx = tvx;
00724       beamvar.btvy = tvy;
00725       beamvar.btvz = tvz;
00726       beamvar.btpx = tpx;
00727       beamvar.btpy = tpy;
00728       beamvar.btpz = tpz;
00729       beamvar.btptype = tptype;
00730       beamvar.btgen = tgen;
00731       beamvar.btgptype = tgptype;
00732       beamvar.btgppx = tgppx;
00733       beamvar.btgppy = tgppy;
00734       beamvar.btgppz = tgppz;
00735       beamvar.btprivx = tprivx;
00736       beamvar.btprivy = tprivy;
00737       beamvar.btprivz = tprivz;
00738       beamvar.bbeamx = beamx;
00739       beamvar.bbeamy = beamy;
00740       beamvar.bbeamz = beamz;
00741       beamvar.bbeampx = beampx;
00742       beamvar.bbeampy = beampy;
00743       beamvar.bbeampz = beampz;
00744 
00745       // beam center (0,0) in detector coord system:
00746       // (148.844,13.97)
00747       // z offset is 104000
00748       // beamdydz = -0.0581
00749 
00750       //re-use neutrinos a few times to spread out importance weights...
00751       for(int j=0;j<10;j++){
00752         //Generate random x,y,z in detector coords:
00753         Double_t r = 1000;
00754         Double_t rand_x = 0;
00755         Double_t rand_y = 0;
00756         Double_t rand_z = gRandom->Uniform(100.0,500.0);
00757         while(r>100.) {
00758           rand_x = gRandom->Uniform(-100,100);
00759           rand_y = gRandom->Uniform(-100,100);
00760           r = TMath::Sqrt(rand_x*rand_x + rand_y*rand_y);
00761         }
00762 
00763         //now translate to beam coords:
00764         Double_t X = rand_x;
00765         Double_t Y = (rand_y*TMath::Cos(3.323155*TMath::DegToRad()) + 
00766                       rand_z*TMath::Sin(3.323155*TMath::DegToRad()));
00767         Double_t Z = 104000.+(rand_z*TMath::Cos(3.323155*TMath::DegToRad()) - 
00768                               rand_y*TMath::Sin(3.323155*TMath::DegToRad()));
00769         ArrAy newvals_nd = nuwte(X,Y,Z,beamvar);
00770 
00771         //get new fd weight and energy:
00772         ArrAy newvals_fd = nuwte(0,0,73200000.,beamvar);
00773 
00774         //recalc xsec
00775         xsec_nd = double(xsec_Graph->Eval(newvals_nd.new_ene,0,""));
00776         if(xsec_nd<=0.0) xsec_nd = 0.0; 
00777 
00778         xsec_fd = double(xsec_Graph->Eval(newvals_fd.new_ene,0,""));
00779         if(xsec_fd<=0.0) xsec_fd = 0.0; 
00780 
00781         //recalc flux weight
00782         flux_weight_near = reWeighter.GetWeight(tptype,tpx,tpy,tpz,reWeightPars);
00783         flux_weight_far = reWeighter.GetWeight(tptype,tpx,tpy,tpz,reWeightPars);
00784         for(int j=0;j<6;j++) {
00785           flux_weight_near *= beamWeighter.GetWeight(1,2,j+1,beamWeighterPars[j],
00786                                                      newvals_nd.new_ene);
00787           flux_weight_far *= beamWeighter.GetWeight(2,2,j+1,beamWeighterPars[j],
00788                                                     newvals_fd.new_ene);
00789         }
00790         if(!doSKZP) {
00791           flux_weight_near = 1.;
00792           flux_weight_far = 1.;
00793         }
00794 
00795         //fill matrix
00796         fFDVsNDMatrixRW->Fill(newvals_nd.new_ene,newvals_fd.new_ene,
00797                               Nimpwt*newvals_fd.new_weight*flux_weight_far);
00798         fTrueEnergyNuFluxRW_ND->Fill(newvals_nd.new_ene,
00799                                      Nimpwt*newvals_nd.new_weight*flux_weight_near);
00800         fTrueEnergyNuFluxRW_FD->Fill(newvals_fd.new_ene,
00801                                      Nimpwt*newvals_fd.new_weight*flux_weight_far);
00802         
00803         fFDVsNDMatrixXSecRW->Fill(newvals_nd.new_ene,newvals_fd.new_ene,
00804                                   Nimpwt*newvals_fd.new_weight*xsec_fd*flux_weight_far);
00805         fTrueEnergyCCFluxRW_ND->Fill(newvals_nd.new_ene,
00806                                      Nimpwt*newvals_nd.new_weight*xsec_nd*flux_weight_near);
00807         fTrueEnergyCCFluxRW_FD->Fill(newvals_fd.new_ene,
00808                                      Nimpwt*newvals_fd.new_weight*xsec_fd*flux_weight_far);
00809         
00810       }
00811     }
00812     //re-norm the histo:
00813     for(int i=1;i<=fNBins;i++){
00814       for(int j=1;j<=fNBins;j++){
00815         if( fTrueEnergyNuFlux_ND->GetBinContent(i)>0 && 
00816             fFDVsNDMatrix->GetBinContent(i,j)>0 ) {
00817           Float_t error = TMath::Power(fFDVsNDMatrix->GetBinError(i,j) / 
00818                                        fFDVsNDMatrix->GetBinContent(i,j),2);
00819           error = TMath::Sqrt(error);
00820           fFDVsNDMatrix->SetBinContent(i,j,
00821                                        fFDVsNDMatrix->GetBinContent(i,j)/
00822                                        fTrueEnergyNuFlux_ND->GetBinContent(i));
00823           fFDVsNDMatrix->SetBinError(i,j,error * 
00824                                      fFDVsNDMatrix->GetBinContent(i,j));
00825         }
00826         else {
00827           fFDVsNDMatrix->SetBinContent(i,j,0.);
00828           fFDVsNDMatrix->SetBinError(i,j,0.);
00829         }
00830 
00831         if( fTrueEnergyNuFluxRW_ND->GetBinContent(i)>0 && 
00832             fFDVsNDMatrixRW->GetBinContent(i,j)>0 ) {
00833           Float_t error = TMath::Power(fFDVsNDMatrixRW->GetBinError(i,j) / 
00834                                        fFDVsNDMatrixRW->GetBinContent(i,j),2);
00835           error = TMath::Sqrt(error);
00836           fFDVsNDMatrixRW->SetBinContent(i,j,
00837                                          fFDVsNDMatrixRW->GetBinContent(i,j)/
00838                                          fTrueEnergyNuFluxRW_ND->GetBinContent(i));
00839           fFDVsNDMatrixRW->SetBinError(i,j,error * 
00840                                        fFDVsNDMatrixRW->GetBinContent(i,j));
00841         }
00842         else {
00843           fFDVsNDMatrix->SetBinContent(i,j,0.);
00844           fFDVsNDMatrix->SetBinError(i,j,0.);
00845         }
00846 
00847         if( fTrueEnergyCCFlux_ND->GetBinContent(i)>0 &&
00848             fFDVsNDMatrixXSec->GetBinContent(i,j)>0 ) {
00849           Float_t error = TMath::Power(fFDVsNDMatrixXSec->GetBinError(i,j) / 
00850                                        fFDVsNDMatrixXSec->GetBinContent(i,j),2);
00851           error = TMath::Sqrt(error);
00852           fFDVsNDMatrixXSec->SetBinContent(i,j,
00853                                            fFDVsNDMatrixXSec->GetBinContent(i,j)/
00854                                            fTrueEnergyCCFlux_ND->GetBinContent(i));
00855           fFDVsNDMatrixXSec->SetBinError(i,j,error * 
00856                                          fFDVsNDMatrixXSec->GetBinContent(i,j));
00857         }
00858         else {
00859           fFDVsNDMatrixXSec->SetBinContent(i,j,0.);
00860           fFDVsNDMatrixXSec->SetBinError(i,j,0.);
00861         }
00862 
00863         if( fTrueEnergyCCFluxRW_ND->GetBinContent(i)>0 &&
00864             fFDVsNDMatrixXSecRW->GetBinContent(i,j)>0 ) {
00865           Float_t error = TMath::Power(fFDVsNDMatrixXSecRW->GetBinError(i,j) / 
00866                                        fFDVsNDMatrixXSecRW->GetBinContent(i,j),2);
00867           error = TMath::Sqrt(error);
00868           fFDVsNDMatrixXSecRW->SetBinContent(i,j,
00869                                              fFDVsNDMatrixXSecRW->GetBinContent(i,j)/
00870                                              fTrueEnergyCCFluxRW_ND->GetBinContent(i));
00871           fFDVsNDMatrixXSecRW->SetBinError(i,j,error *
00872                                            fFDVsNDMatrixXSecRW->GetBinContent(i,j));
00873         }
00874         else {
00875           fFDVsNDMatrixXSecRW->SetBinContent(i,j,0.);
00876           fFDVsNDMatrixXSecRW->SetBinError(i,j,0.);
00877         }
00878       }
00879     }
00880     return true;
00881   }
00882   return false;
00883 }
00884 
00885 Bool_t NueMatrixHelper::ReadXSecPlots(TFile *mikeFile) 
00886 {
00887   cout << "In ReadXSecPlots" << endl;
00888   if(mikeFile && mikeFile->IsOpen() && !mikeFile->IsZombie()) {
00889     TH1F *mike_xsec_cc = (TH1F*) mikeFile->Get("h_numu_cc_tot");
00890     if(!mike_xsec_cc) return false;
00891     for(int i=1;i<=fNBins;i++) {
00892       Float_t xval = fXSec_CC->GetBinCenter(i);
00893       Int_t theBin = mike_xsec_cc->FindBin(xval);
00894       Float_t est = 0;
00895       if(TMath::Abs(mike_xsec_cc->GetBinCenter(theBin)-xval)<1e-6) 
00896         est = mike_xsec_cc->GetBinContent(theBin);
00897       else if(theBin==1 || mike_xsec_cc->GetBinCenter(theBin)>xval) {
00898         est = mike_xsec_cc->GetBinContent(theBin);
00899         est += ( mike_xsec_cc->GetBinContent(theBin+1) - 
00900                  mike_xsec_cc->GetBinContent(theBin) ) * 
00901           ( (xval - mike_xsec_cc->GetBinCenter(theBin)) / 
00902             (mike_xsec_cc->GetBinCenter(theBin+1) - 
00903              mike_xsec_cc->GetBinCenter(theBin)) );
00904       }
00905       else if(theBin==mike_xsec_cc->GetNbinsX() || 
00906               mike_xsec_cc->GetBinCenter(theBin)<xval){
00907         est = mike_xsec_cc->GetBinContent(theBin-1);
00908         est += ( mike_xsec_cc->GetBinContent(theBin) - 
00909                  mike_xsec_cc->GetBinContent(theBin-1) ) * 
00910           ( (xval - mike_xsec_cc->GetBinCenter(theBin-1)) / 
00911             (mike_xsec_cc->GetBinCenter(theBin) - 
00912              mike_xsec_cc->GetBinCenter(theBin-1)) );
00913       }
00914       fXSec_CC->SetBinContent(i,est);
00915     }
00916     mikeFile->Close();
00917     return true;
00918   }
00919   return false;
00920 }
00921 */
00922 
00923 void NueMatrixHelper::WriteFile(std::string tag)
00924 {
00925   std::map<NueSystematic*,
00926     std::map<Background::Background_t,MatrixHists*> >::iterator mapBeg = 
00927     fMatrixHists.begin();
00928   std::map<NueSystematic*,
00929     std::map<Background::Background_t,MatrixHists*> >::iterator mapEnd = 
00930     fMatrixHists.end();
00931   if(mapBeg==mapEnd) return;
00932 
00933   std::string filename = "MatrixHelper_" + tag + ".root";
00934   TFile *file = new TFile(filename.c_str(),"RECREATE");
00935   file->cd();
00936 
00937   char selection[256]; 
00938   sprintf(selection,"%s",Selection::AsString(fCurSel));  
00939   TTree *tree = new TTree("matrixtree","matrixtree");
00940   tree->Branch("Selection",selection,"Selection/C");
00941   tree->Branch("nearPOT",&fNearPOT,"nearPOT/D");
00942   tree->Branch("farPOT",&fFarPOT,"farPOT/D");
00943 
00944   while(mapBeg!=mapEnd){
00945     std::map<Background::Background_t,MatrixHists*>::iterator Matbeg = (mapBeg->second).begin();
00946     std::map<Background::Background_t,MatrixHists*>::iterator Matend = (mapBeg->second).end();
00947     while(Matbeg!=Matend) {
00948       Matbeg->second->fDirectory->Write();
00949       Matbeg++;
00950     }
00951     if((mapBeg->first)) (mapBeg->first)->MakeBranches(tree);
00952     tree->Fill();
00953     mapBeg++;
00954   }
00955   tree->Write();
00956   delete file;
00957 }