12 #include "CaloDetDescr/CaloDetDescrElement.h"
52 #include "HepPDT/ParticleData.hh"
92 m_larEmID=caloIdManager->getEM_ID();
93 if(m_larEmID==
nullptr)
94 throw std::runtime_error(
"ISF_HitAnalysis: Invalid LAr EM ID helper");
95 m_larFcalID=caloIdManager->getFCAL_ID();
96 if(m_larFcalID==
nullptr)
97 throw std::runtime_error(
"ISF_HitAnalysis: Invalid FCAL ID helper");
98 m_larHecID=caloIdManager->getHEC_ID();
99 if(m_larHecID==
nullptr)
100 throw std::runtime_error(
"ISF_HitAnalysis: Invalid HEC ID helper");
101 m_tileID=caloIdManager->getTileID();
102 if(m_tileID==
nullptr)
103 throw std::runtime_error(
"ISF_HitAnalysis: Invalid Tile ID helper");
108 ATH_CHECK( m_tileSamplingFractionKey.initialize() );
110 ATH_CHECK( m_tileCablingSvc.retrieve() );
111 m_tileCabling = m_tileCablingSvc->cablingService();
116 if (!m_extrapolator.empty() && m_extrapolator.retrieve().isFailure()) {
117 return StatusCode::FAILURE;
126 ATH_CHECK (m_FastCaloSimCaloExtrapolation.retrieve());
134 m_particleDataTable = (HepPDT::ParticleDataTable*) m_partPropSvc->PDT();
135 if(m_particleDataTable ==
nullptr) {
137 return StatusCode::FAILURE;
140 std::unique_ptr<TFile> dummyFile = std::unique_ptr<TFile>(TFile::Open(
"dummyFile.root",
"RECREATE"));
141 m_tree =
new TTree(
"FCS_ParametrizationInput",
"FCS_ParametrizationInput");
142 std::string fullNtupleName =
"/"+m_ntupleFileName+
"/"+m_ntupleTreeName;
143 StatusCode sc = m_thistSvc->regTree(fullNtupleName, m_tree);
144 if (
sc.isFailure() || !m_tree )
146 ATH_MSG_ERROR(
"Unable to register TTree: " << fullNtupleName);
147 return StatusCode::FAILURE;
153 ATH_MSG_INFO(
"Successfull registered TTree: " << fullNtupleName);
155 m_hit_x =
new std::vector<float>;
156 m_hit_y =
new std::vector<float>;
157 m_hit_z =
new std::vector<float>;
158 m_hit_energy =
new std::vector<float>;
159 m_hit_time =
new std::vector<float>;
160 m_hit_identifier =
new std::vector<Long64_t>;
161 m_hit_cellidentifier =
new std::vector<Long64_t>;
162 m_islarbarrel =
new std::vector<bool>;
163 m_islarendcap =
new std::vector<bool>;
164 m_islarhec =
new std::vector<bool>;
165 m_islarfcal =
new std::vector<bool>;
166 m_istile =
new std::vector<bool>;
167 m_hit_sampling =
new std::vector<int>;
168 m_hit_samplingfraction =
new std::vector<float>;
170 m_truth_energy =
new std::vector<float>;
171 m_truth_px =
new std::vector<float>;
172 m_truth_py =
new std::vector<float>;
173 m_truth_pz =
new std::vector<float>;
174 m_truth_pdg =
new std::vector<int>;
175 m_truth_barcode =
new std::vector<int>;
176 m_truth_vtxbarcode =
new std::vector<int>;
178 m_cluster_energy =
new std::vector<float>;
179 m_cluster_eta =
new std::vector<float>;
180 m_cluster_phi =
new std::vector<float>;
181 m_cluster_size =
new std::vector<unsigned>;
182 m_cluster_cellID =
new std::vector<std::vector<Long64_t > >;
184 m_cell_identifier =
new std::vector<Long64_t>;
185 m_cell_energy =
new std::vector<float>;
186 m_cell_sampling =
new std::vector<int>;
188 m_g4hit_energy =
new std::vector<float>;
189 m_g4hit_time =
new std::vector<float>;
190 m_g4hit_identifier =
new std::vector<Long64_t>;
191 m_g4hit_cellidentifier =
new std::vector<Long64_t>;
192 m_g4hit_samplingfraction =
new std::vector<float>;
193 m_g4hit_sampling =
new std::vector<int>;
199 m_final_cell_energy =
new std::vector<Float_t>;
200 m_final_hit_energy =
new std::vector<Float_t>;
201 m_final_g4hit_energy =
new std::vector<Float_t>;
203 m_newTTC_entrance_eta =
new std::vector<std::vector<float> >;
204 m_newTTC_entrance_phi =
new std::vector<std::vector<float> >;
205 m_newTTC_entrance_r =
new std::vector<std::vector<float> >;
206 m_newTTC_entrance_z =
new std::vector<std::vector<float> >;
207 m_newTTC_entrance_detaBorder =
new std::vector<std::vector<float> >;
208 m_newTTC_entrance_OK =
new std::vector<std::vector<bool> >;
209 m_newTTC_back_eta =
new std::vector<std::vector<float> >;
210 m_newTTC_back_phi =
new std::vector<std::vector<float> >;
211 m_newTTC_back_r =
new std::vector<std::vector<float> >;
212 m_newTTC_back_z =
new std::vector<std::vector<float> >;
213 m_newTTC_back_detaBorder =
new std::vector<std::vector<float> >;
214 m_newTTC_back_OK =
new std::vector<std::vector<bool> >;
215 m_newTTC_mid_eta =
new std::vector<std::vector<float> >;
216 m_newTTC_mid_phi =
new std::vector<std::vector<float> >;
217 m_newTTC_mid_r =
new std::vector<std::vector<float> >;
218 m_newTTC_mid_z =
new std::vector<std::vector<float> >;
219 m_newTTC_mid_detaBorder =
new std::vector<std::vector<float> >;
220 m_newTTC_mid_OK =
new std::vector<std::vector<bool> >;
221 m_newTTC_IDCaloBoundary_eta =
new std::vector<float>;
222 m_newTTC_IDCaloBoundary_phi =
new std::vector<float>;
223 m_newTTC_IDCaloBoundary_r =
new std::vector<float>;
224 m_newTTC_IDCaloBoundary_z =
new std::vector<float>;
225 m_newTTC_Angle3D =
new std::vector<float>;
226 m_newTTC_AngleEta =
new std::vector<float>;
228 m_MuonEntryLayer_E =
new std::vector<float>;
229 m_MuonEntryLayer_px =
new std::vector<float>;
230 m_MuonEntryLayer_py =
new std::vector<float>;
231 m_MuonEntryLayer_pz =
new std::vector<float>;
232 m_MuonEntryLayer_x =
new std::vector<float>;
233 m_MuonEntryLayer_y =
new std::vector<float>;
234 m_MuonEntryLayer_z =
new std::vector<float>;
235 m_MuonEntryLayer_pdg =
new std::vector<int>;
238 if(m_saveAllBranches){
239 m_tree->Branch(
"HitX", &m_hit_x);
240 m_tree->Branch(
"HitY", &m_hit_y);
241 m_tree->Branch(
"HitZ", &m_hit_z);
242 m_tree->Branch(
"HitE", &m_hit_energy);
243 m_tree->Branch(
"HitT", &m_hit_time);
244 m_tree->Branch(
"HitIdentifier", &m_hit_identifier);
245 m_tree->Branch(
"HitCellIdentifier", &m_hit_cellidentifier);
246 m_tree->Branch(
"HitIsLArBarrel", &m_islarbarrel);
247 m_tree->Branch(
"HitIsLArEndCap", &m_islarendcap);
248 m_tree->Branch(
"HitIsHEC", &m_islarhec);
249 m_tree->Branch(
"HitIsFCAL", &m_islarfcal);
250 m_tree->Branch(
"HitIsTile", &m_istile);
251 m_tree->Branch(
"HitSampling", &m_hit_sampling);
252 m_tree->Branch(
"HitSamplingFraction", &m_hit_samplingfraction);
254 m_tree->Branch(
"CellIdentifier", &m_cell_identifier);
255 m_tree->Branch(
"CellE", &m_cell_energy);
256 m_tree->Branch(
"CellSampling", &m_cell_sampling);
258 m_tree->Branch(
"G4HitE", &m_g4hit_energy);
259 m_tree->Branch(
"G4HitT", &m_g4hit_time);
260 m_tree->Branch(
"G4HitIdentifier", &m_g4hit_identifier);
261 m_tree->Branch(
"G4HitCellIdentifier", &m_g4hit_cellidentifier);
262 m_tree->Branch(
"G4HitSamplingFraction",&m_g4hit_samplingfraction);
263 m_tree->Branch(
"G4HitSampling", &m_g4hit_sampling);
267 m_tree->Branch(
"TruthE", &m_truth_energy);
268 m_tree->Branch(
"TruthPx", &m_truth_px);
269 m_tree->Branch(
"TruthPy", &m_truth_py);
270 m_tree->Branch(
"TruthPz", &m_truth_pz);
271 m_tree->Branch(
"TruthPDG", &m_truth_pdg);
272 m_tree->Branch(
"TruthBarcode", &m_truth_barcode);
273 m_tree->Branch(
"TruthVtxBarcode", &m_truth_vtxbarcode);
276 m_tree->Branch(
"ClusterE", &m_cluster_energy);
277 m_tree->Branch(
"ClusterEta", &m_cluster_eta);
278 m_tree->Branch(
"ClusterPhi", &m_cluster_phi);
279 m_tree->Branch(
"ClusterSize", &m_cluster_size);
280 m_tree->Branch(
"ClusterCellID", &m_cluster_cellID);
285 m_tree->Branch(
"AllCells", &m_oneeventcells);
290 for (Int_t
i = 0;
i < MAX_LAYER;
i++)
301 m_tree->Branch(
"cell_energy", &m_final_cell_energy);
302 m_tree->Branch(
"hit_energy", &m_final_hit_energy);
303 m_tree->Branch(
"g4hit_energy", &m_final_g4hit_energy);
306 m_tree->Branch(
"total_cell_energy", &m_total_cell_e);
307 m_tree->Branch(
"total_hit_energy", &m_total_hit_e);
308 m_tree->Branch(
"total_g4hit_energy", &m_total_g4hit_e);
311 m_tree->Branch(
"newTTC_back_eta",&m_newTTC_back_eta);
312 m_tree->Branch(
"newTTC_back_phi",&m_newTTC_back_phi);
313 m_tree->Branch(
"newTTC_back_r",&m_newTTC_back_r);
314 m_tree->Branch(
"newTTC_back_z",&m_newTTC_back_z);
315 m_tree->Branch(
"newTTC_back_detaBorder",&m_newTTC_back_detaBorder);
316 m_tree->Branch(
"newTTC_back_OK",&m_newTTC_back_OK);
317 m_tree->Branch(
"newTTC_entrance_eta",&m_newTTC_entrance_eta);
318 m_tree->Branch(
"newTTC_entrance_phi",&m_newTTC_entrance_phi);
319 m_tree->Branch(
"newTTC_entrance_r",&m_newTTC_entrance_r);
320 m_tree->Branch(
"newTTC_entrance_z",&m_newTTC_entrance_z);
321 m_tree->Branch(
"newTTC_entrance_detaBorder",&m_newTTC_entrance_detaBorder);
322 m_tree->Branch(
"newTTC_entrance_OK",&m_newTTC_entrance_OK);
323 m_tree->Branch(
"newTTC_mid_eta",&m_newTTC_mid_eta);
324 m_tree->Branch(
"newTTC_mid_phi",&m_newTTC_mid_phi);
325 m_tree->Branch(
"newTTC_mid_r",&m_newTTC_mid_r);
326 m_tree->Branch(
"newTTC_mid_z",&m_newTTC_mid_z);
327 m_tree->Branch(
"newTTC_mid_detaBorder",&m_newTTC_mid_detaBorder);
328 m_tree->Branch(
"newTTC_mid_OK",&m_newTTC_mid_OK);
329 m_tree->Branch(
"newTTC_IDCaloBoundary_eta",&m_newTTC_IDCaloBoundary_eta);
330 m_tree->Branch(
"newTTC_IDCaloBoundary_phi",&m_newTTC_IDCaloBoundary_phi);
331 m_tree->Branch(
"newTTC_IDCaloBoundary_r",&m_newTTC_IDCaloBoundary_r);
332 m_tree->Branch(
"newTTC_IDCaloBoundary_z",&m_newTTC_IDCaloBoundary_z);
333 m_tree->Branch(
"newTTC_Angle3D",&m_newTTC_Angle3D);
334 m_tree->Branch(
"newTTC_AngleEta",&m_newTTC_AngleEta);
336 m_tree->Branch(
"MuonEntryLayer_E",&m_MuonEntryLayer_E);
337 m_tree->Branch(
"MuonEntryLayer_px",&m_MuonEntryLayer_px);
338 m_tree->Branch(
"MuonEntryLayer_py",&m_MuonEntryLayer_py);
339 m_tree->Branch(
"MuonEntryLayer_pz",&m_MuonEntryLayer_pz);
340 m_tree->Branch(
"MuonEntryLayer_x",&m_MuonEntryLayer_x);
341 m_tree->Branch(
"MuonEntryLayer_y",&m_MuonEntryLayer_y);
342 m_tree->Branch(
"MuonEntryLayer_z",&m_MuonEntryLayer_z);
343 m_tree->Branch(
"MuonEntryLayer_pdg",&m_MuonEntryLayer_pdg);
346 return StatusCode::SUCCESS;
358 return StatusCode::FAILURE;
361 for (
auto attrItr = simParam->begin(); attrItr != simParam->end();
363 std::stringstream outstr;
364 attrItr->toOutputStream(outstr);
372 return StatusCode::FAILURE;
375 for (
auto attrItr = digiParam->begin(); attrItr != digiParam->end();
377 std::stringstream outstr;
378 attrItr->toOutputStream(outstr);
379 ATH_MSG_INFO(
"Digitization MetaData: " << outstr.str());
382 std::unique_ptr<TFile> dummyGeoFile = std::unique_ptr<TFile>(TFile::Open(
"dummyGeoFile.root",
"RECREATE"));
383 TTree*
geo =
new TTree( m_geoModel->atlasVersion().c_str() , m_geoModel->atlasVersion().c_str() );
384 std::string fullNtupleName =
"/"+m_geoFileName+
"/"+m_geoModel->atlasVersion();
386 if(
sc.isFailure() || !
geo )
388 ATH_MSG_ERROR(
"Unable to register TTree: " << fullNtupleName);
389 return StatusCode::FAILURE;
394 using GEOCELL =
struct
398 float eta,
phi,
r,eta_raw,phi_raw,r_raw,
x,
y,
z,x_raw,y_raw,z_raw;
399 float deta,dphi,
dr,
dx,
dy,dz;
402 static GEOCELL geocell;
406 ATH_MSG_INFO(
"Successfull registered TTree: " << fullNtupleName);
409 geo->Branch(
"identifier", &geocell.identifier,
"identifier/L");
410 geo->Branch(
"calosample", &geocell.calosample,
"calosample/I");
412 geo->Branch(
"eta", &geocell.eta,
"eta/F");
413 geo->Branch(
"phi", &geocell.phi,
"phi/F");
414 geo->Branch(
"r", &geocell.r,
"r/F");
415 geo->Branch(
"eta_raw", &geocell.eta_raw,
"eta_raw/F");
416 geo->Branch(
"phi_raw", &geocell.phi_raw,
"phi_raw/F");
417 geo->Branch(
"r_raw", &geocell.r_raw,
"r_raw/F");
419 geo->Branch(
"x", &geocell.x,
"x/F");
420 geo->Branch(
"y", &geocell.y,
"y/F");
421 geo->Branch(
"z", &geocell.z,
"z/F");
422 geo->Branch(
"x_raw", &geocell.x_raw,
"x_raw/F");
423 geo->Branch(
"y_raw", &geocell.y_raw,
"y_raw/F");
424 geo->Branch(
"z_raw", &geocell.z_raw,
"z_raw/F");
426 geo->Branch(
"deta", &geocell.deta,
"deta/F");
427 geo->Branch(
"dphi", &geocell.dphi,
"dphi/F");
428 geo->Branch(
"dr", &geocell.dr,
"dr/F");
429 geo->Branch(
"dx", &geocell.dx,
"dx/F");
430 geo->Branch(
"dy", &geocell.dy,
"dy/F");
431 geo->Branch(
"dz", &geocell.dz,
"dz/F");
448 geocell.identifier=theDDE->identify().get_compact();
449 geocell.calosample=
sample;
450 geocell.eta=theDDE->eta();
451 geocell.phi=theDDE->phi();
452 geocell.r=theDDE->r();
453 geocell.eta_raw=theDDE->eta_raw();
454 geocell.phi_raw=theDDE->phi_raw();
455 geocell.r_raw=theDDE->r_raw();
456 geocell.x=theDDE->x();
457 geocell.y=theDDE->y();
458 geocell.z=theDDE->z();
459 geocell.x_raw=theDDE->x_raw();
460 geocell.y_raw=theDDE->y_raw();
461 geocell.z_raw=theDDE->z_raw();
462 geocell.deta=theDDE->deta();
463 geocell.dphi=theDDE->dphi();
464 geocell.dr=theDDE->dr();
465 geocell.dx=theDDE->dx();
466 geocell.dy=theDDE->dy();
467 geocell.dz=theDDE->dz();
476 dummyGeoFile->Close();
477 return StatusCode::SUCCESS;
489 return StatusCode::FAILURE;
501 vectest.SetPtEtaPhi(1.,1.,1.);
539 std::map<Long64_t, FCS_cell>
cells;
540 std::map<Long64_t, std::vector<FCS_g4hit> > g4hits;
541 std::map<Long64_t, std::vector<FCS_hit> >
hits;
601 if (
sc.isFailure()) {
607 m_hit_x->push_back( (*it)->x() );
608 m_hit_y->push_back( (*it)->y() );
609 m_hit_z->push_back( (*it)->z() );
614 bool larbarrel=
false;
615 bool larendcap=
false;
658 Int_t tile_sampling = -1;
662 if(tile_sampling!= -1) sampling = tile_sampling;
684 sc =
evtStore()->retrieve(mcEvent,
"TruthEvent");
690 if(!mcEvent->
empty()) {
693 int particles_size=(*mcEvent->
begin())->particles_size();
695 loopEnd = particles_size;
698 for (
const auto&
part: *(*mcEvent->
begin()))
700 for (
const auto part: *(*mcEvent->
begin()))
704 ATH_MSG_DEBUG(
"Number truth particles="<<particles_size<<
" loopEnd="<<loopEnd);
707 if (particleIndex>loopEnd)
break;
715 moment.SetXYZ(
part->momentum().px(),
part->momentum().py(),
part->momentum().pz());
716 TVector3 direction=moment.Unit();
728 if((
part)->production_vertex()) {
729 truth.
set_vertex((
part)->production_vertex()->position().
x(), (
part)->production_vertex()->position().
y(), (
part)->production_vertex()->position().
z());
731 truth.
set_vertex(direction.X(),direction.Y(),direction.Z());
732 ATH_MSG_WARNING(
"No particle production vetext, use VERTEX from direction: x "<<direction.X()<<
" y "<<direction.Y()<<
" z "<<direction.Z());
735 if( std::abs(direction.X()-truth.
vertex().X())>0.1 || std::abs(direction.Y()-truth.
vertex().Y())>0.1 || std::abs(direction.Z()-truth.
vertex().Z())>0.1 ) {
736 ATH_MSG_WARNING(
"VERTEX from direction: x "<<direction.X()<<
" y "<<direction.Y()<<
" z "<<direction.Z());
759 std::vector<float> eta_vec_ENT;
760 std::vector<float> phi_vec_ENT;
761 std::vector<float> r_vec_ENT;
762 std::vector<float> z_vec_ENT;
763 std::vector<float> detaBorder_vec_ENT;
764 std::vector<bool> OK_vec_ENT;
766 std::vector<float> eta_vec_EXT;
767 std::vector<float> phi_vec_EXT;
768 std::vector<float> r_vec_EXT;
769 std::vector<float> z_vec_EXT;
770 std::vector<float> detaBorder_vec_EXT;
771 std::vector<bool> OK_vec_EXT;
773 std::vector<float> eta_vec_MID;
774 std::vector<float> phi_vec_MID;
775 std::vector<float> r_vec_MID;
776 std::vector<float> z_vec_MID;
777 std::vector<float> detaBorder_vec_MID;
778 std::vector<bool> OK_vec_MID;
841 sc =
evtStore()->retrieve(MuonEntry,
"MuonEntryLayer");
863 std::string clusterContainerName =
"CaloCalTopoClusters";
864 sc =
evtStore()->retrieve(theClusters, clusterContainerName);
865 if (
sc.isFailure()) {
866 ATH_MSG_WARNING(
" Couldn't get cluster container '" << clusterContainerName <<
"'");
867 return StatusCode::SUCCESS;
871 for ( ; itrClus!=itrLastClus; ++itrClus){
889 unsigned cellcount = 0;
890 std::vector<Long64_t> cellIDs_in_cluster;
893 for ( ;cellIter !=cellIterEnd;cellIter++) {
896 cellIDs_in_cluster.push_back(
cell->ID().get_compact());
897 float EnergyCell=
cell->energy();
906 sc =
evtStore()->retrieve(cellColl,
"AllCalo");
917 for ( ; itrCell!=itrLastCell; ++itrCell)
925 else if (calo_dd_man->
get_element((*itrCell)->ID()))
937 std::string lArKey [4] = {
"LArHitEMB",
"LArHitEMEC",
"LArHitFCAL",
"LArHitHEC"};
938 for (
unsigned int i=0;
i<4;
i++)
946 for (hi=(*iter).begin();hi!=(*iter).end();++hi) {
948 const LArHit* larHit = *hi;
956 float larsampfrac=fSampl->
FSAMPL(larhitid);
965 ATH_MSG_INFO(
"Read "<<hitnumber<<
" G4Hits from "<<lArKey[
i]);
993 for (
int tilesubhit_i = 0; tilesubhit_i<(*i_hit).size(); tilesubhit_i++)
1004 ATH_MSG_INFO(
"Read "<<hitnumber<<
" G4Hits from TileHitVec");
1014 one_cell.cell_identifier = (*m_cell_identifier)[cell_i];
1015 one_cell.sampling = (*m_cell_sampling)[cell_i];
1016 one_cell.energy = (*m_cell_energy)[cell_i];
1017 one_cell.center_x = 0.0;
1018 one_cell.center_y = 0.0;
1019 one_cell.center_z = 0.0;
1020 cells.insert(std::pair<Long64_t, FCS_cell>(one_cell.cell_identifier, one_cell));
1042 one_g4hit.identifier = (*m_g4hit_identifier)[g4hit_i];
1043 one_g4hit.cell_identifier = (*m_g4hit_cellidentifier)[g4hit_i];
1044 one_g4hit.sampling = (*m_g4hit_sampling)[g4hit_i];
1045 one_g4hit.hit_time = (*m_g4hit_time)[g4hit_i];
1047 if (one_g4hit.sampling >= 12 && one_g4hit.sampling <= 20)
1051 one_g4hit.hit_energy = (*m_g4hit_energy)[g4hit_i] * (*m_g4hit_samplingfraction)[g4hit_i];
1053 else one_g4hit.hit_energy = 0.;
1057 one_g4hit.hit_energy = (*m_g4hit_energy)[g4hit_i] / (*m_g4hit_samplingfraction)[g4hit_i];
1059 g4hits.insert(std::pair<Long64_t, std::vector<FCS_g4hit> >(one_g4hit.cell_identifier, std::vector<FCS_g4hit>(1, one_g4hit)));
1064 one_g4hit.identifier = (*m_g4hit_identifier)[g4hit_i];
1065 one_g4hit.cell_identifier = (*m_g4hit_cellidentifier)[g4hit_i];
1066 one_g4hit.sampling = (*m_g4hit_sampling)[g4hit_i];
1067 one_g4hit.hit_time = (*m_g4hit_time)[g4hit_i];
1068 if (one_g4hit.sampling >= 12 && one_g4hit.sampling <= 20)
1072 one_g4hit.hit_energy = (*m_g4hit_energy)[g4hit_i] * (*m_g4hit_samplingfraction)[g4hit_i];
1074 else one_g4hit.hit_energy = 0.;
1078 one_g4hit.hit_energy = (*m_g4hit_energy)[g4hit_i] / (*m_g4hit_samplingfraction)[g4hit_i];
1080 g4hits[(*m_g4hit_cellidentifier)[g4hit_i]].push_back(one_g4hit);
1096 one_hit.identifier = (*m_hit_identifier)[hit_i];
1097 one_hit.cell_identifier = (*m_hit_cellidentifier)[hit_i];
1098 one_hit.sampling = (*m_hit_sampling)[hit_i];
1100 if (one_hit.sampling >= 12 && one_hit.sampling <= 20)
1104 one_hit.hit_energy = (*m_hit_energy)[hit_i] * (*m_hit_samplingfraction)[hit_i];
1106 else one_hit.hit_energy = 0.;
1110 one_hit.hit_energy = (*m_hit_energy)[hit_i] / (*m_hit_samplingfraction)[hit_i];
1113 one_hit.hit_time = (*m_hit_time)[hit_i];
1114 one_hit.hit_x = (*m_hit_x)[hit_i];
1115 one_hit.hit_y = (*m_hit_y)[hit_i];
1116 one_hit.hit_z = (*m_hit_z)[hit_i];
1117 hits.insert(std::pair<Long64_t, std::vector<FCS_hit> >(one_hit.cell_identifier, std::vector<FCS_hit>(1, one_hit)));
1122 one_hit.identifier = (*m_hit_identifier)[hit_i];
1123 one_hit.cell_identifier = (*m_hit_cellidentifier)[hit_i];
1124 one_hit.sampling = (*m_hit_sampling)[hit_i];
1126 if (one_hit.sampling >= 12 && one_hit.sampling <= 20)
1130 one_hit.hit_energy = (*m_hit_energy)[hit_i] * (*m_hit_samplingfraction)[hit_i];
1132 else one_hit.hit_energy = 0.;
1136 one_hit.hit_energy = (*m_hit_energy)[hit_i] / (*m_hit_samplingfraction)[hit_i];
1139 one_hit.hit_time = (*m_hit_time)[hit_i];
1140 one_hit.hit_x = (*m_hit_x)[hit_i];
1141 one_hit.hit_y = (*m_hit_y)[hit_i];
1142 one_hit.hit_z = (*m_hit_z)[hit_i];
1143 hits[(*m_hit_cellidentifier)[hit_i]].push_back(one_hit);
1150 one_matchedcell.
clear();
1152 one_matchedcell.
cell =
it->second;
1154 std::map<Long64_t, std::vector<FCS_hit> >
::iterator it2 =
hits.find(
it->first);
1155 if (it2 !=
hits.end())
1158 one_matchedcell.
hit = it2->second;
1164 one_matchedcell.
hit.clear();
1167 std::map<Long64_t, std::vector<FCS_g4hit> >
::iterator it3 = g4hits.find(
it->first);
1168 if (it3 != g4hits.end())
1170 one_matchedcell.
g4hit = it3->second;
1176 one_matchedcell.
g4hit.clear();
1185 ATH_MSG_DEBUG(
"ISF_HitAnalysis Check after cells: " <<
cells.size() <<
" " << g4hits.size() <<
" " <<
hits.size());
1189 one_matchedcell.
clear();
1192 if (!
it->second.empty())
1206 one_matchedcell.
hit =
it->second;
1207 std::map<Long64_t, std::vector<FCS_g4hit> >
::iterator it3 = g4hits.find(
it->first);
1208 if (it3 != g4hits.end())
1210 one_matchedcell.
g4hit = it3->second;
1216 one_matchedcell.
g4hit.clear();
1224 ATH_MSG_DEBUG(
"ISF_HitAnalysis Check after hits: " <<
cells.size() <<
" " << g4hits.size() <<
" " <<
hits.size());
1225 for (std::map<Long64_t, std::vector<FCS_g4hit> >::
iterator it = g4hits.begin();
it != g4hits.end();)
1227 one_matchedcell.
clear();
1229 if (!
it->second.empty())
1243 one_matchedcell.
g4hit =
it->second;
1244 one_matchedcell.
hit.clear();
1267 for (
unsigned int cellindex = 0; cellindex <
m_layercells[
i]->
size(); cellindex++)
1324 return StatusCode::SUCCESS;
1331 ATH_MSG_DEBUG (
"[ fastCaloSim transport ] processing particle "<<
part.pdg_id() );
1333 std::vector<Trk::HitInfo>*
hitVector =
new std::vector<Trk::HitInfo>;
1335 int pdgId =
part.pdg_id();
1346 auto vtx =
part.production_vertex();
1351 pos =
Amg::Vector3D( vtx->position().x(),vtx->position().y(), vtx->position().z());
1355 ATH_MSG_DEBUG(
"[ fastCaloSim transport ] starting transport from position eta="<<
pos.eta()<<
" phi="<<
pos.phi()<<
" d="<<
pos.mag()<<
" pT="<<
mom.perp() );
1362 double freepath = -1.;
1365 double tDec = freepath > 0. ? freepath : -1.;
1399 ATH_MSG_DEBUG(
"[ fastCaloSim transport ] before calo entrance ");
1411 ATH_MSG_DEBUG(
"[ fastCaloSim transport ] after calo entrance ");
1413 std::unique_ptr<const Trk::TrackParameters> caloEntry =
nullptr;
1417 std::vector<Trk::HitInfo>* dummyHitVector =
nullptr;
1442 ATH_MSG_DEBUG(
"[ fastCaloSim transport ] after calo caloEntry ");
1446 std::unique_ptr<const Trk::TrackParameters> eParameters =
nullptr;
1451 ATH_MSG_DEBUG(
"[ fastCaloSim transport ] starting Calo transport from position eta="<<caloEntry->
position().eta()<<
" phi="<<caloEntry->
position().phi()<<
" d="<<caloEntry->
position().mag() );
1463 eParameters =
m_extrapolator->extrapolateWithPathLimit(*caloEntry,
1479 while (it < hitVector->
end() )
1483 ATH_MSG_DEBUG(
" HIT: layer="<<
sample<<
" sample="<<
sample-3000<<
" eta="<<hitPos.eta()<<
" phi="<<hitPos.phi()<<
" d="<<hitPos.mag());