MadMKAnalysis.h

Go to the documentation of this file.

//
// MadMKAnalysis
//
// Mike's MadAnalysis.  Mostly I just make a simple ntuple with it.
//
// Created:  M. Kordosky -- December, 2004
//
// $Author: gmieg $ 
//
// $Revision: 1.27 $
// 
// $Name:  $
//
// $Id: MadMKAnalysis.h,v 1.27 2009/02/28 21:46:13 gmieg Exp $
//
// $Log: MadMKAnalysis.h,v $
// Revision 1.27  2009/02/28 21:46:13  gmieg
// Update for GCC 4.4.0.  Mostly add #includes
// for reorganized C++ headers, such as -
// algorithm for: sort, find
// cassert for:   assert, abort
// climits for:   SHRT_MAX, INT_MAX
// cstdlib for:   abs, getenv, atoi, atof, strtod, exit
// cstring for:   strcmp
//
// Revision 1.26  2007/08/23 15:58:45  kordosky
// small updates, not committed before
//
// Revision 1.25  2007/06/05 17:10:41  kordosky
// blah
//
// Revision 1.24  2007/05/08 12:50:48  kordosky
// support intranuke reweighting
//
// Revision 1.23  2007/01/30 03:17:21  rhatcher
// No longer user BeamType "class" found in MCReweight, but rather the
// BeamType "namespace" found in Conventions package.
// Also covert to using "Detector::" rather than "DetectorType::", though
// for now they are synonyms.
//
// Revision 1.22  2006/10/16 19:53:30  kordosky
// cope with CEDAR strip indices set to -1. allow user to force beam type (useful for 150cm data). Open files with TFile::Open() to allow natural dcache access. IMPORTANT: Big change to pans. Tacking on a spill marker, looks like event=-1 with no muon or shower energy, pass and is_fid set to fail. Beam variables and nevent are correct for these entries. Modification intended to ease POT counting
//
// Revision 1.21  2006/03/09 19:33:30  kordosky
// Properly do energy corrections in MK and DpAnalysis
//
// Revision 1.20  2006/02/16 19:26:09  kordosky
// upgrades
//
// Revision 1.19  2006/01/26 23:41:34  kordosky
// upgrades.
//
// Revision 1.18  2006/01/10 17:40:50  kordosky
// Event display improvements. Now prints kinematic variables, indicates if muon stopped in detector and applies energy corrections.
//
// Revision 1.17  2005/12/14 20:35:59  kordosky
// Updated energy corrections.
//
// Revision 1.16  2005/11/17 15:50:17  vahle
// Added checks for track vertex in fid volume.  Switched goodbeam varible to use BMSpillAna
//
// Revision 1.15  2005/11/17 13:22:05  vahle
// Added a count of the number of reconstructed strips and the ph of reconstructed strips that show up in another event
//
// Revision 1.14  2005/11/14 12:04:52  vahle
// Adding Brians lowph ratio variable for cutting fd junk events to MadMKAnalysis ntuple
//
// Revision 1.13  2005/11/04 15:27:55  vahle
// Remove all vestiges of old beam monitoring ntuple lookup, change the way track containment is decided in the FD
//
// Revision 1.12  2005/10/07 18:32:54  vahle
// Added late activity variables to MadMKAnalysis
//
// Revision 1.11  2005/10/07 18:05:28  kordosky
// a few changes to analysis, refinement of pdfs. MadBase not loading mcHeader correctly when GetEntry() called.  Suprisingly, it had no effect up until now.
//
// Revision 1.10  2005/10/06 17:29:22  vahle
// Added early activity variables to MadMKAnalysis
//
// Revision 1.9  2005/10/06 15:03:26  kordosky
// various improvements, pids finally integrated into pan construction routine, still not fully debugged though
//
// Revision 1.8  2005/09/23 11:48:22  kordosky
//  a few additions for R1.18. code meant to handle R1.16 and R1.18 ntpshowers not working properly. only handles R1.18 correctly.
//
// Revision 1.7  2005/09/14 13:32:45  kordosky
// DpID and NsID classes now fully working and validated with MadTestAnalysis by comparing nannpid vs. annpid and npid0 vs. pid0.  pdf files and weights stored in data subdirectory. FD pdfs needed for DpID.
//
//

#ifndef MADMKANALYSIS
#define MADMKANALYSIS

#include <string>

#include "Mad/MadAnalysis.h"
#include "Conventions/PlaneCoverage.h"
#include "Conventions/PlaneView.h"
#include "Conventions/Detector.h"
#include "Registry/Registry.h"

#include "Conventions/BeamType.h"

#include "Mad/MadNsID.h"
#include "Mad/MadDpID.h"
#include "Mad/MadQEID.h"

class TFile;

class MadMKAnalysis : public MadAnalysis {

 public:  
  void CreatePAN(std::string s, const char* bmonpath="");
  MadMKAnalysis(TChain *chainSR=0,TChain *chainMC=0,
              TChain *chainTH=0,TChain *chainEM=0);
  MadMKAnalysis(JobC *,string,int);
  virtual ~MadMKAnalysis();

  Float_t ClosestPointToLine(const Float_t& x1, 
                             const Float_t& y1,
                             const Float_t& x2,
                             const Float_t& y2,
                             const Float_t& xpt,
                             const Float_t& ypt,
                             Float_t& x, Float_t& y);

