EcalTPGAnalyzer.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Class:      EcalTPGAnalyzer
//
//
// Original Author:  Pascal Paganini
//
//

// system include files
#include <memory>
#include <utility>

// user include files
#include "FWCore/Framework/interface/EDAnalyzer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"

#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
#include "DataFormats/EcalDigi/interface/EcalTriggerPrimitiveDigi.h"

#include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"

#include "CondFormats/L1TObjects/interface/L1GtTriggerMenu.h"
#include "CondFormats/DataRecord/interface/L1GtTriggerMenuRcd.h"
#include "CondFormats/L1TObjects/interface/L1GtTriggerMask.h"
#include "CondFormats/DataRecord/interface/L1GtTriggerMaskAlgoTrigRcd.h"
#include "CondFormats/DataRecord/interface/L1GtTriggerMaskTechTrigRcd.h"
#include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetupFwd.h"
#include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetup.h"
#include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
#include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerRecord.h"

#include "EcalTPGAnalyzer.h"

#include <TMath.h>
#include <sstream>

using namespace edm;
class CaloSubdetectorGeometry;

EcalTPGAnalyzer::EcalTPGAnalyzer(const edm::ParameterSet& iConfig) {
  tpCollection_ = iConfig.getParameter<edm::InputTag>("TPCollection");
  tpEmulatorCollection_ = iConfig.getParameter<edm::InputTag>("TPEmulatorCollection");
  digiCollectionEB_ = iConfig.getParameter<edm::InputTag>("DigiCollectionEB");
  digiCollectionEE_ = iConfig.getParameter<edm::InputTag>("DigiCollectionEE");
  gtRecordCollectionTag_ = iConfig.getParameter<std::string>("GTRecordCollection");

  allowTP_ = iConfig.getParameter<bool>("ReadTriggerPrimitives");
  useEE_ = iConfig.getParameter<bool>("UseEndCap");
  print_ = iConfig.getParameter<bool>("Print");

  // file
  file_ = new TFile("ECALTPGtree.root", "RECREATE");
  file_->cd();

  // tree
  tree_ = new TTree("EcalTPGAnalysis", "EcalTPGAnalysis");

  tree_->Branch("runNb", &treeVariables_.runNb, "runNb/i");                                                //
  tree_->Branch("evtNb", &treeVariables_.evtNb, "evtNb/i");                                                //
  tree_->Branch("bxNb", &treeVariables_.bxNb, "bxNb/i");                                                   //
  tree_->Branch("orbitNb", &treeVariables_.orbitNb, "orbitNb/i");                                          //
  tree_->Branch("nbOfActiveTriggers", &treeVariables_.nbOfActiveTriggers, "nbOfActiveTriggers/i");         //
  tree_->Branch("activeTriggers", treeVariables_.activeTriggers, "activeTriggers[nbOfActiveTriggers]/I");  //

  tree_->Branch("nbOfTowers", &treeVariables_.nbOfTowers, "nbOfTowers/i");             //
  tree_->Branch("ieta", treeVariables_.ieta, "ieta[nbOfTowers]/I");                    //
  tree_->Branch("iphi", treeVariables_.iphi, "iphi[nbOfTowers]/I");                    //
  tree_->Branch("nbOfXtals", treeVariables_.nbOfXtals, "nbOfXtals[nbOfTowers]/I");     //
  tree_->Branch("rawTPData", treeVariables_.rawTPData, "rawTPData[nbOfTowers]/I");     //
  tree_->Branch("rawTPEmul1", treeVariables_.rawTPEmul1, "rawTPEmul1[nbOfTowers]/I");  //
  tree_->Branch("rawTPEmul2", treeVariables_.rawTPEmul2, "rawTPEmul2[nbOfTowers]/I");  //
  tree_->Branch("rawTPEmul3", treeVariables_.rawTPEmul3, "rawTPEmul3[nbOfTowers]/I");  //
  tree_->Branch("rawTPEmul4", treeVariables_.rawTPEmul4, "rawTPEmul4[nbOfTowers]/I");  //
  tree_->Branch("rawTPEmul5", treeVariables_.rawTPEmul5, "rawTPEmul5[nbOfTowers]/I");  //
  tree_->Branch("eRec", treeVariables_.eRec, "eRec[nbOfTowers]/F");                    //
}

EcalTPGAnalyzer::~EcalTPGAnalyzer() {
  file_->cd();
  tree_->Write();
  file_->Close();
}

void EcalTPGAnalyzer::beginRun(edm::Run const&, edm::EventSetup const& evtSetup) {
  // geometry
  ESHandle<CaloGeometry> theGeometry;
  ESHandle<CaloSubdetectorGeometry> theEndcapGeometry_handle, theBarrelGeometry_handle;

  evtSetup.get<CaloGeometryRecord>().get(theGeometry);
  evtSetup.get<EcalEndcapGeometryRecord>().get("EcalEndcap", theEndcapGeometry_handle);
  evtSetup.get<EcalBarrelGeometryRecord>().get("EcalBarrel", theBarrelGeometry_handle);

  evtSetup.get<IdealGeometryRecord>().get(eTTmap_);
  theEndcapGeometry_ = &(*theEndcapGeometry_handle);
  theBarrelGeometry_ = &(*theBarrelGeometry_handle);
}

