CMS 3D CMS Logo

Typedefs | Enumerations | Functions | Variables
WriteMemoryFile.cc File Reference
#include "TFile.h"
#include "TTree.h"
#include <iostream>
#include <list>
#include <unordered_map>
#include "RecoTracker/MkFitCore/standalone/Event.h"
#include "RecoTracker/MkFitCore/interface/Config.h"
#include "RecoTracker/MkFitCore/interface/TrackerInfo.h"
#include "RecoTracker/MkFitCMS/interface/LayerNumberConverter.h"

Go to the source code of this file.

Typedefs

typedef lStr_t::iterator lStr_i
 
typedef std::list< std::string > lStr_t
 
using TrackAlgorithm = TrackBase::TrackAlgorithm
 

Enumerations

enum  HitType {
  HitType::Pixel = 0, HitType::Strip = 1, HitType::Glued = 2, HitType::Invalid = 3,
  HitType::Phase2OT = 4, HitType::Unknown = 99
}
 

Functions

int main (int argc, char *argv[])
 
bool next_arg_option (lStr_t &args, lStr_i &i)
 
void next_arg_or_die (lStr_t &args, lStr_i &i)
 
void printHelp (const char *av0)
 

Variables

constexpr int cleanSimTrack_minRecHits = 2
 
constexpr int cleanSimTrack_minSimHits = 3
 
constexpr bool useMatched = false
 

Typedef Documentation

◆ lStr_i

typedef lStr_t::iterator lStr_i

Definition at line 25 of file WriteMemoryFile.cc.

◆ lStr_t

typedef std::list<std::string> lStr_t

Definition at line 24 of file WriteMemoryFile.cc.

◆ TrackAlgorithm

Definition at line 14 of file WriteMemoryFile.cc.

Enumeration Type Documentation

◆ HitType

enum HitType
strong
Enumerator
Pixel 
Strip 
Glued 
Invalid 
Phase2OT 
Unknown 

Definition at line 22 of file WriteMemoryFile.cc.

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)

===============================================================================================================================================================================================


variant2: for each run define phi-averaged A for normalization channel (Dref,16) and then, divide Rijk on it, i.e. get RRijk




















































































eta=27

eta=25

eta=23

eta=22

eta=21

eta=26

eta=24

eta=19

eta=17

eta=25

eta=23

eta=22

eta=21

eta=26

eta=24

eta=20

eta=19

eta=18

eta=27 L1=1

eta=25 L1=1

eta=23 L1=1

eta=22 L1=1

eta=21 L1=1

eta=29 L1=1

eta=26 L1=1

eta=24 L1=1

eta=20 L1=1

eta=19 L1=1

eta=18 L1=1

eta=17 L1=1

eta=28 L7=1

eta=27 L7=1

eta=25 L7=1

eta=23 L7=1

eta=22 L7=1

eta=21 L7=1

eta=26 L7=1

eta=24 L7=1

eta=20 L7=1

eta=19 L7=1

eta=18 L7=1

eta=17 L7=1

eta=27

eta=25

eta=23

eta=22

eta=21

eta=26

eta=24

eta=19

eta=17

eta=25

eta=23

eta=22

eta=21

eta=26

eta=24

eta=20

eta=19

eta=18

eta=27 L1=1

eta=25 L1=1

eta=23 L1=1

eta=22 L1=1

eta=21 L1=1

eta=26 L1=1

eta=24 L1=1

eta=20 L1=1

eta=19 L1=1

eta=18 L1=1

eta=17 L1=1

eta=28 L7=1

eta=27 L7=1

eta=25 L7=1

eta=23 L7=1

eta=22 L7=1

eta=21 L7=1

eta=26 L7=1

eta=24 L7=1

eta=20 L7=1

eta=19 L7=1

eta=18 L7=1

eta=17 L7=1

eta=27

eta=28

errA with average Amplitudes

errA with average Amplitudes

errA with average Amplitudes

errA with average Amplitudes

Summed Amplitude Plots:





Summed Amplitude Plots:

Summed Amplitude Plots:

Summed Amplitude Plots:

Summed Amplitude Plots:

Summed Amplitude Plots:

RBX:

errA with average Amplitudes

errA with average Amplitudes

errA with average Amplitudes

errA with average Amplitudes

Summed Amplitude Plots:





Summed Amplitude Plots:

Summed Amplitude Plots:

Summed Amplitude Plots:

Summed Amplitude Plots:

Summed Amplitude Plots:

RBX:

errA with average Amplitudes

errA with average Amplitudes

errA with average Amplitudes

errA with average Amplitudes

Summed Amplitude Plots:





Summed Amplitude Plots:

Summed Amplitude Plots:

Summed Amplitude Plots:

Summed Amplitude Plots:

Summed Amplitude Plots:

RBX:

Prepare maps of good/bad channels:

Prepare maps of good/bad channels:

Prepare maps of good/bad channels:

Prepare maps of good/bad channels:


CALO JETS


PF JETS

Definition at line 61 of file WriteMemoryFile.cc.

References qcdUeDQM_cfi::algorithm, dir2webdir::argc, GCPpyPlots::argv, cms::cuda::assert(), cms::cuda::bs, ALCARECOTkAlJpsiMuMu_cff::charge, hltPixelTracks_cff::chi2, cleanSimTrack_minRecHits, cleanSimTrack_minSimHits, mkfit::DataFile::CloseWrite(), mkfit::LayerNumberConverter::convertLayerNumber(), funct::cos(), submitPVResolutionJobs::count, gather_cfg::cout, EE, submitPVResolutionJobs::err, mkfit::DataFile::ES_BeamSpot, mkfit::DataFile::ES_CmsswTracks, mkfit::DataFile::ES_HitIterMasks, mkfit::DataFile::ES_Seeds, edmPickEvents::event, beamvalidation::exit(), f, dqmMemoryStats::float, Glued, hfClusterShapes_cfi::hits, ntupleEnum::hltIter0, trackingPlots::hp, mps_fire::i, cuy::ii, InitialStep_cff::initialStep, InefficientDoubleROC::inputFileName, createfilelist::int, mkfit::TrackBase::InTimePU, Invalid, dqmiolumiharvest::j, label, M_PI_2, SiStripPI::min, mkfit::TrackerInfo::n_layers(), next_arg_option(), next_arg_or_die(), cms::cuda::nh, nHits, mkfit::LayerNumberConverter::nLayers(), EgHLTOffHistBins_cfi::nr, nt, mkfit::DataFile::openWrite(), or, mkfit::TrackBase::OutOfTimePU, makeListRunsInFiles::outputFileName, SiStripOfflineCRack_cfg::p2, mkfit::phase1, mkfit::phase2, Phase2OT, Pixel, funct::pow(), printHelp(), DiDispStaMuonMonitor_cfi::pt, pv::pT, multPhiCorr_741_25nsDY_cfi::px, multPhiCorr_741_25nsDY_cfi::py, mkfit::TrackerInfo::read_bin_file(), mkfit::TrackBase::Signal, funct::sin(), findQualityFiles::size, mathSSE::sqrt(), submitPVResolutionJobs::stderr, AlCaHLTBitMon_QueryRunRegistry::string, Strip, submitPVValidationJobs::t, cmsswSequenceInfo::tp, HLT_2022v15_cff::track, mitigatedMETSequence_cff::U, useMatched, verbosity, and x.

