LaserAlignment.cc

Go to the documentation of this file.

#include "Alignment/LaserAlignment/plugins/LaserAlignment.h"
#include "FWCore/Framework/interface/Run.h"
#include "Geometry/Records/interface/TrackerTopologyRcd.h"
#include "CondFormats/GeometryObjects/interface/PTrackerParameters.h"
#include "Geometry/Records/interface/PTrackerParametersRcd.h"

LaserAlignment::LaserAlignment( edm::ParameterSet const& theConf ) :
  theEvents(0), 
  theDoPedestalSubtraction( theConf.getUntrackedParameter<bool>( "SubtractPedestals", true ) ),
  theUseMinuitAlgorithm( theConf.getUntrackedParameter<bool>( "RunMinuitAlignmentTubeAlgorithm", false ) ),
  theApplyBeamKinkCorrections( theConf.getUntrackedParameter<bool>( "ApplyBeamKinkCorrections", true ) ),
  peakFinderThreshold( theConf.getUntrackedParameter<double>( "PeakFinderThreshold", 10. ) ),
  enableJudgeZeroFilter( theConf.getUntrackedParameter<bool>( "EnableJudgeZeroFilter", true ) ),
  judgeOverdriveThreshold( theConf.getUntrackedParameter<unsigned int>( "JudgeOverdriveThreshold", 220 ) ),
  updateFromInputGeometry( theConf.getUntrackedParameter<bool>( "UpdateFromInputGeometry", false ) ),
  misalignedByRefGeometry( theConf.getUntrackedParameter<bool>( "MisalignedByRefGeometry", false ) ),
  theStoreToDB ( theConf.getUntrackedParameter<bool>( "SaveToDbase", false ) ),
  theDigiProducersList( theConf.getParameter<std::vector<edm::ParameterSet> >( "DigiProducersList" ) ),
  theSaveHistograms( theConf.getUntrackedParameter<bool>( "SaveHistograms", false ) ),
  theCompression( theConf.getUntrackedParameter<int>( "ROOTFileCompression", 1 ) ),
  theFileName( theConf.getUntrackedParameter<std::string>( "ROOTFileName", "test.root" ) ),
  theMaskTecModules( theConf.getUntrackedParameter<std::vector<unsigned int> >( "MaskTECModules" ) ),
  theMaskAtModules( theConf.getUntrackedParameter<std::vector<unsigned int> >( "MaskATModules" ) ),
  theSetNominalStrips( theConf.getUntrackedParameter<bool>( "ForceFitterToNominalStrips", false ) ),
  theLasConstants( theConf.getUntrackedParameter<std::vector<edm::ParameterSet> >( "LaserAlignmentConstants" ) ),
  theFile(),
  theAlignableTracker(),
  theAlignRecordName( "TrackerAlignmentRcd" ),
  theErrorRecordName( "TrackerAlignmentErrorExtendedRcd" ),
  firstEvent_(true)
{


  std::cout << std::endl;
  std::cout <<   "=============================================================="
        << "\n===         LaserAlignment module configuration            ==="
        << "\n"
        << "\n    Write histograms to file       = " << (theSaveHistograms?"true":"false")
        << "\n    Histogram file name            = " << theFileName
        << "\n    Histogram file compression     = " << theCompression
        << "\n    Subtract pedestals             = " << (theDoPedestalSubtraction?"true":"false")
        << "\n    Run Minuit AT algorithm        = " << (theUseMinuitAlgorithm?"true":"false")
        << "\n    Apply beam kink corrections    = " << (theApplyBeamKinkCorrections?"true":"false")
        << "\n    Peak Finder Threshold          = " << peakFinderThreshold
        << "\n    EnableJudgeZeroFilter          = " << (enableJudgeZeroFilter?"true":"false")
        << "\n    JudgeOverdriveThreshold        = " << judgeOverdriveThreshold
        << "\n    Update from input geometry     = " << (updateFromInputGeometry?"true":"false")
        << "\n    Misalignment from ref geometry = " << (misalignedByRefGeometry?"true":"false")
        << "\n    Number of TEC modules masked   = " << theMaskTecModules.size() << " (s. below list if > 0)"
        << "\n    Number of AT modules masked    = " << theMaskAtModules.size()  << " (s. below list if > 0)"
        << "\n    Store to database              = " << (theStoreToDB?"true":"false")
        << "\n    ----------------------------------------------- ----------"
        << (theSetNominalStrips?"\n    Set strips to nominal       =  true":"\n")
        << "\n=============================================================" << std::endl;

  // tell about masked modules
  if( theMaskTecModules.size() ) {
    std::cout << " ===============================================================================================\n" << std::flush;
    std::cout << " The following " << theMaskTecModules.size() << " TEC modules have been masked out and will not be considered by the TEC algorithm:\n " << std::flush;
    for( std::vector<unsigned int>::iterator moduleIt = theMaskTecModules.begin(); moduleIt != theMaskTecModules.end(); ++moduleIt ) {
      std::cout << *moduleIt << (moduleIt!=--theMaskTecModules.end()?", ":"") << std::flush;
    }
    std::cout << std::endl << std::flush;
    std::cout << " ===============================================================================================\n\n" << std::flush;
  }
  if( theMaskAtModules.size() ) {
    std::cout << " ===============================================================================================\n" << std::flush;
    std::cout << " The following " << theMaskAtModules.size() << " AT modules have been masked out and will not be considered by the AT algorithm:\n " << std::flush;
    for( std::vector<unsigned int>::iterator moduleIt = theMaskAtModules.begin(); moduleIt != theMaskAtModules.end(); ++moduleIt ) {
      std::cout << *moduleIt << (moduleIt!=--theMaskAtModules.end()?", ":"") << std::flush;
    }
    std::cout << std::endl << std::flush;
    std::cout << " ===============================================================================================\n\n" << std::flush;
  }
  


  // alias for the Branches in the root files
  std::string alias ( theConf.getParameter<std::string>("@module_label") );  

  // declare the product to produce
  produces<TkLasBeamCollection, edm::Transition::EndRun>( "tkLaserBeams" ).setBranchAlias( alias + "TkLasBeamCollection" );

  // switch judge's zero filter depending on cfg
  judge.EnableZeroFilter( enableJudgeZeroFilter );

  // set the upper threshold for zero suppressed data
  judge.SetOverdriveThreshold( judgeOverdriveThreshold );

}





LaserAlignment::~LaserAlignment() {

  if ( theSaveHistograms ) theFile->Write();
  if ( theFile ) { delete theFile; }
  if ( theAlignableTracker ) { delete theAlignableTracker; }

}





void LaserAlignment::beginJob() {


  // write sumed histograms to file (if selected in cfg)
  if( theSaveHistograms ) {

    // creating a new file
    theFile = new TFile( theFileName.c_str(), "RECREATE", "CMS ROOT file" );
    
    // initialize the histograms
    if ( theFile ) {
      theFile->SetCompressionLevel(theCompression);
      singleModulesDir = theFile->mkdir( "single modules" );
    } else 
      throw cms::Exception( " [LaserAlignment::beginJob]") << " ** ERROR: could not open file:"
                               << theFileName.c_str() << " for writing." << std::endl;

  }

  // detector id maps (hard coded)
  fillDetectorId();

  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  //   PROFILE, HISTOGRAM & FITFUNCTION INITIALIZATION
  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  // object used to build various strings for names and labels
  std::stringstream nameBuilder;

  // loop variables for use with LASGlobalLoop object
  int det, ring, beam, disk, pos;

  // loop TEC modules
  det = 0; ring = 0; beam = 0; disk = 0;
  do { // loop using LASGlobalLoop functionality
    // init the profiles
    pedestalProfiles.GetTECEntry( det, ring, beam, disk ).SetAllValuesTo( 0. );
    currentDataProfiles.GetTECEntry( det, ring, beam, disk ).SetAllValuesTo( 0. );
    collectedDataProfiles.GetTECEntry( det, ring, beam, disk ).SetAllValuesTo( 0. );

    // init the hit maps
    isAcceptedProfile.SetTECEntry( det, ring, beam, disk, 0 );
    numberOfAcceptedProfiles.SetTECEntry( det, ring, beam, disk, 0 );

    // create strings for histo names
    nameBuilder.clear();
    nameBuilder.str( "" );
    nameBuilder << "TEC";
    if( det == 0 ) nameBuilder << "+"; else nameBuilder << "-";
    nameBuilder << "_Ring";
    if( ring == 0 ) nameBuilder << "4"; else nameBuilder << "6";
    nameBuilder << "_Beam" << beam;
    nameBuilder << "_Disk" << disk;
    theProfileNames.SetTECEntry( det, ring, beam, disk, nameBuilder.str() );

    // init the histograms
    if( theSaveHistograms ) {
      nameBuilder << "_Histo";
      summedHistograms.SetTECEntry( det, ring, beam, disk, new TH1D( nameBuilder.str().c_str(), nameBuilder.str().c_str(), 512, 0, 512 ) );
      summedHistograms.GetTECEntry( det, ring, beam, disk )->SetDirectory( singleModulesDir );
    }
    
  } while ( moduleLoop.TECLoop( det, ring, beam, disk ) );


  // TIB & TOB section
  det = 2; beam = 0; pos = 0;
  do { // loop using LASGlobalLoop functionality
    // init the profiles
    pedestalProfiles.GetTIBTOBEntry( det, beam, pos ).SetAllValuesTo( 0. );
    currentDataProfiles.GetTIBTOBEntry( det, beam, pos ).SetAllValuesTo( 0. );
    collectedDataProfiles.GetTIBTOBEntry( det, beam, pos ).SetAllValuesTo( 0. );

    // init the hit maps
    isAcceptedProfile.SetTIBTOBEntry( det, beam, pos, 0 );
    numberOfAcceptedProfiles.SetTIBTOBEntry( det, beam, pos, 0 );

    // create strings for histo names
    nameBuilder.clear();
    nameBuilder.str( "" );
    if( det == 2 ) nameBuilder << "TIB"; else nameBuilder << "TOB";
    nameBuilder << "_Beam" << beam;
    nameBuilder << "_Zpos" << pos;

    theProfileNames.SetTIBTOBEntry( det, beam, pos, nameBuilder.str() );

    // init the histograms
    if( theSaveHistograms ) {
      nameBuilder << "_Histo";
      summedHistograms.SetTIBTOBEntry( det, beam, pos, new TH1D( nameBuilder.str().c_str(), nameBuilder.str().c_str(), 512, 0, 512 ) );
      summedHistograms.GetTIBTOBEntry( det, beam, pos )->SetDirectory( singleModulesDir );
    }
    
  } while( moduleLoop.TIBTOBLoop( det, beam, pos ) );


  // TEC2TEC AT section
  det = 0; beam = 0; disk = 0;
  do { // loop using LASGlobalLoop functionality
    // init the profiles
    pedestalProfiles.GetTEC2TECEntry( det, beam, disk ).SetAllValuesTo( 0. );
    currentDataProfiles.GetTEC2TECEntry( det, beam, disk ).SetAllValuesTo( 0. );
    collectedDataProfiles.GetTEC2TECEntry( det, beam, disk ).SetAllValuesTo( 0. );
    
    // init the hit maps
    isAcceptedProfile.SetTEC2TECEntry( det, beam, disk, 0 );
    numberOfAcceptedProfiles.SetTEC2TECEntry( det, beam, disk, 0 );

    // create strings for histo names
    nameBuilder.clear();
    nameBuilder.str( "" );
    nameBuilder << "TEC(AT)";
    if( det == 0 ) nameBuilder << "+"; else nameBuilder << "-";
    nameBuilder << "_Beam" << beam;
    nameBuilder << "_Disk" << disk;
    theProfileNames.SetTEC2TECEntry( det, beam, disk, nameBuilder.str() );

    // init the histograms
    if( theSaveHistograms ) {
      nameBuilder << "_Histo";
      summedHistograms.SetTEC2TECEntry( det, beam, disk, new TH1D( nameBuilder.str().c_str(), nameBuilder.str().c_str(), 512, 0, 512 ) );
      summedHistograms.GetTEC2TECEntry( det, beam, disk )->SetDirectory( singleModulesDir );
    }
    
  } while( moduleLoop.TEC2TECLoop( det, beam, disk ) );

  firstEvent_ = true;
}






