Reference Manual

Go to the documentation of this file.
00001 #include "L1Trigger/GlobalCaloTrigger/plugins/L1GctValidation.h"
00002 
00003 #include "FWCore/PluginManager/interface/ModuleDef.h"
00004 #include "FWCore/Framework/interface/MakerMacros.h"
00005 #include "FWCore/Framework/interface/ESHandle.h"
00006 
00007 #include "FWCore/ServiceRegistry/interface/Service.h"
00008 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00009 
00010 #include "CondFormats/L1TObjects/interface/L1GctJetFinderParams.h"
00011 #include "CondFormats/DataRecord/interface/L1GctJetFinderParamsRcd.h"
00012 
00013 #include "CondFormats/L1TObjects/interface/L1CaloEtScale.h"
00014 #include "CondFormats/DataRecord/interface/L1HtMissScaleRcd.h"
00015 #include "CondFormats/DataRecord/interface/L1HfRingEtScaleRcd.h"
00016 
00017 #include "DataFormats/L1CaloTrigger/interface/L1CaloCollections.h"
00018 
00019 #include "DataFormats/Math/interface/deltaPhi.h"
00020 
00021 #include <math.h>
00022 
00023 L1GctValidation::L1GctValidation(const edm::ParameterSet& iConfig) :
00024   m_gctinp_tag(iConfig.getUntrackedParameter<edm::InputTag>("rctInputTag", edm::InputTag("rctDigis"))),
00025   m_energy_tag(iConfig.getUntrackedParameter<edm::InputTag>("gctInputTag", edm::InputTag("gctDigis")))
00026 {
00027 }
00028 
00029 
00030 L1GctValidation::~L1GctValidation()
00031 {
00032  
00033    // do anything here that needs to be done at desctruction time
00034    // (e.g. close files, deallocate resources etc.)
00035 
00036 }
00037 
00038 //
00039 // member functions
00040 //
00041 
00042 // ------------ method called to for each event  ------------
00043 void
00044 L1GctValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
00045 {
00046   using namespace edm;
00047 
00048   // Get the scales from the event setup
00049   ESHandle< L1GctJetFinderParams > jfPars ;
00050   iSetup.get< L1GctJetFinderParamsRcd >().get( jfPars ) ; // which record?
00051 
00052   double lsbForEt = jfPars.product()->getRgnEtLsbGeV();
00053   double lsbForHt = jfPars.product()->getHtLsbGeV();
00054 
00055   unsigned httJetThreshold = static_cast<int>(jfPars.product()->getHtJetEtThresholdGeV()/lsbForHt);
00056   unsigned htmJetThreshold = static_cast<int>(jfPars.product()->getMHtJetEtThresholdGeV()/lsbForHt);
00057 
00058   ESHandle< L1CaloEtScale > htMissScale ;
00059   iSetup.get< L1HtMissScaleRcd >().get( htMissScale ) ; // which record?
00060   ESHandle< L1CaloEtScale > hfRingEtScale ;
00061   iSetup.get< L1HfRingEtScaleRcd >().get( hfRingEtScale ) ; // which record?
00062 
00063   // Get the Gct energy sums from the event
00064   Handle< L1GctEtTotalCollection > sumEtColl ;
00065   iEvent.getByLabel( m_energy_tag, sumEtColl ) ;
00066   Handle< L1GctEtHadCollection >   sumHtColl ;
00067   iEvent.getByLabel( m_energy_tag, sumHtColl ) ;
00068   Handle< L1GctEtMissCollection >  missEtColl ;
00069   iEvent.getByLabel( m_energy_tag, missEtColl ) ;
00070   Handle< L1GctHtMissCollection >  missHtColl ;
00071   iEvent.getByLabel( m_energy_tag, missHtColl ) ;
00072 
00073   // Get the input calo regions from the event (for checking MEt)
00074   Handle < L1CaloRegionCollection > inputColl ;
00075   iEvent.getByLabel( m_gctinp_tag, inputColl ) ;
00076 
00077   // Get the internal jet data from the event (for checking Ht)
00078   Handle < L1GctInternJetDataCollection > internalJetsColl;
00079   iEvent.getByLabel( m_energy_tag, internalJetsColl ) ;
00080 
00081   double etTot = 0.0;
00082   for (L1GctEtTotalCollection::const_iterator jbx=sumEtColl->begin(); jbx!=sumEtColl->end(); jbx++) {
00083     if (jbx->bx()==0) { etTot  = static_cast<double>(jbx->et()); }
00084   } 
00085 
00086   double etHad  = 0.0;
00087   for (L1GctEtHadCollection::const_iterator jbx=sumHtColl->begin(); jbx!=sumHtColl->end(); jbx++) {
00088     if (jbx->bx()==0) { etHad  = static_cast<double>(jbx->et()); }
00089   }
00090 
00091   double etMiss = 0.0;
00092   double etMAng = 0.0;
00093   for (L1GctEtMissCollection::const_iterator jbx=missEtColl->begin(); jbx!=missEtColl->end(); jbx++) {
00094     if (jbx->bx()==0) {
00095       etMiss = static_cast<double>(jbx->et());
00096       int phibin = jbx->phi();
00097       if (phibin>=36) phibin -= 72;
00098       double etMPhi = static_cast<double>(phibin);
00099 
00100       etMAng = (etMPhi+0.5)*M_PI/36.;
00101     }
00102   }  
00103 
00104   double etTotFromRegions = 0.0;
00105   double exTotFromRegions = 0.0;
00106   double eyTotFromRegions = 0.0;
00107   for (L1CaloRegionCollection::const_iterator jrg=inputColl->begin(); jrg!=inputColl->end(); jrg++) {
00108     if (jrg->bx()==0) {
00109       double rgEt = static_cast<double>(jrg->et()) * lsbForEt;
00110       double rgPhibin = static_cast<double>(jrg->id().iphi());
00111       double rgPh = (rgPhibin + 0.5)*M_PI/9.;
00112 
00113       etTotFromRegions += rgEt;
00114       exTotFromRegions += rgEt*cos(rgPh);
00115       eyTotFromRegions += rgEt*sin(rgPh);
00116     }
00117   }
00118 
00119   double htMissGct = 0.0;
00120   double htMissAng = 0.0;
00121   double htMissGeV = 0.0;
00122   for (L1GctHtMissCollection::const_iterator jbx=missHtColl->begin(); jbx!=missHtColl->end(); jbx++) {
00123     if (jbx->bx()==0) {
00124       htMissGct = static_cast<double>(jbx->et());
00125       htMissGeV = htMissScale->et(jbx->et());
00126       int phibin = jbx->phi();
00127       if (phibin>=9) phibin -= 18;
00128       double htMPhi = static_cast<double>(phibin);
00129       htMissAng = (htMPhi+0.5)*M_PI/9.;
00130     }
00131   }
00132 
00133   double htFromJets = 0.0;
00134   double hxFromJets = 0.0;
00135   double hyFromJets = 0.0;
00136   for (L1GctInternJetDataCollection::const_iterator jet=internalJetsColl->begin(); jet!=internalJetsColl->end(); jet++) {
00137     if (jet->bx()==0 && !jet->empty()) {
00138       unsigned jetEtGct = jet->et();
00139       double jetEt = static_cast<double>(jetEtGct);
00140       int phibin = jet->regionId().iphi();
00141       if (phibin>=9) phibin -= 18;
00142       // The phi bin centres are at 0, 20, 40, ... degrees
00143       double jetAng = (static_cast<double>(phibin))*M_PI/9.;
00144       if (jetEtGct>httJetThreshold) {
00145         htFromJets += jetEt;
00146       }
00147       if (jetEtGct>htmJetThreshold) {
00148         hxFromJets += jetEt*cos(jetAng);
00149         hyFromJets += jetEt*sin(jetAng);
00150       }
00151     }
00152   }
00153 
00154   double dPhiMetMht = deltaPhi(etMAng,htMissAng);
00155 
00156   theSumEtInLsb->Fill(etTot);
00157   theSumHtInLsb->Fill(etHad);
00158   theMissEtInLsb->Fill(etMiss);
00159   theMissHtInLsb->Fill(htMissGct);
00160   theSumEtInGeV->Fill(etTot*lsbForEt);
00161   theSumHtInGeV->Fill(etHad*lsbForHt);
00162   theMissEtInGeV->Fill(etMiss*lsbForEt);
00163   theMissEtAngle->Fill(etMAng);
00164   theMissEtVector->Fill(etMiss*lsbForEt*cos(etMAng),etMiss*lsbForEt*sin(etMAng));
00165   if (htMissGct<126.5) {
00166     theMissHtInGeV->Fill(htMissGeV);
00167     theMissHtAngle->Fill(htMissAng);
00168     theMissHtVector->Fill(htMissGeV*cos(htMissAng),htMissGeV*sin(htMissAng));
00169   }
00170 
00171   theSumEtVsInputRegions->Fill(etTot*lsbForEt, etTotFromRegions);
00172   theMissEtMagVsInputRegions->Fill(etMiss*lsbForEt, sqrt(exTotFromRegions*exTotFromRegions + eyTotFromRegions*eyTotFromRegions) );
00173   theMissEtAngleVsInputRegions->Fill(etMAng, atan2(-eyTotFromRegions, -exTotFromRegions) );
00174   theMissHtMagVsInputRegions->Fill(htMissGeV, sqrt(exTotFromRegions*exTotFromRegions + eyTotFromRegions*eyTotFromRegions) );
00175 
00176   theMissEtVsMissHt->Fill(etMiss*lsbForEt, htMissGeV);
00177   theMissEtVsMissHtAngle->Fill(etMAng, htMissAng);
00178   theDPhiVsMissEt->Fill(dPhiMetMht,etMiss*lsbForEt);
00179   theDPhiVsMissHt->Fill(dPhiMetMht,htMissGeV);
00180 
00181   theHtVsInternalJetsSum->Fill(etHad*lsbForHt, htFromJets*lsbForHt);
00182   if (htMissGct<126.5) {
00183     theMissHtVsInternalJetsSum->Fill(htMissGeV, sqrt(hxFromJets*hxFromJets + hyFromJets*hyFromJets)*lsbForHt);
00184     theMissHtPhiVsInternalJetsSum->Fill(htMissAng, atan2(-hyFromJets, -hxFromJets));
00185     theMissHxVsInternalJetsSum->Fill(htMissGeV*cos(htMissAng), hxFromJets*lsbForHt);
00186     theMissHyVsInternalJetsSum->Fill(htMissGeV*sin(htMissAng), hyFromJets*lsbForHt);
00187   }
00188 
00189   // Get minbias trigger quantities from HF
00190   Handle<L1GctHFRingEtSumsCollection> HFEtSumsColl;
00191   Handle<L1GctHFBitCountsCollection>  HFCountsColl;
00192   iEvent.getByLabel( m_energy_tag, HFEtSumsColl ) ;
00193   iEvent.getByLabel( m_energy_tag, HFCountsColl ) ;
00194 
00195   for (L1GctHFRingEtSumsCollection::const_iterator es = HFEtSumsColl->begin(); es != HFEtSumsColl->end(); es++) {
00196     if (es->bx()==0) {
00197       theHfRing0EtSumPositiveEta->Fill(hfRingEtScale->et(es->etSum(0)));  
00198       theHfRing0EtSumNegativeEta->Fill(hfRingEtScale->et(es->etSum(1)));  
00199       theHfRing1EtSumPositiveEta->Fill(hfRingEtScale->et(es->etSum(2)));  
00200       theHfRing1EtSumNegativeEta->Fill(hfRingEtScale->et(es->etSum(3)));
00201     }
00202   }  
00203 
00204   for (L1GctHFBitCountsCollection::const_iterator bc = HFCountsColl->begin(); bc != HFCountsColl->end(); bc++) {
00205     if (bc->bx()==0) {
00206       theHfRing0CountPositiveEta->Fill(bc->bitCount(0));  
00207       theHfRing0CountNegativeEta->Fill(bc->bitCount(1));  
00208       theHfRing1CountPositiveEta->Fill(bc->bitCount(2));  
00209       theHfRing1CountNegativeEta->Fill(bc->bitCount(3));
00210     }
00211   }
00212 
00213 }
00214 
00215 // ------------ method called once each job just before starting event loop  ------------
00216 void 
00217 L1GctValidation::beginJob()
00218 {
00219   edm::Service<TFileService> fs;
00220 
00221   TFileDirectory dir0 = fs->mkdir("L1GctEtSums");
00222 
00223   theSumEtInLsb   = dir0.make<TH1F>("SumEtInLsb",   "Total Et (GCT units)",
00224                                     128, 0., 2048.); 
00225   theSumHtInLsb   = dir0.make<TH1F>("SumHtInLsb",   "Total Ht (GCT units)",
00226                                     128, 0., 2048.); 
00227   theMissEtInLsb  = dir0.make<TH1F>("MissEtInLsb",  "Missing Et magnitude (GCT units)",
00228                                     128, 0., 1024.); 
00229   theMissHtInLsb  = dir0.make<TH1F>("MissHtInLsb",  "Missing Ht magnitude (GCT units)",
00230                                     128, 0., 127.); 
00231   theSumEtInGeV   = dir0.make<TH1F>("SumEtInGeV",   "Total Et (in GeV)",
00232                                     100, 0., 1000.); 
00233   theSumHtInGeV   = dir0.make<TH1F>("SumHtInGeV",   "Total Ht (in GeV)",
00234                                     100, 0., 1000.); 
00235   theMissEtInGeV  = dir0.make<TH1F>("MissEtInGeV",  "Missing Et magnitude (in GeV)",
00236                                     100, 0., 500.); 
00237   theMissEtAngle  = dir0.make<TH1F>("MissEtAngle",  "Missing Et angle",
00238                                     72, -M_PI, M_PI);
00239   theMissEtVector = dir0.make<TH2F>("MissEtVector", "Missing Ex vs Missing Ey",
00240                                     100, -100., 100., 100, -100., 100.); 
00241   theMissHtInGeV  = dir0.make<TH1F>("MissHtInGeV",  "Missing Ht magnitude (in GeV)",
00242                                     100, 0., 500.); 
00243   theMissHtAngle  = dir0.make<TH1F>("MissHtAngle",  "Missing Ht angle",
00244                                     72, -M_PI, M_PI);
00245   theMissHtVector = dir0.make<TH2F>("MissHtVector", "Missing Hx vs Missing Hy",
00246                                     100, -100., 100., 100, -100., 100.); 
00247   theSumEtVsInputRegions = dir0.make<TH2F>("SumEtVsInputRegions", "Total Et vs sum of input regions",
00248                                            100, 0., 1000., 100, 0., 1000.);
00249   theMissEtMagVsInputRegions = dir0.make<TH2F>("MissEtMagVsInputRegions", "Missing Et magnitude vs sum of input regions",
00250                                                100, 0., 500., 100, 0., 500.);
00251   theMissEtAngleVsInputRegions = dir0.make<TH2F>("MissEtAngleVsInputRegions", "Missing Et angle vs sum of input regions",
00252                                                 72, -M_PI, M_PI, 72, -M_PI, M_PI);
00253   theMissHtMagVsInputRegions = dir0.make<TH2F>("MissHtMagVsInputRegions", "Missing Ht magnitude vs sum of input regions",
00254                                                100, 0., 500., 100, 0., 500.);
00255   theMissEtVsMissHt = dir0.make<TH2F>("MissEtVsMissHt", "Missing Et vs Missing Ht",
00256                                       100, 0., 500., 100, 0., 500.);
00257   theMissEtVsMissHtAngle = dir0.make<TH2F>("MissEtVsMissHtAngle", "Angle correlation Missing Et vs Missing Ht",
00258                                            72, -M_PI, M_PI, 72, -M_PI, M_PI);
00259   theDPhiVsMissEt = dir0.make<TH2F>("theDPhiVsMissEt", "Angle difference MET-MHT vs MET magnitude",
00260                                     72, -M_PI, M_PI, 100, 0., 500.);
00261   theDPhiVsMissHt = dir0.make<TH2F>("theDPhiVsMissHt", "Angle difference MET-MHT vs MHT magnitude",
00262                                     72, -M_PI, M_PI, 100, 0., 500.);
00263 
00264   theHtVsInternalJetsSum     = dir0.make<TH2F>("HtVsInternalJetsSum", "Ht vs scalar sum of jet Et values (in GeV)",
00265                                                128, 0., 2048., 128, 0., 2048.);
00266   theMissHtVsInternalJetsSum = dir0.make<TH2F>("MissHtVsInternalJetsSum", "Missing Ht vs vector sum of jet Et values (in GeV)",
00267                                                128, 0., 512., 128, 0., 512.);
00268   theMissHtPhiVsInternalJetsSum = dir0.make<TH2F>("MissHtPhiVsInternalJetsSum", "Angle correlation Missing Ht vs vector sum of jet Et values",
00269                                                72, -M_PI, M_PI, 72, -M_PI, M_PI);
00270   theMissHxVsInternalJetsSum = dir0.make<TH2F>("MissHxVsInternalJetsSum", "Missing Ht x component vs sum of jet Et values (in GeV)",
00271                                                128, -256., 256., 128, -256., 256.);
00272   theMissHyVsInternalJetsSum = dir0.make<TH2F>("MissHyVsInternalJetsSum", "Missing Ht y component vs sum of jet Et values (in GeV)",
00273                                                128, -256., 256., 128, -256., 256.);
00274 
00275 
00276   TFileDirectory dir1 = fs->mkdir("L1GctHfSumsAndJetCounts");
00277 
00278   // Minimum bias triggers from Hf inner rings
00279   theHfRing0EtSumPositiveEta = dir1.make<TH1F>("HfRing0EtSumPositiveEta", "Hf Inner Ring0 Et eta+",
00280                                                60, 0., 30.);
00281   theHfRing0EtSumNegativeEta = dir1.make<TH1F>("HfRing0EtSumNegativeEta", "Hf Inner Ring0 Et eta-",
00282                                                60, 0., 30.);
00283   theHfRing1EtSumPositiveEta = dir1.make<TH1F>("HfRing1EtSumPositiveEta", "Hf Inner Ring1 Et eta+",
00284                                                60, 0., 30.);
00285   theHfRing1EtSumNegativeEta = dir1.make<TH1F>("HfRing1EtSumNegativeEta", "Hf Inner Ring1 Et eta-",
00286                                                60, 0., 30.);
00287   theHfRing0CountPositiveEta = dir1.make<TH1F>("HfRing0CountPositiveEta", "Hf Threshold bits Ring0 eta+",
00288                                                20, 0., 20.);
00289   theHfRing0CountNegativeEta = dir1.make<TH1F>("HfRing0CountNegativeEta", "Hf Threshold bits Ring0 eta-",
00290                                                20, 0., 20.);
00291   theHfRing1CountPositiveEta = dir1.make<TH1F>("HfRing1CountPositiveEta", "Hf Threshold bits Ring1 eta+",
00292                                                20, 0., 20.);
00293   theHfRing1CountNegativeEta = dir1.make<TH1F>("HfRing1CountNegativeEta", "Hf Threshold bits Ring1 eta-",
00294                                                20, 0., 20.);
00295 }
00296 
00297 // ------------ method called once each job just after ending the event loop  ------------
00298 void 
00299 L1GctValidation::endJob() {
00300 }
00301 
00302 DEFINE_FWK_MODULE(L1GctValidation);
00303
/data/refman/pasoursint/CMSSW_4_1_8_patch9/src/L1Trigger/GlobalCaloTrigger/plugins/L1GctValidation.cc
