IB-0.2
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
IBAnalyzerEM.h
Go to the documentation of this file.
1 /*////////////////////////////////////////////////////////////////////////////
2  Copyright 2018 Istituto Nazionale di Fisica Nucleare
3 
4  Licensed under the EUPL, Version 1.2 or - as soon they will be approved by
5  the European Commission - subsequent versions of the EUPL (the "Licence").
6  You may not use this work except in compliance with the Licence.
7 
8  You may obtain a copy of the Licence at:
9 
10  https://joinup.ec.europa.eu/software/page/eupl
11 
12  Unless required by applicable law or agreed to in writing, software
13  distributed under the Licence is distributed on an "AS IS" basis, WITHOUT
14  WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
15  Licence for the specific language governing permissions and limitations under
16  the Licence.
18 
19 
20 #ifndef IBANALYZEREM_H
21 #define IBANALYZEREM_H
22 
23 #include "IBAnalyzer.h"
24 #include "IBVoxCollection.h"
25 #include "IBVoxRaytracer.h"
26 #include "IBVoxel.h"
27 #include <string>
28 #include <iomanip>
29 
30 class IBPocaEvaluator;
31 class IBMinimizationVariablesEvaluator;
32 class IBAnalyzerEMAlgorithm;
33 
34 
35 class IBAnalyzerEM : public IBAnalyzer
36 {
37 public:
38  typedef IBAnalyzer BaseClass;
39 
40  struct Event {
41  struct Element {
42  Matrix2f Wij;
43  union {
44  Scalarf lambda;
45  Scalarf Sij;
46  };
47  IBVoxel *voxel;
48  Scalarf pw;
49  };
50 
51  struct {
52  Vector4f Di;
53  Matrix4f E;
54  Scalarf InitialSqrP;
55  Scalarf pTrue; // SV 20160921 variable to store useful quantities to study....
56  } header;
57  std::vector<Element> elements;
58  };
59 
60 
61 public:
62  IBAnalyzerEM(IBVoxCollection &voxels, int nPath=2, double alpha=0., bool doRecoPath=true,
63  bool oldTCalculation=false, float rankLimit=-100., IBVoxCollection* initialSqrPfromVtk=NULL, int pVoxelMean=0);
64  ~IBAnalyzerEM();
65 
66  inline virtual const char *type_name() const { return "IBAnalyzerEM"; }
67 
68  bool AddMuon(const MuonScatterData &muon);//{ return false;}
69  bool AddMuonFullPath(const MuonScatterData &muon, std::vector<HPoint3f>& muonPath);
70 
71  void SetMuonCollection(IBMuonCollection *muons);
72 
73  unsigned int Size();
74 
75  void Run(unsigned int iterations, float muons_ratio);
76 
77  void SetMLAlgorithm(IBAnalyzerEMAlgorithm *MLAlgorithm);
78 
79  inline IBPocaEvaluator *GetPocaAlgorithm() const { return this->m_PocaAlgorithm; }
80  inline void SetPocaAlgorithm(IBPocaEvaluator *algo) { this->m_PocaAlgorithm = algo; }
81 
82  inline IBMinimizationVariablesEvaluator *GetVarAlgorithm() const {
83  return this->m_VarAlgorithm;
84  }
85  inline void SetVarAlgorithm(IBMinimizationVariablesEvaluator *eval) {
86  this->m_VarAlgorithm = eval;
87  }
88 
89  inline IBVoxRaytracer *GetRayAlgorithm() const { return this->m_RayAlgorithm; }
90  inline void SetRayAlgorithm(IBVoxRaytracer *algo) { this->m_RayAlgorithm = algo; }
91 
92  inline IBAbstract::IBVoxCollectionUpdateAlgorithm *GetUpdateAlgorithm() const {
93  return this->m_UpdateAlgorithm;
94  }
95  inline void SetUpdateAlgorithm(IBAbstract::IBVoxCollectionUpdateAlgorithm *algo) {
96  this->m_UpdateAlgorithm = algo;
97  }
98 
99  void filterEventsVoxelMask();
100 
101  void filterEventsLineDistance(float min, float max);
102 
103  void SijCut(float threshold);
104 
105  std::vector<Event > SijCutCount(float threshold_low, float threshold_high);
106 
107  void dumpEventsSijInfo(const char *filename, std::vector<float> N);
108 
109  void SijGuess(std::vector<Vector2f> tpv);
110 
111  void Chi2Cut(float threshold);
112 
113  void SetVoxCollection(IBVoxCollection *voxels);
114 
115  void SetVoxcollectionShift(Vector3f shift);
116 
117  void dumpEventsTTree(const char *filename);
118  void DumpP(const char *filename, float x0 = 0, float x1 = 10);
119  void DumpEvent(Event *evc);
120 
121  std::vector<Event> & Events();
122 
123  float SijMedian(const Event &evc);
124 
125  void SetSijMedianMomentum();
126 
127 protected:
128 
129  IBAnalyzerEMAlgorithm *m_SijAlgorithm;
130  Scalarf nominal_momentum;
131  //Scalarf SijCutEM;
132 
133 private:
134 
135  void Project(Event *evc);
136  void BackProject(Event *evc);
137  void Evaluate(float muons_ratio);
138 
139  IBPocaEvaluator *m_PocaAlgorithm;
140  IBMinimizationVariablesEvaluator *m_VarAlgorithm;
141  IBVoxRaytracer *m_RayAlgorithm;
142  IBAbstract::IBVoxCollectionUpdateAlgorithm *m_UpdateAlgorithm;
143 
144  int m_nPath; //---- Int to indicate whether to build a 1, 2 or 3-path
145  double m_alpha; //---- Relative distance along the trajectory to build the 3-path
146  bool m_useRecoPath; //---- Use the true muon path from MC
147  bool m_oldTCalculation; //---- Use the old method of calculating length parameter T
148  float m_rankLimit;
149 
150  IBVoxCollection* m_initialSqrPfromVtk;
151  int m_pVoxelMean; //---- compute p voxel by hand
152  IBVoxCollection m_imgMC;
153 
154  std::vector<IBAnalyzerEM::Event> m_Events;
155  bool m_firstIteration;
156 
157 };
158 
159 
160 #endif // IBANALYZEREM_H