void LaserAlignment::produce(edm::Event& theEvent, edm::EventSetup const& theSetup)  {

  if (firstEvent_) {

    //Retrieve tracker topology from geometry
    edm::ESHandle<TrackerTopology> tTopoHandle;
    theSetup.get<TrackerTopologyRcd>().get(tTopoHandle);
    const TrackerTopology* const tTopo = tTopoHandle.product();

    // access the tracker
    theSetup.get<TrackerDigiGeometryRecord>().get( theTrackerGeometry );
    theSetup.get<IdealGeometryRecord>().get( gD );
    
    // access pedestals (from db..) if desired
    edm::ESHandle<SiStripPedestals> pedestalsHandle;
    if( theDoPedestalSubtraction ) {
      theSetup.get<SiStripPedestalsRcd>().get( pedestalsHandle );
      fillPedestalProfiles( pedestalsHandle );
    }
    
    // global positions
    //  edm::ESHandle<Alignments> theGlobalPositionRcd;
    theSetup.get<TrackerDigiGeometryRecord>().getRecord<GlobalPositionRcd>().get( theGlobalPositionRcd );

    // select the reference geometry
    if( !updateFromInputGeometry ) {
      // the AlignableTracker object is initialized with the ideal geometry
      edm::ESHandle<GeometricDet> theGeometricDet;
      theSetup.get<IdealGeometryRecord>().get(theGeometricDet);
      edm::ESHandle<PTrackerParameters> ptp;
      theSetup.get<PTrackerParametersRcd>().get( ptp );
      TrackerGeomBuilderFromGeometricDet trackerBuilder;
      TrackerGeometry* theRefTracker = trackerBuilder.build(&*theGeometricDet, *ptp, tTopo );
      
      theAlignableTracker = new AlignableTracker(&(*theRefTracker), tTopo);
    }
    else {
      // the AlignableTracker object is initialized with the input geometry from DB
      theAlignableTracker = new AlignableTracker( &(*theTrackerGeometry), tTopo );
    }
    
    firstEvent_ = false;
  }

  LogDebug("LaserAlignment") << "==========================================================="
                  << "\n   Private analysis of event #"<< theEvent.id().event() 
                  << " in run #" << theEvent.id().run();


  // do the Tracker Statistics to retrieve the current profiles
  fillDataProfiles( theEvent, theSetup );

  // index variables for the LASGlobalLoop object
  int det, ring, beam, disk, pos;

  //
  // first pre-loop on selected entries to find out
  // whether the TEC or the AT beams have fired
  // (pedestal profiles are left empty if false in cfg)
  // 


  // TEC+- (only ring 6)
  ring = 1;
  for( det = 0; det < 2; ++det ) {
    for( beam = 0; beam < 8; ++ beam ) {
      for( disk = 0; disk < 9; ++disk ) {
    if( judge.IsSignalIn( currentDataProfiles.GetTECEntry( det, ring, beam, disk ) - pedestalProfiles.GetTECEntry( det, ring, beam, disk ), 0 ) ) {
      isAcceptedProfile.SetTECEntry( det, ring, beam, disk, 1 );
    }
    else { // assume no initialization
      isAcceptedProfile.SetTECEntry( det, ring, beam, disk, 0 );
    }
      }
    }
  }

  // TIBTOB
  det = 2; beam = 0; pos = 0;
  do {
    // add current event's data and subtract pedestals
    if( judge.IsSignalIn( currentDataProfiles.GetTIBTOBEntry( det, beam, pos ) - pedestalProfiles.GetTIBTOBEntry( det, beam, pos ), getTIBTOBNominalBeamOffset( det, beam, pos ) ) ) {
      isAcceptedProfile.SetTIBTOBEntry( det, beam, pos, 1 );
    }
    else { // dto.
      isAcceptedProfile.SetTIBTOBEntry( det, beam, pos, 0 );
    }

  } while( moduleLoop.TIBTOBLoop( det, beam, pos ) );




  // now come the beam finders
  bool isTECMode = isTECBeam();
  //  LogDebug( " [LaserAlignment::produce]" ) << "LaserAlignment::isTECBeam declares this event " << ( isTECMode ? "" : "NOT " ) << "a TEC event." << std::endl;
  std::cout << " [LaserAlignment::produce] -- LaserAlignment::isTECBeam declares this event " << ( isTECMode ? "" : "NOT " ) << "a TEC event." << std::endl;

  bool isATMode  = isATBeam();
  //  LogDebug( " [LaserAlignment::produce]" ) << "LaserAlignment::isATBeam declares this event "  << ( isATMode ? "" : "NOT " )  << "an AT event." << std::endl;
  std::cout << " [LaserAlignment::produce] -- LaserAlignment::isATBeam declares this event "  << ( isATMode ? "" : "NOT " )  << "an AT event." << std::endl;




  //
  // now pass the pedestal subtracted profiles to the judge
  // if they're accepted, add them on the collectedDataProfiles
  // (pedestal profiles are left empty if false in cfg)
  //


  // loop TEC+- modules
  det = 0; ring = 0; beam = 0; disk = 0;
  do {
    
    LogDebug( "[LaserAlignment::produce]" ) << "Profile is: " << theProfileNames.GetTECEntry( det, ring, beam, disk ) << "." << std::endl;

    // this now depends on the AT/TEC mode, is this a doubly hit module? -> look for it in vector<int> tecDoubleHitDetId
    // (ring == 0 not necessary but makes it a little faster)
    if( ring == 0 && find( tecDoubleHitDetId.begin(), tecDoubleHitDetId.end(), detectorId.GetTECEntry( det, ring, beam, disk ) ) != tecDoubleHitDetId.end() ) {

      if( isTECMode ) { // add profile to TEC collection
    // add current event's data and subtract pedestals
    if( judge.JudgeProfile( currentDataProfiles.GetTECEntry( det, ring, beam, disk ) - pedestalProfiles.GetTECEntry( det, ring, beam, disk ), 0 ) ) {
      collectedDataProfiles.GetTECEntry( det, ring, beam, disk ) += currentDataProfiles.GetTECEntry( det, ring, beam, disk ) - pedestalProfiles.GetTECEntry( det, ring, beam, disk );
      numberOfAcceptedProfiles.GetTECEntry( det, ring, beam, disk )++;
    }
      }
    }

    else { // not a doubly hit module, don't care about the mode
      // add current event's data and subtract pedestals
      if( judge.JudgeProfile( currentDataProfiles.GetTECEntry( det, ring, beam, disk ) - pedestalProfiles.GetTECEntry( det, ring, beam, disk ), 0 ) ) {
    collectedDataProfiles.GetTECEntry( det, ring, beam, disk ) += currentDataProfiles.GetTECEntry( det, ring, beam, disk ) - pedestalProfiles.GetTECEntry( det, ring, beam, disk );
    numberOfAcceptedProfiles.GetTECEntry( det, ring, beam, disk )++;
      }
    }
    
  } while( moduleLoop.TECLoop( det, ring, beam, disk ) );


  


  // loop TIB/TOB modules
  det = 2; beam = 0; pos = 0;
  do {
    
    LogDebug( "[LaserAlignment::produce]" ) << "Profile is: " << theProfileNames.GetTIBTOBEntry( det, beam, pos ) << "." << std::endl;
    
    // add current event's data and subtract pedestals
    if( judge.JudgeProfile( currentDataProfiles.GetTIBTOBEntry( det, beam, pos ) - pedestalProfiles.GetTIBTOBEntry( det, beam, pos ), getTIBTOBNominalBeamOffset( det, beam, pos ) ) ) {
      collectedDataProfiles.GetTIBTOBEntry( det, beam, pos ) += currentDataProfiles.GetTIBTOBEntry( det, beam, pos ) - pedestalProfiles.GetTIBTOBEntry( det, beam, pos );
      numberOfAcceptedProfiles.GetTIBTOBEntry( det, beam, pos )++;
    }

  } while( moduleLoop.TIBTOBLoop( det, beam, pos ) );
  


  // loop TEC2TEC modules
  det = 0; beam = 0; disk = 0;
  do {
    
    LogDebug( "[LaserAlignment::produce]" ) << "Profile is: " << theProfileNames.GetTEC2TECEntry( det, beam, disk ) << "." << std::endl;

    // this again depends on the AT/TEC mode, is this a doubly hit module?
    // (ring == 0 not necessary but makes it a little faster)
    if( ring == 0 && find( tecDoubleHitDetId.begin(), tecDoubleHitDetId.end(), detectorId.GetTECEntry( det, ring, beam, disk ) ) != tecDoubleHitDetId.end() ) {

      if( isATMode ) { // add profile to TEC2TEC collection
    // add current event's data and subtract pedestals
    if( judge.JudgeProfile( currentDataProfiles.GetTEC2TECEntry( det, beam, disk ) - pedestalProfiles.GetTEC2TECEntry( det, beam, disk ), 0 ) ) {
      collectedDataProfiles.GetTEC2TECEntry( det, beam, disk ) += currentDataProfiles.GetTEC2TECEntry( det, beam, disk ) - pedestalProfiles.GetTEC2TECEntry( det, beam, disk );
      numberOfAcceptedProfiles.GetTEC2TECEntry( det, beam, disk )++;
    }
      }

    }     
    
    else { // not a doubly hit module, don't care about the mode
      // add current event's data and subtract pedestals
      if( judge.JudgeProfile( currentDataProfiles.GetTEC2TECEntry( det, beam, disk ) - pedestalProfiles.GetTEC2TECEntry( det, beam, disk ), 0 ) ) {
    collectedDataProfiles.GetTEC2TECEntry( det, beam, disk ) += currentDataProfiles.GetTEC2TECEntry( det, beam, disk ) - pedestalProfiles.GetTEC2TECEntry( det, beam, disk );
    numberOfAcceptedProfiles.GetTEC2TECEntry( det, beam, disk )++;
      }
    }
      

  } while( moduleLoop.TEC2TECLoop( det, beam, disk ) );



  // total event number counter
  theEvents++;

}





