TrackerGeometryCompare.cc

Go to the documentation of this file.

#include "CondFormats/Alignment/interface/Alignments.h"
#include "CondFormats/Alignment/interface/AlignmentErrors.h"
#include "CondFormats/Alignment/interface/AlignmentSurfaceDeformations.h" 
#include "CondFormats/Alignment/interface/Definitions.h" 
#include "CLHEP/Vector/RotationInterfaces.h" 
#include "CondFormats/Alignment/interface/AlignmentSorter.h"
#include "CondFormats/AlignmentRecord/interface/TrackerSurveyRcd.h"
#include "CondFormats/AlignmentRecord/interface/TrackerSurveyErrorRcd.h"
#include "CondFormats/AlignmentRecord/interface/TrackerAlignmentRcd.h"
#include "CondFormats/AlignmentRecord/interface/TrackerAlignmentErrorRcd.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ESTransientHandle.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "Alignment/CommonAlignment/interface/AlignableObjectId.h"
#include "Geometry/CommonTopologies/interface/SurfaceDeformationFactory.h"
#include "Geometry/CommonTopologies/interface/SurfaceDeformation.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeomBuilderFromGeometricDet.h"
#include "Geometry/TrackerNumberingBuilder/interface/GeometricDet.h" 
#include "Geometry/Records/interface/IdealGeometryRecord.h"
#include "Geometry/TrackingGeometryAligner/interface/GeometryAligner.h"
#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 "CondFormats/AlignmentRecord/interface/GlobalPositionRcd.h"
#include "CondFormats/Alignment/interface/DetectorGlobalPosition.h"
#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"

//#include "Alignment/OfflineValidation/interface/ComparisonUtilities.h"
//#include "Alignment/CommonAlignment/interface/AlignTools.h"

//#include "Alignment/OfflineValidation/plugins/TrackerGeometryCompare.h"
#include "TrackerGeometryCompare.h"
#include "TFile.h" 
#include "CLHEP/Vector/ThreeVector.h"

// Database
#include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"
//#include "Geometry/Records/interface/PGeometricDetRcd.h"

#include "DataFormats/SiStripDetId/interface/TECDetId.h"
#include "DataFormats/SiStripDetId/interface/TIBDetId.h"
#include "DataFormats/SiStripDetId/interface/TOBDetId.h"
#include "DataFormats/SiStripDetId/interface/TIDDetId.h"
#include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
#include "DataFormats/SiPixelDetId/interface/PXBDetId.h"

#include <iostream>
#include <fstream>
#include <sstream> 

TrackerGeometryCompare::TrackerGeometryCompare(const edm::ParameterSet& cfg) :
  referenceTracker(0),
  dummyTracker(0),
  currentTracker(0),
  theSurveyIndex(0),
  theSurveyValues(0),
  theSurveyErrors(0),
  _commonTrackerLevel(align::invalid),
  _inputRootFile1(0),
  _inputRootFile2(0),
  _inputTree01(0),
  _inputTree02(0),
  _inputTree11(0),
  _inputTree12(0),
  m_nBins(10000),
  m_rangeLow(-.1),
  m_rangeHigh(.1), 
  firstEvent_(true),
  m_vtkmap(13)   
{
    
    //input is ROOT
    _inputFilename1 = cfg.getUntrackedParameter< std::string > ("inputROOTFile1");
    _inputFilename2 = cfg.getUntrackedParameter< std::string > ("inputROOTFile2");
    _inputTreenameAlign = cfg.getUntrackedParameter< std::string > ("treeNameAlign");
    _inputTreenameDeform = cfg.getUntrackedParameter< std::string > ("treeNameDeform"); 
    
    //output file
    _filename = cfg.getUntrackedParameter< std::string > ("outputFile");
    
    _writeToDB = cfg.getUntrackedParameter< bool > ("writeToDB" );
    
    const std::vector<std::string>& levels = cfg.getUntrackedParameter< std::vector<std::string> > ("levels");
    
    _weightBy = cfg.getUntrackedParameter< std::string > ("weightBy");
    _setCommonTrackerSystem = cfg.getUntrackedParameter< std::string > ("setCommonTrackerSystem");
    _detIdFlag = cfg.getUntrackedParameter< bool > ("detIdFlag");
    _detIdFlagFile = cfg.getUntrackedParameter< std::string > ("detIdFlagFile");
    _weightById  = cfg.getUntrackedParameter< bool > ("weightById");
    _weightByIdFile = cfg.getUntrackedParameter< std::string > ("weightByIdFile");
    
    //setting the levels being used in the geometry comparator
    AlignableObjectId dummy; //DM_534??
    edm::LogInfo("TrackerGeometryCompare") << "levels: " << levels.size();
    for (unsigned int l = 0; l < levels.size(); ++l){
        //DM_61Xm_theLevels.push_back(AlignableObjectId::stringToId(levels[l])) ; 
        m_theLevels.push_back( dummy.nameToType(levels[l])); //DM_534?? 
        edm::LogInfo("TrackerGeometryCompare") << "level: " << levels[l];
        edm::LogInfo("TrackerGeometryCompare") << "structure type: " << dummy.typeToName(m_theLevels.at(l));  //DM_534?? 
    }
    
        
    // if want to use, make id cut list
    if (_detIdFlag){
        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
    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');
            
            _weightByIdVector.push_back( listId );
        }
        inFile.close();
    }
    
    //root configuration
    _theFile = new TFile(_filename.c_str(),"RECREATE");
    _alignTree = new TTree("alignTree","alignTree");//,"id:level:mid:mlevel:sublevel:x:y:z:r:phi:a:b:c:dx:dy:dz:dr:dphi:da:db:dc");
    _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("du", &_duVal, "du/F");
    _alignTree->Branch("dv", &_dvVal, "dv/F");
    _alignTree->Branch("dw", &_dwVal, "dw/F");
    _alignTree->Branch("da", &_daVal, "da/F");
    _alignTree->Branch("db", &_dbVal, "db/F");
    _alignTree->Branch("dg", &_dgVal, "dg/F");
    _alignTree->Branch("useDetId", &_useDetId, "useDetId/I");
    _alignTree->Branch("detDim", &_detDim, "detDim/I"); 
    _alignTree->Branch("surW", &_surWidth, "surW/F");
    _alignTree->Branch("surL", &_surLength, "surL/F");
    _alignTree->Branch("surRot", &_surRot, "surRot[9]/D");
    _alignTree->Branch("identifiers", &_identifiers, "identifiers[6]/I");
    _alignTree->Branch("type", &_type, "type/I");
        _alignTree->Branch("surfDeform", &_surfDeform, "surfDeform[13]/D"); 

    for (std::vector<TrackerMap>::iterator it = m_vtkmap.begin(); it != m_vtkmap.end(); ++it) {
          it->setPalette(1) ;
      it->addPixel(true) ;
    }

        edm::Service<TFileService> fs;
    TFileDirectory subDir_All = fs->mkdir( "AllSubdetectors" );
    TFileDirectory subDir_PXB = fs->mkdir( "PixelBarrel" );
    TFileDirectory subDir_PXF = fs->mkdir( "PixelEndcap" );
        for (int ii = 0; ii < 13; ++ii) { 
          std::stringstream histname0 ;
          histname0 << "SurfDeform_Par_" << ii ; 
          m_h1[histname0.str()] = subDir_All.make<TH1D>((histname0.str()).c_str(),(histname0.str()).c_str(),m_nBins,m_rangeLow,m_rangeHigh); 
        
          std::stringstream histname1 ;
          histname1 << "SurfDeform_PixelBarrel_Par_" << ii ; 
          m_h1[histname1.str()] = subDir_PXB.make<TH1D>((histname1.str()).c_str(),(histname1.str()).c_str(),m_nBins,m_rangeLow,m_rangeHigh); 
        
          std::stringstream histname2 ;
          histname2 << "SurfDeform_PixelEndcap_Par_" << ii ; 
          m_h1[histname2.str()] = subDir_PXF.make<TH1D>((histname2.str()).c_str(),(histname2.str()).c_str(),m_nBins,m_rangeLow,m_rangeHigh); 
        }
    
}

