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 && nClusECalibGt0 > 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()