BremFitinForward

Qualification Requirement

Jira ticket: https://its.cern.ch/jira/browse/ATLIDTRKCP-119

The qualification task for Kaili Zhang, to be carried out in cooperation between ID Tracking and egamma, is defined as follows:

• Implement a calo ROI selection tuned for the forward region > 2.5 , to complement the existing calo seeding for <2.5
• This should then picked up by the pattern recognition to allow bremsstrahlung recovery fits to be performed
• Study the effect of using the above tracks and clusters > 2.5 to seed the Gaussian Sum Fitter
• Quantify the performance in this new region and possible gains for physics.
• Document the selection, code updates, and results

Background Knowledge

Particles, usually electrons, would have bremsstrahlung radiation in the inner detector. With energu loss, paricles passed the track detector and deposite in the calorimeter. Brem recovery method would use the deposition information from the calorimeter cluster and the track information from the inner tracker detector, to refit the particle.

In ATLAS Run-2:

Existing ATLAS calorimeter system covers the pseudorapidity range |η| < 4.9. Sampling calorimeters, based on Liquid Argon (LAr) technology are used for the detection of EM objects such as electrons and photons up to |η| = 3.2. They are LAr electromagnetic barrel(EMB), for |η|<1.475, and The EM End-cap(EMEC), for 1.375<|η|<3.2. The EMEC has 2 co-axial wheels, with the outer wheel (OW) covering the region 1.375 <|η|< 2.5 and the inner wheel (IW) covering the region 2.5< |η|< 3.2. Then we have forward calorimeters (FCal), covering 3.1<|η|<4.9, has 3 layers and FCal1 use copper, and FCal2 and FCal3 use tungsten. With EMB+EMEC+Fcal1, EM calorimeters caver the full range |η| < 4.9.

Existing ATLAS Inner Detector（ID）covers |η|< 2.5. And, in this period, the basic calorimeter reconstruction algorithm is sliding-window. https://twiki.cern.ch/twiki/bin/view/AtlasComputing/SlidingWindowClustering this algorithm only support |η|< 2.5.

The Phase-2 Upgrade, one major change is that ID was replaced by the Inner Tracker(ITK). Besides all silicon design, ITK provides the η coverage to 4.0. Also Calorimeters are upgraded to upgraded small-gap Forward Calorimeter (sFCal), but the coverage stay the same.

Since now ITK covers wider So Kaili's task is to implement the calo-seeded brem fit in the forward region.

Workflow

As it's a joint task for ID tracking and egamma, brem fit has been done for 2 times, in the ITK part and also egamma part. Since 2 parts have all build their own paricle objects. Tools are ready now, and the suitable inputs they need, are the track, the cluster, and the corresponding matching method to match track and cluster together.

All the work here are under athena r21.9.0 for ITK upgrade. Newer r21.9 have different detector layout like Step 3.0 which could have different configuration. And, the version for egamma part in r21.9 is an old version about ~2 years ago. I may consider to migrate the code to master in the following work. The following work didn't totally follow the instrcution from the egamma experts.

For the samples, first we mainly focus on the single particle gun electron sample. Muons do not have large brem loss. For convienence pt in each sample are the same. Here, ITK and Egamma have different electron pt perference. ITK use particle gun electron and the pt could be 5, 10, 20, 50, and 100GeV. While Egamma part enjoy Z->ee sample as benchmark. There is no universal way but currently we use pt=100GeV single electron sample.

in one specific detector version, we generate the corresponding EVNT/HITS files. Then the first reconstruction step generate RDO(Raw Data Object) files. HITS and RDO files do not change in this work. So we can keep the RDO files, try the brem fit and test the performance of the output DAOD_IDTRKVALID files by tools like IDPVM.

The workflow is that，

For ITK, the final trackParticleContainer name is called InDetTrackParticles. The default fit tool is GlobalChi2Fitter.

First, In InnerDetector /InDetExample/InDetRecExample/share/InDetRec_jobOptions.py, the general reconstruction process calls the package InDetCaloClusterROISelector. Here two options would affect the brem fit:

InDetFlags.doBremRecovery.set_Value_and_Lock(True); 
InDetFlags.doCaloSeededBrem.set_Value_and_Lock(True);