void TrackerGeometryCompare::beginJob(){
  firstEvent_ = true;
}

void TrackerGeometryCompare::endJob(){

  int iname(0) ;
  for (std::vector<TrackerMap>::iterator it = m_vtkmap.begin(); it != m_vtkmap.end(); ++it) {
    std::stringstream mapname ;
    mapname << "TkMap_SurfDeform" << iname << ".png" ; 
    it->save(true,0,0,mapname.str());
    mapname.str( std::string() ); 
    mapname.clear() ; 
    mapname << "TkMap_SurfDeform" << iname << ".pdf" ; 
    it->save(true,0,0,mapname.str());
    ++iname ; 
  }

  _theFile->cd();
  _alignTree->Write();
  _theFile->Close();
    
}

void TrackerGeometryCompare::analyze(const edm::Event&, const edm::EventSetup& iSetup){

  if (firstEvent_) {

    //upload the ROOT geometries
    createROOTGeometry(iSetup);

    //set common tracker system first
    // if setting the tracker common system
    if (_setCommonTrackerSystem != "NONE"){
        setCommonTrackerSystem();
    }
    
    //compare the goemetries
    compareGeometries(referenceTracker,currentTracker);
    compareSurfaceDeformations(_inputTree11, _inputTree12); 
    
    //write out ntuple
    //might be better to do within output module
    
    if (_writeToDB){
        Alignments* myAlignments = currentTracker->alignments();
        AlignmentErrors* myAlignmentErrors = currentTracker->alignmentErrors();
        
        // 2. Store alignment[Error]s to DB
        edm::Service<cond::service::PoolDBOutputService> poolDbService;
        // Call service
        if( !poolDbService.isAvailable() ) // Die if not available
            throw cms::Exception("NotAvailable") << "PoolDBOutputService not available";
        
        poolDbService->writeOne<Alignments>(&(*myAlignments), poolDbService->beginOfTime(), "TrackerAlignmentRcd");
        poolDbService->writeOne<AlignmentErrors>(&(*myAlignmentErrors), poolDbService->beginOfTime(), "TrackerAlignmentErrorRcd");
        
    }       

    firstEvent_ = false;
  }
}

