ISF_HitAnalysis.cxx 67.2 KB
Newer Older
1
/*
2
  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
3
4
*/

5
#include "ISF_FastCaloSimParametrization/ISF_HitAnalysis.h"
6
7
#include "ISF_FastCaloSimEvent/TFCSTruthState.h"
#include "ISF_FastCaloSimEvent/TFCSExtrapolationState.h"
8

9
10
11
12
13
// Section of includes for LAr calo tests
#include "LArSimEvent/LArHitContainer.h"
#include "CaloDetDescr/CaloDetDescrElement.h"
#include "CaloDetDescr/CaloDetDescrManager.h"
#include "GeoAdaptors/GeoLArHit.h"
14
15
#include "GeoModelInterfaces/IGeoModelSvc.h"
#include "AthenaPoolUtilities/AthenaAttributeList.h"
16
17
18
19
20
21
22
23
24
25
26
27
28

// Section of includes for tile calo tests
#include "CaloIdentifier/CaloIdManager.h"
#include "CaloIdentifier/LArEM_ID.h"
#include "CaloIdentifier/LArFCAL_ID.h"
#include "CaloIdentifier/LArHEC_ID.h"
#include "TileConditions/TileInfo.h"

#include "TileDetDescr/TileDetDescrManager.h"
#include "CaloIdentifier/TileID.h"
#include "TileSimEvent/TileHit.h"
#include "TileSimEvent/TileHitVector.h"

29
30
31
//Track Record
#include "TrackRecord/TrackRecordCollection.h"

32
33
//CaloCell
#include "CaloEvent/CaloCellContainer.h"
34
35
36
#include "CaloEvent/CaloClusterCellLinkContainer.h"
#include "xAODCaloEvent/CaloClusterContainer.h"
#include "xAODCaloEvent/CaloCluster.h"
37
38
39
40
41
42
43
44

#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/AlgFactory.h"
#include "GaudiKernel/IToolSvc.h"
#include "GaudiKernel/ITHistSvc.h"

#include "StoreGate/StoreGateSvc.h"

45
#include "ISF_FastCaloSimEvent/FCS_StepInfoCollection.h"
46
47
48
49
50

#include "EventInfo/EventInfo.h"
#include "EventInfo/EventID.h"

#include "TTree.h"
51
#include "TFile.h"
52
#include "TString.h"
53
54
#include "TVector3.h"
#include <sstream>
55
56
57
58

// For MC Truth information:
#include "GeneratorObjects/McEventCollection.h"

59
60
61
62
63
64
65
66

//####################
#include "GaudiKernel/ListItem.h"
#include "CaloDetDescr/CaloDepthTool.h"
#include "TrkParameters/TrackParameters.h"
#include "TrkSurfaces/CylinderSurface.h"
#include "TrkSurfaces/DiscSurface.h"
#include "TrkSurfaces/DiscBounds.h"
67
#include "CaloTrackingGeometry/ICaloSurfaceHelper.h"
68
69
70
#include "CaloTrackingGeometry/ICaloSurfaceBuilder.h"
#include "TrkExInterfaces/IExtrapolator.h"
#include "TrkMaterialOnTrack/EnergyLoss.h"
71
#include "TrkGeometry/TrackingGeometry.h"
72
73
74
75
76
77
#include "GaudiKernel/IPartPropSvc.h"
#include "HepPDT/ParticleData.hh"
#include "HepPDT/ParticleDataTable.hh"
//#########################


78
79
80
81
82
83
#include <algorithm>
#include <math.h>
#include <functional>
#include <iostream>

ISF_HitAnalysis::ISF_HitAnalysis(const std::string& name, ISvcLocator* pSvcLocator)
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
   : AthAlgorithm(name, pSvcLocator)
   //, m_storeGate(0)
   , m_geoModel(0)
   , m_tileInfo(0)
   , m_larEmID(0)
   , m_larFcalID(0)
   , m_larHecID(0)
   , m_tileID(0)
   , m_tileMgr(0)
   , m_hit_x(0)
   , m_hit_y(0)
   , m_hit_z(0)
   , m_hit_energy(0)
   , m_hit_time(0)
   , m_hit_identifier(0)
   , m_hit_cellidentifier(0)
   , m_islarbarrel(0)
   , m_islarendcap(0)
   , m_islarhec(0)
   , m_islarfcal(0)
   , m_istile(0)
   , m_hit_sampling(0)
   , m_hit_samplingfraction(0)
   , m_truth_energy(0)
   , m_truth_px(0)
   , m_truth_py(0)
   , m_truth_pz(0)
   , m_truth_pdg(0)
   , m_truth_barcode(0)
   , m_truth_vtxbarcode(0)
114
115
116
117
118
   , m_cluster_energy(0)
   , m_cluster_eta(0)
   , m_cluster_phi(0)
   , m_cluster_size(0)
   , m_cluster_cellID(0)
119
120
121
122
123
124
125
126
127
   , m_cell_identifier(0)
   , m_cell_energy(0)
   , m_cell_sampling(0)
   , m_g4hit_energy(0)
   , m_g4hit_time(0)
   , m_g4hit_identifier(0)
   , m_g4hit_cellidentifier(0)
   , m_g4hit_samplingfraction(0)
   , m_g4hit_sampling(0)
128
129
130
131
132
133
   , m_total_cell_e(0)
   , m_total_hit_e(0)
   , m_total_g4hit_e(0)
   , m_final_cell_energy(0)
   , m_final_hit_energy(0)
   , m_final_g4hit_energy(0)
134
   , m_tree(0)
135
   , m_ntupleFileName("ISF_HitAnalysis")
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
   , m_ntupleTreeName("CaloHitAna")
   , m_metadataTreeName("MetaData")
   , m_geoFileName("ISF_Geometry")
   , m_thistSvc(0)
   , m_calo_dd_man(0)
   //#####################
   , m_eta_calo_surf(0)
   , m_phi_calo_surf(0)
   , m_d_calo_surf(0)
   , m_ptruth_eta(0)
   , m_ptruth_phi(0)
   , m_ptruth_e(0)
   , m_ptruth_et(0)
   , m_ptruth_pt(0)
   , m_ptruth_p(0)
   , m_pdgid(0)
   , m_newTTC_entrance_eta(0)
   , m_newTTC_entrance_phi(0)
   , m_newTTC_entrance_r(0)
   , m_newTTC_entrance_z(0)
156
157
   , m_newTTC_entrance_detaBorder(0)
   , m_newTTC_entrance_OK(0)
158
159
160
161
   , m_newTTC_back_eta(0)
   , m_newTTC_back_phi(0)
   , m_newTTC_back_r(0)
   , m_newTTC_back_z(0)
162
163
   , m_newTTC_back_detaBorder(0)
   , m_newTTC_back_OK(0)
164
165
166
167
   , m_newTTC_mid_eta(0)
   , m_newTTC_mid_phi(0)
   , m_newTTC_mid_r(0)
   , m_newTTC_mid_z(0)
168
169
   , m_newTTC_mid_detaBorder(0)
   , m_newTTC_mid_OK(0)
170
171
172
173
174
175
176
   , m_newTTC_IDCaloBoundary_eta(0)
   , m_newTTC_IDCaloBoundary_phi(0)
   , m_newTTC_IDCaloBoundary_r(0)
   , m_newTTC_IDCaloBoundary_z(0)
   , m_newTTC_Angle3D(0)
   , m_newTTC_AngleEta(0)

177
178
179
180
181
182
183
184
   , m_MuonEntryLayer_E(0)
   , m_MuonEntryLayer_px(0)
   , m_MuonEntryLayer_py(0)
   , m_MuonEntryLayer_pz(0)
   , m_MuonEntryLayer_x(0)
   , m_MuonEntryLayer_y(0)
   , m_MuonEntryLayer_z(0)
   , m_MuonEntryLayer_pdg(0)
185

186
187
188
189
   , m_caloEntrance(0)
   , m_calo_tb_coord(0)
   , m_sample_calo_surf(CaloCell_ID_FCS::noSample)
   , m_particleDataTable(0)
190
191
   , m_CaloBoundaryR(1148)
   , m_CaloBoundaryZ(3550)
192
193
194

   ,m_MC_DIGI_PARAM("/Digitization/Parameters")
   ,m_MC_SIM_PARAM("/Simulation/Parameters")