Default value for these 2 flags are True. (True, True) means we do brem fit and do calorimeter seeded brem fit. For the comparison, (False, False) do not do the brem fit, and (True False) would force all the tracks to do brem fit and do not consider the calorimeter conditions. The force brem fit would affect ~2% events of total electrons. which could be regarded as the upper limit of the calo seeded brem fit.

ROI, means the region of interest, ROISelector would loop all the candidated calorimeter clusters and do the selection. if passed, then it would call package InnerDetector \InDetRecTools\InDetCaloClusterROIBuilder to build the ROI object, containing the local position information of the cluster.

Then, in the ITK track maker InnerDetector /InDetRecTools/SiTrackMakerTool_xk, it loads the ROI information to match the tracks and clusters together.

Finally in the Tracking/TrkTools/TrkAmbiguityProcessor/src/DenseEnvironmentsAmbiguityProcessorTool.cxx, it call the selection, scoring, fitting tools and apply the brem fit. Also it provides a specific output about the fit details.

After the DenseEnvironmentsAmbiguityProcessorTool, the fit in ITK part finished. Now we have InDetTrackParticles container. Default r21.9 already have this track container eta range to 4.0, but the forward region in the output, has no events with brem fit.

For Egamma the container name is called GSFTrackParticles. GSF is GaussianSumFitter, the brem fit algorithm name.

Egamma used GSF fit in package Reconstruction\egamma\egammaTools\EMBremCollectionBuilder. Note that the directory changed in master and r21.9. This package used InDetTrackParticles as input tracker, and called package EMExtrapolationTools to determine whether the cluster and tracker is matched or not. If matched, it called the GSF tools and generate the output GSFTrackParticles track container. Default GSF container only have eta range up to 2.4.

In the tools IDPVM, it could use either InDetTrackParticles or GSFTrackParticles to see the performance, so at this step we could say brem fit in already finished. While from the egamma part, experts have further suggestions, like accompolish the function into the forward electron object. This part is not finished yet.

Modification

From the workflow part we listed the correlated packages about the brem fit both in ITK and Egamma part. So in the following I would list what I have done to make brem fit work in the forward region.

First, the default CaloCluster container is LArClusterEM. As it said, it used the sliding-window algorithm and the eta range is 2.5. It do not contain the forward information. So we need to switch to CaloCalTopoClusters, which means the calibrated topology clusters. The eta range is up to 4. Several validations are done to compare the performance of topoclusters and LArClusterEM, the deviation is accepctable. Actually, in the master release they plan to switch all the clusters to topoclusters.

The second issue is that, many existing egamma functions about calo clusters do not work in the forward region. Like to judge the cluster in in barrel or endcap it has isBarrel() and isEndcap(), for eta it has the second sampling function etaBE(2), and so on.

As eta() method still work in the forward, the simple way is to replace all the etaBE(2) to eta(). phiBE(2) to phi(). And to decide the cluster is forward or not, fabs(cluster->eta())>2.8 is used. This is not an elegant solution, but it work for now.

In the cluster selection part, in the central region the ratio like e237/e277 is used. But it do not work in the forward.

Besides, there are serveral places have the calosampling ID like EMB2 and EME2. In this case, forward CaloID, FCal0 is added.

In the Egamma part, the EMExtrapolationTools used the layer energy to match the cluster and track. But the forward topoclusters do not have layer energy in the beginning. So as suggested,

  CaloClusterKineHelper::calculateKine((xAOD::CaloCluster*)(*clus_iter),true,true);

is used to update the layer information for forward tracks. Currently it is called in the EMBremCollectionBuilder package. With this method, now the deltaeta and deltaphi can successfully return the value for matching. It should move to some other places for conviencence.

But compared to directly matching by cluster->eta() and track->eta(), it seems no difference at this time.

In the ROISelector and EMBremCollectionBuilder package, I add the total forward clusters in summary output. Which could be located with SelectedForwardClusters and SelectedForwardTracks.

The AmbiguityProcessor output is just below the summary of InDetCaloClusterROISelector. Or it can be located with “including number of brem fits”.

There are some selection creteria to select the forward topoclusters, since the topoclusters are not 1 cluster corresponding to 1 track. Usualy 8 or more clusters for one track. But for convienence, we didn't apply complex selection in this matching step. So currently only Et>5GeV are used to select the forward clusters. For pt100GeV electron sample, usually only 1 or 2 topocluster could satisfy Et>5GeV. So it's enough at present.

Latest git branch for merge is located in https://gitlab.cern.ch/kazhang/athena/tree/BremFit_ForMerge