void TrackerGeometryCompare::createROOTGeometry(const edm::EventSetup& iSetup){
    
    int inputRawId1, inputRawId2;
    double inputX1, inputY1, inputZ1, inputX2, inputY2, inputZ2;
    double inputAlpha1, inputBeta1, inputGamma1, inputAlpha2, inputBeta2, inputGamma2;
        
    //declare alignments
    Alignments* alignments1 = new Alignments();
    AlignmentErrors* alignmentErrors1 = new AlignmentErrors();  
    if (_inputFilename1 != "IDEAL"){
        _inputRootFile1 = new TFile(_inputFilename1.c_str());
        TTree* _inputTree01 = (TTree*) _inputRootFile1->Get(_inputTreenameAlign.c_str());
        _inputTree01->SetBranchAddress("rawid", &inputRawId1);
        _inputTree01->SetBranchAddress("x", &inputX1);
        _inputTree01->SetBranchAddress("y", &inputY1);
        _inputTree01->SetBranchAddress("z", &inputZ1);
        _inputTree01->SetBranchAddress("alpha", &inputAlpha1);
        _inputTree01->SetBranchAddress("beta", &inputBeta1);
        _inputTree01->SetBranchAddress("gamma", &inputGamma1);

        int nEntries1 = _inputTree01->GetEntries();
        //fill alignments
        for (int i = 0; i < nEntries1; ++i){
            
            _inputTree01->GetEntry(i);
            CLHEP::Hep3Vector translation1(inputX1, inputY1, inputZ1);
            CLHEP::HepEulerAngles eulerangles1(inputAlpha1,inputBeta1,inputGamma1);
            uint32_t detid1 = inputRawId1;
            AlignTransform transform1(translation1, eulerangles1, detid1);
            alignments1->m_align.push_back(transform1);
            
            //dummy errors
            CLHEP::HepSymMatrix clhepSymMatrix(3,0);
            AlignTransformError transformError(clhepSymMatrix, detid1);
            alignmentErrors1->m_alignError.push_back(transformError);
        }       
        
        // to get the right order
        std::sort( alignments1->m_align.begin(), alignments1->m_align.end(), lessAlignmentDetId<AlignTransform>() );
        std::sort( alignmentErrors1->m_alignError.begin(), alignmentErrors1->m_alignError.end(), lessAlignmentDetId<AlignTransformError>() );
    }
    //------------------
    Alignments* alignments2 = new Alignments();
    AlignmentErrors* alignmentErrors2 = new AlignmentErrors();
    if (_inputFilename2 != "IDEAL"){    
        _inputRootFile2 = new TFile(_inputFilename2.c_str());
        TTree* _inputTree02 = (TTree*) _inputRootFile2->Get(_inputTreenameAlign.c_str());
        _inputTree02->SetBranchAddress("rawid", &inputRawId2);
        _inputTree02->SetBranchAddress("x", &inputX2);
        _inputTree02->SetBranchAddress("y", &inputY2);
        _inputTree02->SetBranchAddress("z", &inputZ2);
        _inputTree02->SetBranchAddress("alpha", &inputAlpha2);
        _inputTree02->SetBranchAddress("beta", &inputBeta2);
        _inputTree02->SetBranchAddress("gamma", &inputGamma2);
        
        int nEntries2 = _inputTree02->GetEntries();
        //fill alignments
        for (int i = 0; i < nEntries2; ++i){
            
            _inputTree02->GetEntry(i);
            CLHEP::Hep3Vector translation2(inputX2, inputY2, inputZ2);
            CLHEP::HepEulerAngles eulerangles2(inputAlpha2,inputBeta2,inputGamma2);
            uint32_t detid2 = inputRawId2;
            AlignTransform transform2(translation2, eulerangles2, detid2);
            alignments2->m_align.push_back(transform2);
            
            //dummy errors
            CLHEP::HepSymMatrix clhepSymMatrix(3,0);
            AlignTransformError transformError(clhepSymMatrix, detid2);
            alignmentErrors2->m_alignError.push_back(transformError); 
        }           
        
        //to get the right order
        std::sort( alignments2->m_align.begin(), alignments2->m_align.end(), lessAlignmentDetId<AlignTransform>() );
        std::sort( alignmentErrors2->m_alignError.begin(), alignmentErrors2->m_alignError.end(), lessAlignmentDetId<AlignTransformError>() );
    }
    
    //accessing the initial geometry
    edm::ESTransientHandle<DDCompactView> cpv;
    iSetup.get<IdealGeometryRecord>().get(cpv);
    edm::ESHandle<GeometricDet> theGeometricDet;
    iSetup.get<IdealGeometryRecord>().get(theGeometricDet);
    TrackerGeomBuilderFromGeometricDet trackerBuilder;
    
    edm::ESHandle<Alignments> globalPositionRcd;
    iSetup.get<TrackerDigiGeometryRecord>().getRecord<GlobalPositionRcd>().get(globalPositionRcd);
    
    //reference tracker
    TrackerGeometry* theRefTracker = trackerBuilder.build(&*theGeometricDet); //DM_534?? 
    if (_inputFilename1 != "IDEAL"){
        GeometryAligner aligner1;
        aligner1.applyAlignments<TrackerGeometry>( &(*theRefTracker), &(*alignments1), &(*alignmentErrors1),
                                                  align::DetectorGlobalPosition(*globalPositionRcd, DetId(DetId::Tracker)));
    }
    referenceTracker = new AlignableTracker(&(*theRefTracker));
    //referenceTracker->setSurfaceDeformation(surfDef1, true) ; 

    int inputRawid1;
    int inputRawid2;
    int inputDtype1, inputDtype2 ; 
    std::vector<double> inputDpar1;
    std::vector<double> inputDpar2 ; 
    std::vector<double>* p_inputDpar1 = &inputDpar1; 
    std::vector<double>* p_inputDpar2 = &inputDpar2; 

    const std::vector<Alignable*> comp1 = referenceTracker->deepComponents(); 

    SurfaceDeformation * surfDef1; 
    if (_inputFilename1 != "IDEAL"){
      TTree* _inputTree11 = (TTree*) _inputRootFile1->Get(_inputTreenameDeform.c_str());
      _inputTree11->SetBranchAddress("irawid", &inputRawid1);
      _inputTree11->SetBranchAddress("dtype", &inputDtype1);
      _inputTree11->SetBranchAddress("dpar", &p_inputDpar1);

      unsigned int nEntries11 = _inputTree11->GetEntries();
      edm::LogInfo("TrackerGeometryCompare") << " nentries11 = " << nEntries11 << std::endl ; 
      for (unsigned int iEntry = 0; iEntry < nEntries11; ++iEntry) {
            _inputTree11->GetEntry(iEntry) ; 

            surfDef1 = SurfaceDeformationFactory::create( inputDtype1, inputDpar1);

        if (int(comp1[iEntry]->id()) == inputRawid1) {
          comp1[iEntry]->setSurfaceDeformation(surfDef1, true) ; 
        }

      }
    }
        
    //currernt tracker
    TrackerGeometry* theCurTracker = trackerBuilder.build(&*theGeometricDet); //DM_534??
    if (_inputFilename2 != "IDEAL"){
        GeometryAligner aligner2;
        aligner2.applyAlignments<TrackerGeometry>( &(*theCurTracker), &(*alignments2), &(*alignmentErrors2),
                                                  align::DetectorGlobalPosition(*globalPositionRcd, DetId(DetId::Tracker)));
    }
    currentTracker = new AlignableTracker(&(*theCurTracker));
    
    const std::vector<Alignable*> comp2 = currentTracker->deepComponents(); 

    SurfaceDeformation * surfDef2 ; 
    if (_inputFilename2 != "IDEAL"){ 
      TTree* _inputTree12 = (TTree*) _inputRootFile2->Get(_inputTreenameDeform.c_str());
      _inputTree12->SetBranchAddress("irawid", &inputRawid2);
      _inputTree12->SetBranchAddress("dtype", &inputDtype2);
      _inputTree12->SetBranchAddress("dpar",  &p_inputDpar2);

      unsigned int nEntries12 = _inputTree12->GetEntries();
      edm::LogInfo("TrackerGeometryCompare") << " nentries12 = " << nEntries12 << std::endl ; 
      for (unsigned int iEntry = 0; iEntry < nEntries12; ++iEntry) {
            _inputTree12->GetEntry(iEntry) ; 
        
            surfDef2 = SurfaceDeformationFactory::create( inputDtype2, inputDpar2);

        if (int(comp2[iEntry]->id()) == inputRawid2) {
          comp2[iEntry]->setSurfaceDeformation(surfDef2, true) ; 
        }

      }
    }
        
    delete alignments1;
    delete alignmentErrors1;
    delete alignments2;
    delete alignmentErrors2;

}

