MuonGeometryArrange.cc

Go to the documentation of this file.

#include "CondFormats/Alignment/interface/Alignments.h"
#include "CondFormats/Alignment/interface/AlignmentErrorsExtended.h"
#include "CLHEP/Vector/RotationInterfaces.h" 
#include "CondFormats/Alignment/interface/AlignmentSorter.h"

#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "Alignment/CommonAlignment/interface/AlignableObjectId.h"

#include "Alignment/MuonAlignment/interface/AlignableMuon.h"
#include "Geometry/TrackerNumberingBuilder/interface/GeometricDet.h" 
#include "Geometry/Records/interface/IdealGeometryRecord.h"
// The following looks generic enough to use
#include "Alignment/CommonAlignment/interface/Utilities.h"
#include "Alignment/CommonAlignment/interface/SurveyDet.h"
#include "Alignment/CommonAlignment/interface/Alignable.h"
#include "DataFormats/DetId/interface/DetId.h"
#include "DataFormats/MuonDetId/interface/CSCDetId.h"
#include "CondFormats/AlignmentRecord/interface/GlobalPositionRcd.h"
#include "CondFormats/Alignment/interface/DetectorGlobalPosition.h"

#include "Alignment/MuonAlignment/interface/MuonAlignmentInputXML.h"
#include "Alignment/MuonAlignment/interface/MuonAlignment.h"

#include "MuonGeometryArrange.h"
#include "TFile.h" 
#include "TLatex.h"
#include "TArrow.h"
#include "TGraph.h"
#include "TH1F.h" 
#include "TH2F.h" 
#include "CLHEP/Vector/ThreeVector.h"

// Database
#include "FWCore/ServiceRegistry/interface/Service.h"

#include <iostream>
#include <fstream>