Repeat my work

For the environment:

mkdir workdir && cd workdirmkdir run source build && cd source asetup 21.9.0,Athena,here,slc6git atlas init-workdir https://:@gitlab.cern.ch:8443/kazhang/athena.gitcd athenagit checkout BremFit_ForMergegit atlas addpkg InDetCaloClusterROIBuilder InDetCaloClusterROISelector egammaTools egammaTrackTools// egammaRec and xAODCaloEvent didn't count in current version.cd ../../buildcmake ../sourcecmake --build .//before you run the reco, source x86_64-centos7-gcc62-opt/setup.sh

For the sample used, now

/eos/user/k/kazhang/QT_Electron_EVNT
/eos/user/k/kazhang/QT_Electron_HITS
/eos/user/k/kazhang/QT_Electron_RDO

is public for all atlas users.

If you want to generate the samples from the beginning,

Sim_tf.py \--inputEVNTFile mc15_14TeV.422030.ParticleGun_single_ele_Pt20.evgen.EVNT.e5286/EVNT.09244564._000001.pool.root.1 \--outputHITSFile HITS_1.root \--skipEvents 0 \--maxEvents 1000000 \--randomSeed 1234 \--geometryVersion ATLAS-P2-ITK-20-03-00_VALIDATION \--conditionsTag OFLCOND-MC15c-SDR-14-03 \--truthStrategy MC15aPlus \--DataRunNumber 242000 \--preInclude all:'InDetSLHC_Example/preInclude.SLHC.py,InDetSLHC_Example/preInclude.NoTRT_NoBCM_NoDBM.py,InDetSLHC_Example/preInclude.SLHC_Setup_InclBrl_4.py,InDetSLHC_Example/preInclude.SLHC_Setup_Strip_GMX.py' \--preExec all:'from AthenaCommon.DetFlags import DetFlags;DetFlags.Muon_setOff();from InDetSLHC_Example.SLHC_JobProperties import SLHC_Flags; SLHC_Flags.doGMX.set_Value_and_Lock(True); SLHC_Flags.LayoutOption="InclinedDuals";' \--postInclude all:'PyJobTransforms/UseFrontier.py,InDetSLHC_Example/postInclude.SLHC_Setup_InclBrl_4.py,InDetSLHC_Example/postInclude.SLHC_Setup.py' \--postExec EVNTtoHITS:'from AthenaCommon.DetFlags import DetFlags;DetFlags.Muon_setOff();ServiceMgr.DetDescrCnvSvc.DoInitNeighbours=False; from AthenaCommon import CfgGetter; CfgGetter.getService("ISF_MC15aPlusTruthService").BeamPipeTruthStrategies+= ["ISF_MCTruthStrategyGroupIDHadInt_MC15"];'

which could generate HITS files.

To generate output DAOD_IDTRKVALIDFile， either use RDO or HITS files as input, and the reco command is