  // fiducial volume and track containment as defined by pitt group
  static Bool_t PittInFidVol(Float_t x, Float_t y, 
                             Float_t z, Float_t u, Float_t v);
  static Int_t PittTrkContained(Float_t x, Float_t y, 
                                Float_t z, Float_t u, Float_t v);

  // is this strip in the partially instrumented region of the ND?
  static Int_t InPartialRegion(UShort_t plane, UShort_t strip);
  
  // pass back an bitfield with bits corresponding to radial cuts
  // Float_t rings[5]={0.1,0.25,0.5,1.0,1.5}; // in meters
  static Int_t NDRadialFidVol(Float_t x,Float_t y, Float_t z);

  void ForceMCBeam(const BeamType::BeamType_t& beam) {fMCBeam=beam;}
  const BeamType::BeamType_t& GetForcedMCBeam() {return fMCBeam;}
  void SetDataUseMCBeam(bool x) {fDataUseMCBeam=x;}
  void OnlyWriteFidEvents(bool x) {fOnlyWriteFidEvents=x;}
  void FillMCHists(TFile* fout);

  static int FarTrkContained(Float_t x, Float_t y, Float_t z);

  void DoInukeReweight(bool t){fDoInukeReweight=t;}
  void SetNRandomSets(int i){fNRandomSets=i;}
  // PID calculator objects
  MadNsID nsid;
  MadDpID dpid;
  MadQEID qeid;
  
 protected:
  // signal inside and outside of the partial plane region
  // for whole event and just for track portion
  // uses currently loaded event and track  
  void SigInOut(Float_t& sigfull, Float_t& sigpart, Float_t& trkfull, Float_t& trkpart);
  void SigInOut(Int_t trkidx, Float_t& sigfull, Float_t& sigpart, Float_t& trkfull, Float_t& trkpart);


  // from base class
//  Float_t RecoAnalysisEnergy(Int_t event);
  Float_t RecoShwEnergy(int ishw, const Detector::Detector_t& det, int mode=1, bool isdata=true);
  Float_t RecoMKMuEnergy(Int_t& opt, Int_t itrk, bool isdata=true, Detector::Detector_t det=Detector::kNear);

  //  virtual Bool_t PassBasicCuts(Int_t event);

  // muon angle with respect to neutrino angle
  // it's detector dependant!
  virtual Float_t RecoMuDCosNeuND(Int_t itr=0, Float_t* vtx=0);
  virtual Float_t RecoMuDCosNeuFD(Int_t itr=0, Float_t* vtx=0);
  
  //  virtual Float_t RecoKinematics(Int_t event);

  virtual bool InFidVol(Int_t evidx);
  virtual bool InFidVol(); // assumes event already loaded,
                           //and that we're interested in event vertex
  bool InFidVol(float vx, float vy, float vz); //looks to see if vx,vy,vz 
                                               //is in fid. vol

  void SetBECFile(const char* f);
  Registry& GetBECRegistry() { return fBECConfig; }

  // number of planes in the reconstructed object
  // possibly with a ph cut in pes
  // return the sigcor, after the ph cut, in sumph
  Int_t NPlanesInObj(TObject *rec, TObject *obj, float phcut, float& sumph);
  Int_t NStripsInObj(TObject *rec, TObject *obj, float phcut, float& sumph);
  //GetEvtTimeChar assumesevent is already loaded
  void GetEvtTimeChar(double &timemin, double &timemax, double &timeln);
  void GetEvtTimeCharPHC(double &timemin, double &timemax, double &timeln);

  float MuonShowerEnergy(const NtpSREvent* evt, const NtpSRTrack* trk);
  
  void EarlyActivity(int evidx, float &earlywph, float &earlywphpw, 
                     float &earlywphsphere, float &ph1000, 
                     float &ph500, float &ph100,
                     float &lateph1000, 
                     float &lateph500, float &lateph100,
                     float &lateph50);

  //a function to count the number of reconstructed strips 
  //that are in this event and in some other event
  //(by reconstructed strips I mean CandStripSR)
  void DupRecoStrips(int evidx, int &nduprs, float &dupphrs);
  
  float ComputeLowPHRatio();


 private:
  
  // Z (m) locations in the ND (for fiducial cuts)
  static const Float_t kVetoEndZ;
  static const Float_t kTargEndZ;
  static const Float_t kCalorEndZ;
  static const Float_t kHalfCalorEndZ;
  static const Float_t kSpecEndZ;
  // xy beam center in ND
  static const Float_t kXcenter;
  static const Float_t kYcenter;

  //constants for EarlyActivity estimation
  static const Double_t fEarlyActivityWindowSize;
  static const Double_t fPMTTimeConstant;
  static const Int_t fNPlaneEarlyAct;
  static const Double_t fEarlyPlaneSphere;


  //  std::string fBECFile; // file for beam energy reweighting
  Registry fBECConfig; // config for beam energy reweighting
  BeamType::BeamType_t fMCBeam;
  bool fDataUseMCBeam;
  bool fDoInukeReweight;
  bool fOnlyWriteFidEvents;//only write out events reco'd inside fid. vol.
  int fNRandomSets; // number of random sets to generate 0= 1sigmas
};

#endif

---