AlignmentRelative/test/ppsTrackBasedAligAnalRes.cc

0001 /****************************************************************************
0002 * Authors:
0003 *  Jan Kašpar (jan.kaspar@gmail.com)
0004 ****************************************************************************/
0005
0006 #include <dirent.h>
0007 #include <map>
0008 #include <string>
0009
0010 #include "TFile.h"
0011 #include "TGraph.h"
0012 #include "TGraphErrors.h"
0013
0014 #include "FWCore/Utilities/interface/Exception.h"
0015
0016 #include "CondFormats/PPSObjects/interface/CTPPSRPAlignmentCorrectionsData.h"
0017 #include "CondFormats/PPSObjects/interface/CTPPSRPAlignmentCorrectionsMethods.h"
0018
0019 using namespace std;
0020
0021 //----------------------------------------------------------------------------------------------------
0022
0023 TDirectory *mkdir(TDirectory *parent, const char *child) {
0024   TDirectory *dir = (TDirectory *)parent->Get(child);
0025   if (!dir)
0026     dir = parent->mkdir(child);
0027   return dir;
0028 }
0029
0030 //----------------------------------------------------------------------------------------------------
0031
0032 struct Stat {
0033   double s1, sx, sxx, sxxx, sxxxx;
0034
0035   Stat() : s1(0.), sx(0.), sxx(0.), sxxx(0.), sxxxx(0.) {}
0036
0037   void fill(double x) {
0038     s1 += 1;
0039     sx += x;
0040     sxx += x * x;
0041     sxxx += x * x * x;
0042     sxxxx += x * x * x * x;
0043   }
0044
0045   double n() const { return s1; }
0046
0047   double mean() const { return sx / s1; }
0048
0049   double meanError() const {
0050     double v = rms() / sqrt(s1);
0051     return v;
0052   }
0053
0054   double rms() const {
0055     double sig_sq = (sxx - sx * sx / s1) / (s1 - 1.);
0056     if (sig_sq < 0)
0057       sig_sq = 0;
0058
0059     return sqrt(sig_sq);
0060   }
0061
0062   double rmsError() const {
0063     // see R.J.Barlow: Statistics, page 78
0064     double mu = mean();
0065     double E2 = rms();
0066     E2 *= E2;
0067     //if (print) printf("\t\t%E, %E, %E, %E, %E\n", sxxxx, -4.*sxxx*mu, 6.*sxx*mu*mu, -4.*sx*mu*mu*mu, s1*mu*mu*mu*mu);
0068     double E4 = (sxxxx - 4. * sxxx * mu + 6. * sxx * mu * mu - 4. * sx * mu * mu * mu + s1 * mu * mu * mu * mu) / s1;
0069     //if (print) printf("\t\tmu = %E, E2 = %E, E4 = %E\n", mu, E2, E4);
0070     double v_sig_sq = (s1 - 1) * ((s1 - 1) * E4 - (s1 - 3) * E2 * E2) / s1 / s1 / s1;
0071     double v_sig = v_sig_sq / 4. / E2;
0072     //if (print) printf("\t\tv_sig_sq = %E\n", v_sig_sq);
0073     if (v_sig < 0 || std::isnan(v_sig) || std::isinf(v_sig)) {
0074       //printf(">> Stat::rmsError > ERROR: v_siq = %E < 0.\n", v_sig);
0075       v_sig = 0;
0076     }
0077     double v = sqrt(v_sig);
0078     return v;
0079   }
0080 };
0081
0082 //----------------------------------------------------------------------------------------------------
0083
0084 struct DetStat {
0085   // _v ... value given by Jan algorithm
0086   // _u ... uncertainty estimated by Jan algorithm
0087   // _e ... difference (Jan - Ideal) algorithm
0088   Stat sh_x_v, sh_x_u, sh_x_e;
0089   Stat sh_y_v, sh_y_u, sh_y_e;
0090   Stat rot_z_v, rot_z_u, rot_z_e;
0091   Stat sh_z_v, sh_z_u, sh_z_e;
0092 };
0093
0094 //----------------------------------------------------------------------------------------------------
0095
0096 struct RPStat {
0097   Stat sh_x_v, sh_x_u, sh_x_e;
0098   Stat sh_y_v, sh_y_u, sh_y_e;
0099   Stat rot_z_v, rot_z_u, rot_z_e;
0100   Stat sh_z_v, sh_z_u, sh_z_e;
0101 };
0102
0103 //----------------------------------------------------------------------------------------------------
0104
0105 struct Desc {
0106   unsigned int N;  // number of events
0107   unsigned int i;  // iteration
0108   unsigned int d;  // detector
0109   Desc(unsigned int _N, unsigned int _i, unsigned int _d) : N(_N), i(_i), d(_d) {}
0110   bool operator<(const Desc &c) const;
0111 };
0112
0113 bool Desc::operator<(const Desc &c) const {
0114   if (this->N < c.N)
0115     return true;
0116   if (this->N > c.N)
0117     return false;
0118   if (this->i < c.i)
0119     return true;
0120   if (this->i > c.i)
0121     return false;
0122   if (this->d < c.d)
0123     return true;
0124   if (this->d > c.d)
0125     return false;
0126   return false;
0127 }
0128
0129 //----------------------------------------------------------------------------------------------------
0130
0131 map<Desc, DetStat> det_stat;
0132 map<Desc, RPStat> rp_stat;
0133
0134 //----------------------------------------------------------------------------------------------------
0135
0136 void resetStatistics() {
0137   det_stat.clear();
0138   rp_stat.clear();
0139 }
0140
0141 //----------------------------------------------------------------------------------------------------
0142
0143 void updateSensorStatistics(unsigned int n_events,
0144                             unsigned iteration,
0145                             const CTPPSRPAlignmentCorrectionsData &r_actual,
0146                             const CTPPSRPAlignmentCorrectionsData &r_ideal) {
0147   for (CTPPSRPAlignmentCorrectionsData::mapType::const_iterator it = r_actual.getSensorMap().begin();
0148        it != r_actual.getSensorMap().end();
0149        ++it) {
0150     unsigned int id = it->first;
0151     const auto &c_actual = r_actual.getFullSensorCorrection(id, false);
0152     const auto &c_ideal = r_ideal.getFullSensorCorrection(id, false);
0153
0154     DetStat &s = det_stat[Desc(n_events, iteration, id)];
0155
0156     s.sh_x_v.fill(c_actual.getShX());
0157     s.sh_x_u.fill(c_actual.getShXUnc());
0158     s.sh_x_e.fill(c_actual.getShX() - c_ideal.getShX());
0159
0160     s.sh_y_v.fill(c_actual.getShY());
0161     s.sh_y_u.fill(c_actual.getShYUnc());
0162     s.sh_y_e.fill(c_actual.getShY() - c_ideal.getShY());
0163
0164     s.sh_z_v.fill(c_actual.getShZ());
0165     s.sh_z_u.fill(c_actual.getShZUnc());
0166     s.sh_z_e.fill(c_actual.getShZ() - c_ideal.getShZ());
0167
0168     s.rot_z_v.fill(c_actual.getRotZ());
0169     s.rot_z_u.fill(c_actual.getRotZUnc());
0170     s.rot_z_e.fill(c_actual.getRotZ() - c_ideal.getRotZ());
0171   }
0172 }
0173
0174 //----------------------------------------------------------------------------------------------------
0175
0176 void updateRPStatistics(unsigned int n_events,
0177                         unsigned iteration,
0178                         const CTPPSRPAlignmentCorrectionsData &r_actual,
0179                         const CTPPSRPAlignmentCorrectionsData &r_ideal) {
0180   for (CTPPSRPAlignmentCorrectionsData::mapType::const_iterator it = r_actual.getRPMap().begin();
0181        it != r_actual.getRPMap().end();
0182        ++it) {
0183     unsigned int id = it->first;
0184     const auto &c_actual = r_actual.getRPCorrection(id);
0185     const auto &c_ideal = r_ideal.getRPCorrection(id);
0186
0187     RPStat &s = rp_stat[Desc(n_events, iteration, id)];
0188
0189     s.sh_x_v.fill(c_actual.getShX());
0190     s.sh_x_u.fill(c_actual.getShXUnc());
0191     s.sh_x_e.fill(c_actual.getShX() - c_ideal.getShX());
0192
0193     s.sh_y_v.fill(c_actual.getShY());
0194     s.sh_y_u.fill(c_actual.getShYUnc());
0195     s.sh_y_e.fill(c_actual.getShY() - c_ideal.getShY());
0196
0197     s.sh_z_v.fill(c_actual.getShZ());
0198     s.sh_z_u.fill(c_actual.getShZUnc());
0199     s.sh_z_e.fill(c_actual.getShZ() - c_ideal.getShZ());
0200
0201     s.rot_z_v.fill(c_actual.getRotZ());
0202     s.rot_z_u.fill(c_actual.getRotZUnc());
0203     s.rot_z_e.fill(c_actual.getRotZ() - c_ideal.getRotZ());
0204   }
0205 }
0206
0207 //----------------------------------------------------------------------------------------------------
0208
0209 struct StatGraphs {
0210   TGraph *n;
0211   TGraphErrors *v_m, *v_v, *u_m, *u_v, *e_m, *e_v, *eR;
0212   StatGraphs()
0213       : n(new TGraph()),
0214         v_m(new TGraphErrors()),
0215         v_v(new TGraphErrors()),
0216         u_m(new TGraphErrors()),
0217         u_v(new TGraphErrors()),
0218         e_m(new TGraphErrors()),
0219         e_v(new TGraphErrors()),
0220         eR(new TGraphErrors()) {}
0221
0222   void write(const char *xLabel);
0223 };
0224
0225 //----------------------------------------------------------------------------------------------------
0226
0227 #define ENTRY(tag, label)             \
0228   tag->SetName(#tag);                 \
0229   sprintf(buf, ";%s;" label, xLabel); \
0230   tag->SetTitle(buf);                 \
0231   tag->Write();
0232
0233 void StatGraphs::write(const char *xLabel) {
0234   char buf[50];
0235   ENTRY(n, "number of repetitions");
0236   ENTRY(v_m, "value mean");
0237   ENTRY(v_v, "value variation");
0238   ENTRY(u_m, "estim. uncertainty mean");
0239   ENTRY(u_v, "estim. uncertainty variation");
0240   ENTRY(e_m, "error mean");
0241   ENTRY(e_v, "error variation");
0242   ENTRY(eR, "error variation / estim. uncertainty");
0243 }
0244
0245 //----------------------------------------------------------------------------------------------------
0246
0247 struct DetGraphs {
0248   StatGraphs sh_x, sh_y, rot_z, sh_z;
0249   void fill(double x, const DetStat &);
0250   void write(const char *xLabel);
0251 };
0252
0253 //----------------------------------------------------------------------------------------------------
0254
0255 double eR_error(const Stat &e, const Stat &u) {
0256   double a = e.rms(), ae = e.rmsError();
0257   double b = u.mean(), be = u.meanError();
0258
0259   return (b <= 0) ? 0. : a / b * sqrt(ae * ae / a / a + be * be / b / b);
0260 }
0261
0262 //----------------------------------------------------------------------------------------------------
0263
0264 void DetGraphs::fill(double x, const DetStat &s) {
0265   int idx = sh_x.n->GetN();
0266
0267   sh_x.n->SetPoint(idx, x, s.sh_x_u.n());
0268   sh_x.v_m->SetPoint(idx, x, s.sh_x_v.mean());
0269   sh_x.v_v->SetPoint(idx, x, s.sh_x_v.rms());
0270   sh_x.u_m->SetPoint(idx, x, s.sh_x_u.mean());
0271   sh_x.u_v->SetPoint(idx, x, s.sh_x_u.rms());
0272   sh_x.e_m->SetPoint(idx, x, s.sh_x_e.mean());
0273   sh_x.e_v->SetPoint(idx, x, s.sh_x_e.rms());
0274   sh_x.eR->SetPoint(idx, x, s.sh_x_e.rms() / s.sh_x_u.mean());
0275
0276   sh_x.v_m->SetPointError(idx, 0., s.sh_x_v.meanError());
0277   sh_x.v_v->SetPointError(idx, 0., s.sh_x_v.rmsError());
0278   sh_x.u_m->SetPointError(idx, 0., s.sh_x_u.meanError());
0279   sh_x.u_v->SetPointError(idx, 0., s.sh_x_u.rmsError());
0280   sh_x.e_m->SetPointError(idx, 0., s.sh_x_e.meanError());
0281   sh_x.e_v->SetPointError(idx, 0., s.sh_x_e.rmsError());
0282   sh_x.eR->SetPointError(idx, 0., eR_error(s.sh_x_e, s.sh_x_u));
0283
0284   sh_y.n->SetPoint(idx, x, s.sh_y_u.n());
0285   sh_y.v_m->SetPoint(idx, x, s.sh_y_v.mean());
0286   sh_y.v_v->SetPoint(idx, x, s.sh_y_v.rms());
0287   sh_y.u_m->SetPoint(idx, x, s.sh_y_u.mean());
0288   sh_y.u_v->SetPoint(idx, x, s.sh_y_u.rms());
0289   sh_y.e_m->SetPoint(idx, x, s.sh_y_e.mean());
0290   sh_y.e_v->SetPoint(idx, x, s.sh_y_e.rms());
0291   sh_y.eR->SetPoint(idx, x, s.sh_y_e.rms() / s.sh_y_u.mean());
0292
0293   sh_y.v_m->SetPointError(idx, 0., s.sh_y_v.meanError());
0294   sh_y.v_v->SetPointError(idx, 0., s.sh_y_v.rmsError());
0295   sh_y.u_m->SetPointError(idx, 0., s.sh_y_u.meanError());
0296   sh_y.u_v->SetPointError(idx, 0., s.sh_y_u.rmsError());
0297   sh_y.e_m->SetPointError(idx, 0., s.sh_y_e.meanError());
0298   sh_y.e_v->SetPointError(idx, 0., s.sh_y_e.rmsError());
0299   sh_y.eR->SetPointError(idx, 0., eR_error(s.sh_y_e, s.sh_y_u));
0300
0301   rot_z.n->SetPoint(idx, x, s.rot_z_u.n());
0302   rot_z.v_m->SetPoint(idx, x, s.rot_z_v.mean());
0303   rot_z.v_v->SetPoint(idx, x, s.rot_z_v.rms());
0304   rot_z.u_m->SetPoint(idx, x, s.rot_z_u.mean());
0305   rot_z.u_v->SetPoint(idx, x, s.rot_z_u.rms());
0306   rot_z.e_m->SetPoint(idx, x, s.rot_z_e.mean());
0307   rot_z.e_v->SetPoint(idx, x, s.rot_z_e.rms());
0308   rot_z.eR->SetPoint(idx, x, s.rot_z_e.rms() / s.rot_z_u.mean());
0309
0310   rot_z.v_m->SetPointError(idx, 0., s.rot_z_v.meanError());
0311   rot_z.v_v->SetPointError(idx, 0., s.rot_z_v.rmsError());
0312   rot_z.u_m->SetPointError(idx, 0., s.rot_z_u.meanError());
0313   rot_z.u_v->SetPointError(idx, 0., s.rot_z_u.rmsError());
0314   rot_z.e_m->SetPointError(idx, 0., s.rot_z_e.meanError());
0315   rot_z.e_v->SetPointError(idx, 0., s.rot_z_e.rmsError());
0316   rot_z.eR->SetPointError(idx, 0., eR_error(s.rot_z_e, s.rot_z_u));
0317
0318   sh_z.n->SetPoint(idx, x, s.sh_z_u.n());
0319   sh_z.v_m->SetPoint(idx, x, s.sh_z_v.mean());
0320   sh_z.v_v->SetPoint(idx, x, s.sh_z_v.rms());
0321   sh_z.u_m->SetPoint(idx, x, s.sh_z_u.mean());
0322   sh_z.u_v->SetPoint(idx, x, s.sh_z_u.rms());
0323   sh_z.e_m->SetPoint(idx, x, s.sh_z_e.mean());
0324   sh_z.e_v->SetPoint(idx, x, s.sh_z_e.rms());
0325   sh_z.eR->SetPoint(idx, x, s.sh_z_e.rms() / s.sh_z_u.mean());
0326
0327   sh_z.v_m->SetPointError(idx, 0., s.sh_z_v.meanError());
0328   sh_z.v_v->SetPointError(idx, 0., s.sh_z_v.rmsError());
0329   sh_z.u_m->SetPointError(idx, 0., s.sh_z_u.meanError());
0330   sh_z.u_v->SetPointError(idx, 0., s.sh_z_u.rmsError());
0331   sh_z.e_m->SetPointError(idx, 0., s.sh_z_e.meanError());
0332   sh_z.e_v->SetPointError(idx, 0., s.sh_z_e.rmsError());
0333   sh_z.eR->SetPointError(idx, 0., eR_error(s.sh_z_e, s.sh_z_u));
0334 }
0335
0336 //----------------------------------------------------------------------------------------------------
0337
0338 void DetGraphs::write(const char *xLabel) {
0339   gDirectory = mkdir(gDirectory, "sh_x");
0340   sh_x.write(xLabel);
0341   gDirectory->cd("..");
0342   gDirectory = mkdir(gDirectory, "sh_y");
0343   sh_y.write(xLabel);
0344   gDirectory->cd("..");
0345   gDirectory = mkdir(gDirectory, "rot_z");
0346   rot_z.write(xLabel);
0347   gDirectory->cd("..");
0348   gDirectory = mkdir(gDirectory, "sh_z");
0349   sh_z.write(xLabel);
0350   gDirectory->cd("..");
0351 }
0352 //----------------------------------------------------------------------------------------------------
0353
0354 struct RPGraphs {
0355   StatGraphs sh_x, sh_y, rot_z, sh_z;
0356   void fill(double x, const RPStat &);
0357   void write(const char *xLabel);
0358 };
0359
0360 //----------------------------------------------------------------------------------------------------
0361
0362 void RPGraphs::fill(double x, const RPStat &s) {
0363   int idx = sh_x.n->GetN();
0364
0365   sh_x.n->SetPoint(idx, x, s.sh_x_u.n());
0366   sh_x.v_m->SetPoint(idx, x, s.sh_x_v.mean());
0367   sh_x.v_v->SetPoint(idx, x, s.sh_x_v.rms());
0368   sh_x.u_m->SetPoint(idx, x, s.sh_x_u.mean());
0369   sh_x.u_v->SetPoint(idx, x, s.sh_x_u.rms());
0370   sh_x.e_m->SetPoint(idx, x, s.sh_x_e.mean());
0371   sh_x.e_v->SetPoint(idx, x, s.sh_x_e.rms());
0372   sh_x.eR->SetPoint(idx, x, s.sh_x_e.rms() / s.sh_x_u.mean());
0373
0374   sh_x.v_m->SetPointError(idx, 0., s.sh_x_v.meanError());
0375   sh_x.v_v->SetPointError(idx, 0., s.sh_x_v.rmsError());
0376   sh_x.u_m->SetPointError(idx, 0., s.sh_x_u.meanError());
0377   sh_x.u_v->SetPointError(idx, 0., s.sh_x_u.rmsError());
0378   sh_x.e_m->SetPointError(idx, 0., s.sh_x_e.meanError());
0379   sh_x.e_v->SetPointError(idx, 0., s.sh_x_e.rmsError());
0380   sh_x.eR->SetPointError(idx, 0., eR_error(s.sh_x_e, s.sh_x_u));
0381
0382   sh_y.n->SetPoint(idx, x, s.sh_y_u.n());
0383   sh_y.v_m->SetPoint(idx, x, s.sh_y_v.mean());
0384   sh_y.v_v->SetPoint(idx, x, s.sh_y_v.rms());
0385   sh_y.u_m->SetPoint(idx, x, s.sh_y_u.mean());
0386   sh_y.u_v->SetPoint(idx, x, s.sh_y_u.rms());
0387   sh_y.e_m->SetPoint(idx, x, s.sh_y_e.mean());
0388   sh_y.e_v->SetPoint(idx, x, s.sh_y_e.rms());
0389   sh_y.eR->SetPoint(idx, x, s.sh_y_e.rms() / s.sh_y_u.mean());
0390
0391   sh_y.v_m->SetPointError(idx, 0., s.sh_y_v.meanError());
0392   sh_y.v_v->SetPointError(idx, 0., s.sh_y_v.rmsError());
0393   sh_y.u_m->SetPointError(idx, 0., s.sh_y_u.meanError());
0394   sh_y.u_v->SetPointError(idx, 0., s.sh_y_u.rmsError());
0395   sh_y.e_m->SetPointError(idx, 0., s.sh_y_e.meanError());
0396   sh_y.e_v->SetPointError(idx, 0., s.sh_y_e.rmsError());
0397   sh_y.eR->SetPointError(idx, 0., eR_error(s.sh_y_e, s.sh_y_u));
0398
0399   rot_z.n->SetPoint(idx, x, s.rot_z_u.n());
0400   rot_z.v_m->SetPoint(idx, x, s.rot_z_v.mean());
0401   rot_z.v_v->SetPoint(idx, x, s.rot_z_v.rms());
0402   rot_z.u_m->SetPoint(idx, x, s.rot_z_u.mean());
0403   rot_z.u_v->SetPoint(idx, x, s.rot_z_u.rms());
0404   rot_z.e_m->SetPoint(idx, x, s.rot_z_e.mean());
0405   rot_z.e_v->SetPoint(idx, x, s.rot_z_e.rms());
0406   rot_z.eR->SetPoint(idx, x, s.rot_z_e.rms() / s.rot_z_u.mean());
0407
0408   rot_z.v_m->SetPointError(idx, 0., s.rot_z_v.meanError());
0409   rot_z.v_v->SetPointError(idx, 0., s.rot_z_v.rmsError());
0410   rot_z.u_m->SetPointError(idx, 0., s.rot_z_u.meanError());
0411   rot_z.u_v->SetPointError(idx, 0., s.rot_z_u.rmsError());
0412   rot_z.e_m->SetPointError(idx, 0., s.rot_z_e.meanError());
0413   rot_z.e_v->SetPointError(idx, 0., s.rot_z_e.rmsError());
0414   rot_z.eR->SetPointError(idx, 0., eR_error(s.rot_z_e, s.rot_z_u));
0415
0416   sh_z.n->SetPoint(idx, x, s.sh_z_u.n());
0417   sh_z.v_m->SetPoint(idx, x, s.sh_z_v.mean());
0418   sh_z.v_v->SetPoint(idx, x, s.sh_z_v.rms());
0419   sh_z.u_m->SetPoint(idx, x, s.sh_z_u.mean());
0420   sh_z.u_v->SetPoint(idx, x, s.sh_z_u.rms());
0421   sh_z.e_m->SetPoint(idx, x, s.sh_z_e.mean());
0422   sh_z.e_v->SetPoint(idx, x, s.sh_z_e.rms());
0423   sh_z.eR->SetPoint(idx, x, s.sh_z_e.rms() / s.sh_z_u.mean());
0424
0425   sh_z.v_m->SetPointError(idx, 0., s.sh_z_v.meanError());
0426   sh_z.v_v->SetPointError(idx, 0., s.sh_z_v.rmsError());
0427   sh_z.u_m->SetPointError(idx, 0., s.sh_z_u.meanError());
0428   sh_z.u_v->SetPointError(idx, 0., s.sh_z_u.rmsError());
0429   sh_z.e_m->SetPointError(idx, 0., s.sh_z_e.meanError());
0430   sh_z.e_v->SetPointError(idx, 0., s.sh_z_e.rmsError());
0431   sh_z.eR->SetPointError(idx, 0., eR_error(s.sh_z_e, s.sh_z_u));
0432 }
0433
0434 //----------------------------------------------------------------------------------------------------
0435
0436 void RPGraphs::write(const char *xLabel) {
0437   gDirectory = mkdir(gDirectory, "sh_x");
0438   sh_x.write(xLabel);
0439   gDirectory->cd("..");
0440   gDirectory = mkdir(gDirectory, "sh_y");
0441   sh_y.write(xLabel);
0442   gDirectory->cd("..");
0443   gDirectory = mkdir(gDirectory, "rot_z");
0444   rot_z.write(xLabel);
0445   gDirectory->cd("..");
0446   gDirectory = mkdir(gDirectory, "sh_z");
0447   sh_z.write(xLabel);
0448   gDirectory->cd("..");
0449 }
0450
0451 //----------------------------------------------------------------------------------------------------
0452
0453 TFile *sf;  // output file
0454
0455 void writeGraphs(string r, string o, string d) {
0456   map<Desc, DetGraphs> dFcnN;   // here Desc::N is 0 by definition
0457   map<Desc, DetGraphs> dFcnIt;  // here Desc::i is 0 by definition
0458   for (map<Desc, DetStat>::iterator it = det_stat.begin(); it != det_stat.end(); ++it) {
0459     dFcnN[Desc(0, it->first.i, it->first.d)].fill(it->first.N, it->second);
0460     dFcnIt[Desc(it->first.N, 0, it->first.d)].fill(it->first.i, it->second);
0461   }
0462
0463   map<Desc, RPGraphs> rFcnN;   // here Desc::N is 0 by definition
0464   map<Desc, RPGraphs> rFcnIt;  // here Desc::i is 0 by definition
0465   for (map<Desc, RPStat>::iterator it = rp_stat.begin(); it != rp_stat.end(); ++it) {
0466     rFcnN[Desc(0, it->first.i, it->first.d)].fill(it->first.N, it->second);
0467     rFcnIt[Desc(it->first.N, 0, it->first.d)].fill(it->first.i, it->second);
0468   }
0469
0470   r = r.replace(r.find(':'), 1, ">");
0471   o = o.replace(o.find(':'), 1, ">");
0472   d = d.replace(d.find(':'), 1, ">");
0473
0474   gDirectory = sf;
0475   gDirectory = mkdir(gDirectory, r.c_str());
0476   gDirectory = mkdir(gDirectory, o.c_str());
0477   gDirectory = mkdir(gDirectory, d.c_str());
0478
0479   char buf[100];
0480   gDirectory = mkdir(gDirectory, "fcn_of_N");
0481   for (map<Desc, DetGraphs>::iterator it = dFcnN.begin(); it != dFcnN.end(); ++it) {
0482     sprintf(buf, "iteration>%u", it->first.i);
0483     gDirectory = mkdir(gDirectory, buf);
0484     sprintf(buf, "%u", it->first.d);
0485     gDirectory = mkdir(gDirectory, buf);
0486     it->second.write("tracks");
0487     gDirectory->cd("../..");
0488   }
0489
0490   for (map<Desc, RPGraphs>::iterator it = rFcnN.begin(); it != rFcnN.end(); ++it) {
0491     sprintf(buf, "iteration>%u", it->first.i);
0492     gDirectory = mkdir(gDirectory, buf);
0493     sprintf(buf, "RP %u", it->first.d);
0494     gDirectory = mkdir(gDirectory, buf);
0495     it->second.write("tracks");
0496     gDirectory->cd("../..");
0497   }
0498   gDirectory->cd("..");
0499
0500   gDirectory = mkdir(gDirectory, "fcn_of_iteration");
0501   for (map<Desc, DetGraphs>::iterator it = dFcnIt.begin(); it != dFcnIt.end(); ++it) {
0502     sprintf(buf, "N>%u", it->first.N);
0503     gDirectory = mkdir(gDirectory, buf);
0504     sprintf(buf, "%u", it->first.d);
0505     gDirectory = mkdir(gDirectory, buf);
0506     it->second.write("iteration");
0507     gDirectory->cd("../..");
0508   }
0509
0510   for (map<Desc, RPGraphs>::iterator it = rFcnIt.begin(); it != rFcnIt.end(); ++it) {
0511     sprintf(buf, "N>%u", it->first.N);
0512     gDirectory = mkdir(gDirectory, buf);
0513     sprintf(buf, "RP %u", it->first.d);
0514     gDirectory = mkdir(gDirectory, buf);
0515     it->second.write("iteration");
0516     gDirectory->cd("../..");
0517   }
0518   gDirectory->cd("..");
0519 }
0520
0521 //----------------------------------------------------------------------------------------------------
0522
0523 bool isRegDir(const dirent *de) {
0524   if (de->d_type != DT_DIR)
0525     return false;
0526
0527   if (!strcmp(de->d_name, "."))
0528     return false;
0529
0530   if (!strcmp(de->d_name, ".."))
0531     return false;
0532
0533   return true;
0534 }
0535
0536 //----------------------------------------------------------------------------------------------------
0537
0538 // Choices for RP input.
0539 // * First of all, the Jan to Ideal comparison is not possible for RP data -
0540 //   the factorizations are different, because of zero errors of the Ideal
0541 //   algorithm.
0542 // * PLAIN results do not contain RP data, do not use.
0543 // * CUMULATIVE contain the factorization from the previous iteration.
0544 // * FACTORED is probably the best choice.
0545 string r_actual_file("./cumulative_factored_results_Jan.xml");
0546 string r_ideal_file("./cumulative_factored_results_Ideal.xml");
0547
0548 // Choices for sensor input.
0549 // * PLAIN results contain only the correction by the last iteration.
0550 // * Expansion is not done, hence CUMULATIVE results contain internal ALIGNMENT
0551 //   only. Factorization inherited from the previous iteration.
0552 // * FACTORED results shall not be used - factorization is different for Jan and
0553 //   Ideal algorithms (Ideal has zero errors).
0554 // * If EXPANDED results are available, they might be used, they contain both
0555 //   internal and RP alignment. Hence comparison between different misalignments
0556 //   is not possible.
0557 string s_actual_file("./cumulative_expanded_results_Jan.xml");
0558 string s_ideal_file("./cumulative_expanded_results_Ideal.xml");
0559
0560 //----------------------------------------------------------------------------------------------------
0561
0562 CTPPSRPAlignmentCorrectionsData LoadAlignment(const string &fn) {
0563   const auto &seq = CTPPSRPAlignmentCorrectionsMethods::loadFromXML(fn);
0564   if (seq.empty())
0565     throw cms::Exception("PPS") << "LoadAlignment: alignment sequence empty, file: " << fn;
0566
0567   return seq[0].second;
0568 }
0569
0570 //----------------------------------------------------------------------------------------------------
0571
0572 void PrintHelp(const char *name) {
0573   printf("USAGE: %s r_actual_file r_ideal_file s_actual_file s_ideal_file\n", name);
0574   printf("       %s --help (to print this help)\n", name);
0575   printf("PARAMETERS:\n");
0576   printf("\tr_actual_file\t file with actual RP results (default %s)\n", r_actual_file.c_str());
0577   printf("\tr_ideal_file\t file with ideal RP results (default %s)\n", r_ideal_file.c_str());
0578   printf("\ts_actual_file\t file with actual sensor results (default %s)\n", s_actual_file.c_str());
0579   printf("\ts_ideal_file\t file with ideal sensor results (default %s)\n", s_ideal_file.c_str());
0580 }
0581
0582 //----------------------------------------------------------------------------------------------------
0583
0584 int main(int argc, const char *argv[]) {
0585   if (argc > 1) {
0586     if (!strcmp(argv[1], "-h") || !strcmp(argv[1], "--help")) {
0587       PrintHelp(argv[0]);
0588       return 0;
0589     }
0590
0591     if (argc > 1)
0592       r_actual_file = argv[1];
0593     if (argc > 2)
0594       r_ideal_file = argv[2];
0595     if (argc > 3)
0596       s_actual_file = argv[3];
0597     if (argc > 4)
0598       s_ideal_file = argv[4];
0599   }
0600
0601   printf("r_actual_file: %s\n", r_actual_file.c_str());
0602   printf("r_ideal_file: %s\n", r_ideal_file.c_str());
0603   printf("s_actual_file: %s\n", s_actual_file.c_str());
0604   printf("s_ideal_file: %s\n", s_ideal_file.c_str());
0605
0606   // open output file
0607   sf = new TFile("result_summary.root", "recreate");
0608
0609   // traverse directory structure
0610   try {
0611     DIR *dp_r = opendir(".");
0612     dirent *de_r;
0613     // traverse RPs directories
0614     while ((de_r = readdir(dp_r))) {
0615       if (!isRegDir(de_r))
0616         continue;
0617
0618       chdir(de_r->d_name);
0619       DIR *dp_o = opendir(".");
0620       // traverse optimized directories
0621       dirent *de_o;
0622       while ((de_o = readdir(dp_o))) {
0623         if (!isRegDir(de_o))
0624           continue;
0625
0626         chdir(de_o->d_name);
0627         DIR *dp_d = opendir(".");
0628         dirent *de_d;
0629         // traverse tr_dist directories
0630         while ((de_d = readdir(dp_d))) {
0631           if (!isRegDir(de_d))
0632             continue;
0633
0634           // sort the tr_N directories by N
0635           chdir(de_d->d_name);
0636           DIR *dp_N = opendir(".");
0637           dirent *de_N;
0638           map<unsigned int, string> nMap;
0639           while ((de_N = readdir(dp_N))) {
0640             if (isRegDir(de_N)) {
0641               string sN = de_N->d_name;
0642               sN = sN.substr(sN.find(':') + 1);
0643               unsigned int N = atof(sN.c_str());
0644               nMap[N] = de_N->d_name;
0645             }
0646           }
0647
0648           resetStatistics();
0649
0650           // traverse tr_N directories
0651           for (map<unsigned int, string>::iterator nit = nMap.begin(); nit != nMap.end(); ++nit) {
0652             chdir(nit->second.c_str());
0653             DIR *dp_i = opendir(".");
0654             dirent *de_i;
0655             // traverse repetitions directories
0656             while ((de_i = readdir(dp_i))) {
0657               if (!isRegDir(de_i))
0658                 continue;
0659
0660               chdir(de_i->d_name);
0661
0662               // sort the iteration directories
0663               DIR *dp_it = opendir(".");
0664               dirent *de_it;
0665               map<unsigned int, string> itMap;
0666               while ((de_it = readdir(dp_it))) {
0667                 if (isRegDir(de_it)) {
0668                   unsigned idx = 9;
0669                   if (de_it->d_name[9] == ':')
0670                     idx = 10;
0671                   unsigned int it = atoi(de_it->d_name + idx);
0672                   itMap[it] = de_it->d_name;
0673                 }
0674               }
0675
0676               // traverse iteration directories
0677               for (map<unsigned int, string>::iterator iit = itMap.begin(); iit != itMap.end(); ++iit) {
0678                 chdir(iit->second.c_str());
0679
0680                 printf("%s|%s|%s|%s|%s|%s\n",
0681                        de_r->d_name,
0682                        de_o->d_name,
0683                        de_d->d_name,
0684                        nit->second.c_str(),
0685                        de_i->d_name,
0686                        iit->second.c_str());
0687
0688                 // load and process alignments
0689                 try {
0690                   const auto &r_actual = LoadAlignment(r_actual_file);
0691                   const auto &r_ideal = LoadAlignment(r_ideal_file);
0692                   updateRPStatistics(nit->first, iit->first, r_actual, r_ideal);
0693
0694                   const auto &s_actual = LoadAlignment(s_actual_file);
0695                   const auto &s_ideal = LoadAlignment(s_ideal_file);
0696                   updateSensorStatistics(nit->first, iit->first, s_actual, s_ideal);
0697                 } catch (cms::Exception &e) {
0698                   printf("ERROR: A CMS exception has been caught:\n%s\nSkipping this directory.\n", e.what());
0699                 }
0700
0701                 chdir("..");
0702               }
0703
0704               closedir(dp_it);
0705               chdir("..");
0706             }
0707
0708             closedir(dp_i);
0709             chdir("..");
0710           }
0711
0712           writeGraphs(de_r->d_name, de_o->d_name, de_d->d_name);
0713
0714           closedir(dp_N);
0715           chdir("..");
0716         }
0717
0718         closedir(dp_d);
0719         chdir("..");
0720       }
0721
0722       closedir(dp_o);
0723       chdir("..");
0724     }
0725
0726     closedir(dp_r);
0727     chdir("..");
0728
0729     delete sf;
0730   }
0731
0732   catch (cms::Exception &e) {
0733     printf("ERROR: A CMS exception has been caught:\n%s\nStopping.\n", e.what());
0734   }
0735
0736   catch (std::exception &e) {
0737     printf("ERROR: A std::exception has been caught:\n%s\nStopping.\n", e.what());
0738   }
0739
0740   catch (...) {
0741     printf("ERROR: An exception has been caught, stopping.\n");
0742   }
0743
0744   return 0;
0745 }