195

196
   //######################
197
198
199


  //Note that m_xxx are pointers to vectors set to 0, not set to empty vector! see note around TBranch
200
{
201
202
203
204
  declareProperty("NtupleFileName", m_ntupleFileName);
  declareProperty("NtupleTreeName", m_ntupleTreeName);
  declareProperty("GeoFileName", m_geoFileName);
  declareProperty("MetadataTreeName", m_metadataTreeName);
205
206
  declareProperty("NTruthParticles", m_NtruthParticles=1, "Number of truth particles saved from the truth collection, -1 to save all");

207
208
  declareProperty("FastCaloSimCaloExtrapolation",   m_FastCaloSimCaloExtrapolation );

209
  //###########################
210
211
212
213
  //declareProperty("ExtrapolatorName",               m_extrapolatorName );
  //declareProperty("ExtrapolatorInstanceName",       m_extrapolatorInstanceName );
  //declareProperty("CalosurfMiddleInstanceName",     m_calosurf_InstanceName);
  //declareProperty("CalosurfEntranceInstanceName",   m_calosurf_entrance_InstanceName);
214
215
216
217
218
219

  declareProperty("CaloBoundaryR", m_CaloBoundaryR);
  declareProperty("CaloBoundaryZ", m_CaloBoundaryZ);
  declareProperty("CaloMargin", m_calomargin);
  //######################

220
221
222
223
224
  declareProperty("Extrapolator",                   m_extrapolator );
  declareProperty("CaloEntrance",                   m_caloEntranceName );
  declareProperty("CaloSurfaceHelper",              m_caloSurfaceHelper );

  declareProperty("CaloGeometryHelper",             m_CaloGeometryHelper );
225

226
227
228
  declareProperty("MetaDataSim", m_MC_SIM_PARAM );
  declareProperty("MetaDataDigi", m_MC_DIGI_PARAM );

229
230
  declareProperty("SaveAllBranches", m_saveAllBranches = false);
  declareProperty("DoAllCells", m_doAllCells = false);
231
  declareProperty("DoClusterInfo", m_doClusterInfo = false);
232
233
234
235
236
237
238
239
  declareProperty("DoLayers", m_doLayers = true);
  declareProperty("DoLayerSums", m_doLayerSums = true);
  declareProperty("DoG4Hits", m_doG4Hits = false);
  declareProperty("TimingCut", m_TimingCut = 999999);




240
241
242
243
244
245
  m_surfacelist.resize(0);
  m_surfacelist.push_back(CaloCell_ID_FCS::PreSamplerB);
  m_surfacelist.push_back(CaloCell_ID_FCS::PreSamplerE);
  m_surfacelist.push_back(CaloCell_ID_FCS::EME1);
  m_surfacelist.push_back(CaloCell_ID_FCS::EME2);
  m_surfacelist.push_back(CaloCell_ID_FCS::FCAL0);
246

247
248
249
}


250
ISF_HitAnalysis::~ISF_HitAnalysis()
251
252
{
}
253

254
255
StatusCode ISF_HitAnalysis::geoInit(IOVSVC_CALLBACK_ARGS)
{
256
 ATH_MSG_INFO("geoInit for " << m_geoModel->atlasVersion() );
257

258
259
260
261
262
263
264
265
266
 StatusCode sc = detStore()->retrieve(m_tileMgr);
 if (sc.isFailure())
 {
  ATH_MSG_ERROR( "Unable to retrieve TileDetDescrManager from DetectorStore" );
  m_tileMgr=0;
 }

 sc = detStore()->retrieve(m_tileID);
 if (sc.isFailure())
267
 {
268
269
270
  ATH_MSG_ERROR( "Unable to retrieve TileID helper from DetectorStore" );
  m_tileID=0;
 }
271

272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
 const DataHandle<CaloIdManager> caloIdManager;
 sc=detStore()->retrieve(caloIdManager);
 if(sc.isSuccess())
  ATH_MSG_DEBUG("CaloIDManager retrieved.");
 else
  throw std::runtime_error("ISF_HitAnalysis: Unable to retrieve CaloIDManeger");
 m_larEmID=caloIdManager->getEM_ID();
 if(m_larEmID==0)
  throw std::runtime_error("ISF_HitAnalysis: Invalid LAr EM ID helper");
 m_larFcalID=caloIdManager->getFCAL_ID();
 if(m_larFcalID==0)
  throw std::runtime_error("ISF_HitAnalysis: Invalid FCAL ID helper");
 m_larHecID=caloIdManager->getHEC_ID();
 if(m_larHecID==0)
  throw std::runtime_error("ISF_HitAnalysis: Invalid HEC ID helper");
 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;
 }
296

297
298
299
300
301
302
303
 detStore()->retrieve(m_tileInfo,"TileInfo");
 if(sc.isFailure())
 {
  ATH_MSG_ERROR("Unable to retrieve TileInfo from DetectorStore");
  return StatusCode::FAILURE;
 }
 m_calo_dd_man  = CaloDetDescrManager::instance();
304

305
306
307
308
 // Retrieve Tools
 IToolSvc* p_toolSvc = 0;
 if ( service("ToolSvc",p_toolSvc).isFailure() )
 {
309
        ATH_MSG_ERROR("Cannot find ToolSvc! ");
310
311
312
313
314
  return StatusCode::FAILURE;
 }
 else
 {
  IAlgTool* algTool;
315

316
317
318
  // Get TimedExtrapolator  ***************************************************************************************************
  if (!m_extrapolator.empty() && m_extrapolator.retrieve().isFailure())
   return StatusCode::FAILURE;
319

320
321
322
323
324
  else ATH_MSG_DEBUG("Extrapolator retrieved "<< m_extrapolator);

  // Get CaloSurfaceHelper
  if(m_caloSurfaceHelper.retrieve().isFailure())
   ATH_MSG_INFO("CaloSurfaceHelper not found ");
325

326
327
328
329
330
331
332
  std::string CaloCoordinateTool_name="TBCaloCoordinate";
  ListItem CaloCoordinateTool(CaloCoordinateTool_name);
  if(p_toolSvc->retrieveTool(CaloCoordinateTool.type(),CaloCoordinateTool.name(), algTool, this).isFailure() )
  {
   ATH_MSG_ERROR("Cannot retrieve " << CaloCoordinateTool_name);
   return StatusCode::FAILURE;
  }
333
  m_calo_tb_coord = dynamic_cast<ICaloCoordinateTool*>(algTool);
334
335
336
337
  if(!m_calo_tb_coord )
  {
   ATH_MSG_ERROR("Cannot retrieve " << CaloCoordinateTool_name);
   return StatusCode::FAILURE;
338
339
  }
  else
340
341
   ATH_MSG_INFO("retrieved " << CaloCoordinateTool_name);
 } //tools
342

343
 return StatusCode::SUCCESS;
344

345
346
}

John Chapman's avatar
John Chapman committed
347
348
StatusCode ISF_HitAnalysis::updateMetaData( IOVSVC_CALLBACK_ARGS_P( I, keys ) )
{
349
 ATH_MSG_INFO( "Updating the Sim+Digi MetaData" );
350

351
352
 // Reset the internal settings:
 bool run_update = false;
353

354
355
 // Check what kind of keys we got. In principle the function should only
 // receive the "/Digitization/Parameters" and "/Simulation/Parameters" key.
356
 ATH_MSG_DEBUG("Update called with " <<I<< " folder " << keys.size() << " keys:");
357
358
359
360
361
362
363
364
365
366
 std::list< std::string >::const_iterator itr = keys.begin();
 std::list< std::string >::const_iterator end = keys.end();
 for( ; itr != end; ++itr )
 {
  if( *itr == m_MC_DIGI_PARAM ) run_update = true;
  if( *itr == m_MC_SIM_PARAM ) run_update = true;
 }
 // If that's not the key that we received after all, let's just return
 // silently...
 if( ! run_update ) return StatusCode::SUCCESS;
367

368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
 const DataHandle< AthenaAttributeList > simParam;
 if( detStore()->retrieve( simParam, m_MC_SIM_PARAM ).isFailure() )
 {
   ATH_MSG_WARNING("Retrieving MC SIM metadata failed");
 }
 else
 {
  AthenaAttributeList::const_iterator attr_itr = simParam->begin();
  AthenaAttributeList::const_iterator attr_end = simParam->end();
  for( ; attr_itr != attr_end; ++attr_itr )
  {
   std::stringstream outstr;
   attr_itr->toOutputStream(outstr);
   ATH_MSG_INFO("MetaData: " << outstr.str());
  }
383
384
 }

385
 return StatusCode::SUCCESS;
386
387
388
}


