MillePedeAlignmentAlgorithm/macros/CompareMillePede.C

0001 // Original Author: Gero Flucke
0002 // last change    : $Date: 2013/03/07 11:27:09 $
0003 // by             : $Author: flucke $
0004
0005 #include "CompareMillePede.h"
0006 #include "PlotMillePede.h"
0007 #include "GFUtils/GFHistManager.h"
0008
0009 #include <TTree.h>
0010 #include <TString.h>
0011 #include <TError.h>
0012 #include <TH1.h>
0013 #include <TH2.h>
0014 #include <TLine.h>
0015
0016 #include <iostream>
0017 const unsigned int CompareMillePede::kNpar = MillePedeTrees::kNpar; // number of parameters we have...
0018
0019 //_________________________________________________________________________________________________
0020 CompareMillePede::CompareMillePede(const char *fileName1, const char *fileName2,
0021                                    Int_t iov1, Int_t iov2, Int_t hieraLevel) :
0022   fPlotMp1(new PlotMillePede(fileName1, iov1, hieraLevel, "first")),
0023   fPlotMp2(new PlotMillePede(fileName2, iov2, hieraLevel, "second")),
0024   fHistManager(new GFHistManager)
0025 {
0026   fHistManager->SetLegendX1Y1X2Y2(0.14, 0.7, 0.45, 0.9);
0027
0028   TTree *tree1 = fPlotMp1->GetMainTree();
0029   TTree *tree2 = fPlotMp2->GetMainTree();
0030   if (!tree1 || !tree2) {
0031     ::Error("CompareMillePede", "Stop here: Previous problems...");
0032   } else if (tree1->GetEntries() != tree2->GetEntries()) {
0033     // '%lld' for Long64_t
0034     ::Error("CompareMillePede", "Stop here: %lld alignables in tree1, %lld in tree2.",
0035         tree1->GetEntries(), tree2->GetEntries());
0036   } else {
0037     tree1->AddFriend(tree2);
0038     if (!this->IsConsistent()) {
0039       ::Error("CompareMillePede", "Alignables in inconsistent order.");
0040     }
0041   }
0042 }
0043
0044 //_________________________________________________________________________________________________
0045 CompareMillePede::~CompareMillePede()
0046 {
0047   delete fHistManager;
0048   delete fPlotMp1;
0049   delete fPlotMp2;
0050 }
0051
0052 //_________________________________________________________________________________________________
0053 void CompareMillePede::DrawPedeParam(Option_t *option, unsigned int nNonRigidParam)
0054 {
0055
0056   const TString opt(option);
0057
0058   const Int_t layer = this->PrepareAdd(opt.Contains("add", TString::kIgnoreCase));
0059   const bool lineVs = opt.Contains("line", TString::kIgnoreCase);
0060
0061   const TString titleAdd = this->TitleAdd();
0062
0063   const PlotMillePede *m = fPlotMp1;
0064   UInt_t nPlot = 0;
0065   for (UInt_t iPar = 0; iPar < kNpar + nNonRigidParam; ++iPar) { //
0066     TString sel("");
0067     this->AddIsFreeSel(sel, opt, iPar);
0068     fPlotMp1->AddBasicSelection(sel);
0069     fPlotMp2->AddBasicSelection(sel);
0070
0071     const TString pedePar1(fPlotMp1->MpT() += m->Par(iPar));
0072     const TString pedePar2(fPlotMp2->MpT() += m->Par(iPar));
0073     const TString deltaPedePar(m->Parenth(pedePar2 + m->Min() += pedePar1)
0074                                += m->ToMumMuRadPede(iPar));
0075     const TString deltaName(m->Unique(Form("deltaPedePar%d", iPar)));
0076     TH1 *h = fPlotMp1->CreateHist(deltaPedePar, sel, deltaName);
0077     if (0. == h->GetEntries()) {
0078       delete h;
0079       continue;
0080     }
0081
0082     const TString name2D(m->Unique(Form("pedePar2D%d", iPar)));
0083     TH2 *hVs = fPlotMp1->CreateHist2D(pedePar1 + m->ToMumMuRadPede(iPar),
0084                       pedePar2 + m->ToMumMuRadPede(iPar), sel, name2D, "BOX");
0085
0086     const TString diff(Form("%s_{2}-%s_{1}", m->NamePede(iPar).Data(), m->NamePede(iPar).Data()));
0087     h->SetTitle(diff + titleAdd + ";" + diff + m->UnitPede(iPar) +=";#parameters");
0088
0089     hVs->SetTitle(m->NamePede(iPar).Data() + titleAdd
0090           + Form(";%s_{1}", m->NamePede(iPar).Data()) + m->UnitPede(iPar)
0091           += Form(";%s_{2}", m->NamePede(iPar).Data()) + m->UnitPede(iPar));
0092     fHistManager->AddHist(h, layer);
0093     fHistManager->AddHist(hVs, layer+1);
0094
0095     if (lineVs) { // add a TLine
0096       const Double_t xMin = hVs->GetXaxis()->GetXmin();
0097       const Double_t xMax = hVs->GetXaxis()->GetXmax();
0098       TLine *line = new TLine(xMin, xMin, xMax, xMax);
0099       line->SetLineColor(kRed); line->SetLineWidth(2);
0100       fHistManager->AddObject(line, layer + 1, iPar);
0101     }
0102
0103     ++nPlot;
0104   }
0105
0106   fHistManager->Draw();
0107 }
0108
0109 //_________________________________________________________________________________________________
0110 void CompareMillePede::DrawPedeParamVsLocation(Option_t *option, unsigned int nNonRigidParam)
0111 {
0112
0113   ::Warning("CompareMillePede::DrawPedeParamVsLocation",
0114         "Not yet fully adjusted to possible bows");
0115   const TString opt(option);
0116   const Int_t layer = this->PrepareAdd(opt.Contains("add", TString::kIgnoreCase));
0117   const TString titleAdd = this->TitleAdd();
0118
0119   const PlotMillePede *m = fPlotMp1;
0120   UInt_t nPlot = 0;
0121   for (UInt_t iPar = 0; iPar < kNpar + nNonRigidParam; ++iPar) { //
0122     TString sel("");
0123     this->AddIsFreeSel(sel, opt, iPar);
0124     fPlotMp1->AddBasicSelection(sel);
0125     fPlotMp2->AddBasicSelection(sel);
0126
0127     const TString deltaPedePar(m->Parenth(fPlotMp2->MpT() += m->Par(iPar) += m->Min()
0128                                           += fPlotMp1->MpT() += m->Par(iPar))
0129                                += m->ToMumMuRadPede(iPar));
0130
0131     const TString deltaNameR(fPlotMp1->Unique(Form("deltaPedePar%dR", iPar))); //
0132     TH2 *hDr = fPlotMp1->CreateHist2D(fPlotMp1->RPos(fPlotMp1->OrgPosT()), deltaPedePar, sel,
0133                                       deltaNameR, "BOX");
0134     if (0. == hDr->GetEntries()) {
0135       delete hDr;
0136       continue;
0137     }
0138
0139     const TString deltaNamePhi(fPlotMp1->Unique(Form("deltaPedePar%dPhi", iPar))); //
0140     TH2 *hDphi = fPlotMp1->CreateHist2D(fPlotMp1->Phi(fPlotMp1->OrgPosT()), deltaPedePar, sel,
0141                                         deltaNamePhi, "BOX");
0142     const TString deltaNameZ(fPlotMp1->Unique(Form("deltaPedePar%dZ", iPar))); //
0143     TH2 *hDz = fPlotMp1->CreateHist2D(fPlotMp1->OrgPosT() += fPlotMp1->ZPos(), deltaPedePar, sel,
0144                                       deltaNameZ, "BOX");
0145
0146     const TString deltaNameH(fPlotMp1->Unique(Form("deltaPedePar%dH", iPar))); //
0147     TH2 *hDhit = fPlotMp1->CreateHist2D(fPlotMp1->HitsX(), deltaPedePar, sel, deltaNameH, "BOX");
0148
0149 //     // now devided by sigma prediction from 1
0150 //     const TString deltaBySi1(this->DeltaParBySigma(iPar, fPlotMp1));
0151 //     TString sel1Ok(fPlotMp1->ParSiOk(iPar));
0152 //     if (!sel.IsNull()) sel1Ok.Append(fPlotMp1->AndL() += sel);
0153 //     const TString deltaBySiName1R(fPlotMp1->Unique(Form("deltaPar%dRRel1", iPar))); //
0154 //     TH2 *hDrS1 = fPlotMp1->CreateHist2D(fPlotMp1->RPos(fPlotMp1->OrgPosT()),
0155 //                                         deltaBySi1, sel1Ok, deltaBySiName1R, "BOX");
0156 //     const TString deltaBySiName1Phi(fPlotMp1->Unique(Form("deltaPar%dPhiRel1", iPar))); //
0157 //     TH2 *hDphiS1 = fPlotMp1->CreateHist2D(fPlotMp1->Phi(fPlotMp1->OrgPosT()),
0158 //                                           deltaBySi1, sel1Ok, deltaBySiName1Phi, "BOX");
0159 //     const TString deltaBySiName1Z(fPlotMp1->Unique(Form("deltaPar%dZRel1", iPar))); //
0160 //     TH2 *hDzS1 = fPlotMp1->CreateHist2D(fPlotMp1->OrgPosT() += fPlotMp1->ZPos(),
0161 //                                         deltaBySi1, sel1Ok, deltaBySiName1Z, "BOX");
0162
0163 //     // same now devided by sigma prediction from 2 (but CreateHist from 1!)
0164 //     const TString deltaBySi2(this->DeltaParBySigma(iPar, fPlotMp2));
0165 //     TString sel2Ok(fPlotMp2->ParSiOk(iPar));
0166 //     if (!sel.IsNull()) sel2Ok.Append(fPlotMp2->AndL() += sel);
0167 //     const TString deltaBySiName2R(fPlotMp2->Unique(Form("deltaPar%dRRel2", iPar))); //
0168 //     TH2 *hDrS2 = fPlotMp1->CreateHist2D(fPlotMp2->RPos(fPlotMp2->OrgPosT()),
0169 //                                         deltaBySi2, sel2Ok, deltaBySiName2R, "BOX");
0170 //     const TString deltaBySiName2Phi(fPlotMp2->Unique(Form("deltaPar%dPhiRel2", iPar))); //
0171 //     TH2 *hDphiS2 = fPlotMp1->CreateHist2D(fPlotMp2->Phi(fPlotMp2->OrgPosT()),
0172 //                                           deltaBySi2, sel2Ok, deltaBySiName2Phi, "BOX");
0173 //     const TString deltaBySiName2Z(fPlotMp2->Unique(Form("deltaPar%dZRel2", iPar))); //
0174 //     TH2 *hDzS2 = fPlotMp1->CreateHist2D(fPlotMp2->OrgPosT() += fPlotMp2->ZPos(),
0175 //                                         deltaBySi2, sel2Ok, deltaBySiName2Z, "BOX");
0176
0177     const TString diff(Form("%s_{2}-%s_{1}", m->NamePede(iPar).Data(), m->NamePede(iPar).Data()));
0178
0179     hDr->SetTitle(m->DelName(iPar) += " vs r" + titleAdd + ";r[cm];" + diff + ' ' + m->UnitPede(iPar));
0180     hDphi->SetTitle(m->DelName(iPar) += " vs #phi" + titleAdd + ";#phi;" + diff + ' '
0181             + m->Unit(iPar));
0182     hDz->SetTitle(m->DelName(iPar) += " vs z" + titleAdd + ";z[cm];" + diff + ' ' + m->UnitPede(iPar));
0183     hDhit->SetTitle(m->DelName(iPar) += " vs #hits" + titleAdd + ";#hits_{x};" + diff + ' '
0184             + m->UnitPede(iPar));
0185
0186     fHistManager->AddHist(hDr, layer + nPlot);
0187     fHistManager->AddHist(hDphi, layer + nPlot);
0188     fHistManager->AddHist(hDz, layer + nPlot);
0189     fHistManager->AddHist(hDhit, layer + nPlot);
0190
0191 //     if (hDrS1->GetEntries()) {
0192 //       hDrS1->SetTitle(m->DelName(iPar) += "/#sigma vs r;r[cm];(" + diff + ")/#sigma");
0193 //       hDphiS1->SetTitle(m->DelName(iPar) += "/#sigma vs #phi;#phi;(" + diff + ")/#sigma");
0194 //       hDzS1->SetTitle(m->DelName(iPar) += "/#sigma vs z;z[cm];(" + diff + ")/#sigma");
0195 //       fHistManager->AddHist(hDrS1, layer + nPlot + 1, "by #sigma_{1}", "f");
0196 //       fHistManager->AddHist(hDphiS1, layer + nPlot + 1, "by #sigma_{1}", "f");
0197 //       fHistManager->AddHist(hDzS1, layer + nPlot + 1, "by #sigma_{1}", "f");
0198 //     } else {
0199 //       delete hDrS1;
0200 //       delete hDphiS1;
0201 //       delete hDzS1;
0202 //     }
0203 //     if (hDrS2->GetEntries()) {
0204 //       hDrS2->SetLineColor(kRed);
0205 //       hDphiS2->SetLineColor(kRed);
0206 //       hDzS2->SetLineColor(kRed);
0207 //       fHistManager->AddHistSame(hDrS2, layer + nPlot + 1, 0, "by #sigma_{2}", "f");
0208 //       fHistManager->AddHistSame(hDphiS2, layer + nPlot + 1, 1, "by #sigma_{2}", "f");
0209 //       fHistManager->AddHistSame(hDzS2, layer + nPlot + 1, 2, "by #sigma_{2}", "f");
0210 //     } else {
0211 //       delete hDrS2;
0212 //       delete hDphiS2;
0213 //       delete hDzS2;
0214 //     }
0215 //     nPlot += 2;
0216     ++nPlot;
0217   }
0218
0219   const bool oldDiffStyle = fHistManager->DrawDiffStyle(false);//avoid automatic hist style changes
0220   fHistManager->Draw();
0221   fHistManager->DrawDiffStyle(oldDiffStyle);
0222 }
0223
0224 //_________________________________________________________________________________________________
0225 void CompareMillePede::DrawParam(Option_t *option)
0226 {
0227
0228   const TString opt(option);
0229
0230   const Int_t layer = this->PrepareAdd(opt.Contains("add", TString::kIgnoreCase));
0231   const TString titleAdd = this->TitleAdd();
0232
0233   const PlotMillePede *m = fPlotMp1;
0234   UInt_t nPlot = 0;
0235   for (UInt_t iPar = 0; iPar < kNpar; ++iPar) { //
0236     TString sel("");
0237     this->AddIsFreeSel(sel, opt, iPar);
0238     fPlotMp1->AddBasicSelection(sel);
0239     fPlotMp2->AddBasicSelection(sel);
0240
0241     const TString deltaPar(m->Parenth(this->DeltaPar(iPar)) += m->ToMumMuRad(iPar));
0242     const TString deltaName(m->Unique(Form("deltaPar%d", iPar))); //
0243     TH1 *hD = fPlotMp1->CreateHist(deltaPar, sel, deltaName);
0244     if (0. == hD->GetEntries()) {
0245       delete hD;
0246       continue;
0247     }
0248
0249     const TString deltaBySi1(this->DeltaParBySigma(iPar, fPlotMp1));
0250     const TString deltaBySiName1(m->Unique(Form("deltaPar%dRel1", iPar))); //
0251     TString sel1Ok(fPlotMp1->ParSiOk(iPar));
0252     if (!sel.IsNull()) sel1Ok.Append(m->AndL() += sel);
0253     TH1 *hDs1 = fPlotMp1->CreateHist(deltaBySi1, sel1Ok, deltaBySiName1);
0254
0255     const TString deltaBySi2(this->DeltaParBySigma(iPar, fPlotMp2));
0256     TString deltaBySiName2(m->Unique(Form("deltaPar%dRel2", iPar))); //
0257     m->CopyAddBinning(deltaBySiName2, hDs1);
0258     TString sel2Ok(fPlotMp2->ParSiOk(iPar));
0259     if (!sel.IsNull()) sel2Ok.Append(m->AndL() += sel);
0260     TH1 *hDs2 = fPlotMp1->CreateHist(deltaBySi2, sel2Ok, deltaBySiName2);
0261
0262     const TString diff(Form("%s_{2}-%s_{1}", m->Name(iPar).Data(), m->Name(iPar).Data()));
0263     hD->SetTitle(m->DelName(iPar) += titleAdd + ";" + diff + ' ' + m->Unit(iPar) +=";#parameters");
0264     hDs1->SetTitle(m->DelName(iPar) += "/#sigma" + titleAdd + ";(" + diff +")/#sigma;#parameters");
0265
0266     fHistManager->AddHist(hD, layer);
0267     if (hDs1->GetEntries()) {
0268       fHistManager->AddHist(hDs1, layer+1, "by #sigma_{1}");
0269     } else delete hDs1;
0270     if (hDs2->GetEntries()) { // does DrawSame work if hDs1 was empty?
0271       fHistManager->AddHistSame(hDs2, layer+1, nPlot, "by #sigma_{2}");
0272     } else delete hDs2;
0273
0274     ++nPlot;
0275   }
0276
0277   fHistManager->Draw();
0278 }
0279
0280 //_________________________________________________________________________________________________
0281 void CompareMillePede::DrawParamVsLocation(Option_t *option)
0282 {
0283
0284   const TString opt(option);
0285   const Int_t layer = this->PrepareAdd(opt.Contains("add", TString::kIgnoreCase));
0286   const TString titleAdd = this->TitleAdd();
0287
0288   const PlotMillePede *m = fPlotMp1;
0289   UInt_t nPlot = 0;
0290   for (UInt_t iPar = 0; iPar < kNpar; ++iPar) { //
0291     TString sel("");
0292     this->AddIsFreeSel(sel, opt, iPar);
0293     fPlotMp1->AddBasicSelection(sel);
0294     fPlotMp2->AddBasicSelection(sel);
0295
0296     const TString deltaPar(fPlotMp1->Parenth(this->DeltaPar(iPar)) += fPlotMp1->ToMumMuRad(iPar));
0297
0298     const TString deltaNameR(fPlotMp1->Unique(Form("deltaPar%dR", iPar))); //
0299     TH2 *hDr = fPlotMp1->CreateHist2D(fPlotMp1->RPos(fPlotMp1->OrgPosT()), deltaPar, sel,
0300                                       deltaNameR, "BOX");
0301     if (0. == hDr->GetEntries()) {
0302       delete hDr;
0303       continue;
0304     }
0305
0306     const TString deltaNamePhi(fPlotMp1->Unique(Form("deltaPar%dPhi", iPar))); //
0307     TH2 *hDphi = fPlotMp1->CreateHist2D(fPlotMp1->Phi(fPlotMp1->OrgPosT()), deltaPar, sel,
0308                                       deltaNamePhi, "BOX");
0309     const TString deltaNameZ(fPlotMp1->Unique(Form("deltaPar%dZ", iPar))); //
0310     TH2 *hDz = fPlotMp1->CreateHist2D(fPlotMp1->OrgPosT() += fPlotMp1->ZPos(), deltaPar, sel,
0311                                       deltaNameZ, "BOX");
0312     // now devided by sigma prediction from 1
0313     const TString deltaBySi1(this->DeltaParBySigma(iPar, fPlotMp1));
0314     TString sel1Ok(fPlotMp1->ParSiOk(iPar));
0315     if (!sel.IsNull()) sel1Ok.Append(fPlotMp1->AndL() += sel);
0316     const TString deltaBySiName1R(fPlotMp1->Unique(Form("deltaPar%dRRel1", iPar))); //
0317     TH2 *hDrS1 = fPlotMp1->CreateHist2D(fPlotMp1->RPos(fPlotMp1->OrgPosT()),
0318                                         deltaBySi1, sel1Ok, deltaBySiName1R, "BOX");
0319     const TString deltaBySiName1Phi(fPlotMp1->Unique(Form("deltaPar%dPhiRel1", iPar))); //
0320     TH2 *hDphiS1 = fPlotMp1->CreateHist2D(fPlotMp1->Phi(fPlotMp1->OrgPosT()),
0321                                           deltaBySi1, sel1Ok, deltaBySiName1Phi, "BOX");
0322     const TString deltaBySiName1Z(fPlotMp1->Unique(Form("deltaPar%dZRel1", iPar))); //
0323     TH2 *hDzS1 = fPlotMp1->CreateHist2D(fPlotMp1->OrgPosT() += fPlotMp1->ZPos(),
0324                                         deltaBySi1, sel1Ok, deltaBySiName1Z, "BOX");
0325
0326     // same now devided by sigma prediction from 2 (but CreateHist from 1!)
0327     const TString deltaBySi2(this->DeltaParBySigma(iPar, fPlotMp2));
0328     TString sel2Ok(fPlotMp2->ParSiOk(iPar));
0329     if (!sel.IsNull()) sel2Ok.Append(fPlotMp2->AndL() += sel);
0330     const TString deltaBySiName2R(fPlotMp2->Unique(Form("deltaPar%dRRel2", iPar))); //
0331     TH2 *hDrS2 = fPlotMp1->CreateHist2D(fPlotMp2->RPos(fPlotMp2->OrgPosT()),
0332                                         deltaBySi2, sel2Ok, deltaBySiName2R, "BOX");
0333     const TString deltaBySiName2Phi(fPlotMp2->Unique(Form("deltaPar%dPhiRel2", iPar))); //
0334     TH2 *hDphiS2 = fPlotMp1->CreateHist2D(fPlotMp2->Phi(fPlotMp2->OrgPosT()),
0335                                           deltaBySi2, sel2Ok, deltaBySiName2Phi, "BOX");
0336     const TString deltaBySiName2Z(fPlotMp2->Unique(Form("deltaPar%dZRel2", iPar))); //
0337     TH2 *hDzS2 = fPlotMp1->CreateHist2D(fPlotMp2->OrgPosT() += fPlotMp2->ZPos(),
0338                                         deltaBySi2, sel2Ok, deltaBySiName2Z, "BOX");
0339
0340     const TString diff(Form("%s_{2}-%s_{1}", m->Name(iPar).Data(), m->Name(iPar).Data()));
0341
0342     hDr->SetTitle(m->DelName(iPar) += " vs r" + titleAdd + ";r[cm];" + diff +' '+ m->Unit(iPar));
0343     hDphi->SetTitle(m->DelName(iPar) += " vs #phi" + titleAdd + ";#phi;"+diff +' '+ m->Unit(iPar));
0344     hDz->SetTitle(m->DelName(iPar) += " vs z" + titleAdd + ";z[cm];" + diff +' '+ m->Unit(iPar));
0345
0346     fHistManager->AddHist(hDr, layer + nPlot);
0347     fHistManager->AddHist(hDphi, layer + nPlot);
0348     fHistManager->AddHist(hDz, layer + nPlot);
0349
0350     if (hDrS1->GetEntries()) {
0351       hDrS1->SetTitle(m->DelName(iPar) += "/#sigma vs r" +titleAdd+ ";r[cm];(" + diff + ")/#sigma");
0352       hDphiS1->SetTitle(m->DelName(iPar) += "/#sigma vs #phi" + titleAdd + ";#phi;("
0353             + diff + ")/#sigma");
0354       hDzS1->SetTitle(m->DelName(iPar) += "/#sigma vs z" +titleAdd+ ";z[cm];(" + diff + ")/#sigma");
0355       fHistManager->AddHist(hDrS1, layer + nPlot + 1, "by #sigma_{1}", "f");
0356       fHistManager->AddHist(hDphiS1, layer + nPlot + 1, "by #sigma_{1}", "f");
0357       fHistManager->AddHist(hDzS1, layer + nPlot + 1, "by #sigma_{1}", "f");
0358     } else {
0359       delete hDrS1;
0360       delete hDphiS1;
0361       delete hDzS1;
0362     }
0363     if (hDrS2->GetEntries()) {
0364       hDrS2->SetLineColor(kRed);
0365       hDphiS2->SetLineColor(kRed);
0366       hDzS2->SetLineColor(kRed);
0367       fHistManager->AddHistSame(hDrS2, layer + nPlot + 1, 0, "by #sigma_{2}", "f");
0368       fHistManager->AddHistSame(hDphiS2, layer + nPlot + 1, 1, "by #sigma_{2}", "f");
0369       fHistManager->AddHistSame(hDzS2, layer + nPlot + 1, 2, "by #sigma_{2}", "f");
0370     } else {
0371       delete hDrS2;
0372       delete hDphiS2;
0373       delete hDzS2;
0374     }
0375     nPlot += 2;
0376   }
0377
0378   const bool oldDiffStyle = fHistManager->DrawDiffStyle(false);//avoid automatic hist style changes
0379   fHistManager->Draw();
0380   fHistManager->DrawDiffStyle(oldDiffStyle);
0381 }
0382
0383 //_________________________________________________________________________________________________
0384 void CompareMillePede::DrawParamDeltaMis(Option_t *option)
0385 {
0386   //"add": keep old canvas
0387   const TString opt(option);
0388
0389   const Int_t layer = this->PrepareAdd(opt.Contains("add", TString::kIgnoreCase));
0390   const TString titleAdd = this->TitleAdd();
0391
0392   const PlotMillePede *m = fPlotMp1;
0393   UInt_t nPlot = 0;
0394   for (UInt_t iPar = 0; iPar < kNpar; ++iPar) { //
0395     TString sel("");
0396     this->AddIsFreeSel(sel, opt, iPar);
0397     fPlotMp1->AddBasicSelection(sel);
0398     fPlotMp2->AddBasicSelection(sel);
0399
0400     const TString deltaMisPar(m->Parenth(this->DeltaMisPar(iPar)) += m->ToMumMuRad(iPar));
0401     const TString deltaName(m->Unique(Form("deltaMisPar%d", iPar))); //
0402     TH1 *hD = fPlotMp1->CreateHist(deltaMisPar, sel, deltaName);
0403     if (0. == hD->GetEntries()) {
0404       delete hD;
0405       continue;
0406     }
0407
0408     const TString deltaMisBySi1(this->DeltaMisParBySigma(iPar, fPlotMp1));
0409     const TString deltaBySiName1(m->Unique(Form("deltaMisPar%dRel1", iPar))); //
0410     TString sel1Ok(fPlotMp1->ParSiOk(iPar));
0411     if (!sel.IsNull()) sel1Ok.Append(m->AndL() += sel);
0412     TH1 *hDs1 = fPlotMp1->CreateHist(deltaMisBySi1, sel1Ok, deltaBySiName1);
0413
0414     const TString deltaMisBySi2(this->DeltaMisParBySigma(iPar, fPlotMp2));
0415     TString deltaBySiName2(m->Unique(Form("deltaMisPar%dRel2", iPar)));
0416     m->CopyAddBinning(deltaBySiName2, hDs1);
0417     TString sel2Ok(fPlotMp2->ParSiOk(iPar));
0418     if (!sel.IsNull()) sel2Ok.Append(m->AndL() += sel);
0419     TH1 *hDs2 = fPlotMp1->CreateHist(deltaMisBySi2, sel2Ok, deltaBySiName2);
0420
0421     const TString dlNm = m->DelName(iPar);
0422     hD->SetTitle(m->Name(iPar) += titleAdd
0423          + Form(";|%s_{2}|-|%s_{1}| ",dlNm.Data(), dlNm.Data())
0424          + m->Unit(iPar) += ";#parameters");
0425
0426     hDs1->SetTitle(m->Name(iPar) += titleAdd
0427            + Form(";(|%s_{2}|-|%s_{1}|)/#sigma%s", dlNm.Data(), dlNm.Data(),
0428               ";#parameters"));
0429
0430     fHistManager->AddHist(hD, layer);
0431     if (hDs1->GetEntries()) {
0432       fHistManager->AddHist(hDs1, layer+1, "by #sigma_{1}");
0433     } else delete hDs1;
0434     if (hDs2->GetEntries()) { // does DrawSame work if hDs1 was empty?
0435       fHistManager->AddHistSame(hDs2, layer+1, nPlot, "by #sigma_{2}");
0436     } else delete hDs2;
0437
0438     ++nPlot;
0439   }
0440
0441   fHistManager->Draw();
0442
0443 }
0444
0445 //_________________________________________________________________________________________________
0446 void CompareMillePede::DrawParamDeltaMisVsLoc(Option_t *option)
0447 {
0448
0449   const TString opt(option);
0450   const Int_t layer = this->PrepareAdd(opt.Contains("add", TString::kIgnoreCase));
0451   const TString titleAdd = this->TitleAdd();
0452
0453   const PlotMillePede *m = fPlotMp1;
0454   UInt_t nPlot = 0;
0455   for (UInt_t iPar = 0; iPar < kNpar; ++iPar) { //
0456     TString sel("");
0457     this->AddIsFreeSel(sel, opt, iPar);
0458     fPlotMp1->AddBasicSelection(sel);
0459     fPlotMp2->AddBasicSelection(sel);
0460
0461     const TString deltaMisPar(m->Parenth(this->DeltaMisPar(iPar)) += m->ToMumMuRad(iPar));
0462
0463     const TString deltaNameR(m->Unique(Form("deltaMisPar%dR", iPar))); //
0464     TH2 *hDr = fPlotMp1->CreateHist2D(fPlotMp1->RPos(fPlotMp1->OrgPosT()), deltaMisPar, sel,
0465                                       deltaNameR, "BOX");
0466     if (0. == hDr->GetEntries()) {
0467       delete hDr;
0468       continue;
0469     }
0470     const TString deltaNamePhi(m->Unique(Form("deltaMisPar%dPhi", iPar))); //
0471     TH2 *hDphi = fPlotMp1->CreateHist2D(fPlotMp1->Phi(fPlotMp1->OrgPosT()), deltaMisPar, sel,
0472                     deltaNamePhi, "BOX");
0473     const TString deltaNameZ(m->Unique(Form("deltaMisPar%dZ", iPar))); //
0474     TH2 *hDz = fPlotMp1->CreateHist2D(fPlotMp1->OrgPosT() += fPlotMp1->ZPos(), deltaMisPar, sel,
0475                                       deltaNameZ, "BOX");
0476     const TString deltaNameNh1(m->Unique(Form("deltaMisPar%dNh1", iPar))); //
0477     TH2 *hDnh1 = fPlotMp1->CreateHist2D(fPlotMp1->HitsX(), deltaMisPar, sel,
0478                     deltaNameNh1, "BOX");
0479     const TString deltaNameNh2(m->Unique(Form("deltaMisPar%dNh2", iPar))); //
0480     TH2 *hDnh2 = fPlotMp1->CreateHist2D(fPlotMp2->HitsX(), deltaMisPar, sel,
0481                     deltaNameNh2, "BOX");
0482
0483     // now devided by sigma prediction from 1
0484     const TString deltaMisBySi1(this->DeltaMisParBySigma(iPar, fPlotMp1));
0485     TString sel1Ok(fPlotMp1->ParSiOk(iPar));
0486     if (!sel.IsNull()) sel1Ok.Append(m->AndL() += sel);
0487     const TString deltaBySiName1R(m->Unique(Form("deltaMisPar%dRRel1", iPar))); //
0488     TH2 *hDrS1 = fPlotMp1->CreateHist2D(fPlotMp1->RPos(fPlotMp1->OrgPosT()),
0489                                         deltaMisBySi1, sel1Ok, deltaBySiName1R, "BOX");
0490     const TString deltaBySiName1Phi(m->Unique(Form("deltaMisPar%dPhiRel1", iPar))); //
0491     TH2 *hDphiS1 = fPlotMp1->CreateHist2D(fPlotMp1->Phi(fPlotMp1->OrgPosT()),
0492                                           deltaMisBySi1, sel1Ok, deltaBySiName1Phi, "BOX");
0493     const TString deltaBySiName1Z(m->Unique(Form("deltaMisPar%dZRel1", iPar))); //
0494     TH2 *hDzS1 = fPlotMp1->CreateHist2D(fPlotMp1->OrgPosT() += fPlotMp1->ZPos(),
0495                                         deltaMisBySi1, sel1Ok, deltaBySiName1Z, "BOX");
0496
0497     // same now devided by sigma prediction from 2 (but CreateHist from 1!)
0498     const TString deltaMisBySi2(this->DeltaMisParBySigma(iPar, fPlotMp2));
0499     TString sel2Ok(fPlotMp2->ParSiOk(iPar));
0500     if (!sel.IsNull()) sel2Ok.Append(m->AndL() += sel);
0501     const TString deltaBySiName2R(m->Unique(Form("deltaMisPar%dRRel2", iPar))); //
0502     TH2 *hDrS2 = fPlotMp1->CreateHist2D(fPlotMp2->RPos(fPlotMp2->OrgPosT()),
0503                                         deltaMisBySi2, sel2Ok, deltaBySiName2R, "BOX");
0504     const TString deltaBySiName2Phi(m->Unique(Form("deltaMisPar%dPhiRel2", iPar))); //
0505     TH2 *hDphiS2 = fPlotMp1->CreateHist2D(fPlotMp2->Phi(fPlotMp2->OrgPosT()),
0506                                           deltaMisBySi2, sel2Ok, deltaBySiName2Phi, "BOX");
0507     const TString deltaBySiName2Z(m->Unique(Form("deltaMisPar%dZRel2", iPar))); //
0508     TH2 *hDzS2 = fPlotMp1->CreateHist2D(fPlotMp2->OrgPosT() += fPlotMp2->ZPos(),
0509                                         deltaMisBySi2, sel2Ok, deltaBySiName2Z, "BOX");
0510
0511
0512     const TString diff(Form("|%s_{2}|-|%s_{1}|",m->DelName(iPar).Data(),m->DelName(iPar).Data()));
0513
0514     hDr->SetTitle(m->DelName(iPar) += " vs r" + titleAdd + ";r[cm];" + diff + ' ' + m->Unit(iPar));
0515     hDphi->SetTitle(m->DelName(iPar) += " vs #phi" + titleAdd + ";#phi;" + diff + ' '
0516             + m->Unit(iPar));
0517     hDz->SetTitle(m->DelName(iPar) += " vs z" + titleAdd + ";z[cm];" + diff + ' ' + m->Unit(iPar));
0518     hDnh1->SetTitle(m->DelName(iPar) += " vs N_{1}(hit" + titleAdd + ";N_{1}(hit);"
0519             + diff + ' ' + m->Unit(iPar));
0520     hDnh2->SetTitle(m->DelName(iPar) += " vs N_{2}(hit" + titleAdd + ";N_{2}(hit);"
0521             + diff + ' ' + m->Unit(iPar));
0522
0523     fHistManager->AddHist(hDr, layer + nPlot);
0524     fHistManager->AddHist(hDphi, layer + nPlot);
0525     fHistManager->AddHist(hDz, layer + nPlot);
0526     fHistManager->AddHist(hDnh1, layer + nPlot);
0527     fHistManager->AddHist(hDnh2, layer + nPlot);
0528
0529     if (hDrS1->GetEntries()) {
0530       hDrS1->SetTitle(m->DelName(iPar) += "/#sigma vs r" + titleAdd + ";r[cm];(" + diff
0531               + ")/#sigma");
0532       hDphiS1->SetTitle(m->DelName(iPar) += "/#sigma vs #phi" + titleAdd + ";#phi;(" + diff
0533             + ")/#sigma");
0534       hDzS1->SetTitle(m->DelName(iPar) += "/#sigma vs z" + titleAdd + ";z[cm];(" + diff
0535               + ")/#sigma");
0536       fHistManager->AddHist(hDrS1, layer + nPlot + 1, "by #sigma_{1}", "f");
0537       fHistManager->AddHist(hDphiS1, layer + nPlot + 1, "by #sigma_{1}", "f");
0538       fHistManager->AddHist(hDzS1, layer + nPlot + 1, "by #sigma_{1}", "f");
0539     } else {
0540       delete hDrS1;
0541       delete hDphiS1;
0542       delete hDzS1;
0543     }
0544     if (hDrS2->GetEntries()) {
0545       hDrS2->SetLineColor(kRed);
0546       hDphiS2->SetLineColor(kRed);
0547       hDzS2->SetLineColor(kRed);
0548       fHistManager->AddHistSame(hDrS2, layer + nPlot + 1, 0, "by #sigma_{2}", "f");
0549       fHistManager->AddHistSame(hDphiS2, layer + nPlot + 1, 1, "by #sigma_{2}", "f");
0550       fHistManager->AddHistSame(hDzS2, layer + nPlot + 1, 2, "by #sigma_{2}", "f");
0551     } else {
0552       delete hDrS2;
0553       delete hDphiS2;
0554       delete hDzS2;
0555     }
0556     nPlot += 2;
0557   }
0558
0559   const bool oldDiffStyle = fHistManager->DrawDiffStyle(false);//avoid automatic hist style changes
0560   fHistManager->Draw();
0561   fHistManager->DrawDiffStyle(oldDiffStyle);
0562 }
0563
0564 //_________________________________________________________________________________________________
0565 void CompareMillePede::DrawNumHits(Option_t *option)
0566 {
0567   const TString opt(option);
0568   const Int_t layer = this->PrepareAdd(opt.Contains("add", TString::kIgnoreCase));
0569   const TString titleAdd = this->TitleAdd();
0570
0571   const PlotMillePede *m = fPlotMp1;
0572
0573   TString sel;
0574   fPlotMp1->AddBasicSelection(sel);
0575   fPlotMp2->AddBasicSelection(sel);
0576
0577   const TString nX1vs2(m->Unique("hitsX1vs2"));
0578   TH2 *hX1vs2 = fPlotMp1->CreateHist2D(fPlotMp1->HitsX(), fPlotMp2->HitsX(), sel, nX1vs2, "BOX");
0579
0580   const TString nY1vs2(m->Unique("hitsY1vs2"));
0581   TH2 *hY1vs2 = fPlotMp1->CreateHist2D(fPlotMp1->HitsY(), fPlotMp2->HitsY(), sel, nY1vs2, "BOX");
0582
0583   const TString nDiffX(m->Unique("deltaHitsX"));
0584   const TString diffX(m->Parenth(fPlotMp2->HitsX() += m->Min() += fPlotMp1->HitsX()));
0585   TH1 *hDiffX = fPlotMp1->CreateHist(diffX, sel, nDiffX);
0586
0587   const TString nDiffY(m->Unique("deltaHitsY"));
0588   const TString diffY(m->Parenth(fPlotMp2->HitsY() += m->Min() += fPlotMp1->HitsY()));
0589   TH1 *hDiffY = fPlotMp1->CreateHist(diffY, sel, nDiffY);
0590
0591   const TString ratioX(m->Parenth(fPlotMp2->HitsX() += m->Div() += fPlotMp1->HitsX()));
0592   const TString selRatioX(fPlotMp1->HitsX() += ">0" + m->AndL() += sel);
0593   const TString nRatioX(m->Unique("ratioHitsX"));
0594   TH1 *hRatioX = fPlotMp1->CreateHist(ratioX, selRatioX, nRatioX);
0595
0596   const TString ratioY(m->Parenth(fPlotMp2->HitsY() += m->Div() += fPlotMp1->HitsY()));
0597   const TString selRatioY(fPlotMp1->HitsY() += ">0" + m->AndL() += sel);
0598   const TString nRatioY(m->Unique("ratioHitsY"));
0599   TH1 *hRatioY = fPlotMp1->CreateHist(ratioY, selRatioY, nRatioY);
0600
0601 //   const TString nDiffYvs1(m->Unique("deltaHitsYvs1"));
0602 //   TH2 *hDiffYvs1 = fPlotMp1->CreateHist2D(diffY, fPlotMp1->HitsY(), sel, nDiffXvs1, "BOX");
0603
0604   // diff vs location/hits
0605   const TString nDiffXvsR(m->Unique("deltaHitsXvsR"));
0606   TH2 *hDiffXvsR = fPlotMp1->CreateHist2D(diffX, fPlotMp1->RPos(fPlotMp1->OrgPosT()),
0607                                           sel, nDiffXvsR, "BOX");
0608   const TString nDiffXvsPhi(m->Unique("deltaHitsXvsPhi"));
0609   TH2 *hDiffXvsPhi = fPlotMp1->CreateHist2D(diffX, fPlotMp1->Phi(fPlotMp1->OrgPosT()),
0610                                             sel, nDiffXvsPhi, "BOX");
0611   const TString nDiffXvsZ(m->Unique("deltaHitsXvsZ"));
0612   TH2 *hDiffXvsZ = fPlotMp1->CreateHist2D(diffX, fPlotMp1->OrgPosT() += fPlotMp1->ZPos(),
0613                                           sel, nDiffXvsZ, "BOX");
0614   const TString nDiffXvsX(m->Unique("deltaHitsXvsX"));
0615   TH2 *hDiffXvsX = fPlotMp1->CreateHist2D(diffX, fPlotMp1->OrgPosT() += fPlotMp1->XPos(),
0616                                           sel, nDiffXvsX, "BOX");
0617   const TString nDiffXvsY(m->Unique("deltaHitsXvsY"));
0618   TH2 *hDiffXvsY = fPlotMp1->CreateHist2D(diffX, fPlotMp1->OrgPosT() += fPlotMp1->YPos(),
0619                                           sel, nDiffXvsY, "BOX");
0620   const TString nDiffXvs1(m->Unique("deltaHitsXvs1"));
0621   TH2 *hDiffXvs1 = fPlotMp1->CreateHist2D(diffX, fPlotMp1->HitsX(), sel, nDiffXvs1, "BOX");
0622
0623   // ratio vs location/hits
0624   const TString nRatioXvsR(m->Unique("ratioHitsXvsR"));
0625   TH2 *hRatioXvsR = fPlotMp1->CreateHist2D(ratioX, fPlotMp1->RPos(fPlotMp1->OrgPosT()),
0626                        selRatioX, nRatioXvsR, "BOX");
0627   const TString nRatioXvsPhi(m->Unique("ratioHitsXvsPhi"));
0628   TH2 *hRatioXvsPhi = fPlotMp1->CreateHist2D(ratioX, fPlotMp1->Phi(fPlotMp1->OrgPosT()),
0629                          selRatioX, nRatioXvsPhi, "BOX");
0630   const TString nRatioXvsZ(m->Unique("ratioHitsXvsZ"));
0631   TH2 *hRatioXvsZ = fPlotMp1->CreateHist2D(ratioX, fPlotMp1->OrgPosT() += fPlotMp1->ZPos(),
0632                        selRatioX, nRatioXvsZ, "BOX");
0633   const TString nRatioXvsX(m->Unique("ratioHitsXvsX"));
0634   TH2 *hRatioXvsX = fPlotMp1->CreateHist2D(ratioX, fPlotMp1->OrgPosT() += fPlotMp1->XPos(),
0635                        selRatioX, nRatioXvsX, "BOX");
0636   const TString nRatioXvsY(m->Unique("ratioHitsXvsY"));
0637   TH2 *hRatioXvsY = fPlotMp1->CreateHist2D(ratioX, fPlotMp1->OrgPosT() += fPlotMp1->YPos(),
0638                        selRatioX, nRatioXvsY, "BOX");
0639   const TString nRatioXvs1(m->Unique("ratioHitsXvs1"));
0640   TH2 *hRatioXvs1 = fPlotMp1->CreateHist2D(ratioX, fPlotMp1->HitsX(), selRatioX, nRatioXvs1, "BOX");
0641
0642
0643   hX1vs2->SetTitle("#hits_{x}" + titleAdd + ";N_{hit,x}^{1};N_{hit,x}^{2}");
0644   hY1vs2->SetTitle("#hits_{y}" + titleAdd + ";N_{hit,y}^{1};N_{hit,y}^{2}");
0645   hDiffX->SetTitle("#Delta#hits_{x}" + titleAdd + ";N_{hit,x}^{2} - N_{hit,x}^{1};#alignables");
0646   hDiffY->SetTitle("#Delta#hits_{y}" + titleAdd + ";N_{hit,y}^{2} - N_{hit,y}^{1};#alignables");
0647 //   hDiffYvs1->SetTitle("#Delta#hits_{y} vs #hits_{y}^{1}" + titleAdd
0648 //            + ";N_{hit,y}^{2} - N_{hit,y}^{1};N_{hit,y}^{1}");
0649   hRatioX->SetTitle("ratio #hits_{x}" + titleAdd + ";N_{hit,x}^{2} / N_{hit,x}^{1};#alignables");
0650   hRatioY->SetTitle("ratio #hits_{y}" + titleAdd + ";N_{hit,y}^{2} / N_{hit,y}^{1};#alignables");
0651
0652   hDiffXvsR->SetTitle("#Delta#hits_{x} vs r" + titleAdd + ";N_{hit,x}^{2} - N_{hit,x}^{1};r[cm]");
0653   hDiffXvsPhi->SetTitle("#Delta#hits_{x} vs #phi" +titleAdd+ ";N_{hit,x}^{2} - N_{hit,x}^{1};#phi");
0654   hDiffXvsZ->SetTitle("#Delta#hits_{x} vs z" + titleAdd + ";N_{hit,x}^{2} - N_{hit,x}^{1};z[cm]");
0655   hDiffXvsX->SetTitle("#Delta#hits_{x} vs x" + titleAdd + ";N_{hit,x}^{2} - N_{hit,x}^{1};x[cm]");
0656   hDiffXvsY->SetTitle("#Delta#hits_{x} vs y" + titleAdd + ";N_{hit,x}^{2} - N_{hit,x}^{1};y[cm]");
0657   hDiffXvs1->SetTitle("#Delta#hits_{x} vs #hits_{x}^{1}" + titleAdd
0658               + ";N_{hit,x}^{2} - N_{hit,x}^{1};N_{hit,x}^{1}");
0659
0660   hRatioXvsR->SetTitle("hit ratio vs r" + titleAdd + ";N_{hit,x}^{2} / N_{hit,x}^{1};r[cm]");
0661   hRatioXvsPhi->SetTitle("hit ratio vs #phi" +titleAdd+ ";N_{hit,x}^{2} / N_{hit,x}^{1};#phi");
0662   hRatioXvsZ->SetTitle("hit ratio vs z" + titleAdd + ";N_{hit,x}^{2} / N_{hit,x}^{1};z[cm]");
0663   hRatioXvsX->SetTitle("hit ratio vs x" + titleAdd + ";N_{hit,x}^{2} / N_{hit,x}^{1};x[cm]");
0664   hRatioXvsY->SetTitle("hit ratio vs y" + titleAdd + ";N_{hit,x}^{2} / N_{hit,x}^{1};y[cm]");
0665   hRatioXvs1->SetTitle("hit ratio vs #hits_{x}^{1}" + titleAdd
0666                + ";N_{hit,x}^{2} / N_{hit,x}^{1};N_{hit,x}^{1}");
0667
0668   fHistManager->AddHist(hX1vs2, layer);
0669   fHistManager->AddHist(hDiffX, layer);
0670   fHistManager->AddHist(hRatioX, layer);
0671   fHistManager->AddHist(hY1vs2, layer);
0672   fHistManager->AddHist(hDiffY, layer);
0673 //   fHistManager->AddHist(hDiffYvs1, layer);
0674   fHistManager->AddHist(hRatioY, layer);
0675
0676   fHistManager->AddHist(hDiffXvsR, layer+1);
0677   fHistManager->AddHist(hDiffXvsPhi, layer+1);
0678   fHistManager->AddHist(hDiffXvsZ, layer+1);
0679   fHistManager->AddHist(hDiffXvsX, layer+1);
0680   fHistManager->AddHist(hDiffXvsY, layer+1);
0681   fHistManager->AddHist(hDiffXvs1, layer+1);
0682
0683   fHistManager->AddHist(hRatioXvsR, layer+2);
0684   fHistManager->AddHist(hRatioXvsPhi, layer+2);
0685   fHistManager->AddHist(hRatioXvsZ, layer+2);
0686   fHistManager->AddHist(hRatioXvsX, layer+2);
0687   fHistManager->AddHist(hRatioXvsY, layer+2);
0688   fHistManager->AddHist(hRatioXvs1, layer+2);
0689
0690   fHistManager->Draw();
0691 }
0692
0693 //_________________________________________________________________________________________________
0694 void CompareMillePede::DrawAbsPos(Option_t *option)
0695 {
0696 //"start": at start (else end), "add": keep old canvas
0697
0698   const TString opt(option);
0699
0700   const Int_t layer = this->PrepareAdd(opt.Contains("add", TString::kIgnoreCase));
0701   const TString titleAdd = this->TitleAdd();
0702
0703   const TString tree1(opt.Contains("start", TString::kIgnoreCase) ? fPlotMp1->OrgPosT() : fPlotMp1->PosT());
0704   const TString tree2(opt.Contains("start", TString::kIgnoreCase) ? fPlotMp2->OrgPosT() : fPlotMp2->PosT());
0705
0706   const PlotMillePede *m = fPlotMp1;
0707   //UInt_t nPlot = 0;
0708
0709   // const TString posNames[] = {"rphi", "r", "z", "phi", "x", "y"};
0710
0711 //   for (UInt_t iPos = 0; iPos < sizeof(posNames)/sizeof(posNames[0]); ++iPos) {
0712 //     const TString &posName = posNames[iPos];
0713   for (UInt_t iPos = 0; iPos < 3; ++iPos) { //
0714     TString sel("");
0715     //this->AddIsFreeSel(sel, opt, iPar);
0716     fPlotMp1->AddBasicSelection(sel);
0717     fPlotMp2->AddBasicSelection(sel);
0718
0719 //   // values: absolute positions and orientations; need a 'position tree' as argument or in front
0720 //   TString Pos() const { return "Pos";}
0721 //   TString Pos(UInt_t ui) const { return Pos() += Bracket(ui);}
0722 //   TString XPos() const { return Pos(0);}
0723 //   TString YPos() const { return Pos(1);}
0724 //   TString ZPos() const { return Pos(2);}
0725
0726 //   TString DelPos(UInt_t ui, const TString &tree1, const TString &tree2) const;
0727
0728 //   // units and names for position strings (r, rphi, phi, x, y, z)
0729 //   TString ToMumMuRad(const TString &pos) const;
0730 //   TString Name(const TString &pos) const;
0731
0732 //   TString OrgPosT() const { return fOrgPos + Dot();} // nominal global positions
0733 //   TString PosT() const { return fPos + Dot();} // aligned global positions
0734
0735
0736     const TString deltaName(m->Unique(Form("deltaPos%d", iPos))); //
0737     //const TString deltaName(m->Unique("deltaPos" + posName)); //
0738     TH1 *h = fPlotMp1->CreateHist(m->DelPos(iPos, tree2, tree1), sel, deltaName);
0739     if (0. == h->GetEntries()) {
0740       delete h;
0741       continue;
0742     }
0743
0744     const TString diff(Form("%s_{2}-%s_{1}", m->NamePos(iPos).Data(), m->NamePos(iPos).Data()));
0745     h->SetTitle(diff + titleAdd + ";" + diff + " [cm];#alignables");
0746
0747     fHistManager->AddHist(h, layer);
0748     //    ++nPlot;
0749   }
0750
0751   fHistManager->Draw();
0752 }
0753
0754 //_______________________________________________________________________________
0755 void CompareMillePede::DrawSurfaceDeformations(Option_t *option, unsigned int firstPar,
0756                            unsigned int lastPar, const TString &whichOne)
0757 {
0758    const TString opt(option);
0759
0760    const Int_t layer = this->PrepareAdd(opt.Contains("add", TString::kIgnoreCase));
0761    const bool limits = opt.Contains("limit", TString::kIgnoreCase);
0762    const bool noVs = opt.Contains("novs", TString::kIgnoreCase);
0763    const bool noDiff = opt.Contains("noDiff", TString::kIgnoreCase);
0764    const TString titleAdd = this->TitleAdd();
0765    const PlotMillePede *m = fPlotMp1;
0766
0767    UInt_t nPlot = 0;
0768    for (UInt_t iPar = firstPar; iPar <= lastPar; ++iPar) { //
0769      TString sel("");
0770      fPlotMp1->AddBasicSelection(sel);
0771      fPlotMp2->AddBasicSelection(sel);
0772      sel += m->AndL() += m->Parenth(fPlotMp1->NumDeformValues(whichOne) += Form(">%u", iPar))
0773        += m->AndL() += m->Parenth(fPlotMp2->NumDeformValues(whichOne) += Form(">%u", iPar));
0774
0775      const TString deform1(fPlotMp1->DeformValue(iPar, whichOne));
0776      const TString deform2(fPlotMp2->DeformValue(iPar, whichOne));
0777      const TString diff(m->Parenth(deform2 + m->Min() += deform1)
0778             += m->ToMumMuRadSurfDef(iPar));
0779      TString hNameDiff(m->Unique(Form("hSurfDiff%s%u", whichOne.Data(), iPar)));
0780      if (limits) hNameDiff += Form("(101,%f,%f)", m->GetMaxDevDown(), m->GetMaxDevUp());
0781
0782      // delta plot
0783      TH1 *hDiff = fPlotMp1->CreateHist(diff, sel, hNameDiff);
0784      if (0. == hDiff->GetEntries()) {
0785        delete hDiff;
0786        continue;
0787      }
0788
0789      if (noDiff) {
0790        delete hDiff;
0791      } else {
0792        // titles
0793        const TString diffTit(m->NameSurfDef(iPar) += "_{2}-"
0794                  + m->NameSurfDef(iPar) += "_{1}");
0795        hDiff->SetTitle(diffTit + titleAdd + ";"
0796                + diffTit + m->UnitSurfDef(iPar) +=";#parameters");
0797        fHistManager->AddHist(hDiff, layer);
0798      }
0799
0800      if (!noVs) { // 2D plot
0801        TString hNameVs(m->Unique(Form("hSurfVs%s%u", whichOne.Data(), iPar)));
0802        if (limits) hNameVs += Form("(101,%f,%f, 101,%f,%f)", m->GetMaxDevDown(),
0803                    m->GetMaxDevUp(), m->GetMaxDevDown(), m->GetMaxDevUp());
0804        TH2 *hVs = fPlotMp1->CreateHist2D(deform1 + m->ToMumMuRadSurfDef(iPar),
0805                      deform2 + m->ToMumMuRadSurfDef(iPar), sel,
0806                      hNameVs, "BOX");
0807        hVs->SetTitle(m->NameSurfDef(iPar) += titleAdd + ";"
0808            + m->NameSurfDef(iPar) += "_{1}" + m->UnitSurfDef(iPar) +=";"
0809            + m->NameSurfDef(iPar) += "_{2}" + m->UnitSurfDef(iPar));
0810
0811        fHistManager->AddHist(hVs, layer + !noDiff);
0812        // add a TLine
0813        const Double_t xMin = hVs->GetXaxis()->GetXmin();
0814        const Double_t xMax = hVs->GetXaxis()->GetXmax();
0815        TLine *line = new TLine(xMin, xMin, xMax, xMax);
0816        line->SetLineColor(kRed); line->SetLineWidth(2);
0817        fHistManager->AddObject(line, layer + !noDiff, nPlot);
0818      }
0819
0820      ++nPlot;
0821    }
0822
0823    fHistManager->Draw();
0824 }
0825 //_________________________________________________________________________________________________
0826 //_________________________________________________________________________________________________
0827 //_________________________________________________________________________________________________
0828 bool CompareMillePede::IsConsistent()
0829 {
0830
0831   // Id and ObjId identify an alignable, so the trees are in same order if the difference
0832   // is always zero. We check this by checking mean and RMS of the difference.
0833
0834   TH1 *hId = fPlotMp1->CreateHist(fPlotMp1->PosT() += "Id - " + fPlotMp2->PosT() += "Id", "");
0835   TH1 *hObjId = fPlotMp1->CreateHist(fPlotMp1->PosT()+="ObjId - " + fPlotMp2->PosT()+="ObjId", "");
0836
0837   if (hId->GetMean() == 0. && hId->GetRMS() == 0.
0838       && hObjId->GetMean() == 0. && hObjId->GetRMS() == 0.) {
0839     delete hId;
0840     delete hObjId;
0841     return true;
0842   } else {
0843     fHistManager->Clear();
0844     fHistManager->AddHist(hId);
0845     fHistManager->AddHist(hObjId);
0846     fHistManager->Draw();
0847
0848     return false;
0849   }
0850 }
0851
0852 //_________________________________________________________________________________________________
0853 TString CompareMillePede::DeltaPar(UInt_t iPar) const
0854 {
0855   // '2' - '1'
0856   return fPlotMp2->FinalMisAlignment(iPar) += fPlotMp1->Min() += fPlotMp1->FinalMisAlignment(iPar);
0857 }
0858
0859 //_________________________________________________________________________________________________
0860 TString CompareMillePede::DeltaParBySigma(UInt_t iPar, const PlotMillePede *sigmaSource) const
0861 {
0862
0863   return
0864     sigmaSource->Parenth(this->DeltaPar(iPar)) += sigmaSource->Div() += sigmaSource->ParSi(iPar);
0865 }
0866
0867 //_________________________________________________________________________________________________
0868 TString CompareMillePede::DeltaMisPar(UInt_t iPar) const
0869 {
0870   // '2' - '1'
0871   const PlotMillePede *m = fPlotMp1;
0872   return            m->Abs(fPlotMp2->FinalMisAlignment(iPar))
0873     += m->Min() += m->Abs(fPlotMp1->FinalMisAlignment(iPar));
0874 }
0875
0876 //_________________________________________________________________________________________________
0877 TString CompareMillePede::DeltaMisParBySigma(UInt_t iPar, const PlotMillePede *sigmaSource) const
0878 {
0879
0880   return sigmaSource->Parenth(this->DeltaMisPar(iPar))
0881     += sigmaSource->Div() += sigmaSource->ParSi(iPar);
0882 }
0883
0884 //_________________________________________________________________________________________________
0885 TString CompareMillePede::DeltaPos(UInt_t iPos) const
0886 {
0887   TString pos1(fPlotMp1->Parenth(fPlotMp1->PosT() += fPlotMp1->Pos(iPos)));
0888   TString pos2(fPlotMp2->Parenth(fPlotMp2->PosT() += fPlotMp2->Pos(iPos)));
0889
0890   return pos1 += fPlotMp1->Min() += pos2;
0891 }
0892
0893 //_________________________________________________________________________________________________
0894 void CompareMillePede::AddIsFreeSel(TString &sel, const TString &option, UInt_t iPar) const
0895 {
0896   if (option.Contains("free1", TString::kIgnoreCase)) {
0897     if (sel.IsNull()) sel = fPlotMp1->Fixed(iPar, false);
0898     else        sel.Prepend(fPlotMp1->Fixed(iPar, false) += fPlotMp1->AndL());
0899   }
0900   if (option.Contains("free2", TString::kIgnoreCase)) {
0901     if (sel.IsNull()) sel = fPlotMp2->Fixed(iPar, false);
0902     else        sel.Prepend(fPlotMp2->Fixed(iPar, false) += fPlotMp2->AndL());
0903   }
0904 }
0905
0906 //_________________________________________________________________________________________________
0907 void CompareMillePede::SetSubDetId(Int_t subDetId)
0908 {
0909   fPlotMp1->SetSubDetId(subDetId);
0910   fPlotMp2->SetSubDetId(subDetId);
0911 }
0912
0913 //_________________________________________________________________________________________________
0914 void CompareMillePede::AddSubDetId(Int_t subDetId)
0915 {
0916   fPlotMp1->AddSubDetId(subDetId);
0917   fPlotMp2->AddSubDetId(subDetId);
0918 }
0919
0920 //_________________________________________________________________________________________________
0921 void CompareMillePede::SetAlignableTypeId(Int_t alignableTypeId)
0922 {
0923   fPlotMp1->SetAlignableTypeId(alignableTypeId);
0924   fPlotMp2->SetAlignableTypeId(alignableTypeId);
0925 }
0926
0927 //_________________________________________________________________________________________________
0928 void CompareMillePede::SetHieraLevel(Int_t hieraLevel)
0929 {
0930   fPlotMp1->SetHieraLevel(hieraLevel);
0931   fPlotMp2->SetHieraLevel(hieraLevel);
0932 }
0933
0934 //_________________________________________________________________________________________________
0935 void CompareMillePede::AddAdditionalSel(const char *selection)
0936 {
0937   fPlotMp1->AddAdditionalSel(selection);
0938   fPlotMp2->AddAdditionalSel(selection);
0939 }
0940
0941 //_________________________________________________________________________________________________
0942 void CompareMillePede::AddAdditionalSel(const TString &xyzrPhiNhit, Float_t min, Float_t max)
0943 {
0944   fPlotMp1->AddAdditionalSel(xyzrPhiNhit, min, max);
0945   fPlotMp2->AddAdditionalSel(xyzrPhiNhit, min, max);
0946 }
0947
0948 //_________________________________________________________________________________________________
0949 void CompareMillePede::ClearAdditionalSel ()
0950 {
0951   fPlotMp1->ClearAdditionalSel();
0952   fPlotMp2->ClearAdditionalSel();
0953 }
0954
0955 //_________________________________________________________________________________________________
0956 void CompareMillePede::SetSurfDefDeltaBows(bool deltaBows)
0957 {
0958   fPlotMp1->SetSurfDefDeltaBows(deltaBows);
0959   fPlotMp2->SetSurfDefDeltaBows(deltaBows);
0960 }
0961
0962 //_________________________________________________________________________________________________
0963 Int_t CompareMillePede::PrepareAdd(bool addPlots)
0964 {
0965   if (addPlots) {
0966     return fHistManager->GetNumLayers();
0967   } else {
0968     fHistManager->Clear();
0969     return 0;
0970   }
0971 }
0972
0973 //_________________________________________________________________________________________________
0974 TString CompareMillePede::TitleAdd() const
0975 {
0976   // try to avoid duplication of cuts in title
0977
0978   // get title add from titleAdd for first
0979   TString titleAdd1 = fPlotMp1->TitleAdd();
0980   const TString title1(fPlotMp1->GetTitle());
0981   TString titleAdd1noTitle(titleAdd1);
0982   const Ssiz_t tit1pos = titleAdd1.First(title1);
0983   if (tit1pos != kNPOS) titleAdd1noTitle.Remove(tit1pos, title1.Length() + 2); //+2: ', '
0984   // now for second
0985   TString titleAdd2 = fPlotMp2->TitleAdd();
0986   const TString title2(fPlotMp2->GetTitle());
0987   TString titleAdd2noTitle(titleAdd2);
0988   const Ssiz_t tit2pos = titleAdd2.First(title2);
0989   if (tit2pos != kNPOS) titleAdd2noTitle.Remove(tit2pos, title2.Length() + 2); //+2: ', '
0990
0991   if (titleAdd1noTitle != titleAdd2noTitle) {
0992     // shortened titleAdds different: keep full length, adding 1/2 info
0993     titleAdd1 +=  " (1), ";
0994     if (titleAdd2(0, 2) == ": ") titleAdd2.Remove(0, 2); // remove ": "
0995     titleAdd1 += titleAdd2 += " (2)";
0996     return titleAdd1;
0997   } else {
0998     // otherwise keep only one, add titles at the end
0999     return (titleAdd1noTitle += ", 1: " + title1 + ", 2: " + title2);
1000   }
1001 }