From 6c594765f806223c21caf9c762c28f0f081f0f1a Mon Sep 17 00:00:00 2001 From: zdemirog Date: Mon, 31 Aug 2026 14:18:15 -0400 Subject: [PATCH] Add macros for background studies --- analysis/background_fill.C | 310 ++++++++++++++++++++++++++++++++++ analysis/background_plotter.C | 117 +++++++++++++ 2 files changed, 427 insertions(+) create mode 100644 analysis/background_fill.C create mode 100644 analysis/background_plotter.C diff --git a/analysis/background_fill.C b/analysis/background_fill.C new file mode 100644 index 000000000..1c6fa2909 --- /dev/null +++ b/analysis/background_fill.C @@ -0,0 +1,310 @@ +// background_fill.C +// Every photon or electron that hits main detector was created somewhere in +// the beamline, in a collimator, shield, quartz tile itself, and so on. This +// macro finds out where each one was born and counts them per component, so we +// can say 'this much of background comes from collimator 2'. +// what this macro does +// 1. loop over every hit on main detector plane and keep the good ones +// (right detector ID, enough energy born outside the target if cut is on, going forwards, +// e+- or gamma); +// 2. take the point where that particle was created, its vertex (vz,vr), and see +// which component box contains it; +// 3. add 1 to that component's counter +// a component here is just a rectangle in (z,r); z-range and r-range that together +// enclose a real object in the geometry. they are listed in the box[][] table below. +// boxes are checked independently so if two of them overlap, a hit inside both is counted in both. +// nothing is normalized here, this macro only counts. it also stores the total number of beam +// events processed, and the cut that was used, so plotting macro (analysis/background_plotter.C) +// can turn counts into a rate later and always know which mode it has. +// how to run +// build/reroot -l -b -q 'background_fill.C("rootfiles.txt","outputs")' +// build/reroot -l -b -q 'background_fill.C("rootfiles.txt","outputs",28, 5, -3440)' +// arg1 = file list +// arg2 = output directory +// arg3 = main detector plane detector ID +// arg4 = signal ring +// arg5 = target cut in mm +// vzTgtCut = -1e9 (default), no cut, keep everything +// vzTgtCut = -3440 drop tracks born in the target +// why target cut? with the blocker out simulation, most of what reaches ring 5 was created +// in the target, that is signal, not background. Dropping tracks born at vz <=vzTgtCut +// leaves only beamline background. + +int sectorOf(double phi){ + const double pi = acos(-1); + if(phi < 0) phi += 2*pi; // atan2 returns [-pi,pi], we want [0,2pi) + double s = fmod(phi, 2*pi/7); // where we are inside one septant + if(s < pi/28) return 0; // 0,1,2 == closed, transition, open + if(s < 3*pi/28) return 1; + if(s < 5*pi/28) return 2; + if(s < 7*pi/28) return 1; + return 0; +} + +//Is the vertex (vz,vr) inside this component's box? +// box = {z_begin, z_end, r_min, r_max}. Lower edge open, upper edge closed, +// so a vertex sitting exactly on shared edge belongs to one box only. +bool inBox(double vz, double vr, double *box){ + return vz > box[0] && vz <= box[1] && vr > box[2] && vr <= box[3]; +} + +void background_fill(TString infile = "rootfiles.txt", + TString outdir = "background_out", // output directory + int MD = 28, // main det plane + int ring = 5, // signal ring for counting + double vzTgtCut = -1e9) //drop tracks born at vz <= this, in mm. +{ + gSystem->mkdir(outdir, kTRUE); + + if(vzTgtCut > -1e8) + cout<<"target cut ON: dropping tracks created at vz <= "< files; + if(infile.EndsWith(".root")) files.push_back(infile); + else {ifstream fl(infile); string line; while (fl>>line) files.push_back(line); } + cout <<"Will process "<IsZombie()){ cout<<"skip, cannot open: " <Get("T"); + if(!T){cout<<"skip, no tree T: "<Close(); continue; } + + std::vector *hit = 0; + T->SetBranchAddress("hit", &hit); + + long nEv = T->GetEntries(); + nPrimary += nEv; // 1 entry = 1 beam event + if(fi%50==0) cout<<"["<GetEntry(i); + + for(int j=0;j<(int)hit->size();j++){ + + remollGenericDetectorHit_t &h = hit->at(j); + + // cut 1: only main detector plane + if(h.det != MD ) continue; + + // cut 2: energy + if(h.e <= EMIN) continue; + + // cut 3: keep electrons/positrons(+-11) and photons (22) + int sp = (abs(h.pid)==11) ? 0 : (h.pid==22 ? 1 : -1); + if (sp < 0) continue; + + // cut 4: drop anything created in target, so only background is left. + if(h.vz <= vzTgtCut) continue; + + // cut5: forward going only + if(h.pz <=0) continue; + + double r = h.r; + + // the vertex, where the particle was created + double vz = h.vz; + double vr = sqrt(h.vx*h.vx + h.vy*h.vy); + + int sec = sectorOf(atan2(h.y,h.x)); // 0 closed, 1 transition, 2 open + + hR[sp]->Fill(r); + hVrVzAll->Fill(vz, vr); + + if(r > rMin[0] && r <= rMax[0]){ + hXY[sp][0]->Fill(h.x,h.y); // all sectors + if(sec==0) hXY[sp][1]->Fill(h.x,h.y); // closed + if(sec==2) hXY[sp][2]->Fill(h.x,h.y); // open + } + + //signal ring only + if(r > rMin[ring] && r <= rMax[ring]){ + + bool inQuartz = false; + for(int c=firstQ; c< nComp; c++) if(inBox(vz,vr,box[c])) inQuartz = true; + + if (!inQuartz) hVz[sp]->Fill(vz); // quartz excluded + hVzAll[sp]->Fill(vz); // everything, quartz included + hVrVz[sp]->Fill(vz,vr); + + // find which components the vertex falls in + bool anyBox = false; + for(int c=0; c< nComp; c++){ + + if(!inBox(vz,vr,box[c])) continue; + anyBox = true; + compCount[sp][c]++; // all sectors + if(sec==0) compCountC[sp][c]++; // closed only + if(sec==2) compCountO[sp][c]++; // open only + } + + // no box at all + if(!anyBox) nUncovered[sp]++; + } + + // z-slice x ring counts, for summary table + for(int s=0; s vzEdge[s] && vz <= vzEdge[s+1]){ + for(int k=0;k<7;k++) + if(r > rMin[k] && r <= rMax[k]) sliceRing[sp][s][k]++; + break; + } + } + } // hits + } // events + + f->Close(); delete f; + + } // files + + //print outs + cout<<"Done. Total beam events processed, nPrimary = " << nPrimary< signal ring hits in no box: " << nUncovered[sp] << " (sum of bars " << tot << ", overlaps counted more than once)" << endl; +} + + // Turn plain counters into histograms, one bin per component + TH1D *hComp[2], *hCompC[2], *hCompO[2]; + TH2D *hSlice[2]; + for(int sp=0; sp<2; sp++){ + hComp[sp] = new TH1D("hComp_"+spName[sp], "component counts;;count",nComp,0,nComp); + hCompC[sp] = new TH1D("hComp_"+spName[sp]+"_closed", "component counts (closed sectors);;count",nComp,0,nComp); + hCompO[sp] = new TH1D("hComp_"+spName[sp]+"_open", "component counts (open sectors);;count",nComp,0,nComp); + for(int c=0; cSetBinContent(c+1,compCount[sp][c]); + hCompC[sp]->SetBinContent(c+1,compCountC[sp][c]); + hCompO[sp]->SetBinContent(c+1,compCountO[sp][c]); + hComp[sp]->GetXaxis()->SetBinLabel(c+1,compName[c]); + hCompC[sp]->GetXaxis()->SetBinLabel(c+1,compName[c]); + hCompO[sp]->GetXaxis()->SetBinLabel(c+1,compName[c]); + } + hSlice[sp] = new TH2D ("hSlice_"+spName[sp],"vz-slice x ring counts; slice; ring",nSlice,0,nSlice,7,0,7); + for(int s=0; sSetBinContent(s+1,k+1, sliceRing[sp][s][k]); + } + + //nPrimary and uncovered counts go in as histograms + TH1D *hNprimary = new TH1D("hNprimary","total beam events processed",1,0,1); + hNprimary->SetBinContent(1, nPrimary); + + //record the cut used + TH1D *hTgtCut = new TH1D("hTgtCut","target vz cut used [mm]; -1e9 means no cut",1,0,1); + hTgtCut->SetBinContent(1, vzTgtCut); + + TH1D *hUncov = new TH1D("hUncovered","ring hits outside every box;species;count",2,0,2); + for(int sp=0;sp<2;sp++){ + hUncov->SetBinContent(sp+1, nUncovered[sp]); + hUncov->GetXaxis()->SetBinLabel(sp+1, spName[sp]); + + } + + //Write everything into ROOT file + TFile *fout = new TFile(outdir+"/background_counts.root","RECREATE"); + hNprimary->Write(); + hTgtCut->Write(); + hUncov->Write(); + hVrVzAll->Write(); + for(int sp=0;sp<2; sp++){ + hR[sp]->Write(); hVz[sp]->Write(); hVzAll[sp]->Write(); hVrVz[sp]->Write(); + hComp[sp]->Write(); hCompC[sp]->Write(); hCompO[sp]->Write(); hSlice[sp]->Write(); + for(int q=0;q<3;q++) hXY[sp][q]->Write(); + } + fout->Close(); + + //cross check plot + gStyle->SetOptStat(0); + TCanvas *cchk = new TCanvas("cchk","boxes check",1100,650); + cchk->SetLogz(); cchk->SetRightMargin(0.13); + hVrVzAll->Draw("colz"); + for(int c=0;c< nComp; c++){ + if(box[c][1] <= -6100 || box[c][0] >= 26100) continue; + TBox *bx = new TBox(box[c][0], box[c][2], box[c][1], box[c][3]); + bx->SetFillStyle(0); bx->SetLineColor(kRed); bx->SetLineWidth(1); bx->Draw("l"); + TText *tt = new TText(0.5*(box[c][0]+box[c][1]),box[c][3], compName[c]); + tt->SetTextSize(0.011); tt->SetTextColor(kRed+1); tt->SetTextAngle(90); tt->Draw(); + + } + cchk->SaveAs(outdir+"/boxes_check.pdf"); + cout<<"Wrote "<SetOptStat(0); + gSystem->mkdir(outdir, kTRUE); + + int cbA = TColor::GetColor("#009E73"); // green = simA + int cbB = TColor::GetColor("#CC79A7"); // purple = simB + TString ssuf = (sector=="") ? "" : "_"+sector; + + //open ROOT files and take component counts + TFile *fA = TFile::Open(fileA), *fB = TFile::Open(fileB); + if(!fA||fA->IsZombie()||!fB||fB->IsZombie()){cout<<"Cannot open one of the inputs"<Get("hComp_"+species+ssuf); + TH1D *hB = (TH1D*)fB->Get("hComp_"+species+ssuf); + if(!hA||!hB){cout<<"Missing hComp_"<Get("hNprimary"); nPrimaryA = n?n->GetBinContent(1):0; } + if(nPrimaryB<=0){ TH1D*n=(TH1D*)fB->Get("hNprimary"); nPrimaryB = n?n->GetBinContent(1):0; } + double NmA = nPrimaryA/beamPerMoller; + double NmB = nPrimaryB/beamPerMoller; + if(NmA<=0||NmB<=0){cout<<"Bad nPrimary (A="<GetNbinsX(); // number of components + + TString spWord = (species=="g") ? "photons" : "e^{#pm}"; + TString secWord = (sector=="") ? "all sectors" : sector+" sectors"; + TString title = "Backgrounds, "+spWord+", E>1 MeV, pz>0, Ring 5, "+secWord; + + //three plots; mode 0 = all, 1 = no quartz, 2 = quartz only + for(int mode=0; mode<3; mode++){ + + //to count how many components this plot will have + std::vectorkeep; + for(int c=0; cGetXaxis()->GetBinLabel(c+1)).BeginsWith("Quartz"); + if(mode==1 && isQuartz) continue; // everything except quartz + if(mode==2 && !isQuartz) continue; // quartz only + keep.push_back(c); // number of bars/components + } + int n = keep.size(); + + // one histogram per simulation + TH1D *bA = new TH1D(Form("bA%d",mode), title+";;Moller Rate %", n,0,n); + TH1D *bB = new TH1D(Form("bB%d",mode), "", n,0,n); + + for(int k=0; kSetBinContent(k+1, 100.0*hA->GetBinContent(c+1)/NmA); + bB->SetBinContent(k+1, 100.0*hB->GetBinContent(c+1)/NmB); + bA->GetXaxis()->SetBinLabel(k+1, hA->GetXaxis()->GetBinLabel(c+1)); + } + + double labSize = 0.022, txtSize=0.012; + if(mode==1){labSize=0.028; txtSize=0.016; } + if(mode==2){labSize=0.040; txtSize=0.026; } + + double ymax = TMath::Max(bA->GetMaximum(), bB->GetMaximum()); + bA->SetMinimum(0); + bA->SetMaximum(1.18*ymax); + bA->GetXaxis()->LabelsOption("v"); + bA->GetXaxis()->SetLabelSize(labSize); + + bA->SetFillColor(cbA); bA->SetBarWidth(0.42); bA->SetBarOffset(0.06); //simA on the left + bB->SetFillColor(cbB); bB->SetBarWidth(0.42); bB->SetBarOffset(0.52); //simB on the right + + //draw + TCanvas *c = new TCanvas(Form("c%d",mode),"compare",1700,850); + c->SetLeftMargin(0.08); c->SetRightMargin(0.03); + c->SetTopMargin(0.08); c->SetBottomMargin(0.22); + bA->Draw("bar"); + bB->Draw("bar same"); + + //print values on top of each bar + TLatex tx; tx.SetTextAngle(90); tx.SetTextAlign(12); tx.SetTextSize(txtSize); + for(int k=0;kGetBinContent(k+1), vb = bB->GetBinContent(k+1); + tx.SetTextColor(cbA); tx.DrawLatex(k+0.27, va+0.01*ymax, Form("%.3f",va)); + tx.SetTextColor(cbB); tx.DrawLatex(k+0.73, vb+0.01*ymax, Form("%.3f",vb)); + } + + TLegend *lg = new TLegend(0.10, 0.80,0.22, 0.917); + lg->SetBorderSize(0); + lg->SetFillStyle(0); + lg->AddEntry(bA,labelA,"f"); + lg->AddEntry(bB,labelB,"f"); + lg->Draw(); + + TString tag = "_compare"; + if(mode==1) tag = "_compare_zoom"; + if(mode==2) tag = "_compare_quartz"; + c->SaveAs(outdir+"/background_"+species+ssuf+tag+".pdf"); + } + cout<<"Wrote "<