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