void LaserAlignment::endRunProduce( edm::Run& theRun, const edm::EventSetup& theSetup ) {


  std::cout << " [LaserAlignment::endRun] -- Total number of events processed: " << theEvents << std::endl;

  // for debugging only..
  DumpHitmaps( numberOfAcceptedProfiles );

  // index variables for the LASGlobalLoop objects
  int det, ring, beam, disk, pos;
    
  // measured positions container for the algorithms
  LASGlobalData<LASCoordinateSet> measuredCoordinates;
  
  // fitted peak positions in units of strips (pair for value,error)
  LASGlobalData<std::pair<float,float> > measuredStripPositions;
  
  // the peak finder, a pair (pos/posErr in units of strips) for its results, and the success confirmation
  LASPeakFinder peakFinder;
  peakFinder.SetAmplitudeThreshold( peakFinderThreshold );
  std::pair<double,double> peakFinderResults;
  bool isGoodFit;
  
  // tracker geom. object for calculating the global beam positions
  const TrackerGeometry& theTracker( *theTrackerGeometry );

  // fill LASGlobalData<LASCoordinateSet> nominalCoordinates
  CalculateNominalCoordinates();
  
  // for determining the phi errors
  //    ErrorFrameTransformer errorTransformer; // later...
  




  // do the fits for TEC+- internal
  det = 0; ring = 0; beam = 0; disk = 0;
  do {
    
    // do the fit
    isGoodFit = peakFinder.FindPeakIn( collectedDataProfiles.GetTECEntry( det, ring, beam, disk ), peakFinderResults,
                       summedHistograms.GetTECEntry( det, ring, beam, disk ), 0 ); // offset is 0 for TEC

    // now we have the measured positions in units of strips. 
    if( !isGoodFit ) std::cout << " [LaserAlignment::endRun] ** WARNING: Fit failed for TEC det: "
                   << det << ", ring: " << ring << ", beam: " << beam << ", disk: " << disk
                   << " (id: " << detectorId.GetTECEntry( det, ring, beam, disk ) << ")." << std::endl;

    

    // <- here we will later implement the kink corrections
      
    // access the tracker geometry for this module
    const DetId theDetId( detectorId.GetTECEntry( det, ring, beam, disk ) );
    const StripGeomDetUnit* const theStripDet = dynamic_cast<const StripGeomDetUnit*>( theTracker.idToDet( theDetId ) );
      
    if (theStripDet) {
      // first, set the measured coordinates to their nominal values
      measuredCoordinates.SetTECEntry( det, ring, beam, disk, nominalCoordinates.GetTECEntry( det, ring, beam, disk ) );
      
      if( isGoodFit ) { // convert strip position to global phi and replace the nominal phi value/error
    
    measuredStripPositions.GetTECEntry( det, ring, beam, disk ) = peakFinderResults;
    const float positionInStrips =  theSetNominalStrips ? 256. : peakFinderResults.first; // implementation of "ForceFitterToNominalStrips" config parameter
    const GlobalPoint& globalPoint = theStripDet->surface().toGlobal( theStripDet->specificTopology().localPosition( positionInStrips ) );
    measuredCoordinates.GetTECEntry( det, ring, beam, disk ).SetPhi( ConvertAngle( globalPoint.barePhi() ) );
    
    // const GlobalError& globalError = errorTransformer.transform( theStripDet->specificTopology().localError( peakFinderResults.first, pow( peakFinderResults.second, 2 ) ), theStripDet->surface() );
    // measuredCoordinates.GetTECEntry( det, ring, beam, disk ).SetPhiError( globalError.phierr( globalPoint ) );
    measuredCoordinates.GetTECEntry( det, ring, beam, disk ).SetPhiError( 0.00046  ); // PRELIMINARY ESTIMATE
    
      }
      else { // keep nominal position (middle-of-module) but set a giant phi error so that the module can be ignored by the alignment algorithm
    measuredStripPositions.GetTECEntry( det, ring, beam, disk ) = std::pair<float,float>( 256., 1000. );
    const GlobalPoint& globalPoint = theStripDet->surface().toGlobal( theStripDet->specificTopology().localPosition( 256. ) );
    measuredCoordinates.GetTECEntry( det, ring, beam, disk ).SetPhi( ConvertAngle( globalPoint.barePhi() ) );
    measuredCoordinates.GetTECEntry( det, ring, beam, disk ).SetPhiError( 1000. );
      }
    }
      
  } while( moduleLoop.TECLoop( det, ring, beam, disk ) );




  // do the fits for TIB/TOB
  det = 2; beam = 0; pos = 0;
  do {

    // do the fit
    isGoodFit = peakFinder.FindPeakIn( collectedDataProfiles.GetTIBTOBEntry( det, beam, pos ), peakFinderResults, 
                       summedHistograms.GetTIBTOBEntry( det, beam, pos ), getTIBTOBNominalBeamOffset( det, beam, pos ) );

    // now we have the measured positions in units of strips.
    if( !isGoodFit ) std::cout << " [LaserAlignment::endJob] ** WARNING: Fit failed for TIB/TOB det: "
                   << det << ", beam: " << beam << ", pos: " << pos 
                   << " (id: " << detectorId.GetTIBTOBEntry( det, beam, pos ) << ")." << std::endl;

      
    // <- here we will later implement the kink corrections
      
    // access the tracker geometry for this module
    const DetId theDetId( detectorId.GetTIBTOBEntry( det, beam, pos ) );
    const StripGeomDetUnit* const theStripDet = dynamic_cast<const StripGeomDetUnit*>( theTracker.idToDet( theDetId ) );
      
    if (theStripDet) {
      // first, set the measured coordinates to their nominal values
      measuredCoordinates.SetTIBTOBEntry( det, beam, pos, nominalCoordinates.GetTIBTOBEntry( det, beam, pos ) );
      
      if( isGoodFit ) { // convert strip position to global phi and replace the nominal phi value/error
    measuredStripPositions.GetTIBTOBEntry( det, beam, pos ) = peakFinderResults;
    const float positionInStrips =  theSetNominalStrips ? 256. + getTIBTOBNominalBeamOffset( det, beam, pos ) : peakFinderResults.first; // implementation of "ForceFitterToNominalStrips" config parameter
    const GlobalPoint& globalPoint = theStripDet->surface().toGlobal( theStripDet->specificTopology().localPosition( positionInStrips ) );
    measuredCoordinates.GetTIBTOBEntry( det, beam, pos ).SetPhi( ConvertAngle( globalPoint.barePhi() ) );
    measuredCoordinates.GetTIBTOBEntry( det, beam, pos ).SetPhiError( 0.00028 ); // PRELIMINARY ESTIMATE
      }
      else { // keep nominal position but set a giant phi error so that the module can be ignored by the alignment algorithm
    measuredStripPositions.GetTIBTOBEntry( det, beam, pos ) = std::pair<float,float>( 256. + getTIBTOBNominalBeamOffset( det, beam, pos ), 1000. );
    const GlobalPoint& globalPoint = theStripDet->surface().toGlobal( theStripDet->specificTopology().localPosition( 256. + getTIBTOBNominalBeamOffset( det, beam, pos ) ) );
    measuredCoordinates.GetTIBTOBEntry( det, beam, pos ).SetPhi( ConvertAngle( globalPoint.barePhi() ) );
    measuredCoordinates.GetTIBTOBEntry( det, beam, pos ).SetPhiError( 1000. );
      }
    }

  } while( moduleLoop.TIBTOBLoop( det, beam, pos ) );




  // do the fits for TEC AT
  det = 0; beam = 0; disk = 0;
  do {

    // do the fit
    isGoodFit = peakFinder.FindPeakIn( collectedDataProfiles.GetTEC2TECEntry( det, beam, disk ), peakFinderResults,
                       summedHistograms.GetTEC2TECEntry( det, beam, disk ), getTEC2TECNominalBeamOffset( det, beam, disk ) );
    // now we have the positions in units of strips.
    if( !isGoodFit ) std::cout << " [LaserAlignment::endRun] ** WARNING: Fit failed for TEC2TEC det: "
                   << det << ", beam: " << beam << ", disk: " << disk
                   << " (id: " << detectorId.GetTEC2TECEntry( det, beam, disk ) << ")." << std::endl;


    // <- here we will later implement the kink corrections
    
    // access the tracker geometry for this module
    const DetId theDetId( detectorId.GetTEC2TECEntry( det, beam, disk ) );
    const StripGeomDetUnit* const theStripDet = dynamic_cast<const StripGeomDetUnit*>( theTracker.idToDet( theDetId ) );

    if (theStripDet) {
      // first, set the measured coordinates to their nominal values
      measuredCoordinates.SetTEC2TECEntry( det, beam, disk, nominalCoordinates.GetTEC2TECEntry( det, beam, disk ) );
      
      if( isGoodFit ) { // convert strip position to global phi and replace the nominal phi value/error
    measuredStripPositions.GetTEC2TECEntry( det, beam, disk ) = peakFinderResults;
    const float positionInStrips =  theSetNominalStrips ? 256. + getTEC2TECNominalBeamOffset( det, beam, disk ) : peakFinderResults.first; // implementation of "ForceFitterToNominalStrips" config parameter
    const GlobalPoint& globalPoint = theStripDet->surface().toGlobal( theStripDet->specificTopology().localPosition( positionInStrips ) );
    measuredCoordinates.GetTEC2TECEntry( det, beam, disk ).SetPhi( ConvertAngle( globalPoint.barePhi() ) );
    measuredCoordinates.GetTEC2TECEntry( det, beam, disk ).SetPhiError( 0.00047 ); // PRELIMINARY ESTIMATE
      }
      else { // keep nominal position but set a giant phi error so that the module can be ignored by the alignment algorithm
    measuredStripPositions.GetTEC2TECEntry( det, beam, disk ) = std::pair<float,float>( 256. + getTEC2TECNominalBeamOffset( det, beam, disk ), 1000. );
    const GlobalPoint& globalPoint = theStripDet->surface().toGlobal( theStripDet->specificTopology().localPosition( 256. + getTEC2TECNominalBeamOffset( det, beam, disk ) ) );
    measuredCoordinates.GetTEC2TECEntry( det, beam, disk ).SetPhi( ConvertAngle( globalPoint.barePhi() ) );
    measuredCoordinates.GetTEC2TECEntry( det, beam, disk ).SetPhiError( 1000. );
      }
    }

  } while( moduleLoop.TEC2TECLoop( det, beam, disk ) );
  






  // see what we got (for debugging)
  //  DumpStripFileSet( measuredStripPositions );
  //  DumpPosFileSet( measuredCoordinates );


  



  // CALCULATE PARAMETERS AND UPDATE DB OBJECT
  // for beam kink corrections, reconstructing the geometry and updating the db object
  LASGeometryUpdater geometryUpdater( nominalCoordinates, theLasConstants );

  // apply all beam corrections
  if( theApplyBeamKinkCorrections ) geometryUpdater.ApplyBeamKinkCorrections( measuredCoordinates );

  // if we start with input geometry instead of IDEAL,
  // reverse the adjustments in the AlignableTracker object
  if( updateFromInputGeometry ) geometryUpdater.SetReverseDirection( true );

  // if we have "virtual" misalignment which is introduced via the reference geometry,
  // tell the LASGeometryUpdater to reverse x & y adjustments
  if( misalignedByRefGeometry ) geometryUpdater.SetMisalignmentFromRefGeometry( true );

  // run the endcap algorithm
  LASEndcapAlgorithm endcapAlgorithm;
  LASEndcapAlignmentParameterSet endcapParameters;


  // this basically sets all the endcap modules to be masked 
  // to their nominal positions (since endcapParameters is overall zero)
  if( theMaskTecModules.size() ) {
    ApplyEndcapMaskingCorrections( measuredCoordinates, nominalCoordinates, endcapParameters );
  }

  // run the algorithm
  endcapParameters = endcapAlgorithm.CalculateParameters( measuredCoordinates, nominalCoordinates );

  // 
  // loop to mask out events
  // DESCRIPTION:
  //

  // do this only if there are modules to be masked..
  if( theMaskTecModules.size() ) {
    
    const unsigned int nIterations = 30;
    for( unsigned int iteration = 0; iteration < nIterations; ++iteration ) {
      
      // set the endcap modules to be masked to their positions
      // according to the reconstructed parameters
      ApplyEndcapMaskingCorrections( measuredCoordinates, nominalCoordinates, endcapParameters );
      
      // modifications applied, so re-run the algorithm
      endcapParameters = endcapAlgorithm.CalculateParameters( measuredCoordinates, nominalCoordinates );
      
    }

  } 

  // these are now final, so:
  endcapParameters.Print();



  
  // do a pre-alignment of the endcaps (TEC2TEC only)
  // so that the alignment tube algorithms finds orderly disks
  geometryUpdater.EndcapUpdate( endcapParameters, measuredCoordinates );


  // the alignment tube algorithms, choose from config
  LASBarrelAlignmentParameterSet alignmentTubeParameters;
  // the MINUIT-BASED alignment tube algorithm
  LASBarrelAlgorithm barrelAlgorithm;
  // the ANALYTICAL alignment tube algorithm
  LASAlignmentTubeAlgorithm alignmentTubeAlgorithm;


  // this basically sets all the modules to be masked 
  // to their nominal positions (since alignmentTubeParameters is overall zero)
  if( theMaskAtModules.size() ) {
    ApplyATMaskingCorrections( measuredCoordinates, nominalCoordinates, alignmentTubeParameters );
  }

  if( theUseMinuitAlgorithm ) {
    // run the MINUIT-BASED alignment tube algorithm
    alignmentTubeParameters = barrelAlgorithm.CalculateParameters( measuredCoordinates, nominalCoordinates );
  }
  else {
    // the ANALYTICAL alignment tube algorithm
    alignmentTubeParameters = alignmentTubeAlgorithm.CalculateParameters( measuredCoordinates, nominalCoordinates );
  }



  // 
  // loop to mask out events
  // DESCRIPTION:
  //

  // do this only if there are modules to be masked..
  if( theMaskAtModules.size() ) {
    
    const unsigned int nIterations = 30;
    for( unsigned int iteration = 0; iteration < nIterations; ++iteration ) {
      
      // set the AT modules to be masked to their positions
      // according to the reconstructed parameters
      ApplyATMaskingCorrections( measuredCoordinates, nominalCoordinates, alignmentTubeParameters );
      
      // modifications applied, so re-run the algorithm
      if( theUseMinuitAlgorithm ) {
    alignmentTubeParameters = barrelAlgorithm.CalculateParameters( measuredCoordinates, nominalCoordinates );
      }
      else {
    alignmentTubeParameters = alignmentTubeAlgorithm.CalculateParameters( measuredCoordinates, nominalCoordinates );
      }
      
    }

  } 


  // these are now final, so:
  alignmentTubeParameters.Print();


  // combine the results and update the db object
  geometryUpdater.TrackerUpdate( endcapParameters, alignmentTubeParameters, *theAlignableTracker );
  




    

    
    
  // the collection container
  auto laserBeams = std::make_unique<TkLasBeamCollection>();

  
  // first for the endcap internal beams
  for( det = 0; det < 2; ++det ) {
    for( ring = 0; ring < 2; ++ring ) {
      for( beam = 0; beam < 8; ++beam ) {
    
    // the beam and its identifier (see TkLasTrackBasedInterface TWiki)
    TkLasBeam currentBeam( 100 * det + 10 * beam + ring );
    
    // order the hits in the beam by increasing z
    const int firstDisk = det==0 ? 0 : 8;
    const int lastDisk  = det==0 ? 8 : 0;
      
    // count upwards or downwards
    for( disk = firstDisk; det==0 ? disk <= lastDisk : disk >= lastDisk; det==0 ? ++disk : --disk ) {
        
      // detId for the SiStripLaserRecHit2D
      const SiStripDetId theDetId( detectorId.GetTECEntry( det, ring, beam, disk ) );
        
      // need this to calculate the localPosition and its error
      const StripGeomDetUnit* const theStripDet = dynamic_cast<const StripGeomDetUnit*>( theTracker.idToDet( theDetId ) );
        
      // the hit container
      const SiStripLaserRecHit2D currentHit(
        theStripDet->specificTopology().localPosition( measuredStripPositions.GetTECEntry( det, ring, beam, disk ).first ),
        theStripDet->specificTopology().localError( measuredStripPositions.GetTECEntry( det, ring, beam, disk ).first, measuredStripPositions.GetTECEntry( det, ring, beam, disk ).second ),
        theDetId
      );

      currentBeam.push_back( currentHit );

    }     
      
    laserBeams->push_back( currentBeam );
      
      }
    }
  }
    
    

  // then, following the convention in TkLasTrackBasedInterface TWiki, the alignment tube beams;
  // they comprise hits in TIBTOB & TEC2TEC

  for( beam = 0; beam < 8; ++beam ) {

    // the beam and its identifier (see TkLasTrackBasedInterface TWiki)
    TkLasBeam currentBeam( 100 * 2 /*beamGroup=AT=2*/ + 10 * beam + 0 /*ring=0*/);


    // first: tec-
    det = 1;
    for( disk = 4; disk >= 0; --disk ) {
    
      // detId for the SiStripLaserRecHit2D
      const SiStripDetId theDetId( detectorId.GetTEC2TECEntry( det, beam, disk ) );
    
      // need this to calculate the localPosition and its error
      const StripGeomDetUnit* const theStripDet = dynamic_cast<const StripGeomDetUnit*>( theTracker.idToDet( theDetId ) );
    
      // the hit container
      const SiStripLaserRecHit2D currentHit(
        theStripDet->specificTopology().localPosition( measuredStripPositions.GetTEC2TECEntry( det, beam, disk ).first ),
    theStripDet->specificTopology().localError( measuredStripPositions.GetTEC2TECEntry( det, beam, disk ).first, measuredStripPositions.GetTEC2TECEntry( det, beam, disk ).second ),
    theDetId
      );

      currentBeam.push_back( currentHit );
    
    }

      
    // now TIB and TOB in one go
    for( det = 2; det < 4; ++det ) {
      for( pos = 5; pos >= 0; --pos ) { // stupidly, pos is defined from +z to -z in LASGlobalLoop
      
    // detId for the SiStripLaserRecHit2D
    const SiStripDetId theDetId( detectorId.GetTIBTOBEntry( det, beam, pos ) );
      
    // need this to calculate the localPosition and its error
    const StripGeomDetUnit* const theStripDet = dynamic_cast<const StripGeomDetUnit*>( theTracker.idToDet( theDetId ) );
      
    // the hit container
    const SiStripLaserRecHit2D currentHit(
      theStripDet->specificTopology().localPosition( measuredStripPositions.GetTIBTOBEntry( det, beam, pos ).first ),
      theStripDet->specificTopology().localError( measuredStripPositions.GetTIBTOBEntry( det, beam, pos ).first, measuredStripPositions.GetTIBTOBEntry( det, beam, pos ).second ),
      theDetId
    );

    currentBeam.push_back( currentHit );
      
      }
    }
      

    // then: tec+
    det = 0;
    for( disk = 0; disk < 5; ++disk ) {
    
      // detId for the SiStripLaserRecHit2D
      const SiStripDetId theDetId( detectorId.GetTEC2TECEntry( det, beam, disk ) );
    
      // need this to calculate the localPosition and its error
      const StripGeomDetUnit* const theStripDet = dynamic_cast<const StripGeomDetUnit*>( theTracker.idToDet( theDetId ) );
    
      // the hit container
      const SiStripLaserRecHit2D currentHit(
        theStripDet->specificTopology().localPosition( measuredStripPositions.GetTEC2TECEntry( det, beam, disk ).first ),
    theStripDet->specificTopology().localError( measuredStripPositions.GetTEC2TECEntry( det, beam, disk ).first, measuredStripPositions.GetTEC2TECEntry( det, beam, disk ).second ),
    theDetId
      );

      currentBeam.push_back( currentHit );
    
    }



    // save this beam to the beamCollection
    laserBeams->push_back( currentBeam );
    
  } // (close beam loop)
  
  
  // now attach the collection to the run
  theRun.put(std::move(laserBeams), "tkLaserBeams" );
  


  

  // store the estimated alignment parameters into the DB
  // first get them
  Alignments* alignments =  theAlignableTracker->alignments();
  AlignmentErrorsExtended* alignmentErrors = theAlignableTracker->alignmentErrors();

  if ( theStoreToDB ) {

    std::cout << " [LaserAlignment::endRun] -- Storing the calculated alignment parameters to the DataBase:" << std::endl;

    // Call service
    edm::Service<cond::service::PoolDBOutputService> poolDbService;
    if( !poolDbService.isAvailable() ) // Die if not available
      throw cms::Exception( "NotAvailable" ) << "PoolDBOutputService not available";
    
    // Store

    //     if ( poolDbService->isNewTagRequest(theAlignRecordName) ) {
    //       poolDbService->createNewIOV<Alignments>( alignments, poolDbService->currentTime(), poolDbService->endOfTime(), theAlignRecordName );
    //     }
    //     else {
    //       poolDbService->appendSinceTime<Alignments>( alignments, poolDbService->currentTime(), theAlignRecordName );
    //     }
    poolDbService->writeOne<Alignments>( alignments, poolDbService->beginOfTime(), theAlignRecordName );

    //     if ( poolDbService->isNewTagRequest(theErrorRecordName) ) {
    //       poolDbService->createNewIOV<AlignmentErrorsExtended>( alignmentErrors, poolDbService->currentTime(), poolDbService->endOfTime(), theErrorRecordName );
    //     }
    //     else {
    //       poolDbService->appendSinceTime<AlignmentErrorsExtended>( alignmentErrors, poolDbService->currentTime(), theErrorRecordName );
    //     }
    poolDbService->writeOne<AlignmentErrorsExtended>( alignmentErrors, poolDbService->beginOfTime(), theErrorRecordName );

    std::cout << " [LaserAlignment::endRun] -- Storing done." << std::endl;
    
  }

}





