 |
ATLAS Offline Software
|
Go to the documentation of this file.
33 #include "CLHEP/Random/Randomize.h"
34 #include "CLHEP/Random/RandomEngine.h"
41 using CLHEP::RandFlat;
42 using CLHEP::RandGaussQ;
43 using CLHEP::RandPoissonT;
44 using CLHEP::RandGeneral;
47 const std::string&
name,
77 ATH_MSG_INFO(
"No photostatistics effect will be simulated");
94 ATH_MSG_INFO(
" In case of pileup, the trigger time subtraction is done in PileUpSvc");
95 ATH_MSG_INFO(
" => TileHitVecToCnt will not apply Trigger Time ");
114 <<
"'; therefore set HitTimeFlag to 2");
127 ATH_MSG_INFO(
"Minimal hit time will be used as trigger time"
128 <<
" with random additional shift between 0 and " << -
m_triggerTime <<
" ns");
131 ATH_MSG_INFO(
"Average time will be calculated in every event");
139 ATH_MSG_INFO(
"Time of all hits will be preserved during copy");
177 ATH_MSG_DEBUG(
"TileHitVecToCntTool initialization completed");
180 return StatusCode::RECOVERABLE;
182 return StatusCode::SUCCESS;
199 for (std::unique_ptr<TileHit>& hit :
m_allHits) {
217 return StatusCode::SUCCESS;
223 ATH_MSG_DEBUG(
"TileHitVecToCntTool prepareEvent initialization started");
232 return StatusCode::SUCCESS;
240 for (; inpItr !=
end; ++inpItr) {
242 const TileHit * cinp = &(*inpItr);
244 eHitTot += cinp->
energy();
247 hits->push_back(pHit);
251 int hitsize = cinp->
size();
254 for (
int i = 0;
i < hitsize; ++
i) {
263 eHitTot += eHit - cinp->
energy();
269 <<
" Copy hit: ener=";
270 for (
int i = 0;
i < hitsize; ++
i)
273 for (
int i = 0;
i < hitsize; ++
i)
282 std::vector<std::unique_ptr<TileHit>>& allHits,
283 std::vector<std::unique_ptr<TileHit>>& allHits_DigiHSTruth,
284 int& nHit,
double& eHitTot,
bool isSignal)
const {
294 for (; inpItr !=
end; ++inpItr) {
296 const TileHit * cinp = &(*inpItr);
311 if (hit_idhash >= allHits.size()) {
316 double ener = cinp->
energy();
317 double time = cinp->
time() + SubEvtTimOffset;
322 std::unique_ptr<TileHit>& pHit = allHits[hit_idhash];
323 std::unique_ptr<TileHit> inValidHit(
nullptr);
324 std::unique_ptr<TileHit>& pHit_DigiHSTruth =
m_doDigiTruth ? allHits_DigiHSTruth[hit_idhash] : inValidHit;
338 if (pHit->
size() > 1 || pHit->
energy() != 0.0)
343 <<
" offs=" << SubEvtTimOffset
344 <<
" double hit" <<
endmsg;
351 <<
" offs=" << SubEvtTimOffset
355 int hitsize = cinp->
size();
356 for (
int ind = 1; ind < hitsize; ++ind) {
358 time = cinp->
time(ind) + SubEvtTimOffset;
378 <<
" double hit from single hit" <<
endmsg;
397 for (; inpItr !=
end; ++inpItr) {
398 const TileHit * cinp = &(*inpItr);
399 eHitTot += cinp->
energy();
406 int hitsize = cinp->
size();
409 for (
int i = 0;
i < hitsize; ++
i) {
418 eHitTot += eHit - cinp->
energy();
424 <<
" Copy hit: ener=";
425 for (
int i = 0;
i < hitsize; ++
i)
428 for (
int i = 0;
i < hitsize; ++
i)
438 for (; inpItr !=
end; ++inpItr) {
439 const TileHit * cinp = &(*inpItr);
442 for (
int i = 0;
i <
size; ++
i) {
455 int hitsize = cinp->
size();
460 <<
" Input hit: ener=";
461 for (
int i = 0;
i < hitsize; ++
i)
464 for (
int i = 0;
i < hitsize; ++
i)
483 avtime = cosTriggerTime->
time();
494 for (; inpItr !=
end; ++inpItr) {
495 const TileHit * cinp = &(*inpItr);
497 for (
int i = 0;
i <
size; ++
i) {
498 if (cinp->
time(
i) < avtime) avtime = cinp->
time(
i);
511 for (; inpItr !=
end; ++inpItr) {
512 const TileHit * cinp = &(*inpItr);
514 for (
int i = 0;
i <
size; ++
i) {
528 inpItr = inputHits->
begin();
531 for (; inpItr !=
end; ++inpItr) {
532 const TileHit * cinp = &(*inpItr);
536 for (
int i = 0;
i <
size; ++
i) {
545 int hitsize = pHit->
size();
548 for (
int i = 0;
i < hitsize; ++
i) {
561 <<
" Output hit: ener=";
562 for (
int i = 0;
i < hitsize; ++
i)
565 for (
int i = 0;
i < hitsize; ++
i)
586 ATH_MSG_DEBUG(
"Inside TileHitVecToCntTool processBunchXing" << bunchXing);
590 CLHEP::HepRandomEngine * engine = rngWrapper->
getEngine(Gaudi::Hive::currentContext());
595 while (iEvt != eSubEvents) {
600 std::vector<std::string>::const_iterator hitVecNamesItr =
m_hitVectorNames.begin();
601 std::vector<std::string>::const_iterator hitVecNamesEnd =
m_hitVectorNames.end();
602 for (; hitVecNamesItr != hitVecNamesEnd; ++hitVecNamesItr) {
604 const std::string hitVectorName(*hitVecNamesItr);
609 if (!(
m_mergeSvc->retrieveSingleSubEvtData(hitVectorName, inputHits, bunchXing, iEvt))){
610 ATH_MSG_ERROR(
" Tile Hit container not found for event key " << hitVectorName);
613 const double SubEvtTimOffset(iEvt->time());
616 if (fabs(SubEvtTimOffset) > 0.1) {
617 ATH_MSG_ERROR(
"Wrong time for in-time event: " << SubEvtTimOffset <<
" Ignoring all hits ");
619 ATH_MSG_DEBUG(
" New HitCont. TimeOffset=" << SubEvtTimOffset <<
", size =" << inputHits->
size());
624 bool isSignal =
false;
625 if(iEvt == bSubEvents) isSignal =
true;
633 if (!(
m_mergeSvc->retrieveSingleSubEvtData(hitVectorName, inputHits, bunchXing, iEvt))){
634 ATH_MSG_ERROR(
" Tile Hit container not found for event key " << hitVectorName);
647 ATH_MSG_DEBUG(
"Exiting processBunchXing in TileHitVecToCntTool");
649 return StatusCode::SUCCESS;
660 ATH_MSG_DEBUG(
"TileHitVecToCntTool processAllSubEvents started");
665 std::vector<std::unique_ptr<TileHit>> allHits;
666 auto hits = std::make_unique<TileHitNonConstContainer>(ownPolicy);
668 std::vector<std::unique_ptr<TileHit>> allHits_DigiHSTruth;
669 std::unique_ptr<TileHitNonConstContainer> hits_DigiHSTruth;
671 hits_DigiHSTruth = std::make_unique<TileHitNonConstContainer>(ownPolicy);
686 CLHEP::HepRandomEngine * engine = rngWrapper->
getEngine(ctx);
690 for (
auto & inputHits : hitVectorHandles) {
691 if (!inputHits.isValid()) {
693 return StatusCode::FAILURE;
695 const double SubEvtTimeOffset(0.0);
697 ATH_MSG_DEBUG(
" New HitCont. TimeOffset=" << SubEvtTimeOffset <<
", size =" << inputHits->size());
702 return StatusCode::SUCCESS;
708 TimedHitContList hitContList;
710 if (!(
m_mergeSvc->retrieveSubEvtsData(hitVectorName, hitContList).isSuccess()) || hitContList.size() == 0) {
711 ATH_MSG_WARNING(
"Could not fill TimedHitContList for hit vector " << hitVectorName);
720 if (iCont != iEndCont) {
722 const double SubEvtTimeOffset(iCont->first.time());
723 if (fabs(SubEvtTimeOffset) > 0.1) {
724 ATH_MSG_ERROR(
"Wrong time for in-time event: " << SubEvtTimeOffset <<
" Ignoring all hits ");
728 ATH_MSG_DEBUG(
" New HitCont. TimeOffset=" << SubEvtTimeOffset <<
", size =" << inputHits->
size());
735 for (; iCont != iEndCont; ++iCont) {
737 const double SubEvtTimeOffset(iCont->first.time());
740 ATH_MSG_VERBOSE(
" New HitCont. TimeOffset=" << SubEvtTimeOffset <<
", size =" << inputHits->
size());
741 bool isSignal =
false;
742 if(iCont == hitContList.begin() ) isSignal =
true;
753 ATH_MSG_WARNING(
"Hit Vector "<< hitVectorName <<
" not found in StoreGate");
768 return StatusCode::SUCCESS;
772 std::vector<std::unique_ptr<TileHit>>& allHits_DigiHSTruth,
773 std::unique_ptr<TileHitNonConstContainer>&
hits,
774 std::unique_ptr<TileHitNonConstContainer>& hits_DigiHSTruth,
777 std::vector<std::unique_ptr<TileHit>>
::iterator iHit = allHits.begin();
778 std::vector<std::unique_ptr<TileHit>>
::iterator lastHit = allHits.end();
779 std::vector<std::unique_ptr<TileHit>>
::iterator iHit_DigiHSTruth = allHits_DigiHSTruth.begin();
782 double eHitInTime = 0.0;
789 for (; iHit != lastHit; ++iHit) {
790 std::unique_ptr<TileHit>& hit = (*iHit);
791 if (hit->
size() > 1 || hit->
energy() != 0.0) {
792 eHitInTime += hit->
energy();
793 hits->push_back(getOrRelease(hit));
800 ATH_MSG_DEBUG(
" nHitUni=" << nHitUni <<
" eHitInTime="<< eHitInTime);
804 std::unique_ptr<TileHitNonConstContainer>&
hits,
805 std::unique_ptr<TileHitNonConstContainer>& hits_DigiHSTruth,
822 for (std::unique_ptr<TileHitCollection>& coll : *
hits ) {
823 int frag_id = coll->identify();
833 for (; collIt != endcollIt; ++collIt) {
834 int frag_id = (*collIt)->identify();
847 CLHEP::HepRandomEngine * engine = rngWrapper->
getEngine(ctx);
855 collIt_DigiHSTruth = hits_DigiHSTruth->
begin();
856 endColl_DigiHSTruth = hits_DigiHSTruth->
end();
859 for (std::unique_ptr<TileHitCollection>& coll : *
hits ) {
864 coll_DigiHSTruth = (*collIt_DigiHSTruth).
get();
865 hitItr_DigiHSTruth = coll_DigiHSTruth->
begin();
866 hitEnd_DigiHSTruth = coll_DigiHSTruth->
end();
876 int hitsize = pHit->size();
877 for (
int i = 0;
i < hitsize; ++
i) {
878 double thit = pHit->time(
i);
890 double scaleFactor = 1.0;
893 pHit->scale(scaleFactor);
897 TileHit *pHit_DigiHSTruth = (*hitItr_DigiHSTruth);
898 pHit_DigiHSTruth->
scale(scaleFactor);
900 ++hitItr_DigiHSTruth;
909 auto hitCont = std::make_unique<TileHitContainer>(
false, ownPolicy);
911 for (std::unique_ptr<TileHitCollection>& coll : *
hits ) {
912 CHECK(hitCont->addCollection (coll.release(), hashId++));
921 auto hitCont_DigiHSTruth = std::make_unique<TileHitContainer>(
false, ownPolicy);
922 size_t hashId_DigiHSTruth = 0;
923 for (std::unique_ptr<TileHitCollection>& coll : *hits_DigiHSTruth ) {
924 ATH_CHECK(hitCont_DigiHSTruth->addCollection (coll.release(), hashId_DigiHSTruth++));
928 ATH_CHECK( hitContainer_DigiHSTruth.
record(std::move(hitCont_DigiHSTruth)) );
932 return StatusCode::SUCCESS;
951 return StatusCode::SUCCESS;
965 nPhotoElectrons =
std::round(nPhotoElectrons * 1000) / 1000;
967 double pe =
energy * nPhotoElectrons;
968 double pe_scale = 1., RndmPois = 1.;
973 RndmPois =
std::max(0.0, RandGaussQ::shoot(engine,
pe, sqrt(
pe)));
974 pe_scale = RndmPois /
pe;
978 double singleMEAN = 1.0;
979 double singleSIGMA = 1.0;
980 RndmPois = RandPoissonT::shoot(engine,
pe);
984 for (
int i = 0;
i < RndmPois;
i++)
985 pe_scale += 1 / (1.08332) *
std::max(0., RandGaussQ::shoot(engine, singleMEAN, singleSIGMA));
987 pe_scale /= RndmPois;
996 RndmPois = RandPoissonT::shoot(engine,
pe);
997 pe_scale = RndmPois /
pe;
1004 RndmPois =
std::max(0.0, RandGaussQ::shoot(engine,
pe, sqrt(
pe)));
1007 double * ProbFunc =
new double[nn];
1008 ProbFunc[0] =
exp(-
pe);
1009 for (
int i = 1;
i < nn; ++
i) {
1010 ProbFunc[
i] = ProbFunc[
i - 1] *
pe /
i;
1012 RandGeneral* RandG =
new RandGeneral(ProbFunc, nn, 0);
1013 RndmPois = RandG->shoot(engine) * nn;
1019 pe_scale = RndmPois /
pe;
1026 <<
", numPhElec=" << nPhotoElectrons
1028 <<
", rndmPoisson=" << RndmPois
1029 <<
", pe_scale=" << pe_scale);
1036 int module = frag_id & 0x3F;
1044 for (; hitIt != endHitIt; ++hitIt) {
1055 if (fromHitIt != coll->
end()) {
1057 <<
"] in module: " <<
module);
1058 bool isToHitNew =
false;
1062 toHit =
new TileHit(to_pmt_id);
1077 toHit->
add(*fromHitIt, 0.1);
1084 coll->
erase(fromHitIt);
1094 int module = frag_id & 0x3F;
1102 for (; hitIt != endHitIt; ++hitIt) {
1113 if (fromHitIt != coll->
end()) {
1115 <<
"] in module: " <<
module);
1116 bool isToHitNew =
false;
1121 toHit =
new TileHit(to_pmt_id);
1136 toHit->
add(*fromHitIt, 0.1);
1144 coll->
erase(fromHitIt);
1153 for (std::unique_ptr<TileHitCollection>& coll : *hitCont) {
1155 int frag_id = coll->identify();
1156 int module = frag_id & 0x3F;
1158 std::vector<TileHit*>
hits(48,
nullptr);
1159 std::vector<std::unique_ptr<TileHit>> otherModuleHits;
1160 coll->erase(std::remove_if(coll->begin(), coll->end(),
1162 Identifier pmt_id = hit->pmt_ID();
1163 int channel = m_tileHWID->channel(hit->pmt_HWID());
1164 TileHit* channelHit = hits[channel];
1166 mergeExtraHitToChannelHit(hit, channelHit);
1172 otherModuleHits.push_back(std::make_unique<TileHit>(*hit));
1175 hits[channel] = hit;
1180 for (std::unique_ptr<TileHit>& hit : otherModuleHits) {
1187 coll->push_back(std::move(hit));
1201 channelHit->
add(extraHit, 0.1);
1211 allHits.reserve(nHits);
1217 allHits.emplace_back(std::make_unique<TileHit>(hit_id, 0., 0.));
1218 allHits.back()->reserve(71);
1227 allHits[mbtsIndex] = std::make_unique<TileHit>(hit_id, 0., 0.);
1228 allHits[mbtsIndex]->reserve(71);
1237 allHits[e4prIndex] = std::make_unique<TileHit>(hit_id, 0., 0.);
1238 allHits[e4prIndex]->reserve(71);
def retrieve(aClass, aKey=None)
JetConstituentVector::iterator iterator
void scale(float coeff)
Scale energy of all sub-hits in a hit
Identifier identify(void) const
Return logical ID of the pmt.
float getSamplingFraction(unsigned int drawerIdx, unsigned int channel) const
Return Tile Calorimeter sampling fraction.
void setSeed(const std::string &algName, const EventContext &ctx)
Set the random seed using a string (e.g.
Identifier channel_id(int type, int module, int channel) const
identifer for one channel of a Tile testbeam detector
std::string to_string(const Identifier &id, int level=0) const
extract all fields from TileTB identifier Identifier get_all_fields ( const Identifier & id,...
Scalar phi() const
phi method
bool is_tiletb(const Identifier &id) const
Test ID if it is TileTBID.
const_pointer_type cptr()
Dereference the pointer.
IdContext channel_context() const
idContext for channels
@ VIEW_ELEMENTS
this data object is a view, it does not own its elmts
AtlasHitsVector< TileHit >::const_iterator TileHitVecConstIterator
Scalar eta() const
pseudorapidity method
float round(const float toRound, const unsigned int decimals)
int side(const Identifier &id) const
int sample(const Identifier &id) const
const T * get(size_type n) const
Access an element, as an rvalue.
int frag(const HWIdentifier &id) const
extract frag field from HW identifier
int tower(const Identifier &id) const
#define ATH_MSG_VERBOSE(x)
const std::string & key() const
Return the StoreGate ID for the referenced object.
bool empty() const
Test if the key is blank.
int type(const Identifier &id) const
extract type field from TileTB identifier
size_type pmt_hash_max() const
const_iterator begin() const
int size(void) const
Return length of energy/time vectors
int module(const Identifier &id) const
extract module field from TileTB identifier
int channel(const HWIdentifier &id) const
extract channel field from HW identifier
int ros(const HWIdentifier &id) const
extract ros field from HW identifier
std::list< value_t > type
type of the collection of timed data object
OwnershipPolicy
describes the possible element ownership policies (see e.g. DataVector)
float getNumberOfPhotoElectrons(unsigned int drawerIdx, unsigned int channel) const
Return number of photoelectrons per 1 GeV in Tile Calorimeter scintilator.
Handle class for recording to StoreGate.
bool isRun2PlusCabling() const
virtual int get_id(const IdentifierHash &hash_id, Identifier &id, const IdContext *context=0) const
create compact id from hash id (return == 0 for OK)
std::vector< std::unique_ptr< TileHitCollection > >::iterator iterator
int module(const Identifier &id) const
Condition object to keep and provide Tile Calorimeter sampling fraction and number of photoelectrons.
(Non-const) Iterator class for DataVector/DataList.
@ OWN_ELEMENTS
this data object owns its elements
::StatusCode StatusCode
StatusCode definition for legacy code.
int add(float energy, float time)
Add sub-hit to a given hit.
#define CHECK(...)
Evaluate an expression and check for errors.
StatusCode initialize(bool used=true)
If this object is used as a property, then this should be called during the initialize phase.
bool isRun2Cabling() const
void initialize(const TileHWID *tileHWID, TYPE type=Default)
int channel(const Identifier &id) const
extract channel field from TileTB identifier
virtual bool isValid() override final
Can the handle be successfully dereferenced?
virtual int get_hash(const Identifier &id, IdentifierHash &hash_id, const IdContext *context=0) const
create hash id from compact id (return == 0 for OK)
void push_back(element_t *rc)
size_type drawer_hash_max() const
drawer hash table max size
HWIdentifier drawer_id(int frag) const
ROS HWIdentifer.
bool is_MBTS_merged_run2plus(int module) const
A wrapper class for event-slot-local random engines.
Helpers for checking error return status codes and reporting errors.
StatusCode initialize(bool used=true)
CLHEP::HepRandomEngine * getEngine(const EventContext &ctx) const
Retrieve the random engine corresponding to the provided EventContext.
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
This is a minimal version of a TileHitContainer in which the saved collections remain non-const.
int drawer(const HWIdentifier &id) const
extract drawer field from HW identifier
StatusCode record(std::unique_ptr< T > data)
Record a const object to the store.
def time(flags, cells_name, *args, **kw)
HWIdentifier pmt_HWID(void) const
Return pmt hardware ID (== channel ID)
std::string to_string(const Identifier &id, int level=0) const
#define ATH_MSG_WARNING(x)
float energy(int ind=0) const
Return energy of ind-th sub-hit
Identifier pmt_id(const Identifier &any_id) const
const_iterator end() const
T * get(TKey *tobj)
get a TObject* from a TKey* (why can't a TObject be a TKey?)
void setTime(float t, int ind=0)
Set time of ind-th sub-hit in a hit
int phi(const Identifier &id) const
extract phi field from MBTS identifier
static unsigned int getDrawerIdx(unsigned int ros, unsigned int drawer)
Returns a drawer hash.
This is a "hash" representation of an Identifier. This encodes a 32 bit index which can be used to lo...
iterator erase(iterator position)
Remove element at a given position.
Identifier cell_id(const Identifier &any_id) const
float time(int ind=0) const
Return time of ind-th sub-hit
This class saves the "context" of an expanded identifier (ExpandedIdentifier) for compact or hash ver...
IdContext pmt_context() const
id for PMTs
int E1_merged_with_run2plus(int ros, int module) const
Identifier pmt_ID(void) const
Return logical ID of the pmt.
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
int side(const Identifier &id) const
define synonyms for minimum bias scintillators
HWIdentifier s2h_channel_id(const Identifier &id) const