MuonGeometryArrange::MuonGeometryArrange(const edm::ParameterSet& cfg) :
  theSurveyIndex(0),
  _levelStrings(cfg.getUntrackedParameter<std::vector<std::string> >("levels")),
  _writeToDB(false), _commonMuonLevel(align::invalid), firstEvent_(true)
{
    referenceMuon=nullptr;
    currentMuon=nullptr;
    // Input is XML
    _inputXMLCurrent = cfg.getUntrackedParameter<std::string> ("inputXMLCurrent");
    _inputXMLReference = cfg.getUntrackedParameter<std::string> ("inputXMLReference");

    //input is ROOT
    _inputFilename1 = cfg.getUntrackedParameter< std::string >
                    ("inputROOTFile1");
    _inputFilename2 = cfg.getUntrackedParameter< std::string >
                    ("inputROOTFile2");
    _inputTreename = cfg.getUntrackedParameter< std::string > ("treeName");
    
    //output file
    _filename = cfg.getUntrackedParameter< std::string > ("outputFile");
    
    
    _weightBy = cfg.getUntrackedParameter< std::string > ("weightBy");
    _detIdFlag = cfg.getUntrackedParameter< bool > ("detIdFlag");
    _detIdFlagFile = cfg.getUntrackedParameter< std::string >
                      ("detIdFlagFile");
    _weightById  = cfg.getUntrackedParameter< bool > ("weightById");
    _weightByIdFile = cfg.getUntrackedParameter< std::string >
             ("weightByIdFile");
    _endcap  = cfg.getUntrackedParameter<int> ("endcapNumber");
    _station = cfg.getUntrackedParameter<int> ("stationNumber");
    _ring    = cfg.getUntrackedParameter<int> ("ringNumber");
    
        
    // if want to use, make id cut list
    if (_detIdFlag){
        std::ifstream fin;
        fin.open( _detIdFlagFile.c_str() );
        
        while (!fin.eof() && fin.good() ){
            
            uint32_t id;
            fin >> id;
            _detIdFlagVector.push_back(id);
        }
        fin.close();
    }       
    
    // turn weightByIdFile into weightByIdVector
    unsigned int lastID=999999999;
    if (_weightById){
        std::ifstream inFile;
        inFile.open( _weightByIdFile.c_str() );
        int ctr = 0;
        while ( !inFile.eof() ){
            ctr++;
            unsigned int listId;
            inFile >> listId;
            inFile.ignore(256, '\n');
            if(listId!=lastID){
              _weightByIdVector.push_back( listId );
            }
            lastID=listId;
        }
        inFile.close();
    }
    
    
    
    //root configuration
    _theFile = new TFile(_filename.c_str(),"RECREATE");
    _alignTree = new TTree("alignTree","alignTree");
    _alignTree->Branch("id", &_id, "id/I");
    _alignTree->Branch("level", &_level, "level/I");
    _alignTree->Branch("mid", &_mid, "mid/I");
    _alignTree->Branch("mlevel", &_mlevel, "mlevel/I");
    _alignTree->Branch("sublevel", &_sublevel, "sublevel/I");
    _alignTree->Branch("x", &_xVal, "x/F");
    _alignTree->Branch("y", &_yVal, "y/F");
    _alignTree->Branch("z", &_zVal, "z/F");
    _alignTree->Branch("r", &_rVal, "r/F");
    _alignTree->Branch("phi", &_phiVal, "phi/F");
    _alignTree->Branch("eta", &_etaVal, "eta/F");
    _alignTree->Branch("alpha", &_alphaVal, "alpha/F");
    _alignTree->Branch("beta", &_betaVal, "beta/F");
    _alignTree->Branch("gamma", &_gammaVal, "gamma/F");
    _alignTree->Branch("dx", &_dxVal, "dx/F");
    _alignTree->Branch("dy", &_dyVal, "dy/F");
    _alignTree->Branch("dz", &_dzVal, "dz/F");
    _alignTree->Branch("dr", &_drVal, "dr/F");
    _alignTree->Branch("dphi", &_dphiVal, "dphi/F");
    _alignTree->Branch("dalpha", &_dalphaVal, "dalpha/F");
    _alignTree->Branch("dbeta", &_dbetaVal, "dbeta/F");
    _alignTree->Branch("dgamma", &_dgammaVal, "dgamma/F");
    _alignTree->Branch("ldx", &_ldxVal, "ldx/F");
    _alignTree->Branch("ldy", &_ldyVal, "ldy/F");
    _alignTree->Branch("ldz", &_ldzVal, "ldz/F");
    _alignTree->Branch("ldr", &_ldrVal, "ldr/F");
    _alignTree->Branch("ldphi", &_ldphiVal, "ldphi/F");
    _alignTree->Branch("useDetId", &_useDetId, "useDetId/I");
    _alignTree->Branch("detDim", &_detDim, "detDim/I");
    _alignTree->Branch("rotx",&_rotxVal, "rotx/F"); 
    _alignTree->Branch("roty",&_rotyVal, "roty/F"); 
    _alignTree->Branch("rotz",&_rotzVal, "rotz/F"); 
    _alignTree->Branch("drotx",&_drotxVal, "drotx/F");  
    _alignTree->Branch("droty",&_drotyVal, "droty/F");  
    _alignTree->Branch("drotz",&_drotzVal, "drotz/F");  
    _alignTree->Branch("surW", &_surWidth, "surW/F");
    _alignTree->Branch("surL", &_surLength, "surL/F");
    _alignTree->Branch("surRot", &_surRot, "surRot[9]/D");

    _mgacollection.clear(); 
}
void MuonGeometryArrange::endHist(){
 // Unpack the list and create ntuples here.

   int size=_mgacollection.size();
   if(size<=0) return;  // nothing to do here.
   std::vector<float> xp(size+1);
   std::vector<float> yp(size+1);
   int i;
   float minV, maxV;
   int minI, maxI;

   minV=99999999.; maxV=-minV;  minI=9999999; maxI=-minI;
   TGraph* grx=nullptr;
   TH2F* dxh=nullptr;

// for position plots:
   for(i=0; i<size; i++){
     if(_mgacollection[i].phipos<minI) minI=_mgacollection[i].phipos;
     if(_mgacollection[i].phipos>maxI) maxI=_mgacollection[i].phipos;
     xp[i]=_mgacollection[i].phipos;
   }
   if(minI>=maxI) return;   // can't do anything?
   xp[size]=xp[size-1]+1;   // wraparound point

   if(1<minI) minI=1;
   if(size>maxI) maxI=size;
   maxI++;  // allow for wraparound to show neighbors
   int sizeI=maxI+1-minI;
   float smi=minI-1;
   float sma=maxI+1;
 

// Dx plot

   for(i=0; i<size; i++){
     if(_mgacollection[i].ldx<minV) minV=_mgacollection[i].ldx;
     if(_mgacollection[i].ldx>maxV) maxV=_mgacollection[i].ldx;
     yp[i]=_mgacollection[i].ldx;
   }
   yp[size]=yp[0];  // wraparound point
   
   makeGraph(sizeI, smi, sma, minV, maxV,
       dxh, grx, "delX_vs_position", "Local #delta X vs position", 
       "GdelX_vs_position","#delta x in cm", xp.data(), yp.data(), size);
// Dy plot
   minV=99999999.; maxV=-minV;
   for(i=0; i<size; i++){
     if(_mgacollection[i].ldy<minV) minV=_mgacollection[i].ldy;
     if(_mgacollection[i].ldy>maxV) maxV=_mgacollection[i].ldy;
     yp[i]=_mgacollection[i].ldy;
   }
   yp[size]=yp[0];  // wraparound point
   
   makeGraph(sizeI, smi, sma, minV, maxV,
       dxh, grx, "delY_vs_position", "Local #delta Y vs position", 
       "GdelY_vs_position","#delta y in cm", xp.data(), yp.data(), size);

// Dz plot
   minV=99999999.; maxV=-minV;
   for(i=0; i<size; i++){
     if(_mgacollection[i].dz<minV) minV=_mgacollection[i].dz;
     if(_mgacollection[i].dz>maxV) maxV=_mgacollection[i].dz;
     yp[i]=_mgacollection[i].dz;
   }
   yp[size]=yp[0];  // wraparound point
   
   makeGraph(sizeI, smi, sma, minV, maxV,
       dxh, grx, "delZ_vs_position", "Local #delta Z vs position", 
       "GdelZ_vs_position","#delta z in cm", xp.data(), yp.data(), size);

// Dphi plot
   minV=99999999.; maxV=-minV;
   for(i=0; i<size; i++){
     if(_mgacollection[i].dphi<minV) minV=_mgacollection[i].dphi;
     if(_mgacollection[i].dphi>maxV) maxV=_mgacollection[i].dphi;
     yp[i]=_mgacollection[i].dphi;
   }
   yp[size]=yp[0];  // wraparound point
   
   makeGraph(sizeI, smi, sma, minV, maxV,
       dxh, grx, "delphi_vs_position", "#delta #phi vs position", 
       "Gdelphi_vs_position","#delta #phi in radians", xp.data(), yp.data(), size);

// Dr plot
   minV=99999999.; maxV=-minV;
   for(i=0; i<size; i++){
     if(_mgacollection[i].dr<minV) minV=_mgacollection[i].dr;
     if(_mgacollection[i].dr>maxV) maxV=_mgacollection[i].dr;
     yp[i]=_mgacollection[i].dr;
   }
   yp[size]=yp[0];  // wraparound point
   
   makeGraph(sizeI, smi, sma, minV, maxV,
       dxh, grx, "delR_vs_position", "#delta R vs position", 
       "GdelR_vs_position","#delta R in cm", xp.data(), yp.data(), size);

// Drphi plot
   minV=99999999.; maxV=-minV;
   for(i=0; i<size; i++){
     float ttemp=_mgacollection[i].r*_mgacollection[i].dphi;
     if(ttemp<minV) minV=ttemp;
     if(ttemp>maxV) maxV=ttemp;
     yp[i]=ttemp;
   }
   yp[size]=yp[0];  // wraparound point
   
   makeGraph(sizeI, smi, sma, minV, maxV,
       dxh, grx, "delRphi_vs_position", "R #delta #phi vs position", 
       "GdelRphi_vs_position","R #delta #phi in cm", xp.data(), yp.data(), size);

// Dalpha plot
   minV=99999999.; maxV=-minV;
   for(i=0; i<size; i++){
     if(_mgacollection[i].dalpha<minV) minV=_mgacollection[i].dalpha;
     if(_mgacollection[i].dalpha>maxV) maxV=_mgacollection[i].dalpha;
     yp[i]=_mgacollection[i].dalpha;
   }
   yp[size]=yp[0];  // wraparound point
   
   makeGraph(sizeI, smi, sma, minV, maxV,
       dxh, grx, "delalpha_vs_position", "#delta #alpha vs position", 
       "Gdelalpha_vs_position","#delta #alpha in rad", xp.data(), yp.data(), size);

// Dbeta plot
   minV=99999999.; maxV=-minV;
   for(i=0; i<size; i++){
     if(_mgacollection[i].dbeta<minV) minV=_mgacollection[i].dbeta;
     if(_mgacollection[i].dbeta>maxV) maxV=_mgacollection[i].dbeta;
     yp[i]=_mgacollection[i].dbeta;
   }
   yp[size]=yp[0];  // wraparound point
   
   makeGraph(sizeI, smi, sma, minV, maxV,
       dxh, grx, "delbeta_vs_position", "#delta #beta vs position", 
       "Gdelbeta_vs_position","#delta #beta in rad", xp.data(), yp.data(), size);

// Dgamma plot
   minV=99999999.; maxV=-minV;
   for(i=0; i<size; i++){
     if(_mgacollection[i].dgamma<minV) minV=_mgacollection[i].dgamma;
     if(_mgacollection[i].dgamma>maxV) maxV=_mgacollection[i].dgamma;
     yp[i]=_mgacollection[i].dgamma;
   }
   yp[size]=yp[0];  // wraparound point
   
   makeGraph(sizeI, smi, sma, minV, maxV,
       dxh, grx, "delgamma_vs_position", "#delta #gamma vs position", 
       "Gdelgamma_vs_position","#delta #gamma in rad", xp.data(), yp.data(), size);

// Drotx plot
   minV=99999999.; maxV=-minV;
   for(i=0; i<size; i++){
     if(_mgacollection[i].drotx<minV) minV=_mgacollection[i].drotx;
     if(_mgacollection[i].drotx>maxV) maxV=_mgacollection[i].drotx;
     yp[i]=_mgacollection[i].drotx;
   }
   yp[size]=yp[0];  // wraparound point
   
   makeGraph(sizeI, smi, sma, minV, maxV,
       dxh, grx, "delrotX_vs_position", "#delta rotX vs position", 
       "GdelrotX_vs_position","#delta rotX in rad", xp.data(), yp.data(), size);

// Droty plot
   minV=99999999.; maxV=-minV;
   for(i=0; i<size; i++){
     if(_mgacollection[i].droty<minV) minV=_mgacollection[i].droty;
     if(_mgacollection[i].droty>maxV) maxV=_mgacollection[i].droty;
     yp[i]=_mgacollection[i].droty;
   }
   yp[size]=yp[0];  // wraparound point
   
   makeGraph(sizeI, smi, sma, minV, maxV,
       dxh, grx, "delrotY_vs_position", "#delta rotY vs position", 
       "GdelrotY_vs_position","#delta rotY in rad", xp.data(), yp.data(), size);

// Drotz plot
   minV=99999999.; maxV=-minV;
   for(i=0; i<size; i++){
     if(_mgacollection[i].drotz<minV) minV=_mgacollection[i].drotz;
     if(_mgacollection[i].drotz>maxV) maxV=_mgacollection[i].drotz;
     yp[i]=_mgacollection[i].drotz;
   }
   yp[size]=yp[0];  // wraparound point
   
   makeGraph(sizeI, smi, sma, minV, maxV,
       dxh, grx, "delrotZ_vs_position", "#delta rotZ vs position", 
       "GdelrotZ_vs_position","#delta rotZ in rad", xp.data(), yp.data(), size);



// Vector plots
// First find the maximum length of sqrt(dx*dx+dy*dy):  we'll have to
// scale these for visibility
   maxV=-99999999.;
   float ttemp, rtemp;
   float maxR=-9999999.;
   for(i=0; i<size; i++){
     ttemp= sqrt(_mgacollection[i].dx*_mgacollection[i].dx+
        _mgacollection[i].dy*_mgacollection[i].dy);
     rtemp= sqrt(_mgacollection[i].x*_mgacollection[i].x+
        _mgacollection[i].y*_mgacollection[i].y);
     if(ttemp>maxV) maxV=ttemp;
     if(rtemp>maxR) maxR=rtemp;
   }
   
  // Don't try to scale rediculously small values 
   float smallestVcm=.001;  // 10 microns
   if(maxV<smallestVcm) maxV=smallestVcm;
   float scale=0.;
   float lside=1.1*maxR;
   if(lside<=0) lside=100.;
   if(maxV>0){scale=.09*lside/maxV;}    // units of pad length!
   char scalename[50];
   int ret=snprintf(scalename,50,"#delta #bar{x}   length =%f cm",maxV);
  // If ret<=0 we don't want to print the scale!
  
   if(ret>0){
     dxh=new TH2F("vecdrplot",scalename,80,-lside,lside,80,-lside,lside);
   }
   else{
     dxh=new TH2F("vecdrplot","delta #bar{x} Bad scale",80,-lside,lside,80,-lside,lside);
   }
   dxh->GetXaxis()->SetTitle("x in cm");
   dxh->GetYaxis()->SetTitle("y in cm");
   dxh->SetStats(kFALSE);
   dxh->Draw();
   TArrow* arrow;
   for(i=0; i<size; i++){
     ttemp= sqrt(_mgacollection[i].dx*_mgacollection[i].dx+
       _mgacollection[i].dy*_mgacollection[i].dy);
//     ttemp=ttemp*scale;
     float nx=_mgacollection[i].x+scale*_mgacollection[i].dx;
     float ny=_mgacollection[i].y+scale*_mgacollection[i].dy;
     arrow = new TArrow(_mgacollection[i].x, 
                   _mgacollection[i].y, nx, ny);// ttemp*.3*.05, "->");
     arrow->SetLineWidth(2); arrow->SetArrowSize(ttemp*.2*.05/maxV);
     arrow->SetLineColor(1); arrow->SetLineStyle(1);
     arrow->Paint();
     dxh->GetListOfFunctions()->Add(static_cast<TObject*>(arrow));
//     arrow->Draw();
//     arrow->Write();
   }
   dxh->Write();
 
   _theFile->Write();   
   _theFile->Close();


}
void MuonGeometryArrange::makeGraph(int sizeI, float smi, float sma, 
    float minV, float maxV,
    TH2F* dxh, TGraph* grx, const char* name, const char* title, 
    const char* titleg, const char* axis,
    const float* xp, const float* yp, int size){

  if(minV>=maxV || smi>=sma || sizeI<=1 || xp==nullptr || yp==nullptr) return;
    // out of bounds, bail
  float diff=maxV-minV;
  float over=.05*diff;
  double ylo=minV-over;
  double yhi=maxV+over;
  double dsmi, dsma;
  dsmi=smi; dsma=sma;
  dxh= new TH2F(name, title,
    sizeI+2, dsmi, dsma, 50, ylo, yhi);
  dxh->GetXaxis()->SetTitle("Position around ring");
  dxh->GetYaxis()->SetTitle(axis);
  dxh->SetStats(kFALSE);
  dxh->Draw();
  grx = new TGraph(size, xp, yp);
  grx->SetName(titleg);
  grx->SetTitle(title);
  grx->SetMarkerColor(2); grx->SetMarkerStyle(3);
  grx->GetXaxis()->SetLimits(dsmi, dsma);
  grx->GetXaxis()->SetTitle("position number");
  grx->GetYaxis()->SetLimits(ylo,yhi);
  grx->GetYaxis()->SetTitle(axis);
  grx->Draw("A*");
  grx->Write();
  return; 
}
void MuonGeometryArrange::beginJob(){
  firstEvent_ = true;
}

