StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StFcsClusterMaker.cxx
1 // $Id: StFcsClusterMaker.cxx,v 1.1 2021/03/30 13:40:02 akio Exp $
2 //
3 // $Log: StFcsClusterMaker.cxx,v $
4 // Revision 1.1 2021/03/30 13:40:02 akio
5 // FCS code after peer review and moved from $CVSROOT/offline/upgrades/akio
6 //
7 // Revision 1.18 2021/02/25 21:52:57 akio
8 // Int_t -> int
9 //
10 // Revision 1.17 2021/02/25 19:24:15 akio
11 // Modified for STAR code review (Hongwei)
12 //
13 // Revision 1.16 2020/12/17 21:18:01 akio
14 // Separate RATIO2SPLIT for ecal/hcal
15 //
16 // Revision 1.15 2020/09/03 19:42:24 akio
17 // moving sum & fit to StFcsWaveformFitMaker
18 //
19 // Revision 1.14 2019/11/22 17:26:17 akio
20 // fix typo
21 //
22 // Revision 1.13 2019/11/22 15:51:16 akio
23 // adding categorization
24 //
25 // Revision 1.12 2019/08/16 16:26:28 akio
26 // adding gaincorr to sum
27 //
28 // Revision 1.11 2019/08/01 18:37:01 akio
29 // calling makesum/makefit for all detectors
30 //
31 // Revision 1.10 2019/07/15 16:58:32 akio
32 // Adding hit->cluster pointer
33 //
34 // Revision 1.9 2019/07/10 06:14:17 akio
35 // adding cluster pointer to hit
36 //
37 // Revision 1.8 2019/07/03 16:13:29 akio
38 // separate neighbor distance for ecal and hcal
39 //
40 // Revision 1.7 2019/06/27 16:12:42 akio
41 // Using getLocalXYinCell for cell/cluster distance measure and in moment analysis
42 //
43 // Revision 1.6 2019/06/26 18:00:34 akio
44 // added some adjustable parameters and way to select how to get energy
45 //
46 // Revision 1.5 2019/06/25 16:46:56 akio
47 // removed debugging comment
48 //
49 // Revision 1.4 2019/06/25 16:38:17 akio
50 // Added TOWER_E_RATIO2SPLIT, neighbor clusters
51 //
52 // Revision 1.3 2019/06/21 16:32:42 akio
53 // Added neighbor cluster and factor threshold for cluster splitting at valley
54 //
55 // Revision 1.2 2019/06/07 18:17:24 akio
56 // added summing adc
57 //
58 // Revision 1.1 2018/11/14 16:50:11 akio
59 // FCS codes in offline/upgrade/akio
60 //
61 
62 #include "StFcsClusterMaker.h"
63 
64 #include "StLorentzVectorF.hh"
65 
66 #include "StMessMgr.h"
67 #include "StEvent/StEventTypes.h"
68 #include "StEvent/StFcsHit.h"
69 #include "StEvent/StFcsCluster.h"
70 #include "StFcsDbMaker/StFcsDb.h"
71 #include "tables/St_g2t_track_Table.h"
72 
73 #include "StMuDSTMaker/COMMON/StMuTypes.hh"
74 #include "StMuDSTMaker/COMMON/StMuDst.h"
75 #include "StMuDSTMaker/COMMON/StMuEvent.h"
76 
77 #include <cmath>
78 #include "TMath.h"
79 #include "TVector2.h"
80 
81 ClassImp(StFcsClusterMaker)
82 
83 StFcsClusterMaker::StFcsClusterMaker(const char* name) : StMaker(name) {}
84 
85 StFcsClusterMaker::~StFcsClusterMaker() { }
86 
87 void StFcsClusterMaker::Clear(Option_t* option) {
88  StMaker::Clear(option);
89 }
90 
91 int StFcsClusterMaker::InitRun(int runNumber) {
92  // Ensure we can access database information
93  LOG_DEBUG << "StFcsClusterMaker initializing run" << endm;
94  //mDb = static_cast<StFcsDbMaker*>(GetMaker("fcsDb"));
95  mDb = static_cast<StFcsDb*>(GetDataSet("fcsDb"));
96  if (!mDb) {
97  LOG_ERROR << "StFcsClusterMaker initializing failed due to no StFcsDb" << endm;
98  return kStErr;
99  }
100  return StMaker::InitRun(runNumber);
101 }
102 
104  LOG_DEBUG << "StFcsClusterMaker Make!!!" << endm;
105 
106  StEvent* event = static_cast<StEvent*>(GetInputDS("StEvent"));
107  mFcsCollection=0;
108  if (event) mFcsCollection = event->fcsCollection();
109  if(!mFcsCollection) {
110  LOG_WARN << "StFcsClusterMaker did not find fcsCollection in StEvent" << endm;
111  return kStWarn;
112  }
113 
114  for(int det=0; det<=kFcsHcalSouthDetId; det++) {
115  if(det==0){
116  mNeighborDistance = mNeighborDistance_Ecal;
117  mDistanceAdvantage = mDistanceAdvantage_Ecal;
118  mTowerEThreSeed = mTowerEThreSeed_Ecal;
119  mTowerEThreshold = mTowerEThreshold_Ecal;
120  mTowerEThreMoment = mTowerEThreMoment_Ecal;
121  mTowerERatio2Split = mTowerERatio2Split_Ecal;
122  }
123  if(det==2){
124  mNeighborDistance = mNeighborDistance_Hcal;
125  mDistanceAdvantage = mDistanceAdvantage_Hcal;
126  mTowerEThreSeed = mTowerEThreSeed_Hcal;
127  mTowerEThreshold = mTowerEThreshold_Hcal;
128  mTowerEThreMoment = mTowerEThreMoment_Hcal;
129  mTowerERatio2Split = mTowerERatio2Split_Hcal;
130  }
131  makeCluster(det);
132  }
133  if(GetDebug()>0) mFcsCollection->print(3);
134  return kStOk;
135 }
136 
137 int StFcsClusterMaker::makeCluster(int det) {
138  StSPtrVecFcsHit& hits = mFcsCollection->hits(det);
139  StSPtrVecFcsCluster& clusters = mFcsCollection->clusters(det);
140  clusters.clear();
141  int nhit=hits.size();
142  for(int i=0; i<nhit; i++) hits[i]->setCluster(0); //reset cluster pointer from hit
143 
144  if(mSortById==0){ //sort by energy
145  std::sort(hits.begin(), hits.end(), [](StFcsHit* a, StFcsHit* b) {
146  return b->energy() < a->energy();
147  });
148  }else{ //sort by Id
149  std::sort(hits.begin(), hits.end(), [](StFcsHit* a, StFcsHit* b) {
150  return b->id() > a->id();
151  });
152  }
153 
154  for(int i=0; i<nhit; i++){ //loop over all hits
155  StFcsHit* hit=hits[i];
156  float e=hit->energy();
157  if(e < mTowerEThreshold && mSortById==0) break;
158  float neighborClusterId=-1;
159  float minDistance=999.0;
160  int ncluster = clusters.size();
161  int nNeighbor= 0;
162  StPtrVecFcsCluster neighbor;
163  for(int j=0; j<ncluster; j++){ //check all existing cluster
164  StFcsCluster* clu=clusters[j];
165  float neighborTowerE = isNeighbor(hit,clu);
166  if(neighborTowerE>0.0) { //found neighbor cluster
167  neighbor.push_back(clu);
168  nNeighbor++;
169  if(neighborTowerE * mTowerERatio2Split > e){ //merge to existing cluster
170  float d = distance(hit,clu);
171  if(d * mDistanceAdvantage < minDistance){
172  neighborClusterId=j;
173  minDistance=d;
174  }
175  }
176  }
177  }
178  StFcsCluster* cluster=0;
179  if(neighborClusterId==-1) {
180  //no neighbor, thus found new cluster seed
181  if(e >= mTowerEThreSeed){
182  cluster = new StFcsCluster();
183  cluster->setId(ncluster);
184  cluster->setDetectorId(det);
185  cluster->hits().push_back(hit);
186  hit->setCluster(cluster);
187  updateCluster(cluster);
188  mFcsCollection->addCluster(det,cluster);
189  neighbor.push_back(cluster);
190  nNeighbor++;
191  }else{
192  //no neighbor and not exceeding seed threshold
193  //what should we do??? I guess nothing...
194  }
195  }else{
196  //found neighbor tower which has higher energy
197  //add to the cluster with closest cluster
198  cluster = clusters[neighborClusterId];
199  cluster->hits().push_back(hit);
200  hit->setCluster(cluster);
201  updateCluster(cluster);
202  }
203  if(nNeighbor>1) { //more than 1 neighbors found
204  for(int k=0; k<nNeighbor-1; k++){
205  StFcsCluster* clu1 = neighbor[k];
206  for(int l=k+1; l<nNeighbor; l++){
207  StFcsCluster* clu2 = neighbor[l];
208  clu1->addNeighbor(clu2);
209  clu2->addNeighbor(clu1);
210  }
211  }
212  }
213  }
214 
215  //loop over all found cluster and fill fourMomentum and do moment analysis
216  int nc = clusters.size();
217  for(int j=0; j<nc; j++){
218  StFcsCluster* clu=clusters[j];
219  StThreeVectorD xyz = mDb->getStarXYZfromColumnRow(det,clu->x(),clu->y());
220  clu->setFourMomentum(mDb->getLorentzVector(xyz,clu->energy(),0.0));
221  const StLorentzVectorD& p = clu->fourMomentum();
222  //LOG_DEBUG << Form("momentum= %lf %lf %lf %lf", p.px(), p.py(), p.pz(), p.e()) << endm;
223  int ret=clusterMomentAnalysis(clu,mTowerEThreMoment); //moment analysis with default threshold
224  if(ret==kStErr) ret=clusterMomentAnalysis(clu,0.0); //Redo with 0 threshold
225  categorization(clu);
226  }
227 
228  //debug MC info
229  if(GetDebug()>=5){
230  g2t_track_st* g2ttrk=0;
231  St_g2t_track* trackTable = static_cast<St_g2t_track*>(GetDataSet("g2t_track"));
232  if(!trackTable) {
233  LOG_INFO << "g2t_track Table not found" << endm;
234  }else{
235  const int nTrk = trackTable->GetNRows();
236  LOG_INFO << Form("g2t_track table has %d tracks",nTrk) << endm;
237  if(nTrk>0){
238  g2ttrk = trackTable->GetTable();
239  if(!g2ttrk){
240  LOG_INFO << "g2t_track GetTable failed" << endm;
241  }
242  }
243  }
244  if(g2ttrk){
245  int ntrk=0;
246  float frc=0;
247  int nh = hits.size();
248  for(int i=0; i<nh; i++){
249  StFcsHit* hit=hits[i];
250  const g2t_track_st* trk = mDb->getParentG2tTrack(hit,g2ttrk,frc,ntrk);
251  //const g2t_track_st* trk=0;
252  //std::tie(trk,frc,ntrk) = mDb->getParentG2tTrack(hit,g2ttrk);
253  LOG_INFO << Form("Det=%1d Id=%3d E=%8.3f Parent Id=%4d Pid=%4d E=%8.3f Frc=%6.3f N=%d",
254  det,hit->id(),hit->energy(),trk->id,trk->ge_pid,trk->e,frc,ntrk)<<endm;
255  const g2t_track_st* ptrk = mDb->getPrimaryG2tTrack(hit,g2ttrk,frc,ntrk);
256  LOG_INFO << Form("Det=%1d Id=%3d E=%8.3f Primary Id=%4d Pid=%4d E=%8.3f Frc=%6.3f N=%d",
257  det,hit->id(),hit->energy(),ptrk->id,ptrk->ge_pid,ptrk->e,frc,ntrk)<<endm;
258  }
259  int nc = clusters.size();
260  for(int j=0; j<nc; j++){
261  StFcsCluster* clu=clusters[j];
262  const g2t_track_st* trk = mDb->getParentG2tTrack(clu,g2ttrk,frc,ntrk);
263  LOG_INFO << Form("Det=%1d C#=%3d E=%8.3f Parent Id=%4d Pid=%4d E=%8.3f Frc=%6.3f N=%d",
264  det,j,clu->energy(),trk->id,trk->ge_pid,trk->e,frc,ntrk)<<endm;
265  const g2t_track_st* ptrk = mDb->getPrimaryG2tTrack(clu,g2ttrk,frc,ntrk);
266  LOG_INFO << Form("Det=%1d C#=%3d E=%8.3f Primary Id=%4d Pid=%4d E=%8.3f Frc=%6.3f N=%d",
267  det,j,clu->energy(),ptrk->id,ptrk->ge_pid,ptrk->e,frc,ntrk)<<endm;
268  }
269  }
270  }
271 
272  return kStOk;
273 }
274 
275 //check if the hit is next to hits in the cluster
276 //returning energy of highest tower in cluster which is neighbor to the hit
277 float StFcsClusterMaker::isNeighbor(StFcsHit* hit, StFcsCluster* clu){
278  int nhit = clu->hits().size();
279  float ne=0.0;
280  for(int i=0; i<nhit; i++){
281  StFcsHit* h=clu->hits()[i];
282  int ehp = mDb->ecalHcalPres(h->detectorId());
283  float thr=1.01;
284  if(ehp==0) thr=mNeighborDistance_Ecal;
285  if(ehp==1) thr=mNeighborDistance_Hcal;
286  float d = distance(hit,h);
287  float e = h->energy();
288  if(d < thr && e>ne) ne=e;
289  }
290  return ne;
291 }
292 
293 // distance between 2 hits in cell unit
294 float StFcsClusterMaker::distance(StFcsHit* hit1, StFcsHit* hit2){
295  int det1=hit1->detectorId();
296  int det2=hit2->detectorId();
297  if(det1 != det2) return 999.0;
298  float x1,x2,y1,y2;
299  mDb->getLocalXYinCell(hit1,x1,y1);
300  mDb->getLocalXYinCell(hit2,x2,y2);
301  float dx=x1-x2;
302  float dy=y1-y2;
303  return sqrt(dx*dx + dy*dy);
304 }
305 
306 // distance between a hit and cluster center
307 float StFcsClusterMaker::distance(StFcsHit* hit, StFcsCluster* clu){
308  int det1=hit->detectorId();
309  int det2=clu->detectorId();
310  if(det1 != det2) return 999.0;
311  float x,y;
312  mDb->getLocalXYinCell(hit,x,y);
313  float dx = x - clu->x();
314  float dy = y - clu->y();
315  return sqrt(dx*dx + dy*dy);
316 }
317 
318 //Fill in number of towers, cluster energy and weighted mean x/y [local grid] based on hits collection
319 void StFcsClusterMaker::updateCluster(StFcsCluster* clu){
320  int nhit=clu->hits().size();
321  double etot=0.0, xe=0.0, ye=0.0;
322  double wtot=0.0, xw=0.0, yw=0.0;
323  for(int i=0; i<nhit; i++){
324  StFcsHit* hit=clu->hits()[i];
325  float x,y;
326  mDb->getLocalXYinCell(hit,x,y);
327  double e= hit->energy();
328  double w=log(e + 1.0 - mTowerEThreMoment);
329  if(w<0.0) w=0.0;
330  etot+= e; xe += x*e; ye += y*e;
331  wtot+= w; xw += x*w; yw += y*w;
332  }
333  clu->setNTowers(nhit);
334  clu->setEnergy(etot);
335  if(wtot>0.0){
336  clu->setX(xw/wtot);
337  clu->setY(yw/wtot);
338  }else if(etot>0.0) {
339  clu->setX(xe/etot);
340  clu->setY(ye/etot);
341  }else{
342  clu->setX(0.0);
343  clu->setY(0.0);
344  }
345 }
346 
347 // perform cluster moment analysis, return kStErr if no tower found above ecut
348 int StFcsClusterMaker::clusterMomentAnalysis(StFcsCluster* clu, float ecut){
349  int nhit=clu->hits().size();
350  double wtot=0.0, xx=0.0, yy=0.0, sx=0.0, sy=0.0, sxy=0.0;
351  for(int i=0; i<nhit; i++){
352  StFcsHit* hit=clu->hits()[i];
353  float xh,yh;
354  mDb->getLocalXYinCell(hit,xh,yh);
355  double e= hit->energy();
356  double w=log(e + 1.0 - ecut);
357  if(w>0.0){
358  wtot+= w;
359  xx += xh*w;
360  yy += yh*w;
361  sx += xh*xh*w;
362  sy += yh*yh*w;
363  sxy+= xh*yh*w;
364  }
365  }
366  if(wtot<=0.0){ //cluster has no tower above ecoff
367  if(ecut>0.0){
368  return kStErr; // return error so one can re-try with lower threshold
369  }else{ //even with ecut=0.0, no energy!?
370  clu->setSigmaMin(0.0);
371  clu->setSigmaMax(0.0);
372  return kStOK;
373  }
374  }else{
375  double x = xx/wtot;
376  double y = yy/wtot;
377  double sigx = sqrt(fabs(sx / wtot - std::pow(x, 2.0)));
378  double sigy = sqrt(fabs(sy / wtot - std::pow(y, 2.0)));
379  double sigxy = sxy/wtot - x*y;
380  double dsig2 = sigx*sigx - sigy*sigy;
381  double aA = sqrt(dsig2 * dsig2 + 4.0 * sigxy * sigxy) + dsig2;
382  double bB = 2.0 * sigxy;
383  if (sigxy < 1e-10 && aA < 1e-10) {
384  bB = sqrt(dsig2 * dsig2 + 4.0 * sigxy * sigxy) - dsig2;
385  aA = 2.0 * sigxy;
386  }
387  double theta = atan2(bB, aA);
388  while (theta > M_PI / 2.0) {
389  theta -= M_PI;
390  }
391  while (theta < -M_PI / 2.0) {
392  theta += M_PI;
393  }
394  clu->setTheta(theta);
395  clu->setSigmaMin(getSigma(clu, theta, ecut));
396  clu->setSigmaMax(getSigma(clu, theta - M_PI/2.0, ecut));
397  }
398  return kStOK;
399 }
400 
401 float StFcsClusterMaker::getSigma(StFcsCluster* clu, double theta, float ecut){
402  double sigma = 0;
403  // 2-d vector vaxis define the axis
404  TVector2 vaxis(cos(theta), sin(theta));
405  double wtot =0.0;
406  int nhit=clu->hits().size();
407  for (int i=0; i<nhit; i++){ // loop over all hits in cluster
408  StFcsHit* hit = clu->hits()[i];
409  int det=hit->detectorId();
410  float x,y;
411  mDb->getLocalXYinCell(hit,x,y);
412  // the 2-d vector from the "center" of cluster to tower
413  TVector2 v1(x - clu->x(), y - clu->y());
414  // perpendicular distance to the axis = length of the component of vector
415  // "v1" that is norm to "vaxis"
416  double dis = (v1.Norm(vaxis)).Mod();
417  // contribution to sigma
418  double w = log(hit->energy() + 1.0 - ecut);
419  if(w>=0.0){
420  wtot += w;
421  sigma += w * dis * dis;
422  }
423  }
424  return wtot > 0.0 ? (float)sqrt(sigma / wtot) : 0.0;
425 }
426 
427 void StFcsClusterMaker::categorization(StFcsCluster* cluster){
428  if(cluster->nTowers() < 5){
429  cluster->setCategory(1);
430  }else{
431  const double sigma=cluster->sigmaMax();
432  const double e =cluster->energy();
433  if(sigma > 1/2.5 + 0.003*e + 7.0/e){
434  cluster->setCategory(2);
435  }else if(sigma < 1/2.1 - 0.001*e + 2.0/e){
436  cluster->setCategory(1);
437  }else{
438  cluster->setCategory(0);
439  }
440  }
441 }
StLorentzVectorD getLorentzVector(const StThreeVectorD &xyz, float energy, float zVertex=0.0)
Get get 4 vector assuing m=0 and taking beamline from DB.
Definition: StFcsDb.cxx:895
void Clear(Option_t *option="")
User defined functions.
virtual void Clear(Option_t *option="")
User defined functions.
Definition: StMaker.cxx:634
const g2t_track_st * getParentG2tTrack(StFcsHit *h, g2t_track_st *g2ttrk, float &fraction, int &ntrk, unsigned int order=0)
Definition: StFcsDb.cxx:2191
StThreeVectorD getStarXYZfromColumnRow(int det, float col, float row, float FcsZ=-1.0) const
get the STAR frame cooridnates from other way
Definition: StFcsDb.cxx:863
void getLocalXYinCell(StFcsHit *hit, float &x, float &y) const
getting XY in local cell coordinate
Definition: StFcsDb.cxx:652
int ecalHcalPres(int det) const
Ecal North=0, Ecal South=1, Hcal North=2, Hcal South=3, Pres=4/5.
Definition: StFcsDb.cxx:423
Definition: Stypes.h:42
Definition: Stypes.h:40
Definition: Stypes.h:44
Definition: Stypes.h:41