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_mcevfNearMCEvent
vector< pan_mcevfFarMCEvent

Detailed Description

Definition at line 39 of file PANAnalysis.h.


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.

00019                                                                      {
00020   fMCChainFar = MCChainFar;
00021   fMCChainNear = MCChainNear;
00022   fDataChainFar = DataChainFar;
00023   fDataChainNear = DataChainNear;  
00024 
00025   fNumberOfBinsX = 10;
00026   fRangeLowerLimitX = 0;
00027   fRangeUpperLimitX = 10;
00028   fNumberOfBinsY = 10;
00029   fRangeLowerLimitY = 0;
00030   fRangeUpperLimitY = 10;
00031 
00032   fOscillationFunction = NULL;
00033   fNumOscPars = 0;
00034   fOscillationParameters = NULL;  
00035   fVaryX = new int[2];
00036   fVaryX[0] = 0; fVaryX[1] = 0;
00037   
00038   fOutFileTag = "test";
00039   
00040   fChi2Surf = NULL;
00041   fReweightTree = NULL;
00042 
00043 }

PANAnalysis::~PANAnalysis (  ) 

Definition at line 45 of file PANAnalysis.cxx.

References EndJob(), fBestFit1D, fBestFit2D, fChi2Surf, fFarHist1D, fFarHist2D, fNearHist1D, fNearHist2D, fOscillationFunction, fReweightTree, and fVaryX.

00045                           {
00046   this->EndJob();
00047 
00048   cout << "Deleting Everything" << endl;
00049   map<string,TH1F*>::iterator beg1D = fFarHist1D.begin();
00050   map<string,TH1F*>::iterator end1D = fFarHist1D.end();
00051   while(beg1D!=end1D) { delete beg1D->second; beg1D++;}
00052   beg1D = fNearHist1D.begin();
00053   end1D = fNearHist1D.end();
00054   while(beg1D!=end1D) { delete beg1D->second; beg1D++;}
00055 
00056   map<string,TH2F*>::iterator beg2D = fFarHist2D.begin();
00057   map<string,TH2F*>::iterator end2D = fFarHist2D.end();
00058   while(beg2D!=end2D) { delete beg2D->second; beg2D++;}
00059   beg2D = fNearHist2D.begin();
00060   end2D = fNearHist2D.end();
00061   while(beg2D!=end2D) { delete beg2D->second; beg2D++;}
00062 
00063   vector<TH1F*>::iterator begBest1D = fBestFit1D.begin();
00064   vector<TH1F*>::iterator endBest1D = fBestFit1D.end();
00065   while(begBest1D!=endBest1D) { delete (*begBest1D); begBest1D++;}
00066   vector<TH2F*>::iterator begBest2D = fBestFit2D.begin();
00067   vector<TH2F*>::iterator endBest2D = fBestFit2D.end();
00068   while(begBest2D!=endBest2D) { delete (*begBest2D); begBest2D++;}
00069 
00070   if(fOscillationFunction) delete fOscillationFunction;
00071   if(fVaryX) delete [] fVaryX;
00072   if(fChi2Surf) delete fChi2Surf;
00073   if(fReweightTree) delete fReweightTree;
00074 }


Member Function Documentation

Bool_t PANAnalysis::Do2DFit ( Int_t  Par1Bins,
Float_t  Par1Min,
Float_t  Par1Max,
Int_t  Par2Bins,
Float_t  Par2Min,
Float_t  Par2Max 
)

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().