void MuonGeometryArrange::createROOTGeometry(const edm::EventSetup& iSetup){}
void MuonGeometryArrange::analyze(const edm::Event&, 
                  const edm::EventSetup& iSetup){
  if (firstEvent_) {
    
    // My stuff
    MuonAlignmentInputXML inputMethod1(_inputXMLCurrent);
    inputAlign1 = new MuonAlignment(iSetup, inputMethod1);
    inputAlign1->fillGapsInSurvey(0, 0);
    MuonAlignmentInputXML inputMethod2(_inputXMLReference);
    inputAlign2 = new MuonAlignment(iSetup, inputMethod2);
    inputAlign2->fillGapsInSurvey(0, 0);
    MuonAlignmentInputXML inputMethod3(_inputXMLReference);
    inputAlign2a = new MuonAlignment(iSetup, inputMethod3);
    inputAlign2a->fillGapsInSurvey(0, 0);
    
    inputGeometry1 = static_cast<Alignable*> (inputAlign1->getAlignableMuon());
    inputGeometry2 = static_cast<Alignable*> (inputAlign2->getAlignableMuon());
    auto inputGeometry2Copy2 = inputAlign2a->getAlignableMuon();

    //setting the levels being used in the geometry comparator
    edm::LogInfo("MuonGeometryArrange") << "levels: " << _levelStrings.size();
    for (const auto& level: _levelStrings){
      theLevels.push_back(inputGeometry2Copy2->objectIdProvider().stringToId(level));
      edm::LogInfo("MuonGeometryArrange") << "level: " << level;
    }
    
    //compare the goemetries
    compare(inputGeometry1, inputGeometry2, inputGeometry2Copy2);
    
    //write out ntuple
    //might be better to do within output module
    _theFile->cd();
    _alignTree->Write();
    endHist();
    //   _theFile->Close();
    
    firstEvent_ = false;
  }
}

