00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include "AtlfastAlgs/CellMaker.h"
00012 #include "AtlfastAlgs/CellsAboveThreshold.h"
00013 #include "AtlfastAlgs/TransportedParticle.h"
00014 #include "AtlfastAlgs/TransportedParticleCollection.h"
00015 #include "AtlfastAlgs/Calorimeter.h"
00016 #include "AtlfastAlgs/MagField.h"
00017 #include "AtlfastAlgs/ICellSelector.h"
00018 #include "AtlfastAlgs/ISmearer.h"
00019 #include "AtlfastAlgs/CellSmearer.h"
00020 #include "AtlfastAlgs/GlobalEventData.h"
00021 #include "AtlfastAlgs/EPileupMap.h"
00022
00023 #include "AtlfastEvent/CollectionDefs.h"
00024
00025 #include "AtlfastUtils/TesIO.h"
00026 #include "AtlfastUtils/HeaderPrinter.h"
00027
00028 #include <cmath>
00029 #include <iomanip>
00030 #include <vector>
00031
00032 #include "GaudiKernel/DataSvc.h"
00033 #include "GaudiKernel/Chrono.h"
00034 #include "StoreGate/DataHandle.h"
00035 #include "GaudiKernel/ISvcLocator.h"
00036 #include "GaudiKernel/MsgStream.h"
00037
00038
00039
00040 namespace Atlfast {
00041 CellMaker:: CellMaker(const std::string& name, ISvcLocator* pSvcLocator):
00042 Algorithm(name,pSvcLocator){
00043
00044
00045 m_inputLocation = "/Event/McEventCollection";
00046 m_outputLocation = "/Event/AtlfastCells" ;
00047
00048
00049 m_etaCoverage = 5.0;
00050 m_minETCell = 0.0;
00051 m_granBarrelEta = 0.1;
00052 m_granBarrelPhi = 0.1;
00053 m_granForwardEta = 0.2;
00054 m_granForwardPhi = 0.2;
00057 m_doSmearing = false;
00058
00059 declareProperty("InputLocation", m_inputLocation);
00060 declareProperty("OutputLocation", m_outputLocation);
00061 declareProperty("EtaCoverage", m_etaCoverage);
00062 declareProperty("MinETCell", m_minETCell);
00063 declareProperty("GranBarrelEta", m_granBarrelEta);
00064 declareProperty("GranBarrelPhi", m_granBarrelPhi);
00065 declareProperty("GranForwardEta", m_granForwardEta);
00066 declareProperty("GranForwardPhi", m_granForwardPhi);
00068 declareProperty("DoSmearing", m_doSmearing);
00069
00070 }
00071
00072
00073 CellMaker::~CellMaker()
00074 {
00075 if (m_calorimeter) {
00076 delete m_calorimeter;
00077 }
00078 if (m_tesIO) {
00079 delete m_tesIO;
00080 }
00081 if (m_cellSelector) {
00082 delete m_cellSelector;
00083 }
00084 if (m_magField) {
00085 delete m_magField;
00086 }
00087 }
00088
00089 StatusCode CellMaker::initialize()
00090 {
00091
00092 MsgStream log( messageService(), name() );
00093
00094
00095
00096 m_cellSelector = new CellsAboveThreshold(m_minETCell);
00097
00098
00099 m_tesIO = new TesIO();
00100
00101
00102 GlobalEventData* ged = GlobalEventData::Instance();
00103 m_fieldOn = ged->fieldOn();
00104 m_mcSelector = ged -> visibleToCal();
00105 m_barrelForwardEta = ged -> barrelForwardEta();
00106 m_magField = new MagField(m_fieldOn);
00107 if (fabs(m_etaCoverage) < fabs(m_barrelForwardEta) ) {
00108 log << MSG::ERROR
00109 << "asked to makecalo with max extent "
00110 << m_etaCoverage
00111 << "but barrel-forward boundary at "
00112 << m_barrelForwardEta << endreq;
00113 return StatusCode::FAILURE;
00114 }
00115 m_calorimeter = new Calorimeter(
00116 log,
00117 m_etaCoverage,
00118 m_barrelForwardEta,
00119 m_granBarrelEta,
00120 m_granBarrelPhi,
00121 m_granForwardEta,
00122 m_granForwardPhi
00123 );
00127 if(m_doSmearing == true){
00128 m_epileupMap = new EPileupMap(ged->lumi());
00129 m_smearer = new CellSmearer(ged->randSeed(),m_barrelForwardEta);
00130 }else{
00131 m_epileupMap = NULL;
00132 m_smearer = NULL;
00133 }
00134
00135 HeaderPrinter hp("Atlfast Cell Maker:", log);
00136
00137 hp.add("min Cell ET ", m_minETCell);
00138 hp.add("Magnetic Field is on? ", m_fieldOn);
00139 hp.add("Input Location ", m_inputLocation);
00140 hp.add("Output Location ", m_outputLocation);
00141 hp.add(" DoSmearing ", m_doSmearing);
00142 hp.add("Calorimeter Geometry: ");
00143 hp.add(" Total Eta Range ", m_etaCoverage);
00144 hp.add(" Barrel Eta Range ", m_barrelForwardEta);
00145 hp.add(" Barrel Eta Granularity ", m_granBarrelEta);
00146 hp.add(" Barrel Phi Granularity ", m_granBarrelPhi);
00147 hp.add(" Forward Eta Granularity", m_granForwardEta);
00148 hp.add(" Forward Phi Granularity", m_granForwardPhi);
00149 hp.print();
00150
00151
00152
00153 return StatusCode::SUCCESS;
00154 }
00155
00156 StatusCode CellMaker::finalize()
00157 {
00158
00159 MsgStream log( messageService(), name() );
00160
00161 return StatusCode::SUCCESS;
00162 }
00163
00164
00165 StatusCode CellMaker::execute()
00166 {
00167
00168
00169 MsgStream log( messageService(), name() );
00170 log << MSG::DEBUG << "CellMaker execute()" << endreq;
00171
00172 std::string mess;
00173
00174
00175
00176 MCparticleCollection mcParticles ;
00177 TesIoStat stat = m_tesIO->getMC( mcParticles, m_mcSelector ) ;
00178 mess = stat? "Retrieved MC from TES ":"Failed MC retrieve from TES";
00179 log << MSG::DEBUG << mess << endreq;
00180
00181
00182 TransportedParticleCollection particlesAtCal;
00183 (*m_magField)(mcParticles, particlesAtCal);
00184 log << MSG::DEBUG << "bend" << endreq;
00185
00186
00187
00188 TransportedParticleCollectionIter fpart = particlesAtCal.begin();
00189 TransportedParticleCollectionIter lpart = particlesAtCal.end();
00190 m_calorimeter->deposit(fpart, lpart);
00191
00192
00193
00194 if(m_doSmearing){
00195
00196 m_calorimeter->smearCells(m_smearer);
00197 log << MSG::DEBUG << "Doing Smearing" << endreq;
00199 m_epileupMap->fillMap();
00200 log << MSG::DEBUG << "Size of EM EPileup vector: " << m_epileupMap->getEMMap()->size() << endreq;
00201 std::vector<Atlfast::EPileupDeposit*>::iterator fdep = m_epileupMap->getEMMap()->begin();
00202 std::vector<Atlfast::EPileupDeposit*>::iterator ldep = m_epileupMap->getEMMap()->end();
00204 m_calorimeter->deposit(fdep, ldep);
00205 log << MSG::DEBUG << "ECAL pileup deposited" << endreq;
00206 log << MSG::DEBUG << "Size of HAD EPileup vector: " << m_epileupMap->getHadMap()->size() << endreq;
00207 fdep = m_epileupMap->getHadMap()->begin();
00208 ldep = m_epileupMap->getHadMap()->end();
00210 m_calorimeter->deposit(fdep, ldep);
00211 log << MSG::DEBUG << "HCAL pileup deposited" << endreq;
00212 }
00213
00214 log << MSG::DEBUG << "deposit" << endreq;
00215
00216
00217
00218
00219 CellCollection* cells = new CellCollection;
00220
00223 m_calorimeter->giveHitCells(m_cellSelector, cells);
00224 log << MSG::DEBUG << "hit cells retrieved, " << cells->size()
00225 << " cells were hit this event " << endreq;
00226
00227
00228
00229
00230 m_calorimeter->reset();
00231 log << MSG::DEBUG << "calorimeter reset" << endreq;
00232
00233
00234
00235
00236 stat = m_tesIO->store(cells, m_outputLocation );
00237 mess = stat? " Stored cells to TES":" Failed store cells to TES";
00238 log << MSG::DEBUG << mess << endreq;
00239
00240 return stat;
00241 }
00242
00243 }
00244
00245
00246
00247
00248
00249
00250