File indexing completed on 2023-10-25 09:56:32
0001
0002
0003
0004
0005 #include <TMath.h>
0006 #include <TROOT.h>
0007 #include <TFile.h>
0008 #include <TTree.h>
0009 #include <TChain.h>
0010 #include <TBranch.h>
0011 #include <TSystem.h>
0012 #include <TH1D.h>
0013 #include <TCanvas.h>
0014 #include <TColor.h>
0015 #include <TStyle.h>
0016 #include <TSystem.h>
0017 #include <TLegend.h>
0018 #include <TLatex.h>
0019
0020 #include <iostream>
0021
0022 #include "../interface/IMATH_TrackletCalculator.h"
0023 #include "../interface/IMATH_TrackletCalculatorDisk.h"
0024 #include "../interface/IMATH_TrackletCalculatorOverlap.h"
0025
0026 #include "../interface/SLHCEvent.h"
0027 #include "../interface/Track.h"
0028 #include "../interface/Settings.h"
0029 #include "../interface/StubStreamData.h"
0030 #include "../interface/TrackletEventProcessor.h"
0031
0032 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0033 #include "DataFormats/Math/interface/deltaPhi.h"
0034
0035 #include <iomanip>
0036 #include <fstream>
0037 #include <string>
0038
0039
0040
0041
0042
0043
0044 #ifdef USEROOT
0045 #include "FPGAEvent.h"
0046 #endif
0047
0048
0049 using namespace trklet;
0050 using namespace std;
0051
0052 int main(const int argc, const char **argv) {
0053 trklet::Settings settings;
0054
0055
0056
0057
0058 settings.setFitPatternFile("../data/fitpattern.txt");
0059 settings.setProcessingModulesFile("../data/processingmodules_" + settings.geomext() + ".dat");
0060 settings.setMemoryModulesFile("../data/memorymodules_" + settings.geomext() + ".dat");
0061 settings.setWiresFile("../data/wires_" + settings.geomext() + ".dat");
0062
0063 if (settings.extended()) {
0064 settings.setTableTEDFile("../data/table_TED/table_TED_Dummy.txt");
0065 settings.setTableTREFile("../data/table_TRE/table_TRE_Dummy.txt");
0066 }
0067
0068 edm::LogVerbatim("Tracklet") << "fit pattern : " << settings.fitPatternFile();
0069 edm::LogVerbatim("Tracklet") << "process modules : " << settings.processingModulesFile();
0070 edm::LogVerbatim("Tracklet") << "memory modules : " << settings.memoryModulesFile();
0071 edm::LogVerbatim("Tracklet") << "wires : " << settings.wiresFile();
0072
0073
0074
0075 TrackletEventProcessor eventProcessor;
0076 eventProcessor.init(settings);
0077
0078 if (argc < 3)
0079 edm::LogVerbatim("Tracklet") << "Need to specify the input ascii file and the number of events to run on!";
0080
0081 int nevents = atoi(argv[2]);
0082
0083 ifstream infile;
0084 istream *in = &cin;
0085 if (strcmp(argv[1], "stdin")) {
0086 infile.open(argv[1]);
0087 in = &infile;
0088 }
0089
0090 ofstream outpars;
0091 if (settings.writeMonitorData("Pars"))
0092 outpars.open("trackpars.txt");
0093
0094
0095
0096 #ifdef USEROOT
0097 TFile *hfile = new TFile("myTest.root", "RECREATE", "Simple ROOT Ntuple");
0098 TTree *trackTree = new TTree("FPGAEvent", "L1Track Tree");
0099 FPGAEvent *fpgaEvent = new FPGAEvent;
0100 fpgaEvent->reset();
0101 trackTree->Branch("Event", &fpgaEvent);
0102 #endif
0103
0104
0105 if (settings.writeMonitorData("Seeds")) {
0106 ofstream fout("seeds.txt", ofstream::out);
0107 fout.close();
0108 }
0109
0110 for (int eventnum = 0; eventnum < nevents && !in->eof(); eventnum++) {
0111 SLHCEvent ev(*in);
0112
0113
0114
0115
0116
0117
0118 #ifdef USEROOT
0119 fpgaEvent->reset();
0120 fpgaEvent->nevt = eventnum;
0121 for (int nst = 0; nst < ev.nsimtracks(); nst++) {
0122 simtrk = ev.simtrack(nst);
0123 FPGAEventMCTrack *mcTrack = new FPGAEventMCTrack(
0124 simtrk.type(), simtrk.pt(), simtrk.eta(), simtrk.phi(), simtrk.vx(), simtrk.vy(), simtrk.vz());
0125 fpgaEvent->mcTracks.push_back(*mcTrack);
0126 }
0127 #endif
0128
0129
0130 if (settings.writeMonitorData("Seeds")) {
0131 ofstream fout("seeds.txt", ofstream::app);
0132 fout << "======== Event " << eventnum << " ========" << endl;
0133 for (unsigned nst = 0; nst < ev.nsimtracks(); nst++) {
0134 const L1SimTrack &simtrk = ev.simtrack(nst);
0135 fout << "SimTrk " << simtrk.pt() << " " << simtrk.eta() << " " << simtrk.phi() << " " << simtrk.d0() << " ";
0136
0137 vector<string> hitPattern;
0138 for (int i = 0; i < ev.nstubs(); i++) {
0139 const L1TStub &stub = ev.stub(i);
0140 if (!stub.tpmatch(simtrk.trackid()))
0141 continue;
0142 if (stub.layer() < 999) {
0143 switch (stub.layer()) {
0144 case 0:
0145 hitPattern.emplace_back("L1");
0146 break;
0147 case 1:
0148 hitPattern.emplace_back("L2");
0149 break;
0150 case 2:
0151 hitPattern.emplace_back("L3");
0152 break;
0153 case 3:
0154 hitPattern.emplace_back("L4");
0155 break;
0156 case 4:
0157 hitPattern.emplace_back("L5");
0158 break;
0159 case 5:
0160 hitPattern.emplace_back("L6");
0161 break;
0162 default:
0163 edm::LogVerbatim("Tracklet") << "Stub layer: " << stub.layer();
0164 assert(0);
0165 }
0166 } else {
0167 string d = (stub.isPSmodule() ? "D" : "d");
0168 switch (abs(stub.disk())) {
0169 case 1:
0170 hitPattern.push_back(d + "1");
0171 break;
0172 case 2:
0173 hitPattern.push_back(d + "2");
0174 break;
0175 case 3:
0176 hitPattern.push_back(d + "3");
0177 break;
0178 case 4:
0179 hitPattern.push_back(d + "4");
0180 break;
0181 case 5:
0182 hitPattern.push_back(d + "5");
0183 break;
0184 default:
0185 edm::LogVerbatim("Tracklet") << "Stub disk: " << stub.disk();
0186 assert(0);
0187 }
0188 }
0189 }
0190 bool (*compare)(const string &, const string &) = [](const string &a, const string &b) -> bool {
0191 if (a.at(0) == 'L' && b.at(0) == 'D')
0192 return true;
0193 else if (a.at(0) == 'D' && b.at(0) == 'L')
0194 return false;
0195 else
0196 return a.at(1) < b.at(1);
0197 };
0198 sort(hitPattern.begin(), hitPattern.end(), compare);
0199 hitPattern.erase(unique(hitPattern.begin(), hitPattern.end()), hitPattern.end());
0200 for (const auto &stub : hitPattern)
0201 fout << stub;
0202 if (hitPattern.empty())
0203 fout << "XX";
0204 fout << endl;
0205 }
0206 fout.close();
0207 }
0208
0209 edm::LogVerbatim("Tracklet") << "Process event: " << eventnum << " with " << ev.nstubs() << " stubs and "
0210 << ev.nsimtracks() << " simtracks";
0211
0212
0213 constexpr unsigned int numStubStreamsPerTrack = settings.extended() ? N_SEED : N_SEED_PROMPT;
0214
0215 constexpr unsigned int maxNumProjectionLayers = 8;
0216 constexpr unsigned int numStreamsTrack = N_SECTOR * numStubStreamsTrack;
0217 constexpr unsigned int numStreamsStub = numStreamsTrack * maxNumProjectionLayers;
0218 std::vector<std::vector<std::string>> tracksStream(numStreamsTrack);
0219 std::vector<std::vector<StubStreamData>> stubsStream(numStreamsStub);
0220
0221 eventProcessor.event(ev, tracksStream, stubsStream);
0222
0223 const std::vector<Track> &tracks = eventProcessor.tracks();
0224
0225
0226
0227 #ifdef USEROOT
0228 #include "FPGATree.icc"
0229 #endif
0230
0231
0232 if (settings.writeMonitorData("MatchEff")) {
0233 static ofstream out("matcheff.txt");
0234 int nsim = 0;
0235 for (unsigned int isimtrack = 0; isimtrack < ev.nsimtracks(); isimtrack++) {
0236 L1SimTrack simtrack = ev.simtrack(isimtrack);
0237 if (simtrack.pt() < 2.0)
0238 continue;
0239 if (std::abs(simtrack.eta()) > 2.4)
0240 continue;
0241 if (std::abs(simtrack.vz()) > 15.0)
0242 continue;
0243 if (hypot(simtrack.vx(), simtrack.vy()) > 0.1)
0244 continue;
0245 bool electron = (abs(simtrack.type()) == 11);
0246 bool muon = (abs(simtrack.type()) == 13);
0247 bool pion = (abs(simtrack.type()) == 211);
0248 bool kaon = (abs(simtrack.type()) == 321);
0249 bool proton = (abs(simtrack.type()) == 2212);
0250 if (!(electron || muon || pion || kaon || proton))
0251 continue;
0252 int nlayers = 0;
0253 int ndisks = 0;
0254 int simeventid = simtrack.eventid();
0255 int simtrackid = simtrack.trackid();
0256 ev.layersHit(simtrackid, nlayers, ndisks);
0257 if (nlayers + ndisks < 4)
0258 continue;
0259 nsim++;
0260 for (int seed = -1; seed < 12; seed++) {
0261 bool eff = false;
0262 bool effloose = false;
0263 int itrackmatch = -1;
0264 for (unsigned int itrack = 0; itrack < tracks.size(); itrack++) {
0265 const std::vector<L1TStub> &stubs = tracks[itrack].stubs();
0266 if (seed == -1) {
0267 if (tracks[itrack].duplicate())
0268 continue;
0269 } else {
0270 if (seed != tracks[itrack].seed())
0271 continue;
0272 }
0273
0274 unsigned int nmatch = 0;
0275 for (auto stub : stubs) {
0276 if (stub.tpmatch(simtrackid)) {
0277 nmatch++;
0278 }
0279 }
0280
0281 if (nmatch == stubs.size()) {
0282 eff = true;
0283 itrackmatch = itrack;
0284 }
0285 if (nmatch >= stubs.size() - 1) {
0286 effloose = true;
0287 if (!eff)
0288 itrackmatch = itrack;
0289 }
0290 }
0291 double dpt = -999;
0292 double dphi = -999;
0293 double deta = -999;
0294 double dz0 = -999;
0295 int q = 1;
0296 if (simtrack.type() == 11 || simtrack.type() == 13 || simtrack.type() == -211 || simtrack.type() == -321 ||
0297 simtrack.type() == -2212) {
0298 q = -1;
0299 }
0300
0301 if (itrackmatch >= 0) {
0302 dpt = tracks[itrackmatch].pt(settings) - q * simtrack.pt();
0303 dphi = tracks[itrackmatch].phi0(settings) - simtrack.phi();
0304 if (dphi > M_PI)
0305 dphi -= 2 * M_PI;
0306 if (dphi < -M_PI)
0307 dphi += 2 * M_PI;
0308 deta = tracks[itrackmatch].eta(settings) - simtrack.eta();
0309 dz0 = tracks[itrackmatch].z0(settings) - simtrack.vz();
0310 }
0311
0312 out << eventnum << " " << simeventid << " " << seed << " " << simtrackid << " " << simtrack.type() << " "
0313 << simtrack.pt() << " " << simtrack.eta() << " " << simtrack.phi() << " " << simtrack.vx() << " "
0314 << simtrack.vy() << " " << simtrack.vz() << " " << eff << " " << effloose << " " << dpt << " " << dphi
0315 << " " << deta << " " << dz0 << endl;
0316 }
0317 }
0318 }
0319
0320 int ntrack = 0;
0321 int i = 0;
0322 for (auto &track : tracks) {
0323 i++;
0324 if (settings.writeMonitorData("Pars")) {
0325 outpars << track.duplicate() << " " << track.eta(settings) << " " << track.phi0(settings) << " "
0326 << track.z0(settings) << " " << angle0to2pi::make0To2pi(track.phi0(settings)) / (2 * M_PI / N_SECTOR)
0327 << " " << track.rinv(settings);
0328 }
0329 if (!track.duplicate()) {
0330 ntrack++;
0331 }
0332 }
0333
0334 edm::LogVerbatim("Tracklet") << "Number of found tracks : " << tracks.size() << " unique " << ntrack << " ";
0335 }
0336
0337 eventProcessor.printSummary();
0338
0339 exit(0);
0340 }
