Back to home page

Project CMSSW displayed by LXR

 
 

    


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

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