John Chapman's avatar
John Chapman committed
389
390
StatusCode ISF_HitAnalysis::initialize()
{
391
392
393
394
395
396
397
398
399
400
 ATH_MSG_INFO( "Initializing ISF_HitAnalysis" );
 //
 // Register the callback(s):
 //
 StatusCode sc = service("GeoModelSvc", m_geoModel);
 if(sc.isFailure())
 {
  ATH_MSG_ERROR( "Could not locate GeoModelSvc" );
  return StatusCode::FAILURE;
 }
401

402
403
404
 // dummy parameters for the callback:
 int dummyInt=0;
 std::list<std::string> dummyList;
405

406
407
408
 if (m_geoModel->geoInitialized())
 {
  sc=geoInit(dummyInt,dummyList);
409
  if(sc.isFailure())
410
411
412
413
414
415
416
417
418
419
420
421
422
  {
   ATH_MSG_ERROR( "Call to geoInit failed" );
   return StatusCode::FAILURE;
  }
 }
 else
 {
  sc = detStore()->regFcn(&IGeoModelSvc::geoInit, m_geoModel, &ISF_HitAnalysis::geoInit,this);
  if(sc.isFailure())
  {
    ATH_MSG_ERROR( "Could not register geoInit callback" );
    return StatusCode::FAILURE;
  }
423
424
 }

425
426
427
428
429
430
431
432
433
434
435
436
437
 if( detStore()->contains< AthenaAttributeList >( m_MC_DIGI_PARAM ) )
 {
  const DataHandle< AthenaAttributeList > aptr;
  if( detStore()->regFcn( &ISF_HitAnalysis::updateMetaData, this, aptr,m_MC_DIGI_PARAM, true ).isFailure() )
  {
   ATH_MSG_ERROR( "Could not register callback for "<< m_MC_DIGI_PARAM );
   return StatusCode::FAILURE;
  }
 }
 else
 {
  ATH_MSG_WARNING( "MetaData not found for "<< m_MC_DIGI_PARAM );
 }
438
439

 if(detStore()->contains< AthenaAttributeList >( m_MC_SIM_PARAM ) )
440
441
442
443
444
445
446
447
448
449
 {
  const DataHandle< AthenaAttributeList > aptr;
  if( detStore()->regFcn( &ISF_HitAnalysis::updateMetaData, this, aptr,m_MC_SIM_PARAM, true ).isFailure() )
  {
   ATH_MSG_ERROR( "Could not register callback for "<< m_MC_SIM_PARAM );
   return StatusCode::FAILURE;
  }
 }
 else
  ATH_MSG_WARNING( "MetaData not found for "<< m_MC_SIM_PARAM );
450

451
452
453
454
455
456
 // Get CaloGeometryHelper
 if (m_CaloGeometryHelper.retrieve().isFailure())
 {
  ATH_MSG_ERROR("CaloGeometryHelper not found ");
  return StatusCode::FAILURE;
 }
457

458
459
460
461
462
463
 // Get FastCaloSimCaloExtrapolation
 if (m_FastCaloSimCaloExtrapolation.retrieve().isFailure())
 {
  ATH_MSG_ERROR("FastCaloSimCaloExtrapolation not found ");
  return StatusCode::FAILURE;
 }
464

465
466
467
468
469
470
471
 // Grab the Ntuple and histogramming service for the tree
 sc = service("THistSvc",m_thistSvc);
 if (sc.isFailure())
 {
  ATH_MSG_ERROR( "Unable to retrieve pointer to THistSvc" );
  return StatusCode::FAILURE;
 }
472

473
474
475
476
477
478
479
 //#########################
 IPartPropSvc* p_PartPropSvc=0;
 if (service("PartPropSvc",p_PartPropSvc).isFailure() || p_PartPropSvc == 0)
 {
  ATH_MSG_ERROR("could not find PartPropService");
  return StatusCode::FAILURE;
 }
480

481
482
483
 m_particleDataTable = (HepPDT::ParticleDataTable*) p_PartPropSvc->PDT();
 if(m_particleDataTable == 0)
 {
484
        ATH_MSG_ERROR("PDG table not found");
485
486
487
  return StatusCode::FAILURE;
 }
 //#########################
488
 std::unique_ptr<TFile> dummyFile = std::unique_ptr<TFile>(TFile::Open("dummyFile.root", "RECREATE")); //This is added to suppress the error messages about memory-resident trees
489
490
 m_tree = new TTree("FCS_ParametrizationInput", "FCS_ParametrizationInput");
 std::string fullNtupleName =  "/"+m_ntupleFileName+"/"+m_ntupleTreeName;
491
 sc = m_thistSvc->regTree(fullNtupleName, m_tree);
492
 if (sc.isFailure() || !m_tree )
493
494
 {
  ATH_MSG_ERROR("Unable to register TTree: " << fullNtupleName);
495
496
  return StatusCode::FAILURE;
 }
497

498
499
500
 /** now add branches and leaves to the tree */
 if (m_tree)
 {
501
  ATH_MSG_INFO("Successfull registered TTree: " << fullNtupleName);
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
  //initialize the variables before creating the branches
  m_hit_x = new std::vector<float>;
  m_hit_y = new std::vector<float>;
  m_hit_z = new std::vector<float>;
  m_hit_energy = new std::vector<float>;
  m_hit_time = new std::vector<float>;
  m_hit_identifier = new std::vector<Long64_t>;
  m_hit_cellidentifier = new std::vector<Long64_t>;
  m_islarbarrel = new std::vector<bool>;
  m_islarendcap = new std::vector<bool>;
  m_islarhec = new std::vector<bool>;
  m_islarfcal = new std::vector<bool>;
  m_istile = new std::vector<bool>;
  m_hit_sampling = new std::vector<int>;
  m_hit_samplingfraction = new std::vector<float>;

  m_truth_energy = new std::vector<float>;
  m_truth_px = new std::vector<float>;
  m_truth_py = new std::vector<float>;
  m_truth_pz = new std::vector<float>;
  m_truth_pdg = new std::vector<int>;
  m_truth_barcode = new std::vector<int>;
  m_truth_vtxbarcode = new std::vector<int>;

526
527
528
529
530
531
532
  m_cluster_energy = new std::vector<float>;
  m_cluster_eta    = new std::vector<float>;
  m_cluster_phi    = new std::vector<float>;
  m_cluster_size   = new std::vector<unsigned>;
  m_cluster_cellID = new std::vector<std::vector<Long64_t > >;


533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
  m_cell_identifier = new std::vector<Long64_t>;
  m_cell_energy = new std::vector<float>;
  m_cell_sampling = new std::vector<int>;

  m_g4hit_energy = new std::vector<float>;
  m_g4hit_time = new std::vector<float>;
  m_g4hit_identifier = new std::vector<Long64_t>;
  m_g4hit_cellidentifier = new std::vector<Long64_t>;
  m_g4hit_samplingfraction = new std::vector<float>;
  m_g4hit_sampling = new std::vector<int>;

  m_total_cell_e = 0;
  m_total_hit_e = 0;
  m_total_g4hit_e = 0;

  m_final_cell_energy = new std::vector<Float_t>;
  m_final_hit_energy = new std::vector<Float_t>;
  m_final_g4hit_energy = new std::vector<Float_t>;

  m_newTTC_entrance_eta = new std::vector<std::vector<float> >;
  m_newTTC_entrance_phi = new std::vector<std::vector<float> >;
  m_newTTC_entrance_r = new std::vector<std::vector<float> >;
  m_newTTC_entrance_z = new std::vector<std::vector<float> >;
556
557
  m_newTTC_entrance_detaBorder = new std::vector<std::vector<float> >;
  m_newTTC_entrance_OK = new std::vector<std::vector<bool> >;
558
559
560
561
  m_newTTC_back_eta = new std::vector<std::vector<float> >;
  m_newTTC_back_phi = new std::vector<std::vector<float> >;
  m_newTTC_back_r = new std::vector<std::vector<float> >;
  m_newTTC_back_z = new std::vector<std::vector<float> >;
562
563
  m_newTTC_back_detaBorder = new std::vector<std::vector<float> >;
  m_newTTC_back_OK = new std::vector<std::vector<bool> >;
564
565
566
567
  m_newTTC_mid_eta = new std::vector<std::vector<float> >;
  m_newTTC_mid_phi = new std::vector<std::vector<float> >;
  m_newTTC_mid_r = new std::vector<std::vector<float> >;
  m_newTTC_mid_z = new std::vector<std::vector<float> >;
568
569
  m_newTTC_mid_detaBorder = new std::vector<std::vector<float> >;
  m_newTTC_mid_OK = new std::vector<std::vector<bool> >;
570
571
572
573
574
575
576
  m_newTTC_IDCaloBoundary_eta = new std::vector<float>;
  m_newTTC_IDCaloBoundary_phi = new std::vector<float>;
  m_newTTC_IDCaloBoundary_r = new std::vector<float>;
  m_newTTC_IDCaloBoundary_z = new std::vector<float>;
  m_newTTC_Angle3D = new std::vector<float>;
  m_newTTC_AngleEta = new std::vector<float>;

577
578
579
580
581
582
583
584
  m_MuonEntryLayer_E = new std::vector<float>;
  m_MuonEntryLayer_px = new std::vector<float>;
  m_MuonEntryLayer_py = new std::vector<float>;
  m_MuonEntryLayer_pz = new std::vector<float>;
  m_MuonEntryLayer_x = new std::vector<float>;
  m_MuonEntryLayer_y = new std::vector<float>;
  m_MuonEntryLayer_z = new std::vector<float>;
  m_MuonEntryLayer_pdg = new std::vector<int>;
585
586


587
  // Optional branches
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
  if(m_saveAllBranches){
    m_tree->Branch("HitX",                 &m_hit_x);
    m_tree->Branch("HitY",                 &m_hit_y);
    m_tree->Branch("HitZ",                 &m_hit_z);
    m_tree->Branch("HitE",                 &m_hit_energy);
    m_tree->Branch("HitT",                 &m_hit_time);
    m_tree->Branch("HitIdentifier",        &m_hit_identifier);
    m_tree->Branch("HitCellIdentifier",    &m_hit_cellidentifier);
    m_tree->Branch("HitIsLArBarrel",       &m_islarbarrel);
    m_tree->Branch("HitIsLArEndCap",       &m_islarendcap);
    m_tree->Branch("HitIsHEC",             &m_islarhec);
    m_tree->Branch("HitIsFCAL",            &m_islarfcal);
    m_tree->Branch("HitIsTile",            &m_istile);
    m_tree->Branch("HitSampling",          &m_hit_sampling);
    m_tree->Branch("HitSamplingFraction",  &m_hit_samplingfraction);

    m_tree->Branch("CellIdentifier",       &m_cell_identifier);
    m_tree->Branch("CellE",                &m_cell_energy);
    m_tree->Branch("CellSampling",         &m_cell_sampling);

    m_tree->Branch("G4HitE",               &m_g4hit_energy);
    m_tree->Branch("G4HitT",               &m_g4hit_time);
    m_tree->Branch("G4HitIdentifier",      &m_g4hit_identifier);
    m_tree->Branch("G4HitCellIdentifier",  &m_g4hit_cellidentifier);
    m_tree->Branch("G4HitSamplingFraction",&m_g4hit_samplingfraction);
    m_tree->Branch("G4HitSampling",        &m_g4hit_sampling);
  }
615

616
  //CaloHitAna output variables
617
618
619
620
621
622
623
  m_tree->Branch("TruthE",               &m_truth_energy);
  m_tree->Branch("TruthPx",              &m_truth_px);
  m_tree->Branch("TruthPy",              &m_truth_py);
  m_tree->Branch("TruthPz",              &m_truth_pz);
  m_tree->Branch("TruthPDG",             &m_truth_pdg);
  m_tree->Branch("TruthBarcode",         &m_truth_barcode);
  m_tree->Branch("TruthVtxBarcode",      &m_truth_vtxbarcode);
624

625
626
627
628
629
630
631
  if(m_doClusterInfo){
    m_tree->Branch("ClusterE",               &m_cluster_energy);
    m_tree->Branch("ClusterEta",             &m_cluster_eta);
    m_tree->Branch("ClusterPhi",             &m_cluster_phi);
    m_tree->Branch("ClusterSize",            &m_cluster_size);
    m_tree->Branch("ClusterCellID",          &m_cluster_cellID);
  }
632
633
634
635
  m_oneeventcells = new FCS_matchedcellvector;
  if(m_doAllCells){
    m_tree->Branch("AllCells", &m_oneeventcells);
  }
636

637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
  //write cells per layer
  if(m_doLayers){
    for (Int_t i = 0; i < MAX_LAYER; i++)
    {
      TString branchname = "Sampling_";
      branchname += i;
      m_layercells[i] = new FCS_matchedcellvector;
      m_tree->Branch(branchname, &m_layercells[i]);
    }
  }

  if(m_doLayerSums){
    //write also energies per layer:
    m_tree->Branch("cell_energy", &m_final_cell_energy);
    m_tree->Branch("hit_energy",  &m_final_hit_energy);
    m_tree->Branch("g4hit_energy", &m_final_g4hit_energy);

    //This is a duplicate of cell_energy[25]
    m_tree->Branch("total_cell_energy", &m_total_cell_e);
    m_tree->Branch("total_hit_energy",  &m_total_hit_e);
    m_tree->Branch("total_g4hit_energy", &m_total_g4hit_e);
  }
659

660
661
662
663
  m_tree->Branch("newTTC_back_eta",&m_newTTC_back_eta);
  m_tree->Branch("newTTC_back_phi",&m_newTTC_back_phi);
  m_tree->Branch("newTTC_back_r",&m_newTTC_back_r);
  m_tree->Branch("newTTC_back_z",&m_newTTC_back_z);
664
665
  m_tree->Branch("newTTC_back_detaBorder",&m_newTTC_back_detaBorder);
  m_tree->Branch("newTTC_back_OK",&m_newTTC_back_OK);
666
667
668
669
  m_tree->Branch("newTTC_entrance_eta",&m_newTTC_entrance_eta);
  m_tree->Branch("newTTC_entrance_phi",&m_newTTC_entrance_phi);
  m_tree->Branch("newTTC_entrance_r",&m_newTTC_entrance_r);
  m_tree->Branch("newTTC_entrance_z",&m_newTTC_entrance_z);
670
671
  m_tree->Branch("newTTC_entrance_detaBorder",&m_newTTC_entrance_detaBorder);
  m_tree->Branch("newTTC_entrance_OK",&m_newTTC_entrance_OK);
672
673
674
675
  m_tree->Branch("newTTC_mid_eta",&m_newTTC_mid_eta);
  m_tree->Branch("newTTC_mid_phi",&m_newTTC_mid_phi);
  m_tree->Branch("newTTC_mid_r",&m_newTTC_mid_r);
  m_tree->Branch("newTTC_mid_z",&m_newTTC_mid_z);
676
677
  m_tree->Branch("newTTC_mid_detaBorder",&m_newTTC_mid_detaBorder);
  m_tree->Branch("newTTC_mid_OK",&m_newTTC_mid_OK);
678
679
680
681
682
683
  m_tree->Branch("newTTC_IDCaloBoundary_eta",&m_newTTC_IDCaloBoundary_eta);
  m_tree->Branch("newTTC_IDCaloBoundary_phi",&m_newTTC_IDCaloBoundary_phi);
  m_tree->Branch("newTTC_IDCaloBoundary_r",&m_newTTC_IDCaloBoundary_r);
  m_tree->Branch("newTTC_IDCaloBoundary_z",&m_newTTC_IDCaloBoundary_z);
  m_tree->Branch("newTTC_Angle3D",&m_newTTC_Angle3D);
  m_tree->Branch("newTTC_AngleEta",&m_newTTC_AngleEta);
684

685
686


687
688
689
690
691
692
693
694
  m_tree->Branch("MuonEntryLayer_E",&m_MuonEntryLayer_E);
  m_tree->Branch("MuonEntryLayer_px",&m_MuonEntryLayer_px);
  m_tree->Branch("MuonEntryLayer_py",&m_MuonEntryLayer_py);
  m_tree->Branch("MuonEntryLayer_pz",&m_MuonEntryLayer_pz);
  m_tree->Branch("MuonEntryLayer_x",&m_MuonEntryLayer_x);
  m_tree->Branch("MuonEntryLayer_y",&m_MuonEntryLayer_y);
  m_tree->Branch("MuonEntryLayer_z",&m_MuonEntryLayer_z);
  m_tree->Branch("MuonEntryLayer_pdg",&m_MuonEntryLayer_pdg);
695
696


697
 }
698
 dummyFile->Close();
699
 return StatusCode::SUCCESS;
700

701
} //initialize
John Chapman's avatar
John Chapman committed
702
703
704