void TrackerGeometryCompare::compareSurfaceDeformations(TTree* refTree, TTree* curTree) {
    
  if (_inputFilename1 != "IDEAL" && _inputFilename2 != "IDEAL") {
             
    int inputRawid1;
    int inputRawid2;
    int inputSubdetid1, inputSubdetid2 ; 
    int inputDtype1, inputDtype2 ; 
    std::vector<double> inputDpar1;
    std::vector<double> inputDpar2 ; 
    std::vector<double>* p_inputDpar1 = &inputDpar1; 
    std::vector<double>* p_inputDpar2 = &inputDpar2; 
  
    TTree* refTree = (TTree*) _inputRootFile1->Get(_inputTreenameDeform.c_str());
    refTree->SetBranchAddress("irawid", &inputRawid1);
    refTree->SetBranchAddress("subdetid", &inputSubdetid1);
    refTree->SetBranchAddress("dtype", &inputDtype1);
    refTree->SetBranchAddress("dpar", &p_inputDpar1);
  
    TTree* curTree = (TTree*) _inputRootFile2->Get(_inputTreenameDeform.c_str());
    curTree->SetBranchAddress("irawid", &inputRawid2);
    curTree->SetBranchAddress("subdetid", &inputSubdetid2);
    curTree->SetBranchAddress("dtype", &inputDtype2);
    curTree->SetBranchAddress("dpar",  &p_inputDpar2);
  
    unsigned int nEntries11 = refTree->GetEntries();
    unsigned int nEntries12 = curTree->GetEntries();

    if (nEntries11 != nEntries12) {
      edm::LogError("TrackerGeometryCompare")   << " Surface deformation parameters in two geometries differ!\n" ;
      return ; 
    }
    
    for (unsigned int iEntry = 0; iEntry < nEntries12; ++iEntry) {
      refTree->GetEntry(iEntry) ;
      curTree->GetEntry(iEntry) ;
      for (int ii = 0; ii < 13; ++ii) { _surfDeform[ii] = -1.0 ; } 
      for (int npar = 0; npar < int(inputDpar2.size()); ++npar ) {
        if (inputRawid1 == inputRawid2) {
        _surfDeform[npar] = inputDpar2.at(npar) - inputDpar1.at(npar) ; 
        std::stringstream histname0 ;
        histname0 << "SurfDeform_Par_" << npar ;  
    if ( TMath::Abs(_surfDeform[npar]) > (m_rangeHigh - m_rangeLow)/(10.*m_nBins) ) m_h1[histname0.str()]->Fill(_surfDeform[npar]) ; 
    if (inputSubdetid1 == 1 && inputSubdetid2 == 1) {
          std::stringstream histname1 ;
          histname1 << "SurfDeform_PixelBarrel_Par_" << npar ;  
          if ( TMath::Abs(_surfDeform[npar]) > (m_rangeHigh - m_rangeLow)/(10.*m_nBins) ) m_h1[histname1.str()]->Fill(_surfDeform[npar]) ; 
    }
    if (inputSubdetid1 == 2 && inputSubdetid2 == 2) {
          std::stringstream histname2 ;
          histname2 << "SurfDeform_PixelEndcap_Par_" << npar ;  
          if ( TMath::Abs(_surfDeform[npar]) > (m_rangeHigh - m_rangeLow)/(10.*m_nBins) ) m_h1[histname2.str()]->Fill(_surfDeform[npar]) ; 
    }
        (m_vtkmap.at(npar)).fill_current_val(inputRawid1,_surfDeform[npar]) ; 
        }
      }
    }
  
  } else if ( _inputFilename1 == "IDEAL" && _inputFilename2 != "IDEAL" ) {
             
    int inputRawid2;
    int inputSubdetid2 ; 
    int inputDtype2 ; 
    std::vector<double> inputDpar2 ; 
    std::vector<double>* p_inputDpar2 = &inputDpar2; 
  
    TTree* curTree = (TTree*) _inputRootFile2->Get(_inputTreenameDeform.c_str());
    curTree->SetBranchAddress("irawid", &inputRawid2);
    curTree->SetBranchAddress("subdetid", &inputSubdetid2);
    curTree->SetBranchAddress("dtype", &inputDtype2);
    curTree->SetBranchAddress("dpar",  &p_inputDpar2);
  
    unsigned int nEntries12 = curTree->GetEntries();
    
    for (unsigned int iEntry = 0; iEntry < nEntries12; ++iEntry) {
      curTree->GetEntry(iEntry) ;
      for (int ii = 0; ii < 12; ++ii) { _surfDeform[ii] = -1.0 ; } 
      for (int npar = 0; npar < int(inputDpar2.size()); ++npar ) {
        _surfDeform[npar] = inputDpar2.at(npar) ; 
        std::stringstream histname0 ;
        histname0 << "SurfDeform_Par_" << npar ;  
        if ( TMath::Abs(_surfDeform[npar]) > (m_rangeHigh - m_rangeLow)/(10.*m_nBins) ) m_h1[histname0.str()]->Fill(_surfDeform[npar]) ; 
    if (inputSubdetid2 == 1) {
          std::stringstream histname1 ;
          histname1 << "SurfDeform_PixelBarrel_Par_" << npar ;  
          if ( TMath::Abs(_surfDeform[npar]) > (m_rangeHigh - m_rangeLow)/(10.*m_nBins) ) m_h1[histname1.str()]->Fill(_surfDeform[npar]) ; 
    }
    if (inputSubdetid2 == 2) {
          std::stringstream histname2 ;
          histname2 << "SurfDeform_PixelEndcap_Par_" << npar ;  
          if ( TMath::Abs(_surfDeform[npar]) > (m_rangeHigh - m_rangeLow)/(10.*m_nBins) ) m_h1[histname2.str()]->Fill(_surfDeform[npar]) ; 
    }
        (m_vtkmap.at(npar)).fill_current_val(inputRawid2,_surfDeform[npar]) ; 
      }
    }
  
  } else if ( _inputFilename1 != "IDEAL" && _inputFilename2 == "IDEAL" ) {
             
    int inputRawid1;
    int inputSubdetid1 ; 
    int inputDtype1 ; 
    std::vector<double> inputDpar1;
    std::vector<double>* p_inputDpar1 = &inputDpar1; 
  
    TTree* refTree = (TTree*) _inputRootFile1->Get(_inputTreenameDeform.c_str());
    refTree->SetBranchAddress("irawid", &inputRawid1);
    refTree->SetBranchAddress("subdetid", &inputSubdetid1);
    refTree->SetBranchAddress("dtype", &inputDtype1);
    refTree->SetBranchAddress("dpar", &p_inputDpar1);
  
    unsigned int nEntries11 = refTree->GetEntries();
    
    for (unsigned int iEntry = 0; iEntry < nEntries11; ++iEntry) {
      refTree->GetEntry(iEntry) ;
      for (int ii = 0; ii < 12; ++ii) { _surfDeform[ii] = -1.0 ; } 
      for (int npar = 0; npar < int(inputDpar1.size()); ++npar ) {
        _surfDeform[npar] = - inputDpar1.at(npar) ; 
        std::stringstream histname0 ;
        histname0 << "SurfDeform_Par_" << npar ;  
        if ( TMath::Abs(_surfDeform[npar]) > (m_rangeHigh - m_rangeLow)/(10.*m_nBins) ) m_h1[histname0.str()]->Fill(_surfDeform[npar]) ; 
    if (inputSubdetid1 == 1) {
          std::stringstream histname1 ;
          histname1 << "SurfDeform_PixelBarrel_Par_" << npar ;  
          if ( TMath::Abs(_surfDeform[npar]) > (m_rangeHigh - m_rangeLow)/(10.*m_nBins) ) m_h1[histname1.str()]->Fill(_surfDeform[npar]) ; 
    }
    if (inputSubdetid1 == 2) {
          std::stringstream histname2 ;
          histname2 << "SurfDeform_PixelEndcap_Par_" << npar ;  
          if ( TMath::Abs(_surfDeform[npar]) > (m_rangeHigh - m_rangeLow)/(10.*m_nBins) ) m_h1[histname2.str()]->Fill(_surfDeform[npar]) ; 
    }
        (m_vtkmap.at(npar)).fill_current_val(inputRawid1,_surfDeform[npar]) ; 
      }
    }
  
  } else if ( _inputFilename1 == "IDEAL" && _inputFilename2 == "IDEAL" ) {

      edm::LogInfo("TrackerGeometryCompare") << ">>>> Comparing IDEAL with IDEAL: nothing to do! <<<<\n" ; 
      
  }

  return ;  
}

