File indexing completed on 2024-04-06 12:19:44
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017 #include "L1Trigger/DTTrigger/interface/DTTrigTest.h"
0018
0019
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
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
0033 #include "TROOT.h"
0034 #include "TTree.h"
0035 #include "TFile.h"
0036
0037
0038 #include "DataFormats/Math/interface/LorentzVector.h"
0039 #include <CLHEP/Vector/LorentzVector.h>
0040
0041
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
0081 void DTTrigTest::beginJob() {
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096
0097
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
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
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
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
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
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
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
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
0212 runn = iEvent.id().run();
0213 eventn = iEvent.id().event();
0214 weight = 1;
0215
0216
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);
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
0261
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
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
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
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
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
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
0428 my_tree->Fill();
0429 }