00524                                                                             {
00525 
00526   cout << "In PANAnalysis::Do2DFit()" << endl;
00527   if(!fOscillationFunction||fFarMCEvent.size()==0||
00528      (fFarHist1D.size()==0&&fFarHist2D.size()==0)) {
00529     cout << "**ERROR** Did not find one of the following: " << endl;
00530     cout << "Oscillation Function," 
00531          << " Far Detector Reconstructed Data Histogram," 
00532          << " Far Detector MC Event Vector" << endl;
00533     cout << "Ensure these are set and try again" << endl;
00534     return false;
00535   }
00536   
00537   MadChi2Calc Chi2Calc; //chi2 calculator object
00538   //Set up chi2 surface:
00539   Float_t Par1Width=(Par1Max-Par1Min)/Float_t(Par1Bins);
00540   Float_t Par2Width=(Par2Max-Par2Min)/Float_t(Par2Bins);
00541   fChi2Surf = new TH2F("Chi2Surf","ChiSquare Surface",
00542                        Par1Bins,Par1Min,Par1Max,
00543                        Par2Bins,Par2Min,Par2Max);
00544 
00545   //for holding MC histos:
00546   map<string,TH1F*> TempNearHist1D;
00547   map<string,TH1F*> TempFarHist1D;
00548   map<string,TH2F*> TempNearHist2D;
00549   map<string,TH2F*> TempFarHist2D;
00550 
00551   //copy dimensions of 1D histos:
00552   for(int i=0;i<2;i++){
00553     map<string,TH1F*>::iterator beg1D;    
00554     map<string,TH1F*>::iterator end1D;
00555     if(i==0){
00556       beg1D = fNearHist1D.begin();
00557       end1D = fNearHist1D.end();
00558     }
00559     else {
00560       beg1D = fFarHist1D.begin();
00561       end1D = fFarHist1D.end();
00562     }
00563     while(beg1D!=end1D){
00564       string histname = "TEMP_";
00565       histname.append(beg1D->second->GetName());
00566       TH1F *temp = (TH1F*) beg1D->second->Clone(histname.c_str());
00567       temp->Reset();
00568       if(i==0) TempNearHist1D[beg1D->first] = temp;
00569       else TempFarHist1D[beg1D->first] = temp;
00570       beg1D++;
00571     }
00572   }
00573   //copy dimensions of 2D histos:
00574   for(int i=0;i<2;i++){
00575     map<string,TH2F*>::iterator beg2D;
00576     map<string,TH2F*>::iterator end2D;
00577     if(i==0){
00578       beg2D = fNearHist2D.begin();
00579       end2D = fNearHist2D.end();
00580     }
00581     else {
00582       beg2D = fFarHist2D.begin();
00583       end2D = fFarHist2D.end();
00584     }
00585     while(beg2D!=end2D){
00586       string histname = "TEMP_";
00587       histname.append(beg2D->second->GetName());
00588       TH2F *temp = (TH2F*) beg2D->second->Clone(histname.c_str());
00589       temp->Reset();
00590       if(i==0) TempNearHist2D[beg2D->first] = temp;
00591       else TempFarHist2D[beg2D->first] = temp;
00592       beg2D++;
00593     }
00594   }
00595 
00596   float ChiSq = 999999;
00597   float bestPar1 = 0;
00598   float bestPar2 = 0;
00599 
00600   if(!fReweightTree){
00601     //get iterators over MC events:
00602     vector<pan_mcev>::iterator begNearEv = fNearMCEvent.begin();
00603     vector<pan_mcev>::iterator endNearEv = fNearMCEvent.end();
00604     vector<pan_mcev>::iterator begFarEv = fFarMCEvent.begin();
00605     vector<pan_mcev>::iterator endFarEv = fFarMCEvent.end();
00606 
00607     //make a map of where each user variable is held in pan_mcev struct
00608     Int_t totNearUserVar = 0;
00609     Int_t *stringMatchNear = MakeVarMapNear(totNearUserVar);
00610     cout << "Total Number of NearDet User Variables = " << totNearUserVar << endl;
00611     Int_t totFarUserVar = 0;
00612     Int_t *stringMatchFar = MakeVarMapFar(totFarUserVar);
00613     cout << "Total Number of FarDet User Variables = " << totFarUserVar << endl;
00614 
00615     //Start loop over oscillation parameters:  
00616     for(int i=0;i<Par1Bins;i++){
00617       float par1 = Par1Min + (float(i+1)*Par1Width) - Par1Width/2.;
00618       if(i%10==0) cout << "Doing Par1 = " << par1 << endl;
00619       for(int j=0;j<Par2Bins;j++){
00620         float par2 = Par2Min + (float(j+1)*Par2Width) - Par2Width/2.;
00621         
00622         //Set oscillation parameters:
00623         fOscillationFunction->SetParameter(fVaryX[0],par1);
00624         fOscillationFunction->SetParameter(fVaryX[1],par2);
00625         
00626         //Reset Temp histos:
00627         for(int k=0;k<2;k++){
00628           map<string,TH1F*>::iterator beg1D;    
00629           map<string,TH1F*>::iterator end1D;
00630           if(k==0){
00631             beg1D = TempNearHist1D.begin();
00632             end1D = TempNearHist1D.end();
00633           }
00634           else {
00635             beg1D = TempFarHist1D.begin();
00636             end1D = TempFarHist1D.end();
00637           }
00638           while(beg1D!=end1D) { beg1D->second->Reset(); beg1D++;}
00639         }
00640         for(int k=0;k<2;k++){
00641           map<string,TH2F*>::iterator beg2D;    
00642           map<string,TH2F*>::iterator end2D;
00643           if(k==0){
00644             beg2D = TempNearHist2D.begin();
00645             end2D = TempNearHist2D.end();
00646           }
00647           else {
00648             beg2D = TempFarHist2D.begin();
00649             end2D = TempFarHist2D.end();
00650           }
00651           while(beg2D!=end2D) { beg2D->second->Reset(); beg2D++;}
00652         }
00653         
00654         //start loop over events:
00655         //Near:
00656         begNearEv = fNearMCEvent.begin();
00657         fOscillationFunction->SetParameter(2,1.0); //set baseline to 1km
00658         while(begNearEv!=endNearEv){
00659           fOscillationFunction->SetParameter(0,begFarEv->flavour);
00660           fOscillationFunction->SetParameter(1,begFarEv->flavour_noosc);
00661           Double_t oscprob = 1;
00662           if(begFarEv->cc_nc==1)
00663             oscprob = fOscillationFunction->Eval(begFarEv->nu[3]);
00664           Int_t counter = 0;
00665           map<string,TH1F*>::iterator beg1D = TempNearHist1D.begin();
00666           map<string,TH1F*>::iterator end1D = TempNearHist1D.end();
00667           map<string,TH2F*>::iterator beg2D = TempNearHist2D.begin();
00668           map<string,TH2F*>::iterator end2D = TempNearHist2D.end();
00669           while(beg1D!=end1D){
00670             beg1D->second->Fill(begNearEv->userVar[stringMatchNear[counter]],
00671                                 oscprob);
00672             counter++;
00673             beg1D++;
00674           }
00675           while(beg2D!=end2D){
00676             beg2D->second->Fill(begNearEv->userVar[stringMatchNear[counter]],
00677                                 begNearEv->userVar[stringMatchNear[counter+1]],
00678                                 oscprob);
00679             counter+=2;
00680             beg2D++;
00681           }
00682           begNearEv++;
00683         }
00684         //Far:
00685         begFarEv = fFarMCEvent.begin();
00686         fOscillationFunction->SetParameter(2,735.0); //set baseline to 735km
00687         while(begFarEv!=endFarEv){        
00688           fOscillationFunction->SetParameter(0,begFarEv->flavour);
00689           fOscillationFunction->SetParameter(1,begFarEv->flavour_noosc);
00690           Double_t oscprob = 1;
00691           if(begFarEv->cc_nc==1)
00692             oscprob = fOscillationFunction->Eval(begFarEv->nu[3]);
00693           Int_t counter = 0;
00694           map<string,TH1F*>::iterator beg1D = TempFarHist1D.begin();
00695           map<string,TH1F*>::iterator end1D = TempFarHist1D.end();
00696           map<string,TH2F*>::iterator beg2D = TempFarHist2D.begin();
00697           map<string,TH2F*>::iterator end2D = TempFarHist2D.end();
00698           while(beg1D!=end1D){
00699             beg1D->second->Fill(begFarEv->userVar[stringMatchFar[counter]],
00700                                 oscprob);
00701             counter++;
00702             beg1D++;
00703           }
00704           while(beg2D!=end2D){
00705             beg2D->second->Fill(begFarEv->userVar[stringMatchFar[counter]],
00706                                 begFarEv->userVar[stringMatchFar[counter+1]],
00707                                 oscprob);
00708             counter+=2;
00709             beg2D++;
00710           }
00711           begFarEv++;
00712         }
00713         //Now calculate deltachi2 for all histos:
00714         Double_t chisquare = 0;
00715         TH1F *hist1 = NULL; TH2F *hist2 = NULL;
00716         //Near:
00717         map<string,TH1F*>::iterator beg1D = fNearHist1D.begin();
00718         map<string,TH1F*>::iterator end1D = fNearHist1D.end();
00719         map<string,TH1F*>::iterator Tempbeg1D = TempNearHist1D.begin();
00720         map<string,TH1F*>::iterator Tempend1D = TempNearHist1D.end();
00721         while(beg1D!=end1D){
00722           chisquare += Chi2Calc.GetChi2(beg1D->second,Tempbeg1D->second,hist1,
00723                                         Tempbeg1D->second->Integral() / 
00724                                         beg1D->second->Integral());
00725           beg1D++;
00726           Tempbeg1D++;
00727         }
00728         map<string,TH2F*>::iterator beg2D = fNearHist2D.begin();
00729         map<string,TH2F*>::iterator end2D = fNearHist2D.end();
00730         map<string,TH2F*>::iterator Tempbeg2D = TempNearHist2D.begin();
00731         map<string,TH2F*>::iterator Tempend2D = TempNearHist2D.end();
00732         while(beg2D!=end2D){
00733           chisquare += Chi2Calc.GetChi2(beg2D->second,Tempbeg2D->second,hist2,
00734                                         Tempbeg2D->second->Integral() /
00735                                         beg2D->second->Integral());
00736           beg2D++;
00737           Tempbeg2D++;
00738         }
00739         //Far:
00740         beg1D = fFarHist1D.begin();
00741         end1D = fFarHist1D.end();
00742         Tempbeg1D = TempFarHist1D.begin();
00743         Tempend1D = TempFarHist1D.end();
00744         while(beg1D!=end1D){
00745           chisquare += Chi2Calc.GetChi2(beg1D->second,Tempbeg1D->second,hist1,
00746                                         Tempbeg1D->second->Integral() / 
00747                                         beg1D->second->Integral());
00748           beg1D++;
00749           Tempbeg1D++;
00750         }
00751         beg2D = fFarHist2D.begin();
00752         end2D = fFarHist2D.end();
00753         Tempbeg2D = TempFarHist2D.begin();
00754         Tempend2D = TempFarHist2D.end();
00755         while(beg2D!=end2D){
00756           chisquare += Chi2Calc.GetChi2(beg2D->second,Tempbeg2D->second,hist2,
00757                                         Tempbeg2D->second->Integral() / 
00758                                         beg2D->second->Integral());
00759           beg2D++;
00760           Tempbeg2D++;
00761         }
00762         //Fill chi2surface
00763         fChi2Surf->Fill(par1,par2,chisquare);
00764         
00765         //if this is best fit so far:
00766         if(chisquare<ChiSq){
00767           vector<TH1F*>::iterator begBest1D = fBestFit1D.begin();
00768           vector<TH1F*>::iterator endBest1D = fBestFit1D.end();
00769           while(begBest1D!=endBest1D) {delete (*begBest1D); begBest1D++;}
00770           fBestFit1D.clear();
00771           vector<TH2F*>::iterator begBest2D = fBestFit2D.begin();
00772           vector<TH2F*>::iterator endBest2D = fBestFit2D.end();
00773           while(begBest2D!=endBest2D) {delete (*begBest2D); begBest2D++;}
00774           fBestFit2D.clear();
00775           Tempbeg1D = TempNearHist1D.begin();
00776           Tempend1D = TempNearHist1D.end();
00777           while(Tempbeg1D!=Tempend1D) {
00778             TString newname(Tempbeg1D->second->GetName());
00779             newname.ReplaceAll("TEMP_","Best_");
00780             TH1F *hist = (TH1F*) Tempbeg1D->second->Clone(newname.Data());
00781             fBestFit1D.push_back(hist);
00782             Tempbeg1D++;
00783           }
00784           Tempbeg2D = TempNearHist2D.begin();
00785           Tempend2D = TempNearHist2D.end();
00786           while(Tempbeg2D!=Tempend2D) {
00787             TString newname(Tempbeg2D->second->GetName());
00788             newname.ReplaceAll("TEMP_","Best_");
00789             TH2F *hist = (TH2F*) Tempbeg2D->second->Clone(newname.Data());
00790             fBestFit2D.push_back(hist);
00791             Tempbeg2D++;
00792           }
00793           Tempbeg1D = TempFarHist1D.begin();
00794           Tempend1D = TempFarHist1D.end();
00795           while(Tempbeg1D!=Tempend1D) {
00796             TString newname(Tempbeg1D->second->GetName());
00797             newname.ReplaceAll("TEMP_","Best_");
00798             TH1F *hist = (TH1F*) Tempbeg1D->second->Clone(newname.Data());
00799             fBestFit1D.push_back(hist);
00800             Tempbeg1D++;
00801           }
00802           Tempbeg2D = TempFarHist2D.begin();
00803           Tempend2D = TempFarHist2D.end();
00804           while(Tempbeg2D!=Tempend2D) {
00805             TString newname(Tempbeg2D->second->GetName());
00806             newname.ReplaceAll("TEMP_","Best_");
00807             TH2F *hist = (TH2F*) Tempbeg2D->second->Clone(newname.Data());
00808             fBestFit2D.push_back(hist);
00809             Tempbeg2D++;
00810           }
00811           ChiSq = chisquare;
00812           bestPar1 = par1;
00813           bestPar2 = par2;
00814         }
00815       }
00816     }
00817     delete [] stringMatchNear;
00818     delete [] stringMatchFar;
00819   }
00820   else {
00821     //Do fits with ReweightTree...
00822     //Make Chi2Surf contents large:
00823     for(int i=1;i<=fChi2Surf->GetNbinsX();i++){
00824       for(int j=1;j<=fChi2Surf->GetNbinsY();j++) {
00825         fChi2Surf->SetBinContent(i,j,999999);
00826       }
00827     }
00828     //Set Branch Addresses:
00829     Float_t par1 = 0; fReweightTree->SetBranchAddress("par1",&par1);
00830     Float_t par2 = 0; fReweightTree->SetBranchAddress("par2",&par2);
00831     map<string,Float_t*> branches;
00832     for(int i=0;i<2;i++){
00833       map<string,TH1F*>::iterator beg1D;    
00834       map<string,TH1F*>::iterator end1D;
00835       if(i==0){
00836         beg1D = TempNearHist1D.begin();
00837         end1D = TempNearHist1D.end();
00838       }
00839       else {
00840         beg1D = TempFarHist1D.begin();
00841         end1D = TempFarHist1D.end();
00842       }
00843       while(beg1D!=end1D) {     
00844         string branchn = beg1D->second->GetName();
00845         Float_t *temp = new Float_t[beg1D->second->GetSize()];
00846         branches[branchn] = temp;
00847         TString branchname = beg1D->second->GetName();
00848         branchname.Remove(0,branchname.First("_")+1);
00849         fReweightTree->SetBranchAddress(branchname.Data(),temp);
00850         //Set correct number of entries in histos based on num MC events:
00851         if(i==0) beg1D->second->SetEntries(fNearMCEvent.size());
00852         else beg1D->second->SetEntries(fFarMCEvent.size());
00853         beg1D++;
00854       }
00855     }
00856     for(int i=0;i<2;i++){
00857       map<string,TH2F*>::iterator beg2D;    
00858       map<string,TH2F*>::iterator end2D;
00859       if(i==0){
00860         beg2D = TempNearHist2D.begin();
00861         end2D = TempNearHist2D.end();
00862       }
00863       else {
00864         beg2D = TempFarHist2D.begin();
00865         end2D = TempFarHist2D.end();
00866       }
00867       while(beg2D!=end2D) {
00868         string branchn = beg2D->second->GetName();
00869         Float_t *temp = new Float_t[beg2D->second->GetSize()];
00870         branches[branchn] = temp;
00871         TString branchname = beg2D->second->GetName();
00872         branchname.Remove(0,branchname.First("_")+1);
00873         fReweightTree->SetBranchAddress(branchname.Data(),temp);
00874         //Set correct number of entries in histos based on num MC events:
00875         if(i==0) beg2D->second->SetEntries(fNearMCEvent.size());
00876         else beg2D->second->SetEntries(fFarMCEvent.size());
00877         beg2D++;
00878       }
00879     }
00880 
00881     Int_t bestEntry = 0;
00882     //loop over entries in reweight tree:
00883     for(int i=0;i<fReweightTree->GetEntries();i++){
00884       fReweightTree->GetEntry(i);
00885       
00886       //calculate chi2 between temp histos and data histos
00887       Double_t chisquare = 0;
00888       TH1F *hist1 = NULL; TH2F *hist2 = NULL;
00889       //Near:
00890       map<string,TH1F*>::iterator beg1D = fNearHist1D.begin();
00891       map<string,TH1F*>::iterator end1D = fNearHist1D.end();
00892       map<string,TH1F*>::iterator Tempbeg1D = TempNearHist1D.begin();
00893       map<string,TH1F*>::iterator Tempend1D = TempNearHist1D.end();
00894       while(beg1D!=end1D){
00895         Tempbeg1D->second->Set(Tempbeg1D->second->GetSize(),
00896                                branches[Tempbeg1D->second->GetName()]);
00897         chisquare += Chi2Calc.GetChi2(beg1D->second,Tempbeg1D->second,hist1,
00898                                       Tempbeg1D->second->Integral() / 
00899                                       beg1D->second->Integral());
00900         beg1D++;
00901         Tempbeg1D++;
00902       }
00903       map<string,TH2F*>::iterator beg2D = fNearHist2D.begin();
00904       map<string,TH2F*>::iterator end2D = fNearHist2D.end();
00905       map<string,TH2F*>::iterator Tempbeg2D = TempNearHist2D.begin();
00906       map<string,TH2F*>::iterator Tempend2D = TempNearHist2D.end();
00907       while(beg2D!=end2D){
00908         Tempbeg2D->second->Set(Tempbeg2D->second->GetSize(),
00909                                branches[Tempbeg2D->second->GetName()]);
00910         chisquare += Chi2Calc.GetChi2(beg2D->second,Tempbeg2D->second,hist2,
00911                                       Tempbeg2D->second->Integral() /
00912                                       beg2D->second->Integral());
00913         beg2D++;
00914         Tempbeg2D++;
00915       }
00916       //Far:
00917       beg1D = fFarHist1D.begin();
00918       end1D = fFarHist1D.end();
00919       Tempbeg1D = TempFarHist1D.begin();
00920       Tempend1D = TempFarHist1D.end();
00921       while(beg1D!=end1D){
00922         Tempbeg1D->second->Set(Tempbeg1D->second->GetSize(),
00923                                branches[Tempbeg1D->second->GetName()]);
00924         chisquare += Chi2Calc.GetChi2(beg1D->second,Tempbeg1D->second,hist1,
00925                                       Tempbeg1D->second->Integral() / 
00926                                       beg1D->second->Integral());
00927         beg1D++;
00928         Tempbeg1D++;
00929       }
00930       beg2D = fFarHist2D.begin();
00931       end2D = fFarHist2D.end();
00932       Tempbeg2D = TempFarHist2D.begin();
00933       Tempend2D = TempFarHist2D.end();
00934       while(beg2D!=end2D){
00935         Tempbeg2D->second->Set(Tempbeg2D->second->GetSize(),
00936                                branches[Tempbeg2D->second->GetName()]);
00937         chisquare += Chi2Calc.GetChi2(beg2D->second,Tempbeg2D->second,hist2,
00938                                       Tempbeg2D->second->Integral() / 
00939                                       beg2D->second->Integral());
00940         beg2D++;
00941         Tempbeg2D++;
00942       }
00943 
00944       //Fill chi2surface
00945       Int_t xbin = fChi2Surf->GetXaxis()->FindBin(par1);
00946       Int_t ybin = fChi2Surf->GetYaxis()->FindBin(par2);
00947       if(chisquare<fChi2Surf->GetBinContent(xbin,ybin))
00948         fChi2Surf->SetBinContent(xbin,ybin,chisquare);
00949 
00950       //if this is best fit so far:
00951       if(chisquare<ChiSq){
00952         ChiSq = chisquare;
00953         bestPar1 = par1;
00954         bestPar2 = par2;
00955         bestEntry = i;
00956       }
00957     }
00958 
00959     fReweightTree->GetEntry(bestEntry);
00960 
00961     map<string,TH1F*>::iterator Tempbeg1D = TempNearHist1D.begin();
00962     map<string,TH1F*>::iterator Tempend1D = TempNearHist1D.end();
00963     while(Tempbeg1D!=Tempend1D) {
00964       Tempbeg1D->second->Set(Tempbeg1D->second->GetSize(),
00965                              branches[Tempbeg1D->second->GetName()]);
00966       TString newname(Tempbeg1D->second->GetName());
00967       newname.ReplaceAll("TEMP_","Best_");
00968       TH1F *hist = (TH1F*) Tempbeg1D->second->Clone(newname.Data());      
00969       fBestFit1D.push_back(hist);
00970       Tempbeg1D++;
00971     }
00972     map<string,TH2F*>::iterator Tempbeg2D = TempNearHist2D.begin();
00973     map<string,TH2F*>::iterator Tempend2D = TempNearHist2D.end();
00974     while(Tempbeg2D!=Tempend2D) {
00975       Tempbeg2D->second->Set(Tempbeg2D->second->GetSize(),
00976                              branches[Tempbeg2D->second->GetName()]);
00977       TString newname(Tempbeg2D->second->GetName());
00978       newname.ReplaceAll("TEMP_","Best_");
00979       TH2F *hist = (TH2F*) Tempbeg2D->second->Clone(newname.Data());
00980       fBestFit2D.push_back(hist);
00981       Tempbeg2D++;
00982     }
00983     Tempbeg1D = TempFarHist1D.begin();
00984     Tempend1D = TempFarHist1D.end();
00985     while(Tempbeg1D!=Tempend1D) {
00986       Tempbeg1D->second->Set(Tempbeg1D->second->GetSize(),
00987                              branches[Tempbeg1D->second->GetName()]);
00988       TString newname(Tempbeg1D->second->GetName());
00989       newname.ReplaceAll("TEMP_","Best_");
00990       TH1F *hist = (TH1F*) Tempbeg1D->second->Clone(newname.Data());
00991       fBestFit1D.push_back(hist);
00992       Tempbeg1D++;
00993     }
00994     Tempbeg2D = TempFarHist2D.begin();
00995     Tempend2D = TempFarHist2D.end();
00996     while(Tempbeg2D!=Tempend2D) {
00997       Tempbeg2D->second->Set(Tempbeg2D->second->GetSize(),
00998                              branches[Tempbeg2D->second->GetName()]);
00999       TString newname(Tempbeg2D->second->GetName());
01000       newname.ReplaceAll("TEMP_","Best_");
01001       TH2F *hist = (TH2F*) Tempbeg2D->second->Clone(newname.Data());
01002       fBestFit2D.push_back(hist);
01003       Tempbeg2D++;
01004     }
01005     
01006     map<string,Float_t*>::iterator branchIter = branches.begin();
01007     map<string,Float_t*>::iterator branchEnd = branches.end();
01008     while(branchIter!=branchEnd) { delete [] branchIter->second; branchIter++; }
01009     branches.clear();
01010     
01011   }
01012 
01013   cout << "Best DeltaChi2 = " << ChiSq 
01014        << " Best Par1 = " << bestPar1 
01015        << " Best Par2 = " << bestPar2 << endl;
01016 
01017   //delete temp histos:
01018   for(int i=0;i<2;i++){
01019     map<string,TH1F*>::iterator beg1D;    
01020     map<string,TH1F*>::iterator end1D;
01021     if(i==0){
01022       beg1D = TempNearHist1D.begin();
01023       end1D = TempNearHist1D.end();
01024     }
01025     else {
01026       beg1D = TempFarHist1D.begin();
01027       end1D = TempFarHist1D.end();
01028     }
01029     while(beg1D!=end1D) { delete beg1D->second; beg1D++;}
01030   }
01031   TempNearHist1D.clear();
01032   TempFarHist1D.clear();
01033   for(int i=0;i<2;i++){
01034     map<string,TH2F*>::iterator beg2D;    
01035     map<string,TH2F*>::iterator end2D;
01036     if(i==0){
01037       beg2D = TempNearHist2D.begin();
01038       end2D = TempNearHist2D.end();
01039     }
01040     else {
01041       beg2D = TempFarHist2D.begin();
01042       end2D = TempFarHist2D.end();
01043     }
01044     while(beg2D!=end2D) { delete beg2D->second; beg2D++;}
01045   }
01046   TempNearHist2D.clear();
01047   TempFarHist2D.clear();
01048   
01049   return true;
01050 }