void LaserAlignment::endJob() {
}





void LaserAlignment::fillDataProfiles( edm::Event const& theEvent, edm::EventSetup const& theSetup ) {

  // two handles for the two different kinds of digis
  edm::Handle< edm::DetSetVector<SiStripRawDigi> > theStripRawDigis;
  edm::Handle< edm::DetSetVector<SiStripDigi> > theStripDigis;

  bool isRawDigi = false;

  // indices for the LASGlobalLoop object
  int det = 0, ring = 0, beam = 0, disk = 0, pos = 0;

  // query config set and loop over all PSets in the VPSet
  for ( std::vector<edm::ParameterSet>::iterator itDigiProducersList = theDigiProducersList.begin(); itDigiProducersList != theDigiProducersList.end(); ++itDigiProducersList ) {

    std::string digiProducer = itDigiProducersList->getParameter<std::string>( "DigiProducer" );
    std::string digiLabel = itDigiProducersList->getParameter<std::string>( "DigiLabel" );
    std::string digiType = itDigiProducersList->getParameter<std::string>( "DigiType" );

    // now branch according to digi type (raw or processed);    
    // first we go for raw digis => SiStripRawDigi
    if( digiType == "Raw" ) {
      theEvent.getByLabel( digiProducer, digiLabel, theStripRawDigis );
      isRawDigi = true;
    }
    else if( digiType == "Processed" ) {
      theEvent.getByLabel( digiProducer, digiLabel, theStripDigis );
      isRawDigi = false;
    }
    else {
      throw cms::Exception( " [LaserAlignment::fillDataProfiles]") << " ** ERROR: Invalid digi type: \"" << digiType << "\" specified in configuration." << std::endl;
    }    



    // loop TEC internal modules
    det = 0; ring = 0; beam = 0; disk = 0;
    do {
      
      // first clear the profile
      currentDataProfiles.GetTECEntry( det, ring, beam, disk ).SetAllValuesTo( 0. );

      // retrieve the raw id of that module
      const int detRawId = detectorId.GetTECEntry( det, ring, beam, disk );
      
      if( isRawDigi ) { // we have raw SiStripRawDigis
    
    // search the digis for the raw id
    edm::DetSetVector<SiStripRawDigi>::const_iterator detSetIter = theStripRawDigis->find( detRawId );
    if( detSetIter == theStripRawDigis->end() ) {
      throw cms::Exception( "[Laser Alignment::fillDataProfiles]" ) << " ** ERROR: No raw DetSet found for det: " << detRawId << "." << std::endl;
    }
      
    // fill the digis to the profiles
    edm::DetSet<SiStripRawDigi>::const_iterator digiRangeIterator = detSetIter->data.begin(); // for the loop
    edm::DetSet<SiStripRawDigi>::const_iterator digiRangeStart = digiRangeIterator; // save starting positions
    
    // loop all digis
    for (; digiRangeIterator != detSetIter->data.end(); ++digiRangeIterator ) {
      const SiStripRawDigi& digi = *digiRangeIterator;
      const int channel = distance( digiRangeStart, digiRangeIterator );
      if ( channel >= 0 && channel < 512 ) currentDataProfiles.GetTECEntry( det, ring, beam, disk ).SetValue( channel, digi.adc() );
      else throw cms::Exception( "[Laser Alignment::fillDataProfiles]" ) << " ** ERROR: raw digi channel: " << channel << " out of range for det: " << detRawId << "." << std::endl;
    }

      }

      else { // we have zero suppressed SiStripDigis

    // search the digis for the raw id
    edm::DetSetVector<SiStripDigi>::const_iterator detSetIter = theStripDigis->find( detRawId );
    
    // processed DetSets may be missing, just skip
    if( detSetIter == theStripDigis->end() ) continue;

    // fill the digis to the profiles
    edm::DetSet<SiStripDigi>::const_iterator digiRangeIterator = detSetIter->data.begin(); // for the loop
    
    for(; digiRangeIterator != detSetIter->data.end(); ++digiRangeIterator ) {
      const SiStripDigi& digi = *digiRangeIterator;
      if ( digi.strip() < 512 ) currentDataProfiles.GetTECEntry( det, ring, beam, disk ).SetValue( digi.strip(), digi.adc() );
      else throw cms::Exception( "[Laser Alignment::fillDataProfiles]" ) << " ** ERROR: digi strip: " << digi.strip() << " out of range for det: " << detRawId << "." << std::endl;
    }

      }

      
    } while( moduleLoop.TECLoop( det, ring, beam, disk ) );
    



    
    // loop TIBTOB modules
    det = 2; beam = 0; pos = 0;
    do {

      // first clear the profile
      currentDataProfiles.GetTIBTOBEntry( det, beam, pos ).SetAllValuesTo( 0. );

      // retrieve the raw id of that module
      const int detRawId = detectorId.GetTIBTOBEntry( det, beam, pos );
      
      if( isRawDigi ) { // we have raw SiStripRawDigis
    
    // search the digis for the raw id
    edm::DetSetVector<SiStripRawDigi>::const_iterator detSetIter = theStripRawDigis->find( detRawId );
    if( detSetIter == theStripRawDigis->end() ) {
      throw cms::Exception( "[Laser Alignment::fillDataProfiles]" ) << " ** ERROR: No raw DetSet found for det: " << detRawId << "." << std::endl;
    }
      
    // fill the digis to the profiles
    edm::DetSet<SiStripRawDigi>::const_iterator digiRangeIterator = detSetIter->data.begin(); // for the loop
    edm::DetSet<SiStripRawDigi>::const_iterator digiRangeStart = digiRangeIterator; // save starting positions
    
    // loop all digis
    for (; digiRangeIterator != detSetIter->data.end(); ++digiRangeIterator ) {
      const SiStripRawDigi& digi = *digiRangeIterator;
      const int channel = distance( digiRangeStart, digiRangeIterator );
      if ( channel >= 0 && channel < 512 ) currentDataProfiles.GetTIBTOBEntry( det, beam, pos ).SetValue( channel, digi.adc() );
      else throw cms::Exception( "[Laser Alignment::fillDataProfiles]" ) << " ** ERROR: raw digi channel: " << channel << " out of range for det: " << detRawId << "." << std::endl;
    }

      }

      else { // we have zero suppressed SiStripDigis

    // search the digis for the raw id
    edm::DetSetVector<SiStripDigi>::const_iterator detSetIter = theStripDigis->find( detRawId );

    // processed DetSets may be missing, just skip
    if( detSetIter == theStripDigis->end() ) continue;

    // fill the digis to the profiles
    edm::DetSet<SiStripDigi>::const_iterator digiRangeIterator = detSetIter->data.begin(); // for the loop
    
    for(; digiRangeIterator != detSetIter->data.end(); ++digiRangeIterator ) {
      const SiStripDigi& digi = *digiRangeIterator;
      if ( digi.strip() < 512 ) currentDataProfiles.GetTIBTOBEntry( det, beam, pos ).SetValue( digi.strip(), digi.adc() );
      else throw cms::Exception( "[Laser Alignment::fillDataProfiles]" ) << " ** ERROR: digi strip: " << digi.strip() << " out of range for det: " << detRawId << "." << std::endl;
    }

      }

    } while( moduleLoop.TIBTOBLoop( det, beam, pos ) );



    // loop TEC AT modules
    det = 0; beam = 0; disk = 0;
    do {

      // first clear the profile
      currentDataProfiles.GetTEC2TECEntry( det, beam, disk ).SetAllValuesTo( 0. );

      // retrieve the raw id of that module
      const int detRawId = detectorId.GetTEC2TECEntry( det, beam, disk );
    
      if( isRawDigi ) { // we have raw SiStripRawDigis
      
    // search the digis for the raw id
    edm::DetSetVector<SiStripRawDigi>::const_iterator detSetIter = theStripRawDigis->find( detRawId );
    if( detSetIter == theStripRawDigis->end() ) {
      throw cms::Exception( "[Laser Alignment::fillDataProfiles]" ) << " ** ERROR: No raw DetSet found for det: " << detRawId << "." << std::endl;
    }
      
    // fill the digis to the profiles
    edm::DetSet<SiStripRawDigi>::const_iterator digiRangeIterator = detSetIter->data.begin(); // for the loop
    edm::DetSet<SiStripRawDigi>::const_iterator digiRangeStart = digiRangeIterator; // save starting positions
      
    // loop all digis
    for (; digiRangeIterator != detSetIter->data.end(); ++digiRangeIterator ) {
      const SiStripRawDigi& digi = *digiRangeIterator;
      const int channel = distance( digiRangeStart, digiRangeIterator );
      if ( channel >= 0 && channel < 512 ) currentDataProfiles.GetTEC2TECEntry( det, beam, disk ).SetValue( channel, digi.adc() );
      else throw cms::Exception( "[Laser Alignment::fillDataProfiles]" ) << " ** ERROR: raw digi channel: " << channel << " out of range for det: " << detRawId << "." << std::endl;
    }
      
      }
    
      else { // we have zero suppressed SiStripDigis
      
    // search the digis for the raw id
    edm::DetSetVector<SiStripDigi>::const_iterator detSetIter = theStripDigis->find( detRawId );
    
    // processed DetSets may be missing, just skip
    if( detSetIter == theStripDigis->end() ) continue;
      
    // fill the digis to the profiles
    edm::DetSet<SiStripDigi>::const_iterator digiRangeIterator = detSetIter->data.begin(); // for the loop
      
    for(; digiRangeIterator != detSetIter->data.end(); ++digiRangeIterator ) {
      const SiStripDigi& digi = *digiRangeIterator;
      if ( digi.strip() < 512 ) currentDataProfiles.GetTEC2TECEntry( det, beam, disk ).SetValue( digi.strip(), digi.adc() );
      else throw cms::Exception( "[Laser Alignment::fillDataProfiles]" ) << " ** ERROR: digi strip: " << digi.strip() << " out of range for det: " << detRawId << "." << std::endl;
    }
      
      }
    
    } while( moduleLoop.TEC2TECLoop( det, beam, disk ) );

  } // theDigiProducersList loop

}