void TrackerGeometryCompare::compareGeometries(Alignable* refAli, Alignable* curAli){

    using namespace align ; 
    
    const std::vector<Alignable*>& refComp = refAli->components();
    const std::vector<Alignable*>& curComp = curAli->components();
    
    unsigned int nComp = refComp.size();
    //only perform for designate levels
    bool useLevel = false;
    for (unsigned int i = 0; i < m_theLevels.size(); ++i){
        if (refAli->alignableObjectId() == m_theLevels[i]) useLevel = true;
    }
    
    //another added level for difference between det and detunit
    //if ((refAli->alignableObjectId()==2)&&(nComp == 1)) useLevel = false;
    
    //coordinate matching, etc etc
    if (useLevel){
        DetId detid(refAli->id());

        CLHEP::Hep3Vector Rtotal, Wtotal, lRtotal, lWtotal;
        Rtotal.set(0.,0.,0.); 
        Wtotal.set(0.,0.,0.);
        lRtotal.set(0.,0.,0.); 
        lWtotal.set(0.,0.,0.);

        for (int i = 0; i < 100; i++){
            AlgebraicVector diff = align::diffAlignables(refAli,curAli, _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());
            // local coordinates
            lRtotal.set(diff[6],diff[7],diff[8]);
            lWtotal.set(diff[9],diff[10],diff[11]);
            align::moveAlignable(curAli, diff);
            float tolerance = 1e-7;
            AlgebraicVector check = align::diffAlignables(refAli,curAli, _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)){
                edm::LogInfo("TrackerGeometryCompare") << "Tolerance Exceeded!(alObjId: " << refAli->alignableObjectId()
                << ", rawId: " << refAli->geomDetId().rawId()
                << ", subdetId: "<< detid.subdetId() << "): " << diff;
                throw cms::Exception("Tolerance in TrackerGeometryCompare exceeded");
            }
            else{
                break;
            }
        }

        AlgebraicVector TRtot(12);
        // global 
        TRtot(1) = Rtotal.x(); TRtot(2) = Rtotal.y(); TRtot(3) = Rtotal.z();
        TRtot(4) = Wtotal.x(); TRtot(5) = Wtotal.y(); TRtot(6) = Wtotal.z();
        // local
        TRtot(7) = lRtotal.x(); TRtot(8) = lRtotal.y(); TRtot(9) = lRtotal.z();
        TRtot(10) = lWtotal.x(); TRtot(11) = lWtotal.y(); TRtot(12) = lWtotal.z();

        fillTree(refAli, TRtot);
    }

    // another added level for difference between det and detunit
    for (unsigned int i = 0; i < nComp; ++i) 
      compareGeometries(refComp[i],curComp[i]); 

}