void PANAnalysis::EndJob (  ) 

Definition at line 1529 of file PANAnalysis.cxx.

References fBestFit1D, fBestFit2D, fChi2Surf, fFarHist1D, fFarHist2D, fNearHist1D, fNearHist2D, fOutFileTag, and fReweightTree.

Referenced by ~PANAnalysis().

01530 {
01531   cout << "In PANAnalysis::EndJob()" << endl;
01532   TString nam("PANAna_");
01533   nam+=fOutFileTag;
01534   nam+=".root";
01535   TFile *fSave = new TFile(nam.Data(),"RECREATE");
01536   fSave->cd();
01537 
01538   map<string,TH1F*>::iterator beg1D = fFarHist1D.begin();
01539   map<string,TH1F*>::iterator end1D = fFarHist1D.end();
01540   while(beg1D!=end1D) { beg1D->second->Write(); beg1D++;}
01541   beg1D = fNearHist1D.begin();
01542   end1D = fNearHist1D.end();
01543   while(beg1D!=end1D) { beg1D->second->Write(); beg1D++;}  
01544 
01545   map<string,TH2F*>::iterator beg2D = fFarHist2D.begin();
01546   map<string,TH2F*>::iterator end2D = fFarHist2D.end();
01547   while(beg2D!=end2D) { beg2D->second->Write(); beg2D++;}
01548   beg2D = fNearHist2D.begin();
01549   end2D = fNearHist2D.end();
01550   while(beg2D!=end2D) { beg2D->second->Write(); beg2D++;}
01551 
01552   vector<TH1F*>::iterator begBest1D = fBestFit1D.begin();
01553   vector<TH1F*>::iterator endBest1D = fBestFit1D.end();
01554   while(begBest1D!=endBest1D) { (*begBest1D)->Write(); begBest1D++;}
01555   vector<TH2F*>::iterator begBest2D = fBestFit2D.begin();
01556   vector<TH2F*>::iterator endBest2D = fBestFit2D.end();
01557   while(begBest2D!=endBest2D) { (*begBest2D)->Write(); begBest2D++;}
01558 
01559   if(fChi2Surf) fChi2Surf->Write();
01560   if(fReweightTree) fReweightTree->Write();
01561   fSave->Close();
01562   cout << "Closed file" << endl;
01563 }

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, fMCChainFar, fMCChainNear, 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.