void LaserAlignment::fillPedestalProfiles( edm::ESHandle<SiStripPedestals>& pedestalsHandle ) {

  int det, ring, beam, disk, pos;

  // loop TEC modules (yet without AT)
  det = 0; ring = 0; beam = 0; disk = 0;
  do { // loop using LASGlobalLoop functionality
    SiStripPedestals::Range pedRange = pedestalsHandle->getRange( detectorId.GetTECEntry( det, ring, beam, disk ) );
    for( int strip = 0; strip < 512; ++strip ) {
      int thePedestal = int( pedestalsHandle->getPed( strip, pedRange ) );
      if( thePedestal > 895 ) thePedestal -= 1024;
      pedestalProfiles.GetTECEntry( det, ring, beam, disk ).SetValue( strip, thePedestal );
    }
  } while ( moduleLoop.TECLoop( det, ring, beam, disk ) );


  // TIB & TOB section
  det = 2; beam = 0; pos = 0;
  do { // loop using LASGlobalLoop functionality
    SiStripPedestals::Range pedRange = pedestalsHandle->getRange( detectorId.GetTIBTOBEntry( det, beam, pos ) );
    for( int strip = 0; strip < 512; ++strip ) {
      int thePedestal = int( pedestalsHandle->getPed( strip, pedRange ) );
      if( thePedestal > 895 ) thePedestal -= 1024;
      pedestalProfiles.GetTIBTOBEntry( det, beam, pos ).SetValue( strip, thePedestal );
    }
  } while( moduleLoop.TIBTOBLoop( det, beam, pos ) );


  // TEC2TEC AT section
  det = 0; beam = 0; disk = 0;
  do { // loop using LASGlobalLoop functionality
    SiStripPedestals::Range pedRange = pedestalsHandle->getRange( detectorId.GetTEC2TECEntry( det, beam, disk ) );
    for( int strip = 0; strip < 512; ++strip ) {
      int thePedestal = int( pedestalsHandle->getPed( strip, pedRange ) );
      if( thePedestal > 895 ) thePedestal -= 1024;
      pedestalProfiles.GetTEC2TECEntry( det, beam, disk ).SetValue( strip, thePedestal );
    }
  } while( moduleLoop.TEC2TECLoop( det, beam, disk ) );

}





