Reference Manual

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    HTrackAssociator
00004 // Class:      HDetIdAssociator
00005 // 
00006 /*
00007 
00008  Description: <one line class summary>
00009 
00010  Implementation:
00011      <Notes on implementation>
00012 */
00013 //
00014 // Original Author:  Dmytro Kovalskyi
00015 // Modified for ECAL+HCAL by:  Michal Szleper
00016 //         Created:  Fri Apr 21 10:59:41 PDT 2006
00017 // $Id: DetIdAssociator.cc,v 1.6 2009/04/08 12:34:29 argiro Exp $
00018 //
00019 //
00020 
00021 
00022 #include "Calibration/Tools/interface/DetIdAssociator.h"
00023 
00024 
00025 // surfaces is a vector of GlobalPoint representing outermost point on a cylinder
00026 std::vector<GlobalPoint> HDetIdAssociator::getTrajectory( const FreeTrajectoryState& ftsStart,
00027                                                     const std::vector<GlobalPoint>& surfaces)
00028 {
00029    check_setup();
00030    std::vector<GlobalPoint> trajectory;
00031    TrajectoryStateOnSurface tSOSDest;
00032    FreeTrajectoryState ftsCurrent = ftsStart;
00033 
00034    for(std::vector<GlobalPoint>::const_iterator surface_iter = surfaces.begin(); 
00035        surface_iter != surfaces.end(); surface_iter++) {
00036       // this stuff is some weird pointer, which destroy itself
00037       Cylinder *cylinder = new Cylinder(Surface::PositionType(0,0,0),
00038                                         Surface::RotationType(), 
00039                                         double (surface_iter->perp()) );
00040       Plane *forwardEndcap = new Plane(Surface::PositionType(0,0,surface_iter->z()),
00041                                        Surface::RotationType());
00042       Plane *backwardEndcap = new Plane(Surface::PositionType(0,0,-surface_iter->z()),
00043                                         Surface::RotationType());
00044 
00045       
00046       LogTrace("StartingPoint")<< "Propagate from "<< "\n"
00047         << "\tx: " << ftsStart.position().x()<< "\n"
00048         << "\ty: " << ftsStart.position().y()<< "\n"
00049         << "\tz: " << ftsStart.position().z()<< "\n"
00050         << "\tmomentum eta: " << ftsStart.momentum().eta()<< "\n"
00051         << "\tmomentum phi: " << ftsStart.momentum().phi()<< "\n"
00052         << "\tmomentum: " << ftsStart.momentum().mag()<< "\n";
00053      
00054          float tanTheta = ftsCurrent.momentum().perp()/ftsCurrent.momentum().z();
00055          float corner = surface_iter->perp()/surface_iter->z();
00056 /*
00057       std::cout<<"Propagate from "<< "\n"
00058         << "\tx: " << ftsCurrent.position().x()<< "\n"
00059         << "\ty: " << ftsCurrent.position().y()<< "\n"
00060         << "\tz: " << ftsCurrent.position().z()<< "\n"
00061         << "\tz: " << ftsCurrent.position().perp()<< "\n"
00062         << "\tz: " << tanTheta<<" "<< corner <<"\n"
00063         << "\tmomentum eta: " << ftsCurrent.momentum().eta()<< "\n"
00064         << "\tmomentum phi: " << ftsCurrent.momentum().phi()<< "\n"
00065         << "\tmomentum: " << ftsCurrent.momentum().mag()<<std::endl;
00066 */ 
00067       // First propage the track to the cylinder if |eta|<1, othewise to the encap
00068       // and correct depending on the result
00069       int ibar = 0;
00070       if (fabs(tanTheta) > corner)
00071         {
00072                    tSOSDest = ivProp_->propagate(ftsCurrent, *cylinder); 
00073  //                  std::cout<<" Propagate to cylinder "<<std::endl;
00074         }
00075       else if(tanTheta > 0.)
00076         {tSOSDest = ivProp_->propagate(ftsCurrent, *forwardEndcap); ibar=1; }
00077       else
00078         {tSOSDest = ivProp_->propagate(ftsCurrent, *backwardEndcap); ibar=-1; }
00079 
00080 //       std::cout<<" Trajectory valid? "<<tSOSDest.isValid()<<" First propagation in "<<ibar<<std::endl;
00081 
00082       if(! tSOSDest.isValid() )
00083       {
00084 // barrel
00085        if(ibar == 0){ 
00086            if (tanTheta < 0 ) tSOSDest = ivProp_->propagate( ftsCurrent,*forwardEndcap);
00087            if (tanTheta >= 0 ) tSOSDest = ivProp_->propagate( ftsCurrent,*backwardEndcap);
00088        }
00089          else
00090          {
00091                 tSOSDest = ivProp_->propagate(ftsCurrent, *cylinder);
00092          }
00093       } else
00094           {
00095 // missed target
00096             if(abs(ibar) > 0)
00097             {
00098               if(tSOSDest.globalPosition().perp() > surface_iter->perp())
00099               {
00100                 tSOSDest = ivProp_->propagate(ftsCurrent, *cylinder);
00101               }           
00102             }
00103               else
00104               {
00105                 if (tanTheta < 0 ) tSOSDest = ivProp_->propagate( ftsCurrent,*forwardEndcap);
00106                 if (tanTheta >= 0 ) tSOSDest = ivProp_->propagate( ftsCurrent,*backwardEndcap);
00107               }
00108           }
00109 
00110 
00111       // If missed the target, propagate to again
00112 //      if ((!tSOSDest.isValid()) && point.perp() > surface_iter->perp())
00113 //      {tSOSDest = ivProp_->propagate(ftsCurrent, *cylinder);std::cout<<" Propagate again 1 "<<std::endl;}
00114 //         std::cout<<" Track is ok after repropagation to cylinder or not? "<<tSOSDest.isValid()<<std::endl;
00115 //      if ((!tSOSDest.isValid()) && ftsStart.momentum().eta()>0. && fabs(ftsStart.momentum().eta())>1.) 
00116 //      {tSOSDest = ivProp_->propagate(ftsStart, *forwardEndcap);std::cout<<" Propagate again 2 "<<std::endl;}
00117 //       std::cout<<" Track is ok after repropagation forward or not? "<<tSOSDest.isValid()<<std::endl;
00118 //      if ((!tSOSDest.isValid()) && ftsStart.momentum().eta()<0.&&fabs(ftsStart.momentum().eta())>1.) 
00119 //      {tSOSDest = ivProp_->propagate(ftsStart, *backwardEndcap);std::cout<<" Propagate again 3 "<<std::endl;}
00120 //       std::cout<<" Track is after repropagation backward ok or not? "<<tSOSDest.isValid()<<std::endl; 
00121 
00122 
00123       if (! tSOSDest.isValid()) return trajectory;
00124       
00125 //      std::cout<<" Propagate reach something"<<std::endl; 
00126       LogTrace("SuccessfullPropagation") << "Great, I reached something." << "\n"
00127         << "\tx: " << tSOSDest.freeState()->position().x() << "\n"
00128         << "\ty: " << tSOSDest.freeState()->position().y() << "\n"
00129         << "\tz: " << tSOSDest.freeState()->position().z() << "\n"
00130         << "\teta: " << tSOSDest.freeState()->position().eta() << "\n"
00131         << "\tphi: " << tSOSDest.freeState()->position().phi() << "\n";
00132 
00133 //      std::cout<<" The position of trajectory "<<tSOSDest.freeState()->position().perp()<<" "<<tSOSDest.freeState()->position().z()<<std::endl;  
00134 
00135       GlobalPoint point = tSOSDest.freeState()->position(); 
00136       point = tSOSDest.freeState()->position();
00137       ftsCurrent = *tSOSDest.freeState();
00138       trajectory.push_back(point);
00139    }
00140    return trajectory;
00141 }
00142 
00143 //------------------------------------------------------------------------------
00144 std::set<DetId> HDetIdAssociator::getDetIdsCloseToAPoint(const GlobalPoint& direction,
00145                                                    const int idR)
00146 {
00147    std::set<DetId> set;
00148    check_setup();
00149    int nDets=0;
00150    if (! theMap_) buildMap();
00151    LogTrace("MatchPoint") << "point (eta,phi): " << direction.eta() << "," << direction.phi() << "\n";
00152    int ieta = iEta(direction);
00153    int iphi = iPhi(direction);
00154 
00155    LogTrace("MatchPoint") << "(ieta,iphi): " << ieta << "," << iphi << "\n";
00156    
00157    if (ieta>=0 && ieta<nEta_ && iphi>=0 && iphi<nPhi_){
00158       set = (*theMap_)[ieta][iphi];
00159       nDets++;
00160       if (idR>0){
00161           LogTrace("MatchPoint") << "Add neighbors (ieta,iphi): " << ieta << "," << iphi << "\n";
00162          //add neighbors
00163          int maxIEta = ieta+idR;
00164          int minIEta = ieta-idR;
00165          if(maxIEta>=nEta_) maxIEta = nEta_-1;
00166          if(minIEta<0) minIEta = 0;
00167          int maxIPhi = iphi+idR;
00168          int minIPhi = iphi-idR;
00169          if(minIPhi<0) {
00170             minIPhi+=nPhi_;
00171             maxIPhi+=nPhi_;
00172          }
00173          LogTrace("MatchPoint") << "\tieta (min,max): " << minIEta << "," << maxIEta<< "\n";
00174          LogTrace("MatchPoint") << "\tiphi (min,max): " << minIPhi << "," << maxIPhi<< "\n";
00175          for (int i=minIEta;i<=maxIEta;i++)
00176            for (int j=minIPhi;j<=maxIPhi;j++) {
00177               if( i==ieta && j==iphi) continue; // already in the set
00178               set.insert((*theMap_)[i][j%nPhi_].begin(),(*theMap_)[i][j%nPhi_].end());
00179               nDets++;
00180            }
00181       }
00182    }
00183 //     if(set.size() > 0) {
00184 //   if (ieta+idR<55 && ieta-idR>14 && set.size() != (2*idR+1)*(2*idR+1)){
00185 //       std::cout<<" RRRA: "<<set.size()<<" DetIds in region "<<ieta<<" "<<iphi<<std::endl;
00186 //       for( std::set<DetId>::const_iterator itr=set.begin(); itr!=set.end(); itr++) {
00187 //         GlobalPoint point = getPosition(*itr);
00188 //         std::cout << "DetId: " <<itr->rawId() <<" (eta,phi): " << point.eta() << "," << point.phi()<<" "<<iEta(point)<<" "<<iPhi(point)<<std::endl;
00189 //       }
00190 //     }
00191 //     else {
00192 //       std::cout <<" HDetIdAssociator::getDetIdsCloseToAPoint::There are strange days "<<std::endl; 
00193 //     }  
00194    return set;
00195 }
00196 
00197 //------------------------------------------------------------------------------
00198 int HDetIdAssociator::iEta (const GlobalPoint& point)
00199 {
00200 // unequal bin sizes for endcap, following HCAL geometry
00201    int iEta1 = int(point.eta()/etaBinSize_ + nEta_/2);
00202    if (point.eta()>1.827 && point.eta()<=1.830) return iEta1-1;
00203    else if (point.eta()>1.914 && point.eta()<=1.930) return iEta1-1;
00204    else if (point.eta()>2.001 && point.eta()<=2.043) return iEta1-1;
00205    else if (point.eta()>2.088 && point.eta()<=2.172) return iEta1-1;
00206    else if (point.eta()>2.175 && point.eta()<=2.262) return iEta1-1;
00207    else if (point.eta()>2.262 && point.eta()<=2.332) return iEta1-2;
00208    else if (point.eta()>2.332 && point.eta()<=2.349) return iEta1-1;
00209    else if (point.eta()>2.349 && point.eta()<=2.436) return iEta1-2;
00210    else if (point.eta()>2.436 && point.eta()<=2.500) return iEta1-3;
00211    else if (point.eta()>2.500 && point.eta()<=2.523) return iEta1-2;
00212    else if (point.eta()>2.523 && point.eta()<=2.610) return iEta1-3;
00213    else if (point.eta()>2.610 && point.eta()<=2.650) return iEta1-4;
00214    else if (point.eta()>2.650 && point.eta()<=2.697) return iEta1-3;
00215    else if (point.eta()>2.697 && point.eta()<=2.784) return iEta1-4;
00216    else if (point.eta()>2.784 && point.eta()<=2.868) return iEta1-5;
00217    else if (point.eta()>2.868 && point.eta()<=2.871) return iEta1-4;
00218    else if (point.eta()>2.871 && point.eta()<=2.958) return iEta1-5;
00219    else if (point.eta()>2.958) return iEta1-6;
00220    else if (point.eta()<-1.827 && point.eta()>=-1.830) return iEta1+1;
00221    else if (point.eta()<-1.914 && point.eta()>=-1.930) return iEta1+1;
00222    else if (point.eta()<-2.001 && point.eta()>=-2.043) return iEta1+1;
00223    else if (point.eta()<-2.088 && point.eta()>=-2.172) return iEta1+1;
00224    else if (point.eta()<-2.175 && point.eta()>=-2.262) return iEta1+1;
00225    else if (point.eta()<-2.262 && point.eta()>=-2.332) return iEta1+2;
00226    else if (point.eta()<-2.332 && point.eta()>=-2.349) return iEta1+1;
00227    else if (point.eta()<-2.349 && point.eta()>=-2.436) return iEta1+2;
00228    else if (point.eta()<-2.436 && point.eta()>=-2.500) return iEta1+3;
00229    else if (point.eta()<-2.500 && point.eta()>=-2.523) return iEta1+2;
00230    else if (point.eta()<-2.523 && point.eta()>=-2.610) return iEta1+3;
00231    else if (point.eta()<-2.610 && point.eta()>=-2.650) return iEta1+4;
00232    else if (point.eta()<-2.650 && point.eta()>=-2.697) return iEta1+3;
00233    else if (point.eta()<-2.697 && point.eta()>=-2.784) return iEta1+4;
00234    else if (point.eta()<-2.784 && point.eta()>=-2.868) return iEta1+5;
00235    else if (point.eta()<-2.868 && point.eta()>=-2.871) return iEta1+4;
00236    else if (point.eta()<-2.871 && point.eta()>=-2.958) return iEta1+5;
00237    else if (point.eta()<-2.349) return iEta1+6;
00238    else return iEta1;
00239 }
00240 
00241 //------------------------------------------------------------------------------
00242 int HDetIdAssociator::iPhi (const GlobalPoint& point)
00243 {
00244    double pi=4*atan(1.);
00245    int iPhi1 = int((double(point.phi())+pi)/(2*pi)*nPhi_);
00246    return iPhi1;
00247 }
00248 
00249 //------------------------------------------------------------------------------
00250 void HDetIdAssociator::buildMap()
00251 {
00252 // modified version: take only detector central position
00253    check_setup();
00254    LogTrace("HDetIdAssociator")<<"building map" << "\n";
00255    if(theMap_) delete theMap_;
00256    theMap_ = new std::vector<std::vector<std::set<DetId> > >(nEta_,std::vector<std::set<DetId> >(nPhi_));
00257    int numberOfDetIdsOutsideEtaRange = 0;
00258    int numberOfDetIdsActive = 0;
00259    std::set<DetId> validIds = getASetOfValidDetIds();
00260    for (std::set<DetId>::const_iterator id_itr = validIds.begin(); id_itr!=validIds.end(); id_itr++) {   
00261 //      std::vector<GlobalPoint> points = getDetIdPoints(*id_itr);
00262       GlobalPoint point = getPosition(*id_itr);
00263 // reject fake DetIds (eta=0 - what are they anyway???)
00264       if(point.eta()==0)continue;
00265 
00266       int ieta = iEta(point);
00267       int iphi = iPhi(point);
00268       int etaMax(-1);
00269       int etaMin(nEta_);
00270       int phiMax(-1);
00271       int phiMin(nPhi_);
00272       if ( iphi >= nPhi_ ) iphi = iphi % nPhi_;
00273       assert (iphi>=0);
00274       if ( etaMin > ieta) etaMin = ieta;
00275       if ( etaMax < ieta) etaMax = ieta;
00276       if ( phiMin > iphi) phiMin = iphi;
00277       if ( phiMax < iphi) phiMax = iphi;
00278 // for abs(eta)>1.8 one tower covers two phi segments
00279       if ((ieta>54||ieta<15) && iphi%2==0) phiMax++;
00280       if ((ieta>54||ieta<15) && iphi%2==1) phiMin--;
00281 
00282       if (etaMax<0||phiMax<0||etaMin>=nEta_||phiMin>=nPhi_) {
00283          LogTrace("HDetIdAssociator")<<"Out of range: DetId:" << id_itr->rawId() <<
00284            "\n\teta (min,max): " << etaMin << "," << etaMax <<
00285            "\n\tphi (min,max): " << phiMin << "," << phiMax <<
00286            "\nTower id: " << id_itr->rawId() << "\n";
00287          numberOfDetIdsOutsideEtaRange++;
00288          continue;
00289       }
00290           
00291       if (phiMax-phiMin > phiMin+nPhi_-phiMax){
00292          phiMin += nPhi_;
00293          std::swap(phiMin,phiMax);
00294       }
00295       for(int ieta = etaMin; ieta <= etaMax; ieta++)
00296         for(int iphi = phiMin; iphi <= phiMax; iphi++)
00297           (*theMap_)[ieta][iphi%nPhi_].insert(*id_itr);
00298       numberOfDetIdsActive++;
00299    }
00300    LogTrace("HDetIdAssociator") << "Number of elements outside the allowed range ( |eta|>"<<
00301      nEta_/2*etaBinSize_ << "): " << numberOfDetIdsOutsideEtaRange << "\n";
00302    LogTrace("HDetIdAssociator") << "Number of active DetId's mapped: " << 
00303      numberOfDetIdsActive << "\n";
00304 }
00305 
00306 //------------------------------------------------------------------------------
00307 std::set<DetId> HDetIdAssociator::getDetIdsInACone(const std::set<DetId>& inset, 
00308                                      const std::vector<GlobalPoint>& trajectory,
00309                                      const double dR)
00310 {
00311 // modified version: if dR<0, returns 3x3 towers around the input one (Michal)
00312    check_setup();
00313    std::set<DetId> outset;
00314   
00315    if(dR>=0) {
00316      for(std::set<DetId>::const_iterator id_iter = inset.begin(); id_iter != inset.end(); id_iter++)
00317        for(std::vector<GlobalPoint>::const_iterator point_iter = trajectory.begin(); point_iter != trajectory.end(); point_iter++)
00318          if (nearElement(*point_iter,*id_iter,dR)) outset.insert(*id_iter);
00319    }
00320    else {
00321      if (inset.size()!=1) return outset;
00322      std::set<DetId>::const_iterator id_inp = inset.begin();
00323      int ieta;
00324      int iphi;
00325      GlobalPoint point = getPosition(*id_inp);
00326      ieta = iEta(point);
00327      iphi = iPhi(point);
00328      for (int i=ieta-1;i<=ieta+1;i++) {
00329        for (int j=iphi-1;j<=iphi+1;j++) {
00330 //         if( i==ieta && j==iphi) continue;
00331          if( i<0 || i>=nEta_) continue;
00332          int j2fill = j%nPhi_;
00333          if(j2fill<0) j2fill+=nPhi_;
00334          if((*theMap_)[i][j2fill].size()==0)continue;
00335          outset.insert((*theMap_)[i][j2fill].begin(),(*theMap_)[i][j2fill].end());
00336        }
00337      }
00338    }
00339 
00340 //   if (outset.size() > 0) {
00341 //     std::cout<<" RRRA: DetIds in cone:"<<std::endl;
00342 //     for( std::set<DetId>::const_iterator itr=outset.begin(); itr!=outset.end(); itr++) {
00343 //       GlobalPoint point = getPosition(*itr);
00344 //       std::cout << "DetId: " <<itr->rawId() <<" (eta,phi): " << point.eta() << "," << point.phi()<<std::endl;
00345 //     }
00346 //   }
00347 
00348    return outset;
00349 }
00350 
00351 //------------------------------------------------------------------------------
00352 std::set<DetId> HDetIdAssociator::getCrossedDetIds(const std::set<DetId>& inset,
00353                                              const std::vector<GlobalPoint>& trajectory)
00354 {
00355    check_setup();
00356    std::set<DetId> outset;
00357    for(std::set<DetId>::const_iterator id_iter = inset.begin(); id_iter != inset.end(); id_iter++)
00358      for(std::vector<GlobalPoint>::const_iterator point_iter = trajectory.begin(); point_iter != trajectory.end(); point_iter++)
00359        if (insideElement(*point_iter, *id_iter))  outset.insert(*id_iter);
00360    return outset;
00361 }
00362 
00363 //------------------------------------------------------------------------------
00364 std::set<DetId> HDetIdAssociator::getMaxEDetId(const std::set<DetId>& inset,
00365                                            edm::Handle<CaloTowerCollection> caloTowers)
00366 {
00367 // returns the most energetic tower in the NxN box (Michal)
00368    check_setup();
00369    std::set<DetId> outset;
00370    std::set<DetId>::const_iterator id_max = inset.begin();
00371    double Ehadmax=0;
00372 
00373    for(std::set<DetId>::const_iterator id_iter = inset.begin(); id_iter != inset.end(); id_iter++) {
00374      DetId id(*id_iter);
00375 //     GlobalPoint point = getPosition(*id_iter);
00376 //     int ieta = iEta(point);
00377 //     int iphi = iPhi(point);
00378      CaloTowerCollection::const_iterator tower = (*caloTowers).find(id);
00379      if(tower != (*caloTowers).end() && tower->hadEnergy()>Ehadmax) {
00380        id_max = id_iter;
00381        Ehadmax = tower->hadEnergy();
00382      }
00383    }
00384 
00385    if (Ehadmax > 0) outset.insert(*id_max);
00386 
00387 //   if (outset.size() > 0) {
00388 //     std::cout<<" RRRA: Most energetic DetId:"<<std::endl;
00389 //     for( std::set<DetId>::const_iterator itr=outset.begin(); itr!=outset.end(); itr++) {
00390 //       GlobalPoint point = getPosition(*itr);
00391 //       std::cout << "DetId: " <<itr->rawId() <<" (eta,phi): " << point.eta() << "," << point.phi()<<std::endl;
00392 //     }
00393 //   }
00394 
00395    return outset;
00396 }
00397 
00398 //------------------------------------------------------------------------------
00399 std::set<DetId> HDetIdAssociator::getMaxEDetId(const std::set<DetId>& inset,
00400                                            edm::Handle<HBHERecHitCollection> recHits)
00401 {
00402 // returns the most energetic tower in the NxN box - from RecHits (Michal)
00403    check_setup();
00404    std::set<DetId> outset;
00405    std::set<DetId>::const_iterator id_max = inset.begin();
00406    double Ehadmax=0;
00407 
00408    for(std::set<DetId>::const_iterator id_iter = inset.begin(); id_iter != inset.end(); id_iter++) {
00409      DetId id(*id_iter);
00410 //     GlobalPoint point = getPosition(*id_iter);
00411 //     int ieta = iEta(point);
00412 //     int iphi = iPhi(point);
00413      HBHERecHitCollection::const_iterator hit = (*recHits).find(id);
00414      if(hit != (*recHits).end() && hit->energy()>Ehadmax) {
00415        id_max = id_iter;
00416        Ehadmax = hit->energy();
00417      }
00418    }
00419 
00420    if (Ehadmax > 0) outset.insert(*id_max);
00421 
00422 //   if (outset.size() > 0) {
00423 //     std::cout<<" RRRA: Most energetic DetId:"<<std::endl;
00424 //     for( std::set<DetId>::const_iterator itr=outset.begin(); itr!=outset.end(); itr++) {
00425 //       GlobalPoint point = getPosition(*itr);
00426 //       std::cout << "DetId: " <<itr->rawId() <<" (eta,phi): " << point.eta() << "," << point.phi()<<std::endl;
00427 //     }
00428 //   }
00429 
00430    return outset;
00431 }
/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/Calibration/Tools/src/DetIdAssociator.cc
