StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StIstScanClusterAlgo.cxx
1 #include "StEvent.h"
2 #include "StEvent/StEnumerations.h"
3 #include "StMessMgr.h"
4 #include "StIstScanClusterAlgo.h"
5 #include "StIstUtil/StIstCollection.h"
6 #include "StIstUtil/StIstRawHitCollection.h"
7 #include "StIstUtil/StIstRawHit.h"
8 #include "StIstUtil/StIstClusterCollection.h"
9 #include "StIstUtil/StIstCluster.h"
10 #include "StIstUtil/StIstConsts.h"
11 
12 #include <math.h>
13 #include <iostream>
14 #include <iterator>
15 #include <vector>
16 #include <algorithm>
17 
18 
19 Int_t StIstScanClusterAlgo::doClustering(const StIstCollection &istCollection, StIstRawHitCollection &rawHitsOriginal, StIstClusterCollection &clusters )
20 {
21  StIstCluster *newCluster = 0;
22  StIstRawHit *rawHitTemp = 0;
23  StIstRawHit *rawHitMaxAdcTemp = 0;
24  Int_t clusterType = kIstScanClusterAlgo;
25  Int_t maxTb = -1, usedTb = -1;
26  Int_t ladder = 0, sensor = 0;
27  Float_t meanRow = 0., meanColumn = 0.;
28  Float_t totCharge = 0., totChargeErr = 0.;
29  Int_t clusterSize = 0, clusterSizeRPhi = 0, clusterSizeZ = 0;
30  unsigned short idTruth = 0;
31 
32  //get number of time bin used in this event
33  Int_t nTimeBins = istCollection.getNumTimeBins();
34 
35  //sort raw hits in increasing order by geometry ID
36  rawHitsOriginal.sortByGeoId();
37 
38  //copy raw hit collection to temporary vectors
39  std::vector<StIstRawHit *> rawHitsVec[kIstNumSensorsPerLadder][kIstNumColumnsPerSensor];
40  std::vector<StIstCluster *> clustersVec[kIstNumSensorsPerLadder][kIstNumColumnsPerSensor];
41  std::vector<StIstRawHit *> rawHitsToMerge;
42 
43  for (int sensorIdx = 0; sensorIdx < kIstNumSensorsPerLadder; sensorIdx++) {
44  for (int columnIdx = 0; columnIdx < kIstNumColumnsPerSensor; columnIdx++) {
45  rawHitsVec[sensorIdx][columnIdx].reserve(kIstNumRowsPerSensor);
46  clustersVec[sensorIdx][columnIdx].reserve(kIstNumRowsPerSensor);
47  }
48  }
49 
50  rawHitsToMerge.reserve(kIstNumRowsPerSensor);
51 
52  for (std::vector< StIstRawHit * >::iterator rawHitPtr = rawHitsOriginal.getRawHitVec().begin(); rawHitPtr != rawHitsOriginal.getRawHitVec().end(); ++rawHitPtr) {
53  int sensorIndex = (int)(*rawHitPtr)->getSensor();
54  int columnIndex = (int)(*rawHitPtr)->getColumn();
55  rawHitsVec[sensorIndex - 1][columnIndex - 1].push_back( new StIstRawHit( *(*rawHitPtr)) );
56  }
57 
58  //do clustering
59  int clusterLabel = 0;
60 
61  for (int sensorIdx = 0; sensorIdx < kIstNumSensorsPerLadder; sensorIdx++)
62  {
63  //step 1: do clustering for each column
64  for (int columnIdx = 0; columnIdx < kIstNumColumnsPerSensor; columnIdx++)
65  {
66  while ( !rawHitsVec[sensorIdx][columnIdx].empty() )
67  {
68  rawHitTemp = rawHitsVec[sensorIdx][columnIdx].back();
69  rawHitsVec[sensorIdx][columnIdx].pop_back();
70  rawHitsToMerge.push_back(rawHitTemp);
71 
72  //count number to merge
73  int nToMerge = 1;
74  //find all raw hits that are neighboring
75  std::vector<StIstRawHit *>::iterator rawHitsToMergePtr = rawHitsToMerge.begin();
76  rawHitMaxAdcTemp = *rawHitsToMergePtr;
77  StIstRawHit *rawHitNext = 0;
78  StIstRawHit *rawHitTempBack = 0;
79 
80  while (rawHitsToMergePtr != rawHitsToMerge.end() && !rawHitsVec[sensorIdx][columnIdx].empty()) {
81  rawHitTemp = rawHitsVec[sensorIdx][columnIdx].back();
82  rawHitTempBack = rawHitsVec[sensorIdx][columnIdx].back();
83 
84  if ( (*rawHitsToMergePtr)->getRow() == rawHitTemp->getRow() + 1 ) {
85  rawHitsVec[sensorIdx][columnIdx].pop_back();
86 
87  if (!rawHitsVec[sensorIdx][columnIdx].empty()) {
88  rawHitNext = rawHitsVec[sensorIdx][columnIdx].back();
89 
90  if ( (rawHitTemp->getRow() == rawHitNext->getRow() + 1) &&
91  (rawHitTemp->getCharge(rawHitTemp->getMaxTimeBin()) < (*rawHitsToMergePtr)->getCharge((*rawHitsToMergePtr)->getMaxTimeBin())) &&
92  (rawHitTemp->getCharge(rawHitTemp->getMaxTimeBin()) < rawHitNext->getCharge(rawHitNext->getMaxTimeBin())) )
93  {
94  float weightBack = rawHitNext->getCharge(rawHitNext->getMaxTimeBin()) / ((*rawHitsToMergePtr)->getCharge((*rawHitsToMergePtr)->getMaxTimeBin()) + rawHitNext->getCharge(rawHitNext->getMaxTimeBin()));
95 
96  for (int iTB = 0; iTB < nTimeBins; iTB++) {
97  rawHitTempBack->setCharge(weightBack * rawHitTemp->getCharge(iTB), iTB);
98  rawHitTemp->setCharge((1.0 - weightBack) * rawHitTemp->getCharge(iTB), iTB);
99  }
100 
101  rawHitsVec[sensorIdx][columnIdx].push_back(rawHitTempBack);
102  }
103  }
104 
105  ++nToMerge;
106  rawHitsToMerge.push_back(rawHitTemp);
107 
108  if ( rawHitTemp->getCharge(rawHitTemp->getMaxTimeBin()) > (*rawHitsToMergePtr)->getCharge((*rawHitsToMergePtr)->getMaxTimeBin()) )
109  rawHitMaxAdcTemp = rawHitTemp;
110  }
111 
112  ++rawHitsToMergePtr;
113  }
114 
115  //used time bin index (raw hits with maximum ADC holds the time-bin priority)
116  maxTb = rawHitMaxAdcTemp->getMaxTimeBin();
117  idTruth = rawHitMaxAdcTemp->getIdTruth();
118 
119  if (maxTb < 0 || maxTb >= nTimeBins) maxTb = rawHitMaxAdcTemp->getDefaultTimeBin();
120 
121  if (mTimeBin < nTimeBins) usedTb = mTimeBin;
122  else usedTb = maxTb;
123 
124  ladder = rawHitMaxAdcTemp->getLadder();
125  sensor = rawHitMaxAdcTemp->getSensor(); // = sensorIdx + 1
126  meanColumn = (float)rawHitMaxAdcTemp->getColumn(); // = columnIdx + 1
127  clusterSize = nToMerge;
128  clusterSizeRPhi = nToMerge;
129  clusterSizeZ = 1;
130 
131  float tempCharge[nToMerge], tempChargeErr[nToMerge], tempRow[nToMerge];
132 
133  for (int i = 0; i < nToMerge; i++) {
134  tempCharge[i] = 0.; tempChargeErr[i] = 0.; tempRow[i] = 0.;
135  }
136 
137  float tempSumCharge = 0, tempSumChargeErrSquare = 0.;
138  int mergeIdx = 0;
139 
140  for (rawHitsToMergePtr = rawHitsToMerge.begin(); rawHitsToMergePtr != rawHitsToMerge.end() && mergeIdx < nToMerge; ++rawHitsToMergePtr, ++mergeIdx) {
141  tempCharge[mergeIdx] = (*rawHitsToMergePtr)->getCharge(usedTb);
142  tempChargeErr[mergeIdx] = (*rawHitsToMergePtr)->getChargeErr(usedTb);
143  tempRow[mergeIdx] = (float)(*rawHitsToMergePtr)->getRow();
144  tempSumCharge += (*rawHitsToMergePtr)->getCharge(usedTb);
145  }
146 
147  meanRow = 0;
148 
149  for (int iRawHit = 0; iRawHit < nToMerge; iRawHit++) {
150  meanRow += tempRow[iRawHit] * tempCharge[iRawHit] / tempSumCharge;
151  tempSumChargeErrSquare += tempChargeErr[iRawHit] * tempChargeErr[iRawHit];
152  }
153 
154  totCharge = tempSumCharge;
155  totChargeErr = sqrt(tempSumChargeErrSquare / nToMerge);
156 
157  newCluster = new StIstCluster((int)ladder * 10000 + clusterLabel, ladder, sensor, meanRow, meanColumn, totCharge, totChargeErr, clusterType);
158  newCluster->setNRawHits(clusterSize);
159  newCluster->setNRawHitsRPhi(clusterSizeRPhi);
160  newCluster->setNRawHitsZ(clusterSizeZ);
161  newCluster->setMaxTimeBin(maxTb);
162  newCluster->setIdTruth(idTruth);
163 
164  clustersVec[sensorIdx][columnIdx].push_back(newCluster);
165  clusterLabel++;
166 
167  rawHitsToMerge.clear();
168  }//end current column raw hits loop
169  }//end current sensor raw hits loop
170 
171  //step 2: do clustering for neighboring rows
172  std::vector<StIstCluster *>::iterator clusterIt1, clusterIt2;
173 
174  for (int columnIdx1 = 0; columnIdx1 < kIstNumColumnsPerSensor - 1; columnIdx1++) {
175  int columnIdx2 = columnIdx1 + 1;
176 
177  if (clustersVec[sensorIdx][columnIdx1].size() > 0 && clustersVec[sensorIdx][columnIdx2].size() > 0) {
178  for (clusterIt1 = clustersVec[sensorIdx][columnIdx1].begin(); clusterIt1 != clustersVec[sensorIdx][columnIdx1].end() && !clustersVec[sensorIdx][columnIdx1].empty(); clusterIt1++) {
179  for (clusterIt2 = clustersVec[sensorIdx][columnIdx2].begin(); clusterIt2 != clustersVec[sensorIdx][columnIdx2].end() && !clustersVec[sensorIdx][columnIdx1].empty(); clusterIt2++) {
180  float rowDistance = (*clusterIt1)->getMeanRow() - (*clusterIt2)->getMeanRow();
181 
182  if (TMath::Abs(rowDistance) < 0.5) { //here 0.5 means the distance between two clusters' weighted centers in row direction smaller than 0.5
183  maxTb = (*clusterIt1)->getMaxTimeBin();
184  idTruth = (*clusterIt1)->getIdTruth();
185  if((*clusterIt1)->getTotCharge() < (*clusterIt2)->getTotCharge()) {
186  maxTb = (*clusterIt2)->getMaxTimeBin();
187  idTruth = (*clusterIt2)->getIdTruth();
188  }
189 
190  totCharge = (*clusterIt1)->getTotCharge() + (*clusterIt2)->getTotCharge();
191  totChargeErr = sqrt(((*clusterIt1)->getTotChargeErr() * (*clusterIt1)->getTotChargeErr() * (*clusterIt1)->getNRawHits() + (*clusterIt2)->getTotChargeErr() * (*clusterIt2)->getTotChargeErr() * (*clusterIt2)->getNRawHits()) / ((*clusterIt1)->getNRawHits() + (*clusterIt2)->getNRawHits()));
192  clusterSize = (*clusterIt1)->getNRawHits() + (*clusterIt2)->getNRawHits();
193  clusterSizeRPhi = (*clusterIt1)->getNRawHitsRPhi() + (*clusterIt2)->getNRawHitsRPhi() - 1;
194  clusterSizeZ = (*clusterIt1)->getNRawHitsZ() + (*clusterIt2)->getNRawHitsZ();
195  meanRow = (*clusterIt1)->getMeanRow() * (*clusterIt1)->getTotCharge() / totCharge + (*clusterIt2)->getMeanRow() * (*clusterIt2)->getTotCharge() / totCharge;
196  meanColumn = (*clusterIt1)->getMeanColumn() * (*clusterIt1)->getTotCharge() / totCharge + (*clusterIt2)->getMeanColumn() * (*clusterIt2)->getTotCharge() / totCharge;
197 
198  (*clusterIt2)->setMeanRow(meanRow);
199  (*clusterIt2)->setMeanColumn(meanColumn);
200  (*clusterIt2)->setTotCharge(totCharge);
201  (*clusterIt2)->setTotChargeErr(totChargeErr);
202  (*clusterIt2)->setNRawHits(clusterSize);
203  (*clusterIt2)->setNRawHitsRPhi(clusterSizeRPhi);
204  (*clusterIt2)->setNRawHitsZ(clusterSizeZ);
205  (*clusterIt2)->setIdTruth(idTruth);
206 
207  int distance1 = std::distance(clustersVec[sensorIdx][columnIdx1].begin(), clusterIt1);
208  clustersVec[sensorIdx][columnIdx1].erase(clusterIt1);
209 
210  if (distance1 == 0)
211  clusterIt1 = clustersVec[sensorIdx][columnIdx1].begin();
212  else
213  --clusterIt1;
214  }//end merge
215  }//end clusters vector 2 loop
216  }//end clusters vector 1 loop
217  }//end column cluster number cut
218  }//end current sensor cluster loop
219  }//end all sensor clustering loop
220 
221  //fill output container
222  std::vector<StIstCluster *>::iterator clusterIt;
223 
224  for (int sensorIdx = 0; sensorIdx < kIstNumSensorsPerLadder; sensorIdx++) {
225  for (int columnIdx = 0; columnIdx < kIstNumColumnsPerSensor; columnIdx++) {
226  if (clustersVec[sensorIdx][columnIdx].size() <= 0) continue;
227 
228  for (clusterIt = clustersVec[sensorIdx][columnIdx].begin(); clusterIt != clustersVec[sensorIdx][columnIdx].end(); ++clusterIt)
229  clusters.getClusterVec().push_back(*clusterIt);
230 
231  rawHitsVec[sensorIdx][columnIdx].clear();
232  clustersVec[sensorIdx][columnIdx].clear();
233  }
234  }
235 
236  return kStOk;
237 }
238 
239 
240 /***************************************************************************
241 *
242 * $Log: StIstScanClusterAlgo.cxx,v $
243 * Revision 1.20 2018/01/04 17:34:37 smirnovd
244 * [Cosmetic] Remove StRoot/ from include path
245 *
246 * $STAR/StRoot is already in the default path search
247 *
248 * Revision 1.19 2016/02/18 17:27:19 smirnovd
249 * StIstScanClusterAlgo: Whitespace corrected in "Add idTruth for clusters."
250 *
251 * Corrected commit 45e825e1
252 *
253 * Revision 1.18 2016/02/17 19:52:50 huangbc
254 * Add idTruth for clusters.
255 *
256 * Revision 1.17 2015/05/20 20:53:57 smirnovd
257 * Removed a priori true condition without changing the logic
258 *
259 * mTimeBin is unsigned char always >= 0
260 *
261 * Revision 1.16 2015/03/04 16:17:16 smirnovd
262 * Added a check for empty container of to-be-merged proto-clusters
263 *
264 * This is done to avoid a situation when we run out of column-wise clusters in the
265 * current (happens when all clusters are merged and removed from the temporary
266 * container) column but the next one still has some cluster left.
267 *
268 * Revision 1.15 2015/03/04 16:17:07 smirnovd
269 * Revert "empty check was added for the to-be-merged column-wise proto-clusters container, to fix the RT #3056 oberserved by Lidia"
270 *
271 * Revision 1.13 2014/09/18 06:27:25 ypwang
272 * remove unneccessary check for raw hit electroincis ID check
273 *
274 * Revision 1.12 2014/09/17 20:39:45 smirnovd
275 * Squashed commit of the following:
276 *
277 * commit 37d3d404a31c9b152811232af55d37177162269d
278 * Author: Dmitri Smirnov <d.s@plexoos.com>
279 * Date: Wed Sep 17 16:11:22 2014 -0400
280 *
281 * Added an author to reflect on contributions
282 *
283 * commit 6ceacb443d2d35bc21295b81a3d25b7433d40260
284 * Author: Dmitri Smirnov <d.s@plexoos.com>
285 * Date: Wed Sep 17 16:09:48 2014 -0400
286 *
287 * [Minor] Reversed the logic and saved one level of intentation
288 *
289 * commit 4bc24031445ecce9f19b940697d13cc8a755aaf1
290 * Author: Dmitri Smirnov <d.s@plexoos.com>
291 * Date: Wed Sep 17 16:06:42 2014 -0400
292 *
293 * Do not use standard ROOT's dictionary macroses since the classes are transient by design
294 *
295 * Revision 1.11 2014/09/17 20:33:32 smirnovd
296 * Squashed commit of the following:
297 *
298 * commit 72dc19a6663ea31c719c1a61f6d2b4752dd766aa
299 * Author: Dmitri Smirnov <d.s@plexoos.com>
300 * Date: Wed Sep 17 12:34:42 2014 -0400
301 *
302 * Minor code refactoring, clean up
303 *
304 * commit e083a10a9fb60b7dcce692ef8043b9227c12768b
305 * Author: Dmitri Smirnov <d.s@plexoos.com>
306 * Date: Wed Sep 17 12:18:16 2014 -0400
307 *
308 * Removed pointless comments
309 *
310 * commit 88d51857362c91c954704cec4a31a0b0fa7fccc5
311 * Author: Dmitri Smirnov <d.s@plexoos.com>
312 * Date: Wed Sep 17 12:17:26 2014 -0400
313 *
314 * Updated description in doxygen comments
315 *
316 * commit eb09527489179fc7dab6aa7f23fd132b25185bb1
317 * Author: Dmitri Smirnov <d.s@plexoos.com>
318 * Date: Tue Sep 9 15:15:56 2014 -0400
319 *
320 * StIstScanClusterAlgo: Removed unused variable
321 *
322 * commit 1a8df63533c71a0e2ba4d8275ebf89f4e3004765
323 * Author: Dmitri Smirnov <d.s@plexoos.com>
324 * Date: Fri Aug 22 16:04:47 2014 -0400
325 *
326 * Neatened headers: Removed unused, spelled paths in includes explicitly as it slightly helps in identifying dependencies
327 *
328 * commit 972e8ed41403bd680ade5ecc509f8bca004e86ee
329 * Author: Dmitri Smirnov <d.s@plexoos.com>
330 * Date: Wed Sep 17 12:34:20 2014 -0400
331 *
332 * Minor stylistic changes
333 *
334 * commit 57daf5a1e0b3246fd12f1dd1c2ca089b62930c83
335 * Author: Dmitri Smirnov <d.s@plexoos.com>
336 * Date: Tue Sep 16 16:29:14 2014 -0400
337 *
338 * Improved doxygen comments
339 *
340 * Revision 1.10 2014/09/09 07:34:07 ypwang
341 * data type was updated from UChar_t to Int_t for several varibales, and the nTimeBins was corrected as a non-static varible
342 *
343 * Revision 1.9 2014/09/07 13:54:45 ypwang
344 * move setUsedTimeBin() and setSplitFlag() setters from inherited classes to their base class StIstIClusterAlgo.h
345 *
346 * Revision 1.8 2014/09/07 11:41:36 ypwang
347 * ClassDef version updated from 1 to 0, and remove Init() function
348 *
349 * Revision 1.7 2014/08/22 15:55:15 smirnovd
350 * Fixed style with astyle -s3 -p -H -A3 -k3 -O -o -y -Y -f
351 *
352 * Revision 1.6 2014/08/22 15:50:00 smirnovd
353 * Moved CVS history to the end of file
354 *
355 * Revision 1.5 2014/03/17 21:51:56 ypwang
356 * minor update due to some IST constants moved to StEnumurations.h
357 *
358 * Revision 1.4 2014/02/16 23:18:34 ypwang
359 * getting number of time bins used in current event by StIstCollection::getNumTimeBins() function
360 *
361 * Revision 1.3 2014/02/08 03:34:16 ypwang
362 * updating scripts
363 *
364 *
365 ****************************************************************************
366 * StIstScanClusterAlgo.cxx,v 1.0
367 * Revision 1.0 2013/11/04 15:55:30 Yaping
368 * Initial version
369 ****************************************************************************/
StIstRawHit::getLadder
unsigned char getLadder() const
1-24
Definition: StIstRawHit.cxx:46
StIstRawHit::getColumn
unsigned char getColumn() const
1-12
Definition: StIstRawHit.cxx:62
StIstRawHit::getRow
unsigned char getRow() const
1-64
Definition: StIstRawHit.cxx:56
StIstConsts::kIstNumColumnsPerSensor
const int kIstNumColumnsPerSensor
12 columns in beam direction per each sensor
Definition: StEnumerations.h:716
StIstRawHit::getIdTruth
unsigned short getIdTruth() const
for embedding, 0 as background
Definition: StIstRawHit.cxx:44
StIstConsts::kIstNumSensorsPerLadder
const int kIstNumSensorsPerLadder
6 sensor per one IST Ladder
Definition: StEnumerations.h:712
StIstRawHit::getSensor
unsigned char getSensor() const
1-6
Definition: StIstRawHit.cxx:51
StIstConsts::kIstNumRowsPerSensor
const int kIstNumRowsPerSensor
64 rows in r-phi direction per each sensor
Definition: StEnumerations.h:715
kStOk
Definition: Stypes.h:41
Generated by   doxygen 1.8.5