Corral/src/corral_finalize.cxx

0001 #include "corral_class.h"
0002 #include <TSystem.h>
0003 #include "fluct_common.h"
0004 #include "CQResonances.h"   // QPT_NOTE for the ir2=2 page
0005 #include "QPTPage.h"        // ir2=2 (Q,pT_pair) summary page, shared with corral_m
0006 #include "PairTypes.h"
0007 #include "EtaLabels.h"      // h+, h- pages in eta
0008 #include "finalize_hists.h"
0009 #include "CrossingCorrect.h"    // same dirty-side rule and correction as CalcRm::Calculate
0010 #include "ZvtxAverage.h"        // same Zvtx average as CalcRm::Calculate
0011 #include "SegLive.h"            // TPC live fraction per segment (the segveto rule)
0012 #include "FinalizeQA.h"         // title page, QA pages, vs-segment pages, info/ of the chunks
0013 #include "corral_version.h"     // CORRAL_GIT, CORRAL_BUILD_TIME of this binary
0014 #include <fstream>
0015 #include <map>
0016 #include <algorithm>
0017 #include <TCanvas.h>
0018 #include <TLegend.h>
0019 #include <TLegendEntry.h>
0020 #include <TStyle.h>
0021 #include <TGraphErrors.h>
0022 #include <TF1.h>
0023 #include <TLine.h>
0024 #include <TBox.h>
0025 #include <TLatex.h>
0026 #include <TFitResult.h>
0027 #include <TParameter.h>
0028 #include <TTree.h>
0029 #include <TObjString.h>
0030 #include <TArrayD.h>
0031 #include <sstream>
0032
0033 //------------------------------------------------------------
0034 // Finalize (README_Finalize.md). Combine the chunk outputs of a chunked production.
0035 //   CF histograms (hrho2_<ir2>, hrho1rho1_<ir2>, hC2_<ir2>, hR2_<ir2>, hrho2C_1, hC2C_1, hR2C_1):
0036 //     step 2, per zvtx bin z: S = sum_c(nevt S)/sum_c(nevt), M = sum_c(nevt df M)/sum_c(nevt df) over the
0037 //     chunks where z was computed (df>0, nevt>=2); then C2/R2 exactly as CalcRm::Calculate, crossing
0038 //     correction (ir2=1) with CrossingCorrect.h, and the Zvtx average with ZvtxAverageValid (weights =
0039 //     combined hzvtx fractions).
0040 //     Zvtx-average validity (option A, user 2026-09-26): zvtx bin z counts for a bin iff rho2(M)>0 and
0041 //     N_exp = nevt_z*rho2(M)_z >= N_min (PairTypes.h 3rd field, per pairtype) at full statistics (edge zvtx bins with O(1) pairs
0042 //     otherwise dominate the noise of the average).
0043 //     step 3, subgroup errors: the SAME chain run on each chunk alone (a subgroup = one chunk), with the
0044 //     FULL-STATS zvtx weights and validity masks, so each subgroup is the same linear combination as the
0045 //     central value. Per bin, error = s/sqrt(N_valid), s = sample std (N-1) over the chunks usable there
0046 //     (every zvtx bin -- and, corrected, every source bin -- the full-stats value uses is defined in the
0047 //     chunk); others are skipped; N_valid < NVALIDMIN -> error 0 (no reliable error, empty-bin rule).
0048 //     The central value stays the full-stats combination.
0049 //   Anything else in the list: plain sum over chunks (step 1).
0050 // Every listed histogram is compared with the single-job reference (root/<set>/<CorralFilePrefix(set)>.root) and written to
0051 // root/<base>.root; one pdf page per name (4x4 pairtypes), plus subgroup-error pages for the CF names.
0052 //
0053 void DrawCQSignposts(int ipid1, int ipid2, double ymin, double ymax, double xmax);  // corral_loop.cxx
0054
0055 namespace {
0056     //---- is this base name a CF histogram Finalize rebuilds?  kind: 0 S, 1 M, 2 C2, 3 R2, 4 rho2C, 5 C2C, 6 R2C
0057     bool ParseCF(const TString& base, int& kind, int& ir2){
0058         const char* pre[4]  = {"hrho2_","hrho1rho1_","hC2_","hR2_"};
0059         if (base=="hrho2C_1"){ kind=4; ir2=1; return true; }
0060         if (base=="hC2C_1")  { kind=5; ir2=1; return true; }
0061         if (base=="hR2C_1")  { kind=6; ir2=1; return true; }
0062         for (int k=0;k<4;k++){
0063             TString p   = pre[k];
0064             if (base.BeginsWith(p) && base.Length()==p.Length()+1 && isdigit(base[p.Length()])){
0065                 kind    = k; ir2 = base[p.Length()]-'0';
0066                 return ir2>=0 && ir2<=2;
0067             }
0068         }
0069         return false;
0070     }
0071     //---- base name of each CF kind (same order as ParseCF), %d = ir2
0072     const char* KINDNAME[7] = {"hrho2_%d","hrho1rho1_%d","hC2_%d","hR2_%d","hrho2C_%d","hC2C_%d","hR2C_%d"};
0073     //---- output of one run of the CF chain: Zvtx-averaged maps per kind, and their validity
0074     //---- (V: some used zvtx bin has rho2(M)>0; VC: same with the crossing-corrected mask; ir2=1 only)
0075     //---- keepMasks: also the per-zvtx-bin validity masks (VZ, uncorrected; MK, corrected) and the Zvtx
0076     //---- weights, so that the subgroups can be run with the full-stats ones (option A)
0077     struct CFSet {
0078         TH2D    *h[7];
0079         TH2D    *V, *VC;
0080         std::vector<TH2D*>  vz, mk;
0081         double  wz[64];
0082         bool    okz[64];
0083         CFSet(){ for (int k=0;k<7;k++) h[k]=0; V=0; VC=0; for (int z=0;z<64;z++){ wz[z]=0.; okz[z]=false; } }
0084         void Delete(){
0085             for (int k=0;k<7;k++){ delete h[k]; h[k]=0; } delete V; V=0; delete VC; VC=0;
0086             for (TH2D* h2 : vz) delete h2; vz.clear();
0087             for (TH2D* h2 : mk) delete h2; mk.clear();
0088         }
0089     };
0090     double Median(std::vector<double> v){
0091         if (v.empty()) return 0.;
0092         std::sort(v.begin(),v.end());
0093         return v[v.size()/2];
0094     }
0095 }
0096
0097 void corral::Finalize(){
0098     //
0099     TString PROJ    = CorralDir();
0100     TString SET     = (Dataset!="") ? Dataset : TString(FINALIZE_SET);  // -d <dataset>
0101     const int NVALIDMIN = 10;       // step 3 pick (a): fewer valid subgroups -> error 0
0102     //---- Zvtx-average validity threshold N_min per pairtype = PairTypes_Info[ipaty][2] (PairTypes.h);
0103     //---- the environment variable FINALIZE_NEXPMIN overrides it for ALL pairtypes (scans)
0104     const bool nexpOverride = (getenv("FINALIZE_NEXPMIN")!=0);
0105     const double nexpGlobal = nexpOverride ? atof(getenv("FINALIZE_NEXPMIN")) : -1.;
0106     //---- README_SplitTracks573 sec 38: acceptance-edge rule. A bin-zvtx whose rho2(M) is below MFRAC x the largest rho2(M) of
0107     //---- that bin over the zvtx slices (full statistics) is a partial-acceptance edge (the central-membrane mask, the y reach
0108     //---- of the outer slices); it is left out of the Zvtx average like N_exp < N_min. FINALIZE_MFRAC overrides (0 = off).
0109     const double mFrac  = getenv("FINALIZE_MFRAC") ? atof(getenv("FINALIZE_MFRAC")) : 0.3;
0110     cout<<"corral::Finalize -- Zvtx-average validity: rho2(M) >= "<<mFrac<<" x the bin's largest rho2(M) over the slices (acceptance edges; FINALIZE_MFRAC)"<<endl;
0111     if (nexpOverride) cout<<"corral::Finalize -- Zvtx-average validity: N_exp = nevt*rho2(M) >= "<<nexpGlobal<<" for ALL pairtypes (FINALIZE_NEXPMIN override)"<<endl;
0112     else              cout<<"corral::Finalize -- Zvtx-average validity: N_exp = nevt*rho2(M) >= N_min per pairtype (PairTypes.h 3rd field)"<<endl;
0113     //
0114     //---- number of chunks from the list manifest
0115     int nlists  = -1;
0116     {
0117         std::ifstream fm((PROJ+"lists/"+SET+"/manifest.txt").Data());
0118         std::string key; std::string val;
0119         while (fm>>key>>val){ if (key=="nlists") nlists = atoi(val.c_str()); }
0120     }
0121     if (nlists<=0){ cout<<"corral::Finalize -- no nlists in "<<PROJ<<"lists/"<<SET<<"/manifest.txt (set CORRAL_DIR?), exit"<<endl; exit(1); }
0122     cout<<"corral::Finalize -- set "<<SET<<": "<<nlists<<" chunks, "<<NFINALIZEHISTS<<" histogram base names x "<<NPairTypes<<" pairtypes"<<endl;
0123     //
0124     //---- open all chunk files (all must exist) and the reference
0125     std::vector<TFile*> fch;
0126     for (int ic=0;ic<nlists;ic++){
0127         TString fn  = PROJ+"root/"+SET+Form("/chunks/%s_%02d.root",CorralFilePrefix(SET.Data()),ic);    // per-dataset convention
0128         TFile *f    = TFile::Open(fn.Data());
0129         if (!f || f->IsZombie()){ cout<<"corral::Finalize -- cannot open "<<fn<<", exit"<<endl; exit(1); }
0130         fch.push_back(f);
0131     }
0132     TString REF = PROJ+"root/"+SET+"/"+CorralFilePrefix(SET.Data())+".root";    // per-dataset convention
0133     //---- a whole-dataset job is optional (most chunked passes have none, user 2026-10-02): without it, no reference pages
0134     TFile *fref = gSystem->AccessPathName(REF.Data()) ? nullptr : TFile::Open(REF.Data());
0135     if (fref && fref->IsZombie()) fref = nullptr;
0136     cout<<"corral::Finalize -- whole-dataset reference "<<REF<<(fref?": found, pages vs the reference follow":": none (optional), no pages vs the reference")<<endl;
0137     //
0138     //---- crossing correction: follow the chunks. A "nocross" job leaves its C maps empty (CalcRm), and Finalize then
0139     //---- leaves the C maps uncorrected too (copies of the plain ones), so a nocross set stays uncorrected end to end.
0140     bool chunkCross = true;
0141     {
0142         TH1 *hp = (TH1*)fch[0]->Get("hrho2_1_0");
0143         TH1 *hc = (TH1*)fch[0]->Get("hrho2C_1_0");
0144         if (hp && hp->Integral()>0 && (!hc || hc->Integral()==0)) chunkCross = false;
0145         cout<<"corral::Finalize -- crossing correction "<<(chunkCross?"ON":"OFF (the chunks ran nocross: C maps = plain maps)")<<endl;
0146     }
0147     bool anyCF  = false;
0148     for (int ih=0;ih<NFINALIZEHISTS;ih++){ int kd,ir; if (ParseCF(TString(FinalizeHists[ih]),kd,ir)) anyCF = true; }
0149     //
0150     //---- results: name -> histogram (all pairtypes), filled pairtype by pairtype; extra = written, not compared
0151     std::map<TString,TH1*> result;
0152     std::vector<TH1*> extra;
0153     std::vector<TObject*> extraW;       // written only (step 1b maps and graphs)
0154     //
0155     //---- README_Finalize step 1b: crossing loss = dphi-mirror asymmetry of C2 on the near side, per pairtype.
0156     //---- Tested half: the pions' dirty side (dirty side sign; dphi<0 here) for ALL pairtypes, so the V0
0157     //---- pairtypes (no correction) are the null test on the same half. Windows: 0 = whole near-side half,
0158     //---- 1 = near side |dy|<XLDYCUT, 2 = away-side half (internal null). Bins: valid at full stats in the
0159     //---- bin AND its mirror. D = sum[C2(mirror)-C2(bin)], f = D / sum rho2(S)(mirror); fM = the same with
0160     //---- rho2(M) (acceptance asymmetry check). D is built from avg(S)-avg(M) (= avg(C2), linear), so every
0161     //---- chunk contributes (no usability cut).
0162     const double XLDYCUT    = 0.2;
0163     struct XLoss {
0164         bool    have;
0165         int     ds, half;
0166         double  f[3], ef[3], fM[3], efM[3], fA[3], efA[3], fC, dycut, xlo, xhi;
0167         int     nuse[3];
0168         std::vector<double> fc[3], efc[3], fMc[3], fAc[3];
0169         std::vector<bool>   usec[3];
0170         std::vector<double> x;      // <n1><n2> per chunk (hmult)
0171         TH2D    *A;                 // full stats C2(mirror)-C2
0172         TH1D    *Ady;               // D per dy column, near half, subgroup errors
0173         double  chi2c, ndfc, slopex, eslopex, slopec, eslopec;  // fits (filled when drawing)
0174         XLoss(){ have=false; ds=0; half=-1; fC=0; dycut=0; xlo=xhi=0; A=0; Ady=0; chi2c=ndfc=slopex=eslopex=slopec=eslopec=0;
0175             for (int w=0;w<3;w++){ f[w]=ef[w]=fM[w]=efM[w]=fA[w]=efA[w]=0; nuse[w]=0; } }
0176     };
0177     std::vector<XLoss> xl(NPairTypes);
0178     //
0179     //---- physics pages (user 2026-09-26): per pairtype, the corral_m pages rebuilt from the Finalize maps,
0180     //---- all projections with SUBGROUP errors (the projection is formed inside every chunk with the
0181     //---- full-stats validity; in a bin where a chunk is not usable it takes the full-stats value -- the
0182     //---- fraction of such fill-ins is logged). dC2 = C2C - C2 = what the crossing correction adds.
0183     const double DYCUTPAGE      = 0.1;      // as corral_m's crossing page
0184     const double DPHICUTPAGE    = 30.;
0185     struct Phys {
0186         TH1D    *R2yy;                          // R2(dy) from R2(y1,y2)
0187         TH1D    *pU[2], *pC[2], *nU[2], *nC[2]; // [0] R2(dy), [1] R2(dphi): whole range / narrow window, uncorrected / corrected
0188         TH1D    *qU[2], *qC[2];                 // [0] C2(dy), [1] C2(dphi): sums over the valid bins, uncorrected / corrected
0189         TH2D    *dC2;                           // C2C - C2 (full stats; bins valid in both)
0190         TH1D    *dC2dy, *dC2dphi;               // its sums over dphi / over dy
0191         double  I, eI, Sfull, frac, efrac;      // integral of dC2 (pairs per event), sum rho2(S), I/S; subgroup errors
0192         std::vector<double> Ic, fracc, efracc;  // per chunk
0193         TH1D    *hy1, *hy2;                     // nevt-weighted averages over chunks
0194         double  fillin;                         // fraction of chunk-bin values taken from full stats in the projections
0195         TH1D    *CQ, *CQ5, *CQ20;               // C(Q) at 1 MeV (written), 5 MeV (display) and 20 MeV (comparisons), ratio of zvtx sums
0196         TH2D    *CQKT, *CQKT10;                 // C(Q,kT) at 2 MeV and 10 MeV (display)
0197         TH2D    *CQvarPT, *QPTsib, *QPTmix;     // ir2=2: C(Q,pT_pair) and its zvtx sums (variable width in Q)
0198         TH1D    *CQvar;                         // ir2=2: C(Q) over all pT_pair (variable width in Q)
0199         TH1D    *MinvS, *MinvM;                 // hMinvFS/FM plain sums (counts)
0200         Phys(){ R2yy=0; for (int k=0;k<2;k++){ pU[k]=pC[k]=nU[k]=nC[k]=qU[k]=qC[k]=0; } dC2=0; dC2dy=dC2dphi=0; I=eI=Sfull=frac=efrac=0; hy1=hy2=0; fillin=0; CQ=CQ5=CQ20=0; CQKT=CQKT10=0; MinvS=MinvM=0; }
0201     };
0202     std::vector<Phys> ph(NPairTypes);
0203     std::vector<long> fillN(NPairTypes,0), fillT(NPairTypes,0);
0204     std::vector<int> chunkRun(nlists,-1);
0205     {
0206         std::ifstream fs((PROJ+"lists/"+SET+"/lists_summary.txt").Data());
0207         std::string line; int nread=0;
0208         while (std::getline(fs,line)){
0209             std::istringstream is(line); int il=-1, run=-1; std::string tok;
0210             if (!(is>>il>>tok>>run) || tok!="run") continue;
0211             if (il>=0 && il<nlists){ chunkRun[il] = run; ++nread; }
0212         }
0213         cout<<"corral::Finalize -- step 1b: runs for "<<nread<<" of "<<nlists<<" chunks from lists_summary.txt"<<endl;
0214     }
0215     //
0216     for (int ipaty=0;ipaty<NPairTypes;ipaty++){
0217         //
0218         //---- plain sums (non-CF names)
0219         for (int ih=0;ih<NFINALIZEHISTS;ih++){
0220             TString base    = TString(FinalizeHists[ih]);
0221             int kd,ir;
0222             if (ParseCF(base,kd,ir)) continue;
0223             TString name    = Form("%s_%d",base.Data(),ipaty);
0224             TH1 *hsum       = 0;
0225             for (int ic=0;ic<nlists;ic++){
0226                 TH1 *h  = (TH1*)fch[ic]->Get(name.Data());
0227                 if (!h){ cout<<"corral::Finalize -- "<<name<<" missing in chunk "<<ic<<", exit"<<endl; exit(1); }
0228                 if (!hsum){ hsum = (TH1*)h->Clone(name.Data()); hsum->SetDirectory(0); }
0229                 else        hsum->Add(h);
0230             }
0231             result[name]    = hsum;
0232         }
0233         //---- hy1, hy2 are per event (dN/dy): average over chunks weighted by the chunk's events (hmult entries)
0234         for (int iyh=0;iyh<2;iyh++){
0235             TH1D *avg   = 0; double sw = 0.;
0236             for (int ic=0;ic<nlists;ic++){
0237                 TH1D *h     = (TH1D*)fch[ic]->Get(Form("hy%d_%d",iyh+1,ipaty));
0238                 TH1 *hm     = (TH1*)fch[ic]->Get(Form("hmult_%d",ipaty));
0239                 if (!h || !hm){ cout<<"corral::Finalize -- hy"<<iyh+1<<"_"<<ipaty<<" or hmult missing in chunk "<<ic<<", exit"<<endl; exit(1); }
0240                 double w    = hm->GetEntries();
0241                 if (!avg){ avg = (TH1D*)h->Clone(Form("hy%d_%d",iyh+1,ipaty)); avg->SetDirectory(0); avg->Reset(); }
0242                 avg->Add(h,w); sw += w;
0243                 delete h; delete hm;
0244             }
0245             if (sw>0.) avg->Scale(1./sw);
0246             if (iyh==0) ph[ipaty].hy1 = avg; else ph[ipaty].hy2 = avg;
0247         }
0248         if (!anyCF) continue;
0249         //
0250         //---- per-chunk, per-zvtx-bin bookkeeping: nevt (hnevtz), rho2(M) normalization (hdenomz), hzvtx
0251         TH1D *hn0   = (TH1D*)fch[0]->Get(Form("hnevtz_%d",ipaty));
0252         if (!hn0){ cout<<"corral::Finalize -- hnevtz_"<<ipaty<<" missing (chunks made before README_Finalize step 1b?), exit"<<endl; exit(1); }
0253         const int NZ    = hn0->GetNbinsX();
0254         if (NZ>64){ cout<<"corral::Finalize -- NZ="<<NZ<<" > 64, exit"<<endl; exit(1); }
0255         std::vector<std::vector<double>> nc(nlists,std::vector<double>(NZ,0.)), dc(nlists,std::vector<double>(NZ,0.));
0256         std::vector<std::vector<double>> zc(nlists,std::vector<double>(NZ,0.));
0257         std::vector<double> zentc(nlists,0.), nevc(nlists,0.);
0258         for (int ic=0;ic<nlists;ic++){
0259             TH1D *hn    = (TH1D*)fch[ic]->Get(Form("hnevtz_%d",ipaty));
0260             TH1D *hd    = (TH1D*)fch[ic]->Get(Form("hdenomz_%d",ipaty));
0261             TH1D *hz    = (TH1D*)fch[ic]->Get(Form("hzvtx_%d",ipaty));
0262             for (int z=0;z<NZ;z++){
0263                 nc[ic][z]   = hn->GetBinContent(z+1);
0264                 dc[ic][z]   = hd->GetBinContent(z+1);
0265                 zc[ic][z]   = hz->GetBinContent(z+1);
0266                 nevc[ic]    += nc[ic][z];
0267             }
0268             zentc[ic]   = hz->GetEntries();
0269             //---- README_SplitTracks573 sec 30: the physics Zvtx range. A slice outside abs(zvtx) < valFinZMax gets no events
0270             //---- and no weight, so every step below (the CF chain, masks, C(Q)) leaves it out; its maps stay in the chunks.
0271             for (int z=0;z<NZ;z++){
0272                 double lo = hz->GetXaxis()->GetBinLowEdge(z+1), hi = hz->GetXaxis()->GetBinUpEdge(z+1);
0273                 if (lo < -valFinZMax-1e-6 || hi > valFinZMax+1e-6){
0274                     zentc[ic]   -= zc[ic][z];
0275                     nevc[ic]    -= nc[ic][z];
0276                     nc[ic][z] = 0.; dc[ic][z] = 0.; zc[ic][z] = 0.;
0277                 }
0278             }
0279             if (ic==0 && ipaty==0){
0280                 int nin=0; for (int z=0;z<NZ;z++) if (zc[0][z]>0.) ++nin;
0281                 cout<<"corral::Finalize -- physics Zvtx range abs(zvtx) < "<<valFinZMax<<" cm: "<<nin<<" of "<<NZ<<" slices (chunk 0)"<<endl;
0282             }
0283             delete hn; delete hd; delete hz;
0284         }
0285         double nevtAll  = 0.;   // events of this pairtype, all chunks (sum of hnevtz)
0286         for (int ic=0;ic<nlists;ic++) nevtAll += nevc[ic];
0287         //---- C(Q) (README_Finalize step 4; user 2026-09-26: RATIO OF SUMS over zvtx bins, standard femtoscopy):
0288         //----   the mixed denominator is still formed within each zvtx bin, in expected-sibling-pair units:
0289         //----   D_z = (sum_c df_c*Qmix_c,z) * N_z/ND_z  (Qmix in the chunk files is already /df_c; N_z = sum nevt,
0290         //----   ND_z = sum nevt*df over the chunks where z was computed), C(Q) = sum_z Qsib_z / sum_z D_z.
0291         //----   Chunk pieces S_c = sum_z Qsib_c,z and D_c = sum_z df_c*Qmix_c,z*N_z/ND_z add up exactly, so the
0292         //----   error is the ratio-estimator subgroup error sigma^2 = n/(n-1) sum_c (S_c - C D_c)^2 / (sum_c D_c)^2
0293         //----   (no per-chunk ratio: sparse chunks are fine). S=0 over D>0: valid zero, one-pair error 1/D.
0294         {
0295             Phys &P = ph[ipaty];
0296             double Nz[64], NDz[64]; bool okq[64];
0297             for (int z=0;z<NZ;z++){
0298                 Nz[z]=0.; NDz[z]=0.;
0299                 for (int ic=0;ic<nlists;ic++) if (nc[ic][z]>=2. && dc[ic][z]>0.){ Nz[z] += nc[ic][z]; NDz[z] += nc[ic][z]*dc[ic][z]; }
0300                 okq[z]  = (Nz[z]>=2. && NDz[z]>0.);
0301             }
0302             //---- ir2=2 (README_meanpT_from_CQ.md sec 11): the same ratio of zvtx sums from the (Q,pT_pair) densities.
0303             //----   perEvent: inputs are rho2(S) = sib/nevt and rho2(M) = mix/nevt/denomfactor (CalcRm::Calculate), so
0304             //----   x nevt (= hnevtz, the same hmult entries) gives back exactly what hQsib, hQmix hold. projQ: sum over
0305             //----   pT_pair (incl. under/overflow) first, for C(Q) over all pT_pair. keepSD names: sname_k, dname_k.
0306             auto buildCQ    = [&](const char* bs, const char* bm, int rebin, const char* outname, const char* title, bool keepSD=false,
0307                                   bool perEvent=false, bool projQ=false, const char* sname="hQsib", const char* dname="hQmix") -> TH1* {
0308                 TH1 *t0 = (TH1*)fch[0]->Get(Form("%s_z00_%d",bs,ipaty));
0309                 if (!t0){ cout<<"corral::Finalize -- "<<bs<<"_z00_"<<ipaty<<" missing, exit"<<endl; exit(1); }
0310                 if (projQ){ TH1 *p0 = ((TH2*)t0)->ProjectionX(Form("%s_t0p",outname),0,-1); p0->SetDirectory(0); delete t0; t0 = p0; }
0311                 TH1 *out    = (TH1*)t0->Clone(outname); out->SetDirectory(0); delete t0;
0312                 if (rebin>1){ if (out->GetDimension()==2) ((TH2*)out)->RebinX(rebin); else out->Rebin(rebin); }
0313                 out->Reset(); out->SetTitle(title);
0314                 int nce = out->GetNcells();
0315                 std::vector<std::vector<double>> Sc(nlists,std::vector<double>(nce,0.)), Dc(nlists,std::vector<double>(nce,0.));
0316                 for (int ic=0;ic<nlists;ic++){
0317                     for (int z=0;z<NZ;z++){
0318                         if (!okq[z] || !(nc[ic][z]>=2. && dc[ic][z]>0.)) continue;
0319                         TH1 *hs = (TH1*)fch[ic]->Get(Form("%s_z%02d_%d",bs,z,ipaty));
0320                         TH1 *hm = (TH1*)fch[ic]->Get(Form("%s_z%02d_%d",bm,z,ipaty));
0321                         if (!hs || !hm){ cout<<"corral::Finalize -- "<<bs<<"/"<<bm<<"_z"<<z<<"_"<<ipaty<<" missing in chunk "<<ic<<", exit"<<endl; exit(1); }
0322                         hs->SetDirectory(0); hm->SetDirectory(0);
0323                         if (perEvent){ hs->Scale(nc[ic][z]); hm->Scale(nc[ic][z]); }
0324                         if (projQ){
0325                             TH1 *ps = ((TH2*)hs)->ProjectionX(Form("%s_ps",outname),0,-1), *pm = ((TH2*)hm)->ProjectionX(Form("%s_pm",outname),0,-1);
0326                             ps->SetDirectory(0); pm->SetDirectory(0); delete hs; delete hm; hs = ps; hm = pm;
0327                         }
0328                         if (rebin>1){
0329                             if (hs->GetDimension()==2){ ((TH2*)hs)->RebinX(rebin); ((TH2*)hm)->RebinX(rebin); }
0330                             else { hs->Rebin(rebin); hm->Rebin(rebin); }
0331                         }
0332                         double w    = dc[ic][z]*Nz[z]/NDz[z];
0333                         for (int ib=0;ib<nce;ib++){ Sc[ic][ib] += hs->GetBinContent(ib); Dc[ic][ib] += w*hm->GetBinContent(ib); }
0334                         delete hs; delete hm;
0335                     }
0336                 }
0337                 //---- keepSD (user 2026-10-06): also keep the numerator and denominator, hQsib_k = sum_z Qsib (sibling pairs,
0338                 //---- sqrt(N) errors) and hQmix_k = sum_z D_z (mixed, expected-sibling-pair units), so peak yields in pairs
0339                 //---- = sum (C - base) * hQmix need no pass over the chunks
0340                 TH1 *hS = 0, *hD = 0;
0341                 if (keepSD){
0342                     TString pn  = TString(out->GetTitle()); pn.Remove(pn.First(','));   // the pair label (title is a Form() buffer, long reused here)
0343                     TString ot  = out->GetTitle();
0344                     TString ax  = (out->GetDimension()==2) ? TString(ot(ot.First(';'),ot.Length())) : TString(";Q_{inv} (GeV);pairs");  // 2D: the x;y axes of C
0345                     hS  = (TH1*)out->Clone(Form("%s_%d",sname,ipaty)); hS->SetDirectory(0);
0346                     hD  = (TH1*)out->Clone(Form("%s_%d",dname,ipaty)); hD->SetDirectory(0);
0347                     double nevQ = 0.; for (int z=0;z<NZ;z++) if (okq[z]) nevQ += Nz[z];
0348                     hS->SetTitle(Form("%s, #Sigma_{z} sibling pairs (N_{evt} = %.0f)%s",pn.Data(),nevQ,ax.Data()));
0349                     hD->SetTitle(Form("%s, #Sigma_{z} mixed pairs, expected-sibling units (N_{evt} = %.0f)%s",pn.Data(),nevQ,ax.Data()));
0350                     if (TString(sname)=="hQsib") extraW.push_back(new TParameter<double>(Form("nevtQ_%d",ipaty),nevQ));
0351                 }
0352                 for (int ib=0;ib<nce;ib++){
0353                     double S=0, D=0;
0354                     for (int ic=0;ic<nlists;ic++){ S += Sc[ic][ib]; D += Dc[ic][ib]; }
0355                     if (hS){ hS->SetBinContent(ib,S); hS->SetBinError(ib,sqrt(S)); hD->SetBinContent(ib,D); hD->SetBinError(ib,0.); }
0356                     if (D<=0.){ out->SetBinContent(ib,0.); out->SetBinError(ib,0.); continue; }
0357                     double C    = S/D;
0358                     double v    = 0.;
0359                     for (int ic=0;ic<nlists;ic++) v += pow(Sc[ic][ib] - C*Dc[ic][ib],2);
0360                     double e    = sqrt(v*nlists/(nlists-1.))/D;
0361                     if (S<=0.) e = 1./D;
0362                     out->SetBinContent(ib,C); out->SetBinError(ib,e);
0363                 }
0364                 extraW.push_back(out);
0365                 if (hS){ extraW.push_back(hS); extraW.push_back(hD); }
0366                 return out;
0367             };
0368             TString pnm = Form("%s%s",ParticleIDNames[PairTypes_Info[ipaty][0]],ParticleIDNames[PairTypes_Info[ipaty][1]]);
0369             P.CQ    = (TH1D*)buildCQ("hQsib","hQmix",1,Form("hFin_CQ_%d",ipaty),    Form("%s, C(Q) = #Sigma_{z}sib / #Sigma_{z}mix (1 MeV bins);Q_{inv} (GeV);C(Q)",pnm.Data()),true);
0370             P.CQ5   = (TH1D*)buildCQ("hQsib","hQmix",5,Form("hFin_CQ5_%d",ipaty),   Form("%s, C(Q) (5 MeV bins, #pm subgroup err);Q_{inv} (GeV);C(Q)",pnm.Data()));
0371             P.CQ20  = (TH1D*)buildCQ("hQsib","hQmix",20,Form("hFin_CQ20_%d",ipaty), Form("%s, C(Q) (20 MeV bins, #pm subgroup err);Q_{inv} (GeV);C(Q)",pnm.Data()));
0372             //---- 2026-10-08: pair-normalized C(Q) at 20 MeV, hCQ20n_k = hCQ20_k / (1+r2), with 1+r2 = total sibling pairs /
0373             //---- total mixed pairs (expected-sibling units) over ALL Q incl. overflow, same zvtx slices = int rho2 / int rho1rho1
0374             //---- of exactly the pairs in C(Q). Removes the event-normalization offset (baseline -> ~1); 1+r2 is in the title
0375             //---- and in work/onePlusR2_k.
0376             {
0377                 TH1 *hS = 0, *hD = 0;
0378                 for (TObject *o : extraW){ TString n = o->GetName(); if (n==Form("hQsib_%d",ipaty)) hS = (TH1*)o; if (n==Form("hQmix_%d",ipaty)) hD = (TH1*)o; }
0379                 double St = hS ? hS->Integral(0,hS->GetNbinsX()+1) : 0., Dt = hD ? hD->Integral(0,hD->GetNbinsX()+1) : 0.;
0380                 double opr  = (St>0. && Dt>0.) ? St/Dt : 0.;
0381                 TH1D *hn    = (TH1D*)P.CQ20->Clone(Form("hFin_CQ20n_%d",ipaty)); hn->SetDirectory(0);
0382                 if (opr>0.) hn->Scale(1./opr); else hn->Reset();
0383                 hn->SetTitle(Form("%s, C(Q)/(1+r_{2}), 1+r_{2} = %.4f (20 MeV bins, #pm subgroup err);Q_{inv} (GeV);C(Q)/(1+r_{2})",pnm.Data(),opr));
0384                 extraW.push_back(hn);
0385                 extraW.push_back(new TParameter<double>(Form("onePlusR2_%d",ipaty),opr));
0386                 cout<<"corral::Finalize -- pairtype "<<ipaty<<" "<<pnm<<": 1+r2 (C(Q) pairs, all Q) = "<<opr<<endl;
0387             }
0388             P.CQKT  = (TH2D*)buildCQ("hQsibKT","hQmixKT",1,Form("hFin_CQKT_%d",ipaty),  Form("%s, C(Q,k_{T});Q_{inv} (GeV);k_{T} (GeV)",pnm.Data()));
0389             //---- ir2=2 (README_meanpT_from_CQ.md sec 11): C(Q,pT_pair) and C(Q) over all pT_pair, variable width in Q;
0390             //---- the sums behind them kept as work/hQPTsib_k (sibling pairs), work/hQPTmix_k (mixed, expected-sibling units)
0391             P.CQvarPT   = (TH2D*)buildCQ("hrho2_2","hrho1rho1_2",1,Form("hFin_CQvarPT_%d",ipaty),
0392                             Form("%s, C(Q,p_{T,pair}) (variable width binning in Q);Q_{inv} (GeV);p_{T,pair} (GeV/c)",pnm.Data()),true,true,false,"hQPTsib","hQPTmix");
0393             P.CQvar     = (TH1D*)buildCQ("hrho2_2","hrho1rho1_2",1,Form("hFin_CQvar_%d",ipaty),
0394                             Form("%s, C(Q) all p_{T,pair} (variable width binning, #pm subgroup err);Q_{inv} (GeV);C(Q)",pnm.Data()),false,true,true);
0395             for (TObject *o : extraW){ TString n = o->GetName(); if (n==Form("hQPTsib_%d",ipaty)) P.QPTsib = (TH2D*)o; if (n==Form("hQPTmix_%d",ipaty)) P.QPTmix = (TH2D*)o; }
0396             P.CQKT10= (TH2D*)buildCQ("hQsibKT","hQmixKT",5,Form("hFin_CQKT10_%d",ipaty),Form("%s, C(Q,k_{T}) (10 MeV bins);Q_{inv} (GeV);k_{T} (GeV)",pnm.Data()));
0397             //---- Minv near threshold: plain sums (counts)
0398             for (int im=0;im<2;im++){
0399                 TH1D *sum   = 0;
0400                 for (int ic=0;ic<nlists;ic++){
0401                     TH1D *h = (TH1D*)fch[ic]->Get(Form("%s_%d",im==0?"hMinvFS":"hMinvFM",ipaty));
0402                     if (!h) continue;
0403                     if (!sum){ sum = (TH1D*)h->Clone(Form("hFin_%s_%d",im==0?"MinvFS":"MinvFM",ipaty)); sum->SetDirectory(0); }
0404                     else sum->Add(h);
0405                     delete h;
0406                 }
0407                 if (im==0) P.MinvS = sum; else P.MinvM = sum;
0408                 if (sum) extraW.push_back(sum);
0409             }
0410             cout<<Form("corral::Finalize -- C(Q) pairtype %2d: ratio of zvtx sums built (%d zvtx bins used)",ipaty,(int)std::count(okq,okq+NZ,true))<<endl;
0411         }
0412         std::vector<int> allc(nlists);
0413         for (int ic=0;ic<nlists;ic++) allc[ic] = ic;
0414         //
0415         const double nexpmin    = nexpOverride ? nexpGlobal : (double)PairTypes_Info[ipaty][2];
0416         cout<<"corral::Finalize -- pairtype "<<ipaty<<": N_min = "<<nexpmin<<endl;
0417         int chg1    = ParticleCharge[PairTypes_Info[ipaty][0]];
0418         int chg2    = ParticleCharge[PairTypes_Info[ipaty][1]];
0419         int ds      = CrossingDirtySide(FINALIZE_FIELD,chg1,chg2);
0420         //
0421         for (int ir2=0;ir2<3;ir2++){
0422             //---- does the list want anything from this ir2?
0423             bool want   = false;
0424             for (int ih=0;ih<NFINALIZEHISTS;ih++){ int kd,ir; if (ParseCF(TString(FinalizeHists[ih]),kd,ir) && ir==ir2) want = true; }
0425             if (!want) continue;
0426             const int NK    = (ir2==1) ? 7 : 4;
0427             //
0428             //---- cache this pairtype's per-zvtx-bin maps: read once, used by the full chain and all subgroups
0429             std::vector<std::vector<TH2D*>> cS(nlists,std::vector<TH2D*>(NZ,(TH2D*)0)), cM(nlists,std::vector<TH2D*>(NZ,(TH2D*)0));
0430             for (int ic=0;ic<nlists;ic++){
0431                 for (int z=0;z<NZ;z++){
0432                     if (!(nc[ic][z]>=2. && dc[ic][z]>0.)) continue;
0433                     cS[ic][z]   = (TH2D*)fch[ic]->Get(Form("hrho2_%d_z%02d_%d"    ,ir2,z,ipaty));
0434                     cM[ic][z]   = (TH2D*)fch[ic]->Get(Form("hrho1rho1_%d_z%02d_%d",ir2,z,ipaty));
0435                     if (!cS[ic][z] || !cM[ic][z]){ cout<<"corral::Finalize -- per-zvtx maps missing in chunk "<<ic<<", exit"<<endl; exit(1); }
0436                     cS[ic][z]->SetDirectory(0); cM[ic][z]->SetDirectory(0);
0437                 }
0438             }
0439             TH2D *tmpl  = (TH2D*)fch[0]->Get(Form("hrho2_%d_z%02d_%d",ir2,0,ipaty));    // binning template
0440             if (!tmpl){ cout<<"corral::Finalize -- hrho2_"<<ir2<<"_z00_"<<ipaty<<" missing in chunk 0, exit"<<endl; exit(1); }
0441             tmpl    = (TH2D*)tmpl->Clone(Form("tmpl_%d_%d",ir2,ipaty)); tmpl->SetDirectory(0); tmpl->Reset();
0442             //---- sec 38: the bin's largest full-stats rho2(M) over the zvtx slices used (reference of the acceptance-edge rule)
0443             TH2D *mRef  = (TH2D*)tmpl->Clone(Form("mRef_%d_%d",ir2,ipaty)); mRef->SetDirectory(0);
0444             for (int z=0;z<NZ;z++){
0445                 double nd = 0.; for (int ic=0;ic<nlists;ic++) if (cM[ic][z]) nd += nc[ic][z]*dc[ic][z];
0446                 if (nd<=0.) continue;
0447                 TH2D *mz    = (TH2D*)tmpl->Clone("mz_tmp"); mz->SetDirectory(0);
0448                 for (int ic=0;ic<nlists;ic++) if (cM[ic][z]) mz->Add(cM[ic][z], nc[ic][z]*dc[ic][z]/nd);
0449                 for (int ib=0;ib<mz->GetNcells();ib++) if (mz->GetBinContent(ib)>mRef->GetBinContent(ib)) mRef->SetBinContent(ib,mz->GetBinContent(ib));
0450                 delete mz;
0451             }
0452             //
0453             //---- the CF chain on a set of chunks (all = full stats, step 2; one = a subgroup, step 3)
0454             //----   own mode (fx==0): zvtx weights, zvtx-bin use and validity masks from these chunks;
0455             //----     VZ[z] = 1 where rho2(M)>0 and N_exp = nevt*rho2(M) >= nmin (nmin=0: CalcRm's M>0 rule)
0456             //----   fixed mode (fx = the full-stats set): the full-stats weights and masks VZ/MK, so a
0457             //----     subgroup is the same linear combination as the full-stats value. out.V/VC then say
0458             //----     where this chunk is USABLE: every zvtx bin (and, corrected, every source bin) that the
0459             //----     full-stats value uses is defined in this chunk (computed, rho2(M)>0).
0460             auto chain  = [&](const std::vector<int>& chs, bool report, const TString& tag, double nmin, const CFSet* fx, bool keepMasks) -> CFSet {
0461                 double N[64], ND[64], wz[64]; bool okz[64], okn[64];
0462                 double zent = 0.;
0463                 for (int z=0;z<NZ;z++){ N[z]=0.; ND[z]=0.; wz[z]=0.; }
0464                 for (int ic : chs){
0465                     zent    += zentc[ic];
0466                     for (int z=0;z<NZ;z++){
0467                         wz[z]   += zc[ic][z];
0468                         if (nc[ic][z]>=2. && dc[ic][z]>0.){ N[z] += nc[ic][z]; ND[z] += nc[ic][z]*dc[ic][z]; }
0469                     }
0470                 }
0471                 for (int z=0;z<NZ;z++){ wz[z] = (zent>0.) ? wz[z]/zent : 0.; okn[z] = (N[z]>=2. && ND[z]>0.); okz[z] = okn[z]; }
0472                 if (fx) for (int z=0;z<NZ;z++){ wz[z] = fx->wz[z]; okz[z] = fx->okz[z]; }
0473                 std::vector<TH2D*> S(NZ,0), M(NZ,0), C2(NZ,0), R2(NZ,0), SC(NZ,0), C2C(NZ,0), R2C(NZ,0), MK(NZ,0), VZ(NZ,0), BADU(NZ,0), BADC(NZ,0);
0474                 long nBelow=0, nBelowS=0;       // bin-zvtx with rho2(M)>0 but N_exp<nmin (own mode)
0475                 long nEdge=0, nEdgeS=0;         // sec 38: bin-zvtx left out as acceptance edges (rho2(M) < mFrac x mRef; own mode, nmin>0)
0476                 for (int z=0;z<NZ;z++){
0477                     S[z]    = (TH2D*)tmpl->Clone(Form("S%d_%d_%d%s",ir2,z,ipaty,tag.Data())); S[z]->SetDirectory(0);
0478                     M[z]    = (TH2D*)tmpl->Clone(Form("M%d_%d_%d%s",ir2,z,ipaty,tag.Data())); M[z]->SetDirectory(0);
0479                     for (int ic : chs){
0480                         if (!cS[ic][z]) continue;
0481                         S[z]->Add(cS[ic][z], nc[ic][z]);
0482                         M[z]->Add(cM[ic][z], nc[ic][z]*dc[ic][z]);
0483                     }
0484                     if (okn[z]){ S[z]->Scale(1./N[z]); M[z]->Scale(1./ND[z]); }
0485                     //---- validity mask of this zvtx bin
0486                     if (fx){
0487                         VZ[z]   = (TH2D*)fx->vz[z]->Clone(Form("VZ_%d_%d_%d%s",ir2,z,ipaty,tag.Data())); VZ[z]->SetDirectory(0);
0488                     } else {
0489                         VZ[z]   = (TH2D*)tmpl->Clone(Form("VZ_%d_%d_%d%s",ir2,z,ipaty,tag.Data())); VZ[z]->SetDirectory(0);
0490                         for (int ib=0;ib<VZ[z]->GetNcells();ib++){
0491                             double m    = M[z]->GetBinContent(ib);
0492                             if (!okn[z] || m<=0.) continue;
0493                             if (nmin>0. && mFrac>0. && m<mFrac*mRef->GetBinContent(ib)){
0494                                 if (!VZ[z]->IsBinUnderflow(ib) && !VZ[z]->IsBinOverflow(ib)){ ++nEdge; if (S[z]->GetBinContent(ib)>0.) ++nEdgeS; }
0495                                 continue;
0496                             }
0497                             if (N[z]*m>=nmin) VZ[z]->SetBinContent(ib,1.);
0498                             else if (!VZ[z]->IsBinUnderflow(ib) && !VZ[z]->IsBinOverflow(ib)){ ++nBelow; if (S[z]->GetBinContent(ib)>0.) ++nBelowS; }
0499                         }
0500                     }
0501                     //---- fixed mode: where the full-stats value uses this zvtx bin but this chunk has nothing
0502                     if (fx){
0503                         BADU[z] = (TH2D*)tmpl->Clone(Form("BADU_%d_%d_%d%s",ir2,z,ipaty,tag.Data())); BADU[z]->SetDirectory(0);
0504                         for (int ib=0;ib<BADU[z]->GetNcells();ib++){
0505                             if (VZ[z]->GetBinContent(ib)<=0.) continue;
0506                             if (!okn[z] || M[z]->GetBinContent(ib)<=0.) BADU[z]->SetBinContent(ib,1.);
0507                         }
0508                     }
0509                     //---- C2 and R2 per zvtx bin: CalcRm::Calculate's rules (defined iff M>0; S=0 is a valid zero)
0510                     C2[z]   = (TH2D*)S[z]->Clone(Form("C2_%d_%d_%d%s",ir2,z,ipaty,tag.Data())); C2[z]->SetDirectory(0); C2[z]->Reset();
0511                     R2[z]   = (TH2D*)S[z]->Clone(Form("R2_%d_%d_%d%s",ir2,z,ipaty,tag.Data())); R2[z]->SetDirectory(0); R2[z]->Reset();
0512                     if (okn[z]){
0513                         for (int ibin=1;ibin<=S[z]->GetNbinsX();ibin++){
0514                             for (int jbin=1;jbin<=S[z]->GetNbinsY();jbin++){
0515                                 double valn = S[z]->GetBinContent(ibin,jbin), valne = S[z]->GetBinError(ibin,jbin);
0516                                 double vald = M[z]->GetBinContent(ibin,jbin), valde = M[z]->GetBinError(ibin,jbin);
0517                                 double valC2,valC2e,valR2,valR2e;
0518                                 if (vald>0. && valn<=0.){
0519                                     double valne1   = 1./N[z];
0520                                     valC2   = -vald;    valC2e  = sqrt(valne1*valne1 + valde*valde);
0521                                     valR2   = -1.;      valR2e  = valne1/vald;
0522                                 } else if (valn>0. && vald>0.){
0523                                     valC2   = valn - vald;  valC2e  = sqrt(valne*valne + valde*valde);
0524                                     valR2   = valn/vald;    valR2e  = valR2*sqrt(pow(valde/vald,2)+pow(valne/valn,2));
0525                                     valR2   -= 1.;
0526                                 } else {
0527                                     valC2   = 0.; valC2e = 0.; valR2 = -1.; valR2e = 0.;
0528                                     //---- fixed mode (a subgroup): C2 = S - M stays linear where this chunk has M=0 but
0529                                     //---- the full-stats mask uses the bin, so S, M, C2, C2C of a chunk are defined
0530                                     //---- everywhere (only R2 needs the chunk to be usable)
0531                                     if (fx){ valC2 = valn - vald; valC2e = sqrt(valne*valne + valde*valde); }
0532                                 }
0533                                 C2[z]->SetBinContent(ibin,jbin,valC2); C2[z]->SetBinError(ibin,jbin,valC2e);
0534                                 R2[z]->SetBinContent(ibin,jbin,valR2); R2[z]->SetBinError(ibin,jbin,valR2e);
0535                             }
0536                         }
0537                     }
0538                     //---- crossing correction per zvtx bin (ir2=1), as in Calculate; validity = the VZ mask
0539                     if (ir2==1){
0540                         TH2D* src[3]    = {S[z], C2[z], R2[z]};
0541                         TH2D** dst[3]   = {&SC[z], &C2C[z], &R2C[z]};
0542                         for (int k=0;k<3;k++){
0543                             *dst[k] = (TH2D*)src[k]->Clone(Form("%s_C",src[k]->GetName())); (*dst[k])->SetDirectory(0);
0544                             if (okz[z] && chunkCross) CrossingCorrect(*dst[k],ds,VZ[z]);
0545                         }
0546                         if (fx){
0547                             MK[z]   = (TH2D*)fx->mk[z]->Clone(Form("MK_%d_%d%s",z,ipaty,tag.Data())); MK[z]->SetDirectory(0);
0548                             //---- a corrected bin draws on itself and/or its dphi mirror exactly as a map does: correct
0549                             //---- the flags with an all-valid mask (near: dirty takes the clean's flag; away: either flags both)
0550                             BADC[z] = (TH2D*)BADU[z]->Clone(Form("BADC_%d_%d%s",z,ipaty,tag.Data())); BADC[z]->SetDirectory(0);
0551                             TH2D *ones  = (TH2D*)tmpl->Clone("ones_tmp"); ones->SetDirectory(0);
0552                             for (int ib=0;ib<ones->GetNcells();ib++) ones->SetBinContent(ib,1.);
0553                             if (okz[z] && chunkCross) CrossingCorrect(BADC[z],ds,ones);
0554                             delete ones;
0555                         } else {
0556                             MK[z]   = (TH2D*)VZ[z]->Clone(Form("MK_%d_%d%s",z,ipaty,tag.Data())); MK[z]->SetDirectory(0);
0557                             if (okz[z] && chunkCross) CrossingCorrect(MK[z],ds,MK[z]);
0558                         }
0559                     }
0560                 }
0561                 //---- empty-denominator report at full statistics (the squawk belongs here, README_Crossing sec 0)
0562                 if (report){
0563                     long nMz=0, nSMz=0; int nbinAny=0;
0564                     for (int ib=0;ib<S[0]->GetNcells();ib++){
0565                         if (S[0]->IsBinUnderflow(ib) || S[0]->IsBinOverflow(ib)) continue;
0566                         int nzb=0;
0567                         for (int z=0;z<NZ;z++){
0568                             if (!okz[z] || wz[z]<=0.) continue;
0569                             if (M[z]->GetBinContent(ib)>0.) continue;
0570                             ++nzb; ++nMz; if (S[z]->GetBinContent(ib)>0.) ++nSMz;
0571                         }
0572                         if (nzb>0) ++nbinAny;
0573                     }
0574                     cout<<Form("corral::Finalize -- pairtype %2d ir2=%d: %d bins have rho2(M)=0 in >=1 used zvtx bin (%ld bin-zvtx, %ld of them with rho2(S)>0); %ld more bin-zvtx below N_exp>=%g (%ld with rho2(S)>0) left out",
0575                         ipaty,ir2,nbinAny,nMz,nSMz,nBelow,nmin,nBelowS)<<endl;
0576                     cout<<Form("corral::Finalize -- pairtype %2d ir2=%d: %ld bin-zvtx left out as acceptance edges (rho2(M) < %.2f x the bin's largest over the slices; %ld with rho2(S)>0)",
0577                         ipaty,ir2,nEdge,mFrac,nEdgeS)<<endl;
0578                 }
0579                 //---- Zvtx averages (masks: VZ uncorrected, MK corrected)
0580                 auto zavg   = [&](std::vector<TH2D*>& v, std::vector<TH2D*>& mk, double sentinel, const char* nm){
0581                     TH2D *avg   = (TH2D*)v[0]->Clone(nm); avg->SetDirectory(0);
0582                     std::vector<TH1*> sl(NZ), ms(NZ);
0583                     for (int z=0;z<NZ;z++){ sl[z] = v[z]; ms[z] = mk[z]; }
0584                     ZvtxAverageValid(avg,sl.data(),ms.data(),wz,okz,NZ,sentinel);
0585                     return avg;
0586                 };
0587                 //---- validity of the average: some used zvtx bin valid; fixed mode also: no used zvtx bin BAD
0588                 auto zvalid = [&](std::vector<TH2D*>& mk, std::vector<TH2D*>& bad, const char* nm){
0589                     TH2D *hv    = (TH2D*)tmpl->Clone(nm); hv->SetDirectory(0); hv->Reset();
0590                     for (int ib=0;ib<hv->GetNcells();ib++){
0591                         bool any=false, isbad=false;
0592                         for (int z=0;z<NZ;z++){
0593                             if (!okz[z] || wz[z]<=0. || mk[z]->GetBinContent(ib)<=0.) continue;
0594                             any = true;
0595                             if (fx && bad[z]->GetBinContent(ib)>0.) isbad = true;
0596                         }
0597                         if (any && !isbad) hv->SetBinContent(ib,1.);
0598                     }
0599                     return hv;
0600                 };
0601                 CFSet out;
0602                 std::vector<TH2D*>* src[7]  = {&S,&M,&C2,&R2,&SC,&C2C,&R2C};
0603                 for (int k=0;k<NK;k++){
0604                     double sentinel = (k==3 || k==6) ? -1. : 0.;
0605                     out.h[k]    = zavg(*src[k], (k<4) ? VZ : MK, sentinel, Form("%s_%d%s",Form(KINDNAME[k],ir2),ipaty,tag.Data()));
0606                 }
0607                 out.V   = zvalid(VZ, BADU, Form("V_%d_%d%s",ir2,ipaty,tag.Data()));
0608                 if (ir2==1) out.VC  = zvalid(MK, BADC, Form("VC_%d_%d%s",ir2,ipaty,tag.Data()));
0609                 if (keepMasks){
0610                     out.vz  = VZ; out.mk = MK;
0611                     for (int z=0;z<NZ;z++){ out.wz[z] = wz[z]; out.okz[z] = okz[z]; }
0612                     VZ.assign(NZ,0); MK.assign(NZ,0);
0613                 }
0614                 for (int z=0;z<NZ;z++){ delete S[z]; delete M[z]; delete C2[z]; delete R2[z]; delete SC[z]; delete C2C[z]; delete R2C[z]; delete MK[z]; delete VZ[z]; delete BADU[z]; delete BADC[z]; }
0615                 return out;
0616             };
0617             //
0618             //---- step 2: full statistics, with the N_exp>=nmin validity; its masks and weights are kept
0619             CFSet full  = chain(allc, true, TString(""), nexpmin, 0, true);
0620             double chk1d[7]={0,0,0,0,0,0,0}, chk1e[7]={0,0,0,0,0,0,0}; int chk1miss[7]={0,0,0,0,0,0,0};
0621             //---- check 1: a one-chunk chain with the chunks' own rule (nmin=0) must reproduce that chunk's
0622             //---- saved Zvtx averages from CalcRm::Calculate
0623             for (int ic=0;ic<nlists;ic++){
0624                 CFSet own   = chain(std::vector<int>(1,ic), false, TString(Form("_own%02d",ic)), 0., 0, false);
0625                 for (int k=0;k<NK;k++){
0626                     int kc      = (!chunkCross && k>=4) ? ((k==4) ? 0 : k-3) : k;   // nocross chunks (rho2C,C2C,R2C -> rho2,C2,R2): C must equal the chunk's plain map
0627                     TString nm  = Form("%s_%d",Form(KINDNAME[kc],ir2),ipaty);
0628                     TH2D *hc    = (TH2D*)fch[ic]->Get(nm.Data());
0629                     if (!hc){ ++chk1miss[k]; continue; }
0630                     for (int ib=0;ib<hc->GetNcells();ib++){
0631                         if (hc->IsBinUnderflow(ib) || hc->IsBinOverflow(ib)) continue;
0632                         chk1d[k]    = std::max(chk1d[k], fabs(own.h[k]->GetBinContent(ib) - hc->GetBinContent(ib)));
0633                         chk1e[k]    = std::max(chk1e[k], fabs(own.h[k]->GetBinError(ib)   - hc->GetBinError(ib)));
0634                     }
0635                     delete hc;
0636                 }
0637                 own.Delete();
0638             }
0639             for (int k=0;k<NK;k++){
0640                 cout<<Form("corral::Finalize -- step 3 check 1 %-16s one-chunk chain (own rule, N_exp>=0) vs the chunk's own map: max |diff| content %.2e uncertainty %.2e (%d chunks missing it)",
0641                     Form("%s_%d",Form(KINDNAME[k],ir2),ipaty),chk1d[k],chk1e[k],chk1miss[k])<<endl;
0642                 chk1d[k]=0.; chk1e[k]=0.; chk1miss[k]=0;
0643             }
0644             //---- check 2: the fixed-mask chain on ALL chunks must reproduce the full-stats result exactly
0645             {
0646                 CFSet fx2   = chain(allc, false, TString("_fx2"), nexpmin, &full, false);
0647                 int nd=0;
0648                 for (int k=0;k<NK;k++){
0649                     for (int ib=0;ib<fx2.h[k]->GetNcells();ib++){
0650                         if (fx2.h[k]->GetBinContent(ib)!=full.h[k]->GetBinContent(ib) || fx2.h[k]->GetBinError(ib)!=full.h[k]->GetBinError(ib)) ++nd;
0651                     }
0652                 }
0653                 for (int ib=0;ib<full.V->GetNcells();ib++){
0654                     if (fx2.V->GetBinContent(ib)!=full.V->GetBinContent(ib)) ++nd;
0655                     if (ir2==1 && fx2.VC->GetBinContent(ib)!=full.VC->GetBinContent(ib)) ++nd;
0656                 }
0657                 cout<<Form("corral::Finalize -- step 3 check 2 pairtype %2d ir2=%d: fixed-mask chain on all chunks vs full: %d differing cells (expect 0)",ipaty,ir2,nd)<<endl;
0658                 fx2.Delete();
0659             }
0660             //---- step 3: each chunk alone, with the full-stats masks and weights
0661             std::vector<CFSet> sub(nlists);
0662             for (int ic=0;ic<nlists;ic++) sub[ic] = chain(std::vector<int>(1,ic), false, TString(Form("_sg%02d",ic)), nexpmin, &full, false);
0663             //
0664             //---- step 3: subgroup errors on the full-stats values
0665             for (int k=0;k<NK;k++){
0666                 TH2D *hf    = full.h[k];
0667                 TH2D *vf    = (k<4) ? full.V : full.VC;
0668                 std::vector<double> rsp, rw;    // sg/per-run-style error, weighted/unweighted error
0669                 int nlow=0, nset=0;
0670                 for (int ib=0;ib<hf->GetNcells();ib++){
0671                     if (hf->IsBinUnderflow(ib) || hf->IsBinOverflow(ib)) continue;
0672                     if (vf->GetBinContent(ib)<=0.) continue;            // not valid at full stats: sentinel, error 0
0673                     double sx=0, sw=0, swx=0; int n=0;
0674                     std::vector<double> xv, wv;
0675                     for (int ic=0;ic<nlists;ic++){
0676                         TH2D *vs    = (k<4) ? sub[ic].V : sub[ic].VC;
0677                         bool linear = (k!=3 && k!=6);               // S, M, C2, rho2C, C2C: defined in every chunk (fixed mode)
0678                         if (!linear && vs->GetBinContent(ib)<=0.) continue; // pick (a), R2/R2C only: undefined in this chunk -> skipped
0679                         double x    = sub[ic].h[k]->GetBinContent(ib);
0680                         xv.push_back(x); wv.push_back(nevc[ic]);
0681                         sx += x; sw += nevc[ic]; swx += nevc[ic]*x; ++n;
0682                     }
0683                     double eold = hf->GetBinError(ib);
0684                     double err  = 0.;
0685                     if (n>=NVALIDMIN){
0686                         double m    = sx/n, mw = swx/sw, s2 = 0., sw2d2 = 0.;
0687                         for (size_t j=0;j<xv.size();j++){ s2 += (xv[j]-m)*(xv[j]-m); sw2d2 += wv[j]*wv[j]*(xv[j]-mw)*(xv[j]-mw); }
0688                         err = sqrt(s2/(n-1.))/sqrt((double)n);
0689                         double errw = sqrt(sw2d2)/sw*sqrt(n/(n-1.));    // pick (b): event-weighted version, printed only
0690                         if (err>0.){ rw.push_back(errw/err); if (eold>0.) rsp.push_back(err/eold); }
0691                         ++nset;
0692                     } else {
0693                         ++nlow;
0694                     }
0695                     hf->SetBinError(ib,err);
0696                 }
0697                 cout<<Form("corral::Finalize -- step 3 %-16s subgroup uncertainties in %d bins, %d valid bins with N_valid<%d -> uncertainty 0; median sg/per-run-style uncertainty %.3f, median weighted/unweighted %.4f",
0698                     hf->GetName(),nset,nlow,NVALIDMIN,Median(rsp),Median(rw))<<endl;
0699             }
0700             //---- N_valid maps (same for kinds 0-3; kinds 4-6 use the corrected validity)
0701             for (int iv=0;iv<((ir2==1)?2:1);iv++){
0702                 TH2D *hnv   = (TH2D*)tmpl->Clone(Form("%s_%d_%d",(iv==0)?"hNvalid":"hNvalidC",ir2,ipaty)); hnv->SetDirectory(0); hnv->Reset();
0703                 hnv->SetTitle(Form("%s: N_{valid} subgroups%s;%s;%s",hnv->GetName(),(iv==0)?"":" (crossing-corrected)",tmpl->GetXaxis()->GetTitle(),tmpl->GetYaxis()->GetTitle()));
0704                 for (int ic=0;ic<nlists;ic++) hnv->Add((iv==0) ? sub[ic].V : sub[ic].VC);
0705                 extra.push_back(hnv);
0706             }
0707             //
0708             //---- physics pages: projections with subgroup errors (ir2=0: R2yy(dy); ir2=1: R2(dy), R2(dphi),
0709             //---- whole and narrow window, uncorrected and corrected) and dC2 = C2C - C2 (ir2=1)
0710             {
0711                 Phys &P     = ph[ipaty];
0712                 TAxis *ax   = tmpl->GetXaxis(), *ay = tmpl->GetYaxis();
0713                 int nx      = ax->GetNbins(), ny = ay->GetNbins();
0714                 //---- mean over the full-stats-valid bins mapped to each output bin; the same inside every chunk
0715                 //---- (a chunk bin that is not usable takes the full-stats value); error = s/sqrt(N_chunks).
0716                 //---- mean=false: the sum over those bins instead (C2 projections, pairs per event)
0717                 auto project    = [&](int kind, bool corr, auto outbin, TH1D* hout, bool mean=true){
0718                     TH2D *vfull = corr ? full.VC : full.V;
0719                     int nb      = hout->GetNbinsX();
0720                     std::vector<double> vf(nb+2,0.), nf(nb+2,0.);
0721                     std::vector<std::vector<double>> vc(nlists,std::vector<double>(nb+2,0.));
0722                     for (int ix=1;ix<=nx;ix++) for (int iy=1;iy<=ny;iy++){
0723                         if (vfull->GetBinContent(ix,iy)<=0.) continue;
0724                         int ob  = outbin(ix,iy);
0725                         if (ob<1 || ob>nb) continue;
0726                         double f    = full.h[kind]->GetBinContent(ix,iy);
0727                         vf[ob] += f; nf[ob] += 1.;
0728                         for (int ic=0;ic<nlists;ic++){
0729                             if (!mean){ vc[ic][ob] += sub[ic].h[kind]->GetBinContent(ix,iy); continue; }    // C2: linear, every chunk counts
0730                             bool u  = ((corr ? sub[ic].VC : sub[ic].V)->GetBinContent(ix,iy)>0.);
0731                             vc[ic][ob]  += u ? sub[ic].h[kind]->GetBinContent(ix,iy) : f;
0732                             ++fillT[ipaty]; if (!u) ++fillN[ipaty];
0733                         }
0734                     }
0735                     for (int ob=1;ob<=nb;ob++){
0736                         if (nf[ob]<=0.) continue;
0737                         double s=0,s2=0;
0738                         double nrm  = mean ? nf[ob] : 1.;
0739                         for (int ic=0;ic<nlists;ic++){ double x = vc[ic][ob]/nrm; s += x; s2 += x*x; }
0740                         double var  = (s2-s*s/nlists)/(nlists-1.);
0741                         hout->SetBinContent(ob, vf[ob]/nrm);
0742                         hout->SetBinError(ob, sqrt(std::max(0.,var))/sqrt((double)nlists));
0743                     }
0744                 };
0745                 //---- README_SplitTracks573 sec 26.1: R2(dy) from (y1,y2) as sum S / sum M - 1 along each dy diagonal (an
0746                 //---- equal-weight mean of R2 cells is pulled up by the acceptance-edge cells); same validity/fill-in as project()
0747                 //---- 2026-09-30: outbins(ix,iy,v) gives the output bins of a cell with their weights (a cell of unequal y widths
0748                 //---- is spread over the dy bins it overlaps; DyOverlap, fluct_common.h)
0749                 auto projectRatio   = [&](auto outbins, TH1D* hout){
0750                     int nb      = hout->GetNbinsX();
0751                     std::vector<double> sf(nb+2,0.), mf(nb+2,0.);
0752                     std::vector<std::vector<double>> sc(nlists,std::vector<double>(nb+2,0.)), mc(nlists,std::vector<double>(nb+2,0.));
0753                     std::vector<std::pair<int,double>> ov;
0754                     for (int ix=1;ix<=nx;ix++) for (int iy=1;iy<=ny;iy++){
0755                         if (full.V->GetBinContent(ix,iy)<=0.) continue;
0756                         outbins(ix,iy,ov);
0757                         double fs   = full.h[0]->GetBinContent(ix,iy), fm = full.h[1]->GetBinContent(ix,iy);
0758                         for (int ic=0;ic<nlists;ic++){ ++fillT[ipaty]; if (sub[ic].V->GetBinContent(ix,iy)<=0.) ++fillN[ipaty]; }
0759                         for (auto &o : ov){
0760                             int ob  = o.first; double w = o.second;
0761                             if (ob<1 || ob>nb) continue;
0762                             sf[ob] += w*fs; mf[ob] += w*fm;
0763                             for (int ic=0;ic<nlists;ic++){
0764                                 bool u  = (sub[ic].V->GetBinContent(ix,iy)>0.);
0765                                 sc[ic][ob]  += w*(u ? sub[ic].h[0]->GetBinContent(ix,iy) : fs);
0766                                 mc[ic][ob]  += w*(u ? sub[ic].h[1]->GetBinContent(ix,iy) : fm);
0767                             }
0768                         }
0769                     }
0770                     for (int ob=1;ob<=nb;ob++){
0771                         if (mf[ob]<=0.) continue;
0772                         double s=0,s2=0; int n=0;
0773                         for (int ic=0;ic<nlists;ic++){ if (mc[ic][ob]<=0.) continue; double x = sc[ic][ob]/mc[ic][ob]-1.; s += x; s2 += x*x; ++n; }
0774                         hout->SetBinContent(ob, sf[ob]/mf[ob]-1.);
0775                         if (n>1) hout->SetBinError(ob, sqrt(std::max(0.,(s2-s*s/n)/(n-1.)))/sqrt((double)n));
0776                     }
0777                 };
0778                 if (ir2==0){
0779                     TH1D *tyy   = (TH1D*)fch[0]->Get(Form("hR2yydy_%d",ipaty));
0780                     if (!tyy){ cout<<"corral::Finalize -- hR2yydy_"<<ipaty<<" missing in chunk 0, exit"<<endl; exit(1); }
0781                     P.R2yy  = (TH1D*)tyy->Clone(Form("hFin_R2yydy_%d",ipaty)); P.R2yy->SetDirectory(0); P.R2yy->Reset(); delete tyy;
0782                     P.R2yy->SetTitle(Form("%s%s, R_{2}(dy) from (y_{1},y_{2}), #SigmaS/#SigmaM (#pm subgroup err);dy;R_{2}",
0783                         ParticleIDNames[PairTypes_Info[ipaty][0]],ParticleIDNames[PairTypes_Info[ipaty][1]]));
0784                     TAxis *ayy  = P.R2yy->GetXaxis();
0785                     double w1 = ax->GetBinWidth(1), w2 = ay->GetBinWidth(1), wd = ayy->GetBinWidth(1);
0786                     bool sameW  = (fabs(w1-w2)<1.e-6 && fabs(w1-wd)<1.e-6);
0787                     projectRatio([&](int ix,int iy,std::vector<std::pair<int,double>>& ov){
0788                         if (sameW){ ov.assign(1,{ayy->FindFixBin(ax->GetBinCenter(ix)-ay->GetBinCenter(iy)),1.0}); return; }
0789                         DyOverlap(ax->GetBinCenter(ix),w1,ay->GetBinCenter(iy),w2,ayy->GetNbins(),ayy->GetXmin(),wd,ov);
0790                     },P.R2yy);
0791                     extraW.push_back(P.R2yy);
0792                 }
0793                 if (ir2==1){
0794                     TString pn  = Form("%s%s",ParticleIDNames[PairTypes_Info[ipaty][0]],ParticleIDNames[PairTypes_Info[ipaty][1]]);
0795                     std::vector<double> cwp(ny+1,0.);
0796                     for (int iy=1;iy<=ny;iy++){ double c = ay->GetBinCenter(iy); while (c>=180.) c -= 360.; while (c<-180.) c += 360.; cwp[iy] = c; }
0797                     double dycut    = std::max(DYCUTPAGE,   0.51*ax->GetBinWidth(1));
0798                     double dphicut  = std::max(DPHICUTPAGE, 0.51*ay->GetBinWidth(1));
0799                     for (int k=0;k<2;k++){
0800                         TString axn = (k==0) ? "dy" : "d#phi";
0801                         TString cut = (k==0) ? TString(Form("|d#phi|<%.0f#circ",dphicut)) : TString(Form("|dy|<%.2g",dycut));
0802                         int nb      = (k==0) ? nx : ny;
0803                         double lo_  = (k==0) ? ax->GetXmin() : ay->GetXmin(), hi_ = (k==0) ? ax->GetXmax() : ay->GetXmax();
0804                         const char* xt  = (k==0) ? "dy" : "d#phi (deg)";
0805                         P.pU[k] = new TH1D(Form("hFin_R2%s_U_%d", k==0?"dy":"dphi",ipaty),Form("%s, R_{2}(%s) from R_{2}(dy,d#phi);%s;R_{2}",pn.Data(),axn.Data(),xt),nb,lo_,hi_);
0806                         P.pC[k] = new TH1D(Form("hFin_R2%s_C_%d", k==0?"dy":"dphi",ipaty),Form("%s, R_{2}(%s) from R_{2}(dy,d#phi);%s;R_{2}",pn.Data(),axn.Data(),xt),nb,lo_,hi_);
0807                         P.nU[k] = new TH1D(Form("hFin_R2%s_nU_%d",k==0?"dy":"dphi",ipaty),Form("%s, R_{2}(%s) from R_{2}(dy,d#phi), %s;%s;R_{2}",pn.Data(),axn.Data(),cut.Data(),xt),nb,lo_,hi_);
0808                         P.nC[k] = new TH1D(Form("hFin_R2%s_nC_%d",k==0?"dy":"dphi",ipaty),Form("%s, R_{2}(%s) from R_{2}(dy,d#phi), %s;%s;R_{2}",pn.Data(),axn.Data(),cut.Data(),xt),nb,lo_,hi_);
0809                         P.qU[k] = new TH1D(Form("hFin_C2%s_U_%d", k==0?"dy":"dphi",ipaty),Form("%s, C_{2}(%s) = #Sigma_{%s} C_{2}(dy,d#phi) (#pm subgroup err);%s;pairs per event",pn.Data(),axn.Data(),k==0?"d#phi":"dy",xt),nb,lo_,hi_);
0810                         P.qC[k] = new TH1D(Form("hFin_C2%s_C_%d", k==0?"dy":"dphi",ipaty),Form("%s, C_{2C}(%s) = #Sigma_{%s} C_{2C}(dy,d#phi), crossing-corrected (#pm subgroup err);%s;pairs per event",pn.Data(),axn.Data(),k==0?"d#phi":"dy",xt),nb,lo_,hi_);
0811                         TH1D *hh[6] = {P.pU[k],P.pC[k],P.nU[k],P.nC[k],P.qU[k],P.qC[k]};
0812                         for (int j=0;j<6;j++){ hh[j]->SetDirectory(0); extraW.push_back(hh[j]); }
0813                     }
0814                     project(2,false,[&](int ix,int iy){ return ix; },P.qU[0],false);
0815                     project(5,true ,[&](int ix,int iy){ return ix; },P.qC[0],false);
0816                     project(2,false,[&](int ix,int iy){ return iy; },P.qU[1],false);
0817                     project(5,true ,[&](int ix,int iy){ return iy; },P.qC[1],false);
0818                     project(3,false,[&](int ix,int iy){ return ix; },P.pU[0]);
0819                     project(6,true ,[&](int ix,int iy){ return ix; },P.pC[0]);
0820                     project(3,false,[&](int ix,int iy){ return iy; },P.pU[1]);
0821                     project(6,true ,[&](int ix,int iy){ return iy; },P.pC[1]);
0822                     project(3,false,[&](int ix,int iy){ return fabs(cwp[iy])<dphicut ? ix : -1; },P.nU[0]);
0823                     project(6,true ,[&](int ix,int iy){ return fabs(cwp[iy])<dphicut ? ix : -1; },P.nC[0]);
0824                     project(3,false,[&](int ix,int iy){ return fabs(ax->GetBinCenter(ix))<dycut ? iy : -1; },P.nU[1]);
0825                     project(6,true ,[&](int ix,int iy){ return fabs(ax->GetBinCenter(ix))<dycut ? iy : -1; },P.nC[1]);
0826                     //---- dC2 = C2C - C2 in bins valid in both; I = its integral = pairs per event added by the
0827                     //---- correction; frac = I / sum rho2(S) over the valid bins. Per chunk: the same bins
0828                     //---- C2 and C2C are linear in a chunk (fixed mode: C2 = S - M everywhere), so no fill-in.
0829                     P.dC2   = (TH2D*)full.h[5]->Clone(Form("hFin_dC2_%d",ipaty)); P.dC2->SetDirectory(0); P.dC2->Reset();
0830                     P.dC2->SetTitle(Form("%s, #DeltaC_{2} = C_{2C} - C_{2} (pairs added by the crossing correction);dy;d#phi (deg)",pn.Data()));
0831                     P.dC2dy     = new TH1D(Form("hFin_dC2dy_%d",ipaty),  Form("%s, #Sigma_{d#phi} #DeltaC_{2} (#pm subgroup err);dy;pairs per event",pn.Data()),nx,ax->GetXmin(),ax->GetXmax());
0832                     P.dC2dphi   = new TH1D(Form("hFin_dC2dphi_%d",ipaty),Form("%s, #Sigma_{dy} #DeltaC_{2} (#pm subgroup err);d#phi (deg);pairs per event",pn.Data()),ny,ay->GetXmin(),ay->GetXmax());
0833                     P.dC2dy->SetDirectory(0); P.dC2dphi->SetDirectory(0);
0834                     P.Ic.assign(nlists,0.); P.fracc.assign(nlists,0.); P.efracc.assign(nlists,0.);
0835                     std::vector<double> Sc(nlists,0.), eI2c(nlists,0.);
0836                     std::vector<std::vector<double>> cdy(nlists,std::vector<double>(nx+2,0.)), cdp(nlists,std::vector<double>(ny+2,0.));
0837                     P.I = 0.; P.Sfull = 0.;
0838                     for (int ix=1;ix<=nx;ix++) for (int iy=1;iy<=ny;iy++){
0839                         if (full.V->GetBinContent(ix,iy)<=0.) continue;
0840                         P.Sfull += full.h[0]->GetBinContent(ix,iy);
0841                         for (int ic=0;ic<nlists;ic++) Sc[ic] += sub[ic].h[0]->GetBinContent(ix,iy);
0842                         if (full.VC->GetBinContent(ix,iy)<=0.) continue;
0843                         double d    = full.h[5]->GetBinContent(ix,iy) - full.h[2]->GetBinContent(ix,iy);
0844                         P.dC2->SetBinContent(ix,iy,d);
0845                         P.I += d;
0846                         P.dC2dy->AddBinContent(ix,d); P.dC2dphi->AddBinContent(iy,d);
0847                         for (int ic=0;ic<nlists;ic++){
0848                             double dc   = sub[ic].h[5]->GetBinContent(ix,iy) - sub[ic].h[2]->GetBinContent(ix,iy);  // linear: no fill-in
0849                             P.Ic[ic] += dc; cdy[ic][ix] += dc; cdp[ic][iy] += dc;
0850                             eI2c[ic] += pow(sub[ic].h[5]->GetBinError(ix,iy),2);    // dirty bins: C2C error = the clean bin's
0851                         }
0852                     }
0853                     P.frac  = (P.Sfull>0.) ? P.I/P.Sfull : 0.;
0854                     double s=0,s2=0,sf=0,sf2=0;
0855                     for (int ic=0;ic<nlists;ic++){
0856                         double fr   = (Sc[ic]>0.) ? P.Ic[ic]/Sc[ic] : 0.;
0857                         P.fracc[ic] = fr; P.efracc[ic] = (Sc[ic]>0.) ? sqrt(eI2c[ic])/Sc[ic] : 0.;
0858                         s += P.Ic[ic]; s2 += P.Ic[ic]*P.Ic[ic]; sf += fr; sf2 += fr*fr;
0859                     }
0860                     P.eI    = sqrt(std::max(0.,(s2-s*s/nlists)/(nlists-1.)))/sqrt((double)nlists);
0861                     P.efrac = sqrt(std::max(0.,(sf2-sf*sf/nlists)/(nlists-1.)))/sqrt((double)nlists);
0862                     for (int ix=1;ix<=nx;ix++){ double a=0,a2=0; for (int ic=0;ic<nlists;ic++){ a+=cdy[ic][ix]; a2+=cdy[ic][ix]*cdy[ic][ix]; }
0863                         P.dC2dy->SetBinError(ix, sqrt(std::max(0.,(a2-a*a/nlists)/(nlists-1.)))/sqrt((double)nlists)); }
0864                     for (int iy=1;iy<=ny;iy++){ double a=0,a2=0; for (int ic=0;ic<nlists;ic++){ a+=cdp[ic][iy]; a2+=cdp[ic][iy]*cdp[ic][iy]; }
0865                         P.dC2dphi->SetBinError(iy, sqrt(std::max(0.,(a2-a*a/nlists)/(nlists-1.)))/sqrt((double)nlists)); }
0866                     extraW.push_back(P.dC2); extraW.push_back(P.dC2dy); extraW.push_back(P.dC2dphi);
0867                     cout<<Form("corral::Finalize -- dC2 pairtype %2d: integral %+.4e +- %.1e pairs/event (x %.0f events = %.4g pairs), sum rho2(S) %.4e, fraction %+.3e +- %.1e",
0868                         ipaty,P.I,P.eI,nevtAll,P.I*nevtAll,P.Sfull,P.frac,P.efrac)<<endl;
0869                 }
0870             }
0871             //---- step 1b: crossing loss / dphi-mirror asymmetry (dy,dphi only), full stats + every chunk
0872             if (ir2==1){
0873                 XLoss &X    = xl[ipaty];
0874                 X.have      = true; X.ds = ds; X.half = (ds!=0) ? ds : -1;
0875                 TAxis *ax   = tmpl->GetXaxis(), *ay = tmpl->GetYaxis();
0876                 int nx      = ax->GetNbins(), ny = ay->GetNbins();
0877                 double lo   = ay->GetXmin();
0878                 std::vector<int> jm(ny+1,0); std::vector<double> cw(ny+1,0.);
0879                 for (int iy=1;iy<=ny;iy++){     // dphi mirror bin, same rule as CrossingCorrect.h
0880                     double c    = ay->GetBinCenter(iy);
0881                     while (c>= 180.) c -= 360.;
0882                     while (c< -180.) c += 360.;
0883                     cw[iy]      = c;
0884                     double cm   = -c;
0885                     while (cm<  lo      ) cm += 360.;
0886                     while (cm>= lo+360. ) cm -= 360.;
0887                     jm[iy]      = ay->FindFixBin(cm);
0888                 }
0889                 double dycut    = std::max(XLDYCUT, 0.51*ax->GetBinWidth(1));
0890                 X.dycut     = dycut;
0891                 auto inW    = [&](int w,int ix,int iy)->bool{
0892                     double c    = cw[iy];
0893                     if (c*X.half<=0.) return false;         // tested half only
0894                     bool near   = fabs(c)<90.;
0895                     if (w==2) return !near;
0896                     if (!near) return false;
0897                     return (w==0) || fabs(ax->GetBinCenter(ix))<dycut;
0898                 };
0899                 auto okF    = [&](int ix,int iy){ return full.V->GetBinContent(ix,iy)>0. && full.V->GetBinContent(ix,jm[iy])>0.; };
0900                 //---- D from the Zvtx-averaged S and M: C2 = S - M in every valid zvtx bin and the average is linear
0901                 //---- with the same masks and weights, so avg(C2) = avg(S) - avg(M) at full stats (checked below),
0902                 //---- and a chunk's S and M are defined wherever the full-stats masks are (no chunk skipped)
0903                 //---- DA = acceptance-corrected deficit: sum[S(mirror)*M(bin)/M(mirror) - S(bin)], M ratio from FULL
0904                 //---- stats (precise, and keeps DA linear in a chunk's S) = sum M(bin)[R2(mirror)-R2(bin)] at full stats
0905                 auto sums   = [&](const CFSet& cs, int w, double& D, double& eD2, double& Sc, double& DM, double& Mc, double& DA){
0906                     D=0; eD2=0; Sc=0; DM=0; Mc=0; DA=0;
0907                     for (int ix=1;ix<=nx;ix++) for (int iy=1;iy<=ny;iy++){
0908                         if (!inW(w,ix,iy) || !okF(ix,iy)) continue;
0909                         int jy  = jm[iy];
0910                         double sj=cs.h[0]->GetBinContent(ix,jy), si=cs.h[0]->GetBinContent(ix,iy);
0911                         double mj=cs.h[1]->GetBinContent(ix,jy), mi=cs.h[1]->GetBinContent(ix,iy);
0912                         D   += (sj-mj) - (si-mi);
0913                         eD2 += pow(cs.h[0]->GetBinError(ix,jy),2) + pow(cs.h[1]->GetBinError(ix,jy),2) + pow(cs.h[0]->GetBinError(ix,iy),2) + pow(cs.h[1]->GetBinError(ix,iy),2);
0914                         Sc  += sj;
0915                         DM  += mj - mi;
0916                         Mc  += mj;
0917                         double mfj=full.h[1]->GetBinContent(ix,jy), mfi=full.h[1]->GetBinContent(ix,iy);
0918                         if (mfj>0.) DA += sj*mfi/mfj - si;
0919                     }
0920                 };
0921                 //---- check: avg(C2) == avg(S) - avg(M) at full stats over all valid bins
0922                 {
0923                     double dmx=0;
0924                     for (int ix=1;ix<=nx;ix++) for (int iy=1;iy<=ny;iy++){
0925                         if (full.V->GetBinContent(ix,iy)<=0.) continue;
0926                         dmx = std::max(dmx, fabs(full.h[2]->GetBinContent(ix,iy) - (full.h[0]->GetBinContent(ix,iy)-full.h[1]->GetBinContent(ix,iy))));
0927                     }
0928                     cout<<Form("corral::Finalize -- step 1b check pairtype %2d: max |avg(C2) - (avg(S)-avg(M))| = %.2e over valid bins",ipaty,dmx)<<endl;
0929                 }
0930                 for (int w=0;w<3;w++){
0931                     double D,eD2,Sc,DM,Mc,DA;
0932                     sums(full,w,D,eD2,Sc,DM,Mc,DA);
0933                     X.f[w]  = (Sc>0.) ? D/Sc : 0.;
0934                     X.fM[w] = (Mc>0.) ? DM/Mc : 0.;
0935                     X.fA[w] = (Sc>0.) ? DA/Sc : 0.;
0936                     X.fc[w].assign(nlists,0.); X.efc[w].assign(nlists,0.); X.fMc[w].assign(nlists,0.); X.fAc[w].assign(nlists,0.); X.usec[w].assign(nlists,false);
0937                     double s=0,s2=0,sm=0,sm2=0,sa=0,sa2=0; int n=0;
0938                     for (int ic=0;ic<nlists;ic++){
0939                         sums(sub[ic],w,D,eD2,Sc,DM,Mc,DA);
0940                         if (Sc<=0. || Mc<=0.) continue;
0941                         X.usec[w][ic]   = true;
0942                         X.fc[w][ic]     = D/Sc;  X.efc[w][ic] = sqrt(eD2)/Sc;  X.fMc[w][ic] = DM/Mc;  X.fAc[w][ic] = DA/Sc;
0943                         s += D/Sc; s2 += pow(D/Sc,2); sm += DM/Mc; sm2 += pow(DM/Mc,2); sa += DA/Sc; sa2 += pow(DA/Sc,2); ++n;
0944                     }
0945                     X.nuse[w]   = n;
0946                     if (n>=NVALIDMIN){
0947                         X.ef[w]     = sqrt(std::max(0.,(s2-s*s/n)/(n-1.)))/sqrt((double)n);
0948                         X.efM[w]    = sqrt(std::max(0.,(sm2-sm*sm/n)/(n-1.)))/sqrt((double)n);
0949                         X.efA[w]    = sqrt(std::max(0.,(sa2-sa*sa/n)/(n-1.)))/sqrt((double)n);
0950                     }
0951                 }
0952                 //---- check: on the near half, D = sum(C2C - C2) over the dirty bins (the correction copies the clean C2)
0953                 {
0954                     double DC=0, Sc=0;
0955                     for (int ix=1;ix<=nx;ix++) for (int iy=1;iy<=ny;iy++){
0956                         if (!inW(0,ix,iy) || !okF(ix,iy)) continue;
0957                         DC  += full.h[5]->GetBinContent(ix,iy) - full.h[2]->GetBinContent(ix,iy);
0958                         Sc  += full.h[0]->GetBinContent(ix,jm[iy]);
0959                     }
0960                     X.fC    = (Sc>0.) ? DC/Sc : 0.;
0961                 }
0962                 //---- maps: A = C2(mirror)-C2 (antisymmetric), and D per dy column over the near half
0963                 X.A     = (TH2D*)full.h[2]->Clone(Form("hC2mirrorAsym_1_%d",ipaty)); X.A->SetDirectory(0); X.A->Reset();
0964                 for (int ix=1;ix<=nx;ix++) for (int iy=1;iy<=ny;iy++){
0965                     if (!okF(ix,iy)) continue;
0966                     int jy  = jm[iy];
0967                     X.A->SetBinContent(ix,iy, full.h[2]->GetBinContent(ix,jy) - full.h[2]->GetBinContent(ix,iy));
0968                     X.A->SetBinError(  ix,iy, sqrt(pow(full.h[2]->GetBinError(ix,jy),2) + pow(full.h[2]->GetBinError(ix,iy),2)));
0969                 }
0970                 X.Ady   = new TH1D(Form("hC2mirrorAsymDy_1_%d",ipaty),"",nx,ax->GetXmin(),ax->GetXmax()); X.Ady->SetDirectory(0);
0971                 for (int ix=1;ix<=nx;ix++){
0972                     auto colD   = [&](const CFSet& cs, bool& ok){
0973                         double D=0; ok = true; int nb=0;
0974                         for (int iy=1;iy<=ny;iy++){
0975                             if (!inW(0,ix,iy) || !okF(ix,iy)) continue;
0976                             D   += (cs.h[0]->GetBinContent(ix,jm[iy]) - cs.h[1]->GetBinContent(ix,jm[iy])) - (cs.h[0]->GetBinContent(ix,iy) - cs.h[1]->GetBinContent(ix,iy)); ++nb;
0977                         }
0978                         if (nb==0) ok = false;
0979                         return D;
0980                     };
0981                     bool ok; double Dfull = colD(full,ok);
0982                     if (!ok) continue;
0983                     double s=0,s2=0; int n=0;
0984                     for (int ic=0;ic<nlists;ic++){ bool okc; double d = colD(sub[ic],okc); if (!okc) continue; s+=d; s2+=d*d; ++n; }
0985                     X.Ady->SetBinContent(ix,Dfull);
0986                     X.Ady->SetBinError(ix, (n>=NVALIDMIN) ? sqrt(std::max(0.,(s2-s*s/n)/(n-1.)))/sqrt((double)n) : 0.);
0987                 }
0988                 //---- occupancy proxy per chunk: <n1><n2> from the chunk's hmult
0989                 X.x.assign(nlists,0.);
0990                 for (int ic=0;ic<nlists;ic++){
0991                     TH2D *hm    = (TH2D*)fch[ic]->Get(Form("hmult_%d",ipaty));
0992                     if (hm){ X.x[ic] = hm->GetMean(1)*hm->GetMean(2); delete hm; }
0993                 }
0994                 { double xs=0; for (double v : X.x) xs += v; xs /= nlists; X.xlo = *std::min_element(X.x.begin(),X.x.end())/xs; X.xhi = *std::max_element(X.x.begin(),X.x.end())/xs; }
0995                 cout<<Form("corral::Finalize -- step 1b pairtype %2d acceptance-corrected fA: near %+.3e +- %.1e | near |dy| %+.3e +- %.1e | away %+.3e +- %.1e ; <n1><n2>/mean over chunks %.3f..%.3f",
0996                     ipaty,X.fA[0],X.efA[0],X.fA[1],X.efA[1],X.fA[2],X.efA[2],X.xlo,X.xhi)<<endl;
0997                 cout<<Form("corral::Finalize -- step 1b pairtype %2d (ds %+d, tested half dphi%s0): f near %+.3e +- %.1e (%d chunks) | near |dy|<%.2g %+.3e +- %.1e | away %+.3e +- %.1e | fM near %+.3e +- %.1e | check sum(C2C-C2)/S %+.3e",
0998                     ipaty,ds,X.half<0?"<":">",X.f[0],X.ef[0],X.nuse[0],dycut,X.f[1],X.ef[1],X.f[2],X.ef[2],X.fM[0],X.efM[0],X.fC)<<endl;
0999             }
1000             for (int k=0;k<NK;k++){
1001                 TString nm  = Form("%s_%d",Form(KINDNAME[k],ir2),ipaty);
1002                 //---- the map was built from a per-zvtx slice clone: take title and axis titles from the chunk's own map
1003                 TH1 *t      = (TH1*)fch[0]->Get(nm.Data());
1004                 if (t){
1005                     full.h[k]->SetTitle(t->GetTitle());
1006                     full.h[k]->GetXaxis()->SetTitle(t->GetXaxis()->GetTitle());
1007                     full.h[k]->GetYaxis()->SetTitle(t->GetYaxis()->GetTitle());
1008                     delete t;
1009                 }
1010                 result[nm]  = full.h[k]; full.h[k] = 0;
1011             }
1012             full.Delete();
1013             for (int ic=0;ic<nlists;ic++) sub[ic].Delete();     // README_Finalize step 1b will use sub[] here
1014             for (int ic=0;ic<nlists;ic++) for (int z=0;z<NZ;z++){ delete cS[ic][z]; delete cM[ic][z]; }
1015             delete tmpl; delete mRef;
1016         }
1017         cout<<"corral::Finalize -- pairtype "<<ipaty<<" done (dirty side "<<ds<<")"<<endl;
1018     }
1019     //
1020     //---- compare with the reference, write, draw
1021     TString OutputFileBase  = (OutputName!="") ? OutputName : TString(Form("corral_%s",RunString.Data()));
1022     TString RootFileName    = TString("./root/") + OutputFileBase + ".root";
1023     TString PdfFileName     = TString("./pdf/")  + OutputFileBase + ".pdf";
1024     TFile *fout = new TFile(RootFileName.Data(),"RECREATE");
1025     gStyle->SetOptStat(0);
1026     TCanvas *cfin   = new TCanvas("cfin","cfin",1400,1000);
1027     cfin->Print((PdfFileName+"[").Data());
1028     //
1029     //==== FRONT PAGES (2026-10-12): title page (how corral ran, from info/ of the chunks; contents), QA pages of the summed
1030     //==== chunk histograms with the cuts drawn, vs-segment pages (segCells). FinalizeQA.h.
1031     FinQA::Info finInfo = FinQA::ReadChunkInfo(fch);
1032     std::vector<TH1*> qaHists;
1033     TTree *segQA    = 0;
1034     {
1035         std::vector<TString> pairNames;
1036         for (int k=0;k<NPairTypes;k++) pairNames.push_back(Form("%s%s",ParticleIDNames[PairTypes_Info[k][0]],ParticleIDNames[PairTypes_Info[k][1]]));
1037         //---- runs: from the chunks' segCells (all runs, in order), else from lists_summary.txt
1038         std::vector<int> runsT;
1039         for (TFile *f : fch){ TTree *t = (TTree*)f->Get("segCells"); if (!t) continue; Int_t r; t->SetBranchStatus("*",0); t->SetBranchStatus("run",1); t->SetBranchAddress("run",&r);
1040             for (Long64_t e=0;e<t->GetEntries();e++){ t->GetEntry(e); if (std::find(runsT.begin(),runsT.end(),r)==runsT.end()) runsT.push_back(r); } delete t; }
1041         if (runsT.empty()) for (int r : chunkRun) if (r>=0 && std::find(runsT.begin(),runsT.end(),r)==runsT.end()) runsT.push_back(r);
1042         //---- events in the physics range: events reaching CalcRm with |zvtx| < fz (hnvtxCR)
1043         double nevtPhys = -1.;
1044         { TH2 *h = (TH2*)FinQA::Sum(fch,"hnvtxCR"); if (h){ nevtPhys = h->Integral(h->GetXaxis()->FindBin(-valFinZMax+1e-6),h->GetXaxis()->FindBin(valFinZMax-1e-6),0,-1); delete h; } }
1045         FinQA::TitlePage(cfin,SET,finInfo,RunString,valFinZMax,Form("git %s, built %s",CORRAL_GIT,CORRAL_BUILD_TIME),runsT,nevtPhys,pairNames,fref!=nullptr);
1046         cfin->Print(PdfFileName.Data());
1047         FinQA::QAPages(cfin,PdfFileName.Data(),fch,finInfo,valFinZMax,qaHists);
1048         segQA   = FinQA::SegmentPages(cfin,PdfFileName.Data(),fch);
1049         cfin->Clear(); cfin->cd();
1050         for (const auto &m : finInfo.mismatch) cout<<"corral::Finalize -- info/: "<<m<<endl;
1051     }
1052     //
1053     //==== PHYSICS PAGES FIRST (user 2026-09-26), per pairtype: 1 overview (as corral_m's first page, but the
1054     //==== (dy,dphi) row crossing-corrected), 2 C(Q) (ratio of zvtx sums), 3 crossing page (3x2, as corral_m),
1055     //==== 4 dC2 = C2C - C2 page. Diagnostic pages follow after the summary.
1056     {
1057         TCanvas *cph    = new TCanvas("cph","cph",1200,800);
1058         //---- V0 mass histograms (global, counts): plain sums over chunks
1059         std::map<TString,TH1*> v0m;
1060         for (const char* n : {"hKSmass","hLAmass","hLAmass_clean","hALmass","hALmass_clean"}){
1061             for (int ic=0;ic<nlists;ic++){
1062                 TH1 *h  = (TH1*)fch[ic]->Get(n);
1063                 if (!h) continue;
1064                 if (!v0m[n]){ v0m[n] = (TH1*)h->Clone(Form("%s_fin",n)); v0m[n]->SetDirectory(0); }
1065                 else v0m[n]->Add(h);
1066                 delete h;
1067             }
1068         }
1069         auto pn = [&](int ipaty){ return TString(Form("%s%s",ParticleIDNames[PairTypes_Info[ipaty][0]],ParticleIDNames[PairTypes_Info[ipaty][1]])); };
1070         auto clone2 = [&](const char* nm, const char* tit){
1071             TH2D *h = (TH2D*)result[nm];
1072             if (!h) return (TH2D*)0;
1073             TH2D *c = (TH2D*)h->Clone(Form("%s_page",nm)); c->SetDirectory(0); c->SetTitle(tit);
1074             return c;
1075         };
1076         auto zeroline   = [&](){ gPad->Update(); TLine *l0 = new TLine(gPad->GetUxmin(),0.,gPad->GetUxmax(),0.); l0->SetLineColor(kGray+1); l0->SetLineStyle(2); l0->Draw(); };
1077         auto runlines   = [&](double a, double b){ for (int k=1;k<nlists;k++){ if (chunkRun[k]==chunkRun[k-1]) continue; TLine *lr = new TLine(k-0.5,a,k-0.5,b); lr->SetLineStyle(2); lr->SetLineColor(kGray+1); lr->Draw(); } };
1078         std::vector<int> runsP;
1079         for (int ic=0;ic<nlists;ic++) if (std::find(runsP.begin(),runsP.end(),chunkRun[ic])==runsP.end()) runsP.push_back(chunkRun[ic]);
1080         const int runcolP[6]    = {kBlack,kRed+1,kBlue+1,kGreen+2,kMagenta+1,kOrange+1};
1081         std::vector<double> slopeX(NPairTypes,0.), eslopeX(NPairTypes,0.);
1082         for (int ipaty=0;ipaty<NPairTypes;ipaty++){
1083             Phys &P = ph[ipaty];
1084             if (!P.pU[0] || !P.R2yy) continue;
1085             TString nm  = pn(ipaty);
1086             int pid1    = PairTypes_Info[ipaty][0], pid2 = PairTypes_Info[ipaty][1];
1087             //---------------- page 1: overview
1088             cph->Clear(); cph->Divide(4,3,0.0001,0.0001);
1089             TH1 *hm1    = result[Form("hmult1_%d",ipaty)], *hm2 = result[Form("hmult2_%d",ipaty)], *hmm = result[Form("hmult_%d",ipaty)];
1090             cph->cd(1); gPad->SetLogy(1);
1091             if (hm1){
1092                 TH1 *d  = (TH1*)hm1->Clone(Form("hm1page_%d",ipaty)); d->SetDirectory(0);
1093                 d->SetTitle(Form("%s, n_{1} per event (Finalize, all chunks);n_{1}",nm.Data()));
1094                 d->Draw();
1095                 TLatex tx; tx.SetNDC(); tx.SetTextFont(42); tx.SetTextSize(0.06); tx.SetTextAlign(31);
1096                 tx.DrawLatex(0.88,0.85,Form("N_{evt}=%.0f",hmm?hmm->GetEntries():0.));
1097                 tx.DrawLatex(0.88,0.78,Form("<N_{1}>=%.3g",hm1->GetMean()));
1098                 if (P.hy1) tx.DrawLatex(0.88,0.71,Form("INT #rho_{1}(1)=%.3g",P.hy1->Integral()));
1099                 if (hm2) tx.DrawLatex(0.88,0.60,Form("<N_{2}>=%.3g",hm2->GetMean()));
1100                 if (P.hy2) tx.DrawLatex(0.88,0.53,Form("INT #rho_{1}(2)=%.3g",P.hy2->Integral()));
1101             }
1102             cph->cd(2); P.R2yy->SetMarkerStyle(20); P.R2yy->SetMarkerSize(0.5); P.R2yy->Draw("E1");
1103             cph->cd(3);
1104             { TH1D *d = (TH1D*)P.pC[1]->Clone(Form("pC1page_%d",ipaty)); d->SetDirectory(0);
1105               d->SetTitle(Form("%s, R_{2C}(d#phi) crossing-corrected (#pm subgroup err);d#phi (deg);R_{2}",nm.Data()));
1106               d->SetMarkerStyle(20); d->SetMarkerSize(0.5); d->Draw("E1"); }
1107             //---- V0 mass pads (as corral_m: one per distinct V0 species of the pair)
1108             {
1109                 std::vector<int> v0s;
1110                 for (int pid : {pid1,pid2}){
1111                     if (std::find(v0s.begin(),v0s.end(),pid)!=v0s.end()) continue;
1112                     if (pid==kParticleIDKshort || pid==kParticleIDLambda || pid==kParticleIDAntiLambda) v0s.push_back(pid);
1113                 }
1114                 for (size_t iv=0;iv<v0s.size() && iv<2;iv++){
1115                     cph->cd(iv==0 ? 4 : 8);
1116                     const char* b   = (v0s[iv]==kParticleIDKshort) ? "hKSmass" : (v0s[iv]==kParticleIDLambda) ? "hLAmass" : "hALmass";
1117                     if (v0m[b]) v0m[b]->Draw();
1118                     TString cl  = TString(b)+"_clean";
1119                     if (v0s[iv]!=kParticleIDKshort && v0m[cl]) v0m[cl]->Draw("same");
1120                 }
1121             }
1122             //---- row 2: (y1,y2), not crossing-correctable
1123             {
1124                 TH2D *a = clone2(Form("hrho2_0_%d",ipaty),    Form("%s, #rho_{2}(y_{1},y_{2});y_{1};y_{2}",nm.Data()));
1125                 TH2D *b = clone2(Form("hrho1rho1_0_%d",ipaty),Form("%s, #rho_{1}#rho_{1}(y_{1},y_{2}) (mixed);y_{1};y_{2}",nm.Data()));
1126                 TH2D *c = clone2(Form("hR2_0_%d",ipaty),      Form("%s, R_{2}(y_{1},y_{2});y_{1};y_{2}",nm.Data()));
1127                 cph->cd(5); if (a) a->Draw("surf3");
1128                 cph->cd(6); if (b) b->Draw("surf3");
1129                 cph->cd(7); if (c) c->Draw("lego2");
1130             }
1131             //---- row 3: (dy,dphi), crossing-corrected
1132             {
1133                 TH2D *a = clone2(Form("hrho2C_1_%d",ipaty),   Form("%s, #rho_{2C}(dy,d#phi) crossing-corrected;dy;d#phi",nm.Data()));
1134                 TH2D *b = clone2(Form("hrho1rho1_1_%d",ipaty),Form("%s, #rho_{1}#rho_{1}(dy,d#phi) (mixed);dy;d#phi",nm.Data()));
1135                 TH2D *c = clone2(Form("hR2C_1_%d",ipaty),     Form("%s, R_{2C}(dy,d#phi) crossing-corrected;dy;d#phi",nm.Data()));
1136                 TH2D *e = clone2(Form("hR2C_1_%d",ipaty),     Form("%s, R_{2C}(dy,d#phi) crossing-corrected;dy;d#phi (deg)",nm.Data()));
1137                 cph->cd(9);  if (a) a->Draw("surf3");
1138                 cph->cd(10); if (b) b->Draw("surf3");
1139                 cph->cd(11); if (c) c->Draw("lego2");
1140                 cph->cd(12); if (e){ gPad->SetRightMargin(0.14); e->Draw("colz"); }
1141             }
1142             cph->cd(); cph->Update(); EtaLabelsIf(cph,ipaty); cph->Print(PdfFileName.Data());
1143             //---------------- page 2: C(Q) (2x2, as corral_m): C(Q) 0-2.0 GeV (5 MeV) with signposts | C(Q) in the
1144             //---------------- 4 kT bins (10 MeV) | C(Q) 0-0.3 GeV (1 MeV) | Minv near threshold, sibling vs mixed
1145             if (P.CQ && P.CQ5 && P.CQKT10){
1146                 cph->Clear(); cph->Divide(2,2,0.0001,0.0001);
1147                 auto yrange = [&](TH1D* h, double qlo, double qhi, double top){
1148                     double ymn=1e30, ymx=-1e30;
1149                     for (int ib=h->FindBin(qlo);ib<=h->FindBin(qhi);ib++){
1150                         if (h->GetBinError(ib)<=0.) continue;
1151                         ymn = std::min(ymn,h->GetBinContent(ib)); ymx = std::max(ymx,h->GetBinContent(ib));
1152                     }
1153                     if (ymx>ymn){ h->SetMinimum(ymn-0.05*(ymx-ymn)); h->SetMaximum(ymx+top*(ymx-ymn)); }
1154                 };
1155                 cph->cd(1);
1156                 { TH1D *h = (TH1D*)P.CQ5->Clone(Form("CQ5page_%d",ipaty)); h->SetDirectory(0);
1157                   yrange(h,0.0201,2.0,0.25); h->Draw(); gPad->Update();
1158                   DrawCQSignposts(pid1,pid2,gPad->GetUymin(),gPad->GetUymax(),h->GetXaxis()->GetXmax()); }
1159                 cph->cd(2);
1160                 {
1161                     const int kcol[4]   = {kBlue,kGreen+2,kOrange+1,kRed};
1162                     TH1D *hk[4]; double ymn=1e30, ymx=-1e30;
1163                     for (int ik=0;ik<4;ik++){
1164                         hk[ik]  = P.CQKT10->ProjectionX(Form("CQKTpage_%d_%d",ipaty,ik),ik+1,ik+1,"e"); hk[ik]->SetDirectory(0);
1165                         hk[ik]->SetLineColor(kcol[ik]); hk[ik]->SetMarkerColor(kcol[ik]);
1166                         for (int ib=hk[ik]->FindBin(0.0201);ib<=hk[ik]->GetNbinsX();ib++){
1167                             if (hk[ik]->GetBinError(ib)<=0.) continue;
1168                             ymn = std::min(ymn,hk[ik]->GetBinContent(ib)); ymx = std::max(ymx,hk[ik]->GetBinContent(ib));
1169                         }
1170                     }
1171                     hk[0]->SetTitle(Form("%s, C(Q) in k_{T} bins (10 MeV bins, #pm subgroup err);Q_{inv} (GeV)",nm.Data()));
1172                     if (ymx>ymn){ hk[0]->SetMinimum(ymn-0.05*(ymx-ymn)); hk[0]->SetMaximum(ymx+0.05*(ymx-ymn)); }
1173                     hk[0]->Draw("hist");
1174                     for (int ik=1;ik<4;ik++) hk[ik]->Draw("hist same");
1175                     TLegend *lgk    = new TLegend(0.55,0.65,0.89,0.89);
1176                     lgk->SetBorderSize(0); lgk->SetFillStyle(0); lgk->SetTextSize(0.035);
1177                     for (int ik=0;ik<4;ik++){
1178                         double k1   = P.CQKT10->GetYaxis()->GetBinLowEdge(ik+1), k2 = P.CQKT10->GetYaxis()->GetBinUpEdge(ik+1);
1179                         TLegendEntry *le    = lgk->AddEntry(hk[ik],Form("%.2f<k_{T}<%.2f GeV",k1,k2),"l");
1180                         le->SetLineColor(kcol[ik]); le->SetLineWidth(2); le->SetTextColor(kcol[ik]);
1181                     }
1182                     lgk->Draw();
1183                 }
1184                 cph->cd(3);
1185                 { TH1D *h = (TH1D*)P.CQ->Clone(Form("CQzoompage_%d",ipaty)); h->SetDirectory(0);
1186                   h->GetXaxis()->SetRangeUser(0.,0.3);
1187                   h->SetTitle(Form("%s, C(Q) (1 MeV bins, zoom, #pm subgroup err);Q_{inv} (GeV);C(Q)",nm.Data()));
1188                   yrange(h,0.0201,0.2999,0.05); h->Draw(); }
1189                 cph->cd(4);
1190                 if (P.MinvS && P.MinvM){
1191                     TH1D *ms    = (TH1D*)P.MinvS->Clone(Form("MinvSpage_%d",ipaty)); ms->SetDirectory(0);
1192                     TH1D *mm    = (TH1D*)P.MinvM->Clone(Form("MinvMpage_%d",ipaty)); mm->SetDirectory(0);
1193                     int b1      = ms->GetXaxis()->FindBin(ms->GetXaxis()->GetXmax()-0.150+1e-6), b2 = ms->GetNbinsX();
1194                     double im   = mm->Integral(b1,b2);
1195                     double fmix = (im>0.) ? ms->Integral(b1,b2)/im : 0.;
1196                     if (fmix>0.) mm->Scale(fmix);
1197                     ms->SetLineColor(kGreen+2); mm->SetLineColor(kBlue);
1198                     ms->SetMinimum(0.); ms->SetMaximum(1.1*std::max(ms->GetMaximum(),mm->GetMaximum()));
1199                     ms->Draw("hist"); mm->Draw("hist same");
1200                     TLegend *lgm    = new TLegend(0.30,0.11,0.89,0.21);
1201                     lgm->SetBorderSize(0); lgm->SetFillColor(kWhite); lgm->SetTextSize(0.030);
1202                     TLegendEntry *l1    = lgm->AddEntry(ms,"sibling (raw counts)","l"); l1->SetLineColor(kGreen+2);
1203                     TLegendEntry *l2    = lgm->AddEntry(mm,Form("mixed #times %.4g (sib/mix, top 150 MeV)",fmix),"l"); l2->SetLineColor(kBlue);
1204                     lgm->Draw();
1205                 }
1206                 cph->cd(); cph->Update(); EtaLabelsIf(cph,ipaty); cph->Print(PdfFileName.Data());
1207             }
1208             //---------------- page 2b: ir2=2 (Q_inv, pT_pair), 2x2, the same page as in corral_m (QPTPage.h); frame 1 with subgroup errors
1209             if (P.QPTsib && P.QPTmix){
1210                 DrawQPTPage(cph,P.QPTsib,P.QPTmix,P.CQvar,pid1,pid2,nm,"#Sigma_{z} sib / #Sigma_{z} mix",Form("fin_%d",ipaty));
1211                 EtaLabelsIf(cph,ipaty); cph->Print(PdfFileName.Data());
1212             }
1213             //---------------- page 3: crossing page (3x2), maps from Finalize, projections with subgroup errors
1214             {
1215                 TH2D *hu    = clone2(Form("hR2_1_%d",ipaty), Form("%s, R_{2}(dy,d#phi) uncorrected;dy;d#phi (deg)",nm.Data()));
1216                 TH2D *hc    = clone2(Form("hR2C_1_%d",ipaty),Form("%s, R_{2}(dy,d#phi) crossing-corrected;dy;d#phi (deg)",nm.Data()));
1217                 TH2D *hmM   = (TH2D*)result[Form("hrho1rho1_1_%d",ipaty)];
1218                 int ds      = xl[ipaty].ds;
1219                 int nx=hu->GetNbinsX(), ny=hu->GetNbinsY(); double lo=hu->GetYaxis()->GetXmin();
1220                 std::vector<int> jmir(ny+1,0); std::vector<double> cwr(ny+1,0.);
1221                 for (int iy=1;iy<=ny;iy++){
1222                     double c    = hu->GetYaxis()->GetBinCenter(iy);
1223                     while (c>= 180.) c -= 360.; while (c< -180.) c += 360.;
1224                     cwr[iy]     = c;
1225                     double cm   = -c; while (cm<lo) cm += 360.; while (cm>=lo+360.) cm -= 360.;
1226                     jmir[iy]    = hu->GetYaxis()->FindFixBin(cm);
1227                 }
1228                 auto okU    = [&](int ix,int iy){ return hmM->GetBinContent(ix,iy)>0.; };
1229                 auto okC    = [&](int ix,int iy){
1230                     if (ds==0) return okU(ix,iy);
1231                     if (fabs(cwr[iy])<90.) return okU(ix, (cwr[iy]*ds>0.) ? jmir[iy] : iy);
1232                     return okU(ix,iy) || okU(ix,jmir[iy]);
1233                 };
1234                 double zmin=1e30, zmax=-1e30;
1235                 for (int ix=1;ix<=nx;ix++) for (int iy=1;iy<=ny;iy++){
1236                     if (okU(ix,iy)){ zmin = std::min(zmin,hu->GetBinContent(ix,iy)); zmax = std::max(zmax,hu->GetBinContent(ix,iy)); }
1237                     if (okC(ix,iy)){ zmin = std::min(zmin,hc->GetBinContent(ix,iy)); zmax = std::max(zmax,hc->GetBinContent(ix,iy)); }
1238                 }
1239                 if (zmin<zmax){ hu->SetMinimum(zmin); hu->SetMaximum(zmax); hc->SetMinimum(zmin); hc->SetMaximum(zmax); }
1240                 std::vector<std::vector<bool>> flag(nx+2,std::vector<bool>(ny+2,false));
1241                 int nflag   = 0;
1242                 if (ds!=0){
1243                     for (int ix=1;ix<=nx;ix++) for (int iy=1;iy<=ny;iy++){
1244                         if (!(fabs(cwr[iy])<90. && cwr[iy]*ds>0.)) continue;
1245                         int jy  = jmir[iy];
1246                         if (!okU(ix,iy) || !okU(ix,jy)) continue;
1247                         double d    = hu->GetBinContent(ix,iy) - hu->GetBinContent(ix,jy);
1248                         double e    = sqrt(pow(hu->GetBinError(ix,iy),2) + pow(hu->GetBinError(ix,jy),2));
1249                         if (e>0. && fabs(d)>3.*e){ flag[ix][iy] = true; ++nflag; }
1250                     }
1251                 }
1252                 cph->Clear(); cph->Divide(3,2,0.0001,0.0001);
1253                 cph->cd(1); gPad->SetRightMargin(0.14); hu->Draw("colz");
1254                 for (int ix=1;ix<=nx;ix++) for (int iy=1;iy<=ny;iy++){
1255                     if (!flag[ix][iy]) continue;
1256                     double x1=hu->GetXaxis()->GetBinLowEdge(ix), x2=hu->GetXaxis()->GetBinUpEdge(ix);
1257                     double y1=hu->GetYaxis()->GetBinLowEdge(iy), y2=hu->GetYaxis()->GetBinUpEdge(iy);
1258                     TLine *ln[4]    = {0,0,0,0};
1259                     if (!flag[ix-1][iy]) ln[0] = new TLine(x1,y1,x1,y2);
1260                     if (!flag[ix+1][iy]) ln[1] = new TLine(x2,y1,x2,y2);
1261                     if (!flag[ix][iy-1]) ln[2] = new TLine(x1,y1,x2,y1);
1262                     if (!flag[ix][iy+1]) ln[3] = new TLine(x1,y2,x2,y2);
1263                     for (int k=0;k<4;k++) if (ln[k]){ ln[k]->SetLineColor(kRed); ln[k]->SetLineWidth(2); ln[k]->Draw(); }
1264                 }
1265                 { TLatex tx; tx.SetNDC(); tx.SetTextFont(42); tx.SetTextSize(0.035); tx.SetTextColor(kRed); tx.SetTextAlign(11);
1266                   tx.DrawLatex(0.02,0.015,Form("red: dirty |R_{2}-mirror|>3#sigma_{sg}, %d bins (d#phi%s0)",nflag,ds>0?">":"<")); }
1267                 cph->cd(4); gPad->SetRightMargin(0.14); hc->Draw("colz");
1268                 TH1D *hyy   = (TH1D*)P.R2yy->Clone(Form("R2yypage_%d",ipaty)); hyy->SetDirectory(0);
1269                 hyy->SetLineColor(kGreen+2); hyy->SetMarkerColor(kGreen+2); hyy->SetMarkerStyle(21); hyy->SetMarkerSize(0.5);
1270                 for (int k=0;k<2;k++){
1271                     TH1D *hs[4] = {P.pU[k],P.pC[k],P.nU[k],P.nC[k]};
1272                     const int col[4]    = {kBlack,kRed,kBlue,kMagenta};
1273                     const int mst[4]    = {20,24,20,24};
1274                     for (int j=0;j<4;j++){ hs[j]->SetLineColor(col[j]); hs[j]->SetMarkerColor(col[j]); hs[j]->SetMarkerStyle(mst[j]); hs[j]->SetMarkerSize(0.5); }
1275                     for (int iw=0;iw<2;iw++){
1276                         cph->cd(2+3*k+iw);
1277                         TH1D *ha    = (iw==0) ? P.pU[k] : P.nU[k];
1278                         TH1D *hb    = (iw==0) ? P.pC[k] : P.nC[k];
1279                         bool doyy   = (iw==0 && k==0);
1280                         double pmin = std::min(ha->GetMinimum(), hb->GetMinimum());
1281                         double pmax = std::max(ha->GetMaximum(), hb->GetMaximum());
1282                         if (doyy){ pmin = std::min(pmin,hyy->GetMinimum()); pmax = std::max(pmax,hyy->GetMaximum()); }
1283                         TH1D *hf    = (TH1D*)ha->Clone(Form("%s_frame",ha->GetName())); hf->SetDirectory(0);
1284                         hf->SetMinimum(pmin-0.1*(pmax-pmin)); hf->SetMaximum(pmax+0.30*(pmax-pmin));
1285                         hf->Draw("E1");
1286                         if (k==1){ TLine *l0 = new TLine(0.,hf->GetMinimum(),0.,hf->GetMaximum()); l0->SetLineColor(kGray); l0->Draw(); hf->Draw("E1 same"); }
1287                         hb->Draw("E1 same");
1288                         if (doyy) hyy->Draw("E1 same");
1289                         TLegend *lgc    = new TLegend(0.25,doyy?0.74:0.78,0.75,0.89);
1290                         lgc->SetBorderSize(1); lgc->SetFillStyle(0); lgc->SetTextSize(0.04);
1291                         lgc->AddEntry(ha,"No Corr","lp");
1292                         lgc->AddEntry(hb,"Cross Corr","lp");
1293                         if (doyy) lgc->AddEntry(hyy,"from R_{2}(y_{1},y_{2}), No Corr","lp");
1294                         lgc->Draw();
1295                     }
1296                 }
1297                 cph->cd(); cph->Update(); EtaLabelsIf(cph,ipaty); cph->Print(PdfFileName.Data());
1298             }
1299             //---------------- page 4: dC2 = C2C - C2
1300             {
1301                 XLoss &X    = xl[ipaty];
1302                 cph->Clear(); cph->Divide(3,2,0.0001,0.0001);
1303                 cph->cd(1); gPad->SetRightMargin(0.14);
1304                 P.dC2->Draw("colz");
1305                 if (X.ds==0){ TLatex tx; tx.SetNDC(); tx.SetTextFont(42); tx.SetTextSize(0.05); tx.DrawLatex(0.15,0.5,"no crossing correction (dirty side 0): C_{2C} = C_{2}"); }
1306                 cph->cd(2); P.dC2dphi->SetMarkerStyle(20); P.dC2dphi->SetMarkerSize(0.5); P.dC2dphi->Draw("E0"); zeroline();    // E0: the clean side is identically 0 (content and error) and must show
1307                 {   // thin gray verticals at dphi = 0 and 90 deg: the clean side between them is identically 0
1308                     gPad->Update();
1309                     for (double xl : {0.,90.}){ TLine *lv = new TLine(xl,gPad->GetUymin(),xl,gPad->GetUymax()); lv->SetLineColor(kGray); lv->SetLineWidth(1); lv->Draw(); }
1310                     P.dC2dphi->Draw("E0 same");
1311                 }
1312                 cph->cd(3); P.dC2dy->SetMarkerStyle(20);   P.dC2dy->SetMarkerSize(0.5);   P.dC2dy->Draw("E0");   zeroline();
1313                 //---- per chunk: fraction of pairs added vs chunk index, and vs occupancy
1314                 double ymn=1e30, ymx=-1e30;
1315                 TGraphErrors *gi    = new TGraphErrors(); gi->SetName(Form("gdC2chunk_%d",ipaty));
1316                 for (int ic=0;ic<nlists;ic++){
1317                     gi->SetPoint(gi->GetN(),ic,P.fracc[ic]); gi->SetPointError(gi->GetN()-1,0.,P.efracc[ic]);
1318                     ymn = std::min(ymn,P.fracc[ic]-P.efracc[ic]); ymx = std::max(ymx,P.fracc[ic]+P.efracc[ic]);
1319                 }
1320                 if (!(ymx>ymn)){ ymn=-1e-3; ymx=1e-3; }
1321                 { double dd=ymx-ymn; ymn -= 0.1*dd; ymx += 0.3*dd; }
1322                 double chi2=0, ndf=0, slc=0, eslc=0, slx=0, eslx=0;
1323                 cph->cd(4);
1324                 gPad->DrawFrame(-1.,ymn,nlists,ymx,Form("%s, pairs added / #Sigma#rho_{2}(S), per chunk;chunk index (time order);fraction added",nm.Data()));
1325                 runlines(ymn,ymx);
1326                 { TBox *bx = new TBox(-1.,P.frac-P.efrac,nlists,P.frac+P.efrac); bx->SetFillColorAlpha(kGray,0.5); bx->SetLineWidth(0); bx->Draw(); }
1327                 gi->SetMarkerStyle(20); gi->SetMarkerSize(0.6);
1328                 if (X.ds!=0 && gi->GetN()>2){
1329                     TFitResultPtr r0    = gi->Fit("pol0","QNS"); if (r0.Get() && r0->IsValid()){ chi2 = r0->Chi2(); ndf = r0->Ndf(); }
1330                     TFitResultPtr r1    = gi->Fit("pol1","QNS"); if (r1.Get() && r1->IsValid()){ slc = r1->Parameter(1); eslc = r1->ParError(1); }
1331                 }
1332                 gi->Draw("P");
1333                 extraW.push_back(gi);
1334                 cph->cd(5);
1335                 {
1336                     double xs=0; for (double v : X.x) xs += v; double xm = (nlists>0 && xs>0.) ? xs/nlists : 1.;
1337                     double xmn=1e30, xmx=-1e30;
1338                     TGraphErrors *ga    = new TGraphErrors(); ga->SetName(Form("gdC2occ_%d",ipaty));
1339                     std::vector<TGraphErrors*> gr(runsP.size(),(TGraphErrors*)0);
1340                     for (int ic=0;ic<nlists;ic++){
1341                         double xx   = X.x[ic]/xm;
1342                         ga->SetPoint(ga->GetN(),xx,P.fracc[ic]); ga->SetPointError(ga->GetN()-1,0.,P.efracc[ic]);
1343                         int ir  = std::find(runsP.begin(),runsP.end(),chunkRun[ic]) - runsP.begin();
1344                         if (!gr[ir]){ gr[ir] = new TGraphErrors(); gr[ir]->SetMarkerStyle(20); gr[ir]->SetMarkerSize(0.7); gr[ir]->SetMarkerColor(runcolP[ir%6]); gr[ir]->SetLineColor(runcolP[ir%6]); }
1345                         gr[ir]->SetPoint(gr[ir]->GetN(),xx,P.fracc[ic]); gr[ir]->SetPointError(gr[ir]->GetN()-1,0.,P.efracc[ic]);
1346                         xmn = std::min(xmn,xx); xmx = std::max(xmx,xx);
1347                     }
1348                     if (xmx>xmn){ double dx=xmx-xmn; xmn -= 0.05*dx; xmx += 0.05*dx; } else { xmn -= 0.05; xmx += 0.05; }
1349                     gPad->DrawFrame(xmn,ymn,xmx,ymx,Form("%s, fraction added per chunk vs occupancy;#LTn_{1}#GT#LTn_{2}#GT / mean;fraction added",nm.Data()));
1350                     if (X.ds!=0 && ga->GetN()>2){
1351                         TF1 *fl = new TF1(Form("fdC2occ_%d",ipaty),"pol1",xmn,xmx);
1352                         TFitResultPtr r = ga->Fit(fl,"QNS");
1353                         if (r.Get() && r->IsValid()){ slx = r->Parameter(1); eslx = r->ParError(1); }
1354                         fl->SetLineColor(kGray+2); fl->SetLineWidth(1); fl->Draw("same");
1355                     }
1356                     TLegend *lg = new TLegend(0.14,0.78,0.89,0.89);
1357                     lg->SetBorderSize(0); lg->SetFillStyle(0); lg->SetTextSize(0.035); lg->SetNColumns(2);
1358                     for (size_t ir=0;ir<runsP.size();ir++) if (gr[ir]){ gr[ir]->Draw("P"); lg->AddEntry(gr[ir],Form("run %d",runsP[ir]),"p"); }
1359                     lg->Draw();
1360                     extraW.push_back(ga);
1361                 }
1362                 slopeX[ipaty] = slx; eslopeX[ipaty] = eslx;
1363                 cph->cd(6);
1364                 {
1365                     TLatex tx; tx.SetNDC(); tx.SetTextFont(42); tx.SetTextSize(0.05); tx.SetTextAlign(11);
1366                     double y    = 0.88;
1367                     tx.DrawLatex(0.06,y,Form("%s, #DeltaC_{2} = C_{2C} - C_{2}, full statistics",nm.Data())); y -= 0.09;
1368                     tx.DrawLatex(0.06,y,Form("integral: %+.4e #pm %.1e pairs/event",P.I,P.eI)); y -= 0.08;
1369                     tx.DrawLatex(0.06,y,Form("#times N_{evt} = %.4g pairs added",P.I*(hmm?hmm->GetEntries():0.))); y -= 0.08;
1370                     tx.DrawLatex(0.06,y,Form("#Sigma#rho_{2}(S) = %.4e pairs/event",P.Sfull)); y -= 0.08;
1371                     tx.DrawLatex(0.06,y,Form("fraction added: %+.3e #pm %.1e",P.frac,P.efrac)); y -= 0.08;
1372                     if (X.ds!=0){
1373                         tx.DrawLatex(0.06,y,Form("const fit over chunks: #chi^{2}/ndf = %.1f/%.0f",chi2,ndf)); y -= 0.08;
1374                         tx.DrawLatex(0.06,y,Form("slope vs chunk: %+.2e #pm %.1e per chunk",slc,eslc)); y -= 0.08;
1375                         tx.DrawLatex(0.06,y,Form("slope vs occupancy: %+.2e #pm %.1e",slx,eslx)); y -= 0.08;
1376                     } else {
1377                         tx.DrawLatex(0.06,y,"not pt-ordered (neutral leg): not corrected"); y -= 0.08;
1378                     }
1379                     tx.SetTextSize(0.035);
1380                     tx.DrawLatex(0.06,0.06,Form("errors: subgroup (%d chunks); per-chunk bars: propagated (approximate)",nlists));
1381                 }
1382                 //cph->cd(); cph->Update(); cph->Print(PdfFileName.Data()); // page 4 (dC2) not printed (2026-10-10, not needed); still drawn, the fit line below stays in the log
1383                 cout<<Form("corral::Finalize -- dC2 fits pairtype %2d %-18s: const chi2/ndf %.1f/%.0f, slope vs chunk %+.2e +- %.1e, slope vs occupancy %+.2e +- %.1e; projection fill-ins %.2e",
1384                     ipaty,nm.Data(),chi2,ndf,slc,eslc,slx,eslx,fillT[ipaty]>0?(double)fillN[ipaty]/fillT[ipaty]:0.)<<endl;
1385             }
1386         }
1387         //---------------- dC2 summary page: fraction of pairs added per pairtype; per-chunk overlay for the pions
1388         {
1389             cph->Clear(); cph->Divide(2,1,0.0001,0.0001);
1390             cph->cd(1); gPad->SetBottomMargin(0.18);
1391             TH1D *hf    = new TH1D("hFin_dC2frac","fraction of pairs added by the crossing correction (full stats #pm subgroup err);;#Sigma#DeltaC_{2} / #Sigma#rho_{2}(S)",NPairTypes,-0.5,NPairTypes-0.5);
1392             hf->SetDirectory(0);
1393             for (int ipaty=0;ipaty<NPairTypes;ipaty++){ hf->SetBinContent(ipaty+1,ph[ipaty].frac); hf->SetBinError(ipaty+1,ph[ipaty].efrac); hf->GetXaxis()->SetBinLabel(ipaty+1,pn(ipaty).Data()); }
1394             hf->SetMarkerStyle(20); hf->SetMarkerSize(0.9); hf->GetXaxis()->SetLabelSize(0.04);
1395             hf->Draw("E1 X0"); zeroline();
1396             extraW.push_back(hf);
1397             cph->cd(2);
1398             double a=1e30, b=-1e30;
1399             for (int ipaty=0;ipaty<4;ipaty++) for (int ic=0;ic<nlists && ic<(int)ph[ipaty].fracc.size();ic++){ a = std::min(a,ph[ipaty].fracc[ic]-ph[ipaty].efracc[ic]); b = std::max(b,ph[ipaty].fracc[ic]+ph[ipaty].efracc[ic]); }
1400             if (!(b>a)){ a=-1e-3; b=1e-3; }
1401             { double dd=b-a; a -= 0.05*dd; b += 0.35*dd; }
1402             gPad->DrawFrame(-1.,a,nlists,b,"fraction of pairs added per chunk, pions;chunk index (time order);fraction added");
1403             runlines(a,b);
1404             TLegend *lg = new TLegend(0.14,0.78,0.89,0.89); lg->SetBorderSize(0); lg->SetFillStyle(0); lg->SetTextSize(0.035); lg->SetNColumns(4);
1405             const int pc[4] = {kBlack,kRed+1,kBlue+1,kGreen+2};
1406             for (int ipaty=0;ipaty<4;ipaty++){
1407                 if (ph[ipaty].fracc.empty()) continue;
1408                 TGraphErrors *g = new TGraphErrors();
1409                 for (int ic=0;ic<nlists;ic++){ g->SetPoint(g->GetN(),ic+0.15*ipaty,ph[ipaty].fracc[ic]); g->SetPointError(g->GetN()-1,0.,ph[ipaty].efracc[ic]); }
1410                 g->SetMarkerStyle(20); g->SetMarkerSize(0.5); g->SetMarkerColor(pc[ipaty]); g->SetLineColor(pc[ipaty]); g->Draw("P");
1411                 lg->AddEntry(g,pn(ipaty).Data(),"p");
1412             }
1413             lg->Draw();
1414             cph->cd(); cph->Update(); cph->Print(PdfFileName.Data());
1415         }
1416         cfin->cd();
1417     }
1418     //==== DIAGNOSTIC PAGES from here on
1419     for (int ih=0;ih<(fref?NFINALIZEHISTS:0);ih++){
1420         TString base    = TString(FinalizeHists[ih]);
1421         int kd=-1,ir;
1422         bool isCF   = ParseCF(base,kd,ir);
1423         for (int ipage=0;ipage<(isCF?2:1);ipage++){     // page 0: vs reference; page 1 (CF): Finalize error / reference error
1424             cfin->Clear(); cfin->Divide(4,4,0.0001,0.0001);
1425             for (int ipaty=0;ipaty<NPairTypes;ipaty++){
1426                 TString name    = Form("%s_%d",base.Data(),ipaty);
1427                 TH1 *hf     = result[name];
1428                 TH1 *href   = (TH1*)fref->Get(name.Data());
1429                 if (!hf || !href){ if (ipage==0) cout<<"corral::Finalize -- "<<name<<" missing (result or reference), skipped"<<endl; continue; }
1430                 cfin->cd(1+ipaty);
1431                 if (ipage==1){
1432                     //---- subgroup error / reference (single-job, per-run-style) error, bins where both are >0
1433                     TH2D *hr    = (TH2D*)hf->Clone(Form("%s_errratio",name.Data())); hr->SetDirectory(0); hr->Reset();
1434                     std::vector<double> rr;
1435                     for (int ib=0;ib<hf->GetNcells();ib++){
1436                         if (hf->IsBinUnderflow(ib) || hf->IsBinOverflow(ib)) continue;  // u/o keep the chain's per-run errors
1437                         double ef=hf->GetBinError(ib), er=href->GetBinError(ib);
1438                         if (ef>0. && er>0.){ hr->SetBinContent(ib,ef/er); rr.push_back(ef/er); }
1439                     }
1440                     cout<<Form("corral::Finalize -- %-17s subgroup err / ref err: median %.3f (%d bins)",name.Data(),Median(rr),(int)rr.size())<<endl;
1441                     hr->SetTitle(Form("%s: #sigma_{Finalize}/#sigma_{ref} (median %.2f)",name.Data(),Median(rr)));
1442                     hr->SetMinimum(0.); hr->SetMaximum(3.);
1443                     hr->Draw("colz");
1444                     continue;
1445                 }
1446                 //---- relative difference (in-range bins with ref!=0) and pull vs the reference error (bins with both errors>0)
1447                 double sd=0,sd2=0,dmax=0,sp=0,sp2=0; int nb=0,np=0;
1448                 for (int ib=0;ib<hf->GetNcells();ib++){
1449                     if (hf->IsBinUnderflow(ib) || hf->IsBinOverflow(ib)) continue;
1450                     double a=hf->GetBinContent(ib), r=href->GetBinContent(ib);
1451                     if (r!=0.){ double d=(a-r)/fabs(r); sd+=d; sd2+=d*d; dmax=std::max(dmax,fabs(d)); ++nb; }
1452                     double er=href->GetBinError(ib);
1453                     if (er>0. && hf->GetBinError(ib)>0.){ double p=(a-r)/er; sp+=p; sp2+=p*p; ++np; }
1454                 }
1455                 double mrel = nb?sd/nb:0., rrel = nb?sqrt(std::max(0.,sd2/nb-mrel*mrel)):0.;
1456                 double mp   = np?sp/np:0., rp = np?sqrt(std::max(0.,sp2/np-mp*mp)):0.;
1457                 cout<<Form("corral::Finalize -- %-17s %s  rel diff mean %+.2e rms %.2e max %.2e (%d bins)   pull vs ref err: mean %+.3f rms %.3f (%d bins)",
1458                     name.Data(), isCF?"[CF] ":"[sum]", mrel,rrel,dmax,nb, mp,rp,np)<<endl;
1459                 if (hf->GetDimension()==1){
1460                     href->SetLineColor(kBlack); hf->SetLineColor(kRed);
1461                     href->Draw("hist"); hf->Draw("hist same");
1462                 } else {
1463                     TH2D *hp    = (TH2D*)hf->Clone(Form("%s_pull",name.Data())); hp->SetDirectory(0); hp->Reset();
1464                     for (int ib=0;ib<hf->GetNcells();ib++){
1465                         double er=href->GetBinError(ib);
1466                         if (er>0.) hp->SetBinContent(ib,(hf->GetBinContent(ib)-href->GetBinContent(ib))/er);
1467                     }
1468                     hp->SetTitle(Form("%s: (Finalize - ref)/#sigma_{ref}",name.Data()));
1469                     hp->SetMinimum(-5.); hp->SetMaximum(5.);
1470                     hp->Draw("colz");
1471                 }
1472             }
1473             cfin->Print(PdfFileName.Data());
1474         }
1475     }
1476     //---- step 1b pages. Per pairtype (2x2): A = C2(mirror)-C2 map | D(dy) over the near tested half |
1477     //---- f per chunk vs chunk index (run boundaries dashed; near half black, |dy| window blue; full-stats
1478     //---- value +- subgroup error as a band) | f per chunk vs <n1><n2>/mean, colored by run, linear fit.
1479     //---- Per-chunk error bars: propagated per-run C2 errors (approximate); pol0 chi2 says whether the
1480     //---- chunk-to-chunk variation exceeds them.
1481     auto pname  = [&](int ipaty){ return TString(Form("%s%s",ParticleIDNames[PairTypes_Info[ipaty][0]],ParticleIDNames[PairTypes_Info[ipaty][1]])); };
1482     std::vector<int> runs;
1483     for (int ic=0;ic<nlists;ic++) if (std::find(runs.begin(),runs.end(),chunkRun[ic])==runs.end()) runs.push_back(chunkRun[ic]);
1484     const int runcol[6] = {kBlack,kRed+1,kBlue+1,kGreen+2,kMagenta+1,kOrange+1};
1485     for (int ipaty=0;ipaty<NPairTypes;ipaty++){
1486         XLoss &X    = xl[ipaty];
1487         if (!X.have) continue;
1488         cfin->Clear(); cfin->Divide(2,2,0.0001,0.0001);
1489         //---- pad 1: A map
1490         cfin->cd(1); gPad->SetRightMargin(0.14);
1491         double amax = 0.;
1492         for (int ib=0;ib<X.A->GetNcells();ib++) amax = std::max(amax,fabs(X.A->GetBinContent(ib)));
1493         X.A->SetTitle(Form("%s, C_{2}(dy,-d#phi) - C_{2}(dy,d#phi) (full stats);dy;d#phi (deg)",pname(ipaty).Data()));
1494         if (amax>0.){ X.A->SetMinimum(-amax); X.A->SetMaximum(amax); }
1495         X.A->Draw("colz");
1496         //---- pad 2: D(dy)
1497         cfin->cd(2);
1498         X.Ady->SetTitle(Form("%s, D(dy) = #Sigma_{near, d#phi%s0}[C_{2}(mirror)-C_{2}], #pm subgroup err;dy;pairs per event",pname(ipaty).Data(),X.half<0?"<":">"));
1499         X.Ady->SetMarkerStyle(20); X.Ady->SetMarkerSize(0.6);
1500         X.Ady->Draw("E1");
1501         gPad->Update();
1502         { TLine *l0 = new TLine(gPad->GetUxmin(),0.,gPad->GetUxmax(),0.); l0->SetLineColor(kGray+1); l0->SetLineStyle(2); l0->Draw(); }
1503         //---- pad 3: f vs chunk index
1504         cfin->cd(3);
1505         TGraphErrors *gc[2];
1506         double ymn=1e30, ymx=-1e30;
1507         for (int w=0;w<2;w++){
1508             gc[w]   = new TGraphErrors();
1509             gc[w]->SetName(Form("gXlossChunk_w%d_%d",w,ipaty));
1510             for (int ic=0;ic<nlists;ic++){
1511                 if (!X.usec[w][ic]) continue;
1512                 int ip  = gc[w]->GetN();
1513                 gc[w]->SetPoint(ip, ic, X.fAc[w][ic]); gc[w]->SetPointError(ip, 0., X.efc[w][ic]);
1514                 ymn = std::min(ymn, X.fAc[w][ic]-X.efc[w][ic]); ymx = std::max(ymx, X.fAc[w][ic]+X.efc[w][ic]);
1515             }
1516             gc[w]->SetMarkerStyle(20); gc[w]->SetMarkerSize(0.6);
1517             gc[w]->SetMarkerColor(w==0?kBlack:kBlue); gc[w]->SetLineColor(w==0?kBlack:kBlue);
1518             extraW.push_back(gc[w]);
1519         }
1520         if (gc[0]->GetN()>1){
1521             gc[0]->Fit("pol0","Q");
1522             TF1 *f0 = gc[0]->GetFunction("pol0");
1523             if (f0){ f0->SetLineColor(kGray+2); f0->SetLineWidth(1); X.chi2c = f0->GetChisquare(); X.ndfc = f0->GetNDF(); }
1524             TFitResultPtr r1    = gc[0]->Fit("pol1","QNS"); // N: not stored, pol0 stays the drawn function
1525             if (r1.Get() && r1->IsValid()){ X.slopec = r1->Parameter(1); X.eslopec = r1->ParError(1); }
1526         }
1527         if (ymx>ymn){ double dd=ymx-ymn; ymn -= 0.1*dd; ymx += 0.35*dd; }
1528         TH1F *fr3   = gPad->DrawFrame(-1., ymn, nlists, ymx, Form("%s, f_{A} (acceptance-corrected mirror deficit / #Sigma#rho_{2}(S)_{mirror}) per chunk;chunk index (time order);f_{A}",pname(ipaty).Data()));
1529         (void)fr3;
1530         for (int k=1;k<nlists;k++){
1531             if (chunkRun[k]==chunkRun[k-1]) continue;
1532             TLine *lr   = new TLine(k-0.5,ymn,k-0.5,ymx); lr->SetLineStyle(2); lr->SetLineColor(kGray+1); lr->Draw();
1533         }
1534         {   // full-stats value +- subgroup error, near half
1535             TBox *bx    = new TBox(-1., X.fA[0]-X.efA[0], nlists, X.fA[0]+X.efA[0]);
1536             bx->SetFillColorAlpha(kGray,0.5); bx->SetLineWidth(0); bx->Draw();
1537             TLine *lz   = new TLine(-1.,0.,nlists,0.); lz->SetLineColor(kGray+1); lz->Draw();
1538         }
1539         gc[0]->Draw("P"); gc[1]->Draw("P");
1540         {
1541             TLegend *lg = new TLegend(0.14,0.72,0.89,0.89);
1542             lg->SetBorderSize(0); lg->SetFillStyle(0); lg->SetTextSize(0.033);
1543             lg->AddEntry(gc[0],Form("near half: full %+.2e #pm %.1e, #chi^{2}/ndf(const) %.1f/%.0f, slope %+.1e #pm %.1e /chunk",X.fA[0],X.efA[0],X.chi2c,X.ndfc,X.slopec,X.eslopec),"p");
1544             lg->AddEntry(gc[1],Form("near |dy|<%.2g: full %+.2e #pm %.1e",X.dycut,X.fA[1],X.efA[1]),"p");
1545             lg->Draw();
1546         }
1547         //---- pad 4: f vs <n1><n2>/mean, by run
1548         cfin->cd(4);
1549         {
1550             double xs=0; int nxs=0;
1551             for (int ic=0;ic<nlists;ic++) if (X.usec[0][ic]){ xs += X.x[ic]; ++nxs; }
1552             double xm   = (nxs>0) ? xs/nxs : 1.;
1553             TGraphErrors *ga    = new TGraphErrors(); ga->SetName(Form("gXlossOcc_%d",ipaty));
1554             std::vector<TGraphErrors*> gr(runs.size(),(TGraphErrors*)0);
1555             double xmn=1e30, xmx=-1e30;
1556             for (int ic=0;ic<nlists;ic++){
1557                 if (!X.usec[0][ic] || xm<=0.) continue;
1558                 double xx   = X.x[ic]/xm;
1559                 ga->SetPoint(ga->GetN(), xx, X.fAc[0][ic]); ga->SetPointError(ga->GetN()-1, 0., X.efc[0][ic]);
1560                 int ir  = std::find(runs.begin(),runs.end(),chunkRun[ic]) - runs.begin();
1561                 if (!gr[ir]){ gr[ir] = new TGraphErrors(); gr[ir]->SetMarkerStyle(20); gr[ir]->SetMarkerSize(0.7); gr[ir]->SetMarkerColor(runcol[ir%6]); gr[ir]->SetLineColor(runcol[ir%6]); }
1562                 gr[ir]->SetPoint(gr[ir]->GetN(), xx, X.fAc[0][ic]); gr[ir]->SetPointError(gr[ir]->GetN()-1, 0., X.efc[0][ic]);
1563                 xmn = std::min(xmn,xx); xmx = std::max(xmx,xx);
1564             }
1565             extraW.push_back(ga);
1566             if (ga->GetN()>1){
1567                 TFitResultPtr r = ga->Fit("pol1","QNS");
1568                 if (r.Get() && r->IsValid()){ X.slopex = r->Parameter(1); X.eslopex = r->ParError(1); }
1569             }
1570             if (xmx>xmn){ double dx=xmx-xmn; xmn -= 0.05*dx; xmx += 0.05*dx; } else { xmn -= 0.05; xmx += 0.05; }
1571             gPad->DrawFrame(xmn, ymn, xmx, ymx, Form("%s, f_{A} per chunk vs occupancy (near half);#LTn_{1}#GT#LTn_{2}#GT / mean;f_{A}",pname(ipaty).Data()));
1572             TLine *lz   = new TLine(xmn,0.,xmx,0.); lz->SetLineColor(kGray+1); lz->Draw();
1573             if (ga->GetN()>1 && X.eslopex>0.){
1574                 TF1 *fl = new TF1(Form("fXlossOcc_%d",ipaty),"pol1",xmn,xmx);
1575                 TFitResultPtr r = ga->Fit(fl,"QNS");
1576                 fl->SetLineColor(kGray+2); fl->SetLineWidth(1); fl->Draw("same");
1577                 (void)r;
1578             }
1579             TLegend *lg = new TLegend(0.14,0.72,0.89,0.89);
1580             lg->SetBorderSize(0); lg->SetFillStyle(0); lg->SetTextSize(0.033); lg->SetNColumns(2);
1581             for (size_t ir=0;ir<runs.size();ir++) if (gr[ir]){ gr[ir]->Draw("P"); lg->AddEntry(gr[ir],Form("run %d",runs[ir]),"p"); }
1582             lg->AddEntry((TObject*)0,Form("slope df/d(x/#bar{x}) %+.2e #pm %.1e",X.slopex,X.eslopex),"");
1583             lg->Draw();
1584         }
1585         EtaLabelsIf(cfin,ipaty); cfin->Print(PdfFileName.Data());
1586         extraW.push_back(X.A); extraW.push_back(X.Ady);
1587     }
1588     //---- step 1b summary pages: (1) full-stats f per pairtype, 3 windows + the rho2(M) asymmetry;
1589     //---- (2) f per chunk vs chunk index, overlaid per pairtype group, and the occupancy slope per pairtype
1590     {
1591         cfin->Clear(); cfin->Divide(1,1);
1592         cfin->cd(1); gPad->SetBottomMargin(0.18);
1593         const int NSM   = 5;
1594         TH1D *hs[NSM];
1595         const char* lab[NSM]    = {"f_{A}, near half","f_{A}, near |dy| window","f_{A}, away half (null)","f (C_{2}, not acc.-corrected), near half","f_{M} = same with #rho_{2}(M), near half"};
1596         const int col[NSM]  = {kBlack,kBlue,kGreen+2,kGray+1,kOrange+1};
1597         double ymn=1e30, ymx=-1e30;
1598         for (int k=0;k<NSM;k++){
1599             hs[k]   = new TH1D(Form("hXlossSummary_%d",k),"",NPairTypes,-0.5,NPairTypes-0.5); hs[k]->SetDirectory(0);
1600             for (int ipaty=0;ipaty<NPairTypes;ipaty++){
1601                 XLoss &X    = xl[ipaty];
1602                 if (!X.have) continue;
1603                 double v    = (k<3) ? X.fA[k] : (k==3) ? X.f[0]  : X.fM[0];
1604                 double e    = (k<3) ? X.efA[k] : (k==3) ? X.ef[0] : X.efM[0];
1605                 hs[k]->SetBinContent(ipaty+1,v); hs[k]->SetBinError(ipaty+1,e);
1606                 hs[k]->GetXaxis()->SetBinLabel(ipaty+1,pname(ipaty).Data());
1607                 ymn = std::min(ymn,v-e); ymx = std::max(ymx,v+e);
1608             }
1609             hs[k]->SetMarkerStyle(20+k); hs[k]->SetMarkerSize(0.9); hs[k]->SetMarkerColor(col[k]); hs[k]->SetLineColor(col[k]);
1610             extraW.push_back(hs[k]);
1611         }
1612         double dd   = ymx-ymn; if (dd<=0.) dd = 1.;
1613         hs[0]->SetMinimum(ymn-0.1*dd); hs[0]->SetMaximum(ymx+0.3*dd);
1614         hs[0]->SetTitle(Form("mirror deficit / #Sigma#rho_{2}(S)_{mirror}, tested half d#phi%s0, full stats #pm subgroup error;;fraction",xl[0].half<0?"<":">"));
1615         hs[0]->GetXaxis()->SetLabelSize(0.035);
1616         hs[0]->Draw("E1 X0");
1617         for (int k=1;k<NSM;k++) hs[k]->Draw("E1 X0 same");
1618         gPad->Update();
1619         { TLine *l0 = new TLine(gPad->GetUxmin(),0.,gPad->GetUxmax(),0.); l0->SetLineColor(kGray+1); l0->SetLineStyle(2); l0->Draw(); }
1620         TLegend *lg = new TLegend(0.55,0.70,0.89,0.89);
1621         lg->SetBorderSize(0); lg->SetFillStyle(0); lg->SetTextSize(0.028);
1622         for (int k=0;k<NSM;k++) lg->AddEntry(hs[k],lab[k],"p");
1623         lg->Draw();
1624         cfin->Print(PdfFileName.Data());
1625         //
1626         cfin->Clear(); cfin->Divide(2,2,0.0001,0.0001);
1627         const char* grp[3]  = {"pions (0-3)","V0-V0 (4-9)","V0-#pi (10-15)"};
1628         const int glo[3]    = {0,4,10}, ghi[3] = {3,9,15};
1629         const int pcol[6]   = {kBlack,kRed+1,kBlue+1,kGreen+2,kMagenta+1,kOrange+1};
1630         for (int g=0;g<3;g++){
1631             cfin->cd(1+g);
1632             double a=1e30, b=-1e30;
1633             for (int ipaty=glo[g];ipaty<=ghi[g];ipaty++){ XLoss &X = xl[ipaty]; if (!X.have) continue;
1634                 for (int ic=0;ic<nlists;ic++) if (X.usec[0][ic]){ a = std::min(a,X.fAc[0][ic]-X.efc[0][ic]); b = std::max(b,X.fAc[0][ic]+X.efc[0][ic]); } }
1635             if (!(b>a)){ a=-1; b=1; }
1636             double d2   = b-a; a -= 0.05*d2; b += 0.35*d2;
1637             gPad->DrawFrame(-1.,a,nlists,b,Form("f_{A} per chunk (near half), %s;chunk index (time order);f_{A}",grp[g]));
1638             for (int k=1;k<nlists;k++){ if (chunkRun[k]==chunkRun[k-1]) continue; TLine *lr = new TLine(k-0.5,a,k-0.5,b); lr->SetLineStyle(2); lr->SetLineColor(kGray+1); lr->Draw(); }
1639             TLine *lz   = new TLine(-1.,0.,nlists,0.); lz->SetLineColor(kGray+1); lz->Draw();
1640             TLegend *lgg    = new TLegend(0.14,0.74,0.89,0.89);
1641             lgg->SetBorderSize(0); lgg->SetFillStyle(0); lgg->SetTextSize(0.035); lgg->SetNColumns(3);
1642             for (int ipaty=glo[g];ipaty<=ghi[g];ipaty++){
1643                 XLoss &X    = xl[ipaty]; if (!X.have) continue;
1644                 TGraphErrors *gg    = new TGraphErrors();
1645                 for (int ic=0;ic<nlists;ic++) if (X.usec[0][ic]){ gg->SetPoint(gg->GetN(), ic+0.12*(ipaty-glo[g]), X.fAc[0][ic]); gg->SetPointError(gg->GetN()-1,0.,X.efc[0][ic]); }
1646                 int c   = pcol[(ipaty-glo[g])%6];
1647                 gg->SetMarkerStyle(20); gg->SetMarkerSize(0.5); gg->SetMarkerColor(c); gg->SetLineColor(c);
1648                 gg->Draw("P");
1649                 lgg->AddEntry(gg,pname(ipaty).Data(),"p");
1650             }
1651             lgg->Draw();
1652         }
1653         cfin->cd(4); gPad->SetBottomMargin(0.18);
1654         TH1D *hsl   = new TH1D("hXlossSlopeOcc","",NPairTypes,-0.5,NPairTypes-0.5); hsl->SetDirectory(0);
1655         for (int ipaty=0;ipaty<NPairTypes;ipaty++){
1656             XLoss &X    = xl[ipaty]; if (!X.have) continue;
1657             hsl->SetBinContent(ipaty+1,X.slopex); hsl->SetBinError(ipaty+1,X.eslopex);
1658             hsl->GetXaxis()->SetBinLabel(ipaty+1,pname(ipaty).Data());
1659         }
1660         hsl->SetTitle("occupancy slope df_{A}/d(#LTn_{1}#GT#LTn_{2}#GT/mean) per pairtype (near half, per-chunk errors);;slope");
1661         hsl->SetMarkerStyle(20); hsl->SetMarkerSize(0.8);
1662         hsl->Draw("E1 X0");
1663         gPad->Update();
1664         { TLine *l0 = new TLine(gPad->GetUxmin(),0.,gPad->GetUxmax(),0.); l0->SetLineColor(kGray+1); l0->SetLineStyle(2); l0->Draw(); }
1665         extraW.push_back(hsl);
1666         cfin->Print(PdfFileName.Data());
1667         for (int ipaty=0;ipaty<NPairTypes;ipaty++){
1668             XLoss &X    = xl[ipaty]; if (!X.have) continue;
1669             cout<<Form("corral::Finalize -- step 1b fits pairtype %2d %-18s: const chi2/ndf %.1f/%.0f; slope vs chunk %+.2e +- %.1e; slope vs occupancy %+.2e +- %.1e",
1670                 ipaty,pname(ipaty).Data(),X.chi2c,X.ndfc,X.slopec,X.eslopec,X.slopex,X.eslopex)<<endl;
1671         }
1672     }
1673     //---- N_valid pages (4x4 pairtypes each): one per map name without the pairtype suffix
1674     {
1675         std::vector<TString> bases;
1676         for (TH1* h : extra){ TString b = h->GetName(); b.Remove(b.Last('_')); if (std::find(bases.begin(),bases.end(),b)==bases.end()) bases.push_back(b); }
1677         for (const TString& b : bases){
1678             cfin->Clear(); cfin->Divide(4,4,0.0001,0.0001);
1679             for (TH1* h : extra){
1680                 TString nm  = h->GetName();
1681                 TString bb  = nm; bb.Remove(bb.Last('_'));
1682                 if (bb!=b) continue;
1683                 int ipaty   = TString(nm(nm.Last('_')+1,nm.Length())).Atoi();
1684                 cfin->cd(1+ipaty);
1685                 h->SetMinimum(0.); h->SetMaximum(nlists);
1686                 h->Draw("colz");
1687             }
1688             cfin->Print(PdfFileName.Data());
1689         }
1690     }
1691     cfin->Print((PdfFileName+"]").Data());
1692     fout->cd();
1693     //---- the full-stats results under their corral names (hmult_k, hC2_1_k, hR2C_1_k, ...; CF maps with subgroup
1694     //---- errors), always -- not only when a single-job reference exists (user 2026-09-30)
1695     //---- top level = the physics set per pairtype, corral names (user 2026-09-30); everything else, incl. all hFin_*, in work/
1696     auto isPhys = [](const TString& n){
1697         if (n.Last('_')<0) return false;
1698         TString t   = n; t.Remove(t.Last('_'));
1699         return (t=="hmult" || t=="hC2C_1" || t=="hR2C_1");
1700     };
1701     //---- top-level copies of the projections and C(Q) under corral-style names (hFin_* originals stay in work/)
1702     auto physCopy   = [](const TString& n) -> TString {
1703         static const char* MAP[][2] = { {"hFin_C2dy_C","hC2Cdy_1"}, {"hFin_C2dphi_C","hC2Cdphi_1"},
1704                                         {"hFin_R2dy_C","hR2Cdy_1"}, {"hFin_R2dphi_C","hR2Cdphi_1"},
1705                                         {"hFin_CQ","hCQ"}, {"hFin_CQ5","hCQ5"}, {"hFin_CQ20","hCQ20"}, {"hFin_CQ20n","hCQ20n"},
1706                                         {"hFin_CQvar","hCQvar"}, {"hFin_CQvarPT","hCQvarPT"} };
1707         if (n.Last('_')<0) return TString("");
1708         TString t   = n, k = n(n.Last('_')+1,n.Length()); t.Remove(t.Last('_'));
1709         for (auto &m : MAP) if (t==m[0]) return TString(m[1])+"_"+k;
1710         return "";
1711     };
1712     TDirectory *dwork   = fout->mkdir("work");
1713     //---- work/chunkStats (2026-10-11, comparison/from_laptop_20261011_chunkStats_spec.md + the f2 request): one entry per
1714     //---- chunk, from the chunk's hmult_z<zz>_<k> (events reaching CalcRm, n = 0 included) for every pair type k. Index
1715     //---- [s]: 0 = |zvtx| < 8 cm (whole zvtx slices), 1 = all slices. Per pair type: 1+r2 = <n1 n2>/(<n1><n2>) (distinct)
1716     //---- or <n(n-1)>/<n>^2 (identical), its error (delta method from the chunk's own hmult), f2 = <n1 n2> - <n1><n2> or
1717     //---- <n(n-1)> - <n>^2 (pairs per event), <n1>, <n2>, P0 of each. Per particle ID: <n>, P0, P(n) (|zvtx| < 8, n >= NPN-1
1718     //---- in the last bin). Totals: N = sum of n over pi+-, K+-, p, pbar ("id", 21 hmults), over pi+- ("pi"), over h+- ("h");
1719     //---- R2tot = (<N(N-1)> - <N>^2)/<N>^2 (each species in its own y window: Species_yu). Occupancy: <nvtx> (raw vertices
1720     //---- in the TF) of the events reaching CalcRm (hnvtxCR; -1 in chunks made before 2026-10-11). Pair type names in
1721     //---- the tree's UserInfo (TObjString "k name"). Whole-sample totals with subgroup errors: log + work/R2tot_* parameters.
1722     {
1723         const int NP    = NPairTypes, NPN = 20;
1724         TH1 *hn0c   = (TH1*)fch[0]->Get("hnevtz_0");
1725         const int NZc   = hn0c ? hn0c->GetNbinsX() : 0;
1726         if (NZc<=0){ cout<<"corral::Finalize -- chunkStats: hnevtz_0 missing, no chunkStats"<<endl; }
1727         struct Mom { double N=0, A=0, B=0, C=0, AA=0, BB=0, CC=0, AB=0, AC=0, BC=0, P01=0, P02=0, Pn1[20]={0}, Pn2[20]={0}; };
1728         auto addMom = [&](Mom &m, TH2 *h, bool ident){
1729             const TArrayD *arr  = dynamic_cast<TArrayD*>(h);        // hmult is a TH2D: scan its cells directly (mostly empty)
1730             const int nxc       = h->GetNbinsX()+2;
1731             for (int ib=0;ib<arr->GetSize();ib++){
1732                 const double w  = arr->At(ib); if (w==0.) continue;
1733                 const int i = ib%nxc, j = ib/nxc;
1734                 if (i<1 || i>h->GetNbinsX() || j<1 || j>h->GetNbinsY()) continue;
1735                 const double b  = h->GetXaxis()->GetBinCenter(i), c = h->GetYaxis()->GetBinCenter(j);
1736                 const double a  = ident ? b*(b-1.) : b*c;
1737                 m.N += w; m.A += w*a; m.B += w*b; m.C += w*c;
1738                 m.AA += w*a*a; m.BB += w*b*b; m.CC += w*c*c; m.AB += w*a*b; m.AC += w*a*c; m.BC += w*b*c;
1739                 if (b<0.5) m.P01 += w;
1740                 if (c<0.5) m.P02 += w;
1741                 m.Pn1[std::min(NPN-1,(int)lround(b))] += w; m.Pn2[std::min(NPN-1,(int)lround(c))] += w;
1742             }
1743         };
1744         //---- 1+r2 and its delta-method error: F = <A>/(<B><C>), gradient (1/BC, -A/B^2C, -A/BC^2), covariance of the means / N
1745         auto opr    = [](const Mom &m, double &e){
1746             e = 0.; if (m.N<2. || m.B<=0. || m.C<=0.) return -1.;
1747             const double a = m.A/m.N, b = m.B/m.N, c = m.C/m.N, F = a/(b*c);
1748             const double g[3] = {1./(b*c), -a/(b*b*c), -a/(b*c*c)};
1749             const double S[3][3] = {{m.AA/m.N-a*a, m.AB/m.N-a*b, m.AC/m.N-a*c},
1750                                     {m.AB/m.N-a*b, m.BB/m.N-b*b, m.BC/m.N-b*c},
1751                                     {m.AC/m.N-a*c, m.BC/m.N-b*c, m.CC/m.N-c*c}};
1752             double v = 0.; for (int i=0;i<3;i++) for (int j=0;j<3;j++) v += g[i]*S[i][j]*g[j];
1753             e = sqrt(std::max(0.,v)/m.N); return F;
1754         };
1755         //---- per particle ID: the pair type whose hmult gives its own n (first match, either side); per particle pair:
1756         //---- the pair type with that pair (either order) -- found by ID, so the pair-type order does not matter
1757         std::vector<int> kOf(NPARTICLENAMES,-1), sideOf(NPARTICLENAMES,0);
1758         for (int ip=0;ip<NPARTICLENAMES;ip++) for (int k=0;k<NP && kOf[ip]<0;k++){
1759             if (PairTypes_Info[k][0]==ip){ kOf[ip] = k; sideOf[ip] = 0; } else if (PairTypes_Info[k][1]==ip){ kOf[ip] = k; sideOf[ip] = 1; }
1760         }
1761         auto kPair  = [&](int a, int b){ for (int k=0;k<NP;k++) if ((PairTypes_Info[k][0]==a && PairTypes_Info[k][1]==b) || (PairTypes_Info[k][0]==b && PairTypes_Info[k][1]==a)) return k; return -1; };
1762         const std::vector<int> setId    = {kParticleIDPionPlus,kParticleIDPionMinus,kParticleIDKaonPlus,kParticleIDKaonMinus,kParticleIDProton,kParticleIDAntiProton};
1763         const std::vector<int> setPi    = {kParticleIDPionPlus,kParticleIDPionMinus};
1764         const std::vector<int> setH     = {kParticleIDhPlus,kParticleIDhMinus};
1765         //---- tree branches
1766         Int_t   b_run, b_chunk;
1767         Double_t    b_nevt[2], b_occ[2];
1768         std::vector<Double_t> b_opr(2*NP), b_eopr(2*NP), b_f2(2*NP), b_n1(2*NP), b_n2(2*NP), b_P01(2*NP), b_P02(2*NP);
1769         std::vector<Double_t> b_nm(2*NPARTICLENAMES), b_P0(2*NPARTICLENAMES), b_Pn(NPARTICLENAMES*NPN);
1770         Double_t    b_Ntot[3][2], b_f2tot[3][2], b_R2tot[3][2];     // [id, pi, h][s]
1771         dwork->cd();
1772         TTree *tcs  = new TTree("chunkStats","per-chunk multiplicity moments (README.md, Finalize: work/chunkStats)");
1773         if (NZc>0){
1774             tcs->Branch("run",&b_run,"run/I");
1775             tcs->Branch("chunk",&b_chunk,"chunk/I");
1776             tcs->Branch("nevt",b_nevt,"nevt[2]/D");
1777             tcs->Branch("occ",b_occ,"occ[2]/D");
1778             tcs->Branch("opr",b_opr.data(),Form("opr[2][%d]/D",NP));
1779             tcs->Branch("eopr",b_eopr.data(),Form("eopr[2][%d]/D",NP));
1780             tcs->Branch("f2",b_f2.data(),Form("f2[2][%d]/D",NP));
1781             tcs->Branch("n1",b_n1.data(),Form("n1[2][%d]/D",NP));
1782             tcs->Branch("n2",b_n2.data(),Form("n2[2][%d]/D",NP));
1783             tcs->Branch("P0_1",b_P01.data(),Form("P0_1[2][%d]/D",NP));
1784             tcs->Branch("P0_2",b_P02.data(),Form("P0_2[2][%d]/D",NP));
1785             tcs->Branch("nmean",b_nm.data(),Form("nmean[2][%d]/D",NPARTICLENAMES));
1786             tcs->Branch("P0",b_P0.data(),Form("P0[2][%d]/D",NPARTICLENAMES));
1787             tcs->Branch("Pn",b_Pn.data(),Form("Pn[%d][%d]/D",NPARTICLENAMES,NPN));
1788             tcs->Branch("Ntot",b_Ntot,"Ntot[3][2]/D");
1789             tcs->Branch("f2tot",b_f2tot,"f2tot[3][2]/D");
1790             tcs->Branch("R2tot",b_R2tot,"R2tot[3][2]/D");
1791             for (int k=0;k<NP;k++) tcs->GetUserInfo()->Add(new TObjString(Form("%d %s%s",k,ParticleIDNames[PairTypes_Info[k][0]],ParticleIDNames[PairTypes_Info[k][1]])));
1792             for (int ip=0;ip<NPARTICLENAMES;ip++) tcs->GetUserInfo()->Add(new TObjString(Form("particle %d %s",ip,ParticleIDNames[ip])));
1793         }
1794         std::vector<double> sumR2tot[3][2];     // per chunk, for the subgroup errors
1795         std::vector<Mom> msum(2*NP);            // whole sample
1796         int nocc    = 0;
1797         for (int ic=0;ic<nlists && NZc>0;ic++){
1798             std::vector<Mom> m(2*NP);
1799             for (int k=0;k<NP;k++){
1800                 const bool ident    = (PairTypes_Info[k][0]==PairTypes_Info[k][1]);
1801                 for (int z=0;z<NZc;z++){
1802                     TH2 *h  = (TH2*)fch[ic]->Get(Form("hmult_z%02d_%d",z,k));
1803                     if (!h){ cout<<"corral::Finalize -- chunkStats: hmult_z"<<z<<"_"<<k<<" missing in chunk "<<ic<<", exit"<<endl; exit(1); }
1804                     const double zlo = -16.+32.*z/NZc, zhi = -16.+32.*(z+1)/NZc;
1805                     if (zlo>=-8.-1e-9 && zhi<=8.+1e-9) addMom(m[k],h,ident);
1806                     addMom(m[NP+k],h,ident);
1807                     delete h;
1808                 }
1809             }
1810             b_run   = chunkRun[ic]; b_chunk = ic;
1811             for (int s=0;s<2;s++) b_nevt[s] = m[s*NP].N;
1812             TH2 *hv = (TH2*)fch[ic]->Get("hnvtxCR");
1813             if (hv){
1814                 for (int s=0;s<2;s++){ TH1 *p = hv->ProjectionY(Form("hnvtxCR_p%d",s), s==0 ? hv->GetXaxis()->FindBin(-7.99) : 0, s==0 ? hv->GetXaxis()->FindBin(7.99) : -1); b_occ[s] = p->GetMean(); delete p; }
1815                 ++nocc; delete hv;
1816             } else { b_occ[0] = b_occ[1] = -1.; }
1817             for (int s=0;s<2;s++) for (int k=0;k<NP;k++){
1818                 const Mom &q    = m[s*NP+k];
1819                 double e; b_opr[s*NP+k] = opr(q,e); b_eopr[s*NP+k] = e;
1820                 const double n1 = q.N>0 ? q.B/q.N : 0., n2 = q.N>0 ? q.C/q.N : 0.;
1821                 b_f2[s*NP+k]    = q.N>0 ? q.A/q.N - n1*n2 : 0.;
1822                 b_n1[s*NP+k] = n1; b_n2[s*NP+k] = n2;
1823                 b_P01[s*NP+k] = q.N>0 ? q.P01/q.N : 0.; b_P02[s*NP+k] = q.N>0 ? q.P02/q.N : 0.;
1824                 Mom &t = msum[s*NP+k]; t.N += q.N; t.A += q.A; t.B += q.B; t.C += q.C;
1825             }
1826             for (int ip=0;ip<NPARTICLENAMES;ip++){
1827                 for (int s=0;s<2;s++){
1828                     b_nm[s*NPARTICLENAMES+ip] = -1.; b_P0[s*NPARTICLENAMES+ip] = -1.;
1829                     if (kOf[ip]<0) continue;
1830                     const Mom &q = m[s*NP+kOf[ip]]; if (q.N<=0.) continue;
1831                     b_nm[s*NPARTICLENAMES+ip]   = (sideOf[ip]==0 ? q.B : q.C)/q.N;
1832                     b_P0[s*NPARTICLENAMES+ip]   = (sideOf[ip]==0 ? q.P01 : q.P02)/q.N;
1833                 }
1834                 for (int n=0;n<NPN;n++){ b_Pn[ip*NPN+n] = -1.; if (kOf[ip]<0) continue; const Mom &q = m[kOf[ip]]; if (q.N>0.) b_Pn[ip*NPN+n] = (sideOf[ip]==0 ? q.Pn1[n] : q.Pn2[n])/q.N; }
1835             }
1836             //---- totals: <N> = sum <n_a>, <N(N-1)> = sum_a <n_a(n_a-1)> + 2 sum_{a<b} <n_a n_b>
1837             const std::vector<int>* sets[3] = {&setId,&setPi,&setH};
1838             for (int is=0;is<3;is++) for (int s=0;s<2;s++){
1839                 const std::vector<int> &S = *sets[is];
1840                 double Nm = 0., NN = 0.; bool ok = true;
1841                 for (size_t i=0;i<S.size() && ok;i++){
1842                     const int ka = kPair(S[i],S[i]); if (ka<0 || m[s*NP+ka].N<=0.){ ok = false; break; }
1843                     Nm += m[s*NP+ka].B/m[s*NP+ka].N; NN += m[s*NP+ka].A/m[s*NP+ka].N;
1844                     for (size_t j=i+1;j<S.size();j++){ const int kb = kPair(S[i],S[j]); if (kb<0){ ok = false; break; } NN += 2.*m[s*NP+kb].A/m[s*NP+kb].N; }
1845                 }
1846                 b_Ntot[is][s] = ok ? Nm : -1.; b_f2tot[is][s] = ok ? NN - Nm*Nm : -1.; b_R2tot[is][s] = (ok && Nm>0.) ? (NN - Nm*Nm)/(Nm*Nm) : -1.;
1847                 if (ok && Nm>0.) sumR2tot[is][s].push_back(b_R2tot[is][s]);
1848             }
1849             tcs->Fill();
1850         }
1851         if (NZc>0){
1852             cout<<"corral::Finalize -- chunkStats: "<<tcs->GetEntries()<<" chunks, occupancy (hnvtxCR) in "<<nocc<<" of them; written to work/chunkStats"<<endl;
1853             //---- whole-sample totals: mean over chunks with the subgroup error (chunk-to-chunk scatter / sqrt(n))
1854             const char* SETNAME[3] = {"id","pi","h"};
1855             for (int is=0;is<3;is++) for (int s=0;s<2;s++){
1856                 const std::vector<double> &v = sumR2tot[is][s]; if (v.size()<2) continue;
1857                 double a = 0., a2 = 0.; for (double x : v){ a += x; a2 += x*x; } const double n = v.size(), mu = a/n;
1858                 const double err = sqrt(std::max(0.,(a2-a*a/n)/(n-1.)))/sqrt(n);
1859                 cout<<Form("corral::Finalize -- chunkStats: R2tot(%s, %s) = %.4f +- %.4f (subgroup, %d chunks)",SETNAME[is],s==0?"|zvtx|<8":"all zvtx",mu,err,(int)n)<<endl;
1860                 extraW.push_back(new TParameter<double>(Form("R2tot_%s_%s",SETNAME[is],s==0?"z8":"zall"),mu));
1861                 extraW.push_back(new TParameter<double>(Form("eR2tot_%s_%s",SETNAME[is],s==0?"z8":"zall"),err));
1862             }
1863             for (int k=0;k<NP;k++){
1864                 const Mom &t = msum[k]; if (t.N<=0. || t.B<=0. || t.C<=0.) continue;
1865                 cout<<Form("corral::Finalize -- chunkStats: pairtype %2d %s%s |zvtx|<8: 1+r2 (hmult) = %.4f, f2 = %.4f pairs/event, <n1> %.3f <n2> %.3f",
1866                     k,ParticleIDNames[PairTypes_Info[k][0]],ParticleIDNames[PairTypes_Info[k][1]],(t.A/t.N)/((t.B/t.N)*(t.C/t.N)),t.A/t.N-(t.B/t.N)*(t.C/t.N),t.B/t.N,t.C/t.N)<<endl;
1867             }
1868         }
1869         if (NZc>0) tcs->Write(); else delete tcs;
1870     }
1871     //---- info/ (README.md): the chunks' record (first chunk's texts, counters summed over the chunks) + this Finalize
1872     {
1873         TDirectory *di  = fout->mkdir("info","how this output was made (README.md, info/)"); di->cd();
1874         TNamed("runstring",finInfo.runstring.Data()).Write();
1875         TNamed("settings",finInfo.settings.Data()).Write();
1876         TNamed("code",finInfo.code.Data()).Write();
1877         TNamed("segveto",finInfo.segveto.Data()).Write();
1878         TNamed("list",finInfo.listdir.Data()).Write();
1879         TString mm; for (auto &m : finInfo.mismatch) mm += m+"\n";
1880         TNamed("chunks",Form("%d chunks of %s, %d with info/%s%s",finInfo.nchunks,SET.Data(),finInfo.nwith,mm!="" ? "; " : "; all agree",mm.Data())).Write();
1881         TNamed("finalize",Form("run string %s, physics |zvtx| < %.0f cm, git %s (commit date %s), built %s",RunString.Data(),valFinZMax,CORRAL_GIT,CORRAL_GIT_DATE,CORRAL_BUILD_TIME)).Write();
1882         for (auto &c : finInfo.counts) TParameter<double>(("count_"+c.first).Data(),c.second).Write();
1883         for (auto &c : finInfo.cuts) TParameter<double>(("cut_"+c.first).Data(),c.second).Write();
1884         TDirectory *dq  = dwork->mkdir("qa","QA histograms summed over the chunks (the QA pages)"); dq->cd();
1885         for (TH1 *h : qaHists) h->Write(TString(h->GetName()).ReplaceAll("_fin","").Data());
1886         if (segQA){ dwork->cd(); segQA->Write(); }
1887     }
1888     for (auto &r : result) if (r.second){ r.second->SetLineColor(kBlack); (isPhys(r.first) ? (TDirectory*)fout : dwork)->cd(); r.second->Write(r.first.Data()); }
1889     dwork->cd();
1890     for (TH1* h : extra) h->Write();
1891     for (TObject* o : extraW){
1892         dwork->cd(); o->Write();
1893         TString nc  = physCopy(o->GetName());
1894         if (nc!=""){ fout->cd(); o->Write(nc.Data()); }
1895     }
1896     fout->cd();
1897     fout->Close();
1898     cout<<"corral::Finalize -- wrote "<<RootFileName<<" and "<<PdfFileName<<endl;
1899     cout<<"done."<<endl;
1900 }