"MarlinReco"  1.16.0
SimDigital.h
1 #ifndef SimDigital_HHH
2 #define SimDigital_HHH
3 
4 #include "marlin/Processor.h"
5 #include "lcio.h"
6 #include <vector>
7 #include <string>
8 #include <map>
9 #include <set>
10 #include <utility>
11 #include <stdlib.h>
12 #include <EVENT/LCCollection.h>
13 #include <EVENT/SimCalorimeterHit.h>
14 #include <IMPL/CalorimeterHitImpl.h>
15 #include <IMPL/LCCollectionVec.h>
16 #include <IMPL/LCFlagImpl.h>
17 #include <UTIL/CellIDEncoder.h>
18 #include <UTIL/CellIDDecoder.h>
19 #include <marlin/Global.h>
20 #include <gear/GearParameters.h>
21 #include <gear/CalorimeterParameters.h>
22 #include <gear/LayerLayout.h>
23 #include <TROOT.h>
24 #include <TFile.h>
25 #include <TF2.h>
26 #include <TH1.h>
27 #include "TH1F.h"
28 
29 #include "CalorimeterHitType.h" //in MarlinUtil
30 #include "marlinutil/LCGeometryTypes.h"
31 
32 class TH1F;
33 class TF1;
34 class TTree;
35 
36 namespace AIDA {
37  class ITuple;
38 }
39 
40 using namespace lcio ;
41 using namespace marlin ;
42 
43 // Gerald Grenier :
44 // start removing the ECAL part whih has no use for SDHCAL
45 // this can be fully activated when there is an ECAL digitizer
46 // that don't do HCAL digitization
47 #define SIMDIGITAL_WITHECAL
48 
49 #ifdef SIMDIGITAL_WITHECAL
50 const int MAX_LAYERS = 200;
51 const int MAX_STAVES = 16;
52 #endif
53 
64 
66 {
67  StepAndCharge() : step(), charge(0) {}
68  StepAndCharge(LCVector3D vec) : step(vec), charge(0) {}
69  LCVector3D step;
70  float charge;
71 };
72 
73 
74 //helper class to manage cellId and local geometry
76 {
77  public:
78  static void bookTuples(const marlin::Processor* proc);
79  SimDigitalGeomCellId(LCCollection *inputCol, LCCollectionVec *outputCol);
81  //return the list of step positions in coordinates corresponding to 'I' ,'J' and 'layer'
82  std::vector<StepAndCharge> decode(SimCalorimeterHit *hit);
83  void encode(CalorimeterHitImpl *hit,int delta_I, int delta_J);
84  void setLayerLayout( CHT::Layout layout);
85  float getCellSize();
86  const LCVector3D& normalToRPCPlane() {return _normal;}
87  const LCVector3D& Iaxis() {return _Iaxis;}
88  const LCVector3D& Jaxis() {return _Jaxis;}
89 
90  int I() {return _Iy;}
91  int J() {return _Jz;}
92  int K() {return _trueLayer;}
93  int stave() {return _stave;}
94  int module() {return _module;}
95 
96  private:
97  enum HCAL_GEOM {VIDEAU,TESLA};
98  HCAL_GEOM _geom;
99  int _trueLayer;
100  int _stave;
101  int _module;
102  int _Iy;
103  int _Jz;
104  const float* _hitPosition;
105  CellIDDecoder<SimCalorimeterHit> _decoder;
106  CellIDEncoder<CalorimeterHitImpl> _encoder;
107  const gear::LayerLayout* _layerLayout;
108  SimDigitalGeomRPCFrame* _normal_I_J_setter;
109  CHT::Layout _currentHCALCollectionCaloLayout;
110  LCVector3D _normal;
111  LCVector3D _Iaxis;
112  LCVector3D _Jaxis;
113  static AIDA::ITuple* _tupleHit;
114  enum {TH_DETECTOR,TH_CHTLAYOUT,TH_MODULE,TH_STAVE,TH_LAYER,TH_I,TH_J,
115  TH_X,TH_Y,TH_Z,
116  TH_NORMALX,TH_NORMALY,TH_NORMALZ,
117  TH_IX,TH_IY,TH_IZ,
118  TH_JX,TH_JY,TH_JZ};
119  static AIDA::ITuple* _tupleStep;
120  enum {TS_DETECTOR,TS_CHTLAYOUT,TS_HITCELLID,TS_NSTEP,
121  TS_HITX,TS_HITY,TS_HITZ,
122  TS_STEPX,TS_STEPY,TS_STEPZ,
123  TS_DELTAI,TS_DELTAJ,TS_DELTALAYER};
124  friend class SimDigitalGeomRPCFrame;
125 };
126 
127 
128 //hierarchy of classes to determine the RPC reference frame
130 {
131  public:
132  SimDigitalGeomRPCFrame(SimDigitalGeomCellId& h) : _layerInfo(h) {;}
133  virtual ~SimDigitalGeomRPCFrame() {}
134  virtual void setRPCFrame()=0;
135  private:
136  SimDigitalGeomCellId& _layerInfo;
137  protected:
138  int stave() {return _layerInfo._stave;}
139  int module() {return _layerInfo._module;}
140  LCVector3D& normal() {return _layerInfo._normal;}
141  LCVector3D& Iaxis() {return _layerInfo._Iaxis;}
142  LCVector3D& Jaxis() {return _layerInfo._Jaxis;}
143 };
144 
146 {
147  public:
149  void setRPCFrame();
150 };
152 {
153  public:
155  void setRPCFrame();
156 };
158 {
159  public:
161  void setRPCFrame();
162 };
164 {
165  public:
167  void setRPCFrame();
168 };
169 
170 
171 
172 //helper class to manage multiplicity
174 {
175  public:
177  virtual ~multiplicityChargeSplitterBase() {}
178  typedef std::pair<int,int> I_J_Coordinates;
179  virtual void addCharge(float charge, float pos_I, float pos_J)=0;
180  void newHit(float cell_Size) {_chargeMap.clear();_cellSize=cell_Size;}
181  const std::map<I_J_Coordinates,float>& chargeMap() {return _chargeMap;}
182  protected:
183  float _cellSize;
184  std::map<I_J_Coordinates,float> _chargeMap;
185 };
186 
187 
189 {
190  public:
192  void addCharge(float charge, float pos_I, float pos_J);
193  private:
194  float _edgeSize;
195  friend class SimDigital;
196 };
197 
198 
200 {
201  public:
204  void init();
205  void addCharge(float charge, float pos_I, float pos_J);
206  private:
207  TF2* _f2;
208  float _range;
209  std::string _function_description;
210  std::vector<float> _functionParameters;
211  float _normalisation;
212  float _RPC_PadSeparation; //distance between cells
213  friend class SimDigital;
214 };
215 
217 {
218  public:
221  void init();
222  void addCharge(float charge, float pos_I, float pos_J);
223  private:
224  float _range;
225  std::vector<float> _erfWidth;
226  std::vector<float> _erfWeigth;
227  float _normalisation;
228  float _RPC_PadSeparation; //distance between cells
229  friend class SimDigital;
230 };
231 
232 
233 //helper class to manage efficiency maps
235 {
236  public:
237  virtual float getEfficiency(int I, int J, int K, int stave, int module)=0;
238 };
239 
241 {
242  public:
243  effMapConstant(float val=1.0) : value(val) {}
244  virtual float getEfficiency(int I, int J, int K, int stave, int module) { return value;}
245  private:
246  float value;
247 };
248 
249 
251 {
252  public:
253  effMapProtoByAsic(std::string fileName);
254  virtual float getEfficiency(int I, int J, int K, int stave, int module)
255  {
256  std::map<int,float>::iterator it=_effMap.find( (I-1)/8+((J-1)/8)*12+K*1000 );
257  return (it != _effMap.end() ? it->second : 1.0);
258  }
259  private:
260  std::map<int,float> _effMap;
261 };
262 
263 
264 class SimDigital : public Processor {
265  public:
266  virtual Processor* newProcessor() { return new SimDigital;}
267  SimDigital();
268 
272  virtual void init() ;
273 
276  virtual void processRunHeader( LCRunHeader* run ) ;
277 
280  virtual void processEvent( LCEvent * evt ) ;
281 
282  virtual void check( LCEvent * evt ) ;
283 
286  virtual void end() ;
287 
288 #ifdef SIMDIGITAL_WITHECAL
289  virtual void fillECALGaps() ;
290 #endif
291 
292  private:
293 #ifdef SIMDIGITAL_WITHECAL
294  std::vector<std::string> _ecalCollections;
295  std::string _outputEcalCollection0;
296  std::string _outputEcalCollection1;
297  std::string _outputEcalCollection2;
298  std::vector<std::string> _outputEcalCollections;
299  std::vector<float> _calibrCoeffEcal;
300  std::vector<int> _ecalLayers;
301  int _digitalEcal;
302  float _thresholdEcal;
303  int _ecalGapCorrection;
304  float _ecalGapCorrectionFactor;
305  float _ecalModuleGapCorrectionFactor;
306  float _ecalEndcapCorrectionFactor;
307  std::vector<CalorimeterHitImpl*> _calHitsByStaveLayer[MAX_STAVES][MAX_LAYERS];
308  std::vector<int> _calHitsByStaveLayerModule[MAX_STAVES][MAX_LAYERS];
309 
310  float _zOfEcalEndcap;
311  float _barrelPixelSizeT[MAX_LAYERS];
312  float _barrelPixelSizeZ[MAX_LAYERS];
313  float _endcapPixelSizeX[MAX_LAYERS];
314  float _endcapPixelSizeY[MAX_LAYERS];
315  float _barrelStaveDir[MAX_STAVES][2];
316 #endif
317  //std::vector<int> _hcalLayers;
318  //std::vector<std::string> _calorimeterHitCollections;
319  std::vector<std::string> _hcalCollections;
320  std::vector<std::string> _outputHcalCollections;
321  std::map<std::string, int> _counters;
322  std::vector<float> _thresholdHcal;
323  std::vector<float>_calibrCoeffHcal;
324  std::string _outputRelCollection;
325  bool _printSimDigital;
326  TF1 * _QPolya;
327  double _polyaAverageCharge, _polyaFunctionWidthParameter;
328 
329  LCCollectionVec* _relcol;
330 #ifdef SIMDIGITAL_WITHECAL
331  void registerECALparameters();
332  void setECALgeometry();
333  void processECAL(LCEvent* evt, LCFlagImpl& flag);
334 #endif
335  void processHCAL(LCEvent* evt, LCFlagImpl& flag);
336 
337  static bool sortStepWithCharge(StepAndCharge s1, StepAndCharge s2){return s1.charge>s2.charge;}
338 
339  //intermediate storage class
340  struct hitMemory
341  {
342  hitMemory() : ahit(0),relatedHits(), maxEnergydueToHit(-1), rawHit(-1) {}
343  CalorimeterHitImpl *ahit;
344  std::set<int> relatedHits;
345  float maxEnergydueToHit;
346  int rawHit;
347  };
348  float depositedEnergyInRPC;
349  multiplicityChargeSplitterUniform _chargeSplitterUniform;
350  multiplicityChargeSplitterFunction _chargeSplitterFunction;
351  multiplicityChargeSplitterErfFunction _chargeSplitterErfFunction;
352  multiplicityChargeSplitterBase* _theChosenSplitter;
353  bool _doThresholds;
354  std::string _chargeSplitterOption;
355  std::string _effMapFileName;
356  float _constEffMapValue;
357  std::string _effMapOption;
358  effMapBase *_effMap;
359  float _absZstepFilter;
360  bool _keepAtLeastOneStep;
361  float _minXYdistanceBetweenStep;
362  AIDA::ITuple* _debugTupleStepFilter;
363  AIDA::ITuple* _tupleStepFilter;
364  AIDA::ITuple* _tupleCollection;
365  multiplicityChargeSplitterBase& getSplitter() { return *_theChosenSplitter; }
366 
367  //predicate class to remove potential hit below threshold
368  class ThresholdIsBelow
369  {
370  float value;
371  public:
372  ThresholdIsBelow(float f) : value(f) {;}
373  bool operator()(std::pair<int,hitMemory> f) { return f.second.ahit->getEnergy()<value;}
374  };
375 
376  //predicate class to remove steps
377  class absZGreaterThan
378  {
379  public:
380  absZGreaterThan(float val) : _value(val) {}
381  bool operator()(StepAndCharge& v) { return fabs( v.step.z() ) >_value;}
382  private:
383  float _value;
384  };
385  class randomGreater
386  {
387  public:
388  randomGreater(float val) : _value(val) {}
389  bool operator()(StepAndCharge& v) { return double(rand())/RAND_MAX>_value;}
390  private:
391  float _value;
392  };
393  //helper function to remove steps too close in I,J
394  void remove_adjacent_step(std::vector<StepAndCharge>& vec);
395  void fillTupleStep(std::vector<StepAndCharge>& vec,int level);
396 
397  LCCollectionVec * processHCALCollection(LCCollection * col ,CHT::Layout layout, LCFlagImpl& flag);
398  void createPotentialOutputHits(std::map<int,hitMemory>& myHitMap, LCCollection *col, SimDigitalGeomCellId& aGeomCellId);
399  void removeHitsBelowThreshold(std::map<int,hitMemory>& myHitMap, float threshold);
400  void applyThresholds(std::map<int,hitMemory>& myHitMap);
401 
402 };
403 
404 #endif
Definition: SimDigital.h:163
Definition: SimDigital.h:173
Definition: SimDigital.h:151
Definition: SimDigital.h:199
Definition: SimDigital.h:75
Definition: SimDigital.h:65
Definition: SimDigital.h:216
Definition: SimDigital.h:129
Definition: SimDigital.h:250
Definition: SimDigital.h:234
Definition: SimDigital.h:264
Definition: SimDigital.h:240
Definition: SimDigital.h:157
Definition: SimDigital.h:145
Definition: SimDigital.h:188