MillePedeAlignmentAlgorithm/macros/createChi2ndfplot.C

0001 // ROOT script to read the output produced by mps_parse_pedechi2hist.pl
0002 // Author : Joerg Behr
0003
0004
0005 // Usage:
0006
0007 // Start ROOT and compile the script:
0008 // root [0] .L createChi2ndfplot.C+
0009
0010 // Call the plotting function and provide the name of inputfile produced by mps_parse_pedechi2hist.pl:
0011 // createChi2ndfplot("chi2pedehis.txt");
0012 //
0013 // Or on the command line:
0014 // root -l -x -b -q 'createChi2ndfplot.C+(\"chi2pedehis.txt\")'
0015
0016 #include <fstream>
0017 #include <vector>
0018 #include <utility>
0019 #include <iostream>
0020 #include <TROOT.h>
0021 #include <TFile.h>
0022 #include <TDirectory.h>
0023 #include <TError.h>
0024 #include <TH1.h>
0025 #include <TH1D.h>
0026 #include <TGraph.h>
0027 #include <TGraphErrors.h>
0028 #include <TCanvas.h>
0029 #include <TMath.h>
0030 #include <TPaveText.h>
0031 #include <map>
0032 #include <list>
0033 #include <string>
0034 #include <TStyle.h>
0035 #include <TGaxis.h>
0036 #include <TLegend.h>
0037
0038 class histdata
0039 {
0040 public:
0041   histdata(
0042            const std::string name2,
0043            const int number2,
0044            const double value2
0045            ) : name(name2),
0046                number(number2),
0047                value(value2)
0048   {
0049   }
0050   //needed for sort.
0051   bool operator < (const histdata& rhs)
0052   {
0053     return number < rhs.number;
0054   }
0055
0056   const std::string name; //name of configuration to which the Mille binary belongs
0057   const int number; //The binary number (milleBinary334.dat)
0058   const double value; // <chi^2/ndf>
0059 };
0060
0061 //Some configurations for the pad in which the plot is drawn.
0062 void pad_cfg(
0063              TH1D *h
0064              )
0065 {
0066   const Int_t maxd = 2;
0067   const Short_t borderMode = 0;
0068   const Double_t leftMargin = 0.16;
0069   const Double_t rightMargin = 0.05;
0070   //const Double_t bottomMargin = 0.12;
0071   gPad->SetFrameBorderMode(0);
0072
0073   gStyle->SetPalette(1,0);
0074   gPad->SetTickx(1);
0075   gPad->SetTicky(1);
0076   gPad->SetBorderMode(borderMode);
0077   gPad->SetLeftMargin(leftMargin);
0078   gPad->SetRightMargin(rightMargin);
0079   gPad->SetBottomMargin(0.15);
0080   gPad->SetTopMargin(0.07);
0081   TGaxis::SetMaxDigits(maxd);
0082   gStyle->SetTitleX(.0);
0083   gStyle->SetTitleW(.95);
0084   gStyle->SetTitleFontSize(.06);
0085   gStyle->SetTitleStyle(0);
0086   gStyle->SetTitleBorderSize(0);
0087   h->GetXaxis()->SetNoExponent();
0088   h->GetYaxis()->SetNoExponent();
0089   h->GetXaxis()->SetLabelSize(0.05);
0090   h->GetXaxis()->SetTitleSize(0.07); //1.0
0091   h->GetXaxis()->SetTitleOffset(0.83); //1.0
0092   h->GetYaxis()->SetLabelSize(0.05);
0093   h->GetYaxis()->SetTitleOffset(0.74); //1.2 //0.7
0094   h->GetYaxis()->SetTitleSize(0.07);
0095   h->SetStats(kFALSE);
0096   h->GetYaxis()->SetNdivisions(505);
0097   h->GetXaxis()->SetNdivisions(505);
0098
0099   //the margins
0100   gPad->SetLeftMargin(0.1162393);
0101   gPad->SetRightMargin(0.006837607);
0102   gPad->SetTopMargin(0.05347594);
0103   gPad->SetBottomMargin(0.1354724);
0104
0105
0106
0107 }
0108
0109 //This methods checks whether a histogram is existing, books a new histogram if nessesary, and fills the histogram.
0110 void fillhisto(std::map<std::string,TH1D*> &histomap, const std::string name, const int value, const double wtx, const int nbins, const double lower, const double upper)
0111 {
0112   std::map<std::string,TH1D*>::iterator it = histomap.find(name); //histogram existing?
0113   if(it != histomap.end())
0114     it->second->Fill(static_cast<double>(value),wtx);
0115   else
0116     {
0117       histomap[name] = new TH1D(TString(name),TString(name),nbins,lower, upper); //book new histogram
0118       histomap[name]->Fill(static_cast<double>(value),wtx);
0119     }
0120 }
0121
0122 //the main (and only) plotting routine.
0123 void createChi2ndfplot(const char *txtFile)
0124 {
0125   std::ifstream theStream(txtFile, ios::in);
0126   if (!theStream.is_open()) {
0127     ::Error("createChi2ndfplot", "file %s could not be opened", txtFile);
0128   } else {
0129     //the collection of histograms.
0130     std::map<std::string,TH1D*> histomap;
0131
0132     //The data read from file txtFile
0133     std::list<histdata> h;
0134
0135     //Variables needed for reading
0136     std::string name("");
0137     int number = 0;
0138     double value = 0.0;
0139
0140     //Read the data
0141     while(
0142           theStream >> name
0143           >> number
0144           >> value
0145           )
0146       {
0147         h.push_back(histdata(name,number,value));
0148       }
0149     theStream.close();
0150
0151     //If data was read
0152     if(h.size() > 0)
0153       {
0154         //Sort the list of data, because the output of mps_parse_pedechi2hist.pl can be arbitrary ordered.
0155         h.sort();
0156
0157         double min = 10000.0, max = 0.0;
0158         int bincounter = 0; //The bin number (If binaries are missing in the pede job, then this quantity differs from the pure binary number)
0159         for(std::list<histdata>::const_iterator it = h.begin(); it != h.end(); it++)
0160           {
0161             bincounter++;
0162
0163             //fill and book histogram
0164             fillhisto(histomap, it->name, bincounter, it->value, static_cast<int>(h.size()), 1.0, static_cast<double>(h.size()+1));
0165
0166             //find minimal and maximal value of <chi2/ndf> used for axis limits.
0167             it->value > max ? max = it->value : 0;
0168             it->value < min ? min = it->value : 0;
0169           }
0170
0171         TCanvas* c = new TCanvas("chi2ndfperbinary","chi2ndfperbinary", 200, 10, 900, 1200);
0172
0173         gStyle->SetPadBorderMode(0);
0174         gStyle->SetOptStat(0);
0175         c->SetFillColor(10);
0176         c->Divide(1,2);
0177         c->cd(1);
0178
0179
0180         //haxis is only used for the axis.
0181         TH1D *haxis = new TH1D("haxis","; binary number; <#chi^{2}/ndf>", 1, 0.0, static_cast<double>(h.size()+1));
0182
0183         //setup the pad.
0184         pad_cfg(haxis);
0185
0186
0187         haxis->SetMaximum(max*1.2);
0188         haxis->SetMinimum(min*0.8);
0189
0190         haxis->DrawCopy("axis");
0191
0192         //The legend which will be drawn in the second pad.
0193         TLegend *leg = new TLegend(0.01826713,0.003922913,0.9976767,0.9979814);
0194         leg->SetBorderSize(0);
0195         leg->SetFillColor(19);
0196         leg->SetFillStyle(0);
0197
0198         //some usable colors
0199         const int ncolors = 11;
0200         const int kcolors[ncolors] = {kRed, kBlue, kBlack, kOrange, kGreen, kMagenta, kBlue+1, kBlue-1, kBlue+2, kBlue-2, kMagenta-2};
0201
0202         int histocounter = 0;
0203         for( std::map<std::string,TH1D*>::iterator it = histomap.begin(); it != histomap.end(); it++)
0204           {
0205             int color = 1;
0206             //assign a color.
0207             if(histocounter < ncolors)
0208               {
0209                 color = kcolors[histocounter];
0210                 histocounter++;
0211               }
0212             else
0213               {
0214                 color = histocounter;
0215                 histocounter++;
0216               }
0217
0218             //optical appearance
0219             it->second->SetMarkerColor(color);
0220             it->second->SetMarkerSize(1);
0221             it->second->SetMarkerStyle(20+histocounter);
0222             it->second->DrawCopy("p same");
0223             leg->AddEntry(it->second, TString(it->first), "p");
0224           }
0225
0226
0227
0228         //go the the second pad where the legend is depicted.
0229         c->cd(2);
0230         pad_cfg(haxis);
0231         leg->Draw();
0232
0233
0234         c->Print("chi2ndfperbinary.pdf");
0235         c->Print("chi2ndfperbinary.C");
0236
0237
0238         //clean up
0239         delete c;
0240         delete leg;
0241         delete haxis;
0242         for( std::map<std::string,TH1D*>::iterator it = histomap.begin(); it != histomap.end(); it++)
0243           {
0244             delete it->second;
0245           }
0246       }
0247   }
0248 }