Reco_tf.py \
--digiSteeringConf StandardInTimeOnlyTruth \
--geometryVersion ATLAS-P2-ITK-20-03-00 \
--conditionsTag OFLCOND-MC15c-SDR-14-03 \
--DataRunNumber 242000 \
--steering doRAWtoALL \
--postInclude all:'InDetSLHC_Example/postInclude.SLHC_Setup_InclBrl_4.py' \
HITtoRDO:'InDetSLHC_Example/postInclude.SLHC_Digitization_lowthresh.py,InDetSLHC_Example/postInclude.RDOAnalysis.py' \
RAWtoALL:'InDetSLHC_Example/postInclude.DigitalClustering.py' \
--preExec all:'from AthenaCommon.GlobalFlags import globalflags; globalflags.DataSource.set_Value_and_Lock("geant4"); from InDetSLHC_Example.SLHC_JobProperties import SLHC_Flags; SLHC_Flags.doGMX.set_Value_and_Lock(True); SLHC_Flags.LayoutOption="InclinedDuals";' \
HITtoRDO:'from Digitization.DigitizationFlags import digitizationFlags; digitizationFlags.doInDetNoise.set_Value_and_Lock(False); digitizationFlags.doBichselSimulation.set_Value_and_Lock(False); digitizationFlags.overrideMetadata+=["SimLayout","PhysicsList"];' \
RAWtoALL:'from InDetRecExample.InDetJobProperties import InDetFlags; from PixelConditionsServices.PixelConditionsServicesConf import PixelCalibSvc;ServiceMgr +=PixelCalibSvc();InDetFlags.useDCS.set_Value_and_Lock(True);ServiceMgr.PixelCalibSvc.DisableDB=True;from InDetPrepRawDataToxAOD.InDetDxAODJobProperties import InDetDxAODFlags;InDetDxAODFlags.DumpLArCollisionTime.set_Value_and_Lock(False);InDetDxAODFlags.DumpSctInfo.set_Value_and_Lock(True);InDetDxAODFlags.ThinHitsOnTrack.set_Value_and_Lock(False); from InDetRecExample.InDetKeys import InDetKeys; InDetKeys.CaloClusterContainer.set_Value_and_Lock("CaloCalTopoClusters");' \
--preInclude all:'InDetSLHC_Example/preInclude.SLHC_Setup_InclBrl_4.py,InDetSLHC_Example/preInclude.SLHC_Setup_Strip_GMX.py,InDetSLHC_Example/preInclude.SLHC_Calorimeter_mu0.py' \
HITtoRDO:'InDetSLHC_Example/preInclude.SLHC.py,InDetSLHC_Example/preInclude.NoTRT_NoBCM_NoDBM.py' \
default:'InDetSLHC_Example/preInclude.SLHC.NoTRT_NoBCM_NoDBM.Reco.py,InDetSLHC_Example/SLHC_Setup_Reco_TrackingGeometry_GMX.py' \
RDOMergeAthenaMP:'InDetSLHC_Example/preInclude.SLHC.py,InDetSLHC_Example/preInclude.NoTRT_NoBCM_NoDBM.py' \
POOLMergeAthenaMPAOD0:'InDetSLHC_Example/preInclude.SLHC.NoTRT_NoBCM_NoDBM.Ana.py,InDetSLHC_Example/SLHC_Setup_Reco_Alpine.py' \
POOLMergeAthenaMPDAODIDTRKVALID0:'InDetSLHC_Example/preInclude.SLHC.NoTRT_NoBCM_NoDBM.Ana.py,InDetSLHC_Example/SLHC_Setup_Reco_Alpine.py' \
--postExec HITtoRDO:'pixeldigi.EnableSpecialPixels=False; CfgMgr.MessageSvc().setError+=["HepMcParticleLink"];' \
RAWtoALL:'ToolSvc.InDetSCT_ClusteringTool.useRowInformation=True; from AthenaCommon.AppMgr import ToolSvc; ToolSvc.InDetTrackSummaryTool.OutputLevel=INFO;from InDetPhysValMonitoring.InDetPhysValMonitoringConf import InDetPhysValDecoratorAlg;decorators = InDetPhysValDecoratorAlg();topSequence += decorators' \
--inputRDOFile=/eos/user/k/kazhang/QT_Electron_RDO/pt100/RDO_${index}.root \
--outputDAOD_IDTRKVALIDFile=/eos/user/k/kazhang/QT_Electron_DAOD/DAOD_IDTRKVALID_${index}.root \
--maxEvents '10000'

change the input and output directory as you need. Notice that in the --preExec part we have one addtional statement to ask ROISelector use CaloCalTopoClusters rather than LArClusterEM.

from InDetRecExample.InDetKeys import InDetKeys; InDetKeys.CaloClusterContainer.set_Value_and_Lock("CaloCalTopoClusters");

For some other places, default values are also changed in the python script or somewhere else.

After the command it should have one DAOD and one log.RAWtoALL file. In log files we could pick the AmbiguityProcessor output and the summary of the InDetCaloClusterROISelector and EMBremCollectionBuilder.

In the ROISelector and EMBremCollectionBuilder package, I add the total forward clusters in summary output. Which could be located with SelectedForwardClusters and SelectedForwardTracks.

The AmbiguityProcessor output is just below the summary of InDetCaloClusterROISelector. Or it can be located with “including number of brem fits”.

Current result for pt100GeV electron is listed below. In the toal 45k events ~2k events are brem fitted while in the forward region there 101 events affected. Affected electrons would be about 0.5% of the total forward tracks.