61  {
62  bool haveInput = false;
64  bool haveOutput = false;
66  std::string geoFileName("CMS-phase1.bin");
67 
68  bool cleanSimTracks = false;
69  bool writeAllEvents = false;
70  bool writeRecTracks = false;
71  bool writeHitIterMasks = false;
72  bool applyCCC = false;
73  bool allSeeds = false;
74 
75  int verbosity = 0;
76  long long maxevt = -1;
77 
78  int cutValueCCC = 1620; //Nominal value (from first iteration of CMSSW) is 1620
79 
80  lStr_t mArgs;
81  for (int i = 1; i < argc; ++i) {
82  mArgs.push_back(argv[i]);
83  }
84 
85  lStr_i i = mArgs.begin();
86  while (i != mArgs.end()) {
87  lStr_i start = i;
88 
89  if (*i == "-h" || *i == "-help" || *i == "--help") {
90  printHelp(argv[0]);
91  } else if (*i == "--input") {
92  next_arg_or_die(mArgs, i);
93  inputFileName = *i;
94  haveInput = true;
95  } else if (*i == "--output") {
96  next_arg_or_die(mArgs, i);
97  outputFileName = *i;
98  haveOutput = true;
99  } else if (*i == "--geo") {
100  next_arg_or_die(mArgs, i);
101  geoFileName = *i;
102  } else if (*i == "--verbosity") {
103  next_arg_or_die(mArgs, i);
104  verbosity = std::atoi(i->c_str());
105  } else if (*i == "--maxevt") {
106  next_arg_or_die(mArgs, i);
107  maxevt = std::atoi(i->c_str());
108  } else if (*i == "--clean-sim-tracks") {
109  cleanSimTracks = true;
110  } else if (*i == "--write-all-events") {
111  writeAllEvents = true;
112  } else if (*i == "--write-rec-tracks") {
113  writeRecTracks = true;
114  } else if (*i == "--write-hit-iter-masks") {
115  writeHitIterMasks = true;
116  } else if (*i == "--apply-ccc") {
117  applyCCC = true;
118  if (next_arg_option(mArgs, i)) {
119  cutValueCCC = std::atoi(i->c_str());
120  }
121  } else if (*i == "--all-seeds") {
122  allSeeds = true;
123  } else {
124  fprintf(stderr, "Error: Unknown option/argument '%s'.\n", i->c_str());
125  printHelp(argv[0]);
126  exit(1);
127  }
128  mArgs.erase(start, ++i);
129  } //while arguments
130 
131  if (not haveOutput or not haveInput) {
132  fprintf(stderr, "Error: both input and output are required\n");
133  printHelp(argv[0]);
134  exit(1);
135  }
136 
137  using namespace std;
138 
139  TFile* f = TFile::Open(inputFileName.c_str());
140  if (f == 0) {
141  fprintf(stderr, "Failed opening input root file '%s'\n", inputFileName.c_str());
142  exit(1);
143  }
144 
145  TTree* t = (TTree*)f->Get("trackingNtuple/tree");
146 
147  bool hasPh2hits = t->GetBranch("ph2_isLower") != nullptr;
148  TrackerInfo tkinfo;
149  tkinfo.read_bin_file(geoFileName);
151  const unsigned int nTotalLayers = lnc.nLayers();
152  assert(nTotalLayers == (unsigned int)tkinfo.n_layers());
153 
154  int nstot = 0;
155  std::vector<int> nhitstot(nTotalLayers, 0);
156 
157  unsigned long long event;
158  t->SetBranchAddress("event", &event);
159 
160  //sim tracks
161  std::vector<float>* sim_eta = 0;
162  std::vector<float>* sim_px = 0;
163  std::vector<float>* sim_py = 0;
164  std::vector<float>* sim_pz = 0;
165  std::vector<int>* sim_parentVtxIdx = 0;
166  std::vector<int>* sim_q = 0;
167  std::vector<int>* sim_event = 0;
168  std::vector<int>* sim_bunchCrossing = 0;
169  std::vector<int>* sim_nValid = 0; //simHit count, actually
170  t->SetBranchAddress("sim_eta", &sim_eta);
171  t->SetBranchAddress("sim_px", &sim_px);
172  t->SetBranchAddress("sim_py", &sim_py);
173  t->SetBranchAddress("sim_pz", &sim_pz);
174  t->SetBranchAddress("sim_parentVtxIdx", &sim_parentVtxIdx);
175  t->SetBranchAddress("sim_q", &sim_q);
176  t->SetBranchAddress("sim_event", &sim_event);
177  t->SetBranchAddress("sim_bunchCrossing", &sim_bunchCrossing);
178  t->SetBranchAddress("sim_nValid", &sim_nValid);
179 
180  std::vector<vector<int>>* sim_trkIdx = 0;
181  t->SetBranchAddress("sim_trkIdx", &sim_trkIdx);
182 
183  //simvtx
184  std::vector<float>* simvtx_x = 0;
185  std::vector<float>* simvtx_y = 0;
186  std::vector<float>* simvtx_z = 0;
187  t->SetBranchAddress("simvtx_x", &simvtx_x);
188  t->SetBranchAddress("simvtx_y", &simvtx_y);
189  t->SetBranchAddress("simvtx_z", &simvtx_z);
190 
191  //simhit
192  std::vector<short>* simhit_process = 0;
193  std::vector<int>* simhit_particle = 0;
194  std::vector<int>* simhit_simTrkIdx = 0;
195  std::vector<float>* simhit_x = 0;
196  std::vector<float>* simhit_y = 0;
197  std::vector<float>* simhit_z = 0;
198  std::vector<float>* simhit_px = 0;
199  std::vector<float>* simhit_py = 0;
200  std::vector<float>* simhit_pz = 0;
201  t->SetBranchAddress("simhit_process", &simhit_process);
202  t->SetBranchAddress("simhit_particle", &simhit_particle);
203  t->SetBranchAddress("simhit_simTrkIdx", &simhit_simTrkIdx);
204  t->SetBranchAddress("simhit_x", &simhit_x);
205  t->SetBranchAddress("simhit_y", &simhit_y);
206  t->SetBranchAddress("simhit_z", &simhit_z);
207  t->SetBranchAddress("simhit_px", &simhit_px);
208  t->SetBranchAddress("simhit_py", &simhit_py);
209  t->SetBranchAddress("simhit_pz", &simhit_pz);
210 
211  std::vector<std::vector<int>>* simhit_hitIdx = 0;
212  t->SetBranchAddress("simhit_hitIdx", &simhit_hitIdx);
213  std::vector<std::vector<int>>* simhit_hitType = 0;
214  t->SetBranchAddress("simhit_hitType", &simhit_hitType);
215 
216  //rec tracks
217  std::vector<int>* trk_q = 0;
218  std::vector<unsigned int>* trk_nValid = 0;
219  std::vector<int>* trk_seedIdx = 0;
220  std::vector<unsigned long long>* trk_algoMask = 0;
221  std::vector<unsigned int>* trk_algo = 0;
222  std::vector<unsigned int>* trk_originalAlgo = 0;
223  std::vector<float>* trk_nChi2 = 0;
224  std::vector<float>* trk_px = 0;
225  std::vector<float>* trk_py = 0;
226  std::vector<float>* trk_pz = 0;
227  std::vector<float>* trk_pt = 0;
228  std::vector<float>* trk_phi = 0;
229  std::vector<float>* trk_lambda = 0;
230  std::vector<float>* trk_refpoint_x = 0;
231  std::vector<float>* trk_refpoint_y = 0;
232  std::vector<float>* trk_refpoint_z = 0;
233  std::vector<float>* trk_dxyErr = 0;
234  std::vector<float>* trk_dzErr = 0;
235  std::vector<float>* trk_ptErr = 0;
236  std::vector<float>* trk_phiErr = 0;
237  std::vector<float>* trk_lambdaErr = 0;
238  t->SetBranchAddress("trk_q", &trk_q);
239  t->SetBranchAddress("trk_nValid", &trk_nValid);
240  t->SetBranchAddress("trk_seedIdx", &trk_seedIdx);
241  t->SetBranchAddress("trk_algoMask", &trk_algoMask);
242  t->SetBranchAddress("trk_algo", &trk_algo);
243  t->SetBranchAddress("trk_originalAlgo", &trk_originalAlgo);
244  t->SetBranchAddress("trk_nChi2", &trk_nChi2);
245  t->SetBranchAddress("trk_px", &trk_px);
246  t->SetBranchAddress("trk_py", &trk_py);
247  t->SetBranchAddress("trk_pz", &trk_pz);
248  t->SetBranchAddress("trk_pt", &trk_pt);
249  t->SetBranchAddress("trk_phi", &trk_phi);
250  t->SetBranchAddress("trk_lambda", &trk_lambda);
251  t->SetBranchAddress("trk_refpoint_x", &trk_refpoint_x);
252  t->SetBranchAddress("trk_refpoint_y", &trk_refpoint_y);
253  t->SetBranchAddress("trk_refpoint_z", &trk_refpoint_z);
254  t->SetBranchAddress("trk_dxyErr", &trk_dxyErr);
255  t->SetBranchAddress("trk_dzErr", &trk_dzErr);
256  t->SetBranchAddress("trk_ptErr", &trk_ptErr);
257  t->SetBranchAddress("trk_phiErr", &trk_phiErr);
258  t->SetBranchAddress("trk_lambdaErr", &trk_lambdaErr);
259 
260  std::vector<std::vector<int>>* trk_hitIdx = 0;
261  t->SetBranchAddress("trk_hitIdx", &trk_hitIdx);
262  std::vector<std::vector<int>>* trk_hitType = 0;
263  t->SetBranchAddress("trk_hitType", &trk_hitType);
264 
265  //seeds
266  std::vector<float>* see_stateTrajGlbX = 0;
267  std::vector<float>* see_stateTrajGlbY = 0;
268  std::vector<float>* see_stateTrajGlbZ = 0;
269  std::vector<float>* see_stateTrajGlbPx = 0;
270  std::vector<float>* see_stateTrajGlbPy = 0;
271  std::vector<float>* see_stateTrajGlbPz = 0;
272  std::vector<float>* see_eta = 0; //PCA parameters
273  std::vector<float>* see_pt = 0; //PCA parameters
274  std::vector<float>* see_stateCcov00 = 0;
275  std::vector<float>* see_stateCcov01 = 0;
276  std::vector<float>* see_stateCcov02 = 0;
277  std::vector<float>* see_stateCcov03 = 0;
278  std::vector<float>* see_stateCcov04 = 0;
279  std::vector<float>* see_stateCcov05 = 0;
280  std::vector<float>* see_stateCcov11 = 0;
281  std::vector<float>* see_stateCcov12 = 0;
282  std::vector<float>* see_stateCcov13 = 0;
283  std::vector<float>* see_stateCcov14 = 0;
284  std::vector<float>* see_stateCcov15 = 0;
285  std::vector<float>* see_stateCcov22 = 0;
286  std::vector<float>* see_stateCcov23 = 0;
287  std::vector<float>* see_stateCcov24 = 0;
288  std::vector<float>* see_stateCcov25 = 0;
289  std::vector<float>* see_stateCcov33 = 0;
290  std::vector<float>* see_stateCcov34 = 0;
291  std::vector<float>* see_stateCcov35 = 0;
292  std::vector<float>* see_stateCcov44 = 0;
293  std::vector<float>* see_stateCcov45 = 0;
294  std::vector<float>* see_stateCcov55 = 0;
295  std::vector<std::vector<float>>* see_stateCurvCov = 0;
296  std::vector<int>* see_q = 0;
297  std::vector<unsigned int>* see_algo = 0;
298  t->SetBranchAddress("see_stateTrajGlbX", &see_stateTrajGlbX);
299  t->SetBranchAddress("see_stateTrajGlbY", &see_stateTrajGlbY);
300  t->SetBranchAddress("see_stateTrajGlbZ", &see_stateTrajGlbZ);
301  t->SetBranchAddress("see_stateTrajGlbPx", &see_stateTrajGlbPx);
302  t->SetBranchAddress("see_stateTrajGlbPy", &see_stateTrajGlbPy);
303  t->SetBranchAddress("see_stateTrajGlbPz", &see_stateTrajGlbPz);
304  t->SetBranchAddress("see_eta", &see_eta);
305  t->SetBranchAddress("see_pt", &see_pt);
306 
307  bool hasCartCov = t->GetBranch("see_stateCcov00") != nullptr;
308  if (hasCartCov) {
309  t->SetBranchAddress("see_stateCcov00", &see_stateCcov00);
310  t->SetBranchAddress("see_stateCcov01", &see_stateCcov01);
311  t->SetBranchAddress("see_stateCcov02", &see_stateCcov02);
312  t->SetBranchAddress("see_stateCcov03", &see_stateCcov03);
313  t->SetBranchAddress("see_stateCcov04", &see_stateCcov04);
314  t->SetBranchAddress("see_stateCcov05", &see_stateCcov05);
315  t->SetBranchAddress("see_stateCcov11", &see_stateCcov11);
316  t->SetBranchAddress("see_stateCcov12", &see_stateCcov12);
317  t->SetBranchAddress("see_stateCcov13", &see_stateCcov13);
318  t->SetBranchAddress("see_stateCcov14", &see_stateCcov14);
319  t->SetBranchAddress("see_stateCcov15", &see_stateCcov15);
320  t->SetBranchAddress("see_stateCcov22", &see_stateCcov22);
321  t->SetBranchAddress("see_stateCcov23", &see_stateCcov23);
322  t->SetBranchAddress("see_stateCcov24", &see_stateCcov24);
323  t->SetBranchAddress("see_stateCcov25", &see_stateCcov25);
324  t->SetBranchAddress("see_stateCcov33", &see_stateCcov33);
325  t->SetBranchAddress("see_stateCcov34", &see_stateCcov34);
326  t->SetBranchAddress("see_stateCcov35", &see_stateCcov35);
327  t->SetBranchAddress("see_stateCcov44", &see_stateCcov44);
328  t->SetBranchAddress("see_stateCcov45", &see_stateCcov45);
329  t->SetBranchAddress("see_stateCcov55", &see_stateCcov55);
330  } else {
331  t->SetBranchAddress("see_stateCurvCov", &see_stateCurvCov);
332  }
333  t->SetBranchAddress("see_q", &see_q);
334  t->SetBranchAddress("see_algo", &see_algo);
335 
336  std::vector<std::vector<int>>* see_hitIdx = 0;
337  t->SetBranchAddress("see_hitIdx", &see_hitIdx);
338  std::vector<std::vector<int>>* see_hitType = 0;
339  t->SetBranchAddress("see_hitType", &see_hitType);
340 
341  //pixel hits
342  vector<unsigned short>* pix_det = 0;
343  vector<unsigned short>* pix_lay = 0;
344  vector<unsigned int>* pix_detId = 0;
345  vector<float>* pix_x = 0;
346  vector<float>* pix_y = 0;
347  vector<float>* pix_z = 0;
348  vector<float>* pix_xx = 0;
349  vector<float>* pix_xy = 0;
350  vector<float>* pix_yy = 0;
351  vector<float>* pix_yz = 0;
352  vector<float>* pix_zz = 0;
353  vector<float>* pix_zx = 0;
354  vector<int>* pix_csize_col = 0;
355  vector<int>* pix_csize_row = 0;
356  vector<uint64_t>* pix_usedMask = 0;
357  //these were renamed in CMSSW_9_1_0: auto-detect
358  bool has910_det_lay = t->GetBranch("pix_det") == nullptr;
359  if (has910_det_lay) {
360  t->SetBranchAddress("pix_subdet", &pix_det);
361  t->SetBranchAddress("pix_layer", &pix_lay);
362  } else {
363  t->SetBranchAddress("pix_det", &pix_det);
364  t->SetBranchAddress("pix_lay", &pix_lay);
365  }
366  t->SetBranchAddress("pix_detId", &pix_detId);
367  t->SetBranchAddress("pix_x", &pix_x);
368  t->SetBranchAddress("pix_y", &pix_y);
369  t->SetBranchAddress("pix_z", &pix_z);
370  t->SetBranchAddress("pix_xx", &pix_xx);
371  t->SetBranchAddress("pix_xy", &pix_xy);
372  t->SetBranchAddress("pix_yy", &pix_yy);
373  t->SetBranchAddress("pix_yz", &pix_yz);
374  t->SetBranchAddress("pix_zz", &pix_zz);
375  t->SetBranchAddress("pix_zx", &pix_zx);
376  t->SetBranchAddress("pix_clustSizeCol", &pix_csize_col);
377  t->SetBranchAddress("pix_clustSizeRow", &pix_csize_row);
378  if (writeHitIterMasks) {
379  t->SetBranchAddress("pix_usedMask", &pix_usedMask);
380  }
381 
382  vector<vector<int>>* pix_simHitIdx = 0;
383  t->SetBranchAddress("pix_simHitIdx", &pix_simHitIdx);
384  vector<vector<float>>* pix_chargeFraction = 0;
385  t->SetBranchAddress("pix_chargeFraction", &pix_chargeFraction);
386 
387  //strip hits
388  vector<short>* glu_isBarrel = 0;
389  vector<unsigned int>* glu_det = 0;
390  vector<unsigned int>* glu_lay = 0;
391  vector<unsigned int>* glu_detId = 0;
392  vector<int>* glu_monoIdx = 0;
393  vector<int>* glu_stereoIdx = 0;
394  vector<float>* glu_x = 0;
395  vector<float>* glu_y = 0;
396  vector<float>* glu_z = 0;
397  vector<float>* glu_xx = 0;
398  vector<float>* glu_xy = 0;
399  vector<float>* glu_yy = 0;
400  vector<float>* glu_yz = 0;
401  vector<float>* glu_zz = 0;
402  vector<float>* glu_zx = 0;
403  if (!hasPh2hits) {
404  t->SetBranchAddress("glu_isBarrel", &glu_isBarrel);
405  if (has910_det_lay) {
406  t->SetBranchAddress("glu_subdet", &glu_det);
407  t->SetBranchAddress("glu_layer", &glu_lay);
408  } else {
409  t->SetBranchAddress("glu_det", &glu_det);
410  t->SetBranchAddress("glu_lay", &glu_lay);
411  }
412  t->SetBranchAddress("glu_detId", &glu_detId);
413  t->SetBranchAddress("glu_monoIdx", &glu_monoIdx);
414  t->SetBranchAddress("glu_stereoIdx", &glu_stereoIdx);
415  t->SetBranchAddress("glu_x", &glu_x);
416  t->SetBranchAddress("glu_y", &glu_y);
417  t->SetBranchAddress("glu_z", &glu_z);
418  t->SetBranchAddress("glu_xx", &glu_xx);
419  t->SetBranchAddress("glu_xy", &glu_xy);
420  t->SetBranchAddress("glu_yy", &glu_yy);
421  t->SetBranchAddress("glu_yz", &glu_yz);
422  t->SetBranchAddress("glu_zz", &glu_zz);
423  t->SetBranchAddress("glu_zx", &glu_zx);
424  }
425 
426  vector<short>* str_isBarrel = 0;
427  vector<short>* str_isStereo = 0;
428  vector<unsigned int>* str_det = 0;
429  vector<unsigned int>* str_lay = 0;
430  vector<unsigned int>* str_detId = 0;
431  vector<unsigned int>* str_simType = 0;
432  vector<float>* str_x = 0;
433  vector<float>* str_y = 0;
434  vector<float>* str_z = 0;
435  vector<float>* str_xx = 0;
436  vector<float>* str_xy = 0;
437  vector<float>* str_yy = 0;
438  vector<float>* str_yz = 0;
439  vector<float>* str_zz = 0;
440  vector<float>* str_zx = 0;
441  vector<float>* str_chargePerCM = 0;
442  vector<int>* str_csize = 0;
443  vector<uint64_t>* str_usedMask = 0;
444  vector<vector<int>>* str_simHitIdx = 0;
445  vector<vector<float>>* str_chargeFraction = 0;
446  if (!hasPh2hits) {
447  t->SetBranchAddress("str_isBarrel", &str_isBarrel);
448  t->SetBranchAddress("str_isStereo", &str_isStereo);
449  if (has910_det_lay) {
450  t->SetBranchAddress("str_subdet", &str_det);
451  t->SetBranchAddress("str_layer", &str_lay);
452  } else {
453  t->SetBranchAddress("str_det", &str_det);
454  t->SetBranchAddress("str_lay", &str_lay);
455  }
456  t->SetBranchAddress("str_detId", &str_detId);
457  t->SetBranchAddress("str_simType", &str_simType);
458  t->SetBranchAddress("str_x", &str_x);
459  t->SetBranchAddress("str_y", &str_y);
460  t->SetBranchAddress("str_z", &str_z);
461  t->SetBranchAddress("str_xx", &str_xx);
462  t->SetBranchAddress("str_xy", &str_xy);
463  t->SetBranchAddress("str_yy", &str_yy);
464  t->SetBranchAddress("str_yz", &str_yz);
465  t->SetBranchAddress("str_zz", &str_zz);
466  t->SetBranchAddress("str_zx", &str_zx);
467  t->SetBranchAddress("str_chargePerCM", &str_chargePerCM);
468  t->SetBranchAddress("str_clustSize", &str_csize);
469  if (writeHitIterMasks) {
470  t->SetBranchAddress("str_usedMask", &str_usedMask);
471  }
472 
473  t->SetBranchAddress("str_simHitIdx", &str_simHitIdx);
474  t->SetBranchAddress("str_chargeFraction", &str_chargeFraction);
475  }
476 
477  vector<unsigned short>* ph2_isLower = 0;
478  vector<unsigned short>* ph2_subdet = 0;
479  vector<unsigned short>* ph2_layer = 0;
480  vector<unsigned int>* ph2_detId = 0;
481  vector<unsigned short>* ph2_simType = 0;
482  vector<float>* ph2_x = 0;
483  vector<float>* ph2_y = 0;
484  vector<float>* ph2_z = 0;
485  vector<float>* ph2_xx = 0;
486  vector<float>* ph2_xy = 0;
487  vector<float>* ph2_yy = 0;
488  vector<float>* ph2_yz = 0;
489  vector<float>* ph2_zz = 0;
490  vector<float>* ph2_zx = 0;
491  vector<uint64_t>* ph2_usedMask = 0;
492  vector<vector<int>>* ph2_simHitIdx = 0;
493  if (hasPh2hits && applyCCC)
494  std::cout << "WARNING: applyCCC is set for Phase2 inputs: applyCCC will be ignored" << std::endl;
495  if (hasPh2hits) {
496  t->SetBranchAddress("ph2_isLower", &ph2_isLower);
497  t->SetBranchAddress("ph2_subdet", &ph2_subdet);
498  t->SetBranchAddress("ph2_layer", &ph2_layer);
499  t->SetBranchAddress("ph2_detId", &ph2_detId);
500  t->SetBranchAddress("ph2_simType", &ph2_simType);
501  t->SetBranchAddress("ph2_x", &ph2_x);
502  t->SetBranchAddress("ph2_y", &ph2_y);
503  t->SetBranchAddress("ph2_z", &ph2_z);
504  t->SetBranchAddress("ph2_xx", &ph2_xx);
505  t->SetBranchAddress("ph2_xy", &ph2_xy);
506  t->SetBranchAddress("ph2_yy", &ph2_yy);
507  t->SetBranchAddress("ph2_yz", &ph2_yz);
508  t->SetBranchAddress("ph2_zz", &ph2_zz);
509  t->SetBranchAddress("ph2_zx", &ph2_zx);
510  if (writeHitIterMasks) {
511  t->SetBranchAddress("ph2_usedMask", &ph2_usedMask);
512  }
513  t->SetBranchAddress("ph2_simHitIdx", &ph2_simHitIdx);
514  }
515  vector<float> ph2_chargeFraction_dummy(16, 0.f);
516 
517  // beam spot
518  float bsp_x;
519  float bsp_y;
520  float bsp_z;
521  float bsp_sigmax;
522  float bsp_sigmay;
523  float bsp_sigmaz;
524  t->SetBranchAddress("bsp_x", &bsp_x);
525  t->SetBranchAddress("bsp_y", &bsp_y);
526  t->SetBranchAddress("bsp_z", &bsp_z);
527  t->SetBranchAddress("bsp_sigmax", &bsp_sigmax);
528  t->SetBranchAddress("bsp_sigmay", &bsp_sigmay);
529  t->SetBranchAddress("bsp_sigmaz", &bsp_sigmaz);
530 
531  long long totentries = t->GetEntries();
532  long long savedEvents = 0;
533 
534  DataFile data_file;
535  int outOptions = DataFile::ES_Seeds;
536  if (writeRecTracks)
537  outOptions |= DataFile::ES_CmsswTracks;
538  if (writeHitIterMasks)
539  outOptions |= DataFile::ES_HitIterMasks;
540  outOptions |= DataFile::ES_BeamSpot;
541 
542  if (maxevt < 0)
543  maxevt = totentries;
544  data_file.openWrite(outputFileName, static_cast<int>(nTotalLayers), std::min(maxevt, totentries), outOptions);
545 
546  Event EE(0, static_cast<int>(nTotalLayers));
547 
548  int numFailCCC = 0;
549  int numTotalStr = 0;
550  // gDebug = 8;
551 
552  for (long long i = 0; savedEvents < maxevt && i < totentries && i < maxevt; ++i) {
553  EE.reset(i);
554 
555  cout << "process entry i=" << i << " out of " << totentries << ", saved so far " << savedEvents
556  << ", with max=" << maxevt << endl;
557 
558  t->GetEntry(i);
559 
560  cout << "edm event=" << event << endl;
561 
562  auto& bs = EE.beamSpot_;
563  bs.x = bsp_x;
564  bs.y = bsp_y;
565  bs.z = bsp_z;
566  bs.sigmaZ = bsp_sigmaz;
567  bs.beamWidthX = bsp_sigmax;
568  bs.beamWidthY = bsp_sigmay;
569  //dxdz and dydz are not in the trackingNtuple at the moment
570 
571  if (!hasPh2hits) {
572  for (unsigned int istr = 0; istr < str_lay->size(); ++istr) {
573  if (str_chargePerCM->at(istr) < cutValueCCC)
574  numFailCCC++;
575  numTotalStr++;
576  }
577  }
578 
579  auto nSims = sim_q->size();
580  if (nSims == 0) {
581  cout << "branches not loaded" << endl;
582  exit(1);
583  }
584  if (verbosity > 0)
585  std::cout << __FILE__ << " " << __LINE__ << " nSims " << nSims << " nSeeds " << see_q->size() << " nRecT "
586  << trk_q->size() << std::endl;
587 
588  //find best matching tkIdx from a list of simhits indices
589  auto bestTkIdx = [&](std::vector<int> const& shs, std::vector<float> const& shfs, int rhIdx, HitType rhType) {
590  //assume that all simhits are associated
591  int ibest = -1;
592  int shbest = -1;
593  float hpbest = -1;
594  float tpbest = -1;
595  float hfbest = -1;
596 
597  float maxfrac = -1;
598  int ish = -1;
599  int nshs = shs.size();
600  for (auto const sh : shs) {
601  ish++;
602  auto tkidx = simhit_simTrkIdx->at(sh);
603  //use only sh with available TP
604  if (tkidx < 0)
605  continue;
606 
607  auto hpx = simhit_px->at(sh);
608  auto hpy = simhit_py->at(sh);
609  auto hpz = simhit_pz->at(sh);
610  auto hp = sqrt(hpx * hpx + hpy * hpy + hpz * hpz);
611 
612  //look only at hits with p> 50 MeV
613  if (hp < 0.05f)
614  continue;
615 
616  auto tpx = sim_px->at(tkidx);
617  auto tpy = sim_py->at(tkidx);
618  auto tpz = sim_pz->at(tkidx);
619  auto tp = sqrt(tpx * tpx + tpy * tpy + tpz * tpz);
620 
621  //take only hits with hp> 0.5*tp
622  if (hp < 0.5 * tp)
623  continue;
624 
625  //pick tkidx corresponding to max hp/tp; .. this is probably redundant
626  if (maxfrac < hp / tp) {
627  maxfrac = hp / tp;
628  ibest = tkidx;
629  shbest = sh;
630  hpbest = hp;
631  tpbest = tp;
632  hfbest = shfs[ish];
633  }
634  }
635 
636  //arbitration: a rechit with one matching sim is matched to sim if it's the first
637  //FIXME: SOME BETTER SELECTION CAN BE DONE (it will require some more correlated knowledge)
638  if (nshs == 1 && ibest >= 0) {
639  auto const& srhIdxV = simhit_hitIdx->at(shbest);
640  auto const& srhTypeV = simhit_hitType->at(shbest);
641  int ih = -1;
642  for (auto itype : srhTypeV) {
643  ih++;
644  if (HitType(itype) == rhType && srhIdxV[ih] != rhIdx) {
645  ibest = -1;
646  break;
647  }
648  }
649  }
650 
651  if (ibest >= 0 && false) {
652  std::cout << " best tkIdx " << ibest << " rh " << rhIdx << " for sh " << shbest << " out of " << shs.size()
653  << " hp " << hpbest << " chF " << hfbest << " tp " << tpbest << " process "
654  << simhit_process->at(shbest) << " particle " << simhit_particle->at(shbest) << std::endl;
655  if (rhType == HitType::Strip) {
656  std::cout << " sh " << simhit_x->at(shbest) << ", " << simhit_y->at(shbest) << ", " << simhit_z->at(shbest)
657  << " rh " << str_x->at(rhIdx) << ", " << str_y->at(rhIdx) << ", " << str_z->at(rhIdx) << std::endl;
658  }
659  }
660  return ibest;
661  };
662 
663  vector<Track>& simTracks_ = EE.simTracks_;
664  vector<int> simTrackIdx_(sim_q->size(), -1); //keep track of original index in ntuple
665  vector<int> seedSimIdx(see_q->size(), -1);
666  for (unsigned int isim = 0; isim < sim_q->size(); ++isim) {
667  //load sim production vertex data
668  auto iVtx = sim_parentVtxIdx->at(isim);
669  constexpr float largeValF = 9999.f;
670  float sim_prodx = iVtx >= 0 ? simvtx_x->at(iVtx) : largeValF;
671  float sim_prody = iVtx >= 0 ? simvtx_y->at(iVtx) : largeValF;
672  float sim_prodz = iVtx >= 0 ? simvtx_z->at(iVtx) : largeValF;
673  //if (fabs(sim_eta->at(isim))>0.8) continue;
674 
675  vector<int> const& trkIdxV = sim_trkIdx->at(isim);
676 
677  //if (trkIdx<0) continue;
678  //FIXME: CHECK IF THE LOOP AND BEST SELECTION IS NEEDED.
679  //Pick the first
680  const int trkIdx = trkIdxV.empty() ? -1 : trkIdxV[0];
681 
682  int nlay = 0;
683  if (trkIdx >= 0) {
684  std::vector<int> hitlay(nTotalLayers, 0);
685  auto const& hits = trk_hitIdx->at(trkIdx);
686  auto const& hitTypes = trk_hitType->at(trkIdx);
687  auto nHits = hits.size();
688  for (auto ihit = 0U; ihit < nHits; ++ihit) {
689  auto ihIdx = hits[ihit];
690  auto const ihType = HitType(hitTypes[ihit]);
691 
692  switch (ihType) {
693  case HitType::Pixel: {
694  int ipix = ihIdx;
695  if (ipix < 0)
696  continue;
697  int cmsswlay =
698  lnc.convertLayerNumber(pix_det->at(ipix), pix_lay->at(ipix), useMatched, -1, pix_z->at(ipix) > 0);
699  if (cmsswlay >= 0 && cmsswlay < static_cast<int>(nTotalLayers))
700  hitlay[cmsswlay]++;
701  break;
702  }
703  case HitType::Strip: {
704  int istr = ihIdx;
705  if (istr < 0)
706  continue;
707  int cmsswlay = lnc.convertLayerNumber(
708  str_det->at(istr), str_lay->at(istr), useMatched, str_isStereo->at(istr), str_z->at(istr) > 0);
709  if (cmsswlay >= 0 && cmsswlay < static_cast<int>(nTotalLayers))
710  hitlay[cmsswlay]++;
711  break;
712  }
713  case HitType::Glued: {
714  if (useMatched) {
715  int iglu = ihIdx;
716  if (iglu < 0)
717  continue;
718  int cmsswlay =
719  lnc.convertLayerNumber(glu_det->at(iglu), glu_lay->at(iglu), useMatched, -1, glu_z->at(iglu) > 0);
720  if (cmsswlay >= 0 && cmsswlay < static_cast<int>(nTotalLayers))
721  hitlay[cmsswlay]++;
722  }
723  break;
724  }
725  case HitType::Phase2OT: {
726  int istr = ihIdx;
727  if (istr < 0)
728  continue;
729  int cmsswlay = lnc.convertLayerNumber(
730  ph2_subdet->at(istr), ph2_layer->at(istr), useMatched, ph2_isLower->at(istr), ph2_z->at(istr) > 0);
731  if (cmsswlay >= 0 && cmsswlay < static_cast<int>(nTotalLayers))
732  hitlay[cmsswlay]++;
733  break;
734  }
735  case HitType::Invalid:
736  break; //FIXME. Skip, really?
737  default:
738  throw std::logic_error("Track type can not be handled");
739  } //hit type
740  } //hits on track
741  for (unsigned int i = 0; i < nTotalLayers; i++)
742  if (hitlay[i] > 0)
743  nlay++;
744  } //count nlay layers on matching reco track
745 
746  //cout << Form("track q=%2i p=(%6.3f, %6.3f, %6.3f) x=(%6.3f, %6.3f, %6.3f) nlay=%i",sim_q->at(isim),sim_px->at(isim),sim_py->at(isim),sim_pz->at(isim),sim_prodx,sim_prody,sim_prodz,nlay) << endl;
747 
748  SVector3 pos(sim_prodx, sim_prody, sim_prodz);
749  SVector3 mom(sim_px->at(isim), sim_py->at(isim), sim_pz->at(isim));
751  err.At(0, 0) = sim_prodx * sim_prodx;
752  err.At(1, 1) = sim_prody * sim_prody;
753  err.At(2, 2) = sim_prodz * sim_prodz;
754  err.At(3, 3) = sim_px->at(isim) * sim_px->at(isim);
755  err.At(4, 4) = sim_py->at(isim) * sim_py->at(isim);
756  err.At(5, 5) = sim_pz->at(isim) * sim_pz->at(isim);
757  TrackState state(sim_q->at(isim), pos, mom, err);
758  state.convertFromCartesianToCCS();
759  //create track: store number of reco hits in place of track chi2; fill hits later
760  // set label to be its own index in the output file
761  Track track(state, float(nlay), simTracks_.size(), 0, nullptr);
762  if (sim_bunchCrossing->at(isim) == 0) { //in time
763  if (sim_event->at(isim) == 0)
764  track.setProdType(Track::ProdType::Signal);
765  else
766  track.setProdType(Track::ProdType::InTimePU);
767  } else {
768  track.setProdType(Track::ProdType::OutOfTimePU);
769  }
770  if (trkIdx >= 0) {
771  int seedIdx = trk_seedIdx->at(trkIdx);
772  // Unused: auto const& shTypes = see_hitType->at(seedIdx);
773  seedSimIdx[seedIdx] = simTracks_.size();
774  }
775  if (cleanSimTracks) {
776  if (sim_nValid->at(isim) < cleanSimTrack_minSimHits)
777  continue;
778  if (cleanSimTrack_minRecHits > 0) {
779  int nRecToSimHit = 0;
780  for (unsigned int ipix = 0; ipix < pix_lay->size() && nRecToSimHit < cleanSimTrack_minRecHits; ++ipix) {
781  int ilay = -1;
782  ilay = lnc.convertLayerNumber(pix_det->at(ipix), pix_lay->at(ipix), useMatched, -1, pix_z->at(ipix) > 0);
783  if (ilay < 0)
784  continue;
785  int simTkIdxNt = bestTkIdx(pix_simHitIdx->at(ipix), pix_chargeFraction->at(ipix), ipix, HitType::Pixel);
786  if (simTkIdxNt >= 0)
787  nRecToSimHit++;
788  }
789  if (hasPh2hits) {
790  for (unsigned int istr = 0; istr < ph2_layer->size() && nRecToSimHit < cleanSimTrack_minRecHits; ++istr) {
791  int ilay = -1;
792  ilay = lnc.convertLayerNumber(
793  ph2_subdet->at(istr), ph2_layer->at(istr), useMatched, ph2_isLower->at(istr), ph2_z->at(istr) > 0);
794  if (useMatched && !ph2_isLower->at(istr))
795  continue;
796  if (ilay == -1)
797  continue;
798  int simTkIdxNt = bestTkIdx(ph2_simHitIdx->at(istr), ph2_chargeFraction_dummy, istr, HitType::Phase2OT);
799  if (simTkIdxNt >= 0)
800  nRecToSimHit++;
801  }
802  } else {
803  if (useMatched) {
804  for (unsigned int iglu = 0; iglu < glu_lay->size() && nRecToSimHit < cleanSimTrack_minRecHits; ++iglu) {
805  if (glu_isBarrel->at(iglu) == 0)
806  continue;
807  int igluMono = glu_monoIdx->at(iglu);
808  int simTkIdxNt =
809  bestTkIdx(str_simHitIdx->at(igluMono), str_chargeFraction->at(igluMono), igluMono, HitType::Strip);
810  if (simTkIdxNt >= 0)
811  nRecToSimHit++;
812  }
813  }
814  for (unsigned int istr = 0; istr < str_lay->size() && nRecToSimHit < cleanSimTrack_minRecHits; ++istr) {
815  int ilay = -1;
816  ilay = lnc.convertLayerNumber(
817  str_det->at(istr), str_lay->at(istr), useMatched, str_isStereo->at(istr), str_z->at(istr) > 0);
818  if (useMatched && str_isBarrel->at(istr) == 1 && str_isStereo->at(istr))
819  continue;
820  if (ilay == -1)
821  continue;
822  int simTkIdxNt = bestTkIdx(str_simHitIdx->at(istr), str_chargeFraction->at(istr), istr, HitType::Strip);
823  if (simTkIdxNt >= 0)
824  nRecToSimHit++;
825  }
826  }
827  if (nRecToSimHit < cleanSimTrack_minRecHits)
828  continue;
829  } //count rec-to-sim hits
830  } //cleanSimTracks
831 
832  simTrackIdx_[isim] = simTracks_.size();
833  simTracks_.push_back(track);
834  }
835 
836  if (simTracks_.empty() and not writeAllEvents)
837  continue;
838 
839  vector<Track>& seedTracks_ = EE.seedTracks_;
840  vector<vector<int>> pixHitSeedIdx(pix_lay->size());
841  vector<vector<int>> strHitSeedIdx(hasPh2hits ? 0 : str_lay->size());
842  vector<vector<int>> gluHitSeedIdx(hasPh2hits ? 0 : glu_lay->size());
843  vector<vector<int>> ph2HitSeedIdx(hasPh2hits ? ph2_layer->size() : 0);
844  for (unsigned int is = 0; is < see_q->size(); ++is) {
845  auto isAlgo = TrackAlgorithm(see_algo->at(is));
846  if (not allSeeds)
847  if (isAlgo != TrackAlgorithm::initialStep && isAlgo != TrackAlgorithm::hltIter0)
848  continue; //select seed in acceptance
849  //if (see_pt->at(is)<0.5 || fabs(see_eta->at(is))>0.8) continue;//select seed in acceptance
850  SVector3 pos = SVector3(see_stateTrajGlbX->at(is), see_stateTrajGlbY->at(is), see_stateTrajGlbZ->at(is));
851  SVector3 mom = SVector3(see_stateTrajGlbPx->at(is), see_stateTrajGlbPy->at(is), see_stateTrajGlbPz->at(is));
853  if (hasCartCov) {
854  err.At(0, 0) = see_stateCcov00->at(is);
855  err.At(0, 1) = see_stateCcov01->at(is);
856  err.At(0, 2) = see_stateCcov02->at(is);
857  err.At(0, 3) = see_stateCcov03->at(is);
858  err.At(0, 4) = see_stateCcov04->at(is);
859  err.At(0, 5) = see_stateCcov05->at(is);
860  err.At(1, 1) = see_stateCcov11->at(is);
861  err.At(1, 2) = see_stateCcov12->at(is);
862  err.At(1, 3) = see_stateCcov13->at(is);
863  err.At(1, 4) = see_stateCcov14->at(is);
864  err.At(1, 5) = see_stateCcov15->at(is);
865  err.At(2, 2) = see_stateCcov22->at(is);
866  err.At(2, 3) = see_stateCcov23->at(is);
867  err.At(2, 4) = see_stateCcov24->at(is);
868  err.At(2, 5) = see_stateCcov25->at(is);
869  err.At(3, 3) = see_stateCcov33->at(is);
870  err.At(3, 4) = see_stateCcov34->at(is);
871  err.At(3, 5) = see_stateCcov35->at(is);
872  err.At(4, 4) = see_stateCcov44->at(is);
873  err.At(4, 5) = see_stateCcov45->at(is);
874  err.At(5, 5) = see_stateCcov55->at(is);
875  } else {
876  auto const& vCov = see_stateCurvCov->at(is);
877  assert(vCov.size() == 15);
878  auto vCovP = vCov.begin();
879  for (int i = 0; i < 5; ++i)
880  for (int j = 0; j <= i; ++j)
881  err.At(i, j) = *(vCovP++);
882  }
883  TrackState state(see_q->at(is), pos, mom, err);
884  if (hasCartCov)
885  state.convertFromCartesianToCCS();
886  else
887  state.convertFromGlbCurvilinearToCCS();
888  Track track(state, 0, seedSimIdx[is], 0, nullptr);
889  track.setAlgorithm(isAlgo);
890  auto const& shTypes = see_hitType->at(is);
891  auto const& shIdxs = see_hitIdx->at(is);
892  if (not allSeeds)
893  if (!((isAlgo == TrackAlgorithm::initialStep || isAlgo == TrackAlgorithm::hltIter0) &&
894  std::count(shTypes.begin(), shTypes.end(), int(HitType::Pixel)) >= 3))
895  continue; //check algo and nhits
896  for (unsigned int ip = 0; ip < shTypes.size(); ip++) {
897  unsigned int hidx = shIdxs[ip];
898  switch (HitType(shTypes[ip])) {
899  case HitType::Pixel: {
900  pixHitSeedIdx[hidx].push_back(seedTracks_.size());
901  break;
902  }
903  case HitType::Strip: {
904  strHitSeedIdx[hidx].push_back(seedTracks_.size());
905  break;
906  }
907  case HitType::Glued: {
908  if (not useMatched) {
909  //decompose
910  int uidx = glu_monoIdx->at(hidx);
911  strHitSeedIdx[uidx].push_back(seedTracks_.size());
912  uidx = glu_stereoIdx->at(hidx);
913  strHitSeedIdx[uidx].push_back(seedTracks_.size());
914  } else {
915  gluHitSeedIdx[hidx].push_back(seedTracks_.size());
916  }
917  break;
918  }
919  case HitType::Phase2OT: {
920  ph2HitSeedIdx[hidx].push_back(seedTracks_.size());
921  break;
922  }
923  case HitType::Invalid:
924  break; //FIXME. Skip, really?
925  default:
926  throw std::logic_error("Track hit type can not be handled");
927  } //switch( HitType
928  }
929  seedTracks_.push_back(track);
930  }
931 
932  if (seedTracks_.empty() and not writeAllEvents)
933  continue;
934 
935  vector<Track>& cmsswTracks_ = EE.cmsswTracks_;
936  vector<vector<int>> pixHitRecIdx(pix_lay->size());
937  vector<vector<int>> strHitRecIdx(hasPh2hits ? 0 : str_lay->size());
938  vector<vector<int>> gluHitRecIdx(hasPh2hits ? 0 : glu_lay->size());
939  vector<vector<int>> ph2HitRecIdx(hasPh2hits ? ph2_layer->size() : 0);
940  for (unsigned int ir = 0; ir < trk_q->size(); ++ir) {
941  //check the origin; redundant for initialStep ntuples
942  if (not allSeeds)
943  if ((trk_algoMask->at(ir) & ((1 << int(TrackAlgorithm::initialStep)) | (1 << int(TrackAlgorithm::hltIter0)))) ==
944  0) {
945  if (verbosity > 1) {
946  std::cout << "track " << ir << " failed algo selection for " << int(TrackAlgorithm::initialStep)
947  << ": mask " << trk_algoMask->at(ir) << " origAlgo " << trk_originalAlgo->at(ir) << " algo "
948  << trk_algo->at(ir) << std::endl;
949  }
950  continue;
951  }
952  //fill the state in CCS upfront
954  /*
955  vx = -dxy*sin(phi) - pt*cos(phi)/p*pz/p*dz;
956  vy = dxy*cos(phi) - pt*sin(phi)/p*pz/p*dz;
957  vz = dz*pt*pt/p/p;
958  //partial: ignores cross-terms
959  c(vx,vx) = c(dxy,dxy)*sin(phi)*sin(phi) + c(dz,dz)*pow(pt*cos(phi)/p*pz/p ,2);
960  c(vx,vy) = -c(dxy,dxy)*cos(phi)*sin(phi) + c(dz,dz)*cos(phi)*sin(phi)*pow(pt/p*pz/p, 2);
961  c(vy,vy) = c(dxy,dxy)*cos(phi)*cos(phi) + c(dz,dz)*pow(pt*sin(phi)/p*pz/p ,2);
962  c(vx,vz) = -c(dz,dz)*pt*pt/p/p*pt/p*pz/p*cos(phi);
963  c(vy,vz) = -c(dz,dz)*pt*pt/p/p*pt/p*pz/p*sin(phi);
964  c(vz,vz) = c(dz,dz)*pow(pt*pt/p/p, 2);
965  */
966  float pt = trk_pt->at(ir);
967  float pz = trk_pz->at(ir);
968  float p2 = pt * pt + pz * pz;
969  float phi = trk_phi->at(ir);
970  float sP = sin(phi);
971  float cP = cos(phi);
972  float dxyErr2 = trk_dxyErr->at(ir);
973  dxyErr2 *= dxyErr2;
974  float dzErr2 = trk_dzErr->at(ir);
975  dzErr2 *= dzErr2;
976  float dzErrF2 = trk_dzErr->at(ir) * (pt * pz / p2);
977  dzErr2 *= dzErr2;
978  err.At(0, 0) = dxyErr2 * sP * sP + dzErrF2 * cP * cP;
979  err.At(0, 1) = -dxyErr2 * cP * sP + dzErrF2 * cP * sP;
980  err.At(1, 1) = dxyErr2 * cP * cP + dzErrF2 * sP * sP;
981  err.At(0, 2) = -dzErrF2 * cP * pt / pz;
982  err.At(1, 2) = -dzErrF2 * sP * pt / pz;
983  err.At(2, 2) = dzErr2 * std::pow((pt * pt / p2), 2);
984  err.At(3, 3) = std::pow(trk_ptErr->at(ir) / pt / pt, 2);
985  err.At(4, 4) = std::pow(trk_phiErr->at(ir), 2);
986  err.At(5, 5) = std::pow(trk_lambdaErr->at(ir), 2);
987  SVector3 pos = SVector3(trk_refpoint_x->at(ir), trk_refpoint_y->at(ir), trk_refpoint_z->at(ir));
988  SVector3 mom = SVector3(1.f / pt, phi, M_PI_2 - trk_lambda->at(ir));
989  TrackState state(trk_q->at(ir), pos, mom, err);
990  Track track(state, trk_nChi2->at(ir), trk_seedIdx->at(ir), 0, nullptr); //hits are filled later
991  track.setAlgorithm(TrackAlgorithm(trk_originalAlgo->at(ir)));
992  auto const& hTypes = trk_hitType->at(ir);
993  auto const& hIdxs = trk_hitIdx->at(ir);
994  for (unsigned int ip = 0; ip < hTypes.size(); ip++) {
995  unsigned int hidx = hIdxs[ip];
996  switch (HitType(hTypes[ip])) {
997  case HitType::Pixel: {
998  //cout << "pix=" << hidx << " track=" << cmsswTracks_.size() << endl;
999  pixHitRecIdx[hidx].push_back(cmsswTracks_.size());
1000  break;
1001  }
1002  case HitType::Strip: {
1003  //cout << "str=" << hidx << " track=" << cmsswTracks_.size() << endl;
1004  strHitRecIdx[hidx].push_back(cmsswTracks_.size());
1005  break;
1006  }
1007  case HitType::Glued: {
1008  if (not useMatched)
1009  throw std::logic_error("Tracks have glued hits, but matchedHit load is not configured");
1010  //cout << "glu=" << hidx << " track=" << cmsswTracks_.size() << endl;
1011  gluHitRecIdx[hidx].push_back(cmsswTracks_.size());
1012  break;
1013  }
1014  case HitType::Phase2OT: {
1015  //cout << "ph2=" << hidx << " track=" << cmsswTracks_.size() << endl;
1016  ph2HitRecIdx[hidx].push_back(cmsswTracks_.size());
1017  break;
1018  }
1019  case HitType::Invalid:
1020  break; //FIXME. Skip, really?
1021  default:
1022  throw std::logic_error("Track hit type can not be handled");
1023  } //switch( HitType
1024  }
1025  cmsswTracks_.push_back(track);
1026  }
1027 
1028  vector<vector<Hit>>& layerHits_ = EE.layerHits_;
1029  vector<vector<uint64_t>>& layerHitMasks_ = EE.layerHitMasks_;
1030  vector<MCHitInfo>& simHitsInfo_ = EE.simHitsInfo_;
1031  int totHits = 0;
1032  layerHits_.resize(nTotalLayers);
1033  layerHitMasks_.resize(nTotalLayers);
1034  for (unsigned int ipix = 0; ipix < pix_lay->size(); ++ipix) {
1035  int ilay = -1;
1036  ilay = lnc.convertLayerNumber(pix_det->at(ipix), pix_lay->at(ipix), useMatched, -1, pix_z->at(ipix) > 0);
1037  if (ilay < 0)
1038  continue;
1039 
1040  unsigned int imoduleid = tkinfo[ilay].short_id(pix_detId->at(ipix));
1041 
1042  int simTkIdxNt = bestTkIdx(pix_simHitIdx->at(ipix), pix_chargeFraction->at(ipix), ipix, HitType::Pixel);
1043  int simTkIdx = simTkIdxNt >= 0 ? simTrackIdx_[simTkIdxNt] : -1; //switch to index in simTracks_
1044  //cout << Form("pix lay=%i det=%i x=(%6.3f, %6.3f, %6.3f)",ilay+1,pix_det->at(ipix),pix_x->at(ipix),pix_y->at(ipix),pix_z->at(ipix)) << endl;
1045  SVector3 pos(pix_x->at(ipix), pix_y->at(ipix), pix_z->at(ipix));
1046  SMatrixSym33 err;
1047  err.At(0, 0) = pix_xx->at(ipix);
1048  err.At(1, 1) = pix_yy->at(ipix);
1049  err.At(2, 2) = pix_zz->at(ipix);
1050  err.At(0, 1) = pix_xy->at(ipix);
1051  err.At(0, 2) = pix_zx->at(ipix);
1052  err.At(1, 2) = pix_yz->at(ipix);
1053  if (simTkIdx >= 0) {
1054  simTracks_[simTkIdx].addHitIdx(layerHits_[ilay].size(), ilay, 0);
1055  }
1056  for (unsigned int is = 0; is < pixHitSeedIdx[ipix].size(); is++) {
1057  //cout << "xxx ipix=" << ipix << " seed=" << pixHitSeedIdx[ipix][is] << endl;
1058  seedTracks_[pixHitSeedIdx[ipix][is]].addHitIdx(layerHits_[ilay].size(), ilay, 0); //per-hit chi2 is not known
1059  }
1060  for (unsigned int ir = 0; ir < pixHitRecIdx[ipix].size(); ir++) {
1061  //cout << "xxx ipix=" << ipix << " recTrack=" << pixHitRecIdx[ipix][ir] << endl;
1062  cmsswTracks_[pixHitRecIdx[ipix][ir]].addHitIdx(layerHits_[ilay].size(), ilay, 0); //per-hit chi2 is not known
1063  }
1064  Hit hit(pos, err, totHits);
1065  hit.setupAsPixel(imoduleid, pix_csize_row->at(ipix), pix_csize_col->at(ipix));
1066  layerHits_[ilay].push_back(hit);
1067  if (writeHitIterMasks)
1068  layerHitMasks_[ilay].push_back(pix_usedMask->at(ipix));
1069  MCHitInfo hitInfo(simTkIdx, ilay, layerHits_[ilay].size() - 1, totHits);
1070  simHitsInfo_.push_back(hitInfo);
1071  totHits++;
1072  }
1073 
1074  if (hasPh2hits) {
1075  vector<int> ph2Idx(ph2_layer->size());
1076  for (unsigned int iph2 = 0; iph2 < ph2_layer->size(); ++iph2) {
1077  int ilay = -1;
1078  ilay = lnc.convertLayerNumber(
1079  ph2_subdet->at(iph2), ph2_layer->at(iph2), useMatched, ph2_isLower->at(iph2), ph2_z->at(iph2) > 0);
1080  if (useMatched && !ph2_isLower->at(iph2))
1081  continue;
1082  if (ilay == -1)
1083  continue;
1084 
1085  unsigned int imoduleid = tkinfo[ilay].short_id(ph2_detId->at(iph2));
1086 
1087  int simTkIdxNt = bestTkIdx(ph2_simHitIdx->at(iph2), ph2_chargeFraction_dummy, iph2, HitType::Phase2OT);
1088  int simTkIdx = simTkIdxNt >= 0 ? simTrackIdx_[simTkIdxNt] : -1; //switch to index in simTracks_
1089 
1090  SVector3 pos(ph2_x->at(iph2), ph2_y->at(iph2), ph2_z->at(iph2));
1091  SMatrixSym33 err;
1092  err.At(0, 0) = ph2_xx->at(iph2);
1093  err.At(1, 1) = ph2_yy->at(iph2);
1094  err.At(2, 2) = ph2_zz->at(iph2);
1095  err.At(0, 1) = ph2_xy->at(iph2);
1096  err.At(0, 2) = ph2_zx->at(iph2);
1097  err.At(1, 2) = ph2_yz->at(iph2);
1098  if (simTkIdx >= 0) {
1099  simTracks_[simTkIdx].addHitIdx(layerHits_[ilay].size(), ilay, 0);
1100  }
1101  for (unsigned int ir = 0; ir < ph2HitSeedIdx[iph2].size(); ir++) {
1102  //cout << "xxx iph2=" << iph2 << " seed=" << ph2HitSeedIdx[iph2][ir] << endl;
1103  seedTracks_[ph2HitSeedIdx[iph2][ir]].addHitIdx(layerHits_[ilay].size(), ilay, 0); //per-hit chi2 is not known
1104  }
1105  for (unsigned int ir = 0; ir < ph2HitRecIdx[iph2].size(); ir++) {
1106  //cout << "xxx iph2=" << iph2 << " recTrack=" << ph2HitRecIdx[iph2][ir] << endl;
1107  cmsswTracks_[ph2HitRecIdx[iph2][ir]].addHitIdx(layerHits_[ilay].size(), ilay, 0); //per-hit chi2 is not known
1108  }
1109  Hit hit(pos, err, totHits);
1110  hit.setupAsStrip(imoduleid, 0, 1);
1111  layerHits_[ilay].push_back(hit);
1112  if (writeHitIterMasks)
1113  layerHitMasks_[ilay].push_back(ph2_usedMask->at(iph2));
1114  MCHitInfo hitInfo(simTkIdx, ilay, layerHits_[ilay].size() - 1, totHits);
1115  simHitsInfo_.push_back(hitInfo);
1116  totHits++;
1117  }
1118  } else {
1119  if (useMatched) {
1120  for (unsigned int iglu = 0; iglu < glu_lay->size(); ++iglu) {
1121  if (glu_isBarrel->at(iglu) == 0)
1122  continue;
1123  int igluMono = glu_monoIdx->at(iglu);
1124  int simTkIdxNt =
1125  bestTkIdx(str_simHitIdx->at(igluMono), str_chargeFraction->at(igluMono), igluMono, HitType::Strip);
1126  int simTkIdx = simTkIdxNt >= 0 ? simTrackIdx_[simTkIdxNt] : -1; //switch to index in simTracks_
1127 
1128  int ilay = lnc.convertLayerNumber(glu_det->at(iglu), glu_lay->at(iglu), useMatched, -1, glu_z->at(iglu) > 0);
1129  // cout << Form("glu lay=%i det=%i bar=%i x=(%6.3f, %6.3f, %6.3f)",ilay+1,glu_det->at(iglu),glu_isBarrel->at(iglu),glu_x->at(iglu),glu_y->at(iglu),glu_z->at(iglu)) << endl;
1130  SVector3 pos(glu_x->at(iglu), glu_y->at(iglu), glu_z->at(iglu));
1131  SMatrixSym33 err;
1132  err.At(0, 0) = glu_xx->at(iglu);
1133  err.At(1, 1) = glu_yy->at(iglu);
1134  err.At(2, 2) = glu_zz->at(iglu);
1135  err.At(0, 1) = glu_xy->at(iglu);
1136  err.At(0, 2) = glu_zx->at(iglu);
1137  err.At(1, 2) = glu_yz->at(iglu);
1138  if (simTkIdx >= 0) {
1139  simTracks_[simTkIdx].addHitIdx(layerHits_[ilay].size(), ilay, 0);
1140  }
1141  for (unsigned int ir = 0; ir < gluHitSeedIdx[iglu].size(); ir++) {
1142  //cout << "xxx iglu=" << iglu << " seed=" << gluHitSeedIdx[iglu][ir] << endl;
1143  seedTracks_[gluHitSeedIdx[iglu][ir]].addHitIdx(
1144  layerHits_[ilay].size(), ilay, 0); //per-hit chi2 is not known
1145  }
1146  for (unsigned int ir = 0; ir < gluHitRecIdx[iglu].size(); ir++) {
1147  //cout << "xxx iglu=" << iglu << " recTrack=" << gluHitRecIdx[iglu][ir] << endl;
1148  cmsswTracks_[gluHitRecIdx[iglu][ir]].addHitIdx(
1149  layerHits_[ilay].size(), ilay, 0); //per-hit chi2 is not known
1150  }
1151 
1152  // QQQQ module-id-in-layer, adc and phi/theta spans are not done for matched hits.
1153  // Will we ever use / need this?
1154  assert(false && "Implement module-ids, cluster adc and spans for matched hits!");
1155 
1156  Hit hit(pos, err, totHits);
1157  layerHits_[ilay].push_back(hit);
1158  MCHitInfo hitInfo(simTkIdx, ilay, layerHits_[ilay].size() - 1, totHits);
1159  simHitsInfo_.push_back(hitInfo);
1160  totHits++;
1161  }
1162  }
1163 
1164  vector<int> strIdx;
1165  strIdx.resize(str_lay->size());
1166  for (unsigned int istr = 0; istr < str_lay->size(); ++istr) {
1167  int ilay = -1;
1168  ilay = lnc.convertLayerNumber(
1169  str_det->at(istr), str_lay->at(istr), useMatched, str_isStereo->at(istr), str_z->at(istr) > 0);
1170  if (useMatched && str_isBarrel->at(istr) == 1 && str_isStereo->at(istr))
1171  continue;
1172  if (ilay == -1)
1173  continue;
1174 
1175  unsigned int imoduleid = tkinfo[ilay].short_id(str_detId->at(istr));
1176 
1177  int simTkIdxNt = bestTkIdx(str_simHitIdx->at(istr), str_chargeFraction->at(istr), istr, HitType::Strip);
1178  int simTkIdx = simTkIdxNt >= 0 ? simTrackIdx_[simTkIdxNt] : -1; //switch to index in simTracks_
1179 
1180  bool passCCC = applyCCC ? (str_chargePerCM->at(istr) > cutValueCCC) : true;
1181 
1182  //if (str_onTrack->at(istr)==0) continue;//do not consider hits that are not on track!
1183  SVector3 pos(str_x->at(istr), str_y->at(istr), str_z->at(istr));
1184  SMatrixSym33 err;
1185  err.At(0, 0) = str_xx->at(istr);
1186  err.At(1, 1) = str_yy->at(istr);
1187  err.At(2, 2) = str_zz->at(istr);
1188  err.At(0, 1) = str_xy->at(istr);
1189  err.At(0, 2) = str_zx->at(istr);
1190  err.At(1, 2) = str_yz->at(istr);
1191  if (simTkIdx >= 0) {
1192  if (passCCC)
1193  simTracks_[simTkIdx].addHitIdx(layerHits_[ilay].size(), ilay, 0);
1194  else
1195  simTracks_[simTkIdx].addHitIdx(-9, ilay, 0);
1196  }
1197  for (unsigned int ir = 0; ir < strHitSeedIdx[istr].size(); ir++) {
1198  //cout << "xxx istr=" << istr << " seed=" << strHitSeedIdx[istr][ir] << endl;
1199  if (passCCC)
1200  seedTracks_[strHitSeedIdx[istr][ir]].addHitIdx(
1201  layerHits_[ilay].size(), ilay, 0); //per-hit chi2 is not known
1202  else
1203  seedTracks_[strHitSeedIdx[istr][ir]].addHitIdx(-9, ilay, 0);
1204  }
1205  for (unsigned int ir = 0; ir < strHitRecIdx[istr].size(); ir++) {
1206  //cout << "xxx istr=" << istr << " recTrack=" << strHitRecIdx[istr][ir] << endl;
1207  if (passCCC)
1208  cmsswTracks_[strHitRecIdx[istr][ir]].addHitIdx(
1209  layerHits_[ilay].size(), ilay, 0); //per-hit chi2 is not known
1210  else
1211  cmsswTracks_[strHitRecIdx[istr][ir]].addHitIdx(-9, ilay, 0);
1212  }
1213  if (passCCC) {
1214  Hit hit(pos, err, totHits);
1215  hit.setupAsStrip(imoduleid, str_chargePerCM->at(istr), str_csize->at(istr));
1216  layerHits_[ilay].push_back(hit);
1217  if (writeHitIterMasks)
1218  layerHitMasks_[ilay].push_back(str_usedMask->at(istr));
1219  MCHitInfo hitInfo(simTkIdx, ilay, layerHits_[ilay].size() - 1, totHits);
1220  simHitsInfo_.push_back(hitInfo);
1221  totHits++;
1222  }
1223  }
1224  }
1225 
1226  // Seed % hit statistics
1227  nstot += seedTracks_.size();
1228  for (unsigned int il = 0; il < layerHits_.size(); ++il) {
1229  int nh = layerHits_[il].size();
1230  nhitstot[il] += nh;
1231  }
1232 
1233  if (verbosity > 0) {
1234  int nt = simTracks_.size();
1235 
1236  int nl = layerHits_.size();
1237 
1238  int nm = simHitsInfo_.size();
1239 
1240  int ns = seedTracks_.size();
1241 
1242  int nr = cmsswTracks_.size();
1243 
1244  printf("number of simTracks %i\n", nt);
1245  printf("number of layerHits %i\n", nl);
1246  printf("number of simHitsInfo %i\n", nm);
1247  printf("number of seedTracks %i\n", ns);
1248  printf("number of recTracks %i\n", nr);
1249 
1250  if (verbosity > 1) {
1251  printf("\n");
1252  for (int il = 0; il < nl; ++il) {
1253  int nh = layerHits_[il].size();
1254  for (int ih = 0; ih < nh; ++ih) {
1255  printf("lay=%i idx=%i mcid=%i x=(%6.3f, %6.3f, %6.3f) r=%6.3f mask=0x%lx\n",
1256  il + 1,
1257  ih,
1258  layerHits_[il][ih].mcHitID(),
1259  layerHits_[il][ih].x(),
1260  layerHits_[il][ih].y(),
1261  layerHits_[il][ih].z(),
1262  sqrt(pow(layerHits_[il][ih].x(), 2) + pow(layerHits_[il][ih].y(), 2)),
1263  writeHitIterMasks ? layerHitMasks_[il][ih] : 0);
1264  }
1265  }
1266 
1267  for (int i = 0; i < nt; ++i) {
1268  float spt = sqrt(pow(simTracks_[i].px(), 2) + pow(simTracks_[i].py(), 2));
1269  printf(
1270  "sim track id=%i q=%2i p=(%6.3f, %6.3f, %6.3f) x=(%6.3f, %6.3f, %6.3f) pT=%7.4f nTotal=%i nFound=%i \n",
1271  i,
1272  simTracks_[i].charge(),
1273  simTracks_[i].px(),
1274  simTracks_[i].py(),
1275  simTracks_[i].pz(),
1276  simTracks_[i].x(),
1277  simTracks_[i].y(),
1278  simTracks_[i].z(),
1279  spt,
1280  simTracks_[i].nTotalHits(),
1281  simTracks_[i].nFoundHits());
1282  int nh = simTracks_[i].nTotalHits();
1283  for (int ih = 0; ih < nh; ++ih) {
1284  int hidx = simTracks_[i].getHitIdx(ih);
1285  int hlay = simTracks_[i].getHitLyr(ih);
1286  float hx = layerHits_[hlay][hidx].x();
1287  float hy = layerHits_[hlay][hidx].y();
1288  float hz = layerHits_[hlay][hidx].z();
1289  printf("track #%4i hit #%2i idx=%4i lay=%2i x=(% 8.3f, % 8.3f, % 8.3f) r=%8.3f\n",
1290  i,
1291  ih,
1292  hidx,
1293  hlay,
1294  hx,
1295  hy,
1296  hz,
1297  sqrt(hx * hx + hy * hy));
1298  }
1299  }
1300 
1301  for (int i = 0; i < ns; ++i) {
1302  printf("seed id=%i label=%i algo=%i q=%2i pT=%6.3f p=(%6.3f, %6.3f, %6.3f) x=(%6.3f, %6.3f, %6.3f)\n",
1303  i,
1304  seedTracks_[i].label(),
1305  (int)seedTracks_[i].algorithm(),
1306  seedTracks_[i].charge(),
1307  seedTracks_[i].pT(),
1308  seedTracks_[i].px(),
1309  seedTracks_[i].py(),
1310  seedTracks_[i].pz(),
1311  seedTracks_[i].x(),
1312  seedTracks_[i].y(),
1313  seedTracks_[i].z());
1314  int nh = seedTracks_[i].nTotalHits();
1315  for (int ih = 0; ih < nh; ++ih)
1316  printf("seed #%i hit #%i idx=%i\n", i, ih, seedTracks_[i].getHitIdx(ih));
1317  }
1318 
1319  if (writeRecTracks) {
1320  for (int i = 0; i < nr; ++i) {
1321  float spt = sqrt(pow(cmsswTracks_[i].px(), 2) + pow(cmsswTracks_[i].py(), 2));
1322  printf(
1323  "rec track id=%i label=%i algo=%i chi2=%6.3f q=%2i p=(%6.3f, %6.3f, %6.3f) x=(%6.3f, %6.3f, %6.3f) "
1324  "pT=%7.4f nTotal=%i nFound=%i \n",
1325  i,
1326  cmsswTracks_[i].label(),
1327  (int)cmsswTracks_[i].algorithm(),
1328  cmsswTracks_[i].chi2(),
1329  cmsswTracks_[i].charge(),
1330  cmsswTracks_[i].px(),
1331  cmsswTracks_[i].py(),
1332  cmsswTracks_[i].pz(),
1333  cmsswTracks_[i].x(),
1334  cmsswTracks_[i].y(),
1335  cmsswTracks_[i].z(),
1336  spt,
1337  cmsswTracks_[i].nTotalHits(),
1338  cmsswTracks_[i].nFoundHits());
1339  int nh = cmsswTracks_[i].nTotalHits();
1340  for (int ih = 0; ih < nh; ++ih) {
1341  int hidx = cmsswTracks_[i].getHitIdx(ih);
1342  int hlay = cmsswTracks_[i].getHitLyr(ih);
1343  float hx = layerHits_[hlay][hidx].x();
1344  float hy = layerHits_[hlay][hidx].y();
1345  float hz = layerHits_[hlay][hidx].z();
1346  printf("track #%4i hit #%2i idx=%4i lay=%2i x=(% 8.3f, % 8.3f, % 8.3f) r=%8.3f\n",
1347  i,
1348  ih,
1349  hidx,
1350  hlay,
1351  hx,
1352  hy,
1353  hz,
1354  sqrt(hx * hx + hy * hy));
1355  }
1356  }
1357  } //if (writeRecTracks){
1358 
1359  } //verbosity>1
1360  } //verbosity>0
1361  EE.write_out(data_file);
1362 
1363  savedEvents++;
1364  printf("end of event %lli\n", savedEvents);
1365  }
1366 
1367  data_file.CloseWrite(savedEvents);
1368  printf("\nSaved %lli events\n\n", savedEvents);
1369 
1370  printf("Average number of seeds per event %f\n", float(nstot) / float(savedEvents));
1371  for (unsigned int il = 0; il < nhitstot.size(); ++il)
1372  printf("Average number of hits in layer %3i = %7.2f\n",
1373  il,
1374  float(nhitstot[il]) / float(savedEvents)); //Includes those that failed the cluster charge cut
1375 
1376  if (!hasPh2hits)
1377  printf("Out of %i hits, %i failed the cut", numTotalStr, numFailCCC);
1378 
1379  //========================================================================
1380 
1381  printf("\n\n================================================================\n");
1382  printf("=== Max module id for %u layers\n", nTotalLayers);
1383  printf("================================================================\n");
1384  for (unsigned int ii = 0; ii < nTotalLayers; ++ii) {
1385  printf("Layer%2d : %d\n", ii, tkinfo[ii].n_modules());
1386  }
1387 }
void CloseWrite(int n_written)
Definition: Event.cc:978
size
Write out results.
Definition: start.py:1
void next_arg_or_die(lStr_t &args, lStr_i &i)
ROOT::Math::SMatrix< float, 6, 6, ROOT::Math::MatRepSym< float, 6 > > SMatrixSym66
Definition: MatrixSTypes.h:8
void read_bin_file(const std::string &fname)
Definition: TrackerInfo.cc:168
bool next_arg_option(lStr_t &args, lStr_i &i)
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
ROOT::Math::SVector< float, 3 > SVector3
Definition: MatrixSTypes.h:14
#define M_PI_2
assert(be >=bs)
float float float z
constexpr bool useMatched
lStr_t::iterator lStr_i
Definition: mkFit.cc:462
char const * label
constexpr int cleanSimTrack_minSimHits
T sqrt(T t)
Definition: SSEVec.h:19
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e< void, edm::EventID const &, edm::Timestamp const & > We also list in braces which AR_WATCH_USING_METHOD_ is used for those or
Definition: Activities.doc:12
int n_layers() const
Definition: TrackerInfo.h:161
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
TrackBase::TrackAlgorithm TrackAlgorithm
double f[11][100]
void printHelp(const char *av0)
uint32_t nh
HitType
int nt
Definition: AMPTWrapper.h:42
ii
Definition: cuy.py:589
const int verbosity
caConstants::TupleMultiplicity const CAHitNtupletGeneratorKernelsGPU::HitToTuple const cms::cuda::AtomicPairCounter GPUCACell const *__restrict__ uint32_t const *__restrict__ gpuPixelDoublets::CellNeighborsVector const gpuPixelDoublets::CellTracksVector const GPUCACell::OuterHitOfCell const int32_t nHits
void openWrite(const std::string &fname, int n_layers, int n_ev, int extra_sections=0)
Definition: Event.cc:930
float x
constexpr int cleanSimTrack_minRecHits
std::list< std::string > lStr_t
Definition: mkFit.cc:461
ROOT::Math::SVector< double, 3 > SVector3
Definition: V0Fitter.cc:47
ROOT::Math::SMatrix< float, 3, 3, ROOT::Math::MatRepSym< float, 3 > > SMatrixSym33
Definition: MatrixSTypes.h:13
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
Definition: event.py:1
def exit(msg="")

