Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-06-06 07:38:46

0001 // Small Program to read Grid Files
0002 // by droll (29/02/04)
0003 // essential files
0004 #include "MagneticField/Interpolation/interface/binary_ifstream.h"
0005 
0006 // used libs
0007 #include <string>
0008 #include <iostream>
0009 #include <fstream>
0010 #include <iomanip>
0011 #include <cmath>
0012 #include <stdlib.h>
0013 #include <optional>
0014 #include "DataFormats/Math/interface/approx_exp.h"
0015 inline int bits(int a) {
0016   unsigned int aa = abs(a);
0017   int b = 0;
0018   if (a == 0)
0019     return 0;
0020   while ((aa /= 2) > 0)
0021     ++b;
0022   return (a > 0) ? b : -b;
0023 }
0024 
0025 using namespace std;
0026 using namespace approx_math;
0027 int main(int argc, char **argv) {
0028   if (argc < 3 or argc > 4) {
0029     cout << "SYNOPSIS:" << endl << " GridFileReader input.bin [fullDump=true|false] [file.index]" << endl;
0030     cout << "Example:" << endl << " GridFileReader grid.217.bin false" << endl;
0031     return 1;
0032   }
0033   const string filename = argv[1];
0034 
0035   bool fullDump = argv[2];
0036 
0037   string indexFile;
0038   if (argc == 4) {
0039     indexFile = argv[3];
0040   }
0041 
0042   std::optional<magneticfield::interpolation::binary_ifstream> tempFile;
0043   if (indexFile.empty()) {
0044     tempFile = magneticfield::interpolation::binary_ifstream(filename);
0045   } else {
0046     auto pos = indexFile.find('.');
0047     tempFile = magneticfield::interpolation::binary_ifstream(indexFile.substr(0, pos) + ".bin");
0048   }
0049   auto inFile = std::move(*tempFile);
0050   if (!inFile) {
0051     cout << "file open failed!" << endl;
0052     return EXIT_FAILURE;
0053   }
0054 
0055   if (not indexFile.empty()) {
0056     std::ifstream index(indexFile);
0057 
0058     bool found = true;
0059     while (index) {
0060       std::string magFile;
0061       unsigned int offset;
0062       index >> magFile >> offset;
0063 
0064       if (magFile == filename) {
0065         inFile.seekg(offset);
0066         break;
0067       }
0068     }
0069     if (not found) {
0070       cout << "failed to find entry " << filename << " in index file " << indexFile << endl;
0071       return EXIT_FAILURE;
0072     }
0073     cout << "Using index file: " << indexFile << endl;
0074   }
0075 
0076   cout << "Data File: " << filename << endl;
0077 
0078   // reading iterator (number of entries)
0079   // declaration of all variable needed for file reading
0080   int type;
0081   int nPnt[3];
0082   double ooPnt[3];
0083   double dist0[3];
0084   double dist1[3][3];
0085   double dist2[3][3];
0086   double rParm[4];
0087   bool easyC[3];
0088   // float Bx, By, Bz;
0089   // reading the type
0090   inFile >> type;
0091   // reading the header
0092   if (type == 1) {
0093     inFile >> nPnt[0] >> nPnt[1] >> nPnt[2];
0094     inFile >> ooPnt[0] >> ooPnt[1] >> ooPnt[2];
0095     inFile >> dist0[0] >> dist0[1] >> dist0[2];
0096   }
0097   if (type == 2) {
0098     inFile >> nPnt[0] >> nPnt[1] >> nPnt[2];
0099     inFile >> ooPnt[0] >> ooPnt[1] >> ooPnt[2];
0100     inFile >> dist0[0] >> dist0[1] >> dist0[2];
0101     inFile >> dist1[0][0] >> dist1[1][0] >> dist1[2][0];
0102     inFile >> dist1[0][1] >> dist1[1][1] >> dist1[2][1];
0103     inFile >> dist1[0][2] >> dist1[1][2] >> dist1[2][2];
0104     inFile >> dist2[0][0] >> dist2[1][0] >> dist2[2][0];
0105     inFile >> dist2[0][1] >> dist2[1][1] >> dist2[2][1];
0106     inFile >> dist2[0][2] >> dist2[1][2] >> dist2[2][2];
0107     inFile >> easyC[0] >> easyC[1] >> easyC[2];
0108   }
0109   if (type == 3) {
0110     inFile >> nPnt[0] >> nPnt[1] >> nPnt[2];
0111     inFile >> ooPnt[0] >> ooPnt[1] >> ooPnt[2];
0112     inFile >> dist0[0] >> dist0[1] >> dist0[2];
0113   }
0114   if (type == 4) {
0115     inFile >> nPnt[0] >> nPnt[1] >> nPnt[2];
0116     inFile >> ooPnt[0] >> ooPnt[1] >> ooPnt[2];
0117     inFile >> dist0[0] >> dist0[1] >> dist0[2];
0118     inFile >> dist1[0][0] >> dist1[1][0] >> dist1[2][0];
0119     inFile >> dist1[0][1] >> dist1[1][1] >> dist1[2][1];
0120     inFile >> dist1[0][2] >> dist1[1][2] >> dist1[2][2];
0121     inFile >> dist2[0][0] >> dist2[1][0] >> dist2[2][0];
0122     inFile >> dist2[0][1] >> dist2[1][1] >> dist2[2][1];
0123     inFile >> dist2[0][2] >> dist2[1][2] >> dist2[2][2];
0124     inFile >> easyC[0] >> easyC[1] >> easyC[2];
0125   }
0126   if (type == 5) {
0127     inFile >> nPnt[0] >> nPnt[1] >> nPnt[2];
0128     inFile >> ooPnt[0] >> ooPnt[1] >> ooPnt[2];
0129     inFile >> dist0[0] >> dist0[1] >> dist0[2];
0130     inFile >> rParm[0] >> rParm[1] >> rParm[2] >> rParm[3];
0131   }
0132 
0133   //reading the field
0134   int nLines = nPnt[0] * nPnt[1] * nPnt[2];
0135 
0136   // print the stuff from above (onlt last line of field values)
0137   cout << "  content of " << filename << endl;
0138   cout << type << endl;
0139   if (type == 1) {
0140     cout << nPnt[0] << " " << nPnt[1] << " " << nPnt[2] << endl;
0141     cout << ooPnt[0] << " " << ooPnt[1] << " " << ooPnt[2] << endl;
0142     cout << dist0[0] << " " << dist0[1] << " " << dist0[2] << endl;
0143   }
0144   if (type == 2) {
0145     cout << nPnt[0] << " " << nPnt[1] << " " << nPnt[2] << endl;
0146     cout << ooPnt[0] << " " << ooPnt[1] << " " << ooPnt[2] << endl;
0147     cout << dist0[0] << " " << dist0[1] << " " << dist0[2] << endl;
0148     cout << dist1[0][0] << " " << dist1[1][0] << " " << dist1[2][0] << endl;
0149     cout << dist1[0][1] << " " << dist1[1][1] << " " << dist1[2][1] << endl;
0150     cout << dist1[0][2] << " " << dist1[1][2] << " " << dist1[2][2] << endl;
0151     cout << dist2[0][0] << " " << dist2[1][0] << " " << dist2[2][0] << endl;
0152     cout << dist2[0][1] << " " << dist2[1][1] << " " << dist2[2][1] << endl;
0153     cout << dist2[0][2] << " " << dist2[1][2] << " " << dist2[2][2] << endl;
0154     cout << easyC[0] << " " << easyC[1] << " " << easyC[2] << endl;
0155   }
0156   if (type == 3) {
0157     cout << nPnt[0] << " " << nPnt[1] << " " << nPnt[2] << endl;
0158     cout << ooPnt[0] << " " << ooPnt[1] << " " << ooPnt[2] << endl;
0159     cout << dist0[0] << " " << dist0[1] << " " << dist0[2] << endl;
0160   }
0161   if (type == 4) {
0162     cout << nPnt[0] << " " << nPnt[1] << " " << nPnt[2] << endl;
0163     cout << ooPnt[0] << " " << ooPnt[1] << " " << ooPnt[2] << endl;
0164     cout << dist0[0] << " " << dist0[1] << " " << dist0[2] << endl;
0165     cout << dist1[0][0] << " " << dist1[1][0] << " " << dist1[2][0] << endl;
0166     cout << dist1[0][1] << " " << dist1[1][1] << " " << dist1[2][1] << endl;
0167     cout << dist1[0][2] << " " << dist1[1][2] << " " << dist1[2][2] << endl;
0168     cout << dist2[0][0] << " " << dist2[1][0] << " " << dist2[2][0] << endl;
0169     cout << dist2[0][1] << " " << dist2[1][1] << " " << dist2[2][1] << endl;
0170     cout << dist2[0][2] << " " << dist2[1][2] << " " << dist2[2][2] << endl;
0171     cout << easyC[0] << " " << easyC[1] << " " << easyC[2] << endl;
0172   }
0173   if (type == 5) {
0174     cout << nPnt[0] << " " << nPnt[1] << " " << nPnt[2] << endl;
0175     cout << ooPnt[0] << " " << ooPnt[1] << " " << ooPnt[2] << endl;
0176     cout << dist0[0] << " " << dist0[1] << " " << dist0[2] << endl;
0177     cout << rParm[0] << " " << rParm[1] << " " << rParm[2] << " " << rParm[3] << endl;
0178   }
0179 
0180   float B[3] = {0, 0, 0}, Bmin[3] = {9999., 9999., 9999.}, Bmax[3] = {0, 0, 0.};
0181   for (int iLine = 0; iLine < nLines; ++iLine) {
0182     inFile >> B[0] >> B[1] >> B[2];
0183     for (int i = 0; i != 3; ++i) {
0184       Bmin[i] = std::min(Bmin[i], std::abs(B[i]));
0185       Bmax[i] = std::max(Bmax[i], std::abs(B[i]));
0186     }
0187     if (fullDump) {
0188       cout << setprecision(12);
0189       cout << "line: " << iLine << " " << B[0] << " " << B[1] << " " << B[2] << endl;
0190     }
0191   }
0192 
0193   if (!fullDump) {
0194     cout << ". . ." << endl;
0195     cout << B[0] << " " << B[1] << " " << B[2] << " (last line of B-field only)" << endl;
0196     cout << Bmin[0] << " " << Bmin[1] << " " << Bmin[2] << " (min B-field abs)" << endl;
0197     cout << Bmax[0] << " " << Bmax[1] << " " << Bmax[2] << " (max B-field abs)" << endl;
0198     for (int i = 0; i != 3; ++i)
0199       std::cout << bits(binary32(Bmax[i]).i32 - binary32(Bmin[i]).i32) << " ";
0200     std::cout << "(max-min in bits)" << std::endl;
0201   }
0202 
0203   // check completeness and close file
0204   string lastEntry;
0205   inFile >> lastEntry;
0206   inFile.close();
0207 
0208   cout << "  file is " << lastEntry;
0209   if (lastEntry == "complete")
0210     cout << "  -->  reading done" << endl;
0211   else
0212     cout << "  -->  reading ERROR" << endl;
0213   cout << endl;
0214 
0215   return (0);
0216 }