void EcalTPGAnalyzer::analyze(edm::Event const& iEvent, edm::EventSetup const& iSetup) {
  using namespace edm;
  using namespace std;

  if (print_)
    std::cout << "===========" << iEvent.id() << std::endl;

  map<EcalTrigTowerDetId, towerEner> mapTower;
  map<EcalTrigTowerDetId, towerEner>::iterator itTT;

  // get Evts info

  treeVariables_.runNb = iEvent.id().run();
  treeVariables_.evtNb = iEvent.id().event();
  treeVariables_.bxNb = iEvent.bunchCrossing();
  treeVariables_.orbitNb = iEvent.orbitNumber();

  // get L1 info

  edm::Handle<L1GlobalTriggerReadoutRecord> gtRecord;
  iEvent.getByLabel(edm::InputTag(gtRecordCollectionTag_), gtRecord);
  DecisionWord dWord = gtRecord->decisionWord();  // this will get the decision word *before* masking disabled bits

  edm::ESHandle<L1GtTriggerMask> l1GtTmAlgo;
  iSetup.get<L1GtTriggerMaskAlgoTrigRcd>().get(l1GtTmAlgo);
  std::vector<unsigned int> triggerMaskAlgoTrig = l1GtTmAlgo.product()->gtTriggerMask();

  // apply masks on algo
  int iDaq = 0;
  int iBit = -1;
  treeVariables_.nbOfActiveTriggers = 0;
  for (std::vector<bool>::iterator itBit = dWord.begin(); itBit != dWord.end(); ++itBit) {
    iBit++;
    int maskBit = triggerMaskAlgoTrig[iBit] & (1 << iDaq);
    if (maskBit)
      *itBit = false;
    if (*itBit) {
      treeVariables_.activeTriggers[treeVariables_.nbOfActiveTriggers] = iBit;
      treeVariables_.nbOfActiveTriggers++;
    }
  }

  // Get TP data

  edm::Handle<EcalTrigPrimDigiCollection> tp;
  iEvent.getByLabel(tpCollection_, tp);
  if (print_)
    std::cout << "TP collection size=" << tp.product()->size() << std::endl;

  for (unsigned int i = 0; i < tp.product()->size(); i++) {
    EcalTriggerPrimitiveDigi d = (*(tp.product()))[i];
    const EcalTrigTowerDetId TPtowid = d.id();
    towerEner tE;
    tE.iphi_ = TPtowid.iphi();
    tE.ieta_ = TPtowid.ieta();
    tE.tpgADC_ = d[0].raw();
    mapTower[TPtowid] = tE;
  }

  // Get Emulators TP

  edm::Handle<EcalTrigPrimDigiCollection> tpEmul;
  iEvent.getByLabel(tpEmulatorCollection_, tpEmul);
  if (print_)
    std::cout << "TPEmulator collection size=" << tpEmul.product()->size() << std::endl;

  for (unsigned int i = 0; i < tpEmul.product()->size(); i++) {
    EcalTriggerPrimitiveDigi d = (*(tpEmul.product()))[i];
    const EcalTrigTowerDetId TPtowid = d.id();
    itTT = mapTower.find(TPtowid);
    if (itTT != mapTower.end())
      for (int j = 0; j < 5; j++)
        (itTT->second).tpgEmul_[j] = d[j].raw();
  }

  // Get nb of crystals read out

  // Get EB xtal digi inputs
  edm::Handle<EBDigiCollection> digiEB;
  iEvent.getByLabel(digiCollectionEB_, digiEB);

  for (unsigned int i = 0; i < digiEB.product()->size(); i++) {
    const EBDataFrame& df = (*(digiEB.product()))[i];
    const EBDetId& id = df.id();
    const EcalTrigTowerDetId towid = id.tower();
    itTT = mapTower.find(towid);
    if (itTT != mapTower.end())
      (itTT->second).nbXtal_++;
  }

  if (useEE_) {
    // Get EE xtal digi inputs
    edm::Handle<EEDigiCollection> digiEE;
    iEvent.getByLabel(digiCollectionEE_, digiEE);
    for (unsigned int i = 0; i < digiEE.product()->size(); i++) {
      const EEDataFrame& df = (*(digiEE.product()))[i];
      const EEDetId& id = df.id();
      const EcalTrigTowerDetId towid = (*eTTmap_).towerOf(id);
      itTT = mapTower.find(towid);
      if (itTT != mapTower.end())
        (itTT->second).nbXtal_++;
    }
  }

  // Get rechits

  //... to be completed

  // fill tree

  treeVariables_.nbOfTowers = mapTower.size();
  int towerNb = 0;
  for (itTT = mapTower.begin(); itTT != mapTower.end(); ++itTT) {
    treeVariables_.ieta[towerNb] = (itTT->second).ieta_;
    treeVariables_.iphi[towerNb] = (itTT->second).iphi_;
    treeVariables_.nbOfXtals[towerNb] = (itTT->second).nbXtal_;
    treeVariables_.rawTPData[towerNb] = (itTT->second).tpgADC_;
    treeVariables_.rawTPEmul1[towerNb] = (itTT->second).tpgEmul_[0];
    treeVariables_.rawTPEmul2[towerNb] = (itTT->second).tpgEmul_[1];
    treeVariables_.rawTPEmul3[towerNb] = (itTT->second).tpgEmul_[2];
    treeVariables_.rawTPEmul4[towerNb] = (itTT->second).tpgEmul_[3];
    treeVariables_.rawTPEmul5[towerNb] = (itTT->second).tpgEmul_[4];
    treeVariables_.eRec[towerNb] = (itTT->second).eRec_;
    towerNb++;
  }

  tree_->Fill();
}