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 63 of file WriteMemoryFile.cc.

References qcdUeDQM_cfi::algorithm, dir2webdir::argc, GCPpyPlots::argv, cms::cuda::assert(), cms::cuda::bs, ALCARECOTkAlJpsiMuMu_cff::charge, nano_mu_local_reco_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, nano_mu_digi_cff::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, conifer::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_2023v12_cff::track, nano_mu_reco_cff::trk_algo, mitigatedMETSequence_cff::U, useMatched, verbosity, and x.

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

◆ 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:454
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  " --help print help and exit\n"
49  " --input <str> input file\n"
50  " --output <str> output file\n"
51  " --geo <file> binary TrackerInfo geometry (def: CMS-phase1.bin)\n"
52  " --verbosity <num> print details (0 quiet, 1 print counts, 2 print all; def: 0)\n"
53  " --maxevt <num> maxevt events to write (-1 for everything in the file def: -1)\n"
54  " --clean-sim-tracks apply sim track cleaning (def: no cleaning)\n"
55  " --write-all-events write all events (def: skip events with 0 simtracks or seeds)\n"
56  " --write-rec-tracks write rec tracks (def: not written)\n"
57  " --apply-ccc apply cluster charge cut to strip hits (def: false)\n"
58  " --all-seeds write all seeds from the input file, not only initialStep and hltIter0 (def: "
59  "false)\n",
60  av0);
61 }

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