MarlinTPC  1.2.0
MMHitFinderProcessor.h
1 #ifndef MMHITFINDERPROCESSOR_H
2 #define MMHITFINDERPROCESSOR_H
3 
4 // MARLIN
5 #include "marlin/Processor.h"
6 #include "marlin/Global.h"
7 
8 // LCIO
9 #include "lcio.h"
10 #include "EVENT/LCCollection.h"
11 #include "IMPL/LCCollectionVec.h"
12 #include "EVENT/TrackerPulse.h"
13 
14 // GEAR
15 #include "gear/TPCParameters.h"
16 
17 // C++
18 #include <string>
19 #include <vector>
20 #include <map>
21 #include <set>
22 #include <deque>
23 
24 // MarlinTPC
25 #include "ChannelCorrectionListener.h"
26 #include "TPCConditionsListener.h"
27 
28 // Aida
29 #include <AIDA/AIDA.h>
30 #include <marlin/AIDAProcessor.h>
31 #include <AIDA/IHistogramFactory.h>
32 #include <AIDA/IHistogram1D.h>
33 
34 // ROOT
35 #include <Minuit2/FCNBase.h>
36 #include <TMath.h>
37 
38 // Minuit2 headers used to find seed track.
39 #include <Minuit2/MnUserParameters.h>
40 #include <Minuit2/MnMigrad.h>
41 #include <Minuit2/FunctionMinimum.h>
42 #include <Minuit2/Minuit2Minimizer.h>
43 
44 #include "Math/Minimizer.h"
45 #include "Math/Factory.h"
46 #include "Math/Functor.h"
47 
48 namespace marlintpc
49 {
50 
51  /* Processor to find Hits from Micromegas data.
52  *
53  * Basic procedure: Sort pulses by module and row, search rows for groups of associated pulses,
54  * then either fit pulses with the Pad Response Function, or use the weighted mean if calibrating.
55  *
56  * Parameters:
57  *
58  * InputTrackerPulses = Input pulses [string]
59  * OutputTrackerHits = Output hits [string]
60  * OutputTrackerHitPulses = Output collection of pulses associated with output hits [string]
61  * WriteOutputToStorage = Set persistency of output collections [string]
62  * MaxEmptyPads = Maximum number of empty pads in a hit [integer]
63  * MinHitSizeOverride = Minimum number of pulses in a hit [integer]
64  * VDriftOverride = Drift velocity (not really important at this stage, but will be later) [float]
65  * ModuleList = A vector of modules to be used in the analysis [IntVec]
66  * MaxTimeSpreadOverride = Maximum time between pulses in a hit [float]
67  * PRFType = Name of the PRF to be used [string]
68  * PRF_Parameters = PRF parameters to be used (can be empty only if it is a calibration run) [FloatVec]
69  * RowDep = Set whether or not the PRF parameters have a row dependency [boolean]
70  * DeadRowList = A list of rows not to be used in the analysis. Format: MODULE ROW (e.g. 0 1 0 2 ... etc.) [IntVec]
71  * CentreColumns = A list of central beam columns (same format as for the rows). Only data in columns +/- 5 of the central column is used [IntVec]
72  * Threshold = Minimum amplitude of the maximum pulse in a hit [integer]
73  * CalibRun = IMPORTANT! Sets whether or not this is a calibration run. Uses the weighted mean instead of the PRF, and related settings [boolean]
74  * ExcludeRows = Whether or not to exclude rows in DeadRowList [boolean]
75  * SelectCentralCols = Whether or not to use the central 10 columns indicated in CentreColumns [boolean]
76  * ModuleDependence = Set wheter or not the PRF parameters have a module dependency [boolean]
77  * PadHits = Whether or not to pad hits with empty pulses to ensure a minimum number of pulses per hit [boolean]
78  *
79  * Author: P. Hayman, Carleton University
80  */
81 
82  class MMHitFinderProcessor : public marlin::Processor
83  {
84 
85  public:
86 
87  virtual Processor* newProcessor(){
88  return new MMHitFinderProcessor;
89  }
90 
92 
93  virtual void init() ;
94  virtual void processRunHeader(EVENT::LCRunHeader* run) ;
95  virtual void processEvent(EVENT::LCEvent * evt) ;
96  virtual void check(EVENT::LCEvent * evt) ;
97  virtual void end() ;
98 
100 
101  //Static function to calculate the row dependency of PRF parameters, if desired. This function can be used from anywhere that needs this information.
102  static double calcRowDep(const double * x, int index, int row)
103  {
104  //Currently a quadratic form: par[i] = a[i]*row^2 + b[i]*row + c[i]
105 
106  return (double)(row*row)*x[3*index] + (double)(row)*x[3*index + 1] + x[3*index + 2];
107  }
108 
109  //Also accessible from anywhere, this function will fill an array compatible with the MMPadResponseFunction, taking into account both row and module
110  //dependencies, if desired.
111  //static const double * calcParArray(const double * x, int nPar, int nMod, int modIndex, int rowIndex)
112  static const std::vector<double> calcParArray(const double * x, int nPar, int nMod, int modIndex, int rowIndex)
113  {
114  //double * retVal = new double[nPar];
115  std::vector<double> retVal(nPar);
116  if(modIndex >= 0)
117  {
118  if(rowIndex >= 0)
119  {
120  for(int i = 0; i < nPar; i++)
121  {
122  retVal[i] = calcRowDep(x, i*nMod + modIndex, rowIndex);
123  }
124  }
125  else
126  {
127  for(int i = 0; i < nPar; i++)
128  {
129  retVal[i] = x[i*nMod + modIndex];
130  }
131  }
132  }
133  else if(rowIndex >= 0)
134  {
135  for(int i = 0; i < nPar; i++)
136  {
137  retVal[i] = calcRowDep(x, i, rowIndex);
138  }
139  }
140  else
141  {
142  //return x;
143  for(int i = 0; i < nPar; i++)
144  retVal[i] = x[i];
145  }
146 
147  return retVal;
148  }
150 
151  // a helper class
153  {
154  public:
155 
156  //Assume the first pulse is the largest. This could be made more dynamic, but hopefully that's not necessary.
157  HitCandidate(EVENT::TrackerPulse* max)
158  {
159  _thePulses.push_back(max);
160  _maxPulse = max;
161  qualityFlag = 0;
162  };
163 
164  void addPulse(const EVENT::TrackerPulse* aPulse, bool left)
165  {
166  if(left)
167  _thePulses.push_front(aPulse);
168  else
169  _thePulses.push_back(aPulse);
170  };
171 
172  const EVENT::TrackerPulse* getMaximumPulse() const{
173  return _maxPulse;
174  };
175 
176  double getMaxCharge() const { return _maxPulse->getCharge(); };
177 
178  int getNumberOfPulses() const { return _thePulses.size(); };
179 
180  const std::deque<const EVENT::TrackerPulse*>& getPulses() const{
181  return _thePulses;
182  };
183 
184  int getQuality() const { return qualityFlag; };
185  void setQuality(int q) { qualityFlag = q; };
186 
187  private:
188 
189  int qualityFlag;
190  std::deque<const EVENT::TrackerPulse*> _thePulses; //Use a double-ended queue. First, the max pulse is added.
191  //Then, the pulses to the left of it are added to the front,
192  //and the pulses to the right of it are added to the end.
193  const EVENT::TrackerPulse* _maxPulse;
194  };
195 
196  private:
197 
198  //Processor Parameters
199  std::string _inputTrackerPulsesCollectionName ;
200  std::string _outputTrackerHitsCollectionName ;
201  std::string _outputTrackerHitPulsesCollectionName ;
202  std::string _inputTPCConditionsCollectionName;
203 
204  // A vector containing the modules which should be used in the analysis
205  std::vector<int> _moduleListParam;
206 
207  bool _outputHitsPersistent;
208  bool _rowDep;
209 
210  int _maxEmptyChannelsOverride;
211  int _maxDeadChannelsOverride;
212  int _minHitSizeOverride;
213  float _minMaxPulseSliceHeightOverride;
214  float _maxTimeSpreadOverride;
215  float _driftVelocity;
216 
217  std::string _rootFileName;
218  std::string _prfType;
219 
220  std::vector<float> _parVec;
221 
222 
223  //Optional Parameters
224  std::vector<int> _centreColListParam;
225  std::vector<int> _deadRowListParam;
226 
227  int _threshold;
228 
229  bool _excludeRows;
230  bool _selectCentralCols;
231  bool _calibrationRun;
232  bool _moduleDep;
233  bool _padHits;
234 
235 
236  //Class Variables
237  int _evt;
238  int _nPar;
239 
240  double * _parArray;
241 
242  std::set<int> _moduleList;
243  std::map<int, std::set<int> > _deadRowList;
244  std::map<int, int > _centreColList;
245  std::map<int, int > _moduleMap;
246 
247  ROOT::Minuit2::Minuit2Minimizer* min;
248 
249  // Some AIDA histograms used for error checking.
250  std::map<int, AIDA::IHistogram1D *> numberOfPulsesPerHit;
251  std::map<int, AIDA::IHistogram1D *> averagePulsesPerHitPerRow;
252  std::map<int, AIDA::IHistogram1D *> numberOfHitsPerEvent;
253  std::map<int, AIDA::IHistogram1D *> chiSquareDistrubtionOfHits;
254  std::map<int, AIDA::IHistogram1D *> hitFitResults;
255  std::map<int, AIDA::IHistogram1D *> hitDistribution;
256  std::map<int, AIDA::IHistogram1D *> hitErrorDistribution;
257  std::map<int, AIDA::IHistogram1D *> hitShapes;
258 
259  //Vectors used for calculating the average number of pulse in a hit per row.
260  std::map<int, std::vector<int> > numberOfHitsInRow;
261  std::map<int, std::vector<int> > numberOfPulsesInRow;
262 
263 
264  //Typedefs
265 
266  //Can't declare this before HitCandidate, for some reason...
267  typedef std::map<int, std::vector<MMHitFinderProcessor::HitCandidate> > SortedHitMap;
268  typedef std::map<int, std::map<int, std::vector<TrackerPulse*> > > SortedPulseMap;
269 
270 
271  /************************************************************************************************/
272  // helper functions:
273 
274  // Sorts pulses into collection corresponding to modules and rows.
275  // i.e a pulse on row 5 of module 1 is put in the same container as the other pulses
276  // from that row and module.
277  void _sortPulses ( EVENT::LCCollection*, const gear::TPCParameters&, SortedPulseMap& );
278 
279  // Collects the pulses in a row into a potential hit. This is the first stage of hit determination.
280  // As long as the pulses are beside each other spatially and close in time they are a potential hit.
281  // Hence the name hit CANDIDATE.
282  void _findHitCandidates ( SortedPulseMap&, const gear::TPCParameters&, SortedHitMap& );
283 
284  // Does all the neseccary checks and calculations to turn a HitCandidate into a full fledged hit!
285  // It then stores this TPCHit in the output hit collection which is saved to the lcio file.
286  void _createHitCollection ( SortedHitMap&, const gear::TPCParameters&, IMPL::LCCollectionVec* );
287 
288  // All these functions, as the names suggest, calculate the position and position error of the hit in 3D.
289  // The "fixed" coordinate is determined by the row position.
290  // The "free" coordinate is determined by fitting the PRF to the pulse amplitudes in the row.
291  // The "Z" coordinate is determined from the time measurement of the main pulse in the hit.
292  // Each has a return value of a pair: first = value, second = error.
293  std::pair<double, double> _calculateFixedCoordinate ( const MMHitFinderProcessor::HitCandidate&, const gear::TPCParameters& );
294  std::pair<double, double> _calculateFreeCoordinate ( const MMHitFinderProcessor::HitCandidate&, const gear::TPCParameters& );
295  std::pair<double, double> _calculateFreeCoordinatePRF ( const MMHitFinderProcessor::HitCandidate&, const gear::TPCParameters& );
296  std::pair<double, double> _calculateFreeCoordinateWeightedMean ( const MMHitFinderProcessor::HitCandidate&, const gear::TPCParameters& );
297  std::pair<double, double> _calculateZCoordinate ( const MMHitFinderProcessor::HitCandidate& );
298 
299  //Useful helper method for finding the largest pulse in a row
300  EVENT::TrackerPulse* _findHighestPulseInRow ( std::vector<EVENT::TrackerPulse*>& );
301  //Useful helper method for finding pulses that are near a max pulse (i.e., forming a HitCandidate)
302  void _addNeighbourPulses ( MMHitFinderProcessor::HitCandidate&, const gear::TPCParameters&,
303  std::vector<EVENT::TrackerPulse*>& );
304 
305  };//End MMHitFinderProcessor class
306 
307 }//End namespace marlintpc
308 
309 #endif
Definition: MMHitFinderProcessor.h:82
Definition: MMHitFinderProcessor.h:152
virtual void processEvent(EVENT::LCEvent *evt)
Definition: MMHitFinderProcessor.cc:319