ATLAS Offline Software
ALFA_PileUpTool.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include "ALFA_Digitization/ALFA_PileUpTool.h"
6 
7 #include "ALFA_SimEv/ALFA_HitCollection.h"
8 #include "ALFA_SimEv/ALFA_Hit.h"
9 #include "ALFA_SimEv/ALFA_ODHitCollection.h"
10 
11 #include "ALFA_RawEv/ALFA_DigitCollection.h"
12 #include "ALFA_SimEv/ALFA_ODHit.h"
13 #include "ALFA_RawEv/ALFA_ODDigitCollection.h"
14 
15 #include "AthenaKernel/RNGWrapper.h"
16 #include "CLHEP/Random/RandomEngine.h"
17 #include "CLHEP/Random/RandFlat.h"
18 #include "CLHEP/Random/RandGaussZiggurat.h"
19 #include "CLHEP/Random/RandPoissonQ.h"
20 
21 
22 
23 #include "Identifier/Identifier.h"
24 #include "AthenaBaseComps/AthMsgStreamMacros.h"
25 
26 // from now old part
27 
28 
29 //for truth particles
30 #include "AtlasHepMC/GenParticle.h"
31 #include "GeneratorObjects/McEventCollection.h"
32 #include "AtlasHepMC/GenEvent.h"
33 #include "AtlasHepMC/GenVertex.h"
34 #include "CLHEP/Vector/LorentzVector.h"
35 #include "AtlasHepMC/SimpleVector.h"
36 
37 #include "TTree.h"
38 
39 #include "GaudiKernel/ITHistSvc.h"
40 
41 #include "PathResolver/PathResolver.h"
42 #include <fstream>
43 #include <algorithm>
44 #include <map>
45 #include <cmath>
46 #include <functional>
47 
48 // *********************************************************************************
49 
50 ALFA_PileUpTool::ALFA_PileUpTool(const std::string& type,
51  const std::string& name,
52  const IInterface* parent) :
53  PileUpToolBase(type, name, parent),
54  m_SimHitCollectionName ("ALFA_HitCollection"),
55  m_SimODHitCollectionName ("ALFA_ODHitCollection"),
56  m_key_DigitCollection ("ALFA_DigitCollection"),
57  m_key_ODDigitCollection ("ALFA_ODDigitCollection"),
58  m_digitCollection (nullptr),
59  m_ODdigitCollection (nullptr),
60  m_sigma0 (0.1),
61  m_sigma1 (1.0),
62  m_meanE_dep (0.0863), //MeV
63  m_meanN_photo (4.11),
64  m_mean (0.0),
65  m_stdDev (1.0),
66  m_AmplitudeCut (0.24),
67  m_mergedALFA_HitList (nullptr),
68  m_mergedALFA_ODHitList (nullptr)
69 {
70  // Properties
71  declareProperty("SimHitCollection", m_SimHitCollectionName, "Name of the input ALFA Sim Hit Collection of simulated hits");
72  declareProperty("SimODHitCollection", m_SimODHitCollectionName, "Name of the input ALFA Sim OD Hit Collection of simulated hits");
73 
74  declareProperty("ALFA_DigitCollection", m_key_DigitCollection, "Name of the Collection to hold the output from the ALFA main detector digitization");
75  declareProperty("ALFA_ODDigitCollection", m_key_ODDigitCollection, "Name of the Collection to hold the output from the ALFA OD digitization");
76 
77 
78 // declareProperty("fillRootTree", m_fillRootTree);
79 
80 
81  declareProperty("sigma0", m_sigma0);
82  declareProperty("sigma1", m_sigma1);
83  declareProperty("meanE_dep", m_meanE_dep);
84  declareProperty("meanN_photo", m_meanN_photo);
85  declareProperty("mean", m_mean);
86  declareProperty("stdDev", m_stdDev);
87  declareProperty("AmplitudeCut", m_AmplitudeCut);
88 // declareProperty("fillRootTree", m_fillRootTree);
89 }
90 
91 
92 
93 
94 StatusCode ALFA_PileUpTool::initialize(){
95 
96  ATH_MSG_DEBUG ("ALFA_digiti initialize() called");
97 
98  ATH_MSG_DEBUG (" Pedestal: " << m_AmplitudeCut);
99 
100  ATH_CHECK (m_randomSvc.retrieve());
101  ATH_MSG_DEBUG ( "Retrieved RandomNumber Service" );
102 
103  ATH_CHECK (m_mergeSvc.retrieve());
104  ATH_MSG_DEBUG("Retrieved PileUpMergeSvc");
105 
106  m_mergedALFA_HitList = new ALFA_HitCollection();
107  m_mergedALFA_ODHitList = new ALFA_ODHitCollection();
108 
109  return StatusCode::SUCCESS;
110 }
111 
112 
113 
114 
115 // *********************************************************************************
116 
117 StatusCode ALFA_PileUpTool::processAllSubEvents(const EventContext& ctx) {
118 
119  ATH_MSG_DEBUG ("ALFA_PileUpTool::processAllSubEvents()");
120 
121  using TimedALFAHitCollList = PileUpMergeSvc::TimedList<ALFA_HitCollection>::type;
122  using TimedALFAODHitCollList = PileUpMergeSvc::TimedList<ALFA_ODHitCollection>::type;
123 
124  TimedALFAHitCollList tHitCollList;
125  TimedALFAODHitCollList tODHitCollList;
126 
127 
128  if (m_mergeSvc->retrieveSubEvtsData(m_SimHitCollectionName, tHitCollList).isFailure()) {
129 
130  ATH_MSG_FATAL ( "Could not fill TimedALFAHitCollList" ); return StatusCode::FAILURE;
131  }
132  else { ATH_MSG_DEBUG ( "Retrieved TimedALFAHitCollList" ); }
133 
134  ATH_MSG_DEBUG ( "PileUp: Merge " << tHitCollList.size() << " ALFA_HitCollection with key" << m_SimHitCollectionName );
135 
136 
137  if (m_mergeSvc->retrieveSubEvtsData(m_SimODHitCollectionName, tODHitCollList).isFailure()) {
138 
139  ATH_MSG_FATAL ( "Could not fill TimedALFAODHitCollList" ); return StatusCode::FAILURE;
140  }
141  else { ATH_MSG_DEBUG ( "Retrieved TimedALFAODHitCollList" ); }
142 
143  ATH_MSG_DEBUG ( "PileUp: Merge " << tODHitCollList.size() << " ALFA_ODHitCollection with key " << m_SimODHitCollectionName );
144 
145 
146  TimedHitCollection<ALFA_Hit> tALFAhit;
147 
148  TimedALFAHitCollList::iterator iHitColl (tHitCollList.begin());
149  TimedALFAHitCollList::iterator eHitColl (tHitCollList.end());
150 
151  while (iHitColl != eHitColl) {
152 
153  const ALFA_HitCollection* tmpColl(iHitColl->second);
154 
155  tALFAhit.insert(iHitColl->first, tmpColl);
156 
157  ATH_MSG_DEBUG ( " ALFA_HitCollection found with " << tmpColl->size() << " hits " << iHitColl->first );
158 
159  ++iHitColl;
160  }
161 
162  TimedHitCollection<ALFA_ODHit> tALFAODhit;
163 
164  TimedALFAODHitCollList::iterator iODHitColl (tODHitCollList.begin());
165  TimedALFAODHitCollList::iterator eODHitColl (tODHitCollList.end());
166 
167  while (iODHitColl != eODHitColl) {
168 
169  const ALFA_ODHitCollection* tmpColl(iODHitColl->second);
170 
171  tALFAODhit.insert(iODHitColl->first, tmpColl);
172 
173  ATH_MSG_DEBUG ( " ALFA_ODHitCollection found with " << tmpColl->size() << " hits " << iODHitColl->first );
174 
175  ++iODHitColl;
176  }
177 
178  // preparing containers
179 
180  if (recordCollection(this->evtStore(), m_key_DigitCollection).isFailure()) {
181 
182  ATH_MSG_FATAL ( " ALFA_PileUpTool::processAllSubEvents(): Could not record the empty ALFA digit container in StoreGate " ); return StatusCode::FAILURE;
183  }
184  else { ATH_MSG_DEBUG ( " ALFA_PileUpTool::processAllSubEvents(): ALFA Digit container is recorded in StoreGate " ); }
185 
186 
187  if (recordODCollection(this->evtStore(), m_key_ODDigitCollection).isFailure()) {
188 
189  ATH_MSG_FATAL ( " ALFA_PileUpTool::processAllSubEvents(): Could not record the empty ALFA OD digit container in StoreGate " ); return StatusCode::FAILURE;
190  }
191  else { ATH_MSG_DEBUG ( " ALFA_PileUpTool::processAllSubEvents(): ALFA OD Digit container is recorded in StoreGate " ); }
192 
193 
194  // filling containers
195 
196  ALFA_MD_info(tALFAhit);
197 
198  // Prepare RNG Service
199  ATHRNG::RNGWrapper* rngWrapper = m_randomSvc->getEngine(this, m_randomStreamName);
200  rngWrapper->setSeed( m_randomStreamName, ctx );
201  CLHEP::HepRandomEngine* rngEngine = rngWrapper->getEngine(ctx);
202  StatusCode sc = fill_MD_DigitCollection(rngEngine);
203  if (sc.isFailure()) {
204  ATH_MSG_WARNING ("ALFA_PileUpTool::fill_MD_DigitCollection() failed");
205  return StatusCode::SUCCESS;
206  }
207 
208 
209  ALFA_OD_info(tALFAODhit);
210 
211  sc = fill_OD_DigitCollection(rngEngine);
212  if (sc.isFailure()) {
213  ATH_MSG_WARNING ("ALFA_PileUpTool::fill_OD_DigitCollection() failed");
214  return StatusCode::SUCCESS;
215  }
216 
217  return StatusCode::SUCCESS;
218 }
219 
220 
221 
222 
223 StatusCode ALFA_PileUpTool::prepareEvent(const EventContext& /*ctx*/, const unsigned int nInputEvents){
224 
225  ATH_MSG_DEBUG ( "ALFA_PileUpTool::prepareEvent() called for " << nInputEvents << " input events" );
226 
227  StatusCode sc = recordCollection(this->evtStore(), m_key_DigitCollection);
228 
229  if (sc.isFailure()) { ATH_MSG_FATAL ( " ALFA_PileUpTool::prepareEvent(): Could not record the empty digit container in StoreGate " ); return sc; }
230  else { ATH_MSG_DEBUG ( " ALFA_PileUpTool::prepareEvent(): Digit container is recorded in StoreGate " ); }
231 
232  sc = recordODCollection(this->evtStore(), m_key_ODDigitCollection);
233 
234  if (sc.isFailure()) { ATH_MSG_FATAL ( " ALFA_PileUpTool::prepareEvent(): Could not record the empty digit OD container in StoreGate " ); return sc; }
235  else { ATH_MSG_DEBUG ( " ALFA_PileUpTool::prepareEvent(): Digit OD container is recorded in StoreGate " ); }
236 
237  m_mergedALFA_HitList->clear();
238  m_mergedALFA_ODHitList->clear();
239 
240  return StatusCode::SUCCESS;
241 }
242 
243 
244 StatusCode ALFA_PileUpTool::processBunchXing(int bunchXing,
245  SubEventIterator bSubEvents,
246  SubEventIterator eSubEvents) {
247  ATH_MSG_DEBUG ( "ALFA_PileUpTool::processBunchXing() " << bunchXing );
248  SubEventIterator iEvt = bSubEvents;
249  for (; iEvt!=eSubEvents; ++iEvt) {
250  StoreGateSvc& seStore = *iEvt->ptr()->evtStore();
251  //PileUpTimeEventIndex thisEventIndex = PileUpTimeEventIndex(static_cast<int>(iEvt->time()),iEvt->index()); not used
252  ATH_MSG_VERBOSE("SubEvt StoreGate " << seStore.name() << " :"
253  << " bunch crossing : " << bunchXing
254  << " time offset : " << iEvt->time()
255  << " event number : " << iEvt->ptr()->eventNumber()
256  << " run number : " << iEvt->ptr()->runNumber()
257  );
258 
259 
260  const ALFA_HitCollection* tmpHitColl = nullptr;
261 
262  if (!seStore.retrieve(tmpHitColl, m_SimHitCollectionName).isSuccess()) {
263 
264  ATH_MSG_ERROR ("SubEvent ALFA_HitCollection not found in StoreGate " << seStore.name() );
265 
266  return StatusCode::FAILURE;
267  }
268 
269  ATH_MSG_DEBUG ("SubEvent, ALFA_HitCollection found with " << tmpHitColl->size() << " hits");
270 
271  ALFA_HitCollection::const_iterator iHitColl = tmpHitColl->begin();
272  ALFA_HitCollection::const_iterator eHitColl = tmpHitColl->end();
273 
274  for (; iHitColl!=eHitColl; ++iHitColl) m_mergedALFA_HitList->push_back((*iHitColl));
275 
276  const ALFA_ODHitCollection* tmpODHitColl = nullptr;
277 
278  if (!seStore.retrieve(tmpODHitColl, m_SimODHitCollectionName).isSuccess()) {
279 
280  ATH_MSG_ERROR ( "SubEvent ALFA_ODHitCollection not found in StoreGate " << seStore.name() );
281 
282  return StatusCode::FAILURE;
283  }
284 
285  ATH_MSG_DEBUG ( "ALFA_ODHitCollection found with " << tmpODHitColl->size() << " hits" );
286 
287  ALFA_ODHitCollection::const_iterator iODHitColl = tmpODHitColl->begin();
288  ALFA_ODHitCollection::const_iterator eODHitColl = tmpODHitColl->end();
289 
290  for (; iODHitColl!=eODHitColl; ++iODHitColl) m_mergedALFA_ODHitList->push_back((*iODHitColl));
291  }
292 
293  return StatusCode::SUCCESS;
294 }
295 
296 
297 
298 
299 StatusCode ALFA_PileUpTool::mergeEvent(const EventContext& ctx){
300 
301 
302  ALFA_MD_info(m_mergedALFA_HitList);
303 
304  // Prepare RNG Service
305  ATHRNG::RNGWrapper* rngWrapper = m_randomSvc->getEngine(this, m_randomStreamName);
306  rngWrapper->setSeed( m_randomStreamName, ctx );
307  CLHEP::HepRandomEngine* rngEngine = rngWrapper->getEngine(ctx);
308  StatusCode sc = fill_MD_DigitCollection(rngEngine);
309  if (sc.isFailure()) {
310  ATH_MSG_WARNING ("ALFA_PileUpTool::fill_MD_DigitCollection() failed");
311  return StatusCode::SUCCESS;
312  }
313 
314 
315  ALFA_OD_info(m_mergedALFA_ODHitList);
316 
317  sc = fill_OD_DigitCollection(rngEngine);
318  if (sc.isFailure()) {
319  ATH_MSG_WARNING ("ALFA_PileUpTool::fill_OD_DigitCollection() failed");
320  return StatusCode::SUCCESS;
321  }
322 
323  return StatusCode::SUCCESS;
324 }
325 
326 StatusCode ALFA_PileUpTool::finalize() { return StatusCode::SUCCESS; }
327 
328 
329 
330 StatusCode ALFA_PileUpTool::recordCollection(ServiceHandle<StoreGateSvc>& evtStore, const std::string& key_digitCnt)
331 {
332  ATH_MSG_DEBUG ("ALFA_Digitization::recordCollection");
333 
334  m_digitCollection = new ALFA_DigitCollection();
335 
336  StatusCode sc = evtStore->record(m_digitCollection, key_digitCnt);
337 
338  return sc;
339 }
340 
341 
342 StatusCode ALFA_PileUpTool::recordODCollection(ServiceHandle<StoreGateSvc>& evtStore, const std::string& key_ODdigitCnt) {
343 
344  ATH_MSG_DEBUG ("ALFA_Digitization::recordODCollection");
345 
346  m_ODdigitCollection = new ALFA_ODDigitCollection();
347 
348  StatusCode sc = evtStore->record(m_ODdigitCollection, key_ODdigitCnt);
349 
350  return sc;
351 }
352 
353 
354 
355 
356 void ALFA_PileUpTool::ALFA_MD_info(TimedHitCollection<ALFA_Hit>& tHitALFA)
357 {
358 
359 // cleaning
360 
361  for ( int l = 0; l < 8; l++ )
362  {
363  for ( int i = 0; i < 20; i++ )
364  {
365  for ( int j = 0; j < 64; j++ )
366  {
367 
368  m_E_fib[l][i][j] = 0.;
369  m_fibres[l][i][j] = 0;
370 
371  }
372  }
373  }
374 
375 // cleaning end
376 
377  int fiber, plate, sign, station;
378 
379  TimedHitCollection<ALFA_Hit> thpc = tHitALFA;
380  TimedHitCollection<ALFA_Hit>::const_iterator i, e, it;
381 
382  while (thpc.nextDetectorElement(i, e)) for (it = i; it != e; ++it) {
383 
384  station = (*it)->GetStationNumber();
385  plate = (*it)->GetPlateNumber();
386  fiber = (*it)->GetFiberNumber();
387  sign = (*it)->GetSignFiber();
388 
389  ATH_MSG_DEBUG ("station=" << station << ", plate= "<< plate << ", fiber=" << fiber << ", sign=" << sign << ", dep energy=" << ((*it)->GetEnergyDeposit()));
390 
391  m_E_fib[station-1][2*(plate-1)+(1-sign)/2][fiber-1] += ((*it)->GetEnergyDeposit());
392 
393  }
394 
395 }
396 
397 
398 
399 
400 void ALFA_PileUpTool::ALFA_MD_info(const ALFA_HitCollection* ALFA_HitCollection)
401 {
402 
403 // cleaning
404 
405  for ( int l = 0; l < 8; l++ )
406  {
407  for ( int i = 0; i < 20; i++ )
408  {
409  for ( int j = 0; j < 64; j++ )
410  {
411 
412  m_E_fib[l][i][j] = 0.;
413  m_fibres[l][i][j] = 0;
414 
415 
416  }
417  }
418  }
419 
420 // cleaning end
421 
422  int fiber, plate, sign, station;
423 
424  ALFA_HitConstIter it = ALFA_HitCollection->begin();
425  ALFA_HitConstIter itend = ALFA_HitCollection->end();
426 
427  for (; it != itend; ++it) {
428 
429  station = it->GetStationNumber();
430  plate = it->GetPlateNumber();
431  fiber = it->GetFiberNumber();
432  sign = it->GetSignFiber();
433 
434  ATH_MSG_DEBUG ("station=" << station << ", plate= "<< plate << ", fiber=" << fiber << ", sign=" << sign << ", dep energy=" << it->GetEnergyDeposit());
435 
436  m_E_fib[station-1][2*(plate-1)+(1-sign)/2][fiber-1] += it->GetEnergyDeposit();
437 
438 
439  }
440 }
441 
442 
443 
444 
445 void ALFA_PileUpTool::ALFA_OD_info(TimedHitCollection<ALFA_ODHit>& tODHitALFA)
446 {
447 
448 // cleaning
449 
450  for ( int l = 0; l < 8; l++ ){
451  for ( int i = 0; i < 2; i++ ){
452  for ( int j = 0; j < 3; j++ ){
453  for ( int k = 0; k < 30; k++ ){
454 
455  m_E_ODfib[l][i][j][k] = 0.;
456 
457  }
458  }
459  }
460  }
461 
462 // cleaning end
463 
464  int fiber, plate, sign, side, station;
465 
466 
467  TimedHitCollection<ALFA_ODHit> thpc = tODHitALFA;
468  TimedHitCollection<ALFA_ODHit>::const_iterator i, e, it;
469 
470  while (thpc.nextDetectorElement(i, e)) for (it = i; it != e; ++it) {
471 
472  station = (*it)->GetStationNumber();
473  side = (*it)->GetODSide();
474  plate = (*it)->GetPlateNumber();
475  fiber = (*it)->GetFiberNumber();
476  sign = (*it)->GetSignFiber();
477 
478 
479  ATH_MSG_DEBUG ("station=" << station << ", side=" << side << ", plate= "<< plate << ", fiber=" << fiber << ", sign=" << sign << ", dep energy=" << ((*it)->GetEnergyDeposit()));
480 
481  if (sign==0) m_E_ODfib[station-1][side-1][plate-1][fiber+15] += ((*it)->GetEnergyDeposit());
482  else m_E_ODfib[station-1][side-1][plate-1][fiber] += ((*it)->GetEnergyDeposit());
483 
484  }
485 
486 }
487 
488 
489 
490 
491 void ALFA_PileUpTool::ALFA_OD_info(const ALFA_ODHitCollection* ALFA_ODHitCollection)
492 {
493 
494  // cleaning
495 
496  for ( int l = 0; l < 8; l++ ){
497  for ( int i = 0; i < 2; i++ ){
498  for ( int j = 0; j < 3; j++ ){
499  for ( int k = 0; k < 30; k++ ){
500 
501  m_E_ODfib[l][i][j][k] = 0.;
502 
503  }
504  }
505  }
506  }
507 
508 // cleaning end
509 
510  int fiber, plate, sign, side, station;
511 
512 
513  ALFA_ODHitConstIter it = ALFA_ODHitCollection->begin();
514  ALFA_ODHitConstIter itend = ALFA_ODHitCollection->end();
515 
516  for (; it != itend; ++it) {
517 
518 
519  station = it->GetStationNumber();
520  side = it->GetODSide();
521  plate = it->GetPlateNumber();
522  fiber = it->GetFiberNumber();
523  sign = it->GetSignFiber();
524 
525 
526  ATH_MSG_DEBUG ("station=" << station << ", side=" << side << ", plate= "<< plate << ", fiber=" << fiber << ", sign=" << sign << ", dep energy=" << it->GetEnergyDeposit());
527 
528  if (sign==0) m_E_ODfib[station-1][side-1][plate-1][fiber+15] += it->GetEnergyDeposit();
529  else m_E_ODfib[station-1][side-1][plate-1][fiber] += it->GetEnergyDeposit();
530 
531  }
532 }
533 
534 
535 
536 
537 StatusCode ALFA_PileUpTool::fill_MD_DigitCollection(CLHEP::HepRandomEngine* rndEngine) {
538 
539  ATH_MSG_DEBUG(" ALFA_PileUpTool::fill_MD_DigitCollection()");
540 
541  double N_photo = 0.;
542  //double N_CTphoto = 0.;
543  //double sigma = 0.;
544  double amplitude = 0.;
545  //double amplitude_CT = 0.;
546 
547 
548  for ( int l = 0; l < 8; l++ )
549  {
550  for ( int i = 0; i < 20; i++ )
551  {
552  for ( int j = 0; j < 64; j++ )
553  {
554 
555  N_photo = CLHEP::RandPoisson::shoot(rndEngine,m_E_fib[l][i][j]*m_meanN_photo/m_meanE_dep);
556 
557  //sigma = m_sigma1 * sqrt(N_photo);
558  //double noise_1 = sigma * CLHEP::RandGaussZiggurat::shoot (rndEngine, m_mean, m_stdDev);
559  //double noise_2 = m_sigma0 * CLHEP::RandGaussZiggurat::shoot (rndEngine, m_mean, m_stdDev);
560 
561  //amplitude = N_photo
562  // + noise_1
563  // + noise_2;
564 
565  amplitude = CLHEP::RandGaussZiggurat::shoot(rndEngine, N_photo, sqrt(pow(m_sigma0,2)+N_photo*pow(m_sigma1,2)));
566 
567  //N_photo_CT = CLHEP::RandPoisson::shoot(rndEngine,0.08*m_E_fib[l][i][j]*m_meanN_photo/m_meanE_dep);
568  //amplitude_CT = CLHEP::RandGaussQ::shoot (rndEngine, N_photo_CT, sqrt(pow(m_sigma0,2)+N_photo_CT*pow(m_sigma1,2)));
569 
570 
571  if (amplitude >= m_AmplitudeCut )
572  {
573  ATH_MSG_DEBUG(" ALFA_Digitization::fillDigitCollection, amplitude " << amplitude);
574  ATH_MSG_DEBUG(" station = " << l << ", plate= " << i << ", fiber=" << j );
575  m_digitCollection->push_back(new ALFA_Digit(l,i,j));
576  m_fibres[l][i][j] = m_fibres[l][i][j] + 1;
577  }
578 
579  }
580  }
581  }
582 
583  double rand_fib;
584 
585  for ( int l = 0; l < 8; l++ )
586  {
587  for ( int i = 0; i < 20; i++ )
588  {
589  for ( int j = 0; j < 64; j++ )
590  {
591  if ( m_fibres[l][i][j] > 0 )
592  {
593  for (int f = j+1; f < 64; f++)
594  {
595  rand_fib = CLHEP::RandFlat::shoot(rndEngine,0.,1.);
596 
597  if ( m_fibres[l][i][f] == 0)
598  {
599  if (rand_fib <= m_cross_talk[l][63+f-j])
600  {
601  m_digitCollection->push_back(new ALFA_Digit(l,i,f));
602  }
603  }
604  }
605 
606  for (int f = j-1; f > -1; f--)
607  {
608  rand_fib = CLHEP::RandFlat::shoot(rndEngine,0.,1.);
609 
610  if ( m_fibres[l][i][f] == 0)
611  {
612  if (rand_fib <= m_cross_talk[l][63-(j-f)])
613  {
614  m_digitCollection->push_back(new ALFA_Digit(l,i,f));
615  }
616  }
617  }
618  }
619  }
620  }
621  }
622  return StatusCode::SUCCESS;
623 
624 }
625 
626 
627 
628 
629 StatusCode ALFA_PileUpTool::fill_OD_DigitCollection(CLHEP::HepRandomEngine* rndEngine) {
630 
631  ATH_MSG_DEBUG(" ALFA_PileUpTool::fill_OD_DigitCollection ");
632 
633  double N_photo = 0.;
634  double sigma = 0.;
635  double amplitude = 0.;
636 
637  for ( int l = 0; l < 8; l++ )
638  {
639  for ( int i = 0; i < 2; i++ )
640  {
641  for ( int j = 0; j < 3; j++ )
642  {
643  for (int k = 0; k < 30; k++)
644  {
645 
646  N_photo = CLHEP::RandPoisson::shoot(rndEngine,m_E_ODfib[l][i][j][k]*m_meanN_photo/m_meanE_dep);
647  sigma = m_sigma1 * sqrt(N_photo);
648 
649  double noise_1 = sigma * CLHEP::RandGaussZiggurat::shoot (rndEngine, m_mean, m_stdDev);
650  double noise_2 = m_sigma0 * CLHEP::RandGaussZiggurat::shoot (rndEngine, m_mean, m_stdDev);
651 
652 
653  amplitude = N_photo
654  + noise_1
655  + noise_2;
656 
657 
658  if (amplitude >= m_AmplitudeCut )
659  {
660  ATH_MSG_DEBUG(" ALFA_Digitization::fill_OD_DigitCollection, amplitude " << amplitude);
661  ATH_MSG_DEBUG(" station = " << l << ", side = " << i << ", plate = " << j << ", fiber = " << k);
662 
663  m_ODdigitCollection->push_back(new ALFA_ODDigit(l,i,j,k));
664  }
665 
666 
667  }
668  }
669  }
670  }
671 
672  return StatusCode::SUCCESS;
673 }
674 
675 
676 StatusCode
677 ALFA_PileUpTool::XTalk() {
678  std::ifstream fXTalk;
679  std::stringstream ss;
680  std::string fname;
681  for (unsigned int j=0; j<8; j++){
682 
683  ss.str("");
684 
685  fname = "Xtalk_station";
686  ss << j+1;
687  fname += ss.str();
688  fname += ".txt";
689 
690  ATH_MSG_DEBUG("file name " << fname);
691 
692  std::string filePath = PathResolver::find_file(fname,"DATAPATH", PathResolver::RecursiveSearch);
693 
694  if(filePath.length() == 0)
695  {
696  ATH_MSG_FATAL(" XTalk file " << fname<< " not found in Datapath");
697  throw std::runtime_error("FATAL: mapping MD maroc-mapmt not found in Datapath.");
698  }
699 
700  else
701  {
702  ATH_MSG_DEBUG("the XTALK file \"" << fname << "\" found in Datapath");
703  ATH_MSG_DEBUG("filePath = " << filePath);
704  }
705 
706  fXTalk.open(filePath.c_str());
707 
708  if (fXTalk.is_open())
709  {
710  for (unsigned int i=0;i<127;i++)
711  {
712  fXTalk >> m_cross_talk[j][i];
713 
714  ATH_MSG_DEBUG("cross_talk[" << j << "][" << i << "] = " << m_cross_talk[j][i]);
715 
716  }
717  fXTalk.close();
718  }
719  else
720  {
721  ATH_MSG_ERROR("the file " << fname<< " was not open");
722  return StatusCode::FAILURE;
723  }
724  }
725  return StatusCode::SUCCESS;
726  }
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
ALFA_PileUpTool::m_fibres
int m_fibres[8][20][64]
Definition: ALFA_PileUpTool.h:111
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
ATHRNG::RNGWrapper::setSeed
void setSeed(const std::string &algName, const EventContext &ctx)
Set the random seed using a string (e.g.
Definition: RNGWrapper.h:169
ALFA_HitCollection
AtlasHitsVector< ALFA_Hit > ALFA_HitCollection
Definition: ALFA_HitCollection.h:13
ALFA_PileUpTool::ALFA_MD_info
void ALFA_MD_info(const ALFA_HitCollection *)
Definition: ALFA_PileUpTool.cxx:400
ALFA_ODDigitCollection.h
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition: AthMsgStreamMacros.h:34
pdg_comparison.sigma
sigma
Definition: pdg_comparison.py:324
GenEvent.h
PathResolver::RecursiveSearch
@ RecursiveSearch
Definition: PathResolver.h:28
ALFA_ODHit.h
ALFA_PileUpTool::processBunchXing
virtual StatusCode processBunchXing(int bunchXing, SubEventIterator bSubEvents, SubEventIterator eSubEvents) override final
called for each active bunch-crossing to process current SubEvents bunchXing is in ns
Definition: ALFA_PileUpTool.cxx:244
PowhegControl_ttHplus_NLO.ss
ss
Definition: PowhegControl_ttHplus_NLO.py:83
ALFA_PileUpTool::m_meanN_photo
double m_meanN_photo
Definition: ALFA_PileUpTool.h:102
PathResolver::find_file
static std::string find_file(const std::string &logical_file_name, const std::string &search_path, SearchType search_type=LocalSearch)
Definition: PathResolver.cxx:251
AthMsgStreamMacros.h
ALFA_PileUpTool::m_key_ODDigitCollection
std::string m_key_ODDigitCollection
Definition: ALFA_PileUpTool.h:93
GenVertex.h
AtlasHitsVector
Definition: AtlasHitsVector.h:33
ALFA_ODHitCollection.h
skel.it
it
Definition: skel.GENtoEVGEN.py:396
ALFA_PileUpTool::ALFA_OD_info
void ALFA_OD_info(const ALFA_ODHitCollection *)
Definition: ALFA_PileUpTool.cxx:491
UploadAMITag.l
list l
Definition: UploadAMITag.larcaf.py:158
ALFA_HitConstIter
AtlasHitsVector< ALFA_Hit >::const_iterator ALFA_HitConstIter
Definition: ALFA_HitCollection.h:15
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
AtlasHitsVector::begin
const_iterator begin() const
Definition: AtlasHitsVector.h:131
ALFA_PileUpTool::m_AmplitudeCut
double m_AmplitudeCut
Definition: ALFA_PileUpTool.h:105
TimedHitCollection::nextDetectorElement
bool nextDetectorElement(const_iterator &b, const_iterator &e)
sets an iterator range with the hits of current detector element returns a bool when done
GenParticle.h
AtlasHitsVector::const_iterator
CONT::const_iterator const_iterator
Definition: AtlasHitsVector.h:43
ALFA_PileUpTool.h
AthenaPoolTestRead.sc
sc
Definition: AthenaPoolTestRead.py:27
AtlasHitsVector::clear
void clear()
Definition: AtlasHitsVector.h:153
ALFA_PileUpTool::initialize
virtual StatusCode initialize() override final
Definition: ALFA_PileUpTool.cxx:94
PileUpMergeSvc::TimedList::type
std::list< value_t > type
type of the collection of timed data object
Definition: PileUpMergeSvc.h:75
StoreGateSvc::retrieve
StatusCode retrieve(const T *&ptr) const
Retrieve the default object into a const T*.
TRT::Hit::side
@ side
Definition: HitInfo.h:83
ALFA_PileUpTool::m_mergedALFA_HitList
ALFA_HitCollection * m_mergedALFA_HitList
Definition: ALFA_PileUpTool.h:107
ALFA_PileUpTool::m_randomStreamName
Gaudi::Property< std::string > m_randomStreamName
Definition: ALFA_PileUpTool.h:84
StoreGateSvc
The Athena Transient Store API.
Definition: StoreGateSvc.h:125
ALFA_PileUpTool::m_ODdigitCollection
ALFA_ODDigitCollection * m_ODdigitCollection
Definition: ALFA_PileUpTool.h:96
ALFA_Digit
Definition: ALFA_Digit.h:8
SimpleVector.h
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
ALFA_PileUpTool::recordCollection
StatusCode recordCollection(ServiceHandle< StoreGateSvc > &evtStore, const std::string &key_digitCnt)
Definition: ALFA_PileUpTool.cxx:330
TimedHitCollection::insert
void insert(const PileUpTimeEventIndex &timeEventIndex, const AtlasHitsVector< HIT > *inputCollection)
McEventCollection.h
lumiFormat.i
int i
Definition: lumiFormat.py:85
ALFA_PileUpTool::m_meanE_dep
double m_meanE_dep
Definition: ALFA_PileUpTool.h:101
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
ALFA_PileUpTool::fill_OD_DigitCollection
StatusCode fill_OD_DigitCollection(CLHEP::HepRandomEngine *)
Definition: ALFA_PileUpTool.cxx:629
ALFA_PileUpTool::m_SimODHitCollectionName
std::string m_SimODHitCollectionName
Definition: ALFA_PileUpTool.h:91
ALFA_PileUpTool::m_E_ODfib
double m_E_ODfib[8][2][3][30]
Definition: ALFA_PileUpTool.h:87
test_pyathena.parent
parent
Definition: test_pyathena.py:15
sign
int sign(int a)
Definition: TRT_StrawNeighbourSvc.h:107
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
ALFA_ODHitCollection
AtlasHitsVector< ALFA_ODHit > ALFA_ODHitCollection
Definition: ALFA_ODHitCollection.h:13
hist_file_dump.f
f
Definition: hist_file_dump.py:135
ALFA_DigitCollection.h
ALFA_PileUpTool::m_E_fib
double m_E_fib[8][20][64]
Definition: ALFA_PileUpTool.h:86
PileUpToolBase
Definition: PileUpToolBase.h:18
ALFA_PileUpTool::processAllSubEvents
virtual StatusCode processAllSubEvents(const EventContext &ctx) override final
Definition: ALFA_PileUpTool.cxx:117
ALFA_PileUpTool::recordODCollection
StatusCode recordODCollection(ServiceHandle< StoreGateSvc > &evtStore, const std::string &key_ODdigitCnt)
Definition: ALFA_PileUpTool.cxx:342
ALFA_PileUpTool::m_mergedALFA_ODHitList
ALFA_ODHitCollection * m_mergedALFA_ODHitList
Definition: ALFA_PileUpTool.h:108
ALFA_PileUpTool::m_sigma0
double m_sigma0
Definition: ALFA_PileUpTool.h:99
PathResolver.h
ALFA_DigitCollection
Definition: ALFA_DigitCollection.h:15
ALFA_PileUpTool::XTalk
StatusCode XTalk()
Definition: ALFA_PileUpTool.cxx:677
ALFA_PileUpTool::m_sigma1
double m_sigma1
Definition: ALFA_PileUpTool.h:100
ALFA_Hit.h
hancool.filePath
string filePath
Definition: hancool.py:28
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:228
ATHRNG::RNGWrapper
A wrapper class for event-slot-local random engines.
Definition: RNGWrapper.h:56
ALFA_PileUpTool::m_SimHitCollectionName
std::string m_SimHitCollectionName
Definition: ALFA_PileUpTool.h:90
ALFA_ODHitConstIter
AtlasHitsVector< ALFA_ODHit >::const_iterator ALFA_ODHitConstIter
Definition: ALFA_ODHitCollection.h:15
DataVector::push_back
value_type push_back(value_type pElem)
Add an element to the end of the collection.
ALFA_PileUpTool::finalize
virtual StatusCode finalize() override final
Definition: ALFA_PileUpTool.cxx:326
ALFA_PileUpTool::m_cross_talk
double m_cross_talk[8][127]
Definition: ALFA_PileUpTool.h:110
ATHRNG::RNGWrapper::getEngine
CLHEP::HepRandomEngine * getEngine(const EventContext &ctx) const
Retrieve the random engine corresponding to the provided EventContext.
Definition: RNGWrapper.h:134
RNGWrapper.h
ALFA_PileUpTool::m_mergeSvc
ServiceHandle< PileUpMergeSvc > m_mergeSvc
Definition: ALFA_PileUpTool.h:82
python.AthDsoLogger.fname
string fname
Definition: AthDsoLogger.py:67
ALFA_PileUpTool::m_key_DigitCollection
std::string m_key_DigitCollection
Definition: ALFA_PileUpTool.h:92
ALFA_PileUpTool::prepareEvent
virtual StatusCode prepareEvent(const EventContext &ctx, const unsigned int nInputEvents) override final
code taken from ZDC; author (highly) probably John Chapman called before the subevts loop.
Definition: ALFA_PileUpTool.cxx:223
ALFA_HitCollection.h
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
AtlasHitsVector::end
const_iterator end() const
Definition: AtlasHitsVector.h:134
ALFA_PileUpTool::ALFA_PileUpTool
ALFA_PileUpTool(const std::string &type, const std::string &name, const IInterface *parent)
Definition: ALFA_PileUpTool.cxx:50
python.CaloScaleNoiseConfig.type
type
Definition: CaloScaleNoiseConfig.py:78
TimedHitCollection::const_iterator
TimedVector::const_iterator const_iterator
Definition: TimedHitCollection.h:20
ALFA_PileUpTool::m_digitCollection
ALFA_DigitCollection * m_digitCollection
Definition: ALFA_PileUpTool.h:95
AtlasHitsVector::size
size_type size() const
Definition: AtlasHitsVector.h:143
SubEventIterator
std::vector< xAOD::EventInfo::SubEvent >::const_iterator SubEventIterator
Definition: IPileUpTool.h:22
ALFA_PileUpTool::m_mean
double m_mean
Definition: ALFA_PileUpTool.h:103
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
ALFA_PileUpTool::mergeEvent
virtual StatusCode mergeEvent(const EventContext &ctx) override final
return false if not interested in certain xing times (in ns) implemented by default in PileUpToolBase...
Definition: ALFA_PileUpTool.cxx:299
ALFA_ODDigit
Definition: ALFA_ODDigit.h:8
ALFA_PileUpTool::fill_MD_DigitCollection
StatusCode fill_MD_DigitCollection(CLHEP::HepRandomEngine *)
Definition: ALFA_PileUpTool.cxx:537
ALFA_PileUpTool::m_randomSvc
ServiceHandle< IAthRNGSvc > m_randomSvc
Definition: ALFA_PileUpTool.h:83
TimedHitCollection
Definition: TimedHitCollection.h:15
AtlasHitsVector::push_back
void push_back(const T &t)
Definition: AtlasHitsVector.h:145
fitman.k
k
Definition: fitman.py:528
ALFA_PileUpTool::m_stdDev
double m_stdDev
Definition: ALFA_PileUpTool.h:104
ServiceHandle< StoreGateSvc >
ALFA_ODDigitCollection
Definition: ALFA_ODDigitCollection.h:15
Generated on Mon Dec 23 2024 21:06:48 for ATLAS Offline Software by   doxygen 1.8.18