bool LaserAlignment::isTECBeam( void ) {
  
  int numberOfProfiles = 0;

  int ring = 1; // search all ring6 modules for signals
  for( int det = 0; det < 2; ++det ) {
    for( int beam = 0; beam < 8; ++ beam ) {
      for( int disk = 0; disk < 9; ++disk ) {
    if( isAcceptedProfile.GetTECEntry( det, ring, beam, disk ) == 1 ) numberOfProfiles++;
      }
    }
  }

  LogDebug( "[LaserAlignment::isTECBeam]" ) << " Found: " << numberOfProfiles << "hits." << std::endl;
  std::cout << " [LaserAlignment::isTECBeam] -- Found: " << numberOfProfiles << " hits." << std::endl; 

  if( numberOfProfiles > 10 ) return( true );
  return( false );
 
}






bool LaserAlignment::isATBeam( void ) {

  int numberOfProfiles = 0;

  int det = 2; int beam = 0; int pos = 0; // search all TIB/TOB for signals
  do {
    if( isAcceptedProfile.GetTIBTOBEntry( det, beam, pos ) == 1 ) numberOfProfiles++;
  } while( moduleLoop.TIBTOBLoop( det, beam, pos ) );

  LogDebug( "[LaserAlignment::isATBeam]" ) << " Found: " << numberOfProfiles << "hits." << std::endl;
  std::cout << " [LaserAlignment::isATBeam] -- Found: " << numberOfProfiles << " hits." << std::endl; 

  if( numberOfProfiles > 10 ) return( true );
  return( false );
    
}