00298 {
00299   cout << "In PANAnalysis::MakeMCVector()" << endl;
00300   TChain *chain = NULL;
00301   if(FN==1) chain = fMCChainNear;
00302   else if(FN==2) chain = fMCChainFar;
00303   else return false;
00304 
00305   Float_t true_enu = 0; chain->SetBranchAddress("true_enu",&true_enu);
00306   Float_t true_pxnu = 0; chain->SetBranchAddress("true_pxnu",&true_pxnu);
00307   Float_t true_pynu = 0; chain->SetBranchAddress("true_pynu",&true_pynu);
00308   Float_t true_pznu = 0; chain->SetBranchAddress("true_pznu",&true_pznu);
00309   Float_t true_etgt = 0; chain->SetBranchAddress("true_etgt",&true_etgt);
00310   Float_t true_pxtgt = 0; chain->SetBranchAddress("true_pxtgt",&true_pxtgt);
00311   Float_t true_pytgt = 0; chain->SetBranchAddress("true_pytgt",&true_pytgt);
00312   Float_t true_pztgt = 0; chain->SetBranchAddress("true_pztgt",&true_pztgt);
00313   Int_t flavour    = 0; chain->SetBranchAddress("flavour",&flavour);
00314   Int_t flavour_noosc = 0; chain->SetBranchAddress("nooscflavour",&flavour_noosc);
00315   Int_t cc_nc      = 0; chain->SetBranchAddress("cc_nc",&cc_nc);
00316   Int_t process    = 0; chain->SetBranchAddress("process",&process);
00317   Int_t init_state = 0; chain->SetBranchAddress("initial_state",&init_state);
00318   Int_t nucleus    = 0; chain->SetBranchAddress("nucleus",&nucleus);
00319   Int_t had_fs     = 0; chain->SetBranchAddress("had_fs",&had_fs);
00320   Float_t x        = 0; chain->SetBranchAddress("true_x",&x);
00321   Float_t y        = 0; chain->SetBranchAddress("true_y",&y);
00322   Float_t q2       = 0; chain->SetBranchAddress("true_q2",&q2);
00323   Float_t w2       = 0; chain->SetBranchAddress("true_w2",&w2);
00324   NuParent *parent = new NuParent();
00325   chain->SetBranchAddress("NuParent",parent);
00326   Int_t pass       = 0; chain->SetBranchAddress("pass",&pass);
00327   Int_t is_fid     = 0; chain->SetBranchAddress("is_fid",&is_fid);
00328   Int_t is_mc      = 0; chain->SetBranchAddress("is_mc",&is_mc);
00329 
00330   //Set branch addresses for reco quantities
00331   //1D first:
00332   map<string,TH1F*>::iterator beg1D;
00333   map<string,TH1F*>::iterator end1D;
00334   Int_t tot1DExpress = 0;
00335   if(FN==1) {
00336     beg1D = fNearHist1D.begin();
00337     end1D = fNearHist1D.end();
00338     tot1DExpress = fNearHist1D.size();
00339   }
00340   else {
00341     beg1D = fFarHist1D.begin();
00342     end1D = fFarHist1D.end();
00343     tot1DExpress = fFarHist1D.size();
00344   }
00345   if(tot1DExpress>10) {
00346     cout << "uh oh... >10 1D histos not supported yet!" << endl;
00347     return false;
00348   }
00349 
00350   TString leaftype_1D[10][10]; Double_t valD_1D[10][10] = {{0}};
00351   Int_t    valI_1D[10][10] = {{0}}; Float_t  valF_1D[10][10] = {{0}};
00352   TF1 *inputExpression1D[10];
00353   Int_t count_1D[10] = {0}; Int_t histCounter1D = 0;
00354   while(beg1D!=end1D){
00355     TString stringExpression1D(beg1D->first);
00356     TString funcName = "FUNC_" + stringExpression1D;
00357     TLeaf *leaf = 0; TIter leafIter(chain->GetListOfLeaves());
00358     while ((leaf = (TLeaf*) leafIter())) {
00359       if(stringExpression1D.Contains(leaf->GetName())) {
00360         leaftype_1D[histCounter1D][count_1D[histCounter1D]] = leaf->GetTypeName();
00361         char tempStr[5]; sprintf(tempStr,"[%i]",count_1D[histCounter1D]);
00362         stringExpression1D.ReplaceAll(leaf->GetName(),tempStr);
00363         if(strcmp(leaf->GetTypeName(),"Int_t")==0) {
00364           chain->SetBranchAddress(leaf->GetName(),
00365                                   &valI_1D[histCounter1D][count_1D[histCounter1D]]);
00366           count_1D[histCounter1D]++;
00367         }
00368         else if(strcmp(leaf->GetTypeName(),"Float_t")==0){
00369           chain->SetBranchAddress(leaf->GetName(),
00370                                   &valF_1D[histCounter1D][count_1D[histCounter1D]]);
00371           count_1D[histCounter1D]++;
00372         }
00373         else if(strcmp(leaf->GetTypeName(),"Double_t")==0){
00374           chain->SetBranchAddress(leaf->GetName(),
00375                                   &valD_1D[histCounter1D][count_1D[histCounter1D]]);
00376           count_1D[histCounter1D]++;
00377         }
00378         else {
00379           cout << "Can't include " << leaf->GetName() 
00380                << " in expressions at the moment..." << endl;
00381           return false;
00382         }
00383       }
00384     }
00385     inputExpression1D[histCounter1D] = 
00386       new TF1(funcName.Data(),stringExpression1D.Data(),0,1);
00387     histCounter1D++;
00388     beg1D++;
00389   }
00390 
00391   //now 2D:
00392   map<string,TH2F*>::iterator beg2D;
00393   map<string,TH2F*>::iterator end2D;
00394   Int_t tot2DExpress = 0;
00395   if(FN==1) {
00396     beg2D = fNearHist2D.begin();
00397     end2D = fNearHist2D.end();
00398     tot2DExpress = fNearHist2D.size();
00399   }
00400   else {
00401     beg2D = fFarHist2D.begin();
00402     end2D = fFarHist2D.end();
00403     tot2DExpress = fFarHist2D.size();
00404   }
00405   if(tot2DExpress>5) {
00406     cout << "uh oh... >5 2D histos not supported yet!" << endl;
00407     return false;
00408   }
00409 
00410   TString leaftype_2D[10][10]; Double_t valD_2D[10][10] = {{}};
00411   Int_t    valI_2D[10][10] = {{}}; Float_t  valF_2D[10][10] = {{}};
00412   TF1 *inputExpression2D[10];
00413   Int_t count_2D[10] = {}; Int_t histCounter2D = 0;
00414   while(beg2D!=end2D){
00415     for(int i=0;i<2;i++){
00416       TString stringExpression2D(beg2D->first);
00417       if(i==0) stringExpression2D.Remove(0,stringExpression2D.First(":")+1);
00418       else stringExpression2D.Remove(stringExpression2D.First(":"));
00419       TString funcName = "FUNC_" + stringExpression2D;
00420       TLeaf *leaf = 0; TIter leafIter(chain->GetListOfLeaves());
00421       while ((leaf = (TLeaf*) leafIter())) {
00422         if(stringExpression2D.Contains(leaf->GetName())) {
00423           leaftype_2D[histCounter2D][count_2D[histCounter2D]] = leaf->GetTypeName();
00424           char tempStr[5]; sprintf(tempStr,"[%i]",count_2D[histCounter2D]);
00425           stringExpression2D.ReplaceAll(leaf->GetName(),tempStr);
00426           if(strcmp(leaf->GetTypeName(),"Int_t")==0) {
00427             chain->SetBranchAddress(leaf->GetName(),
00428                                     &valI_2D[histCounter2D][count_2D[histCounter2D]]);
00429             count_2D[histCounter2D]++;
00430           }
00431           else if(strcmp(leaf->GetTypeName(),"Float_t")==0){
00432             chain->SetBranchAddress(leaf->GetName(),
00433                                     &valF_2D[histCounter2D][count_2D[histCounter2D]]);
00434             count_2D[histCounter2D]++;
00435           }
00436           else if(strcmp(leaf->GetTypeName(),"Double_t")==0){
00437             chain->SetBranchAddress(leaf->GetName(),
00438                                     &valD_2D[histCounter1D][count_2D[histCounter2D]]);
00439             count_2D[histCounter2D]++;
00440           }
00441           else {
00442             cout << "Can't include " << leaf->GetName() 
00443                  << " in expressions at the moment..." << endl;
00444             return false;
00445           }
00446         }
00447       }
00448       inputExpression2D[histCounter2D] = 
00449         new TF1(funcName.Data(),stringExpression2D.Data(),0,1);
00450       histCounter2D++;
00451     }
00452     beg2D++;
00453   }
00454   
00455   Int_t totEnt = chain->GetEntries();
00456   for(int i=entries;i<totEnt;i++){
00457     chain->GetEvent(i);
00458     if(pass&&is_fid&&is_mc){
00459       pan_mcev mmei;
00460       mmei.nu[0] = true_pxnu;
00461       mmei.nu[1] = true_pynu;
00462       mmei.nu[2] = true_pznu;
00463       mmei.nu[3] = true_enu;
00464       mmei.target[0] = true_pxtgt;
00465       mmei.target[1] = true_pytgt;
00466       mmei.target[2] = true_pztgt;
00467       mmei.target[3] = true_etgt;
00468       mmei.flavour = flavour;
00469       mmei.flavour_noosc = flavour_noosc;
00470       mmei.cc_nc = cc_nc;
00471       mmei.process = process;
00472       mmei.init_state = init_state;
00473       mmei.nucleus = nucleus;
00474       mmei.had_fs = had_fs;
00475       mmei.x = x;
00476       mmei.y = y;
00477       mmei.q2 = q2;
00478       mmei.w2 = w2;
00479       mmei.parvtx[0] = parent->GetX();
00480       mmei.parvtx[1] = parent->GetY();
00481       mmei.parvtx[2] = parent->GetZ();
00482       mmei.parmom[0] = parent->GetPx();
00483       mmei.parmom[1] = parent->GetPy();
00484       mmei.parmom[2] = parent->GetPz();
00485       mmei.parpid = parent->GetPID();
00486       mmei.pargen = parent->GetGen();
00487       for(int j=0;j<histCounter1D;j++){
00488         for(int k=0;k<count_1D[j];k++){
00489           if(leaftype_1D[j][k].BeginsWith("I"))
00490             inputExpression1D[j]->SetParameter(k,valI_1D[j][k]);
00491           else if(leaftype_1D[j][k].BeginsWith("F")) 
00492             inputExpression1D[j]->SetParameter(k,valF_1D[j][k]);
00493           else if(leaftype_1D[j][k].BeginsWith("D"))
00494             inputExpression1D[j]->SetParameter(k,valD_1D[j][k]);
00495         }
00496         TString ustring = inputExpression1D[j]->GetName();
00497         ustring.ReplaceAll("FUNC_","");
00498         mmei.userString[j] = ustring;
00499         mmei.userVar[j] = inputExpression1D[j]->Eval(0);
00500       }
00501       for(int j=0;j<histCounter2D;j++){
00502         for(int k=0;k<count_2D[j];k++){
00503           if(leaftype_2D[j][k].BeginsWith("I"))
00504             inputExpression2D[j]->SetParameter(k,valI_2D[j][k]);
00505           else if(leaftype_2D[j][k].BeginsWith("F")) 
00506             inputExpression2D[j]->SetParameter(k,valF_2D[j][k]);
00507           else if(leaftype_2D[j][k].BeginsWith("D"))
00508             inputExpression2D[j]->SetParameter(k,valD_2D[j][k]);
00509         }
00510         TString ustring = inputExpression2D[j]->GetName();
00511         ustring.ReplaceAll("FUNC_","");
00512         mmei.userString[tot1DExpress+j] = ustring;
00513         mmei.userVar[tot1DExpress+j] = inputExpression2D[j]->Eval(0);
00514       }
00515       if(FN==1) fNearMCEvent.push_back(mmei);
00516       else fFarMCEvent.push_back(mmei);
00517     }
00518   }
00519   delete parent;
00520   return true;
00521 }

