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