double LaserAlignment::getTIBTOBNominalBeamOffset( unsigned int det, unsigned int beam, unsigned int pos ) {

  if( det < 2 || det > 3 || beam > 7 || pos > 5 ) {
    throw cms::Exception( "[LaserAlignment::getTIBTOBNominalBeamOffset]" ) << " ERROR ** Called with nonexisting parameter set: det " << det << " beam " << beam << " pos " << pos << "." << std::endl;
  }

  const double nominalOffsetsTIB[8] = { 0.00035, 2.10687, -2.10827, -0.00173446, 2.10072, -0.00135114, 2.10105, -2.10401 };

  // in tob, modules have alternating orientations along the rods.
  // this is described by the following pattern.
  // (even more confusing, this pattern is inversed for beams 0, 5, 6, 7)
  const int orientationPattern[6] = { -1, 1, 1, -1, -1, 1 };
  const double nominalOffsetsTOB[8] = { 0.00217408, 1.58678, 117.733, 119.321, 120.906, 119.328, 117.743, 1.58947 };


  if( det == 2 ) return( -1. * nominalOffsetsTIB[beam] );

  else {
    if( beam == 0 or beam > 4 ) return( nominalOffsetsTOB[beam] * orientationPattern[pos] );
    else return( -1. * nominalOffsetsTOB[beam] * orientationPattern[pos] );
  }

}




double LaserAlignment::getTEC2TECNominalBeamOffset( unsigned int det, unsigned int beam, unsigned int disk ) {

  if( det > 1 || beam > 7 || disk > 5 ) {
    throw cms::Exception( "[LaserAlignment::getTEC2TECNominalBeamOffset]" ) << " ERROR ** Called with nonexisting parameter set: det " << det << " beam " << beam << " disk " << disk << "." << std::endl;
  }

  const double nominalOffsets[8] = { 0., 2.220, -2.221, 0., 2.214, 0., 2.214, -2.217 };
  
  if( det == 0 ) return -1. * nominalOffsets[beam];
  else return nominalOffsets[beam];

}





void LaserAlignment::CalculateNominalCoordinates( void ) {

  //
  // hard coded data yet...
  //

  // nominal phi values of tec beam / alignment tube hits (parameter is beam 0-7)
  const double tecPhiPositions[8]   = { 0.392699, 1.178097, 1.963495, 2.748894, 3.534292, 4.319690, 5.105088, 5.890486 }; // new values calculated by maple
  const double atPhiPositions[8]    = { 0.392699, 1.289799, 1.851794, 2.748894, 3.645995, 4.319690, 5.216791, 5.778784 }; // new values calculated by maple

  // nominal r values (mm) of hits
  const double tobRPosition = 600.;
  const double tibRPosition = 514.;
  const double tecRPosition[2] = { 564., 840. }; // ring 4,6

  // nominal z values (mm) of hits in barrel (parameter is pos 0-6)
  const double tobZPosition[6] = { 1040., 580., 220., -140., -500., -860. };
  const double tibZPosition[6] = { 620., 380., 180., -100., -340., -540. };

  // nominal z values (mm) of hits in tec (parameter is disk 0-8); FOR TEC-: (* -1.)
  const double tecZPosition[9] = { 1322.5, 1462.5, 1602.5, 1742.5, 1882.5, 2057.5, 2247.5, 2452.5, 2667.5 };
  

  //
  // now we fill these into the nominalCoordinates container;
  // errors are zero for nominal values..
  //

  // loop object and its variables
  LASGlobalLoop moduleLoop;
  int det, ring, beam, disk, pos;

  
  // TEC+- section
  det = 0; ring = 0, beam = 0; disk = 0;
  do {
    
    if( det == 0 ) { // this is TEC+
      nominalCoordinates.SetTECEntry( det, ring, beam, disk, LASCoordinateSet( tecPhiPositions[beam], 0., tecRPosition[ring], 0., tecZPosition[disk], 0. ) );
    }
    else { // now TEC-
      nominalCoordinates.SetTECEntry( det, ring, beam, disk, LASCoordinateSet( tecPhiPositions[beam], 0., tecRPosition[ring], 0., -1. * tecZPosition[disk], 0. ) ); // just * -1.
    }
    
  } while( moduleLoop.TECLoop( det, ring, beam, disk ) );



  // TIB & TOB section
  det = 2; beam = 0; pos = 0;
  do {
    if( det == 2 ) { // this is TIB
      nominalCoordinates.SetTIBTOBEntry( det, beam, pos, LASCoordinateSet( atPhiPositions[beam], 0., tibRPosition, 0., tibZPosition[pos], 0. ) );
    }
    else { // now TOB
      nominalCoordinates.SetTIBTOBEntry( det, beam, pos, LASCoordinateSet( atPhiPositions[beam], 0., tobRPosition, 0., tobZPosition[pos], 0. ) );
    }

  } while( moduleLoop.TIBTOBLoop( det, beam, pos ) );




  // TEC2TEC AT section
  det = 0; beam = 0; disk = 0;
  do {
    
    if( det == 0 ) { // this is TEC+, ring4 only
      nominalCoordinates.SetTEC2TECEntry( det, beam, disk, LASCoordinateSet( atPhiPositions[beam], 0., tecRPosition[0], 0., tecZPosition[disk], 0. ) );
    }
    else { // now TEC-
      nominalCoordinates.SetTEC2TECEntry( det, beam, disk, LASCoordinateSet( atPhiPositions[beam], 0., tecRPosition[0], 0., -1. * tecZPosition[disk], 0. ) ); // just * -1.
    }
    
  } while( moduleLoop.TEC2TECLoop( det, beam, disk ) );


}





double LaserAlignment::ConvertAngle( double angle ) {

  if( angle < -1. * M_PI  || angle > M_PI ) {
    throw cms::Exception(" [LaserAlignment::ConvertAngle] ") << "** ERROR: Called with illegal input angle: " << angle << "." << std::endl;
  }

  if( angle >= 0. ) return angle;
  else return( angle + 2. * M_PI );

}





void LaserAlignment::DumpPosFileSet( LASGlobalData<LASCoordinateSet>& coordinates ) {

  LASGlobalLoop loop;
  int det, ring, beam, disk, pos;

  std:: cout << std::endl << " [LaserAlignment::DumpPosFileSet] -- Dump: " << std::endl;

  // TEC INTERNAL
  det = 0; ring = 0; beam = 0; disk = 0;
  do {
    std::cout << "POS " << det << "\t" << beam << "\t" << disk << "\t" << ring << "\t" << coordinates.GetTECEntry( det, ring, beam, disk ).GetPhi() << "\t" << coordinates.GetTECEntry( det, ring, beam, disk ).GetPhiError() << std::endl;
  } while ( loop.TECLoop( det, ring, beam, disk ) );

  // TIBTOB
  det = 2; beam = 0; pos = 0;
  do {
    std::cout << "POS " << det << "\t" << beam << "\t" << pos << "\t" << "-1" << "\t" << coordinates.GetTIBTOBEntry( det, beam, pos ).GetPhi() << "\t" << coordinates.GetTIBTOBEntry( det, beam, pos ).GetPhiError() << std::endl;
  } while( loop.TIBTOBLoop( det, beam, pos ) );

  // TEC2TEC
  det = 0; beam = 0; disk = 0;
  do {
    std::cout << "POS " << det << "\t" << beam << "\t" << disk << "\t" << "-1" << "\t" << coordinates.GetTEC2TECEntry( det, beam, disk ).GetPhi() << "\t" << coordinates.GetTEC2TECEntry( det, beam, disk ).GetPhiError() << std::endl;
  } while( loop.TEC2TECLoop( det, beam, disk ) );

  std:: cout << std::endl << " [LaserAlignment::DumpPosFileSet] -- End dump: " << std::endl;

}





void LaserAlignment::DumpStripFileSet( LASGlobalData<std::pair<float,float> >& measuredStripPositions ) {

  LASGlobalLoop loop;
  int det, ring, beam, disk, pos;

  std:: cout << std::endl << " [LaserAlignment::DumpStripFileSet] -- Dump: " << std::endl;

  // TEC INTERNAL
  det = 0; ring = 0; beam = 0; disk = 0;
  do {
    std::cout << "STRIP " << det << "\t" << beam << "\t" << disk << "\t" << ring << "\t" << measuredStripPositions.GetTECEntry( det, ring, beam, disk ).first
          << "\t" << measuredStripPositions.GetTECEntry( det, ring, beam, disk ).second << std::endl;
  } while ( loop.TECLoop( det, ring, beam, disk ) );

  // TIBTOB
  det = 2; beam = 0; pos = 0;
  do {
    std::cout << "STRIP " << det << "\t" << beam << "\t" << pos << "\t" << "-1" << "\t" << measuredStripPositions.GetTIBTOBEntry( det, beam, pos ).first
          << "\t" << measuredStripPositions.GetTIBTOBEntry( det, beam, pos ).second << std::endl;
  } while( loop.TIBTOBLoop( det, beam, pos ) );

  // TEC2TEC
  det = 0; beam = 0; disk = 0;
  do {
    std::cout << "STRIP " << det << "\t" << beam << "\t" << disk << "\t" << "-1" << "\t" << measuredStripPositions.GetTEC2TECEntry( det, beam, disk ).first
          << "\t" << measuredStripPositions.GetTEC2TECEntry( det, beam, disk ).second << std::endl;
  } while( loop.TEC2TECLoop( det, beam, disk ) );

  std:: cout << std::endl << " [LaserAlignment::DumpStripFileSet] -- End dump: " << std::endl;
  
  
}





