-- FrancescoMariaFollega - 2018-01-30

Introduction

This page is based on my experience on running Zmumu samples with Athena 21.0.20. For setting up release 21, please check the following twiki:

• https://twiki.cern.ch/twiki/bin/viewauth/AtlasComputing/AlignmentGitAthenaRepository

The Zmumu studies are based on the scripts found in /InnerDetector/InDetMonitoring/InDetPerformanceMonitoring/scripts folder.

This work has been carried on in the context of the Qualification Task (https://its.cern.ch/jira/browse/ATLIDTRKCP-104). I produced new maps for the detection of weak modes and I derived a new way for calculating delta sagitta correction using MC samples as a reference. In the following I will explain the script I used in such a way that my studies can be reproduced.

Producing N-tuples:panda/athena

The pre-requisite of this analysis is to have suitable ntuples. For the production and the merging of these samples you can exploit the tools shown in the following tutorial https://twiki.cern.ch/twiki/bin/view/Sandbox/HowtorunZmumumaps.

Analyzing the n-tuple and producing the root files: ZmumuValidationExample

The analysis boils down to three scripts: ZmumuValidationExample.cxx (which contains the methods and the loops on top of the ntuples), run_ZmumuValidationExample.C (which basicaly runs the code)and one script for the setting of the input files. A previous documentation can be found here (https://twiki.cern.ch/twiki/bin/view/Sandbox/HowtorunZmumumaps).

This code is designed to calculate iteratively the correction to the weak modes. At the end of this iteration it will produce maps for the detection and correction of weak modes (eg. sagitta distortion and radial distortion).

Eg. in the sagitta distortion analysis it calculates the correction factor for the pT, then it correct the pTs for each event and finally performs an other loop calculating the residual distortion/correction.

The script run_ZmumuValidationExample.C has different modes that can be set upt:

run_ZmumuValidationExample(Int_t nIterationsUserInput, Int_t sampletype, Int_t analysismode)

1. nIterationUser: the number of iteration expeted to complete the analysis, usually 20 iteration are sufficient;
2. sampleType: the sample on which you are running the code -> {0 = ZSAMPLE, 1 = JPSISAMPLE, 2 = UPSILONSAMPLE, 3 = DIMUONSAMPLE, 4 = KSHORTSAMPLE};
3. analysismode: the type of analysis which you want to carry on -> {ANALYSIS_ALL = 0, ANALYSIS_RADIALDIST = 1, ANALYSIS_SAGITTADIST= 2}.

You can run the code with three steps:

root -l
.L run_ZmumuValidationExample.C++
run_ZmumuValidationExample(20,0,2) #this will run 20 iteration in ANALYSIS_SAGITTADIST alone, on Zmumu samples.

It will give you an output (you can set the output name in the run_ZmumuValidationExample.C) with all the maps and histograms iteration-by-iteration.

Updates on the ZmumuValidationExample package: Part 1.

In this section I will show updates and modifies on the preaviouly listed scripts.

In order to include the markers of the previous section to monitor the weak modes I added new histograms and new variables and methods to fill it. Here below there is a list of them added in the method ZmumuValidationExample::bookHistograms():

// histogram for pT leading study: each eta-phi bin is filled with the corrisponding quantity in the direction of the pT leading muon.This histogram are filled iteration-by-iteration;

 h_3Dhisto_dimuon_etaphi_ptlead = new TH3D ("h_3Dhisto_dimuon_etaphi_ptlead"," 3D histogram of mass_{reco}; leading p_{T} #eta; leading p_{T} #phi; M_{reco}[GeV]", fnEtaBins, -eta_bound, eta_bound, fnPhiBins, -phi_bound, phi_bound,80,70,110);
 h_dimuon_DeltaM_vs_etaphi_iteration = new TH2F ("h_dimuon_DeltaM_vs_etaphi_iteration","(m_{reco} - M_{PDG})/M_{PDG} vs etaphi; leading p_{T} #eta; leading p_{T} #phi", fnEtaBins, -eta_bound, eta_bound, fnPhiBins, -phi_bound, phi_bound);
 h_dimuon_DeltaWidth_vs_etaphi_iteration = new TH2F ("h_dimuon_DeltaWidth_vs_etaphi_iteration","(Width_{reco} - Width_{PDG})/Width_{PDG} vs etaphi; leading p_{T} #eta; leading p_{T} #phi", fnEtaBins-1, -eta_bound, eta_bound, fnPhiBins-1, -phi_bound, phi_bound);
 h_dimuon_DeltaWidth_vs_etaphi_iteration_err = new TH2F ("h_dimuon_DeltaWidth_vs_etaphi_iteration_err","(Width_{reco} - Width_{PDG})/Width_{PDG} vs etaphi; leading p_{T} #eta; leading p_{T} #phi", fnEtaBins-1, -eta_bound, eta_bound, fnPhiBins-1, -phi_bound, phi_bound);

// histogram for pos and neg muons study: each eta-phi bin is filled with the corrisponding quantity in the direction of the pos or neg muon. This histogram are filled iteration-by-iteration;

 h_muonpos_DeltaWidth_vs_etaphi_iteration = new TH2F ("h_muonpos_DeltaWidth_vs_etaphi_iteration","(Width_{reco} - Width_{PDG})/Width_{PDG} vs etaphi; #mu^+ #eta; #mu^+ #phi", fnEtaBins-1, -eta_bound, eta_bound, fnPhiBins-1, -phi_bound, phi_bound);
 h_muonneg_DeltaWidth_vs_etaphi_iteration = new TH2F ("h_muonneg_DeltaWidth_vs_etaphi_iteration","(Width_{reco} - Width_{PDG})/Width_{PDG} vs etaphi; #mu^- #eta; #mu^- #phi", fnEtaBins-1, -eta_bound, eta_bound, fnPhiBins-1, -phi_bound, phi_bound);
 h_muonpos_DeltaM_vs_etaphi_iteration = new TH2F ("h_muonpos_DeltaM_vs_etaphi_iteration","(m_{reco} - M_{PDG})/M_{PDG} vs etaphi; #mu^+ #eta; #mu^+ #phi", fnEtaBins, -eta_bound, eta_bound, fnPhiBins, -phi_bound, phi_bound);
 h_muonneg_DeltaM_vs_etaphi_iteration = new TH2F ("h_muonneg_DeltaM_vs_etaphi_iteration","(m_{reco} - M_{PDG})/M_{PDG} vs etaphi; #mu^+ #eta; #mu^+ #phi", fnEtaBins, -eta_bound, eta_bound, fnPhiBins, -phi_bound, phi_bound);

// histogram for MC subtraction. In this part are declared histograms in which, at the last iteration is perfomed the MC_mass or PDG_mass subratction. They have been designed in order to evaluate marker1 and marker2;

h_dimuon_DeltaM_vs_etaphi = new [[TH2F][TH2F]] ("h_dimuon_DeltaM_vs_etaphi","(m_{reco} - M_{PDG})/M_{PDG} vs etaphi; leading p_{T} #eta; leading p_{T} #phi",fnEtaBins, -eta_bound, eta_bound, fnPhiBins, -phi_bound, phi_bound);
h_dimuon_DeltaWidth_vs_etaphi = new [[TH2F][TH2F]]("h_dimuon_DeltaWidth_vs_etaphi","(Width_{reco} - Width_{PDG})/Width_{PDG} vs etaphi; leading p_{T} #eta; leading p_{T} #phi", fnEtaBins-1, -eta_bound, eta_bound, fnPhiBins-1, -phi_bound, phi_bound);
h_muonpos_DeltaWidth_vs_etaphi = new [[TH2F][TH2F]] ("h_muonpos_DeltaWidth_vs_etaphi","(Width_{reco} - Width_{PDG})/Width_{PDG} vs etaphi; #mu^+ #eta; #mu^+ #phi", fnEtaBins-1, -eta_bound, eta_bound, fnPhiBins-1, -phi_bound, phi_bound);
h_muonneg_DeltaWidth_vs_etaphi = new [[TH2F][TH2F]] ("h_muonneg_DeltaWidth_vs_etaphi","(Width_{reco} - Width_{PDG})/Width_{PDG} vs etaphi; #mu^- #eta; #mu^- #phi", fnEtaBins-1, -eta_bound, eta_bound, fnPhiBins-1, -phi_bound, phi_bound);
h_muonpos_DeltaM_vs_etaphi = new [[TH2F][TH2F]]("h_muonpos_DeltaM_vs_etaphi","(m_{reco} - M_{PDG})/M_{PDG} vs etaphi; #mu^+ #eta; #mu^+ #phi", fnEtaBins, -eta_bound, eta_bound, fnPhiBins, -phi_bound, phi_bound);
h_muonneg_DeltaM_vs_etaphi = new [[TH2F][TH2F]] ("h_muonneg_DeltaM_vs_etaphi","(m_{reco} - M_{PDG})/M_{PDG} vs etaphi; #mu^+ #eta; #mu^+ #phi", fnEtaBins, -eta_bound, eta_bound, fnPhiBins, -phi_bound, phi_bound);

//this histograms are for MC maps (that had to be created previously, running the script over a unbiased MC sample), they will be useful for the MC_mass and MC_width subtraction;

h_MC_mass_Zmumu_etaphi = new TH2F ("h_MC_mass_Zmumu_etaphi","Mass_MC vs etaphi; leading p_{T} #eta; leading p_{T} #phi", fnEtaBins, -eta_bound, eta_bound, fnPhiBins, -phi_bound, phi_bound);
 h_MC_width_Zmumu_etaphi = new TH2F ("h_MC_width_Zmumu_etaphi","Mass_MC vs etaphi; leading p_{T} #eta; leading p_{T} #phi", fnEtaBins, -eta_bound, eta_bound, fnPhiBins, -phi_bound, phi_bound);
 h_MC_mass_Zmumu_muonpos_etaphi = new TH2F ("h_MC_mass_Zmumu_muonpos_etaphi","Mass_MC vs etaphi , #mu^{+}; #mu^{+} #eta; #mu^{+} #phi", fnEtaBins, -eta_bound, eta_bound, fnPhiBins, -phi_bound, phi_bound);
 h_MC_mass_Zmumu_muonneg_etaphi = new TH2F ("h_MC_mass_Zmumu_muonneg_etaphi","Mass_MC vs etaphi , #mu^{-}; #mu^{-} #eta; #mu^{-} #phi",
 fnEtaBins, -eta_bound, eta_bound, fnPhiBins, -phi_bound, phi_bound);
 h_MC_width_Zmumu_muonpos_etaphi = new TH2F ("h_MC_width_Zmumu_muonpos_etaphi","Mass_MC vs etaphi , #mu^{+} ; #mu^{+} #eta; #mu^{+} #phi",
 fnEtaBins, -eta_bound, eta_bound, fnPhiBins, -phi_bound, phi_bound);
 h_MC_width_Zmumu_muonneg_etaphi = new TH2F ("h_MC_width_Zmumu_muonneg_etaphi","Mass_MC vs etaphi , #mu^{-}; #mu^{-} #eta; #mu^{-} #phi", fnEtaBins, -eta_bound, eta_bound, fnPhiBins, -phi_bound, phi_bound);

After the declaration of the histograms I fill the 3D ones in ZmumuValidationExample::fillHistograms and the 2D ones in the method ZmumuValidationExample::loop() using the ZmumuValidationExample::profileZwithIterativeGaussFit. The iterative gaussian fit is a method designed to perform a Z profile of a 3D histogram and then return 2D histogram rapresenting the mean, width and error on the mean. The code is the following:

//find ptlead and fill the 3D histogram

double eta_muonlead,phi_muonlead,pt_muonlead;
double eta_muonsublead,phi_muonsublead,pt_muonsublead; 
if(abs(vec_pos->Pt())>=abs(vec_neg->Pt())){ eta_muonlead = vec_pos->Eta();
phi_muonlead = vec_pos->Phi();
pt_muonlead = abs(vec_pos->Pt()/1000.);
eta_muonsublead = vec_neg->Eta();
phi_muonsublead = vec_neg->Phi();
pt_muonsublead = abs(vec_neg->Pt()/1000.); 
} else{
eta_muonlead = vec_neg->Eta(); phi_muonlead = vec_neg->Phi();
  pt_muonlead = abs(vec_neg->Pt()/1000.); eta_muonsublead = vec_pos->Eta();
 phi_muonsublead = vec_pos->Phi(); pt_muonsublead = abs(vec_pos->Pt()/1000.);
 }
h_3Dhisto_dimuon_etaphi_ptlead->Fill(eta_muonlead,phi_muonlead,vec_dimuon->M()/1000.);

//MC subtraction and normalization

for(int etabin = 1; etabin <= h_dimuon_DeltaM_vs_etaphi->GetNbinsX(); ++etabin){
 for(int phibin = 1; phibin <= h_dimuon_DeltaM_vs_etaphi->GetNbinsY(); ++phibin){
 h_dimuon_DeltaM_vs_etaphi->SetBinContent(etabin, phibin,(h_dimuon_DeltaM_vs_etaphi_iteration->GetBinContent(etabin,phibin) - h_MC_mass_Zmumu_etaphi->GetBinContent(etabin,phibin))/h_MC_mass_Zmumu_etaphi->GetBinContent(etabin,phibin));
 h_dimuon_DeltaWidth_vs_etaphi->SetBinContent(etabin, phibin,(h_dimuon_DeltaWidth_vs_etaphi_iteration->GetBinContent(etabin,phibin) - h_MC_width_Zmumu_etaphi->GetBinContent(etabin,phibin))/h_MC_width_Zmumu_etaphi->GetBinContent(etabin,phibin));
 h_muonpos_DeltaM_vs_etaphi->SetBinContent(etabin, phibin,(h_muonpos_DeltaM_vs_etaphi_iteration->GetBinContent(etabin,phibin) - h_MC_mass_Zmumu_muonpos_etaphi->GetBinContent(etabin,phibin))/h_MC_mass_Zmumu_muonpos_etaphi->GetBinContent(etabin,phibin));
 h_muonpos_DeltaWidth_vs_etaphi->SetBinContent(etabin, phibin,(h_muonpos_DeltaWidth_vs_etaphi_iteration->GetBinContent(etabin,phibin) - h_MC_width_Zmumu_muonpos_etaphi->GetBinContent(etabin,phibin))/h_MC_width_Zmumu_muonpos_etaphi->GetBinContent(etabin,phibin));
 h_muonneg_DeltaM_vs_etaphi->SetBinContent(etabin, phibin,(h_muonneg_DeltaM_vs_etaphi_iteration->GetBinContent(etabin,phibin) - h_MC_mass_Zmumu_muonneg_etaphi->GetBinContent(etabin,phibin))/h_MC_mass_Zmumu_muonneg_etaphi->GetBinContent(etabin,phibin));
 h_muonneg_DeltaWidth_vs_etaphi->SetBinContent(etabin, phibin,(h_muonneg_DeltaWidth_vs_etaphi_iteration->GetBinContent(etabin,phibin) - h_MC_width_Zmumu_muonneg_etaphi->GetBinContent(etabin,phibin))/h_MC_width_Zmumu_muonneg_etaphi->GetBinContent(etabin,phibin));
 }
}

[...]

//fill ptleading muon 2D histograms

profileZwithIterativeGaussFit(h_3Dhisto_dimuon_etaphi_ptlead,h_dimuon_DeltaM_vs_etaphi_iteration,h_dimuon_DeltaWidth_vs_etaphi_iteration, 1, h_dimuon_DeltaWidth_vs_etaphi_iteration_err);
profileZwithIterativeGaussFit(h_dimuon_mass_vs_etaphiAll,h_dimuon_DeltaM_vs_etaphi_iteration,h_dimuon_DeltaWidth_vs_etaphi_iteration, 1, h_dimuon_DeltaWidth_vs_etaphi_iteration_err);

//fill pos and neg muons 2D histograms

profileZwithIterativeGaussFit(h_dimuon_mass_vs_etaphiPos,h_muonpos_DeltaM_vs_etaphi_iteration,h_muonpos_DeltaWidth_vs_etaphi_iteration, 1, 0);
profileZwithIterativeGaussFit(h_dimuon_mass_vs_etaphiNeg,h_muonneg_DeltaM_vs_etaphi_iteration,h_muonneg_DeltaWidth_vs_etaphi_iteration, 1, 0);

These histograms allow you to quickly check if a misalignment is affecting the reconstruction capabilities of the dimuon mass. Here a post you some examples of the maps you will get (WARNING: graphics a bit digested and plotted according to ATLAS format) on 2016 data period D.

As you can notice there are three regions that show a misalignment that has been detected also in this public plots (https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/PLOTS/IDTR-2017-005/).

Updates on the ZmumuValidationExample package: Part 2.

Currently the delta sagitta correction is calculated in the method in the ZmumuValidationExample::fillEtaPhiHistogram() .This method is based on the following formula:

In the currently implemented code Mass_mumu0 is equal to the Zmass from PDG and it is costant over all the the eta-phi coordinates. I investigated an other choice for this reference mass: MC_mass calculated for each eta-phi coordinates.

In order to introduce this new method called ZmumuValidationExample::fillEtaPhiHistogram_usingMCmass() which uses the MC mass replacing the flat PDG one.

void [[ZmumuValidationExample][ZmumuValidationExample]]::fillEtaPhiHistogram_usingMCmass(TH3* hist, [[TLorentzVector][TLorentzVector]]* v_pos, [[TLorentzVector][TLorentzVector]]* v_neg, int use_lambda){
 
 int bin_x_pos = h_MC_mass_Zmumu_muonpos_etaphi->GetXaxis()->FindBin(v_pos->Eta());
 int bin_y_pos = h_MC_mass_Zmumu_muonpos_etaphi->GetYaxis()->FindBin(v_pos->Phi());
 int bin_x_neg = h_MC_mass_Zmumu_muonneg_etaphi->GetXaxis()->FindBin(v_neg->Eta());
 int bin_y_neg = h_MC_mass_Zmumu_muonneg_etaphi->GetYaxis()->FindBin(v_neg->Phi());
 double z_mass_pos = h_MC_mass_Zmumu_muonpos_etaphi->GetBinContent( bin_x_pos , bin_y_pos )*1000.; 
 double z_mass_neg = h_MC_mass_Zmumu_muonneg_etaphi->GetBinContent( bin_x_neg , bin_y_neg )*1000.;
double mass = ((*v_pos) + (*v_neg)).M();
 double delta_M2_pos = (mass*mass - z_mass_pos*z_mass_pos)/(z_mass_pos*z_mass_pos);
 double delta_M2_neg = (mass*mass - z_mass_neg*z_mass_neg)/(z_mass_neg*z_mass_neg);
if(delta_M2_pos > 100 || delta_M2_neg > 100){ 
 std::cout << mass << " " << z_mass_pos << std::endl;
 std::cout << mass << " " << z_mass_neg << std::endl;
 }
if (use_lambda) {
hist->Fill(v_pos->Eta(), v_pos->Phi(), +1*m_factor*delta_M2_pos/v_pos->Pt()*1000000.0);
 hist->Fill(v_neg->Eta(), v_neg->Phi(), -1*m_factor*delta_M2_neg/v_neg->Pt()*1000000.0);
 
 } else {
hist->Fill(v_pos->Eta(), v_pos->Phi(), +1*delta_M2_pos/2);
 hist->Fill(v_neg->Eta(), v_neg->Phi(), +1*delta_M2_neg/2);
 
 }
}

Closure test and convergence plot

After that it could be possible to perform a closure test: comparison between the two methods in order to enstablish which is the more accurete in the estimation of the sagitta bias. What you need is:

1. A perfectly aligned MC16 sample;
2. A biased MC16 sample with a known bias for each eta-phi bin;
3. A tool that calculates the closure between the input bias and the bias found after the loop over all the events;
4. A tool that evaluates the convergence speed and precision after a certain number of iteration.

Therefore I proposed two tools and I called them respectively "closure_test.cxx" and "iterationtrend_overlaied.cxx". The are tools that allows a graphical evaluation of the goodness of the method. I designed two markers ( q1 and q2 ) which help to estimate the closure:

.
To run these tools you need just the output file of run_ZmumuValidationExample.C macro (NB: this file must contain the 2D histogram of the input sagitta bias and the comuted bias correction iteration-by-iteration). If all is positioned in the right place they will return you;

1. either a map of the q1 marker for all the eta-phi bins, either an eta-phi map of the q2 marker (from closure_test.cxx);
2. a plot of the difference between the input bias and the estimated correction depending on the current iteration number(from iterationtrend_overlaied.cxx).

Here there are examples of the two plots.

There are also addictional tools:

1. sagitta_correlation.cxx, in order to estimate correlation between the deviation in the Z reconstucted mass in the direction of positive and negative muons;
2. statisticalerror.cxx, to estimate the statistical error of deviation in the Z reconstructed mass/sigma in the direction of positive and negative muons (bin by bin).

All these tools (as well as a copy of the updated Zmumu scripts of the previus sections) are committed on my branch on git lab with the following path :
(https://gitlab.cern.ch/atlas-idalignment/athena/tree/feature-21.0-ffollega-test/InnerDetector/InDetMonitoring/InDetPerformanceMonitoring/scripts/SagittaMapTools)