From 053c2710ac7d96afc0af2f77494c2bfc172fb37a Mon Sep 17 00:00:00 2001
From: Pavol Strizenec <pavol.strizenec@cern.ch>
Date: Thu, 2 Apr 2020 19:24:07 +0000
Subject: [PATCH] Adding code by A. Kiryunin

---
 .../Analysis/10.1-FTFP_BERT/HEC_SF.txt        |   2 +
 .../samf_eta_fwh-EM_e0100_sc1.dat             |  38 +
 .../samf_eta_rwh-EM_e0100_sc2.dat             |  38 +
 .../Analysis/10.1-FTFP_BERT_ATL/HEC_SF.txt    |   2 +
 .../samf_eta_fwh-EM_e0100_sc1.dat             |  38 +
 .../samf_eta_rwh-EM_e0100_sc2.dat             |  38 +
 .../Analysis/ClCaloHits.C                     |  19 +
 .../Analysis/ClCaloHits.h                     | 857 ++++++++++++++++++
 .../HEC_SamplingFractions/Analysis/env.sh     |  11 +
 .../HEC_SamplingFractions/Analysis/get_SF.C   | 422 +++++++++
 .../Analysis/makeclass_ClCaloHits.C           |  33 +
 .../HEC_SamplingFractions/Analysis/mypage.C   |  32 +
 .../HEC_SamplingFractions/Analysis/myplot.C   |  51 ++
 .../Analysis/read_var_cutX.C                  | 121 +++
 .../Analysis/rootlogon.C                      |  60 ++
 .../Analysis/samfr_etaX.C                     | 605 +++++++++++++
 .../Analysis/store_eta.C                      |  24 +
 .../HEC_SamplingFractions/Analysis/substa.C   |  70 ++
 .../jobOption.ParticleGun_ele_e0100_sc1.py    |  57 ++
 .../jobOption.ParticleGun_ele_e0100_sc2.py    |  57 ++
 .../HEC_SamplingFractions/Generation/env.sh   |   4 +
 .../jobOption.ParticleGun_gam_e0100_sc1.py    |  57 ++
 .../jobOption.ParticleGun_gam_e0100_sc2.py    |  57 ++
 .../jobOption.ParticleGun_pos_e0100_sc1.py    |  57 ++
 .../jobOption.ParticleGun_pos_e0100_sc2.py    |  57 ++
 .../share/HEC_SamplingFractions/ReadMe.txt    | 188 ++++
 .../Simulation/addMCpart.py                   |  13 +
 .../Simulation/dumphit.sh                     |   3 +
 .../HEC_SamplingFractions/Simulation/env.sh   |   2 +
 .../HEC_SamplingFractions/Simulation/mkjob.sh | 115 +++
 .../Simulation/preInclude.VertexPosOff.py     |  15 +
 .../HEC_SamplingFractions/Simulation/simul.sh |   3 +
 32 files changed, 3146 insertions(+)
 create mode 100644 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/10.1-FTFP_BERT/HEC_SF.txt
 create mode 100644 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/10.1-FTFP_BERT/samf_eta_fwh-EM_e0100_sc1.dat
 create mode 100644 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/10.1-FTFP_BERT/samf_eta_rwh-EM_e0100_sc2.dat
 create mode 100644 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/10.1-FTFP_BERT_ATL/HEC_SF.txt
 create mode 100644 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/10.1-FTFP_BERT_ATL/samf_eta_fwh-EM_e0100_sc1.dat
 create mode 100644 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/10.1-FTFP_BERT_ATL/samf_eta_rwh-EM_e0100_sc2.dat
 create mode 100644 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/ClCaloHits.C
 create mode 100644 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/ClCaloHits.h
 create mode 100755 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/env.sh
 create mode 100644 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/get_SF.C
 create mode 100644 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/makeclass_ClCaloHits.C
 create mode 100644 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/mypage.C
 create mode 100644 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/myplot.C
 create mode 100644 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/read_var_cutX.C
 create mode 100644 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/rootlogon.C
 create mode 100644 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/samfr_etaX.C
 create mode 100644 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/store_eta.C
 create mode 100644 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/substa.C
 create mode 100644 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Generation/ele_e0100_sc1/jobOption.ParticleGun_ele_e0100_sc1.py
 create mode 100644 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Generation/ele_e0100_sc2/jobOption.ParticleGun_ele_e0100_sc2.py
 create mode 100644 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Generation/env.sh
 create mode 100644 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Generation/gam_e0100_sc1/jobOption.ParticleGun_gam_e0100_sc1.py
 create mode 100644 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Generation/gam_e0100_sc2/jobOption.ParticleGun_gam_e0100_sc2.py
 create mode 100644 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Generation/pos_e0100_sc1/jobOption.ParticleGun_pos_e0100_sc1.py
 create mode 100644 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Generation/pos_e0100_sc2/jobOption.ParticleGun_pos_e0100_sc2.py
 create mode 100644 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/ReadMe.txt
 create mode 100644 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Simulation/addMCpart.py
 create mode 100644 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Simulation/dumphit.sh
 create mode 100755 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Simulation/env.sh
 create mode 100755 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Simulation/mkjob.sh
 create mode 100644 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Simulation/preInclude.VertexPosOff.py
 create mode 100644 Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Simulation/simul.sh

diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/10.1-FTFP_BERT/HEC_SF.txt b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/10.1-FTFP_BERT/HEC_SF.txt
new file mode 100644
index 000000000000..a8e66aef97f7
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/10.1-FTFP_BERT/HEC_SF.txt
@@ -0,0 +1,2 @@
+ HEC front wheel :  0.044273 +- 0.000040  ->  0.04427 
+ HEC  rear wheel :  0.022708 +- 0.000021  ->  0.02271 
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/10.1-FTFP_BERT/samf_eta_fwh-EM_e0100_sc1.dat b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/10.1-FTFP_BERT/samf_eta_fwh-EM_e0100_sc1.dat
new file mode 100644
index 000000000000..6df658e44cad
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/10.1-FTFP_BERT/samf_eta_fwh-EM_e0100_sc1.dat
@@ -0,0 +1,38 @@
+   36 
+     1.500     3.300     0.050 
+     1.525      402       0.0428589       0.0028954 
+     1.575      435       0.0438664       0.0003374 
+     1.625      411       0.0441752       0.0001435 
+     1.675      396       0.0449313       0.0001836 
+     1.725      447       0.0448591       0.0000912 
+     1.775      423       0.0433535       0.0002865 
+     1.825      405       0.0442279       0.0001449 
+     1.875      414       0.0444864       0.0002880 
+     1.925      405       0.0446255       0.0001457 
+     1.975      420       0.0444718       0.0002480 
+     2.025      408       0.0443906       0.0001448 
+     2.075      423       0.0425058       0.0002507 
+     2.125      417       0.0437704       0.0001140 
+     2.175      438       0.0442831       0.0001958 
+     2.225      429       0.0442608       0.0002194 
+     2.275      390       0.0442893       0.0003061 
+     2.325      426       0.0441252       0.0002982 
+     2.375      408       0.0441959       0.0002692 
+     2.425      387       0.0441293       0.0002478 
+     2.475      408       0.0437028       0.0002502 
+     2.525      459       0.0441487       0.0002711 
+     2.575      477       0.0439979       0.0003421 
+     2.625      363       0.0439880       0.0001778 
+     2.675      384       0.0437530       0.0001590 
+     2.725      402       0.0410577       0.0003695 
+     2.775      408       0.0408122       0.0002791 
+     2.825      462       0.0431493       0.0004233 
+     2.875      423       0.0433043       0.0003165 
+     2.925      429       0.0434001       0.0002532 
+     2.975      393       0.0415793       0.0003663 
+     3.025      387       0.0417385       0.0002838 
+     3.075      399       0.0408984       0.0004739 
+     3.125      426       0.0406061       0.0004078 
+     3.175      390       0.0345971       0.0014713 
+     3.225      459       0.0237174       0.0016906 
+     3.275      447       0.0405070       0.0010892 
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/10.1-FTFP_BERT/samf_eta_rwh-EM_e0100_sc2.dat b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/10.1-FTFP_BERT/samf_eta_rwh-EM_e0100_sc2.dat
new file mode 100644
index 000000000000..1aaf72f7917b
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/10.1-FTFP_BERT/samf_eta_rwh-EM_e0100_sc2.dat
@@ -0,0 +1,38 @@
+   36 
+     1.500     3.300     0.050 
+     1.525        0       0.0000000       0.0000000 
+     1.575        0       0.0000000       0.0000000 
+     1.625      429       0.0250368       0.0022645 
+     1.675      462       0.0200280       0.0019277 
+     1.725      408       0.0227120       0.0001548 
+     1.775      465       0.0223796       0.0001281 
+     1.825      474       0.0229359       0.0000520 
+     1.875      420       0.0227940       0.0001272 
+     1.925      453       0.0220481       0.0001195 
+     1.975      414       0.0223090       0.0001571 
+     2.025      459       0.0228444       0.0000918 
+     2.075      405       0.0227141       0.0001144 
+     2.125      444       0.0226598       0.0000917 
+     2.175      447       0.0228192       0.0000671 
+     2.225      480       0.0224611       0.0000832 
+     2.275      450       0.0209032       0.0001528 
+     2.325      411       0.0225525       0.0001378 
+     2.375      441       0.0226695       0.0001323 
+     2.425      435       0.0226930       0.0001203 
+     2.475      390       0.0226514       0.0001235 
+     2.525      450       0.0224885       0.0001548 
+     2.575      480       0.0228315       0.0000628 
+     2.625      492       0.0225280       0.0001569 
+     2.675      408       0.0227350       0.0000703 
+     2.725      408       0.0224340       0.0001173 
+     2.775      420       0.0226175       0.0000989 
+     2.825      471       0.0226179       0.0001370 
+     2.875      441       0.0218483       0.0001697 
+     2.925      477       0.0195613       0.0002101 
+     2.975      441       0.0212897       0.0001517 
+     3.025      369       0.0221725       0.0000942 
+     3.075      450       0.0213549       0.0001538 
+     3.125      438       0.0195638       0.0003212 
+     3.175      429       0.0179608       0.0016190 
+     3.225      474       0.0249693       0.0012896 
+     3.275      465       0.0265631       0.0012719 
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/10.1-FTFP_BERT_ATL/HEC_SF.txt b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/10.1-FTFP_BERT_ATL/HEC_SF.txt
new file mode 100644
index 000000000000..3f19a131474f
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/10.1-FTFP_BERT_ATL/HEC_SF.txt
@@ -0,0 +1,2 @@
+ HEC front wheel :  0.044269 +- 0.000040  ->  0.04427 
+ HEC  rear wheel :  0.022707 +- 0.000021  ->  0.02271 
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/10.1-FTFP_BERT_ATL/samf_eta_fwh-EM_e0100_sc1.dat b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/10.1-FTFP_BERT_ATL/samf_eta_fwh-EM_e0100_sc1.dat
new file mode 100644
index 000000000000..b48429bc1c37
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/10.1-FTFP_BERT_ATL/samf_eta_fwh-EM_e0100_sc1.dat
@@ -0,0 +1,38 @@
+   36 
+     1.500     3.300     0.050 
+     1.525      402       0.0428600       0.0028955 
+     1.575      435       0.0438666       0.0003375 
+     1.625      411       0.0441752       0.0001435 
+     1.675      396       0.0449419       0.0001911 
+     1.725      447       0.0448591       0.0000912 
+     1.775      423       0.0433522       0.0002865 
+     1.825      405       0.0442234       0.0001448 
+     1.875      414       0.0444921       0.0002872 
+     1.925      405       0.0446148       0.0001464 
+     1.975      420       0.0444718       0.0002480 
+     2.025      408       0.0443906       0.0001448 
+     2.075      423       0.0425058       0.0002507 
+     2.125      417       0.0437689       0.0001139 
+     2.175      438       0.0442831       0.0001958 
+     2.225      429       0.0442582       0.0002212 
+     2.275      390       0.0442789       0.0003062 
+     2.325      426       0.0441242       0.0002983 
+     2.375      408       0.0441984       0.0002692 
+     2.425      387       0.0441306       0.0002477 
+     2.475      408       0.0437026       0.0002503 
+     2.525      459       0.0441372       0.0002710 
+     2.575      477       0.0439975       0.0003421 
+     2.625      363       0.0439806       0.0001776 
+     2.675      384       0.0437530       0.0001590 
+     2.725      402       0.0410577       0.0003695 
+     2.775      408       0.0408137       0.0002791 
+     2.825      462       0.0431455       0.0004233 
+     2.875      423       0.0433051       0.0003165 
+     2.925      429       0.0434012       0.0002531 
+     2.975      393       0.0415793       0.0003663 
+     3.025      387       0.0417326       0.0002846 
+     3.075      399       0.0408962       0.0004740 
+     3.125      426       0.0406061       0.0004078 
+     3.175      390       0.0345971       0.0014713 
+     3.225      459       0.0237441       0.0016943 
+     3.275      447       0.0404975       0.0010890 
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/10.1-FTFP_BERT_ATL/samf_eta_rwh-EM_e0100_sc2.dat b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/10.1-FTFP_BERT_ATL/samf_eta_rwh-EM_e0100_sc2.dat
new file mode 100644
index 000000000000..5148af5071dd
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/10.1-FTFP_BERT_ATL/samf_eta_rwh-EM_e0100_sc2.dat
@@ -0,0 +1,38 @@
+   36 
+     1.500     3.300     0.050 
+     1.525        0       0.0000000       0.0000000 
+     1.575        0       0.0000000       0.0000000 
+     1.625      429       0.0250425       0.0022743 
+     1.675      462       0.0200342       0.0019275 
+     1.725      408       0.0227120       0.0001548 
+     1.775      465       0.0223781       0.0001265 
+     1.825      474       0.0229317       0.0000527 
+     1.875      420       0.0228000       0.0001282 
+     1.925      453       0.0220481       0.0001195 
+     1.975      414       0.0223129       0.0001574 
+     2.025      459       0.0228444       0.0000918 
+     2.075      405       0.0227264       0.0001140 
+     2.125      444       0.0226598       0.0000917 
+     2.175      447       0.0228196       0.0000670 
+     2.225      480       0.0224611       0.0000832 
+     2.275      450       0.0209043       0.0001529 
+     2.325      411       0.0225525       0.0001378 
+     2.375      441       0.0226742       0.0001322 
+     2.425      435       0.0226930       0.0001203 
+     2.475      390       0.0226461       0.0001232 
+     2.525      450       0.0224885       0.0001548 
+     2.575      480       0.0228329       0.0000628 
+     2.625      492       0.0225280       0.0001569 
+     2.675      408       0.0227368       0.0000704 
+     2.725      408       0.0224358       0.0001173 
+     2.775      420       0.0226147       0.0000987 
+     2.825      471       0.0226179       0.0001370 
+     2.875      441       0.0218472       0.0001697 
+     2.925      477       0.0195611       0.0002102 
+     2.975      441       0.0212865       0.0001517 
+     3.025      369       0.0221744       0.0000941 
+     3.075      450       0.0213498       0.0001541 
+     3.125      438       0.0195638       0.0003212 
+     3.175      429       0.0179608       0.0016190 
+     3.225      474       0.0250020       0.0014116 
+     3.275      465       0.0265688       0.0012732 
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/ClCaloHits.C b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/ClCaloHits.C
new file mode 100644
index 000000000000..dbf6c4ac1cab
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/ClCaloHits.C
@@ -0,0 +1,19 @@
+#define ClCaloHits_cxx
+#include "ClCaloHits.h"
+#include <TH2.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+
+void ClCaloHits::Loop()
+{
+   if (fChain == 0) return;
+
+   Long64_t nentries = fChain->GetEntriesFast();
+
+   Long64_t nbytes = 0, nb = 0;
+   for (Long64_t jentry=0; jentry<nentries;jentry++) {
+      Long64_t ientry = LoadTree(jentry);
+      if (ientry < 0) break;
+      nb = fChain->GetEntry(jentry);   nbytes += nb;
+   }
+}
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/ClCaloHits.h b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/ClCaloHits.h
new file mode 100644
index 000000000000..4d1605d7898a
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/ClCaloHits.h
@@ -0,0 +1,857 @@
+//////////////////////////////////////////////////////////
+// This class has been automatically generated on
+// Thu Mar  5 13:04:33 2020 by ROOT version 6.18/04
+// from TTree calohitD3PD/calohitD3PD
+// found on file: /eos/user/k/kiryunin/HEC_SF/10.1-FTFP_BERT_ATL/DumpHits/ele_e0100_sc1/01.root
+//////////////////////////////////////////////////////////
+
+#ifndef ClCaloHits_h
+#define ClCaloHits_h
+
+#include <TROOT.h>
+#include <TChain.h>
+#include <TFile.h>
+
+// Header file for the classes stored in the TTree if any.
+#include "vector"
+#include "vector"
+#include "vector"
+
+class ClCaloHits {
+public :
+   TTree          *fChain;   //!pointer to the analyzed TTree or TChain
+   Int_t           fCurrent; //!current Tree number in a TChain
+
+// Fixed size dimensions of array or collections stored in the TTree if any.
+
+   // Declaration of leaf types
+   UInt_t          RunNumber;
+   ULong64_t       EventNumber;
+   UInt_t          timestamp;
+   UInt_t          timestamp_ns;
+   UInt_t          lbn;
+   UInt_t          bcid;
+   UInt_t          detmask0;
+   UInt_t          detmask1;
+   Float_t         actualIntPerXing;
+   Float_t         averageIntPerXing;
+   UInt_t          pixelFlags;
+   UInt_t          sctFlags;
+   UInt_t          trtFlags;
+   UInt_t          larFlags;
+   UInt_t          tileFlags;
+   UInt_t          muonFlags;
+   UInt_t          fwdFlags;
+   UInt_t          coreFlags;
+   UInt_t          backgroundFlags;
+   UInt_t          lumiFlags;
+   UInt_t          pixelError;
+   UInt_t          sctError;
+   UInt_t          trtError;
+   UInt_t          larError;
+   UInt_t          tileError;
+   UInt_t          muonError;
+   UInt_t          fwdError;
+   UInt_t          coreError;
+   Bool_t          streamDecision_Egamma;
+   Bool_t          streamDecision_Muons;
+   Bool_t          streamDecision_JetTauEtmiss;
+   UInt_t          l1id;
+   Bool_t          isSimulation;
+   Bool_t          isCalibration;
+   Bool_t          isTestBeam;
+   Int_t           hitemb_n;
+   vector<float>   *hitemb_eta;
+   vector<float>   *hitemb_phi;
+   vector<float>   *hitemb_E;
+   vector<float>   *hitemb_Time;
+   vector<unsigned int> *hitemb_ID;
+   Int_t           hitemec_n;
+   vector<float>   *hitemec_eta;
+   vector<float>   *hitemec_phi;
+   vector<float>   *hitemec_E;
+   vector<float>   *hitemec_Time;
+   vector<unsigned int> *hitemec_ID;
+   Int_t           hithec_n;
+   vector<float>   *hithec_eta;
+   vector<float>   *hithec_phi;
+   vector<float>   *hithec_E;
+   vector<float>   *hithec_Time;
+   vector<unsigned int> *hithec_ID;
+   Int_t           hitfcal_n;
+   vector<float>   *hitfcal_eta;
+   vector<float>   *hitfcal_phi;
+   vector<float>   *hitfcal_E;
+   vector<float>   *hitfcal_Time;
+   vector<unsigned int> *hitfcal_ID;
+   Int_t           TileHit_n;
+   vector<vector<float> > *TileHit_energy;
+   vector<vector<float> > *TileHit_time;
+   vector<vector<int> > *TileHit_detector;
+   vector<vector<int> > *TileHit_side;
+   vector<vector<int> > *TileHit_sample;
+   vector<vector<int> > *TileHit_pmt;
+   vector<vector<int> > *TileHit_eta;
+   vector<vector<int> > *TileHit_phi;
+   Int_t           MBTSHit_n;
+   vector<vector<float> > *MBTSHit_energy;
+   vector<vector<float> > *MBTSHit_time;
+   vector<vector<int> > *MBTSHit_detector;
+   vector<vector<int> > *MBTSHit_side;
+   vector<vector<int> > *MBTSHit_sample;
+   vector<vector<int> > *MBTSHit_pmt;
+   vector<vector<int> > *MBTSHit_eta;
+   vector<vector<int> > *MBTSHit_phi;
+   Int_t           laract_n;
+   vector<float>   *laract_eta;
+   vector<float>   *laract_phi;
+   vector<unsigned int> *laract_ID;
+   vector<float>   *laract_EnergyTot;
+   vector<float>   *laract_EnergyVis;
+   vector<unsigned int> *laract_detector;
+   vector<unsigned int> *laract_sample;
+   vector<unsigned int> *laract_side;
+   vector<unsigned int> *laract_reg_sec;
+   vector<unsigned int> *laract_eta_tower;
+   vector<unsigned int> *laract_phi_module;
+   vector<float>   *laract_EnergyEm;
+   vector<float>   *laract_EnergyNonEm;
+   vector<float>   *laract_EnergyInv;
+   vector<float>   *laract_EnergyEsc;
+   Int_t           larinact_n;
+   vector<float>   *larinact_eta;
+   vector<float>   *larinact_phi;
+   vector<unsigned int> *larinact_ID;
+   vector<float>   *larinact_EnergyTot;
+   vector<float>   *larinact_EnergyVis;
+   vector<unsigned int> *larinact_detector;
+   vector<unsigned int> *larinact_sample;
+   vector<unsigned int> *larinact_side;
+   vector<unsigned int> *larinact_reg_sec;
+   vector<unsigned int> *larinact_eta_tower;
+   vector<unsigned int> *larinact_phi_module;
+   vector<float>   *larinact_EnergyEm;
+   vector<float>   *larinact_EnergyNonEm;
+   vector<float>   *larinact_EnergyInv;
+   vector<float>   *larinact_EnergyEsc;
+   Int_t           tileact_n;
+   vector<float>   *tileact_eta;
+   vector<float>   *tileact_phi;
+   vector<unsigned int> *tileact_ID;
+   vector<float>   *tileact_EnergyTot;
+   vector<float>   *tileact_EnergyVis;
+   vector<unsigned int> *tileact_detector;
+   vector<unsigned int> *tileact_sample;
+   vector<unsigned int> *tileact_side;
+   vector<unsigned int> *tileact_reg_sec;
+   vector<unsigned int> *tileact_eta_tower;
+   vector<unsigned int> *tileact_phi_module;
+   vector<float>   *tileact_EnergyEm;
+   vector<float>   *tileact_EnergyNonEm;
+   vector<float>   *tileact_EnergyInv;
+   vector<float>   *tileact_EnergyEsc;
+   Int_t           tileinact_n;
+   vector<float>   *tileinact_eta;
+   vector<float>   *tileinact_phi;
+   vector<unsigned int> *tileinact_ID;
+   vector<float>   *tileinact_EnergyTot;
+   vector<float>   *tileinact_EnergyVis;
+   vector<unsigned int> *tileinact_detector;
+   vector<unsigned int> *tileinact_sample;
+   vector<unsigned int> *tileinact_side;
+   vector<unsigned int> *tileinact_reg_sec;
+   vector<unsigned int> *tileinact_eta_tower;
+   vector<unsigned int> *tileinact_phi_module;
+   vector<float>   *tileinact_EnergyEm;
+   vector<float>   *tileinact_EnergyNonEm;
+   vector<float>   *tileinact_EnergyInv;
+   vector<float>   *tileinact_EnergyEsc;
+   Int_t           lardm_n;
+   vector<float>   *lardm_eta;
+   vector<float>   *lardm_phi;
+   vector<unsigned int> *lardm_ID;
+   vector<float>   *lardm_EnergyTot;
+   vector<float>   *lardm_EnergyVis;
+   vector<unsigned int> *lardm_detzside;
+   vector<unsigned int> *lardm_type;
+   vector<unsigned int> *lardm_sampling;
+   vector<unsigned int> *lardm_region;
+   vector<unsigned int> *lardm_ieta;
+   vector<unsigned int> *lardm_iphi;
+   vector<float>   *lardm_EnergyEm;
+   vector<float>   *lardm_EnergyNonEm;
+   vector<float>   *lardm_EnergyInv;
+   vector<float>   *lardm_EnergyEsc;
+   Int_t           tiledm_n;
+   vector<float>   *tiledm_eta;
+   vector<float>   *tiledm_phi;
+   vector<unsigned int> *tiledm_ID;
+   vector<float>   *tiledm_EnergyTot;
+   vector<float>   *tiledm_EnergyVis;
+   vector<unsigned int> *tiledm_detzside;
+   vector<unsigned int> *tiledm_type;
+   vector<unsigned int> *tiledm_sampling;
+   vector<unsigned int> *tiledm_region;
+   vector<unsigned int> *tiledm_ieta;
+   vector<unsigned int> *tiledm_iphi;
+   vector<float>   *tiledm_EnergyEm;
+   vector<float>   *tiledm_EnergyNonEm;
+   vector<float>   *tiledm_EnergyInv;
+   vector<float>   *tiledm_EnergyEsc;
+   Int_t           mcvtx_n;
+   vector<float>   *mcvtx_x;
+   vector<float>   *mcvtx_y;
+   vector<float>   *mcvtx_z;
+   vector<int>     *mcvtx_barcode;
+   Int_t           mcpart_n;
+   vector<float>   *mcpart_pt;
+   vector<float>   *mcpart_m;
+   vector<float>   *mcpart_eta;
+   vector<float>   *mcpart_phi;
+   vector<int>     *mcpart_type;
+   vector<int>     *mcpart_status;
+   vector<int>     *mcpart_barcode;
+   vector<int>     *mcpart_mothertype;
+   vector<int>     *mcpart_motherbarcode;
+   vector<int>     *mcpart_mcprodvtx_index;
+   vector<int>     *mcpart_mother_n;
+   vector<vector<int> > *mcpart_mother_index;
+   vector<int>     *mcpart_mcdecayvtx_index;
+   vector<int>     *mcpart_child_n;
+   vector<vector<int> > *mcpart_child_index;
+
+   // List of branches
+   TBranch        *b_RunNumber;   //!
+   TBranch        *b_EventNumber;   //!
+   TBranch        *b_timestamp;   //!
+   TBranch        *b_timestamp_ns;   //!
+   TBranch        *b_lbn;   //!
+   TBranch        *b_bcid;   //!
+   TBranch        *b_detmask0;   //!
+   TBranch        *b_detmask1;   //!
+   TBranch        *b_actualIntPerXing;   //!
+   TBranch        *b_averageIntPerXing;   //!
+   TBranch        *b_pixelFlags;   //!
+   TBranch        *b_sctFlags;   //!
+   TBranch        *b_trtFlags;   //!
+   TBranch        *b_larFlags;   //!
+   TBranch        *b_tileFlags;   //!
+   TBranch        *b_muonFlags;   //!
+   TBranch        *b_fwdFlags;   //!
+   TBranch        *b_coreFlags;   //!
+   TBranch        *b_backgroundFlags;   //!
+   TBranch        *b_lumiFlags;   //!
+   TBranch        *b_pixelError;   //!
+   TBranch        *b_sctError;   //!
+   TBranch        *b_trtError;   //!
+   TBranch        *b_larError;   //!
+   TBranch        *b_tileError;   //!
+   TBranch        *b_muonError;   //!
+   TBranch        *b_fwdError;   //!
+   TBranch        *b_coreError;   //!
+   TBranch        *b_streamDecision_Egamma;   //!
+   TBranch        *b_streamDecision_Muons;   //!
+   TBranch        *b_streamDecision_JetTauEtmiss;   //!
+   TBranch        *b_l1id;   //!
+   TBranch        *b_isSimulation;   //!
+   TBranch        *b_isCalibration;   //!
+   TBranch        *b_isTestBeam;   //!
+   TBranch        *b_hitemb_n;   //!
+   TBranch        *b_hitemb_eta;   //!
+   TBranch        *b_hitemb_phi;   //!
+   TBranch        *b_hitemb_E;   //!
+   TBranch        *b_hitemb_Time;   //!
+   TBranch        *b_hitemb_ID;   //!
+   TBranch        *b_hitemec_n;   //!
+   TBranch        *b_hitemec_eta;   //!
+   TBranch        *b_hitemec_phi;   //!
+   TBranch        *b_hitemec_E;   //!
+   TBranch        *b_hitemec_Time;   //!
+   TBranch        *b_hitemec_ID;   //!
+   TBranch        *b_hithec_n;   //!
+   TBranch        *b_hithec_eta;   //!
+   TBranch        *b_hithec_phi;   //!
+   TBranch        *b_hithec_E;   //!
+   TBranch        *b_hithec_Time;   //!
+   TBranch        *b_hithec_ID;   //!
+   TBranch        *b_hitfcal_n;   //!
+   TBranch        *b_hitfcal_eta;   //!
+   TBranch        *b_hitfcal_phi;   //!
+   TBranch        *b_hitfcal_E;   //!
+   TBranch        *b_hitfcal_Time;   //!
+   TBranch        *b_hitfcal_ID;   //!
+   TBranch        *b_TileHit_n;   //!
+   TBranch        *b_TileHit_energy;   //!
+   TBranch        *b_TileHit_time;   //!
+   TBranch        *b_TileHit_detector;   //!
+   TBranch        *b_TileHit_side;   //!
+   TBranch        *b_TileHit_sample;   //!
+   TBranch        *b_TileHit_pmt;   //!
+   TBranch        *b_TileHit_eta;   //!
+   TBranch        *b_TileHit_phi;   //!
+   TBranch        *b_MBTSHit_n;   //!
+   TBranch        *b_MBTSHit_energy;   //!
+   TBranch        *b_MBTSHit_time;   //!
+   TBranch        *b_MBTSHit_detector;   //!
+   TBranch        *b_MBTSHit_side;   //!
+   TBranch        *b_MBTSHit_sample;   //!
+   TBranch        *b_MBTSHit_pmt;   //!
+   TBranch        *b_MBTSHit_eta;   //!
+   TBranch        *b_MBTSHit_phi;   //!
+   TBranch        *b_laract_n;   //!
+   TBranch        *b_laract_eta;   //!
+   TBranch        *b_laract_phi;   //!
+   TBranch        *b_laract_ID;   //!
+   TBranch        *b_laract_EnergyTot;   //!
+   TBranch        *b_laract_EnergyVis;   //!
+   TBranch        *b_laract_detector;   //!
+   TBranch        *b_laract_sample;   //!
+   TBranch        *b_laract_side;   //!
+   TBranch        *b_laract_reg_sec;   //!
+   TBranch        *b_laract_eta_tower;   //!
+   TBranch        *b_laract_phi_module;   //!
+   TBranch        *b_laract_EnergyEm;   //!
+   TBranch        *b_laract_EnergyNonEm;   //!
+   TBranch        *b_laract_EnergyInv;   //!
+   TBranch        *b_laract_EnergyEsc;   //!
+   TBranch        *b_larinact_n;   //!
+   TBranch        *b_larinact_eta;   //!
+   TBranch        *b_larinact_phi;   //!
+   TBranch        *b_larinact_ID;   //!
+   TBranch        *b_larinact_EnergyTot;   //!
+   TBranch        *b_larinact_EnergyVis;   //!
+   TBranch        *b_larinact_detector;   //!
+   TBranch        *b_larinact_sample;   //!
+   TBranch        *b_larinact_side;   //!
+   TBranch        *b_larinact_reg_sec;   //!
+   TBranch        *b_larinact_eta_tower;   //!
+   TBranch        *b_larinact_phi_module;   //!
+   TBranch        *b_larinact_EnergyEm;   //!
+   TBranch        *b_larinact_EnergyNonEm;   //!
+   TBranch        *b_larinact_EnergyInv;   //!
+   TBranch        *b_larinact_EnergyEsc;   //!
+   TBranch        *b_tileact_n;   //!
+   TBranch        *b_tileact_eta;   //!
+   TBranch        *b_tileact_phi;   //!
+   TBranch        *b_tileact_ID;   //!
+   TBranch        *b_tileact_EnergyTot;   //!
+   TBranch        *b_tileact_EnergyVis;   //!
+   TBranch        *b_tileact_detector;   //!
+   TBranch        *b_tileact_sample;   //!
+   TBranch        *b_tileact_side;   //!
+   TBranch        *b_tileact_reg_sec;   //!
+   TBranch        *b_tileact_eta_tower;   //!
+   TBranch        *b_tileact_phi_module;   //!
+   TBranch        *b_tileact_EnergyEm;   //!
+   TBranch        *b_tileact_EnergyNonEm;   //!
+   TBranch        *b_tileact_EnergyInv;   //!
+   TBranch        *b_tileact_EnergyEsc;   //!
+   TBranch        *b_tileinact_n;   //!
+   TBranch        *b_tileinact_eta;   //!
+   TBranch        *b_tileinact_phi;   //!
+   TBranch        *b_tileinact_ID;   //!
+   TBranch        *b_tileinact_EnergyTot;   //!
+   TBranch        *b_tileinact_EnergyVis;   //!
+   TBranch        *b_tileinact_detector;   //!
+   TBranch        *b_tileinact_sample;   //!
+   TBranch        *b_tileinact_side;   //!
+   TBranch        *b_tileinact_reg_sec;   //!
+   TBranch        *b_tileinact_eta_tower;   //!
+   TBranch        *b_tileinact_phi_module;   //!
+   TBranch        *b_tileinact_EnergyEm;   //!
+   TBranch        *b_tileinact_EnergyNonEm;   //!
+   TBranch        *b_tileinact_EnergyInv;   //!
+   TBranch        *b_tileinact_EnergyEsc;   //!
+   TBranch        *b_lardm_n;   //!
+   TBranch        *b_lardm_eta;   //!
+   TBranch        *b_lardm_phi;   //!
+   TBranch        *b_lardm_ID;   //!
+   TBranch        *b_lardm_EnergyTot;   //!
+   TBranch        *b_lardm_EnergyVis;   //!
+   TBranch        *b_lardm_detzside;   //!
+   TBranch        *b_lardm_type;   //!
+   TBranch        *b_lardm_sampling;   //!
+   TBranch        *b_lardm_region;   //!
+   TBranch        *b_lardm_ieta;   //!
+   TBranch        *b_lardm_iphi;   //!
+   TBranch        *b_lardm_EnergyEm;   //!
+   TBranch        *b_lardm_EnergyNonEm;   //!
+   TBranch        *b_lardm_EnergyInv;   //!
+   TBranch        *b_lardm_EnergyEsc;   //!
+   TBranch        *b_tiledm_n;   //!
+   TBranch        *b_tiledm_eta;   //!
+   TBranch        *b_tiledm_phi;   //!
+   TBranch        *b_tiledm_ID;   //!
+   TBranch        *b_tiledm_EnergyTot;   //!
+   TBranch        *b_tiledm_EnergyVis;   //!
+   TBranch        *b_tiledm_detzside;   //!
+   TBranch        *b_tiledm_type;   //!
+   TBranch        *b_tiledm_sampling;   //!
+   TBranch        *b_tiledm_region;   //!
+   TBranch        *b_tiledm_ieta;   //!
+   TBranch        *b_tiledm_iphi;   //!
+   TBranch        *b_tiledm_EnergyEm;   //!
+   TBranch        *b_tiledm_EnergyNonEm;   //!
+   TBranch        *b_tiledm_EnergyInv;   //!
+   TBranch        *b_tiledm_EnergyEsc;   //!
+   TBranch        *b_mcvtx_n;   //!
+   TBranch        *b_mcvtx_x;   //!
+   TBranch        *b_mcvtx_y;   //!
+   TBranch        *b_mcvtx_z;   //!
+   TBranch        *b_mcvtx_barcode;   //!
+   TBranch        *b_mcpart_n;   //!
+   TBranch        *b_mcpart_pt;   //!
+   TBranch        *b_mcpart_m;   //!
+   TBranch        *b_mcpart_eta;   //!
+   TBranch        *b_mcpart_phi;   //!
+   TBranch        *b_mcpart_type;   //!
+   TBranch        *b_mcpart_status;   //!
+   TBranch        *b_mcpart_barcode;   //!
+   TBranch        *b_mcpart_mothertype;   //!
+   TBranch        *b_mcpart_motherbarcode;   //!
+   TBranch        *b_mcpart_mcprodvtx_index;   //!
+   TBranch        *b_mcpart_mother_n;   //!
+   TBranch        *b_mcpart_mother_index;   //!
+   TBranch        *b_mcpart_mcdecayvtx_index;   //!
+   TBranch        *b_mcpart_child_n;   //!
+   TBranch        *b_mcpart_child_index;   //!
+
+   ClCaloHits(TTree *tree=0);
+   virtual ~ClCaloHits();
+   virtual Int_t    Cut(Long64_t entry);
+   virtual Int_t    GetEntry(Long64_t entry);
+   virtual Long64_t LoadTree(Long64_t entry);
+   virtual void     Init(TTree *tree);
+   virtual void     Loop();
+   virtual Bool_t   Notify();
+   virtual void     Show(Long64_t entry = -1);
+};
+
+#endif
+
+#ifdef ClCaloHits_cxx
+ClCaloHits::ClCaloHits(TTree *tree) : fChain(0) 
+{
+// if parameter tree is not specified (or zero), connect the file
+// used to generate this class and read the Tree.
+   if (tree == 0) {
+      TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("/eos/user/k/kiryunin/HEC_SF/10.1-FTFP_BERT_ATL/DumpHits/ele_e0100_sc1/01.root");
+      if (!f || !f->IsOpen()) {
+         f = new TFile("/eos/user/k/kiryunin/HEC_SF/10.1-FTFP_BERT_ATL/DumpHits/ele_e0100_sc1/01.root");
+      }
+      f->GetObject("calohitD3PD",tree);
+
+   }
+   Init(tree);
+}
+
+ClCaloHits::~ClCaloHits()
+{
+   if (!fChain) return;
+   delete fChain->GetCurrentFile();
+}
+
+Int_t ClCaloHits::GetEntry(Long64_t entry)
+{
+// Read contents of entry.
+   if (!fChain) return 0;
+   return fChain->GetEntry(entry);
+}
+Long64_t ClCaloHits::LoadTree(Long64_t entry)
+{
+// Set the environment to read one entry
+   if (!fChain) return -5;
+   Long64_t centry = fChain->LoadTree(entry);
+   if (centry < 0) return centry;
+   if (fChain->GetTreeNumber() != fCurrent) {
+      fCurrent = fChain->GetTreeNumber();
+      Notify();
+   }
+   return centry;
+}
+
+void ClCaloHits::Init(TTree *tree)
+{
+   // The Init() function is called when the selector needs to initialize
+   // a new tree or chain. Typically here the branch addresses and branch
+   // pointers of the tree will be set.
+   // It is normally not necessary to make changes to the generated
+   // code, but the routine can be extended by the user if needed.
+   // Init() will be called many times when running on PROOF
+   // (once per file to be processed).
+
+   // Set object pointer
+   hitemb_eta = 0;
+   hitemb_phi = 0;
+   hitemb_E = 0;
+   hitemb_Time = 0;
+   hitemb_ID = 0;
+   hitemec_eta = 0;
+   hitemec_phi = 0;
+   hitemec_E = 0;
+   hitemec_Time = 0;
+   hitemec_ID = 0;
+   hithec_eta = 0;
+   hithec_phi = 0;
+   hithec_E = 0;
+   hithec_Time = 0;
+   hithec_ID = 0;
+   hitfcal_eta = 0;
+   hitfcal_phi = 0;
+   hitfcal_E = 0;
+   hitfcal_Time = 0;
+   hitfcal_ID = 0;
+   TileHit_energy = 0;
+   TileHit_time = 0;
+   TileHit_detector = 0;
+   TileHit_side = 0;
+   TileHit_sample = 0;
+   TileHit_pmt = 0;
+   TileHit_eta = 0;
+   TileHit_phi = 0;
+   MBTSHit_energy = 0;
+   MBTSHit_time = 0;
+   MBTSHit_detector = 0;
+   MBTSHit_side = 0;
+   MBTSHit_sample = 0;
+   MBTSHit_pmt = 0;
+   MBTSHit_eta = 0;
+   MBTSHit_phi = 0;
+   laract_eta = 0;
+   laract_phi = 0;
+   laract_ID = 0;
+   laract_EnergyTot = 0;
+   laract_EnergyVis = 0;
+   laract_detector = 0;
+   laract_sample = 0;
+   laract_side = 0;
+   laract_reg_sec = 0;
+   laract_eta_tower = 0;
+   laract_phi_module = 0;
+   laract_EnergyEm = 0;
+   laract_EnergyNonEm = 0;
+   laract_EnergyInv = 0;
+   laract_EnergyEsc = 0;
+   larinact_eta = 0;
+   larinact_phi = 0;
+   larinact_ID = 0;
+   larinact_EnergyTot = 0;
+   larinact_EnergyVis = 0;
+   larinact_detector = 0;
+   larinact_sample = 0;
+   larinact_side = 0;
+   larinact_reg_sec = 0;
+   larinact_eta_tower = 0;
+   larinact_phi_module = 0;
+   larinact_EnergyEm = 0;
+   larinact_EnergyNonEm = 0;
+   larinact_EnergyInv = 0;
+   larinact_EnergyEsc = 0;
+   tileact_eta = 0;
+   tileact_phi = 0;
+   tileact_ID = 0;
+   tileact_EnergyTot = 0;
+   tileact_EnergyVis = 0;
+   tileact_detector = 0;
+   tileact_sample = 0;
+   tileact_side = 0;
+   tileact_reg_sec = 0;
+   tileact_eta_tower = 0;
+   tileact_phi_module = 0;
+   tileact_EnergyEm = 0;
+   tileact_EnergyNonEm = 0;
+   tileact_EnergyInv = 0;
+   tileact_EnergyEsc = 0;
+   tileinact_eta = 0;
+   tileinact_phi = 0;
+   tileinact_ID = 0;
+   tileinact_EnergyTot = 0;
+   tileinact_EnergyVis = 0;
+   tileinact_detector = 0;
+   tileinact_sample = 0;
+   tileinact_side = 0;
+   tileinact_reg_sec = 0;
+   tileinact_eta_tower = 0;
+   tileinact_phi_module = 0;
+   tileinact_EnergyEm = 0;
+   tileinact_EnergyNonEm = 0;
+   tileinact_EnergyInv = 0;
+   tileinact_EnergyEsc = 0;
+   lardm_eta = 0;
+   lardm_phi = 0;
+   lardm_ID = 0;
+   lardm_EnergyTot = 0;
+   lardm_EnergyVis = 0;
+   lardm_detzside = 0;
+   lardm_type = 0;
+   lardm_sampling = 0;
+   lardm_region = 0;
+   lardm_ieta = 0;
+   lardm_iphi = 0;
+   lardm_EnergyEm = 0;
+   lardm_EnergyNonEm = 0;
+   lardm_EnergyInv = 0;
+   lardm_EnergyEsc = 0;
+   tiledm_eta = 0;
+   tiledm_phi = 0;
+   tiledm_ID = 0;
+   tiledm_EnergyTot = 0;
+   tiledm_EnergyVis = 0;
+   tiledm_detzside = 0;
+   tiledm_type = 0;
+   tiledm_sampling = 0;
+   tiledm_region = 0;
+   tiledm_ieta = 0;
+   tiledm_iphi = 0;
+   tiledm_EnergyEm = 0;
+   tiledm_EnergyNonEm = 0;
+   tiledm_EnergyInv = 0;
+   tiledm_EnergyEsc = 0;
+   mcvtx_x = 0;
+   mcvtx_y = 0;
+   mcvtx_z = 0;
+   mcvtx_barcode = 0;
+   mcpart_pt = 0;
+   mcpart_m = 0;
+   mcpart_eta = 0;
+   mcpart_phi = 0;
+   mcpart_type = 0;
+   mcpart_status = 0;
+   mcpart_barcode = 0;
+   mcpart_mothertype = 0;
+   mcpart_motherbarcode = 0;
+   mcpart_mcprodvtx_index = 0;
+   mcpart_mother_n = 0;
+   mcpart_mother_index = 0;
+   mcpart_mcdecayvtx_index = 0;
+   mcpart_child_n = 0;
+   mcpart_child_index = 0;
+   // Set branch addresses and branch pointers
+   if (!tree) return;
+   fChain = tree;
+   fCurrent = -1;
+   fChain->SetMakeClass(1);
+
+   fChain->SetBranchAddress("RunNumber", &RunNumber, &b_RunNumber);
+   fChain->SetBranchAddress("EventNumber", &EventNumber, &b_EventNumber);
+   fChain->SetBranchAddress("timestamp", &timestamp, &b_timestamp);
+   fChain->SetBranchAddress("timestamp_ns", &timestamp_ns, &b_timestamp_ns);
+   fChain->SetBranchAddress("lbn", &lbn, &b_lbn);
+   fChain->SetBranchAddress("bcid", &bcid, &b_bcid);
+   fChain->SetBranchAddress("detmask0", &detmask0, &b_detmask0);
+   fChain->SetBranchAddress("detmask1", &detmask1, &b_detmask1);
+   fChain->SetBranchAddress("actualIntPerXing", &actualIntPerXing, &b_actualIntPerXing);
+   fChain->SetBranchAddress("averageIntPerXing", &averageIntPerXing, &b_averageIntPerXing);
+   fChain->SetBranchAddress("pixelFlags", &pixelFlags, &b_pixelFlags);
+   fChain->SetBranchAddress("sctFlags", &sctFlags, &b_sctFlags);
+   fChain->SetBranchAddress("trtFlags", &trtFlags, &b_trtFlags);
+   fChain->SetBranchAddress("larFlags", &larFlags, &b_larFlags);
+   fChain->SetBranchAddress("tileFlags", &tileFlags, &b_tileFlags);
+   fChain->SetBranchAddress("muonFlags", &muonFlags, &b_muonFlags);
+   fChain->SetBranchAddress("fwdFlags", &fwdFlags, &b_fwdFlags);
+   fChain->SetBranchAddress("coreFlags", &coreFlags, &b_coreFlags);
+   fChain->SetBranchAddress("backgroundFlags", &backgroundFlags, &b_backgroundFlags);
+   fChain->SetBranchAddress("lumiFlags", &lumiFlags, &b_lumiFlags);
+   fChain->SetBranchAddress("pixelError", &pixelError, &b_pixelError);
+   fChain->SetBranchAddress("sctError", &sctError, &b_sctError);
+   fChain->SetBranchAddress("trtError", &trtError, &b_trtError);
+   fChain->SetBranchAddress("larError", &larError, &b_larError);
+   fChain->SetBranchAddress("tileError", &tileError, &b_tileError);
+   fChain->SetBranchAddress("muonError", &muonError, &b_muonError);
+   fChain->SetBranchAddress("fwdError", &fwdError, &b_fwdError);
+   fChain->SetBranchAddress("coreError", &coreError, &b_coreError);
+   fChain->SetBranchAddress("streamDecision_Egamma", &streamDecision_Egamma, &b_streamDecision_Egamma);
+   fChain->SetBranchAddress("streamDecision_Muons", &streamDecision_Muons, &b_streamDecision_Muons);
+   fChain->SetBranchAddress("streamDecision_JetTauEtmiss", &streamDecision_JetTauEtmiss, &b_streamDecision_JetTauEtmiss);
+   fChain->SetBranchAddress("l1id", &l1id, &b_l1id);
+   fChain->SetBranchAddress("isSimulation", &isSimulation, &b_isSimulation);
+   fChain->SetBranchAddress("isCalibration", &isCalibration, &b_isCalibration);
+   fChain->SetBranchAddress("isTestBeam", &isTestBeam, &b_isTestBeam);
+   fChain->SetBranchAddress("hitemb_n", &hitemb_n, &b_hitemb_n);
+   fChain->SetBranchAddress("hitemb_eta", &hitemb_eta, &b_hitemb_eta);
+   fChain->SetBranchAddress("hitemb_phi", &hitemb_phi, &b_hitemb_phi);
+   fChain->SetBranchAddress("hitemb_E", &hitemb_E, &b_hitemb_E);
+   fChain->SetBranchAddress("hitemb_Time", &hitemb_Time, &b_hitemb_Time);
+   fChain->SetBranchAddress("hitemb_ID", &hitemb_ID, &b_hitemb_ID);
+   fChain->SetBranchAddress("hitemec_n", &hitemec_n, &b_hitemec_n);
+   fChain->SetBranchAddress("hitemec_eta", &hitemec_eta, &b_hitemec_eta);
+   fChain->SetBranchAddress("hitemec_phi", &hitemec_phi, &b_hitemec_phi);
+   fChain->SetBranchAddress("hitemec_E", &hitemec_E, &b_hitemec_E);
+   fChain->SetBranchAddress("hitemec_Time", &hitemec_Time, &b_hitemec_Time);
+   fChain->SetBranchAddress("hitemec_ID", &hitemec_ID, &b_hitemec_ID);
+   fChain->SetBranchAddress("hithec_n", &hithec_n, &b_hithec_n);
+   fChain->SetBranchAddress("hithec_eta", &hithec_eta, &b_hithec_eta);
+   fChain->SetBranchAddress("hithec_phi", &hithec_phi, &b_hithec_phi);
+   fChain->SetBranchAddress("hithec_E", &hithec_E, &b_hithec_E);
+   fChain->SetBranchAddress("hithec_Time", &hithec_Time, &b_hithec_Time);
+   fChain->SetBranchAddress("hithec_ID", &hithec_ID, &b_hithec_ID);
+   fChain->SetBranchAddress("hitfcal_n", &hitfcal_n, &b_hitfcal_n);
+   fChain->SetBranchAddress("hitfcal_eta", &hitfcal_eta, &b_hitfcal_eta);
+   fChain->SetBranchAddress("hitfcal_phi", &hitfcal_phi, &b_hitfcal_phi);
+   fChain->SetBranchAddress("hitfcal_E", &hitfcal_E, &b_hitfcal_E);
+   fChain->SetBranchAddress("hitfcal_Time", &hitfcal_Time, &b_hitfcal_Time);
+   fChain->SetBranchAddress("hitfcal_ID", &hitfcal_ID, &b_hitfcal_ID);
+   fChain->SetBranchAddress("TileHit_n", &TileHit_n, &b_TileHit_n);
+   fChain->SetBranchAddress("TileHit_energy", &TileHit_energy, &b_TileHit_energy);
+   fChain->SetBranchAddress("TileHit_time", &TileHit_time, &b_TileHit_time);
+   fChain->SetBranchAddress("TileHit_detector", &TileHit_detector, &b_TileHit_detector);
+   fChain->SetBranchAddress("TileHit_side", &TileHit_side, &b_TileHit_side);
+   fChain->SetBranchAddress("TileHit_sample", &TileHit_sample, &b_TileHit_sample);
+   fChain->SetBranchAddress("TileHit_pmt", &TileHit_pmt, &b_TileHit_pmt);
+   fChain->SetBranchAddress("TileHit_eta", &TileHit_eta, &b_TileHit_eta);
+   fChain->SetBranchAddress("TileHit_phi", &TileHit_phi, &b_TileHit_phi);
+   fChain->SetBranchAddress("MBTSHit_n", &MBTSHit_n, &b_MBTSHit_n);
+   fChain->SetBranchAddress("MBTSHit_energy", &MBTSHit_energy, &b_MBTSHit_energy);
+   fChain->SetBranchAddress("MBTSHit_time", &MBTSHit_time, &b_MBTSHit_time);
+   fChain->SetBranchAddress("MBTSHit_detector", &MBTSHit_detector, &b_MBTSHit_detector);
+   fChain->SetBranchAddress("MBTSHit_side", &MBTSHit_side, &b_MBTSHit_side);
+   fChain->SetBranchAddress("MBTSHit_sample", &MBTSHit_sample, &b_MBTSHit_sample);
+   fChain->SetBranchAddress("MBTSHit_pmt", &MBTSHit_pmt, &b_MBTSHit_pmt);
+   fChain->SetBranchAddress("MBTSHit_eta", &MBTSHit_eta, &b_MBTSHit_eta);
+   fChain->SetBranchAddress("MBTSHit_phi", &MBTSHit_phi, &b_MBTSHit_phi);
+   fChain->SetBranchAddress("laract_n", &laract_n, &b_laract_n);
+   fChain->SetBranchAddress("laract_eta", &laract_eta, &b_laract_eta);
+   fChain->SetBranchAddress("laract_phi", &laract_phi, &b_laract_phi);
+   fChain->SetBranchAddress("laract_ID", &laract_ID, &b_laract_ID);
+   fChain->SetBranchAddress("laract_EnergyTot", &laract_EnergyTot, &b_laract_EnergyTot);
+   fChain->SetBranchAddress("laract_EnergyVis", &laract_EnergyVis, &b_laract_EnergyVis);
+   fChain->SetBranchAddress("laract_detector", &laract_detector, &b_laract_detector);
+   fChain->SetBranchAddress("laract_sample", &laract_sample, &b_laract_sample);
+   fChain->SetBranchAddress("laract_side", &laract_side, &b_laract_side);
+   fChain->SetBranchAddress("laract_reg_sec", &laract_reg_sec, &b_laract_reg_sec);
+   fChain->SetBranchAddress("laract_eta_tower", &laract_eta_tower, &b_laract_eta_tower);
+   fChain->SetBranchAddress("laract_phi_module", &laract_phi_module, &b_laract_phi_module);
+   fChain->SetBranchAddress("laract_EnergyEm", &laract_EnergyEm, &b_laract_EnergyEm);
+   fChain->SetBranchAddress("laract_EnergyNonEm", &laract_EnergyNonEm, &b_laract_EnergyNonEm);
+   fChain->SetBranchAddress("laract_EnergyInv", &laract_EnergyInv, &b_laract_EnergyInv);
+   fChain->SetBranchAddress("laract_EnergyEsc", &laract_EnergyEsc, &b_laract_EnergyEsc);
+   fChain->SetBranchAddress("larinact_n", &larinact_n, &b_larinact_n);
+   fChain->SetBranchAddress("larinact_eta", &larinact_eta, &b_larinact_eta);
+   fChain->SetBranchAddress("larinact_phi", &larinact_phi, &b_larinact_phi);
+   fChain->SetBranchAddress("larinact_ID", &larinact_ID, &b_larinact_ID);
+   fChain->SetBranchAddress("larinact_EnergyTot", &larinact_EnergyTot, &b_larinact_EnergyTot);
+   fChain->SetBranchAddress("larinact_EnergyVis", &larinact_EnergyVis, &b_larinact_EnergyVis);
+   fChain->SetBranchAddress("larinact_detector", &larinact_detector, &b_larinact_detector);
+   fChain->SetBranchAddress("larinact_sample", &larinact_sample, &b_larinact_sample);
+   fChain->SetBranchAddress("larinact_side", &larinact_side, &b_larinact_side);
+   fChain->SetBranchAddress("larinact_reg_sec", &larinact_reg_sec, &b_larinact_reg_sec);
+   fChain->SetBranchAddress("larinact_eta_tower", &larinact_eta_tower, &b_larinact_eta_tower);
+   fChain->SetBranchAddress("larinact_phi_module", &larinact_phi_module, &b_larinact_phi_module);
+   fChain->SetBranchAddress("larinact_EnergyEm", &larinact_EnergyEm, &b_larinact_EnergyEm);
+   fChain->SetBranchAddress("larinact_EnergyNonEm", &larinact_EnergyNonEm, &b_larinact_EnergyNonEm);
+   fChain->SetBranchAddress("larinact_EnergyInv", &larinact_EnergyInv, &b_larinact_EnergyInv);
+   fChain->SetBranchAddress("larinact_EnergyEsc", &larinact_EnergyEsc, &b_larinact_EnergyEsc);
+   fChain->SetBranchAddress("tileact_n", &tileact_n, &b_tileact_n);
+   fChain->SetBranchAddress("tileact_eta", &tileact_eta, &b_tileact_eta);
+   fChain->SetBranchAddress("tileact_phi", &tileact_phi, &b_tileact_phi);
+   fChain->SetBranchAddress("tileact_ID", &tileact_ID, &b_tileact_ID);
+   fChain->SetBranchAddress("tileact_EnergyTot", &tileact_EnergyTot, &b_tileact_EnergyTot);
+   fChain->SetBranchAddress("tileact_EnergyVis", &tileact_EnergyVis, &b_tileact_EnergyVis);
+   fChain->SetBranchAddress("tileact_detector", &tileact_detector, &b_tileact_detector);
+   fChain->SetBranchAddress("tileact_sample", &tileact_sample, &b_tileact_sample);
+   fChain->SetBranchAddress("tileact_side", &tileact_side, &b_tileact_side);
+   fChain->SetBranchAddress("tileact_reg_sec", &tileact_reg_sec, &b_tileact_reg_sec);
+   fChain->SetBranchAddress("tileact_eta_tower", &tileact_eta_tower, &b_tileact_eta_tower);
+   fChain->SetBranchAddress("tileact_phi_module", &tileact_phi_module, &b_tileact_phi_module);
+   fChain->SetBranchAddress("tileact_EnergyEm", &tileact_EnergyEm, &b_tileact_EnergyEm);
+   fChain->SetBranchAddress("tileact_EnergyNonEm", &tileact_EnergyNonEm, &b_tileact_EnergyNonEm);
+   fChain->SetBranchAddress("tileact_EnergyInv", &tileact_EnergyInv, &b_tileact_EnergyInv);
+   fChain->SetBranchAddress("tileact_EnergyEsc", &tileact_EnergyEsc, &b_tileact_EnergyEsc);
+   fChain->SetBranchAddress("tileinact_n", &tileinact_n, &b_tileinact_n);
+   fChain->SetBranchAddress("tileinact_eta", &tileinact_eta, &b_tileinact_eta);
+   fChain->SetBranchAddress("tileinact_phi", &tileinact_phi, &b_tileinact_phi);
+   fChain->SetBranchAddress("tileinact_ID", &tileinact_ID, &b_tileinact_ID);
+   fChain->SetBranchAddress("tileinact_EnergyTot", &tileinact_EnergyTot, &b_tileinact_EnergyTot);
+   fChain->SetBranchAddress("tileinact_EnergyVis", &tileinact_EnergyVis, &b_tileinact_EnergyVis);
+   fChain->SetBranchAddress("tileinact_detector", &tileinact_detector, &b_tileinact_detector);
+   fChain->SetBranchAddress("tileinact_sample", &tileinact_sample, &b_tileinact_sample);
+   fChain->SetBranchAddress("tileinact_side", &tileinact_side, &b_tileinact_side);
+   fChain->SetBranchAddress("tileinact_reg_sec", &tileinact_reg_sec, &b_tileinact_reg_sec);
+   fChain->SetBranchAddress("tileinact_eta_tower", &tileinact_eta_tower, &b_tileinact_eta_tower);
+   fChain->SetBranchAddress("tileinact_phi_module", &tileinact_phi_module, &b_tileinact_phi_module);
+   fChain->SetBranchAddress("tileinact_EnergyEm", &tileinact_EnergyEm, &b_tileinact_EnergyEm);
+   fChain->SetBranchAddress("tileinact_EnergyNonEm", &tileinact_EnergyNonEm, &b_tileinact_EnergyNonEm);
+   fChain->SetBranchAddress("tileinact_EnergyInv", &tileinact_EnergyInv, &b_tileinact_EnergyInv);
+   fChain->SetBranchAddress("tileinact_EnergyEsc", &tileinact_EnergyEsc, &b_tileinact_EnergyEsc);
+   fChain->SetBranchAddress("lardm_n", &lardm_n, &b_lardm_n);
+   fChain->SetBranchAddress("lardm_eta", &lardm_eta, &b_lardm_eta);
+   fChain->SetBranchAddress("lardm_phi", &lardm_phi, &b_lardm_phi);
+   fChain->SetBranchAddress("lardm_ID", &lardm_ID, &b_lardm_ID);
+   fChain->SetBranchAddress("lardm_EnergyTot", &lardm_EnergyTot, &b_lardm_EnergyTot);
+   fChain->SetBranchAddress("lardm_EnergyVis", &lardm_EnergyVis, &b_lardm_EnergyVis);
+   fChain->SetBranchAddress("lardm_detzside", &lardm_detzside, &b_lardm_detzside);
+   fChain->SetBranchAddress("lardm_type", &lardm_type, &b_lardm_type);
+   fChain->SetBranchAddress("lardm_sampling", &lardm_sampling, &b_lardm_sampling);
+   fChain->SetBranchAddress("lardm_region", &lardm_region, &b_lardm_region);
+   fChain->SetBranchAddress("lardm_ieta", &lardm_ieta, &b_lardm_ieta);
+   fChain->SetBranchAddress("lardm_iphi", &lardm_iphi, &b_lardm_iphi);
+   fChain->SetBranchAddress("lardm_EnergyEm", &lardm_EnergyEm, &b_lardm_EnergyEm);
+   fChain->SetBranchAddress("lardm_EnergyNonEm", &lardm_EnergyNonEm, &b_lardm_EnergyNonEm);
+   fChain->SetBranchAddress("lardm_EnergyInv", &lardm_EnergyInv, &b_lardm_EnergyInv);
+   fChain->SetBranchAddress("lardm_EnergyEsc", &lardm_EnergyEsc, &b_lardm_EnergyEsc);
+   fChain->SetBranchAddress("tiledm_n", &tiledm_n, &b_tiledm_n);
+   fChain->SetBranchAddress("tiledm_eta", &tiledm_eta, &b_tiledm_eta);
+   fChain->SetBranchAddress("tiledm_phi", &tiledm_phi, &b_tiledm_phi);
+   fChain->SetBranchAddress("tiledm_ID", &tiledm_ID, &b_tiledm_ID);
+   fChain->SetBranchAddress("tiledm_EnergyTot", &tiledm_EnergyTot, &b_tiledm_EnergyTot);
+   fChain->SetBranchAddress("tiledm_EnergyVis", &tiledm_EnergyVis, &b_tiledm_EnergyVis);
+   fChain->SetBranchAddress("tiledm_detzside", &tiledm_detzside, &b_tiledm_detzside);
+   fChain->SetBranchAddress("tiledm_type", &tiledm_type, &b_tiledm_type);
+   fChain->SetBranchAddress("tiledm_sampling", &tiledm_sampling, &b_tiledm_sampling);
+   fChain->SetBranchAddress("tiledm_region", &tiledm_region, &b_tiledm_region);
+   fChain->SetBranchAddress("tiledm_ieta", &tiledm_ieta, &b_tiledm_ieta);
+   fChain->SetBranchAddress("tiledm_iphi", &tiledm_iphi, &b_tiledm_iphi);
+   fChain->SetBranchAddress("tiledm_EnergyEm", &tiledm_EnergyEm, &b_tiledm_EnergyEm);
+   fChain->SetBranchAddress("tiledm_EnergyNonEm", &tiledm_EnergyNonEm, &b_tiledm_EnergyNonEm);
+   fChain->SetBranchAddress("tiledm_EnergyInv", &tiledm_EnergyInv, &b_tiledm_EnergyInv);
+   fChain->SetBranchAddress("tiledm_EnergyEsc", &tiledm_EnergyEsc, &b_tiledm_EnergyEsc);
+   fChain->SetBranchAddress("mcvtx_n", &mcvtx_n, &b_mcvtx_n);
+   fChain->SetBranchAddress("mcvtx_x", &mcvtx_x, &b_mcvtx_x);
+   fChain->SetBranchAddress("mcvtx_y", &mcvtx_y, &b_mcvtx_y);
+   fChain->SetBranchAddress("mcvtx_z", &mcvtx_z, &b_mcvtx_z);
+   fChain->SetBranchAddress("mcvtx_barcode", &mcvtx_barcode, &b_mcvtx_barcode);
+   fChain->SetBranchAddress("mcpart_n", &mcpart_n, &b_mcpart_n);
+   fChain->SetBranchAddress("mcpart_pt", &mcpart_pt, &b_mcpart_pt);
+   fChain->SetBranchAddress("mcpart_m", &mcpart_m, &b_mcpart_m);
+   fChain->SetBranchAddress("mcpart_eta", &mcpart_eta, &b_mcpart_eta);
+   fChain->SetBranchAddress("mcpart_phi", &mcpart_phi, &b_mcpart_phi);
+   fChain->SetBranchAddress("mcpart_type", &mcpart_type, &b_mcpart_type);
+   fChain->SetBranchAddress("mcpart_status", &mcpart_status, &b_mcpart_status);
+   fChain->SetBranchAddress("mcpart_barcode", &mcpart_barcode, &b_mcpart_barcode);
+   fChain->SetBranchAddress("mcpart_mothertype", &mcpart_mothertype, &b_mcpart_mothertype);
+   fChain->SetBranchAddress("mcpart_motherbarcode", &mcpart_motherbarcode, &b_mcpart_motherbarcode);
+   fChain->SetBranchAddress("mcpart_mcprodvtx_index", &mcpart_mcprodvtx_index, &b_mcpart_mcprodvtx_index);
+   fChain->SetBranchAddress("mcpart_mother_n", &mcpart_mother_n, &b_mcpart_mother_n);
+   fChain->SetBranchAddress("mcpart_mother_index", &mcpart_mother_index, &b_mcpart_mother_index);
+   fChain->SetBranchAddress("mcpart_mcdecayvtx_index", &mcpart_mcdecayvtx_index, &b_mcpart_mcdecayvtx_index);
+   fChain->SetBranchAddress("mcpart_child_n", &mcpart_child_n, &b_mcpart_child_n);
+   fChain->SetBranchAddress("mcpart_child_index", &mcpart_child_index, &b_mcpart_child_index);
+   Notify();
+}
+
+Bool_t ClCaloHits::Notify()
+{
+   // The Notify() function is called when a new file is opened. This
+   // can be either for a new TTree in a TChain or when when a new TTree
+   // is started when using PROOF. It is normally not necessary to make changes
+   // to the generated code, but the routine can be extended by the
+   // user if needed. The return value is currently not used.
+
+   return kTRUE;
+}
+
+void ClCaloHits::Show(Long64_t entry)
+{
+// Print contents of entry.
+// If entry is not specified, print current entry
+   if (!fChain) return;
+   fChain->Show(entry);
+}
+Int_t ClCaloHits::Cut(Long64_t entry)
+{
+// This function may be called from Loop.
+// returns  1 if entry is accepted.
+// returns -1 otherwise.
+   return 1;
+}
+#endif // #ifdef ClCaloHits_cxx
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/env.sh b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/env.sh
new file mode 100755
index 000000000000..603e5d1a42bb
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/env.sh
@@ -0,0 +1,11 @@
+#
+export ATLAS_LOCAL_ROOT_BASE=/cvmfs/atlas.cern.ch/repo/ATLASLocalRootBase 
+alias setupATLAS='source ${ATLAS_LOCAL_ROOT_BASE}/user/atlasLocalSetup.sh'
+setupATLAS
+#
+lsetup "root 6.18.04-x86_64-centos7-gcc8-opt"
+#
+echo " " 
+echo " === ROOTSYS: $ROOTSYS === "
+echo " "
+#
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/get_SF.C b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/get_SF.C
new file mode 100644
index 000000000000..0e9dffd3fae1
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/get_SF.C
@@ -0,0 +1,422 @@
+//
+//  Plot and Fit Pseudorapidity dependence of the Sampling Fraction
+//  with P0 and P1
+//  - for whole range and after different Eta selection 
+//
+//  Use P0 fit for selected region without tie-rods to get
+//  HEC sampling fractions (for front and rear wheels)
+//
+
+#include <iomanip>
+
+TCanvas *can;
+TH1     *frame;
+
+void get_SF (
+     Char_t const *G4Ver = "10.1",           //  Geant4 version
+     Char_t const *PhLis = "FTFP_BERT_ATL",  //  Physics list
+     Int_t   iPart = 0,    //  primary particle
+     Int_t   IEner = 100 ) //  particle energy [GeV]
+{
+//
+//
+   gStyle->SetLabelSize(0.040,"Y");
+//   
+   gROOT->ForceStyle();
+//
+//  -------------------------------------------------------------------- 
+//
+//
+   const Int_t NHEC = 2; //  HEC wheels
+   Char_t const *THEC[NHEC] = { "HEC front wheel", "HEC  rear wheel" };
+   Char_t const *CHEC[NHEC] = { "fwh", "rwh" };
+   Char_t const *SHEC[NHEC] = { "sc1", "sc2" };
+   Float_t  Cut1[NHEC] = { 1.65, 1.80 }; 
+   Float_t  Cut2[NHEC] = { 2.95, 3.05 }; 
+   Double_t SamFr[NHEC], eSamFr[NHEC];
+//
+//
+   const Int_t NPart = 4;
+   Char_t const *TPart[NPart] = { 
+      "Electrons+Positrons+Gammas", 
+      "Electron", "Positron", "Gamma" };
+   Char_t const *CPart[NPart] = { "EM", "ele", "pos", "gam" };   
+//
+//
+   const Int_t NVar = 1;
+   const Char_t *TVar[NVar] = { "Sampling Fraction" };
+   const Char_t *UVar[NVar] = { "SF = <E_{ACT}> /  <E_{TOT}>" };
+   const Char_t *CVar[NVar] = { "samf" };
+//
+//  --------------------------------------------------------------------
+// 
+//
+   TF1          *fuc1,  *fuc2;
+   TGraphErrors *graf1, *graf2;
+   TText        *txt;
+//
+//
+   Int_t    *N;
+   Float_t **X, **EX;
+   Float_t **Y, **EY;
+//
+//
+   Int_t    ib,it;
+   Int_t    iHEC, n,nb, ipad;
+   Int_t    inexp, inval;
+   Float_t  Sele1,Sele2, xmin,xmax, ymin,ymax,yb, dmin,dmax,db, xtxt,ytxt;
+   Float_t  A,B, eA,eB, X1,X2, sf,sf1,sf2, esf,esf1,esf2, r1,r2;
+   Double_t expon, powr,mant, round;
+   Char_t   gtit[999], utit[99], ytit[99];   
+   Char_t   Tsel[111],Csel[111], Cfile[999], Text[111];
+   Char_t   comm[111], psfil[111];
+   Char_t   datdir[999], datfil[999];
+//
+   FILE *fout;
+   ifstream fin;
+//
+//
+   Int_t NSlct = 3;
+   Int_t IVar  = 0;
+//
+   Char_t const *Xtit = "#eta";   
+//
+//  -------------------------------------------------------------------- 
+//
+//  Arguments
+//
+   if (iPart < 0 || iPart >= NPart) {
+      cout << " *** Bad iPart=" << iPart << " !!! *** " << endl;
+      exit (-2);
+   }
+//
+//
+   sprintf (ytit, "%s", UVar[IVar]);
+//
+//  Canvas
+//
+   can = new TCanvas("can", "mypage", 980,0,600,900);
+//
+   can->SetFillColor(10);
+   can->SetBorderSize(0);
+   can->SetBorderMode(0);
+//
+//
+   can->Print("fplot.ps[");   // No actual print, just open ps-file
+//
+//  -------------------------------------------------------------------- 
+//
+//  Loop over HEC wheels
+//
+   for (iHEC=0; iHEC<NHEC; iHEC++) {
+       N = new Int_t     [NSlct];
+       X = new Float_t * [NSlct];
+      EX = new Float_t * [NSlct];
+       Y = new Float_t * [NSlct];
+      EY = new Float_t * [NSlct];
+//
+//  Read in sampling fractions without selection
+//
+      it = 0;
+//
+      sprintf (Cfile, 
+               "%s-%s/%s_eta_%s-%s_e%04d_%s.dat",
+               G4Ver,PhLis, CVar[IVar], CHEC[iHEC], 
+	       CPart[iPart],IEner,SHEC[iHEC]);
+      cout << " ===  File:  " << Cfile << endl;		 
+//
+      Sele1 = Sele2 = 0;
+//
+      read_var_cutX ( Cfile, Sele1,Sele2, N[it], X[it],EX[it], Y[it],EY[it]);
+//
+//  Common X-boundaries
+//
+      n = N[it]-1;
+      xmin = X[it][0] - EX[it][0];
+      xmax = X[it][n] + EX[it][n];
+//
+//  Read in sampling fractions after boundary cuts
+//
+      it = 1;
+//   
+      Sele1 = Cut1[iHEC];
+      Sele2 = Cut2[iHEC];
+//
+      read_var_cutX ( Cfile, Sele1,Sele2, N[it], X[it],EX[it], Y[it],EY[it]);
+//
+//  Boundaries of the selected regions (not bin centers)
+//  where to calculate P1-results
+//
+      X1 = Cut1[iHEC];
+      X2 = Cut2[iHEC];
+//
+//  Select bins not affected by tie-rods (using data after boundary cuts)
+//
+      it = 2;
+//   
+      n = N[1];
+       X[it] = new Float_t [n];
+      EX[it] = new Float_t [n];
+       Y[it] = new Float_t [n];
+      EY[it] = new Float_t [n];
+//
+//
+      nb = 0;
+//
+      if (iHEC == 0) {
+         for (ib=0; ib<n; ib++) {
+      	    if (X[1][ib] > 1.55 && X[1][ib] < 1.60) continue;
+      	    if (X[1][ib] > 1.75 && X[1][ib] < 1.80) continue;
+      	    if (X[1][ib] > 2.05 && X[1][ib] < 2.10) continue;
+      	    if (X[1][ib] > 2.70 && X[1][ib] < 2.80) continue;
+//
+      	     X[it][nb] =  X[1][ib];	    
+      	    EX[it][nb] = EX[1][ib];	    
+      	     Y[it][nb] =  Y[1][ib];	    
+      	    EY[it][nb] = EY[1][ib];
+//
+      	    nb ++;	    
+         }
+      }
+//
+      else if (iHEC == 1) {
+         for (ib=0; ib<n; ib++) {
+      	    if (X[1][ib] > 1.90 && X[1][ib] < 2.00) continue;
+      	    if (X[1][ib] > 2.25 && X[1][ib] < 2.30) continue;
+      	    if (X[1][ib] > 2.85 && X[1][ib] < 3.00) continue;
+//
+      	     X[it][nb] =  X[1][ib];	    
+      	    EX[it][nb] = EX[1][ib];	    
+      	     Y[it][nb] =  Y[1][ib];	    
+      	    EY[it][nb] = EY[1][ib];
+//
+      	    nb ++;	    
+         }
+      }
+//
+      N[it] = nb;   
+//
+//  Prepare to draw
+//
+      sprintf (gtit, 
+               "G4.%s, %s:  %s %s VS Pseudorapidity (%d GeV %s)", 
+               G4Ver,PhLis, THEC[iHEC],TVar[IVar],  IEner,TPart[iPart]);
+      cout << " ===  Title: " << gtit << endl;
+//
+      can->Clear();
+//
+      mypage (gtit, 2,NSlct);      
+//
+      sprintf (Tsel,"%.2f <= #eta < %.2f",Cut1[iHEC],Cut2[iHEC]);
+//   
+//  Loop over jobs: no selection and after selection
+//
+      for (it=0; it<NSlct; it++) {
+         if (it == 0) 
+            sprintf (utit, "Whole Range");
+         else if (it == 1)	    
+            sprintf (utit, "Selected Region: %s", Tsel);
+         else
+            sprintf (utit, "Selected Region without Tie-Rods");
+// 
+//     
+         ymin =  1.e10;
+         ymax = -1.e10;
+//
+         for (ib=0; ib<N[it]; ib++) {
+            yb = Y[it][ib] - EY[it][ib];
+            if (yb < ymin)
+               ymin = yb;
+            yb = Y[it][ib] + EY[it][ib];
+            if (yb > ymax)
+               ymax = yb;
+         }	 
+//
+         yb = (ymax-ymin);
+         ymin = ymin - yb*0.2;
+         ymax = ymax + yb*0.5;      
+//
+//  Graph 1: plot and fit with constant
+//
+         ipad = it*2+1;
+         pad->cd(ipad);
+//
+         myplot (xmin,xmax, ymin,ymax, Xtit,ytit, utit);
+//      
+         graf1 = new TGraphErrors (N[it], X[it],Y[it], EX[it],EY[it]);
+         graf1->SetMarkerStyle(20);
+         graf1->SetMarkerColor(4);
+//
+         graf1->Draw("P");	 
+//
+//
+         fuc1 = new TF1 ("fuc1", "[0]", xmin,xmax);
+         fuc1->SetParNames("SF0");	
+         graf1->Fit("fuc1","EQ","",xmin,xmax);   
+//
+          sf = fuc1->GetParameter(0);      
+         esf = fuc1->GetParError(0); 
+//
+         if (it == 2) {
+             SamFr[iHEC] =  sf;
+            eSamFr[iHEC] = esf;		   
+         }
+//
+//  Deviation from the fit
+//      
+         dmin =  1.e10;
+         dmax = -1.e10;
+//
+         for (ib=0; ib<N[it]; ib++) {
+            db = Y[it][ib]/sf;
+            if (db < dmin)
+               dmin = db;
+            if (db > dmax)
+               dmax = db;
+         } 
+//
+         txt = new TText ();
+         txt->SetTextAlign(21);
+//
+         xtxt = 0.5*(xmin+xmax);
+         ytxt = ymin + (ymax-ymin)*0.05;
+         sprintf (Text, "%5.3f < SF/SF0 < %5.3f", dmin,dmax);
+//
+         txt->DrawText(xtxt,ytxt, Text);	
+//
+//  Graph 2: plot and linear fit
+//   
+         ipad = it*2+2;
+         pad->cd(ipad);
+//       
+      	 myplot (xmin,xmax, ymin,ymax, Xtit,ytit, utit);
+//       
+         graf2 = new TGraphErrors (N[it], X[it],Y[it], EX[it],EY[it]);
+         graf2->SetMarkerStyle(20);
+         graf2->SetMarkerColor(4);
+//       
+      	 graf2->Draw("P");	 
+//       
+//
+      	 fuc2 = new TF1 ("fuc2","[0]*x+[1]", xmin,xmax);
+         graf2->Fit("fuc2","EQ","",xmin,xmax);
+//       
+//  Values at boundaries
+//
+          A = fuc2->GetParameter(0);	   
+         eA = fuc2->GetParError(0);	  
+          B = fuc2->GetParameter(1);
+         eB = fuc2->GetParError(1);
+//
+         sf1 = A * X1 + B;
+         esf1 = TMath::Sqrt( X1*X1 * eA*eA + eB*eB );	     
+         sf2 = A * X2 + B;    
+         esf2 = TMath::Sqrt( X2*X2 * eA*eA + eB*eB );	     
+//
+//
+         r1 = sf1/sf;		    
+         r2 = sf2/sf;		    
+//
+         txt = new TText ();
+         txt->SetTextAlign(11);
+//
+         xtxt = xmin + 0.05*(xmax-xmin);
+         ytxt = ymin + (ymax-ymin)*0.15;
+         sprintf (Text, "SF(%.3f) = %8.6f   /SF0 = %5.3f", X1,sf1,r1);
+//
+         txt->DrawText(xtxt,ytxt, Text);	
+//
+//
+         ytxt = ymin + (ymax-ymin)*0.05;
+         sprintf (Text, "SF(%.3f) = %8.6f   /SF0 = %5.3f", X2,sf2,r2);
+//
+         txt->DrawText(xtxt,ytxt, Text);	     
+//
+      } //  over it
+//
+//
+      can->Update();
+      can->Print("fplot.ps");
+//
+//
+      delete     N;
+      delete []  X;	 
+      delete [] EX;	 
+      delete []  Y;	 
+      delete [] EY;	 
+//
+   } //  over iHEC
+//
+//  Save PS-file
+//
+   can->Print("fplot.ps]");   // No actual print, just close ps-file
+//
+   sprintf (psfil, "%s-%s/SF_%s_e%04d.ps",
+     	    G4Ver,PhLis, CPart[iPart],IEner);
+//
+   sprintf (comm, "mv -f fplot.ps %s", psfil);
+   gSystem->Exec( comm );
+//
+//  --------------------------------------------------------------------
+//
+//  Write HEC sampling fractions for both wheels
+//
+   if (iPart == 0 && IEner == 100) 
+      sprintf (datfil, "%s-%s/HEC_SF.txt", G4Ver,PhLis);
+   else {
+      sprintf (datfil, "%s-%s/HEC_SF_%s_e%04d.txt",
+     	       G4Ver,PhLis, CPart[iPart],IEner);
+   }
+//
+   fout = fopen (datfil,"w");  
+//   
+   for (iHEC=0; iHEC<NHEC; iHEC++) {
+//
+//  Make rounding of the sampling fraction 
+//  (keeping number of digits corresponding to its error)
+//
+      if (SamFr[iHEC] <= 0 || eSamFr[iHEC] <= 0) {
+         cerr << " *** Bad SamFr/eSamFr = " << SamFr[iHEC] << "/";
+	 cerr << eSamFr[iHEC] << " for " << THEC[iHEC] << " !!! *** " << endl;
+	 exit (-2);
+      }
+//
+//  Get first meaningful digit of the error
+//
+      expon = log10( eSamFr[iHEC] );	     
+//
+      if (expon >= 0)
+         inexp = int( expon );
+      else
+         inexp = int( expon ) - 1;
+      powr  = pow( 10, inexp );
+      mant  = eSamFr[iHEC]/powr;
+      inval = int( mant + 0.5 );
+      round = inval * powr;
+//
+      if (inval == 10) {
+         inval = 1;
+         inexp ++;
+      }
+//
+//  Round sampling fraction
+//
+      powr  = pow( 10, inexp );
+//
+      mant  = SamFr[iHEC]/powr;
+      inval = int( mant + 0.5 );
+      round = inval * powr;
+//
+//  Write sampling fraction in txt-file
+//
+      fprintf (fout, " %s :  %f +- %f  ->  %g \n", 
+               THEC[iHEC], SamFr[iHEC],eSamFr[iHEC], round);
+   }
+//
+   fclose (fout);
+//
+/* ------------------------------------------------------------------ */
+//
+//      
+}
+////////////////////////////////////////////////////////////////////////
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/makeclass_ClCaloHits.C b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/makeclass_ClCaloHits.C
new file mode 100644
index 000000000000..c1570945cd2e
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/makeclass_ClCaloHits.C
@@ -0,0 +1,33 @@
+void makeclass_ClCaloHits ( )
+{
+   TFile *f;
+   TTree *t;
+//   
+   Char_t const *TreeName  = "calohitD3PD"; //  name of the tree
+   Char_t const *ClassName = "ClCaloHits";  //  class name
+//
+//
+   Char_t const *Dir = "/eos/user/k/kiryunin/HEC_SF";
+   Char_t rootfile[222];
+   sprintf (rootfile, 
+            "%s/10.1-FTFP_BERT_ATL/DumpHits/ele_e0100_sc1/01.root", Dir);
+//
+//   
+   f = TFile::Open(rootfile);
+   if ( !f->IsOpen() ) {
+      printf (" *** File - %s - is not openned *** \n", rootfile);
+      exit (-2);
+   }
+//
+//
+   t = (TTree *) f->Get(TreeName);
+   if ( !t ) {
+      printf (" *** Tree - %s - is not found *** \n", TreeName);
+      exit (-2);
+   }
+//  
+// 
+   t->MakeClass(ClassName);
+//
+//   
+}
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/mypage.C b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/mypage.C
new file mode 100644
index 000000000000..d6f3fdbc6511
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/mypage.C
@@ -0,0 +1,32 @@
+////////////////////////////////////////////////////////////////////////
+#include "TCanvas.h"
+#include "TStyle.h"
+#include "TPaveText.h"
+#include "TPad.h"
+
+TPad *pad;
+
+void mypage (Char_t const *tt="Text to Test", Int_t ix=1, Int_t iy=1) 
+{
+//
+//  Page title
+//
+   TPaveText *title = new TPaveText(0.01,0.955,0.99,0.995,"NDC");
+   title->SetFillColor(10);
+   title->SetTextFont(72);
+   title->AddText(tt);
+   title->SetBorderSize(0);
+   title->Draw();
+//
+//  Main pad
+//
+   pad = new TPad("pad","MainPad", 0.005,0.005, 0.995,0.950, 10,0,0);   
+//
+//  Divide main pad
+//
+   pad->Divide(ix,iy, 0.0001, 0.0001, 10);
+   pad->Draw();
+//
+}
+////////////////////////////////////////////////////////////////////////
+
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/myplot.C b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/myplot.C
new file mode 100644
index 000000000000..3260930c703b
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/myplot.C
@@ -0,0 +1,51 @@
+////////////////////////////////////////////////////////////////////////
+void myplot (Float_t x1=0, Float_t x2=1, Float_t y1=0, Float_t y2=1,
+            Char_t const *titx="X axis", Char_t const *tity="Y axis", 
+            Char_t const *titp=" ")
+{
+   TH1F *frame = gPad->DrawFrame(x1,y1,x2,y2,titp);
+//           
+   frame->GetXaxis()->SetTitleSize(0.05);
+   frame->GetXaxis()->SetTitleOffset(1.2);
+   frame->GetXaxis()->SetTitle(titx);
+//
+   frame->GetYaxis()->SetTitleSize(0.05);
+   frame->GetYaxis()->SetTitleOffset(1.8);
+   frame->GetYaxis()->SetTitle(tity);
+//
+//
+//
+   TPaveText *title = 0;
+   TObject *obj;
+   TIter next(gPad->GetListOfPrimitives());
+   while ((obj = next())) {
+      if (!obj->InheritsFrom(TPaveText::Class())) continue;
+      title = (TPaveText*)obj;
+      if (strcmp(title->GetName(),"title")) {title = 0; continue;}
+      break;
+   }
+//
+//
+   if (title) {
+      title->SetBorderSize(1);
+      title->SetFillColor(10);
+      title->SetTextAlign(22);
+      title->SetTextFont(62);
+//
+      title->SetX1NDC(0.18);      
+      title->SetX2NDC(0.98);      
+//
+      TText *t0 = (TText*)title->GetLine(0);
+//
+      t0->SetTextColor( 1 );
+//      
+      title->Draw();
+   }
+   else {
+      cerr << " *** myplot:  title does not exist !!! *** " << endl;
+      return;
+   }
+//
+//
+}
+////////////////////////////////////////////////////////////////////////
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/read_var_cutX.C b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/read_var_cutX.C
new file mode 100644
index 000000000000..ab217c929e1d
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/read_var_cutX.C
@@ -0,0 +1,121 @@
+#include <iostream>
+#include <fstream>
+
+using namespace std;
+
+////////////////////////////////////////////////////////////////////////
+//
+//  Read average values and apply given cuts (if necessary)
+//
+void read_var_cutX  ( 
+     Char_t  *file, 
+     Float_t X1, Float_t X2,
+     Int_t   &N, 
+     Float_t *&X, Float_t *&EX,
+     Float_t *&Y, Float_t *&EY )
+{
+//
+   Int_t   *nevts;
+   Float_t *x, *ave,*err;
+//
+//
+   Int_t   i;
+   Int_t   n, nx;
+   Float_t xmin,xmax, dx;
+//
+   ifstream fin;
+//   
+/* ------------------------------------------------------------------ */   
+//
+//  Open file
+//
+   fin.open(file);
+   if ( !fin.good() ) {
+      cerr << " *** read_var_cutX:  No file - " << file;
+      cerr << " !!! *** " << endl;
+      exit (-2);
+   }
+//
+//
+   fin >> n;
+   if (n < 1 || !fin.good()) {
+      cerr << " *** read_var_cutX:  Bad n = " << n << " in file - " << file;
+      cerr << " !!! *** " << endl;
+      exit (-2);
+   }
+//
+   fin >> xmin >> xmax >> dx;
+   if ( !fin.good() ) {
+      cerr << " *** read_var_cutX:  Bad 2-nd line in file - " << file;
+      cerr << " !!! *** " << endl;
+      exit (-2);
+   }
+//
+//  Read in data
+//
+   nevts = new Int_t   [n];
+   x     = new Float_t [n];   
+   ave   = new Float_t [n];   
+   err   = new Float_t [n];   
+//
+   N = 0;
+//
+   for (i=0; i<n; i++) {
+      fin >> x[i] >> nevts[i] >> ave[i] >> err[i];
+      if ( !fin.good() ) {
+   	 cerr << " *** read_var_cutX:  Bad reading for i=" << i;
+	 cerr << " in file - " << file << " !!! *** " << endl;
+   	 exit (-2);
+      }
+//
+      if (nevts[i] == 0) continue;
+      if (X1<X2 && ( x[i]<=X1 || x[i]>=X2 )) continue;
+//
+      N ++;
+   }
+//
+   fin.close();
+//
+//
+   if (N < 1) {
+      cerr << " *** read_var_cutX:  Bad N = " << N << " for file - " << file;
+      cerr << " !!! *** " << endl;
+      exit (-2);
+   }
+//
+//  Store non-zero data after cuts
+//
+    X = new Float_t [N]; 
+   EX = new Float_t [N]; 
+    Y = new Float_t [N];   
+   EY = new Float_t [N];   
+//
+   nx = N;
+   N  = 0;
+//
+   for (i=0; i<n; i++) {
+      if (nevts[i] == 0) continue;
+      if (X1<X2 && ( x[i]<=X1 || x[i]>=X2 )) continue;
+//
+       X[N] = x[i];
+      EX[N] = 0.5 * dx;
+       Y[N] = ave[i];
+      EY[N] = err[i];    
+//
+      N ++;        
+   }
+//
+   if (N != nx) {
+      cerr << " *** read_var_cutX:  Bad N != nx : " << N << " " << nx;
+      cerr << " !!! *** " << endl;
+      exit (-2);
+   }
+//
+//
+   delete [] nevts;
+   delete [] x;
+   delete [] ave;
+   delete [] err;
+//      
+}
+////////////////////////////////////////////////////////////////////////
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/rootlogon.C b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/rootlogon.C
new file mode 100644
index 000000000000..8a0eeb8ca752
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/rootlogon.C
@@ -0,0 +1,60 @@
+////////////////////////////////////////////////////////////////////////
+{
+   gStyle->SetHistFillColor(11);   
+   gStyle->SetHistFillStyle(1001); 
+   gStyle->SetOptStat(111110); 
+   gStyle->SetOptFit(111);
+//
+   gStyle->SetPadTickX(1);
+   gStyle->SetPadTickY(1);
+//
+   gStyle->SetNdivisions(507,"X");
+   gStyle->SetNdivisions(507,"Y");
+//
+   gStyle->SetLabelSize(0.045,"X");
+   gStyle->SetLabelSize(0.045,"Y");
+   gStyle->SetLabelOffset(0.01,"X");
+   gStyle->SetLabelOffset(0.01,"Y");
+//
+//
+   gStyle->SetTitleBorderSize(1);  
+//   
+//    
+   gStyle->SetStatBorderSize(1);
+   gStyle->SetStatStyle(0);     
+   gStyle->SetStatX(0.98);   
+   gStyle->SetStatY(0.90);   
+   gStyle->SetStatW(0.30);   
+   gStyle->SetStatH(0.20);   
+//
+   gStyle->SetTitleX(0.18);   
+   gStyle->SetTitleY(0.98);   
+   gStyle->SetTitleW(0.80);   
+   gStyle->SetTitleH(0.08);   
+//
+   gStyle->SetPadBorderSize(0);
+   gStyle->SetPadBorderMode(0);
+//
+   gStyle->SetPadLeftMargin(0.18);
+   gStyle->SetPadRightMargin(0.02);
+   gStyle->SetPadTopMargin(0.10);
+   gStyle->SetPadBottomMargin(0.15);
+//   
+   gStyle->SetMarkerSize(1.0);   
+//
+//
+   gErrorIgnoreLevel = kWarning;
+//
+//
+   gROOT->ForceStyle();
+//
+//
+   gROOT->ProcessLine(".L mypage.C");
+   gROOT->ProcessLine(".L myplot.C");
+//
+   gROOT->ProcessLine(".L samfr_etaX.C+");
+   gROOT->ProcessLine(".L read_var_cutX.C+");
+//
+//
+}
+////////////////////////////////////////////////////////////////////////
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/samfr_etaX.C b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/samfr_etaX.C
new file mode 100644
index 000000000000..09ae4230c95b
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/samfr_etaX.C
@@ -0,0 +1,605 @@
+//
+//  Get average HEC energies and sampling fractions
+//  as functions of Scan parameter, namely:  Pseudorapidity
+//
+//  To execute:
+//  .L samfr_etaX.C+
+//  samfr_etaX()
+//
+#include <iomanip>
+#include <iostream>
+#include <fstream>
+
+#include "TBranch.h"
+#include "TChain.h"
+#include "TMath.h"
+#include "TROOT.h"
+#include "TStyle.h"
+#include "TSystem.h"
+
+#include "ClCaloHits.h"
+#include "ClCaloHits.C"
+
+#include "substa.C"
+
+void samfr_etaX (
+     Int_t   iHEC  = 0,   //  HEC wheel (0 - front, 1 - rear)
+     Int_t   iPart = 0,   //  primary particle
+     Int_t   IEner = 100, //  particle energy [GeV]
+     Char_t const *G4Ver = "10.1",          //  Geant4 version
+     Char_t const *PhLis = "FTFP_BERT_ATL", //  Physics list
+     Char_t const *Dir   = "/eos/user/k/kiryunin/HEC_SF")
+{
+//  
+/* ------------------------------------------------------------------ */
+//
+//
+   const Int_t NHEC = 2; //  HEC wheels
+   Float_t ZHEC[NHEC] = { 431.95, 517.50 };
+   Char_t const *THEC[NHEC] = { "HEC front wheel", "HEC rear wheel" };
+   Char_t const *CHEC[NHEC] = { "fwh", "rwh" };
+   Char_t const *SHEC[NHEC] = { "sc1", "sc2" };
+//
+   const Int_t NPart = 4;
+   Int_t PPart[NPart] = { 0,  11, -11, 22 };
+   Char_t const *TPart[NPart] = { 
+      "EM-particle", "Electron", "Positron", "Gamma" };
+   Char_t const *CPart[NPart] = { "EM", "ele", "pos", "gam" };   
+//
+//
+   const Int_t NEHit = 6;
+   Char_t const *TEHit[NEHit] = { "Total energy", "Visible energy",
+      "EM energy", "Non-EM energy", "Invisible energy", "Escaped energy" };
+   Double_t EHit[NEHit];
+//
+   const Int_t NMat = 2;
+   Char_t const *TMat[NMat] = { "Active material", "Passive material" };
+//
+//
+   const Int_t NVar = 5;
+   Char_t const *TVar[NVar] = { 
+      "Visible energy in active material", 
+      "Visible energy in passive material",
+      "Total visible energy",
+      "Ratio of energies: Active-to-Total",
+      "Sampling fraction" };
+   Char_t const *CVar[NVar] = { "eact", "epas", "etot", "rate", "samf" };
+   Double_t AVE[NVar], ERR[NVar];
+   Double_t *Var[NVar], *eVar[NVar];
+//
+   const Int_t NVar1 = 4; //  NVar-1;
+   Double_t **Escan[NVar1];   
+//    
+//  --------------------------------------------------------------------
+//
+//
+   TChain *chain;
+//
+   ClCaloHits *clcalohits;
+//
+   Int_t *NevtSc;
+   Int_t *NevtZE;
+//   
+//
+   Int_t    ia,ie,ih,im,in,ip,is,iv,iw;
+   Int_t    nentries, nvtx,npart, Nhits, NevtOut, n1,n2, Typ;
+   UInt_t   Det,Sam,Sid;
+   Double_t DScan, Sc;
+   Double_t diff,diffa,diffr, esum, rms, rat,ert, d1,d2;
+   Char_t   JOB1[111],JOB2[111], job[111];
+   Char_t   rootfile[222];
+   Char_t   datdir[222], datfil[222];
+   Char_t   comm[222];
+//
+   FILE *fout;
+//     
+   //  Pseudorapidity -scan
+   Int_t    Nscan = 36;
+   Double_t Scan1 = 1.50; 
+   Double_t Scan2 = 3.30;
+// 
+   Int_t NSamp = 25; //  number of root-files with DumpHits
+//
+   Int_t IVar = 4;
+//
+   Double_t dCut = 1.e-5;
+   Double_t MeV2GeV = 0.001;
+//
+/* ------------------------------------------------------------------ */
+//
+//  Argument(s)
+//
+   if (iHEC<0 || iHEC>=NHEC) {
+      cerr << " *** Bad argument iHEC = " << iHEC << " !!! *** " << endl;
+      exit (-2);
+   }
+//
+   if (iPart<0 || iPart>=NPart) {
+      cerr << " *** Bad argument iPart = " << iPart << " !!! *** " << endl;
+      exit (-2);
+   }
+//
+//
+   sprintf (JOB1, "%s-%s", G4Ver,PhLis); 
+   sprintf (JOB2, "%s_e%04d_%s", CPart[iPart],IEner,SHEC[iHEC]); 
+//
+   cout << endl;
+   cout << " ===  " << JOB1 << " " << JOB2;
+   cout << ":  Pseudorapidity scan over " << THEC[iHEC] << "  === " << endl;     
+//
+/* ------------------------------------------------------------------ */
+//      
+//  Create chain from ROOT files
+//
+   chain = new TChain ("calohitD3PD");
+//
+//  Electrons + Positrons + Gammas
+//   
+   if (iPart == 0) {
+      for (ip=1; ip<=3; ip++) {
+         sprintf (job, "%s_e%04d_%s", CPart[ip],IEner,SHEC[iHEC]); 
+//
+         for (ia=1; ia<=NSamp; ia++) {
+            sprintf (rootfile, 
+	             "%s/%s/DumpHits/%s/%02d.root", 
+		     Dir,JOB1, job,ia);
+            if (ia == 1)
+               cout << " ===  " << rootfile << "  === " << endl;  
+            chain->Add(rootfile);
+	 }   
+      }      
+   }
+//
+//  Single particle
+//
+   else {
+      sprintf (job, "%s_e%04d_%s", CPart[iPart],IEner,SHEC[iHEC]); 
+//
+      for (ia=1; ia<=NSamp; ia++) {
+         sprintf (rootfile, 
+	          "%s/%s/DumpHits/%s/%02d.root", 
+		  Dir,JOB1, job,ia);
+         if (ia == 1)
+            cout << " ===  " << rootfile << "  === " << endl;  
+         chain->Add(rootfile);
+      }      
+   }
+//
+//
+   nentries = chain->GetEntries();
+   cout << " ===  Number of entries - " << nentries << "  === " << endl;
+   if (nentries < 100) {
+      cerr << " *** Two few nentries: " << nentries << " !!! *** "  << endl;
+      exit (-2);
+   }
+//
+   clcalohits = new ClCaloHits(chain);
+//
+//  Initialisation
+//
+   NevtOut = 0; //  number of events with Scan parameter out of boundaries
+//
+   DScan = (Scan2-Scan1) / double( Nscan );
+//      
+   NevtSc = new Int_t [Nscan]; //  number of stored events 
+   NevtZE = new Int_t [Nscan]; //  number of events with zero energy
+//
+   for (is=0; is<Nscan; is++) {
+      NevtSc[is] = 0;
+      NevtZE[is] = 0;
+   }   
+//   
+   for (iv=0; iv<NVar1; iv++) {
+      Escan[iv] = new Double_t * [Nscan];
+      for (is=0; is<Nscan; is++) {
+   	 Escan[iv][is] = new Double_t [nentries];
+   	 for (in=0; in<nentries; in++)
+   	    Escan[iv][is][in] = 0;
+      }  
+   }
+//
+/* ------------------------------------------------------------------ */
+//
+//  Loop over events
+//
+   for (in=0; in<nentries; in++) {
+      chain->GetEntry(in);
+//
+      nvtx = clcalohits->mcvtx_n;  
+      if (nvtx < 1) {
+         cerr << " *** Bad nvtx = " << nvtx;
+	 cerr << " for in=" << in << " !!! *** " << endl;
+         exit (-2);
+      }
+//
+      npart = clcalohits->mcpart_n;  
+      if (npart < 1) {
+         cerr << " *** Bad npart = " << npart;
+         cerr << " for in=" << in << " !!! *** " << endl;
+         exit (-2);
+      }
+//
+      Typ = (*clcalohits->mcpart_type)[0];
+      if (PPart[iPart] != 0 && Typ != PPart[iPart]) {
+         cerr << " *** Bad primary particles Typ = " << Typ;
+         cerr << " for PPart=" << PPart[iPart];
+         cerr << " for in=" << in << " !!! *** " << endl;
+         exit (-2);
+      }
+//
+//
+      Sc = double( (*clcalohits->mcpart_eta)[0] );
+//     
+      is = int( (Sc-Scan1) / DScan );
+//
+      if (is < 0 || is >= Nscan) {
+         NevtOut ++;
+	 continue;
+      }
+//
+//  --------------------------------------------------------------------
+//
+//  Check active LAr hits
+//      
+      iv = im = 0;
+      Nhits = clcalohits->laract_n;
+//
+      if (Nhits > 0) {
+         for (ih=0; ih<Nhits; ih++) {
+            //  any hit
+            EHit[0] = double( (*clcalohits->laract_EnergyTot)[ih] );
+            EHit[1] = double( (*clcalohits->laract_EnergyVis)[ih] );
+            EHit[2] = double( (*clcalohits->laract_EnergyEm)[ih] );
+            EHit[3] = double( (*clcalohits->laract_EnergyNonEm)[ih] );
+            EHit[4] = double( (*clcalohits->laract_EnergyInv)[ih] );
+            EHit[5] = double( (*clcalohits->laract_EnergyEsc)[ih] );
+//
+            //  total energy
+            if ( EHit[0] == 0) {
+	       cout << " xxx  Zero LAr hit in " << TMat[im];
+	       cout << ":  " << TEHit[0];
+	       cout << " = " << EHit[0];
+	       cout << " for in/ih=" << in << "/" << ih << "  xxx " << endl;
+	    }
+//
+            //  no negative energies (except invisible and, hence, total)
+            for (ie=0; ie<NEHit; ie++) {
+               if (ie==0 || ie==4) continue;
+//	       
+               if ( EHit[ie] < 0) {
+	          cerr << " *** Bad LAr hit in " << TMat[im];
+		  cerr << ":  " << TEHit[ie];
+		  cerr << " = " << EHit[ie];
+	          cerr << " for in/ih=" << in << "/" << ih;
+		  cerr << " !!! *** " << endl;
+                  exit (-2);
+	       }
+	    }
+//
+            //  sum for total energy
+            diff = EHit[0] - (EHit[1] + EHit[4]+EHit[5]);
+            if (EHit[0] != 0) 
+               diff = TMath::Abs( diff / EHit[0] );
+	    else   
+               diff = TMath::Abs( diff );
+	    if (diff > dCut) {
+	       cerr << " *** Bad LAr hit in " << TMat[im];
+	       cerr << ":  Sum for Etot - diff = ";
+	       cerr << diff << ", Etot = " << EHit[0];
+	       cerr << " for in/ih=" << in << "/" << ih << " !!! *** " << endl;
+               exit (-2);
+	    }	    
+//
+            //  sum for visible energy
+            diff = EHit[1] - (EHit[2]+EHit[3]);
+	    if (EHit[1] != 0)
+               diff = TMath::Abs( diff ) / EHit[1];
+            else	       
+               diff = TMath::Abs( diff );
+	    if (diff > dCut) {
+	       cerr << " *** Bad LAr hit in " << TMat[im];
+	       cerr << ":  Sum for Evis - diff = ";
+	       cerr << diff << ", Evis = " << EHit[1];
+	       cerr << " for in/ih=" << in << "/" << ih << " !!! *** " << endl;
+               exit (-2);
+	    }	    
+//
+//  HEC active hits
+//
+	    Det = (*clcalohits->laract_detector)[ih];
+	    Sam = (*clcalohits->laract_sample)[ih];
+	    Sid = (*clcalohits->laract_side)[ih];
+//
+            if (Det==1 && Sid==2) {
+	       if (Sam < 8 || Sam > 11) {
+	          cerr << " *** Bad HEC sample = " << Sam;
+	          cerr << " in " << TMat[im];
+	          cerr << " for in/ih=" << in << "/" << ih;
+	          cerr << " !!! *** " << endl;
+                  exit (-2);
+	       }
+//
+               iw = int( Sam/2 ) - 4;
+	       if (iw < 0 || iw >= NHEC) {
+	          cerr << " *** Bad HEC wheel = " << iw;
+		  cerr << " for Sam=" << Sam;
+	          cerr << " in " << TMat[im];
+	          cerr << " for in/ih=" << in << "/" << ih;
+	          cerr << " !!! *** " << endl;
+                  exit (-2);
+	       }
+//
+               if (iw == iHEC)
+                  Escan[iv][is][ NevtSc[is] ] += EHit[1];
+	    }
+//	    	    
+         } //  over hits
+      } //  Nhits > 0
+//
+//  --------------------------------------------------------------------
+//
+//  Check inactive (i.e. passive) LAr hits
+//      
+      iv = im = 1;
+      Nhits = clcalohits->larinact_n;
+//
+      if (Nhits > 0) {
+         for (ih=0; ih<Nhits; ih++) {
+            //  any hit
+            EHit[0] = double( (*clcalohits->larinact_EnergyTot)[ih] );
+            EHit[1] = double( (*clcalohits->larinact_EnergyVis)[ih] );
+            EHit[2] = double( (*clcalohits->larinact_EnergyEm)[ih] );
+            EHit[3] = double( (*clcalohits->larinact_EnergyNonEm)[ih] );
+            EHit[4] = double( (*clcalohits->larinact_EnergyInv)[ih] );
+            EHit[5] = double( (*clcalohits->larinact_EnergyEsc)[ih] );
+//
+            //  total energy
+            if ( EHit[0] == 0) {
+	       cout << " xxx  Zero LAr hit in " << TMat[im];
+	       cout << ":  " << TEHit[0];
+	       cout << " = " << EHit[0];
+	       cout << " for in/ih=" << in << "/" << ih << "  xxx " << endl;
+	    }
+//
+            //  no negative energies (except invisible and, hence, total)
+            for (ie=0; ie<NEHit; ie++) {
+               if (ie==0 || ie==4) continue;
+//	       
+               if ( EHit[ie] < 0) {
+	          cerr << " *** Bad LAr hit in " << TMat[im];
+		  cerr << ":  " << TEHit[ie];
+		  cerr << " = " << EHit[ie];
+	          cerr << " for in/ih=" << in << "/" << ih;
+		  cerr << " !!! *** " << endl;
+                  exit (-2);
+	       }
+	    }
+//
+            //  sum for total energy
+            diffa = TMath::Abs( EHit[0] - (EHit[1] + EHit[4]+EHit[5]) );
+	    if (EHit[0] != 0)
+               diffr = diffa / TMath::Abs( EHit[0] );
+            else
+               diffr = 0;
+	    if (diffa > dCut && diffr > dCut) {
+	       cerr << " *** Bad LAr hit in " << TMat[im];
+	       cerr << ":  Sum for Etot ";
+	       cerr << "- diffa = " << diffa;
+	       cerr << ", diffr = " << diffr;
+	       cerr << ", Etot = " << EHit[0];
+	       cerr << " for in/ih=" << in << "/" << ih << " !!! *** " << endl;
+               exit (-2);
+	    }	    
+//
+            //  sum for visible energy
+            diff = EHit[1] - (EHit[2]+EHit[3]);
+	    if (EHit[1] != 0)
+               diff = TMath::Abs( diff ) / EHit[1];
+            else	       
+               diff = TMath::Abs( diff );
+	    if (diff > dCut) {
+	       cerr << " *** Bad LAr hit in " << TMat[im];
+	       cerr << ":  Sum for Evis - diff = ";
+	       cerr << diff << ", Evis = " << EHit[1];
+	       cerr << " for in/ih=" << in << "/" << ih << " !!! *** " << endl;
+               exit (-2);
+	    }	    
+//
+//  HEC passive hits
+//
+	    Det = (*clcalohits->larinact_detector)[ih];
+	    Sam = (*clcalohits->larinact_sample)[ih];
+	    Sid = (*clcalohits->larinact_side)[ih];
+//
+            if (Det==1 && Sid==2) {
+	       if (Sam < 8 || Sam > 11) {
+	          cerr << " *** Bad HEC sample = " << Sam;
+	          cerr << " in " << TMat[im];
+	          cerr << " for in/ih=" << in << "/" << ih;
+	          cerr << " !!! *** " << endl;
+                  exit (-2);
+	       }
+//
+               iw = int( Sam/2 ) - 4;
+	       if (iw < 0 || iw >= NHEC) {
+	          cerr << " *** Bad HEC wheel = " << iw;
+		  cerr << " for Sam=" << Sam;
+	          cerr << " in " << TMat[im];
+	          cerr << " for in/ih=" << in << "/" << ih;
+	          cerr << " !!! *** " << endl;
+                  exit (-2);
+	       }
+//
+               if (iw == iHEC)
+                  Escan[iv][is][ NevtSc[is] ] += EHit[1];
+	    }
+//	    	    
+         } //  over hits
+      } //  Nhits > 0
+//
+//  --------------------------------------------------------------------
+//
+//  Skip event with zero total energy
+//
+      esum = Escan[0][is][ NevtSc[is] ] + 
+      	     Escan[1][is][ NevtSc[is] ];
+//
+      if (esum == 0) {
+         Escan[0][is][ NevtSc[is] ] = 0;
+         Escan[1][is][ NevtSc[is] ] = 0;
+         NevtZE[is] ++;
+      }
+      else
+         NevtSc[is] ++;
+//
+   } //  over events
+//
+   if (NevtOut > 0) {
+      cout << " xxx  " << NevtOut << " event(s)";
+      cout << " out of scan boundaries  xxx " << endl;
+   }
+//
+//
+   chain->Reset();   
+//
+//  --------------------------------------------------------------------
+//
+//  Get total visible energy and ratio of energies event-by-event
+//
+   n1 = n2 = 0;
+//   
+   for (is=0; is<Nscan; is++) {
+      n1 += NevtSc[is];
+      n2 += NevtZE[is];
+//
+      if (NevtZE[is] != 0) {
+         Sc = Scan1 + 0.5*DScan + is*DScan;
+         cout << " xxx  Non-zero NevtZE = " << NevtZE[is];
+	 cout << " for is/Sc = " << is << "/" << Sc;
+	 cout << " (NevtSc=" << NevtSc[is] << ")  xxx " << endl;
+      }      
+//      
+      if (NevtSc[is] < 1) continue;
+// 	   
+      for (in=0; in<NevtSc[is]; in++) {
+   	 iv = 2;
+   	 Escan[iv][is][in] = Escan[0][is][in] + 
+   			     Escan[1][is][in];
+//
+   	 if (Escan[iv][is][in] == 0) {
+   	    cerr << " *** Bad total visible energy";
+   	    cerr << " = " << Escan[iv][is][in];
+   	    cerr << " for is=" << is;
+   	    cerr << " in non-empty event nr.=" << in << " !!! *** " << endl;
+   	    exit (-2);
+   	 }
+//
+   	 iv = 3;
+   	 Escan[iv][is][in] = Escan[0][is][in] / 
+   			     Escan[2][is][in];
+      }
+   } //  over is      
+//
+//
+   if (n1+n2+NevtOut != nentries) {
+      cerr << " *** Bad number of events n1/n2/NevtOut/nentries = ";
+      cerr << n1 << "/" << n2 << "/" << NevtOut << "/";
+      cerr << nentries << " !!! *** " << endl;
+      exit (-2);
+   }
+//
+/* ------------------------------------------------------------------ */
+//
+//  Get Mean and Error of stored variables
+//
+   for (iv=0; iv<NVar; iv++) {
+       Var[iv] = new Double_t [Nscan];
+      eVar[iv] = new Double_t [Nscan]; 
+//
+      for (is=0; is<Nscan; is++) {
+     	  Var[iv][is] = 0;
+     	 eVar[iv][is] = 0;
+      }           
+   }
+//
+//   
+   for (is=0; is<Nscan; is++) {
+      if (NevtSc[is] < 1) continue;
+//
+      for (iv=0; iv<NVar; iv++) {
+   	 if (iv < NVar1) {   
+   	    substa (Escan[iv][is],NevtSc[is], AVE[iv],ERR[iv],rms);
+     	 }
+     	 else {   
+//
+   	    if (AVE[2] != 0) {
+   	       rat = AVE[0] / AVE[2];
+   	       if (AVE[0] != 0) 
+     		  d1 = ERR[0] / AVE[0];
+   	       else
+     		  d1 = 0;
+   	       d2 = ERR[2] / AVE[2];
+   	       ert = rat * TMath::Sqrt( d1*d1+d2*d2 );
+   	    }  
+   	    else {
+   	       rat = 0;
+   	       ert = 0;
+   	    }
+//     
+   	    AVE[NVar1] = rat;
+     	    ERR[NVar1] = ert;		 
+   	 }
+//
+//
+         if (iv > 2) {
+             Var[iv][is] = AVE[iv];
+	    eVar[iv][is] = ERR[iv];
+         }
+	 else { 
+             Var[iv][is] = AVE[iv] * MeV2GeV;
+	    eVar[iv][is] = ERR[iv] * MeV2GeV;
+	 }   
+//
+      } //  over iv
+   } //  over is
+//
+//  --------------------------------------------------------------------
+//
+//  Write [average energies and] sampling fractions VS scan parameter
+//
+   sprintf (datdir, "%s", JOB1);
+// 
+   sprintf (comm, "mkdir -p %s", datdir);
+   gSystem->Exec(comm);
+//
+//
+   for (iv=0; iv<NVar; iv++) {
+      if (iv != IVar) continue;
+//
+      cout << " ===  Written only:  " << TVar[IVar];
+      cout << " as a function of Pseudorapidity  === " << endl;
+//      
+      sprintf (datfil, "%s/%s_eta_%s-%s.dat", 
+               datdir, CVar[iv], CHEC[iHEC], JOB2);
+      fout = fopen (datfil,"w");  
+//
+      fprintf (fout, "%5d \n", Nscan);
+      fprintf (fout, "%10.3f%10.3f%10.3f \n", Scan1,Scan2,DScan); 
+//
+      for (is=0; is<Nscan; is++) {
+         Sc = Scan1 + 0.5*DScan + is*DScan;
+         fprintf (fout, "%10.3f %8d", Sc, NevtSc[is]);
+         fprintf (fout, " %15.7f %15.7f \n", Var[iv][is],eVar[iv][is]);
+      }
+//
+      fclose (fout);
+//
+   } //  over iv
+//
+/* ------------------------------------------------------------------ */
+//
+//  End
+//
+   return;
+}
+////////////////////////////////////////////////////////////////////////
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/store_eta.C b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/store_eta.C
new file mode 100644
index 000000000000..950379af24e7
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/store_eta.C
@@ -0,0 +1,24 @@
+//
+//  Loop over HEC wheels to store 
+//  Sampling fractions as functions of Pseudorapidity
+//  (for 100 GeV EM particles)
+//
+
+void store_eta (
+     Char_t const *G4Ver = "10.1",          //  Geant4 version
+     Char_t const *PhLis = "FTFP_BERT_ATL", //  Physics list
+     Char_t const *Dir   = "/eos/user/k/kiryunin/HEC_SF")
+{
+//   
+//
+   Int_t NHEC  = 2;
+//
+   Int_t iPart = 0;
+   Int_t IEner = 100;   
+//   
+   for (Int_t iHEC=0; iHEC<NHEC; iHEC++)
+         samfr_etaX ( iHEC, iPart,IEner, G4Ver,PhLis, Dir );
+//   
+//
+}
+////////////////////////////////////////////////////////////////////////
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/substa.C b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/substa.C
new file mode 100644
index 000000000000..6cbc9778a637
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Analysis/substa.C
@@ -0,0 +1,70 @@
+// =====================================================================
+// 
+//  Determine mean value, error of the mean and RMS 
+//  for the given set of values
+//
+//  Inputs     : arr    - array with data
+//               narr   - dimension of the array
+//
+//  Outputs    : ave  - mean value
+//               err  - error of the mean
+//               rms  - RMS
+//              
+// =====================================================================
+
+#include "TMath.h"
+
+void substa (
+     Double_t *arr, 
+     Int_t     narr, 
+     Double_t &ave, 
+     Double_t &err, 
+     Double_t &rms )
+{
+//
+   Int_t i;
+//   
+   Double_t s0=0, s1=0, s2=0; 
+//
+/* ------------------------------------------------------------------ */
+//
+   ave=0; err=0; rms=0;
+//   
+   if (narr < 1) {
+      cout << " *** substa:  Bad narr = " << narr << endl;
+      return;
+   }
+//
+//
+   for (i=0; i<narr; i++) {
+      s0 += 1.0;
+      s1 += arr[i];
+      s2 += arr[i]*arr[i];
+    }
+//
+    if (s0 == 0) {
+      cout << " *** substa:  Bad s0 = " << s0 << endl;
+      return;
+   }
+//
+   s1 = s1/s0;
+//   
+   s2 = s2/s0 - s1*s1;
+   if (s2 < 0) {
+      cout << " *** substa:  Warning s2 = " << s2 << endl;
+      if (TMath::Abs(s2/s1) < 1.e-07)
+         s2 = 0;
+      else {
+         cout << " *** substa:  Bad s2 = " << s2 << endl;
+	 exit (-2);
+      }	 
+   }
+//
+   ave = s1;
+   rms = TMath::Sqrt(s2);
+   err = rms / TMath::Sqrt(s0);
+//
+//
+   return;
+}
+// =====================================================================
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Generation/ele_e0100_sc1/jobOption.ParticleGun_ele_e0100_sc1.py b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Generation/ele_e0100_sc1/jobOption.ParticleGun_ele_e0100_sc1.py
new file mode 100644
index 000000000000..180a9de85d46
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Generation/ele_e0100_sc1/jobOption.ParticleGun_ele_e0100_sc1.py
@@ -0,0 +1,57 @@
+#! -*- python -*-
+import math
+import ROOT
+
+include("GeneratorUtils/StdEvgenSetup.py")
+theApp.EvtMax = 5000
+
+import ParticleGun as PG
+      
+class MyParticleSampler(PG.ParticleSampler):
+
+    "A special sampler to shoot in the HEC wheel 1 or 2"
+	
+    def __init__(self, pid = 11,
+                 energy = 100000.,
+                 phi = PG.UniformSampler(0, 2*math.pi),
+                 eta = PG.UniformSampler(1.5,3.3),
+                 pv  = PG.PosSampler(0, 0, 4319.5)):
+        self.pid = pid
+        self.energy = energy
+        self.eta = eta
+        self.phi = phi
+        self.pv =  pv
+          
+    def shoot(self):
+        mass = 0.510998928
+        momentum = math.sqrt(self.energy**2 - mass**2)
+        v4eta = self.eta.shoot()
+        v4 = ROOT.TLorentzVector()
+        v4.SetPtEtaPhiE(momentum/math.cosh(v4eta), v4eta, self.phi.shoot(), self.energy)
+        vv = self.pv.shoot()
+        v4theta = v4.Theta()
+        v4phi = v4.Phi()
+        z = vv.Z()
+        ## Calculate r at z
+        r = math.tan(v4theta)*z
+        x = r*math.cos(v4phi)
+        y = r*math.sin(v4phi)
+        vv.SetX(x)
+        vv.SetY(y)
+        pdg_id = self.pid.shoot()
+
+        print "pdg_id = ",pdg_id
+        rtn = [PG.SampledParticle(pdg_id,v4,vv)]
+
+        return rtn
+
+
+topSeq += PG.ParticleGun()
+topSeq.ParticleGun.sampler = MyParticleSampler()
+
+include("GeneratorUtils/postJO.WriteHepMC.py")
+include("GeneratorUtils/postJO.Output.py")
+
+print topSeq
+
+
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Generation/ele_e0100_sc2/jobOption.ParticleGun_ele_e0100_sc2.py b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Generation/ele_e0100_sc2/jobOption.ParticleGun_ele_e0100_sc2.py
new file mode 100644
index 000000000000..b95155d6fad4
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Generation/ele_e0100_sc2/jobOption.ParticleGun_ele_e0100_sc2.py
@@ -0,0 +1,57 @@
+#! -*- python -*-
+import math
+import ROOT
+
+include("GeneratorUtils/StdEvgenSetup.py")
+theApp.EvtMax = 5000
+
+import ParticleGun as PG
+      
+class MyParticleSampler(PG.ParticleSampler):
+
+    "A special sampler to shoot in the HEC wheel 1 or 2"
+	
+    def __init__(self, pid = 11,
+                 energy = 100000.,
+                 phi = PG.UniformSampler(0, 2*math.pi),
+                 eta = PG.UniformSampler(1.6,3.3),
+                 pv  = PG.PosSampler(0, 0, 5175.0)):
+        self.pid = pid
+        self.energy = energy
+        self.eta = eta
+        self.phi = phi
+        self.pv =  pv
+          
+    def shoot(self):
+        mass = 0.510998928
+        momentum = math.sqrt(self.energy**2 - mass**2)
+        v4eta = self.eta.shoot()
+        v4 = ROOT.TLorentzVector()
+        v4.SetPtEtaPhiE(momentum/math.cosh(v4eta), v4eta, self.phi.shoot(), self.energy)
+        vv = self.pv.shoot()
+        v4theta = v4.Theta()
+        v4phi = v4.Phi()
+        z = vv.Z()
+        ## Calculate r at z
+        r = math.tan(v4theta)*z
+        x = r*math.cos(v4phi)
+        y = r*math.sin(v4phi)
+        vv.SetX(x)
+        vv.SetY(y)
+        pdg_id = self.pid.shoot()
+
+        print "pdg_id = ",pdg_id
+        rtn = [PG.SampledParticle(pdg_id,v4,vv)]
+
+        return rtn
+
+
+topSeq += PG.ParticleGun()
+topSeq.ParticleGun.sampler = MyParticleSampler()
+
+include("GeneratorUtils/postJO.WriteHepMC.py")
+include("GeneratorUtils/postJO.Output.py")
+
+print topSeq
+
+
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Generation/env.sh b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Generation/env.sh
new file mode 100644
index 000000000000..a8d5b9ceb851
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Generation/env.sh
@@ -0,0 +1,4 @@
+## Example, how the release was set up
+#export AtlasSetup=/afs/cern.ch/atlas/software/releases/17.6.0/AtlasSetup
+#source $AtlasSetup/scripts/asetup.sh --tags=AtlasProduction,slc5,17.6.0.3,setup,opt --testarea $PWD
+##
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Generation/gam_e0100_sc1/jobOption.ParticleGun_gam_e0100_sc1.py b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Generation/gam_e0100_sc1/jobOption.ParticleGun_gam_e0100_sc1.py
new file mode 100644
index 000000000000..754157ab6a62
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Generation/gam_e0100_sc1/jobOption.ParticleGun_gam_e0100_sc1.py
@@ -0,0 +1,57 @@
+#! -*- python -*-
+import math
+import ROOT
+
+include("GeneratorUtils/StdEvgenSetup.py")
+theApp.EvtMax = 5000
+
+import ParticleGun as PG
+      
+class MyParticleSampler(PG.ParticleSampler):
+
+    "A special sampler to shoot in the HEC wheel 1 or 2"
+	
+    def __init__(self, pid = 22,
+                 energy = 100000.,
+                 phi = PG.UniformSampler(0, 2*math.pi),
+                 eta = PG.UniformSampler(1.5,3.3),
+                 pv  = PG.PosSampler(0, 0, 4319.5)):
+        self.pid = pid
+        self.energy = energy
+        self.eta = eta
+        self.phi = phi
+        self.pv =  pv
+          
+    def shoot(self):
+        mass = 0.0
+        momentum = math.sqrt(self.energy**2 - mass**2)
+        v4eta = self.eta.shoot()
+        v4 = ROOT.TLorentzVector()
+        v4.SetPtEtaPhiE(momentum/math.cosh(v4eta), v4eta, self.phi.shoot(), self.energy)
+        vv = self.pv.shoot()
+        v4theta = v4.Theta()
+        v4phi = v4.Phi()
+        z = vv.Z()
+        ## Calculate r at z
+        r = math.tan(v4theta)*z
+        x = r*math.cos(v4phi)
+        y = r*math.sin(v4phi)
+        vv.SetX(x)
+        vv.SetY(y)
+        pdg_id = self.pid.shoot()
+
+        print "pdg_id = ",pdg_id
+        rtn = [PG.SampledParticle(pdg_id,v4,vv)]
+
+        return rtn
+
+
+topSeq += PG.ParticleGun()
+topSeq.ParticleGun.sampler = MyParticleSampler()
+
+include("GeneratorUtils/postJO.WriteHepMC.py")
+include("GeneratorUtils/postJO.Output.py")
+
+print topSeq
+
+
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Generation/gam_e0100_sc2/jobOption.ParticleGun_gam_e0100_sc2.py b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Generation/gam_e0100_sc2/jobOption.ParticleGun_gam_e0100_sc2.py
new file mode 100644
index 000000000000..f298829dfef0
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Generation/gam_e0100_sc2/jobOption.ParticleGun_gam_e0100_sc2.py
@@ -0,0 +1,57 @@
+#! -*- python -*-
+import math
+import ROOT
+
+include("GeneratorUtils/StdEvgenSetup.py")
+theApp.EvtMax = 5000
+
+import ParticleGun as PG
+      
+class MyParticleSampler(PG.ParticleSampler):
+
+    "A special sampler to shoot in the HEC wheel 1 or 2"
+	
+    def __init__(self, pid = 22,
+                 energy = 100000.,
+                 phi = PG.UniformSampler(0, 2*math.pi),
+                 eta = PG.UniformSampler(1.6,3.3),
+                 pv  = PG.PosSampler(0, 0, 5175.0)):
+        self.pid = pid
+        self.energy = energy
+        self.eta = eta
+        self.phi = phi
+        self.pv =  pv
+          
+    def shoot(self):
+        mass = 0.0
+        momentum = math.sqrt(self.energy**2 - mass**2)
+        v4eta = self.eta.shoot()
+        v4 = ROOT.TLorentzVector()
+        v4.SetPtEtaPhiE(momentum/math.cosh(v4eta), v4eta, self.phi.shoot(), self.energy)
+        vv = self.pv.shoot()
+        v4theta = v4.Theta()
+        v4phi = v4.Phi()
+        z = vv.Z()
+        ## Calculate r at z
+        r = math.tan(v4theta)*z
+        x = r*math.cos(v4phi)
+        y = r*math.sin(v4phi)
+        vv.SetX(x)
+        vv.SetY(y)
+        pdg_id = self.pid.shoot()
+
+        print "pdg_id = ",pdg_id
+        rtn = [PG.SampledParticle(pdg_id,v4,vv)]
+
+        return rtn
+
+
+topSeq += PG.ParticleGun()
+topSeq.ParticleGun.sampler = MyParticleSampler()
+
+include("GeneratorUtils/postJO.WriteHepMC.py")
+include("GeneratorUtils/postJO.Output.py")
+
+print topSeq
+
+
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Generation/pos_e0100_sc1/jobOption.ParticleGun_pos_e0100_sc1.py b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Generation/pos_e0100_sc1/jobOption.ParticleGun_pos_e0100_sc1.py
new file mode 100644
index 000000000000..adbe7f150d38
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Generation/pos_e0100_sc1/jobOption.ParticleGun_pos_e0100_sc1.py
@@ -0,0 +1,57 @@
+#! -*- python -*-
+import math
+import ROOT
+
+include("GeneratorUtils/StdEvgenSetup.py")
+theApp.EvtMax = 5000
+
+import ParticleGun as PG
+      
+class MyParticleSampler(PG.ParticleSampler):
+
+    "A special sampler to shoot in the HEC wheel 1 or 2"
+	
+    def __init__(self, pid = -11,
+                 energy = 100000.,
+                 phi = PG.UniformSampler(0, 2*math.pi),
+                 eta = PG.UniformSampler(1.5,3.3),
+                 pv  = PG.PosSampler(0, 0, 4319.5)):
+        self.pid = pid
+        self.energy = energy
+        self.eta = eta
+        self.phi = phi
+        self.pv =  pv
+          
+    def shoot(self):
+        mass = 0.510998928
+        momentum = math.sqrt(self.energy**2 - mass**2)
+        v4eta = self.eta.shoot()
+        v4 = ROOT.TLorentzVector()
+        v4.SetPtEtaPhiE(momentum/math.cosh(v4eta), v4eta, self.phi.shoot(), self.energy)
+        vv = self.pv.shoot()
+        v4theta = v4.Theta()
+        v4phi = v4.Phi()
+        z = vv.Z()
+        ## Calculate r at z
+        r = math.tan(v4theta)*z
+        x = r*math.cos(v4phi)
+        y = r*math.sin(v4phi)
+        vv.SetX(x)
+        vv.SetY(y)
+        pdg_id = self.pid.shoot()
+
+        print "pdg_id = ",pdg_id
+        rtn = [PG.SampledParticle(pdg_id,v4,vv)]
+
+        return rtn
+
+
+topSeq += PG.ParticleGun()
+topSeq.ParticleGun.sampler = MyParticleSampler()
+
+include("GeneratorUtils/postJO.WriteHepMC.py")
+include("GeneratorUtils/postJO.Output.py")
+
+print topSeq
+
+
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Generation/pos_e0100_sc2/jobOption.ParticleGun_pos_e0100_sc2.py b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Generation/pos_e0100_sc2/jobOption.ParticleGun_pos_e0100_sc2.py
new file mode 100644
index 000000000000..85e824c4d465
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Generation/pos_e0100_sc2/jobOption.ParticleGun_pos_e0100_sc2.py
@@ -0,0 +1,57 @@
+#! -*- python -*-
+import math
+import ROOT
+
+include("GeneratorUtils/StdEvgenSetup.py")
+theApp.EvtMax = 5000
+
+import ParticleGun as PG
+      
+class MyParticleSampler(PG.ParticleSampler):
+
+    "A special sampler to shoot in the HEC wheel 1 or 2"
+	
+    def __init__(self, pid = -11,
+                 energy = 100000.,
+                 phi = PG.UniformSampler(0, 2*math.pi),
+                 eta = PG.UniformSampler(1.6,3.3),
+                 pv  = PG.PosSampler(0, 0, 5175.0)):
+        self.pid = pid
+        self.energy = energy
+        self.eta = eta
+        self.phi = phi
+        self.pv =  pv
+          
+    def shoot(self):
+        mass = 0.510998928
+        momentum = math.sqrt(self.energy**2 - mass**2)
+        v4eta = self.eta.shoot()
+        v4 = ROOT.TLorentzVector()
+        v4.SetPtEtaPhiE(momentum/math.cosh(v4eta), v4eta, self.phi.shoot(), self.energy)
+        vv = self.pv.shoot()
+        v4theta = v4.Theta()
+        v4phi = v4.Phi()
+        z = vv.Z()
+        ## Calculate r at z
+        r = math.tan(v4theta)*z
+        x = r*math.cos(v4phi)
+        y = r*math.sin(v4phi)
+        vv.SetX(x)
+        vv.SetY(y)
+        pdg_id = self.pid.shoot()
+
+        print "pdg_id = ",pdg_id
+        rtn = [PG.SampledParticle(pdg_id,v4,vv)]
+
+        return rtn
+
+
+topSeq += PG.ParticleGun()
+topSeq.ParticleGun.sampler = MyParticleSampler()
+
+include("GeneratorUtils/postJO.WriteHepMC.py")
+include("GeneratorUtils/postJO.Output.py")
+
+print topSeq
+
+
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/ReadMe.txt b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/ReadMe.txt
new file mode 100644
index 000000000000..def9ca47f7d9
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/ReadMe.txt
@@ -0,0 +1,188 @@
+
+      HOW TO GET SAMPLING FRACTIONS FOR HEC,
+      USING GEANT4-BASED SIMULATIONS
+  
+      Andrei.Kiryunin@cern.ch
+  
+      March 19-th, 2020
+  
+  
+
+Procedure to obtain sampling fractions of two HEC wheels (front wheel with 25 mm copper plates and rear wheel with 50 mm copper plates) consists of a few steps:
+
+  - Generation of sets of primary particles to scan over HEC wheels
+
+  - Simulation of detector response (keeping calibration hits in calorimeter cells) with Geant4 version under study and selected physics lists
+
+  - Dumping calibration hits and primary particle information to ROOT-tree
+
+  - Analysis of relevant variables in ROOT-tree to calculate sampling fractions
+
+All necessary files to run generation, simulation and analysis can be found in corresponding sub-directories under the main directory:
+
+  MainDir=$PWD/HEC_SamplingFractions
+   
+
+      GENERATION OF PRIMARY PARTICLES
+
+  - Single-particle events
+   
+  - Three types of primary particles: electrons, positrons, and gammas
+
+    + flat-distributed in pseudorapidity (eta) and azimuthal angle
+    + fixed Z-position of vertices (Zvtx) at front face of wheels
+    + back pointed to interaction point
+    + particle energy = 100 GeV
+
+  - Two scans over HEC wheels
+
+    1. Front wheel scan: 1.5 < eta < 3.3, Zvtx = 431.95 cm 
+    2. Rear wheel scan: 1.6 < eta < 3.3, Zvtx = 571.5 cm
+
+  - In total, six generation jobs with names:
+
+    JOB = PPP_eEEEE_scN, where
+
+    + PPP - particle name ("ele" - electrons, "pos" - positrons, "gam" - gammas)
+    + EEEE - particle energy (always 0100)
+    + N - scan number (1 - front wheel, 2 - rear wheel)
+
+  - 5,000 generated events per job
+
+All relevant files can be found in directory '$MainDir/Generation'. Generation was done with Athena release 17.6.0. 
+
+Prepare environment:
+
+  export AtlasSetup=/afs/cern.ch/atlas/software/releases/17.6.0/AtlasSetup
+  source $AtlasSetup/scripts/asetup.sh --tags=AtlasProduction,slc5,17.6.0.3,setup,opt --testarea $PWD
+
+and check out 'ParticleGun' package:
+
+  cmt co -r ParticleGun-01-01-00 Generators/ParticleGun
+  cd Generators/ParticleGun/cmt
+  gmake
+
+Six sub-directories '$MainDir/Generation/$JOB' contain files 'jobOption.ParticleGun_$JOB.py'.
+
+To run generation, issue in corresponding sub-directory:
+
+  . $MainDir/Generation/env.sh
+   
+and then run athena:   
+
+  athena.py  jobOption.ParticleGun_$JOB.py  
+
+Output files: 'evgen.pool.root' - can be found (temporarily) under:
+
+  /eos/user/k/kiryunin/HEC_SF/Generation/$JOB
+   
+   
+      SIMULATION AND CALIBRATION-HITS DUMPING
+
+All relevant files can be found in directory '$MainDir/Simulation'. Simulations were done with Athena release 20.3.7, using Geant4 version 10.1 with patch03 (internal ATLAS build geant4.10.1.patch03.atlas01). 
+
+Prepare environment:
+
+  export AtlasSetup=/afs/cern.ch/atlas/software/releases/20.3.7/AtlasSetup
+  source $AtlasSetup/scripts/asetup.sh --tags=AtlasProduction,20.3.7.1,setup,opt --testarea $PWD
+
+and check out necessary packages:
+
+  cmt co -r G4PhysicsLists-01-00-05 Simulation/G4Utilities/G4PhysicsLists
+  cd Simulation/G4Utilities/G4PhysicsLists/cmt
+  gmake
+  cd $MainDir/Simulation
+
+  cmt co -r G4AtlasApps-00-11-00-18 Simulation/G4Atlas/G4AtlasApps
+  cd Simulation/G4Atlas/G4AtlasApps/cmt
+  gmake
+  cd $MainDir/Simulation
+
+  cmt co -r CaloD3PDMaker-00-03-09-01 PhysicsAnalysis/D3PDMaker/CaloD3PDMaker
+  cd PhysicsAnalysis/D3PDMaker/CaloD3PDMaker/cmt
+  gmake
+  cd $MainDir/Simulation
+
+Each job is splitted to 25 sub-jobs to simulate 200 events. To prepare all necessary working directories and files, run in the '$MainDir/Simulation':
+
+  ./mkjob.sh PhLis
+  
+where PhLis - name of the physics list. Two physics lists: FTFP_BERT and FTFP_BERT_ATL - are foreseen for current work. As the result, 150 working directories: 
+
+  $MainDir/Simulation/10.1_$PhLis/$JOB/$SubJob
+
+are created. '$SubJob' runs from '01' to '25'. Each working directory contains:
+
+  - link to necessary file 'evgen.pool.root' with generated events
+
+  - 'simJob.sh' - script to run simulations, where
+    + physics list
+    + random seed number
+    + number of generated events to be skipped
+    are defined
+  - 'preInclude.VertexPosOff.py' - job config to turn off vertex positioners
+
+  - 'dumpJob.sh' - script to dump calibration hits and primary particle information from simulated files 'myHITS.pool.root' (where name of output ROOT-file is defined)
+  - 'addMCpart.py' - job config to add primary particle information  
+
+To run simulations, issue in a working directory:
+
+  . $MainDir/Simulation/env.sh
+   
+and then run:
+
+  ./simJob.sh
+   
+[Current simulations were done in the LSF batch system.]      
+
+To dump hits (after successful simulations), run in the same working directory:
+
+  ./dumpJob.sh
+   
+Output ROOT-tree files: $SubJob.root - can be found (temporarily) under:
+
+  /eos/user/k/kiryunin/HEC_SF/10.1_$PhLis/DumpHits/$JOB
+
+
+      ANALYSIS OF CALIBRATION HITS TO GET SAMPLING FRACTIONS (SF)
+
+Sampling fraction of a HEC wheel is determined as 
+
+  SF = < E_ACT > / < E_ACT + E_PAS >, where
+
+  - E_ACT - sum of visible energies in active parts of a HEC wheel
+  - E_PAS - sum of visible energies in passive parts of a HEC wheel
+
+To calculate final Sampling Fraction (SF0) for each HEC wheel
+
+  - loop over 75 ROOT-tree files from the appropriate scan (three types of primary particles, 25 sub-jobs -> 15,000 events)
+  - get pseudorapidity dependence of SF
+  - fit this dependence to the constant, excluding edges and regions with tie-rods from the fit
+  - take parameter of this fit as the final Sampling Fraction SF0
+   
+ROOT-based analysis is done in the directory '$MainDir/Analysis', where all necessary C-routines can be found. The workflow of the analysis is the following:
+
+  . $MainDir/Analysis/env.sh  #  use current ATLAS setup
+
+  root -l  #  check calibration hits and get ASCII-files with dependence SF(eta)
+  .L store_eta.C
+  store_eta ("10.1", "FTFP_BERT_ATL", "/eos/user/k/kiryunin/HEC_SF")
+  store_eta ("10.1", "FTFP_BERT",     "/eos/user/k/kiryunin/HEC_SF")
+  .q
+
+  root -l -b  #  get final Sampling Fractions for both HEC wheels
+  .L get_SF.C 
+  get_SF ("10.1", "FTFP_BERT_ATL")
+  get_SF ("10.1", "FTFP_BERT")
+  .q
+
+Output of the analysis is in the directories '$MainDir/Analysis/10.1_$PhLis'. They contain following files:
+
+  - 'samf_eta_fwh-EM_e0100_sc1.dat' and 'samf_eta_rwh-EM_e0100_sc2.dat' with pseudorapidity dependence of SF for HEC front and rear wheels
+    
+  - 'SF_EM_e0100.ps' with plots (for both wheels) of dependence SF(eta) with results of the fit:
+    + fit to constant and linear fit
+    + whole eta-range, cut eta-range, and cut eta-range without tie-rod regions
+    
+  - 'HEC_SF.txt": fit parameter with error and its rounding value as final sampling fraction (for both wheels)
+
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Simulation/addMCpart.py b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Simulation/addMCpart.py
new file mode 100644
index 000000000000..f0b6da492e93
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Simulation/addMCpart.py
@@ -0,0 +1,13 @@
+from D3PDMakerConfig.D3PDMakerFlags import D3PDMakerFlags
+D3PDMakerFlags.TruthSGKey.set_Value_and_Lock('TruthEvent') 
+from TruthD3PDAnalysis.GenObjectsFilterTool import *
+from TruthD3PDMaker.GenVertexD3PDObject import GenVertexD3PDObject
+alg += GenVertexD3PDObject( 1, filter = AllTrackFilterTool() )
+pfilter=TruthTrackFilterTool()
+pfilter.PtMin=-1.
+pfilter.EtaMax=-1.
+pfilter.SelectTruthTracks=False
+from TruthD3PDMaker.GenParticleD3PDObject import GenParticleD3PDObject
+alg += GenParticleD3PDObject( 1, filter = pfilter )
+#from TruthD3PDMaker.GenEventD3PDObject import GenEventD3PDObject
+#alg += GenEventD3PDObject( 0, filter = AllTrackFilterTool() )
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Simulation/dumphit.sh b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Simulation/dumphit.sh
new file mode 100644
index 000000000000..0d1a69da100b
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Simulation/dumphit.sh
@@ -0,0 +1,3 @@
+#
+athena.py -c 'inputFiles=["myHITS.pool.root"];outputNtupleFile="calohitD3PD_from_esd.root"' CaloD3PDMaker/MakeCaloHitD3PD_topOptions.py addMCpart.py
+#
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Simulation/env.sh b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Simulation/env.sh
new file mode 100755
index 000000000000..be339296d171
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Simulation/env.sh
@@ -0,0 +1,2 @@
+export AtlasSetup=/afs/cern.ch/atlas/software/releases/20.3.7/AtlasSetup
+source $AtlasSetup/scripts/asetup.sh --tags=AtlasProduction,20.3.7.1,setup,opt --testarea $PWD
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Simulation/mkjob.sh b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Simulation/mkjob.sh
new file mode 100755
index 000000000000..3432c5df90c9
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Simulation/mkjob.sh
@@ -0,0 +1,115 @@
+#!/bin/sh
+#
+#  Create working directories with necessary files to run simulations:
+#  - different primary particles
+#  - different Z-positions of primary vertices (i.e. scans over HEC)
+#  - 25 sub-jobs with 200 events
+#
+if [ $# != 1 ] ; then
+   echo " *** Argument is missing:  Physics list  !!! *** "
+   exit
+fi
+#
+export PhLis=$1
+#
+export DIRM=$PWD
+#
+export SimDIR=$DIRM/10.1-$PhLis
+mkdir -p $SimDIR
+#
+export GenDIR=/eos/user/k/kiryunin/HEC_SF/Generation
+#
+export ENER=e0100
+#
+#
+cat >dum0.sh <<\end0 
+#
+#  Primary vertices
+#
+   ./dum1.sh $1 sc1  # scan over front HEC wheel
+   ./dum1.sh $1 sc2  # scan over rear  HEC wheel
+# 
+end0
+#
+cat >dum1.sh <<\end1  
+#
+#  Loop over sub-jobs
+#
+   for (( is = 1 ; is <= 25 ; is ++ )) ; do
+      ./dum2.sh $1"_"$ENER"_"$2 $is
+   done   
+#
+end1
+#
+cat >dum2.sh <<\end2  
+#
+   echo " === Job: $1 $2 === "
+#
+#  Create working directory to run simulations and dump calibration hits
+#
+   if [ $2 -lt 10 ] ; then   
+      subjb=0$2
+   else
+      subjb=$2
+   fi
+   njob=$(($2-1))
+#   
+   dir=$SimDIR/$1/$subjb
+#
+   if test -s $dir ; then
+      echo " *** Directory $dir does exist !!! *** "
+      exit
+   fi        
+#
+   mkdir -p $dir
+   cd $dir
+#
+   ln -s $GenDIR/$1/evgen.pool.root .   
+#
+   cp -i $DIRM/preInclude.VertexPosOff.py .
+   cp -i $DIRM/addMCpart.py .
+#
+#  Prepare file to run simulations
+#
+   fsim=$DIRM/simul.sh
+#
+   seed0=5432100
+   seed1=$((${seed0}+$njob))
+   sed -e 's/default:'$seed0'/default:'$seed1'/' $fsim >dumm1
+#
+   skip0=0
+   skip1=$((200*$njob))
+   sed -e 's/skipEvents '$skip0'/skipEvents '$skip1'/' dumm1 >dumm2
+#
+   sed -e 's/PhysicsList/'$PhLis'/' dumm2 >simJob.sh
+   rm -f dumm?
+#   
+   chmod 755 simJob.sh
+#
+#  Prepare file to dump calibration hits
+#
+   fdump=$DIRM/dumphit.sh
+#
+   sed -e 's/calohitD3PD_from_esd/'$subjb'/' $fdump >dumpJob.sh
+#
+   chmod 755 dumpJob.sh
+#
+#
+   cd $DIRM
+   echo " === made === "
+#   
+end2
+#
+chmod 755 dum?.sh
+#
+#  Primary particles 
+#
+./dum0.sh ele  # electrons
+./dum0.sh pos  # positrons
+./dum0.sh gam  # gammas
+#
+#
+rm -f dum?.sh
+#
+exit
+########################################################################
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Simulation/preInclude.VertexPosOff.py b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Simulation/preInclude.VertexPosOff.py
new file mode 100644
index 000000000000..bd279e04e975
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Simulation/preInclude.VertexPosOff.py
@@ -0,0 +1,15 @@
+################################################
+##
+## Job config to turn off vertex positioners
+##
+## Author: Andrea di Simone
+##
+################################################
+
+from G4AtlasApps.SimFlags import simFlags
+if simFlags.EventFilter.statusOn:
+    simFlags.VertexFromCondDB.set_Off()
+    simFlags.EventFilter.switchFilterOff('EtaPhiFilters')
+    simFlags.EventFilter.switchFilterOff('VertexPositioner')
+    simFlags.EventFilter.switchFilterOff('VertexRangeChecker')
+
diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Simulation/simul.sh b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Simulation/simul.sh
new file mode 100644
index 000000000000..088e922313cd
--- /dev/null
+++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/HEC_SamplingFractions/Simulation/simul.sh
@@ -0,0 +1,3 @@
+#
+Sim_tf.py --inputEVNTFile evgen.pool.root --outputHITSFile myHITS.pool.root --conditionsTag "default:OFLCOND-RUN12-SDR-31" --physicsList PhysicsList --randomSeed "default:5432100" --maxEvents 200 --skipEvents 0 --preInclude "EVNTtoHITS:preInclude.VertexPosOff.py,SimulationJobOptions/preInclude.CalHits.py,SimulationJobOptions/preInclude.ParticleID.py" --geometryVersion "default:ATLAS-R2-2016-00-00-00_VALIDATION"
+#
-- 
GitLab