Bool_t PANAnalysis::MakeReweightTree ( ReweightLooper looper,
Int_t  Par1Bins,
Float_t  Par1Min,
Float_t  Par1Max,
Int_t  Par2Bins,
Float_t  Par2Min,
Float_t  Par2Max 
)

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.

01054                                                                                      {
01055 
01056   cout << "In PANAnalysis::MakeReweightTree()" << endl;
01057   Float_t Par1Width=(Par1Max-Par1Min)/Float_t(Par1Bins);
01058   Float_t Par2Width=(Par2Max-Par2Min)/Float_t(Par2Bins);
01059 
01060   fReweightTree = new TTree("ReweightTree","ReweightTree");
01061 
01062   //maps of pointers to histograms to write to tree
01063   map<string,TH1F*> TempNearHist1D;
01064   map<string,TH1F*> TempFarHist1D;
01065   map<string,TH2F*> TempNearHist2D;
01066   map<string,TH2F*> TempFarHist2D;
01067   //copy dimension of 1D histos:
01068   for(int i=0;i<2;i++){
01069     map<string,TH1F*>::iterator beg1D;
01070     map<string,TH1F*>::iterator end1D;
01071     if(i==0){
01072       beg1D = fNearHist1D.begin();
01073       end1D = fNearHist1D.end();
01074     }
01075     else {
01076       beg1D = fFarHist1D.begin();
01077       end1D = fFarHist1D.end();
01078     }
01079     while(beg1D!=end1D){
01080       string histname = "TEMP_";
01081       histname.append(beg1D->second->GetName());
01082       TH1F *temp = (TH1F*) beg1D->second->Clone(histname.c_str());
01083       temp->Reset();
01084       if(i==0) TempNearHist1D[beg1D->first] = temp;
01085       else TempFarHist1D[beg1D->first] = temp;
01086       char leaflist[256]; 
01087       sprintf(leaflist,"%s[%i]/F",beg1D->second->GetName(),
01088               temp->GetSize());
01089       fReweightTree->Branch(beg1D->second->GetName(),
01090                             temp->GetArray(),leaflist);
01091       beg1D++;
01092     }
01093   }
01094   //copy dimensions of 2D histos:
01095   for(int i=0;i<2;i++){
01096     map<string,TH2F*>::iterator beg2D;
01097     map<string,TH2F*>::iterator end2D;
01098     if(i==0){
01099       beg2D = fNearHist2D.begin();
01100       end2D = fNearHist2D.end();
01101     }
01102     else {
01103       beg2D = fFarHist2D.begin();
01104       end2D = fFarHist2D.end();
01105     }
01106     while(beg2D!=end2D){
01107       string histname = "TEMP_";
01108       histname.append(beg2D->second->GetName());
01109       TH2F *temp = (TH2F*) beg2D->second->Clone(histname.c_str());
01110       temp->Reset();
01111       if(i==0) TempNearHist2D[beg2D->first] = temp;
01112       else TempFarHist2D[beg2D->first] = temp;
01113       char leaflist[256]; 
01114       sprintf(leaflist,"%s[%i]/F",beg2D->second->GetName(),
01115               temp->GetSize());
01116       fReweightTree->Branch(beg2D->second->GetName(),
01117                             temp->GetArray(),leaflist);
01118       beg2D++;
01119     }
01120   }
01121 
01122   //make a map of where each user variable is held in pan_mcev struct:
01123   Int_t totNearUserVar = 0;
01124   Int_t *stringMatchNear = MakeVarMapNear(totNearUserVar);
01125   cout << "Total Number of NearDet User Variables = " << totNearUserVar << endl;
01126   Int_t totFarUserVar = 0;
01127   Int_t *stringMatchFar = MakeVarMapFar(totFarUserVar);
01128   cout << "Total Number of FarDet User Variables = " << totFarUserVar << endl;
01129 
01130   //add reweight variable values to tree:
01131   Registry rwtconfig;
01132   map<string,Double_t*> variables;
01133   ReweightLooper *daughter = looper;
01134   if(daughter) {
01135     Double_t *temp = new Double_t[1];
01136     variables[daughter->GetName()] = temp;
01137     TString leafname = daughter->GetName(); leafname.ReplaceAll(":","_");
01138     TString leaflist = leafname; leaflist.Append("/D");
01139     fReweightTree->Branch(leafname.Data(),temp,leaflist.Data());
01140   }
01141   while((daughter = daughter->GetDaughter())) {
01142     Double_t *temp = new Double_t[1];
01143     variables[daughter->GetName()] = temp;
01144     TString leafname = daughter->GetName(); leafname.ReplaceAll(":","_");
01145     TString leaflist = leafname; leaflist.Append("/D");    
01146     fReweightTree->Branch(leafname.Data(),temp,leaflist.Data());
01147   }
01148 
01149   //Set branch address for oscillation parameters:
01150   float par1 = 0;  float par2 = 0;
01151   fReweightTree->Branch("par1",&par1,"par1/F");
01152   fReweightTree->Branch("par2",&par2,"par2/F");
01153 
01154   //Get instance of MCReweight object:
01155   MCReweight &fReweight = MCReweight::Instance();
01156 
01157   //get iterators for event vectors
01158   vector<pan_mcev>::iterator begNearEv = fNearMCEvent.begin();
01159   vector<pan_mcev>::iterator endNearEv = fNearMCEvent.end();
01160   vector<pan_mcev>::iterator begFarEv = fFarMCEvent.begin();
01161   vector<pan_mcev>::iterator endFarEv = fFarMCEvent.end();
01162 
01163   //begin loop over reweight variables:
01164   while(looper->Grind(variables)){
01165     //set reweight registry:
01166     //fReweight.ResetAllReweightConfigs();
01167     map<string,Double_t*>::iterator begVarIter = variables.begin();
01168     map<string,Double_t*>::iterator endVarIter = variables.end();
01169     rwtconfig.UnLockKeys();
01170     rwtconfig.UnLockValues();
01171     while(begVarIter!=endVarIter){
01172       rwtconfig.Set(begVarIter->first.c_str(),(*begVarIter->second));      
01173       begVarIter++;
01174     }
01175     rwtconfig.LockKeys();
01176     rwtconfig.LockValues();
01177     fReweight.ComputeWeight(0,&rwtconfig); //internally sets params
01178 
01179     //loop over oscillation parameters
01180     for(int i=0;i<Par1Bins;i++){
01181       par1 = Par1Min + (float(i+1)*Par1Width) - Par1Width/2.;
01182       if(i%10==0) cout << "Doing Par1 = " << par1 << endl;
01183       for(int j=0;j<Par2Bins;j++){
01184         par2 = Par2Min + (float(j+1)*Par2Width) - Par2Width/2.;
01185         
01186         //Set oscillation parameters:
01187         fOscillationFunction->SetParameter(fVaryX[0],par1);
01188         fOscillationFunction->SetParameter(fVaryX[1],par2);
01189         
01190         //Reset Temp histos:
01191         for(int k=0;k<2;k++){
01192           map<string,TH1F*>::iterator beg1D;    
01193           map<string,TH1F*>::iterator end1D;
01194           if(k==0){
01195             beg1D = TempNearHist1D.begin();
01196             end1D = TempNearHist1D.end();
01197           }
01198           else {
01199             beg1D = TempFarHist1D.begin();
01200             end1D = TempFarHist1D.end();
01201           }
01202           while(beg1D!=end1D) { beg1D->second->Reset(); beg1D++;}
01203         }
01204         for(int k=0;k<2;k++){
01205           map<string,TH2F*>::iterator beg2D;    
01206           map<string,TH2F*>::iterator end2D;
01207           if(k==0){
01208             beg2D = TempNearHist2D.begin();
01209             end2D = TempNearHist2D.end();
01210           }
01211           else {
01212             beg2D = TempFarHist2D.begin();
01213             end2D = TempFarHist2D.end();
01214           }
01215           while(beg2D!=end2D) { beg2D->second->Reset(); beg2D++;}
01216         }
01217         
01218         //start loop over events:
01219         //Registry event;
01220         //event.UnLockKeys();
01221         //event.UnLockValues();
01222         MCEventInfo *eventinfo = new MCEventInfo();
01223         NuParent *parent = new NuParent();
01224         //Near:
01225         begNearEv = fNearMCEvent.begin();
01226         fOscillationFunction->SetParameter(2,1.0); //set baseline to 1km
01227         while(begNearEv!=endNearEv){
01228           //fill event registry:
01229           eventinfo->nuE           = begNearEv->nu[3];
01230           eventinfo->nuPx          = begNearEv->nu[0];
01231           eventinfo->nuPy          = begNearEv->nu[1];
01232           eventinfo->nuPz          = begNearEv->nu[2];
01233           eventinfo->tarE          = begNearEv->target[3];
01234           eventinfo->tarPx         = begNearEv->target[0];
01235           eventinfo->tarPy         = begNearEv->target[1];
01236           eventinfo->tarPz         = begNearEv->target[2];
01237           eventinfo->x             = begNearEv->x;
01238           eventinfo->y             = begNearEv->y;
01239           eventinfo->q2            = begNearEv->q2;
01240           eventinfo->w2            = begNearEv->w2;
01241           eventinfo->iaction       = begNearEv->cc_nc;
01242           eventinfo->inu           = begNearEv->flavour;
01243           eventinfo->iresonance    = begNearEv->process;
01244           eventinfo->initial_state = begNearEv->init_state;
01245           eventinfo->nucleus       = begNearEv->nucleus;
01246           eventinfo->had_fs        = begNearEv->had_fs;
01247 
01248           parent->SetX(begNearEv->parvtx[0]);
01249           parent->SetY(begNearEv->parvtx[1]);
01250           parent->SetZ(begNearEv->parvtx[2]);
01251           parent->SetPx(begNearEv->parmom[0]);
01252           parent->SetPy(begNearEv->parmom[1]);
01253           parent->SetPz(begNearEv->parmom[2]);
01254           parent->SetPID(begNearEv->parpid);
01255           parent->SetGen(begNearEv->pargen);
01256           /*
01257           event.Set("event:nu_en",begNearEv->nu[3]);
01258           event.Set("event:nu_px",begNearEv->nu[0]);
01259           event.Set("event:nu_py",begNearEv->nu[1]);
01260           event.Set("event:nu_pz",begNearEv->nu[2]);
01261           event.Set("event:tar_en",begNearEv->target[3]);
01262           event.Set("event:tar_px",begNearEv->target[0]);
01263           event.Set("event:tar_py",begNearEv->target[1]);
01264           event.Set("event:tar_pz",begNearEv->target[2]);
01265           event.Set("event:x",begNearEv->x);
01266           event.Set("event:y",begNearEv->y);
01267           event.Set("event:q2",begNearEv->q2);
01268           event.Set("event:w2",begNearEv->w2);
01269           event.Set("event:iaction",begNearEv->cc_nc);
01270           event.Set("event:inu",begNearEv->flavour);
01271           event.Set("event:process",begNearEv->process);
01272           event.Set("event:initial_state",begNearEv->init_state);
01273           event.Set("event:nucleus",begNearEv->nucleus);
01274           event.Set("event:had_fs",begNearEv->had_fs);
01275           event.Set("event:nuparent_x",begNearEv->parvtx[0]);
01276           event.Set("event:nuparent_y",begNearEv->parvtx[1]);
01277           event.Set("event:nuparent_z",begNearEv->parvtx[2]);
01278           event.Set("event:nuparent_px",begNearEv->parmom[0]);
01279           event.Set("event:nuparent_py",begNearEv->parmom[1]);
01280           event.Set("event:nuparent_pz",begNearEv->parmom[2]);
01281           event.Set("event:nuparent_pid",begNearEv->parpid);
01282           event.Set("event:nuparent_gen",begNearEv->pargen);
01283           */
01284           //set oscillation dependent params:
01285           fOscillationFunction->SetParameter(0,begNearEv->flavour);
01286           fOscillationFunction->SetParameter(1,begNearEv->flavour_noosc);
01287           Double_t oscprob = 1;
01288           if(begFarEv->cc_nc==1)
01289             oscprob = fOscillationFunction->Eval(begNearEv->nu[3]);
01290           Int_t counter = 0;
01291           map<string,TH1F*>::iterator beg1D = TempNearHist1D.begin();
01292           map<string,TH1F*>::iterator end1D = TempNearHist1D.end();
01293           map<string,TH2F*>::iterator beg2D = TempNearHist2D.begin();
01294           map<string,TH2F*>::iterator end2D = TempNearHist2D.end();
01295           while(beg1D!=end1D){
01296             beg1D->second->Fill(begNearEv->userVar[stringMatchNear[counter]],
01297                                 oscprob*fReweight.ComputeWeight(eventinfo,parent));
01298             counter++;
01299             beg1D++;
01300           }
01301           while(beg2D!=end2D){
01302             beg2D->second->Fill(begNearEv->userVar[stringMatchNear[counter]],
01303                                 begNearEv->userVar[stringMatchNear[counter+1]],
01304                                 oscprob*fReweight.ComputeWeight(eventinfo,parent));
01305             counter+=2;
01306             beg2D++;
01307           }
01308           begNearEv++;
01309         }
01310         //Far:
01311         begFarEv = fFarMCEvent.begin();
01312         fOscillationFunction->SetParameter(2,735.0); //set baseline to 735km
01313         while(begFarEv!=endFarEv){        
01314           //fill event registry:
01315           eventinfo->nuE           = begFarEv->nu[3];
01316           eventinfo->nuPx          = begFarEv->nu[0];
01317           eventinfo->nuPy          = begFarEv->nu[1];
01318           eventinfo->nuPz          = begFarEv->nu[2];
01319           eventinfo->tarE          = begFarEv->target[3];
01320           eventinfo->tarPx         = begFarEv->target[0];
01321           eventinfo->tarPy         = begFarEv->target[1];
01322           eventinfo->tarPz         = begFarEv->target[2];
01323           eventinfo->x             = begFarEv->x;
01324           eventinfo->y             = begFarEv->y;
01325           eventinfo->q2            = begFarEv->q2;
01326           eventinfo->w2            = begFarEv->w2;
01327           eventinfo->iaction       = begFarEv->cc_nc;
01328           eventinfo->inu           = begFarEv->flavour;
01329           eventinfo->iresonance    = begFarEv->process;
01330           eventinfo->initial_state = begFarEv->init_state;
01331           eventinfo->nucleus       = begFarEv->nucleus;
01332           eventinfo->had_fs        = begFarEv->had_fs;
01333 
01334           parent->SetX(begFarEv->parvtx[0]);
01335           parent->SetY(begFarEv->parvtx[1]);
01336           parent->SetZ(begFarEv->parvtx[2]);
01337           parent->SetPx(begFarEv->parmom[0]);
01338           parent->SetPy(begFarEv->parmom[1]);
01339           parent->SetPz(begFarEv->parmom[2]);
01340           parent->SetPID(begFarEv->parpid);
01341           parent->SetGen(begFarEv->pargen);
01342           /*
01343           event.Set("event:nu_en",begFarEv->nu[3]);
01344           event.Set("event:nu_px",begFarEv->nu[0]);
01345           event.Set("event:nu_py",begFarEv->nu[1]);
01346           event.Set("event:nu_pz",begFarEv->nu[2]);
01347           event.Set("event:tar_en",begFarEv->target[3]);
01348           event.Set("event:tar_px",begFarEv->target[0]);
01349           event.Set("event:tar_py",begFarEv->target[1]);
01350           event.Set("event:tar_pz",begFarEv->target[2]);
01351           event.Set("event:x",begFarEv->x);
01352           event.Set("event:y",begFarEv->y);
01353           event.Set("event:q2",begFarEv->q2);
01354           event.Set("event:w2",begFarEv->w2);
01355           event.Set("event:iaction",begFarEv->cc_nc);
01356           event.Set("event:inu",begFarEv->flavour);
01357           event.Set("event:process",begFarEv->process);
01358           event.Set("event:initial_state",begFarEv->init_state);
01359           event.Set("event:nucleus",begFarEv->nucleus);
01360           event.Set("event:had_fs",begFarEv->had_fs);
01361           event.Set("event:nuparent_x",begFarEv->parvtx[0]);
01362           event.Set("event:nuparent_y",begFarEv->parvtx[1]);
01363           event.Set("event:nuparent_z",begFarEv->parvtx[2]);
01364           event.Set("event:nuparent_px",begFarEv->parmom[0]);
01365           event.Set("event:nuparent_py",begFarEv->parmom[1]);
01366           event.Set("event:nuparent_pz",begFarEv->parmom[2]);
01367           event.Set("event:nuparent_pid",begFarEv->parpid);
01368           event.Set("event:nuparent_gen",begFarEv->pargen);
01369           */
01370           //set oscillation dependent params:
01371           fOscillationFunction->SetParameter(0,begFarEv->flavour);
01372           fOscillationFunction->SetParameter(1,begFarEv->flavour_noosc);
01373           Double_t oscprob = 1;
01374           if(begFarEv->cc_nc==1)
01375             oscprob = fOscillationFunction->Eval(begFarEv->nu[3]);
01376           Int_t counter = 0;
01377           map<string,TH1F*>::iterator beg1D = TempFarHist1D.begin();
01378           map<string,TH1F*>::iterator end1D = TempFarHist1D.end();
01379           map<string,TH2F*>::iterator beg2D = TempFarHist2D.begin();
01380           map<string,TH2F*>::iterator end2D = TempFarHist2D.end();
01381           while(beg1D!=end1D){
01382             beg1D->second->Fill(begFarEv->userVar[stringMatchFar[counter]],
01383                                 oscprob*fReweight.ComputeWeight(eventinfo,parent));
01384             counter++;
01385             beg1D++;
01386           }
01387           while(beg2D!=end2D){
01388             beg2D->second->Fill(begFarEv->userVar[stringMatchFar[counter]],
01389                                 begFarEv->userVar[stringMatchFar[counter+1]],
01390                                 oscprob*fReweight.ComputeWeight(eventinfo,parent));
01391             counter+=2;
01392             beg2D++;
01393           }
01394           begFarEv++;
01395         }
01396         delete eventinfo;
01397         delete parent;
01398         fReweightTree->Fill();
01399       }
01400     }
01401   }
01402   looper->ResetCounter();
01403 
01404   //delete temp histos:
01405   for(int i=0;i<2;i++){
01406     map<string,TH1F*>::iterator beg1D;    
01407     map<string,TH1F*>::iterator end1D;
01408     if(i==0){
01409       beg1D = TempNearHist1D.begin();
01410       end1D = TempNearHist1D.end();
01411     }
01412     else {
01413       beg1D = TempFarHist1D.begin();
01414       end1D = TempFarHist1D.end();
01415     }
01416     while(beg1D!=end1D) { delete beg1D->second; beg1D++;}
01417   }
01418   TempNearHist1D.clear();
01419   TempFarHist1D.clear();
01420   for(int i=0;i<2;i++){
01421     map<string,TH2F*>::iterator beg2D;    
01422     map<string,TH2F*>::iterator end2D;
01423     if(i==0){
01424       beg2D = TempNearHist2D.begin();
01425       end2D = TempNearHist2D.end();
01426     }
01427     else {
01428       beg2D = TempFarHist2D.begin();
01429       end2D = TempFarHist2D.end();
01430     }
01431     while(beg2D!=end2D) { delete beg2D->second; beg2D++;}
01432   }
01433   TempNearHist2D.clear();
01434   TempFarHist2D.clear();
01435   //delete string maps
01436   delete [] stringMatchNear;
01437   delete [] stringMatchFar;
01438   //delete vector
01439   map<string,Double_t*>::iterator beg = variables.begin();
01440   map<string,Double_t*>::iterator end = variables.end();
01441   while(beg!=end){ delete [] (beg->second); beg++;}
01442   return true;
01443 }

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

