analysis/plugins/PCCNTupler.cc

0001 // ----------------------------------------------------------------------
0002 // PCCNTupler
0003 // ---------
0004
0005 #include <memory>
0006 #include <string>
0007 #include <iostream>
0008 #include <fstream>
0009 #include <vector>
0010 #include <bitset>
0011
0012 #include "PCCNTupler.h"
0013
0014 #include "CondFormats/Alignment/interface/Definitions.h"
0015 #include "CondFormats/RunInfo/interface/RunInfo.h"
0016
0017 #include "DataFormats/TrackerCommon/interface/PixelBarrelName.h"
0018 #include "DataFormats/TrackerCommon/interface/PixelEndcapName.h"
0019
0020 #include "FWCore/Framework/interface/EventSetup.h"
0021 #include "FWCore/Framework/interface/MakerMacros.h"
0022 #include "FWCore/Framework/interface/ESHandle.h"
0023 #include "FWCore/ServiceRegistry/interface/Service.h"
0024
0025 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0026
0027 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
0028
0029 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0030 DEFINE_FWK_MODULE(PCCNTupler);
0031
0032 #include <TROOT.h>
0033 #include <TSystem.h>
0034 #include <TTree.h>
0035 #include <TFile.h>
0036 #include <TH1F.h>
0037 #include <TH2F.h>
0038
0039 using namespace std;
0040 using namespace edm;
0041 using namespace reco;
0042
0043 // ----------------------------------------------------------------------
0044 PCCNTupler::PCCNTupler(edm::ParameterSet const& iConfig)
0045     : fPrimaryVertexCollectionLabel(
0046           iConfig.getUntrackedParameter<InputTag>("vertexCollLabel", edm::InputTag("offlinePrimaryVertices"))),
0047       fPixelClusterLabel(
0048           iConfig.getUntrackedParameter<InputTag>("pixelClusterLabel", edm::InputTag("siPixelClusters"))),
0049       fPileUpInfoLabel(edm::InputTag("addPileupInfo")),
0050       saveType(iConfig.getUntrackedParameter<string>("saveType")),
0051       sampleType(iConfig.getUntrackedParameter<string>("sampleType")),
0052       includeVertexInformation(iConfig.getUntrackedParameter<bool>("includeVertexInformation", 1)),
0053       includePixels(iConfig.getUntrackedParameter<bool>("includePixels", 1)),
0054       includeJets(iConfig.getUntrackedParameter<bool>("includeJets", 0)),
0055       splitByBX(iConfig.getUntrackedParameter<bool>("splitByBX", 1)),
0056       pixelPhase2Geometry(iConfig.getUntrackedParameter<bool>("pixelPhase2Geometry", 0)) {
0057   cout << "----------------------------------------------------------------------" << endl;
0058   cout << "--- PCCNTupler constructor" << endl;
0059
0060   edm::Service<TFileService> fs;
0061
0062   tree = fs->make<TTree>("tree", "Pixel Cluster Counters");
0063   tree->Branch("run", &run, "run/I");
0064   tree->Branch("LS", &LS, "LS/I");
0065   tree->Branch("LN", &LN, "LN/I");
0066   tree->Branch("timeStamp_begin", &timeStamp_begin, "timeStamp_begin/i");
0067   tree->Branch("timeStamp_end", &timeStamp_end, "timeStamp_end/i");
0068   tree->Branch("eventCounter", &eventCounter, "eventCounter/I");
0069   tree->Branch("BXNo", "map<int,int>", &BXNo);
0070   if (saveType == "Event") {
0071     tree->Branch("event", &event, "event/i");
0072     tree->Branch("orbit", &orbit, "orbit/I");
0073     tree->Branch("bunchCrossing", &bunchCrossing, "bunchCrossing/I");
0074   }
0075
0076   pileup = fs->make<TH1F>("pileup", "pileup", 100, 0, 100);
0077   if (includeVertexInformation) {
0078     tree->Branch("nGoodVtx", "map<int,int>", &nGoodVtx);
0079     tree->Branch("nValidVtx", "map<int,int>", &nValidVtx);
0080     recoVtxToken = consumes<reco::VertexCollection>(fPrimaryVertexCollectionLabel);
0081     if (saveType == "Event") {
0082       tree->Branch("nVtx", &nVtx, "nVtx/I");
0083       tree->Branch("vtx_nTrk", &vtx_nTrk, "vtx_nTrk[nVtx]/I");
0084       tree->Branch("vtx_ndof", &vtx_ndof, "vtx_ndof[nVtx]/I");
0085       tree->Branch("vtx_x", &vtx_x, "vtx_x[nVtx]/F");
0086       tree->Branch("vtx_y", &vtx_y, "vtx_y[nVtx]/F");
0087       tree->Branch("vtx_z", &vtx_z, "vtx_z[nVtx]/F");
0088       tree->Branch("vtx_xError", &vtx_xError, "vtx_xError[nVtx]/F");
0089       tree->Branch("vtx_yError", &vtx_yError, "vtx_yError[nVtx]/F");
0090       tree->Branch("vtx_zError", &vtx_zError, "vtx_zError[nVtx]/F");
0091       tree->Branch("vtx_chi2", &vtx_chi2, "vtx_chi2[nVtx]/F");
0092       tree->Branch("vtx_normchi2", &vtx_normchi2, "vtx_normchi2[nVtx]/F");
0093       tree->Branch("vtx_isValid", &vtx_isValid, "vtx_isValid[nVtx]/O");
0094       tree->Branch("vtx_isFake", &vtx_isFake, "vtx_isFake[nVtx]/O");
0095       tree->Branch("vtx_isGood", &vtx_isGood, "vtx_isGood[nVtx]/O");
0096     }
0097   }
0098
0099   if (includePixels) {
0100     tree->Branch("nPixelClusters", "map<std::pair<int,int>,int>", &nPixelClusters);
0101     tree->Branch("nClusters", "map<std::pair<int,int>,int>", &nClusters);
0102     //tree->Branch("nPixelClusters","map<int,int>",&nPixelClusters);
0103     //tree->Branch("nClusters","map<int,int>",&nClusters);
0104     tree->Branch("layers", "map<int,int>", &layers);
0105     pixelToken = consumes<edmNew::DetSetVector<SiPixelCluster> >(fPixelClusterLabel);
0106   }
0107
0108   if (sampleType == "MC") {
0109     pileUpToken = consumes<std::vector<PileupSummaryInfo> >(fPileUpInfoLabel);
0110     tree->Branch("nPU", &nPU, "nPU/I");
0111   }
0112
0113   if (includeJets) {
0114     hltjetsToken_ = consumes<reco::CaloJetCollection>(edm::InputTag("ak4CaloJets"));
0115     const int kMaxJetCal = 100;
0116     jhcalpt = new float[kMaxJetCal];
0117     jhcalphi = new float[kMaxJetCal];
0118     jhcaleta = new float[kMaxJetCal];
0119     jhcale = new float[kMaxJetCal];
0120     jhcalemf = new float[kMaxJetCal];
0121     jhcaln90 = new float[kMaxJetCal];
0122     jhcaln90hits = new float[kMaxJetCal];
0123
0124     //ccla HLTJETS
0125     tree->Branch("NohJetCal", &nhjetcal, "NohJetCal/I");
0126     tree->Branch("ohJetCalPt", jhcalpt, "ohJetCalPt[NohJetCal]/F");
0127     tree->Branch("ohJetCalPhi", jhcalphi, "ohJetCalPhi[NohJetCal]/F");
0128     tree->Branch("ohJetCalEta", jhcaleta, "ohJetCalEta[NohJetCal]/F");
0129     tree->Branch("ohJetCalE", jhcale, "ohJetCalE[NohJetCal]/F");
0130     tree->Branch("ohJetCalEMF", jhcalemf, "ohJetCalEMF[NohJetCal]/F");
0131     tree->Branch("ohJetCalN90", jhcaln90, "ohJetCalN90[NohJetCal]/F");
0132     tree->Branch("ohJetCalN90hits", jhcaln90hits, "ohJetCalN90hits[NohJetCal]/F");
0133   }
0134 }
0135
0136 // ----------------------------------------------------------------------
0137 PCCNTupler::~PCCNTupler() {}
0138
0139 // ----------------------------------------------------------------------
0140 void PCCNTupler::endJob() { cout << "==>PCCNTupler> Succesfully gracefully ended job" << endl; }
0141
0142 // ----------------------------------------------------------------------
0143 void PCCNTupler::beginJob() {}
0144
0145 void PCCNTupler::beginLuminosityBlock(edm::LuminosityBlock const& lumi, edm::EventSetup const& isetup) {
0146   firstEvent = true;
0147   Reset();
0148 }
0149
0150 void PCCNTupler::endLuminosityBlock(edm::LuminosityBlock const& lumi, edm::EventSetup const& isetup) { tree->Fill(); }
0151
0152 // ----------------------------------------------------------------------
0153 void PCCNTupler::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0154   using namespace edm;
0155   using reco::VertexCollection;
0156
0157   eventCounter++;
0158
0159   saveAndReset = false;
0160   sameEvent = (event == (int)iEvent.id().event());
0161   sameLumiNib = true;  // FIXME where is this info?
0162   sameLumiSect = (LS == (int)iEvent.getLuminosityBlock().luminosityBlock());
0163
0164   // When arriving at the new LS, LN or event the tree
0165   // must be filled and branches must be reset.
0166   // The final entry is saved in the deconstructor.
0167   saveAndReset = (saveType == "LumiSect" && !sameLumiSect) || (saveType == "LumiNib" && !sameLumiNib) ||
0168                  (saveType == "Event" && !sameEvent);
0169
0170   if (!saveAndReset && !sameLumiSect && !sameLumiNib && !sameEvent) {
0171     std::cout << "Diff LS, LN and Event, but not saving/resetting..." << std::endl;
0172   }
0173
0174   if (saveAndReset) {
0175     SaveAndReset();
0176   }
0177
0178   if (sampleType == "MC") {
0179     edm::Handle<std::vector<PileupSummaryInfo> > pileUpInfo;
0180     iEvent.getByToken(pileUpToken, pileUpInfo);
0181     std::vector<PileupSummaryInfo>::const_iterator PVI;
0182     for (PVI = pileUpInfo->begin(); PVI != pileUpInfo->end(); ++PVI) {
0183       int pu_bunchcrossing = PVI->getBunchCrossing();
0184       //std::cout<<"pu_bunchcrossing getPU_NumInteractions getTrueNumInteractions "<<pu_bunchcrossing<<" "<<PVI->getPU_NumInteractions()<<" "<<PVI->getTrueNumInteractions()<<std::endl;
0185       if (pu_bunchcrossing == 0) {
0186         nPU = PVI->getPU_NumInteractions();
0187         pileup->Fill(nPU);
0188       }
0189     }
0190   }
0191
0192   // Get the Run, Lumi Section, and Event numbers, etc.
0193   run = iEvent.id().run();
0194   LS = iEvent.getLuminosityBlock().luminosityBlock();
0195   //LN    = -99; // FIXME need the luminibble
0196   event = iEvent.id().event();
0197   bunchCrossing = iEvent.bunchCrossing();
0198   if (!splitByBX) {  //if no splitting by BX then we can remove info.
0199     bunchCrossing = -10;
0200   }
0201   timeStamp_local = iEvent.time().unixTime();
0202   if (timeStamp_end < timeStamp_local)
0203     timeStamp_end = timeStamp_local;
0204   if (timeStamp_begin > timeStamp_local)
0205     timeStamp_begin = timeStamp_local;
0206   orbit = iEvent.orbitNumber();
0207   //LN    = ((int) (orbit/pow(2,12)) % 64);
0208   //int LN2    = iEvent.nibble;
0209   LN = ((int)(orbit >> 12) % 64);  // FIXME need the luminibble
0210
0211   bxModKey.first = bunchCrossing;
0212   bxModKey.second = -1;
0213
0214   if ((BXNo.count(bunchCrossing) == 0 || nGoodVtx.count(bunchCrossing) == 0) &&
0215       !(BXNo.count(bunchCrossing) == 0 && nGoodVtx.count(bunchCrossing) == 0)) {
0216     std::cout << "BXNo and nGoodVtx should have the same keys but DO NOT!!!" << std::endl;
0217   }
0218
0219   if (BXNo.count(bunchCrossing) == 0) {
0220     BXNo[bunchCrossing] = 0;
0221   }
0222
0223   if (nGoodVtx.count(bunchCrossing) == 0) {
0224     nGoodVtx[bunchCrossing] = 0;
0225     nValidVtx[bunchCrossing] = 0;
0226   }
0227
0228   BXNo[bunchCrossing] = BXNo[bunchCrossing] + 1;
0229   // add the vertex information
0230
0231   if (includeVertexInformation) {
0232     edm::Handle<reco::VertexCollection> recVtxs;
0233     iEvent.getByToken(recoVtxToken, recVtxs);
0234
0235     if (recVtxs.isValid()) {
0236       //nVtx=recVtxs->size();
0237       int ivtx = 0;
0238       for (reco::VertexCollection::const_iterator v = recVtxs->begin(); v != recVtxs->end(); ++v) {
0239         if (v->isFake())
0240           continue;
0241         vtx_isGood[ivtx] = false;
0242         vtx_nTrk[ivtx] = v->tracksSize();
0243         vtx_ndof[ivtx] = (int)v->ndof();
0244         vtx_x[ivtx] = v->x();
0245         vtx_y[ivtx] = v->y();
0246         vtx_z[ivtx] = v->z();
0247         vtx_xError[ivtx] = v->xError();
0248         vtx_yError[ivtx] = v->yError();
0249         vtx_zError[ivtx] = v->zError();
0250         vtx_chi2[ivtx] = v->chi2();
0251         vtx_normchi2[ivtx] = v->normalizedChi2();
0252         vtx_isValid[ivtx] = v->isValid();
0253         vtx_isFake[ivtx] = v->isFake();
0254         if (vtx_isValid[ivtx] && (vtx_isFake[ivtx] == 0)) {
0255           nValidVtx[bunchCrossing] = nValidVtx[bunchCrossing] + 1;
0256         }
0257         if (vtx_ndof[ivtx] > 4 && vtx_isValid[ivtx] && (vtx_isFake[ivtx] == 0)) {
0258           if (vtx_nTrk[ivtx] > 0) {
0259             nGoodVtx[bunchCrossing] = nGoodVtx[bunchCrossing] + 1;
0260             vtx_isGood[ivtx] = true;
0261           }
0262         }
0263         ivtx++;
0264       }
0265       nVtx = ivtx;
0266     }
0267   }
0268
0269   if (includeJets) {
0270     edm::Handle<reco::CaloJetCollection> hltjets;
0271     iEvent.getByToken(hltjetsToken_, hltjets);
0272     bool valid = hltjets.isValid();
0273     if (not valid) {
0274       std::cout << "hltjets not valid " << std::endl;
0275       nhjetcal = -1;
0276     } else {
0277       reco::CaloJetCollection mycalojets;
0278       mycalojets = *hltjets;
0279       //std::sort(mycalojets.begin(),mycalojets.end(),PtGreater());
0280       typedef reco::CaloJetCollection::const_iterator cjiter;
0281       int jhcal = 0;
0282       for (cjiter i = mycalojets.begin(); i != mycalojets.end(); i++) {
0283         if (i->pt() > 5 && i->energy() > 0.) {
0284           jhcalpt[jhcal] = i->pt();
0285           jhcalphi[jhcal] = i->phi();
0286           jhcaleta[jhcal] = i->eta();
0287           jhcale[jhcal] = i->energy();
0288           jhcalemf[jhcal] = i->emEnergyFraction();
0289           jhcaln90[jhcal] = i->n90();
0290           //jetID->calculate( iEvent, *i );
0291           //jhcaln90hits[jhcal] = jetID->n90Hits();
0292           jhcal++;
0293         }
0294       }
0295       nhjetcal = jhcal;
0296     }
0297   }
0298
0299   int NumPixelBarrelLayers = 3;
0300   if (pixelPhase2Geometry) {
0301     NumPixelBarrelLayers = 4;
0302   }
0303   // -- Pixel cluster
0304   if (includePixels) {
0305     edm::Handle<edmNew::DetSetVector<SiPixelCluster> > hClusterColl;
0306     iEvent.getByToken(pixelToken, hClusterColl);
0307     if (!hClusterColl.failedToGet()) {
0308       const edmNew::DetSetVector<SiPixelCluster>& clustColl = *hClusterColl;
0309       // ----------------------------------------------------------------------
0310       // -- Clusters without tracks
0311
0312       for (edmNew::DetSetVector<SiPixelCluster>::const_iterator isearch = clustColl.begin(); isearch != clustColl.end();
0313            ++isearch) {
0314         // these are sorted by modules so we pick the current one
0315         edmNew::DetSet<SiPixelCluster> mod = *isearch;
0316         if (mod.empty()) {
0317           continue;
0318         }  // skip empty modules
0319         DetId detId = mod.id();
0320
0321         bxModKey.second = detId();
0322         for (edmNew::DetSet<SiPixelCluster>::const_iterator di = mod.begin(); di != mod.end(); ++di) {
0323           if (nPixelClusters.count(bxModKey) == 0) {
0324             nPixelClusters[bxModKey] = 0;
0325           }
0326           nPixelClusters[bxModKey] = nPixelClusters[bxModKey] + 1;
0327
0328           int nCluster = isearch->size();
0329           if (nClusters.count(bxModKey) == 0) {
0330             nClusters[bxModKey] = 0;
0331           }
0332           nClusters[bxModKey] += nCluster;
0333
0334           if (detId.subdetId() == PixelSubdetector::PixelBarrel) {
0335             PixelBarrelName detName = PixelBarrelName(detId);
0336             int layer = detName.layerName();
0337             if (layers.count(detId()) == 0) {
0338               layers[detId()] = layer;
0339             }
0340           } else {
0341             assert(detId.subdetId() == PixelSubdetector::PixelEndcap);
0342             PixelEndcapName detName = PixelEndcapName(detId);
0343             int disk = detName.diskName();
0344             if (layers.count(detId()) == 0) {
0345               layers[detId()] = disk + NumPixelBarrelLayers;
0346             }
0347           }
0348           //}
0349         }
0350       }
0351     }
0352   }
0353 }
0354
0355 void PCCNTupler::Reset() {
0356   nVtx = 0;
0357   nPixelClusters.clear();
0358   nClusters.clear();
0359   layers.clear();
0360   BXNo.clear();
0361   nValidVtx.clear();
0362   nGoodVtx.clear();
0363   eventCounter = 1;
0364   timeStamp_end = 0;
0365   timeStamp_begin = -1;
0366 }
0367
0368 void PCCNTupler::SaveAndReset() {
0369   if (!firstEvent)
0370     tree->Fill();
0371   Reset();
0372   firstEvent = false;
0373 }