void LaserAlignment::DumpHitmaps( LASGlobalData<int> &numberOfAcceptedProfiles ) {

  std::cout << " [LaserAlignment::DumpHitmaps] -- Dumping hitmap for TEC+:" << std::endl;
  std::cout << " [LaserAlignment::DumpHitmaps] -- Ring4:" << std::endl;
  std::cout << "     disk0   disk1   disk2   disk3   disk4   disk5   disk6   disk7   disk8" << std::endl;

  for( int beam = 0; beam < 8; ++beam ) {
    std::cout << " beam" << beam << ":";
    for( int disk = 0; disk < 9; ++disk ) {
      std::cout << "\t" << numberOfAcceptedProfiles.GetTECEntry( 0, 0, beam, disk );
    }
    std::cout << std::endl;
  }

  std::cout << " [LaserAlignment::DumpHitmaps] -- Ring6:" << std::endl;
  std::cout << "     disk0   disk1   disk2   disk3   disk4   disk5   disk6   disk7   disk8" << std::endl;

  for( int beam = 0; beam < 8; ++beam ) {
    std::cout << " beam" << beam << ":";
    for( int disk = 0; disk < 9; ++disk ) {
      std::cout << "\t" << numberOfAcceptedProfiles.GetTECEntry( 0, 1, beam, disk );
    }
    std::cout << std::endl;
  }

  std::cout << " [LaserAlignment::DumpHitmaps] -- Dumping hitmap for TEC-:" << std::endl;
  std::cout << " [LaserAlignment::DumpHitmaps] -- Ring4:" << std::endl;
  std::cout << "     disk0   disk1   disk2   disk3   disk4   disk5   disk6   disk7   disk8" << std::endl;

  for( int beam = 0; beam < 8; ++beam ) {
    std::cout << " beam" << beam << ":";
    for( int disk = 0; disk < 9; ++disk ) {
      std::cout << "\t" << numberOfAcceptedProfiles.GetTECEntry( 1, 0, beam, disk );
    }
    std::cout << std::endl;
  }

  std::cout << " [LaserAlignment::DumpHitmaps] -- Ring6:" << std::endl;
  std::cout << "     disk0   disk1   disk2   disk3   disk4   disk5   disk6   disk7   disk8" << std::endl;

  for( int beam = 0; beam < 8; ++beam ) {
    std::cout << " beam" << beam << ":";
    for( int disk = 0; disk < 9; ++disk ) {
      std::cout << "\t" << numberOfAcceptedProfiles.GetTECEntry( 1, 1, beam, disk );
    }
    std::cout << std::endl;
  }

  std::cout << " [LaserAlignment::DumpHitmaps] -- End of dump." << std::endl << std::endl;

}





void LaserAlignment::ApplyEndcapMaskingCorrections( LASGlobalData<LASCoordinateSet>& measuredCoordinates, LASGlobalData<LASCoordinateSet>& nominalCoordinates, LASEndcapAlignmentParameterSet& endcapParameters ) {

  // loop the list of modules to be masked
  for( std::vector<unsigned int>::iterator moduleIt = theMaskTecModules.begin(); moduleIt != theMaskTecModules.end(); ++moduleIt ) {

    // loop variables
    LASGlobalLoop moduleLoop;
    int det, ring, beam, disk;

    // this will calculate the corrections from the alignment parameters
    LASEndcapAlgorithm endcapAlgorithm;

    // find the location of the respective module in the container with this loop
    det = 0; ring = 0; beam = 0; disk = 0;
    do {
      
      // here we got it
      if( detectorId.GetTECEntry( det, ring, beam, disk ) == *moduleIt ) {
    
    // the nominal phi value for this module
    const double nominalPhi = nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetPhi();
    
    // the offset from the alignment parameters
    const double phiCorrection = endcapAlgorithm.GetAlignmentParameterCorrection( det, ring, beam, disk, nominalCoordinates, endcapParameters );
    
    // apply the corrections
    measuredCoordinates.GetTECEntry( det, ring, beam, disk ).SetPhi( nominalPhi - phiCorrection );
    
      }
      
    } while ( moduleLoop.TECLoop( det, ring, beam, disk ) );
    
  }

}





void LaserAlignment::ApplyATMaskingCorrections( LASGlobalData<LASCoordinateSet>& measuredCoordinates, LASGlobalData<LASCoordinateSet>& nominalCoordinates, LASBarrelAlignmentParameterSet& atParameters ) {

  // loop the list of modules to be masked
  for( std::vector<unsigned int>::iterator moduleIt = theMaskAtModules.begin(); moduleIt != theMaskAtModules.end(); ++moduleIt ) {

    // loop variables
    LASGlobalLoop moduleLoop;
    int det, beam, disk, pos;

    // this will calculate the corrections from the alignment parameters
    LASAlignmentTubeAlgorithm atAlgorithm;


    // find the location of the respective module in the container with these loops:

    // first TIB+TOB
    det = 2; beam = 0; pos = 0;
    do {

      // here we got it
      if( detectorId.GetTIBTOBEntry( det, beam, pos ) == *moduleIt ) {

    // the nominal phi value for this module
    const double nominalPhi = nominalCoordinates.GetTIBTOBEntry( det, beam, pos ).GetPhi();

    // the offset from the alignment parameters
    const double phiCorrection = atAlgorithm.GetTIBTOBAlignmentParameterCorrection( det, beam, pos, nominalCoordinates, atParameters );

    // apply the corrections
    measuredCoordinates.GetTIBTOBEntry( det, beam, pos ).SetPhi( nominalPhi - phiCorrection );

      }

    } while ( moduleLoop.TIBTOBLoop( det, beam, pos ) );
      
    
    
    // then TEC(AT)  
    det = 0; beam = 0; disk = 0;
    do {
      
      // here we got it
      if( detectorId.GetTEC2TECEntry( det, beam, disk ) == *moduleIt ) {

    // the nominal phi value for this module
    const double nominalPhi = nominalCoordinates.GetTEC2TECEntry( det, beam, disk ).GetPhi();

    // the offset from the alignment parameters
    const double phiCorrection = atAlgorithm.GetTEC2TECAlignmentParameterCorrection( det, beam, disk, nominalCoordinates, atParameters );

    // apply the corrections
    measuredCoordinates.GetTEC2TECEntry( det, beam, disk ).SetPhi( nominalPhi - phiCorrection );

      }
      
    } while ( moduleLoop.TEC2TECLoop( det, beam, disk ) );
    
  }

}





void LaserAlignment::testRoutine( void ) {


  // tracker geom. object for calculating the global beam positions
  const TrackerGeometry& theTracker( *theTrackerGeometry );

  const double atPhiPositions[8] = { 0.392699, 1.289799, 1.851794, 2.748894, 3.645995, 4.319690, 5.216791, 5.778784 };
  const double tecPhiPositions[8] = { 0.392699, 1.178097, 1.963495, 2.748894, 3.534292, 4.319690, 5.105088, 5.890486 };
  const double zPositions[9] = { 125.0, 139.0, 153.0, 167.0, 181.0, 198.5, 217.5, 238.0, 259.5 };
  const double zPositionsTIB[6] = { 62.0, 38.0, 18.0, -10.0, -34.0, -54.0 };
  const double zPositionsTOB[6] = { 104.0, 58.0, 22.0, -14.0, -50.0, -86.0 };

  int det, beam, disk, pos, ring;
  
  // loop TEC+- internal
  det = 0; ring = 0; beam = 0; disk = 0;
  do {

    const double radius = ring?84.0:56.4;

    // access the tracker geometry for this module
    const DetId theDetId( detectorId.GetTECEntry( det, ring, beam, disk ) );
    const StripGeomDetUnit* const theStripDet = dynamic_cast<const StripGeomDetUnit*>( theTracker.idToDet( theDetId ) );
    
    if (theStripDet) {
      const GlobalPoint gp( GlobalPoint::Cylindrical( radius, tecPhiPositions[beam], zPositions[disk] ) );
      
      const LocalPoint lp( theStripDet->surface().toLocal( gp ) );
      std::cout << "__TEC: " << 256. - theStripDet->specificTopology().strip( lp ) << std::endl; 
    }

  } while( moduleLoop.TECLoop( det, ring, beam, disk ) );


  // loop TIBTOB
  det = 2; beam = 0; pos = 0;
  do {

    const double radius = (det==2?51.4:58.4); 
    const double theZ = (det==2?zPositionsTIB[pos]:zPositionsTOB[pos]);

    // access the tracker geometry for this module
    const DetId theDetId( detectorId.GetTIBTOBEntry( det, beam, pos ) );
    const StripGeomDetUnit* const theStripDet = dynamic_cast<const StripGeomDetUnit*>( theTracker.idToDet( theDetId ) );
    
    if (theStripDet) {
      const GlobalPoint gp( GlobalPoint::Cylindrical( radius, atPhiPositions[beam], theZ ) );
      
      const LocalPoint lp( theStripDet->surface().toLocal( gp ) );
      std::cout << "__TIBTOB det " << det << " beam " << beam << " pos " << pos << "  " << 256. - theStripDet->specificTopology().strip( lp );
      std::cout << "           " << theStripDet->position().perp()<< std::endl; 
    }

  } while( moduleLoop.TIBTOBLoop( det, beam, pos ) );

  
  // loop TEC2TEC
  det = 0; beam = 0; disk = 0;
  do {

    const double radius = 56.4;

    // access the tracker geometry for this module
    const DetId theDetId( detectorId.GetTEC2TECEntry( det, beam, disk ) );
    const StripGeomDetUnit* const theStripDet = dynamic_cast<const StripGeomDetUnit*>( theTracker.idToDet( theDetId ) );
    
    if (theStripDet) {
      const GlobalPoint gp( GlobalPoint::Cylindrical( radius, atPhiPositions[beam], zPositions[disk] ) );
      
      const LocalPoint lp( theStripDet->surface().toLocal( gp ) );
      std::cout << "__TEC2TEC det " << det << " beam " << beam << " disk " << disk << "  " << 256. - theStripDet->specificTopology().strip( lp ) << std::endl; 
    }

  } while( moduleLoop.TEC2TECLoop( det, beam, disk ) );


}









// define the SEAL module
#include "FWCore/Framework/interface/MakerMacros.h"

DEFINE_FWK_MODULE(LaserAlignment);




// the ATTIC