test/Calib_RootScripts/L1RCTCalibrator.h

0001 //////////////////////////////////////////////////////////
0002 // This class has been automatically generated on
0003 // Mon Sep 29 12:19:26 2008 by ROOT version 5.18/00a
0004 // from TTree L1RCTCalibrator/RCT Calibration Tree
0005 // found on file: /scratch/lgray/PGun2pi2gamma_calibration/rctCalibratorFarmout_cfg-merge_cfg-PGun2pi2gamma-0034.root
0006 //////////////////////////////////////////////////////////
0007
0008 #ifndef L1RCTCalibrator_h
0009 #define L1RCTCalibrator_h
0010
0011 #include <TROOT.h>
0012 #include <TChain.h>
0013 #include <TFile.h>
0014 #include <TH2F.h>
0015 #include <TGraphAsymmErrors.h>
0016 #include <TGraph2DErrors.h>
0017 #include <iostream>
0018 #include <vector>
0019 #include <cmath>
0020 #include <TMath.h>
0021 #include <TF1.h>
0022 #include <TF2.h>
0023
0024 class generator
0025 {
0026  public:
0027   int particle_type;
0028   double et, phi, eta;
0029
0030   bool operator==(const generator& r) const { return ( particle_type == r.particle_type && et == r.et &&
0031                                phi == r.phi && eta == r.eta ); }
0032 };
0033
0034 class region
0035 {
0036  public:
0037   int ieta, iphi, linear_et;
0038   double eta,phi;
0039
0040   bool operator==(const region& r) const { return ((ieta == r.ieta) && (iphi == r.iphi)); }
0041 };
0042
0043 class tpg
0044 {
0045  public:
0046   int ieta, iphi;
0047   double ecalEt, hcalEt, ecalE, hcalE;
0048   double eta,phi;
0049
0050   bool operator==(const tpg& r) const { return ((ieta == r.ieta) && (iphi == r.iphi)); }
0051 };
0052
0053 class event_data
0054 {
0055  public:
0056   int event, run;
0057   std::vector<generator> gen_particles;
0058   std::vector<region> regions;
0059   std::vector<tpg> tpgs;
0060 };
0061
0062 class L1RCTCalibrator {
0063 public :
0064
0065   typedef TH1F* TH1Fptr;
0066   typedef TH2F* TH2Fptr;
0067   typedef TGraphAsymmErrors* TGraphptr;
0068   typedef TGraph2DErrors* TGraph2Dptr;
0069
0070    TTree          *fChain;   //!pointer to the analyzed TTree or TChain
0071    Int_t           fCurrent; //!current Tree number in a TChain
0072
0073    // Declaration of leaf types
0074    UInt_t          Event_event;
0075    UInt_t          Event_run;
0076    UInt_t          Generator_nGen;
0077    Int_t           Generator_particle_type[100];
0078    Double_t        Generator_et[100];
0079    Double_t        Generator_eta[100];
0080    Double_t        Generator_phi[100];
0081    UInt_t          Generator_crate[100];
0082    UInt_t          Generator_card[100];
0083    UInt_t          Generator_region[100];
0084    UInt_t          Region_nRegions;
0085    Int_t           Region_linear_et[200];
0086    Int_t           Region_ieta[200];
0087    Int_t           Region_iphi[200];
0088    Double_t        Region_eta[200];
0089    Double_t        Region_phi[200];
0090    UInt_t          Region_crate[200];
0091    UInt_t          Region_card[200];
0092    UInt_t          Region_region[200];
0093    UInt_t          CaloTPG_nTPG;
0094    Int_t           CaloTPG_ieta[3100];
0095    Int_t           CaloTPG_iphi[3100];
0096    Double_t        CaloTPG_eta[3100];
0097    Double_t        CaloTPG_phi[3100];
0098    Double_t        CaloTPG_ecalEt[3100];
0099    Double_t        CaloTPG_hcalEt[3100];
0100    Double_t        CaloTPG_ecalE[3100];
0101    Double_t        CaloTPG_hcalE[3100];
0102    UInt_t          CaloTPG_crate[3100];
0103    UInt_t          CaloTPG_card[3100];
0104    UInt_t          CaloTPG_region[3100];
0105
0106    // List of branches
0107    TBranch        *b_Event;   //!
0108    TBranch        *b_Generator;   //!
0109    TBranch        *b_Region;   //!
0110    TBranch        *b_CaloTPG;   //!
0111
0112    L1RCTCalibrator(TTree *tree=0);
0113    virtual ~L1RCTCalibrator();
0114    virtual Int_t    GetEntry(Long64_t entry);
0115    virtual Long64_t LoadTree(Long64_t entry);
0116    virtual void     Init(TTree *tree);
0117    virtual void     Loop();
0118    virtual Bool_t   Notify();
0119    virtual void     Show(Long64_t entry = -1);
0120    void BookHistos();
0121    void WriteHistos();
0122
0123    // ----------protected member functions---------------
0124   void deltaR(const double& eta1, double phi1,
0125           const double& eta2, double phi2,double& dr) const; // calculates delta R between two coordinates
0126   void etaBin(const double&, int&) const; // calculates Trigger Tower number
0127   void etaValue(const int&, double&) const; // calculates Trigger Tower eta bin center
0128   void phiBin(double, int&) const; // calculates TT phi bin
0129   void phiValue(const int&, double&) const; // calculates TT phi bin center
0130   double uniPhi(const double&) const; // returns phi that is in [0, 2*pi]
0131
0132   // returns -1 if item not present, [index] if it is, T must have an == operator
0133   template<typename T>
0134     int find(const T&, const std::vector<T>&) const;
0135
0136   // returns tpgs within a specified delta r near the point (eta,phi)
0137   std::vector<tpg> tpgsNear(const double& eta, const double& phi, const std::vector<tpg>&, const double& dr = .5) const;
0138   // returns an ordered pair of (Et, deltaR)
0139   std::pair<double,double> showerSize(const std::vector<tpg>&, const double frac = .95, const double& max_dr = .5,
0140                       const bool& ecal = true, const bool& hcal = true) const;
0141   // returns the sum of tpg Et near the point (eta,phi) within a specified delta R,
0142   // can choose to only give ecal or hcal sum through bools
0143   double sumEt(const double& eta, const double& phi, const std::vector<region>&, const double& dr = .5) const;
0144   double sumEt(const double& eta, const double& phi, const std::vector<tpg>&, const double& dr = .5,
0145            const bool& ecal = true, const bool& hcal = true, const bool& apply_corrections = false,
0146            const double& high_low_crossover = 23) const;
0147   // returns energy weighted average of Eta
0148   double avgPhi(const std::vector<tpg>&) const;
0149   // returns energy weighted average of Phi
0150   double avgEta(const std::vector<tpg>&) const;
0151
0152
0153    bool sanityCheck() const;
0154    // prints a CFG or python language configuration file fragment that contains the
0155   // resultant calibration information
0156   void printCfFragment(std::ostream&) const;
0157
0158   std::vector<TObject*> hists_;
0159
0160   // saves pointer to Histogram in a vector, making writing out easier later.
0161   void putHist(TObject* o) { hists_.push_back(o); }
0162
0163   //do the calibration
0164   void makeCalibration();
0165
0166   //finds overlapping particles
0167   std::vector<generator> overlaps(const std::vector<generator>& v) const;
0168
0169   bool python_;
0170   const std::string& fitOpts() const { return fitOpts_; }
0171   const int& debug() const { return debug_; }
0172
0173   int debug_;
0174   std::string fitOpts_;
0175   TFile* output_;
0176   std::vector<event_data> data_;
0177
0178   const double deltaEtaBarrel_, maxEtaBarrel_, deltaPhi_;
0179   std::vector<double> endcapEta_;
0180
0181   //The final output... correction factors
0182   double ecal_[28][3], hcal_[28][3], hcal_high_[28][3], cross_[28][6], he_low_smear_[28], he_high_smear_[28];
0183
0184   // histograms
0185
0186   // diagnostic histograms
0187   TH1Fptr hEvent, hRun, hGenPhi, hGenEta, hGenEt, hGenEtSel,
0188     hRCTRegionEt, hRCTRegionPhi, hRCTRegionEta,
0189     hTpgSumEt, hTpgSumEta, hTpgSumPhi;
0190
0191   TH2Fptr hGenPhivsTpgSumPhi, hGenEtavsTpgSumEta, hGenPhivsRegionPhi, hGenEtavsRegionEta;
0192
0193   TH1Fptr hDeltaEtPeakvsEtaBin_uc[12], hDeltaEtPeakvsEtaBin_c[12], hDeltaEtPeakRatiovsEtaBin[12],
0194     hPhotonDeltaEtPeakvsEtaBinAllEt_uc, hPhotonDeltaEtPeakvsEtaBinAllEt_c, hPhotonDeltaEtPeakRatiovsEtaBinAllEt,
0195     hPionDeltaEtPeakvsEtaBinAllEt_uc, hPionDeltaEtPeakvsEtaBinAllEt_c, hPionDeltaEtPeakRatiovsEtaBinAllEt;
0196
0197   TH1Fptr hPhotonDeltaR95[28], hNIPionDeltaR95[28], hPionDeltaR95[28] ;
0198
0199   // histograms for algorithm
0200   TGraphptr gPhotonEtvsGenEt[28], gNIPionEtvsGenEt[28];
0201   TGraph2Dptr gPionEcalEtvsHcalEtvsGenEt[28];
0202
0203   TH1Fptr hPhotonDeltaEOverE_uncor[28], hPionDeltaEOverE_uncor[28];
0204   TH1Fptr hPhotonDeltaEOverE[28], hPionDeltaEOverE[28];
0205   TH1Fptr hPhotonDeltaEOverE_cor[28], hPionDeltaEOverE_cor[28];
0206 };
0207
0208 #endif
0209
0210 #ifdef L1RCTCalibrator_cxx
0211 L1RCTCalibrator::L1RCTCalibrator(TTree *tree):deltaEtaBarrel_(0.0870), maxEtaBarrel_(20*deltaEtaBarrel_),deltaPhi_(0.0870)
0212 {
0213   python_ = true;
0214   debug_ = 0;
0215   fitOpts_ = ((debug_ > -1) ? "ELMRNF" : "QELMRNF");
0216   output_ = new TFile("L1RCTCalibration_info.root","RECREATE");
0217
0218   endcapEta_.push_back(0.09);
0219   endcapEta_.push_back(0.1);
0220   endcapEta_.push_back(0.113);
0221   endcapEta_.push_back(0.129);
0222   endcapEta_.push_back(0.15);
0223   endcapEta_.push_back(0.178);
0224   endcapEta_.push_back(0.15);
0225   endcapEta_.push_back(0.35);
0226
0227   for(int i = 0; i < 28; ++i)
0228     {
0229       he_low_smear_[i]  = -999;
0230       he_high_smear_[i] = -999;
0231       for(int j = 0; j < 6; ++j)
0232         {
0233           if(j < 3)
0234             {
0235               ecal_[i][j]      = -999;
0236               hcal_[i][j]      = -999;
0237               hcal_high_[i][j] = -999;
0238             }
0239           cross_[i][j] = -999;
0240         }
0241     }
0242
0243
0244 // if parameter tree is not specified (or zero), connect the file
0245 // used to generate this class and read the Tree.
0246    if (tree == 0) {
0247       TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("/scratch/lgray/PGun2pi2gamma_calibration/rctCalibratorFarmout_cfg-merge_cfg-PGun2pi2gamma-0034.root");
0248       if (!f) {
0249          f = new TFile("/scratch/lgray/PGun2pi2gamma_calibration/rctCalibratorFarmout_cfg-merge_cfg-PGun2pi2gamma-0034.root");
0250       }
0251       tree = (TTree*)gDirectory->Get("L1RCTCalibrator");
0252
0253    }
0254    Init(tree);
0255    BookHistos();
0256 }
0257
0258 L1RCTCalibrator::~L1RCTCalibrator()
0259 {
0260    if (!fChain) return;
0261    delete fChain->GetCurrentFile();
0262 }
0263
0264 Int_t L1RCTCalibrator::GetEntry(Long64_t entry)
0265 {
0266 // Read contents of entry.
0267    if (!fChain) return 0;
0268    return fChain->GetEntry(entry);
0269 }
0270 Long64_t L1RCTCalibrator::LoadTree(Long64_t entry)
0271 {
0272 // Set the environment to read one entry
0273    if (!fChain) return -5;
0274    Long64_t centry = fChain->LoadTree(entry);
0275    if (centry < 0) return centry;
0276    if (!fChain->InheritsFrom(TChain::Class()))  return centry;
0277    TChain *chain = (TChain*)fChain;
0278    if (chain->GetTreeNumber() != fCurrent) {
0279       fCurrent = chain->GetTreeNumber();
0280       Notify();
0281    }
0282    return centry;
0283 }
0284
0285 void L1RCTCalibrator::Init(TTree *tree)
0286 {
0287    // The Init() function is called when the selector needs to initialize
0288    // a new tree or chain. Typically here the branch addresses and branch
0289    // pointers of the tree will be set.
0290    // It is normally not necessary to make changes to the generated
0291    // code, but the routine can be extended by the user if needed.
0292    // Init() will be called many times when running on PROOF
0293    // (once per file to be processed).
0294
0295    // Set branch addresses and branch pointers
0296    if (!tree) return;
0297    fChain = tree;
0298    fCurrent = -1;
0299    fChain->SetMakeClass(1);
0300
0301    fChain->SetBranchAddress("Event", &Event_event, &b_Event);
0302    fChain->SetBranchAddress("Generator", &Generator_nGen, &b_Generator);
0303    fChain->SetBranchAddress("Region", &Region_nRegions, &b_Region);
0304    fChain->SetBranchAddress("CaloTPG", &CaloTPG_nTPG, &b_CaloTPG);
0305    Notify();
0306 }
0307
0308 Bool_t L1RCTCalibrator::Notify()
0309 {
0310    // The Notify() function is called when a new file is opened. This
0311    // can be either for a new TTree in a TChain or when when a new TTree
0312    // is started when using PROOF. It is normally not necessary to make changes
0313    // to the generated code, but the routine can be extended by the
0314    // user if needed. The return value is currently not used.
0315
0316    return kTRUE;
0317 }
0318
0319 void L1RCTCalibrator::Show(Long64_t entry)
0320 {
0321 // Print contents of entry.
0322 // If entry is not specified, print current entry
0323    if (!fChain) return;
0324    fChain->Show(entry);
0325 }
0326
0327 template<typename T>
0328 int L1RCTCalibrator::find(const T& item, const std::vector<T>& v) const
0329 {
0330   for(unsigned i = 0; i < v.size(); ++i)
0331     if(item == v[i]) return i;
0332   return -1;
0333 }
0334
0335 #endif // #ifdef L1RCTCalibrator_cxx