32 #include "GaudiKernel/ISvcLocator.h"
33 #include "GaudiKernel/StatusCode.h"
43 ISvcLocator* pSvcLocator)
45 m_outputFile(nullptr),
46 m_clusterCollName(
"CaloTopoClusters"),
47 m_NormalizationType(
"Lin"),
48 m_NormalizationTypeNumber(0),
49 m_ClassificationType(
"None"),
50 m_ClassificationTypeNumber(0)
53 std::vector<Gaudi::Histo1DDef> dims(6);
54 dims[0] = Gaudi::Histo1DDef(
"side",-1.5,1.5,1);
55 dims[1] = Gaudi::Histo1DDef(
"|eta|",0.,5.,25);
56 dims[2] = Gaudi::Histo1DDef(
"phi",-
M_PI,
M_PI,1);
57 dims[3] = Gaudi::Histo1DDef(
"log10(E_clus (MeV))",log10(200),log10(1e6),14);
58 dims[4] = Gaudi::Histo1DDef(
"log10(<rho_cell (MeV/mm^3)>)-log10(E_clus (MeV))",-9.0,-4.0,20);
59 dims[5] = Gaudi::Histo1DDef(
"log10(lambda_clus (mm))",0.0,4.0,20);
97 ATH_MSG_INFO(
"Using weighting proportional to E_calib" );
101 ATH_MSG_INFO(
"Using weighting proportional to log(E_calib)" );
105 ATH_MSG_INFO(
"Using weighting proportional to 1/N_Clus_E_calib>0" );
118 ATH_MSG_INFO(
"Expecting ParticleID simulation as input -- use EM type clusters only" );
122 ATH_MSG_INFO(
"Expecting ParticleID simulation as input -- use HAD type clusters only" );
138 for(
unsigned int idim=0;idim<
m_dimensions.size();idim++) {
155 else if (
m_dimensions[idim].
title() ==
"log10(<rho_cell (MeV/mm^3)>)-log10(E_clus (MeV))" )
160 if ( ilogE < 0 || ilogrho < 0 || iloglambda < 0 ) {
161 ATH_MSG_FATAL(
" Mandatory dimension log10E, log10rho or log10lambda missing ..." );
162 return StatusCode::FAILURE;
173 m_hclus.resize(nside*neta*nphi*nlogE,
nullptr);
174 for (
int jside=0;jside<nside;jside++) {
175 for (
int jeta=0;jeta<neta;jeta++) {
176 for (
int jphi=0;jphi<nphi;jphi++) {
177 for(
int jlogE=0;jlogE<nlogE;jlogE++) {
178 TString
hname(
"hclus");
215 int iH = jlogE*nphi*neta*nside+jphi*neta*nside+jeta*nside+jside;
224 m_hclus[iH]->SetXTitle(
"log10(<#rho_{cell}> (MeV/mm^{3})) - log10(E_{clus} (MeV))");
225 m_hclus[iH]->SetYTitle(
"log10(#lambda_{clus} (mm))");
226 m_hclus[iH]->SetZTitle(
"Number of Clusters");
234 return StatusCode::SUCCESS;
248 return StatusCode::SUCCESS;
258 double eCalibTot(0.);
259 double nClusECalibGt0 = 0.0;
264 ATH_MSG_ERROR(
"Failed to retrieve cluster moment ENG_CALIB_TOT" );
265 return StatusCode::FAILURE;
270 ATH_MSG_ERROR(
"Failed to retrieve cluster moment ENG_CALIB_FRAC_EM" );
271 return StatusCode::FAILURE;
284 if ( eCalibTot > 0 ) {
285 const double inv_eCalibTot = 1. / eCalibTot;
286 const double inv_nClusECalibGt0 = 1. / nClusECalibGt0;
288 double eng = pClus->e();
293 ATH_MSG_ERROR(
"Failed to retrieve cluster moment ENG_CALIB_FRAC_EM" );
294 return StatusCode::FAILURE;
302 double eta = fabs(pClus->eta());
316 iside = (
int)(nside*(((pClus->eta()<0?-1.0:1.0) - hd.lowEdge())
317 /(hd.highEdge()-hd.lowEdge())));
318 if ( iside < 0 || iside > nside-1 ) {
320 iside <<
" not in [0," << nside-1 <<
"]" );
328 ieta = (
int)(neta*((
eta - hd.lowEdge())
329 /(hd.highEdge()-hd.lowEdge())));
330 if ( ieta < 0 || ieta > neta-1 ) {
332 ieta <<
" not in [0," << neta-1 <<
"]" );
339 iphi = (
int)(nphi*((pClus->phi() - hd.lowEdge())
340 /(hd.highEdge()-hd.lowEdge())));
341 if ( iphi < 0 || iphi > nphi-1 ) {
343 iphi <<
" not in [0," << nphi-1 <<
"]" );
350 ilogE = (
int)(nlogE*((log10(eng) - hd.lowEdge())
351 /(hd.highEdge()-hd.lowEdge())));
352 if ( ilogE >= 0 && ilogE < nlogE ) {
353 double dens=0,lamb=0,ecal=0;
355 ATH_MSG_ERROR(
"Failed to retrieve cluster moment ENG_CALIB_TOT" );
356 return StatusCode::FAILURE;
359 ATH_MSG_ERROR(
"Failed to retrieve cluster moment FIRST_ENG_DENS" );
360 return StatusCode::FAILURE;
363 ATH_MSG_ERROR(
"Failed to retrieve cluster moment CENTER_LAMBDA" );
364 return StatusCode::FAILURE;
370 int iH = ilogE*nphi*neta*nside+iphi*neta*nside+ieta*nside+iside;
374 norm = ecal*inv_eCalibTot;
378 norm = log10(ecal*inv_eCalibTot)+2.0;
381 norm = inv_nClusECalibGt0;
387 m_hclus[iH]->Fill(log10(dens)-log10(eng),log10(lamb),
norm);
396 return StatusCode::SUCCESS;
400 for (
unsigned int i=0;
i<dims.size();
i++) {
407 for (std::pair<const std::string, Gaudi::Histo1DDef>&
p :
m_dimensionsmap) {
410 <<
p.second.title() <<
", [" <<
p.second.lowEdge()
411 <<
", " <<
p.second.highEdge()
412 <<
", " <<
p.second.bins()