Configuration/macros/fwliteExample.C

0001 // Usage:
0002 // .x rootlogon.C
0003 // .x fwliteExample.C++
0004
0005 #include <memory>
0006 #include <string>
0007 #include <vector>
0008 #include <iostream>
0009 #include <cmath>
0010
0011 #include <TH1D.h>
0012 #include <TH2D.h>
0013 #include <TNtuple.h>
0014 #include <TFile.h>
0015 #include <TROOT.h>
0016 #include <TSystem.h>
0017 #include <TString.h>
0018
0019 #if !defined(__CINT__) && !defined(__MAKECINT__)
0020 #include "DataFormats/FWLite/interface/Handle.h"
0021 #include "DataFormats/FWLite/interface/Event.h"
0022 #include "DataFormats/FWLite/interface/ChainEvent.h"
0023 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
0024 #include "DataFormats/CaloTowers/interface/CaloTower.h"
0025 #include "DataFormats/CaloTowers/interface/CaloTowerDefs.h"
0026 #include "DataFormats/JetReco/interface/CaloJet.h"
0027 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
0028 #include "DataFormats/TrackReco/interface/Track.h"
0029 #include "DataFormats/VertexReco/interface/Vertex.h"
0030 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0031 #include "FWCore/Common/interface/TriggerNames.h"
0032 #include "DataFormats/Common/interface/TriggerResults.h"
0033 #include "DataFormats/Math/interface/deltaPhi.h"
0034 #endif
0035
0036 void fwliteExample(bool debug=false){
0037
0038   // event cuts
0039   const unsigned int maxEvents = -1;
0040   const double hfEThreshold = 3.0;
0041   const int nTowerThreshold = 1;
0042
0043   // track cuts
0044   const double normD0Cut = 5.0;
0045   const double normDZCut = 5.0;
0046   const double ptDebug = 3.0;
0047
0048   // trigger names
0049   const int nTrigs = 4;
0050   const char *hltNames[nTrigs] = {"HLT_MinBiasBSC","HLT_L1Jet6U","HLT_Jet15U","HLT_Jet30U"};
0051
0052   //----- input files (900 GeV data) -----
0053   vector<string> fileNames;
0054   fileNames.push_back("./hiCommonSkimAOD.root");
0055   //fileNames.push_back("../test/hiCommonSkimAOD.root");
0056   fwlite::ChainEvent event(fileNames);
0057
0058   //----- define output hists/trees in directories of output file -----
0059   TFile *outFile = new TFile("output_fwlite.root","recreate");
0060   TH1D::SetDefaultSumw2();
0061
0062   // evt hists
0063   outFile->cd(); outFile->mkdir("evt"); outFile->cd("evt");
0064   TH1D *hL1TechBits = new TH1D("hL1TechBits","L1 technical trigger bits before mask",64,-0.5,63.5);
0065   TH2D *hHfTowers   = new TH2D("hHfTowers","Number of HF tower above threshold; positive side; negative side",80,-0.5,79.5,80,-0.5,79.5);
0066   TH1D *hHLTPaths   = new TH1D("hHLTPaths","HLT Paths",3,0,3);
0067   hHLTPaths->SetCanExtend(TH1::kAllAxes);
0068
0069   // vtx hists
0070   outFile->cd(); outFile->mkdir("vtx"); outFile->cd("vtx");
0071   TH1D *hVtxTrks    = new TH1D("hVtxTrks","number of tracks used to fit pixel vertex",50,-0.5,49.5);
0072   TH1D *hVtxZ       = new TH1D("hVtxZ","z position of best reconstructed pixel vertex", 80,-20,20);
0073
0074   // track hists
0075   outFile->cd(); outFile->mkdir("trk"); outFile->cd("trk");
0076   TH1D *hTrkPt      = new TH1D("hTrkPt","track p_{T}; p_{T} [GeV/c]", 80, 0.0, 20.0);
0077   TH1D *hTrkEta     = new TH1D("hTrkEta","track #eta; #eta", 60, -3.0, 3.0);
0078   TH1D *hTrkPhi     = new TH1D("hTrkPhi","track #phi; #phi [radians]", 56, -3.5, 3.5);
0079
0080   // correlation hists
0081   outFile->cd(); outFile->mkdir("corr"); outFile->cd("corr");
0082   TH2D *hDetaDphi   = new TH2D("hDetaDphi","raw two-particle correlation; #Delta#eta; #Delta#phi",50,-5.0,5.0,50,-3.1416,3.1416);
0083
0084   // debug ntuple
0085   outFile->cd();
0086   TNtuple *nt=0;
0087   if(debug) nt = new TNtuple("nt","track debug ntuple",
0088                  "pt:eta:phi:hits:pterr:d0:d0err:dz:dzerr:jet6:jet15:jet30");
0089
0090   //----- loop over events -----
0091   unsigned int iEvent=0;
0092   for(event.toBegin(); !event.atEnd(); ++event, ++iEvent){
0093
0094     if( iEvent == maxEvents ) break;
0095     if( iEvent % 1000 == 0 ) cout << "Processing " << iEvent<< "th event: "
0096                   << "run " << event.id().run()
0097                   << ", lumi " << event.luminosityBlock()
0098                   << ", evt " << event.id().event() << endl;
0099
0100     // select on L1 trigger bits
0101     fwlite::Handle<L1GlobalTriggerReadoutRecord> gt;
0102     gt.getByLabel(event, "gtDigis");
0103     const TechnicalTriggerWord&  word = gt->technicalTriggerWord(); //before mask
0104     for(int bit=0; bit<64; bit++) hL1TechBits->Fill(bit,word.at(bit));
0105     if(!word.at(0)) continue;  // BPTX coincidence
0106     if(!(word.at(40) || word.at(41))) continue; // BSC coincidence
0107     if(word.at(36) || word.at(37) || word.at(38) || word.at(39)) continue; // BSC halo
0108
0109     // select on coincidence of HF towers above threshold
0110     fwlite::Handle<CaloTowerCollection> towers;
0111     towers.getByLabel(event, "towerMaker");
0112     int nHfTowersN=0, nHfTowersP=0;
0113     for(CaloTowerCollection::const_iterator calo = towers->begin(); calo != towers->end(); ++calo) {
0114       if(calo->energy() < hfEThreshold) continue;
0115       if(calo->eta()>3) nHfTowersP++;
0116       if(calo->eta()<-3) nHfTowersN++;
0117     }
0118     hHfTowers->Fill(nHfTowersP,nHfTowersN);
0119     if(nHfTowersP < nTowerThreshold || nHfTowersN < nTowerThreshold) continue;
0120
0121
0122     // get hlt bits
0123     bool accept[nTrigs]={};
0124     fwlite::Handle<edm::TriggerResults> triggerResults;
0125     triggerResults.getByLabel(event, "TriggerResults","","HLT");
0126     const edm::TriggerNames triggerNames = event.triggerNames(*triggerResults);
0127     for(int i=0; i<nTrigs; i++) {
0128       accept[i] = triggerResults->accept(triggerNames.triggerIndex(hltNames[i]));
0129       if(accept[i]) hHLTPaths->Fill(hltNames[i],1);
0130     }
0131
0132     // select on requirement of valid vertex
0133     math::XYZPoint vtxpoint(0,0,0);
0134     fwlite::Handle<std::vector<reco::Vertex> > vertices;
0135     vertices.getByLabel(event, "hiSelectedVertex");
0136     if(!vertices->size()) continue;
0137     const reco::Vertex & vtx = (*vertices)[0];
0138     vtxpoint.SetCoordinates(vtx.x(),vtx.y(),vtx.z());
0139     hVtxTrks->Fill(vtx.tracksSize());
0140     hVtxZ->Fill(vtx.z());
0141
0142     // get beamspot
0143     fwlite::Handle<reco::BeamSpot> beamspot;
0144     beamspot.getByLabel(event, "offlineBeamSpot");
0145
0146     //----- loop over tracks -----
0147
0148     fwlite::Handle<std::vector<reco::Track> > tracks;
0149     tracks.getByLabel(event, "hiSelectedTracks");
0150
0151     for(unsigned it=0; it<tracks->size(); ++it){
0152
0153       const reco::Track & trk = (*tracks)[it];
0154
0155       // select tracks based on transverse proximity to beamspot
0156       double dxybeam = trk.dxy(beamspot->position());
0157       if(fabs(dxybeam/trk.d0Error()) > normD0Cut) continue;
0158
0159       // select tracks based on z-proximity to best vertex
0160       double dzvtx = trk.dz(vtxpoint);
0161       if(fabs(dzvtx/trk.dzError()) > normDZCut) continue;
0162
0163       // fill selected track histograms and debug ntuple
0164       hTrkPt->Fill(trk.pt());
0165       hTrkEta->Fill(trk.eta());
0166       hTrkPhi->Fill(trk.phi());
0167       if(debug && trk.pt() > ptDebug) // fill debug ntuple for selection of tracks
0168     nt->Fill(trk.pt(),trk.eta(),trk.phi(),trk.numberOfValidHits(),trk.ptError(),
0169          dxybeam,trk.d0Error(),dzvtx,trk.dzError(),accept[1],accept[2],accept[3]);
0170
0171     }
0172
0173     //----- loop over jets -----
0174
0175     fwlite::Handle<vector<reco::CaloJet> > jets;
0176     jets.getByLabel(event, "iterativeConePu5CaloJets");
0177
0178     //----- loop over muons -----
0179
0180     //----- loop over photons -----
0181
0182     //----- loop over charged candidates -----
0183
0184     fwlite::Handle<vector<reco::RecoChargedCandidate> > candidates;
0185     candidates.getByLabel(event, "allTracks");
0186
0187     for(unsigned it1=0; it1<candidates->size(); ++it1) {
0188       const reco::RecoChargedCandidate & c1 = (*candidates)[it1];
0189       for(unsigned it2=0; it2<candidates->size(); ++it2) {
0190     if(it1==it2) continue;
0191     const reco::RecoChargedCandidate & c2 = (*candidates)[it2];
0192     hDetaDphi->Fill(c1.eta()-c2.eta(),deltaPhi(c1.phi(),c2.phi()));
0193       }
0194     }
0195
0196   } //end loop over events
0197
0198   cout << "Number of events processed : " << iEvent << endl;
0199   cout << "Number passing all event selection cuts: " << hVtxZ->GetEntries() << endl;
0200
0201   // write to output file
0202   hHLTPaths->LabelsDeflate();
0203   outFile->Write();
0204   outFile->ls();
0205   outFile->Close();
0206
0207 }