MillePedeAlignmentAlgorithm/macros/PlotFromMillepedeRes.C

0001 #include "TString.h"
0002 #include "TFile.h"
0003 #include "TTree.h"
0004 #include "TH1F.h"
0005 #include "THStack.h"
0006 #include "TCanvas.h"
0007 #include "TLegend.h"
0008 #include "TStyle.h"
0009 #include "TPad.h"
0010 #include "TGraphErrors.h"
0011 #include <stdio.h>
0012 #include <iostream>
0013 #include <iomanip>
0014
0015 enum {PXB,PXF,TIB,TID,TOB,TEC};
0016 TString subdLabels[6]={"PXB","PXF","TIB","TID","TOB","TEC"};
0017 enum {PULLS, NHITS, PARS, PARSwithERRvsLABEL};
0018
0019 //###########################################################
0020
0021 TString StrPlotType(int plotType){
0022   TString str;
0023   if (plotType==PULLS) return "PULLS";
0024   if (plotType==NHITS) return "NHITS";
0025   if (plotType==PARS) return "PARS";
0026   if (plotType==PARSwithERRvsLABEL) return "PARSwithERRvsLABEL";
0027   return "UNKNOWN";
0028 }//end of StrPlotType
0029
0030 //###########################################################
0031
0032 TString StrPar(int parIndLocal){
0033   TString str;
0034   std::cout<<"parIndLocal="<<parIndLocal<<std::endl;
0035   switch (parIndLocal){
0036     case 1: str="u"; break;
0037     case 2: str="v"; break;
0038     case 3: str="w"; break;
0039     case 4: str="alpha"; break;
0040     case 5: str="beta"; break;
0041     case 6: str="gamma"; break;
0042     case 7: str="def1"; break;
0043     case 8: str="def2"; break;
0044     case 9: str="def3"; break;
0045     default: str="UNKNOWN"; break;
0046   }// end of switch
0047   return str;
0048 }//end of StrPar
0049
0050 //###########################################################
0051
0052 TString StrCutSubd(int isubd){
0053   if (isubd==PXB)  return "label>61 && label<8781";//TPB (PXB)
0054   if (isubd==PXF)  return "label>17541 && label<34561";//TPE (PXF)
0055   if (isubd==TIB)  return "label>37021 && label<79041";//TIB
0056   if (isubd==TID)  return "label>121061 && label<132721";//TID
0057   if (isubd==TOB)  return "label>144401 && label<214301";//TOB
0058   if (isubd==TEC)  return "label>284201 && label<380121";//TEC
0059   return "UNKNOWN";
0060 }//end of StrCutSubd
0061
0062 //###########################################################
0063
0064 void PlotParValVsLabelWithErr(TFile* f, TTree* tr, TString strMillepedeRes, TString strOutdir)
0065 {
0066
0067   f->cd();
0068   TString canvName="c_";
0069   canvName+=strMillepedeRes;
0070   canvName+="_";
0071   canvName+=StrPlotType(PARSwithERRvsLABEL);
0072   canvName.ReplaceAll(".res","");
0073
0074   TCanvas* canv = new TCanvas(canvName,canvName,900,600);
0075   canv->Divide(3,2);
0076
0077   for (int ind=1; ind<=6; ind++){
0078     canv->cd(ind);
0079     TPad* pad = (TPad*)canv->GetPad(ind);
0080     TString strCut="((label%20-1)%9+1)==";
0081     strCut+=ind;
0082     int n = tr->Draw("label%700000:10000*parVal:10000*parErr:0.01*(label%700000)",strCut,"goff");
0083     TGraphErrors *gr = new TGraphErrors(n,tr->GetV1(),tr->GetV2(),tr->GetV4(),tr->GetV3());
0084     gr->SetMarkerStyle(20);
0085     gr->SetLineWidth(2);
0086     for (int i=0; i<n; i++){
0087       std::cout<<tr->GetV1()[i]<<" "<<tr->GetV2()[i]<<"+-"<<tr->GetV3()[i]<<std::endl;
0088     }
0089     gr->SetTitle(StrPar(ind)+TString(", 10000*(par+-err)"));
0090     gr->GetXaxis()->SetTitle("label%700000");
0091     gr->Draw("AP");
0092   }// end of loop over ind
0093   canvName+=".png";
0094   TString saveName=strOutdir+canvName;
0095   canv->SaveAs(saveName);
0096   saveName.ReplaceAll(".png",".pdf");
0097   canv->SaveAs(saveName);
0098
0099 }// end of PlotParValVsLabelWithErr
0100
0101 //###########################################################
0102
0103 void PlotHistsNhitsPerModule(TFile* f, TTree* tr, TString strMillepedeRes, TString strOutdir)
0104 {
0105   TString canvName="c_";
0106   canvName+=strMillepedeRes;
0107   canvName+="_";
0108   canvName+=StrPlotType(NHITS);
0109   canvName.ReplaceAll(".res","");
0110
0111
0112   //enum {PXB,PXF,TIB,TID,TOB,TEC};
0113   int colors[6]={1,2,3,4,6,7};
0114 //  TString labels[6]={"PXB","PXF","TIB","TID","TOB","TEC"};
0115
0116   f->cd();
0117   TCanvas* canv = new TCanvas(canvName,canvName,600,600);
0118   canv->SetLogx();
0119   canv->SetLogy();
0120
0121   for (int ind=1; ind<=1; ind++){
0122     TString strHist = "hNhits_";
0123     strHist+=StrPar(ind);
0124     TString strCut="label<700000 && ((label%20-1)%9+1)==";
0125     strCut+=ind;
0126     TStyle style;
0127     style.SetTitleFontSize(0.2);
0128     THStack *hSt = new THStack("hNhits","# of derivatives (~tracks or hits) per module");
0129     TLegend *leg = new TLegend(0.75,0.65,0.95,0.95);
0130     for (int inv=0; inv<6; inv++){
0131       std::cout<<"- - - - - -"<<std::endl;
0132       std::cout<<subdLabels[inv]<<":"<<std::endl;
0133       std::cout<<StrCutSubd(inv)<<": "<<tr->GetEntries(StrCutSubd(inv))<<" parameters"<<std::endl;
0134       TString strHist1=strHist;
0135       strHist1+=ind;
0136       strHist1+=inv;
0137       TH1F* hValInt = new TH1F(strHist1,strHist1,300,10,15000);
0138       TString strCut1 = strCut+TString(" && ")+StrCutSubd(inv);
0139       tr->Draw(TString("Nhits>>")+strHist1,strCut1,"goff");
0140       std::cout<<"# hits = "<<(int)hValInt->GetMean()<<"+-"<<(int)hValInt->GetRMS()<<std::endl;
0141       hValInt->SetLineColor(1);
0142       hValInt->SetFillColor(colors[inv]);
0143       hValInt->SetLineWidth(2);
0144       hSt->Add(hValInt);
0145       leg->AddEntry(hValInt,subdLabels[inv],"f");
0146       leg->SetFillColor(0);
0147     }
0148     hSt->Draw();
0149     leg->Draw("same");
0150
0151   }//end of loop over ind
0152
0153   canvName+=".png";
0154   TString saveName=strOutdir+canvName;
0155   canv->SaveAs(saveName);
0156   saveName.ReplaceAll(".png",".pdf");
0157   canv->SaveAs(saveName);
0158 }//end of PlotHistsNhitsPerModule
0159
0160 //###########################################################
0161
0162 void PlotPullsDistr(TFile* f, TTree* tr, TString strMillepedeRes, TString strOutdir)
0163 {
0164
0165   TString canvName="c_";
0166   canvName+=strMillepedeRes;
0167   canvName+="_";
0168   canvName+=StrPlotType(PULLS);
0169   canvName.ReplaceAll(".res","");
0170
0171   f->cd();
0172   TCanvas* canv = new TCanvas(canvName,canvName,600,600);
0173 //  tr->Draw("parVal","((parN%20-1)%9+1)==2 && (parVal>-0.4 && parVal<0.4)");// v
0174   canv->Divide(3,3);
0175
0176   for (int parInd=1; parInd<=9; parInd++){
0177     canv->cd(parInd);
0178     TString strCut="((label%20-1)%9+1)==";
0179     strCut+=parInd;
0180     strCut+=" && parErr!=0 && parVal!=0 && label<700000";
0181     TString hName="hPulls_";
0182     hName+= StrPar(parInd);
0183     TH1F* h = new TH1F(hName,hName,100,-3,3);
0184     TString strDraw="parVal/(1.41*parErr)>>";
0185     strDraw+=hName;
0186     tr->Draw(strDraw,strCut,"goff");
0187     h->Draw("EP");
0188   }// end of loop over parInd
0189
0190   canvName+=".png";
0191   TString saveName=strOutdir+canvName;
0192   canv->SaveAs(saveName);
0193   saveName.ReplaceAll(".png",".pdf");
0194   canv->SaveAs(saveName);
0195 }// end of PlotPullsDistr
0196
0197 //###########################################################
0198
0199 void PlotParsDistr(TFile* f, TTree* tr, TString strMillepedeRes, TString strOutdir)
0200 {
0201
0202   for (int isubd=PXB; isubd<=TEC; isubd++){
0203
0204     TString canvName="c_";
0205     canvName+=strMillepedeRes;
0206     canvName+="_";
0207     canvName+=StrPlotType(PARS);
0208     canvName+="_";
0209     canvName+=subdLabels[isubd];
0210     canvName.ReplaceAll(".res","");
0211
0212     f->cd();
0213     TCanvas* canv = new TCanvas(canvName,canvName,600,600);
0214     canv->Divide(3,3);
0215
0216     for (int parInd=1; parInd<=9; parInd++){
0217       canv->cd(parInd);
0218       TString strCut="((label%20-1)%9+1)==";
0219       strCut+=parInd;
0220       strCut+=" && label<700000 && ";
0221       strCut+=StrCutSubd(isubd);
0222       TString hName="hPars_";
0223       hName+=subdLabels[isubd];
0224       hName+="_";
0225       hName+= StrPar(parInd);
0226
0227       TTree* trCut = tr->CopyTree(strCut);
0228       float up = trCut->GetMaximum("parVal");
0229       float low = trCut->GetMinimum("parVal");
0230       std::cout<<"low="<<low<<", up="<<up<<", nent="<<trCut->GetEntries()<<std::endl;
0231       TH1F* h = new TH1F(hName,hName,100,10000*low,10000*up);
0232       TString strDraw="10000*parVal>>";
0233       strDraw+=hName;
0234       trCut->Draw(strDraw,strCut,"goff");
0235       h->SetMarkerStyle(2);
0236       h->Draw("EP");
0237     }// end of loop over parInd
0238
0239     canvName+=".png";
0240     TString saveName=strOutdir+canvName;
0241     canv->SaveAs(saveName);
0242     saveName.ReplaceAll(".png",".pdf");
0243     canv->SaveAs(saveName);
0244
0245   }//end of loop over isubd
0246
0247 }// end of PlotParsDistr
0248
0249 //###########################################################
0250
0251 void PlotFromMillepedeRes(TString strMillepedeRes, TString strOutdir, TString strVars, int plotType)
0252 {
0253   // strPar = "u", "v", "w", "alpha", "beta", "gamma", "def1", "def2", "def3"
0254   TFile* f = new TFile(strOutdir+TString("fOut.root"),"recreate");
0255   TTree* tr = new TTree("tr","tr");
0256   tr->ReadFile(strMillepedeRes,strVars);
0257
0258   if (plotType==PARSwithERRvsLABEL)
0259     PlotParValVsLabelWithErr(f, tr, strMillepedeRes, strOutdir);
0260
0261   if (plotType==NHITS)
0262     PlotHistsNhitsPerModule(f, tr, strMillepedeRes, strOutdir);
0263
0264   if (plotType==PULLS)
0265     PlotPullsDistr(f, tr, strMillepedeRes, strOutdir);
0266
0267   if (plotType==PARS)
0268     PlotParsDistr(f, tr, strMillepedeRes, strOutdir);
0269
0270 }// end of PlotPars