Definition at line 1488 of file PANAnalysis.cxx.

References fFarHist1D, fFarHist2D, and fFarMCEvent.

Referenced by Do2DFit(), and MakeReweightTree().

01488                                                       {
01489   totFarUserVar = fFarHist1D.size() + 2*fFarHist2D.size();
01490   Int_t *stringMatchFar = NULL;
01491   vector<pan_mcev>::iterator begFarEv = fFarMCEvent.begin();
01492   vector<pan_mcev>::iterator endFarEv = fFarMCEvent.end();
01493   if(totFarUserVar>0 && begFarEv!=endFarEv){
01494     Int_t counter = 0;
01495     stringMatchFar = new Int_t[totFarUserVar];    
01496     map<string,TH1F*>::iterator beg1D = fFarHist1D.begin();
01497     map<string,TH1F*>::iterator end1D = fFarHist1D.end();
01498     while(beg1D!=end1D){
01499       for(int i=0;i<totFarUserVar;i++){
01500         if(strcmp((beg1D->first).c_str(),(begFarEv->userString[i]).c_str())==0) {         
01501           stringMatchFar[counter] = i;
01502           break;
01503         }
01504       }
01505       counter++;
01506       beg1D++;
01507     }
01508     map<string,TH2F*>::iterator beg2D = fFarHist2D.begin();
01509     map<string,TH2F*>::iterator end2D = fFarHist2D.end();
01510     while(beg2D!=end2D){
01511       for(int i=0;i<2;i++){
01512         TString histVar = beg2D->first; 
01513         if(i==0) histVar.Remove(0,histVar.First(":")+1);
01514         else histVar.Remove(histVar.First(":"));
01515         for(int j=0;j<totFarUserVar;j++){
01516           if(strcmp(histVar.Data(),(begFarEv->userString[j]).c_str())==0) {
01517             stringMatchFar[counter] = j;
01518             break;
01519           }
01520         }
01521         counter++;
01522       }
01523       beg2D++;
01524     }
01525   }
01526   return stringMatchFar;
01527 }

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

