Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:59:29 

0001 // -*- C++ -*-
0002 //
0003 // Package:    HTrackAssociator
0004 // Class:      HDetIdAssociator
0005 //
0006 /*
0007 
0008  Description: <one line class summary>
0009 
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Dmytro Kovalskyi
0015 // Modified for ECAL+HCAL by:  Michal Szleper
0016 //         Created:  Fri Apr 21 10:59:41 PDT 2006
0017 //
0018 //
0019 
0020 #include "Calibration/Tools/interface/DetIdAssociator.h"
0021 
0022 #include <memory>
0023 
0024 // surfaces is a vector of GlobalPoint representing outermost point on a cylinder
0025 std::vector<GlobalPoint> HDetIdAssociator::getTrajectory(const FreeTrajectoryState& ftsStart,
0026                                                          const std::vector<GlobalPoint>& surfaces) {
0027   check_setup();
0028   std::vector<GlobalPoint> trajectory;
0029   TrajectoryStateOnSurface tSOSDest;
0030   FreeTrajectoryState ftsCurrent = ftsStart;
0031 
0032   for (std::vector<GlobalPoint>::const_iterator surface_iter = surfaces.begin(); surface_iter != surfaces.end();
0033        surface_iter++) {
0034     // this stuff is some weird pointer, which destroy itself
0035     std::unique_ptr<Cylinder> cylinder =
0036         std::make_unique<Cylinder>(surface_iter->perp(), Surface::PositionType(0, 0, 0), Surface::RotationType());
0037     std::unique_ptr<Plane> forwardEndcap =
0038         std::make_unique<Plane>(Surface::PositionType(0, 0, surface_iter->z()), Surface::RotationType());
0039     std::unique_ptr<Plane> backwardEndcap =
0040         std::make_unique<Plane>(Surface::PositionType(0, 0, -surface_iter->z()), Surface::RotationType());
0041 
0042     LogTrace("StartingPoint") << "Propagate from "
0043                               << "\n"
0044                               << "\tx: " << ftsStart.position().x() << "\n"
0045                               << "\ty: " << ftsStart.position().y() << "\n"
0046                               << "\tz: " << ftsStart.position().z() << "\n"
0047                               << "\tmomentum eta: " << ftsStart.momentum().eta() << "\n"
0048                               << "\tmomentum phi: " << ftsStart.momentum().phi() << "\n"
0049                               << "\tmomentum: " << ftsStart.momentum().mag() << "\n";
0050 
0051     float tanTheta = ftsCurrent.momentum().perp() / ftsCurrent.momentum().z();
0052     float corner = surface_iter->perp() / surface_iter->z();
0053     /*
0054       std::cout<<"Propagate from "<< "\n"
0055         << "\tx: " << ftsCurrent.position().x()<< "\n"
0056         << "\ty: " << ftsCurrent.position().y()<< "\n"
0057         << "\tz: " << ftsCurrent.position().z()<< "\n"
0058         << "\tz: " << ftsCurrent.position().perp()<< "\n"
0059         << "\tz: " << tanTheta<<" "<< corner <<"\n"
0060         << "\tmomentum eta: " << ftsCurrent.momentum().eta()<< "\n"
0061         << "\tmomentum phi: " << ftsCurrent.momentum().phi()<< "\n"
0062         << "\tmomentum: " << ftsCurrent.momentum().mag()<<std::endl;
0063 */
0064     // First propage the track to the cylinder if |eta|<1, othewise to the encap
0065     // and correct depending on the result
0066     int ibar = 0;
0067     if (fabs(tanTheta) > corner) {
0068       tSOSDest = ivProp_->propagate(ftsCurrent, *cylinder);
0069       //                  std::cout<<" Propagate to cylinder "<<std::endl;
0070     } else if (tanTheta > 0.) {
0071       tSOSDest = ivProp_->propagate(ftsCurrent, *forwardEndcap);
0072       ibar = 1;
0073     } else {
0074       tSOSDest = ivProp_->propagate(ftsCurrent, *backwardEndcap);
0075       ibar = -1;
0076     }
0077 
0078     //       std::cout<<" Trajectory valid? "<<tSOSDest.isValid()<<" First propagation in "<<ibar<<std::endl;
0079 
0080     if (!tSOSDest.isValid()) {
0081       // barrel
0082       if (ibar == 0) {
0083         if (tanTheta < 0)
0084           tSOSDest = ivProp_->propagate(ftsCurrent, *forwardEndcap);
0085         if (tanTheta >= 0)
0086           tSOSDest = ivProp_->propagate(ftsCurrent, *backwardEndcap);
0087       } else {
0088         tSOSDest = ivProp_->propagate(ftsCurrent, *cylinder);
0089       }
0090     } else {
0091       // missed target
0092       if (abs(ibar) > 0) {
0093         if (tSOSDest.globalPosition().perp() > surface_iter->perp()) {
0094           tSOSDest = ivProp_->propagate(ftsCurrent, *cylinder);
0095         }
0096       } else {
0097         if (tanTheta < 0)
0098           tSOSDest = ivProp_->propagate(ftsCurrent, *forwardEndcap);
0099         if (tanTheta >= 0)
0100           tSOSDest = ivProp_->propagate(ftsCurrent, *backwardEndcap);
0101       }
0102     }
0103 
0104     // If missed the target, propagate to again
0105     //      if ((!tSOSDest.isValid()) && point.perp() > surface_iter->perp())
0106     //  {tSOSDest = ivProp_->propagate(ftsCurrent, *cylinder);std::cout<<" Propagate again 1 "<<std::endl;}
0107     //         std::cout<<" Track is ok after repropagation to cylinder or not? "<<tSOSDest.isValid()<<std::endl;
0108     //      if ((!tSOSDest.isValid()) && ftsStart.momentum().eta()>0. && fabs(ftsStart.momentum().eta())>1.)
0109     //  {tSOSDest = ivProp_->propagate(ftsStart, *forwardEndcap);std::cout<<" Propagate again 2 "<<std::endl;}
0110     //       std::cout<<" Track is ok after repropagation forward or not? "<<tSOSDest.isValid()<<std::endl;
0111     //      if ((!tSOSDest.isValid()) && ftsStart.momentum().eta()<0.&&fabs(ftsStart.momentum().eta())>1.)
0112     //  {tSOSDest = ivProp_->propagate(ftsStart, *backwardEndcap);std::cout<<" Propagate again 3 "<<std::endl;}
0113     //       std::cout<<" Track is after repropagation backward ok or not? "<<tSOSDest.isValid()<<std::endl;
0114 
0115     if (!tSOSDest.isValid())
0116       return trajectory;
0117 
0118     //      std::cout<<" Propagate reach something"<<std::endl;
0119     LogTrace("SuccessfullPropagation") << "Great, I reached something."
0120                                        << "\n"
0121                                        << "\tx: " << tSOSDest.freeState()->position().x() << "\n"
0122                                        << "\ty: " << tSOSDest.freeState()->position().y() << "\n"
0123                                        << "\tz: " << tSOSDest.freeState()->position().z() << "\n"
0124                                        << "\teta: " << tSOSDest.freeState()->position().eta() << "\n"
0125                                        << "\tphi: " << tSOSDest.freeState()->position().phi() << "\n";
0126 
0127     //      std::cout<<" The position of trajectory "<<tSOSDest.freeState()->position().perp()<<" "<<tSOSDest.freeState()->position().z()<<std::endl;
0128 
0129     GlobalPoint point = tSOSDest.freeState()->position();
0130     point = tSOSDest.freeState()->position();
0131     ftsCurrent = *tSOSDest.freeState();
0132     trajectory.push_back(point);
0133   }
0134   return trajectory;
0135 }
0136 
0137 //------------------------------------------------------------------------------
0138 std::set<DetId> HDetIdAssociator::getDetIdsCloseToAPoint(const GlobalPoint& direction, const int idR) {
0139   std::set<DetId> set;
0140   check_setup();
0141   if (!theMap_)
0142     buildMap();
0143   LogTrace("MatchPoint") << "point (eta,phi): " << direction.eta() << "," << direction.phi() << "\n";
0144   int ieta = iEta(direction);
0145   int iphi = iPhi(direction);
0146 
0147   LogTrace("MatchPoint") << "(ieta,iphi): " << ieta << "," << iphi << "\n";
0148 
0149   if (ieta >= 0 && ieta < nEta_ && iphi >= 0 && iphi < nPhi_) {
0150     set = (*theMap_)[ieta][iphi];
0151     if (idR > 0) {
0152       LogTrace("MatchPoint") << "Add neighbors (ieta,iphi): " << ieta << "," << iphi << "\n";
0153       //add neighbors
0154       int maxIEta = ieta + idR;
0155       int minIEta = ieta - idR;
0156       if (maxIEta >= nEta_)
0157         maxIEta = nEta_ - 1;
0158       if (minIEta < 0)
0159         minIEta = 0;
0160       int maxIPhi = iphi + idR;
0161       int minIPhi = iphi - idR;
0162       if (minIPhi < 0) {
0163         minIPhi += nPhi_;
0164         maxIPhi += nPhi_;
0165       }
0166       LogTrace("MatchPoint") << "\tieta (min,max): " << minIEta << "," << maxIEta << "\n";
0167       LogTrace("MatchPoint") << "\tiphi (min,max): " << minIPhi << "," << maxIPhi << "\n";
0168       for (int i = minIEta; i <= maxIEta; i++)
0169         for (int j = minIPhi; j <= maxIPhi; j++) {
0170           if (i == ieta && j == iphi)
0171             continue;  // already in the set
0172           set.insert((*theMap_)[i][j % nPhi_].begin(), (*theMap_)[i][j % nPhi_].end());
0173         }
0174     }
0175   }
0176   //     if(set.size() > 0) {
0177   //   if (ieta+idR<55 && ieta-idR>14 && set.size() != (2*idR+1)*(2*idR+1)){
0178   //       std::cout<<" RRRA: "<<set.size()<<" DetIds in region "<<ieta<<" "<<iphi<<std::endl;
0179   //       for( std::set<DetId>::const_iterator itr=set.begin(); itr!=set.end(); itr++) {
0180   //         GlobalPoint point = getPosition(*itr);
0181   //         std::cout << "DetId: " <<itr->rawId() <<" (eta,phi): " << point.eta() << "," << point.phi()<<" "<<iEta(point)<<" "<<iPhi(point)<<std::endl;
0182   //       }
0183   //     }
0184   //     else {
0185   //       std::cout <<" HDetIdAssociator::getDetIdsCloseToAPoint::There are strange days "<<std::endl;
0186   //     }
0187   return set;
0188 }
0189 
0190 //------------------------------------------------------------------------------
0191 int HDetIdAssociator::iEta(const GlobalPoint& point) {
0192   // unequal bin sizes for endcap, following HCAL geometry
0193   int iEta1 = int(point.eta() / etaBinSize_ + nEta_ / 2);
0194   if (point.eta() > 1.827 && point.eta() <= 1.830)
0195     return iEta1 - 1;
0196   else if (point.eta() > 1.914 && point.eta() <= 1.930)
0197     return iEta1 - 1;
0198   else if (point.eta() > 2.001 && point.eta() <= 2.043)
0199     return iEta1 - 1;
0200   else if (point.eta() > 2.088 && point.eta() <= 2.172)
0201     return iEta1 - 1;
0202   else if (point.eta() > 2.175 && point.eta() <= 2.262)
0203     return iEta1 - 1;
0204   else if (point.eta() > 2.262 && point.eta() <= 2.332)
0205     return iEta1 - 2;
0206   else if (point.eta() > 2.332 && point.eta() <= 2.349)
0207     return iEta1 - 1;
0208   else if (point.eta() > 2.349 && point.eta() <= 2.436)
0209     return iEta1 - 2;
0210   else if (point.eta() > 2.436 && point.eta() <= 2.500)
0211     return iEta1 - 3;
0212   else if (point.eta() > 2.500 && point.eta() <= 2.523)
0213     return iEta1 - 2;
0214   else if (point.eta() > 2.523 && point.eta() <= 2.610)
0215     return iEta1 - 3;
0216   else if (point.eta() > 2.610 && point.eta() <= 2.650)
0217     return iEta1 - 4;
0218   else if (point.eta() > 2.650 && point.eta() <= 2.697)
0219     return iEta1 - 3;
0220   else if (point.eta() > 2.697 && point.eta() <= 2.784)
0221     return iEta1 - 4;
0222   else if (point.eta() > 2.784 && point.eta() <= 2.868)
0223     return iEta1 - 5;
0224   else if (point.eta() > 2.868 && point.eta() <= 2.871)
0225     return iEta1 - 4;
0226   else if (point.eta() > 2.871 && point.eta() <= 2.958)
0227     return iEta1 - 5;
0228   else if (point.eta() > 2.958)
0229     return iEta1 - 6;
0230   else if (point.eta() < -1.827 && point.eta() >= -1.830)
0231     return iEta1 + 1;
0232   else if (point.eta() < -1.914 && point.eta() >= -1.930)
0233     return iEta1 + 1;
0234   else if (point.eta() < -2.001 && point.eta() >= -2.043)
0235     return iEta1 + 1;
0236   else if (point.eta() < -2.088 && point.eta() >= -2.172)
0237     return iEta1 + 1;
0238   else if (point.eta() < -2.175 && point.eta() >= -2.262)
0239     return iEta1 + 1;
0240   else if (point.eta() < -2.262 && point.eta() >= -2.332)
0241     return iEta1 + 2;
0242   else if (point.eta() < -2.332 && point.eta() >= -2.349)
0243     return iEta1 + 1;
0244   else if (point.eta() < -2.349 && point.eta() >= -2.436)
0245     return iEta1 + 2;
0246   else if (point.eta() < -2.436 && point.eta() >= -2.500)
0247     return iEta1 + 3;
0248   else if (point.eta() < -2.500 && point.eta() >= -2.523)
0249     return iEta1 + 2;
0250   else if (point.eta() < -2.523 && point.eta() >= -2.610)
0251     return iEta1 + 3;
0252   else if (point.eta() < -2.610 && point.eta() >= -2.650)
0253     return iEta1 + 4;
0254   else if (point.eta() < -2.650 && point.eta() >= -2.697)
0255     return iEta1 + 3;
0256   else if (point.eta() < -2.697 && point.eta() >= -2.784)
0257     return iEta1 + 4;
0258   else if (point.eta() < -2.784 && point.eta() >= -2.868)
0259     return iEta1 + 5;
0260   else if (point.eta() < -2.868 && point.eta() >= -2.871)
0261     return iEta1 + 4;
0262   else if (point.eta() < -2.871 && point.eta() >= -2.958)
0263     return iEta1 + 5;
0264   else if (point.eta() < -2.349)
0265     return iEta1 + 6;
0266   else
0267     return iEta1;
0268 }
0269 
0270 //------------------------------------------------------------------------------
0271 int HDetIdAssociator::iPhi(const GlobalPoint& point) {
0272   double pi = 4 * atan(1.);
0273   int iPhi1 = int((double(point.phi()) + pi) / (2 * pi) * nPhi_);
0274   return iPhi1;
0275 }
0276 
0277 //------------------------------------------------------------------------------
0278 void HDetIdAssociator::buildMap() {
0279   // modified version: take only detector central position
0280   check_setup();
0281   LogTrace("HDetIdAssociator") << "building map"
0282                                << "\n";
0283   if (theMap_)
0284     delete theMap_;
0285   theMap_ = new std::vector<std::vector<std::set<DetId> > >(nEta_, std::vector<std::set<DetId> >(nPhi_));
0286   int numberOfDetIdsOutsideEtaRange = 0;
0287   int numberOfDetIdsActive = 0;
0288   std::set<DetId> validIds = getASetOfValidDetIds();
0289   for (std::set<DetId>::const_iterator id_itr = validIds.begin(); id_itr != validIds.end(); id_itr++) {
0290     //      std::vector<GlobalPoint> points = getDetIdPoints(*id_itr);
0291     GlobalPoint point = getPosition(*id_itr);
0292     // reject fake DetIds (eta=0 - what are they anyway???)
0293     if (point.eta() == 0)
0294       continue;
0295 
0296     int ieta = iEta(point);
0297     int iphi = iPhi(point);
0298     int etaMax(-1);
0299     int etaMin(nEta_);
0300     int phiMax(-1);
0301     int phiMin(nPhi_);
0302     if (iphi >= nPhi_)
0303       iphi = iphi % nPhi_;
0304     assert(iphi >= 0);
0305     if (etaMin > ieta)
0306       etaMin = ieta;
0307     if (etaMax < ieta)
0308       etaMax = ieta;
0309     if (phiMin > iphi)
0310       phiMin = iphi;
0311     if (phiMax < iphi)
0312       phiMax = iphi;
0313     // for abs(eta)>1.8 one tower covers two phi segments
0314     if ((ieta > 54 || ieta < 15) && iphi % 2 == 0)
0315       phiMax++;
0316     if ((ieta > 54 || ieta < 15) && iphi % 2 == 1)
0317       phiMin--;
0318 
0319     if (etaMax < 0 || phiMax < 0 || etaMin >= nEta_ || phiMin >= nPhi_) {
0320       LogTrace("HDetIdAssociator") << "Out of range: DetId:" << id_itr->rawId() << "\n\teta (min,max): " << etaMin
0321                                    << "," << etaMax << "\n\tphi (min,max): " << phiMin << "," << phiMax
0322                                    << "\nTower id: " << id_itr->rawId() << "\n";
0323       numberOfDetIdsOutsideEtaRange++;
0324       continue;
0325     }
0326 
0327     if (phiMax - phiMin > phiMin + nPhi_ - phiMax) {
0328       phiMin += nPhi_;
0329       std::swap(phiMin, phiMax);
0330     }
0331     for (int ieta = etaMin; ieta <= etaMax; ieta++)
0332       for (int iphi = phiMin; iphi <= phiMax; iphi++)
0333         (*theMap_)[ieta][iphi % nPhi_].insert(*id_itr);
0334     numberOfDetIdsActive++;
0335   }
0336   LogTrace("HDetIdAssociator") << "Number of elements outside the allowed range ( |eta|>" << nEta_ / 2 * etaBinSize_
0337                                << "): " << numberOfDetIdsOutsideEtaRange << "\n";
0338   LogTrace("HDetIdAssociator") << "Number of active DetId's mapped: " << numberOfDetIdsActive << "\n";
0339 }
0340 
0341 //------------------------------------------------------------------------------
0342 std::set<DetId> HDetIdAssociator::getDetIdsInACone(const std::set<DetId>& inset,
0343                                                    const std::vector<GlobalPoint>& trajectory,
0344                                                    const double dR) {
0345   // modified version: if dR<0, returns 3x3 towers around the input one (Michal)
0346   check_setup();
0347   std::set<DetId> outset;
0348 
0349   if (dR >= 0) {
0350     for (std::set<DetId>::const_iterator id_iter = inset.begin(); id_iter != inset.end(); id_iter++)
0351       for (std::vector<GlobalPoint>::const_iterator point_iter = trajectory.begin(); point_iter != trajectory.end();
0352            point_iter++)
0353         if (nearElement(*point_iter, *id_iter, dR))
0354           outset.insert(*id_iter);
0355   } else {
0356     if (inset.size() != 1)
0357       return outset;
0358     std::set<DetId>::const_iterator id_inp = inset.begin();
0359     int ieta;
0360     int iphi;
0361     GlobalPoint point = getPosition(*id_inp);
0362     ieta = iEta(point);
0363     iphi = iPhi(point);
0364     for (int i = ieta - 1; i <= ieta + 1; i++) {
0365       for (int j = iphi - 1; j <= iphi + 1; j++) {
0366         //         if( i==ieta && j==iphi) continue;
0367         if (i < 0 || i >= nEta_)
0368           continue;
0369         int j2fill = j % nPhi_;
0370         if (j2fill < 0)
0371           j2fill += nPhi_;
0372         if ((*theMap_)[i][j2fill].empty())
0373           continue;
0374         outset.insert((*theMap_)[i][j2fill].begin(), (*theMap_)[i][j2fill].end());
0375       }
0376     }
0377   }
0378 
0379   //   if (outset.size() > 0) {
0380   //     std::cout<<" RRRA: DetIds in cone:"<<std::endl;
0381   //     for( std::set<DetId>::const_iterator itr=outset.begin(); itr!=outset.end(); itr++) {
0382   //       GlobalPoint point = getPosition(*itr);
0383   //       std::cout << "DetId: " <<itr->rawId() <<" (eta,phi): " << point.eta() << "," << point.phi()<<std::endl;
0384   //     }
0385   //   }
0386 
0387   return outset;
0388 }
0389 
0390 //------------------------------------------------------------------------------
0391 std::set<DetId> HDetIdAssociator::getCrossedDetIds(const std::set<DetId>& inset,
0392                                                    const std::vector<GlobalPoint>& trajectory) {
0393   check_setup();
0394   std::set<DetId> outset;
0395   for (std::set<DetId>::const_iterator id_iter = inset.begin(); id_iter != inset.end(); id_iter++)
0396     for (std::vector<GlobalPoint>::const_iterator point_iter = trajectory.begin(); point_iter != trajectory.end();
0397          point_iter++)
0398       if (insideElement(*point_iter, *id_iter))
0399         outset.insert(*id_iter);
0400   return outset;
0401 }
0402 
0403 //------------------------------------------------------------------------------
0404 std::set<DetId> HDetIdAssociator::getMaxEDetId(const std::set<DetId>& inset,
0405                                                edm::Handle<CaloTowerCollection> caloTowers) {
0406   // returns the most energetic tower in the NxN box (Michal)
0407   check_setup();
0408   std::set<DetId> outset;
0409   std::set<DetId>::const_iterator id_max = inset.begin();
0410   double Ehadmax = 0;
0411 
0412   for (std::set<DetId>::const_iterator id_iter = inset.begin(); id_iter != inset.end(); id_iter++) {
0413     DetId id(*id_iter);
0414     //     GlobalPoint point = getPosition(*id_iter);
0415     //     int ieta = iEta(point);
0416     //     int iphi = iPhi(point);
0417     CaloTowerCollection::const_iterator tower = (*caloTowers).find(id);
0418     if (tower != (*caloTowers).end() && tower->hadEnergy() > Ehadmax) {
0419       id_max = id_iter;
0420       Ehadmax = tower->hadEnergy();
0421     }
0422   }
0423 
0424   if (Ehadmax > 0)
0425     outset.insert(*id_max);
0426 
0427   //   if (outset.size() > 0) {
0428   //     std::cout<<" RRRA: Most energetic DetId:"<<std::endl;
0429   //     for( std::set<DetId>::const_iterator itr=outset.begin(); itr!=outset.end(); itr++) {
0430   //       GlobalPoint point = getPosition(*itr);
0431   //       std::cout << "DetId: " <<itr->rawId() <<" (eta,phi): " << point.eta() << "," << point.phi()<<std::endl;
0432   //     }
0433   //   }
0434 
0435   return outset;
0436 }
0437 
0438 //------------------------------------------------------------------------------
0439 std::set<DetId> HDetIdAssociator::getMaxEDetId(const std::set<DetId>& inset,
0440                                                edm::Handle<HBHERecHitCollection> recHits) {
0441   // returns the most energetic tower in the NxN box - from RecHits (Michal)
0442   check_setup();
0443   std::set<DetId> outset;
0444   std::set<DetId>::const_iterator id_max = inset.begin();
0445   double Ehadmax = 0;
0446 
0447   for (std::set<DetId>::const_iterator id_iter = inset.begin(); id_iter != inset.end(); id_iter++) {
0448     DetId id(*id_iter);
0449     //     GlobalPoint point = getPosition(*id_iter);
0450     //     int ieta = iEta(point);
0451     //     int iphi = iPhi(point);
0452     HBHERecHitCollection::const_iterator hit = (*recHits).find(id);
0453     if (hit != (*recHits).end() && hit->energy() > Ehadmax) {
0454       id_max = id_iter;
0455       Ehadmax = hit->energy();
0456     }
0457   }
0458 
0459   if (Ehadmax > 0)
0460     outset.insert(*id_max);
0461 
0462   //   if (outset.size() > 0) {
0463   //     std::cout<<" RRRA: Most energetic DetId:"<<std::endl;
0464   //     for( std::set<DetId>::const_iterator itr=outset.begin(); itr!=outset.end(); itr++) {
0465   //       GlobalPoint point = getPosition(*itr);
0466   //       std::cout << "DetId: " <<itr->rawId() <<" (eta,phi): " << point.eta() << "," << point.phi()<<std::endl;
0467   //     }
0468   //   }
0469 
0470   return outset;
0471 }
[ 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!