void TrackerGeometryCompare::setCommonTrackerSystem(){

    edm::LogInfo("TrackerGeometryCompare") << "Setting Common Tracker System....";
    
    AlignableObjectId dummy; // DM_534??
    _commonTrackerLevel = dummy.nameToType(_setCommonTrackerSystem); // DM_534?? 
        
    diffCommonTrackerSystem(referenceTracker, currentTracker);
    
    align::EulerAngles dOmega(3); dOmega[0] = _TrackerCommonR.x() ; dOmega[1] = _TrackerCommonR.y(); dOmega[2] = _TrackerCommonR.z();
    align::RotationType rot = align::toMatrix( dOmega );
    align::GlobalVector theR = _TrackerCommonT;
    
    std::cout << "what we get from overlaying the pixels..." << theR << ", " << rot << std::endl;
    
    //transform to the Tracker System
    align::PositionType trackerCM = currentTracker->globalPosition();
    align::GlobalVector cmDiff( trackerCM.x()-_TrackerCommonCM.x(), trackerCM.y()-_TrackerCommonCM.y(), trackerCM.z()-_TrackerCommonCM.z() );
    
    std::cout << "Pixel CM: " << _TrackerCommonCM << ", tracker CM: " << trackerCM << std::endl;
    
    //adjust translational difference factoring in different rotational CM
    //needed because rotateInGlobalFrame is about CM of alignable, not Tracker
    align::GlobalVector::BasicVectorType lpvgf = cmDiff.basicVector();
    align::GlobalVector moveV( rot.multiplyInverse(lpvgf) - lpvgf);
    align::GlobalVector theRprime(theR + moveV);
    
    AlgebraicVector TrackerCommonTR(6);
    TrackerCommonTR(1) = theRprime.x(); TrackerCommonTR(2) = theRprime.y(); TrackerCommonTR(3) = theRprime.z();
    TrackerCommonTR(4) = _TrackerCommonR.x(); TrackerCommonTR(5) = _TrackerCommonR.y(); TrackerCommonTR(6) = _TrackerCommonR.z();
    
    std::cout << "and after the transformation: " << TrackerCommonTR << std::endl;
    
    align::moveAlignable(currentTracker, TrackerCommonTR );
    
}