void MuonGeometryArrange::compare(Alignable* refAli, Alignable* curAli,
                Alignable* curAliCopy2){

 // First sanity
  if(refAli==nullptr){return;}  
  if(curAli==nullptr){return;}
  
  const std::vector<Alignable*>& refComp = refAli->components();
  const std::vector<Alignable*>& curComp = curAli->components();
  const std::vector<Alignable*>& curComp2 = curAliCopy2->components();
  compareGeometries(refAli, curAli, curAliCopy2);

  int nComp=refComp.size();
  for(int i=0; i<nComp; i++){
    compare(refComp[i], curComp[i], curComp2[i]);
  }
  return;
}
  
void MuonGeometryArrange::compareGeometries(Alignable* refAli, 
        Alignable* curAli, Alignable* curCopy){
 // First sanity
  if(refAli==nullptr){return;}  
  if(curAli==nullptr){return;}
 // Is this the Ring we want to align?  If so it will contain the
 // chambers specified in the configuration file    
  if(!isMother(refAli)) return; // Not the desired alignable object
 // But... There are granddaughters involved--and I don't want to monkey with
 // the layers of the chambers.  So, if the mother of this is also an approved
 // mother, bail.
  if(isMother(refAli->mother() )) return;
  const std::vector<Alignable*>& refComp = refAli->components();
  const std::vector<Alignable*>& curComp = curCopy->components();
  if(refComp.size()!=curComp.size()){
     return;
  }
 // GlobalVectors is a vector of GlobalVector which is a 3D vector
  align::GlobalVectors originalVectors;
  align::GlobalVectors currentVectors;
  align::GlobalVectors originalRelativeVectors;
  align::GlobalVectors currentRelativeVectors;

    
  int nComp = refComp.size();
  int nUsed = 0;
 // Use the total displacements here:
  CLHEP::Hep3Vector TotalX, TotalL;
  TotalX.set(0.,0.,0.);   TotalL.set(0., 0., 0.);
//  CLHEP::Hep3Vector* Rsubtotal, Wsubtotal, DRsubtotal, DWsubtotal;
  std::vector<CLHEP::Hep3Vector> Positions;
  std::vector<CLHEP::Hep3Vector> DelPositions;

  double xrcenter=0.;
  double yrcenter=0.;
  double zrcenter=0.;
  double xccenter=0.;
  double yccenter=0.;
  double zccenter=0.;

  bool useIt;
 // Create the "center" for the reference alignment chambers, and
 // load a vector of their centers
  for(int ich=0; ich<nComp; ich++){
    useIt=true;
    if(_weightById){
      if(!align::readModuleList(curComp[ich]->id(), curComp[ich]->id(), _weightByIdVector))
         useIt=false;
    }
    if(!useIt) continue;
    align::GlobalVectors curVs;
    align::createPoints(&curVs, refComp[ich], 
      _weightBy, _weightById, _weightByIdVector);
    align::GlobalVector pointsCM = align::centerOfMass(curVs);
    originalVectors.push_back(pointsCM);
    nUsed++;
    xrcenter+= pointsCM.x();
    yrcenter+= pointsCM.y();
    zrcenter+= pointsCM.z();
  }
  xrcenter=xrcenter/nUsed;
  yrcenter=yrcenter/nUsed;
  zrcenter=zrcenter/nUsed;

 // Create the "center" for the current alignment chambers, and
 // load a vector of their centers
   for(int ich=0; ich<nComp; ich++){
    useIt=true;
    if(_weightById){
      if(!align::readModuleList(curComp[ich]->id(), curComp[ich]->id(), _weightByIdVector))
         useIt=false;
    }
    if(!useIt)continue;
    align::GlobalVectors curVs;
    align::createPoints(&curVs, curComp[ich], 
      _weightBy, _weightById, _weightByIdVector);
    align::GlobalVector pointsCM = align::centerOfMass(curVs);
    currentVectors.push_back(pointsCM);

    xccenter+= pointsCM.x();
    yccenter+= pointsCM.y();
    zccenter+= pointsCM.z();
  }
  xccenter=xccenter/nUsed;
  yccenter=yccenter/nUsed;
  zccenter=zccenter/nUsed;


 // OK, now load the <very approximate> vectors from the ring "centers"
  align::GlobalVector CCur(xccenter, yccenter, zccenter);
  align::GlobalVector CRef(xrcenter, yrcenter, zrcenter);
  int nCompR = currentVectors.size();
  for(int ich=0; ich<nCompR; ich++){
    originalRelativeVectors.push_back(originalVectors[ich]-CRef);
    currentRelativeVectors.push_back(currentVectors[ich]-CCur);
  }

 // All right.  Now let the hacking begin.
 // First out of the gate let's try using the raw values and see what
 // diffRot does for us.

 
  align::RotationType rtype3=align::diffRot(currentRelativeVectors,
                                            originalRelativeVectors);


  align::EulerAngles angles(3);
  angles = align::toAngles(rtype3);

  for(int ich=0; ich<nComp; ich++){
    if(_weightById){
      if(!align::readModuleList(curComp[ich]->id(), curComp[ich]->id(), _weightByIdVector))
         continue;
    }
      CLHEP::Hep3Vector Rtotal, Wtotal;
      Rtotal.set(0.,0.,0.); Wtotal.set(0.,0.,0.);
      for (int i = 0; i < 100; i++){
        AlgebraicVector diff = align::diffAlignables(refComp[ich],curComp[ich],
                           _weightBy, _weightById, _weightByIdVector);
        CLHEP::Hep3Vector dR(diff[0],diff[1],diff[2]);
        Rtotal+=dR;
        CLHEP::Hep3Vector dW(diff[3],diff[4],diff[5]);
        CLHEP::HepRotation rot(Wtotal.unit(),Wtotal.mag());
        CLHEP::HepRotation drot(dW.unit(),dW.mag());
        rot*=drot;
        Wtotal.set(rot.axis().x()*rot.delta(),
               rot.axis().y()*rot.delta(), rot.axis().z()*rot.delta());
        align::moveAlignable(curComp[ich], diff);
        float tolerance = 1e-7;
        AlgebraicVector check = align::diffAlignables(refComp[ich],curComp[ich],
                        _weightBy, _weightById, _weightByIdVector);
        align::GlobalVector checkR(check[0],check[1],check[2]);
        align::GlobalVector checkW(check[3],check[4],check[5]);
        DetId detid(refComp[ich]->id());
        if ((checkR.mag() > tolerance)||(checkW.mag() > tolerance)){
//   edm::LogInfo("CompareGeoms") << "Tolerance Exceeded!(alObjId: " 
//       << refAli->alignableObjectId()
//   << ", rawId: " << refComp[ich]->geomDetId().rawId()
//   << ", subdetId: "<< detid.subdetId() << "): " << diff;
        }
        else{
         TotalX+=Rtotal;
         break;
        }     // end of else
      }  // end of for on int i
  } // end of for on ich

  // At this point we should have a total displacement and total L
   TotalX=TotalX/nUsed;

  // Now start again!
   AlgebraicVector change(6);
   change(1)=TotalX.x();
   change(2)=TotalX.y();
   change(3)=TotalX.z();

   change(4)=angles[0];
   change(5)=angles[1];
   change(6)=angles[2];
   align::moveAlignable(curAli, change);    // move as a chunk

  // Now get the components again.  They should be in new locations
   const std::vector<Alignable*>& curComp2 = curAli->components();

  for(int ich=0; ich<nComp; ich++){
     CLHEP::Hep3Vector Rtotal, Wtotal;
     Rtotal.set(0.,0.,0.); Wtotal.set(0.,0.,0.);
     if(_weightById){
       if(!align::readModuleList(curComp[ich]->id(), curComp[ich]->id(), _weightByIdVector))
         continue;
     }
    
     for (int i = 0; i < 100; i++){
       AlgebraicVector diff = align::diffAlignables(refComp[ich],curComp2[ich],
                          _weightBy, _weightById, _weightByIdVector);
       CLHEP::Hep3Vector dR(diff[0],diff[1],diff[2]);
       Rtotal+=dR;
       CLHEP::Hep3Vector dW(diff[3],diff[4],diff[5]);
       CLHEP::HepRotation rot(Wtotal.unit(),Wtotal.mag());
       CLHEP::HepRotation drot(dW.unit(),dW.mag());
       rot*=drot;
       Wtotal.set(rot.axis().x()*rot.delta(), rot.axis().y()*rot.delta(),
                  rot.axis().z()*rot.delta());
       align::moveAlignable(curComp2[ich], diff);
       float tolerance = 1e-7;
       AlgebraicVector check = align::diffAlignables(refComp[ich],curComp2[ich],
                       _weightBy, _weightById, _weightByIdVector);
       align::GlobalVector checkR(check[0],check[1],check[2]);
       align::GlobalVector checkW(check[3],check[4],check[5]);
       if ((checkR.mag() > tolerance)||(checkW.mag() > tolerance)){}
       else{break;}
     }   // end of for on int i
     AlgebraicVector TRtot(6);
     TRtot(1) = Rtotal.x(); TRtot(2) = Rtotal.y(); TRtot(3) = Rtotal.z();
     TRtot(4) = Wtotal.x(); TRtot(5) = Wtotal.y(); TRtot(6) = Wtotal.z();
     fillTree(refComp[ich], TRtot);
  }    // end of for on ich



}


