Physics/src/EwkDQM.cc

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  \author Valentina Gori, University of Firenze
0005  */
0006
0007 #include "DQM/Physics/src/EwkDQM.h"
0008
0009 #include <vector>
0010 #include <string>
0011 #include <cmath>
0012
0013 #include "DataFormats/Candidate/interface/Candidate.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/EventSetup.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017
0018 #include "DQMServices/Core/interface/DQMStore.h"
0019 #include "FWCore/ServiceRegistry/interface/Service.h"
0020
0021 #include "DataFormats/Common/interface/Handle.h"
0022
0023 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0024 #include "FWCore/Common/interface/TriggerNames.h"
0025
0026 // Physics Objects
0027 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0028 #include "DataFormats/MuonReco/interface/Muon.h"
0029 #include "DataFormats/JetReco/interface/Jet.h"
0030 #include "DataFormats/METReco/interface/MET.h"
0031 #include "DataFormats/TrackReco/interface/Track.h"
0032 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0033 #include "DataFormats/VertexReco/interface/Vertex.h"
0034 #include "DataFormats/Math/interface/LorentzVector.h"
0035
0036 #include "TLorentzVector.h"
0037
0038 using namespace std;
0039 using namespace edm;
0040 using namespace reco;
0041
0042 // EwkDQM::EwkDQM(const ParameterSet& parameters) {
0043 EwkDQM::EwkDQM(const ParameterSet& parameters) {
0044   eJetMin_ = parameters.getUntrackedParameter<double>("EJetMin", 999999.);
0045
0046   // riguardare questa sintassi
0047   // Get parameters from configuration file
0048   thePFJetCollectionLabel_ = parameters.getParameter<InputTag>("PFJetCollection");
0049   theCaloMETCollectionLabel_ = parameters.getParameter<InputTag>("caloMETCollection");
0050   theTriggerResultsCollection_ = parameters.getParameter<InputTag>("triggerResultsCollection");
0051
0052   theElecTriggerPathToPass_ = parameters.getParameter<std::vector<string> >("elecTriggerPathToPass");
0053   theMuonTriggerPathToPass_ = parameters.getParameter<std::vector<string> >("muonTriggerPathToPass");
0054   theTriggerResultsToken_ =
0055       consumes<edm::TriggerResults>(parameters.getParameter<InputTag>("triggerResultsCollection"));
0056   theMuonCollectionLabel_ = consumes<reco::MuonCollection>(parameters.getParameter<InputTag>("muonCollection"));
0057   theElectronCollectionLabel_ =
0058       consumes<reco::GsfElectronCollection>(parameters.getParameter<InputTag>("electronCollection"));
0059   thePFJetCollectionToken_ = consumes<edm::View<reco::Jet> >(parameters.getParameter<InputTag>("PFJetCollection"));
0060   theCaloMETCollectionToken_ = consumes<edm::View<reco::MET> >(parameters.getParameter<InputTag>("caloMETCollection"));
0061   theVertexToken_ = consumes<reco::VertexCollection>(parameters.getParameter<InputTag>("vertexCollection"));
0062
0063   // just to initialize
0064   isValidHltConfig_ = false;
0065
0066   h_vertex_number = nullptr;
0067   h_vertex_chi2 = nullptr;
0068   h_vertex_numTrks = nullptr;
0069   h_vertex_sumTrks = nullptr;
0070   h_vertex_d0 = nullptr;
0071
0072   h_jet_count = nullptr;
0073   h_jet_et = nullptr;
0074   h_jet_pt = nullptr;
0075   h_jet_eta = nullptr;
0076   h_jet_phi = nullptr;
0077   h_jet2_et = nullptr;
0078   // h_jet2_pt = 0;
0079   h_jet2_eta = nullptr;
0080   h_jet2_phi = nullptr;
0081
0082   h_e1_et = nullptr;
0083   h_e2_et = nullptr;
0084   h_e1_eta = nullptr;
0085   h_e2_eta = nullptr;
0086   h_e1_phi = nullptr;
0087   h_e2_phi = nullptr;
0088
0089   h_m1_pt = nullptr;
0090   h_m2_pt = nullptr;
0091   h_m1_eta = nullptr;
0092   h_m2_eta = nullptr;
0093   h_m1_phi = nullptr;
0094   h_m2_phi = nullptr;
0095
0096   // h_t1_et = 0;
0097   // h_t1_eta = 0;
0098   // h_t1_phi = 0;
0099
0100   h_met = nullptr;
0101   h_met_phi = nullptr;
0102
0103   h_e_invWMass = nullptr;
0104   h_m_invWMass = nullptr;
0105   h_mumu_invMass = nullptr;
0106   h_ee_invMass = nullptr;
0107 }
0108
0109 EwkDQM::~EwkDQM() {}
0110
0111 void EwkDQM::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const&, edm::EventSetup const&) {
0112   ibooker.setCurrentFolder("Physics/EwkDQM");
0113
0114   char chtitle[256] = "";
0115   const size_t title_s = sizeof(chtitle);
0116
0117   logTraceName = "EwkAnalyzer";
0118
0119   LogTrace(logTraceName) << "Parameters initialization";
0120
0121   const float pi = 4 * atan(1);
0122
0123   // Keep the number of plots and number of bins to a minimum!
0124
0125   h_vertex_number = ibooker.book1D("vertex_number", "Number of event vertices in collection", 10, -0.5, 9.5);
0126   h_vertex_chi2 = ibooker.book1D("vertex_chi2", "Event Vertex #chi^{2}/n.d.o.f.", 20, 0.0, 2.0);
0127   h_vertex_numTrks = ibooker.book1D("vertex_numTrks", "Event Vertex, number of tracks", 20, -0.5, 59.5);
0128   h_vertex_sumTrks = ibooker.book1D("vertex_sumTrks", "Event Vertex, sum of track pt", 20, 0.0, 100.0);
0129   h_vertex_d0 = ibooker.book1D("vertex_d0", "Event Vertex d0", 20, 0.0, 0.05);
0130
0131   h_jet_count = ibooker.book1D("jet_count", chtitle, 8, -0.5, 7.5);
0132
0133   snprintf(
0134       chtitle, title_s, "Leading jet E_{T} (from %s);E_{T}(1^{st} jet) (GeV)", thePFJetCollectionLabel_.label().data());
0135   h_jet_et = ibooker.book1D("jet_et", chtitle, 20, 0., 200.0);
0136
0137   snprintf(chtitle,
0138            title_s,
0139            "Leading jet p_{T} (from %s);p_{T}(1^{st} jet) (GeV/c)",
0140            thePFJetCollectionLabel_.label().data());
0141   h_jet_pt = ibooker.book1D("jet_pt", chtitle, 20, 0., 200.0);
0142
0143   snprintf(chtitle, title_s, "Leading jet #eta (from %s); #eta (1^{st} jet)", thePFJetCollectionLabel_.label().data());
0144   h_jet_eta = ibooker.book1D("jet_eta", chtitle, 20, -10., 10.0);
0145
0146   snprintf(chtitle, title_s, "Leading jet #phi (from %s); #phi(1^{st} jet)", thePFJetCollectionLabel_.label().data());
0147   h_jet_phi = ibooker.book1D("jet_phi", chtitle, 22, -1.1 * pi, 1.1 * pi);
0148
0149   snprintf(chtitle,
0150            title_s,
0151            "2^{nd} leading jet E_{T} (from %s);E_{T}(2^{nd} jet) (GeV)",
0152            thePFJetCollectionLabel_.label().data());
0153   h_jet2_et = ibooker.book1D("jet2_et", chtitle, 20, 0., 200.0);
0154
0155   snprintf(chtitle,
0156            title_s,
0157            "2^{nd} leading jet #eta (from %s); #eta (2^{nd} jet)",
0158            thePFJetCollectionLabel_.label().data());
0159   h_jet2_eta = ibooker.book1D("jet2_eta", chtitle, 20, -10., 10.0);
0160
0161   snprintf(
0162       chtitle, title_s, "2^{nd} leading jet #phi (from %s); #phi(2^{nd} jet)", thePFJetCollectionLabel_.label().data());
0163   h_jet2_phi = ibooker.book1D("jet2_phi", chtitle, 22, -1.1 * pi, 1.1 * pi);
0164
0165   h_e1_et = ibooker.book1D("e1_et", "E_{T} of Leading Electron;E_{T} (GeV)", 20, 0.0, 100.0);
0166   h_e2_et = ibooker.book1D("e2_et", "E_{T} of Second Electron;E_{T} (GeV)", 20, 0.0, 100.0);
0167   h_e1_eta = ibooker.book1D("e1_eta", "#eta of Leading Electron;#eta", 20, -4.0, 4.0);
0168
0169   h_e2_eta = ibooker.book1D("e2_eta", "#eta of Second Electron;#eta", 20, -4.0, 4.0);
0170
0171   h_e1_phi = ibooker.book1D("e1_phi", "#phi of Leading Electron;#phi", 22, -1.1 * pi, 1.1 * pi);
0172   h_e2_phi = ibooker.book1D("e2_phi", "#phi of Second Electron;#phi", 22, -1.1 * pi, 1.1 * pi);
0173   h_m1_pt = ibooker.book1D("m1_pt", "p_{T} of Leading Muon;p_{T}(1^{st} #mu) (GeV)", 20, 0.0, 100.0);
0174   h_m2_pt = ibooker.book1D("m2_pt", "p_{T} of Second Muon;p_{T}(2^{nd} #mu) (GeV)", 20, 0.0, 100.0);
0175   h_m1_eta = ibooker.book1D("m1_eta", "#eta of Leading Muon;#eta(1^{st} #mu)", 20, -4.0, 4.0);
0176   h_m2_eta = ibooker.book1D("m2_eta", "#eta of Second Muon;#eta(2^{nd} #mu)", 20, -4.0, 4.0);
0177   h_m1_phi = ibooker.book1D(
0178       "m1_phi", "#phi of Leading Muon;#phi(1^{st} #mu)", 20, (-1. - 1. / 10.) * pi, (1. + 1. / 10.) * pi);
0179   h_m2_phi =
0180       ibooker.book1D("m2_phi", "#phi of Second Muon;#phi(2^{nd} #mu)", 20, (-1. - 1. / 10.) * pi, (1. + 1. / 10.) * pi);
0181
0182   snprintf(chtitle, title_s, "Missing E_{T} (%s); GeV", theCaloMETCollectionLabel_.label().data());
0183   h_met = ibooker.book1D("met", chtitle, 20, 0.0, 100);
0184   h_met_phi =
0185       ibooker.book1D("met_phi", "Missing E_{T} #phi;#phi(MET)", 22, (-1. - 1. / 10.) * pi, (1. + 1. / 10.) * pi);
0186
0187   h_e_invWMass = ibooker.book1D("we_invWMass", "W-> e #nu Transverse Mass;M_{T} (GeV)", 20, 0.0, 140.0);
0188   h_m_invWMass = ibooker.book1D("wm_invWMass", "W-> #mu #nu Transverse Mass;M_{T} (GeV)", 20, 0.0, 140.0);
0189   h_mumu_invMass = ibooker.book1D("z_mm_invMass", "#mu#mu Invariant Mass;InvMass (GeV)", 20, 40.0, 140.0);
0190   h_ee_invMass = ibooker.book1D("z_ee_invMass", "ee Invariant Mass;InvMass (Gev)", 20, 40.0, 140.0);
0191 }
0192
0193 ///
0194 void EwkDQM::dqmBeginRun(const edm::Run& theRun, const edm::EventSetup& theSetup) {
0195   // passed as parameter to HLTConfigProvider::init(), not yet used
0196   bool isConfigChanged = false;
0197
0198   // isValidHltConfig_ used to short-circuit analyze() in case of problems
0199   const std::string hltProcessName(theTriggerResultsCollection_.process());
0200   isValidHltConfig_ = hltConfigProvider_.init(theRun, theSetup, hltProcessName, isConfigChanged);
0201 }
0202
0203 void EwkDQM::analyze(const Event& iEvent, const EventSetup& iSetup) {
0204   // short-circuit if hlt problems
0205   if (!isValidHltConfig_)
0206     return;
0207
0208   LogTrace(logTraceName) << "Analysis of event # ";
0209   // Did it pass certain HLT path?
0210   Handle<TriggerResults> HLTresults;
0211   iEvent.getByToken(theTriggerResultsToken_, HLTresults);
0212   if (!HLTresults.isValid())
0213     return;
0214
0215   const edm::TriggerNames& trigNames = iEvent.triggerNames(*HLTresults);
0216
0217   // a temporary, until we have a list of triggers of interest
0218   std::vector<std::string> eleTrigPathNames;
0219   std::vector<std::string> muTrigPathNames;
0220
0221   // eleTrigPathNames.push_back(theElecTriggerPathToPass_);
0222   // muTrigPathNames.push_back(theMuonTriggerPathToPass_);
0223   // end of temporary
0224
0225   bool passed_electron_HLT = false;
0226   bool passed_muon_HLT = false;
0227   for (unsigned int i = 0; i < HLTresults->size(); i++) {
0228     const std::string& trigName = trigNames.triggerName(i);
0229     // check if triggerName matches electronPath
0230     for (unsigned int index = 0; index < theElecTriggerPathToPass_.size() && !passed_electron_HLT; index++) {
0231       // 0 if found, pos if not
0232       size_t trigPath = trigName.find(theElecTriggerPathToPass_[index]);
0233       if (trigPath == 0) {
0234         //      cout << "MuonTrigger passed (=trigName): " << trigName <<endl;
0235         passed_electron_HLT = HLTresults->accept(i);
0236       }
0237     }
0238     // check if triggerName matches muonPath
0239     for (unsigned int index = 0; index < theMuonTriggerPathToPass_.size() && !passed_muon_HLT; index++) {
0240       // 0 if found, pos if not
0241       size_t trigPath = trigName.find(theMuonTriggerPathToPass_[index]);
0242       if (trigPath == 0) {
0243         //      cout << "MuonTrigger passed (=trigName): " << trigName <<endl;
0244         passed_muon_HLT = HLTresults->accept(i);
0245       }
0246     }
0247   }
0248
0249   // we are interested in events with a valid electron or muon
0250   if (!(passed_electron_HLT || passed_muon_HLT))
0251     return;
0252
0253   ////////////////////////////////////////////////////////////////////////////////
0254   // Vertex information
0255   Handle<VertexCollection> vertexHandle;
0256   iEvent.getByToken(theVertexToken_, vertexHandle);
0257   if (!vertexHandle.isValid())
0258     return;
0259   VertexCollection vertexCollection = *(vertexHandle.product());
0260   VertexCollection::const_iterator v = vertexCollection.begin();
0261   int vertex_number = vertexCollection.size();
0262   double vertex_chi2 = v->normalizedChi2();  // v->chi2();
0263   double vertex_d0 = sqrt(v->x() * v->x() + v->y() * v->y());
0264   double vertex_numTrks = v->tracksSize();
0265   double vertex_sumTrks = 0.0;
0266   // std::cout << "vertex_d0=" << vertex_d0 << "\n";
0267   // double vertex_ndof    = v->ndof();cout << "ndof="<<vertex_ndof<<endl;
0268   for (Vertex::trackRef_iterator vertex_curTrack = v->tracks_begin(); vertex_curTrack != v->tracks_end();
0269        vertex_curTrack++)
0270     vertex_sumTrks += (*vertex_curTrack)->pt();
0271
0272   ////////////////////////////////////////////////////////////////////////////////
0273   // Missing ET
0274   Handle<View<MET> > caloMETCollection;
0275   iEvent.getByToken(theCaloMETCollectionToken_, caloMETCollection);
0276   if (!caloMETCollection.isValid())
0277     return;
0278   float missing_et = caloMETCollection->begin()->et();
0279   float met_phi = caloMETCollection->begin()->phi();
0280
0281   ////////////////////////////////////////////////////////////////////////////////
0282   // grab "gaussian sum fitting" electrons
0283   Handle<GsfElectronCollection> electronCollection;
0284   iEvent.getByToken(theElectronCollectionLabel_, electronCollection);
0285   if (!electronCollection.isValid())
0286     return;
0287
0288   // Find the highest and 2nd highest electron
0289   float electron_et = -8.0;
0290   float electron_eta = -8.0;
0291   float electron_phi = -8.0;
0292   float electron2_et = -9.0;
0293   float electron2_eta = -9.0;
0294   float electron2_phi = -9.0;
0295   float ee_invMass = -9.0;
0296   TLorentzVector e1, e2;
0297
0298   // If it passed electron HLT and the collection was found, find electrons near
0299   // Z mass
0300   if (passed_electron_HLT) {
0301     for (reco::GsfElectronCollection::const_iterator recoElectron = electronCollection->begin();
0302          recoElectron != electronCollection->end();
0303          recoElectron++) {
0304       // Require electron to pass some basic cuts
0305       if (recoElectron->et() < 20 || fabs(recoElectron->eta()) > 2.5)
0306         continue;
0307
0308       // Tighter electron cuts
0309       if (recoElectron->deltaPhiSuperClusterTrackAtVtx() > 0.58 ||
0310           recoElectron->deltaEtaSuperClusterTrackAtVtx() > 0.01 || recoElectron->sigmaIetaIeta() > 0.027)
0311         continue;
0312
0313       if (recoElectron->et() > electron_et) {
0314         electron2_et = electron_et;  // 2nd highest gets values from current highest
0315         electron2_eta = electron_eta;
0316         electron2_phi = electron_phi;
0317         electron_et = recoElectron->et();  // 1st highest gets values from new highest
0318         electron_eta = recoElectron->eta();
0319         electron_phi = recoElectron->phi();
0320         e1 = TLorentzVector(recoElectron->momentum().x(),
0321                             recoElectron->momentum().y(),
0322                             recoElectron->momentum().z(),
0323                             recoElectron->p());
0324       } else if (recoElectron->et() > electron2_et) {
0325         electron2_et = recoElectron->et();
0326         electron2_eta = recoElectron->eta();
0327         electron2_phi = recoElectron->phi();
0328         e2 = TLorentzVector(recoElectron->momentum().x(),
0329                             recoElectron->momentum().y(),
0330                             recoElectron->momentum().z(),
0331                             recoElectron->p());
0332       }
0333     }  // end of loop over electrons
0334     if (electron2_et > 0.0) {
0335       TLorentzVector pair = e1 + e2;
0336       ee_invMass = pair.M();
0337     }
0338   }  // end of "are electrons valid"
0339   ////////////////////////////////////////////////////////////////////////////////
0340
0341   ////////////////////////////////////////////////////////////////////////////////
0342   // Take the STA muon container
0343   Handle<MuonCollection> muonCollection;
0344   iEvent.getByToken(theMuonCollectionLabel_, muonCollection);
0345   if (!muonCollection.isValid())
0346     return;
0347
0348   // Find the highest pt muons
0349   float mm_invMass = -9.0;
0350   float muon_pt = -9.0;
0351   float muon_eta = -9.0;
0352   float muon_phi = -9.0;
0353   float muon2_pt = -9.0;
0354   float muon2_eta = -9.0;
0355   float muon2_phi = -9.0;
0356   TLorentzVector m1, m2;
0357
0358   if (passed_muon_HLT) {
0359     for (reco::MuonCollection::const_iterator recoMuon = muonCollection->begin(); recoMuon != muonCollection->end();
0360          recoMuon++) {
0361       // Require muon to pass some basic cuts
0362       if (recoMuon->pt() < 20 || !recoMuon->isGlobalMuon())
0363         continue;
0364       // Some tighter muon cuts
0365       if (recoMuon->globalTrack()->normalizedChi2() > 10)
0366         continue;
0367
0368       if (recoMuon->pt() > muon_pt) {
0369         muon2_pt = muon_pt;  // 2nd highest gets values from current highest
0370         muon2_eta = muon_eta;
0371         muon2_phi = muon_phi;
0372         muon_pt = recoMuon->pt();  // 1st highest gets values from new highest
0373         muon_eta = recoMuon->eta();
0374         muon_phi = recoMuon->phi();
0375         m1 =
0376             TLorentzVector(recoMuon->momentum().x(), recoMuon->momentum().y(), recoMuon->momentum().z(), recoMuon->p());
0377       } else if (recoMuon->pt() > muon2_pt) {
0378         muon2_pt = recoMuon->pt();
0379         muon2_eta = recoMuon->eta();
0380         muon2_phi = recoMuon->phi();
0381         m2 =
0382             TLorentzVector(recoMuon->momentum().x(), recoMuon->momentum().y(), recoMuon->momentum().z(), recoMuon->p());
0383       }
0384     }
0385   }
0386   if (muon2_pt > 0.0) {
0387     TLorentzVector pair = m1 + m2;
0388     mm_invMass = pair.M();
0389   }
0390   ////////////////////////////////////////////////////////////////////////////////
0391
0392   ////////////////////////////////////////////////////////////////////////////////
0393   // Find the highest et jet
0394
0395   //  Handle<CaloJetCollection> caloJetCollection;
0396   Handle<View<Jet> > PFJetCollection;
0397   iEvent.getByToken(thePFJetCollectionToken_, PFJetCollection);
0398   if (!PFJetCollection.isValid())
0399     return;
0400
0401   unsigned int muonCollectionSize = muonCollection->size();
0402   // unsigned int jetCollectionSize = jetCollection->size();
0403   unsigned int PFJetCollectionSize = PFJetCollection->size();
0404   int jet_count = 0;
0405   // int LEADJET=-1;  double max_pt=0;
0406
0407   float jet_et = -80.0;
0408   float jet_pt = -80.0;   // prova
0409   float jet_eta = -80.0;  // now USED
0410   float jet_phi = -80.0;  // now USED
0411   float jet2_et = -90.0;
0412   float jet2_eta = -90.0;  // now USED
0413   float jet2_phi = -90.0;  // now USED
0414   //  for (CaloJetCollection::const_iterator i_calojet =
0415   // caloJetCollection->begin();
0416   //       i_calojet != caloJetCollection->end(); i_calojet++) {
0417   //  for (PFJetCollection::const_iterator i_pfjet = PFJetCollection->begin();
0418   //       i_pfjet != PFJetCollection->end(); i_pfjet++) {
0419   //  float jet_current_et = i_calojet->et();
0420   //  float jet_current_et = i_pfjet->et();            // e` identico a jet.et()
0421   //    jet_count++;
0422
0423   // cleaning: va messo prima del riempimento dell'istogramma // This is in
0424   // order to use PFJets
0425   for (unsigned int i = 0; i < PFJetCollectionSize; i++) {
0426     const Jet& jet = PFJetCollection->at(i);
0427     // la classe "jet" viene definita qui!!!
0428     double minDistance = 99999;
0429     for (unsigned int j = 0; j < muonCollectionSize; j++) {
0430       const Muon& mu = muonCollection->at(j);
0431       double distance =
0432           sqrt((mu.eta() - jet.eta()) * (mu.eta() - jet.eta()) + (mu.phi() - jet.phi()) * (mu.phi() - jet.phi()));
0433       if (minDistance > distance)
0434         minDistance = distance;
0435     }
0436     if (minDistance < 0.3)
0437       continue;  // 0.3 is the isolation cone around the muon
0438     // se la distanza muone-cono del jet e` minore di 0.3, passo avanti e non
0439     // conteggio il mio jet
0440
0441     // If it overlaps with ELECTRON, it is not a jet
0442     if (electron_et > 0.0 && fabs(jet.eta() - electron_eta) < 0.2 && calcDeltaPhi(jet.phi(), electron_phi) < 0.2)
0443       continue;
0444     if (electron2_et > 0.0 && fabs(jet.eta() - electron2_eta) < 0.2 && calcDeltaPhi(jet.phi(), electron2_phi) < 0.2)
0445       continue;
0446
0447     // provo a cambiare la parte degli elettroni in modo simmetrico alla parte
0448     // per i muoni
0449
0450     // ...
0451     // ...
0452
0453     // if it has too low Et, throw away
0454     if (jet.et() < eJetMin_)
0455       continue;
0456     jet_count++;
0457
0458     // ovvero: incrementa jet_count se:
0459     //   - non c'e un muone entro 0.3 di distanza dal cono del jet;
0460     //   - se il jet non si sovrappone ad un elettrone;
0461     //   - se l'energia trasversa e` maggiore della soglia impostata (15?)
0462
0463     // if(jet.et()>max_pt) { LEADJET=i; max_pt=jet.et();}
0464     // se l'energia del jet e` maggiore di max_pt, diventa "i"
0465     // l'indice del jet piu` energetico e max_pt la sua energia
0466
0467     // riguardare questo!!!
0468     // fino ad ora, jet_et era inizializzato a -8.0
0469     if (jet.et() > jet_et) {
0470       jet2_et = jet_et;  // 2nd highest jet gets et from current highest
0471       // perche` prende l'energia del primo jet??
0472       jet2_eta = jet_eta;  // now USED
0473       jet2_phi = jet_phi;  // now USED
0474       // jet_et = i_calojet->et(); // current highest jet gets
0475       // et from the new highest
0476       jet_et = jet.et();  // current highest jet gets et from the new highest
0477       // ah, ok! lo riaggiorna solo dopo!
0478       jet_pt = jet.pt();                          // e` il pT del leading jet
0479       jet_eta = jet.eta();                        // now USED
0480       jet_phi = jet.phi() * (Geom::pi() / 180.);  // now USED
0481     } else if (jet.et() > jet2_et) {
0482       //      jet2_et  = i_calojet->et();
0483       jet2_et = jet.et();
0484       //      jet2_eta = i_calojet->eta();  // UNUSED
0485       //      jet2_phi = i_calojet->phi();  // UNUSED
0486       jet2_eta = jet.eta();  // now USED
0487       jet2_phi = jet.phi();  // now USED
0488     }
0489     // questo elseif funziona
0490   }
0491   ////////////////////////////////////////////////////////////////////////////////
0492
0493   ////////////////////////////////////////////////////////////////////////////////
0494   //                 Fill Histograms //
0495   ////////////////////////////////////////////////////////////////////////////////
0496
0497   bool fill_e1 = false;
0498   bool fill_e2 = false;
0499   bool fill_m1 = false;
0500   bool fill_m2 = false;
0501   bool fill_met = false;
0502
0503   // Was Z->ee found?
0504   if (ee_invMass > 0.0) {
0505     h_ee_invMass->Fill(ee_invMass);
0506     fill_e1 = true;
0507     fill_e2 = true;
0508   }
0509
0510   // Was Z->mu mu found?
0511   if (mm_invMass > 0.0) {
0512     h_mumu_invMass->Fill(mm_invMass);
0513     fill_m1 = true;
0514     fill_m2 = true;
0515     h_jet2_et->Fill(jet2_et);
0516   }
0517
0518   // Was W->e nu found?
0519   if (electron_et > 0.0 && missing_et > 20.0) {
0520     float dphiW = fabs(met_phi - electron_phi);
0521     float W_mt_e = sqrt(2 * missing_et * electron_et * (1 - cos(dphiW)));
0522     h_e_invWMass->Fill(W_mt_e);
0523     fill_e1 = true;
0524     fill_met = true;
0525   }
0526
0527   // Was W->mu nu found?
0528   if (muon_pt > 0.0 && missing_et > 20.0) {
0529     float dphiW = fabs(met_phi - muon_phi);
0530     float W_mt_m = sqrt(2 * missing_et * muon_pt * (1 - cos(dphiW)));
0531     h_m_invWMass->Fill(W_mt_m);
0532     fill_m1 = true;
0533     fill_met = true;
0534   }
0535
0536   if (jet_et > -10.0) {
0537     h_jet_et->Fill(jet_et);
0538     h_jet_count->Fill(jet_count);
0539   }
0540
0541   if (jet_pt > 0.) {
0542     h_jet_pt->Fill(jet_pt);
0543   }
0544
0545   if (jet_eta > -50.) {
0546     h_jet_eta->Fill(jet_eta);
0547   }
0548
0549   if (jet_phi > -10.) {
0550     h_jet_phi->Fill(jet_phi);
0551   }
0552
0553   if (jet2_et > -10.0) {
0554     h_jet2_et->Fill(jet2_et);
0555   }
0556
0557   // if (jet2_pt>0.) {
0558   //  h_jet2_pt   ->Fill(jet2_pt);
0559   // }
0560
0561   if (jet2_eta > -50.) {
0562     h_jet2_eta->Fill(jet2_eta);
0563   }
0564
0565   if (jet2_phi > -10.) {
0566     h_jet2_phi->Fill(jet2_phi);
0567   }
0568
0569   if (fill_e1 || fill_m1) {
0570     h_vertex_number->Fill(vertex_number);
0571     h_vertex_chi2->Fill(vertex_chi2);
0572     h_vertex_d0->Fill(vertex_d0);
0573     h_vertex_numTrks->Fill(vertex_numTrks);
0574     h_vertex_sumTrks->Fill(vertex_sumTrks);
0575   }
0576
0577   if (fill_e1) {
0578     h_e1_et->Fill(electron_et);
0579     h_e1_eta->Fill(electron_eta);
0580     h_e1_phi->Fill(electron_phi);
0581   }
0582   if (fill_e2) {
0583     h_e2_et->Fill(electron2_et);
0584     h_e2_eta->Fill(electron2_eta);
0585     h_e2_phi->Fill(electron2_phi);
0586   }
0587   if (fill_m1) {
0588     h_m1_pt->Fill(muon_pt);
0589     h_m1_eta->Fill(muon_eta);
0590     h_m1_phi->Fill(muon_phi);
0591   }
0592   if (fill_m2) {
0593     h_m2_pt->Fill(muon2_pt);
0594     h_m2_eta->Fill(muon2_eta);
0595     h_m2_phi->Fill(muon2_phi);
0596   }
0597   if (fill_met) {
0598     h_met->Fill(missing_et);
0599     h_met_phi->Fill(met_phi);
0600   }
0601   ////////////////////////////////////////////////////////////////////////////////
0602 }
0603
0604 // This always returns only a positive deltaPhi
0605 double EwkDQM::calcDeltaPhi(double phi1, double phi2) {
0606   double deltaPhi = phi1 - phi2;
0607
0608   if (deltaPhi < 0)
0609     deltaPhi = -deltaPhi;
0610
0611   if (deltaPhi > 3.1415926)
0612     deltaPhi = 2 * 3.1415926 - deltaPhi;
0613
0614   return deltaPhi;
0615 }
0616
0617 // Local Variables:
0618 // show-trailing-whitespace: t
0619 // truncate-lines: t
0620 // End: