DTTrigTest.cc

Go to the documentation of this file.

//-------------------------------------------------
//
//
//--------------------------------------------------

// This class's header
#include "L1Trigger/DTTrigger/interface/DTTrigTest.h"

// Framework related classes
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ESHandle.h"

#include "L1TriggerConfig/DTTPGConfig/interface/DTConfigManager.h"
#include "L1TriggerConfig/DTTPGConfig/interface/DTConfigManagerRcd.h"

// Trigger and DataFormats headers
#include "L1Trigger/DTSectorCollector/interface/DTSectCollPhSegm.h"
#include "L1Trigger/DTSectorCollector/interface/DTSectCollThSegm.h"
#include "SimDataFormats/Track/interface/SimTrackContainer.h"
#include "SimDataFormats/Vertex/interface/SimVertexContainer.h"

// ROOT headers
#include "TROOT.h"
#include "TTree.h"
#include "TFile.h"

// Collaborating classes
#include "DataFormats/Math/interface/LorentzVector.h"
#include <CLHEP/Vector/LorentzVector.h>

// C++ headers
#include <iostream>
#include <cmath>
#include <ctime>

using namespace std;
using namespace edm;

const double DTTrigTest::my_TtoTDC = 32. / 25.;

DTTrigTest::DTTrigTest(const ParameterSet& pset) : my_trig(nullptr) {
  my_debug = pset.getUntrackedParameter<bool>("debug");
  string outputfile = pset.getUntrackedParameter<string>("outputFileName");
  if (my_debug)
    cout << "[DTTrigTest] Creating rootfile " << outputfile << endl;
  my_rootfile = new TFile(outputfile.c_str(), "RECREATE");
  my_tree = new TTree("h1", "GMT", 0);
  my_params = pset;
  if (my_debug)
    cout << "[DTTrigTest] Constructor executed!!!" << endl;
  my_trig = new DTTrig(my_params, consumesCollector());
}

DTTrigTest::~DTTrigTest() {
  if (my_trig != nullptr)
    delete my_trig;
  delete my_rootfile;
  if (my_debug)
    cout << "[DTTrigTest] Destructor executed!!!" << endl;
}

void DTTrigTest::endJob() {
  if (my_debug)
    cout << "[DTTrigTest] Writing Tree and Closing File" << endl;
  my_tree->Write();
  delete my_tree;
  my_rootfile->Close();
}

