macros/wEEC_3D_unfolding_kinematics/fillData.C

0001 #include "analysisHelper.h"
0002
0003 // ─────────────────────────────────────────────────────────────────────────────
0004 //  fillData.C
0005 //
0006 //  Fills measured (reco-only) histograms from real data.
0007 //  No truth information, no response matrix, no pT-hat cut.
0008 //
0009 //  Output — dijet pT:
0010 //    hJetPt_meas       — measured reco dijet pT (flat index, unfolding input)
0011 //
0012 //  Output — wEEC 2D per reco dijet pT bin fr:
0013 //    hWEEC2D_data_fr   — (ΔΦ, reco pair weight) histogram
0014 //                        pair weight = pT_i * pT_j / <pT>²
0015 //                        using vertex-shifted eta for reco towers
0016 // ─────────────────────────────────────────────────────────────────────────────
0017 void fillData(int runnumber = 47722, int seg = 0, const char* outDir = ".", float towCut = 0.25)
0018 {
0019     TH1D* hJetPt_meas = new TH1D("hJetPt_meas", "", nRecoFlat(), 0, nRecoFlat());
0020     hJetPt_meas->Sumw2();
0021
0022     // Vtx distribution from data — used by makeVtxWeights.C to form the
0023     // data/MC ratio that is stored in analysisHelper.h.
0024     TH1D* hVtxData = new TH1D("hVtxData",
0025         "Data vtx_z;v_{tz} (cm);events",
0026         nVtxBins, vtxZMin, vtxZMax);
0027     // No Sumw2: data events are unweighted counts.
0028
0029     // 3D data histograms: one per ΔΦ bin, in reco-index space
0030     std::vector<TH1D*> hWEEC3D_data(nDphi, nullptr);
0031     for (int k = 0; k < nDphi; ++k) {
0032         hWEEC3D_data[k] = new TH1D(
0033             std::format("hWEEC3D_data_{}", k).c_str(),
0034             std::format("Data wEEC 3D k{};flat(iL,iS,iPw);pairs", k).c_str(),
0035             nRecoFlat3D(), 0, nRecoFlat3D());
0036         hWEEC3D_data[k]->Sumw2();
0037     }
0038
0039     // ── open input ────────────────────────────────────────────────────────
0040     std::string dataFile = std::format(
0041         "/sphenix/tg/tg01/jets/bkimelman/VandyDSTs_Apr4_2026/"
0042         "VandyDSTs_run2pp_ana521_2025p007_v001-{:08}-{}_{}.root",
0043         runnumber, seg, seg + 24);
0044     TFile* fIn = new TFile(dataFile.c_str(), "READ");
0045     if (!fIn || fIn->IsZombie()) {
0046         std::cerr << "Cannot open " << dataFile << "\n";
0047         if (fIn) delete fIn; return;
0048     }
0049     TTree* T = (TTree*)fIn->Get("T");
0050     if (!T) { fIn->Close(); delete fIn; return; }
0051
0052     EventInfo*          eventInfo = nullptr;
0053     std::vector<Tower>* towers    = nullptr;
0054     T->SetBranchAddress("EventInfo", &eventInfo);
0055     T->SetBranchAddress("TowerInfo", &towers);
0056
0057     long nEvents = T->GetEntries();
0058     for (long ev = 0; ev < nEvents; ++ev)
0059     {
0060         std::cout << "working on event " << ev << std::endl;
0061         T->GetEntry(ev);
0062
0063         double vtx_z = eventInfo->get_z_vtx();
0064         if (std::abs(vtx_z) > 10.0) continue;
0065
0066         // Fill vtx distribution for all events passing the vtx cut,
0067         // regardless of whether they pass the dijet selection.
0068         hVtxData->Fill(vtx_z);
0069
0070         if (!eventInfo->is_dijet_event(2)) continue;
0071         double rLeadPT = eventInfo->get_lead_pT(2);
0072         double rSublPT = eventInfo->get_sublead_pT(2);
0073         if (rLeadPT < recoLeadMin || rSublPT < recoSublMin) continue;
0074
0075
0076         int iLreco = FindBin(rLeadPT, recoLeadPtBins);
0077         int iSreco = FindBin(rSublPT, recoSublPtBins);
0078         if (iLreco < 0 || iSreco < 0) continue;
0079
0080         int    fr      = RecoFlatIndex(iLreco, iSreco);
0081         double rFlat   = RecoFlatBinCenter(fr);
0082         hJetPt_meas->Fill(rFlat);   // no MC weight for data
0083
0084         double rPTmean2 = 0.25 * (rLeadPT + rSublPT) * (rLeadPT + rSublPT);
0085
0086         double recoMap[24][64] = {};
0087         shiftEtaEdges(vtx_z);
0088
0089         for (auto& tow : *towers) {
0090             if (tow.get_calo() < 1 || tow.get_calo() > 3) continue;
0091             fastjet::PseudoJet tj(tow.px(), tow.py(), tow.pz(), tow.e());
0092             int etaBin = getEtaBin(tj.pseudorapidity(), vtx_z);
0093             int phiBin = getPhiBin(tj.phi());
0094             if (etaBin >= 0 && phiBin >= 0)
0095                 recoMap[etaBin][phiBin] += tow.e();
0096         }
0097
0098         for (int i = 0; i < 24*64; ++i) {
0099             int ei = i / 64, pi = i % 64;
0100             double phi_i = getPhiCenter(pi);
0101             double rpT_i = recoMap[ei][pi] / cosh(getEtaCenter(ei, vtx_z));
0102             if (rpT_i <= towCut || rpT_i >= 80) continue;
0103
0104             for (int j = i+1; j < 24*64; ++j) {
0105                 int ej = j / 64, pj = j % 64;
0106                 double phi_j = getPhiCenter(pj);
0107                 double rpT_j = recoMap[ej][pj] / cosh(getEtaCenter(ej, vtx_z));
0108                 if (rpT_j <= towCut || rpT_j >= 80) continue;
0109
0110                 double dphi   = DeltaPhi(phi_i, phi_j);
0111                 double rPairW = rpT_i * rpT_j / rPTmean2;
0112                 if (rPairW < pairWeightMin || rPairW > pairWeightMax) continue;
0113                 int iDphi = FindBin(dphi, dPhiBins);
0114                 int iPw   = FindBin(rPairW, pairWeightBins);
0115                 if (iDphi < 0 || iPw < 0) continue;
0116                 hWEEC3D_data[iDphi]->Fill(
0117                     RecoFlat3DBinCenter(RecoFlat3DIndex(iLreco, iSreco, iPw)));
0118             }
0119         }
0120     }
0121
0122     fIn->Close(); delete fIn;
0123
0124     std::string outName = std::format(
0125         "{}/Data/data_measured_{:08}_seg{:06d}_to_{:06d}.root",
0126         outDir, runnumber, seg, seg + 24);
0127
0128     std::cout << "writing file " << outName << std::endl;
0129
0130     TFile* fOut = new TFile(outName.c_str(), "RECREATE");
0131     hJetPt_meas->Write();
0132     for (int k = 0; k < nDphi; ++k)
0133         hWEEC3D_data[k]->Write();
0134     hVtxData->Write();
0135     fOut->Close();
0136     std::cout << "Written: " << outName << "\n";
0137 }