Skip to content

Commit

Permalink
feat: timeline for beam charge asymmetry (#195)
Browse files Browse the repository at this point in the history
  • Loading branch information
@c-dilks
c-dilks authored Mar 20, 2024
1 parent dd87c8e commit 5a776ee
Show file tree
Hide file tree
Showing 2 changed files with 97 additions and 47 deletions.
114 changes: 79 additions & 35 deletions qa-physics/monitorPlot.groovy
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
// reads outmon/monitor* files and generates timeline hipo files
// - to be executed after monitorRead.groovy

import org.jlab.groot.data.IDataSet
import org.jlab.groot.data.TDirectory
import org.jlab.groot.data.GraphErrors
import org.jlab.groot.data.H1F
Expand Down Expand Up @@ -36,7 +37,6 @@ def qaTreeFile = new File(qaTreeFileN)
def qaTree = slurper.parse(qaTreeFile)

// subroutine to recompute defect bitmask
// FIXME: shamefully copied from QA/modifyQaTree.groovy
def recomputeDefMask = { rnum, bnum ->
def defList = []
def defMask = 0
Expand Down Expand Up @@ -184,6 +184,16 @@ def buildAsymGraph = { tObj ->
return gr
}

// calculate beam charge asymmetry
def calculateBeamChargeAsym = { qP, qM ->
if(qP+qM > 0) {
return [
(qP - qM) / (qP + qM), // asymmetry
1 / Math.sqrt( qP + qM ) // error, assuming |beamChargeAsym| << 1
]
}
return ["unknown","unknown"]
}

//-----------------------------------------
// fill the monitors
Expand All @@ -205,7 +215,7 @@ def aveX
def aveXerr
def stddevX
def ent
def helP,helM,helDef,helUndef,helFrac,helFracErr,rellum,rellumErr
def helDef,helUndef,helFrac,helFracErr,rellum,rellumErr

// loop over input hipo files
inList.each { inFile ->
Expand Down Expand Up @@ -303,10 +313,6 @@ inList.each { inFile ->
monTree[runnum]['helic']['dist']['heldef']['heldefDenom'] += denom
}

// cut for rellum from events with FD trigger electrons only
//}
//if(objN.contains("/helic_distGoodOnly_")) {

// relative luminosity
T.addLeaf(monTree,[runnum,'helic','dist','rellum','rellumNumer'],{0})
T.addLeaf(monTree,[runnum,'helic','dist','rellum','rellumDenom'],{0})
Expand All @@ -326,8 +332,8 @@ inList.each { inFile ->
})
if(obj.integral()>0) {
// use values from helic_dist
helM = obj.getBinContent(0) // helicity = -1
helP = obj.getBinContent(2) // helicity = +1
def helM = obj.getBinContent(0) // helicity = -1
def helP = obj.getBinContent(2) // helicity = +1
// use charge from FC (disabled)
//helP = fcTree[runnum][timeBinNum.toInteger()]['fcP']
//helM = fcTree[runnum][timeBinNum.toInteger()]['fcM']
Expand All @@ -339,9 +345,32 @@ inList.each { inFile ->
monTree[runnum]['helic']['dist']['rellum']['rellumNumer'] += helP
monTree[runnum]['helic']['dist']['rellum']['rellumDenom'] += helM
}

}

// beam charge asymmetry
//----------------------------
if(objN.contains("/helic_scaler_chargeWeighted_")) {
['numHelP','numHelM'].each{
T.addLeaf(monTree,[runnum,'helic','beamChargeAsym',it],{0.0})
}
T.addLeaf(monTree,[runnum,'helic','beamChargeAsym','asymGraph'],{
def g = buildMonAveGr(obj)
def gN = g.getName().replaceAll(/_aveGr$/,'_chargeAsymGr')
g.setName(gN)
g.setTitle('beam charge asymmetry vs. time bin number')
return g
})
if(obj.integral()>0) {
def numHelP = obj.getBinContent(2) // helicity = +1
def numHelM = obj.getBinContent(0) // helicity = -1
monTree[runnum]['helic']['beamChargeAsym']['numHelP'] += numHelP
monTree[runnum]['helic']['beamChargeAsym']['numHelM'] += numHelM
def asym = calculateBeamChargeAsym(numHelP, numHelM)
if(!asym.contains("unknown")) {
monTree[runnum]['helic']['beamChargeAsym']['asymGraph'].addPoint(timeBinNum, asym[0], 0, asym[1])
}
}
}

// DIS kinematics monitor
//---------------------------------
Expand Down Expand Up @@ -406,7 +435,7 @@ inList.each { inFile ->
} // eo loop over objects in the file (run)


// fit asymmetry
// fit beam spin asymmetry
T.exeLeaves(monTree[runnum]['helic']['asym'],{
def particle = T.leafPath[0]
def grP = T.getLeaf(monTree,[runnum,'helic','sinPhi',particle,'hp','asymGrid'])
Expand Down Expand Up @@ -457,6 +486,7 @@ T.exeLeaves(monTree,{
if(tlPath.contains('sinPhi')) tlT = "sinPhiH"
else if(T.key=='heldefDist') tlT = "defined helicity fraction"
else if(T.key=='rellumDist') tlT = "n+/n-"
else if(tlPath.contains('beamChargeAsym')) tlT = "beam charge asymmetry"
else if(T.key=='asymGraph') tlT = "beam spin asymmetry: pion sin(phiH) amplitude"
else tlT = "unknown"
}
Expand All @@ -478,7 +508,7 @@ T.exeLeaves(monTree,{

// we also want a few timelines to monitor standard deviations
T.addLeaf(timelineTree,tlPath+'timelineDev',{
if(tlPath.contains('DIS') || tlPath.contains('inclusive') || tlPath.contains('nonMonotonicity')) {
if(tlPath.contains('DIS') || tlPath.contains('inclusive') || T.key=='valGraph') {
def tlN = (tlPath+'timelineDev').join('_')
def tlT
if(tlPath.contains('DIS')) tlT = "DIS kinematics"
Expand Down Expand Up @@ -517,9 +547,7 @@ T.exeLeaves(monTree,{
def devs = vals.collect{ (it-ave)**2 }
def stddev = tot>0 ? Math.sqrt( devs.sum() / tot ) : tot
T.getLeaf(timelineTree,tlPath+'timeline').addPoint(tlRun,ave,0.0,0.0)
if(tlPath.contains('DIS') || tlPath.contains('inclusive') || tlPath.contains('nonMonotonicity')) {
T.getLeaf(timelineTree,tlPath+'timelineDev').addPoint(tlRun,stddev,0.0,0.0)
}
T.getLeaf(timelineTree,tlPath+'timelineDev').addPoint(tlRun,stddev,0.0,0.0)
}
// or if it's a helicity distribution monitor, add the run's overall fractions
if(T.key=='heldefDist' || T.key=='rellumDist') {
Expand All @@ -529,20 +557,34 @@ T.exeLeaves(monTree,{
def frac = denom>0 ? numer/denom : 0
T.getLeaf(timelineTree,tlPath+'timeline').addPoint(tlRun,frac,0.0,0.0)
}
// or if it's an asymmetry graph, add fit results to the timeline
// or if it's an asymmetry graph, add its results to the timeline
if(T.key=='asymGraph') {
def valPath = T.leafPath[0..-2] + 'asymValue'
def errPath = T.leafPath[0..-2] + 'asymError'
def asymVal = T.getLeaf(monTree,valPath)
def asymErr = T.getLeaf(monTree,errPath)
T.getLeaf(timelineTree,tlPath+'timeline').addPoint(tlRun, asymVal, 0.0, asymErr)
// and assign a defect bit for pi+ BSA
if(tlPath.contains('pip')) {
def asymMargin = asymVal.abs() - asymErr
if(asymMargin <= 0) {
addDefectBit(T.bit("BSAUnknown"), tlRun, allBins, allSectors)
} else if(asymVal < 0) {
addDefectBit(T.bit("BSAWrong"), tlRun, allBins, allSectors)
// beam charge asymmetry --------
if(T.leafPath.contains("beamChargeAsym")) {
def numHel = ['numHelP','numHelM'].collect{T.getLeaf(monTree, T.leafPath[0..-2] + it)}
def asym = calculateBeamChargeAsym(*numHel)
if(!asym.contains("unknown")) {
T.getLeaf(timelineTree,tlPath+'timeline').addPoint(tlRun, asym[0], 0.0, asym[1])
}
else {
System.err.println "WARNING: unknown beam charge asymmetry for run $tlRun"
}
}
// beam spin asymmetry --------
else {
def valPath = T.leafPath[0..-2] + 'asymValue'
def errPath = T.leafPath[0..-2] + 'asymError'
def asymVal = T.getLeaf(monTree,valPath)
def asymErr = T.getLeaf(monTree,errPath)
T.getLeaf(timelineTree,tlPath+'timeline').addPoint(tlRun, asymVal, 0.0, asymErr)
// and assign a defect bit for pi+ BSA
if(tlPath.contains('pip')) {
def asymMargin = asymVal.abs() - asymErr
if(asymMargin <= 0) {
addDefectBit(T.bit("BSAUnknown"), tlRun, allBins, allSectors)
} else if(asymVal < 0) {
addDefectBit(T.bit("BSAWrong"), tlRun, allBins, allSectors)
}
}
}
}
Expand All @@ -565,15 +607,16 @@ def hipoWrite = { hipoName, filterList, TLkeys ->
// will be renamed such that the front end plots them together
T.exeLeaves(tree,{
if(checkFilter(T.leafPath,filterList,T.key)) {
if(T.key=='asymValue' || T.key=='asymError' || T.key=='asymGrid') return
def name = T.leaf.getName()
if(name.contains('_hp_')) name = name.replaceAll('_hp_','_')
else if(name.contains('_hm_')) {
name = name.replaceAll('_hm_','_')
name += ":hm"
if(T.leaf instanceof IDataSet) {
def name = T.leaf.getName()
if(name.contains('_hp_')) name = name.replaceAll('_hp_','_')
else if(name.contains('_hm_')) {
name = name.replaceAll('_hm_','_')
name += ":hm"
}
T.leaf.setName(name)
outHipo.addDataSet(T.leaf)
}
T.leaf.setName(name)
outHipo.addDataSet(T.leaf)
}
})
}
Expand All @@ -595,6 +638,7 @@ def hipoWrite = { hipoName, filterList, TLkeys ->
hipoWrite("helicity_sinPhi",['helic','sinPhi'],["timeline"])
hipoWrite("beam_spin_asymmetry",['helic','asym'],["timeline"])
hipoWrite("defined_helicity_fraction",['helic','dist','heldef'],["timeline"])
hipoWrite("beam_charge_asymmetry",['helic','beamChargeAsym'],["timeline"])
hipoWrite("relative_yield",['helic','dist','rellum'],["timeline"])
hipoWrite("q2_W_x_y_means",['DIS'],["timeline"])
hipoWrite("pip_kinematics_means",['inclusive','pip'],["timeline"])
Expand Down
30 changes: 18 additions & 12 deletions qa-physics/monitorRead.groovy
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ def NBINS = 50 // number of bins in some histograms
def SECTORS = 0..<6 // sector range
def ECAL_ID = DetectorType.ECAL.getDetectorId() // ECAL detector ID
// debugging settings
def VERBOSE = true // enable extra log messages, for debugging
def VERBOSE = false // enable extra log messages, for debugging
def LIMITER = 0 // if nonzero, only analyze this many DST files (for quick testing and debugging)
def AUXFILE = false // enable auxfile production, an event-by-event table (a large text file)

Expand Down Expand Up @@ -77,7 +77,7 @@ else {

// limiter: use this to only analyse a few DST files, for quicker testing
if(LIMITER>0) {
inHipoList = inHipoList[0..LIMITER]
inHipoList = inHipoList[0..(LIMITER-1)]
System.err.println("WARNING WARNING WARNING: LIMITER ENABLED, we will only be analyzing ${LIMITER} DST files, and not all of them; this is for testing only!")
}

Expand Down Expand Up @@ -561,19 +561,15 @@ def overallMaxTimestamp = "init"
printDebug "Initialize histograms for each time bin"
timeBins.each{ binNum, timeBin ->
def partList = [ 'pip', 'pim' ]
def helList = [ 'hp', 'hm' ]
def heluList = [ 'hp', 'hm', 'hu' ]

T.buildTree(timeBin.histTree, 'helic', [['sinPhi'],partList,helList], { new H1F() })
T.buildTree(timeBin.histTree, 'helic', [['sinPhi'],partList,['hp','hm']], { new H1F() })
T.buildTree(timeBin.histTree, 'helic', [['dist']], { new H1F() })
T.buildTree(timeBin.histTree, 'helic', [['distGoodOnly']], { new H1F() })
T.buildTree(timeBin.histTree, 'helic', [['scaler'],['chargeWeighted']], { new H1F() })
T.buildTree(timeBin.histTree, 'DIS', [['Q2','W','x','y']], { new H1F() })
T.buildTree(timeBin.histTree, "DIS", [['Q2VsW']], { new H2F() })
T.buildTree(timeBin.histTree, "inclusive", [partList,['p','pT','z','theta','phiH']], { new H1F() })
// if(binNum==0) T.printTree(timeBin.histTree,{T.leaf.getClass()});

timeBin.histTree.helic.dist = buildHist('helic_dist','helicity',[],runnum,3,-1,2)
timeBin.histTree.helic.distGoodOnly = buildHist('helic_distGoodOnly','helicity (with electron cuts)',[],runnum,3,-1,2)
timeBin.histTree.DIS.Q2 = buildHist('DIS_Q2','Q^2',[],runnum,2*NBINS,0,12)
timeBin.histTree.DIS.W = buildHist('DIS_W','W',[],runnum,2*NBINS,0,6)
timeBin.histTree.DIS.x = buildHist('DIS_x','x',[],runnum,2*NBINS,0,1)
Expand All @@ -582,6 +578,7 @@ timeBins.each{ binNum, timeBin ->
T.exeLeaves( timeBin.histTree.helic.sinPhi, {
T.leaf = buildHist('helic_sinPhi','sinPhiH',T.leafPath,runnum,NBINS,-1,1)
})
timeBin.histTree.helic.scaler.chargeWeighted = buildHist('helic_scaler_chargeWeighted','FC-charge-weighted helicity',[],runnum,3,-1,2)
T.exeLeaves( timeBin.histTree.inclusive, {
def lbound=0
def ubound=0
Expand Down Expand Up @@ -635,6 +632,7 @@ inHipoList.each { inHipoFile ->
FTparticleBank = hipoEvent.getBank("RECFT::Particle")
calBank = hipoEvent.getBank("REC::Calorimeter")
scalerBank = hipoEvent.getBank("RUN::scaler")
helScalerBank = hipoEvent.getBank("HEL::scaler")

// get event number
def eventNum
Expand Down Expand Up @@ -739,7 +737,10 @@ inHipoList.each { inHipoFile ->
}

// get helicity and fill helicity distribution
def helicity = hipoEvent.hasBank("REC::Event") ? eventBank.getByte('helicity',0) : 0 // (using "0" for undefined)
def helicity = 0 // if undefined, default to 0
if(hipoEvent.hasBank("REC::Event") && eventBank.rows()>0) {
helicity = eventBank.getByte('helicity',0)
}
def helStr
def helDefined
switch(helicity) {
Expand All @@ -748,6 +749,14 @@ inHipoList.each { inHipoFile ->
default: helDefined = false; helicity = 0; break
}
thisTimeBin.histTree.helic.dist.fill(helicity)
// get scaler helicity from `HEL::scaler`, and fill its charge-weighted distribution
if(hipoEvent.hasBank("HEL::scaler")) {
helScalerBank.rows().times{ row -> // HEL::scaler readouts "pile up", so there are multiple bank rows in an event
def sc_helicity = helScalerBank.getByte("helicity", row)
def sc_fc = helScalerBank.getFloat("fcupgated", row) // helicity-latched FC charge
thisTimeBin.histTree.helic.scaler.chargeWeighted.fill(sc_helicity, sc_fc)
}
}

// get electron list, and increment the number of trigger electrons
// - also finds the DIS electron, and calculates x,Q2,W,y,nu
Expand All @@ -756,9 +765,6 @@ inHipoList.each { inHipoFile ->
// CUT: if a DIS electron was found by `findParticles`
if(disEleFound) {

if(disElectronInTrigger)
thisTimeBin.histTree.helic.distGoodOnly.fill(helicity)

// CUT for pions: Q2 and W and y and helicity
if( Q2>1 && W>2 && y<0.8 && helDefined) {

Expand Down

0 comments on commit 5a776ee

Please sign in to comment.