◆ next_arg_option()

bool next_arg_option ( lStr_t args,
lStr_i i 
)

Definition at line 35 of file WriteMemoryFile.cc.

References writedatasetfile::args, mps_fire::i, and dqmiolumiharvest::j.

Referenced by main().

35  {
36  lStr_i j = i;
37  if (++j == args.end() || ((*j)[0] == '-')) {
38  return false;
39  }
40  i = j;
41  return true;
42 }
lStr_t::iterator lStr_i
Definition: mkFit.cc:462

◆ next_arg_or_die()

void next_arg_or_die ( lStr_t args,
lStr_i i 
)

Definition at line 26 of file WriteMemoryFile.cc.

References writedatasetfile::args, DMR_cfg::cerr, beamvalidation::exit(), mps_fire::i, and dqmiolumiharvest::j.

Referenced by main().

26  {
27  lStr_i j = i;
28  if (++j == args.end() || ((*j)[0] == '-')) {
29  std::cerr << "Error: option " << *i << " requires an argument.\n";
30  exit(1);
31  }
32  i = j;
33 }
lStr_t::iterator lStr_i
Definition: mkFit.cc:462
def exit(msg="")

◆ printHelp()

void printHelp ( const char *  av0)

Definition at line 44 of file WriteMemoryFile.cc.

