Corral/src/QPTPage.h

0001 #ifndef QPTPAGE_H
0002 #define QPTPAGE_H
0003 //---- QPTPage.h -- the ir2=2 (Q_inv, pT_pair) summary page, drawn the same way in the corral and the Finalize PDFs
0004 //----   (README_meanpT_from_CQ.md sec 11). One page per pair type, 2x2:
0005 //----   (1) C(Q) over all pT_pair, with the CQResonances.h signposts and the variable-binning note
0006 //----   (2) C(Q, pT_pair) map
0007 //----   (3) C(Q) in four pT_pair slices
0008 //----   (4) <pT_pair> per Q bin, sibling and mixed (the ingredient of the <pT> moment method)
0009 //---- S, M: sibling and mixed in any common normalisation (corral: zvtx-averaged rho2(S), rho2(M); Finalize: zvtx
0010 //----   sums, pairs and mixed pairs in expected-sibling units); only ratios and means are drawn, so this does not matter.
0011 //---- C1: optional C(Q) over all pT_pair with better errors (Finalize: subgroup errors); 0 = built from S and M.
0012 //---- Needs CQResonances.h (QPT_NOTE) included before it.
0013 #include <TCanvas.h>
0014 #include <TH2.h>
0015 #include <TH1D.h>
0016 #include <TLatex.h>
0017 #include <TLegend.h>
0018 #include <TLegendEntry.h>
0019 #include <vector>
0020 #include <algorithm>
0021
0022 void DrawCQSignposts(int ipid1, int ipid2, double ymin, double ymax, double xmax);  // corral_loop.cxx
0023
0024 //---- C = s/m bin by bin; error from both (relative errors in quadrature); bins with s or m empty stay empty
0025 inline TH1D* QPTRatio(TH1D *s, TH1D *m, const char* name){
0026     TH1D *c = (TH1D*)s->Clone(name); c->SetDirectory(0); c->Reset();
0027     for (int b=1;b<=s->GetNbinsX();b++){
0028         double vs = s->GetBinContent(b), vm = m->GetBinContent(b);
0029         if (vs<=0. || vm<=0.) continue;
0030         double C = vs/vm, es = s->GetBinError(b), em = m->GetBinError(b);
0031         c->SetBinContent(b,C); c->SetBinError(b,C*sqrt(pow(es/vs,2)+pow(em/vm,2)));
0032     }
0033     return c;
0034 }
0035 //---- y range of a C(Q): 1st-99th percentile of the bins above Q = qlo that have an error (the lowest-Q and
0036 //---- the sparse high-pT_pair bins are noisy and would set the range otherwise)
0037 inline void QPTYRange(std::vector<TH1D*> hs, double qlo, double top){
0038     std::vector<double> v;
0039     for (TH1D *h : hs) for (int b=h->FindBin(qlo+1e-6);b<=h->GetNbinsX();b++) if (h->GetBinError(b)>0.) v.push_back(h->GetBinContent(b));
0040     if (v.size()<5) return;
0041     std::sort(v.begin(),v.end());
0042     double ymn = v[(size_t)(0.01*v.size())], ymx = v[(size_t)(0.99*v.size())];
0043     if (ymx>ymn){ hs[0]->SetMinimum(ymn-0.10*(ymx-ymn)); hs[0]->SetMaximum(ymx+top*(ymx-ymn)); }
0044 }
0045
0046 inline void DrawQPTPage(TCanvas *c, TH2 *S, TH2 *M, TH1 *C1, int pid1, int pid2, const TString &pair, const TString &src, const TString &tag){
0047     c->Clear(); c->Divide(2,2,0.0001,0.0001);
0048     const double QLO    = 0.02;     // display ranges from Q > 20 MeV
0049     //---- (1) C(Q) over all pT_pair (under- and overflow included)
0050     c->cd(1);
0051     {
0052         TH1D *h;
0053         if (C1){ h = (TH1D*)C1->Clone(Form("QPT1_%s",tag.Data())); h->SetDirectory(0); }
0054         else {
0055             TH1D *s = S->ProjectionX(Form("QPT1s_%s",tag.Data()),0,-1,"e"), *m = M->ProjectionX(Form("QPT1m_%s",tag.Data()),0,-1,"e");
0056             h = QPTRatio(s,m,Form("QPT1_%s",tag.Data())); delete s; delete m;
0057         }
0058         h->SetTitle(Form("%s, C(Q), all p_{T,pair} (%s);Q_{inv} (GeV);C(Q)",pair.Data(),src.Data()));
0059         h->SetMarkerStyle(20); h->SetMarkerSize(0.3);
0060         QPTYRange({h},QLO,0.30); h->Draw("HIST"); h->Draw("E SAME"); gPad->Update();
0061         DrawCQSignposts(pid1,pid2,gPad->GetUymin(),gPad->GetUymax(),h->GetXaxis()->GetXmax());
0062         TLatex t; t.SetNDC(); t.SetTextFont(42); t.SetTextSize(0.032); t.SetTextColor(kRed+1);
0063         t.DrawLatex(0.15,0.87,QPT_NOTE);
0064     }
0065     //---- (2) C(Q, pT_pair) map, z range from the 2nd-98th percentile of the filled cells
0066     c->cd(2);
0067     {
0068         gPad->SetRightMargin(0.13);
0069         TH2 *h  = (TH2*)S->Clone(Form("QPT2_%s",tag.Data())); h->SetDirectory(0); h->Reset();
0070         std::vector<double> v;
0071         for (int i=1;i<=S->GetNbinsX();i++) for (int j=1;j<=S->GetNbinsY();j++){
0072             double vs = S->GetBinContent(i,j), vm = M->GetBinContent(i,j);
0073             if (vs>0. && vm>0.){ h->SetBinContent(i,j,vs/vm); v.push_back(vs/vm); }
0074         }
0075         h->SetTitle(Form("%s, C(Q,p_{T,pair}) (%s, white: no sibling or no mixed pairs);Q_{inv} (GeV);p_{T,pair} (GeV/c)",pair.Data(),src.Data()));
0076         if (v.size()>10){ std::sort(v.begin(),v.end()); h->SetMinimum(v[(size_t)(0.02*v.size())]); h->SetMaximum(v[(size_t)(0.98*v.size())]); }
0077         h->Draw("colz");
0078     }
0079     //---- (3) C(Q) in four pT_pair slices
0080     c->cd(3);
0081     {
0082         const double edge[5]    = {0.,0.5,1.,2.,5.};
0083         const int kcol[4]       = {kBlue,kGreen+2,kOrange+4,kRed};
0084         std::vector<TH1D*> hk;
0085         for (int k=0;k<4;k++){
0086             int j1 = S->GetYaxis()->FindBin(edge[k]+1e-6), j2 = S->GetYaxis()->FindBin(edge[k+1]-1e-6);
0087             TH1D *s = S->ProjectionX(Form("QPT3s_%s_%d",tag.Data(),k),j1,j2,"e"), *m = M->ProjectionX(Form("QPT3m_%s_%d",tag.Data(),k),j1,j2,"e");
0088             TH1D *h = QPTRatio(s,m,Form("QPT3_%s_%d",tag.Data(),k)); delete s; delete m;
0089             h->SetLineColor(kcol[k]); h->SetMarkerColor(kcol[k]); h->SetMarkerStyle(20); h->SetMarkerSize(0.3);
0090             hk.push_back(h);
0091         }
0092         hk[0]->SetTitle(Form("%s, C(Q) in p_{T,pair} slices (%s);Q_{inv} (GeV);C(Q)",pair.Data(),src.Data()));
0093         QPTYRange(hk,QLO,0.30);
0094         hk[0]->Draw("E1"); for (int k=1;k<4;k++) hk[k]->Draw("E1 same");
0095         TLegend *lg = new TLegend(0.55,0.68,0.89,0.89); lg->SetBorderSize(0); lg->SetFillStyle(0); lg->SetTextSize(0.035);
0096         for (int k=0;k<4;k++){ TLegendEntry *le = lg->AddEntry(hk[k],Form("%.1f < p_{T,pair} < %.1f GeV/c",edge[k],edge[k+1]),"l"); le->SetTextColor(kcol[k]); }
0097         lg->Draw();
0098     }
0099     //---- (4) <pT_pair> per Q bin, sibling and mixed; error = rms / sqrt(effective entries) (= counts where S or M carry no errors)
0100     c->cd(4);
0101     {
0102         TH1D *mp[2];
0103         for (int sm=0;sm<2;sm++){
0104             TH2 *H  = sm ? M : S;
0105             mp[sm]  = H->ProjectionX(Form("QPT4_%s_%d",tag.Data(),sm),0,-1); mp[sm]->SetDirectory(0); mp[sm]->Reset();
0106             for (int i=1;i<=H->GetNbinsX();i++){
0107                 double sw=0, swy=0, swy2=0, se2=0;
0108                 for (int j=1;j<=H->GetNbinsY();j++){
0109                     double w = H->GetBinContent(i,j), y = H->GetYaxis()->GetBinCenter(j);
0110                     sw += w; swy += w*y; swy2 += w*y*y; se2 += pow(H->GetBinError(i,j),2);
0111                 }
0112                 if (sw<=0.) continue;
0113                 double mean = swy/sw, rms = sqrt(std::max(0.,swy2/sw-mean*mean)), neff = (se2>0.) ? sw*sw/se2 : sw; // no errors (Finalize mixed sums): counts
0114                 mp[sm]->SetBinContent(i,mean); mp[sm]->SetBinError(i,rms/sqrt(neff));
0115             }
0116             mp[sm]->SetLineColor(sm ? kRed : kBlack); mp[sm]->SetMarkerColor(sm ? kRed : kBlack); mp[sm]->SetMarkerStyle(20); mp[sm]->SetMarkerSize(0.3);
0117         }
0118         mp[0]->SetTitle(Form("%s, #LTp_{T,pair}#GT per Q bin (%s, p_{T,pair} bin centres);Q_{inv} (GeV);#LTp_{T,pair}#GT (GeV/c)",pair.Data(),src.Data()));
0119         QPTYRange({mp[0],mp[1]},QLO,0.25);
0120         mp[0]->Draw("E1"); mp[1]->Draw("E1 same");
0121         TLegend *lg = new TLegend(0.55,0.75,0.89,0.89); lg->SetBorderSize(0); lg->SetFillStyle(0); lg->SetTextSize(0.035);
0122         lg->AddEntry(mp[0],"sibling","l"); TLegendEntry *le = lg->AddEntry(mp[1],"mixed","l"); le->SetTextColor(kRed);
0123         lg->Draw();
0124     }
0125     c->cd(); c->Update();
0126 }
0127 #endif