Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2024-09-27 07:03:27 

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2023 Christopher Dilks, Connor Pecar, Duane Byer
0003 
0004 #include "PostProcessor.h"
0005 
0006 ClassImp(PostProcessor)
0007 
0008 using std::cout;
0009 using std::cerr;
0010 using std::endl;
0011 
0012 // constructor
0013 PostProcessor::PostProcessor(
0014   TString infileN_,
0015   TString outfileN_
0016 )
0017   : infileN(infileN_)
0018   , outfileN(outfileN_)
0019 {
0020 
0021   // settings
0022   // - summary canvas formatting
0023   summaryColor[0] = kRed;     summaryStyle[0] = kFullCircle;
0024   summaryColor[1] = kGreen+1; summaryStyle[1] = kFullTriangleUp;
0025   summaryColor[2] = kBlue;    summaryStyle[2] = kFullTriangleDown;
0026   // - general
0027   gStyle->SetOptStat(0);
0028   cout << std::fixed; // fixed float precision
0029   cout << std::setprecision(3);
0030 
0031   // set up input and output files
0032   cout << "-- postprocess histograms from " << infileN << endl;
0033   if(outfileN=="") outfileN = infileN; // default
0034   outfileN(TRegexp("\\.root$")) = "";
0035   pngDir = outfileN+".images";
0036   outfileN += ".canvas.root";
0037   outfile = new TFile(outfileN,"RECREATE");
0038   infile = new TFile(infileN,"READ");
0039   gROOT->ProcessLine(".! mkdir -p "+pngDir);
0040 
0041   // rebuild DAGs
0042   HD = new HistosDAG();
0043   HD->BuildFromFile(infile);
0044 
0045   // initialize algorithm-specific vars
0046   this->ResetVars();
0047 
0048 };
0049 
0050 
0051 //=========================================================================
0052 // cleanup and close open files and streams
0053 void PostProcessor::Finish() {
0054   infile->Close();
0055   outfile->Close();
0056   cout << outfileN << " written." << endl;
0057   cout << pngDir << "/ images written." << endl;
0058 };
0059 
0060 
0061 //=========================================================================
0062 // reset global variables, such as summary canvases
0063 void PostProcessor::ResetVars() {
0064   nsum=0;
0065   ndump=0;
0066   summaryCanvMap.clear();
0067   varList.clear();
0068 };
0069 
0070 
0071 //=========================================================================
0072 /* ALGORITHM: dump histogram counts
0073  * - output will be appended to `countsFiles`
0074  * - so far only implemented for 1D
0075  */
0076 void PostProcessor::DumpHist(TString datFile, TString histSet, TString varName) {
0077   cout << "dump " << histSet << " : " << varName << " to " << datFile << endl;
0078   Histos *H = (Histos*) infile->Get(histSet);
0079   TH1 *hist = H->Hist(varName);
0080   if(hist->GetDimension()>1) return;
0081   TString histTformatted = hist->GetTitle();
0082   histTformatted.ReplaceAll("#in"," in ");
0083   histTformatted.ReplaceAll("#pm","+-");
0084 
0085   gSystem->RedirectOutput(datFile,"a");
0086   cout << endl;
0087   cout << "Histogram: " << histTformatted << endl;
0088 
0089   gSystem->RedirectOutput("tempo","w");
0090   cout << "bin"
0091        << " " << varName << "_min"
0092        << " " << varName << "_max"
0093        << " " << "counts"
0094        << " " << "error"
0095        << " " << endl;
0096   for(int b=1; b<=hist->GetNbinsX(); b++) {
0097     cout << b
0098          << " " << hist->GetBinLowEdge(b)
0099          << " " << hist->GetBinLowEdge(b+1)
0100          << " " << hist->GetBinContent(b)
0101          << " " << hist->GetBinError(b)
0102          << endl;
0103   };
0104 
0105   gSystem->RedirectOutput(datFile,"a");
0106   gROOT->ProcessLine(".! cat tempo | column -t");
0107   gSystem->RedirectOutput(0);
0108   gROOT->ProcessLine(".! rm tempo");
0109 
0110 };
0111 
0112 //=========================================================================
0113 
0114 /* ALGORITHM: dump averages of all histograms from `histSet` to the file `datFile`
0115  * - it is best to call this function in a loop, this function will only dump
0116  *   one line of information; if it's the first time you called it, it will
0117  *   dump header information as well
0118  * - use `cutName` to specify a cut defintion you want to print; e.g., loop
0119  *   over x bins and set the `cutName` to the x cut, to include columns for x
0120  *   bin boundaries
0121  * - when done looping, call `FinishDumpAve(datFile)`
0122  */
0123 void PostProcessor::DumpAve(TString datFile, Histos *H, TString cutName) {
0124   cout << "dump averages from " << H->GetSetName()
0125        << " to " << datFile << endl;
0126 
0127   // get associated cut, and print header info if it's the first time
0128   dumpCut = nullptr;
0129   if(ndump==0) {
0130     gSystem->RedirectOutput(datFile,"a");
0131     cout << endl << "Histogram Set:";
0132   };
0133   for(CutDef *cut : H->CutDefList) {
0134     if(cutName.CompareTo(cut->GetVarName(),TString::kIgnoreCase)==0) {
0135       dumpCut = cut;
0136     }
0137     else if(ndump==0) {
0138       TString cutT = cut->GetCutTitle();
0139       cutT.ReplaceAll("#in"," in ");
0140       cutT.ReplaceAll("#pm","+-");
0141       cout << "   " + cutT;
0142     };
0143   };
0144   if(ndump==0) {
0145     cout << endl;
0146     cout << "Averages in bins of " << dumpCut->GetVarTitle() << ":" << endl;
0147     cout << "------------------------" << endl;
0148     gSystem->RedirectOutput(0);
0149   };
0150   if(dumpCut==nullptr) {
0151     cerr << "ERROR: in DumpAve, cannot find cut def" << endl;
0152     return;
0153   };
0154 
0155   // start temporary file for output; this is so you can use `Columnify` later
0156   // for pretty print to full `datFile`
0157   gSystem->RedirectOutput(datFile+".tmp",ndump>0?"a":"w");
0158 
0159   // if this is the first time, populate ordered list of histograms; we keep a
0160   // local list, in case you call this function on a `Histos` with different
0161   // ordering; we want to maintain order for printout
0162   if(ndump==0) {
0163     // header for bin columns
0164     cout << "bin"
0165          << " " << dumpCut->GetVarTitle() << "_min"
0166          << " " << dumpCut->GetVarTitle() << "_max"
0167          << " counts";
0168     // loop over 1D histograms, build `varList`,
0169     // and print header for ave columns
0170     for(TString varName : H->VarNameList) {
0171       if(H->Hist(varName)->GetDimension()==1) {
0172         if(varName.Contains("xsec",TString::kIgnoreCase)) continue;
0173         if(varName.Contains("Res",TString::kIgnoreCase)) continue;
0174         if(varName.Contains("RvG",TString::kIgnoreCase)) continue;
0175         if(H->Hist(varName)->GetEntries()==0) continue;
0176         varList.push_back(varName);
0177         cout << " <" + varName + ">";
0178       };
0179     };
0180     // header newline
0181     cout << endl;
0182   };
0183 
0184   // loop over local list of variables, and print their averages
0185   TH1 *hist;
0186   cout << ndump
0187        << " " << dumpCut->GetMin()
0188        << " " << dumpCut->GetMax();
0189   Bool_t first=true;
0190   Double_t counts;
0191   for(TString varName : varList) {
0192     hist = H->Hist(varName);
0193     if(first) {
0194       counts = hist->GetEntries();
0195       cout << " " << counts;
0196       first=false;
0197     }
0198     else if(hist->GetEntries()!=counts) {
0199       cerr << "WARNING: mismatch counts in " << varName << " histogram" << endl;
0200     };
0201     cout << " " << hist->GetMean();
0202   };
0203   cout << endl;
0204   gSystem->RedirectOutput(0);
0205 
0206 
0207   // iterate row counter
0208   ndump++;
0209 
0210 };
0211 void PostProcessor::FinishDumpAve(TString datFile) {
0212   this->Columnify(datFile+".tmp",datFile);
0213   gROOT->ProcessLine(".! rm "+datFile+".tmp");
0214   this->ResetVars();
0215 };
0216 
0217 
0218 //=========================================================================
0219 /* ALGORITHM: draw a single histogram to a canvas, and write it
0220  * - `histName` is the name of the histogram given in Histos
0221  * - `drawFormat` is the formatting string passed to TH1::Draw()
0222  * - `profileAxis` will draw a TProfile on the specified axis, for 2D dists
0223  *   - 0=disabled(default), 1=x-axis, 2=y-axis
0224  * - `profileOnly` draw only the TProfile, for 2D dists
0225  */
0226 void PostProcessor::DrawSingle(Histos *H, TString histName, TString drawFormat, Int_t profileAxis, Bool_t profileOnly) {
0227   cout << "draw single plot " << histName << "..." << endl;
0228   TH1 *hist = H->Hist(histName);
0229   if(hist==nullptr) {
0230     cerr << "ERROR: cannot find histogram " << histName << endl;
0231     return;
0232   };
0233   TString canvN = "canv_"+histName+"___"+H->GetSetName();
0234   TCanvas *canv = new TCanvas(canvN,canvN,dimx,dimy);
0235   hist->Draw(drawFormat);
0236   if(hist->GetMinimum()>=0 && hist->GetDimension()==1 && H->GetHistConfig(histName)->logy==false)
0237     hist->GetYaxis()->SetRangeUser(0,hist->GetMaximum()*1.1); // do not suppress zero
0238 
0239   canv->SetGrid(1,1);
0240   canv->SetLogx(H->GetHistConfig(histName)->logx);
0241   canv->SetLogy(H->GetHistConfig(histName)->logy);
0242   canv->SetLogz(H->GetHistConfig(histName)->logz);
0243   canv->SetBottomMargin(0.15);
0244   canv->SetLeftMargin(0.15);
0245   canv->SetRightMargin(0.15);
0246 
0247   // profile for 2D plot
0248   if(profileAxis>0 && hist->GetDimension()==2) {
0249     TProfile *prof;
0250     switch(profileAxis) {
0251       case 1: prof = ((TH2*)hist)->ProfileX(); break;
0252       case 2: prof = ((TH2*)hist)->ProfileY(); break;
0253     };
0254     prof->SetLineColor(kBlack);
0255     prof->SetLineWidth(3);
0256     if(profileOnly) {
0257       prof->GetXaxis()->SetRangeUser(hist->GetXaxis()->GetXmin(),hist->GetXaxis()->GetXmax());
0258       prof->GetYaxis()->SetRangeUser(hist->GetYaxis()->GetXmin(),hist->GetYaxis()->GetXmax());
0259       prof->Draw();
0260     } else prof->Draw("same");
0261     outfile->cd("/");
0262     prof->Write();
0263   }
0264 
0265   // profile for 3D plot
0266   if(profileAxis>0 && hist->GetDimension()==3) {
0267     TProfile2D *prof = ((TH3*)hist)->Project3DProfile("yx"); // FIXME: generalize for other axes
0268     prof->GetXaxis()->SetRangeUser(hist->GetXaxis()->GetXmin(),hist->GetXaxis()->GetXmax());
0269     prof->GetYaxis()->SetRangeUser(hist->GetYaxis()->GetXmin(),hist->GetYaxis()->GetXmax());
0270     prof->GetZaxis()->SetRangeUser(hist->GetZaxis()->GetXmin(),hist->GetZaxis()->GetXmax());
0271     prof->GetXaxis()->SetTitle(hist->GetXaxis()->GetTitle());
0272     prof->GetYaxis()->SetTitle(hist->GetYaxis()->GetTitle());
0273     prof->GetZaxis()->SetTitle(hist->GetZaxis()->GetTitle());
0274     prof->Draw(drawFormat);
0275     outfile->cd("/");
0276     prof->Write();
0277   }
0278 
0279   canv->Print(pngDir+"/"+canvN+".png");
0280   outfile->cd("/");
0281   canv->Write();
0282 };
0283 
0284 // OLD VERSION: 
0285 void PostProcessor::DrawSingle(TString histSet, TString histName) {
0286 
0287   Histos *H = (Histos*) infile->Get(histSet);
0288   TH1 *hist = H->Hist(histName,true);
0289   Hist4D *hist4 = H->Hist4(histName,true);
0290 
0291   if (hist != nullptr) {
0292     TString canvN = "canv_"+histName+"_"+H->GetSetName();
0293     TCanvas *canv = new TCanvas(canvN,canvN, dimx, dimy);
0294 
0295     hist->SetLineColor(kBlack);
0296     hist->SetMarkerColor(kBlack);
0297     hist->SetMarkerStyle(kFullCircle);
0298     hist->SetMarkerSize(1.0);
0299     hist->SetLineWidth(2);
0300     hist->GetXaxis()->SetLabelSize(0.06);
0301     hist->GetYaxis()->SetLabelSize(0.06);
0302     hist->GetXaxis()->SetTitleSize(0.06);
0303     hist->GetYaxis()->SetTitleSize(0.06);
0304     hist->GetXaxis()->SetTitleOffset(1.2);
0305     
0306     // determine draw type (TODO: probably could be generalized somehow)
0307     TString drawStr = "";
0308     switch(hist->GetDimension()) {
0309       case 1:
0310         drawStr = "EX0 P";
0311         if(histName=="Q_xsec") {
0312           hist->GetXaxis()->SetRangeUser(1,10);
0313           hist->GetYaxis()->SetRangeUser(1e-6,1);
0314           drawStr = "E P";
0315         };
0316         break;
0317       case 2:
0318         drawStr = "COLZ";
0319         break;
0320       case 3:
0321         drawStr = "BOX";
0322         break;
0323     };
0324 
0325     hist->Draw(drawStr);
0326 
0327     canv->SetGrid(1,1);
0328     canv->SetLogx(H->GetHistConfig(histName)->logx);
0329     canv->SetLogy(H->GetHistConfig(histName)->logy);
0330     canv->SetLogz(H->GetHistConfig(histName)->logz);
0331     canv->SetBottomMargin(0.15);
0332     canv->SetLeftMargin(0.15);
0333     canv->SetRightMargin(0.15);
0334     canv->Print(pngDir+"/"+canvN+".png");
0335     outfile->cd("/");
0336     canv->Write();
0337     outfile->cd("/");
0338   } else if (hist4 != nullptr) {
0339     TString canvN = "canv_"+histName+"_"+H->GetSetName();
0340     TCanvas *canv = new TCanvas(canvN,canvN, dimx, dimy);
0341 
0342     hist4->GetWaxis()->SetLabelSize(0.06);
0343     hist4->GetXaxis()->SetLabelSize(0.06);
0344     hist4->GetWaxis()->SetTitleSize(0.06);
0345     hist4->GetXaxis()->SetTitleSize(0.06);
0346     hist4->GetWaxis()->SetTitleOffset(1.2);
0347     
0348     //canv->SetGrid(1,1);
0349     canv->SetLogx(H->GetHist4Config(histName)->logx);
0350     canv->SetLogy(H->GetHist4Config(histName)->logy);
0351     canv->SetLogz(H->GetHist4Config(histName)->logz);
0352     canv->SetBottomMargin(0.15);
0353     canv->SetLeftMargin(0.15);
0354     canv->SetRightMargin(0.15);
0355     hist4->Draw();
0356 
0357     canv->Print(pngDir+"/"+canvN+".png");
0358     outfile->cd("/");
0359     canv->Write();
0360     outfile->cd("/");
0361   } else {
0362     cerr << "Couldn't find histogram " << histName << std::endl;
0363   }
0364 
0365   /* // deprecated, for combining single plots; TODO if needed, make a separate method
0366   TH1 *histClone = (TH1*) hist->Clone();
0367   histClone->SetLineColor  (nsum<nsumMax?summaryColor[nsum]:kBlack);
0368   histClone->SetMarkerColor(nsum<nsumMax?summaryColor[nsum]:kBlack);
0369   histClone->SetMarkerStyle(nsum<nsumMax?summaryStyle[nsum]:kFullCircle);
0370 
0371   if(nsum==0) {
0372     summaryCanv = new TCanvas(
0373         "summaryCanv_"+histName,
0374         "summaryCanv_"+histName,
0375         dimx,dimy
0376         );
0377     summaryCanv->SetGrid(1,1);
0378     summaryCanv->SetLogx(H->GetHistConfig(histName)->logx);
0379     summaryCanv->SetLogy(H->GetHistConfig(histName)->logy);
0380     summaryCanv->SetLogz(H->GetHistConfig(histName)->logz);
0381     summaryCanv->SetBottomMargin(0.15);
0382     summaryCanv->SetLeftMargin(0.15);
0383   };
0384   summaryCanv->cd();
0385   histClone->Draw(drawStr+(nsum>0?" SAME":""));
0386   nsum++;
0387   */
0388 };
0389 
0390 //=========================================================================
0391 /* ALGORITHM: draw histograms from different bins in their respective bins
0392 on axis of bin variables, e.g. Q2 vs x.
0393 */
0394 // -- process one histArr
0395 void PostProcessor::DrawInBins(
0396     TString outName,    
0397     std::vector<std::vector<Histos*>>& histArr,
0398     TString histName,
0399     TString var1name, int nvar1, double var1low, double var1high, bool var1log,
0400     TString var2name, int nvar2, double var2low, double var2high, bool var2log,
0401     bool intgrid1, bool intgrid2, // grid option for small plots
0402     bool renormalize // if true, normalize histograms
0403 ){
0404   std::vector<std::vector<std::vector<Histos*>>> histArrList;
0405   histArrList.push_back(histArr); // build list of histArrs with the one histArr
0406   this->DrawInBins( // call the main algo below, which uses histArrList
0407       outName,    
0408       histArrList,
0409       histName,
0410       var1name, nvar1, var1low, var1high, var1log,
0411       var2name, nvar2, var2low, var2high, var2log,
0412       intgrid1, intgrid2,
0413       renormalize
0414       );
0415 };
0416 // -- process list of histArrs
0417 void PostProcessor::DrawInBins(
0418     TString outName,    
0419     std::vector<std::vector<std::vector<Histos*>>>& histArrList,
0420     TString histName,
0421     TString var1name, int nvar1, double var1low, double var1high, bool var1log,
0422     TString var2name, int nvar2, double var2low, double var2high, bool var2log,
0423     bool intgrid1, bool intgrid2, // grid option for small plots
0424     bool renormalize // if true, normalize histograms (useful for drawing 1D hists with "SAME")
0425 ){
0426   // default values set for nvar1==nvar2
0427   int canvx = 1400;
0428   int canvy = 1200;
0429   double botmargin = 0.2;
0430   double leftmargin = 0.2;
0431   double xaxisy = 0.04;
0432   double xaxisx1 = 0.08;
0433   double xaxisx2 = 0.97;
0434   double yaxisx = 0.04;
0435   double yaxisy1 = 0.085;
0436   double yaxisy2 = 0.97;
0437   if(nvar1 > nvar2){
0438     // different canvas sizing/axis position for unequal binning
0439     canvx = 2200;
0440     canvy = 1400;
0441     xaxisx1 = 0.075;
0442     xaxisx2 = 0.975;
0443     yaxisy1 = 0.08;
0444   };
0445 
0446   TString canvN = "canv_"+outName+"_"+histName;
0447   TCanvas *canv = new TCanvas(canvN,canvN, canvx, canvy);
0448   TPad *mainpad = new TPad("mainpad", "mainpad", 0.07, 0.07, 0.98, 0.98);
0449 
0450   TLegend *leg;
0451   if(legendLabels.size()>0) {
0452     double legX0 = 0.1;
0453     double legY0 = 0.95;
0454     double legWidth = 0.15;
0455     double legRowHeight = 0.05;
0456     leg = new TLegend(
0457         legX0,
0458         legY0,
0459         legX0 + legWidth,
0460         legY0 - legRowHeight*legendLabels.size()
0461         );
0462   };
0463 
0464   mainpad->SetFillStyle(4000);
0465   mainpad->Divide(nvar1,nvar2,0,0);
0466   mainpad->Draw();
0467   TLine * lDIRC = new TLine(6,-1,6,1);
0468   TLine * lDIRClow = new TLine(0.5,-1,0.5,1);
0469   TLine * lmRICH = new TLine(2,-1,2,-4);
0470   TLine * lDRICH = new TLine(2.5,1,2.5,4);
0471   lDIRC->SetLineColor(kRed);
0472   lDIRClow->SetLineColor(kRed);
0473   lmRICH->SetLineColor(kRed);
0474   lDRICH->SetLineColor(kRed);
0475   // TH1* histArray[nvar1][nvar2]; // TODO: re-enable
0476   int drawpid = 0;
0477   outfile->cd("/");
0478   canv->Write();
0479 
0480   // get histograms from Histos 2D vector
0481   Histos *H;
0482   for(int i = 0; i < nvar1; i++){
0483     for(int j = 0; j < nvar2; j++){
0484       int count = 0;
0485       int dims;
0486       for(auto histArr : histArrList) {
0487         H = histArr[i][j];
0488         TH1 *hist = H->Hist(histName);
0489         dims = hist->GetDimension();
0490         // histArray[i][j] = hist;
0491         hist->SetTitle("");
0492         //hist->GetXaxis()->SetTitle("");
0493         //hist->GetYaxis()->SetTitle("");
0494         //hist->GetXaxis()->SetLabelSize(0);
0495         //hist->GetYaxis()->SetLabelSize(0);
0496 
0497         hist->SetLineWidth(3);
0498         switch(count) {
0499           case 0:
0500             hist->SetLineColor(kBlack);
0501             hist->SetFillColor(kGray);
0502             break;
0503           case 1:
0504             hist->SetLineColor(kAzure);
0505             hist->SetLineStyle(1);
0506             break;
0507           case 2:
0508             hist->SetLineColor(kRed);
0509             hist->SetLineStyle(7);
0510             break;
0511           case 3:
0512             hist->SetLineColor(kGreen+2);
0513             hist->SetLineStyle(9);
0514             break;
0515         }
0516 
0517         if(legendLabels.size()>0 && i==0 && j==0) {
0518           try { leg->AddEntry(hist,legendLabels[count],"LF"); }
0519           catch(const std::out_of_range &e) {
0520             cerr << "ERROR: legendLabels not filled correctly" << endl;
0521           };
0522         };
0523 
0524         if(count==0) {
0525           mainpad->cd((nvar2-j-1)*nvar1 + i + 1);
0526           gPad->SetLogx(H->GetHistConfig(histName)->logx);
0527           gPad->SetLogy(H->GetHistConfig(histName)->logy);
0528           gPad->SetLogz(H->GetHistConfig(histName)->logz);
0529           gPad->SetGridy(intgrid2);
0530           gPad->SetGridx(intgrid1);
0531         }
0532 
0533         // renormalize
0534         if(renormalize) {
0535           // hist->Scale(1/hist->GetMaximum());
0536           // hist->Scale(1/hist->GetEntries());
0537           hist->Scale(1/hist->Integral());
0538         }
0539 
0540         TString drawStr = "";
0541         switch(dims) {
0542           case 1:
0543             drawStr = "HIST MIN0";
0544             break;
0545           case 2:
0546             drawStr = "COLZ";
0547             break;
0548           case 3:
0549             drawStr = "BOX";
0550             break;
0551         };
0552         if(count>0) drawStr += " SAME";
0553 
0554         if( hist->GetEntries() > 0 ) {  
0555           hist->Draw(drawStr);
0556           if(drawpid){
0557             lDIRClow->Draw();
0558             lDIRC->Draw();
0559             lmRICH->Draw();
0560             lDRICH->Draw();
0561           }
0562           hist->GetXaxis()->SetLabelSize(0.04);
0563           hist->GetYaxis()->SetLabelSize(0.04);
0564           hist->GetXaxis()->SetTitleSize(0.05);
0565           hist->GetYaxis()->SetTitleSize(0.05);
0566           hist->GetXaxis()->SetTitleOffset(0.9);
0567           hist->GetXaxis()->SetLabelOffset(0.0005);
0568           if(dims==1) {
0569             hist->GetYaxis()->SetLabelSize(0.00); // suppress y-axis labels (since each subplot has its own scale)
0570           }
0571         }
0572         count++;
0573       }; // end for(histArrList)
0574 
0575       // some formatting for 1D plots drawn with "SAME"
0576       if(dims==1) {
0577         UnzoomVertical(gPad,"",true,H->GetHistConfig(histName)->logy); // allow all "SAME" plots to be visible; 3rd arg forces min to be 0
0578       }
0579     };
0580   };
0581   canv->cd();
0582 
0583   TPad *newpad1 = new TPad("newpad1","full pad",0,0,1,1);
0584   TPad *newpad2 = new TPad("newpad2","full pad",0,0,1,1);
0585   newpad1->SetFillStyle(4000);
0586   newpad1->Draw();
0587   newpad2->SetFillStyle(4000);
0588   newpad2->Draw();
0589 
0590   TString xopt, yopt;
0591   if(var1log) xopt = "GS";
0592   else xopt = "S";
0593   if(var2log) yopt = "GS";
0594   else yopt = "S";
0595 
0596   TGaxis *xaxis = new TGaxis(xaxisx1,xaxisy,xaxisx2,xaxisy,var1low,var1high,510,xopt);
0597   TGaxis *yaxis = new TGaxis(yaxisx,yaxisy1,yaxisx,yaxisy2,var2low,var2high,510,yopt);
0598   xaxis->SetTitle(var1name);
0599   xaxis->SetName("xaxis");
0600   xaxis->SetTitleSize(0.02);
0601   xaxis->SetTextFont(40);
0602   xaxis->SetLabelSize(0.02);
0603   xaxis->SetTickSize(0.02);
0604 
0605   yaxis->SetTitle(var2name);
0606   yaxis->SetTitleSize(0.02);
0607   yaxis->SetName("yaxis");
0608   yaxis->SetTextFont(40);
0609   yaxis->SetLabelSize(0.02);
0610   yaxis->SetTickSize(0.02);
0611 
0612   newpad1->cd();
0613   yaxis->Draw();
0614   newpad2->cd();
0615   xaxis->Draw();
0616 
0617   if(legendLabels.size()>0) leg->Draw();
0618 
0619   //  canv->Write();
0620   canv->Print(pngDir+"/"+canvN+".png");
0621   canv->Print(pngDir+"/"+canvN+".pdf");
0622   outfile->cd("/");
0623   canv->Write();
0624   // for(int i = 0; i <nvar1; i++){
0625   //   for(int j = 0; j < nvar2; j++){
0626   //     histArray[i][j]->Write();
0627   //   }
0628   // }
0629 };
0630 
0631 
0632 //=========================================================================
0633 /* ALGORITHM: draw ratios of histograms from different bins in their respective bins
0634 on axis of bin variables, e.g. Q2 vs x.
0635 */
0636 // -- process one histArr
0637 void PostProcessor::DrawRatioInBins(
0638     TString outName,    
0639     std::vector<std::vector<Histos*>>& histArr,
0640     TString histName,
0641     TString var1name, int nvar1, double var1low, double var1high, bool var1log,
0642     TString var2name, int nvar2, double var2low, double var2high, bool var2log,
0643     bool intgrid1, bool intgrid2, // grid option for small plots
0644     bool renormalize // if true, normalize histograms
0645 ){
0646   std::vector<std::vector<std::vector<Histos*>>> histArrList;
0647   histArrList.push_back(histArr); // build list of histArrs with the one histArr
0648   this->DrawRatioInBins( // call the main algo below, which uses histArrList
0649       outName,    
0650       histArrList,
0651       histName,
0652       var1name, nvar1, var1low, var1high, var1log,
0653       var2name, nvar2, var2low, var2high, var2log,
0654       intgrid1, intgrid2,
0655       renormalize
0656       );
0657 };
0658 // -- process list of histArrs
0659 void PostProcessor::DrawRatioInBins(
0660     TString outName,    
0661     std::vector<std::vector<std::vector<Histos*>>>& histArrList,
0662     TString histName,
0663     TString var1name, int nvar1, double var1low, double var1high, bool var1log,
0664     TString var2name, int nvar2, double var2low, double var2high, bool var2log,
0665     bool intgrid1, bool intgrid2, // grid option for small plots
0666     bool renormalize // if true, normalize histograms (useful for drawing 1D hists with "SAME")
0667 ){
0668   // default values set for nvar1==nvar2
0669   int canvx = 1400;
0670   int canvy = 1200;
0671   double botmargin = 0.2;
0672   double leftmargin = 0.2;
0673   double xaxisy = 0.04;
0674   double xaxisx1 = 0.08;
0675   double xaxisx2 = 0.97;
0676   double yaxisx = 0.04;
0677   double yaxisy1 = 0.085;
0678   double yaxisy2 = 0.97;
0679   if(nvar1 > nvar2){
0680     // different canvas sizing/axis position for unequal binning
0681     canvx = 2200;
0682     canvy = 1400;
0683     xaxisx1 = 0.075;
0684     xaxisx2 = 0.975;
0685     yaxisy1 = 0.08;
0686   };
0687 
0688   TString canvN = "canv_"+outName+"_"+histName;
0689   TCanvas *canv = new TCanvas(canvN,canvN, canvx, canvy);
0690   TPad *mainpad = new TPad("mainpad", "mainpad", 0.07, 0.07, 0.98, 0.98);
0691 
0692   TLegend *leg;
0693   if(legendLabels.size()>0) {
0694     double legX0 = 0.1;
0695     double legY0 = 0.95;
0696     double legWidth = 0.15;
0697     double legRowHeight = 0.05;
0698     leg = new TLegend(
0699         legX0,
0700         legY0,
0701         legX0 + legWidth,
0702         legY0 - legRowHeight*legendLabels.size()
0703         );
0704   };
0705 
0706   mainpad->SetFillStyle(4000);
0707   mainpad->Divide(nvar1,nvar2,0,0);
0708   mainpad->Draw();
0709   TLine * lDIRC = new TLine(6,-1,6,1);
0710   TLine * lDIRClow = new TLine(0.5,-1,0.5,1);
0711   TLine * lmRICH = new TLine(2,-1,2,-4);
0712   TLine * lDRICH = new TLine(2.5,1,2.5,4);
0713   lDIRC->SetLineColor(kRed);
0714   lDIRClow->SetLineColor(kRed);
0715   lmRICH->SetLineColor(kRed);
0716   lDRICH->SetLineColor(kRed);
0717   // TH1* histArray[nvar1][nvar2]; // TODO: re-enable
0718   int drawpid = 0;
0719   outfile->cd("/");
0720   canv->Write();
0721 
0722   // get histograms from Histos 2D vector
0723   Histos *H;
0724   for(int i = 0; i < nvar1; i++){
0725     for(int j = 0; j < nvar2; j++){
0726       int count = 0;
0727       int dims;
0728       TH1 *histref;
0729       for(auto histArr : histArrList) {
0730         H = histArr[i][j];
0731         TH1 *hist = H->Hist(histName);
0732         dims = hist->GetDimension();
0733         // histArray[i][j] = hist;
0734         hist->SetTitle("");
0735         //hist->GetXaxis()->SetTitle("");
0736         //hist->GetYaxis()->SetTitle("");
0737         //hist->GetXaxis()->SetLabelSize(0);
0738         //hist->GetYaxis()->SetLabelSize(0);
0739 
0740         hist->SetLineWidth(3);
0741         switch(count) {
0742           case 0:
0743         histref = (TH1*) hist->Clone();
0744             hist->SetLineColor(kBlack);
0745         //            hist->SetFillColor(kGray);
0746             break;
0747           case 1:
0748             hist->SetLineColor(kAzure);
0749             hist->SetLineStyle(1);
0750             break;
0751           case 2:
0752             hist->SetLineColor(kRed);
0753             hist->SetLineStyle(7);
0754             break;
0755           case 3:
0756             hist->SetLineColor(kGreen+2);
0757             hist->SetLineStyle(9);
0758             break;
0759         }
0760 
0761 
0762     hist->Divide(histref);
0763     
0764         if(legendLabels.size()>0 && i==0 && j==0) {
0765           try { leg->AddEntry(hist,legendLabels[count],"LF"); }
0766           catch(const std::out_of_range &e) {
0767             cerr << "ERROR: legendLabels not filled correctly" << endl;
0768           };
0769         };
0770 
0771         if(count==0) {
0772           mainpad->cd((nvar2-j-1)*nvar1 + i + 1);
0773           gPad->SetLogx(H->GetHistConfig(histName)->logx);
0774           gPad->SetLogy(H->GetHistConfig(histName)->logy);
0775           gPad->SetLogz(H->GetHistConfig(histName)->logz);
0776           gPad->SetGridy(intgrid2);
0777           gPad->SetGridx(intgrid1);
0778         }
0779 
0780         // renormalize
0781         if(renormalize) {
0782           // hist->Scale(1/hist->GetMaximum());
0783           // hist->Scale(1/hist->GetEntries());
0784       //         hist->Scale(1/hist->Integral());
0785       hist->GetYaxis()->SetRangeUser(0.,3.);
0786         }
0787 
0788         TString drawStr = "";
0789         switch(dims) {
0790           case 1:
0791             drawStr = "HIST MIN0";
0792             break;
0793           case 2:
0794             drawStr = "COLZ";
0795             break;
0796           case 3:
0797             drawStr = "BOX";
0798             break;
0799         };
0800         if(count>0) drawStr += " SAME";
0801 
0802         if( hist->GetEntries() > 0 ) {  
0803           hist->Draw(drawStr);
0804           if(drawpid){
0805             lDIRClow->Draw();
0806             lDIRC->Draw();
0807             lmRICH->Draw();
0808             lDRICH->Draw();
0809           }
0810           hist->GetXaxis()->SetLabelSize(0.04);
0811           hist->GetYaxis()->SetLabelSize(0.04);
0812           hist->GetXaxis()->SetTitleSize(0.05);
0813           hist->GetYaxis()->SetTitleSize(0.05);
0814           hist->GetXaxis()->SetTitleOffset(0.9);
0815           hist->GetXaxis()->SetLabelOffset(0.0005);
0816           if(dims==1) {
0817             hist->GetYaxis()->SetLabelSize(0.00); // suppress y-axis labels (since each subplot has its own scale)
0818           }
0819         }
0820         count++;
0821       }; // end for(histArrList)
0822 
0823       // some formatting for 1D plots drawn with "SAME"
0824       if(dims==1) {
0825         UnzoomVertical(gPad,"",true,H->GetHistConfig(histName)->logy); // allow all "SAME" plots to be visible; 3rd arg forces min to be 0
0826       }
0827     };
0828   };
0829   canv->cd();
0830 
0831   TPad *newpad1 = new TPad("newpad1","full pad",0,0,1,1);
0832   TPad *newpad2 = new TPad("newpad2","full pad",0,0,1,1);
0833   newpad1->SetFillStyle(4000);
0834   newpad1->Draw();
0835   newpad2->SetFillStyle(4000);
0836   newpad2->Draw();
0837 
0838   TString xopt, yopt;
0839   if(var1log) xopt = "GS";
0840   else xopt = "S";
0841   if(var2log) yopt = "GS";
0842   else yopt = "S";
0843 
0844   TGaxis *xaxis = new TGaxis(xaxisx1,xaxisy,xaxisx2,xaxisy,var1low,var1high,510,xopt);
0845   TGaxis *yaxis = new TGaxis(yaxisx,yaxisy1,yaxisx,yaxisy2,var2low,var2high,510,yopt);
0846   xaxis->SetTitle(var1name);
0847   xaxis->SetName("xaxis");
0848   xaxis->SetTitleSize(0.02);
0849   xaxis->SetTextFont(40);
0850   xaxis->SetLabelSize(0.02);
0851   xaxis->SetTickSize(0.02);
0852 
0853   yaxis->SetTitle(var2name);
0854   yaxis->SetTitleSize(0.02);
0855   yaxis->SetName("yaxis");
0856   yaxis->SetTextFont(40);
0857   yaxis->SetLabelSize(0.02);
0858   yaxis->SetTickSize(0.02);
0859 
0860   newpad1->cd();
0861   yaxis->Draw();
0862   newpad2->cd();
0863   xaxis->Draw();
0864 
0865   if(legendLabels.size()>0) leg->Draw();
0866 
0867   //  canv->Write();
0868   canv->Print(pngDir+"/"+canvN+".png");
0869   canv->Print(pngDir+"/"+canvN+".pdf");
0870   outfile->cd("/");
0871   canv->Write();
0872   // for(int i = 0; i <nvar1; i++){
0873   //   for(int j = 0; j < nvar2; j++){
0874   //     histArray[i][j]->Write();
0875   //   }
0876   // }
0877 };
0878 
0879 
0880 //=========================================================================
0881 
0882 /* ALGORITHM: draw a ratio of all 1D histograms in the specified histogram set
0883 * - the ratio will be of `numerSet` over `denomSet`
0884 * - canvases will be created, and depending on the setting of `plotRatioOnly`,
0885 *   either just the ratio will be plotted, or the two histograms will also be
0886 *   plotted
0887 * - Use `outName` to specify a names for output canvases
0888 * - summary canvases are also created, which can combine multiple ratio plots;
0889 *   they are accumulated in `summaryCanvMap` and colors are set by `nsum`; call
0890 *   `ResetVars` to reset the accumulation
0891  * - when done looping, call `FinishDrawRatios`
0892 */
0893 void PostProcessor::DrawRatios(
0894     TString outName, Histos *numerSet, Histos *denomSet, Bool_t plotRatioOnly
0895 ) {
0896 
0897   cout << "draw ratios " << outName << "..." << endl;
0898   enum HHenum {num,den};
0899   Histos *HH[2];
0900   HH[num] = numerSet;
0901   HH[den] = denomSet;
0902   TH1 *hist[2];
0903   TH1D *ratio;
0904   Double_t ny,ry,err;
0905   TCanvas *canv;
0906   TString histT[2], tok[2], uniqT[2], sameT[2];
0907   Ssiz_t tf;
0908 
0909   outfile->cd("/");
0910   outfile->mkdir(outName);
0911   outfile->cd(outName);
0912 
0913   // loop over 1D histograms
0914   for(TString varName : HH[num]->VarNameList) {
0915     hist[num] = HH[num]->Hist(varName);
0916     if(hist[num] != nullptr && hist[num]->GetDimension()==1) {
0917       hist[den] = HH[den]->Hist(varName);
0918 
0919       // filter title
0920       // TODO: there is a bug here now, ratio denominators aren't quite right
0921       for(int f=0; f<2; f++)  histT[f] = hist[f]->GetTitle();
0922       for(int f=0; f<2; f++) {
0923         tf = 0;
0924         sameT[f] = "";
0925         uniqT[f] = "";
0926         while(histT[f].Tokenize(tok[f],tf,", ")) {
0927           if(histT[(f+1)%2].Contains(tok[f])) sameT[f] += ", " + tok[f];
0928           else                                uniqT[f] += ", " + tok[f];
0929         };
0930         sameT[f](TRegexp("^, "))="";
0931         uniqT[f](TRegexp("^, "))="";
0932         /*cout <<  "histT = " << histT[f] << endl;
0933           cout <<  "uniqT = " << uniqT[f] << endl;
0934           cout <<  "sameT = " << sameT[f] << endl;*/
0935         hist[f]->SetTitle(sameT[f]);
0936       };
0937       TString ratioStr = "ratio("+uniqT[num]+"/"+uniqT[den]+")";
0938 
0939       // calculate ratio
0940       ratio = (TH1D*) hist[num]->Clone();
0941       ratio->Divide(hist[den]);
0942       TString ratioT = ratio->GetTitle();
0943       TString ratioSummaryT = ratioT;
0944       ratioT(TRegexp("distribution")) = ratioStr;
0945       ratioSummaryT(TRegexp("distribution")) = "ratios";
0946       ratio->SetTitle(ratioT);
0947 
0948       // calculate ratio errors, assuming ratio represents an efficiency
0949       for(int b=0; b<=ratio->GetNbinsX(); b++) {
0950         ny = hist[den]->GetBinContent(b);
0951         ry = ratio->GetBinContent(b);
0952         err = TMath::Abs(ny)>0 ?
0953           TMath::Sqrt( ry*(1-ry) / ny )
0954           : 0;
0955         ratio->SetBinError(b,err);
0956       };
0957 
0958       // canvas
0959       TString canvN = "canv_"+outName+"_"+varName;
0960       Int_t dimx=800;
0961       Int_t dimy=700;
0962       canv = new TCanvas(canvN,canvN, dimx*(plotRatioOnly?1:2), dimy);
0963 
0964       // draw plots
0965       int numPads;
0966       if(!plotRatioOnly) { canv->Divide(2,1); numPads=2; }
0967       else { canv->Divide(1,1); numPads=1; };
0968       if(!plotRatioOnly) {
0969         canv->cd(2);
0970         hist[num]->SetLineColor(kYellow-8);
0971         hist[den]->SetLineColor(kAzure-7);
0972         hist[num]->SetMarkerColor(kYellow-8);
0973         hist[den]->SetMarkerColor(kAzure-7);
0974         hist[num]->SetMarkerStyle(kFullCircle);
0975         hist[den]->SetMarkerStyle(kFullCircle);
0976         for(int f=0; f<2; f++) {
0977           hist[f]->SetMarkerSize(1.0);
0978           hist[f]->SetLineWidth(2);
0979           hist[f]->GetXaxis()->SetLabelSize(0.06);
0980           hist[f]->GetYaxis()->SetLabelSize(0.06);
0981           hist[f]->GetXaxis()->SetTitleSize(0.06);
0982           hist[f]->GetYaxis()->SetTitleSize(0.06);
0983           hist[f]->GetXaxis()->SetTitleOffset(1.2);
0984           //if(plotName=="PperpDistLin") hist[f]->GetXaxis()->SetRangeUser(0,1.7);
0985         };
0986         hist[den]->Draw("EX0 P");
0987         hist[num]->Draw("EX0 P SAME");
0988       };
0989       canv->cd(1);
0990       ratio->GetYaxis()->SetTitle(ratioStr);
0991       ratio->SetLineColor(kBlack);
0992       ratio->SetMarkerColor(kBlack);
0993       ratio->SetMarkerStyle(kFullCircle);
0994       ratio->SetMarkerSize(1.0);
0995       ratio->SetLineWidth(2);
0996       ratio->GetXaxis()->SetLabelSize(0.06);
0997       ratio->GetYaxis()->SetLabelSize(0.06);
0998       ratio->GetXaxis()->SetTitleSize(0.06);
0999       ratio->GetYaxis()->SetTitleSize(0.06);
1000       ratio->GetXaxis()->SetTitleOffset(1.2);
1001       ratio->GetYaxis()->SetRangeUser(0,1.3);
1002       ratio->Draw("EX0 P");
1003       for(int pad=1; pad<=numPads; pad++) {
1004         canv->GetPad(pad)->SetGrid(1,1);
1005         canv->GetPad(pad)->SetLogx(HH[num]->GetHistConfig(varName)->logx);
1006         if(pad>1) {
1007           canv->GetPad(pad)->SetLogy(HH[num]->GetHistConfig(varName)->logy);
1008           canv->GetPad(pad)->SetLogz(HH[num]->GetHistConfig(varName)->logz);
1009         };
1010         canv->GetPad(pad)->SetBottomMargin(0.15);
1011         canv->GetPad(pad)->SetLeftMargin(0.15);
1012       };
1013       canv->Print(pngDir+"/"+canvN+".png");
1014       canv->Write();
1015 
1016       // add to summary canvas
1017       TH1D *ratioSummary = (TH1D*) ratio->Clone();
1018       ratioSummary->SetTitle(ratioSummaryT);
1019       ratioSummary->SetYTitle("ratio");
1020       ratioSummary->SetLineColor  (nsum<nsumMax?summaryColor[nsum]:kBlack);
1021       ratioSummary->SetMarkerColor(nsum<nsumMax?summaryColor[nsum]:kBlack);
1022       ratioSummary->SetMarkerStyle(nsum<nsumMax?summaryStyle[nsum]:kFullCircle);
1023       auto kv = summaryCanvMap.find(varName);
1024       if(kv==summaryCanvMap.end()) {
1025         summaryCanv = new TCanvas(
1026             "summaryCanv_ratio"+varName,
1027             "summaryCanv_ratio"+varName,
1028             dimx,dimy
1029             );
1030         summaryCanv->SetGrid(1,1);
1031         summaryCanv->SetLogx(HH[num]->GetHistConfig(varName)->logx);
1032         summaryCanv->SetBottomMargin(0.15);
1033         summaryCanv->SetLeftMargin(0.15);
1034         summaryCanvMap.insert({varName,summaryCanv});
1035         ratioSummary->Draw();
1036       } else {
1037         kv->second->cd();
1038         ratioSummary->Draw("SAME");
1039       };
1040 
1041     };
1042   };
1043   outfile->cd("/");
1044   nsum++;
1045 };
1046 void PostProcessor::FinishDrawRatios(TString summaryDir) {
1047   std::cout << "CALL PostProcessor::FinishDrawRatios" << endl;
1048   // write summary canvas
1049   outfile->cd("/");
1050   outfile->mkdir(summaryDir);
1051   outfile->cd(summaryDir);
1052   for(auto const& kv : summaryCanvMap) {
1053     kv.second->Write();
1054     kv.second->Print(pngDir+"/"+summaryDir+"_"+kv.first+".png");
1055   };
1056   outfile->cd("/");
1057   this->ResetVars();
1058 };
1059 
1060 
1061 
1062 //=========================================================================
1063 // TEXT FILE: start new text file, append to text file, print text file
1064 void PostProcessor::StartTextFile(TString datFile, TString firstLine) {
1065   gSystem->RedirectOutput(datFile,"w");
1066   if(firstLine!="") cout << firstLine << endl;
1067   gSystem->RedirectOutput(0);
1068 };
1069 void PostProcessor::AppendToTextFile(TString datFile, TString appendText="") {
1070   gSystem->RedirectOutput(datFile,"a");
1071   cout << appendText << endl;
1072   gSystem->RedirectOutput(0);
1073 };
1074 void PostProcessor::PrintTextFile(TString datFile) {
1075   gROOT->ProcessLine(".! cat "+datFile);
1076 };
1077 
1078 
1079 //=========================================================================
1080 // TEXT FILE: pipe `inputFile` through `column -t` and append output to `outputFile`
1081 void PostProcessor::Columnify(TString inputFile, TString outputFile) {
1082   gSystem->RedirectOutput(outputFile,"a");
1083   gROOT->ProcessLine(".! cat "+inputFile+" | column -t");
1084   gSystem->RedirectOutput(0);
1085 };
1086 
1087 
1088 //=========================================================================
1089 // get pointer to BinSet associated with variable
1090 BinSet *PostProcessor::GetBinSet(TString varName) {
1091   return (BinSet*) infile->Get(varName+"_bins");
1092 };
1093 // get pointer to CutDef associated with a particular bin
1094 CutDef *PostProcessor::GetBinCut(TString varName, Int_t binNum) {
1095   return this->GetBinSet(varName)->Cut(binNum);
1096 };
1097 // get list of bin numbers
1098 std::vector<int> PostProcessor::GetBinNums(TString varName) {
1099   BinSet *B = this->GetBinSet(varName);
1100   std::vector<int> retVec;
1101   if(B==nullptr) {
1102     cerr << "ERROR: unknown variable " << varName << " in PostProcessor::GetBinNums" << endl;
1103   } else {
1104     for(int n=0; n<B->GetNumBins(); n++) retVec.push_back(n);
1105   };
1106   return retVec;
1107 };
1108 
1109 // return true if the bin is "full" range, and it's not the only bin
1110 Bool_t PostProcessor::SkipFull(TString varName, Int_t binNum) {
1111   return GetBinSet(varName)->GetNumBins() > 1 &&
1112          GetBinCut(varName,binNum)->GetCutType()=="Full";
1113 };
1114 
1115 
1116 //=========================================================================
1117 // destructor
1118 PostProcessor::~PostProcessor() {
1119 };
1120
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!