Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:19:44

0001 //-------------------------------------------------
0002 //
0003 /**  \class DTTrigTest
0004  *
0005  *   EDAnalyzer that generates a rootfile useful
0006  *   for L1-DTTrigger debugging and performance 
0007  *   studies
0008  *
0009  *
0010  *
0011  *   \author C. Battilana
0012  */
0013 //
0014 //--------------------------------------------------
0015 
0016 // This class's header
0017 #include "L1Trigger/DTTrigger/interface/DTTrigTest.h"
0018 
0019 // Framework related classes
0020 #include "FWCore/Framework/interface/EventSetup.h"
0021 #include "FWCore/Framework/interface/ESHandle.h"
0022 
0023 #include "L1TriggerConfig/DTTPGConfig/interface/DTConfigManager.h"
0024 #include "L1TriggerConfig/DTTPGConfig/interface/DTConfigManagerRcd.h"
0025 
0026 // Trigger and DataFormats headers
0027 #include "L1Trigger/DTSectorCollector/interface/DTSectCollPhSegm.h"
0028 #include "L1Trigger/DTSectorCollector/interface/DTSectCollThSegm.h"
0029 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0030 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0031 
0032 // ROOT headers
0033 #include "TROOT.h"
0034 #include "TTree.h"
0035 #include "TFile.h"
0036 
0037 // Collaborating classes
0038 #include "DataFormats/Math/interface/LorentzVector.h"
0039 #include <CLHEP/Vector/LorentzVector.h>
0040 
0041 // C++ headers
0042 #include <iostream>
0043 #include <cmath>
0044 #include <ctime>
0045 
0046 using namespace std;
0047 using namespace edm;
0048 
0049 const double DTTrigTest::my_TtoTDC = 32. / 25.;
0050 
0051 DTTrigTest::DTTrigTest(const ParameterSet& pset) : my_trig(nullptr) {
0052   my_debug = pset.getUntrackedParameter<bool>("debug");
0053   string outputfile = pset.getUntrackedParameter<string>("outputFileName");
0054   if (my_debug)
0055     cout << "[DTTrigTest] Creating rootfile " << outputfile << endl;
0056   my_rootfile = new TFile(outputfile.c_str(), "RECREATE");
0057   my_tree = new TTree("h1", "GMT", 0);
0058   my_params = pset;
0059   if (my_debug)
0060     cout << "[DTTrigTest] Constructor executed!!!" << endl;
0061   my_trig = new DTTrig(my_params, consumesCollector());
0062 }
0063 
0064 DTTrigTest::~DTTrigTest() {
0065   if (my_trig != nullptr)
0066     delete my_trig;
0067   delete my_rootfile;
0068   if (my_debug)
0069     cout << "[DTTrigTest] Destructor executed!!!" << endl;
0070 }
0071 
0072 void DTTrigTest::endJob() {
0073   if (my_debug)
0074     cout << "[DTTrigTest] Writing Tree and Closing File" << endl;
0075   my_tree->Write();
0076   delete my_tree;
0077   my_rootfile->Close();
0078 }
0079 
0080 //void DTTrigTest::beginJob(const EventSetup & iEventSetup){
0081 void DTTrigTest::beginJob() {
0082   // get DTConfigManager
0083   // ESHandle< DTConfigManager > confManager ;
0084   // iEventSetup.get< DTConfigManagerRcd >().get( confManager ) ;
0085 
0086   //for testing purpose....
0087   //DTBtiId btiid(1,1,1,1,1);
0088   //confManager->getDTConfigBti(btiid)->print();
0089 
0090   //   my_trig = new DTTrig(my_params);
0091 
0092   //   my_trig->createTUs(iEventSetup);
0093   //   if (my_debug)
0094   //     cout << "[DTTrigTest] TU's Created" << endl;
0095 
0096   // BOOKING of the tree's varables
0097   // GENERAL block branches
0098   my_tree->Branch("Run", &runn, "Run/I");
0099   my_tree->Branch("Event", &eventn, "Event/I");
0100   my_tree->Branch("Weight", &weight, "Weight/F");
0101   // GEANT block branches
0102   my_tree->Branch("Ngen", &ngen, "Ngen/I");
0103   my_tree->Branch("Pxgen", pxgen, "Pxgen[Ngen]/F");
0104   my_tree->Branch("Pygen", pygen, "Pygen[Ngen]/F");
0105   my_tree->Branch("Pzgen", pzgen, "Pzgen[Ngen]/F");
0106   my_tree->Branch("Ptgen", ptgen, "Ptgen[Ngen]/F");
0107   my_tree->Branch("Etagen", etagen, "Etagen[Ngen]/F");
0108   my_tree->Branch("Phigen", phigen, "Phigen[Ngen]/F");
0109   my_tree->Branch("Chagen", chagen, "Chagen[Ngen]/I");
0110   my_tree->Branch("Vxgen", vxgen, "Vxgen[Ngen]/F");
0111   my_tree->Branch("Vygen", vygen, "Vygen[Ngen]/F");
0112   my_tree->Branch("Vzgen", vzgen, "Vzgen[Ngen]/F");
0113   // L1MuDTBtiChipS block
0114   my_tree->Branch("Nbti", &nbti, "Nbti/I");
0115   my_tree->Branch("bwh", bwh, "bwh[Nbti]/I");
0116   my_tree->Branch("bstat", bstat, "bstat[Nbti]/I");
0117   my_tree->Branch("bsect", bsect, "bsect[Nbti]/I");
0118   my_tree->Branch("bsl", bsl, "bsl[Nbti]/I");
0119   my_tree->Branch("bnum", bnum, "bnum[Nbti]/I");
0120   my_tree->Branch("bbx", bbx, "bbx[Nbti]/I");
0121   my_tree->Branch("bcod", bcod, "bcod[Nbti]/I");
0122   my_tree->Branch("bk", bk, "bk[Nbti]/I");
0123   my_tree->Branch("bx", bx, "bx[Nbti]/I");
0124   my_tree->Branch("bposx", bposx, "bposx[Nbti]/F");
0125   my_tree->Branch("bposy", bposy, "bposy[Nbti]/F");
0126   my_tree->Branch("bposz", bposz, "bposz[Nbti]/F");
0127   my_tree->Branch("bdirx", bdirx, "bdirx[Nbti]/F");
0128   my_tree->Branch("bdiry", bdiry, "bdiry[Nbti]/F");
0129   my_tree->Branch("bdirz", bdirz, "bdirz[Nbti]/F");
0130   // L1MuDTTracoChipS block
0131   my_tree->Branch("Ntraco", &ntraco, "Ntraco/I");
0132   my_tree->Branch("twh", twh, "twh[Ntraco]/I");
0133   my_tree->Branch("tstat", tstat, "tstat[Ntraco]/I");
0134   my_tree->Branch("tsect", tsect, "tsect[Ntraco]/I");
0135   my_tree->Branch("tnum", tnum, "tnum[Ntraco]/I");
0136   my_tree->Branch("tbx", tbx, "tbx[Ntraco]/I");
0137   my_tree->Branch("tcod", tcod, "tcod[Ntraco]/I");
0138   my_tree->Branch("tk", tk, "tk[Ntraco]/I");
0139   my_tree->Branch("tx", tx, "tx[Ntraco]/I");
0140   my_tree->Branch("tposx", tposx, "tposx[Ntraco]/F");
0141   my_tree->Branch("tposy", tposy, "tposy[Ntraco]/F");
0142   my_tree->Branch("tposz", tposz, "tposz[Ntraco]/F");
0143   my_tree->Branch("tdirx", tdirx, "tdirx[Ntraco]/F");
0144   my_tree->Branch("tdiry", tdiry, "tdiry[Ntraco]/F");
0145   my_tree->Branch("tdirz", tdirz, "tdirz[Ntraco]/F");
0146   // TSPHI block
0147   my_tree->Branch("Ntsphi", &ntsphi, "Ntsphi/I");
0148   my_tree->Branch("swh", swh, "swh[Ntsphi]/I");
0149   my_tree->Branch("sstat", sstat, "sstat[Ntsphi]/I");
0150   my_tree->Branch("ssect", ssect, "ssect[Ntsphi]/I");
0151   my_tree->Branch("sbx", sbx, "sbx[Ntsphi]/I");
0152   my_tree->Branch("scod", scod, "scod[Ntsphi]/I");
0153   my_tree->Branch("sphi", sphi, "sphi[Ntsphi]/I");
0154   my_tree->Branch("sphib", sphib, "sphib[Ntsphi]/I");
0155   my_tree->Branch("sposx", sposx, "sposx[Ntsphi]/F");
0156   my_tree->Branch("sposy", sposy, "sposy[Ntsphi]/F");
0157   my_tree->Branch("sposz", sposz, "sposz[Ntsphi]/F");
0158   my_tree->Branch("sdirx", sdirx, "sdirx[Ntsphi]/F");
0159   my_tree->Branch("sdiry", sdiry, "sdiry[Ntsphi]/F");
0160   my_tree->Branch("sdirz", sdirz, "sdirz[Ntsphi]/F");
0161   // TSTHETA block
0162   my_tree->Branch("Ntstheta", &ntstheta, "Ntstheta/I");
0163   my_tree->Branch("thwh", thwh, "thwh[Ntstheta]/I");
0164   my_tree->Branch("thstat", thstat, "thstat[Ntstheta]/I");
0165   my_tree->Branch("thsect", thsect, "thsect[Ntstheta]/I");
0166   my_tree->Branch("thbx", thbx, "thbx[Ntstheta]/I");
0167   my_tree->Branch("thcode", thcode, "thcode[Ntstheta][7]/I");
0168   my_tree->Branch("thpos", thpos, "thpos[Ntstheta][7]/I");
0169   my_tree->Branch("thqual", thqual, "thqual[Ntstheta][7]/I");
0170   // SC PHI block
0171   my_tree->Branch("Nscphi", &nscphi, "Nscphi/I");
0172   my_tree->Branch("scphwh", scphwh, "scphwh[Nscphi]/I");
0173   my_tree->Branch("scphstat", scphstat, "scphstat[Nscphi]/I");
0174   my_tree->Branch("scphsect", scphsect, "scphsect[Nscphi]/I");
0175   my_tree->Branch("scphbx", scphbx, "scphbx[Nscphi]/I");
0176   my_tree->Branch("scphcod", scphcod, "scphcod[Nscphi]/I");
0177   my_tree->Branch("scphphi", scphphi, "scphphi[Nscphi]/I");
0178   my_tree->Branch("scphphib", scphphib, "scphphib[Nscphi]/I");
0179   my_tree->Branch("scphposx", scphposx, "scphposx[Nscphi]/F");
0180   my_tree->Branch("scphposy", scphposy, "scphposy[Nscphi]/F");
0181   my_tree->Branch("scphposz", scphposz, "scphposz[Nscphi]/F");
0182   my_tree->Branch("scphdirx", scphdirx, "scphdirx[Nscphi]/F");
0183   my_tree->Branch("scphdiry", scphdiry, "scphdiry[Nscphi]/F");
0184   my_tree->Branch("scphdirz", scphdirz, "scphdirz[Nscphi]/F");
0185   // SC THETA block
0186   my_tree->Branch("Nsctheta", &nsctheta, "Nsctheta/I");
0187   my_tree->Branch("scthwh", scthwh, "scthwh[Nsctheta]/I");
0188   my_tree->Branch("scthstat", scthstat, "scthstat[Nsctheta]/I");
0189   my_tree->Branch("scthsect", scthsect, "scthsect[Nsctheta]/I");
0190   my_tree->Branch("scthbx", scthbx, "scthbx[Nsctheta]/I");
0191   my_tree->Branch("scthcode", scthcode, "scthcode[Nsctheta][7]/I");
0192   my_tree->Branch("scthpos", scthpos, "scthpos[Nsctheta][7]/I");
0193   my_tree->Branch("scthqual", scthqual, "scthqual[Nsctheta][7]/I");
0194 }
0195 
0196 void DTTrigTest::beginRun(const edm::Run& iRun, const edm::EventSetup& iEventSetup) {
0197   my_trig->createTUs(iEventSetup);
0198   if (my_debug)
0199     cout << "[DTTrigTest] TU's Created" << endl;
0200 }
0201 
0202 void DTTrigTest::analyze(const Event& iEvent, const EventSetup& iEventSetup) {
0203   const int MAXGEN = 10;
0204   const float ptcut = 1.0;
0205   const float etacut = 2.4;
0206   my_trig->triggerReco(iEvent, iEventSetup);
0207   if (my_debug)
0208     cout << "[DTTrigTest] Trigger algorithm executed for run " << iEvent.id().run() << " event " << iEvent.id().event()
0209          << endl;
0210 
0211   // GENERAL Block
0212   runn = iEvent.id().run();
0213   eventn = iEvent.id().event();
0214   weight = 1;  // FIXME what to do with this varable?
0215 
0216   // GEANT Block
0217   Handle<vector<SimTrack> > MyTracks;
0218   Handle<vector<SimVertex> > MyVertexes;
0219   iEvent.getByLabel("g4SimHits", MyTracks);
0220   iEvent.getByLabel("g4SimHits", MyVertexes);
0221   vector<SimTrack>::const_iterator itrack;
0222   ngen = 0;
0223   if (my_debug)
0224     cout << "[DTTrigTest] Tracks found in the detector (not only muons) " << MyTracks->size() << endl;
0225 
0226   for (itrack = MyTracks->begin(); itrack != MyTracks->end(); itrack++) {
0227     if (abs(itrack->type()) == 13) {
0228       math::XYZTLorentzVectorD momentum = itrack->momentum();
0229       float pt = momentum.Pt();
0230       float eta = momentum.eta();
0231       if (pt > ptcut && fabs(eta) < etacut) {
0232         float phi = momentum.phi();
0233         int charge = static_cast<int>(-itrack->type() / 13);  //static_cast<int> (itrack->charge());
0234         if (phi < 0)
0235           phi = 2 * M_PI + phi;
0236         int vtxindex = itrack->vertIndex();
0237         float gvx = 0, gvy = 0, gvz = 0;
0238         if (vtxindex > -1) {
0239           gvx = MyVertexes->at(vtxindex).position().x();
0240           gvy = MyVertexes->at(vtxindex).position().y();
0241           gvz = MyVertexes->at(vtxindex).position().z();
0242         }
0243         if (ngen < MAXGEN) {
0244           pxgen[ngen] = momentum.x();
0245           pygen[ngen] = momentum.y();
0246           pzgen[ngen] = momentum.z();
0247           ptgen[ngen] = pt;
0248           etagen[ngen] = eta;
0249           phigen[ngen] = phi;
0250           chagen[ngen] = charge;
0251           vxgen[ngen] = gvx;
0252           vygen[ngen] = gvy;
0253           vzgen[ngen] = gvz;
0254           ngen++;
0255         }
0256       }
0257     }
0258   }
0259 
0260   // L1 Local Trigger Block
0261   // BTI
0262   vector<DTBtiTrigData> btitrigs = my_trig->BtiTrigs();
0263   vector<DTBtiTrigData>::const_iterator pbti;
0264   int ibti = 0;
0265   if (my_debug)
0266     cout << "[DTTrigTest] " << btitrigs.size() << " BTI triggers found" << endl;
0267 
0268   for (pbti = btitrigs.begin(); pbti != btitrigs.end(); pbti++) {
0269     if (ibti < 100) {
0270       bwh[ibti] = pbti->wheel();
0271       bstat[ibti] = pbti->station();
0272       bsect[ibti] = pbti->sector();
0273       bsl[ibti] = pbti->btiSL();
0274       bnum[ibti] = pbti->btiNumber();
0275       bbx[ibti] = pbti->step();
0276       bcod[ibti] = pbti->code();
0277       bk[ibti] = pbti->K();
0278       bx[ibti] = pbti->X();
0279       GlobalPoint pos = my_trig->CMSPosition(&(*pbti));
0280       GlobalVector dir = my_trig->CMSDirection(&(*pbti));
0281       bposx[ibti] = pos.x();
0282       bposy[ibti] = pos.y();
0283       bposz[ibti] = pos.z();
0284       bdirx[ibti] = dir.x();
0285       bdiry[ibti] = dir.y();
0286       bdirz[ibti] = dir.z();
0287       ibti++;
0288     }
0289   }
0290   nbti = ibti;
0291 
0292   //TRACO
0293   vector<DTTracoTrigData> tracotrigs = my_trig->TracoTrigs();
0294   vector<DTTracoTrigData>::const_iterator ptc;
0295   int itraco = 0;
0296   if (my_debug)
0297     cout << "[DTTrigTest] " << tracotrigs.size() << " TRACO triggers found" << endl;
0298 
0299   for (ptc = tracotrigs.begin(); ptc != tracotrigs.end(); ptc++) {
0300     if (itraco < 80) {
0301       twh[itraco] = ptc->wheel();
0302       tstat[itraco] = ptc->station();
0303       tsect[itraco] = ptc->sector();
0304       tnum[itraco] = ptc->tracoNumber();
0305       tbx[itraco] = ptc->step();
0306       tcod[itraco] = ptc->code();
0307       tk[itraco] = ptc->K();
0308       tx[itraco] = ptc->X();
0309       GlobalPoint pos = my_trig->CMSPosition(&(*ptc));
0310       GlobalVector dir = my_trig->CMSDirection(&(*ptc));
0311       tposx[itraco] = pos.x();
0312       tposy[itraco] = pos.y();
0313       tposz[itraco] = pos.z();
0314       tdirx[itraco] = dir.x();
0315       tdiry[itraco] = dir.y();
0316       tdirz[itraco] = dir.z();
0317       itraco++;
0318     }
0319   }
0320   ntraco = itraco;
0321 
0322   //TSPHI
0323   vector<DTChambPhSegm> tsphtrigs = my_trig->TSPhTrigs();
0324   vector<DTChambPhSegm>::const_iterator ptsph;
0325   int itsphi = 0;
0326   if (my_debug)
0327     cout << "[DTTrigTest] " << tsphtrigs.size() << " TSPhi triggers found" << endl;
0328 
0329   for (ptsph = tsphtrigs.begin(); ptsph != tsphtrigs.end(); ptsph++) {
0330     if (itsphi < 40) {
0331       swh[itsphi] = ptsph->wheel();
0332       sstat[itsphi] = ptsph->station();
0333       ssect[itsphi] = ptsph->sector();
0334       sbx[itsphi] = ptsph->step();
0335       scod[itsphi] = ptsph->oldCode();
0336       sphi[itsphi] = ptsph->phi();
0337       sphib[itsphi] = ptsph->phiB();
0338       GlobalPoint pos = my_trig->CMSPosition(&(*ptsph));
0339       GlobalVector dir = my_trig->CMSDirection(&(*ptsph));
0340       sposx[itsphi] = pos.x();
0341       sposy[itsphi] = pos.y();
0342       sposz[itsphi] = pos.z();
0343       sdirx[itsphi] = dir.x();
0344       sdiry[itsphi] = dir.y();
0345       sdirz[itsphi] = dir.z();
0346       itsphi++;
0347     }
0348   }
0349   ntsphi = itsphi;
0350 
0351   //TSTHETA
0352   vector<DTChambThSegm> tsthtrigs = my_trig->TSThTrigs();
0353   vector<DTChambThSegm>::const_iterator ptsth;
0354   int itstheta = 0;
0355   if (my_debug)
0356     cout << "[DTTrigTest] " << tsthtrigs.size() << " TSTheta triggers found" << endl;
0357 
0358   for (ptsth = tsthtrigs.begin(); ptsth != tsthtrigs.end(); ptsth++) {
0359     if (itstheta < 40) {
0360       thwh[itstheta] = ptsth->ChamberId().wheel();
0361       thstat[itstheta] = ptsth->ChamberId().station();
0362       thsect[itstheta] = ptsth->ChamberId().sector();
0363       thbx[itstheta] = ptsth->step();
0364       for (int i = 0; i < 7; i++) {
0365         thcode[itstheta][i] = ptsth->code(i);
0366         thpos[itstheta][i] = ptsth->position(i);
0367         thqual[itstheta][i] = ptsth->quality(i);
0368       }
0369       itstheta++;
0370     }
0371   }
0372   ntstheta = itstheta;
0373 
0374   //SCPHI
0375   vector<DTSectCollPhSegm> scphtrigs = my_trig->SCPhTrigs();
0376   vector<DTSectCollPhSegm>::const_iterator pscph;
0377   int iscphi = 0;
0378   if (my_debug)
0379     cout << "[DTTrigTest] " << scphtrigs.size() << " SectCollPhi triggers found" << endl;
0380 
0381   for (pscph = scphtrigs.begin(); pscph != scphtrigs.end(); pscph++) {
0382     if (iscphi < 40) {
0383       const DTChambPhSegm* seg = (*pscph).tsPhiTrig();
0384       scphwh[iscphi] = pscph->wheel();
0385       scphstat[iscphi] = pscph->station();
0386       scphsect[iscphi] = pscph->sector();
0387       scphbx[iscphi] = pscph->step();
0388       scphcod[iscphi] = pscph->oldCode();
0389       scphphi[iscphi] = pscph->phi();
0390       scphphib[iscphi] = pscph->phiB();
0391       GlobalPoint pos = my_trig->CMSPosition(seg);
0392       GlobalVector dir = my_trig->CMSDirection(seg);
0393       scphposx[iscphi] = pos.x();
0394       scphposy[iscphi] = pos.y();
0395       scphposz[iscphi] = pos.z();
0396       scphdirx[iscphi] = dir.x();
0397       scphdiry[iscphi] = dir.y();
0398       scphdirz[iscphi] = dir.z();
0399       iscphi++;
0400     }
0401   }
0402   nscphi = iscphi;
0403 
0404   //SCTHETA
0405   vector<DTSectCollThSegm> scthtrigs = my_trig->SCThTrigs();
0406   vector<DTSectCollThSegm>::const_iterator pscth;
0407   int isctheta = 0;
0408   if (my_debug)
0409     cout << "[DTTrigTest] " << scthtrigs.size() << " SectCollTheta triggers found" << endl;
0410 
0411   for (pscth = scthtrigs.begin(); pscth != scthtrigs.end(); pscth++) {
0412     if (isctheta < 40) {
0413       scthwh[isctheta] = pscth->ChamberId().wheel();
0414       scthstat[isctheta] = pscth->ChamberId().station();
0415       scthsect[isctheta] = pscth->ChamberId().sector();
0416       scthbx[isctheta] = pscth->step();
0417       for (int i = 0; i < 7; i++) {
0418         scthcode[isctheta][i] = pscth->code(i);
0419         scthpos[isctheta][i] = pscth->position(i);
0420         scthqual[isctheta][i] = pscth->quality(i);
0421       }
0422       isctheta++;
0423     }
0424   }
0425   nsctheta = isctheta;
0426 
0427   //Fill the tree
0428   my_tree->Fill();
0429 }