Corral/lists/make_bad_segments.C

0001 //---- make_bad_segments.C -- the segment mask for corral "segveto" (README.md, README_ana573.md), from the segCells trees of a
0002 //---- production's chunk files (the acceptance pass, run WITHOUT segveto).
0003 //----   root -l -b -q 'make_bad_segments.C+("ana573_795xx,ana573_8156x,ana573_8158x")'   (from lists/; reads ../root/<set>/chunks/)
0004 //---- Writes lists/<set>/bad_segments.txt for every set given (and lists/ana573/ when several sets are given): one
0005 //---- "run segment" line per LOW segment; "# band" header lines give each run's band and threshold (SegCellsView draws them).
0006 //---- Rule (user 2026-10-12; the live fraction is computed in src/SegLive.h, shared with Finalize):
0007 //----   x_c = TPC hits in cell c per accepted track in the segment; r_c = 90th percentile of x_c over all segments; cells with
0008 //----   r_c < 0.1 x the median r_c are never live and left out. A cell is dead in a segment if x_c < 0.5 r_c.
0009 //----   live fraction = 1 - (sum of r_c over dead cells) / (sum of r_c over all cells): all, South (side 0), North (side 1).
0010 //----   Per run and view, band = mode of the per-segment live fraction (0.001 bins) +- BANDW, sigma = rms inside the band.
0011 //----   A segment is LOW if, in any view, live < mode - max(LOWMIN, LOWNS sigma).
0012 //---- Cells: per FEE (branch "fee", tpcsector*26 + feeMap) or the older generic grid (branch "hits", 48 per sector-side).
0013 #include "../src/SegLive.h" // the live-fraction rule (shared with Finalize's vs-segment page)
0014 #include <TH1D.h>
0015 #include <TObjString.h>
0016 #include <TObjArray.h>
0017 #include <TSystem.h>
0018 #include <cmath>
0019
0020 //---- write = false: print the per-run counts only, write no bad_segments.txt (a check)
0021 void make_bad_segments(const char* setlist = "ana573_795xx,ana573_8156x,ana573_8158x", bool write = true){
0022     const double BANDW = 0.02, LOWMIN = 0.005, LOWNS = 4.;
0023     TString ldir    = gSystem->DirName(gSystem->UnixPathName(__FILE__));    // lists/
0024     TObjArray *ts   = TString(setlist).Tokenize(",");
0025     const int NS    = ts->GetEntries();
0026     std::vector<TString> setName, files;
0027     for (int s=0;s<NS;s++){
0028         TString sn = ((TObjString*)ts->At(s))->GetName(); setName.push_back(sn);
0029         for (int c=0;c<1000;c++){ TString p = ldir+"/../root/"+sn+Form("/chunks/corral_%02d.root",c); if (gSystem->AccessPathName(p)) break; files.push_back(p); }
0030     }
0031     if (files.empty()){ printf("no chunk files\n"); return; }
0032     SegLiveResult SL    = SegLiveFractions(files);
0033     const bool feeMode  = SL.feeMode;
0034     auto &live          = SL.live;
0035     printf("%zu chunk files, cells %s\n",files.size(),feeMode ? "per FEE" : "generic grid");
0036     printf("%zu segments, %d never-live cells\n",live.size(),SL.nNever);
0037     //---- band per run and view, LOW segments
0038     std::map<int,std::vector<std::pair<int,std::array<float,3>>>> byRun;
0039     for (auto &p : live) byRun[p.first.first].push_back({p.first.second,p.second});
0040     std::vector<TString> header; std::vector<std::pair<int,int>> low;
0041     for (auto &R : byRun){
0042         std::array<double,3> thr;
0043         for (int a=0;a<3;a++){
0044             TH1D hb("hb","",400,0.6,1.0); for (auto &s : R.second) hb.Fill(std::min(0.99999,(double)s.second[a]));
0045             double mode = hb.GetBinCenter(hb.GetMaximumBin()), s1 = 0, s2 = 0; int n = 0;
0046             for (auto &s : R.second){ double v = s.second[a]; if (fabs(v-mode)<BANDW){ s1 += v; s2 += v*v; n++; } }
0047             double sig = n>1 ? sqrt(std::max(0.,s2/n-pow(s1/n,2))) : 0.;
0048             thr[a] = mode - std::max(LOWMIN,LOWNS*sig);
0049             header.push_back(Form("# band run %d view %s mode %.4f sigma %.4f thr %.4f",R.first,a==0?"all":a==1?"S":"N",mode,sig,thr[a]));
0050         }
0051         int nl = 0;
0052         for (auto &s : R.second){ bool l = false; for (int a=0;a<3;a++) if (s.second[a]<thr[a]) l = true; if (l){ low.push_back({R.first,s.first}); nl++; } }
0053         printf("run %d: %d of %zu segments low (%.1f%%)\n",R.first,nl,R.second.size(),100.*nl/R.second.size());
0054     }
0055     //---- write
0056     printf("%zu of %zu segments low\n",low.size(),live.size());
0057     if (!write){ printf("write = false: no bad_segments.txt written\n"); return; }
0058     std::vector<TString> dirs = setName; if (NS>1) dirs.push_back("ana573");
0059     for (auto &d : dirs){
0060         if (gSystem->AccessPathName(ldir+"/"+d+"/manifest.txt")) continue;
0061         TString fn = ldir+"/"+d+"/bad_segments.txt";
0062         FILE *fo = fopen(fn.Data(),"w");
0063         fprintf(fo,"# run segment -- below the run's live-fraction band in any view (all/S/N): live < mode - max(%.3f, %.0f sigma), band = mode +- %.2f\n",LOWMIN,LOWNS,BANDW);
0064         fprintf(fo,"# make_bad_segments.C(\"%s\"), %s cells, %zu of %zu segments; read by corral \"segveto\"\n",setlist,feeMode ? "FEE" : "generic",low.size(),live.size());
0065         for (auto &h : header) fprintf(fo,"%s\n",h.Data());
0066         for (auto &l : low) fprintf(fo,"%d %d\n",l.first,l.second);
0067         fclose(fo); printf("wrote %s\n",fn.Data());
0068     }
0069 }