File indexing completed on 2022-04-07 05:53:41
0001 #include "TFile.h"
0002 #include "TTree.h"
0003
0004 #include <iostream>
0005 #include <list>
0006 #include <unordered_map>
0007 #include "RecoTracker/MkFitCore/standalone/Event.h"
0008 #include "RecoTracker/MkFitCore/interface/Config.h"
0009 #include "RecoTracker/MkFitCore/interface/TrackerInfo.h"
0010 #include "RecoTracker/MkFitCMS/interface/LayerNumberConverter.h"
0011
0012 using namespace mkfit;
0013
0014 using TrackAlgorithm = TrackBase::TrackAlgorithm;
0015
0016 constexpr bool useMatched = false;
0017
0018 constexpr int cleanSimTrack_minSimHits = 3;
0019 constexpr int cleanSimTrack_minRecHits = 2;
0020
0021
0022 enum class HitType { Pixel = 0, Strip = 1, Glued = 2, Invalid = 3, Phase2OT = 4, Unknown = 99 };
0023
0024 typedef std::list<std::string> lStr_t;
0025 typedef lStr_t::iterator lStr_i;
0026 void next_arg_or_die(lStr_t& args, lStr_i& i) {
0027 lStr_i j = i;
0028 if (++j == args.end() || ((*j)[0] == '-')) {
0029 std::cerr << "Error: option " << *i << " requires an argument.\n";
0030 exit(1);
0031 }
0032 i = j;
0033 }
0034
0035 bool next_arg_option(lStr_t& args, lStr_i& i) {
0036 lStr_i j = i;
0037 if (++j == args.end() || ((*j)[0] == '-')) {
0038 return false;
0039 }
0040 i = j;
0041 return true;
0042 }
0043
0044 void printHelp(const char* av0) {
0045 printf(
0046 "Usage: %s [options]\n"
0047 "Options:\n"
0048 " --input <str> input file\n"
0049 " --output <str> output file\n"
0050 " --geo <file> binary TrackerInfo geometry (def: CMS-2017.bin)\n"
0051 " --verbosity <num> print details (0 quiet, 1 print counts, 2 print all; def: 0)\n"
0052 " --maxevt <num> maxevt events to write (-1 for everything in the file def: -1)\n"
0053 " --clean-sim-tracks apply sim track cleaning (def: no cleaning)\n"
0054 " --write-all-events write all events (def: skip events with 0 simtracks or seeds)\n"
0055 " --write-rec-tracks write rec tracks (def: not written)\n"
0056 " --apply-ccc apply cluster charge cut to strip hits (def: false)\n"
0057 " --all-seeds merge all seeds from the input file (def: false)\n",
0058 av0);
0059 }
0060
0061 int main(int argc, char* argv[]) {
0062 bool haveInput = false;
0063 std::string inputFileName;
0064 bool haveOutput = false;
0065 std::string outputFileName;
0066 std::string geoFileName("CMS-2017.bin");
0067
0068 bool cleanSimTracks = false;
0069 bool writeAllEvents = false;
0070 bool writeRecTracks = false;
0071 bool writeHitIterMasks = false;
0072 bool applyCCC = false;
0073 bool allSeeds = false;
0074
0075 int verbosity = 0;
0076 long long maxevt = -1;
0077
0078 int cutValueCCC = 1620;
0079
0080 lStr_t mArgs;
0081 for (int i = 1; i < argc; ++i) {
0082 mArgs.push_back(argv[i]);
0083 }
0084
0085 lStr_i i = mArgs.begin();
0086 while (i != mArgs.end()) {
0087 lStr_i start = i;
0088
0089 if (*i == "-h" || *i == "-help" || *i == "--help") {
0090 printHelp(argv[0]);
0091 } else if (*i == "--input") {
0092 next_arg_or_die(mArgs, i);
0093 inputFileName = *i;
0094 haveInput = true;
0095 } else if (*i == "--output") {
0096 next_arg_or_die(mArgs, i);
0097 outputFileName = *i;
0098 haveOutput = true;
0099 } else if (*i == "--geo") {
0100 next_arg_or_die(mArgs, i);
0101 geoFileName = *i;
0102 } else if (*i == "--verbosity") {
0103 next_arg_or_die(mArgs, i);
0104 verbosity = std::atoi(i->c_str());
0105 } else if (*i == "--maxevt") {
0106 next_arg_or_die(mArgs, i);
0107 maxevt = std::atoi(i->c_str());
0108 } else if (*i == "--clean-sim-tracks") {
0109 cleanSimTracks = true;
0110 } else if (*i == "--write-all-events") {
0111 writeAllEvents = true;
0112 } else if (*i == "--write-rec-tracks") {
0113 writeRecTracks = true;
0114 } else if (*i == "--write-hit-iter-masks") {
0115 writeHitIterMasks = true;
0116 } else if (*i == "--apply-ccc") {
0117 applyCCC = true;
0118 if (next_arg_option(mArgs, i)) {
0119 cutValueCCC = std::atoi(i->c_str());
0120 }
0121 } else if (*i == "--all-seeds") {
0122 allSeeds = true;
0123 } else {
0124 fprintf(stderr, "Error: Unknown option/argument '%s'.\n", i->c_str());
0125 printHelp(argv[0]);
0126 exit(1);
0127 }
0128 mArgs.erase(start, ++i);
0129 }
0130
0131 if (not haveOutput or not haveInput) {
0132 fprintf(stderr, "Error: both input and output are required\n");
0133 printHelp(argv[0]);
0134 exit(1);
0135 }
0136
0137 using namespace std;
0138
0139 TrackerInfo tkinfo;
0140 tkinfo.read_bin_file(geoFileName);
0141 LayerNumberConverter lnc(TkLayout::phase1);
0142 const unsigned int nTotalLayers = lnc.nLayers();
0143 assert(nTotalLayers == (unsigned int)tkinfo.n_layers());
0144
0145 int nstot = 0;
0146 std::vector<int> nhitstot(nTotalLayers, 0);
0147
0148 TFile* f = TFile::Open(inputFileName.c_str());
0149 if (f == 0) {
0150 fprintf(stderr, "Failed opening input root file '%s'\n", inputFileName.c_str());
0151 exit(1);
0152 }
0153
0154 TTree* t = (TTree*)f->Get("trackingNtuple/tree");
0155
0156 unsigned long long event;
0157 t->SetBranchAddress("event", &event);
0158
0159
0160 std::vector<float>* sim_eta = 0;
0161 std::vector<float>* sim_px = 0;
0162 std::vector<float>* sim_py = 0;
0163 std::vector<float>* sim_pz = 0;
0164 std::vector<int>* sim_parentVtxIdx = 0;
0165 std::vector<int>* sim_q = 0;
0166 std::vector<int>* sim_event = 0;
0167 std::vector<int>* sim_bunchCrossing = 0;
0168 std::vector<int>* sim_nValid = 0;
0169 t->SetBranchAddress("sim_eta", &sim_eta);
0170 t->SetBranchAddress("sim_px", &sim_px);
0171 t->SetBranchAddress("sim_py", &sim_py);
0172 t->SetBranchAddress("sim_pz", &sim_pz);
0173 t->SetBranchAddress("sim_parentVtxIdx", &sim_parentVtxIdx);
0174 t->SetBranchAddress("sim_q", &sim_q);
0175 t->SetBranchAddress("sim_event", &sim_event);
0176 t->SetBranchAddress("sim_bunchCrossing", &sim_bunchCrossing);
0177 t->SetBranchAddress("sim_nValid", &sim_nValid);
0178
0179 std::vector<vector<int>>* sim_trkIdx = 0;
0180 t->SetBranchAddress("sim_trkIdx", &sim_trkIdx);
0181
0182
0183 std::vector<float>* simvtx_x = 0;
0184 std::vector<float>* simvtx_y = 0;
0185 std::vector<float>* simvtx_z = 0;
0186 t->SetBranchAddress("simvtx_x", &simvtx_x);
0187 t->SetBranchAddress("simvtx_y", &simvtx_y);
0188 t->SetBranchAddress("simvtx_z", &simvtx_z);
0189
0190
0191 std::vector<short>* simhit_process = 0;
0192 std::vector<int>* simhit_particle = 0;
0193 std::vector<int>* simhit_simTrkIdx = 0;
0194 std::vector<float>* simhit_x = 0;
0195 std::vector<float>* simhit_y = 0;
0196 std::vector<float>* simhit_z = 0;
0197 std::vector<float>* simhit_px = 0;
0198 std::vector<float>* simhit_py = 0;
0199 std::vector<float>* simhit_pz = 0;
0200 t->SetBranchAddress("simhit_process", &simhit_process);
0201 t->SetBranchAddress("simhit_particle", &simhit_particle);
0202 t->SetBranchAddress("simhit_simTrkIdx", &simhit_simTrkIdx);
0203 t->SetBranchAddress("simhit_x", &simhit_x);
0204 t->SetBranchAddress("simhit_y", &simhit_y);
0205 t->SetBranchAddress("simhit_z", &simhit_z);
0206 t->SetBranchAddress("simhit_px", &simhit_px);
0207 t->SetBranchAddress("simhit_py", &simhit_py);
0208 t->SetBranchAddress("simhit_pz", &simhit_pz);
0209
0210 std::vector<std::vector<int>>* simhit_hitIdx = 0;
0211 t->SetBranchAddress("simhit_hitIdx", &simhit_hitIdx);
0212 std::vector<std::vector<int>>* simhit_hitType = 0;
0213 t->SetBranchAddress("simhit_hitType", &simhit_hitType);
0214
0215
0216 std::vector<int>* trk_q = 0;
0217 std::vector<unsigned int>* trk_nValid = 0;
0218 std::vector<int>* trk_seedIdx = 0;
0219 std::vector<unsigned long long>* trk_algoMask = 0;
0220 std::vector<unsigned int>* trk_algo = 0;
0221 std::vector<unsigned int>* trk_originalAlgo = 0;
0222 std::vector<float>* trk_nChi2 = 0;
0223 std::vector<float>* trk_px = 0;
0224 std::vector<float>* trk_py = 0;
0225 std::vector<float>* trk_pz = 0;
0226 std::vector<float>* trk_pt = 0;
0227 std::vector<float>* trk_phi = 0;
0228 std::vector<float>* trk_lambda = 0;
0229 std::vector<float>* trk_refpoint_x = 0;
0230 std::vector<float>* trk_refpoint_y = 0;
0231 std::vector<float>* trk_refpoint_z = 0;
0232 std::vector<float>* trk_dxyErr = 0;
0233 std::vector<float>* trk_dzErr = 0;
0234 std::vector<float>* trk_ptErr = 0;
0235 std::vector<float>* trk_phiErr = 0;
0236 std::vector<float>* trk_lambdaErr = 0;
0237 t->SetBranchAddress("trk_q", &trk_q);
0238 t->SetBranchAddress("trk_nValid", &trk_nValid);
0239 t->SetBranchAddress("trk_seedIdx", &trk_seedIdx);
0240 t->SetBranchAddress("trk_algoMask", &trk_algoMask);
0241 t->SetBranchAddress("trk_algo", &trk_algo);
0242 t->SetBranchAddress("trk_originalAlgo", &trk_originalAlgo);
0243 t->SetBranchAddress("trk_nChi2", &trk_nChi2);
0244 t->SetBranchAddress("trk_px", &trk_px);
0245 t->SetBranchAddress("trk_py", &trk_py);
0246 t->SetBranchAddress("trk_pz", &trk_pz);
0247 t->SetBranchAddress("trk_pt", &trk_pt);
0248 t->SetBranchAddress("trk_phi", &trk_phi);
0249 t->SetBranchAddress("trk_lambda", &trk_lambda);
0250 t->SetBranchAddress("trk_refpoint_x", &trk_refpoint_x);
0251 t->SetBranchAddress("trk_refpoint_y", &trk_refpoint_y);
0252 t->SetBranchAddress("trk_refpoint_z", &trk_refpoint_z);
0253 t->SetBranchAddress("trk_dxyErr", &trk_dxyErr);
0254 t->SetBranchAddress("trk_dzErr", &trk_dzErr);
0255 t->SetBranchAddress("trk_ptErr", &trk_ptErr);
0256 t->SetBranchAddress("trk_phiErr", &trk_phiErr);
0257 t->SetBranchAddress("trk_lambdaErr", &trk_lambdaErr);
0258
0259 std::vector<std::vector<int>>* trk_hitIdx = 0;
0260 t->SetBranchAddress("trk_hitIdx", &trk_hitIdx);
0261 std::vector<std::vector<int>>* trk_hitType = 0;
0262 t->SetBranchAddress("trk_hitType", &trk_hitType);
0263
0264
0265 std::vector<float>* see_stateTrajGlbX = 0;
0266 std::vector<float>* see_stateTrajGlbY = 0;
0267 std::vector<float>* see_stateTrajGlbZ = 0;
0268 std::vector<float>* see_stateTrajGlbPx = 0;
0269 std::vector<float>* see_stateTrajGlbPy = 0;
0270 std::vector<float>* see_stateTrajGlbPz = 0;
0271 std::vector<float>* see_eta = 0;
0272 std::vector<float>* see_pt = 0;
0273 std::vector<float>* see_stateCcov00 = 0;
0274 std::vector<float>* see_stateCcov01 = 0;
0275 std::vector<float>* see_stateCcov02 = 0;
0276 std::vector<float>* see_stateCcov03 = 0;
0277 std::vector<float>* see_stateCcov04 = 0;
0278 std::vector<float>* see_stateCcov05 = 0;
0279 std::vector<float>* see_stateCcov11 = 0;
0280 std::vector<float>* see_stateCcov12 = 0;
0281 std::vector<float>* see_stateCcov13 = 0;
0282 std::vector<float>* see_stateCcov14 = 0;
0283 std::vector<float>* see_stateCcov15 = 0;
0284 std::vector<float>* see_stateCcov22 = 0;
0285 std::vector<float>* see_stateCcov23 = 0;
0286 std::vector<float>* see_stateCcov24 = 0;
0287 std::vector<float>* see_stateCcov25 = 0;
0288 std::vector<float>* see_stateCcov33 = 0;
0289 std::vector<float>* see_stateCcov34 = 0;
0290 std::vector<float>* see_stateCcov35 = 0;
0291 std::vector<float>* see_stateCcov44 = 0;
0292 std::vector<float>* see_stateCcov45 = 0;
0293 std::vector<float>* see_stateCcov55 = 0;
0294 std::vector<std::vector<float>>* see_stateCurvCov = 0;
0295 std::vector<int>* see_q = 0;
0296 std::vector<unsigned int>* see_algo = 0;
0297 t->SetBranchAddress("see_stateTrajGlbX", &see_stateTrajGlbX);
0298 t->SetBranchAddress("see_stateTrajGlbY", &see_stateTrajGlbY);
0299 t->SetBranchAddress("see_stateTrajGlbZ", &see_stateTrajGlbZ);
0300 t->SetBranchAddress("see_stateTrajGlbPx", &see_stateTrajGlbPx);
0301 t->SetBranchAddress("see_stateTrajGlbPy", &see_stateTrajGlbPy);
0302 t->SetBranchAddress("see_stateTrajGlbPz", &see_stateTrajGlbPz);
0303 t->SetBranchAddress("see_eta", &see_eta);
0304 t->SetBranchAddress("see_pt", &see_pt);
0305
0306 bool hasCartCov = t->GetBranch("see_stateCcov00") != nullptr;
0307 if (hasCartCov) {
0308 t->SetBranchAddress("see_stateCcov00", &see_stateCcov00);
0309 t->SetBranchAddress("see_stateCcov01", &see_stateCcov01);
0310 t->SetBranchAddress("see_stateCcov02", &see_stateCcov02);
0311 t->SetBranchAddress("see_stateCcov03", &see_stateCcov03);
0312 t->SetBranchAddress("see_stateCcov04", &see_stateCcov04);
0313 t->SetBranchAddress("see_stateCcov05", &see_stateCcov05);
0314 t->SetBranchAddress("see_stateCcov11", &see_stateCcov11);
0315 t->SetBranchAddress("see_stateCcov12", &see_stateCcov12);
0316 t->SetBranchAddress("see_stateCcov13", &see_stateCcov13);
0317 t->SetBranchAddress("see_stateCcov14", &see_stateCcov14);
0318 t->SetBranchAddress("see_stateCcov15", &see_stateCcov15);
0319 t->SetBranchAddress("see_stateCcov22", &see_stateCcov22);
0320 t->SetBranchAddress("see_stateCcov23", &see_stateCcov23);
0321 t->SetBranchAddress("see_stateCcov24", &see_stateCcov24);
0322 t->SetBranchAddress("see_stateCcov25", &see_stateCcov25);
0323 t->SetBranchAddress("see_stateCcov33", &see_stateCcov33);
0324 t->SetBranchAddress("see_stateCcov34", &see_stateCcov34);
0325 t->SetBranchAddress("see_stateCcov35", &see_stateCcov35);
0326 t->SetBranchAddress("see_stateCcov44", &see_stateCcov44);
0327 t->SetBranchAddress("see_stateCcov45", &see_stateCcov45);
0328 t->SetBranchAddress("see_stateCcov55", &see_stateCcov55);
0329 } else {
0330 t->SetBranchAddress("see_stateCurvCov", &see_stateCurvCov);
0331 }
0332 t->SetBranchAddress("see_q", &see_q);
0333 t->SetBranchAddress("see_algo", &see_algo);
0334
0335 std::vector<std::vector<int>>* see_hitIdx = 0;
0336 t->SetBranchAddress("see_hitIdx", &see_hitIdx);
0337 std::vector<std::vector<int>>* see_hitType = 0;
0338 t->SetBranchAddress("see_hitType", &see_hitType);
0339
0340
0341 vector<unsigned short>* pix_det = 0;
0342 vector<unsigned short>* pix_lay = 0;
0343 vector<unsigned int>* pix_detId = 0;
0344 vector<float>* pix_x = 0;
0345 vector<float>* pix_y = 0;
0346 vector<float>* pix_z = 0;
0347 vector<float>* pix_xx = 0;
0348 vector<float>* pix_xy = 0;
0349 vector<float>* pix_yy = 0;
0350 vector<float>* pix_yz = 0;
0351 vector<float>* pix_zz = 0;
0352 vector<float>* pix_zx = 0;
0353 vector<int>* pix_csize_col = 0;
0354 vector<int>* pix_csize_row = 0;
0355 vector<uint64_t>* pix_usedMask = 0;
0356
0357 bool has910_det_lay = t->GetBranch("pix_det") == nullptr;
0358 if (has910_det_lay) {
0359 t->SetBranchAddress("pix_subdet", &pix_det);
0360 t->SetBranchAddress("pix_layer", &pix_lay);
0361 } else {
0362 t->SetBranchAddress("pix_det", &pix_det);
0363 t->SetBranchAddress("pix_lay", &pix_lay);
0364 }
0365 t->SetBranchAddress("pix_detId", &pix_detId);
0366 t->SetBranchAddress("pix_x", &pix_x);
0367 t->SetBranchAddress("pix_y", &pix_y);
0368 t->SetBranchAddress("pix_z", &pix_z);
0369 t->SetBranchAddress("pix_xx", &pix_xx);
0370 t->SetBranchAddress("pix_xy", &pix_xy);
0371 t->SetBranchAddress("pix_yy", &pix_yy);
0372 t->SetBranchAddress("pix_yz", &pix_yz);
0373 t->SetBranchAddress("pix_zz", &pix_zz);
0374 t->SetBranchAddress("pix_zx", &pix_zx);
0375 t->SetBranchAddress("pix_clustSizeCol", &pix_csize_col);
0376 t->SetBranchAddress("pix_clustSizeRow", &pix_csize_row);
0377 if (writeHitIterMasks) {
0378 t->SetBranchAddress("pix_usedMask", &pix_usedMask);
0379 }
0380
0381 vector<vector<int>>* pix_simHitIdx = 0;
0382 t->SetBranchAddress("pix_simHitIdx", &pix_simHitIdx);
0383 vector<vector<float>>* pix_chargeFraction = 0;
0384 t->SetBranchAddress("pix_chargeFraction", &pix_chargeFraction);
0385
0386
0387 vector<short>* glu_isBarrel = 0;
0388 vector<unsigned int>* glu_det = 0;
0389 vector<unsigned int>* glu_lay = 0;
0390 vector<unsigned int>* glu_detId = 0;
0391 vector<int>* glu_monoIdx = 0;
0392 vector<int>* glu_stereoIdx = 0;
0393 vector<float>* glu_x = 0;
0394 vector<float>* glu_y = 0;
0395 vector<float>* glu_z = 0;
0396 vector<float>* glu_xx = 0;
0397 vector<float>* glu_xy = 0;
0398 vector<float>* glu_yy = 0;
0399 vector<float>* glu_yz = 0;
0400 vector<float>* glu_zz = 0;
0401 vector<float>* glu_zx = 0;
0402 t->SetBranchAddress("glu_isBarrel", &glu_isBarrel);
0403 if (has910_det_lay) {
0404 t->SetBranchAddress("glu_subdet", &glu_det);
0405 t->SetBranchAddress("glu_layer", &glu_lay);
0406 } else {
0407 t->SetBranchAddress("glu_det", &glu_det);
0408 t->SetBranchAddress("glu_lay", &glu_lay);
0409 }
0410 t->SetBranchAddress("glu_detId", &glu_detId);
0411 t->SetBranchAddress("glu_monoIdx", &glu_monoIdx);
0412 t->SetBranchAddress("glu_stereoIdx", &glu_stereoIdx);
0413 t->SetBranchAddress("glu_x", &glu_x);
0414 t->SetBranchAddress("glu_y", &glu_y);
0415 t->SetBranchAddress("glu_z", &glu_z);
0416 t->SetBranchAddress("glu_xx", &glu_xx);
0417 t->SetBranchAddress("glu_xy", &glu_xy);
0418 t->SetBranchAddress("glu_yy", &glu_yy);
0419 t->SetBranchAddress("glu_yz", &glu_yz);
0420 t->SetBranchAddress("glu_zz", &glu_zz);
0421 t->SetBranchAddress("glu_zx", &glu_zx);
0422
0423 vector<short>* str_isBarrel = 0;
0424 vector<short>* str_isStereo = 0;
0425 vector<unsigned int>* str_det = 0;
0426 vector<unsigned int>* str_lay = 0;
0427 vector<unsigned int>* str_detId = 0;
0428 vector<unsigned int>* str_simType = 0;
0429 vector<float>* str_x = 0;
0430 vector<float>* str_y = 0;
0431 vector<float>* str_z = 0;
0432 vector<float>* str_xx = 0;
0433 vector<float>* str_xy = 0;
0434 vector<float>* str_yy = 0;
0435 vector<float>* str_yz = 0;
0436 vector<float>* str_zz = 0;
0437 vector<float>* str_zx = 0;
0438 vector<float>* str_chargePerCM = 0;
0439 vector<int>* str_csize = 0;
0440 vector<uint64_t>* str_usedMask = 0;
0441 t->SetBranchAddress("str_isBarrel", &str_isBarrel);
0442 t->SetBranchAddress("str_isStereo", &str_isStereo);
0443 if (has910_det_lay) {
0444 t->SetBranchAddress("str_subdet", &str_det);
0445 t->SetBranchAddress("str_layer", &str_lay);
0446 } else {
0447 t->SetBranchAddress("str_det", &str_det);
0448 t->SetBranchAddress("str_lay", &str_lay);
0449 }
0450 t->SetBranchAddress("str_detId", &str_detId);
0451 t->SetBranchAddress("str_simType", &str_simType);
0452 t->SetBranchAddress("str_x", &str_x);
0453 t->SetBranchAddress("str_y", &str_y);
0454 t->SetBranchAddress("str_z", &str_z);
0455 t->SetBranchAddress("str_xx", &str_xx);
0456 t->SetBranchAddress("str_xy", &str_xy);
0457 t->SetBranchAddress("str_yy", &str_yy);
0458 t->SetBranchAddress("str_yz", &str_yz);
0459 t->SetBranchAddress("str_zz", &str_zz);
0460 t->SetBranchAddress("str_zx", &str_zx);
0461 t->SetBranchAddress("str_chargePerCM", &str_chargePerCM);
0462 t->SetBranchAddress("str_clustSize", &str_csize);
0463 if (writeHitIterMasks) {
0464 t->SetBranchAddress("str_usedMask", &str_usedMask);
0465 }
0466
0467 vector<vector<int>>* str_simHitIdx = 0;
0468 t->SetBranchAddress("str_simHitIdx", &str_simHitIdx);
0469 vector<vector<float>>* str_chargeFraction = 0;
0470 t->SetBranchAddress("str_chargeFraction", &str_chargeFraction);
0471
0472
0473 float bsp_x;
0474 float bsp_y;
0475 float bsp_z;
0476 float bsp_sigmax;
0477 float bsp_sigmay;
0478 float bsp_sigmaz;
0479 t->SetBranchAddress("bsp_x", &bsp_x);
0480 t->SetBranchAddress("bsp_y", &bsp_y);
0481 t->SetBranchAddress("bsp_z", &bsp_z);
0482 t->SetBranchAddress("bsp_sigmax", &bsp_sigmax);
0483 t->SetBranchAddress("bsp_sigmay", &bsp_sigmay);
0484 t->SetBranchAddress("bsp_sigmaz", &bsp_sigmaz);
0485
0486 long long totentries = t->GetEntries();
0487 long long savedEvents = 0;
0488
0489 DataFile data_file;
0490 int outOptions = DataFile::ES_Seeds;
0491 if (writeRecTracks)
0492 outOptions |= DataFile::ES_CmsswTracks;
0493 if (writeHitIterMasks)
0494 outOptions |= DataFile::ES_HitIterMasks;
0495 outOptions |= DataFile::ES_BeamSpot;
0496
0497 if (maxevt < 0)
0498 maxevt = totentries;
0499 data_file.openWrite(outputFileName, static_cast<int>(nTotalLayers), std::min(maxevt, totentries), outOptions);
0500
0501 Event EE(0, static_cast<int>(nTotalLayers));
0502
0503 int numFailCCC = 0;
0504 int numTotalStr = 0;
0505
0506
0507 for (long long i = 0; savedEvents < maxevt && i < totentries && i < maxevt; ++i) {
0508 EE.reset(i);
0509
0510 cout << "process entry i=" << i << " out of " << totentries << ", saved so far " << savedEvents
0511 << ", with max=" << maxevt << endl;
0512
0513 t->GetEntry(i);
0514
0515 cout << "edm event=" << event << endl;
0516
0517 auto& bs = EE.beamSpot_;
0518 bs.x = bsp_x;
0519 bs.y = bsp_y;
0520 bs.z = bsp_z;
0521 bs.sigmaZ = bsp_sigmaz;
0522 bs.beamWidthX = bsp_sigmax;
0523 bs.beamWidthY = bsp_sigmay;
0524
0525
0526 for (unsigned int istr = 0; istr < str_lay->size(); ++istr) {
0527 if (str_chargePerCM->at(istr) < cutValueCCC)
0528 numFailCCC++;
0529 numTotalStr++;
0530 }
0531
0532 auto nSims = sim_q->size();
0533 if (nSims == 0) {
0534 cout << "branches not loaded" << endl;
0535 exit(1);
0536 }
0537 if (verbosity > 0)
0538 std::cout << __FILE__ << " " << __LINE__ << " nSims " << nSims << " nSeeds " << see_q->size() << " nRecT "
0539 << trk_q->size() << std::endl;
0540
0541
0542 auto bestTkIdx = [&](std::vector<int> const& shs, std::vector<float> const& shfs, int rhIdx, HitType rhType) {
0543
0544 int ibest = -1;
0545 int shbest = -1;
0546 float hpbest = -1;
0547 float tpbest = -1;
0548 float hfbest = -1;
0549
0550 float maxfrac = -1;
0551 int ish = -1;
0552 int nshs = shs.size();
0553 for (auto const sh : shs) {
0554 ish++;
0555 auto tkidx = simhit_simTrkIdx->at(sh);
0556
0557 if (tkidx < 0)
0558 continue;
0559
0560 auto hpx = simhit_px->at(sh);
0561 auto hpy = simhit_py->at(sh);
0562 auto hpz = simhit_pz->at(sh);
0563 auto hp = sqrt(hpx * hpx + hpy * hpy + hpz * hpz);
0564
0565
0566 if (hp < 0.05f)
0567 continue;
0568
0569 auto tpx = sim_px->at(tkidx);
0570 auto tpy = sim_py->at(tkidx);
0571 auto tpz = sim_pz->at(tkidx);
0572 auto tp = sqrt(tpx * tpx + tpy * tpy + tpz * tpz);
0573
0574
0575 if (hp < 0.5 * tp)
0576 continue;
0577
0578
0579 if (maxfrac < hp / tp) {
0580 maxfrac = hp / tp;
0581 ibest = tkidx;
0582 shbest = sh;
0583 hpbest = hp;
0584 tpbest = tp;
0585 hfbest = shfs[ish];
0586 }
0587 }
0588
0589
0590
0591 if (nshs == 1 && ibest >= 0) {
0592 auto const& srhIdxV = simhit_hitIdx->at(shbest);
0593 auto const& srhTypeV = simhit_hitType->at(shbest);
0594 int ih = -1;
0595 for (auto itype : srhTypeV) {
0596 ih++;
0597 if (HitType(itype) == rhType && srhIdxV[ih] != rhIdx) {
0598 ibest = -1;
0599 break;
0600 }
0601 }
0602 }
0603
0604 if (ibest >= 0 && false) {
0605 std::cout << " best tkIdx " << ibest << " rh " << rhIdx << " for sh " << shbest << " out of " << shs.size()
0606 << " hp " << hpbest << " chF " << hfbest << " tp " << tpbest << " process "
0607 << simhit_process->at(shbest) << " particle " << simhit_particle->at(shbest) << std::endl;
0608 if (rhType == HitType::Strip) {
0609 std::cout << " sh " << simhit_x->at(shbest) << ", " << simhit_y->at(shbest) << ", " << simhit_z->at(shbest)
0610 << " rh " << str_x->at(rhIdx) << ", " << str_y->at(rhIdx) << ", " << str_z->at(rhIdx) << std::endl;
0611 }
0612 }
0613 return ibest;
0614 };
0615
0616 vector<Track>& simTracks_ = EE.simTracks_;
0617 vector<int> simTrackIdx_(sim_q->size(), -1);
0618 vector<int> seedSimIdx(see_q->size(), -1);
0619 for (unsigned int isim = 0; isim < sim_q->size(); ++isim) {
0620
0621 auto iVtx = sim_parentVtxIdx->at(isim);
0622 constexpr float largeValF = 9999.f;
0623 float sim_prodx = iVtx >= 0 ? simvtx_x->at(iVtx) : largeValF;
0624 float sim_prody = iVtx >= 0 ? simvtx_y->at(iVtx) : largeValF;
0625 float sim_prodz = iVtx >= 0 ? simvtx_z->at(iVtx) : largeValF;
0626
0627
0628 vector<int> const& trkIdxV = sim_trkIdx->at(isim);
0629
0630
0631
0632
0633 const int trkIdx = trkIdxV.empty() ? -1 : trkIdxV[0];
0634
0635 int nlay = 0;
0636 if (trkIdx >= 0) {
0637 std::vector<int> hitlay(nTotalLayers, 0);
0638 auto const& hits = trk_hitIdx->at(trkIdx);
0639 auto const& hitTypes = trk_hitType->at(trkIdx);
0640 auto nHits = hits.size();
0641 for (auto ihit = 0U; ihit < nHits; ++ihit) {
0642 auto ihIdx = hits[ihit];
0643 auto const ihType = HitType(hitTypes[ihit]);
0644
0645 switch (ihType) {
0646 case HitType::Pixel: {
0647 int ipix = ihIdx;
0648 if (ipix < 0)
0649 continue;
0650 int cmsswlay =
0651 lnc.convertLayerNumber(pix_det->at(ipix), pix_lay->at(ipix), useMatched, -1, pix_z->at(ipix) > 0);
0652 if (cmsswlay >= 0 && cmsswlay < static_cast<int>(nTotalLayers))
0653 hitlay[cmsswlay]++;
0654 break;
0655 }
0656 case HitType::Strip: {
0657 int istr = ihIdx;
0658 if (istr < 0)
0659 continue;
0660 int cmsswlay = lnc.convertLayerNumber(
0661 str_det->at(istr), str_lay->at(istr), useMatched, str_isStereo->at(istr), str_z->at(istr) > 0);
0662 if (cmsswlay >= 0 && cmsswlay < static_cast<int>(nTotalLayers))
0663 hitlay[cmsswlay]++;
0664 break;
0665 }
0666 case HitType::Glued: {
0667 if (useMatched) {
0668 int iglu = ihIdx;
0669 if (iglu < 0)
0670 continue;
0671 int cmsswlay =
0672 lnc.convertLayerNumber(glu_det->at(iglu), glu_lay->at(iglu), useMatched, -1, glu_z->at(iglu) > 0);
0673 if (cmsswlay >= 0 && cmsswlay < static_cast<int>(nTotalLayers))
0674 hitlay[cmsswlay]++;
0675 }
0676 break;
0677 }
0678 case HitType::Invalid:
0679 break;
0680 default:
0681 throw std::logic_error("Track type can not be handled");
0682 }
0683 }
0684 for (unsigned int i = 0; i < nTotalLayers; i++)
0685 if (hitlay[i] > 0)
0686 nlay++;
0687 }
0688
0689
0690
0691 SVector3 pos(sim_prodx, sim_prody, sim_prodz);
0692 SVector3 mom(sim_px->at(isim), sim_py->at(isim), sim_pz->at(isim));
0693 SMatrixSym66 err;
0694 err.At(0, 0) = sim_prodx * sim_prodx;
0695 err.At(1, 1) = sim_prody * sim_prody;
0696 err.At(2, 2) = sim_prodz * sim_prodz;
0697 err.At(3, 3) = sim_px->at(isim) * sim_px->at(isim);
0698 err.At(4, 4) = sim_py->at(isim) * sim_py->at(isim);
0699 err.At(5, 5) = sim_pz->at(isim) * sim_pz->at(isim);
0700 TrackState state(sim_q->at(isim), pos, mom, err);
0701 state.convertFromCartesianToCCS();
0702
0703
0704 Track track(state, float(nlay), simTracks_.size(), 0, nullptr);
0705 if (sim_bunchCrossing->at(isim) == 0) {
0706 if (sim_event->at(isim) == 0)
0707 track.setProdType(Track::ProdType::Signal);
0708 else
0709 track.setProdType(Track::ProdType::InTimePU);
0710 } else {
0711 track.setProdType(Track::ProdType::OutOfTimePU);
0712 }
0713 if (trkIdx >= 0) {
0714 int seedIdx = trk_seedIdx->at(trkIdx);
0715
0716 seedSimIdx[seedIdx] = simTracks_.size();
0717 }
0718 if (cleanSimTracks) {
0719 if (sim_nValid->at(isim) < cleanSimTrack_minSimHits)
0720 continue;
0721 if (cleanSimTrack_minRecHits > 0) {
0722 int nRecToSimHit = 0;
0723 for (unsigned int ipix = 0; ipix < pix_lay->size() && nRecToSimHit < cleanSimTrack_minRecHits; ++ipix) {
0724 int ilay = -1;
0725 ilay = lnc.convertLayerNumber(pix_det->at(ipix), pix_lay->at(ipix), useMatched, -1, pix_z->at(ipix) > 0);
0726 if (ilay < 0)
0727 continue;
0728 int simTkIdxNt = bestTkIdx(pix_simHitIdx->at(ipix), pix_chargeFraction->at(ipix), ipix, HitType::Pixel);
0729 if (simTkIdxNt >= 0)
0730 nRecToSimHit++;
0731 }
0732 if (useMatched) {
0733 for (unsigned int iglu = 0; iglu < glu_lay->size() && nRecToSimHit < cleanSimTrack_minRecHits; ++iglu) {
0734 if (glu_isBarrel->at(iglu) == 0)
0735 continue;
0736 int igluMono = glu_monoIdx->at(iglu);
0737 int simTkIdxNt =
0738 bestTkIdx(str_simHitIdx->at(igluMono), str_chargeFraction->at(igluMono), igluMono, HitType::Strip);
0739 if (simTkIdxNt >= 0)
0740 nRecToSimHit++;
0741 }
0742 }
0743 for (unsigned int istr = 0; istr < str_lay->size() && nRecToSimHit < cleanSimTrack_minRecHits; ++istr) {
0744 int ilay = -1;
0745 ilay = lnc.convertLayerNumber(
0746 str_det->at(istr), str_lay->at(istr), useMatched, str_isStereo->at(istr), str_z->at(istr) > 0);
0747 if (useMatched && str_isBarrel->at(istr) == 1 && str_isStereo->at(istr))
0748 continue;
0749 if (ilay == -1)
0750 continue;
0751 int simTkIdxNt = bestTkIdx(str_simHitIdx->at(istr), str_chargeFraction->at(istr), istr, HitType::Strip);
0752 if (simTkIdxNt >= 0)
0753 nRecToSimHit++;
0754 }
0755 if (nRecToSimHit < cleanSimTrack_minRecHits)
0756 continue;
0757 }
0758 }
0759
0760 simTrackIdx_[isim] = simTracks_.size();
0761 simTracks_.push_back(track);
0762 }
0763
0764 if (simTracks_.empty() and not writeAllEvents)
0765 continue;
0766
0767 vector<Track>& seedTracks_ = EE.seedTracks_;
0768 vector<vector<int>> pixHitSeedIdx(pix_lay->size());
0769 vector<vector<int>> strHitSeedIdx(str_lay->size());
0770 vector<vector<int>> gluHitSeedIdx(glu_lay->size());
0771 for (unsigned int is = 0; is < see_q->size(); ++is) {
0772 auto isAlgo = TrackAlgorithm(see_algo->at(is));
0773 if (not allSeeds)
0774 if (isAlgo != TrackAlgorithm::initialStep && isAlgo != TrackAlgorithm::hltIter0)
0775 continue;
0776
0777 SVector3 pos = SVector3(see_stateTrajGlbX->at(is), see_stateTrajGlbY->at(is), see_stateTrajGlbZ->at(is));
0778 SVector3 mom = SVector3(see_stateTrajGlbPx->at(is), see_stateTrajGlbPy->at(is), see_stateTrajGlbPz->at(is));
0779 SMatrixSym66 err;
0780 if (hasCartCov) {
0781 err.At(0, 0) = see_stateCcov00->at(is);
0782 err.At(0, 1) = see_stateCcov01->at(is);
0783 err.At(0, 2) = see_stateCcov02->at(is);
0784 err.At(0, 3) = see_stateCcov03->at(is);
0785 err.At(0, 4) = see_stateCcov04->at(is);
0786 err.At(0, 5) = see_stateCcov05->at(is);
0787 err.At(1, 1) = see_stateCcov11->at(is);
0788 err.At(1, 2) = see_stateCcov12->at(is);
0789 err.At(1, 3) = see_stateCcov13->at(is);
0790 err.At(1, 4) = see_stateCcov14->at(is);
0791 err.At(1, 5) = see_stateCcov15->at(is);
0792 err.At(2, 2) = see_stateCcov22->at(is);
0793 err.At(2, 3) = see_stateCcov23->at(is);
0794 err.At(2, 4) = see_stateCcov24->at(is);
0795 err.At(2, 5) = see_stateCcov25->at(is);
0796 err.At(3, 3) = see_stateCcov33->at(is);
0797 err.At(3, 4) = see_stateCcov34->at(is);
0798 err.At(3, 5) = see_stateCcov35->at(is);
0799 err.At(4, 4) = see_stateCcov44->at(is);
0800 err.At(4, 5) = see_stateCcov45->at(is);
0801 err.At(5, 5) = see_stateCcov55->at(is);
0802 } else {
0803 auto const& vCov = see_stateCurvCov->at(is);
0804 assert(vCov.size() == 15);
0805 auto vCovP = vCov.begin();
0806 for (int i = 0; i < 5; ++i)
0807 for (int j = 0; j <= i; ++j)
0808 err.At(i, j) = *(vCovP++);
0809 }
0810 TrackState state(see_q->at(is), pos, mom, err);
0811 if (hasCartCov)
0812 state.convertFromCartesianToCCS();
0813 else
0814 state.convertFromGlbCurvilinearToCCS();
0815 Track track(state, 0, seedSimIdx[is], 0, nullptr);
0816 track.setAlgorithm(isAlgo);
0817 auto const& shTypes = see_hitType->at(is);
0818 auto const& shIdxs = see_hitIdx->at(is);
0819 if (not allSeeds)
0820 if (!((isAlgo == TrackAlgorithm::initialStep || isAlgo == TrackAlgorithm::hltIter0) &&
0821 std::count(shTypes.begin(), shTypes.end(), int(HitType::Pixel)) >= 3))
0822 continue;
0823 for (unsigned int ip = 0; ip < shTypes.size(); ip++) {
0824 unsigned int hidx = shIdxs[ip];
0825 switch (HitType(shTypes[ip])) {
0826 case HitType::Pixel: {
0827 pixHitSeedIdx[hidx].push_back(seedTracks_.size());
0828 break;
0829 }
0830 case HitType::Strip: {
0831 strHitSeedIdx[hidx].push_back(seedTracks_.size());
0832 break;
0833 }
0834 case HitType::Glued: {
0835 if (not useMatched) {
0836
0837 int uidx = glu_monoIdx->at(hidx);
0838 strHitSeedIdx[uidx].push_back(seedTracks_.size());
0839 uidx = glu_stereoIdx->at(hidx);
0840 strHitSeedIdx[uidx].push_back(seedTracks_.size());
0841 } else {
0842 gluHitSeedIdx[hidx].push_back(seedTracks_.size());
0843 }
0844 break;
0845 }
0846 case HitType::Invalid:
0847 break;
0848 default:
0849 throw std::logic_error("Track hit type can not be handled");
0850 }
0851 }
0852 seedTracks_.push_back(track);
0853 }
0854
0855 if (seedTracks_.empty() and not writeAllEvents)
0856 continue;
0857
0858 vector<Track>& cmsswTracks_ = EE.cmsswTracks_;
0859 vector<vector<int>> pixHitRecIdx(pix_lay->size());
0860 vector<vector<int>> strHitRecIdx(str_lay->size());
0861 vector<vector<int>> gluHitRecIdx(glu_lay->size());
0862 for (unsigned int ir = 0; ir < trk_q->size(); ++ir) {
0863
0864 if (not allSeeds)
0865 if ((trk_algoMask->at(ir) & ((1 << int(TrackAlgorithm::initialStep)) | (1 << int(TrackAlgorithm::hltIter0)))) ==
0866 0) {
0867 if (verbosity > 1) {
0868 std::cout << "track " << ir << " failed algo selection for " << int(TrackAlgorithm::initialStep)
0869 << ": mask " << trk_algoMask->at(ir) << " origAlgo " << trk_originalAlgo->at(ir) << " algo "
0870 << trk_algo->at(ir) << std::endl;
0871 }
0872 continue;
0873 }
0874
0875 SMatrixSym66 err;
0876
0877
0878
0879
0880
0881
0882
0883
0884
0885
0886
0887
0888 float pt = trk_pt->at(ir);
0889 float pz = trk_pz->at(ir);
0890 float p2 = pt * pt + pz * pz;
0891 float phi = trk_phi->at(ir);
0892 float sP = sin(phi);
0893 float cP = cos(phi);
0894 float dxyErr2 = trk_dxyErr->at(ir);
0895 dxyErr2 *= dxyErr2;
0896 float dzErr2 = trk_dzErr->at(ir);
0897 dzErr2 *= dzErr2;
0898 float dzErrF2 = trk_dzErr->at(ir) * (pt * pz / p2);
0899 dzErr2 *= dzErr2;
0900 err.At(0, 0) = dxyErr2 * sP * sP + dzErrF2 * cP * cP;
0901 err.At(0, 1) = -dxyErr2 * cP * sP + dzErrF2 * cP * sP;
0902 err.At(1, 1) = dxyErr2 * cP * cP + dzErrF2 * sP * sP;
0903 err.At(0, 2) = -dzErrF2 * cP * pt / pz;
0904 err.At(1, 2) = -dzErrF2 * sP * pt / pz;
0905 err.At(2, 2) = dzErr2 * std::pow((pt * pt / p2), 2);
0906 err.At(3, 3) = std::pow(trk_ptErr->at(ir) / pt / pt, 2);
0907 err.At(4, 4) = std::pow(trk_phiErr->at(ir), 2);
0908 err.At(5, 5) = std::pow(trk_lambdaErr->at(ir), 2);
0909 SVector3 pos = SVector3(trk_refpoint_x->at(ir), trk_refpoint_y->at(ir), trk_refpoint_z->at(ir));
0910 SVector3 mom = SVector3(1.f / pt, phi, M_PI_2 - trk_lambda->at(ir));
0911 TrackState state(trk_q->at(ir), pos, mom, err);
0912 Track track(state, trk_nChi2->at(ir), trk_seedIdx->at(ir), 0, nullptr);
0913 track.setAlgorithm(TrackAlgorithm(trk_originalAlgo->at(ir)));
0914 auto const& hTypes = trk_hitType->at(ir);
0915 auto const& hIdxs = trk_hitIdx->at(ir);
0916 for (unsigned int ip = 0; ip < hTypes.size(); ip++) {
0917 unsigned int hidx = hIdxs[ip];
0918 switch (HitType(hTypes[ip])) {
0919 case HitType::Pixel: {
0920
0921 pixHitRecIdx[hidx].push_back(cmsswTracks_.size());
0922 break;
0923 }
0924 case HitType::Strip: {
0925
0926 strHitRecIdx[hidx].push_back(cmsswTracks_.size());
0927 break;
0928 }
0929 case HitType::Glued: {
0930 if (not useMatched)
0931 throw std::logic_error("Tracks have glued hits, but matchedHit load is not configured");
0932
0933 gluHitRecIdx[hidx].push_back(cmsswTracks_.size());
0934 break;
0935 }
0936 case HitType::Invalid:
0937 break;
0938 default:
0939 throw std::logic_error("Track hit type can not be handled");
0940 }
0941 }
0942 cmsswTracks_.push_back(track);
0943 }
0944
0945 vector<vector<Hit>>& layerHits_ = EE.layerHits_;
0946 vector<vector<uint64_t>>& layerHitMasks_ = EE.layerHitMasks_;
0947 vector<MCHitInfo>& simHitsInfo_ = EE.simHitsInfo_;
0948 int totHits = 0;
0949 layerHits_.resize(nTotalLayers);
0950 layerHitMasks_.resize(nTotalLayers);
0951 for (unsigned int ipix = 0; ipix < pix_lay->size(); ++ipix) {
0952 int ilay = -1;
0953 ilay = lnc.convertLayerNumber(pix_det->at(ipix), pix_lay->at(ipix), useMatched, -1, pix_z->at(ipix) > 0);
0954 if (ilay < 0)
0955 continue;
0956
0957 unsigned int imoduleid = tkinfo[ilay].short_id(pix_detId->at(ipix));
0958
0959 int simTkIdxNt = bestTkIdx(pix_simHitIdx->at(ipix), pix_chargeFraction->at(ipix), ipix, HitType::Pixel);
0960 int simTkIdx = simTkIdxNt >= 0 ? simTrackIdx_[simTkIdxNt] : -1;
0961
0962
0963 SVector3 pos(pix_x->at(ipix), pix_y->at(ipix), pix_z->at(ipix));
0964 SMatrixSym33 err;
0965 err.At(0, 0) = pix_xx->at(ipix);
0966 err.At(1, 1) = pix_yy->at(ipix);
0967 err.At(2, 2) = pix_zz->at(ipix);
0968 err.At(0, 1) = pix_xy->at(ipix);
0969 err.At(0, 2) = pix_zx->at(ipix);
0970 err.At(1, 2) = pix_yz->at(ipix);
0971 if (simTkIdx >= 0) {
0972 simTracks_[simTkIdx].addHitIdx(layerHits_[ilay].size(), ilay, 0);
0973 }
0974 for (unsigned int is = 0; is < pixHitSeedIdx[ipix].size(); is++) {
0975
0976 seedTracks_[pixHitSeedIdx[ipix][is]].addHitIdx(layerHits_[ilay].size(), ilay, 0);
0977 }
0978 for (unsigned int ir = 0; ir < pixHitRecIdx[ipix].size(); ir++) {
0979
0980 cmsswTracks_[pixHitRecIdx[ipix][ir]].addHitIdx(layerHits_[ilay].size(), ilay, 0);
0981 }
0982 Hit hit(pos, err, totHits);
0983 hit.setupAsPixel(imoduleid, pix_csize_row->at(ipix), pix_csize_col->at(ipix));
0984 layerHits_[ilay].push_back(hit);
0985 if (writeHitIterMasks)
0986 layerHitMasks_[ilay].push_back(pix_usedMask->at(ipix));
0987 MCHitInfo hitInfo(simTkIdx, ilay, layerHits_[ilay].size() - 1, totHits);
0988 simHitsInfo_.push_back(hitInfo);
0989 totHits++;
0990 }
0991
0992 if (useMatched) {
0993 for (unsigned int iglu = 0; iglu < glu_lay->size(); ++iglu) {
0994 if (glu_isBarrel->at(iglu) == 0)
0995 continue;
0996 int igluMono = glu_monoIdx->at(iglu);
0997 int simTkIdxNt =
0998 bestTkIdx(str_simHitIdx->at(igluMono), str_chargeFraction->at(igluMono), igluMono, HitType::Strip);
0999 int simTkIdx = simTkIdxNt >= 0 ? simTrackIdx_[simTkIdxNt] : -1;
1000
1001 int ilay = lnc.convertLayerNumber(glu_det->at(iglu), glu_lay->at(iglu), useMatched, -1, glu_z->at(iglu) > 0);
1002
1003 SVector3 pos(glu_x->at(iglu), glu_y->at(iglu), glu_z->at(iglu));
1004 SMatrixSym33 err;
1005 err.At(0, 0) = glu_xx->at(iglu);
1006 err.At(1, 1) = glu_yy->at(iglu);
1007 err.At(2, 2) = glu_zz->at(iglu);
1008 err.At(0, 1) = glu_xy->at(iglu);
1009 err.At(0, 2) = glu_zx->at(iglu);
1010 err.At(1, 2) = glu_yz->at(iglu);
1011 if (simTkIdx >= 0) {
1012 simTracks_[simTkIdx].addHitIdx(layerHits_[ilay].size(), ilay, 0);
1013 }
1014 for (unsigned int ir = 0; ir < gluHitSeedIdx[iglu].size(); ir++) {
1015
1016 seedTracks_[gluHitSeedIdx[iglu][ir]].addHitIdx(layerHits_[ilay].size(), ilay, 0);
1017 }
1018 for (unsigned int ir = 0; ir < gluHitRecIdx[iglu].size(); ir++) {
1019
1020 cmsswTracks_[gluHitRecIdx[iglu][ir]].addHitIdx(layerHits_[ilay].size(), ilay, 0);
1021 }
1022
1023
1024
1025 assert(false && "Implement module-ids, cluster adc and spans for matched hits!");
1026
1027 Hit hit(pos, err, totHits);
1028 layerHits_[ilay].push_back(hit);
1029 MCHitInfo hitInfo(simTkIdx, ilay, layerHits_[ilay].size() - 1, totHits);
1030 simHitsInfo_.push_back(hitInfo);
1031 totHits++;
1032 }
1033 }
1034
1035 vector<int> strIdx;
1036 strIdx.resize(str_lay->size());
1037 for (unsigned int istr = 0; istr < str_lay->size(); ++istr) {
1038 int ilay = -1;
1039 ilay = lnc.convertLayerNumber(
1040 str_det->at(istr), str_lay->at(istr), useMatched, str_isStereo->at(istr), str_z->at(istr) > 0);
1041 if (useMatched && str_isBarrel->at(istr) == 1 && str_isStereo->at(istr))
1042 continue;
1043 if (ilay == -1)
1044 continue;
1045
1046 unsigned int imoduleid = tkinfo[ilay].short_id(str_detId->at(istr));
1047
1048 int simTkIdxNt = bestTkIdx(str_simHitIdx->at(istr), str_chargeFraction->at(istr), istr, HitType::Strip);
1049 int simTkIdx = simTkIdxNt >= 0 ? simTrackIdx_[simTkIdxNt] : -1;
1050
1051 bool passCCC = applyCCC ? (str_chargePerCM->at(istr) > cutValueCCC) : true;
1052
1053
1054 SVector3 pos(str_x->at(istr), str_y->at(istr), str_z->at(istr));
1055 SMatrixSym33 err;
1056 err.At(0, 0) = str_xx->at(istr);
1057 err.At(1, 1) = str_yy->at(istr);
1058 err.At(2, 2) = str_zz->at(istr);
1059 err.At(0, 1) = str_xy->at(istr);
1060 err.At(0, 2) = str_zx->at(istr);
1061 err.At(1, 2) = str_yz->at(istr);
1062 if (simTkIdx >= 0) {
1063 if (passCCC)
1064 simTracks_[simTkIdx].addHitIdx(layerHits_[ilay].size(), ilay, 0);
1065 else
1066 simTracks_[simTkIdx].addHitIdx(-9, ilay, 0);
1067 }
1068 for (unsigned int ir = 0; ir < strHitSeedIdx[istr].size(); ir++) {
1069
1070 if (passCCC)
1071 seedTracks_[strHitSeedIdx[istr][ir]].addHitIdx(layerHits_[ilay].size(), ilay, 0);
1072 else
1073 seedTracks_[strHitSeedIdx[istr][ir]].addHitIdx(-9, ilay, 0);
1074 }
1075 for (unsigned int ir = 0; ir < strHitRecIdx[istr].size(); ir++) {
1076
1077 if (passCCC)
1078 cmsswTracks_[strHitRecIdx[istr][ir]].addHitIdx(layerHits_[ilay].size(), ilay, 0);
1079 else
1080 cmsswTracks_[strHitRecIdx[istr][ir]].addHitIdx(-9, ilay, 0);
1081 }
1082 if (passCCC) {
1083 Hit hit(pos, err, totHits);
1084 hit.setupAsStrip(imoduleid, str_chargePerCM->at(istr), str_csize->at(istr));
1085 layerHits_[ilay].push_back(hit);
1086 if (writeHitIterMasks)
1087 layerHitMasks_[ilay].push_back(str_usedMask->at(istr));
1088 MCHitInfo hitInfo(simTkIdx, ilay, layerHits_[ilay].size() - 1, totHits);
1089 simHitsInfo_.push_back(hitInfo);
1090 totHits++;
1091 }
1092 }
1093
1094
1095 nstot += seedTracks_.size();
1096 for (unsigned int il = 0; il < layerHits_.size(); ++il) {
1097 int nh = layerHits_[il].size();
1098 nhitstot[il] += nh;
1099 }
1100
1101 if (verbosity > 0) {
1102 int nt = simTracks_.size();
1103
1104 int nl = layerHits_.size();
1105
1106 int nm = simHitsInfo_.size();
1107
1108 int ns = seedTracks_.size();
1109
1110 int nr = cmsswTracks_.size();
1111
1112 printf("number of simTracks %i\n", nt);
1113 printf("number of layerHits %i\n", nl);
1114 printf("number of simHitsInfo %i\n", nm);
1115 printf("number of seedTracks %i\n", ns);
1116 printf("number of recTracks %i\n", nr);
1117
1118 if (verbosity > 1) {
1119 printf("\n");
1120 for (int il = 0; il < nl; ++il) {
1121 int nh = layerHits_[il].size();
1122 for (int ih = 0; ih < nh; ++ih) {
1123 printf("lay=%i idx=%i mcid=%i x=(%6.3f, %6.3f, %6.3f) r=%6.3f mask=0x%lx\n",
1124 il + 1,
1125 ih,
1126 layerHits_[il][ih].mcHitID(),
1127 layerHits_[il][ih].x(),
1128 layerHits_[il][ih].y(),
1129 layerHits_[il][ih].z(),
1130 sqrt(pow(layerHits_[il][ih].x(), 2) + pow(layerHits_[il][ih].y(), 2)),
1131 writeHitIterMasks ? layerHitMasks_[il][ih] : 0);
1132 }
1133 }
1134
1135 for (int i = 0; i < nt; ++i) {
1136 float spt = sqrt(pow(simTracks_[i].px(), 2) + pow(simTracks_[i].py(), 2));
1137 printf(
1138 "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",
1139 i,
1140 simTracks_[i].charge(),
1141 simTracks_[i].px(),
1142 simTracks_[i].py(),
1143 simTracks_[i].pz(),
1144 simTracks_[i].x(),
1145 simTracks_[i].y(),
1146 simTracks_[i].z(),
1147 spt,
1148 simTracks_[i].nTotalHits(),
1149 simTracks_[i].nFoundHits());
1150 int nh = simTracks_[i].nTotalHits();
1151 for (int ih = 0; ih < nh; ++ih) {
1152 int hidx = simTracks_[i].getHitIdx(ih);
1153 int hlay = simTracks_[i].getHitLyr(ih);
1154 float hx = layerHits_[hlay][hidx].x();
1155 float hy = layerHits_[hlay][hidx].y();
1156 float hz = layerHits_[hlay][hidx].z();
1157 printf("track #%4i hit #%2i idx=%4i lay=%2i x=(% 8.3f, % 8.3f, % 8.3f) r=%8.3f\n",
1158 i,
1159 ih,
1160 hidx,
1161 hlay,
1162 hx,
1163 hy,
1164 hz,
1165 sqrt(hx * hx + hy * hy));
1166 }
1167 }
1168
1169 for (int i = 0; i < ns; ++i) {
1170 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",
1171 i,
1172 seedTracks_[i].label(),
1173 (int)seedTracks_[i].algorithm(),
1174 seedTracks_[i].charge(),
1175 seedTracks_[i].pT(),
1176 seedTracks_[i].px(),
1177 seedTracks_[i].py(),
1178 seedTracks_[i].pz(),
1179 seedTracks_[i].x(),
1180 seedTracks_[i].y(),
1181 seedTracks_[i].z());
1182 int nh = seedTracks_[i].nTotalHits();
1183 for (int ih = 0; ih < nh; ++ih)
1184 printf("seed #%i hit #%i idx=%i\n", i, ih, seedTracks_[i].getHitIdx(ih));
1185 }
1186
1187 if (writeRecTracks) {
1188 for (int i = 0; i < nr; ++i) {
1189 float spt = sqrt(pow(cmsswTracks_[i].px(), 2) + pow(cmsswTracks_[i].py(), 2));
1190 printf(
1191 "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) "
1192 "pT=%7.4f nTotal=%i nFound=%i \n",
1193 i,
1194 cmsswTracks_[i].label(),
1195 (int)cmsswTracks_[i].algorithm(),
1196 cmsswTracks_[i].chi2(),
1197 cmsswTracks_[i].charge(),
1198 cmsswTracks_[i].px(),
1199 cmsswTracks_[i].py(),
1200 cmsswTracks_[i].pz(),
1201 cmsswTracks_[i].x(),
1202 cmsswTracks_[i].y(),
1203 cmsswTracks_[i].z(),
1204 spt,
1205 cmsswTracks_[i].nTotalHits(),
1206 cmsswTracks_[i].nFoundHits());
1207 int nh = cmsswTracks_[i].nTotalHits();
1208 for (int ih = 0; ih < nh; ++ih) {
1209 int hidx = cmsswTracks_[i].getHitIdx(ih);
1210 int hlay = cmsswTracks_[i].getHitLyr(ih);
1211 float hx = layerHits_[hlay][hidx].x();
1212 float hy = layerHits_[hlay][hidx].y();
1213 float hz = layerHits_[hlay][hidx].z();
1214 printf("track #%4i hit #%2i idx=%4i lay=%2i x=(% 8.3f, % 8.3f, % 8.3f) r=%8.3f\n",
1215 i,
1216 ih,
1217 hidx,
1218 hlay,
1219 hx,
1220 hy,
1221 hz,
1222 sqrt(hx * hx + hy * hy));
1223 }
1224 }
1225 }
1226
1227 }
1228 }
1229 EE.write_out(data_file);
1230
1231 savedEvents++;
1232 printf("end of event %lli\n", savedEvents);
1233 }
1234
1235 data_file.CloseWrite(savedEvents);
1236 printf("\nSaved %lli events\n\n", savedEvents);
1237
1238 printf("Average number of seeds per event %f\n", float(nstot) / float(savedEvents));
1239 for (unsigned int il = 0; il < nhitstot.size(); ++il)
1240 printf("Average number of hits in layer %3i = %7.2f\n",
1241 il,
1242 float(nhitstot[il]) / float(savedEvents));
1243
1244 printf("Out of %i hits, %i failed the cut", numTotalStr, numFailCCC);
1245
1246
1247
1248 printf("\n\n================================================================\n");
1249 printf("=== Max module id for %u layers\n", nTotalLayers);
1250 printf("================================================================\n");
1251 for (unsigned int ii = 0; ii < nTotalLayers; ++ii) {
1252 printf("Layer%2d : %d\n", ii, tkinfo[ii].n_modules());
1253 }
1254 }