diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/CMakeLists.txt b/Simulation/Tools/CaloSamplingFractionAnalysis/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..b9a3aef0a54d2c4dae34fe766cc838847ed93043 --- /dev/null +++ b/Simulation/Tools/CaloSamplingFractionAnalysis/CMakeLists.txt @@ -0,0 +1,43 @@ +################################################################################ +# Package: CaloSamplingFractionAnalysis +################################################################################ + +# Declare the package name: +atlas_subdir( CaloSamplingFractionAnalysis ) + +# Declare the package's dependencies: +atlas_depends_on_subdirs( PUBLIC + GaudiKernel + PRIVATE + Calorimeter/CaloDetDescr + Calorimeter/CaloIdentifier + Calorimeter/CaloSimEvent + Calorimeter/CaloEvent + Control/AthenaBaseComps + DetectorDescription/GeoModel/GeoAdaptors + Event/EventInfo + LArCalorimeter/LArSimEvent + LArCalorimeter/LArElecCalib + TileCalorimeter/TileDetDescr + TileCalorimeter/TileSimEvent ) + +# External dependencies: +find_package( CLHEP ) +find_package( HepMC ) +find_package( ROOT COMPONENTS Core Tree MathCore Hist RIO pthread Table MathMore Minuit Minuit2 Matrix Physics HistPainter Rint Graf Graf3d Gpad Html Postscript Gui GX11TTF GX11 ) + +# tag ROOTBasicLibs was not recognized in automatic conversion in cmt2cmake + +# tag ROOTSTLDictLibs was not recognized in automatic conversion in cmt2cmake + +# Component(s) in the package: +atlas_add_component( CaloSamplingFractionAnalysis + src/*.cxx + src/components/*.cxx + INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS} ${HEPMC_INCLUDE_DIRS} + LINK_LIBRARIES ${ROOT_LIBRARIES} ${CLHEP_LIBRARIES} ${HEPMC_LIBRARIES} GaudiKernel CaloDetDescrLib CaloIdentifier CaloSimEvent CaloEvent AthenaBaseComps GeoAdaptors EventInfo LArSimEvent TileDetDescr TileSimEvent ) + +# Install files from the package: +# atlas_install_headers( CaloSamplingFractionAnalysis ) +# atlas_install_joboptions( share/*.py ) + diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/LarEMSamplingFraction_analysis.C b/Simulation/Tools/CaloSamplingFractionAnalysis/share/LarEMSamplingFraction_analysis.C new file mode 100644 index 0000000000000000000000000000000000000000..6b00226def4cbe9e4fae80eb57a9d10237886614 --- /dev/null +++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/LarEMSamplingFraction_analysis.C @@ -0,0 +1,39 @@ +{ + TFile* file = TFile::Open("LArEM_SF.root"); + new TCanvas("Etot","Etot"); + mytree->Draw("Sum$(energy_active_total+energy_inactive_total)"); + + new TCanvas("ELAr_hit","ELAr_hit"); + mytree->Draw("(energy_hit[0]+energy_hit[1]+energy_hit[2]+energy_hit[3])"); + + //new TCanvas("SF","SF"); + //mytree->Draw("Sum$(energy_active_total)/Sum$(energy_active_total+energy_inactive_total):mc_eta>>SF(14,0,1.4,20,0.15,0.25)","","colz"); + + //new TCanvas("LAr_barrel","LAr_barrel"); + //mytree->Draw("(energy_active_total[0]+energy_active_total[1]+energy_active_total[2]+energy_active_total[3]):mc_eta>>LAr_barrel(56,0,1.4,50,7000,12000)","","colz"); + + new TCanvas("SF_LAr_barrel_calibhit","SF_LAr_barrel_calibhit"); + mytree->Draw("(energy_active_total[0]+energy_active_total[1]+energy_active_total[2]+energy_active_total[3])/(energy_active_total[0]+energy_active_total[1]+energy_active_total[2]+energy_active_total[3]+energy_inactive_total[0]+energy_inactive_total[1]+energy_inactive_total[2]+energy_inactive_total[3]):mc_eta>>SF_LAr_barrel_calibhit(56,0,1.4,70,0.18,0.25)","","colz"); + TH2* SF_LAr_barrel_calibhit=(TH2*)gDirectory->Get("SF_LAr_barrel_calibhit"); + //ProfileX(const char* name = "_pfx", Int_t firstybin = 1, Int_t lastybin = -1, Option_t* option = "") + TProfile* SF_prof_calibhit=SF_LAr_barrel_calibhit->ProfileX(); + SF_prof_calibhit->SetLineColor(2); + SF_prof_calibhit->Draw("same"); + //TFitResultPtr* fitres= + SF_prof_calibhit->Fit("pol0","","",0.1,0.75); + SF_prof_calibhit->Fit("pol0","+","",0.85,1.35); + + new TCanvas("SF_LAr_barrel","SF_LAr_barrel"); + mytree->Draw("(energy_hit[0]+energy_hit[1]+energy_hit[2]+energy_hit[3])/(energy_active_total[0]+energy_active_total[1]+energy_active_total[2]+energy_active_total[3]+energy_inactive_total[0]+energy_inactive_total[1]+energy_inactive_total[2]+energy_inactive_total[3]):mc_eta>>SF_LAr_barrel(56,0,1.4,70,0.16,0.23)","","colz"); + TH2* SF_LAr_barrel=(TH2*)gDirectory->Get("SF_LAr_barrel"); + //ProfileX(const char* name = "_pfx", Int_t firstybin = 1, Int_t lastybin = -1, Option_t* option = "") + TProfile* SF_prof=SF_LAr_barrel->ProfileX(); + SF_prof->SetLineColor(2); + SF_prof->Draw("same"); + //TFitResultPtr* fitres= + SF_prof->Fit("pol0","","",0.1,0.75); + SF_prof->Fit("pol0","+","",0.85,1.35); + + //mytree->Print(); + +} diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/share/LarEMSamplingFraction_topOptions.py b/Simulation/Tools/CaloSamplingFractionAnalysis/share/LarEMSamplingFraction_topOptions.py new file mode 100644 index 0000000000000000000000000000000000000000..634e43b458caa261cb483d083ae07dcf9b201ab3 --- /dev/null +++ b/Simulation/Tools/CaloSamplingFractionAnalysis/share/LarEMSamplingFraction_topOptions.py @@ -0,0 +1,79 @@ +import AthenaCommon.Constants as Lvl +from AthenaCommon.AppMgr import ServiceMgr as svcMgr +from AthenaCommon.AppMgr import theApp + +from AthenaCommon.AthenaCommonFlags import athenaCommonFlags +athenaCommonFlags.FilesInput = inFileName + +## load POOL support +import AthenaPoolCnvSvc.ReadAthenaPool + +## general job configuration +from AthenaCommon.AlgSequence import AlgSequence +topSequence = AlgSequence() + +from AthenaCommon.GlobalFlags import jobproperties +jobproperties.Global.DetDescrVersion="ATLAS-R2-2016-01-00-01" + +from RecExConfig.ObjKeyStore import ObjKeyStore, objKeyStore +oks = ObjKeyStore() +#oks.addStreamESD('CaloCellContainer', ['AllCalo'] ) +oks.addStreamESD('CaloCalibrationHitContainer',['LArCalibrationHitActive','LArCalibrationHitInactive','LArCalibrationHitDeadMaterial_DEAD','TileCalibHitActiveCell','TileCalibHitInactiveCell','TileCalibHitDeadMaterial']) + +## re-do topo clusters on EM scale +#from CaloRec.CaloRecFlags import jobproperties +#from CaloRec.CaloTopoClusterFlags import jobproperties +#jobproperties.CaloTopoClusterFlags.doTopoClusterLocalCalib = True +#jobproperties.CaloTopoClusterFlags.doCellWeightCalib = False +#from CaloRec.CaloClusterTopoGetter import CaloClusterTopoGetter, addSnapshot +#CaloClusterTopoGetter() + +from AthenaCommon.DetFlags import DetFlags +DetFlags.all_setOff() +DetFlags.LAr_setOn() +DetFlags.Tile_setOn() +DetFlags.digitize.all_setOff() + +from AthenaCommon.GlobalFlags import globalflags +globalflags.DetGeo.set_Value_and_Lock('atlas') +globalflags.DataSource.set_Value_and_Lock('geant4') + +from AtlasGeoModel import SetGeometryVersion +from AtlasGeoModel import GeoModelInit + +include( "CaloDetMgrDetDescrCnv/CaloDetMgrDetDescrCnv_joboptions.py") +include( "CaloIdCnv/CaloIdCnv_joboptions.py" ) +include( "TileIdCnv/TileIdCnv_jobOptions.py" ) +include( "LArDetDescr/LArDetDescr_joboptions.py" ) +include("TileConditions/TileConditions_jobOptions.py" ) +include("LArConditionsCommon/LArConditionsCommon_MC_jobOptions.py") + +svcMgr.IOVDbSvc.GlobalTag = "OFLCOND-MC16-SDR-16" +svcMgr.AthenaSealSvc.OutputLevel = Lvl.ERROR +svcMgr.IOVDbSvc.forceRunNumber=284500 + +svcMgr += CfgMgr.THistSvc() +svcMgr.THistSvc.Output += ["MYSTREAM DATAFILE='LArEM_SF.root' OPT='RECREATE'"] + +#from CaloLocalHadCalib.CaloLocalHadCalibConf import GetLCWeights +from AthenaCommon.SystemOfUnits import deg, GeV, MeV, TeV +from AthenaCommon.AlgSequence import AlgSequence +#from math import ( pi as m_pi, log10 as m_log10 ) + +LarEMSamplingFraction = CfgMgr.LarEMSamplingFraction() +LarEMSamplingFraction.DoCells = 0 +LarEMSamplingFraction.CalibrationHitContainerNames = ["LArCalibrationHitInactive" + ,"LArCalibrationHitActive" + ,"TileCalibHitInactiveCell" + ,"TileCalibHitActiveCell"] + #,"LArCalibrationHitDeadMaterial_DEAD" + #,"TileCalibHitDeadMaterial"] +topSequence += LarEMSamplingFraction + +theApp.EvtMax = -1 + +MessageSvc = svcMgr.MessageSvc +MessageSvc.OutputLevel = Lvl.INFO #2=DEBUG, 3=INFO, 4=WARNING, 5=ERROR, 6=FATAL +MessageSvc.infoLimit = 100000 + +svcMgr.EventSelector.InputCollections = inFileName diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/src/LarEMSamplingFraction.cxx b/Simulation/Tools/CaloSamplingFractionAnalysis/src/LarEMSamplingFraction.cxx new file mode 100644 index 0000000000000000000000000000000000000000..752b1237f7a09165c30b952a3d97fc0e82dc6210 --- /dev/null +++ b/Simulation/Tools/CaloSamplingFractionAnalysis/src/LarEMSamplingFraction.cxx @@ -0,0 +1,482 @@ +#include "LarEMSamplingFraction.h" + +//#include "xAODCaloEvent/CaloClusterContainer.h" +//#include "CaloEvent/CaloCell.h" +//#include "CaloUtils/CaloSamplingHelper.h" +#include "CaloSimEvent/CaloCalibrationHit.h" +#include "CaloSimEvent/CaloCalibrationHitContainer.h" +#include "CaloEvent/CaloCellContainer.h" +//#include "CaloDetDescr/CaloDetDescrManager.h" +//#include "CaloIdentifier/CaloCell_ID.h" +#include "LArSimEvent/LArHitContainer.h" +#include "TileSimEvent/TileHitVector.h" + +#include "GeoAdaptors/GeoLArHit.h" + +#include "GaudiKernel/ServiceHandle.h" +#include "GaudiKernel/ITHistSvc.h" +#include "GaudiKernel/Property.h" + +#include "GeneratorObjects/McEventCollection.h" + +#include "TString.h" +#include <iterator> +#include <cmath> +#include <map> + +using namespace std; + +//############################################################################### +LarEMSamplingFraction::LarEMSamplingFraction(const std::string& name, ISvcLocator* pSvcLocator) : AthAlgorithm(name, pSvcLocator), + m_docells(true), + m_calo_dd_man(0), + m_calo_id(0) +{ + + // Name of ClusterContainer to use + declareProperty("DoCells",m_docells); + declareProperty("CalibrationHitContainerNames",m_CalibrationHitContainerNames); +} + +//############################################################################### + +LarEMSamplingFraction::~LarEMSamplingFraction() +{ } + +//############################################################################### + +StatusCode LarEMSamplingFraction::initialize() +{ + //---- initialize the StoreGateSvc ptr ---------------- + + cout<<"CHECKME initialize()"<<endl; + + ServiceHandle<ITHistSvc> histSvc("THistSvc",name()); + ATH_CHECK( histSvc.retrieve() ); + + mytree = new TTree("mytree","mytree"); + mytree->Branch("energy_reco", &energy_reco); + mytree->Branch("energy_hit", &energy_hit); + mytree->Branch("energy_inactive_total",&energy_inactive_total); + mytree->Branch("energy_inactive_em", &energy_inactive_em); + mytree->Branch("energy_inactive_nonem",&energy_inactive_nonem); + mytree->Branch("energy_inactive_inv", &energy_inactive_inv); + mytree->Branch("energy_inactive_esc", &energy_inactive_esc); + mytree->Branch("energy_active_total_corrected",&energy_active_total_corrected); + mytree->Branch("energy_active_total",&energy_active_total); + mytree->Branch("energy_active_em", &energy_active_em); + mytree->Branch("energy_active_nonem",&energy_active_nonem); + mytree->Branch("energy_active_inv", &energy_active_inv); + mytree->Branch("energy_active_esc", &energy_active_esc); + mytree->Branch("mc_pdg", &mc_pdg); + mytree->Branch("mc_eta", &mc_eta); + mytree->Branch("mc_phi", &mc_phi); + mytree->Branch("mc_e", &mc_e); + mytree->Branch("mc_pt", &mc_pt); + + if(m_docells) { + mytree->Branch("cell_identifier",&m_cell_identifier); + mytree->Branch("cell_energy_reco",&m_cell_energy_reco); + mytree->Branch("cell_energy_inactive_total",&m_cell_energy_inactive_total); + mytree->Branch("cell_energy_active_total_corrected",&m_cell_energy_active_total_corrected); + mytree->Branch("cell_energy_active_total",&m_cell_energy_active_total); + mytree->Branch("cell_sampling",&m_cell_sampling); + mytree->Branch("cell_eta",&m_cell_eta); + mytree->Branch("cell_phi",&m_cell_phi); + } + + histSvc->regTree("/MYSTREAM/myTree",mytree).ignore(); + + // pointer to detector manager: + m_calo_dd_man = CaloDetDescrManager::instance(); + m_calo_id = m_calo_dd_man->getCaloCell_ID(); + + const DataHandle<CaloIdManager> caloIdManager; + StatusCode sc=detStore()->retrieve(caloIdManager); + if(sc.isSuccess()) + ATH_MSG_DEBUG("CaloIDManager retrieved."); + else + throw std::runtime_error("ISF_HitAnalysis: Unable to retrieve CaloIDManeger"); + + m_tileID=caloIdManager->getTileID(); + if(m_tileID==0) throw std::runtime_error("ISF_HitAnalysis: Invalid Tile ID helper"); + + sc=detStore()->regHandle(m_dd_fSampl,"LArfSampl"); + if(sc.isFailure()) + { + ATH_MSG_ERROR("Unable to register handle for LArfSampl"); + return StatusCode::FAILURE; + } + + detStore()->retrieve(m_tileInfo,"TileInfo"); + if(sc.isFailure()) + { + ATH_MSG_ERROR("Unable to retrieve TileInfo from DetectorStore"); + return StatusCode::FAILURE; + } + + + return StatusCode::SUCCESS; +} + +//############################################################################### + +StatusCode LarEMSamplingFraction::finalize() +{ + + return StatusCode::SUCCESS; +} + +//############################################################################### + +StatusCode LarEMSamplingFraction::execute() +{ + + /* + const DataHandle<xAOD::CaloClusterContainer> cc ; + StatusCode sc = evtStore()->retrieve(cc,m_clusterCollName); + + if(sc != StatusCode::SUCCESS) + { + msg(MSG::ERROR) << "Could not retrieve ClusterContainer " << m_clusterCollName << " from StoreGate" << endmsg; + return sc; + } + */ + + const DataHandle<CaloCalibrationHitContainer> cchc; + std::vector<const CaloCalibrationHitContainer *> v_cchc; + std::vector<std::string>::iterator iter; + for (iter=m_CalibrationHitContainerNames.begin();iter!=m_CalibrationHitContainerNames.end();iter++) + { + if ( !evtStore()->contains<CaloCalibrationHitContainer>(*iter)) + { + msg(MSG::ERROR) << "SG does not contain calibration hit container " << *iter << endmsg; + return StatusCode::FAILURE; + } + else + { + StatusCode sc = evtStore()->retrieve(cchc,*iter); + if (sc.isFailure() ) + { + msg(MSG::ERROR) << "Cannot retrieve calibration hit container " << *iter << endmsg; + return sc; + } + else + { + v_cchc.push_back(cchc); + //cout<<"CHECKME successfully retrieved container "<<*iter<<endl; + } + } + } + + const McEventCollection* truthEvent=0; + StatusCode sc = evtStore()->retrieve(truthEvent, "TruthEvent"); + if (sc.isFailure()||!truthEvent) + { + msg(MSG::ERROR) << "No McEventCollection found"<< endmsg; + return StatusCode::FAILURE; + } + HepMC::GenEvent::particle_const_iterator pit = truthEvent->at(0)->particles_begin(); + const HepMC::GenParticle * gen = *pit; + mc_pdg = gen->pdg_id(); + mc_eta = gen->momentum().pseudoRapidity(); + mc_phi = gen->momentum().phi(); + mc_e = gen->momentum().e(); + mc_pt = sqrt(pow(gen->momentum().px(),2)+pow(gen->momentum().py(),2)); + + //inspiration: + //see https://gitlab.cern.ch/atlas/athena/blob/master/Calorimeter/CaloCalibHitRec/src/CalibHitToCaloCell.cxx + //and https://gitlab.cern.ch/atlas/athena/blob/master/Calorimeter/CaloCalibHitRec/src/CaloDmEnergy.cxx + + const CaloDetDescrManager* caloDDMgr = nullptr; + ATH_CHECK(detStore()->retrieve(caloDDMgr, "CaloMgr")); + + energy_reco =new vector<float>; + energy_hit =new vector<float>; + + energy_inactive_total=new vector<float>; + energy_inactive_em =new vector<float>; + energy_inactive_nonem=new vector<float>; + energy_inactive_inv =new vector<float>; + energy_inactive_esc =new vector<float>; + + energy_active_total_corrected=new vector<float>; + energy_active_total=new vector<float>; + energy_active_em =new vector<float>; + energy_active_nonem=new vector<float>; + energy_active_inv =new vector<float>; + energy_active_esc =new vector<float>; + + if(m_docells) { + m_cell_identifier =new std::vector<Long64_t>; + m_cell_energy_reco =new std::vector<float>; + m_cell_energy_active_total_corrected=new std::vector<float>; + m_cell_energy_active_total =new std::vector<float>; + m_cell_energy_inactive_total =new std::vector<float>; + m_cell_sampling =new std::vector<int>; + m_cell_eta =new std::vector<float>; + m_cell_phi =new std::vector<float>; + } + + struct cell_info { + Long64_t cell_identifier=0; + float cell_energy_reco=0; + float cell_energy_active_total_corrected=0; + float cell_energy_active_total=0; + float cell_energy_inactive_total=0; + int cell_sampling=0; + float cell_eta=0; + float cell_phi=0; + }; + + std::map< Long64_t , cell_info > cell_info_map; + + for(int s=0;s<24;s++) + { + energy_reco->push_back(0.0); + energy_hit->push_back(0.0); + + energy_inactive_total->push_back(0.0); + energy_inactive_em ->push_back(0.0); + energy_inactive_nonem->push_back(0.0); + energy_inactive_inv ->push_back(0.0); + energy_inactive_esc ->push_back(0.0); + + energy_active_total_corrected->push_back(0.0); + energy_active_total->push_back(0.0); + energy_active_em ->push_back(0.0); + energy_active_nonem->push_back(0.0); + energy_active_inv ->push_back(0.0); + energy_active_esc ->push_back(0.0); + } + + std::vector<const CaloCalibrationHitContainer * >::const_iterator it; + int count=0; + for (it=v_cchc.begin();it!=v_cchc.end();it++) + { + //cout<<" this is a loop on "<<m_CalibrationHitContainerNames[count]<<endl; + + CaloCalibrationHitContainer::const_iterator chIter = (*it)->begin(); + CaloCalibrationHitContainer::const_iterator chIterE = (*it)->end(); + + for(;chIter!=chIterE;chIter++) + { + Identifier id=(*chIter)->cellID(); + + double Etot = (*chIter)->energyTotal(); + double Eem = (*chIter)->energyEM(); + double Enonem = (*chIter)->energyNonEM(); + //double Evis = Eem + Enonem; + double Einv = (*chIter)->energyInvisible(); + double Eesc = (*chIter)->energyEscaped(); + + double Efactor=1.0; + + const CaloDetDescrElement* caloDDE = caloDDMgr->get_element(id); + int sampling=-1; + if(caloDDE) { + sampling = caloDDE->getSampling(); + + if((sampling>=0 && sampling<=11) || (sampling>=21 && sampling<=23)) Efactor=1/m_dd_fSampl->FSAMPL(id); + if((sampling>=12 && sampling<=20)) { + Identifier cell_id = m_tileID->cell_id(id); + if(m_calo_dd_man->get_element(cell_id)) { + Efactor = m_tileInfo->HitCalib(cell_id); + id=cell_id; + } else { + Efactor = m_tileInfo->HitCalib(id); + } + } + } + + //cout<<" cellID "<<id<<" layer "<<sampling<<" energyTotal "<<Etot<<" Eem "<<Eem<<" Enonem "<<Enonem<<" Einv "<<Einv<<" Eesc "<<Eesc<<" Efactor="<<Efactor<<endl; + + if(sampling>=0 && sampling<=23) + { + if(m_docells) { + cell_info_map[id.get_compact()].cell_identifier=id.get_compact(); + cell_info_map[id.get_compact()].cell_sampling=sampling; + cell_info_map[id.get_compact()].cell_eta=caloDDE->eta_raw(); + cell_info_map[id.get_compact()].cell_phi=caloDDE->phi_raw(); + } + + if(m_CalibrationHitContainerNames[count]=="LArCalibrationHitInactive" || m_CalibrationHitContainerNames[count]=="TileCalibHitInactiveCell") + { + energy_inactive_total->at(sampling)+=Etot; + energy_inactive_em ->at(sampling)+=Eem; + energy_inactive_nonem->at(sampling)+=Enonem; + energy_inactive_inv ->at(sampling)+=Einv; + energy_inactive_esc ->at(sampling)+=Eesc; + + if(m_docells) cell_info_map[id.get_compact()].cell_energy_inactive_total+=Etot; + } + + if(m_CalibrationHitContainerNames[count]=="LArCalibrationHitActive" || m_CalibrationHitContainerNames[count]=="TileCalibHitActiveCell") + { + energy_active_total_corrected->at(sampling)+=Etot*Efactor; + energy_active_total->at(sampling)+=Etot; + energy_active_em ->at(sampling)+=Eem; + energy_active_nonem->at(sampling)+=Enonem; + energy_active_inv ->at(sampling)+=Einv; + energy_active_esc ->at(sampling)+=Eesc; + + if(m_docells) { + cell_info_map[id.get_compact()].cell_energy_active_total_corrected+=Etot*Efactor; + cell_info_map[id.get_compact()].cell_energy_active_total+=Etot; + } + } + + } + + } + + + + count++; + } + + //Get reco cells if available + const CaloCellContainer *cellColl = 0; + sc = evtStore()->retrieve(cellColl, "AllCalo"); + + if (sc.isFailure()) { + ATH_MSG_WARNING( "Couldn't read AllCalo cells from StoreGate"); + //return NULL; + } else { + ATH_MSG_DEBUG( "Found: "<<cellColl->size()<<" calorimeter cells"); + CaloCellContainer::const_iterator itrCell = cellColl->begin(); + CaloCellContainer::const_iterator itrLastCell = cellColl->end(); + for ( ; itrCell!=itrLastCell; ++itrCell) { + Identifier id=(*itrCell)->ID(); + const CaloDetDescrElement* caloDDE = caloDDMgr->get_element(id); + int sampling=-1; + if(caloDDE) { + sampling = caloDDE->getSampling(); + energy_reco->at(sampling)+=(*itrCell)->energy(); + if((sampling>=12 && sampling<=20)) { + Identifier cell_id = m_tileID->cell_id(id); + if(m_calo_dd_man->get_element(cell_id)) { + id=cell_id; + } + } + } + if(m_docells) { + cell_info_map[id.get_compact()].cell_identifier=id.get_compact(); + cell_info_map[id.get_compact()].cell_sampling=sampling; + cell_info_map[id.get_compact()].cell_eta=caloDDE->eta_raw(); + cell_info_map[id.get_compact()].cell_phi=caloDDE->phi_raw(); + cell_info_map[id.get_compact()].cell_energy_reco+=(*itrCell)->energy(); + } + } + } //calorimeter cells + + + //Get all G4Hits (from CaloHitAnalysis) + std::string lArKey [4] = {"LArHitEMB", "LArHitEMEC", "LArHitFCAL", "LArHitHEC"}; + for (unsigned int i=0;i<4;i++) + { + const DataHandle<LArHitContainer> iter; + ATH_MSG_DEBUG( "Checking G4Hits: "<<lArKey[i]); + if(evtStore()->retrieve(iter,lArKey[i])==StatusCode::SUCCESS) + { + LArHitContainer::const_iterator hi; + int hitnumber = 0; + for (hi=(*iter).begin();hi!=(*iter).end();hi++) + { + hitnumber++; + GeoLArHit ghit(**hi); + if (!ghit) continue; + const CaloDetDescrElement *hitElement = ghit.getDetDescrElement(); + if(!hitElement) continue; + Identifier larhitid = hitElement->identify(); + if(m_calo_dd_man->get_element(larhitid)) + { + CaloCell_ID::CaloSample larlayer = m_calo_dd_man->get_element(larhitid)->getSampling(); + energy_hit->at(larlayer)+=ghit.Energy(); + } + } // End while LAr hits + ATH_MSG_DEBUG( "Read "<<hitnumber<<" G4Hits from "<<lArKey[i]); + } + else + { + ATH_MSG_INFO( "Can't retrieve LAr hits"); + }// End statuscode success upon retrieval of hits + }// End detector type loop + + const TileHitVector * hitVec; + if (evtStore()->retrieve(hitVec,"TileHitVec")==StatusCode::SUCCESS && m_tileID ) + { + int hitnumber = 0; + for(TileHitVecConstIterator i_hit=hitVec->begin() ; i_hit!=hitVec->end() ; ++i_hit) + { + hitnumber++; + Identifier pmt_id = (*i_hit).identify(); + Identifier cell_id = m_tileID->cell_id(pmt_id); + //const CaloDetDescrElement* ddElement = m_tileMgr->get_cell_element(cell_id); + + if (m_calo_dd_man->get_element(cell_id)) + { + CaloCell_ID::CaloSample layer = m_calo_dd_man->get_element(cell_id)->getSampling(); + + //could there be more subhits?? + for (int tilesubhit_i = 0; tilesubhit_i<(*i_hit).size(); tilesubhit_i++) + { + //!! + //std::cout <<"Tile subhit: "<<tilesubhit_i<<"/"<<(*i_hit).size()<< " E: "<<(*i_hit).energy(tilesubhit_i)<<std::endl; + energy_hit->at(layer) += (*i_hit).energy(tilesubhit_i); + } + } + } + ATH_MSG_INFO( "Read "<<hitnumber<<" G4Hits from TileHitVec"); + } + + + + + + + + + + + for(auto& cell:cell_info_map) { + m_cell_identifier ->push_back(cell.second.cell_identifier); + m_cell_sampling ->push_back(cell.second.cell_sampling); + m_cell_eta ->push_back(cell.second.cell_eta); + m_cell_phi ->push_back(cell.second.cell_phi); + m_cell_energy_reco ->push_back(cell.second.cell_energy_reco); + m_cell_energy_active_total_corrected->push_back(cell.second.cell_energy_active_total_corrected); + m_cell_energy_active_total ->push_back(cell.second.cell_energy_active_total); + m_cell_energy_inactive_total ->push_back(cell.second.cell_energy_inactive_total); + } + + mytree->Fill(); + + delete energy_reco; + delete energy_hit; + delete energy_inactive_total; + delete energy_inactive_em; + delete energy_inactive_nonem; + delete energy_inactive_inv; + delete energy_inactive_esc; + + delete energy_active_total_corrected; + delete energy_active_total; + delete energy_active_em; + delete energy_active_nonem; + delete energy_active_inv; + delete energy_active_esc; + + if(m_docells) { + delete m_cell_identifier; + delete m_cell_energy_reco; + delete m_cell_energy_active_total_corrected; + delete m_cell_energy_active_total; + delete m_cell_energy_inactive_total; + delete m_cell_sampling; + delete m_cell_eta; + delete m_cell_phi; + } + + return StatusCode::SUCCESS; +} diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/src/LarEMSamplingFraction.h b/Simulation/Tools/CaloSamplingFractionAnalysis/src/LarEMSamplingFraction.h new file mode 100644 index 0000000000000000000000000000000000000000..359ff38bccaa3844a2bb25d7d4720fbb030434e7 --- /dev/null +++ b/Simulation/Tools/CaloSamplingFractionAnalysis/src/LarEMSamplingFraction.h @@ -0,0 +1,84 @@ + +#ifndef LarEMSamplingFraction_H +#define LarEMSamplingFraction_H 1 + +//#include "GaudiKernel/HistoDef.h" +#include "AthenaBaseComps/AthAlgorithm.h" +#include "LArElecCalib/ILArfSampl.h" +#include "TileConditions/TileInfo.h" +#include "CaloIdentifier/TileID.h" +#include "CaloIdentifier/CaloIdManager.h" + +#include "StoreGate/StoreGateSvc.h" +#include "CaloIdentifier/CaloCell_ID.h" +#include "CaloDetDescr/CaloDetDescrManager.h" + +#include <vector> +#include <string> +#include "TTree.h" + +using namespace std; + +class LarEMSamplingFraction : public ::AthAlgorithm +{ + + public: + + LarEMSamplingFraction(const std::string& name, ISvcLocator* pSvcLocator); + virtual ~LarEMSamplingFraction(); + virtual StatusCode initialize(); + virtual StatusCode execute(); + virtual StatusCode finalize(); + + private: + + /// Default constructor: + LarEMSamplingFraction(); + + bool m_docells; + + int mc_pdg; + double mc_eta; + double mc_phi; + double mc_e; + double mc_pt; + + std::vector<Long64_t>* m_cell_identifier; + std::vector<float>* m_cell_energy_reco; + std::vector<float>* m_cell_energy_active_total_corrected; + std::vector<float>* m_cell_energy_active_total; + std::vector<float>* m_cell_energy_inactive_total; + std::vector<int>* m_cell_sampling; + std::vector<float>* m_cell_eta; + std::vector<float>* m_cell_phi; + + vector<float> *energy_reco; + vector<float> *energy_hit; + + vector<float> *energy_inactive_total; + vector<float> *energy_inactive_em; + vector<float> *energy_inactive_nonem; + vector<float> *energy_inactive_inv; + vector<float> *energy_inactive_esc; + + vector<float> *energy_active_total_corrected; + vector<float> *energy_active_total; + vector<float> *energy_active_em; + vector<float> *energy_active_nonem; + vector<float> *energy_active_inv; + vector<float> *energy_active_esc; + + TTree* mytree =0; + + vector<string> m_CalibrationHitContainerNames; + + const CaloDetDescrManager* m_calo_dd_man; + + const CaloCell_ID* m_calo_id; + + const TileInfo *m_tileInfo; + const TileID * m_tileID; + const DataHandle<ILArfSampl> m_dd_fSampl; +}; + +#endif //> !LarEMSamplingFraction_H diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/src/components/CaloSamplingFractionAnalysis_entries.cxx b/Simulation/Tools/CaloSamplingFractionAnalysis/src/components/CaloSamplingFractionAnalysis_entries.cxx new file mode 100755 index 0000000000000000000000000000000000000000..0e53abcbf6445ab2180841457ab44adc22660110 --- /dev/null +++ b/Simulation/Tools/CaloSamplingFractionAnalysis/src/components/CaloSamplingFractionAnalysis_entries.cxx @@ -0,0 +1,10 @@ +#include "GaudiKernel/DeclareFactoryEntries.h" + +#include "../LarEMSamplingFraction.h" + +DECLARE_ALGORITHM_FACTORY( LarEMSamplingFraction ) + +DECLARE_FACTORY_ENTRIES( CaloSamplingFractionAnalysis ) { + DECLARE_ALGORITHM( LarEMSamplingFraction ) +} + diff --git a/Simulation/Tools/CaloSamplingFractionAnalysis/src/components/CaloSamplingFractionAnalysis_load.cxx b/Simulation/Tools/CaloSamplingFractionAnalysis/src/components/CaloSamplingFractionAnalysis_load.cxx new file mode 100755 index 0000000000000000000000000000000000000000..10397443814f579ca9ef861f890a1ac970a8ec40 --- /dev/null +++ b/Simulation/Tools/CaloSamplingFractionAnalysis/src/components/CaloSamplingFractionAnalysis_load.cxx @@ -0,0 +1,3 @@ +#include "GaudiKernel/LoadFactoryEntries.h" +LOAD_FACTORY_ENTRIES(CaloSamplingFractionAnalysis) +