Engine/test/testMagneticField.cc

0001 /** \file
0002  *
0003  *  A set of tests for regression and validation of the field map.
0004  *
0005  *  outputTable: generate txt file with values to be used for regression. Points are generated
0006  *  according to innerRadius, outerRadius, minZ, maxZ
0007  *
0008  *  inputTable: file with input values to be checked against, format depends on inputFileType:
0009  *  xyz = cartesian coordinates in cm (default)
0010  *  rpz_m = r, phi, Z in m
0011  *  xyz_m = cartesian in m
0012  *  TOSCA = input test tables, searches for the corresponding volume/sector determined from the file name and path.
0013  *  TOSCAFileList = file with a list of TOSCA tables
0014  *  TOSCASecorComparison: compare each if the listed TOSCA txt tables with those of the other sectors
0015  *
0016  *  \author N. Amapane - CERN
0017  */
0018
0019 #include <memory>
0020
0021 #include "FWCore/Framework/interface/Frameworkfwd.h"
0022 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0023
0024 #include "FWCore/Framework/interface/Event.h"
0025 #include "FWCore/Framework/interface/MakerMacros.h"
0026
0027 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0028
0029 #include "MagneticField/Engine/interface/MagneticField.h"
0030 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0031
0032 #include "FWCore/Framework/interface/ESHandle.h"
0033 #include "FWCore/Framework/interface/EventSetup.h"
0034
0035 #include "DataFormats/GeometryVector/interface/Pi.h"
0036 #include "DataFormats/GeometryVector/interface/CoordinateSets.h"
0037 #include "MagneticField/GeomBuilder/test/stubs/GlobalPointProvider.h"
0038 #include "MagneticField/VolumeBasedEngine/interface/MagGeometry.h"
0039 #include "MagneticField/VolumeGeometry/interface/MagVolume6Faces.h"
0040
0041 #include "FWCore/ParameterSet/interface/FileInPath.h"
0042 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0043
0044 #include <iostream>
0045 #include <string>
0046 #include <sstream>
0047 #include <fstream>
0048 #include <iomanip>
0049 #include <libgen.h>
0050
0051 using namespace edm;
0052 using namespace Geom;
0053 using namespace std;
0054
0055 class testMagneticField : public edm::one::EDAnalyzer<> {
0056 public:
0057   testMagneticField(const edm::ParameterSet& pset) : fieldToken(esConsumes()) {
0058     //    verbose::debugOut = true;
0059     outputFile = pset.getUntrackedParameter<string>("outputTable", "");
0060     inputFile = pset.getUntrackedParameter<string>("inputTable", "");
0061     inputFileType = pset.getUntrackedParameter<string>("inputTableType", "xyz");
0062
0063     //    resolution for validation of maps
0064     reso = pset.getUntrackedParameter<double>("resolution", 0.0001);
0065     //    number of random points to try
0066     numberOfPoints = pset.getUntrackedParameter<int>("numberOfPoints", 10000);
0067     //    outer radius of test cylinder
0068     innerRadius = pset.getUntrackedParameter<double>("InnerRadius", 0.);
0069     outerRadius = pset.getUntrackedParameter<double>("OuterRadius", 900);
0070     //    Z extent of test cylinder
0071     minZ = pset.getUntrackedParameter<double>("minZ", -2400);
0072     maxZ = pset.getUntrackedParameter<double>("maxZ", 2400);
0073   }
0074
0075   ~testMagneticField() override {}
0076
0077   void go(GlobalPoint g) { std::cout << "At: " << g << " phi=" << g.phi() << " B= " << field->inTesla(g) << std::endl; }
0078
0079   void analyze(const edm::Event& event, const edm::EventSetup& setup) final {
0080     field = &setup.getData(fieldToken);
0081
0082     std::cout << "Nominal Field " << field->nominalValue() << "\n" << std::endl;
0083
0084     go(GlobalPoint(0, 0, 0));
0085
0086     if (outputFile != "") {
0087       writeValidationTable(numberOfPoints, outputFile);
0088       return;
0089     }
0090
0091     if (inputFileType == "TOSCA") {
0092       validateVsTOSCATable(inputFile);
0093     } else if (inputFileType == "TOSCAFileList") {
0094       ifstream file(inputFile.c_str());
0095       string table;
0096       while (getline(file, table)) {
0097         validateVsTOSCATable(table);
0098       }
0099     } else if (inputFileType == "TOSCASectorComparison") {
0100       ifstream file(inputFile.c_str());
0101       string table;
0102
0103       cout << "Vol.      1       2       3       4       5       6       7       8       9      10      11      12"
0104            << endl;
0105
0106       while (getline(file, table)) {
0107         compareSectorTables(table);
0108       }
0109     } else if (inputFile != "") {
0110       validate(inputFile, inputFileType);
0111     }
0112
0113     // Some ad-hoc test
0114     //    for (float phi = 0; phi<Geom::twoPi(); phi+=Geom::pi()/48.) {
0115     //      go(GlobalPoint(Cylindrical2Cartesian<float>(89.,phi,145.892)), magfield.product());
0116     //   }
0117   }
0118
0119   void writeValidationTable(int npoints, string filename);
0120   void validate(string filename, string type = "xyz");
0121   void validateVsTOSCATable(string filename);
0122
0123   const MagVolume6Faces* findVolume(GlobalPoint& gp);
0124   const MagVolume6Faces* findMasterVolume(int volume, int sector);
0125
0126   void parseTOSCATablePath(string filename, int& volNo, int& sector, string& type);
0127   void fillFromTable(string inputFile, vector<GlobalPoint>& p, vector<GlobalVector>& b, string type);
0128   void compareSectorTables(string file);
0129
0130 private:
0131   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> fieldToken;
0132   const MagneticField* field;
0133   string inputFile;
0134   string inputFileType;
0135   string outputFile;
0136   double reso;
0137   int numberOfPoints;
0138   double outerRadius;
0139   double innerRadius;
0140   double minZ;
0141   double maxZ;
0142 };
0143
0144 void testMagneticField::writeValidationTable(int npoints, string filename) {
0145   GlobalPointProvider p(innerRadius, outerRadius, -Geom::pi(), Geom::pi(), minZ, maxZ);
0146
0147   std::string::size_type ps = filename.rfind(".");
0148   if (ps != std::string::npos && filename.substr(ps) == ".txt") {
0149     ofstream file(filename.c_str());
0150
0151     for (int i = 0; i < npoints; ++i) {
0152       GlobalPoint gp = p.getPoint();
0153       GlobalVector f = field->inTesla(gp);
0154       file << setprecision(9)  //<< i << " "
0155            << gp.x() << " " << gp.y() << " " << gp.z() << " " << f.x() << " " << f.y() << " " << f.z() << endl;
0156     }
0157   } else {
0158     ofstream file(filename.c_str(), ios::binary);
0159     float px, py, pz, bx, by, bz;
0160
0161     for (int i = 0; i < npoints; ++i) {
0162       GlobalPoint gp = p.getPoint();
0163       GlobalVector f = field->inTesla(gp);
0164       px = gp.x();
0165       py = gp.y();
0166       pz = gp.z();
0167       bx = f.x();
0168       by = f.y();
0169       bz = f.z();
0170       file.write((char*)&px, sizeof(float));
0171       file.write((char*)&py, sizeof(float));
0172       file.write((char*)&pz, sizeof(float));
0173       file.write((char*)&bx, sizeof(float));
0174       file.write((char*)&by, sizeof(float));
0175       file.write((char*)&bz, sizeof(float));
0176     }
0177   }
0178 }
0179
0180 void testMagneticField::validate(string filename, string type) {
0181   ifstream file;
0182
0183   bool binary = true;
0184
0185   string fullPath;
0186   edm::FileInPath mydata(filename);
0187   fullPath = mydata.fullPath();
0188
0189   std::string::size_type ps = filename.rfind(".");
0190   if (ps != std::string::npos && filename.substr(filename.rfind(".")) == ".txt") {
0191     binary = false;
0192     file.open(fullPath.c_str());
0193   } else {
0194     file.open(fullPath.c_str(), ios::binary);
0195   }
0196
0197   string line;
0198
0199   int fail = 0;
0200   int count = 0;
0201
0202   float maxdelta = 0.;
0203
0204   float px, py, pz;
0205   float bx, by, bz;
0206   GlobalPoint gp;
0207
0208   do {
0209     if (binary) {
0210       if (!(file.read((char*)&px, sizeof(float)) && file.read((char*)&py, sizeof(float)) &&
0211             file.read((char*)&pz, sizeof(float)) && file.read((char*)&bx, sizeof(float)) &&
0212             file.read((char*)&by, sizeof(float)) && file.read((char*)&bz, sizeof(float))))
0213         break;
0214       gp = GlobalPoint(px, py, pz);
0215     } else {
0216       if (!(getline(file, line)))
0217         break;
0218       if (line == "" || line[0] == '#')
0219         continue;
0220       stringstream linestr;
0221       linestr << line;
0222       linestr >> px >> py >> pz >> bx >> by >> bz;
0223       if (type == "rpz_m") {  // assume rpz file with units in m.
0224         gp = GlobalPoint(GlobalPoint::Cylindrical(px * 100., py, pz * 100.));
0225       } else if (type == "xyz_m") {  // assume xyz file with units in m.
0226         gp = GlobalPoint(px * 100., py * 100., pz * 100.);
0227       } else {  // assume x,y,z with units in cm
0228         gp = GlobalPoint(px, py, pz);
0229       }
0230     }
0231
0232     if (gp.perp() < innerRadius || gp.perp() > outerRadius || gp.z() < minZ || gp.z() > maxZ)
0233       continue;
0234
0235     GlobalVector oldB(bx, by, bz);
0236     GlobalVector newB = field->inTesla(gp);
0237     if ((newB - oldB).mag() > reso) {
0238       ++fail;
0239       float delta = (newB - oldB).mag();
0240       if (delta > maxdelta)
0241         maxdelta = delta;
0242       if (fail < 10) {
0243         cout << " Discrepancy at: # " << count + 1 << " " << gp << " R " << gp.perp() << " Phi " << gp.phi()
0244              << " delta : " << newB - oldB << " " << delta << endl;
0245         const MagVolume6Faces* vol = findVolume(gp);
0246         if (vol)
0247           cout << " volume: " << vol->volumeNo << " " << (int)vol->copyno;
0248         cout << " Old: " << oldB << " New: " << newB << endl;
0249       } else if (fail == 10) {
0250         cout << "..." << endl;
0251       }
0252     }
0253     count++;
0254   } while (count < numberOfPoints);
0255
0256   if (count == 0) {
0257     edm::LogError("MagneticField") << "No input data" << endl;
0258   } else {
0259     cout << endl
0260          << " testMagneticField::validate: tested " << count << " points " << fail
0261          << " failures; max delta = " << maxdelta << endl
0262          << endl;
0263     if (fail != 0)
0264       throw cms::Exception("RegressionFailure") << "MF regression found: " << fail << " failures";
0265     ;
0266   }
0267 }
0268
0269 void testMagneticField::parseTOSCATablePath(string filename, int& volNo, int& sector, string& type) {
0270   // Determine volume number, type, and sector from filename, assumed to be like:
0271   // [path]/s01_1/v-xyz-1156.table
0272
0273   char buf[512];
0274   strcpy(buf, filename.c_str());
0275   string table = basename(buf);
0276   string ssector = basename(dirname(buf));
0277
0278   // Find type
0279   string::size_type ibeg = table.find('-');   // first occurence of "-"
0280   string::size_type iend = table.rfind('-');  // last  occurence of "-"
0281   type = table.substr(ibeg + 1, iend - ibeg - 1);
0282
0283   // Find volume number
0284   string::size_type iext = table.rfind('.');  // last  occurence of "."
0285   volNo = std::stoi(table.substr(iend + 1, iext - iend - 1));
0286
0287   // Find sector number
0288   if (ssector[0] == 's') {
0289     sector = std::stoi(ssector.substr(1, 2));
0290   } else {
0291     cout << "Can not determine sector number, assuming 1" << endl;
0292     sector = 1;
0293   }
0294 }
0295
0296 void testMagneticField::validateVsTOSCATable(string filename) {
0297   // The magic here is that we want to check against the result of the master volume only
0298   // as grid points on the border of volumes can end up in the neighbor volume.
0299
0300   int volNo, sector;
0301   string type;
0302   parseTOSCATablePath(filename, volNo, sector, type);
0303
0304   const MagVolume6Faces* vol = findMasterVolume(volNo, sector);
0305   if (vol == 0) {
0306     cout << "   ERROR: volume " << volNo << ":" << sector << "not found" << endl;
0307     return;
0308   }
0309
0310   cout << "Validate interpolation vs TOSCATable: " << filename << " volume " << volNo << ":[" << sector << "], type "
0311        << type << endl;
0312
0313   ifstream file(filename.c_str());
0314   string line;
0315
0316   int fail = 0;
0317   int count = 0;
0318
0319   float maxdelta = 0.;
0320
0321   // Dump table
0322   //   const MFGrid* interpolator = (const MFGrid*) vol->provider();
0323   //   Dimensions dim = interpolator->dimensions();
0324   //   for (int i=0; i<dim.w; ++i){
0325   //     for (int j=0; j<dim.h; ++j){
0326   //       for (int k=0; k<dim.d; ++k){
0327   //    cout << vol->toGlobal(interpolator->nodePosition(i,j,k)) << " " << vol->toGlobal(interpolator->nodeValue(i,j,k)) <<endl;
0328   //       }
0329   //     }
0330   //   }
0331
0332   while (getline(file, line) && count < numberOfPoints) {
0333     if (line == "" || line[0] == '#')
0334       continue;
0335     stringstream linestr;
0336     linestr << line;
0337     double px, py, pz;
0338     double bx, by, bz;
0339     linestr >> px >> py >> pz >> bx >> by >> bz;
0340     GlobalPoint gp;
0341     if (type == "rpz") {  // rpz file with units in m.
0342       gp = GlobalPoint(GlobalPoint::Cylindrical(px * 100., py, pz * 100.));
0343     } else if (type == "xyz") {  // xyz file with units in m.
0344       gp = GlobalPoint(px * 100., py * 100., pz * 100.);
0345     } else {
0346       cout << "validateVsTOSCATable: type " << type << " unknown " << endl;
0347       return;
0348     }
0349
0350     GlobalVector oldB(bx, by, bz);
0351     if (vol->inside(gp, 0.03)) {
0352       GlobalVector newB = vol->inTesla(gp);
0353       if ((newB - oldB).mag() > reso) {
0354         ++fail;
0355         float delta = (newB - oldB).mag();
0356         if (delta > maxdelta)
0357           maxdelta = delta;
0358         cout << " Discrepancy at: # " << count + 1 << " " << gp << " delta : " << newB - oldB << " " << delta << endl;
0359         cout << " Table: " << oldB << " Map: " << newB << endl;
0360       }
0361     } else {
0362       cout << "ERROR: grid point # " << count + 1 << " " << gp << " is not inside volume " << endl;
0363     }
0364
0365     count++;
0366   }
0367
0368   if (count == 0) {
0369     cout << "ERROR: input table not found" << endl;
0370   } else {
0371     cout << endl
0372          << " testMagneticField::validateVsTOSCATable: tested " << count << " points " << fail
0373          << " failures; max delta = " << maxdelta << endl
0374          << endl;
0375   }
0376 }
0377
0378 // #include <multimap>
0379 // typedef multimap<float, pair<int, int> > VolumesByDiscrepancy ;
0380
0381 #include <sys/types.h>
0382 #include <sys/stat.h>
0383 #include <unistd.h>
0384
0385 // Compare the TOSCA txt table with the corresponding one in other sector.
0386 void testMagneticField::compareSectorTables(string file1) {
0387   bool list = false;  // true: print one line per volume
0388                       // false: print formatted table
0389
0390   int volNo, sector1;
0391   string type;
0392   parseTOSCATablePath(file1, volNo, sector1, type);
0393
0394   //cout << "Comparing tables for all sectors for volume " << volNo << " with " << file1 << endl;
0395   if (!list)
0396     cout << volNo;
0397
0398   float maxmaxdelta = 0.;
0399   for (int sector2 = 1; sector2 <= 12; ++sector2) {
0400     if (sector2 == sector1) {
0401       if (!list)
0402         cout << "   ----";
0403       continue;
0404     }
0405
0406     //     vols.insert(pair<float,int>());
0407
0408     double phi = (sector2 - sector1) * Geom::pi() / 6.;
0409     double cphi = cos(phi);
0410     double sphi = sin(phi);
0411
0412     // Path of file for the other sector
0413     string file2 = file1;
0414     string::size_type ss = file2.rfind("/s");  // last  occurence of "-"
0415     string ssec = "/s";
0416     if (sector2 < 10)
0417       ssec += "0";
0418     ssec += std::to_string(sector2);
0419     file2.replace(ss, 4, ssec);
0420
0421     vector<GlobalPoint> p1, p2;
0422     vector<GlobalVector> b1, b2;
0423
0424     fillFromTable(file1, p1, b1, type);
0425     fillFromTable(file2, p2, b2, type);
0426
0427     struct stat theStat;
0428     //FIXME get table size
0429     string binTable = "/data/n/namapane/MagneticField/120812/grid_120812_3_8t_v7_large";
0430     binTable += ssec;
0431     binTable += "_";
0432     string sVolNo = std::to_string(volNo);
0433     binTable += sVolNo[0];
0434     binTable += "/grid.";
0435     binTable += sVolNo;
0436     binTable += ".bin";
0437     stat(binTable.c_str(), &theStat);
0438     off_t size = theStat.st_size;
0439
0440     if (p1.size() != p2.size() || p1.size() == 0) {
0441       cout << "ERROR: file size: " << p1.size() << " " << p2.size() << endl;
0442     }
0443
0444     float maxdelta = 0;
0445     float avgdelta = 0;
0446     //     int imaxdelta = -1;
0447     vector<GlobalVector> b12;
0448     for (unsigned int i = 0; i < p1.size(); ++i) {
0449       // check positions, need to get appropriate rotation
0450       //Rotate b1 into sector of b2
0451       GlobalPoint p12(cphi * p1[i].x() - sphi * p1[i].y(), sphi * p1[i].x() + cphi * p1[i].y(), p1[i].z());
0452       float pd = (p12 - p2[i]).mag();
0453       if (pd > 0.005) {
0454         cout << "ERROR: " << p12 << " " << p2[i] << " " << (p12 - p2[i]).mag() << endl;
0455       }
0456
0457       b12.push_back(GlobalVector(cphi * b1[i].x() - sphi * b1[i].y(), sphi * b1[i].x() + cphi * b1[i].y(), b1[i].z()));
0458       GlobalVector delta = (b2[i] - b12[i]);
0459       float d = delta.mag();
0460       avgdelta += d;
0461       if (d > maxdelta) {
0462         //  imaxdelta=i;
0463         maxdelta = d;
0464       }
0465     }
0466
0467     if (maxdelta > maxmaxdelta) {
0468       maxmaxdelta = maxdelta;
0469     }
0470
0471     avgdelta /= p1.size();
0472
0473     cout << setprecision(3) << fixed;
0474     if (list)
0475       cout << volNo << " " << sector2 << " " << avgdelta << " " << maxdelta << " " << size << endl;
0476     else {
0477       cout << "   " << maxdelta;
0478       //      cout << "   " << avgdelta;
0479     }
0480
0481     cout.unsetf(ios_base::floatfield);
0482     //     cout << "MAX: " << volNo << " " << sector2 << " " <<  setprecision(3) << maxdelta << endl;
0483     //     cout << imaxdelta << " " << b2[imaxdelta] << " " << b12[imaxdelta] << " " << (b2[imaxdelta]-b12[imaxdelta]) << endl;
0484   }
0485   cout << endl;
0486 }
0487
0488 void testMagneticField::fillFromTable(string inputFile, vector<GlobalPoint>& p, vector<GlobalVector>& b, string type) {
0489   ifstream file(inputFile.c_str());
0490   string line;
0491   while (getline(file, line)) {
0492     stringstream linestr;
0493     linestr << line;
0494     double px, py, pz;
0495     double bx, by, bz;
0496     linestr >> px >> py >> pz >> bx >> by >> bz;
0497     GlobalVector gv(bx, by, bz);
0498     GlobalPoint gp;
0499     if (type == "rpz") {  // rpz file with units in m.
0500       gp = GlobalPoint(GlobalPoint::Cylindrical(px * 100., py, pz * 100.));
0501     } else if (type == "xyz") {  // xyz file with units in m.
0502       gp = GlobalPoint(px * 100., py * 100., pz * 100.);
0503     } else {
0504       cout << "fillFromTable: type " << type << " unknown " << endl;
0505       return;
0506     }
0507     p.push_back(gp);
0508     b.push_back(gv);
0509   }
0510 }
0511
0512 #include "MagneticField/VolumeBasedEngine/interface/VolumeBasedMagneticField.h"
0513
0514 // Get the pointer of the volume containing a point
0515 const MagVolume6Faces* testMagneticField::findVolume(GlobalPoint& gp) {
0516   const VolumeBasedMagneticField* vbffield = dynamic_cast<const VolumeBasedMagneticField*>(field);
0517   if (vbffield) {
0518     return (dynamic_cast<const MagVolume6Faces*>(vbffield->findVolume(gp)));
0519   }
0520   return 0;
0521 }
0522
0523 // Find a specific volume:sector
0524 const MagVolume6Faces* testMagneticField::findMasterVolume(int volume, int sector) {
0525   const MagGeometry* vbffield = (dynamic_cast<const VolumeBasedMagneticField*>(field))->field;
0526
0527   if (vbffield == 0)
0528     return 0;
0529
0530   const vector<MagVolume6Faces const*>& bvol = vbffield->barrelVolumes();
0531   for (vector<MagVolume6Faces const*>::const_iterator i = bvol.begin(); i != bvol.end(); i++) {
0532     if ((*i)->copyno == sector && (*i)->volumeNo == volume) {
0533       return (*i);
0534     }
0535   }
0536
0537   const vector<MagVolume6Faces const*>& evol = vbffield->endcapVolumes();
0538   for (vector<MagVolume6Faces const*>::const_iterator i = evol.begin(); i != evol.end(); i++) {
0539     if ((*i)->copyno == sector && (*i)->volumeNo == volume) {
0540       return (*i);
0541     }
0542   }
0543
0544   return 0;
0545 }
0546
0547 DEFINE_FWK_MODULE(testMagneticField);