Definition at line 1446 of file PANAnalysis.cxx.

References fNearHist1D, fNearHist2D, and fNearMCEvent.

Referenced by Do2DFit(), and MakeReweightTree().

01446                                                         {
01447   totNearUserVar = fNearHist1D.size() + 2*fNearHist2D.size();
01448   Int_t *stringMatchNear = NULL;
01449   vector<pan_mcev>::iterator begNearEv = fNearMCEvent.begin();
01450   vector<pan_mcev>::iterator endNearEv = fNearMCEvent.end();
01451   if(totNearUserVar>0 && begNearEv!=endNearEv){
01452     Int_t counter = 0;
01453     stringMatchNear = new Int_t[totNearUserVar];
01454     map<string,TH1F*>::iterator beg1D = fNearHist1D.begin();
01455     map<string,TH1F*>::iterator end1D = fNearHist1D.end();
01456     while(beg1D!=end1D){
01457       for(int i=0;i<totNearUserVar;i++){
01458         if(strcmp((beg1D->first).c_str(),(begNearEv->userString[i]).c_str())==0) {        
01459           stringMatchNear[counter] = i;
01460           break;
01461         }
01462       }
01463       counter++;
01464       beg1D++;
01465     }
01466     map<string,TH2F*>::iterator beg2D = fNearHist2D.begin();
01467     map<string,TH2F*>::iterator end2D = fNearHist2D.end();
01468     while(beg2D!=end2D){
01469       for(int i=0;i<2;i++){
01470         TString histVar = beg2D->first;
01471         if(i==0) histVar.Remove(0,histVar.First(":")+1);
01472         else histVar.Remove(histVar.First(":"));
01473         for(int j=0;j<totNearUserVar;j++){
01474           if(strcmp(histVar.Data(),(begNearEv->userString[j]).c_str())==0) {
01475             stringMatchNear[counter] = j;
01476             break;
01477           }
01478         }
01479         counter++;
01480       }
01481       beg2D++;
01482     }
01483   }
01484   return stringMatchNear;
01485 }

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, fMCChainFar, fMCChainNear, fNearHist1D, fNearHist2D, fNumberOfBinsX, fNumberOfBinsY, fOscillationFunction, fRangeLowerLimitX, fRangeLowerLimitY, fRangeUpperLimitX, and fRangeUpperLimitY.