void TrackerGeometryCompare::diffCommonTrackerSystem(Alignable *refAli, Alignable *curAli){
    
    const std::vector<Alignable*>& refComp = refAli->components();
    const std::vector<Alignable*>& curComp = curAli->components();
    
    unsigned int nComp = refComp.size();
    //only perform for designate levels
    bool useLevel = false;
    if (refAli->alignableObjectId() == _commonTrackerLevel) useLevel = true;
    
    //useLevel = false;
    if (useLevel){
        CLHEP::Hep3Vector Rtotal, Wtotal;
        Rtotal.set(0.,0.,0.); Wtotal.set(0.,0.,0.);
        
        AlgebraicVector diff = align::diffAlignables(refAli,curAli, _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());
        /*
         //std::cout << "a";
         //if (refAli->alignableObjectId() == 1) std::cout << "DIFF: " << diff << std::endl;
         align::moveAlignable(curAli, diff);
         float tolerance = 1e-7;
         AlgebraicVector check = align::diffAlignables(refAli,curAli, _weightBy, _weightById, _weightByIdVector);
         align::GlobalVector checkR(check[0],check[1],check[2]);
         align::GlobalVector checkW(check[3],check[4],check[5]);
         DetId detid(refAli->id());
         if ((checkR.mag() > tolerance)||(checkW.mag() > tolerance)){
         edm::LogInfo("TrackerGeometryCompare") << "Tolerance Exceeded!(alObjId: " << refAli->alignableObjectId()
         << ", rawId: " << refAli->geomDetId().rawId()
         << ", subdetId: "<< detid.subdetId() << "): " << diff;
         }
         else{
         break;
         }
         }
         */
        
        //_TrackerCommonT.set(Rtotal.x(), Rtotal.y(), Rtotal.z());
        _TrackerCommonT = align::GlobalVector(Rtotal.x(), Rtotal.y(), Rtotal.z());
        _TrackerCommonR = align::GlobalVector(Wtotal.x(), Wtotal.y(), Wtotal.z());
        _TrackerCommonCM = curAli->globalPosition();
        //_TrackerCommonTR(1) = Rtotal.x(); _TrackerCommonTR(2) = Rtotal.y(); _TrackerCommonTR(3) = Rtotal.z();
        //_TrackerCommonTR(4) = Wtotal.x(); _TrackerCommonTR(5) = Wtotal.y(); _TrackerCommonTR(6) = Wtotal.z();
        
        
    }
    else{
        for (unsigned int i = 0; i < nComp; ++i) diffCommonTrackerSystem(refComp[i],curComp[i]);
    }
    
    
}

void TrackerGeometryCompare::fillTree(Alignable *refAli, 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();
    fillIdentifiers( _sublevel, _id );
    _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);
    _alphaVal = eulerAngles[0];
    _betaVal = eulerAngles[1];
    _gammaVal = eulerAngles[2];
    // global
    _dxVal = diff[0];
    _dyVal = diff[1];
    _dzVal = diff[2];
    // local
    _duVal = diff[6];
    _dvVal = diff[7];
    _dwVal = diff[8];
    //...TODO...
    align::GlobalVector g(_dxVal, _dyVal, _dzVal);
    //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();
    // global
    _dalphaVal = diff[3];
    _dbetaVal = diff[4];
    _dgammaVal = diff[5];
    // local
    _daVal = diff[9];
    _dbVal = diff[10];
    _dgVal = diff[11];
    
    //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();
    
    //Fill
    _alignTree->Fill();
    
}

void TrackerGeometryCompare::surveyToTracker(AlignableTracker* ali, Alignments* alignVals, AlignmentErrors* alignErrors){
    
    //getting the right alignables for the alignment record
    std::vector<Alignable*> detPB = ali->pixelHalfBarrelGeomDets();
    std::vector<Alignable*> detPEC = ali->pixelEndcapGeomDets();
    std::vector<Alignable*> detTIB = ali->innerBarrelGeomDets();
    std::vector<Alignable*> detTID = ali->TIDGeomDets();
    std::vector<Alignable*> detTOB = ali->outerBarrelGeomDets();
    std::vector<Alignable*> detTEC = ali->endcapGeomDets();
    
    std::vector<Alignable*> allGeomDets;
    std::copy(detPB.begin(), detPB.end(), std::back_inserter(allGeomDets));
    std::copy(detPEC.begin(), detPEC.end(), std::back_inserter(allGeomDets));
    std::copy(detTIB.begin(), detTIB.end(), std::back_inserter(allGeomDets));
    std::copy(detTID.begin(), detTID.end(), std::back_inserter(allGeomDets));
    std::copy(detTOB.begin(), detTOB.end(), std::back_inserter(allGeomDets));
    std::copy(detTEC.begin(), detTEC.end(), std::back_inserter(allGeomDets));
    
    std::vector<Alignable*> rcdAlis;
    for (std::vector<Alignable*>::iterator i = allGeomDets.begin(); i!= allGeomDets.end(); i++){
        if ((*i)->components().size() == 1){
            rcdAlis.push_back((*i));
        }
        else if ((*i)->components().size() > 1){
            rcdAlis.push_back((*i));
            std::vector<Alignable*> comp = (*i)->components();
            for (std::vector<Alignable*>::iterator j = comp.begin(); j != comp.end(); j++){
                rcdAlis.push_back((*j));
            }
        }
    }
    
    //turning them into alignments
    for(std::vector<Alignable*>::iterator k = rcdAlis.begin(); k != rcdAlis.end(); k++){
        
        const SurveyDet* surveyInfo = (*k)->survey();
        align::PositionType pos(surveyInfo->position());
        align::RotationType rot(surveyInfo->rotation());
        CLHEP::Hep3Vector clhepVector(pos.x(),pos.y(),pos.z());
        CLHEP::HepRotation clhepRotation( CLHEP::HepRep3x3(rot.xx(),rot.xy(),rot.xz(),rot.yx(),rot.yy(),rot.yz(),rot.zx(),rot.zy(),rot.zz()));
        AlignTransform transform(clhepVector, clhepRotation, (*k)->id());
        AlignTransformError transformError(CLHEP::HepSymMatrix(3,1), (*k)->id());
        alignVals->m_align.push_back(transform);
        alignErrors->m_alignError.push_back(transformError);
    }
    
    //to get the right order
    std::sort( alignVals->m_align.begin(), alignVals->m_align.end(), lessAlignmentDetId<AlignTransform>() );
    std::sort( alignErrors->m_alignError.begin(), alignErrors->m_alignError.end(), lessAlignmentDetId<AlignTransformError>() );
    
}

