MillePedeAlignmentAlgorithm/macros/PlotMillePedeIOV.C

0001 // Original Author: Gero Flucke
0002 // last change    : $Date: 2012/03/29 08:42:23 $
0003 // by             : $Author: flucke $
0004
0005 #include <map>
0006 #include <vector>
0007 #include <iostream>
0008
0009 #include "TString.h"
0010 #include "TTree.h"
0011 #include "TGraph.h"
0012 #include "TGraphErrors.h"
0013 #include "TMultiGraph.h"
0014 #include "TCanvas.h"
0015 #include "TROOT.h"
0016 #include "TFile.h"
0017 #include "TList.h"
0018 #include "TLegend.h"
0019
0020 #include "PlotMillePede.h"
0021 #include "PlotMillePedeIOV.h"
0022 #include "GFUtils/GFHistManager.h"
0023
0024 PlotMillePedeIOV::PlotMillePedeIOV(const char *fileName, Int_t maxIov,
0025                    Int_t hieraLevel)
0026   : fHistManager(new GFHistManager), fTitle("")
0027 {
0028   if (maxIov <= 0) { // find maximum IOV in file if not specified
0029     maxIov = 0;
0030     TFile *file = TFile::Open(fileName);
0031     if (!file) return;
0032
0033     // Loop as long as we find MillePedeUser_<n> trees:
0034     TTree *tree = 0;
0035     do {
0036       file->GetObject(Form("MillePedeUser_%d", ++maxIov), tree);
0037     } while (tree);
0038     --maxIov; // one step back since last try did not succeed
0039
0040     delete file;
0041   }
0042
0043   for (Int_t iov = 1; iov <= maxIov; ++iov) {
0044     fIovs.push_back(new PlotMillePede(fileName, iov, hieraLevel));
0045     fIovs.back()->GetHistManager()->SetBatch();
0046     TTree *tree = fIovs.back()->GetMainTree();
0047     tree->SetEstimate(tree->GetEntries()); // for safe use of tree->GetV1() etc.
0048   }
0049
0050   fHistManager->SetLegendX1Y1X2Y2(0.14, 0.7, 0.45, 0.9);
0051 }
0052
0053 PlotMillePedeIOV::~PlotMillePedeIOV()
0054 {
0055   delete fHistManager;
0056   for (unsigned int i = 0; i < fIovs.size(); ++i) {
0057     delete fIovs[i];
0058   }
0059 }
0060
0061 ////////////////////////////////////////////////////////////////////////////////
0062 void PlotMillePedeIOV::DrawPedeParam(Option_t *option, unsigned int nNonRigidParam)
0063 {
0064   const unsigned int nPar = PlotMillePede::kNpar + nNonRigidParam;
0065   const Int_t layer = this->PrepareAdd(TString(option).Contains("add", TString::kIgnoreCase));
0066   const Int_t xInTitle = TString(option).Contains("x", TString::kIgnoreCase);
0067   const Int_t yInTitle = TString(option).Contains("y", TString::kIgnoreCase);
0068   const Int_t zInTitle = TString(option).Contains("z", TString::kIgnoreCase);
0069   // if no position is selected for legend, use DetId:
0070   const Int_t idInTitle = (!xInTitle && !yInTitle && !zInTitle)
0071     || TString(option).Contains("id", TString::kIgnoreCase);
0072   const Int_t doErr = TString(option).Contains("err", TString::kIgnoreCase);
0073   const Int_t noErr = TString(option).Contains("val", TString::kIgnoreCase);
0074
0075   typedef std::map<ParId, TGraphErrors*> GraphMap;
0076   GraphMap graphs; // size will be nPar*numAlis (if same nPar for all alignables)
0077
0078   unsigned int maxPar = 0;
0079   for (unsigned int iIov = 0; iIov < fIovs.size(); ++iIov) { // loop in IOVs
0080     PlotMillePede *iov = fIovs[iIov];
0081     TTree *tree = iov->GetMainTree();
0082     for (unsigned int iPar = 0; iPar < nPar; ++iPar) { // loop on parameters
0083       if (iPar > maxPar) maxPar = iPar;
0084       const TString pedePar(iov->MpT() += iov->Par(iPar) += iov->ToMumMuRadPede(iPar));
0085       const TString parSi(iov->ParSi(iPar) += iov->ToMumMuRadPede(iPar));
0086       const TString par(doErr ? parSi : pedePar); // decide to do param or error
0087       // The selection is tree-name (but not parameter...) dependent:
0088       TString selection;//(iov->Valid(iPar)); // FIXME??
0089       iov->AddBasicSelection(selection);
0090       // main command to get parameter & also (Det)Id and ObjId (e.g. TPBLayer)
0091       const Long64_t numAlis = tree->Draw(par + ":Id:ObjId:" + parSi, selection, "goff");
0092       // copy result of above Draw(..)
0093       const std::vector<double> values(tree->GetV1(), tree->GetV1() + numAlis);
0094       const std::vector<double> ids   (tree->GetV2(), tree->GetV2() + numAlis);
0095       const std::vector<double> objIds(tree->GetV3(), tree->GetV3() + numAlis);
0096       const std::vector<double> sigmas(tree->GetV4(), tree->GetV4() + numAlis);
0097
0098       // now loop on selected alignables and create/fill graphs
0099       for (Long64_t iAli = 0; iAli < numAlis; ++iAli) {
0100     // ParId is Id, ObjId and parameter number - used as key for the map
0101     const ParId id(ids[iAli], objIds[iAli], iPar);
0102     TGraphErrors *&gr = graphs[id]; // pointer by ref (might be created!)
0103     if (!gr) {
0104       // this tree->Draw(..) is only needed if xyz position requested...
0105       tree->Draw(iov->XPos() += ":" + iov->YPos() +=
0106              ":" + iov->ZPos(), selection, "goff");
0107       gr = new TGraphErrors; // Assigns value to map-internal pointer. (!)
0108       // We define title for legend here:
0109       TString title(iov->AlignableObjIdString(id.objId_));
0110       if (idInTitle)title += Form(", %d", id.id_);
0111       if (xInTitle) title += Form(", x=%.1f", tree->GetV1()[iAli]);
0112       if (yInTitle) title += Form(", y=%.1f", tree->GetV2()[iAli]);
0113       if (zInTitle) title += Form(", z=%.f", tree->GetV3()[iAli]);
0114       // if (title.Last(',') != kNPOS) title.Remove(title.Last(','));
0115       gr->SetTitle(title);
0116     }
0117     gr->SetPoint(gr->GetN(), iIov+1, values[iAli]); // add new point for IOV
0118     if (!doErr      // not if we plot error instead of value
0119         && !noErr   // not if error bar is deselected
0120         && sigmas[iAli] > 0.) { // determined by pede (inversion...)
0121       gr->SetPointError(gr->GetN()-1, 0., sigmas[iAli]); // add error in y
0122         }
0123       } // end loop on alignables
0124     } // end loop on parameters
0125   } // end loop on IOVs
0126
0127   //  if (graphs.empty()) return; //
0128
0129   // Now we have all graphs filled - order them such that
0130   // the same parameter is attached to the same multigraph:
0131   std::vector<TMultiGraph*> multis(maxPar+1); // one multigraph per parameter
0132   for (GraphMap::const_iterator iGr = graphs.begin(); iGr != graphs.end();
0133        ++iGr) { // loop on map of graphs
0134     const ParId &id = iGr->first;
0135     if (!multis[id.par_]) multis[id.par_] = new TMultiGraph;
0136     multis[id.par_]->Add(iGr->second, ""); // add an option?
0137   }
0138
0139   // Need to draw the multigraph to get its histogram for the axes.
0140   // This histogram has to be given to the hist manager before the multigraphs.
0141   // Therefore set ROOT to batch and prepare a temporary TCanvas for drawing.
0142   const bool isBatch = gROOT->IsBatch(); gROOT->SetBatch();
0143   TCanvas c; // On stack: later goes out of scope...
0144   // Finally loop on multigraphs
0145   for (unsigned int iMulti = 0; iMulti < multis.size(); ++iMulti) {
0146     if (multis[iMulti] == 0) continue;
0147     TIter iter(multis[iMulti]->GetListOfGraphs());
0148     multis[iMulti]->Draw("ALP"); //'Draw' graph into batch canvas to create hist
0149     TH1 *h = static_cast<TH1*>(multis[iMulti]->GetHistogram() //...but clone it:
0150                    ->Clone(this->Unique(Form("IOV%u", iMulti))));
0151     const PlotMillePede *i0 = fIovs[0]; // IOV does not matter here...
0152     if (doErr) { // title if we draw error instead of value
0153       const TString errPar(Form("#sigma(%s)", i0->NamePede(iMulti).Data()));
0154       h->SetTitle(errPar + " IOVs" += i0->TitleAdd() += ";IOV;"
0155           + errPar + i0->UnitPede(iMulti));
0156     } else if (fTitle != "" ) {
0157       h->SetTitle(i0->NamePede(iMulti) + " IOVs" + i0->TitleAdd() + ", " + fTitle + ";IOV;"
0158           + i0->NamePede(iMulti) += i0->UnitPede(iMulti));
0159     } else {     // 'usual' title for drawing parameter values
0160       h->SetTitle((i0->NamePede(iMulti) += " IOVs") += i0->TitleAdd() += ";IOV;"
0161           + i0->NamePede(iMulti) += i0->UnitPede(iMulti));
0162     }
0163     fHistManager->AddHistSame(h, layer, iMulti); // cloned hist for axes
0164     fHistManager->AddObject(multis[iMulti], layer, iMulti, "LP");
0165     // Create legend refering to graphs and add to manager:
0166     Double_t x1, x2, y1, y2; fHistManager->GetLegendX1Y1X2Y2(x1, y1, x2, y2);
0167     TLegend *legend = new TLegend(x1, y1, x2, y2);
0168     legend->SetFillColor(kWhite);
0169     legend->SetTextFont(42);
0170     legend->SetBorderSize(1);
0171     fHistManager->AddLegend(legend, layer, iMulti);
0172     Int_t nGr = 0;
0173     while (TGraph* graph = static_cast<TGraph*>(iter.Next())){
0174       legend->AddEntry(graph, graph->GetTitle(), "lp"); // title set above
0175       this->SetLineMarkerStyle(*graph, nGr++);
0176     }
0177   }
0178   gROOT->SetBatch(isBatch); // reset batch mode
0179
0180   fHistManager->Draw();
0181 }
0182
0183 //////////////////////////////////////////////////////////////////////////
0184 void PlotMillePedeIOV::SetSubDetId(Int_t subDet)
0185 {
0186   for (unsigned int iIov = 0; iIov < fIovs.size(); ++iIov) {
0187     fIovs[iIov]->SetSubDetId(subDet);
0188   }
0189 }
0190
0191 //////////////////////////////////////////////////////////////////////////
0192 void PlotMillePedeIOV::SetSubDetIds(Int_t id1, Int_t id2, Int_t id3, Int_t id4, Int_t id5)
0193 {
0194   for (unsigned int iIov = 0; iIov < fIovs.size(); ++iIov) {
0195     fIovs[iIov]->SetSubDetIds(id1, id2, id3, id4, id5);
0196   }
0197 }
0198
0199 //////////////////////////////////////////////////////////////////////////
0200 void PlotMillePedeIOV::SetAlignableTypeId(Int_t alignableTypeId)
0201 {
0202   for (unsigned int iIov = 0; iIov < fIovs.size(); ++iIov) {
0203     fIovs[iIov]->SetAlignableTypeId(alignableTypeId);
0204   }
0205 }
0206
0207 //////////////////////////////////////////////////////////////////////////
0208 void PlotMillePedeIOV::SetHieraLevel(Int_t hieraLevel)
0209 {
0210   for (unsigned int iIov = 0; iIov < fIovs.size(); ++iIov) {
0211     fIovs[iIov]->SetHieraLevel(hieraLevel);
0212   }
0213 }
0214
0215 //////////////////////////////////////////////////////////////////////////
0216 void PlotMillePedeIOV::SetBowsParameters(bool use)
0217 {
0218   for (unsigned int iIov = 0; iIov < fIovs.size(); ++iIov) {
0219     fIovs[iIov]->SetBowsParameters(use);
0220   }
0221 }
0222
0223 //////////////////////////////////////////////////////////////////////////
0224 void PlotMillePedeIOV::AddAdditionalSel(const TString &xyzrPhiNhit, Float_t min, Float_t max)
0225 {
0226   for (unsigned int iIov = 0; iIov < fIovs.size(); ++iIov) {
0227     fIovs[iIov]->AddAdditionalSel(xyzrPhiNhit, min, max);
0228   }
0229 }
0230
0231 //////////////////////////////////////////////////////////////////////////
0232 void PlotMillePedeIOV::ClearAdditionalSel()
0233 {
0234   for (unsigned int iIov = 0; iIov < fIovs.size(); ++iIov) {
0235     fIovs[iIov]->ClearAdditionalSel();
0236   }
0237 }
0238
0239 //////////////////////////////////////////////////////////////////////////
0240 TString PlotMillePedeIOV::Unique(const char *name) const
0241 {
0242   if (!gROOT->FindObject(name)) return name;
0243
0244   UInt_t i = 1;
0245   while (gROOT->FindObject(Form("%s_%u", name, i))) ++i;
0246
0247   return Form("%s_%u", name, i);
0248 }
0249
0250 //////////////////////////////////////////////////////////////////////////
0251 Int_t PlotMillePedeIOV::PrepareAdd(bool addPlots)
0252 {
0253   if (addPlots) {
0254     return fHistManager->GetNumLayers();
0255   } else {
0256     fHistManager->Clear();
0257     return 0;
0258   }
0259 }
0260
0261 //////////////////////////////////////////////////////////////////////////
0262 template<class T>
0263 void PlotMillePedeIOV::SetLineMarkerStyle(T &object, Int_t num) const
0264 {
0265   // styles start with num=0
0266   // color: use 1-4, 6-9, 41, 46, 49
0267   Int_t colour = num%11 + 1; // i.e. 1-11
0268   if (colour > 4) ++colour; // skip 5=yellow
0269   if (colour == 10) colour = 41;
0270   else if (colour == 11) colour = 46;
0271   else if (colour == 12) colour = 49;
0272   object.SetLineColor(colour);
0273   object.SetMarkerColor(object.GetLineColor());
0274
0275   // style
0276   Int_t marker = 0; // use 14 styles
0277   switch (num%14) {
0278   case 0: marker = kFullCircle;
0279     break;
0280   case 1: marker = kFullSquare;
0281     break;
0282   case 2: marker = kFullTriangleUp;
0283     break;
0284   case 3: marker = kFullTriangleDown;
0285     break;
0286   case 4: marker = kOpenCircle;
0287     break;
0288   case 5: marker = kOpenSquare;
0289     break;
0290   case 6: marker = kOpenTriangleUp;
0291     break;
0292   case 7: marker = kOpenTriangleDown; // 32
0293     break;
0294   case 8: marker = kOpenStar;
0295     break;
0296   case 9: marker = kFullStar;
0297     break;
0298   case 10: marker = kOpenCross;
0299     break;
0300   case 11: marker = kFullCross;
0301     break;
0302   case 12: marker = kOpenDiamond;
0303     break;
0304   case 13: marker = kFullDiamond;
0305     break;
0306   default: marker = kPlus; // should not be reached...
0307   }
0308   object.SetMarkerStyle(marker);
0309
0310   //object.SetLineStyle(num%10 + 1); // use (default) line styles 1-10
0311   object.SetLineStyle(num%4 + 1); // use only line styles 1-4
0312
0313 }
0314
0315 bool PlotMillePedeIOV::ParId::operator< (const ParId& other) const
0316 {
0317   // Sorting needed for use as key in std::map:
0318   // first sort by id_, second objId_, finally par_.
0319   if (id_ < other.id_) return true;
0320   else if (id_ > other.id_) return false;
0321   else { // id the same
0322     if (objId_ < other.objId_) return true;
0323     else if (objId_ > other.objId_) return false;
0324     else { // id and objId the same
0325       if (par_ < other.par_) return true;
0326       // redundant: else if (par_ > other.par_) return false;
0327       else return false; // all are the same!
0328     }
0329   }
0330 }