00124 {
00125   cout << "In PANAnalysis::RecoDist()" << endl;
00126   
00127   Double_t baseLine = 735.0;
00128   if(FN==1) baseLine = 1.0;
00129   fOscillationFunction->SetParameter(2,baseLine);
00130 
00131   Bool_t is2D = false;
00132   TString stringExpression1D(expression);
00133   if(stringExpression1D.Contains(":")) is2D = true;
00134 
00135   TH1F *hist1D = NULL;
00136   TH2F *hist2D = NULL;
00137   if(!is2D) hist1D = new TH1F(name.c_str(),name.c_str(),fNumberOfBinsX,
00138                               fRangeLowerLimitX,fRangeUpperLimitX);
00139   
00140   else hist2D = new TH2F(name.c_str(),name.c_str(),fNumberOfBinsX,
00141                          fRangeLowerLimitX,fRangeUpperLimitX,
00142                          fNumberOfBinsY,fRangeLowerLimitY,fRangeUpperLimitY);
00143 
00144   string toHisto = ">>";
00145   if(FN==2&&fDataChainFar){
00146     fDataChainFar->Draw((expression+toHisto+name).c_str(),
00147                         "pass&&is_fid");
00148     if(!is2D) fFarHist1D[expression] = hist1D;
00149     else fFarHist2D[expression] = hist2D;
00150     return true;
00151   }
00152   else if(FN==1&&fDataChainNear){
00153     fDataChainNear->Draw((expression+toHisto+name).c_str(),
00154                          "pass&&is_fid");
00155     if(!is2D) fNearHist1D[expression] = hist1D;
00156     else fNearHist2D[expression] = hist2D;
00157     return true;
00158   }
00159   else {
00160     TChain *chain = NULL;
00161     if(FN==2&&fMCChainFar) chain = fMCChainFar;
00162     else if(FN==1&&fMCChainNear) chain = fMCChainNear;
00163     else return false;
00164 
00165     //Set branch address for true energy
00166     Float_t trueEnergy = 0;  chain->SetBranchAddress("true_enu",&trueEnergy);
00167     Int_t flavour = 0;       chain->SetBranchAddress("flavour",&flavour);
00168     Int_t flavour_noosc = 0; chain->SetBranchAddress("nooscflavour",&flavour_noosc);
00169     Int_t cc_nc = 0;         chain->SetBranchAddress("cc_nc",&cc_nc);
00170     Int_t pass = 0;          chain->SetBranchAddress("pass",&pass);
00171     Int_t is_fid = 0;        chain->SetBranchAddress("is_fid",&is_fid);
00172     Int_t is_mc = 0;         chain->SetBranchAddress("is_mc",&is_mc);
00173 
00174     //Set branch addresses for reco quantities    
00175     TString stringExpression2D;
00176     if(is2D) {
00177       stringExpression1D.Remove(0,stringExpression1D.First(":")+1);
00178       stringExpression2D = expression;
00179       stringExpression2D.Remove(stringExpression2D.First(":"));
00180     }
00181 
00182     TLeaf *leaf = 0; TIter leafIter(chain->GetListOfLeaves());
00183     
00184     TString leaftype_1D[10]; Double_t valD_1D[10] = {0};
00185     Int_t valI_1D[10] = {0}; Float_t valF_1D[10] = {0};
00186 
00187     TString leaftype_2D[10]; Double_t valD_2D[10] = {0};
00188     Int_t valI_2D[10] = {0}; Float_t valF_2D[10] = {0};
00189 
00190     Int_t count_1D = 0;  Int_t count_2D = 0;    
00191     
00192     while ((leaf = (TLeaf*) leafIter())) {
00193       if(stringExpression1D.Contains(leaf->GetName())) {
00194         leaftype_1D[count_1D] = leaf->GetTypeName();
00195         char tempStr[5]; sprintf(tempStr,"[%i]",count_1D);      
00196         stringExpression1D.ReplaceAll(leaf->GetName(),tempStr);
00197         if(leaftype_1D[count_1D]=="Int_t"){
00198           chain->SetBranchAddress(leaf->GetName(),&valI_1D[count_1D]);
00199           count_1D++;
00200         }
00201         else if(leaftype_1D[count_1D]=="Float_t"){
00202           chain->SetBranchAddress(leaf->GetName(),&valF_1D[count_1D]);
00203           count_1D++;
00204         }
00205         else if(leaftype_1D[count_1D]=="Double_t"){
00206           chain->SetBranchAddress(leaf->GetName(),&valD_1D[count_1D]);
00207           count_1D++;
00208         }
00209         else {
00210           cout << "Can't include " << leaf->GetName() 
00211                << " in expressions at the moment..." << endl;
00212           return false;
00213         }
00214       }
00215       if(is2D){
00216         if(stringExpression2D.Contains(leaf->GetName())) {       
00217           leaftype_2D[count_2D] = leaf->GetTypeName();
00218           char tempStr[5]; sprintf(tempStr,"[%i]",count_2D);
00219           stringExpression2D.ReplaceAll(leaf->GetName(),tempStr);
00220           if(leaftype_2D[count_2D]=="Int_t") {
00221             chain->SetBranchAddress(leaf->GetName(),&valI_2D[count_2D]);
00222             count_2D++;
00223           }
00224           else if(leaftype_2D[count_2D]=="Float_t"){
00225             chain->SetBranchAddress(leaf->GetName(),&valF_2D[count_2D]);
00226             count_2D++;
00227           }
00228           else if(leaftype_2D[count_2D]=="Double_t"){
00229             chain->SetBranchAddress(leaf->GetName(),&valD_2D[count_2D]);
00230             count_2D++;
00231           }
00232           else {
00233             cout << "Can't include " << leaf->GetName() 
00234                  << " in expressions at the moment..." << endl;
00235             return false;
00236           }
00237         }
00238       }
00239     }
00240 
00241     TF1 *inputExpression1D = 
00242       new TF1("inputExpression1D",stringExpression1D.Data(),0,1);
00243     TF1 *inputExpression2D = NULL;
00244     if(is2D) inputExpression2D = 
00245       new TF1("inputExpression2D",stringExpression2D.Data(),0,1);
00246 
00247     if(entries<0) entries = chain->GetEntries();
00248     for(int i=0;i<entries;i++){
00249       chain->GetEntry(i);
00250       if(pass&&is_fid&&is_mc){
00251         //get survival probability:
00252         fOscillationFunction->SetParameter(0,flavour);
00253         fOscillationFunction->SetParameter(1,flavour_noosc);
00254         Double_t oscprob = 1;
00255         if(cc_nc==1) oscprob = fOscillationFunction->Eval(trueEnergy);
00256         //Fill histos:
00257         for(int j=0;j<count_1D;j++){      
00258           if(leaftype_1D[j].BeginsWith("I"))
00259             inputExpression1D->SetParameter(j,double(valI_1D[j]));
00260           else if(leaftype_1D[j].BeginsWith("F")) 
00261             inputExpression1D->SetParameter(j,double(valF_1D[j]));
00262           else if(leaftype_1D[j].BeginsWith("D"))
00263             inputExpression1D->SetParameter(j,double(valD_1D[j]));
00264         }
00265         if(is2D){
00266           for(int j=0;j<count_2D;j++){
00267             if(leaftype_2D[j].BeginsWith("I"))
00268               inputExpression2D->SetParameter(j,double(valI_2D[j]));
00269             else if(leaftype_2D[j].BeginsWith("F")) 
00270               inputExpression2D->SetParameter(j,double(valF_2D[j]));
00271             else if(leaftype_2D[j].BeginsWith("D"))
00272               inputExpression2D->SetParameter(j,double(valD_2D[j]));
00273           }
00274           hist2D->Fill(inputExpression1D->Eval(0),
00275                        inputExpression2D->Eval(0),oscprob);
00276         }
00277         else hist1D->Fill(inputExpression1D->Eval(0),oscprob);
00278       }
00279     }
00280     delete inputExpression1D;
00281     if(is2D) delete inputExpression2D;
00282 
00283     if(FN==2&&fMCChainFar){
00284       if(!is2D) fFarHist1D[expression] = hist1D;
00285       else fFarHist2D[expression] = hist2D;
00286       return true;
00287     }
00288     else if(FN==1&&fMCChainNear){
00289       if(!is2D) fNearHist1D[expression] = hist1D;
00290       else fNearHist2D[expression] = hist2D;
00291       return true;
00292     }
00293   }
00294   return false;
00295 }

void PANAnalysis::SetFileNameTag ( std::string  tag  ) 

Definition at line 76 of file PANAnalysis.cxx.

References fOutFileTag.

00077 {
00078   fOutFileTag = tag;
00079 }

void PANAnalysis::SetHistBookInfo ( int  binsX,
float  lowX,
float  upX,
int  binsY,
float  lowY,
float  upY 
)

Definition at line 88 of file PANAnalysis.cxx.

References fNumberOfBinsX, fNumberOfBinsY, fRangeLowerLimitX, fRangeLowerLimitY, fRangeUpperLimitX, and fRangeUpperLimitY.

00090 {
00091   fNumberOfBinsX = binsX;
00092   fRangeLowerLimitX = lowX;
00093   fRangeUpperLimitX = upX;
00094   fNumberOfBinsY = binsY;
00095   fRangeLowerLimitY = lowY;
00096   fRangeUpperLimitY = upY;
00097 }

void PANAnalysis::SetHistBookInfo ( int  bins,
float  low,
float  up 
)

Definition at line 81 of file PANAnalysis.cxx.

References fNumberOfBinsX, fRangeLowerLimitX, and fRangeUpperLimitX.

00082 {
00083   fNumberOfBinsX = bins;
00084   fRangeLowerLimitX = low;
00085   fRangeUpperLimitX = up;
00086 }

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

Definition at line 112 of file PANAnalysis.cxx.

References fNumOscPars, fOscillationFunction, and fOscillationParameters.

00113 {
00114   fOscillationParameters = params;
00115   fOscillationFunction = new TF1(*form);
00116   fNumOscPars = fOscillationFunction->GetNpar();
00117   for(int i=0;i<fNumOscPars;i++){
00118     fOscillationFunction->SetParameter(i,fOscillationParameters[i]);
00119   }
00120 }

void PANAnalysis::SetOscillationFunction ( Double_t(*)(Double_t *, Double_t *)  fcn,
Float_t  lowend,
Float_t  highend,
int  numpars,
double *  params 
)

Definition at line 99 of file PANAnalysis.cxx.

References fcn(), fNumOscPars, fOscillationFunction, and fOscillationParameters.

00102 {
00103   fOscillationParameters = params;
00104   fOscillationFunction = new TF1("OscFunc",fcn,lowend,
00105                                  highend,numpars);  
00106   fNumOscPars = numpars; 
00107   for(int i=0;i<fNumOscPars;i++){
00108     fOscillationFunction->SetParameter(i,fOscillationParameters[i]);
00109   } 
00110 }

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

Definition at line 101 of file PANAnalysis.h.

References fReweightTree.

00101 {fReweightTree = tree;}

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

Definition at line 94 of file PANAnalysis.h.

References fVaryX.

00094 {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(), PANAnalysis(), 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]
map<string,TH2F*> PANAnalysis::fFarHist2D [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 MakeMCVector(), PANAnalysis(), and RecoDist().

TChain* PANAnalysis::fMCChainNear [protected]

Definition at line 44 of file PANAnalysis.h.

Referenced by MakeMCVector(), PANAnalysis(), and RecoDist().

map<string,TH1F*> PANAnalysis::fNearHist1D [protected]
map<string,TH2F*> PANAnalysis::fNearHist2D [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(), RecoDist(), and SetHistBookInfo().

int PANAnalysis::fNumberOfBinsY [protected]

Definition at line 52 of file PANAnalysis.h.

Referenced by PANAnalysis(), RecoDist(), and SetHistBookInfo().

int PANAnalysis::fNumOscPars [protected]

Definition at line 57 of file PANAnalysis.h.

Referenced by PANAnalysis(), and SetOscillationFunction().

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 EndJob(), PANAnalysis(), and SetFileNameTag().

float PANAnalysis::fRangeLowerLimitX [protected]

Definition at line 50 of file PANAnalysis.h.

Referenced by PANAnalysis(), RecoDist(), and SetHistBookInfo().

float PANAnalysis::fRangeLowerLimitY [protected]

Definition at line 53 of file PANAnalysis.h.

Referenced by PANAnalysis(), RecoDist(), and SetHistBookInfo().

float PANAnalysis::fRangeUpperLimitX [protected]

Definition at line 51 of file PANAnalysis.h.

Referenced by PANAnalysis(), RecoDist(), and SetHistBookInfo().

float PANAnalysis::fRangeUpperLimitY [protected]

Definition at line 54 of file PANAnalysis.h.

Referenced by PANAnalysis(), RecoDist(), and SetHistBookInfo().

TTree* PANAnalysis::fReweightTree [protected]
int* PANAnalysis::fVaryX [protected]

Definition at line 59 of file PANAnalysis.h.

Referenced by Do2DFit(), MakeReweightTree(), PANAnalysis(), VaryFitParam(), and ~PANAnalysis().


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

Generated on 22 Nov 2017 for loon by  doxygen 1.6.1