Referenced by main().

44  {
45  printf(
46  "Usage: %s [options]\n"
47  "Options:\n"
48  " --input <str> input file\n"
49  " --output <str> output file\n"
50  " --geo <file> binary TrackerInfo geometry (def: CMS-phase1.bin)\n"
51  " --verbosity <num> print details (0 quiet, 1 print counts, 2 print all; def: 0)\n"
52  " --maxevt <num> maxevt events to write (-1 for everything in the file def: -1)\n"
53  " --clean-sim-tracks apply sim track cleaning (def: no cleaning)\n"
54  " --write-all-events write all events (def: skip events with 0 simtracks or seeds)\n"
55  " --write-rec-tracks write rec tracks (def: not written)\n"
56  " --apply-ccc apply cluster charge cut to strip hits (def: false)\n"
57  " --all-seeds merge all seeds from the input file (def: false)\n",
58  av0);
59 }

Variable Documentation

◆ cleanSimTrack_minRecHits

constexpr int cleanSimTrack_minRecHits = 2

Definition at line 19 of file WriteMemoryFile.cc.

Referenced by main().

◆ cleanSimTrack_minSimHits

constexpr int cleanSimTrack_minSimHits = 3

Definition at line 18 of file WriteMemoryFile.cc.

Referenced by main().

◆ useMatched

constexpr bool useMatched = false