void MuonGeometryArrange::fillTree(Alignable *refAli, const AlgebraicVector& diff){
    
    
    _id = refAli->id();
    _level = refAli->alignableObjectId();
    //need if ali has no mother
    if (refAli->mother()){
        _mid = refAli->mother()->geomDetId().rawId();
        _mlevel = refAli->mother()->alignableObjectId();
    }
    else{
        _mid = -1;
        _mlevel = -1;
    }
    DetId detid(_id);
    _sublevel = detid.subdetId();
    int ringPhiPos=-99;
    if(detid.det()==DetId::Muon && detid.subdetId()== MuonSubdetId::CSC){
       CSCDetId cscId(refAli->geomDetId());
       ringPhiPos = cscId.chamber();
    }
    _xVal = refAli->globalPosition().x();
    _yVal = refAli->globalPosition().y();
    _zVal = refAli->globalPosition().z();
    align::GlobalVector vec(_xVal,_yVal,_zVal);
    _rVal = vec.perp();
    _phiVal = vec.phi();
    _etaVal = vec.eta();
    align::RotationType rot = refAli->globalRotation();
    align::EulerAngles eulerAngles = align::toAngles(rot);
    _rotxVal = atan2(rot.yz(), rot.zz());
    float ttt=-rot.xz();
    if(ttt>1.) ttt=1.;
    if(ttt<-1.) ttt=-1.;
    _rotyVal = asin(ttt);
    _rotzVal = atan2(rot.xy(), rot.xx());
    _alphaVal = eulerAngles[0];
    _betaVal = eulerAngles[1];
    _gammaVal = eulerAngles[2];
    _dxVal = diff[0];
    _dyVal = diff[1];
    _dzVal = diff[2];
    //getting dR and dPhi
    align::GlobalVector vRef(_xVal,_yVal,_zVal);
    align::GlobalVector vCur(_xVal - _dxVal, _yVal-_dyVal, _zVal - _dzVal);
    _drVal = vCur.perp() - vRef.perp();
    _dphiVal = vCur.phi() - vRef.phi();
    
    _dalphaVal = diff[3];
    _dbetaVal = diff[4];
    _dgammaVal = diff[5];
    _drotxVal=-999.; _drotyVal=-999.; _drotzVal=-999.; 

    align::EulerAngles deuler(3);
    deuler(1)=_dalphaVal;
    deuler(2)= _dbetaVal;
    deuler(3)= _dgammaVal;
    align::RotationType drot = align::toMatrix(deuler);
        double xx=rot.xx();
        double xy=rot.xy();
        double xz=rot.xz();
        double yx=rot.yx();
        double yy=rot.yy();
        double yz=rot.yz();
        double zx=rot.zx();
        double zy=rot.zy();
        double zz=rot.zz();
        double detrot=(zz*yy - zy*yz)*xx + (-zz*yx + zx*yz)*xy + (zy*yx - zx*yy)*xz;
        detrot=1/detrot;
        double ixx=(zz*yy - zy*yz)*detrot;
        double ixy=(-zz*xy + zy*xz)*detrot;
        double ixz=(yz*xy - yy*xz)*detrot;
        double iyx=(-zz*yx + zx*yz)*detrot;
        double iyy=(zz*xx - zx*xz)*detrot;
        double iyz=(-yz*xx + yx*xz)*detrot;
        double izx=(zy*yx - zx*yy)*detrot;
        double izy=(-zy*xx + zx*xy)*detrot;
        double izz=(yy*xx - yx*xy)*detrot;
        align::RotationType invrot(ixx,ixy,ixz, iyx,iyy,iyz, izx,izy,izz);
    align::RotationType prot = rot*drot*invrot;
//  align::RotationType prot = rot*drot;
    float protx; //, proty, protz;
    protx = atan2(prot.yz(), prot.zz());
    _drotxVal = protx;//_rotxVal-protx; //atan2(drot.yz(), drot.zz());
    ttt=-prot.xz();
    if(ttt>1.) ttt=1.;
    if(ttt<-1.) ttt=-1.;
    _drotyVal = asin(ttt);// -_rotyVal;
    _drotzVal = atan2(prot.xy(), prot.xx());// - _rotzVal;
// Above does not account for 2Pi wraparounds!
// Prior knowledge:  these are supposed to be small rotations.  Therefore:
    if(_drotxVal>3.141592656) _drotxVal=-6.2831853072+_drotxVal;
    if(_drotxVal<-3.141592656) _drotxVal=6.2831853072+_drotxVal;
    if(_drotyVal>3.141592656) _drotyVal=-6.2831853072+_drotyVal;
    if(_drotyVal<-3.141592656) _drotyVal=6.2831853072+_drotyVal;
    if(_drotzVal>3.141592656) _drotzVal=-6.2831853072+_drotzVal;
    if(_drotzVal<-3.141592656) _drotzVal=6.2831853072+_drotzVal;
    
    _ldxVal=-999.; _ldyVal=-999.; _ldxVal=-999.; 
    _ldrVal=-999.; _ldphiVal=-999; // set fake

//  if(refAli->alignableObjectId() == align::AlignableDetUnit){
     align::GlobalVector dV(_dxVal, _dyVal, _dzVal); 
     align::LocalVector pointL = refAli->surface().toLocal(dV);
         //align::LocalVector pointL = (refAli->mother())->surface().toLocal(dV);
     _ldxVal=pointL.x(); _ldyVal=pointL.y(); _ldzVal=pointL.z();
     _ldphiVal=pointL.phi(); _ldrVal=pointL.perp();
//  }
    //detIdFlag
    if (refAli->alignableObjectId() == align::AlignableDetUnit){
      if (_detIdFlag){
       if ((passIdCut(refAli->id()))||(passIdCut(refAli->mother()->id()))){
        _useDetId = 1;
       }
       else{
        _useDetId = 0;
       }
      }
    }
    // det module dimension
    if (refAli->alignableObjectId() == align::AlignableDetUnit){
       if (refAli->mother()->alignableObjectId() != align::AlignableDet){
         _detDim = 1;}
       else if (refAli->mother()->alignableObjectId() == 
                                   align::AlignableDet) {_detDim = 2;}
    }
    else _detDim = 0;
    
    
    
    _surWidth = refAli->surface().width();
    _surLength = refAli->surface().length();
    align::RotationType rt = refAli->globalRotation();
    _surRot[0] = rt.xx(); _surRot[1] = rt.xy(); _surRot[2] = rt.xz();
    _surRot[3] = rt.yx(); _surRot[4] = rt.yy(); _surRot[5] = rt.yz();
    _surRot[6] = rt.zx(); _surRot[7] = rt.zy(); _surRot[8] = rt.zz();

        MGACollection holdit;
    holdit.id=_id;  holdit.level=_level;  holdit.mid=_mid;
    holdit.mlevel=_mlevel;
    holdit.sublevel=_sublevel;
    holdit.x=_xVal;  holdit.y=_yVal;  holdit.z=_zVal;   
    holdit.r=_rVal;  holdit.phi=_phiVal;  holdit.eta=_etaVal;
    holdit.alpha=_alphaVal;  holdit.beta=_betaVal;  holdit.gamma=_gammaVal;
    holdit.dx=_dxVal;  holdit.dy=_dyVal;  holdit.dz=_dzVal; 
    holdit.dr=_drVal;  holdit.dphi=_dphiVal; 
    holdit.dalpha=_dalphaVal;  holdit.dbeta=_dbetaVal;  
    holdit.dgamma=_dgammaVal;
    holdit.useDetId=_useDetId; holdit.detDim=_detDim;
    holdit.surW=_surWidth; holdit.surL=_surLength;
    holdit.ldx=_ldxVal; holdit.ldy=_ldyVal; holdit.ldz=_ldzVal;
    holdit.ldr=_ldrVal; holdit.ldphi=_ldphiVal;
    holdit.rotx=_rotxVal; holdit.roty=_rotyVal; holdit.rotz=_rotzVal;
    holdit.drotx=_drotxVal; holdit.droty=_drotyVal; holdit.drotz=_drotzVal;
    for(int i=0; i<9; i++){holdit.surRot[i]=_surRot[i];}
    holdit.phipos=ringPhiPos;
    _mgacollection.push_back(holdit);

    
    //Fill
    _alignTree->Fill();
    
}

