PrintGeomInfo/test/SimFileCompare.cpp

0001 ////////////////////////////////////////////////////////////////////////////////
0002 //
0003 //    Compares output files from PrintGeomInfo created using DDD and DD4hep
0004 //    inputs. Usage:
0005 //
0006 //    SimFileCompare infile1 infile2 type files debug
0007 //    infile1  (const char*)   First file name
0008 //    infile2  (const char*)   Second file name
0009 //    type     (int)           Type of file: material (0), solid (1),
0010 //                             LogicalVolume (2), PhysicalVolume (3);
0011 //                             Region (4)
0012 //    files    (int)           Double digits each inidicating the file source
0013 //                             (0 for DDD, 1 for DD4hep). So if first file is
0014 //                             DDD and second is DD4hep, it will be 10
0015 //    debug    (int)           Single digit number (0 minimum printout)
0016 //
0017 ////////////////////////////////////////////////////////////////////////////////
0018
0019 #include <algorithm>
0020 #include <cstdlib>
0021 #include <fstream>
0022 #include <iomanip>
0023 #include <iostream>
0024 #include <map>
0025 #include <string>
0026 #include <vector>
0027 #include "SimG4Core/Geometry/interface/DD4hep2DDDName.h"
0028
0029 struct materials {
0030   int occ;
0031   double radl, intl;
0032   materials(int oc = 1, double rd = 0, double in = 0) : occ(oc), radl(rd), intl(in) {}
0033 };
0034
0035 struct solids {
0036   int occ;
0037   double volume;
0038   solids(int oc = 1, double vol = 0) : occ(oc), volume(vol) {}
0039 };
0040
0041 struct lvs {
0042   int occ;
0043   double mass;
0044   lvs(int oc = 1, double m = 0) : occ(oc), mass(m) {}
0045 };
0046
0047 struct pvs {
0048   int occ;
0049   double xx, yy, zz;
0050   pvs(int oc = 1, double x = 0, double y = 0, double z = 0) : occ(oc), xx(x), yy(y), zz(z) {}
0051 };
0052
0053 struct regions {
0054   int occ;
0055   double nmat, nvol;
0056   regions(int oc = 1, double mat = 0, double vol = 0) : occ(oc), nmat(mat), nvol(vol) {}
0057 };
0058
0059 std::string removeExtraName(const std::string& name, int debug) {
0060   std::string nam(name);
0061   std::string nam1 = name.substr(0, 2);
0062   if (((nam1 == "GE") || (nam1 == "GH") || (nam1 == "MB") || (nam1 == "ME") || (nam1 == "RE") || (nam1 == "RR") ||
0063        (nam1 == "RT")) &&
0064       (name.size() > 5)) {
0065     uint32_t loc = name.size() - 5;
0066     if ((name.substr(0, 15) != "MBCables_Wheels") && (name.substr(loc, 1) == "_")) {
0067       std::string nam2 = (name.substr(loc + 3, 1) == "0") ? name.substr(loc + 4, 1) : name.substr(loc + 3, 2);
0068       nam = name.substr(0, loc + 1) + nam2;
0069     }
0070   }
0071   if (debug)
0072     std::cout << name << " : " << nam1 << " " << nam << std::endl;
0073   return nam;
0074 }
0075
0076 std::string reducedName(const std::string& name, int debug) {
0077   std::string nam(name);
0078   uint32_t first = ((name.find(":") == std::string::npos) ? 0 : (name.find(":") + 1));
0079   uint32_t last(name.size() + 1);
0080   uint32_t loc(first);
0081   while (1) {
0082     if (name.find("_", loc) == std::string::npos)
0083       break;
0084     if (((loc + 5) < name.size()) && (name.substr(loc, 5) == "shape")) {
0085       last = loc;
0086       break;
0087     }
0088     loc = name.find("_", loc) + 1;
0089     if (loc > name.size())
0090       break;
0091   }
0092   nam = name.substr(first, last - first - 1);
0093   if ((last < name.size()) && (name.substr(name.size() - 5, 5) == "_refl"))
0094     nam += "_refl";
0095   if (debug > 0)
0096     std::cout << name << " col " << first << ":" << last << " " << nam << std::endl;
0097   return nam;
0098 }
0099
0100 std::vector<std::string> splitString(const std::string& fLine) {
0101   std::vector<std::string> result;
0102   int start = 0;
0103   bool empty = true;
0104   for (unsigned i = 0; i <= fLine.size(); i++) {
0105     if (fLine[i] == ' ' || i == fLine.size()) {
0106       if (!empty) {
0107         std::string item(fLine, start, i - start);
0108         result.push_back(item);
0109         empty = true;
0110       }
0111       start = i + 1;
0112     } else {
0113       if (empty)
0114         empty = false;
0115     }
0116   }
0117   return result;
0118 }
0119
0120 template <typename T>
0121 void myPrint1(std::map<std::string, T> const& obj) {
0122   for (auto it : obj) {
0123     if (it.second.occ > 1)
0124       std::cout << it.first << " : " << it.second.occ << std::endl;
0125   }
0126 }
0127
0128 template <typename T>
0129 void myPrint2(std::map<std::string, T> const& obj1, std::map<std::string, T> const& obj2) {
0130   for (auto it : obj1) {
0131     if (obj2.find(it.first) == obj2.end())
0132       std::cout << it.first << " appearing " << it.second.occ << " times" << std::endl;
0133   }
0134 }
0135
0136 void CompareFiles(const char* fileFile1, const char* fileFile2, int type, int files, int debug) {
0137   std::map<std::string, materials> matFile1, matFile2;
0138   std::map<std::string, solids> solidFile1, solidFile2;
0139   std::map<std::string, lvs> lvFile1, lvFile2;
0140   std::map<std::string, pvs> pvFile1, pvFile2;
0141   std::map<std::string, regions> regFile1, regFile2;
0142   bool typeFile1 = ((files % 10) == 0);
0143   bool typeFile2 = (((files / 10) % 10) == 0);
0144   char buffer[100];
0145   std::string name;
0146   std::ifstream fInput1(fileFile1);
0147   unsigned int sizeFile1(0), sizeFile2(0);
0148   if (!fInput1.good()) {
0149     std::cout << "Cannot open file " << fileFile1 << std::endl;
0150   } else {
0151     while (fInput1.getline(buffer, 100)) {
0152       std::vector<std::string> items = splitString(std::string(buffer));
0153       if ((type == 0) || (type == 2))
0154         name = DD4hep2DDDName::nameMatterLV(items[0], !typeFile1);
0155       else if (type == 1)
0156         name = DD4hep2DDDName::nameSolid(items[0], !typeFile1);
0157       else if (type == 3)
0158         name = DD4hep2DDDName::namePV(items[0], !typeFile1);
0159       else
0160         name = items[0];
0161       double r1 = (items.size() > 1) ? atof(items[1].c_str()) : 0;
0162       double r2 = (items.size() > 2) ? atof(items[2].c_str()) : 0;
0163       double r3 = (items.size() > 3) ? atof(items[3].c_str()) : 0;
0164       if (type == 0) {
0165         auto it = matFile1.find(name);
0166         if (it == matFile1.end())
0167           matFile1[name] = materials(1, r1, r2);
0168         else
0169           ++((it->second).occ);
0170       } else if (type == 1) {
0171         auto it = solidFile1.find(name);
0172         if (it == solidFile1.end())
0173           solidFile1[name] = solids(1, r1);
0174         else
0175           ++((it->second).occ);
0176       } else if (type == 2) {
0177         auto it = lvFile1.find(name);
0178         if (it == lvFile1.end())
0179           lvFile1[name] = lvs(1, r1);
0180         else
0181           ++((it->second).occ);
0182       } else if (type == 3) {
0183         auto it = pvFile1.find(name);
0184         if (it == pvFile1.end())
0185           pvFile1[name] = pvs(1, r1, r2, r3);
0186         else
0187           ++((it->second).occ);
0188       } else {
0189         auto it = regFile1.find(name);
0190         if (it == regFile1.end())
0191           regFile1[name] = regions(1, r1, r2);
0192         else
0193           ++((it->second).occ);
0194       }
0195     }
0196     fInput1.close();
0197     sizeFile1 = ((type == 0) ? matFile1.size()
0198                              : ((type == 1) ? solidFile1.size() : ((type == 2) ? lvFile1.size() : pvFile1.size())));
0199   }
0200   std::ifstream fInput2(fileFile2);
0201   if (!fInput2.good()) {
0202     std::cout << "Cannot open file " << fileFile2 << std::endl;
0203   } else {
0204     while (fInput2.getline(buffer, 100)) {
0205       std::vector<std::string> items = splitString(std::string(buffer));
0206       if ((type == 0) || (type == 2))
0207         name = DD4hep2DDDName::nameMatterLV(items[0], !typeFile2);
0208       else if (type == 1)
0209         name = DD4hep2DDDName::nameSolid(items[0], !typeFile2);
0210       else if (type == 3)
0211         name = DD4hep2DDDName::namePV(items[0], !typeFile2);
0212       else
0213         name = items[0];
0214       double r1 = (items.size() > 1) ? atof(items[1].c_str()) : 0;
0215       double r2 = (items.size() > 2) ? atof(items[2].c_str()) : 0;
0216       double r3 = (items.size() > 3) ? atof(items[3].c_str()) : 0;
0217       if (type == 0) {
0218         auto it = matFile2.find(name);
0219         if (it == matFile2.end())
0220           matFile2[name] = materials(1, r1, r2);
0221         else
0222           ++((it->second).occ);
0223       } else if (type == 1) {
0224         auto it = solidFile2.find(name);
0225         if (it == solidFile2.end())
0226           solidFile2[name] = solids(1, r1);
0227         else
0228           ++((it->second).occ);
0229       } else if (type == 2) {
0230         auto it = lvFile2.find(name);
0231         if (it == lvFile2.end())
0232           lvFile2[name] = lvs(1, r1);
0233         else
0234           ++((it->second).occ);
0235       } else if (type == 3) {
0236         auto it = pvFile2.find(name);
0237         if (it == pvFile2.end())
0238           pvFile2[name] = pvs(1, r1, r2, r3);
0239         else
0240           ++((it->second).occ);
0241       } else {
0242         auto it = regFile2.find(name);
0243         if (it == regFile2.end())
0244           regFile2[name] = regions(1, r1, r2);
0245         else
0246           ++((it->second).occ);
0247       }
0248     }
0249     fInput2.close();
0250     sizeFile2 = ((type == 0) ? matFile2.size()
0251                              : ((type == 1) ? solidFile2.size() : ((type == 2) ? lvFile2.size() : pvFile2.size())));
0252   }
0253   std::cout << "Reads " << sizeFile1 << " names from " << fileFile1 << " and " << sizeFile2 << " names from "
0254             << fileFile2 << std::endl;
0255
0256   std::cout << "\nMore than one entry for a given name in " << fileFile1 << std::endl;
0257   if (type == 0) {
0258     myPrint1(matFile1);
0259   } else if (type == 1) {
0260     myPrint1(solidFile1);
0261   } else if (type == 2) {
0262     myPrint1(lvFile1);
0263   } else if (type == 3) {
0264     myPrint1(pvFile1);
0265   } else {
0266     myPrint1(regFile1);
0267   }
0268
0269   std::cout << "\nMore than one entry for a given name in " << fileFile2 << std::endl;
0270   if (type == 0) {
0271     myPrint1(matFile2);
0272   } else if (type == 1) {
0273     myPrint1(solidFile2);
0274   } else if (type == 2) {
0275     myPrint1(lvFile2);
0276   } else if (type == 3) {
0277     myPrint1(pvFile2);
0278   } else {
0279     myPrint1(regFile2);
0280   }
0281
0282   std::cout << "\nEntry in " << fileFile1 << " not in " << fileFile2 << std::endl;
0283   if (type == 0) {
0284     myPrint2(matFile1, matFile2);
0285   } else if (type == 1) {
0286     myPrint2(solidFile1, solidFile2);
0287   } else if (type == 2) {
0288     myPrint2(lvFile1, lvFile2);
0289   } else if (type == 3) {
0290     myPrint2(pvFile1, pvFile2);
0291   } else {
0292     myPrint2(regFile1, regFile2);
0293   }
0294
0295   std::cout << "\nEntry in " << fileFile2 << " not in " << fileFile1 << std::endl;
0296   if (type == 0) {
0297     myPrint2(matFile2, matFile1);
0298   } else if (type == 1) {
0299     myPrint2(solidFile2, solidFile1);
0300   } else if (type == 2) {
0301     myPrint2(lvFile2, lvFile1);
0302   } else if (type == 2) {
0303     myPrint2(pvFile2, pvFile1);
0304   } else {
0305     myPrint2(regFile2, regFile1);
0306   }
0307
0308   //Now type specific changes
0309   std::cout << "\nEntries in " << fileFile1 << " and " << fileFile2 << " do not match in the content\n";
0310   const double denmin = 0.0001;
0311   int kount1(0), kount2(0);
0312   double difmax1(0), difmax2(0);
0313   std::string nameMax("");
0314   if (type == 0) {
0315     const double tol1 = 0.00001;
0316     for (auto it1 : matFile1) {
0317       auto it2 = matFile2.find(it1.first);
0318       if (it2 != matFile2.end()) {
0319         ++kount1;
0320         double rdif =
0321             0.5 * (it1.second.radl - it2->second.radl) / std::max(denmin, (it1.second.radl + it2->second.radl));
0322         double idif =
0323             0.5 * (it1.second.intl - it2->second.intl) / std::max(denmin, (it1.second.intl + it2->second.intl));
0324         if (std::abs(rdif) > difmax1) {
0325           difmax1 = std::abs(rdif);
0326           difmax2 = std::abs(idif);
0327           nameMax = it1.first;
0328         }
0329         if ((std::abs(rdif) > tol1) || (std::abs(idif) > tol1)) {
0330           ++kount2;
0331           std::cout << it1.first << " X0 " << it1.second.radl << ":" << it2->second.radl << ":" << rdif << " #L "
0332                     << it1.second.intl << ":" << it2->second.intl << ":" << idif << std::endl;
0333         }
0334       }
0335     }
0336     std::cout << "\n " << kount2 << " out of " << kount1 << " entries having discrpancies at the level of " << tol1
0337               << " or more; the maximum happens for " << nameMax << " with " << difmax1 << ":" << difmax2 << "\n";
0338   } else if (type == 1) {
0339     const double tol2 = 0.0001;
0340     for (auto it1 : solidFile1) {
0341       auto it2 = solidFile2.find(it1.first);
0342       if (it2 != solidFile2.end()) {
0343         ++kount1;
0344         double vdif =
0345             0.5 * (it1.second.volume - it2->second.volume) / std::max(denmin, (it1.second.volume + it2->second.volume));
0346         if (std::abs(vdif) > difmax1) {
0347           difmax1 = std::abs(vdif);
0348           nameMax = it1.first;
0349         }
0350         if (std::abs(vdif) > tol2) {
0351           ++kount2;
0352           std::cout << it1.first << " Volume " << it1.second.volume << ":" << it2->second.volume << ":" << vdif
0353                     << std::endl;
0354         }
0355       }
0356     }
0357     std::cout << "\n " << kount2 << " out of " << kount1 << " entries having discrpancies at the level of " << tol2
0358               << " or more; the maximum happens for " << nameMax << " with " << difmax1 << "\n";
0359   } else if (type == 2) {
0360     const double tol3 = 0.0001;
0361     for (auto it1 : lvFile1) {
0362       auto it2 = lvFile2.find(it1.first);
0363       if (it2 != lvFile2.end()) {
0364         ++kount1;
0365         double vdif =
0366             0.5 * (it1.second.mass - it2->second.mass) / std::max(denmin, (it1.second.mass + it2->second.mass));
0367         if (std::abs(vdif) > difmax1) {
0368           difmax1 = std::abs(vdif);
0369           nameMax = it1.first;
0370         }
0371         if (std::abs(vdif) > tol3) {
0372           ++kount2;
0373           std::cout << it1.first << " Mass " << it1.second.mass << ":" << it2->second.mass << ":" << vdif << std::endl;
0374         }
0375       }
0376     }
0377     std::cout << "\n " << kount2 << " out of " << kount1 << " entries having discrpancies at the level of " << tol3
0378               << " or more; the maximum happens for " << nameMax << " with " << difmax1 << "\n";
0379   } else if (type == 3) {
0380     const double tol4 = 0.0001;
0381     for (auto it1 : pvFile1) {
0382       auto it2 = pvFile2.find(it1.first);
0383       if (it2 != pvFile2.end()) {
0384         ++kount1;
0385         double xdif = (it1.second.xx - it2->second.xx);
0386         double ydif = (it1.second.yy - it2->second.yy);
0387         double zdif = (it1.second.zz - it2->second.zz);
0388         double vdif = std::max(std::abs(xdif), std::abs(ydif));
0389         vdif = std::max(vdif, std::abs(zdif));
0390         if (vdif > difmax1) {
0391           difmax1 = vdif;
0392           nameMax = it1.first;
0393         }
0394         if ((std::abs(xdif) > tol4) || (std::abs(ydif) > tol4) || (std::abs(zdif) > tol4)) {
0395           ++kount2;
0396           std::cout << it1.first << " x " << it1.second.xx << ":" << it2->second.xx << ":" << xdif << " y "
0397                     << it1.second.yy << ":" << it2->second.yy << ":" << ydif << " z " << it1.second.zz << ":"
0398                     << it2->second.zz << ":" << zdif << std::endl;
0399         }
0400       }
0401     }
0402     std::cout << "\n " << kount2 << " out of " << kount1 << " entries having discrpancies at the level of " << tol4
0403               << " or more; the maximum happens for " << nameMax << " with " << difmax1 << "\n";
0404   } else {
0405     const double tol5 = 0.0001;
0406     for (auto it1 : regFile1) {
0407       auto it2 = regFile2.find(it1.first);
0408       if (it2 != regFile2.end()) {
0409         ++kount1;
0410         double matdif = (it1.second.nmat - it2->second.nmat);
0411         double voldif = (it1.second.nvol - it2->second.nvol);
0412         if (std::abs(matdif) > difmax1) {
0413           difmax1 = std::abs(matdif);
0414           nameMax = it1.first;
0415         }
0416         if (std::abs(voldif) > difmax2) {
0417           difmax2 = std::abs(voldif);
0418           nameMax = it1.first;
0419         }
0420         if ((std::abs(matdif) > tol5) || (std::abs(voldif) > tol5)) {
0421           ++kount2;
0422           std::cout << it1.first << " Material " << it1.second.nmat << ":" << it2->second.nmat << ":" << matdif
0423                     << " Volume " << it1.second.nvol << ":" << it2->second.nvol << ":" << voldif << std::endl;
0424         }
0425       }
0426     }
0427     std::cout << "\n " << kount2 << " out of " << kount1 << " entries having discrpancies at the level of " << tol5
0428               << " or more; the maximum happens for " << nameMax << " with " << difmax1 << ":" << difmax2 << "\n";
0429   }
0430 }
0431
0432 int main(int argc, char* argv[]) {
0433   if (argc <= 5) {
0434     std::cout << "Please give a minimum of 2 arguments \n"
0435               << "name of the first input file\n"
0436               << "name of the second input file\n"
0437               << "type (Material:0, Solid:1, LV:2, PV:3, Region:4)\n"
0438               << "files (10 if first file from DDD and second from DD4hep)\n"
0439               << "debug flag (0 for minimum printout)\n"
0440               << std::endl;
0441     return 0;
0442   }
0443
0444   const char* infile1 = argv[1];
0445   const char* infile2 = argv[2];
0446   int type = ((argc > 3) ? atoi(argv[3]) : 0);
0447   if (type < 0 || type > 3)
0448     type = 0;
0449   int files = ((argc > 4) ? atoi(argv[4]) : 10);
0450   int debug = ((argc > 6) ? atoi(argv[6]) : 0);
0451   CompareFiles(infile1, infile2, type, files, debug);
0452   return 0;
0453 }