File indexing completed on 2024-04-06 12:26:02
0001
0002
0003
0004
0005
0006
0007
0008 #include <RecoLocalMuon/CSCSegment/test/CSCSegmentReader.h>
0009
0010 #include <FWCore/Framework/interface/MakerMacros.h>
0011 #include <Geometry/Records/interface/MuonGeometryRecord.h>
0012
0013 #include <DataFormats/CSCRecHit/interface/CSCRecHit2DCollection.h>
0014
0015 #include <DataFormats/CSCRecHit/interface/CSCRangeMapAccessor.h>
0016 #include <Geometry/CSCGeometry/interface/CSCChamber.h>
0017
0018
0019 CSCSegmentReader::CSCSegmentReader(const edm::ParameterSet& pset) {
0020 simhit = 0;
0021 near_segment = 0;
0022 filename = pset.getUntrackedParameter<std::string>("RootFileName");
0023 minLayerWithRechitChamber = pset.getUntrackedParameter<int>("minLayerWithRechitPerChamber");
0024 minLayerWithSimhitChamber = pset.getUntrackedParameter<int>("minLayerWithSimhitPerChamber");
0025 maxNhits = pset.getUntrackedParameter<int>("maxNhits");
0026 minNhits = pset.getUntrackedParameter<int>("minNhits");
0027 minRechitSegment = pset.getUntrackedParameter<int>("minRechitPerSegment");
0028 maxPhi = pset.getUntrackedParameter<double>("maxPhiSeparation");
0029 maxTheta = pset.getUntrackedParameter<double>("maxThetaSeparation");
0030
0031 geomToken_ = esConsumes();
0032 simTracksToken_ = consumes(edm::InputTag("g4SimHits"));
0033 simHitsToken_ = consumes(edm::InputTag("g4SimHits", "MuonCSCHits"));
0034 recHitsToken_ = consumes(edm::InputTag("csc2DRecHits"));
0035 cscSegmentsToken_ = consumes(edm::InputTag("cscSegments"));
0036
0037 file = new TFile(filename.c_str(), "RECREATE");
0038
0039 if (file->IsOpen())
0040 std::cout << "file open!" << std::endl;
0041 else
0042 std::cout << "*** Error in opening file ***" << std::endl;
0043
0044 hchi2 = new TH1F("h4", "chi2", 200, 0, 400);
0045 hrechit = new TH1I("h5", "nrechit", 6, 2, 8);
0046 hsegment = new TH1I("h6", "segments multiplicity", 20, 0, 20);
0047 heta = new TH1F("h7", "eta sim muons", 50, -2.5, 2.5);
0048 hpt = new TH1F("h8", "pT sim muons", 120, 0, 60);
0049 hx = new TH1F("h9", "deltaX", 400, -100, +100);
0050 hy = new TH1F("h10", "deltaY", 400, -100, +100);
0051
0052 char a[4];
0053 for (unsigned char i = 0; i < 9; i++) {
0054 snprintf(a, sizeof(a), "h2%hhu", i);
0055 hdxOri[i] = new TH1F(a, "#Delta X", 101, -5.05, 5.05);
0056 snprintf(a, sizeof(a), "h3%hhu", i);
0057 hdyOri[i] = new TH1F(a, "#Delta Y", 101, -5.05, 5.05);
0058 snprintf(a, sizeof(a), "h4%hhu", i);
0059 hphiDir[i] = new TH1F(a, "#Delta #phi", 101, -0.505, 0.505);
0060 snprintf(a, sizeof(a), "h5%hhu", i);
0061 hthetaDir[i] = new TH1F(a, "#Delta #theta", 101, -0.505, 0.505);
0062 }
0063 }
0064
0065
0066 CSCSegmentReader::~CSCSegmentReader() {
0067 int ibin = 0;
0068 heff0 = new TH1F("h0", "efficiency", segMap1.size() * 2 + 2, 0, segMap1.size() * 2 + 2);
0069 heff1 = new TH1F("h1", "efficiency", segMap1.size() * 2 + 2, 0, segMap1.size() * 2 + 2);
0070 heff2 = new TH1F("h2", "efficiency", segMap1.size() * 2 + 2, 0, segMap1.size() * 2 + 2);
0071 heff3 = new TH1F("h3", "efficiency", segMap1.size() * 2 + 2, 0, segMap1.size() * 2 + 2);
0072
0073 std::cout << "Raw reco efficiency for 6-hit simulated segment (nhits > 3)" << std::endl;
0074 for (std::map<std::string, int>::const_iterator it = segMap1.begin(); it != segMap1.end(); it++) {
0075 ibin++;
0076 float eff = (float)it->second / (float)chaMap1[it->first];
0077 heff0->SetBinContent(ibin * 2, eff);
0078 heff0->GetXaxis()->SetBinLabel(ibin * 2, (it->first).c_str());
0079 std::cout << it->first << ": " << it->second << " " << chaMap1[it->first] << " " << eff << std::endl;
0080 }
0081 ibin = 0;
0082 std::cout << "Raw reco efficiency for chamber with 6 layers with rechits (nhits > 3)" << std::endl;
0083 for (std::map<std::string, int>::const_iterator it = segMap2.begin(); it != segMap2.end(); it++) {
0084 ibin++;
0085 float eff = (float)it->second / (float)chaMap2[it->first];
0086 heff1->SetBinContent(ibin * 2, eff);
0087 heff1->GetXaxis()->SetBinLabel(ibin * 2, (it->first).c_str());
0088 std::cout << it->first << ": " << it->second << " " << chaMap2[it->first] << " " << eff << std::endl;
0089 }
0090 ibin = 0;
0091 std::cout << "Reco efficiency for building 6-hit segment for 6-hit simulated segment" << std::endl;
0092 for (std::map<std::string, int>::const_iterator it = segMap3.begin(); it != segMap3.end(); it++) {
0093 ibin++;
0094 float eff = (float)it->second / (float)chaMap1[it->first];
0095 heff2->SetBinContent(ibin * 2, eff);
0096 heff2->GetXaxis()->SetBinLabel(ibin * 2, (it->first).c_str());
0097 std::cout << it->first << ": " << it->second << " " << chaMap1[it->first] << " " << eff << std::endl;
0098 }
0099 ibin = 0;
0100 std::cout << "Reco efficiency for building 6-hit segment for chamber with 6 layers with rechits" << std::endl;
0101 for (std::map<std::string, int>::const_iterator it = segMap3.begin(); it != segMap3.end(); it++) {
0102 ibin++;
0103 float eff = (float)it->second / (float)chaMap2[it->first];
0104 heff3->SetBinContent(ibin * 2, eff);
0105 heff3->GetXaxis()->SetBinLabel(ibin * 2, (it->first).c_str());
0106 std::cout << it->first << ": " << it->second << " " << chaMap2[it->first] << " " << eff << std::endl;
0107 }
0108
0109 file->cd();
0110 heff1->Write();
0111 heff2->Write();
0112 heff3->Write();
0113 hchi2->Write();
0114 hrechit->Write();
0115 hsegment->Write();
0116 hpt->Write();
0117 heta->Write();
0118 hx->Write();
0119 hy->Write();
0120
0121 for (int i = 0; i < 9; i++) {
0122 hdxOri[i]->Write();
0123 hdyOri[i]->Write();
0124 hphiDir[i]->Write();
0125 hthetaDir[i]->Write();
0126 }
0127 file->Close();
0128 }
0129
0130
0131 void CSCSegmentReader::analyze(const edm::Event& event, const edm::EventSetup& eventSetup) {
0132 const CSCGeometry* geom = &eventSetup.getData(geomToken_);
0133
0134 edm::Handle<edm::SimTrackContainer> simTracks = event.getHandle(simTracksToken_);
0135
0136 edm::Handle<edm::PSimHitContainer> simHits = event.getHandle(simHitsToken_);
0137
0138 edm::Handle<CSCRecHit2DCollection> recHits = event.getHandle(recHitsToken_);
0139
0140 edm::Handle<CSCSegmentCollection> cscSegments = event.getHandle(cscSegmentsToken_);
0141
0142 simInfo(simTracks);
0143 resolution(simHits, recHits, cscSegments, geom);
0144 recInfo(simHits, recHits, cscSegments, geom);
0145 }
0146
0147 void CSCSegmentReader::recInfo(const edm::Handle<edm::PSimHitContainer> simHits,
0148 const edm::Handle<CSCRecHit2DCollection> recHits,
0149 const edm::Handle<CSCSegmentCollection> cscSegments,
0150 const CSCGeometry* geom) {
0151 hsegment->Fill(cscSegments->end() - cscSegments->begin());
0152
0153 std::vector<CSCDetId> cscChambers;
0154 for (edm::PSimHitContainer::const_iterator simIt = simHits->begin(); simIt != simHits->end(); simIt++) {
0155 bool usedChamber = false;
0156
0157 CSCDetId simId = (CSCDetId)(*simIt).detUnitId();
0158
0159 if (abs((*simIt).particleType()) != 13)
0160 continue;
0161
0162 unsigned sizeCh = cscChambers.size();
0163 if (sizeCh > 0) {
0164 for (unsigned i = 0; i < sizeCh; ++i) {
0165 if (simId == cscChambers[i]) {
0166 usedChamber = true;
0167 } else {
0168 cscChambers.push_back(simId);
0169 }
0170 }
0171 }
0172
0173 if (usedChamber)
0174 continue;
0175
0176 double g6 = geom->chamber(simId)->layer(6)->surface().toGlobal(LocalPoint(0, 0, 0)).z();
0177 double g1 = geom->chamber(simId)->layer(1)->surface().toGlobal(LocalPoint(0, 0, 0)).z();
0178
0179 int firstLayer = 1;
0180 if (fabs(g1) > fabs(g6))
0181 firstLayer = 6;
0182
0183 const CSCChamber* chamber = geom->chamber(simId);
0184
0185 if (simId.layer() == firstLayer) {
0186
0187 int ith_layer = 0;
0188 int nLayerWithSimhitsInChamber = 0;
0189
0190 for (edm::PSimHitContainer::const_iterator it2 = simHits->begin(); it2 != simHits->end(); it2++) {
0191 if (abs((*it2).particleType()) != 13)
0192 continue;
0193
0194 CSCDetId simId2 = (CSCDetId)(*it2).detUnitId();
0195 if ((simId2.chamber() == simId.chamber()) && (simId2.station() == simId.station()) &&
0196 (simId2.ring() == simId.ring()) && (simId2.endcap() == simId.endcap()) && (simId2.layer() != ith_layer)) {
0197 nLayerWithSimhitsInChamber++;
0198 ith_layer = simId2.layer();
0199 }
0200 }
0201
0202 if (nLayerWithSimhitsInChamber < minLayerWithSimhitChamber)
0203 continue;
0204
0205 bool satisfied0 = false;
0206 ith_layer = 0;
0207 int nLayerWithRechitsInChamber = 0;
0208 int nRecHitChamber = 0;
0209
0210 for (CSCRecHit2DCollection::const_iterator recIt = recHits->begin(); recIt != recHits->end(); recIt++) {
0211 CSCDetId idrec = (CSCDetId)(*recIt).cscDetId();
0212 if ((idrec.endcap() == simId.endcap()) && (idrec.ring() == simId.ring()) &&
0213 (idrec.station() == simId.station()) && (idrec.chamber() == simId.chamber()) &&
0214 (idrec.layer() != ith_layer)) {
0215 nLayerWithRechitsInChamber++;
0216 ith_layer = idrec.layer();
0217 } else if ((idrec.endcap() == simId.endcap()) && (idrec.ring() == simId.ring()) &&
0218 (idrec.station() == simId.station()) && (idrec.chamber() == simId.chamber())) {
0219 nRecHitChamber++;
0220 }
0221 }
0222
0223 if (simId.ring() < 4) {
0224 if (nRecHitChamber > maxNhits)
0225 continue;
0226 if (nRecHitChamber < minNhits)
0227 continue;
0228 } else {
0229 if (nRecHitChamber > 3 * maxNhits)
0230 continue;
0231 if (nRecHitChamber < 3 * minNhits)
0232 continue;
0233 }
0234
0235 std::string s = chamber->specs()->chamberTypeName();
0236
0237 chaMap1[s]++;
0238
0239 if (nLayerWithRechitsInChamber >= minLayerWithRechitChamber) {
0240 chaMap2[chamber->specs()->chamberTypeName()]++;
0241 satisfied0 = true;
0242 }
0243
0244 bool satisfied1 = false;
0245 bool satisfied2 = false;
0246
0247 int bestMatchIdx = bestMatch(simId, simHits, cscSegments, geom);
0248
0249 int idx = -1;
0250
0251 for (CSCSegmentCollection::const_iterator segIt = cscSegments->begin(); segIt != cscSegments->end(); ++segIt) {
0252 idx++;
0253
0254 CSCDetId id = (*segIt).cscDetId();
0255 if ((simId.endcap() == id.endcap()) && (simId.ring() == id.ring()) && (simId.station() == id.station()) &&
0256 (simId.chamber() == id.chamber())) {
0257 satisfied1 = true;
0258
0259 if (simId.ring() < 4)
0260 hrechit->Fill((*segIt).nRecHits());
0261
0262 if ((*segIt).nRecHits() >= minRechitSegment) {
0263 satisfied2 = true;
0264
0265 if (bestMatchIdx == idx && simId.ring() < 4) {
0266 hchi2->Fill(((*segIt).chi2() / (2 * (*segIt).nRecHits() - 4)));
0267
0268 }
0269 }
0270 }
0271 }
0272
0273 if (satisfied1)
0274 segMap1[s]++;
0275 if (satisfied1 && satisfied0)
0276 segMap2[s]++;
0277 if (satisfied2 && satisfied0)
0278 segMap3[s]++;
0279 }
0280 }
0281 }
0282
0283 void CSCSegmentReader::simInfo(const edm::Handle<edm::SimTrackContainer> simTracks) {
0284 for (edm::SimTrackContainer::const_iterator it = simTracks->begin(); it != simTracks->end(); it++) {
0285 if (abs((*it).type()) == 13) {
0286 hpt->Fill((*it).momentum().pt());
0287 heta->Fill((*it).momentum().eta());
0288 }
0289 }
0290 }
0291
0292 void CSCSegmentReader::resolution(const edm::Handle<edm::PSimHitContainer> simHits,
0293 const edm::Handle<CSCRecHit2DCollection> recHits,
0294 const edm::Handle<CSCSegmentCollection> cscSegments,
0295 const CSCGeometry* geom) {
0296 int idx = -1;
0297
0298 for (CSCSegmentCollection::const_iterator its = cscSegments->begin(); its != cscSegments->end(); its++) {
0299 idx++;
0300
0301 if ((*its).nRecHits() < minRechitSegment)
0302 continue;
0303
0304
0305 CSCDetId recId = (*its).cscDetId();
0306
0307 if (recId.ring() > 3)
0308 continue;
0309
0310 int ith_layer = 0;
0311 int nLayerWithRechitsInChamber = 0;
0312 int nRecHitChamber = 0;
0313
0314 for (CSCRecHit2DCollection::const_iterator recIt = recHits->begin(); recIt != recHits->end(); recIt++) {
0315 CSCDetId idrec = (CSCDetId)(*recIt).cscDetId();
0316 if ((idrec.endcap() == recId.endcap()) && (idrec.ring() == recId.ring()) &&
0317 (idrec.station() == recId.station()) && (idrec.chamber() == recId.chamber()) &&
0318 (idrec.layer() != ith_layer)) {
0319 nLayerWithRechitsInChamber++;
0320 ith_layer = idrec.layer();
0321 } else if ((idrec.endcap() == recId.endcap()) && (idrec.ring() == recId.ring()) &&
0322 (idrec.station() == recId.station()) && (idrec.chamber() == recId.chamber())) {
0323 nRecHitChamber++;
0324 }
0325 }
0326
0327 if (nLayerWithRechitsInChamber < minLayerWithRechitChamber)
0328 continue;
0329
0330 if (recId.ring() < 4) {
0331 if (nRecHitChamber > maxNhits)
0332 continue;
0333 if (nRecHitChamber < minNhits)
0334 continue;
0335 } else {
0336 if (nRecHitChamber > 3 * maxNhits)
0337 continue;
0338 if (nRecHitChamber < 3 * minNhits)
0339 continue;
0340 }
0341
0342 int bestMatchIdx = bestMatch(recId, simHits, cscSegments, geom);
0343
0344 if (bestMatchIdx != idx)
0345 continue;
0346
0347 double segX = -99999., segY = -99999.;
0348 unsigned int simTrack = 0;
0349
0350 const CSCChamber* ch = geom->chamber((*its).cscDetId());
0351
0352 LocalVector segDir = (*its).localDirection();
0353 LocalPoint segLoc = (*its).localPosition();
0354
0355 segX = segLoc.x();
0356 segY = segLoc.y();
0357
0358 GlobalPoint gpn = ch->layer(6)->surface().toGlobal(LocalPoint(0, 0, 0));
0359 GlobalPoint gpf = ch->layer(1)->surface().toGlobal(LocalPoint(0, 0, 0));
0360
0361 int firstLayer = 6;
0362 int lastLayer = 1;
0363
0364 if (fabs(gpn.z()) > fabs(gpf.z())) {
0365 firstLayer = 1;
0366 lastLayer = 6;
0367 }
0368
0369 bool check1 = false;
0370
0371 int counter = 0;
0372
0373 float sim1X = 0.;
0374 float sim1Y = 0.;
0375 double sim1Phi = 0.;
0376 double sim1Theta = 0.;
0377
0378 for (edm::PSimHitContainer::const_iterator ith = simHits->begin(); ith != simHits->end(); ith++) {
0379 if (abs((*ith).particleType()) != 13)
0380 continue;
0381
0382 CSCDetId simId = (CSCDetId)(*ith).detUnitId();
0383
0384 if (ch != geom->chamber(simId))
0385 continue;
0386
0387 if ((simId.layer() == firstLayer)) {
0388 check1 = true;
0389
0390 LocalVector simDir1_temp = (*ith).momentumAtEntry().unit();
0391
0392 sim1Phi += simDir1_temp.phi();
0393 sim1Theta += simDir1_temp.theta();
0394
0395 sim1X += (*ith).localPosition().x();
0396 sim1Y += (*ith).localPosition().y();
0397
0398 counter++;
0399 }
0400 }
0401
0402 if (!check1)
0403 continue;
0404
0405 sim1X = sim1X / counter;
0406 sim1Y = sim1Y / counter;
0407 sim1Phi = sim1Phi / counter;
0408 sim1Theta = sim1Theta / counter;
0409
0410 double dPhi1 = sim1Phi - segDir.phi();
0411 double dTheta1 = sim1Theta - segDir.theta();
0412
0413 bool check6 = false;
0414
0415 float sim2X = 0.;
0416 float sim2Y = 0.;
0417 double sim2Phi = 0.;
0418 double sim2Theta = 0.;
0419
0420 counter = 0;
0421
0422 for (edm::PSimHitContainer::const_iterator ith = simHits->begin(); ith != simHits->end(); ith++) {
0423 if (abs((*ith).particleType()) != 13)
0424 continue;
0425
0426 CSCDetId simId = (CSCDetId)(*ith).detUnitId();
0427
0428 if (ch != geom->chamber(simId))
0429 continue;
0430
0431 if ((simId.layer() == lastLayer) && (simTrack = (*ith).trackId())) {
0432 check6 = true;
0433
0434 LocalVector simDir6_temp = (*ith).momentumAtEntry().unit();
0435
0436 sim2Phi += simDir6_temp.phi();
0437 sim2Theta += simDir6_temp.theta();
0438
0439 sim2X += (*ith).localPosition().x();
0440 sim2Y += (*ith).localPosition().y();
0441 counter++;
0442 }
0443 }
0444
0445 if (!check6)
0446 continue;
0447
0448 sim2X = sim2X / counter;
0449 sim2Y = sim2Y / counter;
0450 sim2Phi = sim2Phi / counter;
0451 sim2Theta = sim2Theta / counter;
0452
0453 double dPhi2 = sim2Phi - segDir.phi();
0454 double dTheta2 = sim2Theta - segDir.theta();
0455
0456 double deltaTheta = (dTheta1 + dTheta2) / 2.;
0457 double deltaPhi = (dPhi1 + dPhi2) / 2.;
0458
0459 float simX = (sim1X + sim2X) / 2.;
0460 float simY = (sim1Y + sim2Y) / 2.;
0461
0462 float deltaX = segX - simX;
0463 float deltaY = segY - simY;
0464
0465 int indice = 0;
0466
0467 hdxOri[indice]->Fill(deltaX);
0468 hdyOri[indice]->Fill(deltaY);
0469 hphiDir[indice]->Fill(deltaPhi);
0470 hthetaDir[indice]->Fill(deltaTheta);
0471 }
0472 }
0473
0474 int CSCSegmentReader::bestMatch(CSCDetId id0,
0475 const edm::Handle<edm::PSimHitContainer> simHits,
0476 const edm::Handle<CSCSegmentCollection> cscSegments,
0477 const CSCGeometry* geom) {
0478 int bestIndex = -1;
0479
0480 const CSCChamber* ch = geom->chamber(id0);
0481
0482 bool check1 = false;
0483 bool check6 = false;
0484
0485 GlobalPoint gpn = ch->layer(6)->surface().toGlobal(LocalPoint(0, 0, 0));
0486 GlobalPoint gpf = ch->layer(1)->surface().toGlobal(LocalPoint(0, 0, 0));
0487
0488 int firstLayer = 6;
0489 int lastLayer = 1;
0490
0491 if (fabs(gpn.z()) > fabs(gpf.z())) {
0492 firstLayer = 1;
0493 lastLayer = 6;
0494 }
0495
0496 LocalVector simDir1;
0497
0498 float sim1X = 0.;
0499 float sim1Y = 0.;
0500 float sim1Z = 0.;
0501 int counter = 0;
0502
0503 for (edm::PSimHitContainer::const_iterator ith = simHits->begin(); ith != simHits->end(); ith++) {
0504 if (abs((*ith).particleType()) != 13)
0505 continue;
0506
0507 CSCDetId simId = (CSCDetId)(*ith).detUnitId();
0508
0509 if (simId.layer() == firstLayer && ch == geom->chamber(simId)) {
0510 check1 = true;
0511
0512 sim1X += (*ith).localPosition().x();
0513 sim1Y += (*ith).localPosition().y();
0514 sim1Z = (*ith).localPosition().z();
0515
0516 counter++;
0517 }
0518 }
0519
0520 if (!check1)
0521 return bestIndex;
0522 sim1X = sim1X / counter;
0523 sim1Y = sim1Y / counter;
0524
0525 float sim2X = 0.;
0526 float sim2Y = 0.;
0527 float sim2Z = 0.;
0528 counter = 0;
0529
0530 for (edm::PSimHitContainer::const_iterator ith = simHits->begin(); ith != simHits->end(); ith++) {
0531 if (abs((*ith).particleType()) != 13)
0532 continue;
0533
0534 CSCDetId simId = (CSCDetId)(*ith).detUnitId();
0535 if (simId.layer() == lastLayer && ch == geom->chamber(simId)) {
0536 check6 = true;
0537
0538 sim2X += (*ith).localPosition().x();
0539 sim2Y += (*ith).localPosition().y();
0540 sim2Z = (*ith).localPosition().z();
0541 counter++;
0542 }
0543 }
0544
0545 if (!check6)
0546 return bestIndex;
0547
0548 sim2X = sim2X / counter;
0549 sim2Y = sim2Y / counter;
0550
0551 float x = (sim2X + sim1X) / 2.;
0552 float y = (sim2Y + sim1Y) / 2.;
0553
0554 float dx = (sim2X - sim1X);
0555 float dy = (sim2Y - sim1Y);
0556 float dz = (sim2Z - sim1Z);
0557 float magSim = sqrt(dx * dx + dy * dy + dz * dz);
0558
0559 int idxCounter = -1;
0560 float bestCosTheta = 0.;
0561
0562 for (CSCSegmentCollection::const_iterator its = cscSegments->begin(); its != cscSegments->end(); its++) {
0563 CSCDetId id = (*its).cscDetId();
0564
0565 if ((id0.endcap() == id.endcap()) && (id0.ring() == id.ring()) && (id0.station() == id.station()) &&
0566 (id0.chamber() == id.chamber())) {
0567 idxCounter++;
0568 } else {
0569 idxCounter++;
0570 continue;
0571 }
0572
0573 if ((*its).nRecHits() < minRechitSegment)
0574 continue;
0575
0576 float xreco = (*its).localPosition().x();
0577 float yreco = (*its).localPosition().y();
0578
0579 float dR = (xreco - x) * (xreco - x) + (yreco - y) * (yreco - y);
0580 dR = sqrt(dR);
0581
0582 if (dR > 10.)
0583 continue;
0584
0585 LocalVector segDir = (*its).localDirection();
0586
0587 float xdir = segDir.x();
0588 float ydir = segDir.y();
0589 float zdir = segDir.z();
0590 float magReco = sqrt(xdir * xdir + ydir * ydir + zdir * zdir);
0591
0592
0593
0594
0595 double costheta = (xdir * dx + ydir * dy + zdir * dz) / (magSim * magReco);
0596 if (costheta < 0.)
0597 costheta = -costheta;
0598
0599 if (costheta > bestCosTheta) {
0600 bestCosTheta = costheta;
0601 bestIndex = idxCounter;
0602 }
0603 }
0604
0605 return bestIndex;
0606 }
0607
0608 DEFINE_FWK_MODULE(CSCSegmentReader);
