MadAnalysis.h

Go to the documentation of this file.

#ifndef madanalysis_h
#define madanalysis_h
#include <vector>
#include <string>
#include "TROOT.h"
#include "TChain.h"
#include "TH1.h"
#include "TH2.h"
#include "TF1.h"
#include "TVector3.h"
#include "Mad/MadQuantities.h"
#include "Mad/MadChi2Calc.h"

struct mcev{
  double TrueEnergy;
  double RecoEnergy;
  double MuonEnergy;
  double ShowerEnergy;
  double Weight;
  int Inu;
  int Iaction;
};

struct mcev_reweight{
  float TrueNuEnergy;
  float TrueMuEnergy;
  float TrueShwEnergy;
  float TrueMuDCosZ;
  float RecoNuEnergy;
  float RecoMuEnergy;
  float RecoShwEnergy;
  float RecoMuDCosZ;
  float Weight;
  bool PassTruth;
  int INu;
  int IAction;
  int Process;
  int InitState;
  int Nucleus;
  float KinX;
  float KinY;
  float KinQ2;
  float KinW2;
};

using namespace std;

class MadAnalysis : public MadQuantities
{

 protected:

  int NumberOfBins;        //number of bins to use in Reco/Bkgd spectra
  float RangeLowerLimit;   //lower range limit for Reco/Bkgd spectra
  float RangeUpperLimit;   //upper range limit for Reco/Bkgd spectra
  
  //two indices for two detectors: 0 = Near, 1 = Far.
  TH1F *MCSignalHist[2];   //reconstructed energy for MC sample signal events
  TH1F *MCBkgdHist[2];     //reconstructed energy for MC sample bkgd events
  TH1F *DataHist[2];       //reconstructed energy for data events passing cuts
  TH1F *DataSignalHist[2]; //reconstructed energy for data sample signal events
  float DataPOT[2];          //in units of e20
  int numDataEvts[2];        //actual number of signal snarls processed
  float signalFracInData[2]; //fraction of signal events passing analysis cuts
  TH1F *DataBkgdHist[2];     //reco energy for data sample background events
  
  vector<mcev_reweight> MCEvents[2]; //vector holding info for each MC event
  float MCPOT[2];                    //in units of e20
  int numMCEvts[2];         //actual number of MC snarls processed

  Int_t NgenBins[3];        // for doing systematics loops
  Double_t NgenMin[3];
  Double_t NgenMax[3];
  Double_t NgenMean[3];
  Double_t NgenErr[3];
  std::string NgenName[3];

  TF1 *OscillationFunction;
  int NumOscPars;
  double *OscillationParameters;
  int *varyX;  //tells the DoFit functions which params to vary (x,y,...)
  TH1F *BestFit[2];

  std::string tag; //for inserting into histo filename

  MadChi2Calc *Chi2Calc; //chi2 calculator object
  TH2F *Chi2Surf;
  float ChiMin;
  float NDOF;
  float Par1MinVal;
  float Par2MinVal;

  //Used for calculating systematic errors due to energy shifts
  float EShiftErr;  //error in energy scale (i.e. absolute calibration error)

  //define truth signal criteria here = events that will be oscillated
  virtual Bool_t PassTruthSignal(Int_t mcevent=0) { 
                                       if(!LoadTruth(mcevent)) return false;
                                       return true; }
  //put basic initial cuts in here for analysis
  virtual Bool_t PassBasicCuts(Int_t event=0) { if(event>=0) return true; 
                                                return false; }
  //put your favourite cuts in here for analysis
  virtual Bool_t PassAnalysisCuts(Int_t event=0) { if(event>=0) return true; 
                                                   return false; }
  //return a pid quantity
  virtual Float_t PID(Int_t event=0,Int_t method=0) { 
                                       if(method>=0 && event>=0) return 1.0;
                                       return 0.0; }
  //reconstruct nu energy using favourite method
  virtual Float_t RecoAnalysisEnergy(Int_t event=0) { 
                                       if(!LoadEvent(event)) return 0;
                                       return ntpEvent->ph.gev; }
  //allow event weighting:
  virtual Float_t GetWeight(Int_t event=0) {if(event>=0) return 1.0; return 0;}

  virtual Bool_t AddUserBranches(TTree *tree) {if(tree) return true; return false;}
  virtual void CallUserFunctions(Int_t) {return;}

  Bool_t writeOut; //write output file <=> do EndJob

 public:
 
  MadAnalysis(TChain *chainSR=0,TChain *chainMC=0,
              TChain *chainTH=0,TChain *chainEM=0);
  MadAnalysis(JobC *,string,int);
  virtual ~MadAnalysis();

  void ReInit(TChain *chainSR=0,TChain *chainMC=0,
              TChain *chainTH=0,TChain *chainEM=0);
  void ReInit(JobC *,string,int);

  void CreatePAN(std::string);
  void CreateANtpPAN(std::string);

  void SetFileNameTag(std::string);
  void SetHistBookInfo(int,float,float); //set Reco/Bkgd bins, hist range
  
  void RecoMCExperiment(int startnum,double NearPOT,double FarPOT);
  void RecoExperiment(int startnum,double NearExpectedNormToPOT,
                      double FarExpectedNormToPOT);
  void RecoMC(int startnum,double NearPOT,double FarPOT);
  
  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) {varyX[ix] = vx;}

  Bool_t Do2DFit(float n_a=90.,float min_a=0.5,float max_a=1.0,
                 float n_b=100.,float min_b=0.0,float max_b=0.005,
                 float n_c=1.,float f_c=0.);

  Bool_t Do2DFitNearFar(float n_a=90.,float min_a=0.5,float max_a=1.0,
                        float n_b=100.,float min_b=0.0,float max_b=0.005,
                        float n_c=1.,float f_c=0.);
  
  void SetNeugenReweightInfo(Int_t,Int_t*,Float_t*,Float_t*,
                             Float_t*,Float_t*,std::string*);
  
  MadChi2Calc *GetChi2Calc() {return Chi2Calc;} //so user can set options
  void SetEShiftErr(float err) {EShiftErr = err;}  

  int SetDataInfo(TH1F *hist=0, int numev=0, float pot=0, float sigfrac=0, 
                  int det=1);  
  TH1F *GetDataHist(int det) { if(det!=0&&det!=1) return 0; 
                               return DataHist[det]; }  
  int GetNumDataEvents(int det=1) { return numDataEvts[det]; }
  float GetDataPOT(int det=1) { return DataPOT[det]; }
  float GetSignalFracInData(int det=1) { return signalFracInData[det]; }
  
  void NoOutFile() { writeOut = false;}  //don't write output file
  void EndJob();

};
#endif // #ifdef madanalysis_h

---