File indexing completed on 2024-04-06 12:27:13
0001
0002 #include "SegSelector.h"
0003
0004
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006
0007 #include "TFile.h"
0008 #include "TVector3.h"
0009
0010 #include <iostream>
0011 #include <fstream>
0012 #include <map>
0013 #include <utility>
0014 #include <string>
0015 #include <stdio.h>
0016 #include <algorithm>
0017
0018
0019 using namespace std;
0020 using namespace edm;
0021
0022
0023 SegSelector::SegSelector(const ParameterSet& pset) {
0024
0025
0026 debug = pset.getUntrackedParameter<bool>("debug");
0027 recHitLabel = pset.getUntrackedParameter<string>("recHitLabel");
0028 cscSegmentLabel = pset.getUntrackedParameter<string>("cscSegmentLabel");
0029 dtrecHitLabel = pset.getUntrackedParameter<string>("dtrecHitLabel");
0030 dtSegmentLabel = pset.getUntrackedParameter<string>("dtSegmentLabel");
0031 simHitLabel = pset.getUntrackedParameter<string>("simHitLabel");
0032 simTrackLabel = pset.getUntrackedParameter<string>("simTrackLabel");
0033 }
0034
0035
0036 SegSelector::~SegSelector() {
0037 if (debug)
0038 cout << "[SeedQualityAnalysis] Destructor called" << endl;
0039 }
0040
0041
0042
0043
0044 std::vector<SimSegment> SegSelector::Sim_DTSegments(int trkId,
0045 Handle<edm::PSimHitContainer> dsimHits,
0046 ESHandle<DTGeometry> dtGeom) {
0047 SimSegment sim_dseg;
0048
0049 hit_V1.clear();
0050 sDT_v.clear();
0051 int d1 = 0;
0052 for (PSimHitContainer::const_iterator sh_i = dsimHits->begin(); sh_i != dsimHits->end(); ++sh_i) {
0053 if (static_cast<int>((*sh_i).trackId()) != trkId)
0054 continue;
0055 if (abs((*sh_i).particleType()) != 13)
0056 continue;
0057
0058 DTLayerId detId = DTLayerId((*sh_i).detUnitId());
0059
0060 int d2 = (100000 * detId.station()) + (1000 * detId.sector()) + (100 + detId.wheel());
0061
0062 if (d1 == 0) {
0063 d1 = d2;
0064 hit_V1.clear();
0065 hit_V1.push_back(*sh_i);
0066 } else if (d1 == d2) {
0067 DTLayerId detId1 = DTLayerId(hit_V1[hit_V1.size() - 1].detUnitId());
0068 if (detId1.layer() != detId.layer()) {
0069 hit_V1.push_back(*sh_i);
0070 }
0071 if ((hit_V1.size() == 12) || ((hit_V1.size() == 8) && (d2 / 100000 == 4))) {
0072 DTSimHitFit(dtGeom);
0073 DTChamberId det_seg(detId.wheel(), detId.station(), detId.sector());
0074 sim_dseg.chamber_type = 2;
0075 sim_dseg.dt_DetId = det_seg;
0076 sim_dseg.sLocalOrg = LSimOrg1;
0077 sim_dseg.sGlobalVec = GSimVec1;
0078 sim_dseg.sGlobalOrg = GSimOrg1;
0079 sim_dseg.simhit_v = hit_V1;
0080 sDT_v.push_back(sim_dseg);
0081 }
0082 } else {
0083 if ((hit_V1.size() < 12) && (hit_V1.size() > 5)) {
0084 DTSimHitFit(dtGeom);
0085 int old_st = d1 / 100000;
0086 int old_se = (d1 - (old_st * 100000)) / 1000;
0087 int old_wl = (d1 - (old_st * 100000) - (old_se * 1000)) - 100;
0088 DTChamberId det_seg(old_wl, old_st, old_se);
0089 sim_dseg.chamber_type = 2;
0090 sim_dseg.dt_DetId = det_seg;
0091 sim_dseg.sLocalOrg = LSimOrg1;
0092 sim_dseg.sGlobalVec = GSimVec1;
0093 sim_dseg.sGlobalOrg = GSimOrg1;
0094 sim_dseg.simhit_v = hit_V1;
0095 sDT_v.push_back(sim_dseg);
0096 }
0097 d1 = d2;
0098 hit_V1.clear();
0099 hit_V1.push_back(*sh_i);
0100 }
0101 }
0102 return sDT_v;
0103 }
0104
0105
0106 std::vector<SimSegment> SegSelector::Sim_CSCSegments(int trkId,
0107 Handle<edm::PSimHitContainer> csimHits,
0108 ESHandle<CSCGeometry> cscGeom) {
0109 SimSegment sim_cseg;
0110
0111
0112 hit_V.clear();
0113 sCSC_v.clear();
0114 int d1 = 0;
0115 for (PSimHitContainer::const_iterator sh_i = csimHits->begin(); sh_i != csimHits->end(); ++sh_i) {
0116 if (static_cast<int>((*sh_i).trackId()) != trkId)
0117 continue;
0118 if (abs((*sh_i).particleType()) != 13)
0119 continue;
0120
0121 CSCDetId detId = CSCDetId((*sh_i).detUnitId());
0122
0123 int d2 = (1000 * detId.station()) + (100 * detId.ring()) + detId.chamber();
0124 if (d1 == 0) {
0125 d1 = d2;
0126 hit_V.clear();
0127 hit_V.push_back(*sh_i);
0128 } else if (d1 == d2) {
0129 hit_V.push_back(*sh_i);
0130 if (hit_V.size() == 6) {
0131 CSCSimHitFit(cscGeom);
0132 int old_st = d1 / 1000;
0133 int old_rg = (d1 - (old_st * 1000)) / 100;
0134 int old_cb = d1 - (old_st * 1000) - (old_rg * 100);
0135 CSCDetId det_seg(detId.endcap(), old_st, old_rg, old_cb, 0);
0136 sim_cseg.chamber_type = 1;
0137 sim_cseg.csc_DetId = det_seg;
0138 sim_cseg.sLocalOrg = LSimOrg;
0139 sim_cseg.sGlobalVec = GSimVec;
0140 sim_cseg.sGlobalOrg = GSimOrg;
0141 sim_cseg.simhit_v = hit_V;
0142 sCSC_v.push_back(sim_cseg);
0143 }
0144 } else {
0145 if ((hit_V.size() < 6) && (hit_V.size() > 2)) {
0146 CSCSimHitFit(cscGeom);
0147 CSCDetId det_seg(detId.endcap(), detId.station(), detId.ring(), detId.chamber(), 0);
0148 sim_cseg.chamber_type = 1;
0149 sim_cseg.csc_DetId = det_seg;
0150 sim_cseg.sLocalOrg = LSimOrg;
0151 sim_cseg.sGlobalVec = GSimVec;
0152 sim_cseg.sGlobalOrg = GSimOrg;
0153 sim_cseg.simhit_v = hit_V;
0154 sCSC_v.push_back(sim_cseg);
0155 }
0156 d1 = d2;
0157 hit_V.clear();
0158 hit_V.push_back(*sh_i);
0159 }
0160 }
0161 return sCSC_v;
0162 }
0163
0164
0165 std::vector<DTRecSegment4D> SegSelector::Select_DTSeg(Handle<DTRecSegment4DCollection> dtSeg,
0166 ESHandle<DTGeometry> dtGeom,
0167 std::vector<SimSegment> sDT_v1) {
0168 dtseg_V.clear();
0169 for (std::vector<SimSegment>::const_iterator it1 = sDT_v1.begin(); it1 != sDT_v1.end(); it1++) {
0170
0171 std::vector<DTRecSegment4D> segtemp;
0172 segtemp.clear();
0173 for (DTRecSegment4DCollection::const_iterator it2 = dtSeg->begin(); it2 != dtSeg->end(); it2++) {
0174
0175 if (!(*it2).hasPhi())
0176 continue;
0177 if ((*it2).chamberId().station() < 4 && !(*it2).hasZed())
0178 continue;
0179
0180 if ((*it1).dt_DetId != (*it2).chamberId())
0181 continue;
0182
0183 segtemp.push_back(*it2);
0184 }
0185 if (segtemp.size() == 1) {
0186 dtseg_V.push_back(segtemp[0]);
0187 }
0188 if (segtemp.size() > 1) {
0189
0190 LongDTSegment(segtemp);
0191 if (longsegV1.size() == 1) {
0192 dtseg_V.push_back(longsegV1[0]);
0193 }
0194 if (longsegV1.size() > 1) {
0195
0196 std::vector<double> dgv1(7, 999.);
0197 double k = 0.;
0198 DTRecSegment4D closestseg;
0199 for (std::vector<DTRecSegment4D>::const_iterator it3 = longsegV1.begin(); it3 != longsegV1.end(); it3++) {
0200 k += 1.0;
0201 DTChamberId rDetId = (*it3).chamberId();
0202 LocalVector rec_v = (*it3).localDirection();
0203 LocalPoint rec_o = (*it3).localPosition();
0204 const DTChamber* dtchamber = dtGeom->chamber(rDetId);
0205 GlobalVector g_rec_v = dtchamber->toGlobal(rec_v);
0206 GlobalPoint g_rec_o = dtchamber->toGlobal(rec_o);
0207
0208 std::vector<double> dgv2(7, 999.);
0209 dgv2[0] = ((*it1).sGlobalVec).x() - g_rec_v.x();
0210 dgv2[1] = ((*it1).sGlobalVec).y() - g_rec_v.y();
0211 dgv2[2] = ((*it1).sGlobalVec).z() - g_rec_v.z();
0212 dgv2[3] = fabs(((*it1).sGlobalOrg).x() - g_rec_o.x());
0213 dgv2[4] = fabs(((*it1).sGlobalOrg).y() - g_rec_o.y());
0214 dgv2[5] = sqrt((dgv2[0] * dgv2[0]) + (dgv2[1] * dgv2[1]) + (dgv2[2] * dgv2[2]));
0215 dgv2[6] = k;
0216
0217 if (k == 1.0) {
0218 dgv1 = dgv2;
0219 closestseg = *it3;
0220 } else {
0221 closestseg = (dgv2[3] < dgv1[3]) ? (*it3) : closestseg;
0222 dgv1 = (dgv2[3] < dgv1[3]) ? dgv2 : dgv1;
0223 }
0224 }
0225 dtseg_V.push_back(closestseg);
0226 }
0227 }
0228 }
0229 return dtseg_V;
0230 }
0231
0232
0233 std::vector<CSCSegment> SegSelector::Select_CSCSeg(Handle<CSCSegmentCollection> cscSeg,
0234 ESHandle<CSCGeometry> cscGeom,
0235 std::vector<SimSegment> sCSC_v1) {
0236 cscseg_V.clear();
0237 for (std::vector<SimSegment>::const_iterator it1 = sCSC_v1.begin(); it1 != sCSC_v1.end(); it1++) {
0238
0239 std::vector<CSCSegment> segtemp;
0240 segtemp.clear();
0241 for (CSCSegmentCollection::const_iterator seg_i = cscSeg->begin(); seg_i != cscSeg->end(); seg_i++) {
0242 CSCDetId rDetId = (CSCDetId)(*seg_i).cscDetId();
0243 if ((*it1).csc_DetId == rDetId) {
0244 segtemp.push_back(*seg_i);
0245 }
0246 }
0247 if (segtemp.size() == 1) {
0248 cscseg_V.push_back(segtemp[0]);
0249 }
0250 if (segtemp.size() > 1) {
0251
0252 LongCSCSegment(segtemp);
0253 if (longsegV.size() == 1) {
0254 cscseg_V.push_back(longsegV[0]);
0255 }
0256 if (longsegV.size() > 1) {
0257
0258 std::vector<double> dgv1(7, 999.);
0259 double k = 0.;
0260 CSCSegment closestseg;
0261 for (std::vector<CSCSegment>::const_iterator i1 = longsegV.begin(); i1 != longsegV.end(); i1++) {
0262 k += 1.0;
0263 CSCDetId rDetId = (CSCDetId)(*i1).cscDetId();
0264 LocalVector rec_v = (*i1).localDirection();
0265 LocalPoint rec_o = (*i1).localPosition();
0266 const CSCChamber* cscchamber = cscGeom->chamber(rDetId);
0267 GlobalVector g_rec_v = cscchamber->toGlobal(rec_v);
0268 GlobalPoint g_rec_o = cscchamber->toGlobal(rec_o);
0269
0270 std::vector<double> dgv2(7, 999.);
0271 dgv2[0] = ((*it1).sGlobalVec).x() - g_rec_v.x();
0272 dgv2[1] = ((*it1).sGlobalVec).y() - g_rec_v.y();
0273 dgv2[2] = ((*it1).sGlobalVec).z() - g_rec_v.z();
0274 dgv2[3] = fabs(((*it1).sGlobalOrg).x() - g_rec_o.x());
0275 dgv2[4] = fabs(((*it1).sGlobalOrg).y() - g_rec_o.y());
0276 dgv2[5] = sqrt((dgv2[0] * dgv2[0]) + (dgv2[1] * dgv2[1]) + (dgv2[2] * dgv2[2]));
0277 dgv2[6] = k;
0278
0279 if (k == 1.0) {
0280 dgv1 = dgv2;
0281 closestseg = *i1;
0282 } else {
0283 closestseg = (dgv2[3] < dgv1[3]) ? (*i1) : closestseg;
0284 dgv1 = (dgv2[3] < dgv1[3]) ? dgv2 : dgv1;
0285 }
0286 }
0287 cscseg_V.push_back(closestseg);
0288 }
0289 }
0290 }
0291 return cscseg_V;
0292 }
0293
0294
0295 void SegSelector::LongCSCSegment(std::vector<CSCSegment> cscsegs) {
0296 int n = 0;
0297 CSCSegment longseg;
0298 longsegV.clear();
0299 for (std::vector<CSCSegment>::const_iterator i1 = cscsegs.begin(); i1 != cscsegs.end(); i1++) {
0300 n++;
0301 if (n == 1) {
0302 longseg = *i1;
0303 longsegV.push_back(*i1);
0304 } else if ((*i1).nRecHits() == longseg.nRecHits()) {
0305 longsegV.push_back(*i1);
0306 } else {
0307 longseg = ((*i1).nRecHits() > longseg.nRecHits()) ? (*i1) : longseg;
0308 longsegV.clear();
0309 longsegV.push_back(longseg);
0310 }
0311 }
0312 }
0313
0314
0315 void SegSelector::LongDTSegment(std::vector<DTRecSegment4D> dtsegs) {
0316 int n = 0;
0317 int nh = 0;
0318 DTRecSegment4D longseg;
0319 longsegV1.clear();
0320 for (std::vector<DTRecSegment4D>::const_iterator i1 = dtsegs.begin(); i1 != dtsegs.end(); i1++) {
0321 n++;
0322 int n_phi = 0;
0323 int n_zed = 0;
0324 if ((*i1).hasPhi())
0325 n_phi = ((*i1).phiSegment())->specificRecHits().size();
0326 if ((*i1).hasZed())
0327 n_zed = ((*i1).zSegment())->specificRecHits().size();
0328 int n_hits = n_phi + n_zed;
0329
0330 if (n == 1) {
0331 longseg = *i1;
0332 longsegV1.push_back(*i1);
0333 nh = n_hits;
0334 } else if (n_hits == nh) {
0335 longsegV1.push_back(*i1);
0336 } else {
0337 longseg = (n_hits > nh) ? (*i1) : longseg;
0338 longsegV1.clear();
0339 longsegV1.push_back(longseg);
0340 }
0341 }
0342 }
0343
0344
0345 void SegSelector::DTSimHitFit(ESHandle<DTGeometry> dtGeom) {
0346 std::vector<PSimHit> sp;
0347 std::vector<PSimHit> sp1;
0348 std::vector<PSimHit> sp2;
0349 sp1.clear();
0350 sp2.clear();
0351 DTLayerId DT_Id;
0352 for (std::vector<PSimHit>::const_iterator sh_i = hit_V1.begin(); sh_i != hit_V1.end(); ++sh_i) {
0353 DT_Id = DTLayerId((*sh_i).detUnitId());
0354 if ((DT_Id.superLayer() == 1) || (DT_Id.superLayer() == 3)) {
0355 sp1.push_back(*sh_i);
0356 }
0357 if (DT_Id.superLayer() == 2) {
0358 sp2.push_back(*sh_i);
0359 }
0360 }
0361 DTLayerId DT_Id0(DT_Id.wheel(), DT_Id.station(), DT_Id.sector(), 0, 0);
0362
0363
0364
0365
0366
0367 double par1[2] = {0.0};
0368 double par2[2] = {0.0};
0369 for (int i = 1; i < 4; i++) {
0370 double sum[4] = {0.0};
0371 double par[2] = {0.0};
0372 double N = 0.0;
0373 if (i == 1) {
0374 sp = sp1;
0375 }
0376 if (i == 2) {
0377 sp = sp2;
0378 }
0379 if ((i == 3) && (sp2.size() > 0))
0380 continue;
0381 if (i == 3) {
0382 sp = sp1;
0383 }
0384 for (std::vector<PSimHit>::const_iterator sh_i = sp.begin(); sh_i != sp.end(); ++sh_i) {
0385 DTWireId DT_Id = DTWireId((*sh_i).detUnitId());
0386 LocalPoint dt_xyz = (*sh_i).localPosition();
0387
0388 const DTLayer* dtlayer = dtGeom->layer(DT_Id);
0389 const DTChamber* dtchamber = dtGeom->chamber(DT_Id);
0390 GlobalPoint gdt_xyz = dtlayer->toGlobal(dt_xyz);
0391 LocalPoint ldt_xyz = dtchamber->toLocal(gdt_xyz);
0392
0393 N += 1.0;
0394 if (i == 1) {
0395 sum[0] += ldt_xyz.x();
0396 sum[1] += ldt_xyz.z();
0397 sum[2] += ldt_xyz.x() * ldt_xyz.z();
0398 sum[3] += ldt_xyz.z() * ldt_xyz.z();
0399 }
0400 if ((i == 2) || (i == 3)) {
0401 sum[0] += ldt_xyz.y();
0402 sum[1] += ldt_xyz.z();
0403 sum[2] += ldt_xyz.y() * ldt_xyz.z();
0404 sum[3] += ldt_xyz.z() * ldt_xyz.z();
0405 }
0406 }
0407 par[0] = ((sum[0] * sum[3]) - (sum[2] * sum[1])) / ((N * sum[3]) - (sum[1] * sum[1]));
0408 par[1] = ((N * sum[2]) - (sum[0] * sum[1])) / ((N * sum[3]) - (sum[1] * sum[1]));
0409 if (i == 1) {
0410 par1[0] = par[0];
0411 par1[1] = par[1];
0412 }
0413 if ((i == 2) || (i == 3)) {
0414 par2[0] = par[0];
0415 par2[1] = par[1];
0416 }
0417 }
0418
0419 double v_zz = -1. / sqrt((par1[1] * par1[1]) + (par2[1] * par2[1]) + 1.0);
0420 double v_xx = v_zz * par1[1];
0421 double v_yy = v_zz * par2[1];
0422 LSimVec1 = LocalVector(v_xx, v_yy, v_zz);
0423 LSimOrg1 = LocalPoint(par1[0], par2[0], 0.);
0424
0425 const DTLayer* DT_layer = dtGeom->layer(DT_Id0);
0426 GSimOrg1 = DT_layer->toGlobal(LSimOrg1);
0427 GSimVec1 = DT_layer->toGlobal(LSimVec1);
0428 }
0429
0430
0431 void SegSelector::CSCSimHitFit(ESHandle<CSCGeometry> cscGeom) {
0432 bool rv_flag = false;
0433 for (std::vector<PSimHit>::const_iterator sh_i = hit_V.begin(); sh_i != hit_V.end(); ++sh_i) {
0434 CSCDetId Sim_Id1 = (CSCDetId)(*sh_i).detUnitId();
0435 const CSCChamber* cscchamber = cscGeom->chamber(Sim_Id1);
0436
0437 double z1 = (cscchamber->layer(1))->position().z();
0438 double z6 = (cscchamber->layer(6))->position().z();
0439
0440 if (((z1 > z6) && (z1 > 0.)) || ((z1 < z6) && (z1 < 0.))) {
0441 rv_flag = true;
0442 }
0443 }
0444
0445
0446
0447
0448
0449
0450 double sum1[4] = {0.0};
0451 double sum2[4] = {0.0};
0452 par1[0] = 0.0;
0453 par1[1] = 0.0;
0454 par2[0] = 0.0;
0455 par2[1] = 0.0;
0456 double N = 0.0;
0457 CSCDetId Sim_Id;
0458 for (std::vector<PSimHit>::const_iterator sh_i = hit_V.begin(); sh_i != hit_V.end(); ++sh_i) {
0459 Sim_Id = (CSCDetId)(*sh_i).detUnitId();
0460 LocalPoint Sim_xyz = (*sh_i).localPosition();
0461
0462 const CSCLayer* csclayer = cscGeom->layer(Sim_Id);
0463 const CSCChamber* cscchamber = cscGeom->chamber(Sim_Id);
0464 GlobalPoint gsh_xyz = csclayer->toGlobal(Sim_xyz);
0465 LocalPoint lsh_xyz = cscchamber->toLocal(gsh_xyz);
0466 double zz = lsh_xyz.z();
0467
0468 if (rv_flag) {
0469 zz = (-1.0) * lsh_xyz.z();
0470 }
0471
0472 N += 1.0;
0473 sum1[0] += lsh_xyz.x();
0474 sum1[1] += zz;
0475 sum1[2] += lsh_xyz.x() * zz;
0476 sum1[3] += zz * zz;
0477
0478 sum2[0] += lsh_xyz.y();
0479 sum2[1] += zz;
0480 sum2[2] += lsh_xyz.y() * zz;
0481 sum2[3] += zz * zz;
0482 }
0483 par1[0] = ((sum1[0] * sum1[3]) - (sum1[2] * sum1[1])) / ((N * sum1[3]) - (sum1[1] * sum1[1]));
0484 par1[1] = ((N * sum1[2]) - (sum1[0] * sum1[1])) / ((N * sum1[3]) - (sum1[1] * sum1[1]));
0485
0486 par2[0] = ((sum2[0] * sum2[3]) - (sum2[2] * sum2[1])) / ((N * sum2[3]) - (sum2[1] * sum2[1]));
0487 par2[1] = ((N * sum2[2]) - (sum2[0] * sum2[1])) / ((N * sum2[3]) - (sum2[1] * sum2[1]));
0488
0489 double dzz = 1. / sqrt((par1[1] * par1[1]) + (par2[1] * par2[1]) + 1.0);
0490 double dxx = dzz * par1[1];
0491 double dyy = dzz * par2[1];
0492 LocalVector LV = LocalVector(dxx, dyy, dzz);
0493 LSimOrg = LocalPoint(par1[0], par2[0], 0.);
0494
0495 const CSCChamber* cscchamber = cscGeom->chamber(Sim_Id);
0496 GlobalVector GV = cscchamber->toGlobal(LV);
0497 GlobalPoint GP = cscchamber->toGlobal(LSimOrg);
0498 LSimVec = LV;
0499 GSimOrg = GP;
0500
0501 double directionSign = GP.z() * GV.z();
0502 LV = (directionSign * LV).unit();
0503 GV = cscchamber->toGlobal(LV);
0504 GSimVec = GV;
0505
0506 if ((Sim_Id.station() == 3) || (Sim_Id.station() == 4)) {
0507 GSimVec = GlobalVector(-1. * GV.x(), -1. * GV.y(), GV.z());
0508 } else {
0509 GSimVec = GV;
0510 }
0511
0512
0513
0514
0515
0516
0517
0518
0519
0520
0521
0522
0523
0524
0525
0526
0527
0528
0529
0530
0531 }