//void DTTrigTest::beginJob(const EventSetup & iEventSetup){
void DTTrigTest::beginJob() {
  // get DTConfigManager
  // ESHandle< DTConfigManager > confManager ;
  // iEventSetup.get< DTConfigManagerRcd >().get( confManager ) ;

  //for testing purpose....
  //DTBtiId btiid(1,1,1,1,1);
  //confManager->getDTConfigBti(btiid)->print();

  //   my_trig = new DTTrig(my_params);

  //   my_trig->createTUs(iEventSetup);
  //   if (my_debug)
  //     cout << "[DTTrigTest] TU's Created" << endl;

  // BOOKING of the tree's varables
  // GENERAL block branches
  my_tree->Branch("Run", &runn, "Run/I");
  my_tree->Branch("Event", &eventn, "Event/I");
  my_tree->Branch("Weight", &weight, "Weight/F");
  // GEANT block branches
  my_tree->Branch("Ngen", &ngen, "Ngen/I");
  my_tree->Branch("Pxgen", pxgen, "Pxgen[Ngen]/F");
  my_tree->Branch("Pygen", pygen, "Pygen[Ngen]/F");
  my_tree->Branch("Pzgen", pzgen, "Pzgen[Ngen]/F");
  my_tree->Branch("Ptgen", ptgen, "Ptgen[Ngen]/F");
  my_tree->Branch("Etagen", etagen, "Etagen[Ngen]/F");
  my_tree->Branch("Phigen", phigen, "Phigen[Ngen]/F");
  my_tree->Branch("Chagen", chagen, "Chagen[Ngen]/I");
  my_tree->Branch("Vxgen", vxgen, "Vxgen[Ngen]/F");
  my_tree->Branch("Vygen", vygen, "Vygen[Ngen]/F");
  my_tree->Branch("Vzgen", vzgen, "Vzgen[Ngen]/F");
  // L1MuDTBtiChipS block
  my_tree->Branch("Nbti", &nbti, "Nbti/I");
  my_tree->Branch("bwh", bwh, "bwh[Nbti]/I");
  my_tree->Branch("bstat", bstat, "bstat[Nbti]/I");
  my_tree->Branch("bsect", bsect, "bsect[Nbti]/I");
  my_tree->Branch("bsl", bsl, "bsl[Nbti]/I");
  my_tree->Branch("bnum", bnum, "bnum[Nbti]/I");
  my_tree->Branch("bbx", bbx, "bbx[Nbti]/I");
  my_tree->Branch("bcod", bcod, "bcod[Nbti]/I");
  my_tree->Branch("bk", bk, "bk[Nbti]/I");
  my_tree->Branch("bx", bx, "bx[Nbti]/I");
  my_tree->Branch("bposx", bposx, "bposx[Nbti]/F");
  my_tree->Branch("bposy", bposy, "bposy[Nbti]/F");
  my_tree->Branch("bposz", bposz, "bposz[Nbti]/F");
  my_tree->Branch("bdirx", bdirx, "bdirx[Nbti]/F");
  my_tree->Branch("bdiry", bdiry, "bdiry[Nbti]/F");
  my_tree->Branch("bdirz", bdirz, "bdirz[Nbti]/F");
  // L1MuDTTracoChipS block
  my_tree->Branch("Ntraco", &ntraco, "Ntraco/I");
  my_tree->Branch("twh", twh, "twh[Ntraco]/I");
  my_tree->Branch("tstat", tstat, "tstat[Ntraco]/I");
  my_tree->Branch("tsect", tsect, "tsect[Ntraco]/I");
  my_tree->Branch("tnum", tnum, "tnum[Ntraco]/I");
  my_tree->Branch("tbx", tbx, "tbx[Ntraco]/I");
  my_tree->Branch("tcod", tcod, "tcod[Ntraco]/I");
  my_tree->Branch("tk", tk, "tk[Ntraco]/I");
  my_tree->Branch("tx", tx, "tx[Ntraco]/I");
  my_tree->Branch("tposx", tposx, "tposx[Ntraco]/F");
  my_tree->Branch("tposy", tposy, "tposy[Ntraco]/F");
  my_tree->Branch("tposz", tposz, "tposz[Ntraco]/F");
  my_tree->Branch("tdirx", tdirx, "tdirx[Ntraco]/F");
  my_tree->Branch("tdiry", tdiry, "tdiry[Ntraco]/F");
  my_tree->Branch("tdirz", tdirz, "tdirz[Ntraco]/F");
  // TSPHI block
  my_tree->Branch("Ntsphi", &ntsphi, "Ntsphi/I");
  my_tree->Branch("swh", swh, "swh[Ntsphi]/I");
  my_tree->Branch("sstat", sstat, "sstat[Ntsphi]/I");
  my_tree->Branch("ssect", ssect, "ssect[Ntsphi]/I");
  my_tree->Branch("sbx", sbx, "sbx[Ntsphi]/I");
  my_tree->Branch("scod", scod, "scod[Ntsphi]/I");
  my_tree->Branch("sphi", sphi, "sphi[Ntsphi]/I");
  my_tree->Branch("sphib", sphib, "sphib[Ntsphi]/I");
  my_tree->Branch("sposx", sposx, "sposx[Ntsphi]/F");
  my_tree->Branch("sposy", sposy, "sposy[Ntsphi]/F");
  my_tree->Branch("sposz", sposz, "sposz[Ntsphi]/F");
  my_tree->Branch("sdirx", sdirx, "sdirx[Ntsphi]/F");
  my_tree->Branch("sdiry", sdiry, "sdiry[Ntsphi]/F");
  my_tree->Branch("sdirz", sdirz, "sdirz[Ntsphi]/F");
  // TSTHETA block
  my_tree->Branch("Ntstheta", &ntstheta, "Ntstheta/I");
  my_tree->Branch("thwh", thwh, "thwh[Ntstheta]/I");
  my_tree->Branch("thstat", thstat, "thstat[Ntstheta]/I");
  my_tree->Branch("thsect", thsect, "thsect[Ntstheta]/I");
  my_tree->Branch("thbx", thbx, "thbx[Ntstheta]/I");
  my_tree->Branch("thcode", thcode, "thcode[Ntstheta][7]/I");
  my_tree->Branch("thpos", thpos, "thpos[Ntstheta][7]/I");
  my_tree->Branch("thqual", thqual, "thqual[Ntstheta][7]/I");
  // SC PHI block
  my_tree->Branch("Nscphi", &nscphi, "Nscphi/I");
  my_tree->Branch("scphwh", scphwh, "scphwh[Nscphi]/I");
  my_tree->Branch("scphstat", scphstat, "scphstat[Nscphi]/I");
  my_tree->Branch("scphsect", scphsect, "scphsect[Nscphi]/I");
  my_tree->Branch("scphbx", scphbx, "scphbx[Nscphi]/I");
  my_tree->Branch("scphcod", scphcod, "scphcod[Nscphi]/I");
  my_tree->Branch("scphphi", scphphi, "scphphi[Nscphi]/I");
  my_tree->Branch("scphphib", scphphib, "scphphib[Nscphi]/I");
  my_tree->Branch("scphposx", scphposx, "scphposx[Nscphi]/F");
  my_tree->Branch("scphposy", scphposy, "scphposy[Nscphi]/F");
  my_tree->Branch("scphposz", scphposz, "scphposz[Nscphi]/F");
  my_tree->Branch("scphdirx", scphdirx, "scphdirx[Nscphi]/F");
  my_tree->Branch("scphdiry", scphdiry, "scphdiry[Nscphi]/F");
  my_tree->Branch("scphdirz", scphdirz, "scphdirz[Nscphi]/F");
  // SC THETA block
  my_tree->Branch("Nsctheta", &nsctheta, "Nsctheta/I");
  my_tree->Branch("scthwh", scthwh, "scthwh[Nsctheta]/I");
  my_tree->Branch("scthstat", scthstat, "scthstat[Nsctheta]/I");
  my_tree->Branch("scthsect", scthsect, "scthsect[Nsctheta]/I");
  my_tree->Branch("scthbx", scthbx, "scthbx[Nsctheta]/I");
  my_tree->Branch("scthcode", scthcode, "scthcode[Nsctheta][7]/I");
  my_tree->Branch("scthpos", scthpos, "scthpos[Nsctheta][7]/I");
  my_tree->Branch("scthqual", scthqual, "scthqual[Nsctheta][7]/I");
}

