Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:22:27 

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 //check if this is the same as in the release
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       "  --help                    print help and exit\n"
0049       "  --input          <str>    input file\n"
0050       "  --output         <str>    output file\n"
0051       "  --geo            <file>   binary TrackerInfo geometry (def: CMS-phase1.bin)\n"
0052       "  --verbosity      <num>    print details (0 quiet, 1 print counts, 2 print all; def: 0)\n"
0053       "  --maxevt         <num>    maxevt events to write (-1 for everything in the file def: -1)\n"
0054       "  --clean-sim-tracks        apply sim track cleaning (def: no cleaning)\n"
0055       "  --write-all-events        write all events (def: skip events with 0 simtracks or seeds)\n"
0056       "  --write-rec-tracks        write rec tracks (def: not written)\n"
0057       "  --apply-ccc               apply cluster charge cut to strip hits (def: false)\n"
0058       "  --all-seeds               write all seeds from the input file, not only initialStep and hltIter0 (def: "
0059       "false)\n",
0060       av0);
0061 }
0062 
0063 int main(int argc, char* argv[]) {
0064   bool haveInput = false;
0065   std::string inputFileName;
0066   bool haveOutput = false;
0067   std::string outputFileName;
0068   std::string geoFileName("CMS-phase1.bin");
0069 
0070   bool cleanSimTracks = false;
0071   bool writeAllEvents = false;
0072   bool writeRecTracks = false;
0073   bool writeHitIterMasks = false;
0074   bool applyCCC = false;
0075   bool allSeeds = false;
0076 
0077   int verbosity = 0;
0078   long long maxevt = -1;
0079 
0080   int cutValueCCC = 1620;  //Nominal value (from first iteration of CMSSW) is 1620
0081 
0082   lStr_t mArgs;
0083   for (int i = 1; i < argc; ++i) {
0084     mArgs.push_back(argv[i]);
0085   }
0086 
0087   lStr_i i = mArgs.begin();
0088   while (i != mArgs.end()) {
0089     lStr_i start = i;
0090 
0091     if (*i == "-h" || *i == "-help" || *i == "--help") {
0092       printHelp(argv[0]);
0093       exit(0);
0094     } else if (*i == "--input") {
0095       next_arg_or_die(mArgs, i);
0096       inputFileName = *i;
0097       haveInput = true;
0098     } else if (*i == "--output") {
0099       next_arg_or_die(mArgs, i);
0100       outputFileName = *i;
0101       haveOutput = true;
0102     } else if (*i == "--geo") {
0103       next_arg_or_die(mArgs, i);
0104       geoFileName = *i;
0105     } else if (*i == "--verbosity") {
0106       next_arg_or_die(mArgs, i);
0107       verbosity = std::atoi(i->c_str());
0108     } else if (*i == "--maxevt") {
0109       next_arg_or_die(mArgs, i);
0110       maxevt = std::atoi(i->c_str());
0111     } else if (*i == "--clean-sim-tracks") {
0112       cleanSimTracks = true;
0113     } else if (*i == "--write-all-events") {
0114       writeAllEvents = true;
0115     } else if (*i == "--write-rec-tracks") {
0116       writeRecTracks = true;
0117     } else if (*i == "--write-hit-iter-masks") {
0118       writeHitIterMasks = true;
0119     } else if (*i == "--apply-ccc") {
0120       applyCCC = true;
0121       if (next_arg_option(mArgs, i)) {
0122         cutValueCCC = std::atoi(i->c_str());
0123       }
0124     } else if (*i == "--all-seeds") {
0125       allSeeds = true;
0126     } else {
0127       fprintf(stderr, "Error: Unknown option/argument '%s'.\n", i->c_str());
0128       printHelp(argv[0]);
0129       exit(1);
0130     }
0131     mArgs.erase(start, ++i);
0132   }  //while arguments
0133 
0134   if (not haveOutput or not haveInput) {
0135     fprintf(stderr, "Error: both input and output are required\n");
0136     printHelp(argv[0]);
0137     exit(1);
0138   }
0139 
0140   using namespace std;
0141 
0142   TFile* f = TFile::Open(inputFileName.c_str());
0143   if (f == 0) {
0144     fprintf(stderr, "Failed opening input root file '%s'\n", inputFileName.c_str());
0145     exit(1);
0146   }
0147 
0148   TTree* t = (TTree*)f->Get("trackingNtuple/tree");
0149 
0150   bool hasPh2hits = t->GetBranch("ph2_isLower") != nullptr;
0151   TrackerInfo tkinfo;
0152   tkinfo.read_bin_file(geoFileName);
0153   LayerNumberConverter lnc(hasPh2hits ? TkLayout::phase2 : TkLayout::phase1);
0154   const unsigned int nTotalLayers = lnc.nLayers();
0155   assert(nTotalLayers == (unsigned int)tkinfo.n_layers());
0156 
0157   int nstot = 0;
0158   std::vector<int> nhitstot(nTotalLayers, 0);
0159 
0160   unsigned long long event;
0161   t->SetBranchAddress("event", &event);
0162 
0163   //sim tracks
0164   std::vector<float>* sim_eta = 0;
0165   std::vector<float>* sim_px = 0;
0166   std::vector<float>* sim_py = 0;
0167   std::vector<float>* sim_pz = 0;
0168   std::vector<int>* sim_parentVtxIdx = 0;
0169   std::vector<int>* sim_q = 0;
0170   std::vector<int>* sim_event = 0;
0171   std::vector<int>* sim_bunchCrossing = 0;
0172   std::vector<int>* sim_nValid = 0;  //simHit count, actually
0173   t->SetBranchAddress("sim_eta", &sim_eta);
0174   t->SetBranchAddress("sim_px", &sim_px);
0175   t->SetBranchAddress("sim_py", &sim_py);
0176   t->SetBranchAddress("sim_pz", &sim_pz);
0177   t->SetBranchAddress("sim_parentVtxIdx", &sim_parentVtxIdx);
0178   t->SetBranchAddress("sim_q", &sim_q);
0179   t->SetBranchAddress("sim_event", &sim_event);
0180   t->SetBranchAddress("sim_bunchCrossing", &sim_bunchCrossing);
0181   t->SetBranchAddress("sim_nValid", &sim_nValid);
0182 
0183   std::vector<vector<int>>* sim_trkIdx = 0;
0184   t->SetBranchAddress("sim_trkIdx", &sim_trkIdx);
0185 
0186   //simvtx
0187   std::vector<float>* simvtx_x = 0;
0188   std::vector<float>* simvtx_y = 0;
0189   std::vector<float>* simvtx_z = 0;
0190   t->SetBranchAddress("simvtx_x", &simvtx_x);
0191   t->SetBranchAddress("simvtx_y", &simvtx_y);
0192   t->SetBranchAddress("simvtx_z", &simvtx_z);
0193 
0194   //simhit
0195   std::vector<short>* simhit_process = 0;
0196   std::vector<int>* simhit_particle = 0;
0197   std::vector<int>* simhit_simTrkIdx = 0;
0198   std::vector<float>* simhit_x = 0;
0199   std::vector<float>* simhit_y = 0;
0200   std::vector<float>* simhit_z = 0;
0201   std::vector<float>* simhit_px = 0;
0202   std::vector<float>* simhit_py = 0;
0203   std::vector<float>* simhit_pz = 0;
0204   t->SetBranchAddress("simhit_process", &simhit_process);
0205   t->SetBranchAddress("simhit_particle", &simhit_particle);
0206   t->SetBranchAddress("simhit_simTrkIdx", &simhit_simTrkIdx);
0207   t->SetBranchAddress("simhit_x", &simhit_x);
0208   t->SetBranchAddress("simhit_y", &simhit_y);
0209   t->SetBranchAddress("simhit_z", &simhit_z);
0210   t->SetBranchAddress("simhit_px", &simhit_px);
0211   t->SetBranchAddress("simhit_py", &simhit_py);
0212   t->SetBranchAddress("simhit_pz", &simhit_pz);
0213 
0214   std::vector<std::vector<int>>* simhit_hitIdx = 0;
0215   t->SetBranchAddress("simhit_hitIdx", &simhit_hitIdx);
0216   std::vector<std::vector<int>>* simhit_hitType = 0;
0217   t->SetBranchAddress("simhit_hitType", &simhit_hitType);
0218 
0219   //rec tracks
0220   std::vector<int>* trk_q = 0;
0221   std::vector<unsigned int>* trk_nValid = 0;
0222   std::vector<int>* trk_seedIdx = 0;
0223   std::vector<unsigned long long>* trk_algoMask = 0;
0224   std::vector<unsigned int>* trk_algo = 0;
0225   std::vector<unsigned int>* trk_originalAlgo = 0;
0226   std::vector<float>* trk_nChi2 = 0;
0227   std::vector<float>* trk_px = 0;
0228   std::vector<float>* trk_py = 0;
0229   std::vector<float>* trk_pz = 0;
0230   std::vector<float>* trk_pt = 0;
0231   std::vector<float>* trk_phi = 0;
0232   std::vector<float>* trk_lambda = 0;
0233   std::vector<float>* trk_refpoint_x = 0;
0234   std::vector<float>* trk_refpoint_y = 0;
0235   std::vector<float>* trk_refpoint_z = 0;
0236   std::vector<float>* trk_dxyErr = 0;
0237   std::vector<float>* trk_dzErr = 0;
0238   std::vector<float>* trk_ptErr = 0;
0239   std::vector<float>* trk_phiErr = 0;
0240   std::vector<float>* trk_lambdaErr = 0;
0241   t->SetBranchAddress("trk_q", &trk_q);
0242   t->SetBranchAddress("trk_nValid", &trk_nValid);
0243   t->SetBranchAddress("trk_seedIdx", &trk_seedIdx);
0244   t->SetBranchAddress("trk_algoMask", &trk_algoMask);
0245   t->SetBranchAddress("trk_algo", &trk_algo);
0246   t->SetBranchAddress("trk_originalAlgo", &trk_originalAlgo);
0247   t->SetBranchAddress("trk_nChi2", &trk_nChi2);
0248   t->SetBranchAddress("trk_px", &trk_px);
0249   t->SetBranchAddress("trk_py", &trk_py);
0250   t->SetBranchAddress("trk_pz", &trk_pz);
0251   t->SetBranchAddress("trk_pt", &trk_pt);
0252   t->SetBranchAddress("trk_phi", &trk_phi);
0253   t->SetBranchAddress("trk_lambda", &trk_lambda);
0254   t->SetBranchAddress("trk_refpoint_x", &trk_refpoint_x);
0255   t->SetBranchAddress("trk_refpoint_y", &trk_refpoint_y);
0256   t->SetBranchAddress("trk_refpoint_z", &trk_refpoint_z);
0257   t->SetBranchAddress("trk_dxyErr", &trk_dxyErr);
0258   t->SetBranchAddress("trk_dzErr", &trk_dzErr);
0259   t->SetBranchAddress("trk_ptErr", &trk_ptErr);
0260   t->SetBranchAddress("trk_phiErr", &trk_phiErr);
0261   t->SetBranchAddress("trk_lambdaErr", &trk_lambdaErr);
0262 
0263   std::vector<std::vector<int>>* trk_hitIdx = 0;
0264   t->SetBranchAddress("trk_hitIdx", &trk_hitIdx);
0265   std::vector<std::vector<int>>* trk_hitType = 0;
0266   t->SetBranchAddress("trk_hitType", &trk_hitType);
0267 
0268   //seeds
0269   std::vector<float>* see_stateTrajGlbX = 0;
0270   std::vector<float>* see_stateTrajGlbY = 0;
0271   std::vector<float>* see_stateTrajGlbZ = 0;
0272   std::vector<float>* see_stateTrajGlbPx = 0;
0273   std::vector<float>* see_stateTrajGlbPy = 0;
0274   std::vector<float>* see_stateTrajGlbPz = 0;
0275   std::vector<float>* see_eta = 0;  //PCA parameters
0276   std::vector<float>* see_pt = 0;   //PCA parameters
0277   std::vector<float>* see_stateCcov00 = 0;
0278   std::vector<float>* see_stateCcov01 = 0;
0279   std::vector<float>* see_stateCcov02 = 0;
0280   std::vector<float>* see_stateCcov03 = 0;
0281   std::vector<float>* see_stateCcov04 = 0;
0282   std::vector<float>* see_stateCcov05 = 0;
0283   std::vector<float>* see_stateCcov11 = 0;
0284   std::vector<float>* see_stateCcov12 = 0;
0285   std::vector<float>* see_stateCcov13 = 0;
0286   std::vector<float>* see_stateCcov14 = 0;
0287   std::vector<float>* see_stateCcov15 = 0;
0288   std::vector<float>* see_stateCcov22 = 0;
0289   std::vector<float>* see_stateCcov23 = 0;
0290   std::vector<float>* see_stateCcov24 = 0;
0291   std::vector<float>* see_stateCcov25 = 0;
0292   std::vector<float>* see_stateCcov33 = 0;
0293   std::vector<float>* see_stateCcov34 = 0;
0294   std::vector<float>* see_stateCcov35 = 0;
0295   std::vector<float>* see_stateCcov44 = 0;
0296   std::vector<float>* see_stateCcov45 = 0;
0297   std::vector<float>* see_stateCcov55 = 0;
0298   std::vector<std::vector<float>>* see_stateCurvCov = 0;
0299   std::vector<int>* see_q = 0;
0300   std::vector<unsigned int>* see_algo = 0;
0301   t->SetBranchAddress("see_stateTrajGlbX", &see_stateTrajGlbX);
0302   t->SetBranchAddress("see_stateTrajGlbY", &see_stateTrajGlbY);
0303   t->SetBranchAddress("see_stateTrajGlbZ", &see_stateTrajGlbZ);
0304   t->SetBranchAddress("see_stateTrajGlbPx", &see_stateTrajGlbPx);
0305   t->SetBranchAddress("see_stateTrajGlbPy", &see_stateTrajGlbPy);
0306   t->SetBranchAddress("see_stateTrajGlbPz", &see_stateTrajGlbPz);
0307   t->SetBranchAddress("see_eta", &see_eta);
0308   t->SetBranchAddress("see_pt", &see_pt);
0309 
0310   bool hasCartCov = t->GetBranch("see_stateCcov00") != nullptr;
0311   if (hasCartCov) {
0312     t->SetBranchAddress("see_stateCcov00", &see_stateCcov00);
0313     t->SetBranchAddress("see_stateCcov01", &see_stateCcov01);
0314     t->SetBranchAddress("see_stateCcov02", &see_stateCcov02);
0315     t->SetBranchAddress("see_stateCcov03", &see_stateCcov03);
0316     t->SetBranchAddress("see_stateCcov04", &see_stateCcov04);
0317     t->SetBranchAddress("see_stateCcov05", &see_stateCcov05);
0318     t->SetBranchAddress("see_stateCcov11", &see_stateCcov11);
0319     t->SetBranchAddress("see_stateCcov12", &see_stateCcov12);
0320     t->SetBranchAddress("see_stateCcov13", &see_stateCcov13);
0321     t->SetBranchAddress("see_stateCcov14", &see_stateCcov14);
0322     t->SetBranchAddress("see_stateCcov15", &see_stateCcov15);
0323     t->SetBranchAddress("see_stateCcov22", &see_stateCcov22);
0324     t->SetBranchAddress("see_stateCcov23", &see_stateCcov23);
0325     t->SetBranchAddress("see_stateCcov24", &see_stateCcov24);
0326     t->SetBranchAddress("see_stateCcov25", &see_stateCcov25);
0327     t->SetBranchAddress("see_stateCcov33", &see_stateCcov33);
0328     t->SetBranchAddress("see_stateCcov34", &see_stateCcov34);
0329     t->SetBranchAddress("see_stateCcov35", &see_stateCcov35);
0330     t->SetBranchAddress("see_stateCcov44", &see_stateCcov44);
0331     t->SetBranchAddress("see_stateCcov45", &see_stateCcov45);
0332     t->SetBranchAddress("see_stateCcov55", &see_stateCcov55);
0333   } else {
0334     t->SetBranchAddress("see_stateCurvCov", &see_stateCurvCov);
0335   }
0336   t->SetBranchAddress("see_q", &see_q);
0337   t->SetBranchAddress("see_algo", &see_algo);
0338 
0339   std::vector<std::vector<int>>* see_hitIdx = 0;
0340   t->SetBranchAddress("see_hitIdx", &see_hitIdx);
0341   std::vector<std::vector<int>>* see_hitType = 0;
0342   t->SetBranchAddress("see_hitType", &see_hitType);
0343 
0344   //pixel hits
0345   vector<unsigned short>* pix_det = 0;
0346   vector<unsigned short>* pix_lay = 0;
0347   vector<unsigned int>* pix_detId = 0;
0348   vector<float>* pix_x = 0;
0349   vector<float>* pix_y = 0;
0350   vector<float>* pix_z = 0;
0351   vector<float>* pix_xx = 0;
0352   vector<float>* pix_xy = 0;
0353   vector<float>* pix_yy = 0;
0354   vector<float>* pix_yz = 0;
0355   vector<float>* pix_zz = 0;
0356   vector<float>* pix_zx = 0;
0357   vector<int>* pix_csize_col = 0;
0358   vector<int>* pix_csize_row = 0;
0359   vector<uint64_t>* pix_usedMask = 0;
0360   //these were renamed in CMSSW_9_1_0: auto-detect
0361   bool has910_det_lay = t->GetBranch("pix_det") == nullptr;
0362   if (has910_det_lay) {
0363     t->SetBranchAddress("pix_subdet", &pix_det);
0364     t->SetBranchAddress("pix_layer", &pix_lay);
0365   } else {
0366     t->SetBranchAddress("pix_det", &pix_det);
0367     t->SetBranchAddress("pix_lay", &pix_lay);
0368   }
0369   t->SetBranchAddress("pix_detId", &pix_detId);
0370   t->SetBranchAddress("pix_x", &pix_x);
0371   t->SetBranchAddress("pix_y", &pix_y);
0372   t->SetBranchAddress("pix_z", &pix_z);
0373   t->SetBranchAddress("pix_xx", &pix_xx);
0374   t->SetBranchAddress("pix_xy", &pix_xy);
0375   t->SetBranchAddress("pix_yy", &pix_yy);
0376   t->SetBranchAddress("pix_yz", &pix_yz);
0377   t->SetBranchAddress("pix_zz", &pix_zz);
0378   t->SetBranchAddress("pix_zx", &pix_zx);
0379   t->SetBranchAddress("pix_clustSizeCol", &pix_csize_col);
0380   t->SetBranchAddress("pix_clustSizeRow", &pix_csize_row);
0381   if (writeHitIterMasks) {
0382     t->SetBranchAddress("pix_usedMask", &pix_usedMask);
0383   }
0384 
0385   vector<vector<int>>* pix_simHitIdx = 0;
0386   t->SetBranchAddress("pix_simHitIdx", &pix_simHitIdx);
0387   vector<vector<float>>* pix_chargeFraction = 0;
0388   t->SetBranchAddress("pix_chargeFraction", &pix_chargeFraction);
0389 
0390   //strip hits
0391   vector<short>* glu_isBarrel = 0;
0392   vector<unsigned int>* glu_det = 0;
0393   vector<unsigned int>* glu_lay = 0;
0394   vector<unsigned int>* glu_detId = 0;
0395   vector<int>* glu_monoIdx = 0;
0396   vector<int>* glu_stereoIdx = 0;
0397   vector<float>* glu_x = 0;
0398   vector<float>* glu_y = 0;
0399   vector<float>* glu_z = 0;
0400   vector<float>* glu_xx = 0;
0401   vector<float>* glu_xy = 0;
0402   vector<float>* glu_yy = 0;
0403   vector<float>* glu_yz = 0;
0404   vector<float>* glu_zz = 0;
0405   vector<float>* glu_zx = 0;
0406   if (!hasPh2hits) {
0407     t->SetBranchAddress("glu_isBarrel", &glu_isBarrel);
0408     if (has910_det_lay) {
0409       t->SetBranchAddress("glu_subdet", &glu_det);
0410       t->SetBranchAddress("glu_layer", &glu_lay);
0411     } else {
0412       t->SetBranchAddress("glu_det", &glu_det);
0413       t->SetBranchAddress("glu_lay", &glu_lay);
0414     }
0415     t->SetBranchAddress("glu_detId", &glu_detId);
0416     t->SetBranchAddress("glu_monoIdx", &glu_monoIdx);
0417     t->SetBranchAddress("glu_stereoIdx", &glu_stereoIdx);
0418     t->SetBranchAddress("glu_x", &glu_x);
0419     t->SetBranchAddress("glu_y", &glu_y);
0420     t->SetBranchAddress("glu_z", &glu_z);
0421     t->SetBranchAddress("glu_xx", &glu_xx);
0422     t->SetBranchAddress("glu_xy", &glu_xy);
0423     t->SetBranchAddress("glu_yy", &glu_yy);
0424     t->SetBranchAddress("glu_yz", &glu_yz);
0425     t->SetBranchAddress("glu_zz", &glu_zz);
0426     t->SetBranchAddress("glu_zx", &glu_zx);
0427   }
0428 
0429   vector<short>* str_isBarrel = 0;
0430   vector<short>* str_isStereo = 0;
0431   vector<unsigned int>* str_det = 0;
0432   vector<unsigned int>* str_lay = 0;
0433   vector<unsigned int>* str_detId = 0;
0434   vector<unsigned int>* str_simType = 0;
0435   vector<float>* str_x = 0;
0436   vector<float>* str_y = 0;
0437   vector<float>* str_z = 0;
0438   vector<float>* str_xx = 0;
0439   vector<float>* str_xy = 0;
0440   vector<float>* str_yy = 0;
0441   vector<float>* str_yz = 0;
0442   vector<float>* str_zz = 0;
0443   vector<float>* str_zx = 0;
0444   vector<float>* str_chargePerCM = 0;
0445   vector<int>* str_csize = 0;
0446   vector<uint64_t>* str_usedMask = 0;
0447   vector<vector<int>>* str_simHitIdx = 0;
0448   vector<vector<float>>* str_chargeFraction = 0;
0449   if (!hasPh2hits) {
0450     t->SetBranchAddress("str_isBarrel", &str_isBarrel);
0451     t->SetBranchAddress("str_isStereo", &str_isStereo);
0452     if (has910_det_lay) {
0453       t->SetBranchAddress("str_subdet", &str_det);
0454       t->SetBranchAddress("str_layer", &str_lay);
0455     } else {
0456       t->SetBranchAddress("str_det", &str_det);
0457       t->SetBranchAddress("str_lay", &str_lay);
0458     }
0459     t->SetBranchAddress("str_detId", &str_detId);
0460     t->SetBranchAddress("str_simType", &str_simType);
0461     t->SetBranchAddress("str_x", &str_x);
0462     t->SetBranchAddress("str_y", &str_y);
0463     t->SetBranchAddress("str_z", &str_z);
0464     t->SetBranchAddress("str_xx", &str_xx);
0465     t->SetBranchAddress("str_xy", &str_xy);
0466     t->SetBranchAddress("str_yy", &str_yy);
0467     t->SetBranchAddress("str_yz", &str_yz);
0468     t->SetBranchAddress("str_zz", &str_zz);
0469     t->SetBranchAddress("str_zx", &str_zx);
0470     t->SetBranchAddress("str_chargePerCM", &str_chargePerCM);
0471     t->SetBranchAddress("str_clustSize", &str_csize);
0472     if (writeHitIterMasks) {
0473       t->SetBranchAddress("str_usedMask", &str_usedMask);
0474     }
0475 
0476     t->SetBranchAddress("str_simHitIdx", &str_simHitIdx);
0477     t->SetBranchAddress("str_chargeFraction", &str_chargeFraction);
0478   }
0479 
0480   vector<unsigned short>* ph2_isLower = 0;
0481   vector<unsigned short>* ph2_subdet = 0;
0482   vector<unsigned short>* ph2_layer = 0;
0483   vector<unsigned int>* ph2_detId = 0;
0484   vector<unsigned short>* ph2_simType = 0;
0485   vector<float>* ph2_x = 0;
0486   vector<float>* ph2_y = 0;
0487   vector<float>* ph2_z = 0;
0488   vector<float>* ph2_xx = 0;
0489   vector<float>* ph2_xy = 0;
0490   vector<float>* ph2_yy = 0;
0491   vector<float>* ph2_yz = 0;
0492   vector<float>* ph2_zz = 0;
0493   vector<float>* ph2_zx = 0;
0494   vector<uint64_t>* ph2_usedMask = 0;
0495   vector<vector<int>>* ph2_simHitIdx = 0;
0496   if (hasPh2hits && applyCCC)
0497     std::cout << "WARNING: applyCCC is set for Phase2 inputs: applyCCC will be ignored" << std::endl;
0498   if (hasPh2hits) {
0499     t->SetBranchAddress("ph2_isLower", &ph2_isLower);
0500     t->SetBranchAddress("ph2_subdet", &ph2_subdet);
0501     t->SetBranchAddress("ph2_layer", &ph2_layer);
0502     t->SetBranchAddress("ph2_detId", &ph2_detId);
0503     t->SetBranchAddress("ph2_simType", &ph2_simType);
0504     t->SetBranchAddress("ph2_x", &ph2_x);
0505     t->SetBranchAddress("ph2_y", &ph2_y);
0506     t->SetBranchAddress("ph2_z", &ph2_z);
0507     t->SetBranchAddress("ph2_xx", &ph2_xx);
0508     t->SetBranchAddress("ph2_xy", &ph2_xy);
0509     t->SetBranchAddress("ph2_yy", &ph2_yy);
0510     t->SetBranchAddress("ph2_yz", &ph2_yz);
0511     t->SetBranchAddress("ph2_zz", &ph2_zz);
0512     t->SetBranchAddress("ph2_zx", &ph2_zx);
0513     if (writeHitIterMasks) {
0514       t->SetBranchAddress("ph2_usedMask", &ph2_usedMask);
0515     }
0516     t->SetBranchAddress("ph2_simHitIdx", &ph2_simHitIdx);
0517   }
0518   vector<float> ph2_chargeFraction_dummy(16, 0.f);
0519 
0520   // beam spot
0521   float bsp_x;
0522   float bsp_y;
0523   float bsp_z;
0524   float bsp_sigmax;
0525   float bsp_sigmay;
0526   float bsp_sigmaz;
0527   t->SetBranchAddress("bsp_x", &bsp_x);
0528   t->SetBranchAddress("bsp_y", &bsp_y);
0529   t->SetBranchAddress("bsp_z", &bsp_z);
0530   t->SetBranchAddress("bsp_sigmax", &bsp_sigmax);
0531   t->SetBranchAddress("bsp_sigmay", &bsp_sigmay);
0532   t->SetBranchAddress("bsp_sigmaz", &bsp_sigmaz);
0533 
0534   long long totentries = t->GetEntries();
0535   long long savedEvents = 0;
0536 
0537   DataFile data_file;
0538   int outOptions = DataFile::ES_Seeds;
0539   if (writeRecTracks)
0540     outOptions |= DataFile::ES_CmsswTracks;
0541   if (writeHitIterMasks)
0542     outOptions |= DataFile::ES_HitIterMasks;
0543   outOptions |= DataFile::ES_BeamSpot;
0544 
0545   if (maxevt < 0)
0546     maxevt = totentries;
0547   data_file.openWrite(outputFileName, static_cast<int>(nTotalLayers), std::min(maxevt, totentries), outOptions);
0548 
0549   Event EE(0, static_cast<int>(nTotalLayers));
0550 
0551   int numFailCCC = 0;
0552   int numTotalStr = 0;
0553   // gDebug = 8;
0554 
0555   for (long long i = 0; savedEvents < maxevt && i < totentries && i < maxevt; ++i) {
0556     EE.reset(i);
0557 
0558     cout << "process entry i=" << i << " out of " << totentries << ", saved so far " << savedEvents
0559          << ", with max=" << maxevt << endl;
0560 
0561     t->GetEntry(i);
0562 
0563     cout << "edm event=" << event << endl;
0564 
0565     auto& bs = EE.beamSpot_;
0566     bs.x = bsp_x;
0567     bs.y = bsp_y;
0568     bs.z = bsp_z;
0569     bs.sigmaZ = bsp_sigmaz;
0570     bs.beamWidthX = bsp_sigmax;
0571     bs.beamWidthY = bsp_sigmay;
0572     //dxdz and dydz are not in the trackingNtuple at the moment
0573 
0574     if (!hasPh2hits) {
0575       for (unsigned int istr = 0; istr < str_lay->size(); ++istr) {
0576         if (str_chargePerCM->at(istr) < cutValueCCC)
0577           numFailCCC++;
0578         numTotalStr++;
0579       }
0580     }
0581 
0582     auto nSims = sim_q->size();
0583     if (nSims == 0) {
0584       cout << "branches not loaded" << endl;
0585       exit(1);
0586     }
0587     if (verbosity > 0)
0588       std::cout << __FILE__ << " " << __LINE__ << " nSims " << nSims << " nSeeds " << see_q->size() << " nRecT "
0589                 << trk_q->size() << std::endl;
0590 
0591     //find best matching tkIdx from a list of simhits indices
0592     auto bestTkIdx = [&](std::vector<int> const& shs, std::vector<float> const& shfs, int rhIdx, HitType rhType) {
0593       //assume that all simhits are associated
0594       int ibest = -1;
0595       int shbest = -1;
0596       float hpbest = -1;
0597       float tpbest = -1;
0598       float hfbest = -1;
0599 
0600       float maxfrac = -1;
0601       int ish = -1;
0602       int nshs = shs.size();
0603       for (auto const sh : shs) {
0604         ish++;
0605         auto tkidx = simhit_simTrkIdx->at(sh);
0606         //use only sh with available TP
0607         if (tkidx < 0)
0608           continue;
0609 
0610         auto hpx = simhit_px->at(sh);
0611         auto hpy = simhit_py->at(sh);
0612         auto hpz = simhit_pz->at(sh);
0613         auto hp = sqrt(hpx * hpx + hpy * hpy + hpz * hpz);
0614 
0615         //look only at hits with p> 50 MeV
0616         if (hp < 0.05f)
0617           continue;
0618 
0619         auto tpx = sim_px->at(tkidx);
0620         auto tpy = sim_py->at(tkidx);
0621         auto tpz = sim_pz->at(tkidx);
0622         auto tp = sqrt(tpx * tpx + tpy * tpy + tpz * tpz);
0623 
0624         //take only hits with hp> 0.5*tp
0625         if (hp < 0.5 * tp)
0626           continue;
0627 
0628         //pick tkidx corresponding to max hp/tp; .. this is probably redundant
0629         if (maxfrac < hp / tp) {
0630           maxfrac = hp / tp;
0631           ibest = tkidx;
0632           shbest = sh;
0633           hpbest = hp;
0634           tpbest = tp;
0635           hfbest = shfs[ish];
0636         }
0637       }
0638 
0639       //arbitration: a rechit with one matching sim is matched to sim if it's the first
0640       //FIXME: SOME BETTER SELECTION CAN BE DONE (it will require some more correlated knowledge)
0641       if (nshs == 1 && ibest >= 0) {
0642         auto const& srhIdxV = simhit_hitIdx->at(shbest);
0643         auto const& srhTypeV = simhit_hitType->at(shbest);
0644         int ih = -1;
0645         for (auto itype : srhTypeV) {
0646           ih++;
0647           if (HitType(itype) == rhType && srhIdxV[ih] != rhIdx) {
0648             ibest = -1;
0649             break;
0650           }
0651         }
0652       }
0653 
0654       if (ibest >= 0 && false) {
0655         std::cout << " best tkIdx " << ibest << " rh " << rhIdx << " for sh " << shbest << " out of " << shs.size()
0656                   << " hp " << hpbest << " chF " << hfbest << " tp " << tpbest << " process "
0657                   << simhit_process->at(shbest) << " particle " << simhit_particle->at(shbest) << std::endl;
0658         if (rhType == HitType::Strip) {
0659           std::cout << "    sh " << simhit_x->at(shbest) << ", " << simhit_y->at(shbest) << ", " << simhit_z->at(shbest)
0660                     << "  rh " << str_x->at(rhIdx) << ", " << str_y->at(rhIdx) << ", " << str_z->at(rhIdx) << std::endl;
0661         }
0662       }
0663       return ibest;
0664     };
0665 
0666     vector<Track>& simTracks_ = EE.simTracks_;
0667     vector<int> simTrackIdx_(sim_q->size(), -1);  //keep track of original index in ntuple
0668     vector<int> seedSimIdx(see_q->size(), -1);
0669     for (unsigned int isim = 0; isim < sim_q->size(); ++isim) {
0670       //load sim production vertex data
0671       auto iVtx = sim_parentVtxIdx->at(isim);
0672       constexpr float largeValF = 9999.f;
0673       float sim_prodx = iVtx >= 0 ? simvtx_x->at(iVtx) : largeValF;
0674       float sim_prody = iVtx >= 0 ? simvtx_y->at(iVtx) : largeValF;
0675       float sim_prodz = iVtx >= 0 ? simvtx_z->at(iVtx) : largeValF;
0676       //if (fabs(sim_eta->at(isim))>0.8) continue;
0677 
0678       vector<int> const& trkIdxV = sim_trkIdx->at(isim);
0679 
0680       //if (trkIdx<0) continue;
0681       //FIXME: CHECK IF THE LOOP AND BEST SELECTION IS NEEDED.
0682       //Pick the first
0683       const int trkIdx = trkIdxV.empty() ? -1 : trkIdxV[0];
0684 
0685       int nlay = 0;
0686       if (trkIdx >= 0) {
0687         std::vector<int> hitlay(nTotalLayers, 0);
0688         auto const& hits = trk_hitIdx->at(trkIdx);
0689         auto const& hitTypes = trk_hitType->at(trkIdx);
0690         auto nHits = hits.size();
0691         for (auto ihit = 0U; ihit < nHits; ++ihit) {
0692           auto ihIdx = hits[ihit];
0693           auto const ihType = HitType(hitTypes[ihit]);
0694 
0695           switch (ihType) {
0696             case HitType::Pixel: {
0697               int ipix = ihIdx;
0698               if (ipix < 0)
0699                 continue;
0700               int cmsswlay =
0701                   lnc.convertLayerNumber(pix_det->at(ipix), pix_lay->at(ipix), useMatched, -1, pix_z->at(ipix) > 0);
0702               if (cmsswlay >= 0 && cmsswlay < static_cast<int>(nTotalLayers))
0703                 hitlay[cmsswlay]++;
0704               break;
0705             }
0706             case HitType::Strip: {
0707               int istr = ihIdx;
0708               if (istr < 0)
0709                 continue;
0710               int cmsswlay = lnc.convertLayerNumber(
0711                   str_det->at(istr), str_lay->at(istr), useMatched, str_isStereo->at(istr), str_z->at(istr) > 0);
0712               if (cmsswlay >= 0 && cmsswlay < static_cast<int>(nTotalLayers))
0713                 hitlay[cmsswlay]++;
0714               break;
0715             }
0716             case HitType::Glued: {
0717               if (useMatched) {
0718                 int iglu = ihIdx;
0719                 if (iglu < 0)
0720                   continue;
0721                 int cmsswlay =
0722                     lnc.convertLayerNumber(glu_det->at(iglu), glu_lay->at(iglu), useMatched, -1, glu_z->at(iglu) > 0);
0723                 if (cmsswlay >= 0 && cmsswlay < static_cast<int>(nTotalLayers))
0724                   hitlay[cmsswlay]++;
0725               }
0726               break;
0727             }
0728             case HitType::Phase2OT: {
0729               int istr = ihIdx;
0730               if (istr < 0)
0731                 continue;
0732               int cmsswlay = lnc.convertLayerNumber(
0733                   ph2_subdet->at(istr), ph2_layer->at(istr), useMatched, ph2_isLower->at(istr), ph2_z->at(istr) > 0);
0734               if (cmsswlay >= 0 && cmsswlay < static_cast<int>(nTotalLayers))
0735                 hitlay[cmsswlay]++;
0736               break;
0737             }
0738             case HitType::Invalid:
0739               break;  //FIXME. Skip, really?
0740             default:
0741               throw std::logic_error("Track type can not be handled");
0742           }  //hit type
0743         }    //hits on track
0744         for (unsigned int i = 0; i < nTotalLayers; i++)
0745           if (hitlay[i] > 0)
0746             nlay++;
0747       }  //count nlay layers on matching reco track
0748 
0749       //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;
0750 
0751       SVector3 pos(sim_prodx, sim_prody, sim_prodz);
0752       SVector3 mom(sim_px->at(isim), sim_py->at(isim), sim_pz->at(isim));
0753       SMatrixSym66 err;
0754       err.At(0, 0) = sim_prodx * sim_prodx;
0755       err.At(1, 1) = sim_prody * sim_prody;
0756       err.At(2, 2) = sim_prodz * sim_prodz;
0757       err.At(3, 3) = sim_px->at(isim) * sim_px->at(isim);
0758       err.At(4, 4) = sim_py->at(isim) * sim_py->at(isim);
0759       err.At(5, 5) = sim_pz->at(isim) * sim_pz->at(isim);
0760       TrackState state(sim_q->at(isim), pos, mom, err);
0761       state.convertFromCartesianToCCS();
0762       //create track: store number of reco hits in place of track chi2; fill hits later
0763       //              set label to be its own index in the output file
0764       Track track(state, float(nlay), simTracks_.size(), 0, nullptr);
0765       if (sim_bunchCrossing->at(isim) == 0) {  //in time
0766         if (sim_event->at(isim) == 0)
0767           track.setProdType(Track::ProdType::Signal);
0768         else
0769           track.setProdType(Track::ProdType::InTimePU);
0770       } else {
0771         track.setProdType(Track::ProdType::OutOfTimePU);
0772       }
0773       if (trkIdx >= 0) {
0774         int seedIdx = trk_seedIdx->at(trkIdx);
0775         // Unused: auto const& shTypes = see_hitType->at(seedIdx);
0776         seedSimIdx[seedIdx] = simTracks_.size();
0777       }
0778       if (cleanSimTracks) {
0779         if (sim_nValid->at(isim) < cleanSimTrack_minSimHits)
0780           continue;
0781         if (cleanSimTrack_minRecHits > 0) {
0782           int nRecToSimHit = 0;
0783           for (unsigned int ipix = 0; ipix < pix_lay->size() && nRecToSimHit < cleanSimTrack_minRecHits; ++ipix) {
0784             int ilay = -1;
0785             ilay = lnc.convertLayerNumber(pix_det->at(ipix), pix_lay->at(ipix), useMatched, -1, pix_z->at(ipix) > 0);
0786             if (ilay < 0)
0787               continue;
0788             int simTkIdxNt = bestTkIdx(pix_simHitIdx->at(ipix), pix_chargeFraction->at(ipix), ipix, HitType::Pixel);
0789             if (simTkIdxNt >= 0)
0790               nRecToSimHit++;
0791           }
0792           if (hasPh2hits) {
0793             for (unsigned int istr = 0; istr < ph2_layer->size() && nRecToSimHit < cleanSimTrack_minRecHits; ++istr) {
0794               int ilay = -1;
0795               ilay = lnc.convertLayerNumber(
0796                   ph2_subdet->at(istr), ph2_layer->at(istr), useMatched, ph2_isLower->at(istr), ph2_z->at(istr) > 0);
0797               if (useMatched && !ph2_isLower->at(istr))
0798                 continue;
0799               if (ilay == -1)
0800                 continue;
0801               int simTkIdxNt = bestTkIdx(ph2_simHitIdx->at(istr), ph2_chargeFraction_dummy, istr, HitType::Phase2OT);
0802               if (simTkIdxNt >= 0)
0803                 nRecToSimHit++;
0804             }
0805           } else {
0806             if (useMatched) {
0807               for (unsigned int iglu = 0; iglu < glu_lay->size() && nRecToSimHit < cleanSimTrack_minRecHits; ++iglu) {
0808                 if (glu_isBarrel->at(iglu) == 0)
0809                   continue;
0810                 int igluMono = glu_monoIdx->at(iglu);
0811                 int simTkIdxNt =
0812                     bestTkIdx(str_simHitIdx->at(igluMono), str_chargeFraction->at(igluMono), igluMono, HitType::Strip);
0813                 if (simTkIdxNt >= 0)
0814                   nRecToSimHit++;
0815               }
0816             }
0817             for (unsigned int istr = 0; istr < str_lay->size() && nRecToSimHit < cleanSimTrack_minRecHits; ++istr) {
0818               int ilay = -1;
0819               ilay = lnc.convertLayerNumber(
0820                   str_det->at(istr), str_lay->at(istr), useMatched, str_isStereo->at(istr), str_z->at(istr) > 0);
0821               if (useMatched && str_isBarrel->at(istr) == 1 && str_isStereo->at(istr))
0822                 continue;
0823               if (ilay == -1)
0824                 continue;
0825               int simTkIdxNt = bestTkIdx(str_simHitIdx->at(istr), str_chargeFraction->at(istr), istr, HitType::Strip);
0826               if (simTkIdxNt >= 0)
0827                 nRecToSimHit++;
0828             }
0829           }
0830           if (nRecToSimHit < cleanSimTrack_minRecHits)
0831             continue;
0832         }  //count rec-to-sim hits
0833       }    //cleanSimTracks
0834 
0835       simTrackIdx_[isim] = simTracks_.size();
0836       simTracks_.push_back(track);
0837     }
0838 
0839     if (simTracks_.empty() and not writeAllEvents)
0840       continue;
0841 
0842     vector<Track>& seedTracks_ = EE.seedTracks_;
0843     vector<vector<int>> pixHitSeedIdx(pix_lay->size());
0844     vector<vector<int>> strHitSeedIdx(hasPh2hits ? 0 : str_lay->size());
0845     vector<vector<int>> gluHitSeedIdx(hasPh2hits ? 0 : glu_lay->size());
0846     vector<vector<int>> ph2HitSeedIdx(hasPh2hits ? ph2_layer->size() : 0);
0847     for (unsigned int is = 0; is < see_q->size(); ++is) {
0848       auto isAlgo = TrackAlgorithm(see_algo->at(is));
0849       if (not allSeeds)
0850         if (isAlgo != TrackAlgorithm::initialStep && isAlgo != TrackAlgorithm::hltIter0)
0851           continue;  //select seed in acceptance
0852       //if (see_pt->at(is)<0.5 || fabs(see_eta->at(is))>0.8) continue;//select seed in acceptance
0853       SVector3 pos = SVector3(see_stateTrajGlbX->at(is), see_stateTrajGlbY->at(is), see_stateTrajGlbZ->at(is));
0854       SVector3 mom = SVector3(see_stateTrajGlbPx->at(is), see_stateTrajGlbPy->at(is), see_stateTrajGlbPz->at(is));
0855       SMatrixSym66 err;
0856       if (hasCartCov) {
0857         err.At(0, 0) = see_stateCcov00->at(is);
0858         err.At(0, 1) = see_stateCcov01->at(is);
0859         err.At(0, 2) = see_stateCcov02->at(is);
0860         err.At(0, 3) = see_stateCcov03->at(is);
0861         err.At(0, 4) = see_stateCcov04->at(is);
0862         err.At(0, 5) = see_stateCcov05->at(is);
0863         err.At(1, 1) = see_stateCcov11->at(is);
0864         err.At(1, 2) = see_stateCcov12->at(is);
0865         err.At(1, 3) = see_stateCcov13->at(is);
0866         err.At(1, 4) = see_stateCcov14->at(is);
0867         err.At(1, 5) = see_stateCcov15->at(is);
0868         err.At(2, 2) = see_stateCcov22->at(is);
0869         err.At(2, 3) = see_stateCcov23->at(is);
0870         err.At(2, 4) = see_stateCcov24->at(is);
0871         err.At(2, 5) = see_stateCcov25->at(is);
0872         err.At(3, 3) = see_stateCcov33->at(is);
0873         err.At(3, 4) = see_stateCcov34->at(is);
0874         err.At(3, 5) = see_stateCcov35->at(is);
0875         err.At(4, 4) = see_stateCcov44->at(is);
0876         err.At(4, 5) = see_stateCcov45->at(is);
0877         err.At(5, 5) = see_stateCcov55->at(is);
0878       } else {
0879         auto const& vCov = see_stateCurvCov->at(is);
0880         assert(vCov.size() == 15);
0881         auto vCovP = vCov.begin();
0882         for (int i = 0; i < 5; ++i)
0883           for (int j = 0; j <= i; ++j)
0884             err.At(i, j) = *(vCovP++);
0885       }
0886       TrackState state(see_q->at(is), pos, mom, err);
0887       if (hasCartCov)
0888         state.convertFromCartesianToCCS();
0889       else
0890         state.convertFromGlbCurvilinearToCCS();
0891       Track track(state, 0, seedSimIdx[is], 0, nullptr);
0892       track.setAlgorithm(isAlgo);
0893       auto const& shTypes = see_hitType->at(is);
0894       auto const& shIdxs = see_hitIdx->at(is);
0895       if (not allSeeds)
0896         if (!((isAlgo == TrackAlgorithm::initialStep || isAlgo == TrackAlgorithm::hltIter0) &&
0897               std::count(shTypes.begin(), shTypes.end(), int(HitType::Pixel)) >= 3))
0898           continue;  //check algo and nhits
0899       for (unsigned int ip = 0; ip < shTypes.size(); ip++) {
0900         unsigned int hidx = shIdxs[ip];
0901         switch (HitType(shTypes[ip])) {
0902           case HitType::Pixel: {
0903             pixHitSeedIdx[hidx].push_back(seedTracks_.size());
0904             break;
0905           }
0906           case HitType::Strip: {
0907             strHitSeedIdx[hidx].push_back(seedTracks_.size());
0908             break;
0909           }
0910           case HitType::Glued: {
0911             if (not useMatched) {
0912               //decompose
0913               int uidx = glu_monoIdx->at(hidx);
0914               strHitSeedIdx[uidx].push_back(seedTracks_.size());
0915               uidx = glu_stereoIdx->at(hidx);
0916               strHitSeedIdx[uidx].push_back(seedTracks_.size());
0917             } else {
0918               gluHitSeedIdx[hidx].push_back(seedTracks_.size());
0919             }
0920             break;
0921           }
0922           case HitType::Phase2OT: {
0923             ph2HitSeedIdx[hidx].push_back(seedTracks_.size());
0924             break;
0925           }
0926           case HitType::Invalid:
0927             break;  //FIXME. Skip, really?
0928           default:
0929             throw std::logic_error("Track hit type can not be handled");
0930         }  //switch( HitType
0931       }
0932       seedTracks_.push_back(track);
0933     }
0934 
0935     if (seedTracks_.empty() and not writeAllEvents)
0936       continue;
0937 
0938     vector<Track>& cmsswTracks_ = EE.cmsswTracks_;
0939     vector<vector<int>> pixHitRecIdx(pix_lay->size());
0940     vector<vector<int>> strHitRecIdx(hasPh2hits ? 0 : str_lay->size());
0941     vector<vector<int>> gluHitRecIdx(hasPh2hits ? 0 : glu_lay->size());
0942     vector<vector<int>> ph2HitRecIdx(hasPh2hits ? ph2_layer->size() : 0);
0943     for (unsigned int ir = 0; ir < trk_q->size(); ++ir) {
0944       //check the origin; redundant for initialStep ntuples
0945       if (not allSeeds)
0946         if ((trk_algoMask->at(ir) & ((1 << int(TrackAlgorithm::initialStep)) | (1 << int(TrackAlgorithm::hltIter0)))) ==
0947             0) {
0948           if (verbosity > 1) {
0949             std::cout << "track " << ir << " failed algo selection for " << int(TrackAlgorithm::initialStep)
0950                       << ": mask " << trk_algoMask->at(ir) << " origAlgo " << trk_originalAlgo->at(ir) << " algo "
0951                       << trk_algo->at(ir) << std::endl;
0952           }
0953           continue;
0954         }
0955       //fill the state in CCS upfront
0956       SMatrixSym66 err;
0957       /*    
0958     vx = -dxy*sin(phi) - pt*cos(phi)/p*pz/p*dz;
0959     vy =  dxy*cos(phi) - pt*sin(phi)/p*pz/p*dz;
0960     vz = dz*pt*pt/p/p;
0961     //partial: ignores cross-terms
0962     c(vx,vx) = c(dxy,dxy)*sin(phi)*sin(phi) + c(dz,dz)*pow(pt*cos(phi)/p*pz/p ,2);
0963     c(vx,vy) = -c(dxy,dxy)*cos(phi)*sin(phi) + c(dz,dz)*cos(phi)*sin(phi)*pow(pt/p*pz/p, 2);
0964     c(vy,vy) = c(dxy,dxy)*cos(phi)*cos(phi) + c(dz,dz)*pow(pt*sin(phi)/p*pz/p ,2);
0965     c(vx,vz) = -c(dz,dz)*pt*pt/p/p*pt/p*pz/p*cos(phi);
0966     c(vy,vz) = -c(dz,dz)*pt*pt/p/p*pt/p*pz/p*sin(phi);
0967     c(vz,vz) = c(dz,dz)*pow(pt*pt/p/p, 2);
0968       */
0969       float pt = trk_pt->at(ir);
0970       float pz = trk_pz->at(ir);
0971       float p2 = pt * pt + pz * pz;
0972       float phi = trk_phi->at(ir);
0973       float sP = sin(phi);
0974       float cP = cos(phi);
0975       float dxyErr2 = trk_dxyErr->at(ir);
0976       dxyErr2 *= dxyErr2;
0977       float dzErr2 = trk_dzErr->at(ir);
0978       dzErr2 *= dzErr2;
0979       float dzErrF2 = trk_dzErr->at(ir) * (pt * pz / p2);
0980       dzErr2 *= dzErr2;
0981       err.At(0, 0) = dxyErr2 * sP * sP + dzErrF2 * cP * cP;
0982       err.At(0, 1) = -dxyErr2 * cP * sP + dzErrF2 * cP * sP;
0983       err.At(1, 1) = dxyErr2 * cP * cP + dzErrF2 * sP * sP;
0984       err.At(0, 2) = -dzErrF2 * cP * pt / pz;
0985       err.At(1, 2) = -dzErrF2 * sP * pt / pz;
0986       err.At(2, 2) = dzErr2 * std::pow((pt * pt / p2), 2);
0987       err.At(3, 3) = std::pow(trk_ptErr->at(ir) / pt / pt, 2);
0988       err.At(4, 4) = std::pow(trk_phiErr->at(ir), 2);
0989       err.At(5, 5) = std::pow(trk_lambdaErr->at(ir), 2);
0990       SVector3 pos = SVector3(trk_refpoint_x->at(ir), trk_refpoint_y->at(ir), trk_refpoint_z->at(ir));
0991       SVector3 mom = SVector3(1.f / pt, phi, M_PI_2 - trk_lambda->at(ir));
0992       TrackState state(trk_q->at(ir), pos, mom, err);
0993       Track track(state, trk_nChi2->at(ir), trk_seedIdx->at(ir), 0, nullptr);  //hits are filled later
0994       track.setAlgorithm(TrackAlgorithm(trk_originalAlgo->at(ir)));
0995       auto const& hTypes = trk_hitType->at(ir);
0996       auto const& hIdxs = trk_hitIdx->at(ir);
0997       for (unsigned int ip = 0; ip < hTypes.size(); ip++) {
0998         unsigned int hidx = hIdxs[ip];
0999         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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!