|------------------------------------------------------------------------------------
|  Number of events processed      :  45000
|  statistics by eta range          ------All---Barrel---Trans.-- Endcap-- Forwrd--
|------------------------------------------------------------------------------------
|  Number of candidates at input   :     94720  9501    11288   29749   44173
|  - candidates rejected score 0   :     2015   1   21  1719    265
|  - candidates rejected as double :        0        0        0        0        0
|------------------------------------------------------------------------------------
|  candidates with good score      :   96146    9022    10184   28042   48898
|  + recovered after brem refit    :       37        1        5       28        3
|  candidates rejected score 0     :    3015    11  71  1047    1886
|  + rejected failed brem refit    :        28  0   1   19  8
|  + rejected brem refit score 0   :     7525   101 495 5336    1593
|------------------------------------------------------------------------------------
|  number of normal fits           :    78332   8815    9727    23774   36016
|  + 2nd brem fit for failed fit   :      292   6   40  226 20
|  normal brem fits for electrons  :    2917    251 580 1869    217
|------------------------------------------------------------------------------------
|  sum of succesful fits           :   88880    9169    10826   31061   37824
|  sum of failed fits              :      212   3   23  153 33
|------------------------------------------------------------------------------------
|  Number of subtracks created     :     25939  78  514 9145    16202
|  Number of candidates excluded   :    29302   676 1701    7321    19604
|------------------------------------------------------------------------------------
|  Number of tracks accepted       :    52651   8708    8977    14176   20790
|  including number of brem fits   :     1975   219 458 1197    101
|------------------------------------------------------------------------------------
|    definition: ( 0.0 < Barrel < 0.80 < Transition < 1.60 < Endcap < 2.50 < Forward < 4.00 )
|------------------------------------------------------------------------------------

Compared to the force brem fit (True, False), 45k events would have 295 forward brem fitted events. Though in the central region the affected events could be also smaller, but in the forward region it reduced even more. It could be bad scoring or some thing like that. Could be dig further.

Then we want to see the performance of the output DAOD. Get one copy of InnerDetector /InDetValidation/InDetPhysValMonitoring/run/PhysValITk_jobOptions.py

Change the FNAME as you need. And change the following configuration

runDAOD=True
TrackTruthSelectionTool.maxEta = 4.2

Then, run

athena PhysValITk _jobOptions.py

to get the IDPVM output.

The default command for ParticleContainer is that

InDetPhysValMonitoringTool.TruthParticleContainerName = "TruthParticles"
InDetPhysValMonitoringTool.TrackParticleContainerName = "

InDetTrackParticles

"
//So we can change the track container to GSF:
InDetPhysValMonitoringTool.TrackParticleContainerName= "GSFTrackParticles"

Now we have 2 output for both ITK and Egamma.This process can be repeat successfully in Sept. 24, 2019.

IDPVM Performance

Several plots from IDPVM are interesting.

This plot shows the track reconstruction effciency versus eta. First we can see now both ITK and Egamma has full eta range up to 4.0. While the default GSF particle only work in the central region. Now we see it works in the forward region.

Secondly, both ITK and Egamma has high effciency in the forward region. while at the second bin GSF has lower effciency corrsponding to ITK. Specificly, in the IDPVM output eff_vs_mu2 (2.7<|eta|<3.5) is 94.9% for GSF, and for ITK it's 98.5%.

The reason for this drop is still unknown. Using layer energy to match the cluster and track helps a little.

The Eta_reco plot also shows that in most regions GSF and ITK share the similar result.

One another plot is that the 2d pt vs eta plot.

As all the events with truth pt 100GeV, the reco pt vs reco eta distrubted like this. First thing is that though total number similar, GSF obvious has larger mean pt value 71 over 63, and more close to truth value 100. It could prove that GSF fit surely improves the final track pt and make it closer to the truth value.

Secondly, both ITK and GSF there are large amount of forward events,have reconstructed with a small pt. Most likely it is due to that track spilt to 2 smaller track, or 2 pieces of topoclusters generated. Current algorithm works with one track mathing one cluster. Here we should note that brem fit do not help to improve this situation since brem fit in the ITK part only affect ~0.5% forward events. While all the events in the GSF are brem fitted.

Other useful performance benchmarks usually take the following combination.

These plots are packed in Performance.zip.

Another plot to see the pt difference can found in the ptmeanRMS_vs_eta. Both 2 tracks suffered in the forward region while GSF is a bit better.

Also we can see qoverppullresolutionRMS_vs_eta to see better performance GSF over than ITK. (Small pull width the better.)

Major updates:
-- KailiZhang - 2019-09-24