Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-12-17 23:35:49 

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