APEEstimation/macros/writeAPEsInASCII.C

0001 ////////////////////////////////////////////////////////////////////////////////////////
0002 ///
0003 ///  Write APEs in ASCII-file format
0004 ///
0005 ///  The ASCII file contains one row per module, where the first column
0006 ///  lists the module id and the following 21 columns the diagonal and
0007 ///  lower elements x11,x21,x22,x31,x32,x33,0... of the 6x6 covariance
0008 ///  matrix, where the upper 3x3 sub-matrix contains the position APEs.
0009 ///  The elements are stored in units of cm^2.
0010 ///
0011 ///
0012 ///  Before first usage, create 'TrackerTree' that maps DetIDs with tracker
0013 ///  structures by doing:
0014 ///    cd $CMSSW_BASE/src/Alignment/TrackerAlignment
0015 ///    mkdir hists
0016 ///    cd test
0017 ///    cmsRun trackerTreeGenerator_cfg.py
0018 ///
0019 ///
0020 ///  This format is understood as input by
0021 ///  Alignment/CommonAlignmentAlgorithm/python/ApeSettingAlgorithm_cfi.py
0022 ///
0023 ////////////////////////////////////////////////////////////////////////////////////////
0024
0025 #include <exception>
0026 #include <fstream>
0027 #include <iostream>
0028 #include <vector>
0029
0030 #include "TFile.h"
0031 #include "TString.h"
0032 #include "TTree.h"
0033
0034
0035 // APEs in mum, local [x,y,z]
0036 std::vector<double> apesBPXLayer1(3,0.);
0037 std::vector<double> apesBPXLayer2(3,0.);
0038 std::vector<double> apesBPXLayer3(3,0.);
0039 std::vector<double> apesBPXLayer4(3,0.);
0040 std::vector<double> apesFPX(3,0.);
0041 std::vector<double> apesTIB(3,0.);
0042 std::vector<double> apesTOB(3,0.);
0043 std::vector<double> apesTID(3,0.);
0044 std::vector<double> apesTEC(3,0.);
0045
0046 const TString trackerTreeFileName = "../../TrackerAlignment/hists/TrackerTree.root";
0047
0048
0049 // Transform APE into covariance elements in cm^{2} units
0050 std::vector<double> transformIntoCovarianceMatrixElements(const std::vector<double>& apes) {
0051   std::vector<double> cov(21,0.);
0052   cov[0] = apes[0]*apes[0]*1E-8;
0053   cov[2] = apes[1]*apes[1]*1E-8;
0054   cov[5] = apes[2]*apes[2]*1E-8;
0055
0056   return cov;
0057 }
0058
0059
0060 void scaleAPEs(std::vector<double>& apes, const double scale) {
0061   for(std::vector<double>::iterator it = apes.begin();
0062       it != apes.end(); ++it) {
0063     *it = (*it)*scale;
0064   }
0065 }
0066
0067
0068 void writeAPEsInASCII(const TString& outName="ape.txt") {
0069   // set APEs (in mum) for different subdetectors
0070   apesBPXLayer1[0] = 500;
0071   apesBPXLayer1[1] = 500;
0072   apesBPXLayer1[2] = 500;
0073   apesBPXLayer2[0] = 10;
0074   apesBPXLayer2[1] = 40;
0075   apesBPXLayer2[2] = 10;
0076   apesBPXLayer3[0] = 10;
0077   apesBPXLayer3[1] = 10;
0078   apesBPXLayer3[2] = 10;
0079   apesBPXLayer4[0] = 10;
0080   apesBPXLayer4[1] = 10;
0081   apesBPXLayer4[2] = 10;
0082
0083   apesFPX[0] = 10;
0084   apesFPX[1] = 10;
0085   apesFPX[2] = 10;
0086
0087   apesTIB[0] = 10;
0088   apesTIB[1] = 10;
0089   apesTIB[2] = 10;
0090
0091   apesTOB[0] = 10;
0092   apesTOB[1] = 10;
0093   apesTOB[2] = 10;
0094
0095   apesTID[0] = 20;
0096   apesTID[1] = 20;
0097   apesTID[2] = 20;
0098
0099   apesTEC[0] = 20;
0100   apesTEC[1] = 20;
0101   apesTEC[2] = 20;
0102
0103   // scale APEs by
0104   const double scale = 1.;
0105   std::cout << "Scaling APEs by " << scale << std::endl;
0106   scaleAPEs(apesBPXLayer1,scale);
0107   scaleAPEs(apesBPXLayer2,scale);
0108   scaleAPEs(apesBPXLayer3,scale);
0109   scaleAPEs(apesBPXLayer4,scale);
0110   scaleAPEs(apesFPX,scale);
0111   scaleAPEs(apesTIB,scale);
0112   scaleAPEs(apesTOB,scale);
0113   scaleAPEs(apesTID,scale);
0114   scaleAPEs(apesTEC,scale);
0115
0116   // transform into covariance elements
0117   const std::vector<double> covBPXLayer1 = transformIntoCovarianceMatrixElements(apesBPXLayer1);
0118   const std::vector<double> covBPXLayer2 = transformIntoCovarianceMatrixElements(apesBPXLayer2);
0119   const std::vector<double> covBPXLayer3 = transformIntoCovarianceMatrixElements(apesBPXLayer3);
0120   const std::vector<double> covBPXLayer4 = transformIntoCovarianceMatrixElements(apesBPXLayer4);
0121   const std::vector<double> covFPX = transformIntoCovarianceMatrixElements(apesFPX);
0122   const std::vector<double> covTIB = transformIntoCovarianceMatrixElements(apesTIB);
0123   const std::vector<double> covTOB = transformIntoCovarianceMatrixElements(apesTOB);
0124   const std::vector<double> covTID = transformIntoCovarianceMatrixElements(apesTID);
0125   const std::vector<double> covTEC = transformIntoCovarianceMatrixElements(apesTEC);
0126
0127   // open file with tracker-geometry info
0128   TFile file(trackerTreeFileName,"READ");
0129   if( !file.IsOpen() ) {
0130     std::cerr << "\n\nERROR opening file '" << trackerTreeFileName << "'\n" << std::endl;
0131     throw std::exception();
0132   }
0133
0134   // get tree with geometry info
0135   TTree* tree = 0;
0136   const TString treeName = "TrackerTreeGenerator/TrackerTree/TrackerTree";
0137   file.GetObject(treeName,tree);
0138   if( tree == 0 ) {
0139     std::cerr << "\n\nERROR reading tree '" << treeName << "' from file '" << trackerTreeFileName << "'\n" << std::endl;
0140     throw std::exception();
0141   }
0142
0143   // tree variables
0144   unsigned int theRawId = 0;
0145   unsigned int theSubdetId = 0;
0146   unsigned int theLayerId = 0;
0147   tree->SetBranchAddress("RawId",&theRawId);
0148   tree->SetBranchAddress("SubdetId",&theSubdetId);
0149   tree->SetBranchAddress("Layer",&theLayerId);
0150
0151   // open the output file
0152   std::ofstream apeSaveFile(outName.Data());
0153
0154   for(int iE = 0; iE < tree->GetEntries(); ++iE) {
0155     tree->GetEntry(iE);
0156
0157     // Set the APE according to the subdetector.
0158     // The subdetector encoding in tree
0159     // BPIX: 1
0160     // FPIX: 2
0161     // TIB:  3
0162     // TID:  4
0163     // TOB:  5
0164     // TEC:  6
0165     const std::vector<double>* cov = 0;
0166     if(      theSubdetId == 1 ) {
0167       if(      theLayerId == 1 ) cov = &covBPXLayer1;
0168       else if( theLayerId == 2 ) cov = &covBPXLayer2;
0169       else if( theLayerId == 3 ) cov = &covBPXLayer3;
0170       else                       cov = &covBPXLayer4;
0171     }
0172     else if( theSubdetId == 2 ) cov = &covFPX;
0173     else if( theSubdetId == 3 ) cov = &covTIB;
0174     else if( theSubdetId == 4 ) cov = &covTID;
0175     else if( theSubdetId == 5 ) cov = &covTOB;
0176     else if( theSubdetId == 6 ) cov = &covTEC;
0177
0178     // write APE to ASCII file
0179     apeSaveFile << theRawId;
0180     for(std::vector<double>::const_iterator it = cov->begin();
0181     it != cov->end(); ++it) {
0182       apeSaveFile << "  " << *it;
0183     }
0184     apeSaveFile << std::endl;
0185
0186   } // end of loop over tree (=modules)
0187
0188   apeSaveFile.close();
0189   delete tree;
0190   file.Close();
0191
0192   std::cout << "Wrote APEs to '" << outName << "'" << std::endl;
0193 }
0194
0195
0196