StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StFtpcClusterFinder.cc
1 // $Id: StFtpcClusterFinder.cc,v 1.83 2015/07/20 17:37:01 jeromel Exp $
2 //
3 // $Log: StFtpcClusterFinder.cc,v $
4 // Revision 1.83 2015/07/20 17:37:01 jeromel
5 // Use std::isnan to satisfy C++11
6 //
7 // Revision 1.82 2011/09/09 21:32:42 jcs
8 // Undo rev 1.81 since different cluster finding parameters are used for laser and t0 calibration
9 // Comment out the different parameters when studying cluster finding and/or tracking
10 //
11 // Revision 1.81 2011/09/09 12:53:19 jcs
12 // use the same cluster finding parameters for debug (fdbg) that are used for production
13 //
14 // Revision 1.80 2011/05/21 15:33:30 jcs
15 // make corrections to save correct ChargeSum for clusters (one line of code
16 // was deleted inadvertently with update 1.79 and the wrong conversion between
17 // integer and float was performed when splitting merged clusters)
18 //
19 // Revision 1.79 2010/04/08 16:46:16 jcs
20 // swap data for RDO6,RDO7 FTPC East when Calibrations_ftpc/ftpcElectronics->swapRDO6RDO7East=1
21 //
22 // Revision 1.78 2009/12/11 15:41:06 jcs
23 // For laser run use laserTZero and no inner cathode correction
24 //
25 // Revision 1.77 2009/11/25 19:50:15 jcs
26 // remove all references to StFtpcSoftwareMonitor
27 //
28 // Revision 1.76 2009/11/14 12:51:08 jcs
29 // added suggested parentheses to avoid warnings which appeared with system upgrade
30 //
31 // Revision 1.75 2009/10/14 15:52:43 jcs
32 // write out all gas temperature, air pressure info to Run branch of FTPC debug root file
33 //
34 // Revision 1.74 2008/07/14 19:47:02 jcs
35 // The inner cathode correction is an Ftpc geometry correction; it should be applied to both data and laser runs
36 //
37 // Revision 1.73 2007/02/01 11:57:00 jcs
38 // move unessential output from INFO to DEBUG
39 //
40 // Revision 1.72 2007/01/15 07:49:21 jcs
41 // replace printf, cout and gMesMgr with Logger
42 //
43 // Revision 1.71 2006/05/04 06:18:48 jcs
44 // remove old debug statement
45 //
46 // Revision 1.70 2006/03/19 19:29:44 jcs
47 // Move cluster struct definitions to StFtpcClustersStructures.hh
48 // Create DEBUGFILE with bfc option 'fdbg'
49 //
50 // Revision 1.69 2005/12/02 09:03:10 jcs
51 // Delete pradius,pdeflection before error exits to avoid memory leaks
52 //
53 // Revision 1.68 2005/07/12 15:18:01 jcs
54 // save raw (not corrected) timeposition
55 //
56 // Revision 1.67 2005/07/11 13:33:03 jcs
57 // remove duplicate code
58 //
59 // Revision 1.66 2005/06/30 09:19:06 jcs
60 // remove commented out values which were moved to Geometry_ftpc/ftpcClusterGeom a long time ago
61 //
62 // Revision 1.65 2004/12/10 23:07:46 jcs
63 // Only fill FTPC software monitor if it exists
64 //
65 // Revision 1.64 2004/09/07 14:08:17 jcs
66 // use the IAttr(".histos") to control histogramming
67 // remove obsolete clfradius histogram
68 //
69 // Revision 1.63 2004/08/09 12:50:14 jcs
70 // move initialization or westHits and eastHits outside of loop over iftpc
71 //
72 // Revision 1.62 2004/08/03 00:48:32 oldi
73 // Printout of found clusters in FTPC east and west fixed.
74 //
75 // Revision 1.61 2004/06/18 12:04:57 jcs
76 // replace #ifdef...#elif...#endif conditional compiler directives with #ifdef...#endif #ifdef...#endif
77 //
78 // Revision 1.60 2004/06/18 09:04:39 jcs
79 // replace obsolete DEBUGFILE code with code to write out a root file for cluster/laser analysis
80 //
81 // Revision 1.59 2004/05/24 13:38:47 jcs
82 // save number of clusters found in StFtpcSoftwareMonitor
83 //
84 // Revision 1.58 2004/05/19 09:53:38 jcs
85 // print out number of clusters found for both Ftpcsy
86 //
87 // Revision 1.57 2004/04/06 18:36:51 oldi
88 // New data mebers for pad and time position and pad and time sigma filled.
89 //
90 // Revision 1.56 2004/01/28 02:04:43 jcs
91 // replace all instances of StFtpcReducedPoint and StFtpcPoint with StFtpcConfMapPoint
92 //
93 // Revision 1.55 2004/01/28 01:41:13 jeromel
94 // Change OST to OS everywhere since defaultoption is now not to print
95 // the date.
96 //
97 // Revision 1.54 2004/01/02 10:27:38 jcs
98 // remove incorrect correction to lower limit in fastlog filling loop
99 //
100 // Revision 1.53 2003/12/09 10:48:57 jcs
101 // correct lower limit in fastlog filling loop
102 //
103 // Revision 1.52 2003/11/27 03:13:13 perev
104 // protection sqrt(-1) added
105 //
106 // Revision 1.51 2003/10/08 13:49:53 jcs
107 // initialize pointers and arrays
108 //
109 // Revision 1.50 2003/09/02 17:58:14 perev
110 // gcc 3.2 updates + WarnOff
111 //
112 // Revision 1.49 2003/06/11 12:06:03 jcs
113 // get inner cathode and cluster geometry parameters from database
114 //
115 // Revision 1.48 2003/05/15 20:29:50 putschke
116 // bug fixed
117 //
118 // Revision 1.47 2003/05/07 15:09:22 putschke
119 // improvements for cathode offset corretions
120 //
121 // Revision 1.46 2003/04/23 15:13:13 putschke
122 // minor change in padtrans
123 //
124 // Revision 1.45 2003/04/15 11:35:39 putschke
125 // Include corrections for inner cathode offset and move some parameter to database
126 //
127 // Revision 1.44 2003/01/21 09:34:31 jcs
128 // initialize variables to eliminate compiler warnings for NODEBUG=yes
129 //
130 // Revision 1.43 2003/01/16 18:03:52 oldi
131 // Bug eliminated. Now it compiles on Solaris again.
132 //
133 // Revision 1.42 2003/01/14 12:58:00 jcs
134 // use Geometry_ftpc/ftpcAsicMap to control corrections for error in Y2001-2002
135 // FTPC asic mapping
136 //
137 // Revision 1.41 2002/11/15 11:52:05 jcs
138 // correct error in merging CUCs for statement condition
139 // (CurrentCUC->Sequence[159] was not filled for the case of >= MAXNUMSEQUENCES)
140 //
141 // Revision 1.40 2002/11/07 13:27:30 oldi
142 // Eliminated a very dumb mistake.
143 //
144 // Revision 1.39 2002/11/06 13:43:59 oldi
145 // Flag for clusters not to be used for tracking introduced.
146 //
147 // Revision 1.38 2002/08/02 13:03:55 oldi
148 // Wrong pad ordering was corrected in previous version already
149 // (see code changes there).
150 //
151 // Revision 1.37 2002/08/02 11:26:33 oldi
152 // Chargestep corrected (it was looping over the sequences twice).
153 // Pad ordering corrected in this version already.
154 //
155 // Revision 1.36 2002/07/15 13:30:35 jcs
156 // incorporate charge step histos into cluster finder and remove StFtpcChargeStep
157 //
158 // Revision 1.35 2002/06/04 12:33:09 putschke
159 // new 2-dimenisional hitfinding algorithm
160 // correct error in padposition numbering
161 //
162 // Revision 1.34 2002/04/22 09:53:45 jcs
163 // correct errors in calculation of cluster phi angle for ftpc east (Frank Simon)
164 //
165 // Revision 1.33 2002/04/05 16:45:54 oldi
166 // Small code clean ups, to be sure that this part is recompiled. It relies
167 // on StFtpcTracker/StFtpcPoint.* which were changed.
168 //
169 // Revision 1.32 2002/03/22 08:58:46 jcs
170 // invert cluster pad number in FTPC East
171 // convert cluster cartestian coordinates for FTPC west into STAR global coordinate system
172 //
173 // Revision 1.31 2002/03/22 08:52:52 jcs
174 // correct memory leaks found by Insure
175 //
176 // Revision 1.30 2002/03/08 16:32:34 jcs
177 // initialize bSlope to prevent nan's
178 //
179 // Revision 1.29 2002/03/05 16:51:35 jcs
180 // force data type definitions to avoid compiler warnings (this is a correct
181 // but inelegant fix which must be changed)
182 //
183 // Revision 1.28 2002/03/01 14:22:20 jcs
184 // add additional histograms to monitor cluster finding
185 //
186 // Revision 1.26 2002/02/10 21:10:44 jcs
187 // allow for individual west/east Ftpc temperature/pressure corrections
188 //
189 // Revision 1.25 2002/01/21 22:09:29 jcs
190 // use average FTPC gas temperature adjusted air pressure
191 //
192 // Revision 1.24 2001/12/12 16:01:17 jcs
193 // Remove old PhiDeflect definitions
194 // PhiDeflect has correct sign, add instead of subtracting in calculation of phi angle
195 //
196 // Revision 1.23 2001/11/21 12:51:27 jcs
197 // replace malloc/free with new/delete; prevent overstepping table boundaries
198 //
199 // Revision 1.22 2001/08/20 00:35:17 jcs
200 // check index (j) before using it
201 //
202 // Revision 1.21 2001/07/12 10:42:05 jcs
203 // reject clusters outside FTPC sensitive volume
204 // change to linear interpolation method in padtrans
205 // create and fill new histogram (currently inactive)
206 //
207 // Revision 1.20 2001/07/05 13:47:03 jcs
208 // move debug statement to correct location
209 //
210 // Revision 1.19 2001/06/25 00:45:35 jcs
211 // change DEBUG print statement to print out both soft and hard row/sector info
212 //
213 // Revision 1.18 2001/06/24 21:08:22 jcs
214 // Change Hit rejected message since this also occurs for FTPC daq data
215 //
216 // Revision 1.17 2001/04/23 19:37:40 oldi
217 // Output will be sent to StMessMgr.
218 //
219 // Revision 1.16 2001/04/19 11:32:39 oldi
220 // Reject clusters with x = 0. and y = 0.
221 // This is just to avoid a crash in StFtpcTrackMaker but the origin of the problem
222 // is in the FTPC slow simulator.
223 //
224 // Revision 1.15 2001/04/02 12:10:11 jcs
225 // get FTPC calibrations,geometry from MySQL database and code parameters
226 // from StarDb/ftpc
227 //
228 // Revision 1.14 2001/03/19 15:52:47 jcs
229 // use ftpcDimensions from database
230 //
231 // Revision 1.13 2001/03/06 23:33:38 jcs
232 // use database instead of params
233 //
234 // Revision 1.12 2001/01/25 15:25:21 oldi
235 // Fix of several bugs which caused memory leaks:
236 // - Some arrays were not allocated and/or deleted properly.
237 // - TClonesArray seems to have a problem (it could be that I used it in a
238 // wrong way in StFtpcTrackMaker form where Holm cut and pasted it).
239 // I changed all occurences to TObjArray which makes the program slightly
240 // slower but much more save (in terms of memory usage).
241 //
242 // Revision 1.11 2001/01/15 18:18:00 jcs
243 // replace hidden constants with parameters
244 //
245 // Revision 1.10 2000/11/27 14:09:07 hummler
246 // implement tzero and lorentz angle correction factor
247 //
248 // Revision 1.9 2000/11/14 13:08:04 hummler
249 // add charge step calculation, minor cleanup
250 //
251 // Revision 1.8 2000/08/03 14:39:00 hummler
252 // Create param reader to keep parameter tables away from cluster finder and
253 // fast simulator. StFtpcClusterFinder now knows nothing about tables anymore!
254 //
255 // Revision 1.6 2000/04/13 18:08:21 fine
256 // Adjusted for ROOT 2.24
257 //
258 // Revision 1.5 2000/01/27 09:47:16 hummler
259 // implement raw data reader, remove type ambiguities that bothered kcc
260 //
261 // Revision 1.4 2000/01/03 12:48:44 jcs
262 // Add CVS Id strings
263 //
264 
265 #include <Stiostream.h>
266 #include <stdlib.h>
267 #include "StMessMgr.h"
268 #include "StFtpcClusterFinder.hh"
269 #include "StFtpcTrackMaker/StFtpcConfMapPoint.hh"
270 #include "math_constants.h"
271 #include <math.h>
272 #include "TH1.h"
273 
274 #include "PhysicalConstants.h"
275 
276 #include "asic_map_correction.h"
277 
278 
279 // Constructor used for production
280 
281 StFtpcClusterFinder::StFtpcClusterFinder(StFTPCReader *reader,
282  StFtpcParamReader *paramReader,
283  StFtpcDbReader *dbReader,
284  TObjArray *pointarray,
285  TH2F *hpad,
286  TH2F *htime,
287  TH2F *histo,
288  TH1F *histoW,
289  TH1F *histoE)
290 {
291  LOG_DEBUG << "StFtpcClusterFinder constructed for production" << endm;
292  mReader = reader;
293  mParam = paramReader;
294  mDb = dbReader;
295  mPoint = pointarray;
296  mHisto=histo;
297  mHistoW=histoW;
298  mHistoE=histoE;
299 
300  mcldebug = NULL;
301 
302  MAXSEQPEAKS = mParam->maxNumSeqPeaks();
303  MAXPEAKS = mParam->maxNumPeaks();
304  MAXLOOPS = mParam->maxLoops();
305  MAXFASTLOOPS = mParam->maxFastLoops();
306  UNFOLDLIMIT = mParam->unfoldLimit();
307  UNFOLDFAILEDLIMIT = mParam->unfoldFailedLimit();
308 
309  mMinTimeBin = mDb->minTimeBin();
310  mMinTimeBinMed = mDb->minTimeBinMed();
311  mMinTimeBinOut = mDb->minTimeBinOut();
312 
313  mMaxPadlength = mDb->maxPadLength();
314  mMaxTimelength = mDb->maxTimeLength();
315  mMaxPadlengthMed = mDb->maxPadLengthMed();
316  mMaxTimelengthMed = mDb->maxTimeLengthMed();
317  mMaxPadlengthOut = mDb->maxPadLengthOut();
318  mMaxTimelengthOut = mDb->maxTimeLengthOut();
319 
320  DeltaTime = mDb->deltaTime();
321  DeltaPad = mDb->deltaPad();
322 
323  mMinChargeWindow = mDb->minChargeWindow();
324 
325  mOffsetCathodeWest = mDb->offsetCathodeWest();
326  mOffsetCathodeEast = mDb->offsetCathodeEast();
327  mAngleOffsetWest = mDb->angleOffsetWest();
328  mAngleOffsetEast = mDb->angleOffsetEast();
329 
330  mhpad = hpad;
331  mhtime = htime;
332 }
333 
334 // Constructor used to produce the FTPC special root file used for FTPC
335 // laser, t0 and cluster calibration
336 
337 StFtpcClusterFinder::StFtpcClusterFinder(StFTPCReader *reader,
338  StFtpcParamReader *paramReader,
339  StFtpcDbReader *dbReader,
340  TObjArray *pointarray,
341  TH2F *hpad,
342  TH2F *htime,
343  TH2F *histo,
344  TH1F *histoW,
345  TH1F *histoE,
346  StFtpcClusterDebug *cldebug)
347 {
348  LOG_DEBUG << "StFtpcClusterFinder constructed for calibration" << endm;
349  mReader = reader;
350  mParam = paramReader;
351  mDb = dbReader;
352  mPoint = pointarray;
353  mHisto=histo;
354  mHistoW=histoW;
355  mHistoE=histoE;
356 
357  mcldebug = cldebug;
358 
359  MAXSEQPEAKS = mParam->maxNumSeqPeaks();
360  MAXPEAKS = mParam->maxNumPeaks();
361  MAXLOOPS = mParam->maxLoops();
362  MAXFASTLOOPS = mParam->maxFastLoops();
363  UNFOLDLIMIT = mParam->unfoldLimit();
364  UNFOLDFAILEDLIMIT = mParam->unfoldFailedLimit();
365 
366  mMinTimeBin = mDb->minTimeBin();
367  mMinTimeBinMed = mDb->minTimeBinMed();
368  mMinTimeBinOut = mDb->minTimeBinOut();
369 
370  mMaxPadlength = mDb->maxPadLength();
371  mMaxTimelength = mDb->maxTimeLength();
372  mMaxPadlengthMed = mDb->maxPadLengthMed();
373  mMaxTimelengthMed = mDb->maxTimeLengthMed();
374  mMaxPadlengthOut = mDb->maxPadLengthOut();
375  mMaxTimelengthOut = mDb->maxTimeLengthOut();
376 
377  DeltaTime = mDb->deltaTime();
378  DeltaPad = mDb->deltaPad();
379 
380  mMinChargeWindow = mDb->minChargeWindow();
381 
382 // For studying cluster finding and tracking comment out the following lines of code
383 // from here ->
384 // Set FTPC inner cathode offsets = 0
385  mOffsetCathodeWest = 0.0;
386  mOffsetCathodeEast = 0.0;
387  mAngleOffsetWest = 0.0;
388  mAngleOffsetEast = 0.0;
389 
390 
391 // Set ftpcClusterPars
392  MAXPEAKS = 160;
393 
394 // Set FTPC Cluster Geometry values for cluster/laser analysis
395  mMinChargeWindow = 30;
396  mMaxPadlengthOut = 30;
397  mMaxTimelengthOut = 30;
398 // <- to here
399 
400  mhpad = hpad;
401  mhtime = htime;
402 
403 }
404 
405 StFtpcClusterFinder::~StFtpcClusterFinder()
406 {
407 // LOG_DEBUG << "StFtpcClusterFinder destructed" << endm;
408 }
409 
410 int StFtpcClusterFinder::search()
411 {
412 
413  Double_t *pradius = 0;
414  Double_t *pdeflection = 0;
415  int iRow, iSec, iPad, iPadBuf;
416  int iRowBuf, iSecBuf;
417  int firstPadrowToSearch;
418  int bNewSec;
419  int clusters;
420  int iNowSeqIndex, iNewSeqIndex, iOldSeqNumber, iOldSeqIndex;
421  int iCUCSequence, iMoveSequence, iOldSeqBuf;
422  int bOldSequenceUsed, bLastSequence;
423  TClusterUC *FirstCUC, *CurrentCUC, *LastCUC, *DeleteCUC;
424  TClusterUC *SequenceInCUC;
425  TPCSequence *OldSequences, *NewSequences, *(SequencePointer[3]);
426  int iNewSeqNumber;
427 
428  int iIndex;
429  float fastlog[256];
430 
431  double deltaAirPressure;
432 
433  /* variables for dynamic CUC memory handling */
434  TClusterUC CUCMemory[MAXNUMCUC];
435  int CUCMemoryArray[MAXNUMCUC];
436  int CUCMemoryPtr;
437 
438  /* allocate memory for padtrans table */
439  pradius = new Double_t[mParam->numberOfDriftSteps()
440  *mDb->numberOfPadrowsPerSide()];
441  memset(pradius, 0, (mParam->numberOfDriftSteps()*mDb->numberOfPadrowsPerSide())*sizeof(Double_t));
442  pdeflection = new Double_t[mParam->numberOfDriftSteps()
443  *mDb->numberOfPadrowsPerSide()];
444  memset(pdeflection, 0, (mParam->numberOfDriftSteps()*mDb->numberOfPadrowsPerSide())*sizeof(Double_t));
445 
446  if(pradius == 0 || pdeflection == 0)
447  {
448  LOG_ERROR << "Padtrans memory allocation failed, exiting!" << endm;
449  if (pradius != 0) delete[] pradius; // release the pradius array
450  if (pdeflection != 0) delete[] pdeflection; // release the pdeflection array
451  return kStERR;
452  }
453 
454 // Loop over FTPC West and East individually
455  // initialize to elimate compiler warning when setenv NODEBUG yes
456  firstPadrowToSearch = 0;
457  deltaAirPressure = 0;
458 
459 Int_t westHits = 0;
460 Int_t eastHits = 0;
461 for ( int iftpc=0; iftpc<2; iftpc++) {
462  if ( iftpc == 0 ) {
463  deltaAirPressure = mParam->adjustedAirPressureWest() - mParam->standardPressure();
464  firstPadrowToSearch = mDb->firstPadrowToSearch() - 1;
465  LOG_DEBUG << "Ftpc West: deltaAirPressure = adjustedAirPressureWest ("<<mParam->adjustedAirPressureWest()<<") - standardPressure ("<<mParam->standardPressure()<<") = "<<deltaAirPressure<<endm;
466  }
467  if ( iftpc == 1 ) {
468  deltaAirPressure = mParam->adjustedAirPressureEast() - mParam->standardPressure();
469  firstPadrowToSearch = mDb->firstPadrowToSearch() - 1 + mDb->numberOfPadrowsPerSide();
470  LOG_DEBUG <<"Ftpc East: deltaAirPressure = adjustedAirPressureEast ("<<mParam->adjustedAirPressureEast()<<") - standardPressure ("<<mParam->standardPressure()<<") = "<<deltaAirPressure<<endm;
471  }
472 
473  /* integrate padtrans table from magboltz database */
474  if(!calcpadtrans(pradius, pdeflection,deltaAirPressure))
475  {
476  LOG_ERROR << "Couldn't calculate padtrans table, exiting!" << endm;
477  delete[] pradius; // release the pradius array
478  delete[] pdeflection; // release the pdeflection array
479  return kStERR;
480  }
481 
482  /* initialize CUC memory handling */
483  if(!cucInit(CUCMemory, CUCMemoryArray, &CUCMemoryPtr))
484  {
485  LOG_ERROR << "Couldn't initialize CUC memory, exiting!" << endm;
486  delete[] pradius; // release the pradius array
487  delete[] pdeflection; // release the pdeflection array
488  return kStERR;
489  }
490 
491  /* calculate fastlog lookup */
492  for(iIndex=1; iIndex<256; iIndex++)
493  {
494  fastlog[iIndex] = ::log((double) iIndex);
495  }
496 
497  /* reset counter for found clusters */
498  clusters = 0;
499 
500  /* initialize sequence and cluster lists */
501  NewSequences = 0;
502  FirstCUC = 0;
503  LastCUC = 0;
504  iOldSeqNumber = 0;
505  iNewSeqIndex = 0;
506 
507  /* loop over raw data sequences */
508  iNowSeqIndex = 0;
509  iPad=0;
510  iSec=0;
511  iRow=0;
512  bLastSequence=0;
513 
514  for(iRow=firstPadrowToSearch,iRowBuf=firstPadrowToSearch;iRow<firstPadrowToSearch+mDb->numberOfPadrowsPerSide(); iRow++)
515  {
516  for(iSec=mDb->firstSectorToSearch()-1,
517  iSecBuf=mDb->firstSectorToSearch()-1;
518  iSec<mDb->lastSectorToSearch(); iSec++)
519  {
520  // new sector, set pad buffer so there can be no matches
521  iPadBuf=-2;
522  bNewSec=TRUE;
523 
524  // calculate hardware (daq) sectors from software position
525  iHardSec = mDb->numberOfSectors()*(int)(iRow/2) + iSec + 1;
526  iHardRow = iRow%2 + 1;
527 
528  if (DebugOn) {
529  LOG_DEBUG<<"Cluster Finder: Now on Sector "<<iSec<<", Row "<<iRow<<" (iHardSec "<<iHardSec<<", iHardRow "<<iHardRow<<")"<<endm;
530  }
531 
532  // get list of occupied pads in sector
533  unsigned char *(padlist[2]);
534  int iOccPads=mReader->getPadList(iHardSec, iHardRow,
535  padlist[iHardRow-1]);
536 
537  // loop over occupied pads
538  int iThPad;
539 //===============================================================
540  int newpadlist[160];
541 
542  for (iThPad=0; iThPad<160; iThPad++)
543  newpadlist[iThPad]=0;
544 
545  for(iThPad=0; iThPad<iOccPads; iThPad++)
546  {
547  iPad=padlist[iHardRow-1][iThPad];
548  if ( mDb->Asic2EastNotInverted() && iRow>=10 && (iPad>=65 && iPad<=96))
549  newpadlist[padkey[iPad-1]-1] = iPad;
550  else
551  newpadlist[iPad-1] = iPad;
552  }
553 
554 //===============================================================
555  for(iThPad=0; iThPad<160; iThPad++)
556  {
557  if (newpadlist[iThPad] == 0 ) continue;
558  iPad=iThPad+1;
559 
560  // search, fit and remove finished CUCs
561  for(CurrentCUC = FirstCUC; CurrentCUC!=0;
562  CurrentCUC = CurrentCUC->NextClusterUC)
563  {
564  // check if CurrentCUC is adjacent to this pad
565  if(iPad > CurrentCUC->EndPad + 1 || bNewSec)
566  {
567  // CurrentCUC is finished
568  // if CUC has not been lost by merging, process cluster
569  if(CurrentCUC->EndPad > CurrentCUC->StartPad)
570  {
571  // cluster processing: increment cluster counter
572  clusters ++;
573 
574  // cluster processing: call hitfinder
575  if (geometryCut(CurrentCUC))
576  if(!findHits(CurrentCUC, iRowBuf, iSecBuf,
577  pradius, pdeflection,
578  fastlog)
579  )
580  {
581  if (DebugOn) {
582  LOG_DEBUG<<"Hitfinder failed! Cluster is lost"<<endm;
583  }
584  }
585  }
586  DeleteCUC=CurrentCUC;
587  // bypass CurrentCUC in CUC list
588  if(CurrentCUC==FirstCUC)
589  {
590  FirstCUC=CurrentCUC->NextClusterUC;
591  }
592  else
593  {
594  LastCUC->NextClusterUC=CurrentCUC->NextClusterUC;
595  CurrentCUC=LastCUC;
596  }
597  // free CurrentCUC memory
598  if(!cucFree(CUCMemory, CUCMemoryArray,
599  &CUCMemoryPtr, DeleteCUC))
600  {
601  LOG_ERROR << "Fatal memory management error." << endm;
602  delete[] pradius; // release the pradius array
603  delete[] pdeflection; // release the pdeflection array
604  return kStERR;
605  }
606  }
607  LastCUC=CurrentCUC;
608  }
609  iRowBuf=iRow;
610  iSecBuf=iSec;
611 
612  // initialize sequence lists:
613  // new-array is moved to old
614  // and first element initialized
615  OldSequences=NewSequences;
616  iOldSeqNumber=iNewSeqIndex;
617  iNewSeqIndex=0;
618 
619  if(iPad!=iPadBuf+1)
620  {
621  iOldSeqNumber=0;
622  }
623  iPadBuf=iPad;
624 
625  // reset beginning of sequence comparison
626  iOldSeqBuf=0;
627 
628  // get sequences on this pad
629  mReader->getSequences(iHardSec, iHardRow, newpadlist[iThPad], &iNewSeqNumber,
630  SequencePointer[iHardRow]);
631  NewSequences=SequencePointer[iHardRow];
632 
633  // loop over sequences
634  for(iNewSeqIndex=0; iNewSeqIndex < iNewSeqNumber;
635  iNewSeqIndex++)
636  {
637 if (mcldebug){
638  if (mcldebug->drawclhisto!=0)
639  {
640  mcldebug->drawhisto(iHardSec,iHardRow,iPad,NewSequences[iNewSeqIndex]);
641  mcldebug->drawgainhisto(iHardSec,iHardRow,iPad,((float)(mDb->amplitudeSlope((iSec*mDb->numberOfPads()+iPad),iRow))),NewSequences[iNewSeqIndex]);
642  }
643 }
644 //+++++++++++++++++++ fill charge step histograms +++++++++++++++
645  // This loop is running already, but the running variable is called iNewSeqIndex instead of iSeqIndex .
646  //int iSeqIndex;
647  //for(iSeqIndex=0; iSeqIndex < iNewSeqNumber; iSeqIndex++)
648  //{
649  int entry;
650  for(entry=0; entry<NewSequences[iNewSeqIndex].Length; entry++)
651  {
652  if (mHisto) mHisto->Fill(iHardSec-1, // sector
653  entry+NewSequences[iNewSeqIndex].startTimeBin, //bin
654  NewSequences[iNewSeqIndex].FirstAdc[entry]); // weight
655  if (iHardSec >= 1 && iHardSec <= 30 ) {
656  mHistoW->Fill( entry+NewSequences[iNewSeqIndex].startTimeBin, //bin
657  NewSequences[iNewSeqIndex].FirstAdc[entry]); // weight
658  }
659  if (iHardSec >= 31 && iHardSec <= 60 ) {
660  mHistoE->Fill( entry+NewSequences[iNewSeqIndex].startTimeBin, //bin
661 
662  NewSequences[iNewSeqIndex].FirstAdc[entry]); // weight
663  }
664  }
665  //}
666 
667 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
668 
669  // mark sequence as unused
670  SequenceInCUC = 0;
671 
672  // compare this sequence to old sequences
673  for(iOldSeqIndex=iOldSeqBuf; iOldSeqIndex < iOldSeqNumber;
674  iOldSeqIndex++)
675  {
676  // are beginning or end of new sequence between
677  // beginning and end of old sequence?
678  if(((NewSequences[iNewSeqIndex].startTimeBin >=
679  OldSequences[iOldSeqIndex].startTimeBin) &&
680  (NewSequences[iNewSeqIndex].startTimeBin <=
681  OldSequences[iOldSeqIndex].startTimeBin +
682  OldSequences[iOldSeqIndex].Length-1)) ||
683  ((NewSequences[iNewSeqIndex].startTimeBin +
684  NewSequences[iNewSeqIndex].Length-1 >=
685  OldSequences[iOldSeqIndex].startTimeBin) &&
686  (NewSequences[iNewSeqIndex].startTimeBin +
687  NewSequences[iNewSeqIndex].Length-1 <=
688  OldSequences[iOldSeqIndex].startTimeBin +
689  OldSequences[iOldSeqIndex].Length-1)) ||
690  ((OldSequences[iOldSeqIndex].startTimeBin >=
691  NewSequences[iNewSeqIndex].startTimeBin) &&
692  (OldSequences[iOldSeqIndex].startTimeBin <=
693  NewSequences[iNewSeqIndex].startTimeBin +
694  NewSequences[iNewSeqIndex].Length-1)))
695  {
696  // yes, matching sequence found
697  // set beginning of search for next sequence
698  iOldSeqBuf=iOldSeqIndex;
699  bOldSequenceUsed=0;
700 
701  // compare matching sequences to old CUCs
702  // loop over all active CUCs
703  for(CurrentCUC = FirstCUC; CurrentCUC != 0;
704  CurrentCUC = CurrentCUC->NextClusterUC)
705  {
706  LastCUC=CurrentCUC;
707  // loop over all sequences in CUC
708  for(iCUCSequence=1;
709  iCUCSequence<CurrentCUC->NumSequences;
710  iCUCSequence++)
711  {
712  // is cuc sequence identical to
713  // matching old sequence?
714  if((OldSequences[iOldSeqIndex].startTimeBin ==
715  CurrentCUC->Sequence[iCUCSequence].startTimeBin)
716  && (CurrentCUC->SequencePad[iCUCSequence] ==
717  iPad-1))
718  {
719  bOldSequenceUsed=1;
720  // matching old sequence is in CUC
721  // check if new sequence is already used
722  if(SequenceInCUC!=CurrentCUC)
723  {
724  if(SequenceInCUC != 0 && SequenceInCUC!=CurrentCUC)
725  {
726  // yes, already used, merge CUCs
727  // mark old CUC for removal
728  SequenceInCUC->EndPad =
729  SequenceInCUC->StartPad;
730  // set StartPad to the smaller StartPad
731  if(SequenceInCUC->StartPad <
732  CurrentCUC->StartPad)
733  {
734  CurrentCUC->StartPad =
735  SequenceInCUC->StartPad;
736  }
737  CurrentCUC->EndPad=iPad;
738  // append SequenceInCUC to CurrentCUC
739  // copy all sequences to CurrentCUC
740  for(iMoveSequence=0;
741  (iMoveSequence <SequenceInCUC->NumSequences) &&
742  ( (CurrentCUC->NumSequences+iMoveSequence) < MAXNUMSEQUENCES);
743  iMoveSequence++)
744  {
745  CurrentCUC->
746  Sequence[iMoveSequence +
747  CurrentCUC->NumSequences] =
748  SequenceInCUC->
749  Sequence[iMoveSequence];
750  CurrentCUC->
751  SequencePad[iMoveSequence +
752  CurrentCUC->NumSequences] =
753  SequenceInCUC->
754  SequencePad[iMoveSequence];
755  }
756  // add up number of sequences
757  CurrentCUC->NumSequences +=
758  SequenceInCUC->NumSequences;
759  if(CurrentCUC->NumSequences > MAXNUMSEQUENCES)
760  {
761  CurrentCUC->NumSequences = MAXNUMSEQUENCES;
762  }
763 
764  // sequence is now in CurrentCUC
765  SequenceInCUC=CurrentCUC;
766  }
767  else // to: if(SequenceInCUC!=0)
768  {
769  // sequence did not belong to any CUC before
770  // add sequence to CurrentCUC
771  if(CurrentCUC->NumSequences<MAXNUMSEQUENCES)
772  {
773  CurrentCUC->Sequence[CurrentCUC->NumSequences]
774  = NewSequences[iNewSeqIndex];
775  CurrentCUC->SequencePad[CurrentCUC->NumSequences]
776  =iPad;
777  CurrentCUC->NumSequences++;
778  }
779  CurrentCUC->EndPad=iPad;
780  SequenceInCUC=CurrentCUC;
781 
782  // check if new sequence touches sector limit
783  if(NewSequences[iNewSeqIndex].startTimeBin==0 ||
784  NewSequences[iNewSeqIndex].startTimeBin
785  +NewSequences[iNewSeqIndex].Length
786  ==mDb->numberOfTimebins()-1 ||
787  iPad==mDb->numberOfPads())
788  {
789  CurrentCUC->CutOff=1;
790  }
791  } // end of: if(SequenceInCUC!=0) ... else
792  } // end of: if(SequenceInCUC...)
793  } // end of: if((OldSequences...))
794  } // end of: for(ICUCSequence...)
795  } // end of: for(CurrentCUC...)
796  if(SequenceInCUC==0 && NewSequences[iNewSeqIndex].Length>1)
797  {
798  // no matching CUC was found: create new CUC
799  // allocate memory
800  CurrentCUC=cucAlloc(CUCMemory, CUCMemoryArray,
801  &CUCMemoryPtr);
802  if(CurrentCUC == 0)
803  {
804  // no free memory, overwrite last CUC
805  if (DebugOn) {
806  LOG_DEBUG<<"Previous cluster is now lost"<<endm;
807  }
808  CurrentCUC=LastCUC;
809  delete[] pradius; // release the pradius array
810  delete[] pdeflection; // release the pdeflection array
811  return kStWarn;
812  }
813  else
814  {
815  // set pointers to this CUC
816  if(FirstCUC == 0)
817  {
818  FirstCUC = CurrentCUC;
819  }
820  else
821  {
822  LastCUC->NextClusterUC=CurrentCUC;
823  }
824 
825  // this is the newest CUC
826  CurrentCUC->NextClusterUC=0;
827  }
828 
829  // fill new CUC structure
830  CurrentCUC->StartPad=iPad-1;
831  CurrentCUC->EndPad=iPad;
832  CurrentCUC->NumSequences=2;
833 
834  // copy sequences to CUC
835  CurrentCUC->Sequence[0] =
836  OldSequences[iOldSeqIndex];
837  // struct assignment operator used
838  // created copy of sequence in CUC
839  // pointer operation avoided because data in CUC
840  // will be calibrated
841  CurrentCUC->SequencePad[0] = iPad-1;
842 
843  CurrentCUC->Sequence[1] =
844  NewSequences[iNewSeqIndex];
845  CurrentCUC->SequencePad[1]=iPad;
846  SequenceInCUC=CurrentCUC;
847 
848  // check if new CUC touches sector limits
849  if(iPad==1 || iPad==mDb->numberOfPads() ||
850  CurrentCUC->Sequence[0].startTimeBin==0 ||
851  CurrentCUC->Sequence[0].startTimeBin
852  +CurrentCUC->Sequence[0].Length
853  ==mDb->numberOfTimebins()-1 ||
854  CurrentCUC->Sequence[1].startTimeBin==0 ||
855  CurrentCUC->Sequence[1].startTimeBin
856  +CurrentCUC->Sequence[1].Length
857  ==mDb->numberOfTimebins()-1)
858  {
859  CurrentCUC->CutOff=1;
860  }
861  else
862  {
863  CurrentCUC->CutOff=0;
864  }
865  } // end of: if(SequenceInCUC==0)
866  else
867  {
868  if(bOldSequenceUsed==0 && SequenceInCUC!=0)
869  {
870  // new sequence has been used but old one hasn't
871  // append that to cluster, too
872  if(SequenceInCUC->NumSequences<MAXNUMSEQUENCES)
873  {
874  SequenceInCUC->Sequence[SequenceInCUC->NumSequences]
875  = OldSequences[iOldSeqIndex];
876  SequenceInCUC->SequencePad[SequenceInCUC->NumSequences]
877  =iPad-1;
878  SequenceInCUC->NumSequences++;
879  }
880 
881  // check if Old sequence touches sector limit
882  if(OldSequences[iOldSeqIndex].startTimeBin==0 ||
883  OldSequences[iOldSeqIndex].startTimeBin
884  +OldSequences[iOldSeqIndex].Length
885  ==mDb->numberOfTimebins()-1 ||
886  iPad>=mDb->numberOfPads()-1) //change == to >=
887  {
888  SequenceInCUC->CutOff=1;
889  }
890  }
891  }
892  } // end of: if(sequence matching)
893  } // end of: for(iOldSeqIndex...)
894  } // end of: for(iNewSeqIndex...)
895  bNewSec=FALSE;
896  } // end of: for(iThPad...)
897 
898  // sector done, fit and remove all remaining CUCs
899  for(CurrentCUC = FirstCUC; CurrentCUC!=0;
900  CurrentCUC = CurrentCUC->NextClusterUC)
901  {
902  // CurrentCUC is finished
903  // if CUC has not been lost by merging, process cluster
904  if(CurrentCUC->EndPad > CurrentCUC->StartPad)
905  {
906  // cluster processing: increment cluster counter
907  clusters ++;
908 
909  // cluster processing: call hitfinder
910  if (geometryCut(CurrentCUC))
911  if(!findHits(CurrentCUC, iRowBuf, iSecBuf,
912  pradius, pdeflection,
913  fastlog)
914  )
915  {
916  if (DebugOn) {
917  LOG_DEBUG<<"Hitfinder failed! Cluster is lost"<<endm;
918  }
919  }
920  }
921  DeleteCUC=CurrentCUC;
922  // bypass CurrentCUC in CUC list
923  if(CurrentCUC==FirstCUC)
924  {
925  FirstCUC=CurrentCUC->NextClusterUC;
926  }
927  else
928  {
929  LastCUC->NextClusterUC=CurrentCUC->NextClusterUC;
930  CurrentCUC=LastCUC;
931  }
932  // free CurrentCUC memory
933  if(!cucFree(CUCMemory, CUCMemoryArray,
934  &CUCMemoryPtr, DeleteCUC))
935  {
936  LOG_ERROR << "Fatal memory management error." << endm;
937  delete[] pradius; // release the pradius array
938  delete[] pdeflection; // release the pdeflection array
939  return kStERR;
940  }
941  LastCUC=CurrentCUC;
942  }
943 
944  } // end of: for(iSec...)
945  } // end of: for(iRow...)
946 
947  if (iftpc == 0 ) {
948  westHits = mPoint->GetEntriesFast();
949  LOG_INFO << "StFtpcClusterFinder found " << clusters << " clusters and processed to " << westHits << " hits in Ftpc West." << endm;
950  }
951  if (iftpc == 1 ) {
952  eastHits = mPoint->GetEntriesFast() - westHits;
953  LOG_INFO << "StFtpcClusterFinder found " << clusters << " clusters and processed to " << eastHits << " hits in Ftpc East." << endm;
954  }
955 } // end of: for(iftpc
956 
957  delete[] pradius; // release the pradius array
958  delete[] pdeflection; // release the pdeflection array
959 
960  return kStOK;
961 }
962 
963 bool StFtpcClusterFinder::geometryCut(TClusterUC *Cluster)
964 {
965 
966  int seqlength=0;
967  int minTimebin=256;
968  int maxTimebin=0;
969 
970  for(int iPad=Cluster->StartPad; iPad<=Cluster->EndPad; iPad++)
971  {
972  for(int iSequence = 0; iSequence < Cluster->NumSequences; iSequence++)
973  {
974  if(Cluster->SequencePad[iSequence] == iPad)
975  {
976  if (Cluster->Sequence[iSequence].Length > seqlength)
977  seqlength=Cluster->Sequence[iSequence].Length;
978  if (Cluster->Sequence[iSequence].startTimeBin<minTimebin)
979  minTimebin=Cluster->Sequence[iSequence].startTimeBin;
980  if ((Cluster->Sequence[iSequence].startTimeBin+Cluster->Sequence[iSequence].Length) > maxTimebin)
981  maxTimebin=(Cluster->Sequence[iSequence].startTimeBin+Cluster->Sequence[iSequence].Length);
982  }
983  }
984  }
985 
986 if (mcldebug) {
987  if (abs(Cluster->EndPad-Cluster->StartPad)<mMaxPadlengthOut && seqlength<mMaxTimelengthOut)
988  return true;
989  else
990  return false;
991 }
992  if (minTimebin>mMinTimeBin) return true;
993  else if ((minTimebin>mMinTimeBinMed && minTimebin<=mMinTimeBin) && abs(Cluster->EndPad-Cluster->StartPad)<mMaxPadlengthMed && seqlength<mMaxTimelengthMed)
994  return true;
995  else if (minTimebin>mMinTimeBinOut && minTimebin<=mMinTimeBinMed && abs(Cluster->EndPad-Cluster->StartPad)<mMaxPadlengthOut && seqlength<mMaxTimelengthOut)
996  return true;
997  else
998  return false;
999 
1000 }
1001 
1002 int StFtpcClusterFinder::findHits(TClusterUC *Cluster,
1003  int iRow,
1004  int iSec,
1005  double *pRadius,
1006  double *pDeflection,
1007  float fastlog[256])
1008 {
1009 
1010  int iNumPeaks;
1011  int i,k;
1012 
1013  bool PeakFound;
1014 
1015  TPeak *Peaks = 0;
1016  Peaks = new TPeak[MAXPEAKS];
1017  float TClSearch [162][258];
1018 
1019  for (int i=0;i<162;i++)
1020  for (int j=0;j<258;j++)
1021  TClSearch[i][j]=0;
1022 
1023  float cTemp;
1024 
1025  iNumPeaks=0;
1026 
1027  // read out cluster and fill window with gain correction
1028  // =====================================================
1029 
1030  for(int iPad=Cluster->StartPad; iPad<=Cluster->EndPad; iPad++)
1031  {
1032  for(int iSequence = 0; iSequence < Cluster->NumSequences; iSequence++)
1033  {
1034  if(Cluster->SequencePad[iSequence] == iPad)
1035  {
1036  for(int iIndex=0; iIndex< Cluster->Sequence[iSequence].Length; iIndex++)
1037  {
1038  {
1039  cTemp=((float)(unsigned int)(Cluster->Sequence[iSequence].FirstAdc[iIndex])
1040  * mDb->amplitudeSlope(iSec*mDb->numberOfPads()+iPad,iRow)
1041  + mDb->amplitudeOffset(iSec*mDb->numberOfPads()+iPad,iRow));
1042  Cluster->Sequence[iSequence].FirstAdc[iIndex]=(unsigned char) cTemp;
1043  TClSearch[iPad][Cluster->Sequence[iSequence].startTimeBin+iIndex]=cTemp;
1044  }
1045  }
1046  }
1047  }
1048  }
1049 
1050  // check each cluster pixel for local maximum
1051  // ==========================================
1052 
1053  float cl_charge;
1054 
1055  for(int iPad=Cluster->StartPad; iPad<=Cluster->EndPad; iPad++)
1056  {
1057  for(int iSequence = 0; iSequence < Cluster->NumSequences; iSequence++)
1058  {
1059  if(Cluster->SequencePad[iSequence] == iPad)
1060  {
1061  for(int iTime=Cluster->Sequence[iSequence].startTimeBin; iTime <=Cluster->Sequence[iSequence].startTimeBin+Cluster->Sequence[iSequence].Length; iTime++)
1062  {
1063 
1064  // ------------------------------------------------------------------
1065 
1066  PeakFound=true;
1067  cl_charge=0;
1068 
1069  // search window
1070  // =============
1071 
1072  // check if > MaximumPad/Time
1073  for (i=iPad-DeltaPad;i<=(iPad+DeltaPad) && (iPad+DeltaPad)<=160 && (iPad-DeltaPad)>=0;i++)
1074  {
1075  for (k=iTime-DeltaTime;k<=(iTime+DeltaTime) && (iTime+DeltaTime)<=256 && (iTime-DeltaTime)>=0;k++)
1076  {
1077  if (i<iPad+DeltaPad && i>iPad-DeltaPad && k<iTime+DeltaTime && k>iTime-DeltaTime)
1078  cl_charge=cl_charge+TClSearch[i][k];
1079 
1080  if (TClSearch[iPad][iTime]>=TClSearch[i][k]+1.25 && TClSearch[iPad][iTime]>0
1081  && TClSearch[i][k]>0)
1082  {
1083  }
1084  else if (TClSearch[iPad][iTime]<TClSearch[i][k] && TClSearch[i][k]>0 && TClSearch[i][k]>0)
1085  {
1086  PeakFound=false;
1087  }
1088  }
1089  }
1090 
1091  // if local maximum Found make Peak
1092  // ================================
1093 
1094  if (PeakFound && cl_charge>mMinChargeWindow && iNumPeaks<MAXPEAKS)
1095  {
1096 
1097  if (iNumPeaks>0)
1098  {
1099  if ((Peaks[iNumPeaks-1].Timebin!=iTime && Peaks[iNumPeaks-1].Timebin!=iTime-1 && Peaks[iNumPeaks-1].Timebin!=iTime+1) || Peaks[iNumPeaks-1].pad!=iPad )
1100  {
1101  Peaks[iNumPeaks].pad=iPad;
1102  Peaks[iNumPeaks].Timebin=iTime;
1103  Peaks[iNumPeaks].Sequence=Cluster->Sequence[iSequence];
1104  iNumPeaks++;
1105  }
1106  }
1107  else
1108  {
1109  Peaks[iNumPeaks].pad=iPad;
1110  Peaks[iNumPeaks].Timebin=iTime;
1111  Peaks[iNumPeaks].Sequence=Cluster->Sequence[iSequence];
1112  iNumPeaks++;
1113  }
1114  }
1115 
1116  // ------------------------------------------------------------------
1117 
1118  }
1119  }
1120 
1121  }
1122  }
1123 
1124  if(!fitPoints(Cluster, iRow, iSec, pRadius, pDeflection, Peaks,
1125  iNumPeaks, fastlog))
1126  {
1127  if (DebugOn) {
1128  LOG_DEBUG<<"Point fitting failed! Cluster is lost."<<endm;
1129  }
1130  delete[] Peaks;
1131  return FALSE;
1132  }
1133 
1134  delete[] Peaks;
1135  return TRUE;
1136 }
1137 
1138 int StFtpcClusterFinder::fitPoints(TClusterUC* Cluster,
1139  int iRow,
1140  int iSec,
1141  double *pRadius,
1142  double *pDeflection,
1143  TPeak *Peak,
1144  int iNumPeaks,
1145  float fastlog[256])
1146 {
1147  int iRowSAVE, iSecSAVE;
1148  int iADCValue, iADCPlus, iADCMinus;
1149  int iADCTimePlus, iADCTimeMinus;
1150  int iSequence, iBin;
1151  int ChargeSum, PadSum;
1152  float TimeSum;
1153  int iUseGauss;
1154  int iPeakIndex, iInnerIndex;
1155  int iNumUnfoldLoops;
1156  int BadFit, PeakShifted, FailedFit;
1157  int PadtransPerTimebin, PadtransBin;
1158  float SumDeltaPos, SumDeltaAmp;
1159  float NewTimePosition, NewPadPosition;
1160  float NewPeakHeight, PeakHeightSum;
1161  float fDeltaADC, fDeltaADCPlus, fDeltaADCMinus;
1162  float fDeltaADCTimePlus, fDeltaADCTimeMinus;
1163  float fDriftLength, fRadError, fPhiError;
1164 
1165  PadtransPerTimebin=(int) mParam->numberOfDriftSteps()
1166  / mDb->numberOfTimebins();
1167 
1168  if(iNumPeaks == 0)
1169  {
1170  if (DebugOn) {
1171  LOG_DEBUG<<"Cluster starting "<<Cluster->StartPad<<", "<<Cluster->Sequence->startTimeBin<<" has no peak!"<<endm;
1172  }
1173  return FALSE;
1174  }
1175 
1176  /* sum up cluster charge */
1177  ChargeSum=0;
1178  PadSum=0;
1179  TimeSum=0;
1180  iADCValue=0;
1181  for(iSequence = 0; iSequence < Cluster->NumSequences; iSequence++)
1182  {
1183  for(iBin = 0; iBin < Cluster->Sequence[iSequence].Length; iBin++)
1184  {
1185  iADCValue = Cluster->Sequence[iSequence].FirstAdc[iBin];
1186  ChargeSum += iADCValue;
1187  PadSum += Cluster->SequencePad[iSequence] * iADCValue;
1188  TimeSum += (Cluster->Sequence[iSequence].startTimeBin + iBin +
1189  mDb->timeOffset(iSec*mDb->numberOfPads()
1190  +Cluster->SequencePad[iSequence],iRow))
1191  * iADCValue;
1192  }
1193  }
1194 
1195  if(iNumPeaks == 1)
1196  {
1197  Peak->PeakHeight =
1198  Peak->Sequence.FirstAdc[Peak->Timebin - Peak->Sequence.startTimeBin];
1199 
1200  /* one-peak cluster, fit according to preference from det */
1201  iUseGauss=mParam->gaussFittingFlags();
1202 
1203 
1204  /* calculate pad position first */
1205  if((iUseGauss & 1) == 1)
1206  {
1207  /* get values for gaussfit */
1208  iADCValue = (int) Peak->PeakHeight;
1209  iADCPlus=0;
1210  iADCMinus=0;
1211  /* find sequences with neighbor pixels */
1212  for(iSequence = 0; iSequence < Cluster->NumSequences; iSequence++)
1213  {
1214  if((Cluster->SequencePad[iSequence] == Peak->pad - 1) &&
1215  (Cluster->Sequence[iSequence].startTimeBin <= Peak->Timebin) &&
1216  (Cluster->Sequence[iSequence].startTimeBin +
1217  Cluster->Sequence[iSequence].Length > Peak->Timebin))
1218  {
1219  iADCMinus = Cluster->Sequence[iSequence]
1220  .FirstAdc[Peak->Timebin
1221  - Cluster->Sequence[iSequence].startTimeBin];
1222  }
1223  if((Cluster->SequencePad[iSequence] == Peak->pad + 1) &&
1224  (Cluster->Sequence[iSequence].startTimeBin < Peak->Timebin) &&
1225  (Cluster->Sequence[iSequence].startTimeBin +
1226  Cluster->Sequence[iSequence].Length > Peak->Timebin))
1227  {
1228  iADCPlus = Cluster->Sequence[iSequence]
1229  .FirstAdc[Peak->Timebin
1230  - Cluster->Sequence[iSequence].startTimeBin];
1231  }
1232  }
1233 
1234  /* check if Gaussfit is possible */
1235  if((iADCValue == 0) || (iADCPlus == 0) || (iADCMinus == 0) ||
1236  ((iADCValue <= iADCPlus) && (iADCValue <= iADCMinus)))
1237  {
1238  /* use weighted mean instead */
1239  iUseGauss = iUseGauss & 2;
1240  }
1241  else
1242  {
1243  /* do gaussfit */
1244  Peak->PadSigma = sqrt (1 /
1245  ((2 * fastlog[iADCValue]) -
1246  (fastlog[iADCPlus] + fastlog[iADCMinus])));
1247  Peak->PadPosition = (float) Peak->pad +
1248  sqr(Peak->PadSigma) * (fastlog[iADCPlus] - fastlog[iADCMinus]);
1249  }
1250  }
1251 
1252  if((iUseGauss & 1) == 0)
1253  {
1254  /* calculate pad position by weighted mean */
1255  Peak->PadPosition = (float) PadSum / (float) ChargeSum;
1256 
1257  /* calculate pad direction cluster width by weighted mean */
1258  Peak->PadSigma=0;
1259  for(iSequence = 0; iSequence < Cluster->NumSequences; iSequence++)
1260  {
1261  for(iBin = 0; iBin < Cluster->Sequence[iSequence].Length; iBin++)
1262  {
1263  Peak->PadSigma+=(sqr((float) Cluster->SequencePad[iSequence]
1264  - Peak->PadPosition)
1265  *Cluster->Sequence[iSequence].FirstAdc[iBin]);
1266  }
1267  }
1268  Peak->PadSigma /= (float) ChargeSum;
1269  }
1270 
1271  /* now calculate the time position */
1272  if((iUseGauss & 2) > 0)
1273  {
1274  /* get values for gaussfit */
1275  iADCValue = (int) Peak->PeakHeight;
1276  iADCPlus=0;
1277  iADCMinus=0;
1278 
1279  if(Peak->Timebin > Peak->Sequence.startTimeBin)
1280  {
1281  iADCMinus = Peak->Sequence.FirstAdc[Peak->Timebin -
1282  Peak->Sequence.startTimeBin - 1];
1283  }
1284 
1285  if(Peak->Timebin + 1 <
1286  Peak->Sequence.startTimeBin + Peak->Sequence.Length)
1287  {
1288  iADCPlus = Peak->Sequence.FirstAdc[Peak->Timebin -
1289  Peak->Sequence.startTimeBin + 1];
1290  }
1291 
1292  /* check if Gaussfit will fail */
1293  if((iADCValue == 0) || (iADCPlus == 0) || (iADCMinus == 0) ||
1294  ((iADCValue <= iADCPlus) && (iADCValue <= iADCMinus)))
1295  {
1296  /* use weighted mean instead */
1297  iUseGauss = 0;
1298  }
1299  else
1300  {
1301  /* do gaussfit */
1302  Peak->TimeSigma = sqrt (1 /
1303  ((2 * fastlog[iADCValue]) -
1304  (fastlog[iADCPlus] + fastlog[iADCMinus])));
1305  Peak->TimePosition = (float) Peak->Timebin +
1306  mDb->timeOffset(iSec*mDb->numberOfPads()+Peak->pad,iRow) +
1307  sqr(Peak->TimeSigma) * (fastlog[iADCPlus] - fastlog[iADCMinus]);
1308  }
1309  }
1310 
1311  if((iUseGauss & 2) == 0)
1312  {
1313  /* calculate time position by weighted mean */
1314  Peak->TimePosition = (float) TimeSum / (float) ChargeSum;
1315 
1316  /* calculate pad direction cluster width by weighted mean */
1317  Peak->TimeSigma=0;
1318  for(iSequence = 0; iSequence < Cluster->NumSequences; iSequence++)
1319  {
1320  for(iBin = 0; iBin < Cluster->Sequence[iSequence].Length; iBin++)
1321  {
1322  Peak->TimeSigma+=(sqr((float)
1323  Cluster->Sequence[iSequence].startTimeBin
1324  + (float) iBin - Peak->TimePosition)
1325  *Cluster->Sequence[iSequence].FirstAdc[iBin]);
1326  }
1327  }
1328  Peak->TimeSigma /= (float) ChargeSum;
1329  }
1330 
1331  /* transform from pad/time to x/y/z */
1332  if(padtrans(Peak, iRow, iSec,
1333  pRadius, pDeflection))
1334  {
1335  if (Peak->x == 0. && Peak->y == 0.) {
1336  // This if-statement can be deleted as soon as the slow simulator is fixed. This also occurs for FTPC DAQ data.
1337  LOG_WARN << "Hit rejected because of an error in the FTPC data. (x, y, z) = (0. ,0., z)" << endm;
1338  }
1339 
1340 /*
1341  if(!std::isnan(Peak->x) && !std::isnan(Peak->y) && !std::isnan(Peak->PadSigma) &&
1342 !std::isnan(Peak->TimeSigma) // && Peak->PeakHeight>=mParam->minimumClusterMaxADC())
1343  && Peak->Rad <= mDb->sensitiveVolumeOuterRadius() && Peak->Rad >= mDb->sensitiveVolumeInnerRadius() )
1344 */
1345  if(!std::isnan(Peak->x) && !std::isnan(Peak->y) && !std::isnan(Peak->PadSigma) && !std::isnan(Peak->TimeSigma))
1346  //&& (Cluster->EndPad +1 - Cluster->StartPad)<=MAXPADLENGTH
1347  //&& Peak->Sequence.Length<=MAXTIMELENGTH )
1348 
1349  {
1350  // create new point
1351 
1352  // fill QA histograms - unless BFC option "noHists" selected
1353  if (mhpad) mhpad->Fill(Cluster->EndPad +1 - Cluster->StartPad,1);
1354  if (mhtime) mhtime->Fill(Peak->Sequence.Length,1);
1355 
1356  Int_t numPoint = mPoint->GetEntriesFast();
1357  if (numPoint >= mPoint->GetSize()) mPoint->Expand(mPoint->GetSize()+5000);
1358 
1359  mPoint->AddAt(new StFtpcConfMapPoint(), numPoint);
1360  StFtpcConfMapPoint *thispoint = (StFtpcConfMapPoint *) mPoint->At(numPoint);
1361 
1362 
1363  if (iRow>=10 && mDb->SwapRDO6RDO7East()) {
1364  iSecSAVE = iSec;
1365  iRowSAVE = iRow;
1366  if(iSec==3&&iRow==10){
1367  iRow = 16;
1368  }
1369  else if(iSec==3&&iRow==11){
1370  iRow = 17;
1371  }
1372  else if(iSec==3&&iRow==12){
1373  iRow = 18;
1374  }
1375  else if(iSec==3&&iRow==13){
1376  iRow = 19;
1377  }
1378  else if(iSec==3&&iRow==14){
1379  iRow = 18;
1380  iSec = 4;
1381  }
1382  else if(iSec==3&&iRow==15){
1383  iRow = 19;
1384  iSec = 4;
1385  }
1386  else if(iSec==3&&iRow==16){
1387  iRow = 10;
1388  }
1389  else if(iSec==3&&iRow==17) {
1390  iRow = 11;
1391  }
1392  else if(iSec==3&&iRow==18){
1393  iRow = 12;
1394  }
1395  else if(iSec==3&&iRow==19) {
1396  iRow = 13;
1397  }
1398  else if(iSec==4&&iRow==18){
1399  iRow = 14;
1400  iSec = 3;
1401  }
1402  else if(iSec==4&&iRow==19){
1403  iRow = 15;
1404  iSec = 3;
1405  }
1406  thispoint->SetPadRow(iRow+1);
1407  thispoint->SetSector(iSec+1);
1408  iSec = iSecSAVE;
1409  iRow = iRowSAVE;
1410  }
1411  else {
1412  thispoint->SetPadRow(iRow+1);
1413  thispoint->SetSector(iSec+1);
1414  }
1415 
1416  thispoint->SetNumberPads(Cluster->EndPad +1 - Cluster->StartPad);
1417  thispoint->SetNumberBins(Peak->Sequence.Length);
1418  thispoint->SetMaxADC((long)Peak->PeakHeight);
1419  thispoint->SetCharge(ChargeSum);
1420  thispoint->SetPadPos(Peak->PadPosition);
1421  thispoint->SetTimePos(Peak->TimePosition);
1422  thispoint->SetPadPosSigma(Peak->PadSigma);
1423  thispoint->SetTimePosSigma(Peak->TimeSigma);
1424  thispoint->SetX(Peak->x);
1425  thispoint->SetY(Peak->y);
1426  thispoint->SetZ(Peak->z);
1427  thispoint->SetSigmaPhi(Peak->PadSigma*mDb->radiansPerPad());
1428  thispoint->SetSigmaR(Peak->TimeSigma*Peak->Rad/Peak->TimePosition);
1429  fDriftLength = mDb->sensitiveVolumeOuterRadius() - Peak->Rad;
1430  fPhiError = mParam->padDiffusionErrors(0)
1431  + fDriftLength*mParam->padDiffusionErrors(1)
1432  + fDriftLength*fDriftLength*mParam->padDiffusionErrors(2);
1433  fRadError = mParam->timeDiffusionErrors(0)
1434  + fDriftLength*mParam->timeDiffusionErrors(1)
1435  + fDriftLength*fDriftLength*mParam->timeDiffusionErrors(2);
1436  if(thispoint->GetNumberPads()==2)
1437  {
1438  fPhiError = ::sqrt(fPhiError * fPhiError
1439  + sqr(mParam->twoPadWeightedError()));
1440  }
1441  if((thispoint->GetNumberPads()==3) && ((iUseGauss & 1) == 1))
1442  {
1443  fPhiError = ::sqrt(fPhiError * fPhiError
1444  + sqr(mParam->threePadGaussError()));
1445  }
1446  if((thispoint->GetNumberPads()==3) && ((iUseGauss & 1) == 0))
1447  {
1448  fPhiError = ::sqrt(fPhiError * fPhiError
1449  + sqr(mParam->threePadWeightedError()));
1450  }
1451 
1452  thispoint->SetFlags(0);
1453 
1454  if(iADCValue>254)
1455  {
1456  thispoint->SetFlags(4);
1457  fPhiError = ::sqrt(fPhiError * fPhiError
1458  + sqr(mParam->padSaturatedClusterError()));
1459  fRadError = ::sqrt(fRadError * fRadError
1460  + sqr(mParam->timeSaturatedClusterError()));
1461  }
1462  if(Cluster->CutOff==1)
1463  {
1464  thispoint->SetFlags(thispoint->GetFlags() | 16);
1465  fPhiError = ::sqrt(fPhiError * fPhiError
1466  + sqr(mParam->padCutoffClusterError()));
1467  fRadError = ::sqrt(fRadError * fRadError
1468  + sqr(mParam->timeCutoffClusterError()));
1469  }
1470 
1471  if (Peak->Rad > mDb->sensitiveVolumeOuterRadius() ||
1472  Peak->Rad < mDb->sensitiveVolumeInnerRadius()
1473  /* please add additional criteria here*/) {
1474  // don't use this point for tracking
1475 
1476  thispoint->SetFlags(thispoint->GetFlags() | 128);
1477  }
1478 
1479  /* transform errors to actual hit position */
1480  PadtransBin=(int) ((Peak->TimePosition+0.5)*PadtransPerTimebin);
1481  fPhiError *= Peak->Rad / mDb->sensitiveVolumeOuterRadius();
1482  fRadError *= (pRadius[10*PadtransBin]-pRadius[10*PadtransBin+10])
1483  / (pRadius[10]-pRadius[20]);
1484 
1485 if (mcldebug){
1486  mcldebug->fillclustertree(Peak,Cluster,ChargeSum,iHardSec,iHardRow,fRadError,fPhiError,0,0,iNumPeaks);
1487 }
1488 
1489  thispoint->SetXerr(::sqrt(fRadError*cos(Peak->Phi)
1490  *fRadError*cos(Peak->Phi)
1491  + fPhiError*sin(Peak->Phi)
1492  *fPhiError*sin(Peak->Phi)));
1493  thispoint->SetYerr(::sqrt(fRadError*sin(Peak->Phi)
1494  *fRadError*sin(Peak->Phi)
1495  + fPhiError*cos(Peak->Phi)
1496  *fPhiError*cos(Peak->Phi)));
1497  thispoint->SetZerr(mParam->zDirectionError());
1498 
1499  }
1500  else
1501  {
1502  if (DebugOn) {
1503  LOG_DEBUG<<"Cluster fitting error. Point not stored"<<endm;
1504  }
1505  }
1506  }
1507  else{
1508  if (DebugOn) {
1509  LOG_DEBUG<<"Peak position can't be transformed"<<endm;
1510  }
1511  }
1512  } /* end of: if(iNumPeaks == 1) */
1513  else
1514  {
1515  /* multi-peak cluster, unfold */
1516  /* initialize all peaks for unfolding */
1517  BadFit=0;
1518  FailedFit=0;
1519  for(iPeakIndex=0; iPeakIndex < iNumPeaks; iPeakIndex++)
1520  {
1521  Peak[iPeakIndex].TimePosition = (float) Peak[iPeakIndex].Timebin;
1522  Peak[iPeakIndex].Timebin_saved = Peak[iPeakIndex].Timebin;
1523  Peak[iPeakIndex].PadPosition = (float) Peak[iPeakIndex].pad;
1524  Peak[iPeakIndex].pad_saved = Peak[iPeakIndex].pad;
1525  Peak[iPeakIndex].OldTimePosition = 0;
1526  Peak[iPeakIndex].OldPadPosition = 0;
1527  Peak[iPeakIndex].OldPeakHeight = 0;
1528  Peak[iPeakIndex].TimeSigma = sigmat((float) Peak[iPeakIndex].Timebin);
1529  Peak[iPeakIndex].PadSigma = sigmax((float) Peak[iPeakIndex].Timebin);
1530  Peak[iPeakIndex].PeakHeight =
1531  (float)Peak[iPeakIndex].Sequence
1532  .FirstAdc[Peak[iPeakIndex].Timebin
1533  - Peak[iPeakIndex].Sequence.startTimeBin];
1534  }
1535 
1536  /* approximation loop */
1537  iNumUnfoldLoops=0;
1538  do
1539  { /* beginning of: while((SumDeltaPos > 0.01... */
1540  SumDeltaPos=0;
1541  SumDeltaAmp=0;
1542  PeakHeightSum=0;
1543 
1544  /* loop over all peaks in cluster */
1545  for(iPeakIndex=0; iPeakIndex < iNumPeaks; iPeakIndex++)
1546  {
1547 
1548  do
1549  {
1550 
1551  PeakShifted=0;
1552  /* get values for pad direction gaussfit */
1553  iADCValue=0;
1554  iADCPlus=0;
1555  iADCMinus=0;
1556  iADCTimePlus=0;
1557  iADCTimeMinus=0;
1558  /* find sequences on the right pads */
1559  for(iSequence = 0; iSequence < Cluster->NumSequences; iSequence++)
1560  {
1561  if((Cluster->SequencePad[iSequence]
1562  == Peak[iPeakIndex].pad) &&
1563  (Cluster->Sequence[iSequence].startTimeBin
1564  <= Peak[iPeakIndex].Timebin) &&
1565  (Cluster->Sequence[iSequence].startTimeBin +
1566  Cluster->Sequence[iSequence].Length
1567  > Peak[iPeakIndex].Timebin))
1568  {
1569  iADCValue = Cluster->Sequence[iSequence]
1570  .FirstAdc[Peak[iPeakIndex].Timebin
1571  - Cluster->Sequence[iSequence].startTimeBin];
1572  /* get values for time direction gaussfit */
1573  if(Peak[iPeakIndex].Timebin >
1574  Cluster->Sequence[iSequence].startTimeBin)
1575  {
1576  iADCTimeMinus = Cluster->Sequence[iSequence]
1577  .FirstAdc[Peak[iPeakIndex].Timebin
1578  - Cluster->Sequence[iSequence].startTimeBin - 1];
1579  }
1580 
1581  if(Peak[iPeakIndex].Timebin + 1 <
1582  Cluster->Sequence[iSequence].startTimeBin
1583  + Cluster->Sequence[iSequence].Length)
1584  {
1585  iADCTimePlus = Cluster->Sequence[iSequence]
1586  .FirstAdc[Peak[iPeakIndex].Timebin
1587  - Cluster->Sequence[iSequence].startTimeBin + 1];
1588  }
1589 
1590  }
1591  if((Cluster->SequencePad[iSequence]
1592  == Peak[iPeakIndex].pad -1) &&
1593  (Cluster->Sequence[iSequence].startTimeBin
1594  <= Peak[iPeakIndex].Timebin) &&
1595  (Cluster->Sequence[iSequence].startTimeBin +
1596  Cluster->Sequence[iSequence].Length
1597  > Peak[iPeakIndex].Timebin))
1598  {
1599  iADCMinus = Cluster->Sequence[iSequence]
1600  .FirstAdc[Peak[iPeakIndex].Timebin
1601  - Cluster->Sequence[iSequence].startTimeBin];
1602  }
1603  if((Cluster->SequencePad[iSequence]
1604  == Peak[iPeakIndex].pad + 1) &&
1605  (Cluster->Sequence[iSequence].startTimeBin
1606  <= Peak[iPeakIndex].Timebin) &&
1607  (Cluster->Sequence[iSequence].startTimeBin +
1608  Cluster->Sequence[iSequence].Length
1609  > Peak[iPeakIndex].Timebin))
1610  {
1611  iADCPlus = Cluster->Sequence[iSequence]
1612  .FirstAdc[Peak[iPeakIndex].Timebin
1613  - Cluster->Sequence[iSequence].startTimeBin];
1614  }
1615  }
1616 
1617  /* calculate influence of other peaks */
1618  fDeltaADC=0;
1619 
1620  fDeltaADCPlus=0;
1621  fDeltaADCMinus=0;
1622  fDeltaADCTimePlus=0;
1623  fDeltaADCTimeMinus=0;
1624  for(iInnerIndex=0; iInnerIndex < iNumPeaks; iInnerIndex++)
1625  {
1626  if(iInnerIndex != iPeakIndex)
1627  {
1628  /* check if peaks are overlaid */
1629  if((Peak[iPeakIndex].pad == Peak[iInnerIndex].pad)
1630  && (Peak[iPeakIndex].Timebin == Peak[iInnerIndex].Timebin))
1631  {
1632  /* one peak has been shifted on top of the other */
1633  /* that kills the unfold procedure */
1634  FailedFit=1;
1635  break;
1636  }
1637 
1638  fDeltaADC+=gauss_2d(Peak[iPeakIndex].pad,
1639  Peak[iPeakIndex].Timebin,
1640  Peak[iInnerIndex].PadPosition,
1641  Peak[iInnerIndex].PadSigma,
1642  Peak[iInnerIndex].TimePosition,
1643  Peak[iInnerIndex].TimeSigma,
1644  Peak[iInnerIndex].PeakHeight);
1645  fDeltaADCPlus+=gauss_2d(Peak[iPeakIndex].pad+1,
1646  Peak[iPeakIndex].Timebin,
1647  Peak[iInnerIndex].PadPosition,
1648  Peak[iInnerIndex].PadSigma,
1649  Peak[iInnerIndex].TimePosition,
1650  Peak[iInnerIndex].TimeSigma,
1651  Peak[iInnerIndex].PeakHeight);
1652  fDeltaADCMinus+=gauss_2d(Peak[iPeakIndex].pad-1,
1653  Peak[iPeakIndex].Timebin,
1654  Peak[iInnerIndex].PadPosition,
1655  Peak[iInnerIndex].PadSigma,
1656  Peak[iInnerIndex].TimePosition,
1657  Peak[iInnerIndex].TimeSigma,
1658  Peak[iInnerIndex].PeakHeight);
1659  fDeltaADCTimePlus+=gauss_2d(Peak[iPeakIndex].pad,
1660  Peak[iPeakIndex].Timebin+1,
1661  Peak[iInnerIndex].PadPosition,
1662  Peak[iInnerIndex].PadSigma,
1663  Peak[iInnerIndex].TimePosition,
1664  Peak[iInnerIndex].TimeSigma,
1665  Peak[iInnerIndex].PeakHeight);
1666  fDeltaADCTimeMinus+=gauss_2d(Peak[iPeakIndex].pad,
1667  Peak[iPeakIndex].Timebin-1,
1668  Peak[iInnerIndex].PadPosition,
1669  Peak[iInnerIndex].PadSigma,
1670  Peak[iInnerIndex].TimePosition,
1671  Peak[iInnerIndex].TimeSigma,
1672  Peak[iInnerIndex].PeakHeight);
1673  }
1674  }
1675 
1676  /* subtract influence of other peaks */
1677  iADCValue-=(int) (fDeltaADC+0.5);
1678  iADCPlus-=(int) (fDeltaADCPlus+0.5);
1679  iADCMinus-=(int) (fDeltaADCMinus+0.5);
1680  iADCTimePlus-=(int) (fDeltaADCTimePlus+0.5);
1681  iADCTimeMinus-=(int) (fDeltaADCTimeMinus+0.5);
1682 
1683  /* check if Gaussfit is possible */
1684 
1685  /* introduce minimal gradient if none is there */
1686  if(iADCValue == iADCPlus)
1687  {
1688  iADCPlus--;
1689  }
1690  if(iADCValue == iADCMinus)
1691  {
1692  iADCMinus--;
1693  }
1694  if(iADCValue == iADCTimePlus)
1695  {
1696  iADCTimePlus--;
1697  }
1698  if(iADCValue == iADCTimeMinus)
1699  {
1700  iADCTimeMinus--;
1701  }
1702 
1703 
1704  /* set all 0s to 1 to keep gauss from crashing */
1705  /* and iADCValue=2 if iADCValue <= 1 (09/28/98) */
1706  if(iADCValue <= 1)
1707  {
1708  iADCValue=2;
1709  BadFit=1;
1710  }
1711  if(iADCPlus <= 0)
1712  {
1713  iADCPlus=1;
1714  BadFit=1;
1715  }
1716  if(iADCMinus <= 0)
1717  {
1718  iADCMinus=1;
1719  BadFit=1;
1720  }
1721  if(iADCTimePlus <= 0)
1722  {
1723  iADCTimePlus=1;
1724  BadFit=1;
1725  }
1726  if(iADCTimeMinus <= 0)
1727  {
1728  iADCTimeMinus=1;
1729  BadFit=1;
1730  }
1731 
1732 /* move here to protect against iADCValue=iADCPlus=iADCMinus which
1733  causes Peak[..].PadSigma=Inf 09.30.98 HH,JCS */
1734  if(FailedFit == 1)
1735  {
1736  break;
1737  }
1738 
1739  /* shift if peak has moved after correction */
1740  if(iADCValue < iADCPlus)
1741  {
1742  Peak[iPeakIndex].pad++;
1743  PeakShifted=1;
1744  }
1745  if(iADCValue < iADCMinus && PeakShifted==0)
1746  {
1747  Peak[iPeakIndex].pad--;
1748  PeakShifted=1;
1749  }
1750  if(iADCValue < iADCTimePlus && PeakShifted==0)
1751  {
1752  Peak[iPeakIndex].Timebin++;
1753  PeakShifted=1;
1754  }
1755  if(iADCValue < iADCTimeMinus && PeakShifted==0)
1756  {
1757  Peak[iPeakIndex].Timebin--;
1758  PeakShifted=1;
1759  }
1760  }
1761  while(PeakShifted>0 && PeakShifted<5);
1762 
1763  /* recalculate peak position */
1764  NewPeakHeight=iADCValue;
1765  double qwe;
1766  qwe = (2 * fastlog[iADCValue]) -(fastlog[iADCPlus] + fastlog[iADCMinus]);
1767  if (qwe<=0.) continue; //skip case for non physical value (VP)
1768  Peak[iPeakIndex].PadSigma = sqrt (1 / qwe);
1769  NewPadPosition = (float) Peak[iPeakIndex].pad +
1770  sqr(Peak[iPeakIndex].PadSigma) *
1771  (fastlog[iADCPlus] - fastlog[iADCMinus]);
1772  qwe = (2 * fastlog[iADCValue]) -(fastlog[iADCTimePlus] + fastlog[iADCTimeMinus]);
1773  if (qwe<=0.) continue; //skip case for non physical value (VP)
1774  Peak[iPeakIndex].TimeSigma = sqrt (1 / qwe);
1775 
1776  NewTimePosition = (float) Peak[iPeakIndex].Timebin +
1777  mDb->timeOffset(iSec*mDb->numberOfPads()
1778  +Peak[iPeakIndex].pad,iRow) +
1779  sqr(Peak[iPeakIndex].TimeSigma) *
1780  (fastlog[iADCTimePlus] - fastlog[iADCTimeMinus]);
1781 
1782  /* introduce inertia to break infinite loops */
1783  if(iNumUnfoldLoops > MAXFASTLOOPS)
1784  {
1785  NewPeakHeight=(NewPeakHeight+Peak[iPeakIndex].PeakHeight)/2;
1786  NewPadPosition=(NewPadPosition+Peak[iPeakIndex].PadPosition)/2;
1787  NewTimePosition=(NewTimePosition+Peak[iPeakIndex].TimePosition)/2;
1788  }
1789 
1790  SumDeltaAmp+=fabs(NewPeakHeight - Peak[iPeakIndex].PeakHeight);
1791  Peak[iPeakIndex].OldPeakHeight=Peak[iPeakIndex].PeakHeight;
1792  Peak[iPeakIndex].PeakHeight=NewPeakHeight;
1793  SumDeltaPos+=fabs(NewPadPosition - Peak[iPeakIndex].PadPosition);
1794  Peak[iPeakIndex].OldPadPosition=Peak[iPeakIndex].PadPosition;
1795  Peak[iPeakIndex].PadPosition=NewPadPosition;
1796  SumDeltaPos+=fabs(NewTimePosition - Peak[iPeakIndex].TimePosition);
1797  Peak[iPeakIndex].OldTimePosition=Peak[iPeakIndex].TimePosition;
1798  Peak[iPeakIndex].TimePosition=NewTimePosition;
1799 
1800  PeakHeightSum+=NewPeakHeight;
1801 
1802  } /* end of: for(iPeakIndex=0;... */
1803  iNumUnfoldLoops++;
1804  }
1805  while((SumDeltaPos > UNFOLDLIMIT || SumDeltaAmp > 1)
1806  && (iNumUnfoldLoops < MAXLOOPS || SumDeltaPos > UNFOLDFAILEDLIMIT)
1807  && iNumUnfoldLoops < MAXLOOPS+1 && FailedFit == 0);
1808  /* try unfold loop once more if error > failedlimit */
1809 
1810  if(SumDeltaPos > UNFOLDFAILEDLIMIT || FailedFit == 1)
1811  {
1812  /* unfolding has failed, use original peak positions */
1813 
1814  for(iPeakIndex=0; iPeakIndex < iNumPeaks; iPeakIndex++)
1815  {
1816  Peak[iPeakIndex].TimePosition = (float) Peak[iPeakIndex].Timebin_saved;
1817  Peak[iPeakIndex].Timebin = Peak[iPeakIndex].Timebin_saved;
1818  Peak[iPeakIndex].PadPosition = (float) Peak[iPeakIndex].pad_saved;
1819  Peak[iPeakIndex].pad = Peak[iPeakIndex].pad_saved;
1820  for(iSequence = 0; iSequence < Cluster->NumSequences; iSequence++)
1821  {
1822  if((Cluster->SequencePad[iSequence]
1823  == Peak[iPeakIndex].pad) &&
1824  (Cluster->Sequence[iSequence].startTimeBin
1825  <= Peak[iPeakIndex].Timebin) &&
1826  (Cluster->Sequence[iSequence].startTimeBin +
1827  Cluster->Sequence[iSequence].Length
1828  > Peak[iPeakIndex].Timebin))
1829  {
1830  Peak[iPeakIndex].PeakHeight = Cluster->Sequence[iSequence]
1831  .FirstAdc[Peak[iPeakIndex].Timebin
1832  - Cluster->Sequence[iSequence].startTimeBin];
1833  }
1834  }
1835  }
1836  if (DebugOn) {
1837  LOG_DEBUG<<"unfold failed!"<<endm;
1838  }
1839  } /* end of: if(SumDeltaPos > UNFOLDFAILEDLIMIT ... */
1840 
1841  /* write unfolded peaks to point table */
1842  for(iPeakIndex=0; iPeakIndex < iNumPeaks; iPeakIndex++)
1843  {
1844  /* transform from pad/time to x/y/z */
1845  if(padtrans(&(Peak[iPeakIndex]), iRow, iSec,
1846  pRadius, pDeflection))
1847  {
1848  if (Peak[iPeakIndex].x == 0. && Peak[iPeakIndex].y == 0.) {
1849  // This if-statement can be deleted as soon as the slow simulator is fixed. This also occurs for FTPC DAQ data.
1850  LOG_WARN << "Hit rejected because of an error in the FTPC data. (x, y, z) = (0. ,0., z)" << endm;
1851  }
1852 
1853  /* in very complicated clusters some hits may have been unfolded
1854  with errors while the rest of the cluster is okay, don't fill
1855  these hits into array: */
1856 /*
1857  if(!std::isnan(Peak[iPeakIndex].x) && !std::isnan(Peak[iPeakIndex].y) &&
1858 !std::isnan(Peak[iPeakIndex].PadSigma) && !std::isnan(Peak[iPeakIndex].TimeSigma) // && Peak[iPeakIndex].PeakHeight>=mParam->minimumClusterMaxADC())
1859  && Peak[iPeakIndex].Rad <= mDb->sensitiveVolumeOuterRadius()
1860  && Peak[iPeakIndex].Rad >= mDb->sensitiveVolumeInnerRadius() )
1861 */
1862  if(!std::isnan(Peak[iPeakIndex].x) && !std::isnan(Peak[iPeakIndex].y)
1863  && !std::isnan(Peak[iPeakIndex].PadSigma) && !std::isnan(Peak[iPeakIndex].TimeSigma))
1864  //&& (Cluster->EndPad +1 - Cluster->StartPad)<=MAXPADLENGTH
1865  //&& Peak[iPeakIndex].Sequence.Length<=MAXTIMELENGTH)
1866  {
1867 
1868  // fill QA histograms - unless BFC option "noHists" selected
1869  if (mhpad) mhpad->Fill(Cluster->EndPad +1 - Cluster->StartPad,iNumPeaks);
1870  if (mhtime) mhtime->Fill(Peak[iPeakIndex].Sequence.Length,iNumPeaks);
1871 
1872  // create new point
1873  Int_t numPoint = mPoint->GetEntriesFast();
1874  if (numPoint >= mPoint->GetSize()) mPoint->Expand(mPoint->GetSize()+5000);
1875 
1876  mPoint->AddAt(new StFtpcConfMapPoint(), numPoint);
1877  StFtpcConfMapPoint *thispoint = (StFtpcConfMapPoint *) mPoint->At(numPoint);
1878 
1879  if (iRow>=10 && mDb->SwapRDO6RDO7East()) {
1880  iSecSAVE = iSec;
1881  iRowSAVE = iRow;
1882  if(iSec==3&&iRow==10){
1883  iRow = 16;
1884  }
1885  else if(iSec==3&&iRow==11){
1886  iRow = 17;
1887  }
1888  else if(iSec==3&&iRow==12){
1889  iRow = 18;
1890  }
1891  else if(iSec==3&&iRow==13){
1892  iRow = 19;
1893  }
1894  else if(iSec==3&&iRow==14){
1895  iRow = 18;
1896  iSec = 4;
1897  }
1898  else if(iSec==3&&iRow==15){
1899  iRow = 19;
1900  iSec = 4;
1901  }
1902  else if(iSec==3&&iRow==16){
1903  iRow = 10;
1904  }
1905  else if(iSec==3&&iRow==17){
1906  iRow = 11;
1907  }
1908  else if(iSec==3&&iRow==18){
1909  iRow = 12;
1910  }
1911  else if(iSec==3&&iRow==19){
1912  iRow = 13;
1913  }
1914  else if(iSec==4&&iRow==18){
1915  iRow = 14;
1916  iSec = 3;
1917  }
1918  else if(iSec==4&&iRow==19){
1919  iRow = 15;
1920  iSec = 3;
1921  }
1922  thispoint->SetPadRow(iRow+1);
1923  thispoint->SetSector(iSec+1);
1924  iSec = iSecSAVE;
1925  iRow = iRowSAVE;
1926  }
1927  else {
1928  thispoint->SetPadRow(iRow+1);
1929  thispoint->SetSector(iSec+1);
1930  }
1931 
1932  thispoint->SetNumberPads(Cluster->EndPad +1 - Cluster->StartPad);
1933  thispoint->SetNumberBins(Peak[iPeakIndex].Sequence.Length);
1934  thispoint->SetMaxADC((long)Peak[iPeakIndex].PeakHeight);
1935  thispoint->SetCharge((Long_t)(ChargeSum*(Peak[iPeakIndex].PeakHeight
1936  /PeakHeightSum)));
1937  thispoint->SetPadPos(Peak[iPeakIndex].PadPosition);
1938  thispoint->SetTimePos(Peak[iPeakIndex].TimePosition);
1939  thispoint->SetPadPosSigma(Peak[iPeakIndex].PadSigma);
1940  thispoint->SetTimePosSigma(Peak[iPeakIndex].TimeSigma);
1941  thispoint->SetX(Peak[iPeakIndex].x);
1942  thispoint->SetY(Peak[iPeakIndex].y);
1943  thispoint->SetZ(Peak[iPeakIndex].z);
1944  thispoint->SetSigmaPhi(Peak[iPeakIndex].PadSigma
1945  *mDb->radiansPerPad());
1946  thispoint->SetSigmaR(Peak[iPeakIndex].TimeSigma
1947  * Peak[iPeakIndex].Rad
1948  /Peak[iPeakIndex].TimePosition);
1949 
1950  fDriftLength = mDb->sensitiveVolumeOuterRadius()
1951  - Peak[iPeakIndex].Rad;
1952  fPhiError = mParam->padDiffusionErrors(0)
1953  + fDriftLength*mParam->padDiffusionErrors(1)
1954  + fDriftLength*fDriftLength*mParam->padDiffusionErrors(2);
1955  fRadError = mParam->timeDiffusionErrors(0)
1956  + fDriftLength*mParam->timeDiffusionErrors(1)
1957  + fDriftLength*fDriftLength*mParam->timeDiffusionErrors(2);
1958 
1959  /* clusters are unfolded */
1960  thispoint->SetFlags(1);
1961  fPhiError = ::sqrt(fPhiError * fPhiError
1962  + sqr(mParam->padUnfoldError()));
1963  fRadError = ::sqrt(fRadError * fRadError
1964  + sqr(mParam->timeUnfoldError()));
1965 
1966  if(iADCValue>254)
1967  {
1968  thispoint->SetFlags(5);
1969  fPhiError = ::sqrt(fPhiError * fPhiError
1970  + sqr(mParam->padSaturatedClusterError()));
1971  fRadError = ::sqrt(fRadError * fRadError
1972  + sqr(mParam->timeSaturatedClusterError()));
1973  }
1974  if(BadFit==1)
1975  {
1976  thispoint->SetFlags(thispoint->GetFlags() | 8);
1977  fPhiError = ::sqrt(fPhiError * fPhiError
1978  + sqr(mParam->padBadFitError()));
1979  fRadError = ::sqrt(fRadError * fRadError
1980  + sqr(mParam->timeBadFitError()));
1981  }
1982  if(iNumUnfoldLoops == MAXLOOPS)
1983  {
1984  thispoint->SetFlags(thispoint->GetFlags() | 10);
1985  fPhiError = ::sqrt(fPhiError * fPhiError
1986  + sqr(mParam->padFailedFitError()));
1987  fRadError = ::sqrt(fRadError * fRadError
1988  + sqr(mParam->timeFailedFitError()));
1989  }
1990  if(Cluster->CutOff==1)
1991  {
1992  thispoint->SetFlags(thispoint->GetFlags() | 16);
1993  fPhiError = ::sqrt(fPhiError * fPhiError
1994  + sqr(mParam->padCutoffClusterError()));
1995  fRadError = ::sqrt(fRadError * fRadError
1996  + sqr(mParam->timeCutoffClusterError()));
1997  }
1998 
1999  if (Peak[iPeakIndex].Rad > mDb->sensitiveVolumeOuterRadius() ||
2000  Peak[iPeakIndex].Rad < mDb->sensitiveVolumeInnerRadius()
2001  /* please add additional criteria here*/) {
2002  // don't use this point for tracking
2003 
2004  thispoint->SetFlags(thispoint->GetFlags() | 128);
2005  }
2006 
2007  /* transform errors to actual hit position */
2008  PadtransBin=(int) ((Peak->TimePosition+0.5)*PadtransPerTimebin);
2009  fPhiError *= Peak->Rad / mDb->sensitiveVolumeOuterRadius();
2010  fRadError *= (pRadius[10*PadtransBin]-pRadius[10*PadtransBin+10])
2011  / (pRadius[10]-pRadius[20]);
2012 
2013 if (mcldebug){
2014  mcldebug->fillclustertree(Peak[iPeakIndex],Cluster,ChargeSum*Peak[iPeakIndex].PeakHeight/PeakHeightSum,iHardSec,iHardRow,fRadError,fPhiError,10,0,iNumPeaks);
2015 }
2016 
2017  thispoint->SetXerr(::sqrt(fRadError*cos(Peak[iPeakIndex].Phi)
2018  *fRadError*cos(Peak[iPeakIndex].Phi)
2019  + fPhiError*sin(Peak[iPeakIndex].Phi)
2020  *fPhiError*sin(Peak[iPeakIndex].Phi)));
2021  thispoint->SetYerr(::sqrt(fRadError*sin(Peak[iPeakIndex].Phi)
2022  *fRadError*sin(Peak[iPeakIndex].Phi)
2023  + fPhiError*cos(Peak[iPeakIndex].Phi)
2024  *fPhiError*cos(Peak[iPeakIndex].Phi)));
2025  thispoint->SetZerr(mParam->zDirectionError());
2026 
2027  }
2028  }
2029  else
2030  {
2031  if (DebugOn) {
2032  LOG_DEBUG<<"Peak position can't be transformed!"<<endm;
2033  }
2034  }
2035  } /* end of: for(iPeakIndex=0;... */
2036  } /*end of: if(iNumPeaks==1) ... else { */
2037  return TRUE;
2038 }
2039 
2040 int StFtpcClusterFinder::padtrans(TPeak *Peak,
2041  int iRow,
2042  int iSec,
2043  double *pRadius,
2044  double *pDeflection)
2045 {
2046  int PadtransPerTimebin;
2047  int PadtransLower;
2048  float PhiDeflect, TimeCoordinate;
2049 
2050  if (iRow>=10 && mDb->SwapRDO6RDO7East()) {
2051  if(iSec==3&&iRow==10) {
2052  iRow=16;
2053  }
2054  else if(iSec==3&&iRow==11) {
2055  iRow=17;
2056  }
2057  else if(iSec==3&&iRow==12) {
2058  iRow=18;
2059  }
2060  else if(iSec==3&&iRow==13) {
2061  iRow=19;
2062  }
2063  else if(iSec==3&&iRow==14) {
2064  iSec=4;
2065  iRow=18;
2066  }
2067  else if(iSec==3&&iRow==15) {
2068  iSec=4;
2069  iRow=19;
2070  }
2071  else if(iSec==3&&iRow==16) {
2072  iRow=10;
2073  }
2074  else if(iSec==3&&iRow==17) {
2075  iRow=11;
2076  }
2077  else if(iSec==3&&iRow==18) {
2078  iRow=12;
2079  }
2080  else if(iSec==3&&iRow==19) {
2081  iRow=13;
2082  }
2083  else if(iSec==4&&iRow==18) {
2084  iSec=3;
2085  iRow=14;
2086  }
2087  else if(iSec==4&&iRow==19) {
2088  iSec=3;
2089  iRow=15;
2090  }
2091  }
2092 
2093  /* preparatory calculations */
2094  TimeCoordinate = Peak->TimePosition + 0.5; /*time start at beginning of bin 0*/
2095  // include tZero = time from collision to beginning of bin 0
2096  TimeCoordinate += mDb->tZero()/mDb->microsecondsPerTimebin();
2097  PadtransPerTimebin = (int) mParam->numberOfDriftSteps() / mDb->numberOfTimebins();
2098  PadtransLower= (int) (TimeCoordinate*PadtransPerTimebin);
2099 
2100  if ( TimeCoordinate > mDb->numberOfTimebins() || TimeCoordinate <= 0 )
2101  {
2102  // exceeds table dimensions
2103  return FALSE;
2104  }
2105 
2106  /* linear interpolation in radius table */
2107  Peak->Rad=pRadius[iRow + mDb->numberOfPadrowsPerSide() * PadtransLower]
2108  -(pRadius[iRow + mDb->numberOfPadrowsPerSide()*(PadtransLower)]
2109  -pRadius[iRow + mDb->numberOfPadrowsPerSide() * (PadtransLower+1)])/2;
2110 
2111  if ( pRadius[iRow + mDb->numberOfPadrowsPerSide() * PadtransLower] == 0 )
2112  {
2113  // outside FTPC sensitive region
2114  return FALSE;
2115  }
2116 
2117  /* linear interpolation in deflection table */
2118  PhiDeflect=pDeflection[iRow + mDb->numberOfPadrowsPerSide() * PadtransLower]
2119  +(pDeflection[iRow + mDb->numberOfPadrowsPerSide() * (PadtransLower+1)]
2120  -pDeflection[iRow + mDb->numberOfPadrowsPerSide() * PadtransLower])/2;
2121 
2122 
2123  /* calculate phi angle from pad position */
2124  // for FTPC West
2125  if (iRow <10) {
2126  Peak->Phi = mDb->radiansPerBoundary() / 2
2127  + ((Peak->PadPosition-1) + 0.5) * mDb->radiansPerPad()
2128  + PhiDeflect + iSec * (mDb->numberOfPads() * mDb->radiansPerPad()
2129  + mDb->radiansPerBoundary())+halfpi;
2130  }
2131  // for FTPC East
2132  /* Invert pad number (== Peak->PadPosition) for FTPC East */
2133  /* (not yet understood where and why pad numbers were inverted) */
2134  if (iRow >= 10) {
2135  Peak->Phi = mDb->radiansPerBoundary() / 2
2136  + (159.5 - (Peak->PadPosition-1))* mDb->radiansPerPad()
2137  - PhiDeflect + iSec * (mDb->numberOfPads() * mDb->radiansPerPad()
2138  + mDb->radiansPerBoundary())+halfpi;
2139  }
2140 
2141  // ===================================================================
2142 
2143  // shift time => radius to correct for the offset of the inner cathode :
2144 
2145  if (fabs(mOffsetCathodeWest)>0 || fabs(mOffsetCathodeEast)>0)
2146  {
2147 
2148  if (iRow<10) // correct for west chamber
2149  TimeCoordinate=(0.999997-0.09739494018294076*mOffsetCathodeWest*cos(Peak->Phi-mAngleOffsetWest))*TimeCoordinate;
2150  //TimeCoordinate=(1.0-0.09739494018294076*mOffsetCathodeWest*cos(Peak->Phi))*TimeCoordinate;
2151  else // correct for east chamber
2152  TimeCoordinate=(0.999997-0.09739494018294076*mOffsetCathodeEast*sin(Peak->Phi+mAngleOffsetEast))*TimeCoordinate;
2153  //TimeCoordinate=(1.0-0.09739494018294076*mOffsetCathodeEast*sin(Peak->Phi))*TimeCoordinate;
2154 
2155  // calculate new corrected radius and deflection angle :
2156 
2157  //Peak->TimePosition=TimeCoordinate; // DEBUG only !
2158 
2159  PadtransLower= (int) (TimeCoordinate*PadtransPerTimebin);
2160 
2161  if ( TimeCoordinate > mDb->numberOfTimebins() || TimeCoordinate <= 0 )
2162  {
2163  // exceeds table dimensions
2164  return FALSE;
2165  }
2166 
2167  /* linear interpolation in radius table */
2168 
2169  Peak->Rad=pRadius[iRow + mDb->numberOfPadrowsPerSide() * PadtransLower]
2170  -(pRadius[iRow + mDb->numberOfPadrowsPerSide()*(PadtransLower)]
2171  -pRadius[iRow + mDb->numberOfPadrowsPerSide() * (PadtransLower+1)])/2;
2172 
2173  if ( pRadius[iRow + mDb->numberOfPadrowsPerSide() * PadtransLower] == 0 )
2174  {
2175  // outside FTPC sensitive region
2176  return FALSE;
2177  }
2178 
2179  // calculate new corrected deflection angle and phi
2180 
2181  /* linear interpolation in deflection table */
2182  PhiDeflect=pDeflection[iRow + mDb->numberOfPadrowsPerSide() * PadtransLower]
2183  +(pDeflection[iRow + mDb->numberOfPadrowsPerSide() * (PadtransLower+1)]
2184  -pDeflection[iRow + mDb->numberOfPadrowsPerSide() * PadtransLower])/2;
2185 
2186 
2187  /* calculate phi angle from pad position */
2188  // for FTPC West
2189  if (iRow < 10) {
2190  Peak->Phi = mDb->radiansPerBoundary() / 2
2191  + ((Peak->PadPosition-1) + 0.5) * mDb->radiansPerPad()
2192  + PhiDeflect + iSec * (mDb->numberOfPads() * mDb->radiansPerPad()
2193  + mDb->radiansPerBoundary())+halfpi;
2194  }
2195  // for FTPC East
2196  /* Invert pad number (== Peak->PadPosition) for FTPC East */
2197  /* (not yet understood where and why pad numbers were inverted) */
2198  if (iRow >= 10) {
2199  Peak->Phi = mDb->radiansPerBoundary() / 2
2200  + (159.5 - (Peak->PadPosition-1))* mDb->radiansPerPad()
2201  - PhiDeflect + iSec * (mDb->numberOfPads() * mDb->radiansPerPad()
2202  + mDb->radiansPerBoundary())+halfpi;
2203  }
2204  }
2205 
2206  // ===================================================================
2207 
2208  /* transform to cartesian */
2209  Peak->x = Peak->Rad*cos(Peak->Phi);
2210  /* Peak->x = -Peak->x to transform FTPC West into STAR global coordinate system */
2211  if (iRow <10) {
2212  Peak->x = -Peak->x;
2213  }
2214  Peak->y = Peak->Rad*sin(Peak->Phi);
2215  Peak->z = mDb->padrowZPosition(iRow);
2216  return TRUE;
2217 }
2218 
2219 float StFtpcClusterFinder::gauss_2d(int x, int y, float x1, float sigma_x1, float y1, float sigma_y1, float amp1)
2220 {
2221  float result;
2222  float g1, g2;
2223 
2224  if(sigma_x1==0 || sigma_y1==0)
2225  {
2226  return 0.0;
2227  }
2228 
2229  g1=exp(-sqr((float)x-x1)/(2*sqr(sigma_x1)));
2230  g2=exp(-sqr((float)y-y1)/(2*sqr(sigma_y1)));
2231 
2232  result=amp1*g1*g2;
2233 
2234  return result;
2235 }
2236 
2237 float StFtpcClusterFinder::sigmax(float timebin)
2238 {
2239  return 0.9;
2240 }
2241 
2242 float StFtpcClusterFinder::sigmat(float timebin)
2243 {
2244  return 0.6;
2245 }
2246 
2247 int StFtpcClusterFinder::calcpadtrans(double *pradius,
2248  double *pdeflection, double deltap)
2249 {
2250  int i, j, v_buf, padrow;
2251  double t_last, t_next, r_last, r_next, e_now, v_now, psi_now;
2252  double step_size;
2253 
2254  step_size=((float) mDb->numberOfTimebins()
2255  / (float) mParam->numberOfDriftSteps());
2256 
2257  LOG_DEBUG <<"StFtpcClusterFinder::calcpadtrans deltap = "<<deltap<<endm;
2258 
2259  for (padrow=0; padrow<mDb->numberOfPadrowsPerSide(); padrow++)
2260  {
2261  /* determine starting values */
2262  t_last=0;
2263  v_buf=0;
2264  r_last=mDb->sensitiveVolumeOuterRadius();
2265  pradius[padrow]=mDb->sensitiveVolumeOuterRadius();
2266  pdeflection[padrow]=0;
2267  e_now = mDb->radiusTimesField() / (0.5*r_last);
2268  for(j=v_buf; j<mDb->numberOfMagboltzBins()-1
2269  && mDb->magboltzEField(j) < e_now; j++);
2270  if(j<1 || j>mDb->numberOfMagboltzBins())
2271  {
2272  LOG_ERROR << "calcpadtrans error 1: j=" << j << ", v_buf=" << v_buf << " e_drift=" << mDb->magboltzEField(j) << ", e_now=" << e_now << endm;
2273  return FALSE;
2274  }
2275  v_buf=j-1;
2276  v_now=((mDb->magboltzVDrift(v_buf, padrow)
2277  +deltap*mDb->magboltzdVDriftdP(v_buf, padrow))
2278  *(mDb->magboltzEField(j)-e_now)
2279  +(mDb->magboltzVDrift(j, padrow)
2280  +deltap*mDb->magboltzdVDriftdP(j, padrow))
2281  *(e_now-mDb->magboltzEField(v_buf)))
2282  /(mDb->magboltzEField(j)-mDb->magboltzEField(v_buf));
2283  psi_now=((mDb->magboltzDeflection(v_buf,padrow)
2284  +deltap*mDb->magboltzdDeflectiondP(v_buf,padrow))
2285  *mParam->lorentzAngleFactor()
2286  *(mDb->magboltzEField(j)-e_now)
2287  +(mDb->magboltzDeflection(j,padrow)
2288  +deltap*mDb->magboltzdDeflectiondP(j,padrow))
2289  *mParam->lorentzAngleFactor()
2290  *(e_now-mDb->magboltzEField(v_buf)))
2291  /(mDb->magboltzEField(j)-mDb->magboltzEField(v_buf));
2292  for (i=0; i<mParam->numberOfDriftSteps()-1
2293  && e_now < mDb->magboltzEField(mDb->numberOfMagboltzBins()-2)
2294  ; i++)
2295  {
2296  t_next = t_last + step_size;
2297  /* first guess for r_next: */
2298  r_next = r_last - v_now * step_size * mDb->microsecondsPerTimebin();
2299  e_now = mDb->radiusTimesField() / (0.5*(r_last+r_next));
2300 
2301  for(j=v_buf; j<mDb->numberOfMagboltzBins()-1
2302  && mDb->magboltzEField(j) < e_now; j++);
2303 
2304  if(j<1 || j>mDb->numberOfMagboltzBins())
2305  {
2306  LOG_ERROR << "calcpadtrans error 2: j=" << j << ", v_buf=" << v_buf << " e_drift=" << mDb->magboltzEField(j) << ", e_now=" << e_now << endm;
2307  return FALSE;
2308  }
2309 
2310  v_buf=j-1;
2311  v_now=((mDb->magboltzVDrift(v_buf, padrow)
2312  +deltap*mDb->magboltzdVDriftdP(v_buf, padrow))
2313  *(mDb->magboltzEField(j)-e_now)
2314  +(mDb->magboltzVDrift(j, padrow)
2315  +deltap*mDb->magboltzdVDriftdP(j, padrow))
2316  *(e_now-mDb->magboltzEField(v_buf)))
2317  /(mDb->magboltzEField(j)-mDb->magboltzEField(v_buf));
2318  psi_now=((mDb->magboltzDeflection(v_buf,padrow)
2319  +deltap*mDb->magboltzdDeflectiondP(v_buf,padrow))
2320  *mParam->lorentzAngleFactor()
2321  *(mDb->magboltzEField(j)-e_now)
2322  +(mDb->magboltzDeflection(j,padrow)
2323  +deltap*mDb->magboltzdDeflectiondP(j,padrow))
2324  *mParam->lorentzAngleFactor()
2325  *(e_now-mDb->magboltzEField(v_buf)))
2326  /(mDb->magboltzEField(j)-mDb->magboltzEField(v_buf));
2327 
2328  /* correct r_next: */
2329  r_next = r_last - v_now * step_size *mDb->microsecondsPerTimebin();
2330  pradius[padrow+mDb->numberOfPadrowsPerSide()*(i+1)]=r_next;
2331  pdeflection[padrow+mDb->numberOfPadrowsPerSide()*(i+1)]
2332  =pdeflection[padrow+mDb->numberOfPadrowsPerSide()*i]
2333  +((r_last-r_next)*tan(degree * psi_now)/r_last);
2334  t_last=t_next;
2335  r_last=r_next;
2336  }
2337  if (DebugOn) {
2338  LOG_DEBUG<<i<<" steps calculated, padrow "<<padrow<<endm;
2339  }
2340  }
2341  return TRUE;
2342 }
2343 
2344 int StFtpcClusterFinder::cucInit(TClusterUC memory[MAXNUMCUC], int RealMemory[MAXNUMCUC], int *pointer)
2345 {
2346  int i;
2347 
2348  for(i=0; i<MAXNUMCUC; i++)
2349  {
2350  RealMemory[i]=i;
2351  }
2352  *pointer = -1;
2353  return TRUE;
2354 }
2355 
2356 TClusterUC *StFtpcClusterFinder::cucAlloc(TClusterUC memory[MAXNUMCUC], int RealMemory[MAXNUMCUC], int *pointer)
2357 {
2358 
2359  if((*pointer)+1 < MAXNUMCUC)
2360  {
2361  (*pointer)++;
2362  memory[RealMemory[*pointer]].MemoryPtr = RealMemory[*pointer];
2363  return &memory[RealMemory[*pointer]];
2364  }
2365  else
2366  {
2367  if (DebugOn) {
2368  LOG_DEBUG<<"CUC memory exhausted! requested "<<*pointer<<" CUCs"<<endm;
2369  }
2370  return 0;
2371  }
2372 }
2373 
2374 int StFtpcClusterFinder::cucFree(TClusterUC memory[MAXNUMCUC], int RealMemory[MAXNUMCUC], int *pointer, TClusterUC *OldCUC)
2375 {
2376  if(*pointer >= 0)
2377  {
2378  RealMemory[*pointer] = OldCUC->MemoryPtr;
2379  (*pointer)--;
2380  return TRUE;
2381  }
2382  else
2383  {
2384  if (DebugOn) {
2385  LOG_DEBUG<<"CUC memory management confused!"<<endm;
2386  }
2387  return FALSE;
2388  }
2389 }
kStWarn
Definition: Stypes.h:42
kStOK
Definition: Stypes.h:40
kStERR
Definition: Stypes.h:45
Generated by   doxygen 1.8.5