Back to home page

Project CMSSW displayed by LXR
 
 

    

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 //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       "  --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;  //Nominal value (from first iteration of CMSSW) is 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   }  //while arguments
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   //sim tracks
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;  //simHit count, actually
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   //simvtx
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   //simhit
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   //rec tracks
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   //seeds
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;  //PCA parameters
0272   std::vector<float>* see_pt = 0;   //PCA parameters
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   //pixel hits
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   //these were renamed in CMSSW_9_1_0: auto-detect
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   //strip hits
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   // beam spot
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   // gDebug = 8;
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     //dxdz and dydz are not in the trackingNtuple at the moment
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     //find best matching tkIdx from a list of simhits indices
0542     auto bestTkIdx = [&](std::vector<int> const& shs, std::vector<float> const& shfs, int rhIdx, HitType rhType) {
0543       //assume that all simhits are associated
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         //use only sh with available TP
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         //look only at hits with p> 50 MeV
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         //take only hits with hp> 0.5*tp
0575         if (hp < 0.5 * tp)
0576           continue;
0577 
0578         //pick tkidx corresponding to max hp/tp; .. this is probably redundant
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       //arbitration: a rechit with one matching sim is matched to sim if it's the first
0590       //FIXME: SOME BETTER SELECTION CAN BE DONE (it will require some more correlated knowledge)
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);  //keep track of original index in ntuple
0618     vector<int> seedSimIdx(see_q->size(), -1);
0619     for (unsigned int isim = 0; isim < sim_q->size(); ++isim) {
0620       //load sim production vertex data
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       //if (fabs(sim_eta->at(isim))>0.8) continue;
0627 
0628       vector<int> const& trkIdxV = sim_trkIdx->at(isim);
0629 
0630       //if (trkIdx<0) continue;
0631       //FIXME: CHECK IF THE LOOP AND BEST SELECTION IS NEEDED.
0632       //Pick the first
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;  //FIXME. Skip, really?
0680             default:
0681               throw std::logic_error("Track type can not be handled");
0682           }  //hit type
0683         }    //hits on track
0684         for (unsigned int i = 0; i < nTotalLayers; i++)
0685           if (hitlay[i] > 0)
0686             nlay++;
0687       }  //count nlay layers on matching reco track
0688 
0689       //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;
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       //create track: store number of reco hits in place of track chi2; fill hits later
0703       //              set label to be its own index in the output file
0704       Track track(state, float(nlay), simTracks_.size(), 0, nullptr);
0705       if (sim_bunchCrossing->at(isim) == 0) {  //in time
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         // Unused: auto const& shTypes = see_hitType->at(seedIdx);
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         }  //count rec-to-sim hits
0758       }    //cleanSimTracks
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;  //select seed in acceptance
0776       //if (see_pt->at(is)<0.5 || fabs(see_eta->at(is))>0.8) continue;//select seed in acceptance
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;  //check algo and nhits
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               //decompose
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;  //FIXME. Skip, really?
0848           default:
0849             throw std::logic_error("Track hit type can not be handled");
0850         }  //switch( HitType
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       //check the origin; redundant for initialStep ntuples
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       //fill the state in CCS upfront
0875       SMatrixSym66 err;
0876       /*    
0877     vx = -dxy*sin(phi) - pt*cos(phi)/p*pz/p*dz;
0878     vy =  dxy*cos(phi) - pt*sin(phi)/p*pz/p*dz;
0879     vz = dz*pt*pt/p/p;
0880     //partial: ignores cross-terms
0881     c(vx,vx) = c(dxy,dxy)*sin(phi)*sin(phi) + c(dz,dz)*pow(pt*cos(phi)/p*pz/p ,2);
0882     c(vx,vy) = -c(dxy,dxy)*cos(phi)*sin(phi) + c(dz,dz)*cos(phi)*sin(phi)*pow(pt/p*pz/p, 2);
0883     c(vy,vy) = c(dxy,dxy)*cos(phi)*cos(phi) + c(dz,dz)*pow(pt*sin(phi)/p*pz/p ,2);
0884     c(vx,vz) = -c(dz,dz)*pt*pt/p/p*pt/p*pz/p*cos(phi);
0885     c(vy,vz) = -c(dz,dz)*pt*pt/p/p*pt/p*pz/p*sin(phi);
0886     c(vz,vz) = c(dz,dz)*pow(pt*pt/p/p, 2);
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);  //hits are filled later
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             //cout << "pix=" << hidx << " track=" << cmsswTracks_.size() << endl;
0921             pixHitRecIdx[hidx].push_back(cmsswTracks_.size());
0922             break;
0923           }
0924           case HitType::Strip: {
0925             //cout << "pix=" << hidx << " track=" << cmsswTracks_.size() << endl;
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             //cout << "pix=" << hidx << " track=" << cmsswTracks_.size() << endl;
0933             gluHitRecIdx[hidx].push_back(cmsswTracks_.size());
0934             break;
0935           }
0936           case HitType::Invalid:
0937             break;  //FIXME. Skip, really?
0938           default:
0939             throw std::logic_error("Track hit type can not be handled");
0940         }  //switch( HitType
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;  //switch to index in simTracks_
0961 
0962       //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;
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         //cout << "xxx ipix=" << ipix << " seed=" << pixHitSeedIdx[ipix][is] << endl;
0976         seedTracks_[pixHitSeedIdx[ipix][is]].addHitIdx(layerHits_[ilay].size(), ilay, 0);  //per-hit chi2 is not known
0977       }
0978       for (unsigned int ir = 0; ir < pixHitRecIdx[ipix].size(); ir++) {
0979         //cout << "xxx ipix=" << ipix << " recTrack=" << pixHitRecIdx[ipix][ir] << endl;
0980         cmsswTracks_[pixHitRecIdx[ipix][ir]].addHitIdx(layerHits_[ilay].size(), ilay, 0);  //per-hit chi2 is not known
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;  //switch to index in simTracks_
1000 
1001         int ilay = lnc.convertLayerNumber(glu_det->at(iglu), glu_lay->at(iglu), useMatched, -1, glu_z->at(iglu) > 0);
1002         // 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;
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           //cout << "xxx iglu=" << iglu << " seed=" << gluHitSeedIdx[iglu][ir] << endl;
1016           seedTracks_[gluHitSeedIdx[iglu][ir]].addHitIdx(layerHits_[ilay].size(), ilay, 0);  //per-hit chi2 is not known
1017         }
1018         for (unsigned int ir = 0; ir < gluHitRecIdx[iglu].size(); ir++) {
1019           //cout << "xxx iglu=" << iglu << " recTrack=" << gluHitRecIdx[iglu][ir] << endl;
1020           cmsswTracks_[gluHitRecIdx[iglu][ir]].addHitIdx(layerHits_[ilay].size(), ilay, 0);  //per-hit chi2 is not known
1021         }
1022 
1023         // QQQQ module-id-in-layer, adc and phi/theta spans are not done for matched hits.
1024         // Will we ever use / need this?
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;  //switch to index in simTracks_
1050 
1051       bool passCCC = applyCCC ? (str_chargePerCM->at(istr) > cutValueCCC) : true;
1052 
1053       //if (str_onTrack->at(istr)==0) continue;//do not consider hits that are not on track!
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         //cout << "xxx istr=" << istr << " seed=" << strHitSeedIdx[istr][ir] << endl;
1070         if (passCCC)
1071           seedTracks_[strHitSeedIdx[istr][ir]].addHitIdx(layerHits_[ilay].size(), ilay, 0);  //per-hit chi2 is not known
1072         else
1073           seedTracks_[strHitSeedIdx[istr][ir]].addHitIdx(-9, ilay, 0);
1074       }
1075       for (unsigned int ir = 0; ir < strHitRecIdx[istr].size(); ir++) {
1076         //cout << "xxx istr=" << istr << " recTrack=" << strHitRecIdx[istr][ir] << endl;
1077         if (passCCC)
1078           cmsswTracks_[strHitRecIdx[istr][ir]].addHitIdx(layerHits_[ilay].size(), ilay, 0);  //per-hit chi2 is not known
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     // Seed % hit statistics
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         }  //if (writeRecTracks){
1226 
1227       }  //verbosity>1
1228     }    //verbosity>0
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));  //Includes those that failed the cluster charge cut
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 }
[ 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!