PANAnalysis Class Reference

#include <PANAnalysis.h>

List of all members.

Public Member Functions
	PANAnalysis (TChain *MCChainFar, TChain *MCChainNear=0, TChain *DataChainFar=0, TChain *DataChainNear=0)
	~PANAnalysis ()
void	SetFileNameTag (std::string)
void	SetHistBookInfo (int, float, float)
void	SetHistBookInfo (int, float, float, int, float, float)
void	SetOscillationFunction (Double_t(*fcn)(Double_t *, Double_t *), Float_t, Float_t, int, double *)
void	SetOscillationFunction (TF1 *, double *)
void	VaryFitParam (int ix, int vx)
Bool_t	RecoDist (string name, string expression, int FN, int entries=-1)
Bool_t	MakeMCVector (int FN, int entries=0)
Bool_t	MakeReweightTree (ReweightLooper *, Int_t, Float_t, Float_t, Int_t, Float_t, Float_t)
void	SetReweightTree (TTree *tree)
Bool_t	Do2DFit (Int_t, Float_t, Float_t, Int_t, Float_t, Float_t)
void	EndJob ()
Protected Member Functions
Int_t *	MakeVarMapNear (Int_t &)
Int_t *	MakeVarMapFar (Int_t &)
Protected Attributes
TChain *	fMCChainNear
TChain *	fDataChainNear
TChain *	fMCChainFar
TChain *	fDataChainFar
int	fNumberOfBinsX
float	fRangeLowerLimitX
float	fRangeUpperLimitX
int	fNumberOfBinsY
float	fRangeLowerLimitY
float	fRangeUpperLimitY
TF1 *	fOscillationFunction
int	fNumOscPars
double *	fOscillationParameters
int *	fVaryX
std::string	fOutFileTag
TH2F *	fChi2Surf
TTree *	fReweightTree
map< string, TH1F * >	fNearHist1D
map< string, TH1F * >	fFarHist1D
map< string, TH2F * >	fNearHist2D
map< string, TH2F * >	fFarHist2D
vector< TH1F * >	fBestFit1D
vector< TH2F * >	fBestFit2D
vector< pan_mcev >	fNearMCEvent
vector< pan_mcev >	fFarMCEvent

---

Constructor & Destructor Documentation

PANAnalysis::PANAnalysis (  TChain *  MCChainFar, TChain *  MCChainNear = 0, TChain *  DataChainFar = 0, TChain *  DataChainNear = 0 )

Definition at line 18 of file PANAnalysis.cxx. References fChi2Surf, fDataChainFar, fDataChainNear, fMCChainFar, fMCChainNear, fNumberOfBinsX, fNumberOfBinsY, fNumOscPars, fOscillationFunction, fOscillationParameters, fOutFileTag, fRangeLowerLimitX, fRangeLowerLimitY, fRangeUpperLimitX, fRangeUpperLimitY, fReweightTree, and fVaryX. { fMCChainFar = MCChainFar; fMCChainNear = MCChainNear; fDataChainFar = DataChainFar; fDataChainNear = DataChainNear; fNumberOfBinsX = 10; fRangeLowerLimitX = 0; fRangeUpperLimitX = 10; fNumberOfBinsY = 10; fRangeLowerLimitY = 0; fRangeUpperLimitY = 10; fOscillationFunction = NULL; fNumOscPars = 0; fOscillationParameters = NULL; fVaryX = new int[2]; fVaryX[0] = 0; fVaryX[1] = 0; fOutFileTag = "test"; fChi2Surf = NULL; fReweightTree = NULL; }

PANAnalysis::~PANAnalysis (   )

Definition at line 45 of file PANAnalysis.cxx. References EndJob(), fBestFit1D, fBestFit2D, fFarHist1D, fFarHist2D, fNearHist1D, and fNearHist2D. { this->EndJob(); cout << "Deleting Everything" << endl; map<string,TH1F*>::iterator beg1D = fFarHist1D.begin(); map<string,TH1F*>::iterator end1D = fFarHist1D.end(); while(beg1D!=end1D) { delete beg1D->second; beg1D++;} beg1D = fNearHist1D.begin(); end1D = fNearHist1D.end(); while(beg1D!=end1D) { delete beg1D->second; beg1D++;} map<string,TH2F*>::iterator beg2D = fFarHist2D.begin(); map<string,TH2F*>::iterator end2D = fFarHist2D.end(); while(beg2D!=end2D) { delete beg2D->second; beg2D++;} beg2D = fNearHist2D.begin(); end2D = fNearHist2D.end(); while(beg2D!=end2D) { delete beg2D->second; beg2D++;} vector<TH1F*>::iterator begBest1D = fBestFit1D.begin(); vector<TH1F*>::iterator endBest1D = fBestFit1D.end(); while(begBest1D!=endBest1D) { delete (*begBest1D); begBest1D++;} vector<TH2F*>::iterator begBest2D = fBestFit2D.begin(); vector<TH2F*>::iterator endBest2D = fBestFit2D.end(); while(begBest2D!=endBest2D) { delete (*begBest2D); begBest2D++;} if(fOscillationFunction) delete fOscillationFunction; if(fVaryX) delete [] fVaryX; if(fChi2Surf) delete fChi2Surf; if(fReweightTree) delete fReweightTree; }

---

Member Function Documentation

Bool_t PANAnalysis::Do2DFit (  Int_t  , Float_t  , Float_t  , Int_t  , Float_t  , Float_t  )

