From e2eebfd37e29a75b2b7176d552846f1ceac708ed Mon Sep 17 00:00:00 2001 From: Christopher Dilks Date: Tue, 26 Mar 2024 16:09:21 -0400 Subject: [PATCH] feat: calculate gap and overlapping charge from Pass 1 QADB --- util/analyzeGapCharge.C | 50 ++++++++++++++ util/calculateGapCharge.groovy | 118 +++++++++++++++++++++++++++++++++ 2 files changed, 168 insertions(+) create mode 100644 util/analyzeGapCharge.C create mode 100644 util/calculateGapCharge.groovy diff --git a/util/analyzeGapCharge.C b/util/analyzeGapCharge.C new file mode 100644 index 0000000..33d619d --- /dev/null +++ b/util/analyzeGapCharge.C @@ -0,0 +1,50 @@ +// run `calculateGapCharge.groovy` first, then this macro +void analyzeGapCharge(TString datFileN = "charge_gaps.dat") { + auto tr = new TTree("tr", "tr"); + tr->ReadFile(datFileN); + tr->Print(); + + Int_t runset; + Char_t comment[1000]; + Double_t chargeGap; + Double_t chargeGapMin=1e6; + Double_t chargeGapMax=-1e6; + tr->SetBranchAddress("runset", &runset); + tr->SetBranchAddress("comment", comment); + tr->SetBranchAddress("fcChargeGapToNextBin", &chargeGap); + std::map runsets; + runsets[-1] = "full_run_list"; + for(Long64_t e = 0; e < tr->GetEntries(); e++) { + tr->GetEntry(e); + runsets[runset] = comment; + chargeGapMin = std::min(chargeGapMin, chargeGap); + chargeGapMax = std::max(chargeGapMax, chargeGap); + } + std::cout << "chargeGapMin = " << chargeGapMin << std::endl; + std::cout << "chargeGapMax = " << chargeGapMax << std::endl; + + std::map gapDists; + for(auto const& [r, c] : runsets) { + gapDists.insert({r, new TH1D( + Form("gapDist_%d", r), + TString("Gap charge distribution: ") + c + ";gap [nC]", + 100000, + chargeGapMin, + chargeGapMax + )}); + } + + auto canv = new TCanvas(); + canv->Divide(2,2); + int padnum = 1; + for(auto const& [r, c] : runsets) { + TString cut = "golden == 1"; + if(r >= 0) cut += Form(" && runset == %d", r); + tr->Project(gapDists.at(r)->GetName(), "fcChargeGapToNextBin", cut); + auto pad = canv->GetPad(padnum++); + // pad->SetLogy(); + pad->cd(); + gapDists.at(r)->Draw(); + // gapDists.at(r)->GetXaxis()->SetRangeUser(-1000, 1000); + }; +} diff --git a/util/calculateGapCharge.groovy b/util/calculateGapCharge.groovy new file mode 100644 index 0000000..344bba0 --- /dev/null +++ b/util/calculateGapCharge.groovy @@ -0,0 +1,118 @@ +// determine the amount of charge *between* subsequent QA bins +// - for QA using DST 5-files, this may be nonzero +// - for QA using time bins, this should *always* be zero +// - since this is for estimating systematic for cross sections based on RG-A +// Pass 1 data (Valerii, Sangbaek), we include a list of runs +// +// Run this script first, followed by `analyzeGapCharge.C` to make plots + +import org.jlab.io.hipo.HipoDataSource +import clasqa.QADB + +// arguments +if(args.length<1) { + System.err.println "USAGE: run-groovy ${this.class.getSimpleName()}.groovy [RUN NUMBER]" + System.err.println "- set [RUN NUMBER] to 0 to run this for hard-coded list of runs" + System.exit(101) +} +def arg_runnum = args[0].toInteger() +QADB qa = new QADB() + +// define output file +def datfileName = "charge_gaps.dat" +def datfile = new File(datfileName) +def datfileWriter = datfile.newWriter(false) +datfileWriter << '# Runs together with FC charge and gaps between bins\n' +datfileWriter << [ + 'runnum/I', + 'binnum/I', + 'golden/I', + 'fcChargeMin/D', + 'fcChargeMax/D', + 'fcChargeGapToNextBin/D', + 'runset/I', + 'comment/C', +].join(':') << '\n' + +// set run list(s) +def runListHash = [] +if(arg_runnum == 0) { + // Valerii's run lists + runListHash = [ + [ + comment: 'Runs 45 nA', + runs: [ + 5032, 5036, 5038, 5039, 5040, 5041, 5043, 5045, 5052, 5053, 5116, 5117, 5119, 5120, + 5124, 5125, 5126, 5127, 5139, 5153, 5158, 5162, 5163, 5164, 5181, 5191, 5193, 5195, + 5196, 5197, 5198, 5199, 5200, 5201, 5202, 5203, 5204, 5205, 5206, 5208, 5211, 5212, + 5215, 5216, 5219, 5220, 5221, 5222, 5223, 5230, 5231, 5232, 5233, 5234, 5235, 5237, + 5238, 5248, 5345, 5346, 5347, 5349, 5351, 5354, 5355, 5367, + ] + ], + [ + comment: 'Runs 50 - 55 nA', + runs: [ + 5342, 5343, 5344, 5356, 5357, 5358, 5359, 5360, 5361, 5362, 5366, 5368, 5369, 5372, + 5373, 5374, 5375, 5376, 5377, 5378, 5379, 5380, 5381, 5383, 5386, 5390, 5391, 5392, + 5393, 5398, 5401, 5403, 5404, 5406, 5407, + ] + ], + [ + comment: 'Runs with wrong CCDB beam blocker constant 45 nA', + runs: [ + 5249, 5252, 5253, 5257, 5258, 5259, 5261, 5262, 5303, 5304, 5305, 5306, 5307, 5310, + 5311, 5315, 5317, 5318, 5319, 5320, 5323, 5324, 5335, 5339, 5340, + ], + ], + ] +} +else { + runListHash = [[comment: 'Single run', runs: [arg_runnum]]] +} + +// loop over runs +runListHash.eachWithIndex{ hash, runset -> + println "================================================================" + println hash.comment + println "================================================================" + hash.runs.each{ runnum -> + System.out.println runnum + def qaTree = qa.getQaTree()["$runnum"] + def chargeTree = qa.getChargeTree()["$runnum"] + qaTree.each{ binnum, bintree -> + System.out.print " $binnum" + + // the last bin has no subsequent bin, so skip it + def nextBinnum = "${binnum.toInteger()+5}" + if(qaTree[nextBinnum] == null) { + System.out.println " => last bin => skipping" + return + } + System.out.print '\n' + + // check GOLDEN bins + def golden = bintree.defect == 0 + + // calculate the gap charge + // def gapCharge = nextBintree.fcChargeMin - bintree.fcChargeMax + def thisBinChargeTree = chargeTree["$binnum"] + def nextBinChargeTree = chargeTree["$nextBinnum"] + def gapCharge = nextBinChargeTree.fcChargeMin - thisBinChargeTree.fcChargeMax + datfileWriter << [ + runnum, + binnum, + golden ? 1 : 0, + thisBinChargeTree.fcChargeMin, + thisBinChargeTree.fcChargeMax, + gapCharge, + runset, + "\"${hash.comment.replace(' ','_')}\"", + ].join(' ') << '\n' + } + } +} + +// close +datfileWriter.flush() +datfileWriter.close() +System.out.println "Wrote $datfileName"