Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-03-26 02:43:47

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/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 //Uncomment if you want root output
0039 //#define USEROOT
0040 
0041 // Include file to define ROOT-Tree
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   // these are options that are read from python configuration files for the CMSSW running, set manually for the standalone version
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     // Open file to hold ROOT-Tree
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     //auto simtrk = ev.simtrack(0);
0113     //if (std::abs(std::abs(simtrk.eta())-1.3)>0.1) continue;
0114 
0115     // -----------------------------------------------------------------
0116     // setup ROOT Tree and Add Monte Carlo tracks to the ROOT-Tree Event
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     // Block for producing ROOT-Tree
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 }