bool MuonGeometryArrange::isMother(Alignable* ali){
  // Is this the mother ring?
 if(ali==nullptr) return false; // elementary sanity
 const std::vector<Alignable*>& aliComp = ali->components();

 int size=aliComp.size();
 if(size<=0) return false;  // no subcomponents
 
 for(int i=0; i<size; i++){
   if(checkChosen(aliComp[i])) return true; // A ring has CSC chambers
 }                  // as subcomponents
 return false;      // 1'st layer of subcomponents weren't CSC chambers
} 

bool MuonGeometryArrange::checkChosen(Alignable* ali){
 // Check whether the item passed satisfies the criteria given.
  if(ali==nullptr) return false;    // elementary sanity
 // Is this in the CSC section?  If not, bail.  Later may extend.
  if(ali->geomDetId().det()!=DetId::Muon ||
     ali->geomDetId().subdetId()!=MuonSubdetId::CSC) return false;
 // If it is a CSC alignable, then check that the station, etc are
 // those requested.
 // One might think of aligning more than a single ring at a time,
 // by using a vector of ring numbers.  I don't see the sense in
 // trying to align more than one station at a time for comparison.
  CSCDetId cscId(ali->geomDetId());
#ifdef jnbdebug
std::cout<<"JNB "<<ali->id()<<" "<<cscId.endcap()<<" "
<<cscId.station()<<" "<<cscId.ring()<<" "<<cscId.chamber()<<"   "
<<_endcap<<" "<<_station<<" "<<_ring
<<"\n"<<std::flush;
#endif
  if(cscId.endcap()==_endcap && cscId.station()==_station &&
     cscId.ring()==_ring) {
      return true;
  }
  return false;
}
  
