File indexing completed on 2024-04-06 11:57:23
0001 #include <algorithm>
0002 #include <cmath>
0003 #include <iomanip>
0004 #include <iostream>
0005 #include <map>
0006 #include <sstream>
0007 #include <string>
0008 #include <vector>
0009
0010 #include "TCanvas.h"
0011 #include "TF1.h"
0012 #include "TFile.h"
0013 #include "TGaxis.h"
0014 #include "TGraph.h"
0015 #include "TH1F.h"
0016 #include "TLegend.h"
0017 #include "TList.h"
0018 #include "TMatrixD.h"
0019 #include "TNtuple.h"
0020 #include "TStyle.h"
0021 #include "TText.h"
0022 #include "TTree.h"
0023 #include "TVectorD.h"
0024
0025 #include "NumbersAndNames.C"
0026
0027 #if !defined(__CINT__) && !defined(__MAKECINT__)
0028
0029 #include "Alignment/CommonAlignment/interface/AlignableObjectId.h"
0030 #endif
0031
0032 const int nPar = 6;
0033
0034 void setBranch(TTree* t, double p[], double r[])
0035 {
0036 t->SetBranchStatus("*", 0);
0037 t->SetBranchStatus("Pos");
0038 t->SetBranchStatus("Rot");
0039
0040 t->SetBranchAddress("Pos", p);
0041 t->SetBranchAddress("Rot", r);
0042 }
0043
0044 void writeShifts(std::string path)
0045 {
0046 TFile ft((path + "IOTruePositions.root" ).c_str());
0047 TFile f0((path + "IOMisalignedPositions.root").c_str());
0048 TFile fp((path + "IOAlignedPositions.root" ).c_str());
0049
0050 TTree* tt = (TTree*)ft.Get("AlignablesOrgPos_1");
0051
0052 const TList* keys = fp.GetListOfKeys();
0053
0054 const std::string key0Name = keys->At(0)->GetName();
0055
0056 const bool iter0Exist = ('0' == *(key0Name.end() - 1));
0057
0058 const int nIteration = iter0Exist ? keys->GetSize() - 1 : keys->GetSize();
0059 const int nAlignable = tt->GetEntries();
0060
0061 std::vector<TVectorD> post(nAlignable, TVectorD(3));
0062 std::vector<TMatrixD> rott(nAlignable, TMatrixD(3, 3));
0063
0064 double pos[3];
0065 double rot[9];
0066
0067 setBranch(tt, pos, rot);
0068
0069 for (int n = 0; n < nAlignable; ++n)
0070 {
0071 tt->GetEntry(n);
0072
0073 for (int p = 0; p < 3; ++p) post[n][p] = pos[p];
0074 for (int p = 0; p < 9; ++p) rott[n](p / 3, p % 3) = rot[p];
0075 }
0076
0077 std::vector<TTree*> trees(nIteration + 1);
0078
0079 trees[0] = (TTree*)f0.Get("AlignablesAbsPos_1");
0080
0081 for (int i = 1; i <= nIteration; ++i)
0082 {
0083 trees[i] = (TTree*)fp.Get(keys->At(iter0Exist ? i : i - 1)->GetName());
0084 }
0085
0086 TFile fout((path + "shifts.root").c_str(), "RECREATE");
0087
0088 std::vector<TNtuple*> tuples(nIteration + 1);
0089
0090 for (int i = 0; i <= nIteration; ++i)
0091 {
0092
0093
0094
0095
0096
0097
0098
0099 std::ostringstream o; o << 't' << i;
0100
0101 tuples[i] = new TNtuple(o.str().c_str(), "", "u:v:w:a:b:g");
0102
0103 setBranch(trees[i], pos, rot);
0104
0105 for (int n = 0; n < nAlignable; ++n)
0106 {
0107 trees[i]->GetEntry(n);
0108
0109 TVectorD dr(3, pos);
0110 TMatrixD dR(3, 3, rot);
0111
0112
0113
0114
0115 dr -= post[n];
0116 dr *= 1e4;
0117 dr = rott[n] * dr;
0118 dR = dR * TMatrixD(TMatrixD::kTransposed, rott[n]);
0119
0120 tuples[i]->Fill(dr[0], dr[1], dr[2],
0121 1e3 * -std::atan2(dR(2, 1), dR(2, 2)),
0122 1e3 * std::asin(dR(2, 0)),
0123 1e3 * -std::atan2(dR(1, 0), dR(0, 0)));
0124 }
0125 }
0126
0127 fout.Write();
0128
0129 for (int i = 0; i <= nIteration; ++i) delete tuples[i];
0130 }
0131
0132 void writePars(std::string path)
0133 {
0134 TFile ft((path + "IOTruePositions.root" ).c_str());
0135 TFile f0((path + "IOMisalignedPositions.root").c_str());
0136 TFile fp((path + "IOAlignmentParameters.root").c_str());
0137
0138 TTree* tt = (TTree*)ft.Get("AlignablesOrgPos_1");
0139 TTree* t0 = (TTree*)f0.Get("AlignablesAbsPos_1");
0140
0141 const TList* keys = fp.GetListOfKeys();
0142
0143 const std::string key0Name = keys->At(0)->GetName();
0144
0145 const bool iter0Exist = ('0' == *(key0Name.end() - 1));
0146
0147 const int nIteration = iter0Exist ? keys->GetSize() - 1 : keys->GetSize();
0148 const int nAlignable = tt->GetEntries();
0149
0150 std::vector<TVectorD> post(nAlignable, TVectorD(3));
0151 std::vector<TMatrixD> rott(nAlignable, TMatrixD(3, 3));
0152
0153 double pos[3];
0154 double rot[9];
0155 double par[nPar];
0156
0157 setBranch(tt, pos, rot);
0158
0159 for (int n = 0; n < nAlignable; ++n)
0160 {
0161 tt->GetEntry(n);
0162
0163 for (int p = 0; p < 3; ++p) post[n][p] = pos[p];
0164 for (int p = 0; p < 9; ++p) rott[n](p / 3, p % 3) = rot[p];
0165 }
0166
0167 TFile fout((path + "pars.root").c_str(), "RECREATE");
0168
0169 std::vector<TNtuple*> tuples(nIteration + 1);
0170
0171 tuples[0] = new TNtuple("t0", "", "u:v:w:a:b:g");
0172
0173 setBranch(t0, pos, rot);
0174
0175 for (int n = 0; n < nAlignable; ++n)
0176 {
0177 t0->GetEntry(n);
0178
0179 TVectorD dr(3, pos);
0180 TMatrixD dR(3, 3, rot);
0181
0182 dr -= post[n];
0183 dr *= 1e4;
0184 dr = rott[n] * dr;
0185 dR = dR * TMatrixD(TMatrixD::kTransposed, rott[n]);
0186
0187 tuples[0]->Fill(dr[0], dr[1], dr[2],
0188 1e3 * -std::atan2(dR(2, 1), dR(2, 2)),
0189 1e3 * std::asin(dR(2, 0)),
0190 1e3 * -std::atan2(dR(1, 0), dR(0, 0)));
0191 }
0192
0193 for (int i = 1; i <= nIteration; ++i)
0194 {
0195 TTree* tp = (TTree*)fp.Get(keys->At(iter0Exist ? i : i - 1)->GetName());
0196
0197 tp->SetBranchStatus("*", 0);
0198 tp->SetBranchStatus("Par");
0199 tp->SetBranchAddress("Par", par);
0200
0201 std::ostringstream o; o << 't' << i;
0202
0203 tuples[i] = new TNtuple(o.str().c_str(), "", "u:v:w:a:b:g");
0204
0205 const int nAlignable = tp->GetEntries();
0206
0207 for (int n = 0; n < nAlignable; ++n)
0208 {
0209 tp->GetEntry(n);
0210
0211 tuples[i]->Fill(par[0] * 1e4, par[1] * 1e4, par[2] * 1e4,
0212 par[3] * 1e3, par[4] * 1e3, par[5] * 1e3);
0213 }
0214 }
0215
0216 fout.Write();
0217
0218 for (int i = 0; i <= nIteration; ++i) delete tuples[i];
0219 }
0220
0221 void writeSurvey(std::string path, std::string type)
0222 {
0223 static AlignableObjectId objId;
0224
0225 TFile f0((path + "histograms.root").c_str());
0226
0227 const TList* keys = f0.GetListOfKeys();
0228
0229 const int nIteration = keys->GetSize();
0230
0231 std::ostringstream o; o << path << "survey_" << type << ".root";
0232
0233 TFile fout(o.str().c_str(), "RECREATE");
0234
0235 std::vector<TNtuple*> tuples(nIteration);
0236
0237 int level;
0238 double par[nPar];
0239
0240 for (int i = 0; i < nIteration; ++i)
0241 {
0242 std::string name = keys->At(i)->GetName();
0243 name += "/survey";
0244
0245 TTree* t0 = (TTree*)f0.Get(name.c_str());
0246
0247 t0->SetBranchAddress("level", &level);
0248 t0->SetBranchAddress("par", par);
0249
0250 std::ostringstream o; o << 't' << i + 1;
0251
0252 tuples[i] = new TNtuple(o.str().c_str(), "", "u:v:w:a:b:g");
0253
0254 const int nEntry = t0->GetEntries();
0255
0256 for (int n = 0; n < nEntry; ++n)
0257 {
0258 t0->GetEntry(n);
0259
0260 if (objId.nameToType(type) == level)
0261 tuples[i]->Fill(par[0] * 1e4, par[1] * 1e4, par[2] * 1e4,
0262 par[3] * 1e3, par[4] * 1e3, par[5] * 1e3);
0263 }
0264 }
0265
0266 fout.Write();
0267
0268 for (int i = 0; i < nIteration; ++i) delete tuples[i];
0269 }
0270
0271 class AlignPlots
0272 {
0273 typedef std::vector<float> AlignSet;
0274
0275 static const char* const titles_[nPar];
0276
0277 public:
0278
0279 AlignPlots(std::string file);
0280 AlignPlots(std::string file, std::vector<unsigned int> levels, int minHit = 0);
0281
0282 void iters() const;
0283
0284 void iter(int iter) const;
0285
0286 void dump(int index, int iterN = 0) const;
0287
0288 private:
0289
0290 static float sum(const AlignSet&);
0291 static float sum2(const AlignSet&);
0292 static float width(const AlignSet&);
0293
0294 std::string file_;
0295
0296 std::vector<AlignSet> alignSets_[nPar];
0297 };
0298
0299 const char* const AlignPlots::titles_[nPar] =
0300 {"#Deltau (#mum)", "#Deltav (#mum)", "#Deltaw (#mum)",
0301 "#Delta#alpha (mrad)", "#Delta#beta (mrad)", "#Delta#gamma (mrad)"};
0302
0303 AlignPlots::AlignPlots(std::string file):
0304 file_(file)
0305 {
0306 TFile fin((file_ + ".root").c_str());
0307
0308 const TList* keys = fin.GetListOfKeys();
0309
0310 TNtuple* t = (TNtuple*)fin.Get(keys->At(0)->GetName());
0311
0312 const int nIteration = keys->GetSize();
0313 const int nAlignable = t->GetEntries();
0314
0315 for (int p = 0; p < nPar; ++p)
0316 alignSets_[p].resize(nIteration, AlignSet(nAlignable));
0317
0318 for (int i = 0; i < nIteration; ++i)
0319 {
0320 t = (TNtuple*)fin.Get(keys->At(i)->GetName());
0321
0322 for (int n = 0; n < nAlignable; ++n)
0323 {
0324 t->GetEntry(n);
0325
0326 const float* pars = t->GetArgs();
0327
0328 for (int p = 0; p < nPar; ++p) alignSets_[p][i][n] = pars[p];
0329 }
0330 }
0331 }
0332
0333 AlignPlots::AlignPlots(std::string file, std::vector<unsigned int> levels, int minHit):
0334 file_(file)
0335 {
0336 const unsigned int nLevel = levels.size();
0337
0338 if (nLevel == 0)
0339 {
0340 std::cout << "No input levels." << std::endl; return;
0341 }
0342
0343 const CounterNames& cn = Helper::counterNames(levels[0]);
0344
0345 if (cn.size() < nLevel)
0346 {
0347 std::cout << "Too many input levels for " << cn[0].second
0348 << ". Max " << cn.size()
0349 << std::endl;
0350 return;
0351 }
0352
0353 std::string path = file.substr(0, file.find_last_of('/'));
0354
0355 TFile fu((path + "/IOUserVariables.root").c_str());
0356
0357 TTree* tu = (TTree*)fu.Get("T9_1");
0358
0359 unsigned int id(0);
0360 int nHit(0);
0361
0362 tu->SetBranchStatus("*", 0);
0363 tu->SetBranchStatus("Id");
0364 tu->SetBranchStatus("Nhit");
0365 tu->SetBranchAddress("Id" , &id);
0366 tu->SetBranchAddress("Nhit", &nHit);
0367
0368 TFile fin((file + ".root").c_str());
0369
0370 const TList* keys = fin.GetListOfKeys();
0371
0372 TNtuple* t = (TNtuple*)fin.Get(keys->At(0)->GetName());
0373
0374 const int nIteration = keys->GetSize();
0375 const int nAlignable = t->GetEntries();
0376
0377 for (int p = 0; p < nPar; ++p) alignSets_[p].resize(nIteration);
0378
0379 for (int i = 0; i < nIteration; ++i)
0380 {
0381 for (int p = 0; p < nPar; ++p) alignSets_[p][i].reserve(nAlignable);
0382
0383 t = (TNtuple*)fin.Get(keys->At(i)->GetName());
0384
0385 for (int n = 0; n < nAlignable; ++n)
0386 {
0387 t->GetEntry(n);
0388 tu->GetEntry(n);
0389
0390 const float* pars = t->GetArgs();
0391
0392 bool selected = (nHit >= minHit);
0393
0394 for (unsigned int l = 0; selected && l < nLevel; ++l)
0395 selected = (cn[l].first(id) == levels[l]);
0396
0397 if (selected)
0398 for (int p = 0; p < nPar; ++p)
0399 alignSets_[p][i].push_back(pars[p]);
0400 }
0401 }
0402
0403 std::ostringstream o;
0404
0405 for (unsigned int l = 0; l < nLevel; ++l)
0406 o << '_' << cn[l].second << levels[l];
0407
0408 o << "_minHit" << minHit;
0409 file_ += o.str();
0410 }
0411
0412 void AlignPlots::iters() const
0413 {
0414 gStyle->SetOptTitle(0);
0415 gStyle->SetOptStat(0);
0416
0417 TGaxis::SetMaxDigits(3);
0418
0419 const int nIteration = alignSets_[0].size();
0420 const int nAlignable = alignSets_[0][0].size();
0421
0422 if (0 == nAlignable)
0423 {
0424 std::cout << "0 Alignables selected." << std::endl; return;
0425 }
0426
0427 TCanvas c("c", "c", 1200, 800);
0428
0429 c.Divide(3, 2);
0430
0431 std::vector<TGraph> graphs[nPar];
0432
0433 for (int p = 0; p < nPar; ++p)
0434 {
0435 c.cd(p + 1);
0436
0437
0438
0439 std::vector<float> ylimits(nIteration);
0440
0441 for (int i = 0; i < nIteration; ++i)
0442 {
0443 const AlignSet& aSet = alignSets_[p][i];
0444
0445 float ymin = *std::min_element(aSet.begin(), aSet.end());
0446 float ymax = *std::max_element(aSet.begin(), aSet.end());
0447
0448 ylimits[i] = std::max(std::abs(ymin), std::abs(ymax));
0449
0450 }
0451
0452 float ylimit = *std::max_element(ylimits.begin(), ylimits.end());
0453
0454 graphs[p].resize(nAlignable, nIteration);
0455
0456 TGraph& g = graphs[p][0];
0457
0458 g.SetMinimum(-ylimit);
0459 g.SetMaximum( ylimit);
0460 g.GetXaxis()->SetLimits(0., nIteration - 1.);
0461 g.GetXaxis()->SetTitle("iteration");
0462 g.GetYaxis()->SetTitle(titles_[p]);
0463 g.GetYaxis()->SetTitleSize(.04);
0464 g.Draw("AP");
0465
0466 for (int n = 0; n < nAlignable; ++n)
0467 {
0468 TGraph& g = graphs[p][n];
0469
0470 for (int i = 0; i < nIteration; ++i)
0471 {
0472 float y = alignSets_[p][i][n];
0473
0474 g.SetPoint(i, i, isnan(y) ? 0.f : y);
0475 }
0476
0477 g.Draw("L");
0478 }
0479 }
0480
0481 std::ostringstream o;
0482 o << file_ << "_vs_iter";
0483
0484 c.SaveAs((o.str() + ".png").c_str());
0485 c.SaveAs((o.str() + ".eps").c_str());
0486 }
0487
0488 void AlignPlots::dump(int index, int iterN) const
0489 {
0490 gStyle->SetOptTitle(0);
0491 gStyle->SetOptStat(0);
0492
0493 TGaxis::SetMaxDigits(3);
0494
0495 const int nIteration = alignSets_[0].size();
0496 const int nAlignable = alignSets_[0][0].size();
0497
0498 if (0 >= iterN || nIteration <= iterN) iterN = nIteration;
0499 else ++iterN;
0500
0501 if (index >= nAlignable)
0502 {
0503 std::cout << "Alignable index too big. "
0504 << "Number of Alignables is " << nAlignable
0505 << std::endl;
0506 return;
0507 }
0508
0509 TCanvas c("c", "c", 1200, 800);
0510
0511 c.Divide(3, 2);
0512
0513 TGraph graphs[nPar];
0514
0515 for (int p = 0; p < nPar; ++p)
0516 {
0517 c.cd(p + 1);
0518
0519
0520
0521 std::vector<float> ylimits(iterN);
0522
0523 for (int i = 0; i < iterN; ++i)
0524 {
0525 ylimits[i] = std::abs(alignSets_[p][i][index]);
0526 }
0527
0528 float ylimit = *std::max_element(ylimits.begin(), ylimits.end());
0529
0530 TGraph& g = graphs[p];
0531
0532 g.Set(iterN);
0533 g.SetMinimum(-ylimit);
0534 g.SetMaximum( ylimit);
0535 g.GetXaxis()->SetLimits(0., iterN - 1.);
0536 g.GetXaxis()->SetTitle("iteration");
0537 g.GetYaxis()->SetTitle(titles_[p]);
0538 g.GetYaxis()->SetTitleSize(.04);
0539 g.Draw("AP");
0540
0541 for (int i = 0; i < iterN; ++i)
0542 {
0543 g.SetPoint(i, i, alignSets_[p][i][index]);
0544 }
0545
0546 g.Draw("L");
0547 }
0548
0549 std::ostringstream o;
0550 o << file_ << "_vs_iter_Alignable" << index;
0551
0552 c.SaveAs((o.str() + ".png").c_str());
0553 c.SaveAs((o.str() + ".eps").c_str());
0554 }
0555
0556 float AlignPlots::sum(const AlignSet& aSet)
0557 {
0558 float sum(0.);
0559
0560 for (unsigned int i = 0; i < aSet.size(); ++i)
0561 if (std::abs(aSet[i]) < 1e8) sum += aSet[i];
0562
0563 return sum;
0564 }
0565
0566 float AlignPlots::sum2(const AlignSet& aSet)
0567 {
0568 float sum(0.);
0569
0570 for (unsigned int i = 0; i < aSet.size(); ++i)
0571 if (std::abs(aSet[i]) < 1e8) sum += aSet[i] * aSet[i];
0572
0573 return sum;
0574 }
0575
0576 float AlignPlots::width(const AlignSet& aSet)
0577 {
0578 float invN = 1. / aSet.size();
0579 float mean = sum(aSet) * invN;
0580 float rms2 = sum2(aSet) * invN - mean * mean;
0581
0582 return rms2 > 0. ? 3. * std::sqrt(rms2) : 1e-3;
0583 }
0584
0585 void AlignPlots::iter(int iter) const
0586 {
0587 gStyle->SetOptTitle(0);
0588 gStyle->SetOptStat(10);
0589 gStyle->SetOptFit(0);
0590 gStyle->SetStatH(.2);
0591 gStyle->SetStatW(.3);
0592
0593 const int nIteration = alignSets_[0].size();
0594 const int nAlignable = alignSets_[0][0].size();
0595
0596 if (nIteration <= iter)
0597 {
0598 std::cout << "Iteration number too big. "
0599 << "Setting to max " << nIteration - 1
0600 << std::endl;
0601
0602 iter = nIteration - 1;
0603 }
0604
0605 if (0 == nAlignable)
0606 {
0607 std::cout << "0 Alignables selected." << std::endl; return;
0608 }
0609
0610 TCanvas c("c", "c", 1200, 800);
0611
0612 c.Divide(3, 2);
0613
0614 TH1F hists[nPar];
0615
0616 for (int p = 0; p < nPar; ++p)
0617 {
0618 const AlignSet& setI = alignSets_[p][iter];
0619
0620 float mean = sum(setI) / setI.size();
0621 float rms3 = std::min(width(alignSets_[p][0]), width(setI));
0622
0623 TH1F& h = hists[p];
0624
0625 h.SetBins(50, mean - rms3, mean + rms3);
0626
0627 for (int n = 0; n < nAlignable; ++n) h.Fill(setI[n]);
0628
0629 c.cd(p + 1);
0630 h.SetXTitle(titles_[p]);
0631 h.SetTitleSize(.04);
0632 h.Fit("gaus", "LQ");
0633
0634 TF1* f = h.GetFunction("gaus");
0635
0636 f->SetLineWidth(1);
0637 f->SetLineColor(kRed);
0638
0639 std::ostringstream o;
0640
0641 o << std::fixed << std::setprecision(2)
0642 << "width = " << f->GetParameter(2);
0643
0644 TText width; width.DrawTextNDC(.6, .8, o.str().c_str());
0645 }
0646
0647 std::ostringstream o;
0648 o << file_ << "_iter" << iter;
0649
0650 c.SaveAs((o.str() + ".png").c_str());
0651 c.SaveAs((o.str() + ".eps").c_str());
0652 }
0653
0654 void compareShifts(std::string tree)
0655 {
0656 gStyle->SetOptStat(0);
0657
0658 TFile fm("merged/shifts.root");
0659 TFile ft("tracks/shifts.root");
0660 TFile fs("survey/shifts.root");
0661
0662 TTree* tm = (TTree*)fm.Get(tree.c_str());
0663 TTree* tt = (TTree*)ft.Get(tree.c_str());
0664 TTree* ts = (TTree*)fs.Get(tree.c_str());
0665
0666 tm->SetMarkerStyle(kCircle);
0667 tt->SetMarkerStyle(kPlus);
0668
0669
0670
0671
0672
0673
0674
0675
0676 TCanvas c("c", "c", 1200, 800);
0677
0678 c.Divide(3, 2);
0679
0680 const char* const vars[6] = {"x", "y", "z", "a", "b", "g"};
0681
0682 for (int i = 0; i < 6; ++i)
0683 {
0684 c.cd(i + 1);
0685 ts->Draw(vars[i]);
0686 tt->Draw(vars[i], "", "sameP");
0687 tm->Draw(vars[i], "", "sameP");
0688 }
0689
0690 TLegend leg(.7, .7, .9, .9);
0691
0692 leg.AddEntry(ts, "survey");
0693 leg.AddEntry(tt, "tracks");
0694 leg.AddEntry(tm, "merged");
0695 leg.Draw();
0696
0697 c.SaveAs(("compareShifts_" + tree + ".png").c_str());
0698 }