Definition at line 523 of file PANAnalysis.cxx. References fBestFit1D, fBestFit2D, fChi2Surf, fFarHist1D, fFarHist2D, fFarMCEvent, fNearHist1D, fNearHist2D, fNearMCEvent, fOscillationFunction, fReweightTree, fVaryX, MadChi2Calc::GetChi2(), MakeVarMapFar(), and MakeVarMapNear(). { cout << "In PANAnalysis::Do2DFit()" << endl; if(!fOscillationFunction||fFarMCEvent.size()==0|| (fFarHist1D.size()==0&&fFarHist2D.size()==0)) { cout << "**ERROR** Did not find one of the following: " << endl; cout << "Oscillation Function," << " Far Detector Reconstructed Data Histogram," << " Far Detector MC Event Vector" << endl; cout << "Ensure these are set and try again" << endl; return false; } MadChi2Calc Chi2Calc; //chi2 calculator object //Set up chi2 surface: Float_t Par1Width=(Par1Max-Par1Min)/Float_t(Par1Bins); Float_t Par2Width=(Par2Max-Par2Min)/Float_t(Par2Bins); fChi2Surf = new TH2F("Chi2Surf","ChiSquare Surface", Par1Bins,Par1Min,Par1Max, Par2Bins,Par2Min,Par2Max); //for holding MC histos: map<string,TH1F*> TempNearHist1D; map<string,TH1F*> TempFarHist1D; map<string,TH2F*> TempNearHist2D; map<string,TH2F*> TempFarHist2D; //copy dimensions of 1D histos: for(int i=0;i<2;i++){ map<string,TH1F*>::iterator beg1D; map<string,TH1F*>::iterator end1D; if(i==0){ beg1D = fNearHist1D.begin(); end1D = fNearHist1D.end(); } else { beg1D = fFarHist1D.begin(); end1D = fFarHist1D.end(); } while(beg1D!=end1D){ string histname = "TEMP_"; histname.append(beg1D->second->GetName()); TH1F *temp = (TH1F*) beg1D->second->Clone(histname.c_str()); temp->Reset(); if(i==0) TempNearHist1D[beg1D->first] = temp; else TempFarHist1D[beg1D->first] = temp; beg1D++; } } //copy dimensions of 2D histos: for(int i=0;i<2;i++){ map<string,TH2F*>::iterator beg2D; map<string,TH2F*>::iterator end2D; if(i==0){ beg2D = fNearHist2D.begin(); end2D = fNearHist2D.end(); } else { beg2D = fFarHist2D.begin(); end2D = fFarHist2D.end(); } while(beg2D!=end2D){ string histname = "TEMP_"; histname.append(beg2D->second->GetName()); TH2F *temp = (TH2F*) beg2D->second->Clone(histname.c_str()); temp->Reset(); if(i==0) TempNearHist2D[beg2D->first] = temp; else TempFarHist2D[beg2D->first] = temp; beg2D++; } } float ChiSq = 999999; float bestPar1 = 0; float bestPar2 = 0; if(!fReweightTree){ //get iterators over MC events: vector<pan_mcev>::iterator begNearEv = fNearMCEvent.begin(); vector<pan_mcev>::iterator endNearEv = fNearMCEvent.end(); vector<pan_mcev>::iterator begFarEv = fFarMCEvent.begin(); vector<pan_mcev>::iterator endFarEv = fFarMCEvent.end(); //make a map of where each user variable is held in pan_mcev struct Int_t totNearUserVar = 0; Int_t *stringMatchNear = MakeVarMapNear(totNearUserVar); cout << "Total Number of NearDet User Variables = " << totNearUserVar << endl; Int_t totFarUserVar = 0; Int_t *stringMatchFar = MakeVarMapFar(totFarUserVar); cout << "Total Number of FarDet User Variables = " << totFarUserVar << endl; //Start loop over oscillation parameters: for(int i=0;i<Par1Bins;i++){ float par1 = Par1Min + (float(i+1)*Par1Width) - Par1Width/2.; if(i%10==0) cout << "Doing Par1 = " << par1 << endl; for(int j=0;j<Par2Bins;j++){ float par2 = Par2Min + (float(j+1)*Par2Width) - Par2Width/2.; //Set oscillation parameters: fOscillationFunction->SetParameter(fVaryX[0],par1); fOscillationFunction->SetParameter(fVaryX[1],par2); //Reset Temp histos: for(int k=0;k<2;k++){ map<string,TH1F*>::iterator beg1D; map<string,TH1F*>::iterator end1D; if(k==0){ beg1D = TempNearHist1D.begin(); end1D = TempNearHist1D.end(); } else { beg1D = TempFarHist1D.begin(); end1D = TempFarHist1D.end(); } while(beg1D!=end1D) { beg1D->second->Reset(); beg1D++;} } for(int k=0;k<2;k++){ map<string,TH2F*>::iterator beg2D; map<string,TH2F*>::iterator end2D; if(k==0){ beg2D = TempNearHist2D.begin(); end2D = TempNearHist2D.end(); } else { beg2D = TempFarHist2D.begin(); end2D = TempFarHist2D.end(); } while(beg2D!=end2D) { beg2D->second->Reset(); beg2D++;} } //start loop over events: //Near: begNearEv = fNearMCEvent.begin(); fOscillationFunction->SetParameter(2,1.0); //set baseline to 1km while(begNearEv!=endNearEv){ fOscillationFunction->SetParameter(0,begFarEv->flavour); fOscillationFunction->SetParameter(1,begFarEv->flavour_noosc); Double_t oscprob = 1; if(begFarEv->cc_nc==1) oscprob = fOscillationFunction->Eval(begFarEv->nu[3]); Int_t counter = 0; map<string,TH1F*>::iterator beg1D = TempNearHist1D.begin(); map<string,TH1F*>::iterator end1D = TempNearHist1D.end(); map<string,TH2F*>::iterator beg2D = TempNearHist2D.begin(); map<string,TH2F*>::iterator end2D = TempNearHist2D.end(); while(beg1D!=end1D){ beg1D->second->Fill(begNearEv->userVar[stringMatchNear[counter]], oscprob); counter++; beg1D++; } while(beg2D!=end2D){ beg2D->second->Fill(begNearEv->userVar[stringMatchNear[counter]], begNearEv->userVar[stringMatchNear[counter+1]], oscprob); counter+=2; beg2D++; } begNearEv++; } //Far: begFarEv = fFarMCEvent.begin(); fOscillationFunction->SetParameter(2,735.0); //set baseline to 735km while(begFarEv!=endFarEv){ fOscillationFunction->SetParameter(0,begFarEv->flavour); fOscillationFunction->SetParameter(1,begFarEv->flavour_noosc); Double_t oscprob = 1; if(begFarEv->cc_nc==1) oscprob = fOscillationFunction->Eval(begFarEv->nu[3]); Int_t counter = 0; map<string,TH1F*>::iterator beg1D = TempFarHist1D.begin(); map<string,TH1F*>::iterator end1D = TempFarHist1D.end(); map<string,TH2F*>::iterator beg2D = TempFarHist2D.begin(); map<string,TH2F*>::iterator end2D = TempFarHist2D.end(); while(beg1D!=end1D){ beg1D->second->Fill(begFarEv->userVar[stringMatchFar[counter]], oscprob); counter++; beg1D++; } while(beg2D!=end2D){ beg2D->second->Fill(begFarEv->userVar[stringMatchFar[counter]], begFarEv->userVar[stringMatchFar[counter+1]], oscprob); counter+=2; beg2D++; } begFarEv++; } //Now calculate deltachi2 for all histos: Double_t chisquare = 0; TH1F *hist1 = NULL; TH2F *hist2 = NULL; //Near: map<string,TH1F*>::iterator beg1D = fNearHist1D.begin(); map<string,TH1F*>::iterator end1D = fNearHist1D.end(); map<string,TH1F*>::iterator Tempbeg1D = TempNearHist1D.begin(); map<string,TH1F*>::iterator Tempend1D = TempNearHist1D.end(); while(beg1D!=end1D){ chisquare += Chi2Calc.GetChi2(beg1D->second,Tempbeg1D->second,hist1, Tempbeg1D->second->Integral() / beg1D->second->Integral()); beg1D++; Tempbeg1D++; } map<string,TH2F*>::iterator beg2D = fNearHist2D.begin(); map<string,TH2F*>::iterator end2D = fNearHist2D.end(); map<string,TH2F*>::iterator Tempbeg2D = TempNearHist2D.begin(); map<string,TH2F*>::iterator Tempend2D = TempNearHist2D.end(); while(beg2D!=end2D){ chisquare += Chi2Calc.GetChi2(beg2D->second,Tempbeg2D->second,hist2, Tempbeg2D->second->Integral() / beg2D->second->Integral()); beg2D++; Tempbeg2D++; } //Far: beg1D = fFarHist1D.begin(); end1D = fFarHist1D.end(); Tempbeg1D = TempFarHist1D.begin(); Tempend1D = TempFarHist1D.end(); while(beg1D!=end1D){ chisquare += Chi2Calc.GetChi2(beg1D->second,Tempbeg1D->second,hist1, Tempbeg1D->second->Integral() / beg1D->second->Integral()); beg1D++; Tempbeg1D++; } beg2D = fFarHist2D.begin(); end2D = fFarHist2D.end(); Tempbeg2D = TempFarHist2D.begin(); Tempend2D = TempFarHist2D.end(); while(beg2D!=end2D){ chisquare += Chi2Calc.GetChi2(beg2D->second,Tempbeg2D->second,hist2, Tempbeg2D->second->Integral() / beg2D->second->Integral()); beg2D++; Tempbeg2D++; } //Fill chi2surface fChi2Surf->Fill(par1,par2,chisquare); //if this is best fit so far: if(chisquare<ChiSq){ vector<TH1F*>::iterator begBest1D = fBestFit1D.begin(); vector<TH1F*>::iterator endBest1D = fBestFit1D.end(); while(begBest1D!=endBest1D) {delete (*begBest1D); begBest1D++;} fBestFit1D.clear(); vector<TH2F*>::iterator begBest2D = fBestFit2D.begin(); vector<TH2F*>::iterator endBest2D = fBestFit2D.end(); while(begBest2D!=endBest2D) {delete (*begBest2D); begBest2D++;} fBestFit2D.clear(); Tempbeg1D = TempNearHist1D.begin(); Tempend1D = TempNearHist1D.end(); while(Tempbeg1D!=Tempend1D) { TString newname(Tempbeg1D->second->GetName()); newname.ReplaceAll("TEMP_","Best_"); TH1F *hist = (TH1F*) Tempbeg1D->second->Clone(newname.Data()); fBestFit1D.push_back(hist); Tempbeg1D++; } Tempbeg2D = TempNearHist2D.begin(); Tempend2D = TempNearHist2D.end(); while(Tempbeg2D!=Tempend2D) { TString newname(Tempbeg2D->second->GetName()); newname.ReplaceAll("TEMP_","Best_"); TH2F *hist = (TH2F*) Tempbeg2D->second->Clone(newname.Data()); fBestFit2D.push_back(hist); Tempbeg2D++; } Tempbeg1D = TempFarHist1D.begin(); Tempend1D = TempFarHist1D.end(); while(Tempbeg1D!=Tempend1D) { TString newname(Tempbeg1D->second->GetName()); newname.ReplaceAll("TEMP_","Best_"); TH1F *hist = (TH1F*) Tempbeg1D->second->Clone(newname.Data()); fBestFit1D.push_back(hist); Tempbeg1D++; } Tempbeg2D = TempFarHist2D.begin(); Tempend2D = TempFarHist2D.end(); while(Tempbeg2D!=Tempend2D) { TString newname(Tempbeg2D->second->GetName()); newname.ReplaceAll("TEMP_","Best_"); TH2F *hist = (TH2F*) Tempbeg2D->second->Clone(newname.Data()); fBestFit2D.push_back(hist); Tempbeg2D++; } ChiSq = chisquare; bestPar1 = par1; bestPar2 = par2; } } } delete [] stringMatchNear; delete [] stringMatchFar; } else { //Do fits with ReweightTree... //Make Chi2Surf contents large: for(int i=1;i<=fChi2Surf->GetNbinsX();i++){ for(int j=1;j<=fChi2Surf->GetNbinsY();j++) { fChi2Surf->SetBinContent(i,j,999999); } } //Set Branch Addresses: Float_t par1 = 0; fReweightTree->SetBranchAddress("par1",&par1); Float_t par2 = 0; fReweightTree->SetBranchAddress("par2",&par2); map<string,Float_t*> branches; for(int i=0;i<2;i++){ map<string,TH1F*>::iterator beg1D; map<string,TH1F*>::iterator end1D; if(i==0){ beg1D = TempNearHist1D.begin(); end1D = TempNearHist1D.end(); } else { beg1D = TempFarHist1D.begin(); end1D = TempFarHist1D.end(); } while(beg1D!=end1D) { string branchn = beg1D->second->GetName(); Float_t *temp = new Float_t[beg1D->second->GetSize()]; branches[branchn] = temp; TString branchname = beg1D->second->GetName(); branchname.Remove(0,branchname.First("_")+1); fReweightTree->SetBranchAddress(branchname.Data(),temp); //Set correct number of entries in histos based on num MC events: if(i==0) beg1D->second->SetEntries(fNearMCEvent.size()); else beg1D->second->SetEntries(fFarMCEvent.size()); beg1D++; } } for(int i=0;i<2;i++){ map<string,TH2F*>::iterator beg2D; map<string,TH2F*>::iterator end2D; if(i==0){ beg2D = TempNearHist2D.begin(); end2D = TempNearHist2D.end(); } else { beg2D = TempFarHist2D.begin(); end2D = TempFarHist2D.end(); } while(beg2D!=end2D) { string branchn = beg2D->second->GetName(); Float_t *temp = new Float_t[beg2D->second->GetSize()]; branches[branchn] = temp; TString branchname = beg2D->second->GetName(); branchname.Remove(0,branchname.First("_")+1); fReweightTree->SetBranchAddress(branchname.Data(),temp); //Set correct number of entries in histos based on num MC events: if(i==0) beg2D->second->SetEntries(fNearMCEvent.size()); else beg2D->second->SetEntries(fFarMCEvent.size()); beg2D++; } } Int_t bestEntry = 0; //loop over entries in reweight tree: for(int i=0;i<fReweightTree->GetEntries();i++){ fReweightTree->GetEntry(i); //calculate chi2 between temp histos and data histos Double_t chisquare = 0; TH1F *hist1 = NULL; TH2F *hist2 = NULL; //Near: map<string,TH1F*>::iterator beg1D = fNearHist1D.begin(); map<string,TH1F*>::iterator end1D = fNearHist1D.end(); map<string,TH1F*>::iterator Tempbeg1D = TempNearHist1D.begin(); map<string,TH1F*>::iterator Tempend1D = TempNearHist1D.end(); while(beg1D!=end1D){ Tempbeg1D->second->Set(Tempbeg1D->second->GetSize(), branches[Tempbeg1D->second->GetName()]); chisquare += Chi2Calc.GetChi2(beg1D->second,Tempbeg1D->second,hist1, Tempbeg1D->second->Integral() / beg1D->second->Integral()); beg1D++; Tempbeg1D++; } map<string,TH2F*>::iterator beg2D = fNearHist2D.begin(); map<string,TH2F*>::iterator end2D = fNearHist2D.end(); map<string,TH2F*>::iterator Tempbeg2D = TempNearHist2D.begin(); map<string,TH2F*>::iterator Tempend2D = TempNearHist2D.end(); while(beg2D!=end2D){ Tempbeg2D->second->Set(Tempbeg2D->second->GetSize(), branches[Tempbeg2D->second->GetName()]); chisquare += Chi2Calc.GetChi2(beg2D->second,Tempbeg2D->second,hist2, Tempbeg2D->second->Integral() / beg2D->second->Integral()); beg2D++; Tempbeg2D++; } //Far: beg1D = fFarHist1D.begin(); end1D = fFarHist1D.end(); Tempbeg1D = TempFarHist1D.begin(); Tempend1D = TempFarHist1D.end(); while(beg1D!=end1D){ Tempbeg1D->second->Set(Tempbeg1D->second->GetSize(), branches[Tempbeg1D->second->GetName()]); chisquare += Chi2Calc.GetChi2(beg1D->second,Tempbeg1D->second,hist1, Tempbeg1D->second->Integral() / beg1D->second->Integral()); beg1D++; Tempbeg1D++; } beg2D = fFarHist2D.begin(); end2D = fFarHist2D.end(); Tempbeg2D = TempFarHist2D.begin(); Tempend2D = TempFarHist2D.end(); while(beg2D!=end2D){ Tempbeg2D->second->Set(Tempbeg2D->second->GetSize(), branches[Tempbeg2D->second->GetName()]); chisquare += Chi2Calc.GetChi2(beg2D->second,Tempbeg2D->second,hist2, Tempbeg2D->second->Integral() / beg2D->second->Integral()); beg2D++; Tempbeg2D++; } //Fill chi2surface Int_t xbin = fChi2Surf->GetXaxis()->FindBin(par1); Int_t ybin = fChi2Surf->GetYaxis()->FindBin(par2); if(chisquare<fChi2Surf->GetBinContent(xbin,ybin)) fChi2Surf->SetBinContent(xbin,ybin,chisquare); //if this is best fit so far: if(chisquare<ChiSq){ ChiSq = chisquare; bestPar1 = par1; bestPar2 = par2; bestEntry = i; } } fReweightTree->GetEntry(bestEntry); map<string,TH1F*>::iterator Tempbeg1D = TempNearHist1D.begin(); map<string,TH1F*>::iterator Tempend1D = TempNearHist1D.end(); while(Tempbeg1D!=Tempend1D) { Tempbeg1D->second->Set(Tempbeg1D->second->GetSize(), branches[Tempbeg1D->second->GetName()]); TString newname(Tempbeg1D->second->GetName()); newname.ReplaceAll("TEMP_","Best_"); TH1F *hist = (TH1F*) Tempbeg1D->second->Clone(newname.Data()); fBestFit1D.push_back(hist); Tempbeg1D++; } map<string,TH2F*>::iterator Tempbeg2D = TempNearHist2D.begin(); map<string,TH2F*>::iterator Tempend2D = TempNearHist2D.end(); while(Tempbeg2D!=Tempend2D) { Tempbeg2D->second->Set(Tempbeg2D->second->GetSize(), branches[Tempbeg2D->second->GetName()]); TString newname(Tempbeg2D->second->GetName()); newname.ReplaceAll("TEMP_","Best_"); TH2F *hist = (TH2F*) Tempbeg2D->second->Clone(newname.Data()); fBestFit2D.push_back(hist); Tempbeg2D++; } Tempbeg1D = TempFarHist1D.begin(); Tempend1D = TempFarHist1D.end(); while(Tempbeg1D!=Tempend1D) { Tempbeg1D->second->Set(Tempbeg1D->second->GetSize(), branches[Tempbeg1D->second->GetName()]); TString newname(Tempbeg1D->second->GetName()); newname.ReplaceAll("TEMP_","Best_"); TH1F *hist = (TH1F*) Tempbeg1D->second->Clone(newname.Data()); fBestFit1D.push_back(hist); Tempbeg1D++; } Tempbeg2D = TempFarHist2D.begin(); Tempend2D = TempFarHist2D.end(); while(Tempbeg2D!=Tempend2D) { Tempbeg2D->second->Set(Tempbeg2D->second->GetSize(), branches[Tempbeg2D->second->GetName()]); TString newname(Tempbeg2D->second->GetName()); newname.ReplaceAll("TEMP_","Best_"); TH2F *hist = (TH2F*) Tempbeg2D->second->Clone(newname.Data()); fBestFit2D.push_back(hist); Tempbeg2D++; } map<string,Float_t*>::iterator branchIter = branches.begin(); map<string,Float_t*>::iterator branchEnd = branches.end(); while(branchIter!=branchEnd) { delete [] branchIter->second; branchIter++; } branches.clear(); } cout << "Best DeltaChi2 = " << ChiSq << " Best Par1 = " << bestPar1 << " Best Par2 = " << bestPar2 << endl; //delete temp histos: for(int i=0;i<2;i++){ map<string,TH1F*>::iterator beg1D; map<string,TH1F*>::iterator end1D; if(i==0){ beg1D = TempNearHist1D.begin(); end1D = TempNearHist1D.end(); } else { beg1D = TempFarHist1D.begin(); end1D = TempFarHist1D.end(); } while(beg1D!=end1D) { delete beg1D->second; beg1D++;} } TempNearHist1D.clear(); TempFarHist1D.clear(); for(int i=0;i<2;i++){ map<string,TH2F*>::iterator beg2D; map<string,TH2F*>::iterator end2D; if(i==0){ beg2D = TempNearHist2D.begin(); end2D = TempNearHist2D.end(); } else { beg2D = TempFarHist2D.begin(); end2D = TempFarHist2D.end(); } while(beg2D!=end2D) { delete beg2D->second; beg2D++;} } TempNearHist2D.clear(); TempFarHist2D.clear(); return true; }

void PANAnalysis::EndJob (   )

Definition at line 1529 of file PANAnalysis.cxx. References fBestFit1D, fBestFit2D, fChi2Surf, fFarHist1D, fFarHist2D, fNearHist1D, fNearHist2D, and fReweightTree. Referenced by ~PANAnalysis(). { cout << "In PANAnalysis::EndJob()" << endl; TString nam("PANAna_"); nam+=fOutFileTag; nam+=".root"; TFile *fSave = new TFile(nam.Data(),"RECREATE"); fSave->cd(); map<string,TH1F*>::iterator beg1D = fFarHist1D.begin(); map<string,TH1F*>::iterator end1D = fFarHist1D.end(); while(beg1D!=end1D) { beg1D->second->Write(); beg1D++;} beg1D = fNearHist1D.begin(); end1D = fNearHist1D.end(); while(beg1D!=end1D) { beg1D->second->Write(); beg1D++;} map<string,TH2F*>::iterator beg2D = fFarHist2D.begin(); map<string,TH2F*>::iterator end2D = fFarHist2D.end(); while(beg2D!=end2D) { beg2D->second->Write(); beg2D++;} beg2D = fNearHist2D.begin(); end2D = fNearHist2D.end(); while(beg2D!=end2D) { beg2D->second->Write(); beg2D++;} vector<TH1F*>::iterator begBest1D = fBestFit1D.begin(); vector<TH1F*>::iterator endBest1D = fBestFit1D.end(); while(begBest1D!=endBest1D) { (*begBest1D)->Write(); begBest1D++;} vector<TH2F*>::iterator begBest2D = fBestFit2D.begin(); vector<TH2F*>::iterator endBest2D = fBestFit2D.end(); while(begBest2D!=endBest2D) { (*begBest2D)->Write(); begBest2D++;} if(fChi2Surf) fChi2Surf->Write(); if(fReweightTree) fReweightTree->Write(); fSave->Close(); cout << "Closed file" << endl; }

Bool_t PANAnalysis::MakeMCVector (  int  FN, int  entries = 0 )

Definition at line 297 of file PANAnalysis.cxx. References pan_mcev::cc_nc, fFarHist1D, fFarHist2D, fFarMCEvent, pan_mcev::flavour, pan_mcev::flavour_noosc, fNearHist1D, fNearHist2D, fNearMCEvent, NuParent::GetGen(), NuParent::GetPID(), NuParent::GetPx(), NuParent::GetPy(), NuParent::GetPz(), NuParent::GetX(), NuParent::GetY(), NuParent::GetZ(), pan_mcev::had_fs, pan_mcev::init_state, pan_mcev::nu, pan_mcev::nucleus, pan_mcev::pargen, pan_mcev::parmom, pan_mcev::parpid, pan_mcev::parvtx, pan_mcev::process, pan_mcev::q2, pan_mcev::target, pan_mcev::userString, pan_mcev::userVar, pan_mcev::w2, pan_mcev::x, and pan_mcev::y. { cout << "In PANAnalysis::MakeMCVector()" << endl; TChain *chain = NULL; if(FN==1) chain = fMCChainNear; else if(FN==2) chain = fMCChainFar; else return false; Float_t true_enu = 0; chain->SetBranchAddress("true_enu",&true_enu); Float_t true_pxnu = 0; chain->SetBranchAddress("true_pxnu",&true_pxnu); Float_t true_pynu = 0; chain->SetBranchAddress("true_pynu",&true_pynu); Float_t true_pznu = 0; chain->SetBranchAddress("true_pznu",&true_pznu); Float_t true_etgt = 0; chain->SetBranchAddress("true_etgt",&true_etgt); Float_t true_pxtgt = 0; chain->SetBranchAddress("true_pxtgt",&true_pxtgt); Float_t true_pytgt = 0; chain->SetBranchAddress("true_pytgt",&true_pytgt); Float_t true_pztgt = 0; chain->SetBranchAddress("true_pztgt",&true_pztgt); Int_t flavour = 0; chain->SetBranchAddress("flavour",&flavour); Int_t flavour_noosc = 0; chain->SetBranchAddress("nooscflavour",&flavour_noosc); Int_t cc_nc = 0; chain->SetBranchAddress("cc_nc",&cc_nc); Int_t process = 0; chain->SetBranchAddress("process",&process); Int_t init_state = 0; chain->SetBranchAddress("initial_state",&init_state); Int_t nucleus = 0; chain->SetBranchAddress("nucleus",&nucleus); Int_t had_fs = 0; chain->SetBranchAddress("had_fs",&had_fs); Float_t x = 0; chain->SetBranchAddress("true_x",&x); Float_t y = 0; chain->SetBranchAddress("true_y",&y); Float_t q2 = 0; chain->SetBranchAddress("true_q2",&q2); Float_t w2 = 0; chain->SetBranchAddress("true_w2",&w2); NuParent *parent = new NuParent(); chain->SetBranchAddress("NuParent",parent); Int_t pass = 0; chain->SetBranchAddress("pass",&pass); Int_t is_fid = 0; chain->SetBranchAddress("is_fid",&is_fid); Int_t is_mc = 0; chain->SetBranchAddress("is_mc",&is_mc); //Set branch addresses for reco quantities //1D first: map<string,TH1F*>::iterator beg1D; map<string,TH1F*>::iterator end1D; Int_t tot1DExpress = 0; if(FN==1) { beg1D = fNearHist1D.begin(); end1D = fNearHist1D.end(); tot1DExpress = fNearHist1D.size(); } else { beg1D = fFarHist1D.begin(); end1D = fFarHist1D.end(); tot1DExpress = fFarHist1D.size(); } if(tot1DExpress>10) { cout << "uh oh... >10 1D histos not supported yet!" << endl; return false; } TString leaftype_1D[10][10]; Double_t valD_1D[10][10] = {{0}}; Int_t valI_1D[10][10] = {{0}}; Float_t valF_1D[10][10] = {{0}}; TF1 *inputExpression1D[10]; Int_t count_1D[10] = {0}; Int_t histCounter1D = 0; while(beg1D!=end1D){ TString stringExpression1D(beg1D->first); TString funcName = "FUNC_" + stringExpression1D; TLeaf *leaf = 0; TIter leafIter(chain->GetListOfLeaves()); while ((leaf = (TLeaf*) leafIter())) { if(stringExpression1D.Contains(leaf->GetName())) { leaftype_1D[histCounter1D][count_1D[histCounter1D]] = leaf->GetTypeName(); char tempStr[5]; sprintf(tempStr,"[%i]",count_1D[histCounter1D]); stringExpression1D.ReplaceAll(leaf->GetName(),tempStr); if(strcmp(leaf->GetTypeName(),"Int_t")==0) { chain->SetBranchAddress(leaf->GetName(), &valI_1D[histCounter1D][count_1D[histCounter1D]]); count_1D[histCounter1D]++; } else if(strcmp(leaf->GetTypeName(),"Float_t")==0){ chain->SetBranchAddress(leaf->GetName(), &valF_1D[histCounter1D][count_1D[histCounter1D]]); count_1D[histCounter1D]++; } else if(strcmp(leaf->GetTypeName(),"Double_t")==0){ chain->SetBranchAddress(leaf->GetName(), &valD_1D[histCounter1D][count_1D[histCounter1D]]); count_1D[histCounter1D]++; } else { cout << "Can't include " << leaf->GetName() << " in expressions at the moment..." << endl; return false; } } } inputExpression1D[histCounter1D] = new TF1(funcName.Data(),stringExpression1D.Data(),0,1); histCounter1D++; beg1D++; } //now 2D: map<string,TH2F*>::iterator beg2D; map<string,TH2F*>::iterator end2D; Int_t tot2DExpress = 0; if(FN==1) { beg2D = fNearHist2D.begin(); end2D = fNearHist2D.end(); tot2DExpress = fNearHist2D.size(); } else { beg2D = fFarHist2D.begin(); end2D = fFarHist2D.end(); tot2DExpress = fFarHist2D.size(); } if(tot2DExpress>5) { cout << "uh oh... >5 2D histos not supported yet!" << endl; return false; } TString leaftype_2D[10][10]; Double_t valD_2D[10][10] = {{}}; Int_t valI_2D[10][10] = {{}}; Float_t valF_2D[10][10] = {{}}; TF1 *inputExpression2D[10]; Int_t count_2D[10] = {}; Int_t histCounter2D = 0; while(beg2D!=end2D){ for(int i=0;i<2;i++){ TString stringExpression2D(beg2D->first); if(i==0) stringExpression2D.Remove(0,stringExpression2D.First(":")+1); else stringExpression2D.Remove(stringExpression2D.First(":")); TString funcName = "FUNC_" + stringExpression2D; TLeaf *leaf = 0; TIter leafIter(chain->GetListOfLeaves()); while ((leaf = (TLeaf*) leafIter())) { if(stringExpression2D.Contains(leaf->GetName())) { leaftype_2D[histCounter2D][count_2D[histCounter2D]] = leaf->GetTypeName(); char tempStr[5]; sprintf(tempStr,"[%i]",count_2D[histCounter2D]); stringExpression2D.ReplaceAll(leaf->GetName(),tempStr); if(strcmp(leaf->GetTypeName(),"Int_t")==0) { chain->SetBranchAddress(leaf->GetName(), &valI_2D[histCounter2D][count_2D[histCounter2D]]); count_2D[histCounter2D]++; } else if(strcmp(leaf->GetTypeName(),"Float_t")==0){ chain->SetBranchAddress(leaf->GetName(), &valF_2D[histCounter2D][count_2D[histCounter2D]]); count_2D[histCounter2D]++; } else if(strcmp(leaf->GetTypeName(),"Double_t")==0){ chain->SetBranchAddress(leaf->GetName(), &valD_2D[histCounter1D][count_2D[histCounter2D]]); count_2D[histCounter2D]++; } else { cout << "Can't include " << leaf->GetName() << " in expressions at the moment..." << endl; return false; } } } inputExpression2D[histCounter2D] = new TF1(funcName.Data(),stringExpression2D.Data(),0,1); histCounter2D++; } beg2D++; } Int_t totEnt = chain->GetEntries(); for(int i=entries;i<totEnt;i++){ chain->GetEvent(i); if(pass&&is_fid&&is_mc){ pan_mcev mmei; mmei.nu[0] = true_pxnu; mmei.nu[1] = true_pynu; mmei.nu[2] = true_pznu; mmei.nu[3] = true_enu; mmei.target[0] = true_pxtgt; mmei.target[1] = true_pytgt; mmei.target[2] = true_pztgt; mmei.target[3] = true_etgt; mmei.flavour = flavour; mmei.flavour_noosc = flavour_noosc; mmei.cc_nc = cc_nc; mmei.process = process; mmei.init_state = init_state; mmei.nucleus = nucleus; mmei.had_fs = had_fs; mmei.x = x; mmei.y = y; mmei.q2 = q2; mmei.w2 = w2; mmei.parvtx[0] = parent->GetX(); mmei.parvtx[1] = parent->GetY(); mmei.parvtx[2] = parent->GetZ(); mmei.parmom[0] = parent->GetPx(); mmei.parmom[1] = parent->GetPy(); mmei.parmom[2] = parent->GetPz(); mmei.parpid = parent->GetPID(); mmei.pargen = parent->GetGen(); for(int j=0;j<histCounter1D;j++){ for(int k=0;k<count_1D[j];k++){ if(leaftype_1D[j][k].BeginsWith("I")) inputExpression1D[j]->SetParameter(k,valI_1D[j][k]); else if(leaftype_1D[j][k].BeginsWith("F")) inputExpression1D[j]->SetParameter(k,valF_1D[j][k]); else if(leaftype_1D[j][k].BeginsWith("D")) inputExpression1D[j]->SetParameter(k,valD_1D[j][k]); } TString ustring = inputExpression1D[j]->GetName(); ustring.ReplaceAll("FUNC_",""); mmei.userString[j] = ustring; mmei.userVar[j] = inputExpression1D[j]->Eval(0); } for(int j=0;j<histCounter2D;j++){ for(int k=0;k<count_2D[j];k++){ if(leaftype_2D[j][k].BeginsWith("I")) inputExpression2D[j]->SetParameter(k,valI_2D[j][k]); else if(leaftype_2D[j][k].BeginsWith("F")) inputExpression2D[j]->SetParameter(k,valF_2D[j][k]); else if(leaftype_2D[j][k].BeginsWith("D")) inputExpression2D[j]->SetParameter(k,valD_2D[j][k]); } TString ustring = inputExpression2D[j]->GetName(); ustring.ReplaceAll("FUNC_",""); mmei.userString[tot1DExpress+j] = ustring; mmei.userVar[tot1DExpress+j] = inputExpression2D[j]->Eval(0); } if(FN==1) fNearMCEvent.push_back(mmei); else fFarMCEvent.push_back(mmei); } } delete parent; return true; }

Bool_t PANAnalysis::MakeReweightTree (  ReweightLooper *  , Int_t  , Float_t  , Float_t  , Int_t  , Float_t  , Float_t  )

Definition at line 1052 of file PANAnalysis.cxx. References MCReweight::ComputeWeight(), fFarHist1D, fFarHist2D, fFarMCEvent, fNearHist1D, fNearHist2D, fNearMCEvent, fOscillationFunction, fReweightTree, fVaryX, ReweightLooper::GetDaughter(), ReweightLooper::GetName(), ReweightLooper::Grind(), MCEventInfo::had_fs, MCEventInfo::iaction, MCEventInfo::initial_state, MCReweight::Instance(), MCEventInfo::inu, MCEventInfo::iresonance, Registry::LockKeys(), Registry::LockValues(), MakeVarMapFar(), MakeVarMapNear(), MCEventInfo::nucleus, MCEventInfo::nuE, MCEventInfo::nuPx, MCEventInfo::nuPy, MCEventInfo::nuPz, MCEventInfo::q2, ReweightLooper::ResetCounter(), Registry::Set(), NuParent::SetGen(), NuParent::SetPID(), NuParent::SetPx(), NuParent::SetPy(), NuParent::SetPz(), NuParent::SetX(), NuParent::SetY(), NuParent::SetZ(), MCEventInfo::tarE, MCEventInfo::tarPx, MCEventInfo::tarPy, MCEventInfo::tarPz, Registry::UnLockKeys(), Registry::UnLockValues(), MCEventInfo::w2, MCEventInfo::x, and MCEventInfo::y. { cout << "In PANAnalysis::MakeReweightTree()" << endl; Float_t Par1Width=(Par1Max-Par1Min)/Float_t(Par1Bins); Float_t Par2Width=(Par2Max-Par2Min)/Float_t(Par2Bins); fReweightTree = new TTree("ReweightTree","ReweightTree"); //maps of pointers to histograms to write to tree map<string,TH1F*> TempNearHist1D; map<string,TH1F*> TempFarHist1D; map<string,TH2F*> TempNearHist2D; map<string,TH2F*> TempFarHist2D; //copy dimension of 1D histos: for(int i=0;i<2;i++){ map<string,TH1F*>::iterator beg1D; map<string,TH1F*>::iterator end1D; if(i==0){ beg1D = fNearHist1D.begin(); end1D = fNearHist1D.end(); } else { beg1D = fFarHist1D.begin(); end1D = fFarHist1D.end(); } while(beg1D!=end1D){ string histname = "TEMP_"; histname.append(beg1D->second->GetName()); TH1F *temp = (TH1F*) beg1D->second->Clone(histname.c_str()); temp->Reset(); if(i==0) TempNearHist1D[beg1D->first] = temp; else TempFarHist1D[beg1D->first] = temp; char leaflist[256]; sprintf(leaflist,"%s[%i]/F",beg1D->second->GetName(), temp->GetSize()); fReweightTree->Branch(beg1D->second->GetName(), temp->GetArray(),leaflist); beg1D++; } } //copy dimensions of 2D histos: for(int i=0;i<2;i++){ map<string,TH2F*>::iterator beg2D; map<string,TH2F*>::iterator end2D; if(i==0){ beg2D = fNearHist2D.begin(); end2D = fNearHist2D.end(); } else { beg2D = fFarHist2D.begin(); end2D = fFarHist2D.end(); } while(beg2D!=end2D){ string histname = "TEMP_"; histname.append(beg2D->second->GetName()); TH2F *temp = (TH2F*) beg2D->second->Clone(histname.c_str()); temp->Reset(); if(i==0) TempNearHist2D[beg2D->first] = temp; else TempFarHist2D[beg2D->first] = temp; char leaflist[256]; sprintf(leaflist,"%s[%i]/F",beg2D->second->GetName(), temp->GetSize()); fReweightTree->Branch(beg2D->second->GetName(), temp->GetArray(),leaflist); beg2D++; } } //make a map of where each user variable is held in pan_mcev struct: Int_t totNearUserVar = 0; Int_t *stringMatchNear = MakeVarMapNear(totNearUserVar); cout << "Total Number of NearDet User Variables = " << totNearUserVar << endl; Int_t totFarUserVar = 0; Int_t *stringMatchFar = MakeVarMapFar(totFarUserVar); cout << "Total Number of FarDet User Variables = " << totFarUserVar << endl; //add reweight variable values to tree: Registry rwtconfig; map<string,Double_t*> variables; ReweightLooper *daughter = looper; if(daughter) { Double_t *temp = new Double_t[1]; variables[daughter->GetName()] = temp; TString leafname = daughter->GetName(); leafname.ReplaceAll(":","_"); TString leaflist = leafname; leaflist.Append("/D"); fReweightTree->Branch(leafname.Data(),temp,leaflist.Data()); } while((daughter = daughter->GetDaughter())) { Double_t *temp = new Double_t[1]; variables[daughter->GetName()] = temp; TString leafname = daughter->GetName(); leafname.ReplaceAll(":","_"); TString leaflist = leafname; leaflist.Append("/D"); fReweightTree->Branch(leafname.Data(),temp,leaflist.Data()); } //Set branch address for oscillation parameters: float par1 = 0; float par2 = 0; fReweightTree->Branch("par1",&par1,"par1/F"); fReweightTree->Branch("par2",&par2,"par2/F"); //Get instance of MCReweight object: MCReweight &fReweight = MCReweight::Instance(); //get iterators for event vectors vector<pan_mcev>::iterator begNearEv = fNearMCEvent.begin(); vector<pan_mcev>::iterator endNearEv = fNearMCEvent.end(); vector<pan_mcev>::iterator begFarEv = fFarMCEvent.begin(); vector<pan_mcev>::iterator endFarEv = fFarMCEvent.end(); //begin loop over reweight variables: while(looper->Grind(variables)){ //set reweight registry: //fReweight.ResetAllReweightConfigs(); map<string,Double_t*>::iterator begVarIter = variables.begin(); map<string,Double_t*>::iterator endVarIter = variables.end(); rwtconfig.UnLockKeys(); rwtconfig.UnLockValues(); while(begVarIter!=endVarIter){ rwtconfig.Set(begVarIter->first.c_str(),(*begVarIter->second)); begVarIter++; } rwtconfig.LockKeys(); rwtconfig.LockValues(); fReweight.ComputeWeight(0,&rwtconfig); //internally sets params //loop over oscillation parameters for(int i=0;i<Par1Bins;i++){ par1 = Par1Min + (float(i+1)*Par1Width) - Par1Width/2.; if(i%10==0) cout << "Doing Par1 = " << par1 << endl; for(int j=0;j<Par2Bins;j++){ par2 = Par2Min + (float(j+1)*Par2Width) - Par2Width/2.; //Set oscillation parameters: fOscillationFunction->SetParameter(fVaryX[0],par1); fOscillationFunction->SetParameter(fVaryX[1],par2); //Reset Temp histos: for(int k=0;k<2;k++){ map<string,TH1F*>::iterator beg1D; map<string,TH1F*>::iterator end1D; if(k==0){ beg1D = TempNearHist1D.begin(); end1D = TempNearHist1D.end(); } else { beg1D = TempFarHist1D.begin(); end1D = TempFarHist1D.end(); } while(beg1D!=end1D) { beg1D->second->Reset(); beg1D++;} } for(int k=0;k<2;k++){ map<string,TH2F*>::iterator beg2D; map<string,TH2F*>::iterator end2D; if(k==0){ beg2D = TempNearHist2D.begin(); end2D = TempNearHist2D.end(); } else { beg2D = TempFarHist2D.begin(); end2D = TempFarHist2D.end(); } while(beg2D!=end2D) { beg2D->second->Reset(); beg2D++;} } //start loop over events: //Registry event; //event.UnLockKeys(); //event.UnLockValues(); MCEventInfo *eventinfo = new MCEventInfo(); NuParent *parent = new NuParent(); //Near: begNearEv = fNearMCEvent.begin(); fOscillationFunction->SetParameter(2,1.0); //set baseline to 1km while(begNearEv!=endNearEv){ //fill event registry: eventinfo->nuE = begNearEv->nu[3]; eventinfo->nuPx = begNearEv->nu[0]; eventinfo->nuPy = begNearEv->nu[1]; eventinfo->nuPz = begNearEv->nu[2]; eventinfo->tarE = begNearEv->target[3]; eventinfo->tarPx = begNearEv->target[0]; eventinfo->tarPy = begNearEv->target[1]; eventinfo->tarPz = begNearEv->target[2]; eventinfo->x = begNearEv->x; eventinfo->y = begNearEv->y; eventinfo->q2 = begNearEv->q2; eventinfo->w2 = begNearEv->w2; eventinfo->iaction = begNearEv->cc_nc; eventinfo->inu = begNearEv->flavour; eventinfo->iresonance = begNearEv->process; eventinfo->initial_state = begNearEv->init_state; eventinfo->nucleus = begNearEv->nucleus; eventinfo->had_fs = begNearEv->had_fs; parent->SetX(begNearEv->parvtx[0]); parent->SetY(begNearEv->parvtx[1]); parent->SetZ(begNearEv->parvtx[2]); parent->SetPx(begNearEv->parmom[0]); parent->SetPy(begNearEv->parmom[1]); parent->SetPz(begNearEv->parmom[2]); parent->SetPID(begNearEv->parpid); parent->SetGen(begNearEv->pargen); /* event.Set("event:nu_en",begNearEv->nu[3]); event.Set("event:nu_px",begNearEv->nu[0]); event.Set("event:nu_py",begNearEv->nu[1]); event.Set("event:nu_pz",begNearEv->nu[2]); event.Set("event:tar_en",begNearEv->target[3]); event.Set("event:tar_px",begNearEv->target[0]); event.Set("event:tar_py",begNearEv->target[1]); event.Set("event:tar_pz",begNearEv->target[2]); event.Set("event:x",begNearEv->x); event.Set("event:y",begNearEv->y); event.Set("event:q2",begNearEv->q2); event.Set("event:w2",begNearEv->w2); event.Set("event:iaction",begNearEv->cc_nc); event.Set("event:inu",begNearEv->flavour); event.Set("event:process",begNearEv->process); event.Set("event:initial_state",begNearEv->init_state); event.Set("event:nucleus",begNearEv->nucleus); event.Set("event:had_fs",begNearEv->had_fs); event.Set("event:nuparent_x",begNearEv->parvtx[0]); event.Set("event:nuparent_y",begNearEv->parvtx[1]); event.Set("event:nuparent_z",begNearEv->parvtx[2]); event.Set("event:nuparent_px",begNearEv->parmom[0]); event.Set("event:nuparent_py",begNearEv->parmom[1]); event.Set("event:nuparent_pz",begNearEv->parmom[2]); event.Set("event:nuparent_pid",begNearEv->parpid); event.Set("event:nuparent_gen",begNearEv->pargen); */ //set oscillation dependent params: fOscillationFunction->SetParameter(0,begNearEv->flavour); fOscillationFunction->SetParameter(1,begNearEv->flavour_noosc); Double_t oscprob = 1; if(begFarEv->cc_nc==1) oscprob = fOscillationFunction->Eval(begNearEv->nu[3]); Int_t counter = 0; map<string,TH1F*>::iterator beg1D = TempNearHist1D.begin(); map<string,TH1F*>::iterator end1D = TempNearHist1D.end(); map<string,TH2F*>::iterator beg2D = TempNearHist2D.begin(); map<string,TH2F*>::iterator end2D = TempNearHist2D.end(); while(beg1D!=end1D){ beg1D->second->Fill(begNearEv->userVar[stringMatchNear[counter]], oscprob*fReweight.ComputeWeight(eventinfo,parent)); counter++; beg1D++; } while(beg2D!=end2D){ beg2D->second->Fill(begNearEv->userVar[stringMatchNear[counter]], begNearEv->userVar[stringMatchNear[counter+1]], oscprob*fReweight.ComputeWeight(eventinfo,parent)); counter+=2; beg2D++; } begNearEv++; } //Far: begFarEv = fFarMCEvent.begin(); fOscillationFunction->SetParameter(2,735.0); //set baseline to 735km while(begFarEv!=endFarEv){ //fill event registry: eventinfo->nuE = begFarEv->nu[3]; eventinfo->nuPx = begFarEv->nu[0]; eventinfo->nuPy = begFarEv->nu[1]; eventinfo->nuPz = begFarEv->nu[2]; eventinfo->tarE = begFarEv->target[3]; eventinfo->tarPx = begFarEv->target[0]; eventinfo->tarPy = begFarEv->target[1]; eventinfo->tarPz = begFarEv->target[2]; eventinfo->x = begFarEv->x; eventinfo->y = begFarEv->y; eventinfo->q2 = begFarEv->q2; eventinfo->w2 = begFarEv->w2; eventinfo->iaction = begFarEv->cc_nc; eventinfo->inu = begFarEv->flavour; eventinfo->iresonance = begFarEv->process; eventinfo->initial_state = begFarEv->init_state; eventinfo->nucleus = begFarEv->nucleus; eventinfo->had_fs = begFarEv->had_fs; parent->SetX(begFarEv->parvtx[0]); parent->SetY(begFarEv->parvtx[1]); parent->SetZ(begFarEv->parvtx[2]); parent->SetPx(begFarEv->parmom[0]); parent->SetPy(begFarEv->parmom[1]); parent->SetPz(begFarEv->parmom[2]); parent->SetPID(begFarEv->parpid); parent->SetGen(begFarEv->pargen); /* event.Set("event:nu_en",begFarEv->nu[3]); event.Set("event:nu_px",begFarEv->nu[0]); event.Set("event:nu_py",begFarEv->nu[1]); event.Set("event:nu_pz",begFarEv->nu[2]); event.Set("event:tar_en",begFarEv->target[3]); event.Set("event:tar_px",begFarEv->target[0]); event.Set("event:tar_py",begFarEv->target[1]); event.Set("event:tar_pz",begFarEv->target[2]); event.Set("event:x",begFarEv->x); event.Set("event:y",begFarEv->y); event.Set("event:q2",begFarEv->q2); event.Set("event:w2",begFarEv->w2); event.Set("event:iaction",begFarEv->cc_nc); event.Set("event:inu",begFarEv->flavour); event.Set("event:process",begFarEv->process); event.Set("event:initial_state",begFarEv->init_state); event.Set("event:nucleus",begFarEv->nucleus); event.Set("event:had_fs",begFarEv->had_fs); event.Set("event:nuparent_x",begFarEv->parvtx[0]); event.Set("event:nuparent_y",begFarEv->parvtx[1]); event.Set("event:nuparent_z",begFarEv->parvtx[2]); event.Set("event:nuparent_px",begFarEv->parmom[0]); event.Set("event:nuparent_py",begFarEv->parmom[1]); event.Set("event:nuparent_pz",begFarEv->parmom[2]); event.Set("event:nuparent_pid",begFarEv->parpid); event.Set("event:nuparent_gen",begFarEv->pargen); */ //set oscillation dependent params: fOscillationFunction->SetParameter(0,begFarEv->flavour); fOscillationFunction->SetParameter(1,begFarEv->flavour_noosc); Double_t oscprob = 1; if(begFarEv->cc_nc==1) oscprob = fOscillationFunction->Eval(begFarEv->nu[3]); Int_t counter = 0; map<string,TH1F*>::iterator beg1D = TempFarHist1D.begin(); map<string,TH1F*>::iterator end1D = TempFarHist1D.end(); map<string,TH2F*>::iterator beg2D = TempFarHist2D.begin(); map<string,TH2F*>::iterator end2D = TempFarHist2D.end(); while(beg1D!=end1D){ beg1D->second->Fill(begFarEv->userVar[stringMatchFar[counter]], oscprob*fReweight.ComputeWeight(eventinfo,parent)); counter++; beg1D++; } while(beg2D!=end2D){ beg2D->second->Fill(begFarEv->userVar[stringMatchFar[counter]], begFarEv->userVar[stringMatchFar[counter+1]], oscprob*fReweight.ComputeWeight(eventinfo,parent)); counter+=2; beg2D++; } begFarEv++; } delete eventinfo; delete parent; fReweightTree->Fill(); } } } looper->ResetCounter(); //delete temp histos: for(int i=0;i<2;i++){ map<string,TH1F*>::iterator beg1D; map<string,TH1F*>::iterator end1D; if(i==0){ beg1D = TempNearHist1D.begin(); end1D = TempNearHist1D.end(); } else { beg1D = TempFarHist1D.begin(); end1D = TempFarHist1D.end(); } while(beg1D!=end1D) { delete beg1D->second; beg1D++;} } TempNearHist1D.clear(); TempFarHist1D.clear(); for(int i=0;i<2;i++){ map<string,TH2F*>::iterator beg2D; map<string,TH2F*>::iterator end2D; if(i==0){ beg2D = TempNearHist2D.begin(); end2D = TempNearHist2D.end(); } else { beg2D = TempFarHist2D.begin(); end2D = TempFarHist2D.end(); } while(beg2D!=end2D) { delete beg2D->second; beg2D++;} } TempNearHist2D.clear(); TempFarHist2D.clear(); //delete string maps delete [] stringMatchNear; delete [] stringMatchFar; //delete vector map<string,Double_t*>::iterator beg = variables.begin(); map<string,Double_t*>::iterator end = variables.end(); while(beg!=end){ delete [] (beg->second); beg++;} return true; }

Int_t * PANAnalysis::MakeVarMapFar (  Int_t &   )  [protected]

Definition at line 1488 of file PANAnalysis.cxx. References fFarHist1D, fFarHist2D, and fFarMCEvent. Referenced by Do2DFit(), and MakeReweightTree(). { totFarUserVar = fFarHist1D.size() + 2*fFarHist2D.size(); Int_t *stringMatchFar = NULL; vector<pan_mcev>::iterator begFarEv = fFarMCEvent.begin(); vector<pan_mcev>::iterator endFarEv = fFarMCEvent.end(); if(totFarUserVar>0 && begFarEv!=endFarEv){ Int_t counter = 0; stringMatchFar = new Int_t[totFarUserVar]; map<string,TH1F*>::iterator beg1D = fFarHist1D.begin(); map<string,TH1F*>::iterator end1D = fFarHist1D.end(); while(beg1D!=end1D){ for(int i=0;i<totFarUserVar;i++){ if(strcmp((beg1D->first).c_str(),(begFarEv->userString[i]).c_str())==0) { stringMatchFar[counter] = i; break; } } counter++; beg1D++; } map<string,TH2F*>::iterator beg2D = fFarHist2D.begin(); map<string,TH2F*>::iterator end2D = fFarHist2D.end(); while(beg2D!=end2D){ for(int i=0;i<2;i++){ TString histVar = beg2D->first; if(i==0) histVar.Remove(0,histVar.First(":")+1); else histVar.Remove(histVar.First(":")); for(int j=0;j<totFarUserVar;j++){ if(strcmp(histVar.Data(),(begFarEv->userString[j]).c_str())==0) { stringMatchFar[counter] = j; break; } } counter++; } beg2D++; } } return stringMatchFar; }

Int_t * PANAnalysis::MakeVarMapNear (  Int_t &   )  [protected]

Definition at line 1446 of file PANAnalysis.cxx. References fNearHist1D, fNearHist2D, and fNearMCEvent. Referenced by Do2DFit(), and MakeReweightTree(). { totNearUserVar = fNearHist1D.size() + 2*fNearHist2D.size(); Int_t *stringMatchNear = NULL; vector<pan_mcev>::iterator begNearEv = fNearMCEvent.begin(); vector<pan_mcev>::iterator endNearEv = fNearMCEvent.end(); if(totNearUserVar>0 && begNearEv!=endNearEv){ Int_t counter = 0; stringMatchNear = new Int_t[totNearUserVar]; map<string,TH1F*>::iterator beg1D = fNearHist1D.begin(); map<string,TH1F*>::iterator end1D = fNearHist1D.end(); while(beg1D!=end1D){ for(int i=0;i<totNearUserVar;i++){ if(strcmp((beg1D->first).c_str(),(begNearEv->userString[i]).c_str())==0) { stringMatchNear[counter] = i; break; } } counter++; beg1D++; } map<string,TH2F*>::iterator beg2D = fNearHist2D.begin(); map<string,TH2F*>::iterator end2D = fNearHist2D.end(); while(beg2D!=end2D){ for(int i=0;i<2;i++){ TString histVar = beg2D->first; if(i==0) histVar.Remove(0,histVar.First(":")+1); else histVar.Remove(histVar.First(":")); for(int j=0;j<totNearUserVar;j++){ if(strcmp(histVar.Data(),(begNearEv->userString[j]).c_str())==0) { stringMatchNear[counter] = j; break; } } counter++; } beg2D++; } } return stringMatchNear; }

Bool_t PANAnalysis::RecoDist (  string  name, string  expression, int  FN, int  entries = -1 )

Definition at line 122 of file PANAnalysis.cxx. References fDataChainFar, fDataChainNear, fFarHist1D, fFarHist2D, fNearHist1D, fNearHist2D, and fOscillationFunction. { cout << "In PANAnalysis::RecoDist()" << endl; Double_t baseLine = 735.0; if(FN==1) baseLine = 1.0; fOscillationFunction->SetParameter(2,baseLine); Bool_t is2D = false; TString stringExpression1D(expression); if(stringExpression1D.Contains(":")) is2D = true; TH1F *hist1D = NULL; TH2F *hist2D = NULL; if(!is2D) hist1D = new TH1F(name.c_str(),name.c_str(),fNumberOfBinsX, fRangeLowerLimitX,fRangeUpperLimitX); else hist2D = new TH2F(name.c_str(),name.c_str(),fNumberOfBinsX, fRangeLowerLimitX,fRangeUpperLimitX, fNumberOfBinsY,fRangeLowerLimitY,fRangeUpperLimitY); string toHisto = ">>"; if(FN==2&&fDataChainFar){ fDataChainFar->Draw((expression+toHisto+name).c_str(), "pass&&is_fid"); if(!is2D) fFarHist1D[expression] = hist1D; else fFarHist2D[expression] = hist2D; return true; } else if(FN==1&&fDataChainNear){ fDataChainNear->Draw((expression+toHisto+name).c_str(), "pass&&is_fid"); if(!is2D) fNearHist1D[expression] = hist1D; else fNearHist2D[expression] = hist2D; return true; } else { TChain *chain = NULL; if(FN==2&&fMCChainFar) chain = fMCChainFar; else if(FN==1&&fMCChainNear) chain = fMCChainNear; else return false; //Set branch address for true energy Float_t trueEnergy = 0; chain->SetBranchAddress("true_enu",&trueEnergy); Int_t flavour = 0; chain->SetBranchAddress("flavour",&flavour); Int_t flavour_noosc = 0; chain->SetBranchAddress("nooscflavour",&flavour_noosc); Int_t cc_nc = 0; chain->SetBranchAddress("cc_nc",&cc_nc); Int_t pass = 0; chain->SetBranchAddress("pass",&pass); Int_t is_fid = 0; chain->SetBranchAddress("is_fid",&is_fid); Int_t is_mc = 0; chain->SetBranchAddress("is_mc",&is_mc); //Set branch addresses for reco quantities TString stringExpression2D; if(is2D) { stringExpression1D.Remove(0,stringExpression1D.First(":")+1); stringExpression2D = expression; stringExpression2D.Remove(stringExpression2D.First(":")); } TLeaf *leaf = 0; TIter leafIter(chain->GetListOfLeaves()); TString leaftype_1D[10]; Double_t valD_1D[10] = {0}; Int_t valI_1D[10] = {0}; Float_t valF_1D[10] = {0}; TString leaftype_2D[10]; Double_t valD_2D[10] = {0}; Int_t valI_2D[10] = {0}; Float_t valF_2D[10] = {0}; Int_t count_1D = 0; Int_t count_2D = 0; while ((leaf = (TLeaf*) leafIter())) { if(stringExpression1D.Contains(leaf->GetName())) { leaftype_1D[count_1D] = leaf->GetTypeName(); char tempStr[5]; sprintf(tempStr,"[%i]",count_1D); stringExpression1D.ReplaceAll(leaf->GetName(),tempStr); if(leaftype_1D[count_1D]=="Int_t"){ chain->SetBranchAddress(leaf->GetName(),&valI_1D[count_1D]); count_1D++; } else if(leaftype_1D[count_1D]=="Float_t"){ chain->SetBranchAddress(leaf->GetName(),&valF_1D[count_1D]); count_1D++; } else if(leaftype_1D[count_1D]=="Double_t"){ chain->SetBranchAddress(leaf->GetName(),&valD_1D[count_1D]); count_1D++; } else { cout << "Can't include " << leaf->GetName() << " in expressions at the moment..." << endl; return false; } } if(is2D){ if(stringExpression2D.Contains(leaf->GetName())) { leaftype_2D[count_2D] = leaf->GetTypeName(); char tempStr[5]; sprintf(tempStr,"[%i]",count_2D); stringExpression2D.ReplaceAll(leaf->GetName(),tempStr); if(leaftype_2D[count_2D]=="Int_t") { chain->SetBranchAddress(leaf->GetName(),&valI_2D[count_2D]); count_2D++; } else if(leaftype_2D[count_2D]=="Float_t"){ chain->SetBranchAddress(leaf->GetName(),&valF_2D[count_2D]); count_2D++; } else if(leaftype_2D[count_2D]=="Double_t"){ chain->SetBranchAddress(leaf->GetName(),&valD_2D[count_2D]); count_2D++; } else { cout << "Can't include " << leaf->GetName() << " in expressions at the moment..." << endl; return false; } } } } TF1 *inputExpression1D = new TF1("inputExpression1D",stringExpression1D.Data(),0,1); TF1 *inputExpression2D = NULL; if(is2D) inputExpression2D = new TF1("inputExpression2D",stringExpression2D.Data(),0,1); if(entries<0) entries = chain->GetEntries(); for(int i=0;i<entries;i++){ chain->GetEntry(i); if(pass&&is_fid&&is_mc){ //get survival probability: fOscillationFunction->SetParameter(0,flavour); fOscillationFunction->SetParameter(1,flavour_noosc); Double_t oscprob = 1; if(cc_nc==1) oscprob = fOscillationFunction->Eval(trueEnergy); //Fill histos: for(int j=0;j<count_1D;j++){ if(leaftype_1D[j].BeginsWith("I")) inputExpression1D->SetParameter(j,double(valI_1D[j])); else if(leaftype_1D[j].BeginsWith("F")) inputExpression1D->SetParameter(j,double(valF_1D[j])); else if(leaftype_1D[j].BeginsWith("D")) inputExpression1D->SetParameter(j,double(valD_1D[j])); } if(is2D){ for(int j=0;j<count_2D;j++){ if(leaftype_2D[j].BeginsWith("I")) inputExpression2D->SetParameter(j,double(valI_2D[j])); else if(leaftype_2D[j].BeginsWith("F")) inputExpression2D->SetParameter(j,double(valF_2D[j])); else if(leaftype_2D[j].BeginsWith("D")) inputExpression2D->SetParameter(j,double(valD_2D[j])); } hist2D->Fill(inputExpression1D->Eval(0), inputExpression2D->Eval(0),oscprob); } else hist1D->Fill(inputExpression1D->Eval(0),oscprob); } } delete inputExpression1D; if(is2D) delete inputExpression2D; if(FN==2&&fMCChainFar){ if(!is2D) fFarHist1D[expression] = hist1D; else fFarHist2D[expression] = hist2D; return true; } else if(FN==1&&fMCChainNear){ if(!is2D) fNearHist1D[expression] = hist1D; else fNearHist2D[expression] = hist2D; return true; } } return false; }

void PANAnalysis::SetFileNameTag (  std::string   )

Definition at line 76 of file PANAnalysis.cxx. References fOutFileTag. { fOutFileTag = tag; }

void PANAnalysis::SetHistBookInfo (  int  , float  , float  , int  , float  , float  )

Definition at line 88 of file PANAnalysis.cxx. References fNumberOfBinsX, fNumberOfBinsY, fRangeLowerLimitX, fRangeLowerLimitY, fRangeUpperLimitX, and fRangeUpperLimitY. { fNumberOfBinsX = binsX; fRangeLowerLimitX = lowX; fRangeUpperLimitX = upX; fNumberOfBinsY = binsY; fRangeLowerLimitY = lowY; fRangeUpperLimitY = upY; }

void PANAnalysis::SetHistBookInfo (  int  , float  , float  )

Definition at line 81 of file PANAnalysis.cxx. References fNumberOfBinsX, fRangeLowerLimitX, and fRangeUpperLimitX. { fNumberOfBinsX = bins; fRangeLowerLimitX = low; fRangeUpperLimitX = up; }

void PANAnalysis::SetOscillationFunction (  TF1 *  , double *  )

Definition at line 112 of file PANAnalysis.cxx. References fNumOscPars, fOscillationFunction, and fOscillationParameters. { fOscillationParameters = params; fOscillationFunction = new TF1(*form); fNumOscPars = fOscillationFunction->GetNpar(); for(int i=0;i<fNumOscPars;i++){ fOscillationFunction->SetParameter(i,fOscillationParameters[i]); } }

void PANAnalysis::SetOscillationFunction (  Double_t(*)(Double_t *, Double_t *)  fcn, Float_t  , Float_t  , int  , double *  )

Definition at line 99 of file PANAnalysis.cxx. References fNumOscPars, fOscillationFunction, and fOscillationParameters. { fOscillationParameters = params; fOscillationFunction = new TF1("OscFunc",fcn,lowend, highend,numpars); fNumOscPars = numpars; for(int i=0;i<fNumOscPars;i++){ fOscillationFunction->SetParameter(i,fOscillationParameters[i]); } }

void PANAnalysis::SetReweightTree (  TTree *  tree  )  [inline]

Definition at line 101 of file PANAnalysis.h. {fReweightTree = tree;}

void PANAnalysis::VaryFitParam (  int  ix, int  vx )  [inline]

Definition at line 94 of file PANAnalysis.h. {fVaryX[ix] = vx;}

---

Member Data Documentation

vector<TH1F*> PANAnalysis::fBestFit1D [protected]

Definition at line 70 of file PANAnalysis.h. Referenced by Do2DFit(), EndJob(), and ~PANAnalysis().

vector<TH2F*> PANAnalysis::fBestFit2D [protected]

Definition at line 71 of file PANAnalysis.h. Referenced by Do2DFit(), EndJob(), and ~PANAnalysis().

TH2F* PANAnalysis::fChi2Surf [protected]

Definition at line 63 of file PANAnalysis.h. Referenced by Do2DFit(), EndJob(), and PANAnalysis().

TChain* PANAnalysis::fDataChainFar [protected]

Definition at line 47 of file PANAnalysis.h. Referenced by PANAnalysis(), and RecoDist().

TChain* PANAnalysis::fDataChainNear [protected]

Definition at line 45 of file PANAnalysis.h. Referenced by PANAnalysis(), and RecoDist().

map<string,TH1F*> PANAnalysis::fFarHist1D [protected]

Definition at line 67 of file PANAnalysis.h. Referenced by Do2DFit(), EndJob(), MakeMCVector(), MakeReweightTree(), MakeVarMapFar(), RecoDist(), and ~PANAnalysis().

map<string,TH2F*> PANAnalysis::fFarHist2D [protected]

Definition at line 69 of file PANAnalysis.h. Referenced by Do2DFit(), EndJob(), MakeMCVector(), MakeReweightTree(), MakeVarMapFar(), RecoDist(), and ~PANAnalysis().

vector<pan_mcev> PANAnalysis::fFarMCEvent [protected]

Definition at line 74 of file PANAnalysis.h. Referenced by Do2DFit(), MakeMCVector(), MakeReweightTree(), and MakeVarMapFar().

TChain* PANAnalysis::fMCChainFar [protected]

Definition at line 46 of file PANAnalysis.h. Referenced by PANAnalysis().

TChain* PANAnalysis::fMCChainNear [protected]

Definition at line 44 of file PANAnalysis.h. Referenced by PANAnalysis().

map<string,TH1F*> PANAnalysis::fNearHist1D [protected]

Definition at line 66 of file PANAnalysis.h. Referenced by Do2DFit(), EndJob(), MakeMCVector(), MakeReweightTree(), MakeVarMapNear(), RecoDist(), and ~PANAnalysis().

map<string,TH2F*> PANAnalysis::fNearHist2D [protected]

Definition at line 68 of file PANAnalysis.h. Referenced by Do2DFit(), EndJob(), MakeMCVector(), MakeReweightTree(), MakeVarMapNear(), RecoDist(), and ~PANAnalysis().

vector<pan_mcev> PANAnalysis::fNearMCEvent [protected]

Definition at line 73 of file PANAnalysis.h. Referenced by Do2DFit(), MakeMCVector(), MakeReweightTree(), and MakeVarMapNear().

int PANAnalysis::fNumberOfBinsX [protected]

Definition at line 49 of file PANAnalysis.h. Referenced by PANAnalysis(), and SetHistBookInfo().

int PANAnalysis::fNumberOfBinsY [protected]

Definition at line 52 of file PANAnalysis.h. Referenced by PANAnalysis(), and SetHistBookInfo().

int PANAnalysis::fNumOscPars [protected]

Definition at line 57 of file PANAnalysis.h. Referenced by PANAnalysis(), and SetOscillationFunction().

TF1* PANAnalysis::fOscillationFunction [protected]

Definition at line 56 of file PANAnalysis.h. Referenced by Do2DFit(), MakeReweightTree(), PANAnalysis(), RecoDist(), and SetOscillationFunction().

double* PANAnalysis::fOscillationParameters [protected]

Definition at line 58 of file PANAnalysis.h. Referenced by PANAnalysis(), and SetOscillationFunction().

std::string PANAnalysis::fOutFileTag [protected]

Definition at line 61 of file PANAnalysis.h. Referenced by PANAnalysis(), and SetFileNameTag().

float PANAnalysis::fRangeLowerLimitX [protected]

Definition at line 50 of file PANAnalysis.h. Referenced by PANAnalysis(), and SetHistBookInfo().

float PANAnalysis::fRangeLowerLimitY [protected]

Definition at line 53 of file PANAnalysis.h. Referenced by PANAnalysis(), and SetHistBookInfo().

float PANAnalysis::fRangeUpperLimitX [protected]

Definition at line 51 of file PANAnalysis.h. Referenced by PANAnalysis(), and SetHistBookInfo().

float PANAnalysis::fRangeUpperLimitY [protected]

Definition at line 54 of file PANAnalysis.h. Referenced by PANAnalysis(), and SetHistBookInfo().

TTree* PANAnalysis::fReweightTree [protected]

Definition at line 64 of file PANAnalysis.h. Referenced by Do2DFit(), EndJob(), MakeReweightTree(), and PANAnalysis().

int* PANAnalysis::fVaryX [protected]

Definition at line 59 of file PANAnalysis.h. Referenced by Do2DFit(), MakeReweightTree(), and PANAnalysis().

---

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

• PANAnalysis.h
• PANAnalysis.cxx

---