Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:27:12 

0001 // Class Header
0002 #include "MuonSeedValidator.h"
0003 #include "SegSelector.h"
0004 
0005 // for MuonSeedBuilder
0006 #include "RecoMuon/DetLayers/interface/MuonDetLayerGeometry.h"
0007 #include "RecoMuon/Records/interface/MuonRecoGeometryRecord.h"
0008 //#include "MagneticField/Engine/interface/MagneticField.h"
0009 //#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0010 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0011 
0012 // Framework
0013 #include "FWCore/Framework/interface/ConsumesCollector.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/EventSetup.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018 #include "FWCore/Utilities/interface/Exception.h"
0019 
0020 //#include "PluginManager/ModuleDef.h"
0021 #include "FWCore/Framework/interface/MakerMacros.h"
0022 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0023 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0024 
0025 #include "TFile.h"
0026 #include "TVector3.h"
0027 
0028 #include <iostream>
0029 #include <fstream>
0030 #include <map>
0031 #include <utility>
0032 #include <string>
0033 #include <stdio.h>
0034 #include <algorithm>
0035 
0036 DEFINE_FWK_MODULE(MuonSeedValidator);
0037 using namespace std;
0038 using namespace edm;
0039 using namespace reco;
0040 
0041 // constructors
0042 MuonSeedValidator::MuonSeedValidator(const ParameterSet& pset) : cscGeomToken(esConsumes()), dtGeomToken(esConsumes()) {
0043   rootFileName = pset.getUntrackedParameter<string>("rootFileName");
0044   recHitLabel = pset.getUntrackedParameter<string>("recHitLabel");
0045   cscSegmentLabel = pset.getUntrackedParameter<string>("cscSegmentLabel");
0046   dtrecHitLabel = pset.getUntrackedParameter<string>("dtrecHitLabel");
0047   dtSegmentLabel = pset.getUntrackedParameter<string>("dtSegmentLabel");
0048   simHitLabel = pset.getUntrackedParameter<string>("simHitLabel");
0049   simTrackLabel = pset.getUntrackedParameter<string>("simTrackLabel");
0050   muonseedLabel = pset.getUntrackedParameter<string>("muonseedLabel");
0051   staTrackLabel = pset.getParameter<InputTag>("staTrackLabel");
0052   glbTrackLabel = pset.getParameter<InputTag>("glbTrackLabel");
0053 
0054   debug = pset.getUntrackedParameter<bool>("debug");
0055   dtMax = pset.getUntrackedParameter<double>("dtMax");
0056   dfMax = pset.getUntrackedParameter<double>("dfMax");
0057   scope = pset.getUntrackedParameter<bool>("scope");
0058   pTCutMax = pset.getUntrackedParameter<double>("pTCutMax");
0059   pTCutMin = pset.getUntrackedParameter<double>("pTCutMin");
0060   eta_Low = pset.getUntrackedParameter<double>("eta_Low");
0061   eta_High = pset.getUntrackedParameter<double>("eta_High");
0062 
0063   ParameterSet serviceParameters = pset.getParameter<ParameterSet>("ServiceParameters");
0064   theService = new MuonServiceProxy(serviceParameters, consumesCollector());
0065 
0066   recsegSelector = new SegSelector(pset);
0067 
0068   // Create the root file
0069   theFile = new TFile(rootFileName.c_str(), "RECREATE");
0070   theFile->mkdir("AllMuonSys");
0071   theFile->cd();
0072   theFile->mkdir("No_Seed");
0073   theFile->cd();
0074   theFile->mkdir("No_STA");
0075   theFile->cd();
0076   theFile->mkdir("EventScope");
0077   theFile->cd();
0078   theFile->mkdir("UnRelated");
0079   theFile->cd();
0080   // TTree test
0081 
0082   h_all = new H2DRecHit1("AllMu_");
0083   h_NoSeed = new H2DRecHit2("NoSeed");
0084   h_NoSta = new H2DRecHit3("NoSta");
0085   h_Scope = new H2DRecHit4();
0086   h_UnRel = new H2DRecHit5();
0087 }
0088 
0089 // destructor
0090 MuonSeedValidator::~MuonSeedValidator() {
0091   if (debug)
0092     cout << "[Seed Validation] Destructor called" << endl;
0093   delete recsegSelector;
0094   //delete muonSeedBuilder_;
0095   // Write the histos to file
0096   theFile->cd();
0097   theFile->cd("AllMuonSys");
0098   h_all->Write();
0099 
0100   theFile->cd();
0101   theFile->cd("No_Seed");
0102   h_NoSeed->Write();
0103 
0104   theFile->cd();
0105   theFile->cd("No_STA");
0106   h_NoSta->Write();
0107 
0108   theFile->cd();
0109   theFile->cd("EventScope");
0110   h_Scope->Write();
0111 
0112   theFile->cd();
0113   theFile->cd("UnRelated");
0114   h_UnRel->Write();
0115   // for tree
0116   theFile->cd();
0117 
0118   // Release the memory ...
0119   delete h_all;
0120   delete h_NoSeed;
0121   delete h_NoSta;
0122   delete h_Scope;
0123   delete h_UnRel;
0124 
0125   theFile->Close();
0126   if (debug)
0127     cout << "************* Finished writing histograms to file" << endl;
0128   if (theService)
0129     delete theService;
0130 }
0131 
0132 // The Main...Aanlysis...
0133 
0134 void MuonSeedValidator::analyze(const Event& event, const EventSetup& eventSetup) {
0135   //Get the CSC Geometry :
0136   theService->update(eventSetup);
0137 
0138   //ESHandle<GlobalTrackingGeometry> globalGeometry;
0139   //eventSetup.get<GlobalTrackingGeometryRecord>().get(globalGeometry);
0140 
0141   ESHandle<CSCGeometry> cscGeom = eventSetup.getHandle(cscGeomToken);
0142   ESHandle<DTGeometry> dtGeom = eventSetup.getHandle(dtGeomToken);
0143 
0144   // Get the RecHits collection :
0145   Handle<CSCRecHit2DCollection> csc2DRecHits;
0146   event.getByLabel(recHitLabel, csc2DRecHits);
0147 
0148   // Get the CSC Segments collection :
0149   Handle<CSCSegmentCollection> cscSegments;
0150   event.getByLabel(cscSegmentLabel, cscSegments);
0151 
0152   // Get the DT RecHits collection :
0153   Handle<DTRecHitCollection> dt1DRecHits;
0154   event.getByLabel(dtrecHitLabel, dt1DRecHits);
0155 
0156   // Get the DT Segments collection :
0157   Handle<DTRecSegment4DCollection> dt4DSegments;
0158   event.getByLabel(dtSegmentLabel, dt4DSegments);
0159 
0160   // Get the SimHit collection :
0161   Handle<PSimHitContainer> csimHits;
0162   event.getByLabel(simHitLabel, "MuonCSCHits", csimHits);
0163   Handle<PSimHitContainer> dsimHits;
0164   event.getByLabel(simHitLabel, "MuonDTHits", dsimHits);
0165 
0166   // Get the SimTrack
0167   Handle<SimTrackContainer> simTracks;
0168   event.getByLabel(simTrackLabel, simTracks);
0169 
0170   // Get the muon seeds
0171   Handle<TrajectorySeedCollection> muonseeds;
0172   event.getByLabel(muonseedLabel, muonseeds);
0173   // Get sta muon tracks
0174   Handle<TrackCollection> standAloneMuons;
0175   event.getByLabel(staTrackLabel, standAloneMuons);
0176   // Get global muon tracks
0177   Handle<TrackCollection> globalMuons;
0178   event.getByLabel(glbTrackLabel, globalMuons);
0179 
0180   // Magnetic field
0181   //ESHandle<MagneticField> field;
0182   //eventSetup.get<IdealMagneticFieldRecord>().get(field);
0183 
0184   H2DRecHit1* histo1 = 0;
0185   H2DRecHit2* histo2 = 0;
0186   H2DRecHit3* histo3 = 0;
0187   H2DRecHit4* histo4 = 0;
0188   H2DRecHit5* histo5 = 0;
0189 
0190   // Get sim track information
0191   // return  theta_v, theta_p, phi_v, phi_p :  theta and phi of position and momentum
0192   //                                           of first simhit from CSC and DT
0193   //         theQ, eta_trk, phi_trk, pt_trk, trackId : Q,eta,phi,pt and trackID of simtrack
0194   //         pt_layer,palayer : the pt or pa in each layer
0195   SimInfo(simTracks, dsimHits, csimHits, dtGeom, cscGeom);
0196 
0197   /// A. statistic information for segment, seed and sta muon w.r.t eta
0198   // 2. Reading the damn seeds
0199   RecSeedReader(muonseeds);
0200 
0201   // 3. Get muon track information ; sta=0, glb=1
0202   StaTrackReader(standAloneMuons, 0);
0203 
0204   // 4. Associate seeds with sim tracks
0205   //    seed_trk => the simtrack number j associate with seed i
0206   std::vector<int> seeds_simtrk(seed_gp.size(), -1);
0207   for (unsigned int i = 0; i < seed_gp.size(); i++) {
0208     double dR1 = 99.0;
0209     int preferTrack = -1;
0210     for (unsigned int j = 0; j < theta_p.size(); j++) {
0211       double dt = fabs(seed_gp[i].theta() - theta_p[j]);
0212       double df = fabs(seed_gp[i].phi() - phi_p[j]);
0213       if (df > (6.283 - dfMax))
0214         df = 6.283 - df;
0215       double dR2 = sqrt(dt * dt + df * df);
0216 
0217       if (dR2 < dR1 && dt < dtMax && df < dfMax) {
0218         preferTrack = static_cast<int>(j);
0219         dR1 = dR2;
0220       }
0221     }
0222     seeds_simtrk[i] = preferTrack;
0223   }
0224 
0225   // 5. Associate sta with sim tracks
0226   //    sta_simtrk => the simtrack number j associate with sta  i
0227   std::vector<int> sta_simtrk(sta_thetaV.size(), -1);
0228   for (unsigned int i = 0; i < sta_thetaV.size(); i++) {
0229     double dR1 = 99.0;
0230     int preferTrack = -1;
0231     for (unsigned int j = 0; j < theta_p.size(); j++) {
0232       double dt = fabs(sta_thetaP[i] - theta_p[j]);
0233       double df = fabs(sta_phiP[i] - phi_p[j]);
0234       if (df > (6.283 - dfMax))
0235         df = 6.283 - df;
0236       double dR2 = sqrt(dt * dt + df * df);
0237       if (dR2 < dR1 && dt < dtMax && df < dfMax) {
0238         preferTrack = static_cast<int>(j);
0239         dR1 = dR2;
0240       }
0241     }
0242     sta_simtrk[i] = preferTrack;
0243   }
0244 
0245   // look at the number of reconstructed seed/sta for each sim-track
0246   histo1 = h_all;
0247   int sta_Evt = 0;
0248   int seed_Evt = 0;
0249   std::vector<int> nuSTA;
0250   std::vector<int> nuSEED;
0251   for (size_t i = 0; i < theta_p.size(); i++) {
0252     // 1. statistic for seeds
0253     int nu_seeds_trk = 0;
0254     for (size_t j = 0; j < seeds_simtrk.size(); j++) {
0255       if (seeds_simtrk[j] == static_cast<int>(i))
0256         nu_seeds_trk++;
0257       if (nu_seeds_trk > 19)
0258         nu_seeds_trk = 19;
0259     }
0260 
0261     // 2. statistic for sta
0262     int nu_sta_trk = 0;
0263     for (size_t j = 0; j < sta_simtrk.size(); j++) {
0264       if (sta_simtrk[j] == static_cast<int>(i))
0265         nu_sta_trk++;
0266     }
0267     std::vector<double> pa_tmp = palayer[i];
0268     histo1->Fill1b(getEta(theta_p[i]), nu_seeds_trk, nu_sta_trk, pt_trk[i], theQ[i] / pt_trk[i], theQ[i] / pa_tmp[0]);
0269 
0270     nuSTA.push_back(nu_sta_trk);
0271     nuSEED.push_back(nu_seeds_trk);
0272 
0273     if (nu_seeds_trk > 0)
0274       seed_Evt++;
0275     if (nu_sta_trk > 0)
0276       sta_Evt++;
0277 
0278     // 3. statistics for segment
0279     std::vector<SimSegment> sCSC_v = recsegSelector->Sim_CSCSegments(trackID[i], csimHits, cscGeom);
0280     std::vector<SimSegment> sDT_v = recsegSelector->Sim_DTSegments(trackID[i], dsimHits, dtGeom);
0281     std::vector<CSCSegment> cscseg_V = recsegSelector->Select_CSCSeg(cscSegments, cscGeom, sCSC_v);
0282     std::vector<DTRecSegment4D> dtseg_V = recsegSelector->Select_DTSeg(dt4DSegments, dtGeom, sDT_v);
0283 
0284     // 4. check # of segments and rechits in each chambers for this track
0285     CSCsegment_stat(cscSegments, cscGeom, theta_p[i], phi_p[i]);
0286     DTsegment_stat(dt4DSegments, dtGeom, theta_p[i], phi_p[i]);
0287     Simsegment_stat(sCSC_v, sDT_v);
0288     // 5. if less than 1 segment, check the # of rechits in each chambers for this track
0289     if ((cscseg_stat[5] < 2) && (dtseg_stat[5] < 2)) {
0290       CSCRecHit_Stat(csc2DRecHits, cscGeom, getEta(theta_p[i]), phi_p[i]);
0291       DTRecHit_Stat(dt1DRecHits, dtGeom, getEta(theta_p[i]), phi_p[i]);
0292     }
0293 
0294     // 6. look at all information
0295     // seg_stat[0] = total segments in all stations
0296     // seg_stat[5] = the number of stations which have segments
0297     // seg_stat1[x] = the number of stations/segments which only count segments w/ more than 4 rechits
0298     int layer_sum = cscseg_stat[5] + dtseg_stat[5];
0299     int goodlayer_sum = cscseg_stat1[5] + dtseg_stat1[5];
0300     int seg_sum = cscseg_stat[0] + dtseg_stat[0];
0301     int leftrh_sum = cscrh_sum[5] + dtrh_sum[5];
0302     std::vector<double> pa1 = palayer[i];
0303     histo1->Fill1(layer_sum, goodlayer_sum, seg_sum, simseg_sum, getEta(theta_p[i]));
0304     histo1->Fill1c(pt_trk[i], pa1[0], cscseg_stat[0] + dtseg_stat[0]);
0305 
0306     // 7. look at those events without any reco segments => no seed can be produced!!
0307     if ((cscseg_stat[5] == 0) && (dtseg_stat[5] == 0)) {
0308       histo1->Fill1a(leftrh_sum, getEta(theta_p[i]));
0309       if (cscrh_sum[0] != 0) {
0310         histo2 = h_NoSeed;
0311         histo2->Fill2b(getEta(theta_p[i]), cscrh_sum[0]);
0312       }
0313       if (dtrh_sum[0] != 0) {
0314         histo2 = h_NoSeed;
0315         histo2->Fill2c(getEta(theta_p[i]), dtrh_sum[0]);
0316       }
0317     }
0318     bool hasSeed = false;
0319     for (unsigned int j = 0; j < seeds_simtrk.size(); j++) {
0320       if (seeds_simtrk[j] == static_cast<int>(i))
0321         hasSeed = true;
0322     }
0323     if (!hasSeed) {
0324       histo2 = h_NoSeed;
0325       histo2->Fill2a(getEta(theta_p[i]), layer_sum, seg_sum, simseg_sum, simseg_sum - seg_sum, simseg_sum - layer_sum);
0326     }
0327     // 8. Read out reco segments
0328     //    return : the ave_phi, ave_eta, phi_resid, eta_resid, dx_error, dy_error, x_error, y_error
0329     int types = RecSegReader(cscSegments, dt4DSegments, cscGeom, dtGeom, theta_p[i], phi_p[i]);
0330     if (debug)
0331       cout << " seed type = " << types << endl;
0332   }
0333   histo1->Fill1o(sta_mT.size(), seed_mT.size());
0334 
0335   /// showering study
0336   for (unsigned int i = 0; i < theta_p.size(); i++) {
0337     if (theQ[i] < 0.0)
0338       continue;
0339     std::vector<int> showers = IdentifyShowering(cscSegments, cscGeom, dt4DSegments, dtGeom, theta_p[i], phi_p[i]);
0340 
0341     double maxR = 0.;
0342     for (size_t j = 0; j < muCone.size(); j++) {
0343       if (muCone[j] > maxR)
0344         maxR = muCone[j];
0345     }
0346 
0347     double aveShower = 0;
0348     if (maxR > 0) {
0349       aveShower = static_cast<double>(showers[0]) / maxR;
0350     }
0351 
0352     if (fabs(theta_p[i]) > 0.78)
0353       histo1->Fill1d1(pt_trk[i], showers[0], showers[1], showers[2], aveShower, maxR);
0354     if (fabs(theta_p[i]) <= 0.78 && fabs(theta_p[i]) > 0.36)
0355       histo1->Fill1d2(pt_trk[i], showers[0], showers[1], showers[2], aveShower, maxR);
0356     if (fabs(theta_p[i]) <= 0.36)
0357       histo1->Fill1d3(pt_trk[i], showers[0], showers[1], showers[2], aveShower, maxR);
0358   }
0359 
0360   // look at those un-associated seeds and sta
0361   histo5 = h_UnRel;
0362   for (size_t i = 0; i < seeds_simtrk.size(); i++) {
0363     if (seeds_simtrk[i] != -1)
0364       continue;
0365     // un-related !
0366     if (seed_Evt != 0) {
0367       histo5->Fill5a(seed_gp[i].eta(), seed_mT[i]);
0368     }
0369     // Orphan
0370     else {
0371       for (size_t j = 0; j < theta_p.size(); j++) {
0372         histo5->Fill5b(seed_gp[i].eta(), getEta(theta_p[j]), seed_mT[i]);
0373       }
0374     }
0375   }
0376   for (size_t i = 0; i < sta_simtrk.size(); i++) {
0377     if (sta_simtrk[i] != -1)
0378       continue;
0379     // un-related !
0380     if (sta_Evt != 0) {
0381       histo5->Fill5c(getEta(sta_thetaV[i]), sta_mT[i]);
0382     }
0383     // Orphan
0384     else {
0385       for (unsigned int j = 0; j < theta_p.size(); j++) {
0386         histo5->Fill5d(getEta(sta_thetaV[i]), getEta(theta_p[j]), sta_mT[i], sta_phiV[i], phi_v[j]);
0387       }
0388     }
0389   }
0390 
0391   /// B. seeds  Q/PT pulls
0392   if (nu_seed > 0) {
0393     // look the seed for every sta tracks
0394     for (unsigned int i = 0; i < eta_trk.size(); i++) {
0395       int bestSeed = -1;
0396       double dSeedPT = 99999.9;
0397 
0398       // get the correponding sim p and sim pt
0399       std::vector<double> pa1 = palayer[i];
0400       std::vector<double> pt1 = ptlayer[i];
0401 
0402       // find the best seed whose pt is closest to a simTrack which associate with
0403       for (unsigned int j = 0; j < seeds_simtrk.size(); j++) {
0404         if ((seeds_simtrk[j] == static_cast<int>(i)) && (fabs(seed_mT[j] - pt1[1]) < dSeedPT)) {
0405           dSeedPT = fabs(seed_mT[j] - pt1[1]);
0406           bestSeed = static_cast<int>(j);
0407         }
0408       }
0409 
0410       for (unsigned int j = 0; j < seeds_simtrk.size(); j++) {
0411         if (seeds_simtrk[j] != static_cast<int>(i))
0412           continue;
0413 
0414         // put the rest seeds in the minus side
0415         double bestSeedPt = 99999.9;
0416         if (bestSeed == static_cast<int>(j)) {
0417           bestSeedPt = seed_mT[j];
0418 
0419           // look at pull and resolution
0420           double pull_qbp = (qbp[j] - (theQ[i] / pa1[seed_layer[j]])) / err_qbp[j];
0421           double pull_qbpt = (qbpt[j] - (theQ[i] / pt1[seed_layer[j]])) / err_qbpt[j];
0422           //double resol_qbpt = ( qbpt[j] - (theQ[j]/pt_trk[j]) ) / (theQ[j]/pt_trk[j]) ;
0423           double resol_qbpt = (qbpt[j] - (theQ[i] / pt1[1])) / (theQ[i] / pt1[1]);
0424           //cout<<"resol = "<<resol_qbpt <<" = "<<qbpt[j] <<" - "<< theQ[j]/pt1[1]
0425           histo1->Fill1i(pull_qbp, seed_gp[j].eta(), qbpt[j], pull_qbpt, err_qbp[j], err_qbpt[j], resol_qbpt);
0426 
0427         } else {
0428           bestSeedPt = -1.0 * seed_mT[j];
0429         }
0430 
0431         double ptLoss = pt1[seed_layer[j]] / pt1[0];
0432         histo1->Fill1g(seed_mT[j], seed_mA[j], bestSeedPt, seed_gp[j].eta(), ptLoss, pt1[0]);
0433         histo1->Fill1f(seed_gp[j].eta(), err_dx[j], err_dy[j], err_x[j], err_y[j]);
0434       }
0435     }
0436   }
0437 
0438   /// C. sta track information
0439   if (nu_sta > 0) {
0440     for (unsigned int i = 0; i < eta_trk.size(); i++) {
0441       double dPT = 999999.;
0442       int best = 0;
0443       double expectPT = -1.0;
0444       // find the best sta whose pt is closest to a simTrack which associate with
0445       for (unsigned int j = 0; j < sta_simtrk.size(); j++) {
0446         if ((sta_simtrk[j] == static_cast<int>(i)) && (fabs(sta_mT[j] - pt_trk[i]) < dPT)) {
0447           std::vector<double> pt1 = ptlayer[i];
0448           dPT = fabs(sta_mT[j] - pt_trk[i]);
0449           //dPT = fabs( sta_mT[j] - pt1[1] );
0450           best = j;
0451           //expectPT = pt1[1];
0452           expectPT = pt_trk[i];
0453         }
0454       }
0455       // looking for the sta muon which is closest to simTrack pt
0456       if (expectPT > 0) {
0457         double sim_qbpt = theQ[i] / expectPT;
0458         double resol_qbpt = (sta_qbpt[best] - sim_qbpt) / sim_qbpt;
0459         /*
0460           double sta_q = 1.0;
0461           if ( sta_qbpt[best] < 0 ) sta_q = -1.0;
0462           double sta_qbmt = sta_q / sta_mT[best];
0463           double resol_qbpt = (sta_qbmt - sim_qbpt ) / sim_qbpt ;
0464           */
0465         histo1->Fill1j(
0466             getEta(theta_p[i]), sta_qbp[best], sta_qbpt[best], sta_mT[best], sta_mA[best], pt_trk[i], resol_qbpt);
0467       }
0468     }
0469 
0470     // look at the dPhi, dEta, dx, dy from sim-segment and reco-segment of a seed
0471     for (unsigned int i = 0; i < nSegInSeed.size(); i++) {
0472       int ii = seeds_simtrk[i];
0473       std::vector<SimSegment> sCSC_v = recsegSelector->Sim_CSCSegments(trackID[ii], csimHits, cscGeom);
0474       std::vector<SimSegment> sDT_v = recsegSelector->Sim_DTSegments(trackID[ii], dsimHits, dtGeom);
0475 
0476       SegOfRecSeed(muonseeds, i, sCSC_v, sDT_v);
0477       for (unsigned int j = 0; j < geoID.size(); j++) {
0478         if (geoID[j].subdetId() == MuonSubdetId::DT) {
0479           DTChamberId MB_Id = DTChamberId(geoID[j]);
0480           histo1->Fill1e(-1 * MB_Id.station(), d_h[j], d_f[j], d_x[j], d_y[j]);
0481         }
0482         if (geoID[j].subdetId() == MuonSubdetId::CSC) {
0483           CSCDetId ME_Id = CSCDetId(geoID[j]);
0484           histo1->Fill1e(ME_Id.station(), d_h[j], d_f[j], d_x[j], d_y[j]);
0485         }
0486       }
0487     }
0488   }
0489 
0490   /// open a scope to check all information
0491   if (scope) {
0492     cout << " ************ pT-eta scope *************** " << endl;
0493   }
0494   for (size_t i = 0; i < theta_p.size(); i++) {
0495     if (fabs(getEta(theta_p[i])) < eta_Low || fabs(getEta(theta_p[i])) > eta_High)
0496       continue;
0497 
0498     histo4 = h_Scope;
0499     // check the seed phi eta dirstribition
0500     for (size_t j = 0; j < seeds_simtrk.size(); j++) {
0501       if (seeds_simtrk[j] != static_cast<int>(i))
0502         continue;
0503       if (seed_mT[j] < pTCutMin || seed_mT[j] > pTCutMax)
0504         continue;
0505 
0506       histo4->Fill4b(seed_gp[j].phi(), seed_gp[j].eta(), getEta(theta_p[i]));
0507 
0508       // look at the dPhi, dEta, dx, dy from sim-segment and reco-segment of a seed
0509 
0510       if (scope) {
0511         cout << " " << j << "st seed w/ " << nSegInSeed[j] << " segs";
0512         cout << " & pt= " << seed_mT[j] << " +/- " << err_qbp[j] * seed_mT[j] * seed_mA[j];
0513         cout << " q/p= " << qbp[j] << " @ " << seed_layer[j];
0514         cout << " dx= " << err_dx[j] << " dy= " << err_dy[j] << " x= " << err_x[j] << " y= " << err_y[j] << endl;
0515       }
0516 
0517       // check the segments of the seeds in the scope
0518       bool flip_debug = false;
0519       if (!debug) {
0520         debug = true;
0521         flip_debug = true;
0522       }
0523 
0524       std::vector<SimSegment> sCSC_v = recsegSelector->Sim_CSCSegments(trackID[i], csimHits, cscGeom);
0525       std::vector<SimSegment> sDT_v = recsegSelector->Sim_DTSegments(trackID[i], dsimHits, dtGeom);
0526       SegOfRecSeed(muonseeds, j, sCSC_v, sDT_v);
0527 
0528       if (flip_debug) {
0529         debug = false;
0530       }
0531       // check the dEta, dPhi, dx, dy resolution
0532       for (size_t k = 0; k < geoID.size(); k++) {
0533         if (geoID[k].subdetId() == MuonSubdetId::DT) {
0534           DTChamberId MB_Id = DTChamberId(geoID[k]);
0535           histo4->Fill4c(-1 * MB_Id.station(), d_h[k], d_f[k], d_x[k], d_y[k]);
0536         }
0537         if (geoID[k].subdetId() == MuonSubdetId::CSC) {
0538           CSCDetId ME_Id = CSCDetId(geoID[k]);
0539           histo4->Fill4c(ME_Id.station(), d_h[k], d_f[k], d_x[k], d_y[k]);
0540         }
0541       }
0542     }
0543     // check the sta information
0544     for (size_t j = 0; j < sta_simtrk.size(); j++) {
0545       if (sta_simtrk[j] != static_cast<int>(i))
0546         continue;
0547       if (sta_mT[j] < pTCutMin || sta_mT[j] > pTCutMax)
0548         continue;
0549 
0550       histo4->Fill4a(sta_phiV[j], getEta(sta_thetaV[j]), getEta(theta_p[i]));
0551 
0552       // look at the sta pt vs. #_of_hits and chi2
0553       double ndf = static_cast<double>(2 * sta_nHits[j] - 4);
0554       double chi2_ndf = sta_chi2[j] / ndf;
0555       histo4->Fill4d(sta_mT[j], sta_nHits[j], chi2_ndf);
0556 
0557       if (scope) {
0558         cout << " sta_pt= " << sta_mT[j] << " q/p= " << sta_qbp[j] << " w/" << sta_nHits[j] << endl;
0559       }
0560       //cout <<"************************************************"<<endl;
0561       //cout <<"  "<<endl;
0562     }
0563   }
0564 
0565   // D. fail sta cases
0566   for (size_t i = 0; i < nuSEED.size(); i++) {
0567     if (nuSEED[i] == 0)
0568       continue;
0569     if (nuSTA[i] != 0)
0570       continue;
0571 
0572     histo3 = h_NoSta;
0573     double sim_eta = getEta(theta_p[i]);
0574 
0575     for (size_t j = 0; j < seeds_simtrk.size(); j++) {
0576       if (seeds_simtrk[j] != static_cast<int>(i))
0577         continue;
0578       seeds_simtrk[j] = static_cast<int>(i);
0579       double pull_qbpt = (qbpt[j] - (theQ[i] / pt_trk[i])) / err_qbpt[j];
0580       double pt_err = err_qbpt[j] * seed_mT[j] * seed_mT[j];
0581 
0582       histo3->Fill3a(sim_eta, seed_gp[j].eta(), seed_mT[j], pt_err, pull_qbpt);
0583 
0584       CSCsegment_stat(cscSegments, cscGeom, seed_gp[j].theta(), seed_gp[j].phi());
0585       DTsegment_stat(dt4DSegments, dtGeom, seed_gp[j].theta(), seed_gp[j].phi());
0586 
0587       int allseg1 = cscseg_stat1[5] + dtseg_stat1[5];  // # of stations which have good segments
0588       int allseg = cscseg_stat[0] + dtseg_stat[0];     // # of segments
0589       histo3->Fill3b(sim_eta, allseg1, allseg);
0590 
0591       // compare segments of seed with seed itself
0592       SegOfRecSeed(muonseeds, j);
0593       for (unsigned k = 0; k < d_h.size(); k++) {
0594         histo3->Fill3d(sim_eta, d_h[k], d_f[k]);
0595       }
0596 
0597       int types = RecSegReader(cscSegments, dt4DSegments, cscGeom, dtGeom, seed_gp[j].eta(), seed_gp[j].phi());
0598       if (debug)
0599         cout << " seed type for fail STA= " << types << endl;
0600       for (unsigned k = 0; k < phi_resid.size(); k++) {
0601         histo3->Fill3c(phi_resid[k], eta_resid[k]);
0602       }
0603     }
0604 
0605     int nu_seed_trk = 0;
0606     for (size_t j = 0; j < seeds_simtrk.size(); j++) {
0607       if (seeds_simtrk[j] == static_cast<int>(i))
0608         nu_seed_trk++;
0609     }
0610     histo3->Fill3f(sim_eta, nu_seed_trk);
0611   }
0612 }
0613 
0614 // ********************************************
0615 // ***********  Utility functions  ************
0616 // ********************************************
0617 
0618 // number of csc segments in one chamber for each station
0619 // cscseg_stat[0] = total segments in all stations
0620 // cscseg_stat[5] = the number of stations which have segments
0621 void MuonSeedValidator::CSCsegment_stat(Handle<CSCSegmentCollection> cscSeg,
0622                                         ESHandle<CSCGeometry> cscGeom,
0623                                         double trkTheta,
0624                                         double trkPhi) {
0625   for (int i = 0; i < 6; i++) {
0626     cscseg_stat[i] = 0;
0627     cscseg_stat1[i] = 0;
0628   }
0629   for (CSCSegmentCollection::const_iterator seg_It = cscSeg->begin(); seg_It != cscSeg->end(); seg_It++) {
0630     CSCDetId DetId = (CSCDetId)(*seg_It).cscDetId();
0631     const CSCChamber* cscchamber = cscGeom->chamber(DetId);
0632     GlobalPoint gp = cscchamber->toGlobal((*seg_It).localPosition());
0633     GlobalVector gv = cscchamber->toGlobal((*seg_It).localDirection());
0634     if ((fabs(gp.theta() - trkTheta) > dtMax) || (fabs(gv.phi() - trkPhi) > dfMax))
0635       continue;
0636 
0637     double dof = static_cast<double>((*seg_It).degreesOfFreedom());
0638     double chi2_dof = (*seg_It).chi2() / dof;
0639 
0640     cscseg_stat[DetId.station()] += 1;
0641     if ((*seg_It).nRecHits() > 4 && chi2_dof < 200.0) {
0642       cscseg_stat1[DetId.station()] += 1;
0643     }
0644   }
0645   cscseg_stat[0] = cscseg_stat[1] + cscseg_stat[2] + cscseg_stat[3] + cscseg_stat[4];
0646   cscseg_stat1[0] = cscseg_stat1[1] + cscseg_stat1[2] + cscseg_stat1[3] + cscseg_stat1[4];
0647   for (int i = 1; i < 5; i++) {
0648     if (cscseg_stat[i] != 0) {
0649       cscseg_stat[5]++;
0650     }
0651     if (cscseg_stat1[i] != 0) {
0652       cscseg_stat1[5]++;
0653     }
0654   }
0655 }
0656 // number of dt segments in one chamber for each station
0657 void MuonSeedValidator::DTsegment_stat(Handle<DTRecSegment4DCollection> dtSeg,
0658                                        ESHandle<DTGeometry> dtGeom,
0659                                        double trkTheta,
0660                                        double trkPhi) {
0661   for (int i = 0; i < 6; i++) {
0662     dtseg_stat[i] = 0;
0663     dtseg_stat1[i] = 0;
0664   }
0665   for (DTRecSegment4DCollection::const_iterator seg_It = dtSeg->begin(); seg_It != dtSeg->end(); seg_It++) {
0666     //if ( !(*seg_It).hasPhi() || !(*seg_It).hasZed()  ) continue;
0667     if (!(*seg_It).hasPhi())
0668       continue;
0669 
0670     DTChamberId DetId = (*seg_It).chamberId();
0671     const DTChamber* dtchamber = dtGeom->chamber(DetId);
0672     GlobalPoint gp = dtchamber->toGlobal((*seg_It).localPosition());
0673     GlobalVector gv = dtchamber->toGlobal((*seg_It).localDirection());
0674 
0675     if ((fabs(gp.theta() - trkTheta) > dtMax) || (fabs(gv.phi() - trkPhi) > dfMax))
0676       continue;
0677 
0678     dtseg_stat[DetId.station()] += 1;
0679 
0680     int n_phiHits = ((*seg_It).phiSegment())->specificRecHits().size();
0681     //if ( DetId.station() < 4 && (*seg_It).hasZed() && (n_phiHits > 4) ) {
0682     if (DetId.station() < 4 && (*seg_It).dimension() == 4) {
0683       dtseg_stat1[DetId.station()] += 1;
0684     }
0685     if (DetId.station() == 4 && (n_phiHits > 4)) {
0686       dtseg_stat1[DetId.station()] += 1;
0687     }
0688   }
0689   dtseg_stat[0] = dtseg_stat[1] + dtseg_stat[2] + dtseg_stat[3] + dtseg_stat[4];
0690   dtseg_stat1[0] = dtseg_stat1[1] + dtseg_stat1[2] + dtseg_stat1[3] + dtseg_stat1[4];
0691   for (int i = 1; i < 5; i++) {
0692     if (dtseg_stat[i] != 0) {
0693       dtseg_stat[5]++;
0694     }
0695     if (dtseg_stat1[i] != 0) {
0696       dtseg_stat1[5]++;
0697     }
0698   }
0699 }
0700 
0701 // number of sim segments in one chamber for each station
0702 void MuonSeedValidator::Simsegment_stat(std::vector<SimSegment> sCSC, std::vector<SimSegment> sDT) {
0703   for (int i = 0; i < 6; i++) {
0704     simcscseg[i] = 0;
0705     simdtseg[i] = 0;
0706   }
0707 
0708   double ns1 = 0.0;
0709   double eta_sim1 = 0;
0710   for (std::vector<SimSegment>::const_iterator it = sCSC.begin(); it != sCSC.end(); it++) {
0711     int st = ((*it).csc_DetId).station();
0712     eta_sim1 += ((*it).sGlobalOrg).eta();
0713     simcscseg[st]++;
0714     ns1++;
0715   }
0716   simcscseg[0] = simcscseg[1] + simcscseg[2] + simcscseg[3] + simcscseg[4];
0717   for (int i = 1; i < 5; i++) {
0718     if (simcscseg[i] != 0)
0719       simcscseg[5]++;
0720   }
0721   eta_sim1 = eta_sim1 / ns1;
0722 
0723   double ns2 = 0.0;
0724   double eta_sim2 = 0;
0725   for (std::vector<SimSegment>::const_iterator it = sDT.begin(); it != sDT.end(); it++) {
0726     int st = ((*it).dt_DetId).station();
0727     eta_sim2 += ((*it).sGlobalOrg).eta();
0728     simdtseg[st]++;
0729     ns2++;
0730   }
0731   simdtseg[0] = simdtseg[1] + simdtseg[2] + simdtseg[3] + simdtseg[4];
0732   for (int i = 1; i < 5; i++) {
0733     if (simdtseg[i] != 0)
0734       simdtseg[5]++;
0735   }
0736   eta_sim2 = eta_sim2 / ns2;
0737 
0738   simseg_sum = simcscseg[5] + simdtseg[5];
0739   simseg_eta = -9.0;
0740   if ((simcscseg[0] == 0) && (simdtseg[0] != 0)) {
0741     simseg_eta = eta_sim2;
0742   } else if ((simdtseg[0] == 0) && (simcscseg[0] != 0)) {
0743     simseg_eta = eta_sim1;
0744   } else {
0745     simseg_eta = (eta_sim1 + eta_sim2) / 2.0;
0746   }
0747 }
0748 
0749 void MuonSeedValidator::CSCRecHit_Stat(Handle<CSCRecHit2DCollection> cscrechit,
0750                                        ESHandle<CSCGeometry> cscGeom,
0751                                        double trkEta,
0752                                        double trkPhi) {
0753   for (int i = 0; i < 6; i++) {
0754     cscrh_sum[i] = 0;
0755   }
0756   for (CSCRecHit2DCollection::const_iterator r_it = cscrechit->begin(); r_it != cscrechit->end(); r_it++) {
0757     CSCDetId det_id = (CSCDetId)(*r_it).cscDetId();
0758     const CSCChamber* cscchamber = cscGeom->chamber(det_id);
0759     GlobalPoint gp = cscchamber->toGlobal((*r_it).localPosition());
0760     if ((fabs(gp.eta() - trkEta) > dtMax) || (fabs(gp.phi() - trkPhi) > dfMax))
0761       continue;
0762 
0763     cscrh_sum[det_id.station()]++;
0764   }
0765   cscrh_sum[0] = cscrh_sum[1] + cscrh_sum[2] + cscrh_sum[3] + cscrh_sum[4];
0766   for (int i = 1; i < 5; i++) {
0767     if (cscrh_sum[i] != 0) {
0768       cscrh_sum[5]++;
0769     }
0770   }
0771 }
0772 
0773 void MuonSeedValidator::DTRecHit_Stat(Handle<DTRecHitCollection> dtrechit,
0774                                       ESHandle<DTGeometry> dtGeom,
0775                                       double trkEta,
0776                                       double trkPhi) {
0777   //double phi[4]={999.0};
0778   for (int i = 0; i < 6; i++) {
0779     dtrh_sum[i] = 0;
0780   }
0781 
0782   double eta = -9.0;
0783   double nn = 0.0;
0784   for (DTRecHitCollection::const_iterator r_it = dtrechit->begin(); r_it != dtrechit->end(); r_it++) {
0785     DTWireId det_id = (*r_it).wireId();
0786     const DTChamber* dtchamber = dtGeom->chamber(det_id);
0787     LocalPoint lrh = (*r_it).localPosition();
0788     GlobalPoint grh = dtchamber->toGlobal(lrh);
0789     if ((fabs(grh.eta() - trkEta) > dtMax) || (fabs(grh.phi() - trkPhi) > dfMax))
0790       continue;
0791 
0792     dtrh_sum[det_id.station()]++;
0793     eta += grh.eta();
0794     nn += 1.0;
0795   }
0796   eta = eta / nn;
0797 
0798   dtrh_sum[0] = dtrh_sum[1] + dtrh_sum[2] + dtrh_sum[3] + dtrh_sum[4];
0799   for (int i = 1; i < 5; i++) {
0800     if (dtrh_sum[i] != 0) {
0801       dtrh_sum[5]++;
0802     }
0803   }
0804 }
0805 
0806 int MuonSeedValidator::ChargeAssignment(GlobalVector Va, GlobalVector Vb) {
0807   int charge = 0;
0808   float axb = (Va.x() * Vb.y()) - (Vb.x() * Va.y());
0809   if (axb != 0.0) {
0810     charge = ((axb > 0.0) ? 1 : -1);
0811   }
0812   return charge;
0813 }
0814 
0815 void MuonSeedValidator::RecSeedReader(Handle<TrajectorySeedCollection> rec_seeds) {
0816   nu_seed = 0;
0817   seed_gp.clear();
0818   seed_gm.clear();
0819   seed_lp.clear();
0820   seed_lv.clear();
0821   qbp.clear();
0822   qbpt.clear();
0823   err_qbp.clear();
0824   err_qbpt.clear();
0825   err_dx.clear();
0826   err_dy.clear();
0827   err_x.clear();
0828   err_y.clear();
0829   seed_mT.clear();
0830   seed_mA.clear();
0831   seed_layer.clear();
0832   nSegInSeed.clear();
0833 
0834   TrajectorySeedCollection::const_iterator seed_it;
0835   for (seed_it = rec_seeds->begin(); seed_it != rec_seeds->end(); seed_it++) {
0836     PTrajectoryStateOnDet pTSOD = (*seed_it).startingState();
0837 
0838     nSegInSeed.push_back((*seed_it).nHits());
0839     // Get the tsos of the seed
0840     DetId seedDetId(pTSOD.detId());
0841     const GeomDet* gdet = theService->trackingGeometry()->idToDet(seedDetId);
0842     TrajectoryStateOnSurface seedTSOS =
0843         trajectoryStateTransform::transientState(pTSOD, &(gdet->surface()), &*theService->magneticField());
0844     // seed global position and momentum(direction)
0845     seed_gp.push_back(seedTSOS.globalPosition());
0846     seed_gm.push_back(seedTSOS.globalMomentum());
0847     seed_lp.push_back(seedTSOS.localPosition());
0848     seed_lv.push_back(seedTSOS.localDirection());
0849 
0850     LocalTrajectoryParameters seed_para = pTSOD.parameters();
0851 
0852     // the error_vector v[15]-> [0] , [2] , [5] , [9] , [14]
0853     //                          q/p   dx    dy     x      y
0854     err_qbp.push_back(sqrt(pTSOD.error(0)));
0855     err_dx.push_back(sqrt(pTSOD.error(2)));
0856     err_dy.push_back(sqrt(pTSOD.error(5)));
0857     err_x.push_back(sqrt(pTSOD.error(9)));
0858     err_y.push_back(sqrt(pTSOD.error(14)));
0859 
0860     // seed layer
0861     DetId pdid(pTSOD.detId());
0862     if (pdid.subdetId() == MuonSubdetId::DT) {
0863       DTChamberId MB_Id = DTChamberId(pTSOD.detId());
0864       seed_layer.push_back(MB_Id.station());
0865     }
0866     if (pdid.subdetId() == MuonSubdetId::CSC) {
0867       CSCDetId ME_Id = CSCDetId(pTSOD.detId());
0868       seed_layer.push_back(ME_Id.station());
0869     }
0870     double seed_mx = seed_gm[nu_seed].x();
0871     double seed_my = seed_gm[nu_seed].y();
0872     double seed_mz = seed_gm[nu_seed].z();
0873     double seed_q = 1.0 * seed_para.charge();
0874     double seed_sin = sqrt((seed_mx * seed_mx) + (seed_my * seed_my)) /
0875                       sqrt((seed_mx * seed_mx) + (seed_my * seed_my) + (seed_mz * seed_mz));
0876     // seed pt, pa, and q/pt or q/pa
0877     seed_mT.push_back(sqrt((seed_mx * seed_mx) + (seed_my * seed_my)));
0878     seed_mA.push_back(sqrt((seed_mx * seed_mx) + (seed_my * seed_my) + (seed_mz * seed_mz)));
0879     qbpt.push_back(seed_q / sqrt((seed_mx * seed_mx) + (seed_my * seed_my)));
0880     qbp.push_back(seed_para.signedInverseMomentum());
0881     err_qbpt.push_back(err_qbp[nu_seed] / seed_sin);
0882     nu_seed++;
0883     if (debug) {
0884       cout << " seed pt: " << sqrt((seed_mx * seed_mx) + (seed_my * seed_my)) << endl;
0885     }
0886   }
0887 }
0888 
0889 // read the segments associated with the seed and compare with seed
0890 void MuonSeedValidator::SegOfRecSeed(Handle<TrajectorySeedCollection> rec_seeds, int seed_idx) {
0891   int idx = 0;
0892   d_h.clear();
0893   d_f.clear();
0894   d_x.clear();
0895   d_y.clear();
0896   geoID.clear();
0897 
0898   TrajectorySeedCollection::const_iterator seed_it;
0899   for (seed_it = rec_seeds->begin(); seed_it != rec_seeds->end(); seed_it++) {
0900     idx++;
0901     if (seed_idx != (idx - 1))
0902       continue;
0903 
0904     PTrajectoryStateOnDet pTSOD = (*seed_it).startingState();
0905     DetId seedDetId(pTSOD.detId());
0906     const GeomDet* seedDet = theService->trackingGeometry()->idToDet(seedDetId);
0907     TrajectoryStateOnSurface seedTSOS =
0908         trajectoryStateTransform::transientState(pTSOD, &(seedDet->surface()), &*theService->magneticField());
0909     GlobalPoint seedgp = seedTSOS.globalPosition();
0910 
0911     for (auto const& recHit : seed_it->recHits()) {
0912       const GeomDet* gdet = theService->trackingGeometry()->idToDet(recHit.geographicalId());
0913       LocalPoint lp = recHit.localPosition();
0914       GlobalPoint gp = gdet->toGlobal(lp);
0915       LocalPoint slp = gdet->toLocal(gp);
0916       DetId pdid = recHit.geographicalId();
0917 
0918       geoID.push_back(pdid);
0919       d_h.push_back(gp.eta() - seedgp.eta());
0920       d_f.push_back(gp.phi() - seedgp.phi());
0921       d_x.push_back(lp.x() - slp.x());
0922       d_y.push_back(lp.y() - slp.y());
0923     }
0924   }
0925 }
0926 
0927 // read the segments associated with the seed and compare with sim-segment
0928 void MuonSeedValidator::SegOfRecSeed(Handle<TrajectorySeedCollection> rec_seeds,
0929                                      int seed_idx,
0930                                      std::vector<SimSegment> sCSC,
0931                                      std::vector<SimSegment> sDT) {
0932   int idx = 0;
0933   d_h.clear();
0934   d_f.clear();
0935   d_x.clear();
0936   d_y.clear();
0937   geoID.clear();
0938 
0939   TrajectorySeedCollection::const_iterator seed_it;
0940   for (seed_it = rec_seeds->begin(); seed_it != rec_seeds->end(); seed_it++) {
0941     idx++;
0942     if (seed_idx != (idx - 1))
0943       continue;
0944     if (debug && scope) {
0945       cout << " " << endl;
0946     }
0947 
0948     for (auto const& recHit : seed_it->recHits()) {
0949       const GeomDet* gdet = theService->trackingGeometry()->idToDet(recHit.geographicalId());
0950       GlobalPoint gp = gdet->toGlobal(recHit.localPosition());
0951       LocalPoint lp = recHit.localPosition();
0952       DetId pdid = recHit.geographicalId();
0953 
0954       // for parameters [1]:dx/dz, [2]:dy/dz, [3]:x, [4]:y
0955       double dxz = recHit.parameters()[0];
0956       double dyz = recHit.parameters()[1];
0957       double dz = 1.0 / sqrt((dxz * dxz) + (dyz * dyz) + 1.0);
0958       if (pdid.subdetId() == MuonSubdetId::DT) {
0959         dz = -1.0 * dz;
0960       }
0961       double dx = dxz * dz;
0962       double dy = dyz * dz;
0963       LocalVector lv = LocalVector(dx, dy, dz);
0964       GlobalVector gv = gdet->toGlobal(lv);
0965 
0966       if (pdid.subdetId() == MuonSubdetId::CSC) {
0967         double directionSign = gp.z() * gv.z();
0968         lv = (directionSign * lv).unit();
0969         gv = gdet->toGlobal(lv);
0970       }
0971 
0972       if (debug && scope) {
0973         cout << "============= segs from seed ============== " << endl;
0974       }
0975 
0976       if (pdid.subdetId() == MuonSubdetId::DT) {
0977         DTChamberId MB_Id = DTChamberId(pdid);
0978 
0979         if (debug && scope) {
0980           cout << "DId: " << MB_Id << endl;
0981         }
0982 
0983         // find the sim-reco match case and store the difference
0984         if (sDT.size() == 0) {
0985           geoID.push_back(pdid);
0986           d_h.push_back(999.0);
0987           d_f.push_back(999.0);
0988           d_x.push_back(999.0);
0989           d_y.push_back(999.0);
0990         } else {
0991           bool match = false;
0992           for (std::vector<SimSegment>::const_iterator it = sDT.begin(); it != sDT.end(); it++) {
0993             if ((*it).dt_DetId == MB_Id) {
0994               geoID.push_back(pdid);
0995               d_h.push_back(gv.eta() - (*it).sGlobalVec.eta());
0996               d_f.push_back(gv.phi() - (*it).sGlobalVec.phi());
0997               d_x.push_back(lp.x() - (*it).sLocalOrg.x());
0998               d_y.push_back(lp.y() - (*it).sLocalOrg.y());
0999               match = true;
1000             }
1001           }
1002           if (!match) {
1003             geoID.push_back(pdid);
1004             d_h.push_back(999.0);
1005             d_f.push_back(999.0);
1006             d_x.push_back(999.0);
1007             d_y.push_back(999.0);
1008           }
1009         }
1010       }
1011       if (pdid.subdetId() == MuonSubdetId::CSC) {
1012         CSCDetId ME_Id = CSCDetId(pdid);
1013         if (debug && scope) {
1014           cout << "DId: " << ME_Id << endl;
1015         }
1016 
1017         /*/ flip the z-sign for station 1 and 2 because of the chamber orientation
1018                 if ( (ME_Id.station()==1)||(ME_Id.station()==2) ) {
1019                    gv = GlobalVector( gv.x(), gv.y(), (-1.*gv.z()) );
1020                 }*/
1021 
1022         // find the sim-reco match case and store the difference
1023         if (sCSC.size() == 0) {
1024           geoID.push_back(pdid);
1025           d_h.push_back(999.0);
1026           d_f.push_back(999.0);
1027           d_x.push_back(999.0);
1028           d_y.push_back(999.0);
1029         } else {
1030           bool match = false;
1031           for (std::vector<SimSegment>::const_iterator it = sCSC.begin(); it != sCSC.end(); it++) {
1032             if ((*it).csc_DetId == ME_Id) {
1033               geoID.push_back(pdid);
1034               d_h.push_back(gv.eta() - (*it).sGlobalVec.eta());
1035               d_f.push_back(gv.phi() - (*it).sGlobalVec.phi());
1036               d_x.push_back(lp.x() - (*it).sLocalOrg.x());
1037               d_y.push_back(lp.y() - (*it).sLocalOrg.y());
1038               match = true;
1039             }
1040           }
1041           if (!match) {
1042             geoID.push_back(pdid);
1043             d_h.push_back(999.0);
1044             d_f.push_back(999.0);
1045             d_x.push_back(999.0);
1046             d_y.push_back(999.0);
1047           }
1048         }
1049       }
1050 
1051       if (debug && scope) {
1052         cout << "  h0= " << gp.eta() << "  f0= " << gp.phi() << "   gp= " << gp << endl;
1053         cout << "  h1= " << gv.eta() << "  f1= " << gv.phi() << "   gv= " << gv << endl;
1054       }
1055     }
1056   }
1057 }
1058 
1059 void MuonSeedValidator::StaTrackReader(Handle<reco::TrackCollection> sta_trk, int sta_glb) {
1060   // look at the inner most momentum and position
1061   nu_sta = 0;
1062 
1063   sta_phiP.clear();
1064   sta_thetaP.clear();
1065 
1066   sta_mT.clear();
1067   sta_mA.clear();
1068   sta_thetaV.clear();
1069   sta_phiV.clear();
1070   sta_qbp.clear();
1071   sta_qbpt.clear();
1072   sta_chi2.clear();
1073   sta_nHits.clear();
1074 
1075   TrackCollection::const_iterator iTrk;
1076   for (iTrk = sta_trk->begin(); iTrk != sta_trk->end(); iTrk++) {
1077     nu_sta++;
1078 
1079     // get the inner poistion( eta & phi )
1080     math::XYZPoint staPos = (*iTrk).innerPosition();
1081     double posMag = sqrt(staPos.x() * staPos.x() + staPos.y() * staPos.y() + staPos.z() * staPos.z());
1082     //math::XYZVector staMom = (*iTrk).innerMomentum();
1083     //double innerMt = sqrt( (staMom.x()*staMom.x()) + (staMom.y()*staMom.y()) );
1084 
1085     sta_phiP.push_back(atan2(staPos.y(), staPos.x()));
1086     sta_thetaP.push_back(acos(staPos.z() / posMag));
1087 
1088     //cout<<" momentum= "<<(*iTrk).momentum()<<"  pt= "<<(*iTrk).pt()<<endl;
1089     //cout<<" innerMom= "<<(*iTrk).innerMomentum()<<" iMt= "<<innerMt<<endl;
1090 
1091     sta_mA.push_back((*iTrk).p());
1092     sta_mT.push_back((*iTrk).pt());
1093     //sta_mT.push_back( innerMt );
1094     sta_thetaV.push_back((*iTrk).theta());
1095     sta_phiV.push_back((*iTrk).phi());
1096     sta_qbp.push_back((*iTrk).qoverp());
1097     sta_qbpt.push_back(((*iTrk).qoverp() / (*iTrk).pt()) * (*iTrk).p());
1098     sta_chi2.push_back((*iTrk).chi2());
1099     sta_nHits.push_back((*iTrk).recHitsSize());
1100 
1101     if (debug) {
1102       if (sta_glb == 0)
1103         cout << "sta";
1104       if (sta_glb == 1)
1105         cout << "glb";
1106 
1107       cout << " track  pt: " << (*iTrk).pt() << endl;
1108     }
1109   }
1110 }
1111 
1112 void MuonSeedValidator::SimInfo(Handle<edm::SimTrackContainer> simTracks,
1113                                 Handle<edm::PSimHitContainer> dsimHits,
1114                                 Handle<edm::PSimHitContainer> csimHits,
1115                                 ESHandle<DTGeometry> dtGeom,
1116                                 ESHandle<CSCGeometry> cscGeom) {
1117   // theta and phi at inner-most layer of Muon System
1118   theta_p.clear();
1119   theta_v.clear();
1120   phi_p.clear();
1121   phi_v.clear();
1122   // basic sim track infomation
1123   eta_trk.clear();
1124   theta_trk.clear();
1125   phi_trk.clear();
1126   theQ.clear();
1127   pt_trk.clear();
1128   ptlayer.clear();
1129   palayer.clear();
1130   trackID.clear();
1131 
1132   for (SimTrackContainer::const_iterator simTk_It = simTracks->begin(); simTk_It != simTracks->end(); simTk_It++) {
1133     //if (abs((*simTk_It).type())!=13 || (*simTk_It).vertIndex() != 0 ) continue;
1134     bool rechitSize = (dsimHits->size() < 8 && csimHits->size() < 4) ? true : false;
1135     if (abs((*simTk_It).type()) != 13 || rechitSize || (*simTk_It).vertIndex() != 0)
1136       continue;
1137 
1138     trackID.push_back(static_cast<int>((*simTk_It).trackId()));
1139     if ((*simTk_It).type() == 13) {
1140       theQ.push_back(-1.0);
1141     } else {
1142       theQ.push_back(1.0);
1143     }
1144 
1145     std::vector<double> pt1(5, 0.0);
1146     std::vector<double> pa1(5, 0.0);
1147 
1148     double px = ((*simTk_It).momentum()).x();
1149     double py = ((*simTk_It).momentum()).y();
1150     double pz = ((*simTk_It).momentum()).z();
1151     pa1[0] = sqrt(px * px + py * py + pz * pz);
1152     pt1[0] = sqrt(px * px + py * py);
1153 
1154     eta_trk.push_back(getEta(px, py, pz));
1155     theta_trk.push_back(acos(pz / pa1[0]));
1156     phi_trk.push_back(atan2(py, px));
1157     pt_trk.push_back(pt1[0]);
1158 
1159     double enu2 = 0.0;
1160     for (PSimHitContainer::const_iterator ds_It = dsimHits->begin(); ds_It != dsimHits->end(); ds_It++) {
1161       Local3DPoint lp = (*ds_It).localPosition();
1162 
1163       DTLayerId D_Id = DTLayerId((*ds_It).detUnitId());
1164       const DTLayer* dtlayer = dtGeom->layer(D_Id);
1165       GlobalVector m2 = dtlayer->toGlobal((*ds_It).momentumAtEntry());
1166       GlobalPoint gp = dtlayer->toGlobal(lp);
1167 
1168       if ((abs((*ds_It).particleType()) == 13) && ((*ds_It).trackId() == (*simTk_It).trackId())) {
1169         pt1[D_Id.station()] = sqrt((m2.x() * m2.x()) + (m2.y() * m2.y()));
1170         pa1[D_Id.station()] = sqrt((m2.x() * m2.x()) + (m2.y() * m2.y()) + (m2.z() * m2.z()));
1171 
1172         if (enu2 == 0) {
1173           theta_p.push_back(gp.theta());
1174           theta_v.push_back(m2.theta());
1175           phi_p.push_back(gp.phi());
1176           phi_v.push_back(m2.phi());
1177         }
1178         enu2 += 1.0;
1179       }
1180     }
1181 
1182     double enu1 = 0.0;
1183     for (PSimHitContainer::const_iterator cs_It = csimHits->begin(); cs_It != csimHits->end(); cs_It++) {
1184       CSCDetId C_Id = CSCDetId((*cs_It).detUnitId());
1185       const CSCChamber* cscchamber = cscGeom->chamber(C_Id);
1186       GlobalVector m1 = cscchamber->toGlobal((*cs_It).momentumAtEntry());
1187       Local3DPoint lp = (*cs_It).localPosition();
1188       GlobalPoint gp = cscchamber->toGlobal(lp);
1189 
1190       if ((abs((*cs_It).particleType()) == 13) && ((*cs_It).trackId() == (*simTk_It).trackId())) {
1191         if (enu2 == 0.0) {
1192           pt1[C_Id.station()] = sqrt((m1.x() * m1.x()) + (m1.y() * m1.y()));
1193           pa1[C_Id.station()] = sqrt((m1.x() * m1.x()) + (m1.y() * m1.y()) + (m1.z() * m1.z()));
1194 
1195           if (enu1 == 0) {
1196             theta_p.push_back(gp.theta());
1197             theta_v.push_back(m1.theta());
1198             phi_p.push_back(gp.phi());
1199             phi_v.push_back(m1.phi());
1200           }
1201         }
1202         enu1 += 1.0;
1203       }
1204     }
1205 
1206     ptlayer.push_back(pt1);
1207     palayer.push_back(pa1);
1208     //cout<<" simTrk momentum= "<<(*simTk_It).momentum()<<" pa= "<<pa1[0]<<" pt= "<<pt1[0]<<endl;
1209     //cout<<" simhit momentum= "<< pa1[1] <<" pt= "<<pt1[1]<<endl;
1210   }
1211 }
1212 
1213 // Look up what segments we have in a event
1214 int MuonSeedValidator::RecSegReader(Handle<CSCSegmentCollection> cscSeg,
1215                                     Handle<DTRecSegment4DCollection> dtSeg,
1216                                     ESHandle<CSCGeometry> cscGeom,
1217                                     ESHandle<DTGeometry> dtGeom,
1218                                     double trkTheta,
1219                                     double trkPhi) {
1220   // Calculate the ave. eta & phi
1221   ave_phi = 0.0;
1222   ave_eta = 0.0;
1223   phi_resid.clear();
1224   eta_resid.clear();
1225 
1226   double n = 0.0;
1227   double m = 0.0;
1228   for (CSCSegmentCollection::const_iterator it = cscSeg->begin(); it != cscSeg->end(); it++) {
1229     if ((*it).nRecHits() < 4)
1230       continue;
1231     CSCDetId DetId = (CSCDetId)(*it).cscDetId();
1232     const CSCChamber* cscchamber = cscGeom->chamber(DetId);
1233     GlobalPoint gp = cscchamber->toGlobal((*it).localPosition());
1234     if ((fabs(gp.theta() - trkTheta) > 0.5) || (fabs(gp.phi() - trkPhi) > 0.5))
1235       continue;
1236     GlobalVector gv = cscchamber->toGlobal((*it).localDirection());
1237 
1238     ave_phi += gp.phi();
1239     ave_eta += gp.eta();
1240     dx_error.push_back((*it).localDirectionError().xx());
1241     dy_error.push_back((*it).localDirectionError().yy());
1242     x_error.push_back((*it).localPositionError().xx());
1243     y_error.push_back((*it).localPositionError().yy());
1244     n++;
1245     if (debug) {
1246       cout << "~~~~~~~~~~~~~~~~  reco segs  ~~~~~~~~~~~~~~~~  " << endl;
1247       cout << "DId: " << DetId << endl;
1248       cout << "  h0= " << gp.eta() << "  f0= " << gp.phi() << "   gp= " << gp << endl;
1249       cout << "  h1= " << gv.eta() << "  f1= " << gv.phi() << "   gv= " << gv << endl;
1250     }
1251   }
1252   for (DTRecSegment4DCollection::const_iterator it = dtSeg->begin(); it != dtSeg->end(); it++) {
1253     if (!(*it).hasPhi() || !(*it).hasZed())
1254       continue;
1255     DTChamberId DetId = (*it).chamberId();
1256     const DTChamber* dtchamber = dtGeom->chamber(DetId);
1257     GlobalPoint gp = dtchamber->toGlobal((*it).localPosition());
1258     if ((fabs(gp.eta() - trkTheta) > 0.5) || (fabs(gp.phi() - trkPhi) > 0.5))
1259       continue;
1260     GlobalVector gv = dtchamber->toGlobal((*it).localDirection());
1261 
1262     ave_phi += gp.phi();
1263     ave_eta += gp.eta();
1264     m++;
1265     if (debug) {
1266       cout << "~~~~~~~~~~~~~~~~  reco segs  ~~~~~~~~~~~~~~~~  " << endl;
1267       cout << "DId: " << DetId << endl;
1268       cout << "  h0= " << gp.eta() << "  f0= " << gp.phi() << "   gp= " << gp << endl;
1269       cout << "  h1= " << gv.eta() << "  f1= " << gv.phi() << "   gv= " << gv << endl;
1270     }
1271   }
1272   ave_phi = ave_phi / (n + m);
1273   ave_eta = ave_eta / (n + m);
1274 
1275   // Calculate the residual of phi and eta
1276   for (CSCSegmentCollection::const_iterator it = cscSeg->begin(); it != cscSeg->end(); it++) {
1277     if ((*it).nRecHits() < 4)
1278       continue;
1279     CSCDetId DetId = (CSCDetId)(*it).cscDetId();
1280     const CSCChamber* cscchamber = cscGeom->chamber(DetId);
1281     GlobalPoint gp = cscchamber->toGlobal((*it).localPosition());
1282     if ((fabs(gp.theta() - trkTheta) > 0.5) || (fabs(gp.phi() - trkPhi) > 0.5))
1283       continue;
1284     phi_resid.push_back(gp.phi() - ave_phi);
1285     eta_resid.push_back(gp.eta() - ave_eta);
1286   }
1287   for (DTRecSegment4DCollection::const_iterator it = dtSeg->begin(); it != dtSeg->end(); it++) {
1288     if (!(*it).hasPhi() || !(*it).hasZed())
1289       continue;
1290     DTChamberId DetId = (*it).chamberId();
1291     const DTChamber* dtchamber = dtGeom->chamber(DetId);
1292     GlobalPoint gp = dtchamber->toGlobal((*it).localPosition());
1293     if ((fabs(gp.eta() - trkTheta) > 0.5) || (fabs(gp.phi() - trkPhi) > 0.5))
1294       continue;
1295     phi_resid.push_back(gp.phi() - ave_phi);
1296     eta_resid.push_back(gp.eta() - ave_eta);
1297   }
1298 
1299   if (n != 0 && m == 0) {
1300     return 1;  // csc
1301   } else if (n == 0 && m != 0) {
1302     return 3;  // dt
1303   } else if (n != 0 && m != 0) {
1304     return 2;  // overlap
1305   } else {
1306     return 0;  // fail
1307   }
1308 }
1309 
1310 double MuonSeedValidator::getEta(double vx, double vy, double vz) {
1311   double va = sqrt(vx * vx + vy * vy + vz * vz);
1312 
1313   double theta = acos(vz / va);
1314   double eta = (-1.0) * log(tan(theta / 2.0));
1315   return eta;
1316 }
1317 double MuonSeedValidator::getEta(double theta) {
1318   double eta = (-1.0) * log(tan(theta / 2.0));
1319   return eta;
1320 }
1321 
1322 std::vector<int> MuonSeedValidator::IdentifyShowering(Handle<CSCSegmentCollection> cscSeg,
1323                                                       ESHandle<CSCGeometry> cscGeom,
1324                                                       Handle<DTRecSegment4DCollection> dtSeg,
1325                                                       ESHandle<DTGeometry> dtGeom,
1326                                                       double trkTheta,
1327                                                       double trkPhi) {
1328   muCone.clear();
1329 
1330   std::vector<double> csch1;
1331   std::vector<double> csch2;
1332   std::vector<double> csch3;
1333   std::vector<double> csch4;
1334 
1335   std::vector<double> cscf1;
1336   std::vector<double> cscf2;
1337   std::vector<double> cscf3;
1338   std::vector<double> cscf4;
1339 
1340   std::vector<int> showers(4, 0);
1341   int CSC1 = 0;
1342   int CSC2 = 0;
1343   int CSC3 = 0;
1344   int CSC4 = 0;
1345   for (CSCSegmentCollection::const_iterator i1 = cscSeg->begin(); i1 != cscSeg->end(); i1++) {
1346     if (trkTheta > 0.87)
1347       continue;
1348     CSCDetId DId = (CSCDetId)(*i1).cscDetId();
1349     const CSCChamber* cscchamber = cscGeom->chamber(DId);
1350     GlobalPoint gp = cscchamber->toGlobal((*i1).localPosition());
1351     double dh = gp.eta() - getEta(trkTheta);
1352     double df = gp.phi() - trkPhi;
1353     double dR = sqrt((dh * dh) + (df * df));
1354 
1355     if (DId.station() == 1 && dR < 0.5) {
1356       CSC1++;
1357       csch1.push_back(gp.eta());
1358       cscf1.push_back(gp.phi());
1359     }
1360     if (DId.station() == 2 && dR < 0.5) {
1361       CSC2++;
1362       csch2.push_back(gp.eta());
1363       cscf2.push_back(gp.phi());
1364     }
1365     if (DId.station() == 3 && dR < 0.5) {
1366       CSC3++;
1367       csch3.push_back(gp.eta());
1368       cscf3.push_back(gp.phi());
1369     }
1370     if (DId.station() == 4 && dR < 0.5) {
1371       CSC4++;
1372       csch4.push_back(gp.eta());
1373       cscf4.push_back(gp.phi());
1374     }
1375   }
1376 
1377   if (CSC1 > 2) {
1378     showers[0] += CSC1;
1379     showers[3] += 1;
1380     muCone.push_back(getdR(csch1, cscf1));
1381   }
1382   if (CSC2 > 2) {
1383     showers[0] += CSC2;
1384     showers[1] += CSC2;
1385     showers[2] += 1;
1386     showers[3] += 1;
1387     muCone.push_back(getdR(csch1, cscf1));
1388   }
1389   if (CSC3 > 2) {
1390     showers[0] += CSC3;
1391     showers[1] += CSC3;
1392     showers[2] += 1;
1393     showers[3] += 1;
1394     muCone.push_back(getdR(csch1, cscf1));
1395   }
1396   if (CSC4 > 2) {
1397     showers[0] += CSC4;
1398     showers[1] += CSC4;
1399     showers[2] += 1;
1400     showers[3] += 1;
1401     muCone.push_back(getdR(csch1, cscf1));
1402   }
1403 
1404   std::vector<double> dth1;
1405   std::vector<double> dth2;
1406   std::vector<double> dth3;
1407 
1408   std::vector<double> dtf1;
1409   std::vector<double> dtf2;
1410   std::vector<double> dtf3;
1411 
1412   int DT1 = 0;
1413   int DT2 = 0;
1414   int DT3 = 0;
1415   for (DTRecSegment4DCollection::const_iterator i1 = dtSeg->begin(); i1 != dtSeg->end(); i1++) {
1416     if (trkTheta < 0.52)
1417       continue;
1418     DTChamberId DId = (*i1).chamberId();
1419     const DTChamber* dtchamber = dtGeom->chamber(DId);
1420     GlobalPoint gp = dtchamber->toGlobal((*i1).localPosition());
1421 
1422     double dh = gp.eta() - getEta(trkTheta);
1423     double df = gp.phi() - trkPhi;
1424     double dR = sqrt((dh * dh) + (df * df));
1425 
1426     if (DId.station() == 1 && dR < 0.5) {
1427       DT1++;
1428       dth1.push_back(gp.eta());
1429       dtf1.push_back(gp.phi());
1430     }
1431     if (DId.station() == 2 && dR < 0.5) {
1432       DT2++;
1433       dth2.push_back(gp.eta());
1434       dtf2.push_back(gp.phi());
1435     }
1436     if (DId.station() == 3 && dR < 0.5) {
1437       DT3++;
1438       dth3.push_back(gp.eta());
1439       dtf3.push_back(gp.phi());
1440     }
1441   }
1442   if (DT1 > 2) {
1443     showers[0] += DT1;
1444     muCone.push_back(getdR(dth1, dtf1));
1445     showers[3] += 1;
1446   }
1447   if (DT2 > 2) {
1448     showers[0] += DT2;
1449     showers[1] += DT2;
1450     showers[2] += 1;
1451     showers[3] += 1;
1452     muCone.push_back(getdR(dth2, dtf2));
1453   }
1454   if (DT3 > 2) {
1455     showers[0] += DT3;
1456     showers[1] += DT3;
1457     showers[2] += 1;
1458     showers[3] += 1;
1459     muCone.push_back(getdR(dth3, dtf3));
1460   }
1461 
1462   return showers;
1463 }
1464 
1465 double MuonSeedValidator::getdR(std::vector<double> etav, std::vector<double> phiv) {
1466   if (etav.size() < 3)
1467     return -1.0;
1468 
1469   double avh = 0.;
1470   double avf = 0.;
1471   double sz = 0.;
1472   for (size_t i = 0; i < etav.size(); i++) {
1473     avh += etav[i];
1474     avf += phiv[i];
1475     sz += 1.;
1476   }
1477   avh = avh / sz;
1478   avf = avf / sz;
1479 
1480   double maxdh = 0.;
1481   double maxdf = 0.;
1482   for (size_t i = 0; i < etav.size(); i++) {
1483     double dh = fabs(etav[i] - avh);
1484     double df = fabs(phiv[i] - avf);
1485     if (dh > maxdh)
1486       maxdh = dh;
1487     if (df > maxdf)
1488       maxdf = df;
1489     //cout <<" df:"<<df<<"  dh:"<<dh<<endl;
1490   }
1491   //double maxdR = sqrt( (maxdh*maxdh) + (maxdf*maxdf)  ) ;
1492   //return maxdR ;
1493   return maxdf;
1494 }
[ 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!