StatusCode ISF_HitAnalysis::finalize()
{
705

706
 ATH_MSG_INFO( "doing finalize()" );
707
 std::unique_ptr<TFile> dummyGeoFile = std::unique_ptr<TFile>(TFile::Open("dummyGeoFile.root", "RECREATE")); //This is added to suppress the error messages about memory-resident trees
708
 TTree* geo = new TTree( m_geoModel->atlasVersion().c_str() , m_geoModel->atlasVersion().c_str() );
709
710
 std::string fullNtupleName =  "/"+m_geoFileName+"/"+m_geoModel->atlasVersion();
 StatusCode sc = m_thistSvc->regTree(fullNtupleName, geo);
711
 if(sc.isFailure() || !geo )
712
713
 {
  ATH_MSG_ERROR("Unable to register TTree: " << fullNtupleName);
714
715
  return StatusCode::FAILURE;
 }
716

717
 /** now add branches and leaves to the tree */
718

719
720
721
722
723
724
725
 typedef struct
 {
  Long64_t identifier;
  Int_t calosample;
  float eta,phi,r,eta_raw,phi_raw,r_raw,x,y,z,x_raw,y_raw,z_raw;
  float deta,dphi,dr,dx,dy,dz;
 } GEOCELL;
726
727
728

 static GEOCELL geocell;

729
730
 if(geo)
 {
731
  ATH_MSG_INFO("Successfull registered TTree: " << fullNtupleName);
732
733
734
735
  //this actually creates the vector itself! And only if it succeeds! Note that the result is not checked! And the code is probably leaking memory in the end
  //geo->Branch("cells", &geocell,"identifier/L:eta,phi,r,eta_raw,phi_raw,r_raw,x,y,z,x_raw,y_raw,z_raw/F:Deta,Dphi,Dr,Dx,Dy,Dz/F");
  geo->Branch("identifier", &geocell.identifier,"identifier/L");
  geo->Branch("calosample", &geocell.calosample,"calosample/I");
736

737
738
739
740
741
742
  geo->Branch("eta", &geocell.eta,"eta/F");
  geo->Branch("phi", &geocell.phi,"phi/F");
  geo->Branch("r", &geocell.r,"r/F");
  geo->Branch("eta_raw", &geocell.eta_raw,"eta_raw/F");
  geo->Branch("phi_raw", &geocell.phi_raw,"phi_raw/F");
  geo->Branch("r_raw", &geocell.r_raw,"r_raw/F");
743

744
745
746
747
748
749
  geo->Branch("x", &geocell.x,"x/F");
  geo->Branch("y", &geocell.y,"y/F");
  geo->Branch("z", &geocell.z,"z/F");
  geo->Branch("x_raw", &geocell.x_raw,"x_raw/F");
  geo->Branch("y_raw", &geocell.y_raw,"y_raw/F");
  geo->Branch("z_raw", &geocell.z_raw,"z_raw/F");
750

751
752
753
754
755
756
757
  geo->Branch("deta", &geocell.deta,"deta/F");
  geo->Branch("dphi", &geocell.dphi,"dphi/F");
  geo->Branch("dr", &geocell.dr,"dr/F");
  geo->Branch("dx", &geocell.dx,"dx/F");
  geo->Branch("dy", &geocell.dy,"dy/F");
  geo->Branch("dz", &geocell.dz,"dz/F");
 }
758

759
760
761
762
 if(m_calo_dd_man)
 {
  int ncells=0;
  for(CaloDetDescrManager::calo_element_const_iterator calo_iter=m_calo_dd_man->element_begin();calo_iter<m_calo_dd_man->element_end();++calo_iter)
John Chapman's avatar
John Chapman committed
763
  {
764
765
766
767
768
769
770
   const CaloDetDescrElement* theDDE=*calo_iter;
   if(theDDE)
   {
    CaloCell_ID::CaloSample sample=theDDE->getSampling();
    //CaloCell_ID::SUBCALO calo=theDDE->getSubCalo();
    ++ncells;
    if(geo)
John Chapman's avatar
John Chapman committed
771
    {
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
     geocell.identifier=theDDE->identify().get_compact();
     geocell.calosample=sample;
     geocell.eta=theDDE->eta();
     geocell.phi=theDDE->phi();
     geocell.r=theDDE->r();
     geocell.eta_raw=theDDE->eta_raw();
     geocell.phi_raw=theDDE->phi_raw();
     geocell.r_raw=theDDE->r_raw();
     geocell.x=theDDE->x();
     geocell.y=theDDE->y();
     geocell.z=theDDE->z();
     geocell.x_raw=theDDE->x_raw();
     geocell.y_raw=theDDE->y_raw();
     geocell.z_raw=theDDE->z_raw();
     geocell.deta=theDDE->deta();
     geocell.dphi=theDDE->dphi();
     geocell.dr=theDDE->dr();
     geocell.dx=theDDE->dx();
     geocell.dy=theDDE->dy();
     geocell.dz=theDDE->dz();
792

793
     geo->Fill();
794
    }
795
   }
796
  }
797

798
799
  ATH_MSG_INFO( ncells<<" cells found" );
 }
800
 dummyGeoFile->Close();
801
 return StatusCode::SUCCESS;
John Chapman's avatar
John Chapman committed
802
803
804
805
806
} //finalize