void DTTrigTest::beginRun(const edm::Run& iRun, const edm::EventSetup& iEventSetup) {
  my_trig->createTUs(iEventSetup);
  if (my_debug)
    cout << "[DTTrigTest] TU's Created" << endl;
}

void DTTrigTest::analyze(const Event& iEvent, const EventSetup& iEventSetup) {
  const int MAXGEN = 10;
  const float ptcut = 1.0;
  const float etacut = 2.4;
  my_trig->triggerReco(iEvent, iEventSetup);
  if (my_debug)
    cout << "[DTTrigTest] Trigger algorithm executed for run " << iEvent.id().run() << " event " << iEvent.id().event()
         << endl;

  // GENERAL Block
  runn = iEvent.id().run();
  eventn = iEvent.id().event();
  weight = 1;  // FIXME what to do with this varable?

  // GEANT Block
  Handle<vector<SimTrack> > MyTracks;
  Handle<vector<SimVertex> > MyVertexes;
  iEvent.getByLabel("g4SimHits", MyTracks);
  iEvent.getByLabel("g4SimHits", MyVertexes);
  vector<SimTrack>::const_iterator itrack;
  ngen = 0;
  if (my_debug)
    cout << "[DTTrigTest] Tracks found in the detector (not only muons) " << MyTracks->size() << endl;

  for (itrack = MyTracks->begin(); itrack != MyTracks->end(); itrack++) {
    if (abs(itrack->type()) == 13) {
      math::XYZTLorentzVectorD momentum = itrack->momentum();
      float pt = momentum.Pt();
      float eta = momentum.eta();
      if (pt > ptcut && fabs(eta) < etacut) {
        float phi = momentum.phi();
        int charge = static_cast<int>(-itrack->type() / 13);  //static_cast<int> (itrack->charge());
        if (phi < 0)
          phi = 2 * M_PI + phi;
        int vtxindex = itrack->vertIndex();
        float gvx = 0, gvy = 0, gvz = 0;
        if (vtxindex > -1) {
          gvx = MyVertexes->at(vtxindex).position().x();
          gvy = MyVertexes->at(vtxindex).position().y();
          gvz = MyVertexes->at(vtxindex).position().z();
        }
        if (ngen < MAXGEN) {
          pxgen[ngen] = momentum.x();
          pygen[ngen] = momentum.y();
          pzgen[ngen] = momentum.z();
          ptgen[ngen] = pt;
          etagen[ngen] = eta;
          phigen[ngen] = phi;
          chagen[ngen] = charge;
          vxgen[ngen] = gvx;
          vygen[ngen] = gvy;
          vzgen[ngen] = gvz;
          ngen++;
        }
      }
    }
  }

  // L1 Local Trigger Block
  // BTI
  vector<DTBtiTrigData> btitrigs = my_trig->BtiTrigs();
  vector<DTBtiTrigData>::const_iterator pbti;
  int ibti = 0;
  if (my_debug)
    cout << "[DTTrigTest] " << btitrigs.size() << " BTI triggers found" << endl;

  for (pbti = btitrigs.begin(); pbti != btitrigs.end(); pbti++) {
    if (ibti < 100) {
      bwh[ibti] = pbti->wheel();
      bstat[ibti] = pbti->station();
      bsect[ibti] = pbti->sector();
      bsl[ibti] = pbti->btiSL();
      bnum[ibti] = pbti->btiNumber();
      bbx[ibti] = pbti->step();
      bcod[ibti] = pbti->code();
      bk[ibti] = pbti->K();
      bx[ibti] = pbti->X();
      GlobalPoint pos = my_trig->CMSPosition(&(*pbti));
      GlobalVector dir = my_trig->CMSDirection(&(*pbti));
      bposx[ibti] = pos.x();
      bposy[ibti] = pos.y();
      bposz[ibti] = pos.z();
      bdirx[ibti] = dir.x();
      bdiry[ibti] = dir.y();
      bdirz[ibti] = dir.z();
      ibti++;
    }
  }
  nbti = ibti;

  //TRACO
  vector<DTTracoTrigData> tracotrigs = my_trig->TracoTrigs();
  vector<DTTracoTrigData>::const_iterator ptc;
  int itraco = 0;
  if (my_debug)
    cout << "[DTTrigTest] " << tracotrigs.size() << " TRACO triggers found" << endl;

  for (ptc = tracotrigs.begin(); ptc != tracotrigs.end(); ptc++) {
    if (itraco < 80) {
      twh[itraco] = ptc->wheel();
      tstat[itraco] = ptc->station();
      tsect[itraco] = ptc->sector();
      tnum[itraco] = ptc->tracoNumber();
      tbx[itraco] = ptc->step();
      tcod[itraco] = ptc->code();
      tk[itraco] = ptc->K();
      tx[itraco] = ptc->X();
      GlobalPoint pos = my_trig->CMSPosition(&(*ptc));
      GlobalVector dir = my_trig->CMSDirection(&(*ptc));
      tposx[itraco] = pos.x();
      tposy[itraco] = pos.y();
      tposz[itraco] = pos.z();
      tdirx[itraco] = dir.x();
      tdiry[itraco] = dir.y();
      tdirz[itraco] = dir.z();
      itraco++;
    }
  }
  ntraco = itraco;

  //TSPHI
  vector<DTChambPhSegm> tsphtrigs = my_trig->TSPhTrigs();
  vector<DTChambPhSegm>::const_iterator ptsph;
  int itsphi = 0;
  if (my_debug)
    cout << "[DTTrigTest] " << tsphtrigs.size() << " TSPhi triggers found" << endl;

  for (ptsph = tsphtrigs.begin(); ptsph != tsphtrigs.end(); ptsph++) {
    if (itsphi < 40) {
      swh[itsphi] = ptsph->wheel();
      sstat[itsphi] = ptsph->station();
      ssect[itsphi] = ptsph->sector();
      sbx[itsphi] = ptsph->step();
      scod[itsphi] = ptsph->oldCode();
      sphi[itsphi] = ptsph->phi();
      sphib[itsphi] = ptsph->phiB();
      GlobalPoint pos = my_trig->CMSPosition(&(*ptsph));
      GlobalVector dir = my_trig->CMSDirection(&(*ptsph));
      sposx[itsphi] = pos.x();
      sposy[itsphi] = pos.y();
      sposz[itsphi] = pos.z();
      sdirx[itsphi] = dir.x();
      sdiry[itsphi] = dir.y();
      sdirz[itsphi] = dir.z();
      itsphi++;
    }
  }
  ntsphi = itsphi;

  //TSTHETA
  vector<DTChambThSegm> tsthtrigs = my_trig->TSThTrigs();
  vector<DTChambThSegm>::const_iterator ptsth;
  int itstheta = 0;
  if (my_debug)
    cout << "[DTTrigTest] " << tsthtrigs.size() << " TSTheta triggers found" << endl;

  for (ptsth = tsthtrigs.begin(); ptsth != tsthtrigs.end(); ptsth++) {
    if (itstheta < 40) {
      thwh[itstheta] = ptsth->ChamberId().wheel();
      thstat[itstheta] = ptsth->ChamberId().station();
      thsect[itstheta] = ptsth->ChamberId().sector();
      thbx[itstheta] = ptsth->step();
      for (int i = 0; i < 7; i++) {
        thcode[itstheta][i] = ptsth->code(i);
        thpos[itstheta][i] = ptsth->position(i);
        thqual[itstheta][i] = ptsth->quality(i);
      }
      itstheta++;
    }
  }
  ntstheta = itstheta;

  //SCPHI
  vector<DTSectCollPhSegm> scphtrigs = my_trig->SCPhTrigs();
  vector<DTSectCollPhSegm>::const_iterator pscph;
  int iscphi = 0;
  if (my_debug)
    cout << "[DTTrigTest] " << scphtrigs.size() << " SectCollPhi triggers found" << endl;

  for (pscph = scphtrigs.begin(); pscph != scphtrigs.end(); pscph++) {
    if (iscphi < 40) {
      const DTChambPhSegm* seg = (*pscph).tsPhiTrig();
      scphwh[iscphi] = pscph->wheel();
      scphstat[iscphi] = pscph->station();
      scphsect[iscphi] = pscph->sector();
      scphbx[iscphi] = pscph->step();
      scphcod[iscphi] = pscph->oldCode();
      scphphi[iscphi] = pscph->phi();
      scphphib[iscphi] = pscph->phiB();
      GlobalPoint pos = my_trig->CMSPosition(seg);
      GlobalVector dir = my_trig->CMSDirection(seg);
      scphposx[iscphi] = pos.x();
      scphposy[iscphi] = pos.y();
      scphposz[iscphi] = pos.z();
      scphdirx[iscphi] = dir.x();
      scphdiry[iscphi] = dir.y();
      scphdirz[iscphi] = dir.z();
      iscphi++;
    }
  }
  nscphi = iscphi;

  //SCTHETA
  vector<DTSectCollThSegm> scthtrigs = my_trig->SCThTrigs();
  vector<DTSectCollThSegm>::const_iterator pscth;
  int isctheta = 0;
  if (my_debug)
    cout << "[DTTrigTest] " << scthtrigs.size() << " SectCollTheta triggers found" << endl;

  for (pscth = scthtrigs.begin(); pscth != scthtrigs.end(); pscth++) {
    if (isctheta < 40) {
      scthwh[isctheta] = pscth->ChamberId().wheel();
      scthstat[isctheta] = pscth->ChamberId().station();
      scthsect[isctheta] = pscth->ChamberId().sector();
      scthbx[isctheta] = pscth->step();
      for (int i = 0; i < 7; i++) {
        scthcode[isctheta][i] = pscth->code(i);
        scthpos[isctheta][i] = pscth->position(i);
        scthqual[isctheta][i] = pscth->quality(i);
      }
      isctheta++;
    }
  }
  nsctheta = isctheta;

  //Fill the tree
  my_tree->Fill();
}