File indexing completed on 2024-04-06 12:21:37
0001
0002
0003
0004
0005
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;
0071 Int_t fCurrent;
0072
0073
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
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
0124 void deltaR(const double& eta1, double phi1,
0125 const double& eta2, double phi2,double& dr) const;
0126 void etaBin(const double&, int&) const;
0127 void etaValue(const int&, double&) const;
0128 void phiBin(double, int&) const;
0129 void phiValue(const int&, double&) const;
0130 double uniPhi(const double&) const;
0131
0132
0133 template<typename T>
0134 int find(const T&, const std::vector<T>&) const;
0135
0136
0137 std::vector<tpg> tpgsNear(const double& eta, const double& phi, const std::vector<tpg>&, const double& dr = .5) const;
0138
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
0142
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
0148 double avgPhi(const std::vector<tpg>&) const;
0149
0150 double avgEta(const std::vector<tpg>&) const;
0151
0152
0153 bool sanityCheck() const;
0154
0155
0156 void printCfFragment(std::ostream&) const;
0157
0158 std::vector<TObject*> hists_;
0159
0160
0161 void putHist(TObject* o) { hists_.push_back(o); }
0162
0163
0164 void makeCalibration();
0165
0166
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
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
0185
0186
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
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
0245
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
0267 if (!fChain) return 0;
0268 return fChain->GetEntry(entry);
0269 }
0270 Long64_t L1RCTCalibrator::LoadTree(Long64_t entry)
0271 {
0272
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
0288
0289
0290
0291
0292
0293
0294
0295
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
0311
0312
0313
0314
0315
0316 return kTRUE;
0317 }
0318
0319 void L1RCTCalibrator::Show(Long64_t entry)
0320 {
0321
0322
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