Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-05-03 02:57:00

0001 /** \file CSCSegmentReader.cc
0002  *
0003  *  \author M. Sani
0004  *
0005  *  Modified by D. Fortin - UC Riverside
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 // Constructor
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 // Destructor
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 // The real analysis
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;  // Chamber already used to determine efficiency
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       // Test if have 6 layers with simhits in chamber
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       // Test if have 6 layers with rechits in chamber
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               //          hchi2->Fill((*segIt).chi2());
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;  // Only plot resolution for certain segments
0303 
0304     // Look if have enough rechits for resolution plot
0305     CSCDetId recId = (*its).cscDetId();
0306 
0307     if (recId.ring() > 3)
0308       continue;  // ignore ME1/a
0309 
0310     int ith_layer = 0;
0311     int nLayerWithRechitsInChamber = 0;
0312     int nRecHitChamber = 0;
0313     // Test if have 6 layers with rechits in chamber
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;  // 10 centimeter radius
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     // Find angular difference between two segments.
0593     // Use dot product:  cos(theta_12) = v1 . v2 / [ |v1|*|v2| ]
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);