StatusCode ISF_HitAnalysis::execute()
{
807

808
 ATH_MSG_DEBUG( "In ISF_HitAnalysis::execute()" );
809

810
811
812
813
814
 if (! m_tree)
 {
  ATH_MSG_ERROR( "tree not registered" );
  return StatusCode::FAILURE;
 }
815
816

 //now if the branches were created correctly, the pointers point to something and it is possible to clear the vectors
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
 TVector3 vectest;
 vectest.SetPtEtaPhi(1.,1.,1.);
 m_hit_x->clear();
 m_hit_y->clear();
 m_hit_z->clear();
 m_hit_energy->clear();
 m_hit_time->clear();
 m_hit_identifier->clear();
 m_hit_cellidentifier->clear();
 m_islarbarrel->clear();
 m_islarendcap->clear();
 m_islarhec->clear();
 m_islarfcal->clear();
 m_istile->clear();
 m_hit_sampling->clear();
 m_hit_samplingfraction->clear();
 m_truth_energy->clear();
 m_truth_px->clear();
 m_truth_py->clear();
 m_truth_pz->clear();
 m_truth_pdg->clear();
 m_truth_barcode->clear();
 m_truth_vtxbarcode->clear();
840
841
842
843
844
 m_cluster_energy->clear();
 m_cluster_eta->clear();
 m_cluster_phi->clear();
 m_cluster_size->clear();
 m_cluster_cellID->clear();
845
846
847
848
849
850
851
852
853
854
855
 m_cell_identifier->clear();
 m_cell_energy->clear();
 m_cell_sampling->clear();
 m_g4hit_energy->clear();
 m_g4hit_time->clear();
 m_g4hit_identifier->clear();
 m_g4hit_cellidentifier->clear();
 m_g4hit_sampling->clear();
 m_g4hit_samplingfraction->clear();
 //which fails for this one!!
 //m_matched_cells->clear();
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
 std::map<Long64_t, FCS_cell> cells; //read all objects and collect them by identifier (Long64_t)
 std::map<Long64_t, std::vector<FCS_g4hit> > g4hits;
 std::map<Long64_t, std::vector<FCS_hit> > hits;

 cells.clear();
 g4hits.clear();
 hits.clear();

 FCS_cell   one_cell; //note that this is not extra safe if I don't have a clear method!
 FCS_g4hit  one_g4hit;
 FCS_hit    one_hit;
 FCS_matchedcell one_matchedcell;

 m_oneeventcells->m_vector.clear();
 m_final_g4hit_energy->clear();
 m_final_hit_energy->clear();
 m_final_cell_energy->clear();
873

874
875
876
877
 m_newTTC_back_eta->clear();
 m_newTTC_back_phi->clear();
 m_newTTC_back_r->clear();
 m_newTTC_back_z->clear();
878
879
 m_newTTC_back_detaBorder->clear();
 m_newTTC_back_OK->clear();
880
881
882
883
 m_newTTC_entrance_eta->clear();
 m_newTTC_entrance_phi->clear();
 m_newTTC_entrance_r->clear();
 m_newTTC_entrance_z->clear();
884
885
 m_newTTC_entrance_detaBorder->clear();
 m_newTTC_entrance_OK->clear();
886
887
888
889
 m_newTTC_mid_eta->clear();
 m_newTTC_mid_phi->clear();
 m_newTTC_mid_r->clear();
 m_newTTC_mid_z->clear();
890
891
 m_newTTC_mid_detaBorder->clear();
 m_newTTC_mid_OK->clear();
892
893
894
895
896
897
 m_newTTC_IDCaloBoundary_eta->clear();
 m_newTTC_IDCaloBoundary_phi->clear();
 m_newTTC_IDCaloBoundary_r->clear();
 m_newTTC_IDCaloBoundary_z->clear();
 m_newTTC_Angle3D->clear();
 m_newTTC_AngleEta->clear();
898
899


900
901
902
903
904
905
906
 m_MuonEntryLayer_E->clear();
 m_MuonEntryLayer_x->clear();
 m_MuonEntryLayer_y->clear();
 m_MuonEntryLayer_z->clear();
 m_MuonEntryLayer_px->clear();
 m_MuonEntryLayer_py->clear();
 m_MuonEntryLayer_pz->clear();
907

908
 //##########################
909

910
911
 //Get the FastCaloSim step info collection from store
 const ISF_FCS_Parametrization::FCS_StepInfoCollection* eventStepsES;
John Chapman's avatar
John Chapman committed
912
 StatusCode sc = evtStore()->retrieve(eventStepsES, "MergedEventSteps");
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
 if (sc.isFailure()) {
   ATH_MSG_WARNING( "No FastCaloSim steps read from StoreGate?" );
   //return StatusCode::FAILURE;
 } else {
   ATH_MSG_INFO("Read: "<<eventStepsES->size()<<" position hits");
   for (ISF_FCS_Parametrization::FCS_StepInfoCollection::const_iterator it = eventStepsES->begin(); it != eventStepsES->end(); ++it) {
     m_hit_x->push_back( (*it)->x() );
     m_hit_y->push_back( (*it)->y() );
     m_hit_z->push_back( (*it)->z() );
     m_hit_energy->push_back( (*it)->energy() );
     m_hit_time->push_back( (*it)->time());

     //Try to get the samplings, sampling fractions from identifiers
     bool larbarrel=false;
     bool larendcap=false;
     bool larhec=false;
     bool larfcal=false;
     bool tile=false;
     int sampling=-1;
     double sampfrac=0.0;

     Identifier id = (*it)->identify();
     Identifier cell_id = (*it)->identify(); //to be replaced by cell_id in tile

     if(m_calo_dd_man->get_element(id)) {
       CaloCell_ID::CaloSample layer = m_calo_dd_man->get_element(id)->getSampling();
       sampling = layer; //use CaloCell layer immediately
     } else {
       ATH_MSG_WARNING( "Warning no sampling info for "<<id.getString());
942
     }
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973

     if(m_larEmID->is_lar_em(id) || m_larHecID->is_lar_hec(id) || m_larFcalID->is_lar_fcal(id)) sampfrac=m_dd_fSampl->FSAMPL(id);

     if(m_larEmID->is_lar_em(id)) {
       //LAr EM cells
       if (m_larEmID->is_em_barrel(id)) larbarrel=true;
        else if(m_larEmID->is_em_endcap(id)) larendcap=true;
     } else if(m_larHecID->is_lar_hec(id)) {
       //LAr HEC cells
       larhec = true;
     } else if(m_larFcalID->is_lar_fcal(id)) {
       //LAr FCal cells
       larfcal = true;
     } else if (m_tileID->is_tile_aux(id)) {
       // special case for E4'
       tile = true;
       cell_id = m_tileID->cell_id(id);
       sampling = CaloCell_ID::TileGap3;
       sampfrac = m_tileInfo->HitCalib(id);
     } else if(m_tileID->is_tile_barrel(id) || m_tileID->is_tile_extbarrel(id) || m_tileID->is_tile_gap(id)) {
       // all other Tile cells
       tile = true;
       cell_id = m_tileID->cell_id(id);
       Int_t tile_sampling = -1;
       if(m_calo_dd_man->get_element(cell_id)) {
         tile_sampling = m_calo_dd_man->get_element(cell_id)->getSampling();
         sampfrac = m_tileInfo->HitCalib(cell_id);
       }
       if(tile_sampling!= -1) sampling = tile_sampling; //m_calo_dd_man needs to be called with cell_id not pmt_id!!
     } else {
       ATH_MSG_WARNING( "This hit is somewhere. Please check!");
974
     }
975
976
977
978
979
980
981
982
983
984
985
986
987

     m_hit_identifier->push_back(id.get_compact());
     m_hit_cellidentifier->push_back(cell_id.get_compact());
     //push things into vectors:
     m_islarbarrel->push_back(larbarrel);
     m_islarendcap->push_back(larendcap);
     m_islarhec->push_back(larhec);
     m_islarfcal->push_back(larfcal);
     m_istile->push_back(tile);
     m_hit_sampling->push_back(sampling);
     m_hit_samplingfraction->push_back(sampfrac);

   } //event steps
988
 }//event steps read correctly
989

990
991
992
993
 //Get truth particle info
 //Note that there can be more truth particles, the first one is usually the one we need.
 const DataHandle<McEventCollection> mcEvent;
 sc = evtStore()->retrieve(mcEvent,"TruthEvent");
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
 if(sc.isFailure()) {
   ATH_MSG_WARNING( "No truth event!");
 } else {
   if(mcEvent) {
     //std::cout<<"ISF_HitAnalysis: MC event size: "<<mcEvent->size()<<std::endl;
     if(mcEvent->size()) {
       int particleIndex=0;
       int loopEnd = m_NtruthParticles;
       if(loopEnd==-1) {
         loopEnd = (*mcEvent->begin())->particles_size(); //is this the correct thing?
       }
       //std::cout <<"ISF_HitAnalysis: MC first truth event size: "<<(*mcEvent->begin())->particles_size()<<std::endl;
       for (HepMC::GenEvent::particle_const_iterator it = (*mcEvent->begin())->particles_begin(); it != (*mcEvent->begin())->particles_end(); ++it) {
         ATH_MSG_DEBUG("Number truth particles="<<(*mcEvent->begin())->particles_size()<<" loopEnd="<<loopEnd);
         particleIndex++;

         if (particleIndex>loopEnd) break; //enough particles

         //UPDATE EXTRAPOLATION WITH ALGTOOL***********************************************

         TFCSTruthState truth((*it)->momentum().px(),(*it)->momentum().py(),(*it)->momentum().pz(),(*it)->momentum().e(),(*it)->pdg_id());

         //calculate the vertex
         TVector3 moment;
         moment.SetXYZ((*it)->momentum().px(),(*it)->momentum().py(),(*it)->momentum().pz());
         TVector3 direction=moment.Unit();

         //does it hit the barrel or the EC?
1022
1023

         if(abs(direction.Z())/m_CaloBoundaryZ < direction.Perp()/m_CaloBoundaryR) {
1024
           //BARREL
1025
           direction*=m_CaloBoundaryR/direction.Perp();
1026
1027
         } else {
           //EC
1028
           direction*=m_CaloBoundaryZ/abs(direction.Z());
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
         }  

         if((*it)->production_vertex()) {
           truth.set_vertex((*it)->production_vertex()->point3d().x(), (*it)->production_vertex()->point3d().y(), (*it)->production_vertex()->point3d().z());
         } else {
           truth.set_vertex(direction.X(),direction.Y(),direction.Z());
           ATH_MSG_WARNING("No particle production vetext, use VERTEX from direction: x "<<direction.X()<<" y "<<direction.Y()<<" z "<<direction.Z());
         }  
         
         if( fabs(direction.X()-truth.vertex().X())>0.1 || fabs(direction.Y()-truth.vertex().Y())>0.1 || fabs(direction.Z()-truth.vertex().Z())>0.1 ) {
           ATH_MSG_WARNING("VERTEX from direction: x "<<direction.X()<<" y "<<direction.Y()<<" z "<<direction.Z());
           ATH_MSG_WARNING("but VERTEX from hepmc: x "<<truth.vertex().X()<<" y "<<truth.vertex().Y()<<" z "<<truth.vertex().Z());
         }  

         TFCSExtrapolationState result;
         m_FastCaloSimCaloExtrapolation->extrapolate(result,&truth);

         //write the result into the ntuple variables:

         ATH_MSG_DEBUG("IDCaloBoundary_eta() "<<result.IDCaloBoundary_eta());
         ATH_MSG_DEBUG("IDCaloBoundary_phi() "<<result.IDCaloBoundary_phi());
         ATH_MSG_DEBUG("IDCaloBoundary_r() "<<result.IDCaloBoundary_r());
         ATH_MSG_DEBUG("IDCaloBoundary_z() "<<result.IDCaloBoundary_z());
         ATH_MSG_DEBUG("AngleEta "<<result.IDCaloBoundary_AngleEta());
         ATH_MSG_DEBUG("Angle3D "<<result.IDCaloBoundary_Angle3D());

         m_newTTC_IDCaloBoundary_eta->push_back(float(result.IDCaloBoundary_eta()));
         m_newTTC_IDCaloBoundary_phi->push_back(float(result.IDCaloBoundary_phi()));
         m_newTTC_IDCaloBoundary_r->push_back(float(result.IDCaloBoundary_r()));
         m_newTTC_IDCaloBoundary_z->push_back(float(result.IDCaloBoundary_z()));
         m_newTTC_Angle3D ->push_back(float(result.IDCaloBoundary_Angle3D()));
         m_newTTC_AngleEta->push_back(float(result.IDCaloBoundary_AngleEta()));

         std::vector<float> eta_vec_ENT;
         std::vector<float> phi_vec_ENT;
         std::vector<float> r_vec_ENT;
         std::vector<float> z_vec_ENT;
         std::vector<float> detaBorder_vec_ENT;
         std::vector<bool>  OK_vec_ENT;

         std::vector<float> eta_vec_EXT;
         std::vector<float> phi_vec_EXT;
         std::vector<float> r_vec_EXT;
         std::vector<float> z_vec_EXT;
         std::vector<float> detaBorder_vec_EXT;
         std::vector<bool>  OK_vec_EXT;

         std::vector<float> eta_vec_MID;
         std::vector<float> phi_vec_MID;
         std::vector<float> r_vec_MID;
         std::vector<float> z_vec_MID;
         std::vector<float> detaBorder_vec_MID;
         std::vector<bool>  OK_vec_MID;

         for(int sample=CaloCell_ID_FCS::FirstSample;sample<CaloCell_ID_FCS::MaxSample;++sample) {
           ATH_MSG_DEBUG("sample "<<sample);
           ATH_MSG_DEBUG(" eta ENT "<<result.eta(sample,1)<<" eta EXT "<<result.eta(sample,2));
           ATH_MSG_DEBUG(" phi ENT "<<result.phi(sample,1)<<" phi EXT "<<result.phi(sample,2));
           ATH_MSG_DEBUG(" r   ENT "<<result.r(sample,1)  <<" r   EXT "<<result.r(sample,2)  );
           ATH_MSG_DEBUG(" z   ENT "<<result.z(sample,1)  <<" z   EXT "<<result.z(sample,2)  );
           ATH_MSG_DEBUG(" detaBorder   ENT "<<result.detaBorder(sample,1)  <<" detaBorder   EXT "<<result.detaBorder(sample,2)  );
           ATH_MSG_DEBUG(" OK  ENT "<<result.OK(sample,1) <<" OK  EXT "<<result.OK(sample,2)  );
           eta_vec_ENT.push_back(float(result.eta(sample,TFCSExtrapolationState::SUBPOS_ENT)));
           eta_vec_EXT.push_back(float(result.eta(sample,TFCSExtrapolationState::SUBPOS_EXT)));
           eta_vec_MID.push_back(float(result.eta(sample,TFCSExtrapolationState::SUBPOS_MID)));
           phi_vec_ENT.push_back(float(result.phi(sample,TFCSExtrapolationState::SUBPOS_ENT)));
           phi_vec_EXT.push_back(float(result.phi(sample,TFCSExtrapolationState::SUBPOS_EXT)));
           phi_vec_MID.push_back(float(result.phi(sample,TFCSExtrapolationState::SUBPOS_MID)));
           r_vec_ENT.push_back(float(result.r(sample,TFCSExtrapolationState::SUBPOS_ENT)));
           r_vec_EXT.push_back(float(result.r(sample,TFCSExtrapolationState::SUBPOS_EXT)));
           r_vec_MID.push_back(float(result.r(sample,TFCSExtrapolationState::SUBPOS_MID)));
           z_vec_ENT.push_back(float(result.z(sample,TFCSExtrapolationState::SUBPOS_ENT)));
           z_vec_EXT.push_back(float(result.z(sample,TFCSExtrapolationState::SUBPOS_EXT)));
           z_vec_MID.push_back(float(result.z(sample,TFCSExtrapolationState::SUBPOS_MID)));
           detaBorder_vec_ENT.push_back(float(result.detaBorder(sample,TFCSExtrapolationState::SUBPOS_ENT)));
           detaBorder_vec_EXT.push_back(float(result.detaBorder(sample,TFCSExtrapolationState::SUBPOS_EXT)));
           detaBorder_vec_MID.push_back(float(result.detaBorder(sample,TFCSExtrapolationState::SUBPOS_MID)));
           OK_vec_ENT.push_back(result.OK(sample,TFCSExtrapolationState::SUBPOS_ENT));
           OK_vec_EXT.push_back(result.OK(sample,TFCSExtrapolationState::SUBPOS_EXT));
           OK_vec_MID.push_back(result.OK(sample,TFCSExtrapolationState::SUBPOS_MID));
         }

         m_newTTC_back_eta->push_back(eta_vec_EXT);
         m_newTTC_back_phi->push_back(phi_vec_EXT);
         m_newTTC_back_r  ->push_back(r_vec_EXT);
         m_newTTC_back_z  ->push_back(z_vec_EXT);
         m_newTTC_back_detaBorder  ->push_back(detaBorder_vec_EXT);
         m_newTTC_back_OK  ->push_back(OK_vec_EXT);
         m_newTTC_entrance_eta->push_back(eta_vec_ENT);
         m_newTTC_entrance_phi->push_back(phi_vec_ENT);
         m_newTTC_entrance_r  ->push_back(r_vec_ENT);
         m_newTTC_entrance_z  ->push_back(z_vec_ENT);
         m_newTTC_entrance_detaBorder  ->push_back(detaBorder_vec_ENT);
         m_newTTC_entrance_OK  ->push_back(OK_vec_ENT);
         m_newTTC_mid_eta->push_back(eta_vec_MID);
         m_newTTC_mid_phi->push_back(phi_vec_MID);
         m_newTTC_mid_r  ->push_back(r_vec_MID);
         m_newTTC_mid_z  ->push_back(z_vec_MID);
         m_newTTC_mid_detaBorder  ->push_back(detaBorder_vec_MID);
         m_newTTC_mid_OK  ->push_back(OK_vec_MID);

         //*****************************

         //OLD EXTRAPOLATION
         /*
         std::vector<Trk::HitInfo>* hitVector = caloHits(*(*it));
         for(std::vector<Trk::HitInfo>::iterator it = hitVector->begin();it < hitVector->end();++it) {
           if((*it).trackParms) {
             delete (*it).trackParms;
             (*it).trackParms=0;
           }
         }
         delete hitVector;
         */

         //Amg::Vector3D mom((*it)->momentum().x(),(*it)->momentum().y(),(*it)->momentum().z());

         m_truth_energy->push_back((*it)->momentum().e());
         m_truth_px->push_back((*it)->momentum().px());
         m_truth_py->push_back((*it)->momentum().py());
         m_truth_pz->push_back((*it)->momentum().pz());
         m_truth_pdg->push_back((*it)->pdg_id());
         m_truth_barcode->push_back((*it)->barcode());

       } //for mcevent
     } //mcevent size
   } //mcEvent
1156
 }//truth event
1157

1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
// Get the reco clusters if available
// retreiving cluster container
  const DataHandle<xAOD::CaloClusterContainer > theClusters;
  std::string clusterContainerName = "CaloCalTopoClusters";
  sc = evtStore()->retrieve(theClusters, clusterContainerName);
  if (sc.isFailure()) {
    ATH_MSG_WARNING(" Couldn't get cluster container '" << clusterContainerName << "'");
    return 0;
  }
  xAOD::CaloClusterContainer::const_iterator itrClus = theClusters->begin();
  xAOD::CaloClusterContainer::const_iterator itrLastClus = theClusters->end();
  for ( ; itrClus!=itrLastClus; ++itrClus){
    const xAOD::CaloCluster *cluster =(*itrClus);
    m_cluster_energy->push_back(cluster->e());
    m_cluster_eta->push_back(cluster->eta());
    m_cluster_phi->push_back(cluster->phi());
    ATH_MSG_VERBOSE("Cluster energy: " << cluster->e() << " cells: " << " links: " << cluster->getCellLinks());
    //cluster->getCellLinks();
    const CaloClusterCellLink* cellLinks = cluster->getCellLinks();
    if (!cellLinks) {
      ATH_MSG_DEBUG( "No cell links for this cluster"  );
      continue;
    }

    const CaloCellContainer* cellCont=cellLinks->getCellContainer();
    if (!cellCont) {
      ATH_MSG_DEBUG( "DataLink to cell container is broken"  );
      continue;
    }
    unsigned cellcount = 0;
    std::vector<Long64_t> cellIDs_in_cluster;
    xAOD::CaloCluster::const_cell_iterator cellIter =cluster->cell_begin();
    xAOD::CaloCluster::const_cell_iterator cellIterEnd =cluster->cell_end();
    for ( ;cellIter !=cellIterEnd;cellIter++) {
      ++cellcount;
      const CaloCell* cell= (*cellIter);
      cellIDs_in_cluster.push_back(cell->ID().get_compact());
      float EnergyCell=cell->energy(); //ID, time, phi, eta
      ATH_MSG_DEBUG("   Cell energy: " << EnergyCell);
    }// end of cells inside cluster loop
    m_cluster_size->push_back(cellcount);
    m_cluster_cellID->push_back(cellIDs_in_cluster);
  }
1201
1202
1203


 //Retrieve and save MuonEntryLayer information 
1204
 const TrackRecordCollection *MuonEntry = nullptr;
1205
 sc = evtStore()->retrieve(MuonEntry, "MuonEntryLayer");
1206
1207
1208
1209
1210
 if (sc.isFailure())
 {
 ATH_MSG_WARNING( "Couldn't read MuonEntry from StoreGate");
 //return NULL;
 }
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
 else{
  for ( const TrackRecord &record : *MuonEntry){
    m_MuonEntryLayer_E->push_back((record).GetEnergy());
    m_MuonEntryLayer_px->push_back((record).GetMomentum().getX());
    m_MuonEntryLayer_py->push_back((record).GetMomentum().getY());
    m_MuonEntryLayer_pz->push_back((record).GetMomentum().getZ());
    m_MuonEntryLayer_x->push_back((record).GetPosition().getX());
    m_MuonEntryLayer_y->push_back((record).GetPosition().getY());
    m_MuonEntryLayer_z->push_back((record).GetPosition().getZ());
    m_MuonEntryLayer_pdg->push_back((record).GetPDGCode());
  }
1222
1223
 }

1224
1225
1226
 //Get reco cells if available
 const CaloCellContainer *cellColl = 0;
 sc = evtStore()->retrieve(cellColl, "AllCalo");
1227

1228
1229
1230
1231
1232
1233
1234
1235
1236
 if (sc.isFailure())
 {
  ATH_MSG_WARNING( "Couldn't read AllCalo cells from StoreGate");
  //return NULL;
 }
 else
 {
  ATH_MSG_INFO( "Found: "<<cellColl->size()<<" calorimeter cells");
  CaloCellContainer::const_iterator itrCell = cellColl->begin();
1237
  CaloCellContainer::const_iterator itrLastCell = cellColl->end();
1238
1239
  for ( ; itrCell!=itrLastCell; ++itrCell)
  {
1240
1241
         m_cell_energy->push_back((*itrCell)->energy());
         m_cell_identifier->push_back((*itrCell)->ID().get_compact());
1242
1243
1244
1245
1246
         if (m_tileID->is_tile_aux((*itrCell)->ID())) {
           // special case for E4'
           m_cell_sampling->push_back(CaloCell_ID::TileGap3);
         }
         else if (m_calo_dd_man->get_element((*itrCell)->ID()))
1247
         {
1248
          // all other Tile cells
1249
1250
1251
1252
1253
          CaloCell_ID::CaloSample layer = m_calo_dd_man->get_element((*itrCell)->ID())->getSampling();
          m_cell_sampling->push_back(layer);
         }
         else
          m_cell_sampling->push_back(-1);
1254
1255
  }
 } //calorimeter cells
1256

1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
 //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++)
   {
1269
1270
          hitnumber++;
          GeoLArHit ghit(**hi);
1271
1272
1273
    if (!ghit)
     continue;
    const CaloDetDescrElement *hitElement = ghit.getDetDescrElement();
1274
1275
          if(!hitElement)
           continue;
1276
    Identifier larhitid = hitElement->identify();
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
          if(m_calo_dd_man->get_element(larhitid))
          {
           CaloCell_ID::CaloSample larlayer = m_calo_dd_man->get_element(larhitid)->getSampling();

           float larsampfrac=m_dd_fSampl->FSAMPL(larhitid);
           m_g4hit_energy->push_back( ghit.Energy() );
           m_g4hit_time->push_back( ghit.Time() );
           m_g4hit_identifier->push_back( larhitid.get_compact() );
           m_g4hit_cellidentifier->push_back( larhitid.get_compact() );
           m_g4hit_sampling->push_back( larlayer);
           m_g4hit_samplingfraction->push_back( larsampfrac );
          }
1289
1290
1291
   } // End while LAr hits
   ATH_MSG_INFO( "Read "<<hitnumber<<" G4Hits from "<<lArKey[i]);
  }
