File indexing completed on 2024-04-06 11:59:11
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 #include <algorithm>
0021 #include <cmath>
0022 #include <iostream>
0023 #include <vector>
0024 #include <cstring>
0025 #include <TCanvas.h>
0026 #include <TChain.h>
0027 #include <TFile.h>
0028 #include <TH1.h>
0029 #include <TH1D.h>
0030 #include <TH2.h>
0031 #include <TH2D.h>
0032 #include <TProfile.h>
0033 #include <TROOT.h>
0034 #include <TStyle.h>
0035 #include <TTree.h>
0036
0037 class HBHEMuonOfflineSimAnalyzer {
0038 private:
0039 TTree *fChain;
0040 Int_t fCurrent;
0041
0042 UInt_t Run_No;
0043 UInt_t Event_No;
0044 UInt_t LumiNumber;
0045 UInt_t BXNumber;
0046 double pt_of_muon;
0047 double eta_of_muon;
0048 double phi_of_muon;
0049 double p_of_muon;
0050 double ecal_3x3;
0051 unsigned int ecal_detID;
0052 double hcal_1x1;
0053 double matchedId;
0054 unsigned int hcal_detID;
0055 unsigned int hcal_cellHot;
0056 double activeLength;
0057 double hcal_edepth1;
0058 double hcal_edepth2;
0059 double hcal_edepth3;
0060 double hcal_edepth4;
0061 double hcal_activeL1;
0062 double hcal_activeL2;
0063 double hcal_activeL3;
0064 double hcal_activeL4;
0065 double activeLengthHot;
0066 double hcal_edepthHot1;
0067 double hcal_edepthHot2;
0068 double hcal_edepthHot3;
0069 double hcal_edepthHot4;
0070 double hcal_activeHotL1;
0071 double hcal_activeHotL2;
0072 double hcal_activeHotL3;
0073 double hcal_activeHotL4;
0074 double hcal_edepth5;
0075 double hcal_activeL5;
0076 double hcal_edepthHot5;
0077 double hcal_activeHotL5;
0078 double hcal_edepth6;
0079 double hcal_activeL6;
0080 double hcal_edepthHot6;
0081 double hcal_activeHotL6;
0082 double hcal_edepth7;
0083 double hcal_activeL7;
0084 double hcal_edepthHot7;
0085 double hcal_activeHotL7;
0086
0087 TBranch *b_Run_No;
0088 TBranch *b_Event_No;
0089 TBranch *b_LumiNumber;
0090 TBranch *b_BXNumber;
0091 TBranch *b_pt_of_muon;
0092 TBranch *b_eta_of_muon;
0093 TBranch *b_phi_of_muon;
0094 TBranch *b_p_of_muon;
0095 TBranch *b_ecal_3x3;
0096 TBranch *b_ecal_detID;
0097 TBranch *b_hcal_1x1;
0098 TBranch *b_hcal_detID;
0099 TBranch *b_hcal_cellHot;
0100 TBranch *b_activeLength;
0101 TBranch *b_hcal_edepth1;
0102 TBranch *b_hcal_edepth2;
0103 TBranch *b_hcal_edepth3;
0104 TBranch *b_hcal_edepth4;
0105 TBranch *b_hcal_activeL1;
0106 TBranch *b_hcal_activeL2;
0107 TBranch *b_hcal_activeL3;
0108 TBranch *b_hcal_activeL4;
0109 TBranch *b_activeLengthHot;
0110 TBranch *b_hcal_edepthHot1;
0111 TBranch *b_hcal_edepthHot2;
0112 TBranch *b_hcal_edepthHot3;
0113 TBranch *b_hcal_edepthHot4;
0114 TBranch *b_hcal_activeHotL1;
0115 TBranch *b_hcal_activeHotL2;
0116 TBranch *b_hcal_activeHotL3;
0117 TBranch *b_hcal_activeHotL4;
0118 TBranch *b_hcal_edepth5;
0119 TBranch *b_hcal_activeL5;
0120 TBranch *b_hcal_edepthHot5;
0121 TBranch *b_hcal_activeHotL5;
0122 TBranch *b_hcal_edepth6;
0123 TBranch *b_hcal_activeL6;
0124 TBranch *b_hcal_edepthHot6;
0125 TBranch *b_hcal_activeHotL6;
0126 TBranch *b_hcal_edepth7;
0127 TBranch *b_hcal_activeL7;
0128 TBranch *b_hcal_edepthHot7;
0129 TBranch *b_hcal_activeHotL7;
0130 TBranch *b_matchedId;
0131
0132 public:
0133 HBHEMuonOfflineSimAnalyzer(const char *infile,
0134 const char *outfile = "dyll_PU20_25_output.root",
0135 const int mode = 0,
0136 const int maxDHB = 4,
0137 const int maxDHE = 7);
0138 virtual ~HBHEMuonOfflineSimAnalyzer();
0139 virtual Int_t Cut(Long64_t entry);
0140 virtual Int_t GetEntry(Long64_t entry);
0141 virtual Long64_t LoadTree(Long64_t entry);
0142 virtual void Init(TTree *tree);
0143 virtual void Loop();
0144 virtual Bool_t Notify();
0145 virtual void Show(Long64_t entry = -1);
0146
0147 std::vector<std::string> firedTriggers;
0148 void BookHistograms(const char *);
0149 void WriteHistograms();
0150 bool LooseMuon();
0151 bool tightMuon();
0152 bool SoftMuon();
0153 void etaPhiHcal(unsigned int detId, int &eta, int &phi, int &depth);
0154 void etaPhiEcal(unsigned int detId, int &type, int &zside, int &etaX, int &phiY, int &plane, int &strip);
0155 void calculateP(double pt, double eta, double &pM);
0156 void close();
0157 int NDepthBins(int ieta, int iphi);
0158 int NPhiBins(int ieta);
0159
0160 private:
0161 static const bool debug_ = false;
0162 static const int maxDep = 7;
0163 static const int maxEta = 29;
0164 static const int maxPhi = 72;
0165
0166 static const int maxHist = 20000;
0167 int modeLHC_, maxDepthHB_, maxDepthHE_, maxDepth_;
0168 int nHist, nDepths[maxEta], nDepthsPhi[maxEta], indxEta[maxEta][maxDep][maxPhi];
0169 TFile *output_file;
0170
0171 TH1D *h_Pt_Muon[3], *h_Eta_Muon[3], *h_Phi_Muon[3], *h_P_Muon[3];
0172 TH1D *h_PF_Muon[3], *h_GlobTrack_Chi[3], *h_Global_Muon_Hits[3];
0173 TH1D *h_MatchedStations[3], *h_Tight_TransImpactparameter[3];
0174 TH1D *h_Tight_LongitudinalImpactparameter[3], *h_InnerTrackPixelHits[3];
0175 TH1D *h_TrackerLayer[3], *h_IsolationR04[3], *h_Global_Muon[3];
0176 TH1D *h_LongImpactParameter[3], *h_LongImpactParameterBin1[3], *h_LongImpactParameterBin2[3];
0177
0178 TH1D *h_TransImpactParameter[3], *h_TransImpactParameterBin1[3], *h_TransImpactParameterBin2[3];
0179 TH1D *h_Hot_MuonEnergy_hcal_ClosestCell[3][maxHist], *h_Hot_MuonEnergy_hcal_HotCell[3][maxHist],
0180 *h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[3][maxHist], *h_HotCell_MuonEnergy_phi[3][maxHist];
0181 TH2D *h_2D_Bin1[3], *h_2D_Bin2[3];
0182 TH1D *h_ecal_energy[3], *h_hcal_energy[3], *h_3x3_ecal[3], *h_1x1_hcal[3];
0183 TH1D *h_MuonHittingEcal[3], *h_HotCell[3], *h_MuonEnergy_hcal[3][maxHist];
0184 TH1D *h_Hot_MuonEnergy_hcal[3][maxHist];
0185 TH2D *hcal_ietaVsEnergy[3];
0186 TProfile *h_EtaX_hcal[3], *h_PhiY_hcal[3], *h_EtaX_ecal[3], *h_PhiY_ecal[3];
0187 TProfile *h_Eta_ecal[3], *h_Phi_ecal[3];
0188 TProfile *h_MuonEnergy_eta[3][maxDep], *h_MuonEnergy_phi[3][maxDep], *h_MuonEnergy_muon_eta[3][maxDep];
0189 TProfile *h_Hot_MuonEnergy_eta[3][maxDep], *h_Hot_MuonEnergy_phi[3][maxDep], *h_Hot_MuonEnergy_muon_eta[3][maxDep];
0190 TProfile *h_IsoHot_MuonEnergy_eta[3][maxDep], *h_IsoHot_MuonEnergy_phi[3][maxDep],
0191 *h_IsoHot_MuonEnergy_muon_eta[3][maxDep];
0192 TProfile *h_IsoWithoutHot_MuonEnergy_eta[3][maxDep], *h_IsoWithoutHot_MuonEnergy_phi[3][maxDep],
0193 *h_IsoWithoutHot_MuonEnergy_muon_eta[3][maxDep];
0194 TProfile *h_HotWithoutIso_MuonEnergy_eta[3][maxDep], *h_HotWithoutIso_MuonEnergy_phi[3][maxDep],
0195 *h_HotWithoutIso_MuonEnergy_muon_eta[3][maxDep];
0196 };
0197
0198 HBHEMuonOfflineSimAnalyzer::HBHEMuonOfflineSimAnalyzer(
0199 const char *infile, const char *outFileName, const int mode, const int maxDHB, const int maxDHE) {
0200 modeLHC_ = mode;
0201 maxDepthHB_ = maxDHB;
0202 maxDepthHE_ = maxDHE;
0203 maxDepth_ = (maxDepthHB_ > maxDepthHE_) ? maxDepthHB_ : maxDepthHE_;
0204
0205
0206 TFile *f = new TFile(infile);
0207 TDirectory *dir = (TDirectory *)f->Get("HcalHBHEMuonAnalyzer");
0208 TTree *tree(0);
0209 dir->GetObject("TREE", tree);
0210 Init(tree);
0211
0212
0213 BookHistograms(outFileName);
0214 }
0215
0216 HBHEMuonOfflineSimAnalyzer::~HBHEMuonOfflineSimAnalyzer() {
0217 if (!fChain)
0218 return;
0219 delete fChain->GetCurrentFile();
0220 }
0221
0222 Int_t HBHEMuonOfflineSimAnalyzer::Cut(Long64_t) {
0223
0224
0225
0226 return 1;
0227 }
0228
0229 Int_t HBHEMuonOfflineSimAnalyzer::GetEntry(Long64_t entry) {
0230
0231 if (!fChain)
0232 return 0;
0233 return fChain->GetEntry(entry);
0234 }
0235
0236 Long64_t HBHEMuonOfflineSimAnalyzer::LoadTree(Long64_t entry) {
0237
0238 if (!fChain)
0239 return -5;
0240 Long64_t centry = fChain->LoadTree(entry);
0241 if (centry < 0)
0242 return centry;
0243 if (fChain->GetTreeNumber() != fCurrent) {
0244 fCurrent = fChain->GetTreeNumber();
0245 Notify();
0246 }
0247 return centry;
0248 }
0249
0250 void HBHEMuonOfflineSimAnalyzer::Init(TTree *tree) {
0251
0252
0253
0254
0255
0256
0257
0258
0259
0260 pt_of_muon = 0;
0261 eta_of_muon = 0;
0262 phi_of_muon = 0;
0263 p_of_muon = 0;
0264 ecal_3x3 = 0;
0265 ecal_detID = 0;
0266 hcal_1x1 = 0;
0267 hcal_detID = 0;
0268 hcal_cellHot = 0;
0269 activeLength = 0;
0270 hcal_edepth1 = 0;
0271 hcal_edepth2 = 0;
0272 hcal_edepth3 = 0;
0273 hcal_edepth4 = 0;
0274 hcal_activeL1 = 0;
0275 hcal_activeL2 = 0;
0276 hcal_activeL3 = 0;
0277 hcal_activeL4 = 0;
0278 activeLengthHot = 0;
0279 hcal_edepthHot1 = 0;
0280 hcal_edepthHot2 = 0;
0281 hcal_edepthHot3 = 0;
0282 hcal_edepthHot4 = 0;
0283 hcal_activeHotL1 = 0;
0284 hcal_activeHotL2 = 0;
0285 hcal_activeHotL3 = 0;
0286 hcal_activeHotL4 = 0;
0287 hcal_edepth5 = 0;
0288 hcal_activeL5 = 0;
0289 hcal_edepthHot5 = 0;
0290 hcal_activeHotL5 = 0;
0291 hcal_edepth6 = 0;
0292 hcal_activeL6 = 0;
0293 hcal_edepthHot6 = 0;
0294 hcal_activeHotL6 = 0;
0295 hcal_edepth7 = 0;
0296 hcal_activeL7 = 0;
0297 hcal_edepthHot7 = 0;
0298 hcal_activeHotL7 = 0;
0299 matchedId = 0;
0300 if (!tree)
0301 return;
0302 fChain = tree;
0303 fCurrent = -1;
0304 fChain->SetMakeClass(1);
0305
0306 fChain->SetBranchAddress("Run_No", &Run_No, &b_Run_No);
0307 fChain->SetBranchAddress("Event_No", &Event_No, &b_Event_No);
0308 fChain->SetBranchAddress("LumiNumber", &LumiNumber, &b_LumiNumber);
0309 fChain->SetBranchAddress("BXNumber", &BXNumber, &b_BXNumber);
0310 fChain->SetBranchAddress("pt_of_muon", &pt_of_muon, &b_pt_of_muon);
0311 fChain->SetBranchAddress("eta_of_muon", &eta_of_muon, &b_eta_of_muon);
0312 fChain->SetBranchAddress("phi_of_muon", &phi_of_muon, &b_phi_of_muon);
0313 fChain->SetBranchAddress("p_of_muon", &p_of_muon, &b_p_of_muon);
0314 fChain->SetBranchAddress("ecal_3x3", &ecal_3x3, &b_ecal_3x3);
0315 fChain->SetBranchAddress("ecal_detID", &ecal_detID, &b_ecal_detID);
0316 fChain->SetBranchAddress("hcal_1x1", &hcal_1x1, &b_hcal_1x1);
0317 fChain->SetBranchAddress("matchedId", &matchedId, &b_matchedId);
0318 fChain->SetBranchAddress("hcal_detID", &hcal_detID, &b_hcal_detID);
0319 fChain->SetBranchAddress("hcal_cellHot", &hcal_cellHot, &b_hcal_cellHot);
0320 fChain->SetBranchAddress("activeLength", &activeLength, &b_activeLength);
0321 fChain->SetBranchAddress("hcal_edepth1", &hcal_edepth1, &b_hcal_edepth1);
0322 fChain->SetBranchAddress("hcal_edepth2", &hcal_edepth2, &b_hcal_edepth2);
0323 fChain->SetBranchAddress("hcal_edepth3", &hcal_edepth3, &b_hcal_edepth3);
0324 fChain->SetBranchAddress("hcal_edepth4", &hcal_edepth4, &b_hcal_edepth4);
0325 fChain->SetBranchAddress("hcal_edepth5", &hcal_edepth5, &b_hcal_edepth5);
0326 fChain->SetBranchAddress("hcal_edepth6", &hcal_edepth6, &b_hcal_edepth6);
0327 fChain->SetBranchAddress("hcal_edepth7", &hcal_edepth7, &b_hcal_edepth7);
0328 fChain->SetBranchAddress("hcal_activeL1", &hcal_activeL1, &b_hcal_activeL1);
0329 fChain->SetBranchAddress("hcal_activeL2", &hcal_activeL2, &b_hcal_activeL2);
0330 fChain->SetBranchAddress("hcal_activeL3", &hcal_activeL3, &b_hcal_activeL3);
0331 fChain->SetBranchAddress("hcal_activeL4", &hcal_activeL4, &b_hcal_activeL4);
0332 fChain->SetBranchAddress("hcal_activeL5", &hcal_activeL5, &b_hcal_activeL5);
0333 fChain->SetBranchAddress("hcal_activeL6", &hcal_activeL6, &b_hcal_activeL6);
0334 fChain->SetBranchAddress("hcal_activeL7", &hcal_activeL7, &b_hcal_activeL7);
0335 fChain->SetBranchAddress("activeLengthHot", &activeLengthHot, &b_activeLengthHot);
0336 fChain->SetBranchAddress("hcal_edepthHot1", &hcal_edepthHot1, &b_hcal_edepthHot1);
0337 fChain->SetBranchAddress("hcal_edepthHot2", &hcal_edepthHot2, &b_hcal_edepthHot2);
0338 fChain->SetBranchAddress("hcal_edepthHot3", &hcal_edepthHot3, &b_hcal_edepthHot3);
0339 fChain->SetBranchAddress("hcal_edepthHot4", &hcal_edepthHot4, &b_hcal_edepthHot4);
0340 fChain->SetBranchAddress("hcal_edepthHot5", &hcal_edepthHot5, &b_hcal_edepthHot5);
0341 fChain->SetBranchAddress("hcal_edepthHot6", &hcal_edepthHot6, &b_hcal_edepthHot6);
0342 fChain->SetBranchAddress("hcal_edepthHot7", &hcal_edepthHot7, &b_hcal_edepthHot7);
0343 fChain->SetBranchAddress("hcal_activeHotL1", &hcal_activeHotL1, &b_hcal_activeHotL1);
0344 fChain->SetBranchAddress("hcal_activeHotL2", &hcal_activeHotL2, &b_hcal_activeHotL2);
0345 fChain->SetBranchAddress("hcal_activeHotL3", &hcal_activeHotL3, &b_hcal_activeHotL3);
0346 fChain->SetBranchAddress("hcal_activeHotL4", &hcal_activeHotL4, &b_hcal_activeHotL4);
0347 fChain->SetBranchAddress("hcal_activeHotL5", &hcal_activeHotL5, &b_hcal_activeHotL5);
0348 fChain->SetBranchAddress("hcal_activeHotL6", &hcal_activeHotL6, &b_hcal_activeHotL6);
0349 fChain->SetBranchAddress("hcal_activeHotL7", &hcal_activeHotL7, &b_hcal_activeHotL7);
0350 Notify();
0351 }
0352
0353 void HBHEMuonOfflineSimAnalyzer::Loop() {
0354
0355 if (fChain == 0)
0356 return;
0357
0358 Long64_t nentries = fChain->GetEntriesFast();
0359
0360 if (debug_)
0361 std::cout << "nevent = " << nentries << std::endl;
0362
0363 Long64_t nbytes = 0, nb = 0;
0364
0365 for (Long64_t jentry = 0; jentry < nentries; jentry++) {
0366 Long64_t ientry = LoadTree(jentry);
0367 if (ientry < 0)
0368 break;
0369 nb = fChain->GetEntry(jentry);
0370 nbytes += nb;
0371
0372 if (debug_) {
0373 std::cout << "ecal_det_id " << ecal_detID << std::endl;
0374 std::cout << "hcal_det_id " << std::hex << hcal_detID << std::dec;
0375 }
0376 int etaHcal, phiHcal, depthHcal;
0377 etaPhiHcal(hcal_detID, etaHcal, phiHcal, depthHcal);
0378
0379 int eta = (etaHcal > 0) ? (etaHcal - 1) : -(1 + etaHcal);
0380 int nDepth = NDepthBins(eta + 1, phiHcal);
0381 int nPhi = NPhiBins(eta + 1);
0382
0383 double phiYHcal = (phiHcal - 0.5);
0384 if (debug_)
0385 std::cout << "phiHcal" << phiHcal << " phiYHcal" << phiYHcal << std::endl;
0386
0387 for (int cut = 0; cut < 3; ++cut) {
0388 bool select(false);
0389 if (cut == 0)
0390 select = tightMuon();
0391 else if (cut == 1)
0392 select = SoftMuon();
0393 else
0394 select = LooseMuon();
0395
0396 if (select) {
0397
0398 h_P_Muon[cut]->Fill(p_of_muon);
0399 h_Pt_Muon[cut]->Fill(pt_of_muon);
0400 h_Eta_Muon[cut]->Fill(eta_of_muon);
0401
0402 double energyFill;
0403 for (int dep = 0; dep < nDepth; ++dep) {
0404 if (debug_) {
0405 std::cout << "why on 15/2 only" << std::endl;
0406 std::cout << "dep:" << dep << std::endl;
0407 }
0408 int PHI = (nPhi > 36) ? (phiHcal - 1) : (phiHcal - 1) / 2;
0409 double en1(-9999), en2(-9999);
0410 if (dep == 0) {
0411 en1 = hcal_edepth1;
0412 en2 = hcal_edepthHot1;
0413 energyFill = (hcal_activeHotL1 > 0) ? hcal_activeHotL1 : 999;
0414 } else if (dep == 1) {
0415 en1 = hcal_edepth2;
0416 en2 = hcal_edepthHot2;
0417 energyFill = (hcal_activeHotL2 > 0) ? hcal_activeHotL2 : 999;
0418 if (debug_)
0419 std::cout << "problem here.. lets see if it got printed\n";
0420 } else if (dep == 2) {
0421 en1 = hcal_edepth3;
0422 en2 = hcal_edepthHot3;
0423 energyFill = (hcal_activeHotL3 > 0) ? hcal_activeHotL3 : 999;
0424 } else if (dep == 3) {
0425 en1 = hcal_edepth4;
0426 en2 = hcal_edepthHot4;
0427 if (debug_)
0428 std::cout << "Hello in 4" << std::endl;
0429 energyFill = (hcal_activeHotL4 > 0) ? hcal_activeHotL4 : 999;
0430 } else if (dep == 4) {
0431 en1 = hcal_edepth5;
0432 en2 = hcal_edepthHot5;
0433 energyFill = (hcal_activeHotL5 > 0) ? hcal_activeHotL5 : 999;
0434 } else if (dep == 5) {
0435 if (debug_)
0436 std::cout << "Energy in depth 6 " << maxDepth_ << ":" << hcal_edepth6 << ":" << hcal_edepthHot6
0437 << std::endl;
0438 en1 = (maxDepth_ > 5) ? hcal_edepth6 : 0;
0439 en2 = (maxDepth_ > 5) ? hcal_edepthHot6 : 0;
0440 energyFill = (hcal_activeHotL6 > 0) ? hcal_activeHotL6 : 999;
0441 } else if (dep == 6) {
0442 if (debug_)
0443 std::cout << "Energy in depth 7 " << maxDepth_ << ":" << hcal_edepth7 << ":" << hcal_edepthHot7
0444 << std::endl;
0445 en1 = (maxDepth_ > 6) ? hcal_edepth7 : 0;
0446 en2 = (maxDepth_ > 6) ? hcal_edepthHot7 : 0;
0447 energyFill = (hcal_activeHotL7 > 0) ? hcal_activeHotL7 : 999;
0448 }
0449
0450 if (debug_) {
0451 std::cout << " Debug2" << std::endl;
0452 std::cout << "ok1" << en1 << std::endl;
0453 std::cout << "ok2" << en2 << std::endl;
0454 }
0455 bool ok1 = (en1 > -9999);
0456 bool ok2 = (en2 > -9999);
0457
0458 if (debug_)
0459 std::cout << "Before Index" << std::endl;
0460
0461 int ind = (etaHcal > 0) ? indxEta[eta][dep][PHI] : 1 + indxEta[eta][dep][PHI];
0462
0463 if (debug_) {
0464 std::cout << "ieta " << eta << "depth " << dep << "indxEta[eta][dep]:" << indxEta[eta][dep][PHI]
0465 << std::endl;
0466 std::cout << "index showing eta,depth:" << ind << std::endl;
0467 std::cout << "etaHcal: " << etaHcal << " eta " << eta << " dep " << dep << " indx " << ind << std::endl;
0468 }
0469 if (!(matchedId))
0470 continue;
0471 if (ok1) {
0472 if (debug_)
0473 std::cout << "enter ok1" << std::endl;
0474 if (hcal_cellHot == 1) {
0475 if (en2 > 0) {
0476 h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[cut][ind]->Fill(en2 / energyFill);
0477 }
0478 if (debug_)
0479 std::cout << "enter hot cell" << std::endl;
0480 }
0481 }
0482
0483 if (ok2) {
0484 if (debug_)
0485 std::cout << "enter ok2" << std::endl;
0486 if (hcal_cellHot != 1) {
0487 }
0488 }
0489
0490 if (debug_)
0491 std::cout << "ETA \t" << eta << "DEPTH \t" << dep << std::endl;
0492 }
0493 }
0494 }
0495 }
0496 close();
0497 }
0498
0499 Bool_t HBHEMuonOfflineSimAnalyzer::Notify() {
0500
0501
0502
0503
0504
0505
0506 return kTRUE;
0507 }
0508
0509 void HBHEMuonOfflineSimAnalyzer::Show(Long64_t entry) {
0510
0511
0512 if (!fChain)
0513 return;
0514 fChain->Show(entry);
0515 }
0516
0517 void HBHEMuonOfflineSimAnalyzer::BookHistograms(const char *fname) {
0518 output_file = TFile::Open(fname, "RECREATE");
0519
0520 std::string type[] = {"tight", "soft", "loose"};
0521 char name[128], title[500];
0522
0523 std::cout << "BookHistograms" << std::endl;
0524
0525 nHist = 0;
0526 for (int eta = 0; eta < 29; ++eta) {
0527 int nDepth = NDepthBins(eta + 1, -1);
0528 int nPhi = NPhiBins(eta + 1);
0529 for (int depth = 0; depth < nDepth; depth++) {
0530 for (int PHI = 0; PHI < nPhi; ++PHI) {
0531 indxEta[eta][depth][PHI] = nHist;
0532 nHist += 2;
0533 }
0534 }
0535 }
0536
0537 for (int i = 0; i < 3; ++i) {
0538 sprintf(name, "h_Pt_Muon_%s", type[i].c_str());
0539 sprintf(title, "p_{T} of %s muons (GeV)", type[i].c_str());
0540 h_Pt_Muon[i] = new TH1D(name, title, 100, 0, 200);
0541
0542 sprintf(name, "h_Eta_Muon_%s", type[i].c_str());
0543 sprintf(title, "#eta of %s muons", type[i].c_str());
0544 h_Eta_Muon[i] = new TH1D(name, title, 50, -2.5, 2.5);
0545
0546 sprintf(name, "h_Phi_Muon_%s", type[i].c_str());
0547 sprintf(title, "#phi of %s muons", type[i].c_str());
0548 h_Phi_Muon[i] = new TH1D(name, title, 100, -3.1415926, 3.1415926);
0549
0550 sprintf(name, "h_P_Muon_%s", type[i].c_str());
0551 sprintf(title, "p of %s muons (GeV)", type[i].c_str());
0552 h_P_Muon[i] = new TH1D(name, title, 100, 0, 200);
0553
0554 sprintf(name, "h_PF_Muon_%s", type[i].c_str());
0555 sprintf(title, "PF %s muons (GeV)", type[i].c_str());
0556 h_PF_Muon[i] = new TH1D(name, title, 2, 0, 2);
0557
0558 sprintf(name, "h_Global_Muon_Chi2_%s", type[i].c_str());
0559 sprintf(title, "Chi2 Global %s muons (GeV)", type[i].c_str());
0560 h_GlobTrack_Chi[i] = new TH1D(name, title, 15, 0, 15);
0561
0562 sprintf(name, "h_Global_Muon_Hits_%s", type[i].c_str());
0563 sprintf(title, "Global Hits %s muons (GeV)", type[i].c_str());
0564 h_Global_Muon_Hits[i] = new TH1D(name, title, 10, 0, 10);
0565
0566 sprintf(name, "h_Matched_Stations_%s", type[i].c_str());
0567 sprintf(title, "Matched Stations %s muons (GeV)", type[i].c_str());
0568 h_MatchedStations[i] = new TH1D(name, title, 10, 0, 10);
0569
0570 sprintf(name, "h_Transverse_ImpactParameter_%s", type[i].c_str());
0571 sprintf(title, "Transverse_ImpactParameter of %s muons (GeV)", type[i].c_str());
0572 h_Tight_TransImpactparameter[i] = new TH1D(name, title, 50, 0, 10);
0573
0574 sprintf(name, "h_Longitudinal_ImpactParameter_%s", type[i].c_str());
0575 sprintf(title, "Longitudinal_ImpactParameter of %s muons (GeV)", type[i].c_str());
0576 h_Tight_LongitudinalImpactparameter[i] = new TH1D(name, title, 20, 0, 10);
0577
0578 sprintf(name, "h_InnerTrack_PixelHits_%s", type[i].c_str());
0579 sprintf(title, "InnerTrack_PixelHits of %s muons (GeV)", type[i].c_str());
0580 h_InnerTrackPixelHits[i] = new TH1D(name, title, 20, 0, 20);
0581
0582 sprintf(name, "h_TrackLayers_%s", type[i].c_str());
0583 sprintf(title, "No. of Tracker Layers of %s muons (GeV)", type[i].c_str());
0584 h_TrackerLayer[i] = new TH1D(name, title, 20, 0, 20);
0585 ;
0586
0587 sprintf(name, "h_IsolationR04_%s", type[i].c_str());
0588 sprintf(title, "IsolationR04 %s muons (GeV)", type[i].c_str());
0589 h_IsolationR04[i] = new TH1D(name, title, 45, 0, 5);
0590 ;
0591
0592 sprintf(name, "h_Global_Muon_%s", type[i].c_str());
0593 sprintf(title, "Global %s muons (GeV)", type[i].c_str());
0594 h_Global_Muon[i] = new TH1D(name, title, 2, 0, 2);
0595
0596 sprintf(name, "h_TransImpactParameter_%s", type[i].c_str());
0597 sprintf(title, "TransImpactParameter of %s muons (GeV)", type[i].c_str());
0598 h_TransImpactParameter[i] = new TH1D(name, title, 100, 0, 0.5);
0599
0600 sprintf(name, "h_TransImpactParameterBin1_%s", type[i].c_str());
0601 sprintf(title, "TransImpactParameter of %s muons (GeV) in -1.5 <= #phi <= 0.5", type[i].c_str());
0602 h_TransImpactParameterBin1[i] = new TH1D(name, title, 100, 0, 0.5);
0603
0604 sprintf(name, "h_TransImpactParameterBin2_%s", type[i].c_str());
0605 sprintf(title, "TransImpactParameter of %s muons (GeV) in #phi> 0.5 and #phi< -1.5 ", type[i].c_str());
0606 h_TransImpactParameterBin2[i] = new TH1D(name, title, 100, 0, 0.5);
0607
0608 sprintf(name, "h_LongImpactParameter_%s", type[i].c_str());
0609 sprintf(title, "LongImpactParameter of %s muons (GeV)", type[i].c_str());
0610 h_LongImpactParameter[i] = new TH1D(name, title, 100, 0, 30);
0611
0612 sprintf(name, "h_LongImpactParameterBin1_%s", type[i].c_str());
0613 sprintf(title, "LongImpactParameter of %s muons (GeV) in -1.5 <= #phi <= 0.5", type[i].c_str());
0614 h_LongImpactParameterBin1[i] = new TH1D(name, title, 100, 0, 30);
0615
0616 sprintf(name, "h_LongImpactParameterBin2_%s", type[i].c_str());
0617 sprintf(title, "LongImpactParameter of %s muons (GeV) in #phi> 0.5 and #phi< -1.5 ", type[i].c_str());
0618 h_LongImpactParameterBin2[i] = new TH1D(name, title, 100, 0, 30);
0619
0620 sprintf(name, "h_2D_Bin1_%s", type[i].c_str());
0621 sprintf(title, "Trans/Long ImpactParameter of %s muons (GeV) in -1.5 <= #phi< 0.5 ", type[i].c_str());
0622 h_2D_Bin1[i] = new TH2D(name, title, 100, 0, 0.5, 100, 0, 30);
0623
0624 sprintf(name, "h_2D_Bin2_%s", type[i].c_str());
0625 sprintf(title, "Trans/Long ImpactParameter of %s muons (GeV) in #phi> 0.5 and #phi< -1.5 ", type[i].c_str());
0626 h_2D_Bin2[i] = new TH2D(name, title, 100, 0, 0.5, 100, 0, 30);
0627
0628 sprintf(name, "h_ecal_energy_%s", type[i].c_str());
0629 sprintf(title, "ECAL energy for %s muons", type[i].c_str());
0630 h_ecal_energy[i] = new TH1D(name, title, 1000, -10.0, 90.0);
0631
0632 sprintf(name, "h_hcal_energy_%s", type[i].c_str());
0633 sprintf(title, "HCAL energy for %s muons", type[i].c_str());
0634 h_hcal_energy[i] = new TH1D(name, title, 500, -10.0, 90.0);
0635
0636 sprintf(name, "h_3x3_ecal_%s", type[i].c_str());
0637 sprintf(title, "ECAL energy in 3x3 for %s muons", type[i].c_str());
0638 h_3x3_ecal[i] = new TH1D(name, title, 1000, -10.0, 90.0);
0639
0640 sprintf(name, "h_1x1_hcal_%s", type[i].c_str());
0641 sprintf(title, "HCAL energy in 1x1 for %s muons", type[i].c_str());
0642 h_1x1_hcal[i] = new TH1D(name, title, 500, -10.0, 90.0);
0643
0644 sprintf(name, "h_EtaX_hcal_%s", type[i].c_str());
0645 sprintf(title, "HCAL energy as a function of i#eta for %s muons", type[i].c_str());
0646 h_EtaX_hcal[i] = new TProfile(name, title, 60, -30.0, 30.0);
0647
0648 sprintf(name, "h_PhiY_hcal_%s", type[i].c_str());
0649 sprintf(title, "HCAL energy as a function of i#phi for %s muons", type[i].c_str());
0650 h_PhiY_hcal[i] = new TProfile(name, title, 72, 0, 72);
0651
0652 sprintf(name, "h_EtaX_ecal_%s", type[i].c_str());
0653 sprintf(title, "EB energy as a function of i#eta for %s muons", type[i].c_str());
0654 h_EtaX_ecal[i] = new TProfile(name, title, 170, -85.0, 85.0);
0655
0656 sprintf(name, "h_PhiY_ecal_%s", type[i].c_str());
0657 sprintf(title, "EB energy as a function of i#phi for %s muons", type[i].c_str());
0658 h_PhiY_ecal[i] = new TProfile(name, title, 360, 0, 360);
0659
0660 sprintf(name, "h_Eta_ecal_%s", type[i].c_str());
0661 sprintf(title, "ECAL energy as a function of #eta for %s muons", type[i].c_str());
0662 h_Eta_ecal[i] = new TProfile(name, title, 100, -2.5, 2.5);
0663
0664 sprintf(name, "h_Phi_ecal_%s", type[i].c_str());
0665 sprintf(title, "ECAL energy as a function of #phi for %s muons", type[i].c_str());
0666 h_Phi_ecal[i] = new TProfile(name, title, 100, -3.1415926, 3.1415926);
0667
0668 sprintf(name, "h_MuonHittingEcal_%s", type[i].c_str());
0669 sprintf(title, "%s muons hitting ECAL", type[i].c_str());
0670 h_MuonHittingEcal[i] = new TH1D(name, title, 100, 0, 5.0);
0671
0672 sprintf(name, "h_HotCell_%s", type[i].c_str());
0673 sprintf(title, "Hot cell for %s muons", type[i].c_str());
0674 h_HotCell[i] = new TH1D(name, title, 100, 0, 2);
0675
0676 std::cout << "problem here" << std::endl;
0677 for (int eta = 0; eta < 29; ++eta) {
0678 int nDepth = NDepthBins(eta + 1, -1);
0679 int nPhi = NPhiBins(eta + 1);
0680 for (int depth = 0; depth < nDepth; ++depth) {
0681 for (int PHI = 0; PHI < nPhi; ++PHI) {
0682 int PHI0 = (nPhi == 72) ? PHI + 1 : 2 * PHI + 1;
0683 int ih = indxEta[eta][depth][PHI];
0684 std::cout << "eta:" << eta << " depth:" << depth << " PHI:" << PHI << ":" << PHI0 << " ih:" << ih
0685 << std::endl;
0686
0687 sprintf(name,
0688 "h_Hot_MuonEnergy_hc_%d_%d_%d_%s_HotCell_ByActiveLength",
0689 (eta + 1),
0690 (depth + 1),
0691 PHI0,
0692 type[i].c_str());
0693 sprintf(title,
0694 "HCAL energy in hot tower (i#eta=%d, depth=%d, i#phi = %d) for extrapolated %s muons (Hot Cell) "
0695 "divided by Active Length",
0696 (eta + 1),
0697 (depth + 1),
0698 PHI0,
0699 type[i].c_str());
0700 h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[i][ih] = new TH1D(name, title, 4000, 0.0, 1.0);
0701 h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[i][ih]->Sumw2();
0702
0703 ih++;
0704 sprintf(name,
0705 "h_Hot_MuonEnergy_hc_%d_%d_%d_%s_HotCell_ByActiveLength",
0706 -(eta + 1),
0707 (depth + 1),
0708 PHI0,
0709 type[i].c_str());
0710 sprintf(title,
0711 "HCAL energy in hot tower (i#eta=%d, depth=%d, i#phi=%d) for extrapolated %s muons (Hot Cell) "
0712 "divided by Active Length",
0713 -(eta + 1),
0714 (depth + 1),
0715 PHI0,
0716 type[i].c_str());
0717 h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[i][ih] = new TH1D(name, title, 4000, 0.0, 1.0);
0718 h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[i][ih]->Sumw2();
0719 }
0720 }
0721
0722 }
0723 }
0724
0725 }
0726
0727 bool HBHEMuonOfflineSimAnalyzer::LooseMuon() {
0728 if (pt_of_muon > 20.) {
0729 if (fabs(eta_of_muon) <= 5.0) {
0730 return true;
0731 }
0732 }
0733
0734 return false;
0735 }
0736
0737 bool HBHEMuonOfflineSimAnalyzer::SoftMuon() {
0738 if (pt_of_muon > 20.) {
0739 if (fabs(eta_of_muon) <= 5.0) {
0740 return true;
0741 }
0742 }
0743
0744 return false;
0745 }
0746
0747 bool HBHEMuonOfflineSimAnalyzer::tightMuon() {
0748 if (pt_of_muon > 20.) {
0749 if (fabs(eta_of_muon) <= 5.0) {
0750 return true;
0751 }
0752 }
0753
0754 return false;
0755 }
0756
0757 void HBHEMuonOfflineSimAnalyzer::etaPhiHcal(unsigned int detId, int &eta, int &phi, int &depth) {
0758 int zside, etaAbs;
0759 if ((detId & 0x1000000) == 0) {
0760 zside = (detId & 0x2000) ? (1) : (-1);
0761 etaAbs = (detId >> 7) & 0x3F;
0762 phi = detId & 0x7F;
0763 depth = (detId >> 14) & 0x1F;
0764 } else {
0765 zside = (detId & 0x80000) ? (1) : (-1);
0766 etaAbs = (detId >> 10) & 0x1FF;
0767 phi = detId & 0x3FF;
0768 depth = (detId >> 20) & 0xF;
0769 }
0770 eta = etaAbs * zside;
0771 }
0772
0773 void HBHEMuonOfflineSimAnalyzer::etaPhiEcal(
0774 unsigned int detId, int &type, int &zside, int &etaX, int &phiY, int &plane, int &strip) {
0775 type = ((detId >> 25) & 0x7);
0776 plane = strip = 0;
0777 if (type == 1) {
0778
0779 zside = (detId & 0x10000) ? (1) : (-1);
0780 etaX = (detId >> 9) & 0x7F;
0781 phiY = detId & 0x1FF;
0782 } else if (type == 2) {
0783 zside = (detId & 0x4000) ? (1) : (-1);
0784 etaX = (detId >> 7) & 0x7F;
0785 phiY = (detId & 0x7F);
0786 } else if (type == 3) {
0787 zside = (detId & 0x80000) ? (1) : (-1);
0788 etaX = (detId >> 6) & 0x3F;
0789
0790 phiY = (detId >> 12) & 0x3F;
0791
0792 plane = ((detId >> 18) & 0x1) + 1;
0793 strip = detId & 0x3F;
0794 } else {
0795 zside = etaX = phiY = 0;
0796 }
0797 }
0798
0799 void HBHEMuonOfflineSimAnalyzer::calculateP(double pt, double eta, double &pM) {
0800 pM = (pt * cos(2 * (1 / atan(exp(eta)))));
0801 }
0802
0803 void HBHEMuonOfflineSimAnalyzer::close() {
0804 output_file->cd();
0805 std::cout << "file yet to be Written" << std::endl;
0806 WriteHistograms();
0807
0808 std::cout << "file Written" << std::endl;
0809 output_file->Close();
0810 std::cout << "now doing return" << std::endl;
0811 }
0812
0813 void HBHEMuonOfflineSimAnalyzer::WriteHistograms() {
0814 std::string type[] = {"tight", "soft", "loose"};
0815 char name[128];
0816
0817 std::cout << "WriteHistograms" << std::endl;
0818 nHist = 0;
0819 for (int eta = 0; eta < 29; ++eta) {
0820 int nDepth = NDepthBins(eta + 1, -1);
0821 int nPhi = NPhiBins(eta + 1);
0822 if (debug_)
0823 std::cout << "Eta:" << eta << " nDepths " << nDepth << " nPhis " << nPhi << std::endl;
0824 for (int depth = 0; depth < nDepth; ++depth) {
0825 if (debug_)
0826 std::cout << "Eta:" << eta << "Depth:" << depth << std::endl;
0827 for (int PHI = 0; PHI < nPhi; ++PHI) {
0828 indxEta[eta][depth][PHI] = nHist;
0829 nHist += 2;
0830 }
0831 }
0832 }
0833
0834 TDirectory *d_output_file[3][29];
0835 for (int i = 0; i < 3; ++i) {
0836 h_Pt_Muon[i]->Write();
0837 h_Eta_Muon[i]->Write();
0838 h_Phi_Muon[i]->Write();
0839 h_P_Muon[i]->Write();
0840 h_PF_Muon[i]->Write();
0841
0842 h_GlobTrack_Chi[i]->Write();
0843 h_Global_Muon_Hits[i]->Write();
0844 h_MatchedStations[i]->Write();
0845
0846 h_Tight_TransImpactparameter[i]->Write();
0847 h_Tight_LongitudinalImpactparameter[i]->Write();
0848
0849 h_InnerTrackPixelHits[i]->Write();
0850 h_TrackerLayer[i]->Write();
0851 h_IsolationR04[i]->Write();
0852
0853 h_Global_Muon[i]->Write();
0854 h_TransImpactParameter[i]->Write();
0855 ;
0856 h_TransImpactParameterBin1[i]->Write();
0857 h_TransImpactParameterBin2[i]->Write();
0858
0859 h_LongImpactParameter[i]->Write();
0860 h_LongImpactParameterBin1[i]->Write();
0861 h_LongImpactParameterBin2[i]->Write();
0862
0863 h_ecal_energy[i]->Write();
0864 h_hcal_energy[i]->Write();
0865 ;
0866 h_3x3_ecal[i]->Write();
0867 h_1x1_hcal[i]->Write();
0868 ;
0869
0870 h_EtaX_hcal[i]->Write();
0871 h_PhiY_hcal[i]->Write();
0872 ;
0873
0874 h_EtaX_ecal[i]->Write();
0875 ;
0876 h_PhiY_ecal[i]->Write();
0877 ;
0878 h_Eta_ecal[i]->Write();
0879 ;
0880 h_Phi_ecal[i]->Write();
0881 ;
0882 h_MuonHittingEcal[i]->Write();
0883 ;
0884 h_HotCell[i]->Write();
0885 ;
0886
0887 output_file->cd();
0888 for (int eta = 0; eta < 29; ++eta) {
0889 int nDepth = NDepthBins(eta + 1, -1);
0890 int nPhi = NPhiBins(eta + 1);
0891 sprintf(name, "Dir_muon_type_%s_ieta%d", type[i].c_str(), eta + 1);
0892 d_output_file[i][eta] = output_file->mkdir(name);
0893 d_output_file[i][eta]->cd();
0894 for (int depth = 0; depth < nDepth; ++depth) {
0895 for (int PHI = 0; PHI < nPhi; ++PHI) {
0896 int ih = indxEta[eta][depth][PHI];
0897 h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[i][ih]->Write();
0898 ih++;
0899 h_Hot_MuonEnergy_hcal_HotCell_VsActiveLength[i][ih]->Write();
0900 }
0901 }
0902 output_file->cd();
0903 }
0904 }
0905 output_file->cd();
0906 }
0907
0908 int HBHEMuonOfflineSimAnalyzer::NDepthBins(int eta, int phi) {
0909
0910 int nDepthR1[29] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 3, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 2};
0911
0912 int nDepthR2[29] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 4, 3, 5, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 3};
0913
0914 int nDepthR3[29] = {4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 5, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 3};
0915
0916 int nDepthR4[29] = {4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7};
0917
0918
0919
0920
0921
0922
0923
0924
0925 int nbin(0);
0926 if (modeLHC_ == 0) {
0927 if (eta <= 15) {
0928 nbin = maxDepthHB_;
0929 } else if (eta == 16) {
0930 nbin = 4;
0931 } else {
0932 nbin = maxDepthHE_;
0933 }
0934 } else if (modeLHC_ == 1) {
0935 nbin = nDepthR1[eta - 1];
0936 } else if (modeLHC_ == 2) {
0937 nbin = nDepthR2[eta - 1];
0938 } else if (modeLHC_ == 3) {
0939 nbin = nDepthR3[eta - 1];
0940 } else if (modeLHC_ == 4) {
0941 if (phi > 0) {
0942 if (eta >= 16 && phi >= 63 && phi <= 66) {
0943 nbin = nDepthR2[eta - 1];
0944 } else {
0945 nbin = nDepthR1[eta - 1];
0946 }
0947 } else {
0948 if (eta >= 16) {
0949 nbin = (nDepthR2[eta - 1] > nDepthR1[eta - 1]) ? nDepthR2[eta - 1] : nDepthR1[eta - 1];
0950 } else {
0951 nbin = nDepthR1[eta - 1];
0952 }
0953 }
0954 } else {
0955 if (eta > 0 && eta < 30) {
0956 nbin = nDepthR4[eta - 1];
0957 } else {
0958 nbin = nDepthR4[28];
0959 }
0960 }
0961 return nbin;
0962 }
0963
0964 int HBHEMuonOfflineSimAnalyzer::NPhiBins(int eta) {
0965 int nphi = (eta <= 20) ? 72 : 36;
0966 return nphi;
0967 }