bool MuonGeometryArrange::passChosen( Alignable* ali ){

 // Check to see if this contains CSC components of the appropriate ring
 // Ring will contain N Alignables which represent chambers, each of which
 // in turn contains M planes.  For our purposes we don't care about the
 // planes.
 // Hmm.  Interesting question:  Do I want to try to fit the chamber as
 // such, or use the geometry?
 // I want to fit the chamber, so I'll try to use its presence as the marker.
 // What specifically identifies a chamber as a chamber, and not as a layer?
 // The fact that it has layers as sub components, or the fact that it is
 // the first item with a non-zero ID breakdown?  Pick the latter.
 //
 if(ali==nullptr) return false;
 if(checkChosen(ali)) return true;  // If this is one of the desired
                    // CSC chambers, accept it
 const std::vector<Alignable*>& aliComp = ali->components();

 int size=aliComp.size();
 if(size<=0) return false;  // no subcomponents
 
 for(int i=0; i<size; i++){
   if(checkChosen(aliComp[i])) return true; // A ring has CSC chambers
 }                  // as subcomponents
 return false;      // 1'st layer of subcomponents weren't CSC chambers
}
bool MuonGeometryArrange::passIdCut( uint32_t id ){
    
    bool pass = false;
    DetId detid(id);
//  if(detid.det()==DetId::Muon && detid.subdetId()== MuonSubdetId::CSC){
//     CSCDetId cscId(refAli->geomDetId());
//     if(cscId.layer()!=1) return false;       // ONLY FIRST LAYER!
//  }
    int nEntries = _detIdFlagVector.size();
    
    for (int i = 0; i < nEntries; i++){
        if (_detIdFlagVector[i] == id) pass = true;
    }
    
    return pass;
    
}


DEFINE_FWK_MODULE(MuonGeometryArrange);