1292
  else
1293
  {
1294
         ATH_MSG_INFO( "Can't retrieve LAr hits");
1295
1296
1297
  }// End statuscode success upon retrieval of hits
  //std::cout <<"ZH G4Hit size: "<<m_g4hit_e->size()<<std::endl;
 }// End detector type loop
1298

1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
 const TileHitVector * hitVec;
 if (evtStore()->retrieve(hitVec,"TileHitVec")==StatusCode::SUCCESS &&  m_tileMgr &&  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);
1309

1310
   if (m_calo_dd_man->get_element(cell_id))
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
         {
          CaloCell_ID::CaloSample layer = m_calo_dd_man->get_element(cell_id)->getSampling();

          float tilesampfrac = m_tileInfo->HitCalib(cell_id);

          //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;
           m_g4hit_energy->push_back( (*i_hit).energy(tilesubhit_i) );
           m_g4hit_time->push_back(   (*i_hit).time(tilesubhit_i)   );
           m_g4hit_identifier->push_back( pmt_id.get_compact() );
           m_g4hit_cellidentifier->push_back( cell_id.get_compact() );
           m_g4hit_sampling->push_back( layer );
           m_g4hit_samplingfraction->push_back( tilesampfrac );
          }
         }
1329
1330
1331
  }
  ATH_MSG_INFO( "Read "<<hitnumber<<" G4Hits from TileHitVec");
 }
1332

1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403