void TrackerGeometryCompare::addSurveyInfo(Alignable* ali){
    
    const std::vector<Alignable*>& comp = ali->components();
    
    unsigned int nComp = comp.size();
    
    for (unsigned int i = 0; i < nComp; ++i) addSurveyInfo(comp[i]);
    
    const SurveyError& error = theSurveyErrors->m_surveyErrors[theSurveyIndex];
    
    if ( ali->geomDetId().rawId() != error.rawId() ||
        ali->alignableObjectId() != error.structureType() )
    {
        throw cms::Exception("DatabaseError")
        << "Error reading survey info from DB. Mismatched id!";
    }
    
    const CLHEP::Hep3Vector&  pos = theSurveyValues->m_align[theSurveyIndex].translation();
    const CLHEP::HepRotation& rot = theSurveyValues->m_align[theSurveyIndex].rotation();
    
    AlignableSurface surf( align::PositionType( pos.x(), pos.y(), pos.z() ),
                          align::RotationType( rot.xx(), rot.xy(), rot.xz(),
                                              rot.yx(), rot.yy(), rot.yz(),
                                              rot.zx(), rot.zy(), rot.zz() ) );
    
    surf.setWidth( ali->surface().width() );
    surf.setLength( ali->surface().length() );
    
    ali->setSurvey( new SurveyDet( surf, error.matrix() ) );
    
    ++theSurveyIndex;
    
}

bool TrackerGeometryCompare::passIdCut( uint32_t id ){
    
    bool pass = false;
    int nEntries = _detIdFlagVector.size();
    
    for (int i = 0; i < nEntries; i++){
        if (_detIdFlagVector[i] == id) pass = true;
    }
    
    return pass;
    
}

void TrackerGeometryCompare::fillIdentifiers( int subdetlevel, int rawid ){
    
    
    switch( subdetlevel ){
            
                case 1:
            {
            PXBDetId pxbid( rawid );
            _identifiers[0] = pxbid.module();
            _identifiers[1] = pxbid.ladder();
            _identifiers[2] = pxbid.layer();
            _identifiers[3] = 999;
            _identifiers[4] = 999;
            _identifiers[5] = 999;
            break;
        }
        case 2:
        {
            PXFDetId pxfid( rawid );
            _identifiers[0] = pxfid.module();
            _identifiers[1] = pxfid.panel();
            _identifiers[2] = pxfid.blade();
            _identifiers[3] = pxfid.disk();
            _identifiers[4] = pxfid.side();
            _identifiers[5] = 999;
            break;
        }
        case 3:
        {
            TIBDetId tibid( rawid );
            _identifiers[0] = tibid.module();
            _identifiers[1] = tibid.string()[0];
            _identifiers[2] = tibid.string()[1];
            _identifiers[3] = tibid.string()[2];
            _identifiers[4] = tibid.layer();
            _identifiers[5] = 999;
            break;
        }
        case 4: 
        {
            TIDDetId tidid( rawid );
            _identifiers[0] = tidid.module()[0];
            _identifiers[1] = tidid.module()[1];
            _identifiers[2] = tidid.ring();
            _identifiers[3] = tidid.wheel();
            _identifiers[4] = tidid.side();
            _identifiers[5] = 999;
            break;
        }
        case 5: 
        {
            TOBDetId tobid( rawid );
            _identifiers[0] = tobid.module();
            _identifiers[1] = tobid.rod()[0];
            _identifiers[2] = tobid.rod()[1];
            _identifiers[3] = tobid.layer();
            _identifiers[4] = 999;
            _identifiers[5] = 999;
            break;
        }
        case 6: 
        {
            TECDetId tecid( rawid );
            _identifiers[0] = tecid.module();
            _identifiers[1] = tecid.ring();
            _identifiers[2] = tecid.petal()[0];
            _identifiers[3] = tecid.petal()[1];
            _identifiers[4] = tecid.wheel();
            _identifiers[5] = tecid.side();
            break;
        }
        default:
        {
            std::cout << "Error: bad subdetid!!" << std::endl;
            break;
        }
            
    }
}


DEFINE_FWK_MODULE(TrackerGeometryCompare);