Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-01-21 00:19:05 

0001 /** \file LaserAlignment.cc
0002  *  LAS reconstruction module
0003  *
0004  *  $Date: 2013/01/07 20:26:37 $
0005  *  $Revision: 1.47 $
0006  *  \author Maarten Thomas
0007  *  \author Jan Olzem
0008  */
0009 
0010 #include "Alignment/LaserAlignment/plugins/LaserAlignment.h"
0011 #include "FWCore/Framework/interface/Run.h"
0012 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0013 #include "CondFormats/GeometryObjects/interface/PTrackerParameters.h"
0014 #include "CondFormats/GeometryObjects/interface/PTrackerAdditionalParametersPerDet.h"
0015 #include "Geometry/Records/interface/PTrackerParametersRcd.h"
0016 #include "Geometry/Records/interface/PTrackerAdditionalParametersPerDetRcd.h"
0017 
0018 ///
0019 ///
0020 ///
0021 LaserAlignment::LaserAlignment(edm::ParameterSet const& theConf)
0022     : topoToken_(esConsumes()),
0023       geomToken_(esConsumes()),
0024       geomDetToken_(esConsumes()),
0025       ptpToken_(esConsumes()),
0026       gprToken_(esConsumes()),
0027       stripPedestalsToken_(esConsumes()),
0028       theEvents(0),
0029       theDoPedestalSubtraction(theConf.getUntrackedParameter<bool>("SubtractPedestals", true)),
0030       theUseMinuitAlgorithm(theConf.getUntrackedParameter<bool>("RunMinuitAlignmentTubeAlgorithm", false)),
0031       theApplyBeamKinkCorrections(theConf.getUntrackedParameter<bool>("ApplyBeamKinkCorrections", true)),
0032       peakFinderThreshold(theConf.getUntrackedParameter<double>("PeakFinderThreshold", 10.)),
0033       enableJudgeZeroFilter(theConf.getUntrackedParameter<bool>("EnableJudgeZeroFilter", true)),
0034       judgeOverdriveThreshold(theConf.getUntrackedParameter<unsigned int>("JudgeOverdriveThreshold", 220)),
0035       updateFromInputGeometry(theConf.getUntrackedParameter<bool>("UpdateFromInputGeometry", false)),
0036       misalignedByRefGeometry(theConf.getUntrackedParameter<bool>("MisalignedByRefGeometry", false)),
0037       theStoreToDB(theConf.getUntrackedParameter<bool>("SaveToDbase", false)),
0038       theDigiProducersList(theConf.getParameter<std::vector<edm::ParameterSet> >("DigiProducersList")),
0039       theSaveHistograms(theConf.getUntrackedParameter<bool>("SaveHistograms", false)),
0040       theCompression(theConf.getUntrackedParameter<int>("ROOTFileCompression", 1)),
0041       theFileName(theConf.getUntrackedParameter<std::string>("ROOTFileName", "test.root")),
0042       theMaskTecModules(theConf.getUntrackedParameter<std::vector<unsigned int> >("MaskTECModules")),
0043       theMaskAtModules(theConf.getUntrackedParameter<std::vector<unsigned int> >("MaskATModules")),
0044       theSetNominalStrips(theConf.getUntrackedParameter<bool>("ForceFitterToNominalStrips", false)),
0045       theLasConstants(theConf.getUntrackedParameter<std::vector<edm::ParameterSet> >("LaserAlignmentConstants")),
0046       theFile(),
0047       theAlignableTracker(),
0048       theAlignRecordName("TrackerAlignmentRcd"),
0049       theErrorRecordName("TrackerAlignmentErrorExtendedRcd"),
0050       firstEvent_(true) {
0051   std::cout << std::endl;
0052   std::cout << "=============================================================="
0053             << "\n===         LaserAlignment module configuration            ==="
0054             << "\n"
0055             << "\n    Write histograms to file       = " << (theSaveHistograms ? "true" : "false")
0056             << "\n    Histogram file name            = " << theFileName
0057             << "\n    Histogram file compression     = " << theCompression
0058             << "\n    Subtract pedestals             = " << (theDoPedestalSubtraction ? "true" : "false")
0059             << "\n    Run Minuit AT algorithm        = " << (theUseMinuitAlgorithm ? "true" : "false")
0060             << "\n    Apply beam kink corrections    = " << (theApplyBeamKinkCorrections ? "true" : "false")
0061             << "\n    Peak Finder Threshold          = " << peakFinderThreshold
0062             << "\n    EnableJudgeZeroFilter          = " << (enableJudgeZeroFilter ? "true" : "false")
0063             << "\n    JudgeOverdriveThreshold        = " << judgeOverdriveThreshold
0064             << "\n    Update from input geometry     = " << (updateFromInputGeometry ? "true" : "false")
0065             << "\n    Misalignment from ref geometry = " << (misalignedByRefGeometry ? "true" : "false")
0066             << "\n    Number of TEC modules masked   = " << theMaskTecModules.size() << " (s. below list if > 0)"
0067             << "\n    Number of AT modules masked    = " << theMaskAtModules.size() << " (s. below list if > 0)"
0068             << "\n    Store to database              = " << (theStoreToDB ? "true" : "false")
0069             << "\n    ----------------------------------------------- ----------"
0070             << (theSetNominalStrips ? "\n    Set strips to nominal       =  true" : "\n")
0071             << "\n=============================================================" << std::endl;
0072 
0073   // tell about masked modules
0074   if (!theMaskTecModules.empty()) {
0075     std::cout << " ===============================================================================================\n"
0076               << std::flush;
0077     std::cout << " The following " << theMaskTecModules.size()
0078               << " TEC modules have been masked out and will not be considered by the TEC algorithm:\n " << std::flush;
0079     for (std::vector<unsigned int>::iterator moduleIt = theMaskTecModules.begin(); moduleIt != theMaskTecModules.end();
0080          ++moduleIt) {
0081       std::cout << *moduleIt << (moduleIt != --theMaskTecModules.end() ? ", " : "") << std::flush;
0082     }
0083     std::cout << std::endl << std::flush;
0084     std::cout << " ===============================================================================================\n\n"
0085               << std::flush;
0086   }
0087   if (!theMaskAtModules.empty()) {
0088     std::cout << " ===============================================================================================\n"
0089               << std::flush;
0090     std::cout << " The following " << theMaskAtModules.size()
0091               << " AT modules have been masked out and will not be considered by the AT algorithm:\n " << std::flush;
0092     for (std::vector<unsigned int>::iterator moduleIt = theMaskAtModules.begin(); moduleIt != theMaskAtModules.end();
0093          ++moduleIt) {
0094       std::cout << *moduleIt << (moduleIt != --theMaskAtModules.end() ? ", " : "") << std::flush;
0095     }
0096     std::cout << std::endl << std::flush;
0097     std::cout << " ===============================================================================================\n\n"
0098               << std::flush;
0099   }
0100 
0101   // alias for the Branches in the root files
0102   std::string alias(theConf.getParameter<std::string>("@module_label"));
0103 
0104   // declare the product to produce
0105   produces<TkLasBeamCollection, edm::Transition::EndRun>("tkLaserBeams").setBranchAlias(alias + "TkLasBeamCollection");
0106 
0107   // switch judge's zero filter depending on cfg
0108   judge.EnableZeroFilter(enableJudgeZeroFilter);
0109 
0110   // set the upper threshold for zero suppressed data
0111   judge.SetOverdriveThreshold(judgeOverdriveThreshold);
0112 }
0113 
0114 ///
0115 ///
0116 ///
0117 LaserAlignment::~LaserAlignment() {
0118   if (theSaveHistograms)
0119     theFile->Write();
0120   if (theFile) {
0121     delete theFile;
0122   }
0123   if (theAlignableTracker) {
0124     delete theAlignableTracker;
0125   }
0126 }
0127 
0128 ///
0129 ///
0130 ///
0131 void LaserAlignment::beginJob() {
0132   // write sumed histograms to file (if selected in cfg)
0133   if (theSaveHistograms) {
0134     // creating a new file
0135     theFile = new TFile(theFileName.c_str(), "RECREATE", "CMS ROOT file");
0136 
0137     // initialize the histograms
0138     if (theFile) {
0139       theFile->SetCompressionLevel(theCompression);
0140       singleModulesDir = theFile->mkdir("single modules");
0141     } else
0142       throw cms::Exception(" [LaserAlignment::beginJob]")
0143           << " ** ERROR: could not open file:" << theFileName.c_str() << " for writing." << std::endl;
0144   }
0145 
0146   // detector id maps (hard coded)
0147   fillDetectorId();
0148 
0149   // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0150   //   PROFILE, HISTOGRAM & FITFUNCTION INITIALIZATION
0151   // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0152 
0153   // object used to build various strings for names and labels
0154   std::stringstream nameBuilder;
0155 
0156   // loop variables for use with LASGlobalLoop object
0157   int det, ring, beam, disk, pos;
0158 
0159   // loop TEC modules
0160   det = 0;
0161   ring = 0;
0162   beam = 0;
0163   disk = 0;
0164   do {  // loop using LASGlobalLoop functionality
0165     // init the profiles
0166     pedestalProfiles.GetTECEntry(det, ring, beam, disk).SetAllValuesTo(0.);
0167     currentDataProfiles.GetTECEntry(det, ring, beam, disk).SetAllValuesTo(0.);
0168     collectedDataProfiles.GetTECEntry(det, ring, beam, disk).SetAllValuesTo(0.);
0169 
0170     // init the hit maps
0171     isAcceptedProfile.SetTECEntry(det, ring, beam, disk, 0);
0172     numberOfAcceptedProfiles.SetTECEntry(det, ring, beam, disk, 0);
0173 
0174     // create strings for histo names
0175     nameBuilder.clear();
0176     nameBuilder.str("");
0177     nameBuilder << "TEC";
0178     if (det == 0)
0179       nameBuilder << "+";
0180     else
0181       nameBuilder << "-";
0182     nameBuilder << "_Ring";
0183     if (ring == 0)
0184       nameBuilder << "4";
0185     else
0186       nameBuilder << "6";
0187     nameBuilder << "_Beam" << beam;
0188     nameBuilder << "_Disk" << disk;
0189     theProfileNames.SetTECEntry(det, ring, beam, disk, nameBuilder.str());
0190 
0191     // init the histograms
0192     if (theSaveHistograms) {
0193       nameBuilder << "_Histo";
0194       summedHistograms.SetTECEntry(
0195           det, ring, beam, disk, new TH1D(nameBuilder.str().c_str(), nameBuilder.str().c_str(), 512, 0, 512));
0196       summedHistograms.GetTECEntry(det, ring, beam, disk)->SetDirectory(singleModulesDir);
0197     }
0198 
0199   } while (moduleLoop.TECLoop(det, ring, beam, disk));
0200 
0201   // TIB & TOB section
0202   det = 2;
0203   beam = 0;
0204   pos = 0;
0205   do {  // loop using LASGlobalLoop functionality
0206     // init the profiles
0207     pedestalProfiles.GetTIBTOBEntry(det, beam, pos).SetAllValuesTo(0.);
0208     currentDataProfiles.GetTIBTOBEntry(det, beam, pos).SetAllValuesTo(0.);
0209     collectedDataProfiles.GetTIBTOBEntry(det, beam, pos).SetAllValuesTo(0.);
0210 
0211     // init the hit maps
0212     isAcceptedProfile.SetTIBTOBEntry(det, beam, pos, 0);
0213     numberOfAcceptedProfiles.SetTIBTOBEntry(det, beam, pos, 0);
0214 
0215     // create strings for histo names
0216     nameBuilder.clear();
0217     nameBuilder.str("");
0218     if (det == 2)
0219       nameBuilder << "TIB";
0220     else
0221       nameBuilder << "TOB";
0222     nameBuilder << "_Beam" << beam;
0223     nameBuilder << "_Zpos" << pos;
0224 
0225     theProfileNames.SetTIBTOBEntry(det, beam, pos, nameBuilder.str());
0226 
0227     // init the histograms
0228     if (theSaveHistograms) {
0229       nameBuilder << "_Histo";
0230       summedHistograms.SetTIBTOBEntry(
0231           det, beam, pos, new TH1D(nameBuilder.str().c_str(), nameBuilder.str().c_str(), 512, 0, 512));
0232       summedHistograms.GetTIBTOBEntry(det, beam, pos)->SetDirectory(singleModulesDir);
0233     }
0234 
0235   } while (moduleLoop.TIBTOBLoop(det, beam, pos));
0236 
0237   // TEC2TEC AT section
0238   det = 0;
0239   beam = 0;
0240   disk = 0;
0241   do {  // loop using LASGlobalLoop functionality
0242     // init the profiles
0243     pedestalProfiles.GetTEC2TECEntry(det, beam, disk).SetAllValuesTo(0.);
0244     currentDataProfiles.GetTEC2TECEntry(det, beam, disk).SetAllValuesTo(0.);
0245     collectedDataProfiles.GetTEC2TECEntry(det, beam, disk).SetAllValuesTo(0.);
0246 
0247     // init the hit maps
0248     isAcceptedProfile.SetTEC2TECEntry(det, beam, disk, 0);
0249     numberOfAcceptedProfiles.SetTEC2TECEntry(det, beam, disk, 0);
0250 
0251     // create strings for histo names
0252     nameBuilder.clear();
0253     nameBuilder.str("");
0254     nameBuilder << "TEC(AT)";
0255     if (det == 0)
0256       nameBuilder << "+";
0257     else
0258       nameBuilder << "-";
0259     nameBuilder << "_Beam" << beam;
0260     nameBuilder << "_Disk" << disk;
0261     theProfileNames.SetTEC2TECEntry(det, beam, disk, nameBuilder.str());
0262 
0263     // init the histograms
0264     if (theSaveHistograms) {
0265       nameBuilder << "_Histo";
0266       summedHistograms.SetTEC2TECEntry(
0267           det, beam, disk, new TH1D(nameBuilder.str().c_str(), nameBuilder.str().c_str(), 512, 0, 512));
0268       summedHistograms.GetTEC2TECEntry(det, beam, disk)->SetDirectory(singleModulesDir);
0269     }
0270 
0271   } while (moduleLoop.TEC2TECLoop(det, beam, disk));
0272 
0273   firstEvent_ = true;
0274 }
0275 
0276 ///
0277 ///
0278 ///
0279 void LaserAlignment::produce(edm::Event& theEvent, edm::EventSetup const& theSetup) {
0280   if (firstEvent_) {
0281     //Retrieve tracker topology from geometry
0282     const TrackerTopology* const tTopo = &theSetup.getData(topoToken_);
0283 
0284     // access the tracker
0285     gD = theSetup.getHandle(geomDetToken_);
0286     theTrackerGeometry = theSetup.getHandle(geomToken_);
0287 
0288     // access pedestals (from db..) if desired
0289     edm::ESHandle<SiStripPedestals> pedestalsHandle;
0290     if (theDoPedestalSubtraction) {
0291       pedestalsHandle = theSetup.getHandle(stripPedestalsToken_);
0292       fillPedestalProfiles(pedestalsHandle);
0293     }
0294 
0295     // global positions
0296     theGlobalPositionRcd = &theSetup.getData(gprToken_);
0297 
0298     // select the reference geometry
0299     if (!updateFromInputGeometry) {
0300       // the AlignableTracker object is initialized with the ideal geometry
0301       const GeometricDet* theGeometricDet = &theSetup.getData(geomDetToken_);
0302       const PTrackerParameters* ptp = &theSetup.getData(ptpToken_);
0303 
0304       TrackerGeomBuilderFromGeometricDet trackerBuilder;
0305       TrackerGeometry* theRefTracker = trackerBuilder.build(&*theGeometricDet, *ptp, tTopo);
0306 
0307       theAlignableTracker = new AlignableTracker(&(*theRefTracker), tTopo);
0308     } else {
0309       // the AlignableTracker object is initialized with the input geometry from DB
0310       theAlignableTracker = new AlignableTracker(&(*theTrackerGeometry), tTopo);
0311     }
0312 
0313     firstEvent_ = false;
0314   }
0315 
0316   LogDebug("LaserAlignment") << "==========================================================="
0317                              << "\n   Private analysis of event #" << theEvent.id().event() << " in run #"
0318                              << theEvent.id().run();
0319 
0320   // do the Tracker Statistics to retrieve the current profiles
0321   fillDataProfiles(theEvent, theSetup);
0322 
0323   // index variables for the LASGlobalLoop object
0324   int det, ring, beam, disk, pos;
0325 
0326   //
0327   // first pre-loop on selected entries to find out
0328   // whether the TEC or the AT beams have fired
0329   // (pedestal profiles are left empty if false in cfg)
0330   //
0331 
0332   // TEC+- (only ring 6)
0333   ring = 1;
0334   for (det = 0; det < 2; ++det) {
0335     for (beam = 0; beam < 8; ++beam) {
0336       for (disk = 0; disk < 9; ++disk) {
0337         if (judge.IsSignalIn(currentDataProfiles.GetTECEntry(det, ring, beam, disk) -
0338                                  pedestalProfiles.GetTECEntry(det, ring, beam, disk),
0339                              0)) {
0340           isAcceptedProfile.SetTECEntry(det, ring, beam, disk, 1);
0341         } else {  // assume no initialization
0342           isAcceptedProfile.SetTECEntry(det, ring, beam, disk, 0);
0343         }
0344       }
0345     }
0346   }
0347 
0348   // TIBTOB
0349   det = 2;
0350   beam = 0;
0351   pos = 0;
0352   do {
0353     // add current event's data and subtract pedestals
0354     if (judge.IsSignalIn(
0355             currentDataProfiles.GetTIBTOBEntry(det, beam, pos) - pedestalProfiles.GetTIBTOBEntry(det, beam, pos),
0356             getTIBTOBNominalBeamOffset(det, beam, pos))) {
0357       isAcceptedProfile.SetTIBTOBEntry(det, beam, pos, 1);
0358     } else {  // dto.
0359       isAcceptedProfile.SetTIBTOBEntry(det, beam, pos, 0);
0360     }
0361 
0362   } while (moduleLoop.TIBTOBLoop(det, beam, pos));
0363 
0364   // now come the beam finders
0365   bool isTECMode = isTECBeam();
0366   //  LogDebug( " [LaserAlignment::produce]" ) << "LaserAlignment::isTECBeam declares this event " << ( isTECMode ? "" : "NOT " ) << "a TEC event." << std::endl;
0367   std::cout << " [LaserAlignment::produce] -- LaserAlignment::isTECBeam declares this event "
0368             << (isTECMode ? "" : "NOT ") << "a TEC event." << std::endl;
0369 
0370   bool isATMode = isATBeam();
0371   //  LogDebug( " [LaserAlignment::produce]" ) << "LaserAlignment::isATBeam declares this event "  << ( isATMode ? "" : "NOT " )  << "an AT event." << std::endl;
0372   std::cout << " [LaserAlignment::produce] -- LaserAlignment::isATBeam declares this event " << (isATMode ? "" : "NOT ")
0373             << "an AT event." << std::endl;
0374 
0375   //
0376   // now pass the pedestal subtracted profiles to the judge
0377   // if they're accepted, add them on the collectedDataProfiles
0378   // (pedestal profiles are left empty if false in cfg)
0379   //
0380 
0381   // loop TEC+- modules
0382   det = 0;
0383   ring = 0;
0384   beam = 0;
0385   disk = 0;
0386   do {
0387     LogDebug("[LaserAlignment::produce]")
0388         << "Profile is: " << theProfileNames.GetTECEntry(det, ring, beam, disk) << "." << std::endl;
0389 
0390     // this now depends on the AT/TEC mode, is this a doubly hit module? -> look for it in vector<int> tecDoubleHitDetId
0391     // (ring == 0 not necessary but makes it a little faster)
0392     if (ring == 0 &&
0393         find(tecDoubleHitDetId.begin(), tecDoubleHitDetId.end(), detectorId.GetTECEntry(det, ring, beam, disk)) !=
0394             tecDoubleHitDetId.end()) {
0395       if (isTECMode) {  // add profile to TEC collection
0396         // add current event's data and subtract pedestals
0397         if (judge.JudgeProfile(currentDataProfiles.GetTECEntry(det, ring, beam, disk) -
0398                                    pedestalProfiles.GetTECEntry(det, ring, beam, disk),
0399                                0)) {
0400           collectedDataProfiles.GetTECEntry(det, ring, beam, disk) +=
0401               currentDataProfiles.GetTECEntry(det, ring, beam, disk) -
0402               pedestalProfiles.GetTECEntry(det, ring, beam, disk);
0403           numberOfAcceptedProfiles.GetTECEntry(det, ring, beam, disk)++;
0404         }
0405       }
0406     }
0407 
0408     else {  // not a doubly hit module, don't care about the mode
0409       // add current event's data and subtract pedestals
0410       if (judge.JudgeProfile(currentDataProfiles.GetTECEntry(det, ring, beam, disk) -
0411                                  pedestalProfiles.GetTECEntry(det, ring, beam, disk),
0412                              0)) {
0413         collectedDataProfiles.GetTECEntry(det, ring, beam, disk) +=
0414             currentDataProfiles.GetTECEntry(det, ring, beam, disk) -
0415             pedestalProfiles.GetTECEntry(det, ring, beam, disk);
0416         numberOfAcceptedProfiles.GetTECEntry(det, ring, beam, disk)++;
0417       }
0418     }
0419 
0420   } while (moduleLoop.TECLoop(det, ring, beam, disk));
0421 
0422   // loop TIB/TOB modules
0423   det = 2;
0424   beam = 0;
0425   pos = 0;
0426   do {
0427     LogDebug("[LaserAlignment::produce]")
0428         << "Profile is: " << theProfileNames.GetTIBTOBEntry(det, beam, pos) << "." << std::endl;
0429 
0430     // add current event's data and subtract pedestals
0431     if (judge.JudgeProfile(
0432             currentDataProfiles.GetTIBTOBEntry(det, beam, pos) - pedestalProfiles.GetTIBTOBEntry(det, beam, pos),
0433             getTIBTOBNominalBeamOffset(det, beam, pos))) {
0434       collectedDataProfiles.GetTIBTOBEntry(det, beam, pos) +=
0435           currentDataProfiles.GetTIBTOBEntry(det, beam, pos) - pedestalProfiles.GetTIBTOBEntry(det, beam, pos);
0436       numberOfAcceptedProfiles.GetTIBTOBEntry(det, beam, pos)++;
0437     }
0438 
0439   } while (moduleLoop.TIBTOBLoop(det, beam, pos));
0440 
0441   // loop TEC2TEC modules
0442   det = 0;
0443   beam = 0;
0444   disk = 0;
0445   do {
0446     LogDebug("[LaserAlignment::produce]")
0447         << "Profile is: " << theProfileNames.GetTEC2TECEntry(det, beam, disk) << "." << std::endl;
0448 
0449     // this again depends on the AT/TEC mode, is this a doubly hit module?
0450     // (ring == 0 not necessary but makes it a little faster)
0451     if (ring == 0 &&
0452         find(tecDoubleHitDetId.begin(), tecDoubleHitDetId.end(), detectorId.GetTECEntry(det, ring, beam, disk)) !=
0453             tecDoubleHitDetId.end()) {
0454       if (isATMode) {  // add profile to TEC2TEC collection
0455         // add current event's data and subtract pedestals
0456         if (judge.JudgeProfile(currentDataProfiles.GetTEC2TECEntry(det, beam, disk) -
0457                                    pedestalProfiles.GetTEC2TECEntry(det, beam, disk),
0458                                0)) {
0459           collectedDataProfiles.GetTEC2TECEntry(det, beam, disk) +=
0460               currentDataProfiles.GetTEC2TECEntry(det, beam, disk) - pedestalProfiles.GetTEC2TECEntry(det, beam, disk);
0461           numberOfAcceptedProfiles.GetTEC2TECEntry(det, beam, disk)++;
0462         }
0463       }
0464 
0465     }
0466 
0467     else {  // not a doubly hit module, don't care about the mode
0468       // add current event's data and subtract pedestals
0469       if (judge.JudgeProfile(
0470               currentDataProfiles.GetTEC2TECEntry(det, beam, disk) - pedestalProfiles.GetTEC2TECEntry(det, beam, disk),
0471               0)) {
0472         collectedDataProfiles.GetTEC2TECEntry(det, beam, disk) +=
0473             currentDataProfiles.GetTEC2TECEntry(det, beam, disk) - pedestalProfiles.GetTEC2TECEntry(det, beam, disk);
0474         numberOfAcceptedProfiles.GetTEC2TECEntry(det, beam, disk)++;
0475       }
0476     }
0477 
0478   } while (moduleLoop.TEC2TECLoop(det, beam, disk));
0479 
0480   // total event number counter
0481   theEvents++;
0482 }
0483 
0484 ///
0485 ///
0486 ///
0487 void LaserAlignment::endRunProduce(edm::Run& theRun, const edm::EventSetup& theSetup) {
0488   std::cout << " [LaserAlignment::endRun] -- Total number of events processed: " << theEvents << std::endl;
0489 
0490   // for debugging only..
0491   DumpHitmaps(numberOfAcceptedProfiles);
0492 
0493   // index variables for the LASGlobalLoop objects
0494   int det, ring, beam, disk, pos;
0495 
0496   // measured positions container for the algorithms
0497   LASGlobalData<LASCoordinateSet> measuredCoordinates;
0498 
0499   // fitted peak positions in units of strips (pair for value,error)
0500   LASGlobalData<std::pair<float, float> > measuredStripPositions;
0501 
0502   // the peak finder, a pair (pos/posErr in units of strips) for its results, and the success confirmation
0503   LASPeakFinder peakFinder;
0504   peakFinder.SetAmplitudeThreshold(peakFinderThreshold);
0505   std::pair<double, double> peakFinderResults;
0506   bool isGoodFit;
0507 
0508   // tracker geom. object for calculating the global beam positions
0509   const TrackerGeometry& theTracker(*theTrackerGeometry);
0510 
0511   // fill LASGlobalData<LASCoordinateSet> nominalCoordinates
0512   CalculateNominalCoordinates();
0513 
0514   // for determining the phi errors
0515   //    ErrorFrameTransformer errorTransformer; // later...
0516 
0517   // do the fits for TEC+- internal
0518   det = 0;
0519   ring = 0;
0520   beam = 0;
0521   disk = 0;
0522   do {
0523     // do the fit
0524     isGoodFit = peakFinder.FindPeakIn(collectedDataProfiles.GetTECEntry(det, ring, beam, disk),
0525                                       peakFinderResults,
0526                                       summedHistograms.GetTECEntry(det, ring, beam, disk),
0527                                       0);  // offset is 0 for TEC
0528 
0529     // now we have the measured positions in units of strips.
0530     if (!isGoodFit)
0531       std::cout << " [LaserAlignment::endRun] ** WARNING: Fit failed for TEC det: " << det << ", ring: " << ring
0532                 << ", beam: " << beam << ", disk: " << disk << " (id: " << detectorId.GetTECEntry(det, ring, beam, disk)
0533                 << ")." << std::endl;
0534 
0535     // <- here we will later implement the kink corrections
0536 
0537     // access the tracker geometry for this module
0538     const DetId theDetId(detectorId.GetTECEntry(det, ring, beam, disk));
0539     const StripGeomDetUnit* const theStripDet = dynamic_cast<const StripGeomDetUnit*>(theTracker.idToDet(theDetId));
0540 
0541     if (theStripDet) {
0542       // first, set the measured coordinates to their nominal values
0543       measuredCoordinates.SetTECEntry(det, ring, beam, disk, nominalCoordinates.GetTECEntry(det, ring, beam, disk));
0544 
0545       if (isGoodFit) {  // convert strip position to global phi and replace the nominal phi value/error
0546 
0547         measuredStripPositions.GetTECEntry(det, ring, beam, disk) = peakFinderResults;
0548         const float positionInStrips =
0549             theSetNominalStrips
0550                 ? 256.
0551                 : peakFinderResults.first;  // implementation of "ForceFitterToNominalStrips" config parameter
0552         const GlobalPoint& globalPoint =
0553             theStripDet->surface().toGlobal(theStripDet->specificTopology().localPosition(positionInStrips));
0554         measuredCoordinates.GetTECEntry(det, ring, beam, disk).SetPhi(ConvertAngle(globalPoint.barePhi()));
0555 
0556         // const GlobalError& globalError = errorTransformer.transform( theStripDet->specificTopology().localError( peakFinderResults.first, pow( peakFinderResults.second, 2 ) ), theStripDet->surface() );
0557         // measuredCoordinates.GetTECEntry( det, ring, beam, disk ).SetPhiError( globalError.phierr( globalPoint ) );
0558         measuredCoordinates.GetTECEntry(det, ring, beam, disk).SetPhiError(0.00046);  // PRELIMINARY ESTIMATE
0559 
0560       } else {  // keep nominal position (middle-of-module) but set a giant phi error so that the module can be ignored by the alignment algorithm
0561         measuredStripPositions.GetTECEntry(det, ring, beam, disk) = std::pair<float, float>(256., 1000.);
0562         const GlobalPoint& globalPoint =
0563             theStripDet->surface().toGlobal(theStripDet->specificTopology().localPosition(256.));
0564         measuredCoordinates.GetTECEntry(det, ring, beam, disk).SetPhi(ConvertAngle(globalPoint.barePhi()));
0565         measuredCoordinates.GetTECEntry(det, ring, beam, disk).SetPhiError(1000.);
0566       }
0567     }
0568 
0569   } while (moduleLoop.TECLoop(det, ring, beam, disk));
0570 
0571   // do the fits for TIB/TOB
0572   det = 2;
0573   beam = 0;
0574   pos = 0;
0575   do {
0576     // do the fit
0577     isGoodFit = peakFinder.FindPeakIn(collectedDataProfiles.GetTIBTOBEntry(det, beam, pos),
0578                                       peakFinderResults,
0579                                       summedHistograms.GetTIBTOBEntry(det, beam, pos),
0580                                       getTIBTOBNominalBeamOffset(det, beam, pos));
0581 
0582     // now we have the measured positions in units of strips.
0583     if (!isGoodFit)
0584       std::cout << " [LaserAlignment::endJob] ** WARNING: Fit failed for TIB/TOB det: " << det << ", beam: " << beam
0585                 << ", pos: " << pos << " (id: " << detectorId.GetTIBTOBEntry(det, beam, pos) << ")." << std::endl;
0586 
0587     // <- here we will later implement the kink corrections
0588 
0589     // access the tracker geometry for this module
0590     const DetId theDetId(detectorId.GetTIBTOBEntry(det, beam, pos));
0591     const StripGeomDetUnit* const theStripDet = dynamic_cast<const StripGeomDetUnit*>(theTracker.idToDet(theDetId));
0592 
0593     if (theStripDet) {
0594       // first, set the measured coordinates to their nominal values
0595       measuredCoordinates.SetTIBTOBEntry(det, beam, pos, nominalCoordinates.GetTIBTOBEntry(det, beam, pos));
0596 
0597       if (isGoodFit) {  // convert strip position to global phi and replace the nominal phi value/error
0598         measuredStripPositions.GetTIBTOBEntry(det, beam, pos) = peakFinderResults;
0599         const float positionInStrips =
0600             theSetNominalStrips
0601                 ? 256. + getTIBTOBNominalBeamOffset(det, beam, pos)
0602                 : peakFinderResults.first;  // implementation of "ForceFitterToNominalStrips" config parameter
0603         const GlobalPoint& globalPoint =
0604             theStripDet->surface().toGlobal(theStripDet->specificTopology().localPosition(positionInStrips));
0605         measuredCoordinates.GetTIBTOBEntry(det, beam, pos).SetPhi(ConvertAngle(globalPoint.barePhi()));
0606         measuredCoordinates.GetTIBTOBEntry(det, beam, pos).SetPhiError(0.00028);  // PRELIMINARY ESTIMATE
0607       } else {  // keep nominal position but set a giant phi error so that the module can be ignored by the alignment algorithm
0608         measuredStripPositions.GetTIBTOBEntry(det, beam, pos) =
0609             std::pair<float, float>(256. + getTIBTOBNominalBeamOffset(det, beam, pos), 1000.);
0610         const GlobalPoint& globalPoint = theStripDet->surface().toGlobal(
0611             theStripDet->specificTopology().localPosition(256. + getTIBTOBNominalBeamOffset(det, beam, pos)));
0612         measuredCoordinates.GetTIBTOBEntry(det, beam, pos).SetPhi(ConvertAngle(globalPoint.barePhi()));
0613         measuredCoordinates.GetTIBTOBEntry(det, beam, pos).SetPhiError(1000.);
0614       }
0615     }
0616 
0617   } while (moduleLoop.TIBTOBLoop(det, beam, pos));
0618 
0619   // do the fits for TEC AT
0620   det = 0;
0621   beam = 0;
0622   disk = 0;
0623   do {
0624     // do the fit
0625     isGoodFit = peakFinder.FindPeakIn(collectedDataProfiles.GetTEC2TECEntry(det, beam, disk),
0626                                       peakFinderResults,
0627                                       summedHistograms.GetTEC2TECEntry(det, beam, disk),
0628                                       getTEC2TECNominalBeamOffset(det, beam, disk));
0629     // now we have the positions in units of strips.
0630     if (!isGoodFit)
0631       std::cout << " [LaserAlignment::endRun] ** WARNING: Fit failed for TEC2TEC det: " << det << ", beam: " << beam
0632                 << ", disk: " << disk << " (id: " << detectorId.GetTEC2TECEntry(det, beam, disk) << ")." << std::endl;
0633 
0634     // <- here we will later implement the kink corrections
0635 
0636     // access the tracker geometry for this module
0637     const DetId theDetId(detectorId.GetTEC2TECEntry(det, beam, disk));
0638     const StripGeomDetUnit* const theStripDet = dynamic_cast<const StripGeomDetUnit*>(theTracker.idToDet(theDetId));
0639 
0640     if (theStripDet) {
0641       // first, set the measured coordinates to their nominal values
0642       measuredCoordinates.SetTEC2TECEntry(det, beam, disk, nominalCoordinates.GetTEC2TECEntry(det, beam, disk));
0643 
0644       if (isGoodFit) {  // convert strip position to global phi and replace the nominal phi value/error
0645         measuredStripPositions.GetTEC2TECEntry(det, beam, disk) = peakFinderResults;
0646         const float positionInStrips =
0647             theSetNominalStrips
0648                 ? 256. + getTEC2TECNominalBeamOffset(det, beam, disk)
0649                 : peakFinderResults.first;  // implementation of "ForceFitterToNominalStrips" config parameter
0650         const GlobalPoint& globalPoint =
0651             theStripDet->surface().toGlobal(theStripDet->specificTopology().localPosition(positionInStrips));
0652         measuredCoordinates.GetTEC2TECEntry(det, beam, disk).SetPhi(ConvertAngle(globalPoint.barePhi()));
0653         measuredCoordinates.GetTEC2TECEntry(det, beam, disk).SetPhiError(0.00047);  // PRELIMINARY ESTIMATE
0654       } else {  // keep nominal position but set a giant phi error so that the module can be ignored by the alignment algorithm
0655         measuredStripPositions.GetTEC2TECEntry(det, beam, disk) =
0656             std::pair<float, float>(256. + getTEC2TECNominalBeamOffset(det, beam, disk), 1000.);
0657         const GlobalPoint& globalPoint = theStripDet->surface().toGlobal(
0658             theStripDet->specificTopology().localPosition(256. + getTEC2TECNominalBeamOffset(det, beam, disk)));
0659         measuredCoordinates.GetTEC2TECEntry(det, beam, disk).SetPhi(ConvertAngle(globalPoint.barePhi()));
0660         measuredCoordinates.GetTEC2TECEntry(det, beam, disk).SetPhiError(1000.);
0661       }
0662     }
0663 
0664   } while (moduleLoop.TEC2TECLoop(det, beam, disk));
0665 
0666   // see what we got (for debugging)
0667   //  DumpStripFileSet( measuredStripPositions );
0668   //  DumpPosFileSet( measuredCoordinates );
0669 
0670   // CALCULATE PARAMETERS AND UPDATE DB OBJECT
0671   // for beam kink corrections, reconstructing the geometry and updating the db object
0672   LASGeometryUpdater geometryUpdater(nominalCoordinates, theLasConstants);
0673 
0674   // apply all beam corrections
0675   if (theApplyBeamKinkCorrections)
0676     geometryUpdater.ApplyBeamKinkCorrections(measuredCoordinates);
0677 
0678   // if we start with input geometry instead of IDEAL,
0679   // reverse the adjustments in the AlignableTracker object
0680   if (updateFromInputGeometry)
0681     geometryUpdater.SetReverseDirection(true);
0682 
0683   // if we have "virtual" misalignment which is introduced via the reference geometry,
0684   // tell the LASGeometryUpdater to reverse x & y adjustments
0685   if (misalignedByRefGeometry)
0686     geometryUpdater.SetMisalignmentFromRefGeometry(true);
0687 
0688   // run the endcap algorithm
0689   LASEndcapAlgorithm endcapAlgorithm;
0690   LASEndcapAlignmentParameterSet endcapParameters;
0691 
0692   // this basically sets all the endcap modules to be masked
0693   // to their nominal positions (since endcapParameters is overall zero)
0694   if (!theMaskTecModules.empty()) {
0695     ApplyEndcapMaskingCorrections(measuredCoordinates, nominalCoordinates, endcapParameters);
0696   }
0697 
0698   // run the algorithm
0699   endcapParameters = endcapAlgorithm.CalculateParameters(measuredCoordinates, nominalCoordinates);
0700 
0701   //
0702   // loop to mask out events
0703   // DESCRIPTION:
0704   //
0705 
0706   // do this only if there are modules to be masked..
0707   if (!theMaskTecModules.empty()) {
0708     const unsigned int nIterations = 30;
0709     for (unsigned int iteration = 0; iteration < nIterations; ++iteration) {
0710       // set the endcap modules to be masked to their positions
0711       // according to the reconstructed parameters
0712       ApplyEndcapMaskingCorrections(measuredCoordinates, nominalCoordinates, endcapParameters);
0713 
0714       // modifications applied, so re-run the algorithm
0715       endcapParameters = endcapAlgorithm.CalculateParameters(measuredCoordinates, nominalCoordinates);
0716     }
0717   }
0718 
0719   // these are now final, so:
0720   endcapParameters.Print();
0721 
0722   // do a pre-alignment of the endcaps (TEC2TEC only)
0723   // so that the alignment tube algorithms finds orderly disks
0724   geometryUpdater.EndcapUpdate(endcapParameters, measuredCoordinates);
0725 
0726   // the alignment tube algorithms, choose from config
0727   LASBarrelAlignmentParameterSet alignmentTubeParameters;
0728   // the MINUIT-BASED alignment tube algorithm
0729   LASBarrelAlgorithm barrelAlgorithm;
0730   // the ANALYTICAL alignment tube algorithm
0731   LASAlignmentTubeAlgorithm alignmentTubeAlgorithm;
0732 
0733   // this basically sets all the modules to be masked
0734   // to their nominal positions (since alignmentTubeParameters is overall zero)
0735   if (!theMaskAtModules.empty()) {
0736     ApplyATMaskingCorrections(measuredCoordinates, nominalCoordinates, alignmentTubeParameters);
0737   }
0738 
0739   if (theUseMinuitAlgorithm) {
0740     // run the MINUIT-BASED alignment tube algorithm
0741     alignmentTubeParameters = barrelAlgorithm.CalculateParameters(measuredCoordinates, nominalCoordinates);
0742   } else {
0743     // the ANALYTICAL alignment tube algorithm
0744     alignmentTubeParameters = alignmentTubeAlgorithm.CalculateParameters(measuredCoordinates, nominalCoordinates);
0745   }
0746 
0747   //
0748   // loop to mask out events
0749   // DESCRIPTION:
0750   //
0751 
0752   // do this only if there are modules to be masked..
0753   if (!theMaskAtModules.empty()) {
0754     const unsigned int nIterations = 30;
0755     for (unsigned int iteration = 0; iteration < nIterations; ++iteration) {
0756       // set the AT modules to be masked to their positions
0757       // according to the reconstructed parameters
0758       ApplyATMaskingCorrections(measuredCoordinates, nominalCoordinates, alignmentTubeParameters);
0759 
0760       // modifications applied, so re-run the algorithm
0761       if (theUseMinuitAlgorithm) {
0762         alignmentTubeParameters = barrelAlgorithm.CalculateParameters(measuredCoordinates, nominalCoordinates);
0763       } else {
0764         alignmentTubeParameters = alignmentTubeAlgorithm.CalculateParameters(measuredCoordinates, nominalCoordinates);
0765       }
0766     }
0767   }
0768 
0769   // these are now final, so:
0770   alignmentTubeParameters.Print();
0771 
0772   // combine the results and update the db object
0773   geometryUpdater.TrackerUpdate(endcapParameters, alignmentTubeParameters, *theAlignableTracker);
0774 
0775   /// laser hit section for trackbased interface
0776   ///
0777   /// due to the peculiar order of beams in TkLasBeamCollection,
0778   /// we cannot use the LASGlobalLoop object here
0779 
0780   // the collection container
0781   auto laserBeams = std::make_unique<TkLasBeamCollection>();
0782 
0783   // first for the endcap internal beams
0784   for (det = 0; det < 2; ++det) {
0785     for (ring = 0; ring < 2; ++ring) {
0786       for (beam = 0; beam < 8; ++beam) {
0787         // the beam and its identifier (see TkLasTrackBasedInterface TWiki)
0788         TkLasBeam currentBeam(100 * det + 10 * beam + ring);
0789 
0790         // order the hits in the beam by increasing z
0791         const int firstDisk = det == 0 ? 0 : 8;
0792         const int lastDisk = det == 0 ? 8 : 0;
0793 
0794         // count upwards or downwards
0795         for (disk = firstDisk; det == 0 ? disk <= lastDisk : disk >= lastDisk; det == 0 ? ++disk : --disk) {
0796           // detId for the SiStripLaserRecHit2D
0797           const SiStripDetId theDetId(detectorId.GetTECEntry(det, ring, beam, disk));
0798 
0799           // need this to calculate the localPosition and its error
0800           const StripGeomDetUnit* const theStripDet =
0801               dynamic_cast<const StripGeomDetUnit*>(theTracker.idToDet(theDetId));
0802 
0803           // the hit container
0804           const SiStripLaserRecHit2D currentHit(theStripDet->specificTopology().localPosition(
0805                                                     measuredStripPositions.GetTECEntry(det, ring, beam, disk).first),
0806                                                 theStripDet->specificTopology().localError(
0807                                                     measuredStripPositions.GetTECEntry(det, ring, beam, disk).first,
0808                                                     measuredStripPositions.GetTECEntry(det, ring, beam, disk).second),
0809                                                 theDetId);
0810 
0811           currentBeam.push_back(currentHit);
0812         }
0813 
0814         laserBeams->push_back(currentBeam);
0815       }
0816     }
0817   }
0818 
0819   // then, following the convention in TkLasTrackBasedInterface TWiki, the alignment tube beams;
0820   // they comprise hits in TIBTOB & TEC2TEC
0821 
0822   for (beam = 0; beam < 8; ++beam) {
0823     // the beam and its identifier (see TkLasTrackBasedInterface TWiki)
0824     TkLasBeam currentBeam(100 * 2 /*beamGroup=AT=2*/ + 10 * beam + 0 /*ring=0*/);
0825 
0826     // first: tec-
0827     det = 1;
0828     for (disk = 4; disk >= 0; --disk) {
0829       // detId for the SiStripLaserRecHit2D
0830       const SiStripDetId theDetId(detectorId.GetTEC2TECEntry(det, beam, disk));
0831 
0832       // need this to calculate the localPosition and its error
0833       const StripGeomDetUnit* const theStripDet = dynamic_cast<const StripGeomDetUnit*>(theTracker.idToDet(theDetId));
0834 
0835       // the hit container
0836       const SiStripLaserRecHit2D currentHit(
0837           theStripDet->specificTopology().localPosition(measuredStripPositions.GetTEC2TECEntry(det, beam, disk).first),
0838           theStripDet->specificTopology().localError(measuredStripPositions.GetTEC2TECEntry(det, beam, disk).first,
0839                                                      measuredStripPositions.GetTEC2TECEntry(det, beam, disk).second),
0840           theDetId);
0841 
0842       currentBeam.push_back(currentHit);
0843     }
0844 
0845     // now TIB and TOB in one go
0846     for (det = 2; det < 4; ++det) {
0847       for (pos = 5; pos >= 0; --pos) {  // stupidly, pos is defined from +z to -z in LASGlobalLoop
0848 
0849         // detId for the SiStripLaserRecHit2D
0850         const SiStripDetId theDetId(detectorId.GetTIBTOBEntry(det, beam, pos));
0851 
0852         // need this to calculate the localPosition and its error
0853         const StripGeomDetUnit* const theStripDet = dynamic_cast<const StripGeomDetUnit*>(theTracker.idToDet(theDetId));
0854 
0855         // the hit container
0856         const SiStripLaserRecHit2D currentHit(
0857             theStripDet->specificTopology().localPosition(measuredStripPositions.GetTIBTOBEntry(det, beam, pos).first),
0858             theStripDet->specificTopology().localError(measuredStripPositions.GetTIBTOBEntry(det, beam, pos).first,
0859                                                        measuredStripPositions.GetTIBTOBEntry(det, beam, pos).second),
0860             theDetId);
0861 
0862         currentBeam.push_back(currentHit);
0863       }
0864     }
0865 
0866     // then: tec+
0867     det = 0;
0868     for (disk = 0; disk < 5; ++disk) {
0869       // detId for the SiStripLaserRecHit2D
0870       const SiStripDetId theDetId(detectorId.GetTEC2TECEntry(det, beam, disk));
0871 
0872       // need this to calculate the localPosition and its error
0873       const StripGeomDetUnit* const theStripDet = dynamic_cast<const StripGeomDetUnit*>(theTracker.idToDet(theDetId));
0874 
0875       // the hit container
0876       const SiStripLaserRecHit2D currentHit(
0877           theStripDet->specificTopology().localPosition(measuredStripPositions.GetTEC2TECEntry(det, beam, disk).first),
0878           theStripDet->specificTopology().localError(measuredStripPositions.GetTEC2TECEntry(det, beam, disk).first,
0879                                                      measuredStripPositions.GetTEC2TECEntry(det, beam, disk).second),
0880           theDetId);
0881 
0882       currentBeam.push_back(currentHit);
0883     }
0884 
0885     // save this beam to the beamCollection
0886     laserBeams->push_back(currentBeam);
0887 
0888   }  // (close beam loop)
0889 
0890   // now attach the collection to the run
0891   theRun.put(std::move(laserBeams), "tkLaserBeams");
0892 
0893   // store the estimated alignment parameters into the DB
0894   // first get them
0895   Alignments alignments = *(theAlignableTracker->alignments());
0896   AlignmentErrorsExtended alignmentErrors = *(theAlignableTracker->alignmentErrors());
0897 
0898   if (theStoreToDB) {
0899     std::cout << " [LaserAlignment::endRun] -- Storing the calculated alignment parameters to the DataBase:"
0900               << std::endl;
0901 
0902     // Call service
0903     edm::Service<cond::service::PoolDBOutputService> poolDbService;
0904     if (!poolDbService.isAvailable())  // Die if not available
0905       throw cms::Exception("NotAvailable") << "PoolDBOutputService not available";
0906 
0907     // Store
0908 
0909     //     if ( poolDbService->isNewTagRequest(theAlignRecordName) ) {
0910     //       poolDbService->createOneIOV<Alignments>( alignments, poolDbService->currentTime(), theAlignRecordName );
0911     //     }
0912     //     else {
0913     //       poolDbService->appendOneIOV<Alignments>( alignments, poolDbService->currentTime(), theAlignRecordName );
0914     //     }
0915     poolDbService->writeOneIOV<Alignments>(alignments, poolDbService->beginOfTime(), theAlignRecordName);
0916 
0917     //     if ( poolDbService->isNewTagRequest(theErrorRecordName) ) {
0918     //       poolDbService->createOneIOV<AlignmentErrorsExtended>( alignmentErrors, poolDbService->currentTime(), poolDbService->endOfTime(), theErrorRecordName );
0919     //     }
0920     //     else {
0921     //       poolDbService->appendOneIOV<AlignmentErrorsExtended>( alignmentErrors, poolDbService->currentTime(), theErrorRecordName );
0922     //     }
0923     poolDbService->writeOneIOV<AlignmentErrorsExtended>(
0924         alignmentErrors, poolDbService->beginOfTime(), theErrorRecordName);
0925 
0926     std::cout << " [LaserAlignment::endRun] -- Storing done." << std::endl;
0927   }
0928 }
0929 
0930 ///
0931 ///
0932 ///
0933 void LaserAlignment::endJob() {}
0934 
0935 ///
0936 /// fills the module profiles (LASGlobalLoop<LASModuleProfile> currentDataProfiles)
0937 /// from the event digi containers, distinguishing between SiStripDigi or SiStripRawDigi.
0938 ///
0939 void LaserAlignment::fillDataProfiles(edm::Event const& theEvent, edm::EventSetup const& theSetup) {
0940   // two handles for the two different kinds of digis
0941   edm::Handle<edm::DetSetVector<SiStripRawDigi> > theStripRawDigis;
0942   edm::Handle<edm::DetSetVector<SiStripDigi> > theStripDigis;
0943 
0944   bool isRawDigi = false;
0945 
0946   // indices for the LASGlobalLoop object
0947   int det = 0, ring = 0, beam = 0, disk = 0, pos = 0;
0948 
0949   // query config set and loop over all PSets in the VPSet
0950   for (std::vector<edm::ParameterSet>::iterator itDigiProducersList = theDigiProducersList.begin();
0951        itDigiProducersList != theDigiProducersList.end();
0952        ++itDigiProducersList) {
0953     std::string digiProducer = itDigiProducersList->getParameter<std::string>("DigiProducer");
0954     std::string digiLabel = itDigiProducersList->getParameter<std::string>("DigiLabel");
0955     std::string digiType = itDigiProducersList->getParameter<std::string>("DigiType");
0956 
0957     // now branch according to digi type (raw or processed);
0958     // first we go for raw digis => SiStripRawDigi
0959     if (digiType == "Raw") {
0960       theEvent.getByLabel(digiProducer, digiLabel, theStripRawDigis);
0961       isRawDigi = true;
0962     } else if (digiType == "Processed") {
0963       theEvent.getByLabel(digiProducer, digiLabel, theStripDigis);
0964       isRawDigi = false;
0965     } else {
0966       throw cms::Exception(" [LaserAlignment::fillDataProfiles]")
0967           << " ** ERROR: Invalid digi type: \"" << digiType << "\" specified in configuration." << std::endl;
0968     }
0969 
0970     // loop TEC internal modules
0971     det = 0;
0972     ring = 0;
0973     beam = 0;
0974     disk = 0;
0975     do {
0976       // first clear the profile
0977       currentDataProfiles.GetTECEntry(det, ring, beam, disk).SetAllValuesTo(0.);
0978 
0979       // retrieve the raw id of that module
0980       const int detRawId = detectorId.GetTECEntry(det, ring, beam, disk);
0981 
0982       if (isRawDigi) {  // we have raw SiStripRawDigis
0983 
0984         // search the digis for the raw id
0985         edm::DetSetVector<SiStripRawDigi>::const_iterator detSetIter = theStripRawDigis->find(detRawId);
0986         if (detSetIter == theStripRawDigis->end()) {
0987           throw cms::Exception("[Laser Alignment::fillDataProfiles]")
0988               << " ** ERROR: No raw DetSet found for det: " << detRawId << "." << std::endl;
0989         }
0990 
0991         // fill the digis to the profiles
0992         edm::DetSet<SiStripRawDigi>::const_iterator digiRangeIterator = detSetIter->data.begin();  // for the loop
0993         edm::DetSet<SiStripRawDigi>::const_iterator digiRangeStart = digiRangeIterator;  // save starting positions
0994 
0995         // loop all digis
0996         for (; digiRangeIterator != detSetIter->data.end(); ++digiRangeIterator) {
0997           const SiStripRawDigi& digi = *digiRangeIterator;
0998           const int channel = distance(digiRangeStart, digiRangeIterator);
0999           if (channel >= 0 && channel < 512)
1000             currentDataProfiles.GetTECEntry(det, ring, beam, disk).SetValue(channel, digi.adc());
1001           else
1002             throw cms::Exception("[Laser Alignment::fillDataProfiles]")
1003                 << " ** ERROR: raw digi channel: " << channel << " out of range for det: " << detRawId << "."
1004                 << std::endl;
1005         }
1006 
1007       }
1008 
1009       else {  // we have zero suppressed SiStripDigis
1010 
1011         // search the digis for the raw id
1012         edm::DetSetVector<SiStripDigi>::const_iterator detSetIter = theStripDigis->find(detRawId);
1013 
1014         // processed DetSets may be missing, just skip
1015         if (detSetIter == theStripDigis->end())
1016           continue;
1017 
1018         // fill the digis to the profiles
1019         edm::DetSet<SiStripDigi>::const_iterator digiRangeIterator = detSetIter->data.begin();  // for the loop
1020 
1021         for (; digiRangeIterator != detSetIter->data.end(); ++digiRangeIterator) {
1022           const SiStripDigi& digi = *digiRangeIterator;
1023           if (digi.strip() < 512)
1024             currentDataProfiles.GetTECEntry(det, ring, beam, disk).SetValue(digi.strip(), digi.adc());
1025           else
1026             throw cms::Exception("[Laser Alignment::fillDataProfiles]")
1027                 << " ** ERROR: digi strip: " << digi.strip() << " out of range for det: " << detRawId << "."
1028                 << std::endl;
1029         }
1030       }
1031 
1032     } while (moduleLoop.TECLoop(det, ring, beam, disk));
1033 
1034     // loop TIBTOB modules
1035     det = 2;
1036     beam = 0;
1037     pos = 0;
1038     do {
1039       // first clear the profile
1040       currentDataProfiles.GetTIBTOBEntry(det, beam, pos).SetAllValuesTo(0.);
1041 
1042       // retrieve the raw id of that module
1043       const int detRawId = detectorId.GetTIBTOBEntry(det, beam, pos);
1044 
1045       if (isRawDigi) {  // we have raw SiStripRawDigis
1046 
1047         // search the digis for the raw id
1048         edm::DetSetVector<SiStripRawDigi>::const_iterator detSetIter = theStripRawDigis->find(detRawId);
1049         if (detSetIter == theStripRawDigis->end()) {
1050           throw cms::Exception("[Laser Alignment::fillDataProfiles]")
1051               << " ** ERROR: No raw DetSet found for det: " << detRawId << "." << std::endl;
1052         }
1053 
1054         // fill the digis to the profiles
1055         edm::DetSet<SiStripRawDigi>::const_iterator digiRangeIterator = detSetIter->data.begin();  // for the loop
1056         edm::DetSet<SiStripRawDigi>::const_iterator digiRangeStart = digiRangeIterator;  // save starting positions
1057 
1058         // loop all digis
1059         for (; digiRangeIterator != detSetIter->data.end(); ++digiRangeIterator) {
1060           const SiStripRawDigi& digi = *digiRangeIterator;
1061           const int channel = distance(digiRangeStart, digiRangeIterator);
1062           if (channel >= 0 && channel < 512)
1063             currentDataProfiles.GetTIBTOBEntry(det, beam, pos).SetValue(channel, digi.adc());
1064           else
1065             throw cms::Exception("[Laser Alignment::fillDataProfiles]")
1066                 << " ** ERROR: raw digi channel: " << channel << " out of range for det: " << detRawId << "."
1067                 << std::endl;
1068         }
1069 
1070       }
1071 
1072       else {  // we have zero suppressed SiStripDigis
1073 
1074         // search the digis for the raw id
1075         edm::DetSetVector<SiStripDigi>::const_iterator detSetIter = theStripDigis->find(detRawId);
1076 
1077         // processed DetSets may be missing, just skip
1078         if (detSetIter == theStripDigis->end())
1079           continue;
1080 
1081         // fill the digis to the profiles
1082         edm::DetSet<SiStripDigi>::const_iterator digiRangeIterator = detSetIter->data.begin();  // for the loop
1083 
1084         for (; digiRangeIterator != detSetIter->data.end(); ++digiRangeIterator) {
1085           const SiStripDigi& digi = *digiRangeIterator;
1086           if (digi.strip() < 512)
1087             currentDataProfiles.GetTIBTOBEntry(det, beam, pos).SetValue(digi.strip(), digi.adc());
1088           else
1089             throw cms::Exception("[Laser Alignment::fillDataProfiles]")
1090                 << " ** ERROR: digi strip: " << digi.strip() << " out of range for det: " << detRawId << "."
1091                 << std::endl;
1092         }
1093       }
1094 
1095     } while (moduleLoop.TIBTOBLoop(det, beam, pos));
1096 
1097     // loop TEC AT modules
1098     det = 0;
1099     beam = 0;
1100     disk = 0;
1101     do {
1102       // first clear the profile
1103       currentDataProfiles.GetTEC2TECEntry(det, beam, disk).SetAllValuesTo(0.);
1104 
1105       // retrieve the raw id of that module
1106       const int detRawId = detectorId.GetTEC2TECEntry(det, beam, disk);
1107 
1108       if (isRawDigi) {  // we have raw SiStripRawDigis
1109 
1110         // search the digis for the raw id
1111         edm::DetSetVector<SiStripRawDigi>::const_iterator detSetIter = theStripRawDigis->find(detRawId);
1112         if (detSetIter == theStripRawDigis->end()) {
1113           throw cms::Exception("[Laser Alignment::fillDataProfiles]")
1114               << " ** ERROR: No raw DetSet found for det: " << detRawId << "." << std::endl;
1115         }
1116 
1117         // fill the digis to the profiles
1118         edm::DetSet<SiStripRawDigi>::const_iterator digiRangeIterator = detSetIter->data.begin();  // for the loop
1119         edm::DetSet<SiStripRawDigi>::const_iterator digiRangeStart = digiRangeIterator;  // save starting positions
1120 
1121         // loop all digis
1122         for (; digiRangeIterator != detSetIter->data.end(); ++digiRangeIterator) {
1123           const SiStripRawDigi& digi = *digiRangeIterator;
1124           const int channel = distance(digiRangeStart, digiRangeIterator);
1125           if (channel >= 0 && channel < 512)
1126             currentDataProfiles.GetTEC2TECEntry(det, beam, disk).SetValue(channel, digi.adc());
1127           else
1128             throw cms::Exception("[Laser Alignment::fillDataProfiles]")
1129                 << " ** ERROR: raw digi channel: " << channel << " out of range for det: " << detRawId << "."
1130                 << std::endl;
1131         }
1132 
1133       }
1134 
1135       else {  // we have zero suppressed SiStripDigis
1136 
1137         // search the digis for the raw id
1138         edm::DetSetVector<SiStripDigi>::const_iterator detSetIter = theStripDigis->find(detRawId);
1139 
1140         // processed DetSets may be missing, just skip
1141         if (detSetIter == theStripDigis->end())
1142           continue;
1143 
1144         // fill the digis to the profiles
1145         edm::DetSet<SiStripDigi>::const_iterator digiRangeIterator = detSetIter->data.begin();  // for the loop
1146 
1147         for (; digiRangeIterator != detSetIter->data.end(); ++digiRangeIterator) {
1148           const SiStripDigi& digi = *digiRangeIterator;
1149           if (digi.strip() < 512)
1150             currentDataProfiles.GetTEC2TECEntry(det, beam, disk).SetValue(digi.strip(), digi.adc());
1151           else
1152             throw cms::Exception("[Laser Alignment::fillDataProfiles]")
1153                 << " ** ERROR: digi strip: " << digi.strip() << " out of range for det: " << detRawId << "."
1154                 << std::endl;
1155         }
1156       }
1157 
1158     } while (moduleLoop.TEC2TECLoop(det, beam, disk));
1159 
1160   }  // theDigiProducersList loop
1161 }
1162 
1163 ///
1164 /// This function fills the pedestal profiles (LASGlobalData<LASModuleProfiles> pedestalProfiles)
1165 /// from the ESHandle (from file or DB)
1166 ///
1167 /// Argument: readily connected SiStripPedestals object (get() alredy called)
1168 /// The functionality inside the loops is basically taken from:
1169 /// CommonTools/SiStripZeroSuppression/src/SiStripPedestalsSubtractor.cc
1170 ///
1171 void LaserAlignment::fillPedestalProfiles(edm::ESHandle<SiStripPedestals>& pedestalsHandle) {
1172   int det, ring, beam, disk, pos;
1173 
1174   // loop TEC modules (yet without AT)
1175   det = 0;
1176   ring = 0;
1177   beam = 0;
1178   disk = 0;
1179   do {  // loop using LASGlobalLoop functionality
1180     SiStripPedestals::Range pedRange = pedestalsHandle->getRange(detectorId.GetTECEntry(det, ring, beam, disk));
1181     for (int strip = 0; strip < 512; ++strip) {
1182       int thePedestal = int(pedestalsHandle->getPed(strip, pedRange));
1183       if (thePedestal > 895)
1184         thePedestal -= 1024;
1185       pedestalProfiles.GetTECEntry(det, ring, beam, disk).SetValue(strip, thePedestal);
1186     }
1187   } while (moduleLoop.TECLoop(det, ring, beam, disk));
1188 
1189   // TIB & TOB section
1190   det = 2;
1191   beam = 0;
1192   pos = 0;
1193   do {  // loop using LASGlobalLoop functionality
1194     SiStripPedestals::Range pedRange = pedestalsHandle->getRange(detectorId.GetTIBTOBEntry(det, beam, pos));
1195     for (int strip = 0; strip < 512; ++strip) {
1196       int thePedestal = int(pedestalsHandle->getPed(strip, pedRange));
1197       if (thePedestal > 895)
1198         thePedestal -= 1024;
1199       pedestalProfiles.GetTIBTOBEntry(det, beam, pos).SetValue(strip, thePedestal);
1200     }
1201   } while (moduleLoop.TIBTOBLoop(det, beam, pos));
1202 
1203   // TEC2TEC AT section
1204   det = 0;
1205   beam = 0;
1206   disk = 0;
1207   do {  // loop using LASGlobalLoop functionality
1208     SiStripPedestals::Range pedRange = pedestalsHandle->getRange(detectorId.GetTEC2TECEntry(det, beam, disk));
1209     for (int strip = 0; strip < 512; ++strip) {
1210       int thePedestal = int(pedestalsHandle->getPed(strip, pedRange));
1211       if (thePedestal > 895)
1212         thePedestal -= 1024;
1213       pedestalProfiles.GetTEC2TECEntry(det, beam, disk).SetValue(strip, thePedestal);
1214     }
1215   } while (moduleLoop.TEC2TECLoop(det, beam, disk));
1216 }
1217 
1218 ///
1219 /// count useable profiles in TEC,
1220 /// operates on LASGlobalData<int> LaserAlignment::isAcceptedProfile
1221 /// to allow for more elaborate patterns in the future
1222 ///
1223 bool LaserAlignment::isTECBeam(void) {
1224   int numberOfProfiles = 0;
1225 
1226   int ring = 1;  // search all ring6 modules for signals
1227   for (int det = 0; det < 2; ++det) {
1228     for (int beam = 0; beam < 8; ++beam) {
1229       for (int disk = 0; disk < 9; ++disk) {
1230         if (isAcceptedProfile.GetTECEntry(det, ring, beam, disk) == 1)
1231           numberOfProfiles++;
1232       }
1233     }
1234   }
1235 
1236   LogDebug("[LaserAlignment::isTECBeam]") << " Found: " << numberOfProfiles << "hits." << std::endl;
1237   std::cout << " [LaserAlignment::isTECBeam] -- Found: " << numberOfProfiles << " hits." << std::endl;  ////
1238 
1239   if (numberOfProfiles > 10)
1240     return (true);
1241   return (false);
1242 }
1243 
1244 ///
1245 /// count useable profiles in TIBTOB,
1246 /// operates on LASGlobalData<bool> LaserAlignment::isAcceptedProfile
1247 /// to allow for more elaborate patterns in the future
1248 ///
1249 
1250 bool LaserAlignment::isATBeam(void) {
1251   int numberOfProfiles = 0;
1252 
1253   int det = 2;
1254   int beam = 0;
1255   int pos = 0;  // search all TIB/TOB for signals
1256   do {
1257     if (isAcceptedProfile.GetTIBTOBEntry(det, beam, pos) == 1)
1258       numberOfProfiles++;
1259   } while (moduleLoop.TIBTOBLoop(det, beam, pos));
1260 
1261   LogDebug("[LaserAlignment::isATBeam]") << " Found: " << numberOfProfiles << "hits." << std::endl;
1262   std::cout << " [LaserAlignment::isATBeam] -- Found: " << numberOfProfiles << " hits." << std::endl;  /////
1263 
1264   if (numberOfProfiles > 10)
1265     return (true);
1266   return (false);
1267 }
1268 
1269 ///
1270 /// not all TIB & TOB modules are hit in the center;
1271 /// this func returns the nominal beam offset locally on a module (in strips)
1272 /// for the ProfileJudge and the LASPeakFinder in strips.
1273 /// (offset = middle of module - nominal position)
1274 ///
1275 /// the hard coded numbers will later be supplied by a special geometry class..
1276 ///
1277 double LaserAlignment::getTIBTOBNominalBeamOffset(unsigned int det, unsigned int beam, unsigned int pos) {
1278   if (det < 2 || det > 3 || beam > 7 || pos > 5) {
1279     throw cms::Exception("[LaserAlignment::getTIBTOBNominalBeamOffset]")
1280         << " ERROR ** Called with nonexisting parameter set: det " << det << " beam " << beam << " pos " << pos << "."
1281         << std::endl;
1282   }
1283 
1284   const double nominalOffsetsTIB[8] = {
1285       0.00035, 2.10687, -2.10827, -0.00173446, 2.10072, -0.00135114, 2.10105, -2.10401};
1286 
1287   // in tob, modules have alternating orientations along the rods.
1288   // this is described by the following pattern.
1289   // (even more confusing, this pattern is inversed for beams 0, 5, 6, 7)
1290   const int orientationPattern[6] = {-1, 1, 1, -1, -1, 1};
1291   const double nominalOffsetsTOB[8] = {0.00217408, 1.58678, 117.733, 119.321, 120.906, 119.328, 117.743, 1.58947};
1292 
1293   if (det == 2)
1294     return (-1. * nominalOffsetsTIB[beam]);
1295 
1296   else {
1297     if (beam == 0 or beam > 4)
1298       return (nominalOffsetsTOB[beam] * orientationPattern[pos]);
1299     else
1300       return (-1. * nominalOffsetsTOB[beam] * orientationPattern[pos]);
1301   }
1302 }
1303 
1304 ///
1305 /// not all TEC-AT modules are hit in the center;
1306 /// this func returns the nominal beam offset locally on a module (in strips)
1307 /// for the ProfileJudge and the LASPeakFinder in strips.
1308 /// (offset = middle of module - nominal position)
1309 ///
1310 /// the hard coded numbers will later be supplied by a special geometry class..
1311 ///
1312 double LaserAlignment::getTEC2TECNominalBeamOffset(unsigned int det, unsigned int beam, unsigned int disk) {
1313   if (det > 1 || beam > 7 || disk > 5) {
1314     throw cms::Exception("[LaserAlignment::getTEC2TECNominalBeamOffset]")
1315         << " ERROR ** Called with nonexisting parameter set: det " << det << " beam " << beam << " disk " << disk << "."
1316         << std::endl;
1317   }
1318 
1319   const double nominalOffsets[8] = {0., 2.220, -2.221, 0., 2.214, 0., 2.214, -2.217};
1320 
1321   if (det == 0)
1322     return -1. * nominalOffsets[beam];
1323   else
1324     return nominalOffsets[beam];
1325 }
1326 
1327 ///
1328 ///
1329 ///
1330 void LaserAlignment::CalculateNominalCoordinates(void) {
1331   //
1332   // hard coded data yet...
1333   //
1334 
1335   // nominal phi values of tec beam / alignment tube hits (parameter is beam 0-7)
1336   const double tecPhiPositions[8] = {
1337       0.392699, 1.178097, 1.963495, 2.748894, 3.534292, 4.319690, 5.105088, 5.890486};  // new values calculated by maple
1338   const double atPhiPositions[8] = {
1339       0.392699, 1.289799, 1.851794, 2.748894, 3.645995, 4.319690, 5.216791, 5.778784};  // new values calculated by maple
1340 
1341   // nominal r values (mm) of hits
1342   const double tobRPosition = 600.;
1343   const double tibRPosition = 514.;
1344   const double tecRPosition[2] = {564., 840.};  // ring 4,6
1345 
1346   // nominal z values (mm) of hits in barrel (parameter is pos 0-6)
1347   const double tobZPosition[6] = {1040., 580., 220., -140., -500., -860.};
1348   const double tibZPosition[6] = {620., 380., 180., -100., -340., -540.};
1349 
1350   // nominal z values (mm) of hits in tec (parameter is disk 0-8); FOR TEC-: (* -1.)
1351   const double tecZPosition[9] = {1322.5, 1462.5, 1602.5, 1742.5, 1882.5, 2057.5, 2247.5, 2452.5, 2667.5};
1352 
1353   //
1354   // now we fill these into the nominalCoordinates container;
1355   // errors are zero for nominal values..
1356   //
1357 
1358   // loop object and its variables
1359   LASGlobalLoop moduleLoop;
1360   int det, ring, beam, disk, pos;
1361 
1362   // TEC+- section
1363   det = 0;
1364   ring = 0, beam = 0;
1365   disk = 0;
1366   do {
1367     if (det == 0) {  // this is TEC+
1368       nominalCoordinates.SetTECEntry(
1369           det,
1370           ring,
1371           beam,
1372           disk,
1373           LASCoordinateSet(tecPhiPositions[beam], 0., tecRPosition[ring], 0., tecZPosition[disk], 0.));
1374     } else {  // now TEC-
1375       nominalCoordinates.SetTECEntry(
1376           det,
1377           ring,
1378           beam,
1379           disk,
1380           LASCoordinateSet(
1381               tecPhiPositions[beam], 0., tecRPosition[ring], 0., -1. * tecZPosition[disk], 0.));  // just * -1.
1382     }
1383 
1384   } while (moduleLoop.TECLoop(det, ring, beam, disk));
1385 
1386   // TIB & TOB section
1387   det = 2;
1388   beam = 0;
1389   pos = 0;
1390   do {
1391     if (det == 2) {  // this is TIB
1392       nominalCoordinates.SetTIBTOBEntry(
1393           det, beam, pos, LASCoordinateSet(atPhiPositions[beam], 0., tibRPosition, 0., tibZPosition[pos], 0.));
1394     } else {  // now TOB
1395       nominalCoordinates.SetTIBTOBEntry(
1396           det, beam, pos, LASCoordinateSet(atPhiPositions[beam], 0., tobRPosition, 0., tobZPosition[pos], 0.));
1397     }
1398 
1399   } while (moduleLoop.TIBTOBLoop(det, beam, pos));
1400 
1401   // TEC2TEC AT section
1402   det = 0;
1403   beam = 0;
1404   disk = 0;
1405   do {
1406     if (det == 0) {  // this is TEC+, ring4 only
1407       nominalCoordinates.SetTEC2TECEntry(
1408           det, beam, disk, LASCoordinateSet(atPhiPositions[beam], 0., tecRPosition[0], 0., tecZPosition[disk], 0.));
1409     } else {  // now TEC-
1410       nominalCoordinates.SetTEC2TECEntry(
1411           det,
1412           beam,
1413           disk,
1414           LASCoordinateSet(atPhiPositions[beam], 0., tecRPosition[0], 0., -1. * tecZPosition[disk], 0.));  // just * -1.
1415     }
1416 
1417   } while (moduleLoop.TEC2TECLoop(det, beam, disk));
1418 }
1419 
1420 ///
1421 /// convert an angle in the [-pi,pi] range
1422 /// to the [0,2*pi] range
1423 ///
1424 double LaserAlignment::ConvertAngle(double angle) {
1425   if (angle < -1. * M_PI || angle > M_PI) {
1426     throw cms::Exception(" [LaserAlignment::ConvertAngle] ")
1427         << "** ERROR: Called with illegal input angle: " << angle << "." << std::endl;
1428   }
1429 
1430   if (angle >= 0.)
1431     return angle;
1432   else
1433     return (angle + 2. * M_PI);
1434 }
1435 
1436 ///
1437 /// debug only, will disappear
1438 ///
1439 void LaserAlignment::DumpPosFileSet(LASGlobalData<LASCoordinateSet>& coordinates) {
1440   LASGlobalLoop loop;
1441   int det, ring, beam, disk, pos;
1442 
1443   std::cout << std::endl << " [LaserAlignment::DumpPosFileSet] -- Dump: " << std::endl;
1444 
1445   // TEC INTERNAL
1446   det = 0;
1447   ring = 0;
1448   beam = 0;
1449   disk = 0;
1450   do {
1451     std::cout << "POS " << det << "\t" << beam << "\t" << disk << "\t" << ring << "\t"
1452               << coordinates.GetTECEntry(det, ring, beam, disk).GetPhi() << "\t"
1453               << coordinates.GetTECEntry(det, ring, beam, disk).GetPhiError() << std::endl;
1454   } while (loop.TECLoop(det, ring, beam, disk));
1455 
1456   // TIBTOB
1457   det = 2;
1458   beam = 0;
1459   pos = 0;
1460   do {
1461     std::cout << "POS " << det << "\t" << beam << "\t" << pos << "\t"
1462               << "-1"
1463               << "\t" << coordinates.GetTIBTOBEntry(det, beam, pos).GetPhi() << "\t"
1464               << coordinates.GetTIBTOBEntry(det, beam, pos).GetPhiError() << std::endl;
1465   } while (loop.TIBTOBLoop(det, beam, pos));
1466 
1467   // TEC2TEC
1468   det = 0;
1469   beam = 0;
1470   disk = 0;
1471   do {
1472     std::cout << "POS " << det << "\t" << beam << "\t" << disk << "\t"
1473               << "-1"
1474               << "\t" << coordinates.GetTEC2TECEntry(det, beam, disk).GetPhi() << "\t"
1475               << coordinates.GetTEC2TECEntry(det, beam, disk).GetPhiError() << std::endl;
1476   } while (loop.TEC2TECLoop(det, beam, disk));
1477 
1478   std::cout << std::endl << " [LaserAlignment::DumpPosFileSet] -- End dump: " << std::endl;
1479 }
1480 
1481 ///
1482 ///
1483 ///
1484 void LaserAlignment::DumpStripFileSet(LASGlobalData<std::pair<float, float> >& measuredStripPositions) {
1485   LASGlobalLoop loop;
1486   int det, ring, beam, disk, pos;
1487 
1488   std::cout << std::endl << " [LaserAlignment::DumpStripFileSet] -- Dump: " << std::endl;
1489 
1490   // TEC INTERNAL
1491   det = 0;
1492   ring = 0;
1493   beam = 0;
1494   disk = 0;
1495   do {
1496     std::cout << "STRIP " << det << "\t" << beam << "\t" << disk << "\t" << ring << "\t"
1497               << measuredStripPositions.GetTECEntry(det, ring, beam, disk).first << "\t"
1498               << measuredStripPositions.GetTECEntry(det, ring, beam, disk).second << std::endl;
1499   } while (loop.TECLoop(det, ring, beam, disk));
1500 
1501   // TIBTOB
1502   det = 2;
1503   beam = 0;
1504   pos = 0;
1505   do {
1506     std::cout << "STRIP " << det << "\t" << beam << "\t" << pos << "\t"
1507               << "-1"
1508               << "\t" << measuredStripPositions.GetTIBTOBEntry(det, beam, pos).first << "\t"
1509               << measuredStripPositions.GetTIBTOBEntry(det, beam, pos).second << std::endl;
1510   } while (loop.TIBTOBLoop(det, beam, pos));
1511 
1512   // TEC2TEC
1513   det = 0;
1514   beam = 0;
1515   disk = 0;
1516   do {
1517     std::cout << "STRIP " << det << "\t" << beam << "\t" << disk << "\t"
1518               << "-1"
1519               << "\t" << measuredStripPositions.GetTEC2TECEntry(det, beam, disk).first << "\t"
1520               << measuredStripPositions.GetTEC2TECEntry(det, beam, disk).second << std::endl;
1521   } while (loop.TEC2TECLoop(det, beam, disk));
1522 
1523   std::cout << std::endl << " [LaserAlignment::DumpStripFileSet] -- End dump: " << std::endl;
1524 }
1525 
1526 ///
1527 ///
1528 ///
1529 void LaserAlignment::DumpHitmaps(LASGlobalData<int>& numberOfAcceptedProfiles) {
1530   std::cout << " [LaserAlignment::DumpHitmaps] -- Dumping hitmap for TEC+:" << std::endl;
1531   std::cout << " [LaserAlignment::DumpHitmaps] -- Ring4:" << std::endl;
1532   std::cout << "     disk0   disk1   disk2   disk3   disk4   disk5   disk6   disk7   disk8" << std::endl;
1533 
1534   for (int beam = 0; beam < 8; ++beam) {
1535     std::cout << " beam" << beam << ":";
1536     for (int disk = 0; disk < 9; ++disk) {
1537       std::cout << "\t" << numberOfAcceptedProfiles.GetTECEntry(0, 0, beam, disk);
1538     }
1539     std::cout << std::endl;
1540   }
1541 
1542   std::cout << " [LaserAlignment::DumpHitmaps] -- Ring6:" << std::endl;
1543   std::cout << "     disk0   disk1   disk2   disk3   disk4   disk5   disk6   disk7   disk8" << std::endl;
1544 
1545   for (int beam = 0; beam < 8; ++beam) {
1546     std::cout << " beam" << beam << ":";
1547     for (int disk = 0; disk < 9; ++disk) {
1548       std::cout << "\t" << numberOfAcceptedProfiles.GetTECEntry(0, 1, beam, disk);
1549     }
1550     std::cout << std::endl;
1551   }
1552 
1553   std::cout << " [LaserAlignment::DumpHitmaps] -- Dumping hitmap for TEC-:" << std::endl;
1554   std::cout << " [LaserAlignment::DumpHitmaps] -- Ring4:" << std::endl;
1555   std::cout << "     disk0   disk1   disk2   disk3   disk4   disk5   disk6   disk7   disk8" << std::endl;
1556 
1557   for (int beam = 0; beam < 8; ++beam) {
1558     std::cout << " beam" << beam << ":";
1559     for (int disk = 0; disk < 9; ++disk) {
1560       std::cout << "\t" << numberOfAcceptedProfiles.GetTECEntry(1, 0, beam, disk);
1561     }
1562     std::cout << std::endl;
1563   }
1564 
1565   std::cout << " [LaserAlignment::DumpHitmaps] -- Ring6:" << std::endl;
1566   std::cout << "     disk0   disk1   disk2   disk3   disk4   disk5   disk6   disk7   disk8" << std::endl;
1567 
1568   for (int beam = 0; beam < 8; ++beam) {
1569     std::cout << " beam" << beam << ":";
1570     for (int disk = 0; disk < 9; ++disk) {
1571       std::cout << "\t" << numberOfAcceptedProfiles.GetTECEntry(1, 1, beam, disk);
1572     }
1573     std::cout << std::endl;
1574   }
1575 
1576   std::cout << " [LaserAlignment::DumpHitmaps] -- End of dump." << std::endl << std::endl;
1577 }
1578 
1579 ///
1580 /// loop the list of endcap modules to be masked and
1581 /// apply the corrections from the "endcapParameters" to them
1582 ///
1583 void LaserAlignment::ApplyEndcapMaskingCorrections(LASGlobalData<LASCoordinateSet>& measuredCoordinates,
1584                                                    LASGlobalData<LASCoordinateSet>& nominalCoordinates,
1585                                                    LASEndcapAlignmentParameterSet& endcapParameters) {
1586   // loop the list of modules to be masked
1587   for (std::vector<unsigned int>::iterator moduleIt = theMaskTecModules.begin(); moduleIt != theMaskTecModules.end();
1588        ++moduleIt) {
1589     // loop variables
1590     LASGlobalLoop moduleLoop;
1591     int det, ring, beam, disk;
1592 
1593     // this will calculate the corrections from the alignment parameters
1594     LASEndcapAlgorithm endcapAlgorithm;
1595 
1596     // find the location of the respective module in the container with this loop
1597     det = 0;
1598     ring = 0;
1599     beam = 0;
1600     disk = 0;
1601     do {
1602       // here we got it
1603       if (detectorId.GetTECEntry(det, ring, beam, disk) == *moduleIt) {
1604         // the nominal phi value for this module
1605         const double nominalPhi = nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetPhi();
1606 
1607         // the offset from the alignment parameters
1608         const double phiCorrection = endcapAlgorithm.GetAlignmentParameterCorrection(
1609             det, ring, beam, disk, nominalCoordinates, endcapParameters);
1610 
1611         // apply the corrections
1612         measuredCoordinates.GetTECEntry(det, ring, beam, disk).SetPhi(nominalPhi - phiCorrection);
1613       }
1614 
1615     } while (moduleLoop.TECLoop(det, ring, beam, disk));
1616   }
1617 }
1618 
1619 ///
1620 /// loop the list of alignment tube modules to be masked and
1621 /// apply the corrections from the "barrelParameters" to them
1622 ///
1623 void LaserAlignment::ApplyATMaskingCorrections(LASGlobalData<LASCoordinateSet>& measuredCoordinates,
1624                                                LASGlobalData<LASCoordinateSet>& nominalCoordinates,
1625                                                LASBarrelAlignmentParameterSet& atParameters) {
1626   // loop the list of modules to be masked
1627   for (std::vector<unsigned int>::iterator moduleIt = theMaskAtModules.begin(); moduleIt != theMaskAtModules.end();
1628        ++moduleIt) {
1629     // loop variables
1630     LASGlobalLoop moduleLoop;
1631     int det, beam, disk, pos;
1632 
1633     // this will calculate the corrections from the alignment parameters
1634     LASAlignmentTubeAlgorithm atAlgorithm;
1635 
1636     // find the location of the respective module in the container with these loops:
1637 
1638     // first TIB+TOB
1639     det = 2;
1640     beam = 0;
1641     pos = 0;
1642     do {
1643       // here we got it
1644       if (detectorId.GetTIBTOBEntry(det, beam, pos) == *moduleIt) {
1645         // the nominal phi value for this module
1646         const double nominalPhi = nominalCoordinates.GetTIBTOBEntry(det, beam, pos).GetPhi();
1647 
1648         // the offset from the alignment parameters
1649         const double phiCorrection =
1650             atAlgorithm.GetTIBTOBAlignmentParameterCorrection(det, beam, pos, nominalCoordinates, atParameters);
1651 
1652         // apply the corrections
1653         measuredCoordinates.GetTIBTOBEntry(det, beam, pos).SetPhi(nominalPhi - phiCorrection);
1654       }
1655 
1656     } while (moduleLoop.TIBTOBLoop(det, beam, pos));
1657 
1658     // then TEC(AT)
1659     det = 0;
1660     beam = 0;
1661     disk = 0;
1662     do {
1663       // here we got it
1664       if (detectorId.GetTEC2TECEntry(det, beam, disk) == *moduleIt) {
1665         // the nominal phi value for this module
1666         const double nominalPhi = nominalCoordinates.GetTEC2TECEntry(det, beam, disk).GetPhi();
1667 
1668         // the offset from the alignment parameters
1669         const double phiCorrection =
1670             atAlgorithm.GetTEC2TECAlignmentParameterCorrection(det, beam, disk, nominalCoordinates, atParameters);
1671 
1672         // apply the corrections
1673         measuredCoordinates.GetTEC2TECEntry(det, beam, disk).SetPhi(nominalPhi - phiCorrection);
1674       }
1675 
1676     } while (moduleLoop.TEC2TECLoop(det, beam, disk));
1677   }
1678 }
1679 
1680 ///
1681 /// this function is for debugging and testing only
1682 /// and will disappear..
1683 ///
1684 void LaserAlignment::testRoutine(void) {
1685   // tracker geom. object for calculating the global beam positions
1686   const TrackerGeometry& theTracker(*theTrackerGeometry);
1687 
1688   const double atPhiPositions[8] = {0.392699, 1.289799, 1.851794, 2.748894, 3.645995, 4.319690, 5.216791, 5.778784};
1689   const double tecPhiPositions[8] = {0.392699, 1.178097, 1.963495, 2.748894, 3.534292, 4.319690, 5.105088, 5.890486};
1690   const double zPositions[9] = {125.0, 139.0, 153.0, 167.0, 181.0, 198.5, 217.5, 238.0, 259.5};
1691   const double zPositionsTIB[6] = {62.0, 38.0, 18.0, -10.0, -34.0, -54.0};
1692   const double zPositionsTOB[6] = {104.0, 58.0, 22.0, -14.0, -50.0, -86.0};
1693 
1694   int det, beam, disk, pos, ring;
1695 
1696   // loop TEC+- internal
1697   det = 0;
1698   ring = 0;
1699   beam = 0;
1700   disk = 0;
1701   do {
1702     const double radius = ring ? 84.0 : 56.4;
1703 
1704     // access the tracker geometry for this module
1705     const DetId theDetId(detectorId.GetTECEntry(det, ring, beam, disk));
1706     const StripGeomDetUnit* const theStripDet = dynamic_cast<const StripGeomDetUnit*>(theTracker.idToDet(theDetId));
1707 
1708     if (theStripDet) {
1709       const GlobalPoint gp(GlobalPoint::Cylindrical(radius, tecPhiPositions[beam], zPositions[disk]));
1710 
1711       const LocalPoint lp(theStripDet->surface().toLocal(gp));
1712       std::cout << "__TEC: " << 256. - theStripDet->specificTopology().strip(lp)
1713                 << std::endl;  /////////////////////////////////
1714     }
1715 
1716   } while (moduleLoop.TECLoop(det, ring, beam, disk));
1717 
1718   // loop TIBTOB
1719   det = 2;
1720   beam = 0;
1721   pos = 0;
1722   do {
1723     const double radius =
1724         (det == 2 ? 51.4 : 58.4);  /////////////////////////////////////////////////////////////////////////////
1725     const double theZ = (det == 2 ? zPositionsTIB[pos] : zPositionsTOB[pos]);
1726 
1727     // access the tracker geometry for this module
1728     const DetId theDetId(detectorId.GetTIBTOBEntry(det, beam, pos));
1729     const StripGeomDetUnit* const theStripDet = dynamic_cast<const StripGeomDetUnit*>(theTracker.idToDet(theDetId));
1730 
1731     if (theStripDet) {
1732       const GlobalPoint gp(GlobalPoint::Cylindrical(radius, atPhiPositions[beam], theZ));
1733 
1734       const LocalPoint lp(theStripDet->surface().toLocal(gp));
1735       std::cout << "__TIBTOB det " << det << " beam " << beam << " pos " << pos << "  "
1736                 << 256. - theStripDet->specificTopology().strip(lp);
1737       std::cout << "           " << theStripDet->position().perp() << std::endl;  /////////////////////////////////
1738     }
1739 
1740   } while (moduleLoop.TIBTOBLoop(det, beam, pos));
1741 
1742   // loop TEC2TEC
1743   det = 0;
1744   beam = 0;
1745   disk = 0;
1746   do {
1747     const double radius = 56.4;
1748 
1749     // access the tracker geometry for this module
1750     const DetId theDetId(detectorId.GetTEC2TECEntry(det, beam, disk));
1751     const StripGeomDetUnit* const theStripDet = dynamic_cast<const StripGeomDetUnit*>(theTracker.idToDet(theDetId));
1752 
1753     if (theStripDet) {
1754       const GlobalPoint gp(GlobalPoint::Cylindrical(radius, atPhiPositions[beam], zPositions[disk]));
1755 
1756       const LocalPoint lp(theStripDet->surface().toLocal(gp));
1757       std::cout << "__TEC2TEC det " << det << " beam " << beam << " disk " << disk << "  "
1758                 << 256. - theStripDet->specificTopology().strip(lp) << std::endl;  /////////////////////////////////
1759     }
1760 
1761   } while (moduleLoop.TEC2TECLoop(det, beam, disk));
1762 }
1763 
1764 // define the SEAL module
1765 #include "FWCore/Framework/interface/MakerMacros.h"
1766 
1767 DEFINE_FWK_MODULE(LaserAlignment);
1768 
1769 // the ATTIC
[ 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!