Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-10-14 01:44:15

0001 // ----------------------------------------------------------------
0002 // STANDALONE (non-CMSSW) producer for L1 tracking
0003 // ----------------------------------------------------------------
0004 // ROOT includes
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 //Uncomment if you want root output
0040 //#define USEROOT
0041 
0042 // Include file to define ROOT-Tree
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   // these are options that are read from python configuration files for the CMSSW running, set manually for the standalone version
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     // Open file to hold ROOT-Tree
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     //auto simtrk = ev.simtrack(0);
0114     //if (std::abs(std::abs(simtrk.eta())-1.3)>0.1) continue;
0115 
0116     // -----------------------------------------------------------------
0117     // setup ROOT Tree and Add Monte Carlo tracks to the ROOT-Tree Event
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     // Output track streams per nonant from track-builders.
0213     constexpr unsigned int numStubStreamsPerTrack = settings.extended() ? N_SEED : N_SEED_PROMPT;
0214     // Max. number of projection layers for any tracklet seed.
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     // Block for producing ROOT-Tree
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 }