GlobalCaloTrigger/plugins/L1GctValidation.cc

0001 #include "L1Trigger/GlobalCaloTrigger/plugins/L1GctValidation.h"
0002
0003 #include "FWCore/PluginManager/interface/ModuleDef.h"
0004 #include "FWCore/Framework/interface/MakerMacros.h"
0005 #include "FWCore/Framework/interface/ESHandle.h"
0006
0007 #include "FWCore/ServiceRegistry/interface/Service.h"
0008 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0009
0010 #include "DataFormats/L1CaloTrigger/interface/L1CaloCollections.h"
0011
0012 #include "DataFormats/Math/interface/deltaPhi.h"
0013
0014 #include <cmath>
0015
0016 L1GctValidation::L1GctValidation(const edm::ParameterSet& iConfig)
0017     : m_gctinp_tag(iConfig.getUntrackedParameter<edm::InputTag>("rctInputTag", edm::InputTag("rctDigis"))),
0018       m_energy_tag(iConfig.getUntrackedParameter<edm::InputTag>("gctInputTag", edm::InputTag("gctDigis"))),
0019       m_jfParsToken(esConsumes<L1GctJetFinderParams, L1GctJetFinderParamsRcd>()),
0020       m_htMissScaleToken(esConsumes<L1CaloEtScale, L1HtMissScaleRcd>()),
0021       m_hfRingEtScaleToken(esConsumes<L1CaloEtScale, L1HfRingEtScaleRcd>()) {
0022   usesResource(TFileService::kSharedResource);
0023 }
0024
0025 L1GctValidation::~L1GctValidation() {
0026   // do anything here that needs to be done at desctruction time
0027   // (e.g. close files, deallocate resources etc.)
0028 }
0029
0030 //
0031 // member functions
0032 //
0033
0034 // ------------ method called to for each event  ------------
0035 void L1GctValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0036   using namespace edm;
0037
0038   // Get the scales from the event setup
0039   ESHandle<L1GctJetFinderParams> jfPars = iSetup.getHandle(m_jfParsToken);
0040
0041   double lsbForEt = jfPars.product()->getRgnEtLsbGeV();
0042   double lsbForHt = jfPars.product()->getHtLsbGeV();
0043
0044   unsigned httJetThreshold = static_cast<int>(jfPars.product()->getHtJetEtThresholdGeV() / lsbForHt);
0045   unsigned htmJetThreshold = static_cast<int>(jfPars.product()->getMHtJetEtThresholdGeV() / lsbForHt);
0046
0047   ESHandle<L1CaloEtScale> htMissScale = iSetup.getHandle(m_htMissScaleToken);
0048   ESHandle<L1CaloEtScale> hfRingEtScale = iSetup.getHandle(m_hfRingEtScaleToken);
0049
0050   // Get the Gct energy sums from the event
0051   Handle<L1GctEtTotalCollection> sumEtColl;
0052   iEvent.getByLabel(m_energy_tag, sumEtColl);
0053   Handle<L1GctEtHadCollection> sumHtColl;
0054   iEvent.getByLabel(m_energy_tag, sumHtColl);
0055   Handle<L1GctEtMissCollection> missEtColl;
0056   iEvent.getByLabel(m_energy_tag, missEtColl);
0057   Handle<L1GctHtMissCollection> missHtColl;
0058   iEvent.getByLabel(m_energy_tag, missHtColl);
0059
0060   // Get the input calo regions from the event (for checking MEt)
0061   Handle<L1CaloRegionCollection> inputColl;
0062   iEvent.getByLabel(m_gctinp_tag, inputColl);
0063
0064   // Get the internal jet data from the event (for checking Ht)
0065   Handle<L1GctInternJetDataCollection> internalJetsColl;
0066   iEvent.getByLabel(m_energy_tag, internalJetsColl);
0067
0068   double etTot = 0.0;
0069   for (L1GctEtTotalCollection::const_iterator jbx = sumEtColl->begin(); jbx != sumEtColl->end(); jbx++) {
0070     if (jbx->bx() == 0) {
0071       etTot = static_cast<double>(jbx->et());
0072     }
0073   }
0074
0075   double etHad = 0.0;
0076   for (L1GctEtHadCollection::const_iterator jbx = sumHtColl->begin(); jbx != sumHtColl->end(); jbx++) {
0077     if (jbx->bx() == 0) {
0078       etHad = static_cast<double>(jbx->et());
0079     }
0080   }
0081
0082   double etMiss = 0.0;
0083   double etMAng = 0.0;
0084   for (L1GctEtMissCollection::const_iterator jbx = missEtColl->begin(); jbx != missEtColl->end(); jbx++) {
0085     if (jbx->bx() == 0) {
0086       etMiss = static_cast<double>(jbx->et());
0087       int phibin = jbx->phi();
0088       if (phibin >= 36)
0089         phibin -= 72;
0090       double etMPhi = static_cast<double>(phibin);
0091
0092       etMAng = (etMPhi + 0.5) * M_PI / 36.;
0093     }
0094   }
0095
0096   double etTotFromRegions = 0.0;
0097   double exTotFromRegions = 0.0;
0098   double eyTotFromRegions = 0.0;
0099   for (L1CaloRegionCollection::const_iterator jrg = inputColl->begin(); jrg != inputColl->end(); jrg++) {
0100     if (jrg->bx() == 0) {
0101       double rgEt = static_cast<double>(jrg->et()) * lsbForEt;
0102       double rgPhibin = static_cast<double>(jrg->id().iphi());
0103       double rgPh = (rgPhibin + 0.5) * M_PI / 9.;
0104
0105       etTotFromRegions += rgEt;
0106       exTotFromRegions += rgEt * cos(rgPh);
0107       eyTotFromRegions += rgEt * sin(rgPh);
0108     }
0109   }
0110
0111   double htMissGct = 0.0;
0112   double htMissAng = 0.0;
0113   double htMissGeV = 0.0;
0114   for (L1GctHtMissCollection::const_iterator jbx = missHtColl->begin(); jbx != missHtColl->end(); jbx++) {
0115     if (jbx->bx() == 0) {
0116       htMissGct = static_cast<double>(jbx->et());
0117       htMissGeV = htMissScale->et(jbx->et());
0118       int phibin = jbx->phi();
0119       if (phibin >= 9)
0120         phibin -= 18;
0121       double htMPhi = static_cast<double>(phibin);
0122       htMissAng = (htMPhi + 0.5) * M_PI / 9.;
0123     }
0124   }
0125
0126   double htFromJets = 0.0;
0127   double hxFromJets = 0.0;
0128   double hyFromJets = 0.0;
0129   for (L1GctInternJetDataCollection::const_iterator jet = internalJetsColl->begin(); jet != internalJetsColl->end();
0130        jet++) {
0131     if (jet->bx() == 0 && !jet->empty()) {
0132       unsigned jetEtGct = jet->et();
0133       double jetEt = static_cast<double>(jetEtGct);
0134       int phibin = jet->regionId().iphi();
0135       if (phibin >= 9)
0136         phibin -= 18;
0137       // The phi bin centres are at 0, 20, 40, ... degrees
0138       double jetAng = (static_cast<double>(phibin)) * M_PI / 9.;
0139       if (jetEtGct > httJetThreshold) {
0140         htFromJets += jetEt;
0141       }
0142       if (jetEtGct > htmJetThreshold) {
0143         hxFromJets += jetEt * cos(jetAng);
0144         hyFromJets += jetEt * sin(jetAng);
0145       }
0146     }
0147   }
0148
0149   double dPhiMetMht = deltaPhi(etMAng, htMissAng);
0150
0151   theSumEtInLsb->Fill(etTot);
0152   theSumHtInLsb->Fill(etHad);
0153   theMissEtInLsb->Fill(etMiss);
0154   theMissHtInLsb->Fill(htMissGct);
0155   theSumEtInGeV->Fill(etTot * lsbForEt);
0156   theSumHtInGeV->Fill(etHad * lsbForHt);
0157   theMissEtInGeV->Fill(etMiss * lsbForEt);
0158   theMissEtAngle->Fill(etMAng);
0159   theMissEtVector->Fill(etMiss * lsbForEt * cos(etMAng), etMiss * lsbForEt * sin(etMAng));
0160   if (htMissGct < 126.5) {
0161     theMissHtInGeV->Fill(htMissGeV);
0162     theMissHtAngle->Fill(htMissAng);
0163     theMissHtVector->Fill(htMissGeV * cos(htMissAng), htMissGeV * sin(htMissAng));
0164   }
0165
0166   theSumEtVsInputRegions->Fill(etTot * lsbForEt, etTotFromRegions);
0167   theMissEtMagVsInputRegions->Fill(etMiss * lsbForEt,
0168                                    sqrt(exTotFromRegions * exTotFromRegions + eyTotFromRegions * eyTotFromRegions));
0169   theMissEtAngleVsInputRegions->Fill(etMAng, atan2(-eyTotFromRegions, -exTotFromRegions));
0170   theMissHtMagVsInputRegions->Fill(htMissGeV,
0171                                    sqrt(exTotFromRegions * exTotFromRegions + eyTotFromRegions * eyTotFromRegions));
0172
0173   theMissEtVsMissHt->Fill(etMiss * lsbForEt, htMissGeV);
0174   theMissEtVsMissHtAngle->Fill(etMAng, htMissAng);
0175   theDPhiVsMissEt->Fill(dPhiMetMht, etMiss * lsbForEt);
0176   theDPhiVsMissHt->Fill(dPhiMetMht, htMissGeV);
0177
0178   theHtVsInternalJetsSum->Fill(etHad * lsbForHt, htFromJets * lsbForHt);
0179   if (htMissGct < 126.5) {
0180     theMissHtVsInternalJetsSum->Fill(htMissGeV, sqrt(hxFromJets * hxFromJets + hyFromJets * hyFromJets) * lsbForHt);
0181     theMissHtPhiVsInternalJetsSum->Fill(htMissAng, atan2(-hyFromJets, -hxFromJets));
0182     theMissHxVsInternalJetsSum->Fill(htMissGeV * cos(htMissAng), hxFromJets * lsbForHt);
0183     theMissHyVsInternalJetsSum->Fill(htMissGeV * sin(htMissAng), hyFromJets * lsbForHt);
0184   }
0185
0186   // Get minbias trigger quantities from HF
0187   Handle<L1GctHFRingEtSumsCollection> HFEtSumsColl;
0188   Handle<L1GctHFBitCountsCollection> HFCountsColl;
0189   iEvent.getByLabel(m_energy_tag, HFEtSumsColl);
0190   iEvent.getByLabel(m_energy_tag, HFCountsColl);
0191
0192   for (L1GctHFRingEtSumsCollection::const_iterator es = HFEtSumsColl->begin(); es != HFEtSumsColl->end(); es++) {
0193     if (es->bx() == 0) {
0194       theHfRing0EtSumPositiveEta->Fill(hfRingEtScale->et(es->etSum(0)));
0195       theHfRing0EtSumNegativeEta->Fill(hfRingEtScale->et(es->etSum(1)));
0196       theHfRing1EtSumPositiveEta->Fill(hfRingEtScale->et(es->etSum(2)));
0197       theHfRing1EtSumNegativeEta->Fill(hfRingEtScale->et(es->etSum(3)));
0198     }
0199   }
0200
0201   for (L1GctHFBitCountsCollection::const_iterator bc = HFCountsColl->begin(); bc != HFCountsColl->end(); bc++) {
0202     if (bc->bx() == 0) {
0203       theHfRing0CountPositiveEta->Fill(bc->bitCount(0));
0204       theHfRing0CountNegativeEta->Fill(bc->bitCount(1));
0205       theHfRing1CountPositiveEta->Fill(bc->bitCount(2));
0206       theHfRing1CountNegativeEta->Fill(bc->bitCount(3));
0207     }
0208   }
0209 }
0210
0211 // ------------ method called once each job just before starting event loop  ------------
0212 void L1GctValidation::beginJob() {
0213   edm::Service<TFileService> fs;
0214
0215   TFileDirectory dir0 = fs->mkdir("L1GctEtSums");
0216
0217   theSumEtInLsb = dir0.make<TH1F>("SumEtInLsb", "Total Et (GCT units)", 128, 0., 2048.);
0218   theSumHtInLsb = dir0.make<TH1F>("SumHtInLsb", "Total Ht (GCT units)", 128, 0., 2048.);
0219   theMissEtInLsb = dir0.make<TH1F>("MissEtInLsb", "Missing Et magnitude (GCT units)", 128, 0., 1024.);
0220   theMissHtInLsb = dir0.make<TH1F>("MissHtInLsb", "Missing Ht magnitude (GCT units)", 128, 0., 127.);
0221   theSumEtInGeV = dir0.make<TH1F>("SumEtInGeV", "Total Et (in GeV)", 100, 0., 1000.);
0222   theSumHtInGeV = dir0.make<TH1F>("SumHtInGeV", "Total Ht (in GeV)", 100, 0., 1000.);
0223   theMissEtInGeV = dir0.make<TH1F>("MissEtInGeV", "Missing Et magnitude (in GeV)", 100, 0., 500.);
0224   theMissEtAngle = dir0.make<TH1F>("MissEtAngle", "Missing Et angle", 72, -M_PI, M_PI);
0225   theMissEtVector = dir0.make<TH2F>("MissEtVector", "Missing Ex vs Missing Ey", 100, -100., 100., 100, -100., 100.);
0226   theMissHtInGeV = dir0.make<TH1F>("MissHtInGeV", "Missing Ht magnitude (in GeV)", 100, 0., 500.);
0227   theMissHtAngle = dir0.make<TH1F>("MissHtAngle", "Missing Ht angle", 72, -M_PI, M_PI);
0228   theMissHtVector = dir0.make<TH2F>("MissHtVector", "Missing Hx vs Missing Hy", 100, -100., 100., 100, -100., 100.);
0229   theSumEtVsInputRegions =
0230       dir0.make<TH2F>("SumEtVsInputRegions", "Total Et vs sum of input regions", 100, 0., 1000., 100, 0., 1000.);
0231   theMissEtMagVsInputRegions = dir0.make<TH2F>(
0232       "MissEtMagVsInputRegions", "Missing Et magnitude vs sum of input regions", 100, 0., 500., 100, 0., 500.);
0233   theMissEtAngleVsInputRegions = dir0.make<TH2F>(
0234       "MissEtAngleVsInputRegions", "Missing Et angle vs sum of input regions", 72, -M_PI, M_PI, 72, -M_PI, M_PI);
0235   theMissHtMagVsInputRegions = dir0.make<TH2F>(
0236       "MissHtMagVsInputRegions", "Missing Ht magnitude vs sum of input regions", 100, 0., 500., 100, 0., 500.);
0237   theMissEtVsMissHt = dir0.make<TH2F>("MissEtVsMissHt", "Missing Et vs Missing Ht", 100, 0., 500., 100, 0., 500.);
0238   theMissEtVsMissHtAngle = dir0.make<TH2F>(
0239       "MissEtVsMissHtAngle", "Angle correlation Missing Et vs Missing Ht", 72, -M_PI, M_PI, 72, -M_PI, M_PI);
0240   theDPhiVsMissEt =
0241       dir0.make<TH2F>("theDPhiVsMissEt", "Angle difference MET-MHT vs MET magnitude", 72, -M_PI, M_PI, 100, 0., 500.);
0242   theDPhiVsMissHt =
0243       dir0.make<TH2F>("theDPhiVsMissHt", "Angle difference MET-MHT vs MHT magnitude", 72, -M_PI, M_PI, 100, 0., 500.);
0244
0245   theHtVsInternalJetsSum = dir0.make<TH2F>(
0246       "HtVsInternalJetsSum", "Ht vs scalar sum of jet Et values (in GeV)", 128, 0., 2048., 128, 0., 2048.);
0247   theMissHtVsInternalJetsSum = dir0.make<TH2F>(
0248       "MissHtVsInternalJetsSum", "Missing Ht vs vector sum of jet Et values (in GeV)", 128, 0., 512., 128, 0., 512.);
0249   theMissHtPhiVsInternalJetsSum = dir0.make<TH2F>("MissHtPhiVsInternalJetsSum",
0250                                                   "Angle correlation Missing Ht vs vector sum of jet Et values",
0251                                                   72,
0252                                                   -M_PI,
0253                                                   M_PI,
0254                                                   72,
0255                                                   -M_PI,
0256                                                   M_PI);
0257   theMissHxVsInternalJetsSum = dir0.make<TH2F>("MissHxVsInternalJetsSum",
0258                                                "Missing Ht x component vs sum of jet Et values (in GeV)",
0259                                                128,
0260                                                -256.,
0261                                                256.,
0262                                                128,
0263                                                -256.,
0264                                                256.);
0265   theMissHyVsInternalJetsSum = dir0.make<TH2F>("MissHyVsInternalJetsSum",
0266                                                "Missing Ht y component vs sum of jet Et values (in GeV)",
0267                                                128,
0268                                                -256.,
0269                                                256.,
0270                                                128,
0271                                                -256.,
0272                                                256.);
0273
0274   TFileDirectory dir1 = fs->mkdir("L1GctHfSumsAndJetCounts");
0275
0276   // Minimum bias triggers from Hf inner rings
0277   theHfRing0EtSumPositiveEta = dir1.make<TH1F>("HfRing0EtSumPositiveEta", "Hf Inner Ring0 Et eta+", 60, 0., 30.);
0278   theHfRing0EtSumNegativeEta = dir1.make<TH1F>("HfRing0EtSumNegativeEta", "Hf Inner Ring0 Et eta-", 60, 0., 30.);
0279   theHfRing1EtSumPositiveEta = dir1.make<TH1F>("HfRing1EtSumPositiveEta", "Hf Inner Ring1 Et eta+", 60, 0., 30.);
0280   theHfRing1EtSumNegativeEta = dir1.make<TH1F>("HfRing1EtSumNegativeEta", "Hf Inner Ring1 Et eta-", 60, 0., 30.);
0281   theHfRing0CountPositiveEta = dir1.make<TH1F>("HfRing0CountPositiveEta", "Hf Threshold bits Ring0 eta+", 20, 0., 20.);
0282   theHfRing0CountNegativeEta = dir1.make<TH1F>("HfRing0CountNegativeEta", "Hf Threshold bits Ring0 eta-", 20, 0., 20.);
0283   theHfRing1CountPositiveEta = dir1.make<TH1F>("HfRing1CountPositiveEta", "Hf Threshold bits Ring1 eta+", 20, 0., 20.);
0284   theHfRing1CountNegativeEta = dir1.make<TH1F>("HfRing1CountNegativeEta", "Hf Threshold bits Ring1 eta-", 20, 0., 20.);
0285 }
0286
0287 // ------------ method called once each job just after ending the event loop  ------------
0288 void L1GctValidation::endJob() {}
0289
0290 DEFINE_FWK_MODULE(L1GctValidation);