Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 // Class Header
0002 #include "SegSelector.h"
0003 #include "MuonSeedParametrization.h"
0004 
0005 // Framework
0006 #include "FWCore/Framework/interface/Event.h"
0007 #include "FWCore/Framework/interface/EventSetup.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "FWCore/Utilities/interface/Exception.h"
0011 
0012 //#include "PluginManager/ModuleDef.h"
0013 #include "FWCore/Framework/interface/MakerMacros.h"
0014 
0015 #include "TFile.h"
0016 #include "TVector3.h"
0017 
0018 #include <iostream>
0019 #include <fstream>
0020 #include <map>
0021 #include <utility>
0022 #include <string>
0023 #include <stdio.h>
0024 #include <algorithm>
0025 
0026 DEFINE_FWK_MODULE(MuonSeedParametrization);
0027 using namespace std;
0028 using namespace edm;
0029 
0030 // constructors
0031 MuonSeedParametrization::MuonSeedParametrization(const ParameterSet& pset)
0032     : cscGeomToken(esConsumes()), dtGeomToken(esConsumes()) {
0033   debug = pset.getUntrackedParameter<bool>("debug");
0034   scale = pset.getUntrackedParameter<bool>("scale");
0035   rootFileName = pset.getUntrackedParameter<string>("rootFileName");
0036   recHitLabel = pset.getUntrackedParameter<string>("recHitLabel");
0037   cscSegmentLabel = pset.getUntrackedParameter<string>("cscSegmentLabel");
0038   dtrecHitLabel = pset.getUntrackedParameter<string>("dtrecHitLabel");
0039   dtSegmentLabel = pset.getUntrackedParameter<string>("dtSegmentLabel");
0040   //dt2DSegmentLabel  = pset.getUntrackedParameter<string>("dt2DSegmentLabel");
0041   simHitLabel = pset.getUntrackedParameter<string>("simHitLabel");
0042   simTrackLabel = pset.getUntrackedParameter<string>("simTrackLabel");
0043 
0044   recsegSelector = new SegSelector(pset);
0045   HistoFill = new MuonSeedParaFillHisto();
0046   if (scale)
0047     ScaledPhi = new MuonSeeddPhiScale(pset);
0048 
0049   // Create the root file
0050   theFile = new TFile(rootFileName.c_str(), "RECREATE");
0051   theFile->mkdir("AllMuonSys");
0052   theFile->cd();
0053   theFile->mkdir("CSC_All");
0054   theFile->cd();
0055   theFile->mkdir("DT_All");
0056   theFile->cd();
0057   theFile->mkdir("ME_All");
0058   theFile->cd();
0059   theFile->mkdir("MB_All");
0060   theFile->cd();
0061   theFile->mkdir("OL_All");
0062   // TTree test
0063 
0064   // The dphi between chamber by chamber
0065   // All possible segment pair in CSC
0066   //                  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14
0067   int csc1[2][15] = {{11, 11, 12, 12, 13, 11, 11, 12, 13, 11, 21, 21, 22, 21, 31},
0068                      {12, 21, 21, 22, 22, 31, 32, 32, 32, 41, 31, 32, 32, 41, 41}};
0069   char ME_nu1[8];
0070   for (int i = 0; i < 15; i++) {
0071     sprintf(ME_nu1, "ME_%d-%d", csc1[0][i], csc1[1][i]);
0072     hME1[i] = new H2DRecHit4(ME_nu1);
0073     cout << "hME1_" << i << " = " << ME_nu1 << endl;
0074   }
0075   // All possible segment pair in DT
0076   //                0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
0077   int dt1[2][26] = {
0078       {10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 11, 12, 12, 20, 20, 20, 21, 21, 21, 21, 22, 22, 30, 31, 31, 32},
0079       {20, 30, 31, 40, 41, 21, 22, 31, 32, 41, 42, 22, 32, 30, 40, 41, 31, 32, 41, 42, 32, 42, 40, 41, 42, 42}};
0080   char MB_nu1[8];
0081   for (int i = 0; i < 26; i++) {
0082     sprintf(MB_nu1, "MB_%d-%d", dt1[0][i], dt1[1][i]);
0083     hMB1[i] = new H2DRecHit5(MB_nu1);
0084     cout << "hMB1_" << i << " = " << MB_nu1 << endl;
0085   }
0086   // All possible segment pair in Overlap region between DT and CSC
0087   //               0  1  2  3  4  5
0088   int olp[2][6] = {{12, 12, 12, 22, 22, 32}, {13, 22, 32, 13, 22, 13}};
0089   char OL_nu1[7];
0090   for (int i = 0; i < 6; i++) {
0091     sprintf(OL_nu1, "OL_%d%d", olp[0][i], olp[1][i]);
0092     hOL1[i] = new H2DRecHit10(OL_nu1);
0093     cout << "hOL_" << i << " = " << OL_nu1 << endl;
0094   }
0095 
0096   // All single chamber segment in CSC
0097   //             0  1  2  3  4  5  6  7
0098   int csc2[8] = {11, 12, 13, 21, 22, 31, 32, 41};
0099   char ME_nu2[6];
0100   for (int i = 0; i < 8; i++) {
0101     sprintf(ME_nu2, "SME_%d", csc2[i]);
0102     hME2[i] = new H2DRecHit6(ME_nu2);
0103     cout << "hME2_" << i << " = " << ME_nu2 << endl;
0104   }
0105 
0106   // All single chamber segment in DT
0107   //           0  1  2  3  4  5  6  7  8  9 10 11
0108   int dt2[12] = {10, 11, 12, 20, 21, 22, 30, 31, 32, 40, 41, 42};
0109   char MB_nu2[6];
0110   for (int i = 0; i < 12; i++) {
0111     sprintf(MB_nu2, "SMB_%d", dt2[i]);
0112     hMB2[i] = new H2DRecHit7(MB_nu2);
0113     cout << "hMB2_" << i << " = " << MB_nu2 << endl;
0114   }
0115   h_all = new H2DRecHit1("AllMu_");
0116   h_csc = new H2DRecHit2("CSC_");
0117   h_dt = new H2DRecHit3("DT_");
0118 }
0119 
0120 // destructor
0121 MuonSeedParametrization::~MuonSeedParametrization() {
0122   if (debug)
0123     cout << "[SeedQualityAnalysis] Destructor called" << endl;
0124   delete recsegSelector;
0125   delete HistoFill;
0126   if (scale)
0127     delete ScaledPhi;
0128   // Write the histos to file
0129   theFile->cd();
0130   theFile->cd("AllMuonSys");
0131   h_all->Write();
0132   theFile->cd();
0133 
0134   theFile->cd("CSC_All");
0135   h_csc->Write();
0136   theFile->cd();
0137 
0138   theFile->cd("DT_All");
0139   h_dt->Write();
0140   theFile->cd();
0141 
0142   theFile->cd("ME_All");
0143   for (int i = 0; i < 15; i++) {
0144     hME1[i]->Write();
0145   }
0146   for (int i = 0; i < 8; i++) {
0147     hME2[i]->Write();
0148   }
0149   theFile->cd();
0150 
0151   theFile->cd("MB_All");
0152   for (int i = 0; i < 26; i++) {
0153     hMB1[i]->Write();
0154   }
0155   for (int i = 0; i < 12; i++) {
0156     hMB2[i]->Write();
0157   }
0158   theFile->cd();
0159 
0160   theFile->cd("OL_All");
0161   for (int i = 0; i < 6; i++) {
0162     hOL1[i]->Write();
0163   }
0164   // for tree
0165   theFile->cd();
0166 
0167   // Release the memory ...
0168   for (int i = 0; i < 15; i++) {
0169     delete hME1[i];
0170   }
0171   for (int i = 0; i < 26; i++) {
0172     delete hMB1[i];
0173   }
0174   for (int i = 0; i < 6; i++) {
0175     delete hOL1[i];
0176   }
0177   for (int i = 0; i < 8; i++) {
0178     delete hME2[i];
0179   }
0180   for (int i = 0; i < 12; i++) {
0181     delete hMB2[i];
0182   }
0183   delete h_all;
0184   delete h_csc;
0185   delete h_dt;
0186 
0187   theFile->Close();
0188   if (debug)
0189     cout << "************* Finished writing histograms to file" << endl;
0190 }
0191 
0192 // The Main...Aanlysis...
0193 
0194 void MuonSeedParametrization::analyze(const Event& event, const EventSetup& eventSetup) {
0195   //Get the CSC Geometry :
0196   ESHandle<CSCGeometry> cscGeom = eventSetup.getHandle(cscGeomToken);
0197 
0198   //Get the DT Geometry :
0199   ESHandle<DTGeometry> dtGeom = eventSetup.getHandle(dtGeomToken);
0200 
0201   // Get the RecHits collection :
0202   Handle<CSCRecHit2DCollection> csc2DRecHits;
0203   event.getByLabel(recHitLabel, csc2DRecHits);
0204 
0205   // Get the CSC Segments collection :
0206   Handle<CSCSegmentCollection> cscSegments;
0207   event.getByLabel(cscSegmentLabel, cscSegments);
0208 
0209   // Get the DT RecHits collection :
0210   Handle<DTRecHitCollection> dt1DRecHits;
0211   event.getByLabel(dtrecHitLabel, dt1DRecHits);
0212 
0213   // Get the DT 4D Segments collection :
0214   Handle<DTRecSegment4DCollection> dt4DSegments;
0215   event.getByLabel(dtSegmentLabel, dt4DSegments);
0216 
0217   // Get the DT 2D Segments collection :
0218   //Handle<DTRecSegment2DCollection> dt2DSegments;
0219   //event.getByLabel(dt2DSegmentLabel, dt2DSegments);
0220 
0221   // Get the SimHit collection :
0222   Handle<PSimHitContainer> csimHits;
0223   event.getByLabel(simHitLabel, "MuonCSCHits", csimHits);
0224   Handle<PSimHitContainer> dsimHits;
0225   event.getByLabel(simHitLabel, "MuonDTHits", dsimHits);
0226 
0227   // Get the SimTrack
0228   Handle<SimTrackContainer> simTracks;
0229   event.getByLabel(simTrackLabel, simTracks);
0230 
0231   H2DRecHit1* histo1 = 0;
0232   H2DRecHit2* histo2 = 0;
0233   H2DRecHit3* histo3 = 0;
0234 
0235   // 0. Run the class SegSelector
0236   //SegSelector recsegSelector();
0237   std::vector<SimSegment> sCSC_v = recsegSelector->Sim_CSCSegments(1, csimHits, cscGeom);
0238   std::vector<SimSegment> sDT_v = recsegSelector->Sim_DTSegments(1, dsimHits, dtGeom);
0239   std::vector<CSCSegment> cscseg_V = recsegSelector->Select_CSCSeg(cscSegments, cscGeom, sCSC_v);
0240   std::vector<DTRecSegment4D> dtseg_V = recsegSelector->Select_DTSeg(dt4DSegments, dtGeom, sDT_v);
0241 
0242   // 1. Get the track information
0243   SimInfo(simTracks, dsimHits, csimHits, dtGeom, cscGeom);
0244 
0245   // 2. Check # of segments in each chambers
0246   CSCsegment_stat(cscSegments);
0247   DTsegment_stat(dt4DSegments);
0248   if ((cscseg_stat[5] < 2) && (dtseg_stat[5] < 2)) {
0249     /// check # of rechit for those 1 or 0 segment event
0250     CSCRecHit_stat(csc2DRecHits, cscGeom);
0251     DTRecHit_stat(dt1DRecHits, dtGeom);
0252   }
0253   // 3. Event statistic
0254   // number of segments or rechits w.r.t. eta in each event
0255   int allmu_stat = cscseg_stat[5] + dtseg_stat[5];
0256   int allrh_stat = cscrh_sum[5] + dtrh_sum[5];
0257   double allmu_eta = -9.0;
0258   if ((cscseg_stat[5] == 0) && (dtseg_stat[5] != 0)) {
0259     allmu_eta = eta_d;
0260   } else if ((cscseg_stat[5] != 0) && (dtseg_stat[5] == 0)) {
0261     allmu_eta = eta_c;
0262   } else if ((cscseg_stat[5] != 0) && (dtseg_stat[5] != 0)) {
0263     allmu_eta = (eta_c + eta_d) / 2;
0264   } else {
0265     // look up the rechits vs eta in events with no segments
0266     if ((eta_d == -9.0) && (eta_c != -9.0)) {
0267       allmu_eta = eta_c;
0268     } else if ((eta_d != -9.0) && (eta_c == -9.0)) {
0269       allmu_eta = eta_d;
0270     } else {
0271       allmu_eta = (eta_c + eta_d) / 2;
0272     }
0273 
0274     histo1 = h_all;
0275     histo1->Fill1a(allmu_eta, allrh_stat, cscrh_sum[0], dtrh_sum[0]);
0276   }
0277   histo1 = h_all;
0278   histo1->Fill1(cscseg_stat[5], dtseg_stat[5], allmu_stat, eta_c, eta_d, allmu_eta, eta_trk);
0279   // check the energy/pt loss in each layer
0280 
0281   // look up the # of segments in each station
0282   if (cscseg_stat[0] != 0) {
0283     histo2 = h_csc;
0284     histo2->Fill3(pt1[0], pa[0], cscseg_stat[0]);
0285     for (int k = 1; k < 5; k++) {
0286       histo2->Fill4(k, cscseg_stat[k], cscseg_stat1[k]);
0287     }
0288     //if (pt1[0] < 50) {
0289     if (etaLc[1] != 0.0) {
0290       histo1->Fill1c1(etaLc[1], ptLossC[1]);
0291     }
0292     if (etaLc[2] != 0.0) {
0293       histo1->Fill1c2(etaLc[2], ptLossC[2], pt1[0]);
0294     }
0295     if (etaLc[3] != 0.0) {
0296       histo1->Fill1c3(etaLc[3], ptLossC[3]);
0297     }
0298     if (etaLc[4] != 0.0) {
0299       histo1->Fill1c4(etaLc[4], ptLossC[4]);
0300     }
0301     //}
0302   }
0303   if (dtseg_stat[0] != 0) {
0304     histo3 = h_dt;
0305     histo3->Fill3a(pt1[0], pa[0], dtseg_stat[0]);
0306     for (int k = 1; k < 5; k++) {
0307       histo3->Fill4a(k, dtseg_stat[k], dtseg_stat1[k]);
0308     }
0309     //if (pt1[0] < 50) {
0310     if (etaLd[1] != 0.0) {
0311       histo1->Fill1d1(etaLd[1], ptLossD[1], pt1[0]);
0312     }
0313     if (etaLd[2] != 0.0) {
0314       histo1->Fill1d2(etaLd[2], ptLossD[2]);
0315     }
0316     if (etaLd[3] != 0.0) {
0317       histo1->Fill1d3(etaLd[3], ptLossD[3]);
0318     }
0319     if (etaLd[4] != 0.0) {
0320       histo1->Fill1d4(etaLd[4], ptLossD[4]);
0321     }
0322     //}
0323   }
0324 
0325   /*
0326     NOT USED
0327     // flag the overlap case
0328     bool overlap = false;
0329     if ( (dtseg_stat[0] != 0)&&(cscseg_stat[0] != 0) ){
0330     overlap = true;
0331     }
0332   */
0333 
0334   // 4. the simulated segments statistic
0335   int simcscseg[6] = {0};
0336   double ns1 = 0.0;
0337   double eta_sim1 = 0;
0338   for (std::vector<SimSegment>::const_iterator it = sCSC_v.begin(); it != sCSC_v.end(); it++) {
0339     int st = ((*it).csc_DetId).station();
0340     eta_sim1 += fabs(((*it).sGlobalOrg).eta());
0341     simcscseg[st]++;
0342     ns1++;
0343   }
0344   simcscseg[0] = simcscseg[1] + simcscseg[2] + simcscseg[3] + simcscseg[4];
0345   for (int i = 1; i < 5; i++) {
0346     if (simcscseg[i] != 0) {
0347       simcscseg[5]++;
0348     }
0349   }
0350   eta_sim1 = eta_sim1 / ns1;
0351 
0352   int simdtseg[6] = {0};
0353   double ns2 = 0.0;
0354   double eta_sim2 = 0;
0355   for (std::vector<SimSegment>::const_iterator it = sDT_v.begin(); it != sDT_v.end(); it++) {
0356     int st = ((*it).dt_DetId).station();
0357     eta_sim2 += fabs(((*it).sGlobalOrg).eta());
0358     simdtseg[st]++;
0359     ns2++;
0360   }
0361   simdtseg[0] = simdtseg[1] + simdtseg[2] + simdtseg[3] + simdtseg[4];
0362   for (int i = 1; i < 5; i++) {
0363     if (simdtseg[i] != 0) {
0364       simdtseg[5]++;
0365     }
0366   }
0367   eta_sim2 = eta_sim2 / ns2;
0368 
0369   int allmu_stat1 = simcscseg[5] + simdtseg[5];
0370   double allmu_eta1 = -9.0;
0371   if ((simcscseg[0] == 0) && (simdtseg[0] != 0)) {
0372     allmu_eta1 = eta_sim2;
0373   } else if ((simdtseg[0] == 0) && (simcscseg[0] != 0)) {
0374     allmu_eta1 = eta_sim1;
0375   } else {
0376     allmu_eta1 = (eta_sim1 + eta_sim2) / 2;
0377   }
0378 
0379   histo1 = h_all;
0380   histo1->Fill1b(allmu_stat1, allmu_eta1);
0381 
0382   //5.  look at different Bxdl btw. stations
0383   FromCSCSeg(cscseg_V, cscGeom, sCSC_v);
0384   FromDTSeg(dtseg_V, dtGeom, sDT_v);
0385   FromOverlap();
0386   FromCSCSingleSeg(cscseg_V, cscGeom, sCSC_v);
0387   FromDTSingleSeg(dtseg_V, dtGeom, sDT_v);
0388 
0389   /// chi2 distribution
0390   for (int i = 0; i < 5; i++) {
0391     if (chi2_dof1[i] < 0.0)
0392       continue;
0393     histo2 = h_csc;
0394     histo2->Fill3b(chi2_dof1[i]);
0395   }
0396   for (int i = 1; i < 5; i++) {
0397     if (chi2_dof3[i] < 0.0)
0398       continue;
0399     histo3 = h_dt;
0400     histo3->Fill3c(chi2_dof3[i]);
0401   }
0402 
0403   // Scale the dphi
0404   if (scale) {
0405     ScaledPhi->ScaleCSCdPhi(dPhiP1, EtaP1);
0406     ScaledPhi->ScaleDTdPhi(dPhiP3, EtaP3);
0407     ScaledPhi->ScaleOLdPhi(dPhiP2, MBPath, MEPath);
0408     ScaledPhi->ScaleMBSingle(MB_phi, MBPath);
0409     ScaledPhi->ScaleMESingle(ME_phi, MEPath);
0410   }
0411   // fill the information for CSC pT parameterization from segment pair
0412   histo2 = h_csc;
0413   HistoFill->FillCSCSegmentPair(histo2, pt1, chi2_dof1, dPhiP1, EtaP1);
0414 
0415   /// For DT
0416   //  fill the information for DT pT parameterization from segment pair
0417   histo3 = h_dt;
0418   HistoFill->FillDTSegmentPair(histo3, pt1, chi2_dof3, dPhiP3, EtaP3);
0419 
0420   //  Look at different Bxdl btw. stations & rings
0421   HistoFill->FillCSCSegmentPairByChamber(hME1, pt1, dPhiP1, EtaP1, MEPath, dEtaP1);
0422 
0423   HistoFill->FillDTSegmentPairByChamber(hMB1, pt1, dPhiP3, EtaP3, MBPath, dEtaP3);
0424 
0425   HistoFill->FillOLSegmentPairByChamber(hOL1, pt1, dPhiP2, EtaP3, MBPath, MEPath, dEtaP2);
0426 
0427   HistoFill->FillCSCSegmentSingle(hME2, pt1, ME_phi, ME_eta, MEPath);
0428 
0429   HistoFill->FillDTSegmentSingle(hMB2, pt1, MB_phi, MB_eta, MBPath);
0430 
0431   /*
0432   /// For reco-segment treea
0433   tt = tr_muon;
0434   if ( MEPath[1][1] && MEPath[1][2] && MEPath[1][3] ) {
0435       tt->Fill_b1(fabs(EtaP1[1][1]),fabs(EtaP1[1][2]),fabs(EtaP1[1][3]),fabs(EtaP1[1][4]), 
0436                   EtaP1[1][1], EtaP1[1][2], EtaP1[1][3], EtaP1[1][4], pt1[0]);
0437       tt->Fill_l1(pa[0]);
0438   }
0439   if ( MBPath[1][1] && MBPath[1][2] && MBPath[1][3] ) {
0440       tt->Fill_b2(fabs(EtaP3[1][1]),fabs(EtaP3[1][2]),fabs(EtaP3[1][3]),fabs(EtaP3[1][4]), 
0441                   EtaP3[1][1], EtaP3[1][2], EtaP3[1][3], EtaP3[1][4], pt1[0]);
0442       tt->Fill_l1(pa[0]);
0443   }
0444   tt->FillTree();
0445   */
0446 }
0447 
0448 // ********************************************
0449 // ***********  Utility functions  ************
0450 // ********************************************
0451 
0452 // number of csc segments in one chamber for each station
0453 // cscseg_stat[0] = total segments in all stations
0454 // cscseg_stat[i] = the number of segments in station i
0455 // cscseg_stat[5] = the number of stations which have segments
0456 // cscseg_stat1 is the statistic for segments w/ more than 4 rechits
0457 void MuonSeedParametrization::CSCsegment_stat(Handle<CSCSegmentCollection> cscSeg) {
0458   for (int i = 0; i < 6; i++) {
0459     cscseg_stat[i] = 0;
0460     cscseg_stat1[i] = 0;
0461   }
0462   for (CSCSegmentCollection::const_iterator seg_It = cscSeg->begin(); seg_It != cscSeg->end(); seg_It++) {
0463     CSCDetId DetId = (CSCDetId)(*seg_It).cscDetId();
0464     cscseg_stat[DetId.station()] += 1;
0465     if ((*seg_It).nRecHits() > 4) {
0466       cscseg_stat1[DetId.station()] += 1;
0467     }
0468   }
0469   cscseg_stat[0] = cscseg_stat[1] + cscseg_stat[2] + cscseg_stat[3] + cscseg_stat[4];
0470   cscseg_stat1[0] = cscseg_stat1[1] + cscseg_stat1[2] + cscseg_stat1[3] + cscseg_stat1[4];
0471   for (int i = 1; i < 5; i++) {
0472     if (cscseg_stat[i] != 0) {
0473       cscseg_stat[5]++;
0474     }
0475     if (cscseg_stat1[i] != 0) {
0476       cscseg_stat1[5]++;
0477     }
0478   }
0479 }
0480 // the same as CSCsegment_stat
0481 void MuonSeedParametrization::DTsegment_stat(Handle<DTRecSegment4DCollection> dtSeg) {
0482   for (int i = 0; i < 6; i++) {
0483     dtseg_stat[i] = 0;
0484     dtseg_stat1[i] = 0;
0485     dt2Dseg_stat[i] = 0;
0486   }
0487   for (DTRecSegment4DCollection::const_iterator seg_It = dtSeg->begin(); seg_It != dtSeg->end(); seg_It++) {
0488     if ((*seg_It).hasPhi() && (*seg_It).hasZed()) {
0489       DTChamberId DId1 = (*seg_It).chamberId();
0490       dtseg_stat[DId1.station()] += 1;
0491       int n_phiHits = ((*seg_It).phiSegment())->specificRecHits().size();
0492       if (n_phiHits > 4) {
0493         dtseg_stat1[DId1.station()] += 1;
0494       }
0495     }
0496     if ((*seg_It).hasPhi() && !(*seg_It).hasZed()) {
0497       const DTChamberRecSegment2D* phiSeg = (*seg_It).phiSegment();
0498       DetId geoId = (phiSeg)->geographicalId();
0499       DTChamberId DId2 = DTChamberId(geoId);
0500       dt2Dseg_stat[DId2.station()] += 1;
0501     }
0502   }
0503 
0504   dtseg_stat[0] = dtseg_stat[1] + dtseg_stat[2] + dtseg_stat[3] + dt2Dseg_stat[4];
0505   dtseg_stat1[0] = dtseg_stat1[1] + dtseg_stat1[2] + dtseg_stat1[3] + dtseg_stat1[4];
0506 
0507   for (int i = 1; i < 5; i++) {
0508     if (dtseg_stat[i] != 0 || dt2Dseg_stat[4] != 0) {
0509       dtseg_stat[5]++;
0510     }
0511     if (dtseg_stat1[i] != 0 || dt2Dseg_stat[4] != 0) {
0512       dtseg_stat1[5]++;
0513     }
0514   }
0515 }
0516 // number rechit in each station, basically only for those station cannot form a segment
0517 void MuonSeedParametrization::CSCRecHit_stat(Handle<CSCRecHit2DCollection> cscrechit, ESHandle<CSCGeometry> cscGeom) {
0518   for (int i = 0; i < 6; i++) {
0519     cscrh_sum[i] = 0;
0520   }
0521   for (CSCRecHit2DCollection::const_iterator r_it = cscrechit->begin(); r_it != cscrechit->end(); r_it++) {
0522     CSCDetId det_id = (CSCDetId)(*r_it).cscDetId();
0523     //const CSCLayer* csclayer = cscGeom->layer( det_id );
0524     //const CSCChamber* cscchamber = cscGeom->chamber( det_id );
0525     //LocalPoint lrh = (*r_it).localPosition();
0526     //GlobalPoint grh = csclayer->toGlobal(lrh);
0527     cscrh_sum[det_id.station()]++;
0528   }
0529   cscrh_sum[0] = cscrh_sum[1] + cscrh_sum[2] + cscrh_sum[3] + cscrh_sum[4];
0530   for (int i = 1; i < 5; i++) {
0531     if (cscrh_sum[i] != 0) {
0532       cscrh_sum[5]++;
0533     }
0534   }
0535 }
0536 
0537 void MuonSeedParametrization::DTRecHit_stat(Handle<DTRecHitCollection> dtrechit, ESHandle<DTGeometry> dtGeom) {
0538   //double phi[4]={999.0};
0539   for (int i = 0; i < 6; i++) {
0540     dtrh_sum[i] = 0;
0541   }
0542 
0543   double eta = -9.0;
0544   double nn = 0.0;
0545   for (DTRecHitCollection::const_iterator r_it = dtrechit->begin(); r_it != dtrechit->end(); r_it++) {
0546     DTWireId det_id = (*r_it).wireId();
0547     const DTChamber* dtchamber = dtGeom->chamber(det_id);
0548     LocalPoint lrh = (*r_it).localPosition();
0549     GlobalPoint grh = dtchamber->toGlobal(lrh);
0550     dtrh_sum[det_id.station()]++;
0551     eta += grh.eta();
0552     nn += 1.0;
0553   }
0554   eta = eta / nn;
0555 
0556   dtrh_sum[0] = dtrh_sum[1] + dtrh_sum[2] + dtrh_sum[3] + dtrh_sum[4];
0557   for (int i = 1; i < 5; i++) {
0558     if (dtrh_sum[i] != 0) {
0559       dtrh_sum[5]++;
0560     }
0561   }
0562 }
0563 
0564 // find the simHits which is corresponding to the segment
0565 bool MuonSeedParametrization::SameChamber(CSCDetId SimDetId, CSCDetId SegDetId) {
0566   if (SimDetId.endcap() == SegDetId.endcap() && SimDetId.station() == SegDetId.station() &&
0567       SimDetId.ring() == SegDetId.ring() && SimDetId.chamber() == SegDetId.chamber()) {
0568     return true;
0569   } else {
0570     return false;
0571   }
0572 }
0573 
0574 void MuonSeedParametrization::SimInfo(Handle<edm::SimTrackContainer> simTracks,
0575                                       Handle<edm::PSimHitContainer> dsimHits,
0576                                       Handle<edm::PSimHitContainer> csimHits,
0577                                       ESHandle<DTGeometry> dtGeom,
0578                                       ESHandle<CSCGeometry> cscGeom) {
0579   // pt1 -> pt in different station pt1[0] is the track pt
0580   for (int i = 0; i < 5; i++) {
0581     pt1[i] = 0.0;
0582     pa[i] = 0.0;
0583     etaLc[i] = 0.0;
0584     etaLd[i] = 0.0;
0585     ptLossC[i] = 0.0;
0586     ptLossD[i] = 0.0;
0587   }
0588 
0589   // eta_c -> ave.eta from all csc stations ; cta_d -> ave.eta from dt stations
0590   eta_c = -9.0;
0591   eta_d = -9.0;
0592   eta_trk = -9.0;
0593 
0594   for (SimTrackContainer::const_iterator simTk_It = simTracks->begin(); simTk_It != simTracks->end(); simTk_It++) {
0595     if (abs((*simTk_It).type()) != 13)
0596       continue;
0597 
0598     if ((*simTk_It).type() == 13) {
0599       theQ = -1.0;
0600     } else {
0601       theQ = 1.0;
0602     }
0603 
0604     float px = ((*simTk_It).momentum()).x();
0605     float py = ((*simTk_It).momentum()).y();
0606     float pz = ((*simTk_It).momentum()).z();
0607     pt1[0] = sqrt((px * px) + (py * py));
0608     pa[0] = sqrt((px * px) + (py * py) + (pz * pz));
0609     double theta = acos(pz / pa[0]);
0610     eta_trk = fabs((-1.0) * log(tan(theta / 2.0)));
0611 
0612     double eta_c1 = 0.0;
0613     double enu1 = 0.0;
0614     for (PSimHitContainer::const_iterator cs_It = csimHits->begin(); cs_It != csimHits->end(); cs_It++) {
0615       CSCDetId C_Id = CSCDetId((*cs_It).detUnitId());
0616       const CSCChamber* cscchamber = cscGeom->chamber(C_Id);
0617       GlobalVector m1 = cscchamber->toGlobal((*cs_It).momentumAtEntry());
0618       Local3DPoint lp = (*cs_It).localPosition();
0619       GlobalPoint gp = cscchamber->toGlobal(lp);
0620       if ((abs((*cs_It).particleType()) == 13) && ((*cs_It).trackId() == 1)) {
0621         pt1[C_Id.station()] = sqrt((m1.x() * m1.x()) + (m1.y() * m1.y()));
0622         pa[C_Id.station()] = sqrt((m1.x() * m1.x()) + (m1.y() * m1.y()) + (m1.z() * m1.z()));
0623         etaLc[C_Id.station()] = fabs(gp.eta());
0624         ptLossC[C_Id.station()] = pt1[C_Id.station()] / pt1[0];
0625         eta_c1 += fabs(gp.eta());
0626         enu1 += 1.0;
0627       }
0628     }
0629     if (enu1 != 0.0) {
0630       eta_c = eta_c1 / enu1;
0631     } else {
0632       eta_c = -9.0;
0633     }
0634 
0635     double eta_d1 = 0.0;
0636     double enu2 = 0.0;
0637     for (PSimHitContainer::const_iterator ds_It = dsimHits->begin(); ds_It != dsimHits->end(); ds_It++) {
0638       Local3DPoint lp = (*ds_It).localPosition();
0639 
0640       DTLayerId D_Id = DTLayerId((*ds_It).detUnitId());
0641       const DTLayer* dtlayer = dtGeom->layer(D_Id);
0642       GlobalVector m2 = dtlayer->toGlobal((*ds_It).momentumAtEntry());
0643       GlobalPoint gp = dtlayer->toGlobal(lp);
0644 
0645       if ((abs((*ds_It).particleType()) == 13) && ((*ds_It).trackId() == 1)) {
0646         pt1[D_Id.station()] = sqrt((m2.x() * m2.x()) + (m2.y() * m2.y()));
0647         pa[D_Id.station()] = sqrt((m2.x() * m2.x()) + (m2.y() * m2.y()) + (m2.z() * m2.z()));
0648         etaLd[D_Id.station()] = fabs(gp.eta());
0649         ptLossD[D_Id.station()] = pt1[D_Id.station()] / pt1[0];
0650         eta_d1 += fabs(gp.eta());
0651         enu2 += 1.0;
0652       }
0653     }
0654     if (enu2 != 0.0) {
0655       eta_d = eta_d1 / enu2;
0656     } else {
0657       eta_d = -9.0;
0658     }
0659   }
0660 }
0661 
0662 // Fill the phi and eta information for CSC
0663 void MuonSeedParametrization::FromCSCSeg(std::vector<CSCSegment> cscSeg,
0664                                          ESHandle<CSCGeometry> cscGeom,
0665                                          std::vector<SimSegment> seg) {
0666   /// get the global dphi from recHits for CSC
0667   //unused     double Theta[2][5] ;
0668   for (int l = 0; l < 10; l++) {
0669     int i = l / 2;  // 0 0 1 1 2 2 3 3 4 4
0670     int k = l % 2;  // 0 1 0 1 0 1 0 1 0 1
0671     PhiV1[k][i] = 99.;
0672     EtaV1[k][i] = 99.;
0673     PhiP1[k][i] = 99.;
0674     EtaP1[k][i] = 99.;
0675     //unused        Theta[k][i]=99.;
0676     chi2_dof1[i] = -1.0;
0677     for (int j = 0; j < 5; j++) {
0678       dPhiV1[k][i][j] = 99.;
0679       dEtaV1[k][i][j] = 99.;
0680       dPhiP1[k][i][j] = 99.;
0681       dEtaP1[k][i][j] = 99.;
0682     }
0683   }
0684   bool layer[5] = {false};
0685   // Fill the phi and eta of segment direction in different station
0686   for (std::vector<CSCSegment>::const_iterator it = cscSeg.begin(); it != cscSeg.end(); it++) {
0687     CSCDetId DetId = (CSCDetId)(*it).cscDetId();
0688     const CSCChamber* cscchamber = cscGeom->chamber(DetId);
0689     GlobalPoint gp = cscchamber->toGlobal((*it).localPosition());
0690     GlobalVector gv = cscchamber->toGlobal((*it).localDirection());
0691     int st = DetId.station();
0692     int rg = DetId.ring();
0693 
0694     if (st == 1 && (rg == 1 || rg == 4)) {
0695       PhiV1[1][0] = gv.phi();
0696       EtaV1[1][0] = gv.eta();
0697       PhiP1[1][0] = gp.phi();
0698       EtaP1[1][0] = gp.eta();
0699       //unused            Theta[1][0] = gp.theta();
0700       layer[0] = true;
0701       chi2_dof1[st] = (*it).chi2() / static_cast<double>((*it).degreesOfFreedom());
0702     } else {
0703       PhiV1[1][st] = gv.phi();
0704       EtaV1[1][st] = gv.eta();
0705       PhiP1[1][st] = gp.phi();
0706       EtaP1[1][st] = gp.eta();
0707       //unused            Theta[1][st] = gp.theta();
0708       layer[st] = true;
0709       chi2_dof1[st] = (*it).chi2() / static_cast<double>((*it).degreesOfFreedom());
0710     }
0711   }
0712   for (std::vector<SimSegment>::const_iterator it = seg.begin(); it != seg.end(); it++) {
0713     if ((*it).chamber_type != 1)
0714       continue;
0715     GlobalPoint gp = (*it).sGlobalOrg;
0716     GlobalVector gv = (*it).sGlobalVec;
0717     int st = (*it).csc_DetId.station();
0718     int rg = (*it).csc_DetId.ring();
0719 
0720     if (st == 1 && (rg == 1 || rg == 4)) {
0721       PhiV1[0][0] = gv.phi();
0722       EtaV1[0][0] = gv.eta();
0723       PhiP1[0][0] = gp.phi();
0724       EtaP1[0][0] = gp.eta();
0725       //unused            Theta[0][0] = gp.theta();
0726     } else {
0727       PhiV1[0][st] = gv.phi();
0728       EtaV1[0][st] = gv.eta();
0729       PhiP1[0][st] = gp.phi();
0730       EtaP1[0][st] = gp.eta();
0731       //unused            Theta[0][st] = gp.theta();
0732     }
0733   }
0734 
0735   // Get the dPhi and dEta for different combination
0736   for (int m = 0; m < 2; m++) {
0737     for (int l = 0; l < 16; l++) {
0738       int s1 = (l % 4);      // 0  0,1  0,1,2  0,1,2,3
0739       int s2 = (l / 4) + 1;  // 1  2,2  3,3,3  4,4,4,4
0740       if (layer[s1] && layer[s2] && (s1 < s2)) {
0741         dPhiV1[m][s1][s2] = (PhiV1[m][s1] - PhiV1[m][s2]);
0742         dEtaV1[m][s1][s2] = EtaV1[m][s1] - EtaV1[m][s2];
0743         dPhiP1[m][s1][s2] = (PhiP1[m][s1] - PhiP1[m][s2]);
0744         dEtaP1[m][s1][s2] = EtaP1[m][s1] - EtaP1[m][s2];
0745       }
0746     }
0747   }
0748 }
0749 
0750 // Fill the phi and eta information for DT
0751 void MuonSeedParametrization::FromDTSeg(std::vector<DTRecSegment4D> dtSeg,
0752                                         ESHandle<DTGeometry> dtGeom,
0753                                         std::vector<SimSegment> seg) {
0754   /// get the global dphi from recHits for DT
0755   for (int l = 0; l < 10; l++) {
0756     int i = l / 2;  // 0 0 1 1 2 2 3 3 4 4
0757     int k = l % 2;  // 0 1 0 1 0 1 0 1 0 1
0758     PhiV3[k][i] = 99.;
0759     EtaV3[k][i] = 99.;
0760     PhiP3[k][i] = 99.;
0761     EtaP3[k][i] = 99.;
0762     chi2_dof3[i] = -1.0;
0763     for (int j = 0; j < 5; j++) {
0764       dPhiV3[k][i][j] = 99.;
0765       dEtaV3[k][i][j] = 99.;
0766       dPhiP3[k][i][j] = 99.;
0767       dEtaP3[k][i][j] = 99.;
0768     }
0769   }
0770   bool layer[5] = {false};
0771   // Fill the phi and eta of segment direction in different station
0772   for (std::vector<DTRecSegment4D>::const_iterator it = dtSeg.begin(); it != dtSeg.end(); it++) {
0773     if (!(*it).hasPhi())
0774       continue;
0775     DetId geoId = (*it).geographicalId();
0776 
0777     if ((*it).hasPhi() && !(*it).hasZed()) {
0778       const DTChamberRecSegment2D* phiSeg = (*it).phiSegment();
0779       geoId = (phiSeg)->geographicalId();
0780     }
0781 
0782     DTChamberId cbId = DTChamberId(geoId);
0783     const DTChamber* dtchamber = dtGeom->chamber(cbId);
0784 
0785     GlobalPoint gp = dtchamber->toGlobal((*it).localPosition());
0786     GlobalVector gv = dtchamber->toGlobal((*it).localDirection());
0787     int st = cbId.station();
0788     PhiV3[1][st] = gv.phi();
0789     EtaV3[1][st] = gv.eta();
0790     PhiP3[1][st] = gp.phi();
0791     EtaP3[1][st] = gp.eta();
0792     layer[st] = true;
0793     chi2_dof3[st] = (*it).chi2() / static_cast<double>((*it).degreesOfFreedom());
0794   }
0795   for (std::vector<SimSegment>::const_iterator it = seg.begin(); it != seg.end(); it++) {
0796     if ((*it).chamber_type != 2)
0797       continue;
0798     GlobalPoint gp = (*it).sGlobalOrg;
0799     GlobalVector gv = (*it).sGlobalVec;
0800     int st = (*it).dt_DetId.station();
0801     PhiV3[0][st] = gv.phi();
0802     EtaV3[0][st] = gv.eta();
0803     PhiP3[0][st] = gp.phi();
0804     EtaP3[0][st] = gp.eta();
0805   }
0806   // Get the dPhi and dEta for different combination
0807   for (int m = 0; m < 2; m++) {
0808     for (int l = 0; l < 9; l++) {
0809       int s1 = (l % 3) + 1;  // 1  1,2  1,2,3
0810       int s2 = (l / 3) + 2;  // 2  3,3  4,4,4
0811       if (layer[s1] && layer[s2] && (s1 < s2)) {
0812         dPhiV3[m][s1][s2] = PhiV3[m][s1] - PhiV3[m][s2];
0813         dEtaV3[m][s1][s2] = EtaV3[m][s1] - EtaV3[m][s2];
0814         dPhiP3[m][s1][s2] = PhiP3[m][s1] - PhiP3[m][s2];
0815         dEtaP3[m][s1][s2] = EtaP3[m][s1] - EtaP3[m][s2];
0816         //cout <<" ---------------------- DT "<<s1<<"_"<<s2<<"-------------------------------------"<<endl;
0817         //cout <<"dPhi from P = "<< PhiP3[m][s1] <<" - "<<PhiP3[m][s2]<<" = "<<dPhiP3[m][s1][s2]<<endl;
0818         //cout <<"dPhi from V = "<< PhiV3[m][s1] <<" - "<<PhiV3[m][s2]<<" = "<<dPhiV3[m][s1][s2]<<endl;
0819       }
0820     }
0821   }
0822 }
0823 
0824 void MuonSeedParametrization::FromOverlap() {
0825   for (int l = 0; l < 10; l++) {
0826     int i = l / 2;
0827     int k = l % 2;
0828     for (int j = 0; j < 5; j++) {
0829       dPhiV2[k][i][j] = 99.;
0830       dEtaV2[k][i][j] = 99.;
0831       dPhiP2[k][i][j] = 99.;
0832       dEtaP2[k][i][j] = 99.;
0833     }
0834   }
0835   for (int m = 0; m < 2; m++) {
0836     for (int l = 0; l < 9; l++) {
0837       int s1 = (l % 3) + 1;  // 1,2,3  1,2,3  1,2,3
0838       int s2 = (l / 3) + 1;  // 1,1,1  2,2,2  3,3,3
0839 
0840       if ((PhiV3[m][s1] == 99.0) || (PhiV1[m][s2] == 99.0))
0841         continue;
0842       dPhiV2[m][s1][s2] = PhiV3[m][s1] - PhiV1[m][s2];
0843       dEtaV2[m][s1][s2] = EtaV3[m][s1] - EtaV1[m][s2];
0844       dPhiP2[m][s1][s2] = PhiP3[m][s1] - PhiP1[m][s2];
0845       dEtaP2[m][s1][s2] = EtaP3[m][s1] - EtaP1[m][s2];
0846     }
0847   }
0848 }
0849 
0850 void MuonSeedParametrization::FromCSCSingleSeg(std::vector<CSCSegment> cscSeg,
0851                                                ESHandle<CSCGeometry> cscGeom,
0852                                                std::vector<SimSegment> seg) {
0853   for (int i = 0; i < 2; i++) {
0854     for (int j = 0; j < 5; j++) {
0855       for (int k = 0; k < 4; k++) {
0856         MEPath[i][j][k] = false;
0857         ME_phi[i][j][k] = 99.;
0858         ME_eta[i][j][k] = 99.;
0859       }
0860     }
0861   }
0862   for (std::vector<CSCSegment>::const_iterator it = cscSeg.begin(); it != cscSeg.end(); it++) {
0863     CSCDetId DetId = (CSCDetId)(*it).cscDetId();
0864     const CSCChamber* cscchamber = cscGeom->chamber(DetId);
0865     GlobalPoint gp = cscchamber->toGlobal((*it).localPosition());
0866     GlobalVector gv = cscchamber->toGlobal((*it).localDirection());
0867     int st = DetId.station();
0868     int rg = DetId.ring();
0869     if (rg == 4) {
0870       rg = 1;
0871     }
0872     if (st == 1 && rg == 1) {
0873       st = 0;
0874     }
0875 
0876     // for single-chamber segment in csc
0877     double ab = (gv.x() * gp.x()) + (gv.y() * gp.y());
0878     double al = sqrt((gp.x() * gp.x()) + (gp.y() * gp.y()));
0879     double bl = sqrt((gv.x() * gv.x()) + (gv.y() * gv.y()));
0880     double axb = (gp.x() * gv.y()) - (gp.y() * gv.x());
0881     double cc = (axb < 0.) ? 1.0 : -1.0;
0882 
0883     ME_phi[1][st][rg] = cc * acos(ab / (al * bl));
0884     if (ME_phi[1][st][rg] > 1.570796) {
0885       ME_phi[1][st][rg] = 3.141592 - ME_phi[1][st][rg];
0886     }
0887     ME_eta[1][st][rg] = fabs(gp.eta());
0888     MEPath[1][st][rg] = true;
0889   }
0890   for (std::vector<SimSegment>::const_iterator it = seg.begin(); it != seg.end(); it++) {
0891     int st = ((*it).csc_DetId).station();
0892     int rg = ((*it).csc_DetId).ring();
0893     if (rg == 4) {
0894       rg = 1;
0895     }
0896     if (st == 1 && rg == 1) {
0897       st = 0;
0898     }
0899 
0900     // for single-chamber segment in csc
0901     double ab = (((*it).sGlobalVec).x() * ((*it).sGlobalOrg).x()) + (((*it).sGlobalVec).y() * ((*it).sGlobalOrg).y());
0902     double al =
0903         sqrt((((*it).sGlobalOrg).x() * ((*it).sGlobalOrg).x()) + (((*it).sGlobalOrg).y() * ((*it).sGlobalOrg).y()));
0904     double bl =
0905         sqrt((((*it).sGlobalVec).x() * ((*it).sGlobalVec).x()) + (((*it).sGlobalVec).y() * ((*it).sGlobalVec).y()));
0906     double axb = (((*it).sGlobalOrg).x() * ((*it).sGlobalVec).y()) - (((*it).sGlobalOrg).y() * ((*it).sGlobalVec).x());
0907     double cc = (axb < 0.) ? 1.0 : -1.0;
0908     ME_phi[0][st][rg] = cc * acos(ab / (al * bl));
0909 
0910     if (ME_phi[0][st][rg] > 1.570796) {
0911       ME_phi[0][st][rg] = 3.141592 - ME_phi[0][st][rg];
0912     }
0913     ME_eta[0][st][rg] = fabs(((*it).sGlobalOrg).eta());
0914     MEPath[0][st][rg] = true;
0915   }
0916 }
0917 
0918 void MuonSeedParametrization::FromDTSingleSeg(std::vector<DTRecSegment4D> dtSeg,
0919                                               ESHandle<DTGeometry> dtGeom,
0920                                               std::vector<SimSegment> seg) {
0921   for (int i = 0; i < 2; i++) {
0922     for (int j = 0; j < 5; j++) {
0923       for (int k = 0; k < 3; k++) {
0924         MBPath[i][j][k] = false;
0925         MB_phi[i][j][k] = 99.;
0926         MB_eta[i][j][k] = 99.;
0927       }
0928     }
0929   }
0930   for (std::vector<DTRecSegment4D>::const_iterator it = dtSeg.begin(); it != dtSeg.end(); it++) {
0931     DetId geoId = (*it).geographicalId();
0932 
0933     if ((*it).hasPhi() && !(*it).hasZed()) {
0934       const DTChamberRecSegment2D* phiSeg = (*it).phiSegment();
0935       geoId = (phiSeg)->geographicalId();
0936     }
0937     if (!(*it).hasPhi())
0938       continue;
0939 
0940     DTChamberId DId = DTChamberId(geoId);
0941     const DTChamber* dtchamber = dtGeom->chamber(DId);
0942 
0943     int st = DId.station();
0944     int wl = abs(DId.wheel());
0945 
0946     MBPath[1][st][wl] = true;
0947 
0948     if ((*it).dimension() != 4)
0949       continue;
0950 
0951     GlobalPoint g_seg_o = dtchamber->toGlobal((*it).localPosition());
0952     GlobalVector g_seg_v = dtchamber->toGlobal((*it).localDirection());
0953 
0954     // for single-chamber segment in dt
0955     double ab = (g_seg_v.x() * g_seg_o.x()) + (g_seg_v.y() * g_seg_o.y());
0956     double al = sqrt((g_seg_o.x() * g_seg_o.x()) + (g_seg_o.y() * g_seg_o.y()));
0957     double bl = sqrt((g_seg_v.x() * g_seg_v.x()) + (g_seg_v.y() * g_seg_v.y()));
0958     double axb = (g_seg_o.x() * g_seg_v.y()) - (g_seg_o.y() * g_seg_v.x());
0959     double cc = (axb < 0.) ? 1.0 : -1.0;
0960 
0961     MB_phi[1][st][wl] = cc * acos(ab / (al * bl));
0962     MB_eta[1][st][wl] = fabs(g_seg_o.eta());
0963   }
0964   /// Single segments dphi
0965   for (std::vector<SimSegment>::const_iterator it = seg.begin(); it != seg.end(); it++) {
0966     int st = ((*it).dt_DetId).station();
0967     int wl = abs(((*it).dt_DetId).wheel());
0968 
0969     MBPath[0][st][wl] = true;
0970     if (st == 4)
0971       continue;
0972 
0973     // for single-chamber segment in dt
0974     double ab = (((*it).sGlobalVec).x() * ((*it).sGlobalOrg).x()) + (((*it).sGlobalVec).y() * ((*it).sGlobalOrg).y());
0975     double al =
0976         sqrt((((*it).sGlobalOrg).x() * ((*it).sGlobalOrg).x()) + (((*it).sGlobalOrg).y() * ((*it).sGlobalOrg).y()));
0977     double bl =
0978         sqrt((((*it).sGlobalVec).x() * ((*it).sGlobalVec).x()) + (((*it).sGlobalVec).y() * ((*it).sGlobalVec).y()));
0979     double axb = (((*it).sGlobalOrg).x() * ((*it).sGlobalVec).y()) - (((*it).sGlobalOrg).y() * ((*it).sGlobalVec).x());
0980     double cc = (axb < 0.) ? 1.0 : -1.0;
0981     MB_phi[0][st][wl] = cc * acos(ab / (al * bl));
0982     MB_eta[0][st][wl] = fabs(((*it).sGlobalOrg).eta());
0983   }
0984 }
[ 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!