CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_14/src/Alignment/OfflineValidation/plugins/TrackerGeometryCompare.cc

Go to the documentation of this file.
00001 #include "CondFormats/Alignment/interface/Alignments.h"
00002 #include "CondFormats/Alignment/interface/AlignmentErrors.h"
00003 #include "CondFormats/Alignment/interface/AlignmentSurfaceDeformations.h" 
00004 #include "CondFormats/Alignment/interface/Definitions.h" 
00005 #include "CLHEP/Vector/RotationInterfaces.h" 
00006 #include "CondFormats/Alignment/interface/AlignmentSorter.h"
00007 #include "CondFormats/AlignmentRecord/interface/TrackerSurveyRcd.h"
00008 #include "CondFormats/AlignmentRecord/interface/TrackerSurveyErrorRcd.h"
00009 #include "CondFormats/AlignmentRecord/interface/TrackerAlignmentRcd.h"
00010 #include "CondFormats/AlignmentRecord/interface/TrackerAlignmentErrorRcd.h"
00011 #include "FWCore/Framework/interface/MakerMacros.h"
00012 #include "FWCore/Framework/interface/ESHandle.h"
00013 #include "FWCore/Framework/interface/EventSetup.h"
00014 #include "FWCore/Framework/interface/ESTransientHandle.h"
00015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00016 
00017 #include "Alignment/CommonAlignment/interface/AlignableObjectId.h"
00018 #include "Geometry/CommonTopologies/interface/SurfaceDeformationFactory.h"
00019 #include "Geometry/CommonTopologies/interface/SurfaceDeformation.h"
00020 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00021 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeomBuilderFromGeometricDet.h"
00022 #include "Geometry/TrackerNumberingBuilder/interface/GeometricDet.h" 
00023 #include "Geometry/Records/interface/IdealGeometryRecord.h"
00024 #include "Geometry/TrackingGeometryAligner/interface/GeometryAligner.h"
00025 #include "Alignment/CommonAlignment/interface/Utilities.h"
00026 #include "Alignment/CommonAlignment/interface/SurveyDet.h"
00027 #include "Alignment/CommonAlignment/interface/Alignable.h"
00028 #include "DataFormats/DetId/interface/DetId.h"
00029 #include "CondFormats/AlignmentRecord/interface/GlobalPositionRcd.h"
00030 #include "CondFormats/Alignment/interface/DetectorGlobalPosition.h"
00031 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00032 
00033 //#include "Alignment/OfflineValidation/interface/ComparisonUtilities.h"
00034 //#include "Alignment/CommonAlignment/interface/AlignTools.h"
00035 
00036 //#include "Alignment/OfflineValidation/plugins/TrackerGeometryCompare.h"
00037 #include "TrackerGeometryCompare.h"
00038 #include "TFile.h" 
00039 #include "CLHEP/Vector/ThreeVector.h"
00040 
00041 // Database
00042 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
00043 #include "FWCore/ServiceRegistry/interface/Service.h"
00044 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00045 //#include "Geometry/Records/interface/PGeometricDetRcd.h"
00046 
00047 #include "DataFormats/SiStripDetId/interface/TECDetId.h"
00048 #include "DataFormats/SiStripDetId/interface/TIBDetId.h"
00049 #include "DataFormats/SiStripDetId/interface/TOBDetId.h"
00050 #include "DataFormats/SiStripDetId/interface/TIDDetId.h"
00051 #include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
00052 #include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
00053 
00054 #include <iostream>
00055 #include <fstream>
00056 #include <sstream> 
00057 
00058 TrackerGeometryCompare::TrackerGeometryCompare(const edm::ParameterSet& cfg) :
00059   referenceTracker(0),
00060   dummyTracker(0),
00061   currentTracker(0),
00062   theSurveyIndex(0),
00063   theSurveyValues(0),
00064   theSurveyErrors(0),
00065   _commonTrackerLevel(align::invalid),
00066   _inputRootFile1(0),
00067   _inputRootFile2(0),
00068   _inputTree01(0),
00069   _inputTree02(0),
00070   _inputTree11(0),
00071   _inputTree12(0),
00072   m_nBins(10000),
00073   m_rangeLow(-.1),
00074   m_rangeHigh(.1), 
00075   firstEvent_(true),
00076   m_vtkmap(13)   
00077 {
00078         
00079         //input is ROOT
00080         _inputFilename1 = cfg.getUntrackedParameter< std::string > ("inputROOTFile1");
00081         _inputFilename2 = cfg.getUntrackedParameter< std::string > ("inputROOTFile2");
00082         _inputTreenameAlign = cfg.getUntrackedParameter< std::string > ("treeNameAlign");
00083         _inputTreenameDeform = cfg.getUntrackedParameter< std::string > ("treeNameDeform"); 
00084         
00085         //output file
00086         _filename = cfg.getUntrackedParameter< std::string > ("outputFile");
00087         
00088         _writeToDB = cfg.getUntrackedParameter< bool > ("writeToDB" );
00089         
00090         const std::vector<std::string>& levels = cfg.getUntrackedParameter< std::vector<std::string> > ("levels");
00091         
00092         _weightBy = cfg.getUntrackedParameter< std::string > ("weightBy");
00093         _setCommonTrackerSystem = cfg.getUntrackedParameter< std::string > ("setCommonTrackerSystem");
00094         _detIdFlag = cfg.getUntrackedParameter< bool > ("detIdFlag");
00095         _detIdFlagFile = cfg.getUntrackedParameter< std::string > ("detIdFlagFile");
00096         _weightById  = cfg.getUntrackedParameter< bool > ("weightById");
00097         _weightByIdFile = cfg.getUntrackedParameter< std::string > ("weightByIdFile");
00098         
00099         //setting the levels being used in the geometry comparator
00100         AlignableObjectId dummy; //DM_534??
00101         edm::LogInfo("TrackerGeometryCompare") << "levels: " << levels.size();
00102         for (unsigned int l = 0; l < levels.size(); ++l){
00103                 //DM_61Xm_theLevels.push_back(AlignableObjectId::stringToId(levels[l])) ; 
00104                 m_theLevels.push_back( dummy.nameToType(levels[l])); //DM_534?? 
00105                 edm::LogInfo("TrackerGeometryCompare") << "level: " << levels[l];
00106                 edm::LogInfo("TrackerGeometryCompare") << "structure type: " << dummy.typeToName(m_theLevels.at(l));  //DM_534?? 
00107         }
00108         
00109                 
00110         // if want to use, make id cut list
00111         if (_detIdFlag){
00112         ifstream fin;
00113         fin.open( _detIdFlagFile.c_str() );
00114         
00115         while (!fin.eof() && fin.good() ){
00116                         
00117                         uint32_t id;
00118                         fin >> id;
00119                         _detIdFlagVector.push_back(id);
00120         }
00121         fin.close();
00122         }               
00123         
00124         // turn weightByIdFile into weightByIdVector
00125         if (_weightById){
00126                 std::ifstream inFile;
00127                 inFile.open( _weightByIdFile.c_str() );
00128                 int ctr = 0;
00129                 while ( !inFile.eof() ){
00130                         ctr++;
00131                         unsigned int listId;
00132                         inFile >> listId;
00133                         inFile.ignore(256, '\n');
00134                         
00135                         _weightByIdVector.push_back( listId );
00136                 }
00137                 inFile.close();
00138         }
00139         
00140         //root configuration
00141         _theFile = new TFile(_filename.c_str(),"RECREATE");
00142         _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");
00143         _alignTree->Branch("id", &_id, "id/I");
00144         _alignTree->Branch("level", &_level, "level/I");
00145         _alignTree->Branch("mid", &_mid, "mid/I");
00146         _alignTree->Branch("mlevel", &_mlevel, "mlevel/I");
00147         _alignTree->Branch("sublevel", &_sublevel, "sublevel/I");
00148         _alignTree->Branch("x", &_xVal, "x/F");
00149         _alignTree->Branch("y", &_yVal, "y/F");
00150         _alignTree->Branch("z", &_zVal, "z/F");
00151         _alignTree->Branch("r", &_rVal, "r/F");
00152         _alignTree->Branch("phi", &_phiVal, "phi/F");
00153         _alignTree->Branch("eta", &_etaVal, "eta/F");
00154         _alignTree->Branch("alpha", &_alphaVal, "alpha/F");
00155         _alignTree->Branch("beta", &_betaVal, "beta/F");
00156         _alignTree->Branch("gamma", &_gammaVal, "gamma/F");
00157         _alignTree->Branch("dx", &_dxVal, "dx/F");
00158         _alignTree->Branch("dy", &_dyVal, "dy/F");
00159         _alignTree->Branch("dz", &_dzVal, "dz/F");
00160         _alignTree->Branch("dr", &_drVal, "dr/F");
00161         _alignTree->Branch("dphi", &_dphiVal, "dphi/F");
00162         _alignTree->Branch("dalpha", &_dalphaVal, "dalpha/F");
00163         _alignTree->Branch("dbeta", &_dbetaVal, "dbeta/F");
00164         _alignTree->Branch("dgamma", &_dgammaVal, "dgamma/F");
00165         _alignTree->Branch("du", &_duVal, "du/F");
00166         _alignTree->Branch("dv", &_dvVal, "dv/F");
00167         _alignTree->Branch("dw", &_dwVal, "dw/F");
00168         _alignTree->Branch("da", &_daVal, "da/F");
00169         _alignTree->Branch("db", &_dbVal, "db/F");
00170         _alignTree->Branch("dg", &_dgVal, "dg/F");
00171         _alignTree->Branch("useDetId", &_useDetId, "useDetId/I");
00172         _alignTree->Branch("detDim", &_detDim, "detDim/I");     
00173         _alignTree->Branch("surW", &_surWidth, "surW/F");
00174         _alignTree->Branch("surL", &_surLength, "surL/F");
00175         _alignTree->Branch("surRot", &_surRot, "surRot[9]/D");
00176         _alignTree->Branch("identifiers", &_identifiers, "identifiers[6]/I");
00177         _alignTree->Branch("type", &_type, "type/I");
00178         _alignTree->Branch("surfDeform", &_surfDeform, "surfDeform[13]/D"); 
00179 
00180         for (std::vector<TrackerMap>::iterator it = m_vtkmap.begin(); it != m_vtkmap.end(); ++it) {
00181           it->setPalette(1) ;
00182           it->addPixel(true) ;
00183         }
00184 
00185         edm::Service<TFileService> fs;
00186         TFileDirectory subDir_All = fs->mkdir( "AllSubdetectors" );
00187         TFileDirectory subDir_PXB = fs->mkdir( "PixelBarrel" );
00188         TFileDirectory subDir_PXF = fs->mkdir( "PixelEndcap" );
00189         for (int ii = 0; ii < 13; ++ii) { 
00190           std::stringstream histname0 ;
00191           histname0 << "SurfDeform_Par_" << ii ; 
00192           m_h1[histname0.str()] = subDir_All.make<TH1D>((histname0.str()).c_str(),(histname0.str()).c_str(),m_nBins,m_rangeLow,m_rangeHigh); 
00193         
00194           std::stringstream histname1 ;
00195           histname1 << "SurfDeform_PixelBarrel_Par_" << ii ; 
00196           m_h1[histname1.str()] = subDir_PXB.make<TH1D>((histname1.str()).c_str(),(histname1.str()).c_str(),m_nBins,m_rangeLow,m_rangeHigh); 
00197         
00198           std::stringstream histname2 ;
00199           histname2 << "SurfDeform_PixelEndcap_Par_" << ii ; 
00200           m_h1[histname2.str()] = subDir_PXF.make<TH1D>((histname2.str()).c_str(),(histname2.str()).c_str(),m_nBins,m_rangeLow,m_rangeHigh); 
00201         }
00202         
00203 }
00204 
00205 void TrackerGeometryCompare::beginJob(){
00206   firstEvent_ = true;
00207 }
00208 
00209 void TrackerGeometryCompare::endJob(){
00210 
00211   int iname(0) ;
00212   for (std::vector<TrackerMap>::iterator it = m_vtkmap.begin(); it != m_vtkmap.end(); ++it) {
00213     std::stringstream mapname ;
00214     mapname << "TkMap_SurfDeform" << iname << ".png" ; 
00215     it->save(true,0,0,mapname.str());
00216     mapname.str( std::string() ); 
00217     mapname.clear() ; 
00218     mapname << "TkMap_SurfDeform" << iname << ".pdf" ; 
00219     it->save(true,0,0,mapname.str());
00220     ++iname ; 
00221   }
00222 
00223   _theFile->cd();
00224   _alignTree->Write();
00225   _theFile->Close();
00226         
00227 }
00228 
00229 void TrackerGeometryCompare::analyze(const edm::Event&, const edm::EventSetup& iSetup){
00230 
00231   if (firstEvent_) {
00232 
00233         //upload the ROOT geometries
00234         createROOTGeometry(iSetup);
00235 
00236         //set common tracker system first
00237         // if setting the tracker common system
00238         if (_setCommonTrackerSystem != "NONE"){
00239                 setCommonTrackerSystem();
00240         }
00241         
00242         //compare the goemetries
00243         compareGeometries(referenceTracker,currentTracker);
00244         compareSurfaceDeformations(_inputTree11, _inputTree12); 
00245         
00246         //write out ntuple
00247         //might be better to do within output module
00248         
00249         if (_writeToDB){
00250                 Alignments* myAlignments = currentTracker->alignments();
00251                 AlignmentErrors* myAlignmentErrors = currentTracker->alignmentErrors();
00252                 
00253                 // 2. Store alignment[Error]s to DB
00254                 edm::Service<cond::service::PoolDBOutputService> poolDbService;
00255                 // Call service
00256                 if( !poolDbService.isAvailable() ) // Die if not available
00257                         throw cms::Exception("NotAvailable") << "PoolDBOutputService not available";
00258                 
00259                 poolDbService->writeOne<Alignments>(&(*myAlignments), poolDbService->beginOfTime(), "TrackerAlignmentRcd");
00260                 poolDbService->writeOne<AlignmentErrors>(&(*myAlignmentErrors), poolDbService->beginOfTime(), "TrackerAlignmentErrorRcd");
00261                 
00262         }               
00263 
00264         firstEvent_ = false;
00265   }
00266 }
00267 
00268 void TrackerGeometryCompare::createROOTGeometry(const edm::EventSetup& iSetup){
00269         
00270         int inputRawId1, inputRawId2;
00271         double inputX1, inputY1, inputZ1, inputX2, inputY2, inputZ2;
00272         double inputAlpha1, inputBeta1, inputGamma1, inputAlpha2, inputBeta2, inputGamma2;
00273                 
00274         //declare alignments
00275         Alignments* alignments1 = new Alignments();
00276         AlignmentErrors* alignmentErrors1 = new AlignmentErrors();      
00277         if (_inputFilename1 != "IDEAL"){
00278                 _inputRootFile1 = new TFile(_inputFilename1.c_str());
00279                 TTree* _inputTree01 = (TTree*) _inputRootFile1->Get(_inputTreenameAlign.c_str());
00280                 _inputTree01->SetBranchAddress("rawid", &inputRawId1);
00281                 _inputTree01->SetBranchAddress("x", &inputX1);
00282                 _inputTree01->SetBranchAddress("y", &inputY1);
00283                 _inputTree01->SetBranchAddress("z", &inputZ1);
00284                 _inputTree01->SetBranchAddress("alpha", &inputAlpha1);
00285                 _inputTree01->SetBranchAddress("beta", &inputBeta1);
00286                 _inputTree01->SetBranchAddress("gamma", &inputGamma1);
00287 
00288                 int nEntries1 = _inputTree01->GetEntries();
00289                 //fill alignments
00290                 for (int i = 0; i < nEntries1; ++i){
00291                         
00292                         _inputTree01->GetEntry(i);
00293                         CLHEP::Hep3Vector translation1(inputX1, inputY1, inputZ1);
00294                         CLHEP::HepEulerAngles eulerangles1(inputAlpha1,inputBeta1,inputGamma1);
00295                         uint32_t detid1 = inputRawId1;
00296                         AlignTransform transform1(translation1, eulerangles1, detid1);
00297                         alignments1->m_align.push_back(transform1);
00298                         
00299                         //dummy errors
00300                         CLHEP::HepSymMatrix clhepSymMatrix(3,0);
00301                         AlignTransformError transformError(clhepSymMatrix, detid1);
00302                         alignmentErrors1->m_alignError.push_back(transformError);
00303                 }               
00304                 
00305                 // to get the right order
00306                 std::sort( alignments1->m_align.begin(), alignments1->m_align.end(), lessAlignmentDetId<AlignTransform>() );
00307                 std::sort( alignmentErrors1->m_alignError.begin(), alignmentErrors1->m_alignError.end(), lessAlignmentDetId<AlignTransformError>() );
00308         }
00309         //------------------
00310         Alignments* alignments2 = new Alignments();
00311         AlignmentErrors* alignmentErrors2 = new AlignmentErrors();
00312         if (_inputFilename2 != "IDEAL"){        
00313                 _inputRootFile2 = new TFile(_inputFilename2.c_str());
00314                 TTree* _inputTree02 = (TTree*) _inputRootFile2->Get(_inputTreenameAlign.c_str());
00315                 _inputTree02->SetBranchAddress("rawid", &inputRawId2);
00316                 _inputTree02->SetBranchAddress("x", &inputX2);
00317                 _inputTree02->SetBranchAddress("y", &inputY2);
00318                 _inputTree02->SetBranchAddress("z", &inputZ2);
00319                 _inputTree02->SetBranchAddress("alpha", &inputAlpha2);
00320                 _inputTree02->SetBranchAddress("beta", &inputBeta2);
00321                 _inputTree02->SetBranchAddress("gamma", &inputGamma2);
00322                 
00323                 int nEntries2 = _inputTree02->GetEntries();
00324                 //fill alignments
00325                 for (int i = 0; i < nEntries2; ++i){
00326                         
00327                         _inputTree02->GetEntry(i);
00328                         CLHEP::Hep3Vector translation2(inputX2, inputY2, inputZ2);
00329                         CLHEP::HepEulerAngles eulerangles2(inputAlpha2,inputBeta2,inputGamma2);
00330                         uint32_t detid2 = inputRawId2;
00331                         AlignTransform transform2(translation2, eulerangles2, detid2);
00332                         alignments2->m_align.push_back(transform2);
00333                         
00334                         //dummy errors
00335                         CLHEP::HepSymMatrix clhepSymMatrix(3,0);
00336                         AlignTransformError transformError(clhepSymMatrix, detid2);
00337                         alignmentErrors2->m_alignError.push_back(transformError); 
00338                 }                       
00339                 
00340                 //to get the right order
00341                 std::sort( alignments2->m_align.begin(), alignments2->m_align.end(), lessAlignmentDetId<AlignTransform>() );
00342                 std::sort( alignmentErrors2->m_alignError.begin(), alignmentErrors2->m_alignError.end(), lessAlignmentDetId<AlignTransformError>() );
00343         }
00344         
00345         //accessing the initial geometry
00346         edm::ESTransientHandle<DDCompactView> cpv;
00347         iSetup.get<IdealGeometryRecord>().get(cpv);
00348         edm::ESHandle<GeometricDet> theGeometricDet;
00349         iSetup.get<IdealGeometryRecord>().get(theGeometricDet);
00350         TrackerGeomBuilderFromGeometricDet trackerBuilder;
00351         
00352         edm::ESHandle<Alignments> globalPositionRcd;
00353         iSetup.get<TrackerDigiGeometryRecord>().getRecord<GlobalPositionRcd>().get(globalPositionRcd);
00354         
00355         //reference tracker
00356         TrackerGeometry* theRefTracker = trackerBuilder.build(&*theGeometricDet); //DM_534?? 
00357         if (_inputFilename1 != "IDEAL"){
00358                 GeometryAligner aligner1;
00359                 aligner1.applyAlignments<TrackerGeometry>( &(*theRefTracker), &(*alignments1), &(*alignmentErrors1),
00360                                                                                                   align::DetectorGlobalPosition(*globalPositionRcd, DetId(DetId::Tracker)));
00361         }
00362         referenceTracker = new AlignableTracker(&(*theRefTracker));
00363         //referenceTracker->setSurfaceDeformation(surfDef1, true) ; 
00364 
00365         int inputRawid1;
00366         int inputRawid2;
00367         int inputDtype1, inputDtype2 ; 
00368         std::vector<double> inputDpar1;
00369         std::vector<double> inputDpar2 ; 
00370         std::vector<double>* p_inputDpar1 = &inputDpar1; 
00371         std::vector<double>* p_inputDpar2 = &inputDpar2; 
00372 
00373         const std::vector<Alignable*> comp1 = referenceTracker->deepComponents(); 
00374 
00375         SurfaceDeformation * surfDef1; 
00376         if (_inputFilename1 != "IDEAL"){
00377           TTree* _inputTree11 = (TTree*) _inputRootFile1->Get(_inputTreenameDeform.c_str());
00378           _inputTree11->SetBranchAddress("irawid", &inputRawid1);
00379           _inputTree11->SetBranchAddress("dtype", &inputDtype1);
00380           _inputTree11->SetBranchAddress("dpar", &p_inputDpar1);
00381 
00382           unsigned int nEntries11 = _inputTree11->GetEntries();
00383           edm::LogInfo("TrackerGeometryCompare") << " nentries11 = " << nEntries11 << std::endl ; 
00384           for (unsigned int iEntry = 0; iEntry < nEntries11; ++iEntry) {
00385             _inputTree11->GetEntry(iEntry) ; 
00386 
00387             surfDef1 = SurfaceDeformationFactory::create( inputDtype1, inputDpar1);
00388 
00389             if (int(comp1[iEntry]->id()) == inputRawid1) {
00390               comp1[iEntry]->setSurfaceDeformation(surfDef1, true) ; 
00391             }
00392 
00393           }
00394         }
00395                 
00396         //currernt tracker
00397         TrackerGeometry* theCurTracker = trackerBuilder.build(&*theGeometricDet); //DM_534??
00398         if (_inputFilename2 != "IDEAL"){
00399                 GeometryAligner aligner2;
00400                 aligner2.applyAlignments<TrackerGeometry>( &(*theCurTracker), &(*alignments2), &(*alignmentErrors2),
00401                                                                                                   align::DetectorGlobalPosition(*globalPositionRcd, DetId(DetId::Tracker)));
00402         }
00403         currentTracker = new AlignableTracker(&(*theCurTracker));
00404         
00405         const std::vector<Alignable*> comp2 = currentTracker->deepComponents(); 
00406 
00407         SurfaceDeformation * surfDef2 ; 
00408         if (_inputFilename2 != "IDEAL"){ 
00409           TTree* _inputTree12 = (TTree*) _inputRootFile2->Get(_inputTreenameDeform.c_str());
00410           _inputTree12->SetBranchAddress("irawid", &inputRawid2);
00411           _inputTree12->SetBranchAddress("dtype", &inputDtype2);
00412           _inputTree12->SetBranchAddress("dpar",  &p_inputDpar2);
00413 
00414           unsigned int nEntries12 = _inputTree12->GetEntries();
00415           edm::LogInfo("TrackerGeometryCompare") << " nentries12 = " << nEntries12 << std::endl ; 
00416           for (unsigned int iEntry = 0; iEntry < nEntries12; ++iEntry) {
00417             _inputTree12->GetEntry(iEntry) ; 
00418             
00419             surfDef2 = SurfaceDeformationFactory::create( inputDtype2, inputDpar2);
00420 
00421             if (int(comp2[iEntry]->id()) == inputRawid2) {
00422               comp2[iEntry]->setSurfaceDeformation(surfDef2, true) ; 
00423             }
00424 
00425           }
00426         }
00427                 
00428         delete alignments1;
00429         delete alignmentErrors1;
00430         delete alignments2;
00431         delete alignmentErrors2;
00432 
00433 }
00434 
00435 void TrackerGeometryCompare::compareSurfaceDeformations(TTree* refTree, TTree* curTree) {
00436         
00437   if (_inputFilename1 != "IDEAL" && _inputFilename2 != "IDEAL") {
00438              
00439     int inputRawid1;
00440     int inputRawid2;
00441     int inputSubdetid1, inputSubdetid2 ; 
00442     int inputDtype1, inputDtype2 ; 
00443     std::vector<double> inputDpar1;
00444     std::vector<double> inputDpar2 ; 
00445     std::vector<double>* p_inputDpar1 = &inputDpar1; 
00446     std::vector<double>* p_inputDpar2 = &inputDpar2; 
00447   
00448     TTree* refTree = (TTree*) _inputRootFile1->Get(_inputTreenameDeform.c_str());
00449     refTree->SetBranchAddress("irawid", &inputRawid1);
00450     refTree->SetBranchAddress("subdetid", &inputSubdetid1);
00451     refTree->SetBranchAddress("dtype", &inputDtype1);
00452     refTree->SetBranchAddress("dpar", &p_inputDpar1);
00453   
00454     TTree* curTree = (TTree*) _inputRootFile2->Get(_inputTreenameDeform.c_str());
00455     curTree->SetBranchAddress("irawid", &inputRawid2);
00456     curTree->SetBranchAddress("subdetid", &inputSubdetid2);
00457     curTree->SetBranchAddress("dtype", &inputDtype2);
00458     curTree->SetBranchAddress("dpar",  &p_inputDpar2);
00459   
00460     unsigned int nEntries11 = refTree->GetEntries();
00461     unsigned int nEntries12 = curTree->GetEntries();
00462 
00463     if (nEntries11 != nEntries12) {
00464       edm::LogError("TrackerGeometryCompare")   << " Surface deformation parameters in two geometries differ!\n" ;
00465       return ; 
00466     }
00467     
00468     for (unsigned int iEntry = 0; iEntry < nEntries12; ++iEntry) {
00469       refTree->GetEntry(iEntry) ;
00470       curTree->GetEntry(iEntry) ;
00471       for (int ii = 0; ii < 13; ++ii) { _surfDeform[ii] = -1.0 ; } 
00472       for (int npar = 0; npar < int(inputDpar2.size()); ++npar ) {
00473             if (inputRawid1 == inputRawid2) {
00474         _surfDeform[npar] = inputDpar2.at(npar) - inputDpar1.at(npar) ; 
00475         std::stringstream histname0 ;
00476         histname0 << "SurfDeform_Par_" << npar ;  
00477         if ( TMath::Abs(_surfDeform[npar]) > (m_rangeHigh - m_rangeLow)/(10.*m_nBins) ) m_h1[histname0.str()]->Fill(_surfDeform[npar]) ; 
00478         if (inputSubdetid1 == 1 && inputSubdetid2 == 1) {
00479           std::stringstream histname1 ;
00480           histname1 << "SurfDeform_PixelBarrel_Par_" << npar ;  
00481           if ( TMath::Abs(_surfDeform[npar]) > (m_rangeHigh - m_rangeLow)/(10.*m_nBins) ) m_h1[histname1.str()]->Fill(_surfDeform[npar]) ; 
00482         }
00483         if (inputSubdetid1 == 2 && inputSubdetid2 == 2) {
00484           std::stringstream histname2 ;
00485           histname2 << "SurfDeform_PixelEndcap_Par_" << npar ;  
00486           if ( TMath::Abs(_surfDeform[npar]) > (m_rangeHigh - m_rangeLow)/(10.*m_nBins) ) m_h1[histname2.str()]->Fill(_surfDeform[npar]) ; 
00487         }
00488         (m_vtkmap.at(npar)).fill_current_val(inputRawid1,_surfDeform[npar]) ; 
00489             }
00490       }
00491     }
00492   
00493   } else if ( _inputFilename1 == "IDEAL" && _inputFilename2 != "IDEAL" ) {
00494              
00495     int inputRawid2;
00496     int inputSubdetid2 ; 
00497     int inputDtype2 ; 
00498     std::vector<double> inputDpar2 ; 
00499     std::vector<double>* p_inputDpar2 = &inputDpar2; 
00500   
00501     TTree* curTree = (TTree*) _inputRootFile2->Get(_inputTreenameDeform.c_str());
00502     curTree->SetBranchAddress("irawid", &inputRawid2);
00503     curTree->SetBranchAddress("subdetid", &inputSubdetid2);
00504     curTree->SetBranchAddress("dtype", &inputDtype2);
00505     curTree->SetBranchAddress("dpar",  &p_inputDpar2);
00506   
00507     unsigned int nEntries12 = curTree->GetEntries();
00508     
00509     for (unsigned int iEntry = 0; iEntry < nEntries12; ++iEntry) {
00510       curTree->GetEntry(iEntry) ;
00511       for (int ii = 0; ii < 12; ++ii) { _surfDeform[ii] = -1.0 ; } 
00512       for (int npar = 0; npar < int(inputDpar2.size()); ++npar ) {
00513         _surfDeform[npar] = inputDpar2.at(npar) ; 
00514         std::stringstream histname0 ;
00515         histname0 << "SurfDeform_Par_" << npar ;  
00516         if ( TMath::Abs(_surfDeform[npar]) > (m_rangeHigh - m_rangeLow)/(10.*m_nBins) ) m_h1[histname0.str()]->Fill(_surfDeform[npar]) ; 
00517         if (inputSubdetid2 == 1) {
00518           std::stringstream histname1 ;
00519           histname1 << "SurfDeform_PixelBarrel_Par_" << npar ;  
00520           if ( TMath::Abs(_surfDeform[npar]) > (m_rangeHigh - m_rangeLow)/(10.*m_nBins) ) m_h1[histname1.str()]->Fill(_surfDeform[npar]) ; 
00521         }
00522         if (inputSubdetid2 == 2) {
00523           std::stringstream histname2 ;
00524           histname2 << "SurfDeform_PixelEndcap_Par_" << npar ;  
00525           if ( TMath::Abs(_surfDeform[npar]) > (m_rangeHigh - m_rangeLow)/(10.*m_nBins) ) m_h1[histname2.str()]->Fill(_surfDeform[npar]) ; 
00526         }
00527         (m_vtkmap.at(npar)).fill_current_val(inputRawid2,_surfDeform[npar]) ; 
00528       }
00529     }
00530   
00531   } else if ( _inputFilename1 != "IDEAL" && _inputFilename2 == "IDEAL" ) {
00532              
00533     int inputRawid1;
00534     int inputSubdetid1 ; 
00535     int inputDtype1 ; 
00536     std::vector<double> inputDpar1;
00537     std::vector<double>* p_inputDpar1 = &inputDpar1; 
00538   
00539     TTree* refTree = (TTree*) _inputRootFile1->Get(_inputTreenameDeform.c_str());
00540     refTree->SetBranchAddress("irawid", &inputRawid1);
00541     refTree->SetBranchAddress("subdetid", &inputSubdetid1);
00542     refTree->SetBranchAddress("dtype", &inputDtype1);
00543     refTree->SetBranchAddress("dpar", &p_inputDpar1);
00544   
00545     unsigned int nEntries11 = refTree->GetEntries();
00546     
00547     for (unsigned int iEntry = 0; iEntry < nEntries11; ++iEntry) {
00548       refTree->GetEntry(iEntry) ;
00549       for (int ii = 0; ii < 12; ++ii) { _surfDeform[ii] = -1.0 ; } 
00550       for (int npar = 0; npar < int(inputDpar1.size()); ++npar ) {
00551         _surfDeform[npar] = - inputDpar1.at(npar) ; 
00552         std::stringstream histname0 ;
00553         histname0 << "SurfDeform_Par_" << npar ;  
00554         if ( TMath::Abs(_surfDeform[npar]) > (m_rangeHigh - m_rangeLow)/(10.*m_nBins) ) m_h1[histname0.str()]->Fill(_surfDeform[npar]) ; 
00555         if (inputSubdetid1 == 1) {
00556           std::stringstream histname1 ;
00557           histname1 << "SurfDeform_PixelBarrel_Par_" << npar ;  
00558           if ( TMath::Abs(_surfDeform[npar]) > (m_rangeHigh - m_rangeLow)/(10.*m_nBins) ) m_h1[histname1.str()]->Fill(_surfDeform[npar]) ; 
00559         }
00560         if (inputSubdetid1 == 2) {
00561           std::stringstream histname2 ;
00562           histname2 << "SurfDeform_PixelEndcap_Par_" << npar ;  
00563           if ( TMath::Abs(_surfDeform[npar]) > (m_rangeHigh - m_rangeLow)/(10.*m_nBins) ) m_h1[histname2.str()]->Fill(_surfDeform[npar]) ; 
00564         }
00565         (m_vtkmap.at(npar)).fill_current_val(inputRawid1,_surfDeform[npar]) ; 
00566       }
00567     }
00568   
00569   } else if ( _inputFilename1 == "IDEAL" && _inputFilename2 == "IDEAL" ) {
00570 
00571           edm::LogInfo("TrackerGeometryCompare") << ">>>> Comparing IDEAL with IDEAL: nothing to do! <<<<\n" ; 
00572           
00573   }
00574 
00575   return ;      
00576 }
00577 
00578 void TrackerGeometryCompare::compareGeometries(Alignable* refAli, Alignable* curAli){
00579 
00580         using namespace align ; 
00581         
00582         const std::vector<Alignable*>& refComp = refAli->components();
00583         const std::vector<Alignable*>& curComp = curAli->components();
00584         
00585         unsigned int nComp = refComp.size();
00586         //only perform for designate levels
00587         bool useLevel = false;
00588         for (unsigned int i = 0; i < m_theLevels.size(); ++i){
00589                 if (refAli->alignableObjectId() == m_theLevels[i]) useLevel = true;
00590         }
00591         
00592         //another added level for difference between det and detunit
00593         //if ((refAli->alignableObjectId()==2)&&(nComp == 1)) useLevel = false;
00594         
00595         //coordinate matching, etc etc
00596         if (useLevel){
00597                 DetId detid(refAli->id());
00598 
00599                 CLHEP::Hep3Vector Rtotal, Wtotal, lRtotal, lWtotal;
00600                 Rtotal.set(0.,0.,0.); 
00601                 Wtotal.set(0.,0.,0.);
00602                 lRtotal.set(0.,0.,0.); 
00603                 lWtotal.set(0.,0.,0.);
00604 
00605                 for (int i = 0; i < 100; i++){
00606                         AlgebraicVector diff = align::diffAlignables(refAli,curAli, _weightBy, _weightById, _weightByIdVector);
00607                         CLHEP::Hep3Vector dR(diff[0],diff[1],diff[2]);
00608                         Rtotal+=dR;
00609                         CLHEP::Hep3Vector dW(diff[3],diff[4],diff[5]);
00610                         CLHEP::HepRotation rot(Wtotal.unit(),Wtotal.mag());
00611                         CLHEP::HepRotation drot(dW.unit(),dW.mag());
00612                         rot*=drot;
00613                         Wtotal.set(rot.axis().x()*rot.delta(), rot.axis().y()*rot.delta(), rot.axis().z()*rot.delta());
00614                         // local coordinates
00615                         lRtotal.set(diff[6],diff[7],diff[8]);
00616                         lWtotal.set(diff[9],diff[10],diff[11]);
00617                         align::moveAlignable(curAli, diff);
00618                         float tolerance = 1e-7;
00619                         AlgebraicVector check = align::diffAlignables(refAli,curAli, _weightBy, _weightById, _weightByIdVector);
00620                         align::GlobalVector checkR(check[0],check[1],check[2]);
00621                         align::GlobalVector checkW(check[3],check[4],check[5]);
00622                         if ((checkR.mag() > tolerance)||(checkW.mag() > tolerance)){
00623                                 edm::LogInfo("TrackerGeometryCompare") << "Tolerance Exceeded!(alObjId: " << refAli->alignableObjectId()
00624                                 << ", rawId: " << refAli->geomDetId().rawId()
00625                                 << ", subdetId: "<< detid.subdetId() << "): " << diff;
00626                                 throw cms::Exception("Tolerance in TrackerGeometryCompare exceeded");
00627                         }
00628                         else{
00629                                 break;
00630                         }
00631                 }
00632 
00633                 AlgebraicVector TRtot(12);
00634                 // global 
00635                 TRtot(1) = Rtotal.x(); TRtot(2) = Rtotal.y(); TRtot(3) = Rtotal.z();
00636                 TRtot(4) = Wtotal.x(); TRtot(5) = Wtotal.y(); TRtot(6) = Wtotal.z();
00637                 // local
00638                 TRtot(7) = lRtotal.x(); TRtot(8) = lRtotal.y(); TRtot(9) = lRtotal.z();
00639                 TRtot(10) = lWtotal.x(); TRtot(11) = lWtotal.y(); TRtot(12) = lWtotal.z();
00640 
00641                 fillTree(refAli, TRtot);
00642         }
00643 
00644         // another added level for difference between det and detunit
00645         for (unsigned int i = 0; i < nComp; ++i) 
00646           compareGeometries(refComp[i],curComp[i]);     
00647 
00648 }
00649 
00650 void TrackerGeometryCompare::setCommonTrackerSystem(){
00651 
00652         edm::LogInfo("TrackerGeometryCompare") << "Setting Common Tracker System....";
00653         
00654         AlignableObjectId dummy; // DM_534??
00655         _commonTrackerLevel = dummy.nameToType(_setCommonTrackerSystem); // DM_534?? 
00656                 
00657         diffCommonTrackerSystem(referenceTracker, currentTracker);
00658         
00659         align::EulerAngles dOmega(3); dOmega[0] = _TrackerCommonR.x() ; dOmega[1] = _TrackerCommonR.y(); dOmega[2] = _TrackerCommonR.z();
00660         align::RotationType rot = align::toMatrix( dOmega );
00661         align::GlobalVector theR = _TrackerCommonT;
00662         
00663         std::cout << "what we get from overlaying the pixels..." << theR << ", " << rot << std::endl;
00664         
00665         //transform to the Tracker System
00666         align::PositionType trackerCM = currentTracker->globalPosition();
00667         align::GlobalVector cmDiff( trackerCM.x()-_TrackerCommonCM.x(), trackerCM.y()-_TrackerCommonCM.y(), trackerCM.z()-_TrackerCommonCM.z() );
00668         
00669         std::cout << "Pixel CM: " << _TrackerCommonCM << ", tracker CM: " << trackerCM << std::endl;
00670         
00671         //adjust translational difference factoring in different rotational CM
00672         //needed because rotateInGlobalFrame is about CM of alignable, not Tracker
00673         align::GlobalVector::BasicVectorType lpvgf = cmDiff.basicVector();
00674         align::GlobalVector moveV( rot.multiplyInverse(lpvgf) - lpvgf);
00675         align::GlobalVector theRprime(theR + moveV);
00676         
00677         AlgebraicVector TrackerCommonTR(6);
00678         TrackerCommonTR(1) = theRprime.x(); TrackerCommonTR(2) = theRprime.y(); TrackerCommonTR(3) = theRprime.z();
00679         TrackerCommonTR(4) = _TrackerCommonR.x(); TrackerCommonTR(5) = _TrackerCommonR.y(); TrackerCommonTR(6) = _TrackerCommonR.z();
00680         
00681         std::cout << "and after the transformation: " << TrackerCommonTR << std::endl;
00682         
00683         align::moveAlignable(currentTracker, TrackerCommonTR );
00684         
00685 }
00686 
00687 void TrackerGeometryCompare::diffCommonTrackerSystem(Alignable *refAli, Alignable *curAli){
00688         
00689         const std::vector<Alignable*>& refComp = refAli->components();
00690         const std::vector<Alignable*>& curComp = curAli->components();
00691         
00692         unsigned int nComp = refComp.size();
00693         //only perform for designate levels
00694         bool useLevel = false;
00695         if (refAli->alignableObjectId() == _commonTrackerLevel) useLevel = true;
00696         
00697         //useLevel = false;
00698         if (useLevel){
00699                 CLHEP::Hep3Vector Rtotal, Wtotal;
00700                 Rtotal.set(0.,0.,0.); Wtotal.set(0.,0.,0.);
00701                 
00702                 AlgebraicVector diff = align::diffAlignables(refAli,curAli, _weightBy, _weightById, _weightByIdVector);
00703                 CLHEP::Hep3Vector dR(diff[0],diff[1],diff[2]);
00704                 Rtotal+=dR;
00705                 CLHEP::Hep3Vector dW(diff[3],diff[4],diff[5]);
00706                 CLHEP::HepRotation rot(Wtotal.unit(),Wtotal.mag());
00707                 CLHEP::HepRotation drot(dW.unit(),dW.mag());
00708                 rot*=drot;
00709                 Wtotal.set(rot.axis().x()*rot.delta(), rot.axis().y()*rot.delta(), rot.axis().z()*rot.delta());
00710                 /*
00711                  //std::cout << "a";
00712                  //if (refAli->alignableObjectId() == 1) std::cout << "DIFF: " << diff << std::endl;
00713                  align::moveAlignable(curAli, diff);
00714                  float tolerance = 1e-7;
00715                  AlgebraicVector check = align::diffAlignables(refAli,curAli, _weightBy, _weightById, _weightByIdVector);
00716                  align::GlobalVector checkR(check[0],check[1],check[2]);
00717                  align::GlobalVector checkW(check[3],check[4],check[5]);
00718                  DetId detid(refAli->id());
00719                  if ((checkR.mag() > tolerance)||(checkW.mag() > tolerance)){
00720                  edm::LogInfo("TrackerGeometryCompare") << "Tolerance Exceeded!(alObjId: " << refAli->alignableObjectId()
00721                  << ", rawId: " << refAli->geomDetId().rawId()
00722                  << ", subdetId: "<< detid.subdetId() << "): " << diff;
00723                  }
00724                  else{
00725                  break;
00726                  }
00727                  }
00728                  */
00729                 
00730                 //_TrackerCommonT.set(Rtotal.x(), Rtotal.y(), Rtotal.z());
00731                 _TrackerCommonT = align::GlobalVector(Rtotal.x(), Rtotal.y(), Rtotal.z());
00732                 _TrackerCommonR = align::GlobalVector(Wtotal.x(), Wtotal.y(), Wtotal.z());
00733                 _TrackerCommonCM = curAli->globalPosition();
00734                 //_TrackerCommonTR(1) = Rtotal.x(); _TrackerCommonTR(2) = Rtotal.y(); _TrackerCommonTR(3) = Rtotal.z();
00735                 //_TrackerCommonTR(4) = Wtotal.x(); _TrackerCommonTR(5) = Wtotal.y(); _TrackerCommonTR(6) = Wtotal.z();
00736                 
00737                 
00738         }
00739         else{
00740                 for (unsigned int i = 0; i < nComp; ++i) diffCommonTrackerSystem(refComp[i],curComp[i]);
00741         }
00742         
00743         
00744 }
00745 
00746 void TrackerGeometryCompare::fillTree(Alignable *refAli, AlgebraicVector diff){
00747         
00748         _id = refAli->id();
00749         _level = refAli->alignableObjectId();
00750         //need if ali has no mother
00751         if (refAli->mother()){
00752                 _mid = refAli->mother()->geomDetId().rawId();
00753                 _mlevel = refAli->mother()->alignableObjectId();
00754         }
00755         else{
00756                 _mid = -1;
00757                 _mlevel = -1;
00758         }
00759         DetId detid(_id);
00760         _sublevel = detid.subdetId();
00761         fillIdentifiers( _sublevel, _id );
00762         _xVal = refAli->globalPosition().x();
00763         _yVal = refAli->globalPosition().y();
00764         _zVal = refAli->globalPosition().z();
00765         align::GlobalVector vec(_xVal,_yVal,_zVal);
00766         _rVal = vec.perp();
00767         _phiVal = vec.phi();
00768         _etaVal = vec.eta();
00769         align::RotationType rot = refAli->globalRotation();
00770         align::EulerAngles eulerAngles = align::toAngles(rot);
00771         _alphaVal = eulerAngles[0];
00772         _betaVal = eulerAngles[1];
00773         _gammaVal = eulerAngles[2];
00774         // global
00775         _dxVal = diff[0];
00776         _dyVal = diff[1];
00777         _dzVal = diff[2];
00778         // local
00779         _duVal = diff[6];
00780         _dvVal = diff[7];
00781         _dwVal = diff[8];
00782         //...TODO...
00783         align::GlobalVector g(_dxVal, _dyVal, _dzVal);
00784         //getting dR and dPhi
00785         align::GlobalVector vRef(_xVal,_yVal,_zVal);
00786         align::GlobalVector vCur(_xVal + _dxVal, _yVal + _dyVal, _zVal + _dzVal);
00787         _drVal = vCur.perp() - vRef.perp();
00788         _dphiVal = vCur.phi() - vRef.phi();
00789         // global
00790         _dalphaVal = diff[3];
00791         _dbetaVal = diff[4];
00792         _dgammaVal = diff[5];
00793         // local
00794         _daVal = diff[9];
00795         _dbVal = diff[10];
00796         _dgVal = diff[11];
00797         
00798         //detIdFlag
00799         if (refAli->alignableObjectId() == align::AlignableDetUnit){
00800                 if (_detIdFlag){
00801                         if ((passIdCut(refAli->id()))||(passIdCut(refAli->mother()->id()))){
00802                                 _useDetId = 1;
00803                         }
00804                         else{
00805                                 _useDetId = 0;
00806                         }
00807                 }
00808         }
00809         // det module dimension
00810         if (refAli->alignableObjectId() == align::AlignableDetUnit){
00811                 if (refAli->mother()->alignableObjectId() != align::AlignableDet) _detDim = 1;
00812                 else if (refAli->mother()->alignableObjectId() == align::AlignableDet) _detDim = 2;
00813         }
00814         else _detDim = 0;
00815         
00816         _surWidth = refAli->surface().width();
00817         _surLength = refAli->surface().length();
00818         align::RotationType rt = refAli->globalRotation();
00819         _surRot[0] = rt.xx(); _surRot[1] = rt.xy(); _surRot[2] = rt.xz();
00820         _surRot[3] = rt.yx(); _surRot[4] = rt.yy(); _surRot[5] = rt.yz();
00821         _surRot[6] = rt.zx(); _surRot[7] = rt.zy(); _surRot[8] = rt.zz();
00822         
00823         //Fill
00824         _alignTree->Fill();
00825         
00826 }
00827 
00828 void TrackerGeometryCompare::surveyToTracker(AlignableTracker* ali, Alignments* alignVals, AlignmentErrors* alignErrors){
00829         
00830         //getting the right alignables for the alignment record
00831         std::vector<Alignable*> detPB = ali->pixelHalfBarrelGeomDets();
00832         std::vector<Alignable*> detPEC = ali->pixelEndcapGeomDets();
00833         std::vector<Alignable*> detTIB = ali->innerBarrelGeomDets();
00834         std::vector<Alignable*> detTID = ali->TIDGeomDets();
00835         std::vector<Alignable*> detTOB = ali->outerBarrelGeomDets();
00836         std::vector<Alignable*> detTEC = ali->endcapGeomDets();
00837         
00838         std::vector<Alignable*> allGeomDets;
00839         std::copy(detPB.begin(), detPB.end(), std::back_inserter(allGeomDets));
00840         std::copy(detPEC.begin(), detPEC.end(), std::back_inserter(allGeomDets));
00841         std::copy(detTIB.begin(), detTIB.end(), std::back_inserter(allGeomDets));
00842         std::copy(detTID.begin(), detTID.end(), std::back_inserter(allGeomDets));
00843         std::copy(detTOB.begin(), detTOB.end(), std::back_inserter(allGeomDets));
00844         std::copy(detTEC.begin(), detTEC.end(), std::back_inserter(allGeomDets));
00845         
00846         std::vector<Alignable*> rcdAlis;
00847         for (std::vector<Alignable*>::iterator i = allGeomDets.begin(); i!= allGeomDets.end(); i++){
00848                 if ((*i)->components().size() == 1){
00849                         rcdAlis.push_back((*i));
00850                 }
00851                 else if ((*i)->components().size() > 1){
00852                         rcdAlis.push_back((*i));
00853                         std::vector<Alignable*> comp = (*i)->components();
00854                         for (std::vector<Alignable*>::iterator j = comp.begin(); j != comp.end(); j++){
00855                                 rcdAlis.push_back((*j));
00856                         }
00857                 }
00858         }
00859         
00860         //turning them into alignments
00861         for(std::vector<Alignable*>::iterator k = rcdAlis.begin(); k != rcdAlis.end(); k++){
00862                 
00863                 const SurveyDet* surveyInfo = (*k)->survey();
00864                 align::PositionType pos(surveyInfo->position());
00865                 align::RotationType rot(surveyInfo->rotation());
00866                 CLHEP::Hep3Vector clhepVector(pos.x(),pos.y(),pos.z());
00867                 CLHEP::HepRotation clhepRotation( CLHEP::HepRep3x3(rot.xx(),rot.xy(),rot.xz(),rot.yx(),rot.yy(),rot.yz(),rot.zx(),rot.zy(),rot.zz()));
00868                 AlignTransform transform(clhepVector, clhepRotation, (*k)->id());
00869                 AlignTransformError transformError(CLHEP::HepSymMatrix(3,1), (*k)->id());
00870                 alignVals->m_align.push_back(transform);
00871                 alignErrors->m_alignError.push_back(transformError);
00872         }
00873         
00874         //to get the right order
00875         std::sort( alignVals->m_align.begin(), alignVals->m_align.end(), lessAlignmentDetId<AlignTransform>() );
00876         std::sort( alignErrors->m_alignError.begin(), alignErrors->m_alignError.end(), lessAlignmentDetId<AlignTransformError>() );
00877         
00878 }
00879 
00880 void TrackerGeometryCompare::addSurveyInfo(Alignable* ali){
00881         
00882         const std::vector<Alignable*>& comp = ali->components();
00883         
00884         unsigned int nComp = comp.size();
00885         
00886         for (unsigned int i = 0; i < nComp; ++i) addSurveyInfo(comp[i]);
00887         
00888         const SurveyError& error = theSurveyErrors->m_surveyErrors[theSurveyIndex];
00889         
00890         if ( ali->geomDetId().rawId() != error.rawId() ||
00891                 ali->alignableObjectId() != error.structureType() )
00892         {
00893                 throw cms::Exception("DatabaseError")
00894                 << "Error reading survey info from DB. Mismatched id!";
00895         }
00896         
00897         const CLHEP::Hep3Vector&  pos = theSurveyValues->m_align[theSurveyIndex].translation();
00898         const CLHEP::HepRotation& rot = theSurveyValues->m_align[theSurveyIndex].rotation();
00899         
00900         AlignableSurface surf( align::PositionType( pos.x(), pos.y(), pos.z() ),
00901                                                   align::RotationType( rot.xx(), rot.xy(), rot.xz(),
00902                                                                                           rot.yx(), rot.yy(), rot.yz(),
00903                                                                                           rot.zx(), rot.zy(), rot.zz() ) );
00904         
00905         surf.setWidth( ali->surface().width() );
00906         surf.setLength( ali->surface().length() );
00907         
00908         ali->setSurvey( new SurveyDet( surf, error.matrix() ) );
00909         
00910         ++theSurveyIndex;
00911         
00912 }
00913 
00914 bool TrackerGeometryCompare::passIdCut( uint32_t id ){
00915         
00916         bool pass = false;
00917         int nEntries = _detIdFlagVector.size();
00918         
00919         for (int i = 0; i < nEntries; i++){
00920                 if (_detIdFlagVector[i] == id) pass = true;
00921         }
00922         
00923         return pass;
00924         
00925 }
00926 
00927 void TrackerGeometryCompare::fillIdentifiers( int subdetlevel, int rawid ){
00928         
00929         
00930         switch( subdetlevel ){
00931                         
00932                 case 1:
00933                 {
00934                         PXBDetId pxbid( rawid );
00935                         _identifiers[0] = pxbid.module();
00936                         _identifiers[1] = pxbid.ladder();
00937                         _identifiers[2] = pxbid.layer();
00938                         _identifiers[3] = 999;
00939                         _identifiers[4] = 999;
00940                         _identifiers[5] = 999;
00941                         break;
00942                 }
00943                 case 2:
00944                 {
00945                         PXFDetId pxfid( rawid );
00946                         _identifiers[0] = pxfid.module();
00947                         _identifiers[1] = pxfid.panel();
00948                         _identifiers[2] = pxfid.blade();
00949                         _identifiers[3] = pxfid.disk();
00950                         _identifiers[4] = pxfid.side();
00951                         _identifiers[5] = 999;
00952                         break;
00953                 }
00954                 case 3:
00955                 {
00956                         TIBDetId tibid( rawid );
00957                         _identifiers[0] = tibid.module();
00958                         _identifiers[1] = tibid.string()[0];
00959                         _identifiers[2] = tibid.string()[1];
00960                         _identifiers[3] = tibid.string()[2];
00961                         _identifiers[4] = tibid.layer();
00962                         _identifiers[5] = 999;
00963                         break;
00964                 }
00965                 case 4: 
00966                 {
00967                         TIDDetId tidid( rawid );
00968                         _identifiers[0] = tidid.module()[0];
00969                         _identifiers[1] = tidid.module()[1];
00970                         _identifiers[2] = tidid.ring();
00971                         _identifiers[3] = tidid.wheel();
00972                         _identifiers[4] = tidid.side();
00973                         _identifiers[5] = 999;
00974                         break;
00975                 }
00976                 case 5: 
00977                 {
00978                         TOBDetId tobid( rawid );
00979                         _identifiers[0] = tobid.module();
00980                         _identifiers[1] = tobid.rod()[0];
00981                         _identifiers[2] = tobid.rod()[1];
00982                         _identifiers[3] = tobid.layer();
00983                         _identifiers[4] = 999;
00984                         _identifiers[5] = 999;
00985                         break;
00986                 }
00987                 case 6: 
00988                 {
00989                         TECDetId tecid( rawid );
00990                         _identifiers[0] = tecid.module();
00991                         _identifiers[1] = tecid.ring();
00992                         _identifiers[2] = tecid.petal()[0];
00993                         _identifiers[3] = tecid.petal()[1];
00994                         _identifiers[4] = tecid.wheel();
00995                         _identifiers[5] = tecid.side();
00996                         break;
00997                 }
00998                 default:
00999                 {
01000                         std::cout << "Error: bad subdetid!!" << std::endl;
01001                         break;
01002                 }
01003                         
01004         }
01005 }
01006 
01007 
01008 DEFINE_FWK_MODULE(TrackerGeometryCompare);