StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StPmdClustering.cxx
1 /***********************************************************
2  *
3  * $Id: StPmdClustering.cxx,v 1.30 2010/05/31 22:13:02 rashmi Exp $
4  *
5  * Author: based on original routine written by S. C. Phatak.
6  *
7  ***********************************************************
8  *
9  * Description: Base class for PMD clusters
10  *
11  ***********************************************************
12  * ***
13  * 2004/06/24 : Dipak : findCpvClusters() function removed as because
14  * we are doing clustering same way for both the planes
15  * 2004/05/05 : Dipak: New Centroid calculation algorithm has been
16  * implemented as corrected by Prof. Phatak. A new function
17  * 'CentroidCal()' has been put in place of 'gaussfit()'.
18  **
19  * $Log: StPmdClustering.cxx,v $
20  * Revision 1.30 2010/05/31 22:13:02 rashmi
21  * New Clustering routine uses Edep() instead of Adc() of PmdHit; mOptCalibrate Flag set to kTRUE
22  *
23  * Revision 1.29 2010/05/29 00:47:10 rashmi
24  * Added a call to new clustering routines in StPmdClustering
25  *
26  * Revision 1.28 2010/04/20 23:08:51 rashmi
27  * removed refinedcluidet2.dat; no refined clusters in 2010 data
28  *
29  * Revision 1.27 2010/04/15 06:52:28 rashmi
30  * Clustering with option to turn calibration refineclustering on/off
31  *
32  * Revision 1.26 2007/11/02 11:00:06 rashmi
33  * Applying hitcalibration; eta,phi wrt primary vertex
34  *
35  * Revision 1.25 2007/08/31 10:52:30 rashmi
36  * Setting cutoff from SetAdcCutOff(); default cutoff=7
37  *
38  * Revision 1.24 2007/04/26 04:11:11 perev
39  * Remove StBFChain dependency
40  *
41  * Revision 1.23 2005/06/09 19:43:54 perev
42  * Avoid FloatPointException
43  *
44  * Revision 1.22 2004/11/15 23:35:28 subhasis
45  * Refs in centroidCal initialised to solve valgrind error
46  *
47  * Revision 1.20 2004/09/22 19:24:55 perev
48  * Leak fixed + mess with i,j indexes
49  *
50  * Revision 1.19 2004/09/03 14:31:22 subhasis
51  * memset, memcpy used (Gene's suggesstion and order() changed
52  *
53  * Revision 1.18 2004/08/01 06:41:31 subhasis
54  * nclust limit put to <200
55  *
56  * Revision 1.17 2004/07/26 12:01:31 subhasis
57  * sigmaL, sigmaS stored in one place till StPhmdCluster.h modified
58  *
59  * Revision 1.16 2004/07/21 13:02:31 subhasis
60  * refclust called only when incr <2000
61  *
62  * Revision 1.15 2004/07/19 13:23:34 subhasis
63  * checks applied on clust_cell dimension
64  *
65  * Revision 1.14 2004/07/15 13:32:41 subhasis
66  * clust dimension changed
67  *
68  * Revision 1.13 2004/06/29 07:35:25 subhasis
69  * const.h dropped
70  *
71  * Revision 1.12 2004/06/29 07:13:02 subhasis
72  * limit of clust_ fixed
73  *
74  * Revision 1.11 2004/06/24 13:46:07 subhasis
75  * several changes in clustering code
76  *
77  * Revision 1.7 2003/10/23 04:24:14 perev
78  * Stiostream again
79  *
80  * Revision 1.6 2003/10/14 07:23:18 subhasis
81  * Dipak's changes on centroid, correct insure warning
82  *
83  *
84  * Revision 1.5 2003/06/24 17:50:10 Dipak,Tapan
85  * 'Gaussfit()' function corrected as done by Prof. Phatak,
86  * for better cluster centroid determination and
87  * initialization of ncl[i] has been corrected.
88  * New condtions introduced in refclust() function.
89  *
90  *
91  * Revision 1.4 2003/05/29 13:11:51 subhasis
92  * lev1, lev2 dimension increased from 20 to 50
93  *
94  * Revision 1.3 2003/05/14 10:49:12 subhasis
95  * CPV clustering added
96  *
97  * Revision 1.2 2003/05/14 10:21:05 Dipak
98  * Clustering for CPV plane implemented same as PMD plane
99  ***********************************************************/
100 
101 #include<Stiostream.h>
102 #include<assert.h>
103 #include<math.h>
104 #include"TROOT.h"
105 #include<TRandom.h>
106 #include<TBrowser.h>
107 #include<TPad.h>
108 #include<StMessMgr.h>
109 #include<TFile.h>
110 //#include "StConstants.hh"
111 
112 #include <TTableSorter.h>
113 
114 #include "StPmdUtil/StPmdGeom.h"
115 #include "StPmdUtil/StPmdHit.h"
116 #include "StPmdClustering.h"
117 #include "StPmdClusterMaker.h"
118 #include "StPmdUtil/StPmdClusterCollection.h"
119 #include "StPmdUtil/StPmdCluster.h"
120 #include "StPmdUtil/StPmdModule.h"
121 #include "StPmdUtil/StPmdDetector.h"
122 
123 #include "StThreeVectorD.hh"
124 #include "StHelixD.hh"
125 #include "StPhysicalHelixD.hh"
126 #include "StThreeVector.hh"
127 #include "StHelix.hh"
128 #include "SystemOfUnits.h"
129 #include "StEventTypes.h"
130 #include "StEvent.h"
131 
132 
133 
134 ClassImp(StPmdClustering)
135 
136  Double_t d1[96][72],clusters[6][6912], coord[2][96][72];
137 Double_t crd_org[2][96][72];
138 
139 Int_t iord[2][6912], infocl[2][96][72], inford[3][6912], clno;
140 
141 const Double_t pi=3.141592653, sqrth=sqrt(3.)/2.;
142 const Int_t nmx = 6912;
143 Float_t cell_frac[200][2000];
144 //ofstream fout1("cluster.dat");
145 
146 //ofstream fout("refinedcluidet2.dat");
147 
148 StPmdGeom *geom=new StPmdGeom();
149 //-------------------------------------------------
151 
152 StPmdClustering::StPmdClustering(StPmdDetector *pmd_det, StPmdDetector *cpv_det):StPmdAbsClustering(pmd_det,cpv_det){
153  m_pmd_det=pmd_det;
154  m_cpv_det=cpv_det;
155  SetAdcCutOff(7.0);
156  mVertexPos = StThreeVectorF(0.,0.,0.);
157  mOptCalibrate = kTRUE; // calibration is on as default . Set it to kFALSE if calibration not desired
158  mOptSimulate = kFALSE; // working with real data in default. Set it to kTRUE if working with simulated data.
159  mOptRefineCluster = kTRUE; // refine clustering on as default
160  // cout<<" Calibration is "<<mOptCalibrate<<" Simulated data ="<<mOptSimulate<<" Refined Clustering ="<<mOptRefineCluster<<endl;
161 }
162 
163 //------------------------------
164 StPmdClustering::~StPmdClustering()
165 {
166 }
167 //---------------------------------
169 void StPmdClustering::findPmdClusters(StPmdDetector *mdet)
170 {
171  cout<<"cutoff="<<cutoff<<endl;
172  if(mdet)
173  {
174  StPmdClusterCollection * pmdclus = new StPmdClusterCollection();
175  mdet->setCluster(pmdclus);
176 
177  Int_t i, i1, i2, j,xpad, ypad,nmx1, incr,idet;
178  Int_t gsuper;
179  Double_t edep, ave;
180 
181  for(i=0; i<96; i++){
182  for(j=0;j<72;j++){
183  coord[0][i][j]=i+j/2.; coord[1][i][j]=sqrth*j;
184  crd_org[0][i][j]=i; crd_org[1][i][j]=j;
185  }
186  }
187  i=0;
188  for(Int_t id=1;id<=12;id++)
189  {
190 
192  memset(d1[0],0,96*72*sizeof(Double_t));
193 
194  StPmdModule * pmd_mod=mdet->module(id);
195 
196 
197  if(mdet->module_hit(id)>0)
198  {
199 
200  Int_t nmh=mdet->module_hit(id);
201  // cout<<" number of hits in module "<<id<<" is ="<<nmh<<endl;
202  TIter next(pmd_mod->Hits());
203  StPmdHit *spmcl;
204  for(Int_t im=0; im<nmh; im++)
205  {
206 
207  spmcl = (StPmdHit*)next();
208  if(spmcl)
209  {
210  ypad=spmcl->Row();
211  xpad=spmcl->Column();
212  //edep=spmcl->Edep(); //! Cell Edep in keV
213  edep=spmcl->Adc();
214  if(mOptSimulate==kTRUE){edep = spmcl->Edep();}
215  gsuper = spmcl->Gsuper();
216  idet=spmcl->SubDetector();
217 
218  xpad = xpad -1; // for using as array parameter make xpad = 0
219  ypad = ypad -1; // for using as array parameter make ypad = 0
220  // 11th Oct'07 : application of gain factors in ClusterMaker
221  if(mOptCalibrate==kTRUE){
222 
223  Float_t cellgain=spmcl->GainCell();
224  Float_t smchaingain=spmcl->GainSmChain();
225  Float_t cellstatus=spmcl->CellStatus();
226  Float_t finalfactor = cellgain*smchaingain*cellstatus;
227  if(finalfactor>0)edep/=finalfactor;
228  if(finalfactor<=0)edep=0;
229  }
230  //if(idet==1)fout1<<gsuper-1<<" "<<xpad<<" "<<ypad<<" "<<cellgain<<" "<<smchaingain<<" "<<cellstatus<<endl;
231 
232  // fout1<<"gain factors "<<gsuper<<" "<<ypad+1<<" "<<xpad+1<<" "<<edep<<" "<<cellgain<<" "<<smchaingain<<" "<<cellstatus<<" "<<finalfactor<<endl;
233  // end of application of gain factors
234 
235  d1[xpad][ypad]=d1[xpad][ypad]+edep;
236  }
237  }
238 
239  order(idet);
240  /*
241  ave=0.; nmx1=-1;
242  for(Int_t jj=0;jj<nmx; jj++)
243  {
244  i1=iord[0][jj]; i2=iord[1][jj];
245  if(d1[i1][i2] > 0.){nmx1=nmx1+1;ave=ave+d1[i1][i2];}
246  }
247  ave=ave/nmx1;
248  */
249  ave=0.; nmx1=0;
250  for(Int_t jj=0;jj<nmx; jj++)
251  {
252  i1=iord[0][jj]; i2=iord[1][jj];
253  if(d1[i1][i2] > 0.){nmx1=nmx1+1;ave=ave+d1[i1][i2];}
254  }
255 
256  if(nmx1>0){
257  ave=ave/nmx1; //Earlier averaging was being done using no of max one less
258  nmx1--; // this is done, assuming later he needs one less in nmx1 : Sub
259  }
260 
262 
263  // This is now set in Init() or
264  //explicitly by hand from StPmdClusterMaker
265  //AFTER instantiating.
266  // cutoff=0.4; //! cutoff(in KeV) is the threshold above which value the data is analysed.
267 
268  /* crude clusters. superclusters are formed. These are separated from each other by cells having edep smaller than cutoff. */
269  incr=crclust(ave, cutoff, nmx1, idet);
270 
271  arrange(incr);
272 
273  // cout<<"incr="<<incr<<endl;
274  if(incr<2000)refclust(mdet,incr, id, idet,pmdclus);
275  }
276  }
277  }
278 }
279 //-----------------------------------------
280 
281 void StPmdClustering::printclust(Int_t i,Int_t m, StPmdCluster* pclust)
282 {
283 
286  Float_t x0, y0, x, y,xc,yc,zc, cells;
287 
288  xc = clusters[0][m];yc = clusters[1][m];
289  zc = clusters[2][m];
290  cells = clusters[3][m];
291 
292  Float_t clusigmaL = clusters[4][m];
293  Float_t clusigmaS = clusters[5][m];
294 
295  Float_t cluedep=zc; y0 = yc/sqrth; x0 = xc - y0/2.;
296  Float_t clueta,cluphi;
297 
298  // modified by RR on 1 Nov 2007 to get eta & phi wrt to vertex
299  // cout<<" StPmdClustering: vertex="<<mVertexPos.x()<<","<<mVertexPos.y()<<","<<mVertexPos.z()<<endl;
300  if(mVertexPos.x()==0 && mVertexPos.y()==0 && mVertexPos.z()==0) {
301  // This loop is for backcompatibility
302  geom->DetCell_xy(i,y0+1,x0+1,x,y,clueta,cluphi);
303  // cout<<"clueta,phi="<<clueta<<","<<cluphi<<endl;
304  }else{
305  // Gives x,y of cluster in STAR Coordinates using x0,y0 wrt SM corner
306  geom->DetCell_xy(i,y0+1,x0+1,x,y,clueta,cluphi);
307  Float_t xwrtv = x-mVertexPos.x();
308  Float_t ywrtv = y-mVertexPos.y();
309  // StPmdGeom returns PMDZ as a positive number
310  // But in STAR coordinated the PMDZ in negative
311  Float_t zreal = -geom->GetPmdZ();
312  Float_t zwrtv = zreal-mVertexPos.z();
313 
314  // cout<<"x,y,z real="<<x<<","<<y<<","<<zreal<<endl;
315  // cout<<"x,y,z wrtvertex="<<xwrtv<<","<<ywrtv<<","<<zwrtv<<endl;
316  Cluster_Eta_Phi(xwrtv,ywrtv,zwrtv,clueta,cluphi);
317  // cout<<"clueta,phi="<<clueta<<","<<cluphi<<endl;
318  }
319 
321 
322  pclust->setCluX(x);
323  pclust->setCluY(y);
324  pclust->setModule(i);
325  pclust->setCluEdep(cluedep);
326  pclust->setCluEta(clueta);
327  pclust->setCluPhi(cluphi);
328  pclust->setCluSigmaL(clusigmaL);
329  pclust->setCluSigmaS(clusigmaS);
330  pclust->setNumofMems(cells);
331 
332 }
333 //------------------------------------------------------
334 
336 void StPmdClustering::order(Int_t idet)
337 {
338  Double_t d[nmx];
339  Int_t curl=0;
340  Int_t curh=nmx-1;
341  for (int i1=0; i1<96; i1++) {
342  for(int i2=0; i2 < 72; i2++){
343  if (d1[i1][i2] > 0.) {
344  // Insert where it belongs in descending order
345  int j = 0;
346  while ((j<curl) && (d[j]>=d1[i1][i2])) j++; // Find insertion point
347  for (int k=curl; k>j; k--) { // Shift other points forward
348  int l = k-1;
349  d[k] = d[l];
350  iord[0][k] = iord[0][l];
351  iord[1][k] = iord[1][l];
352  }
353  // Insert data
354  d[j] = d1[i1][i2];
355  iord[0][j] = i1;
356  iord[1][j] = i2;
357  curl++;
358  } else {
359  // Insert at back end
360  iord[0][curh] = i1;
361  iord[1][curh] = i2;
362  curh--;
363  }
364  }
365  }
366  }
367 
368 //---------------------------------------------
369 
371 void StPmdClustering::arrange(Int_t incr)
372 {
373  Int_t i, j, k, i1, itest, ihld1, ihld2, ihld3;
374  i1=-1;
375  for(i=0; i<96; i++){
376  for(j=0; j<72; j++){
377  if(infocl[0][i][j] == 1){i1=i1+1;
378  inford[0][i1]=infocl[1][i][j]; inford[1][i1]=i; inford[2][i1]=j;
379  }
380  if(infocl[0][i][j] == 2){i1=i1+1;
381  inford[0][i1]=infocl[1][i][j]; inford[1][i1]=i; inford[2][i1]=j;
382  }
383  }
384  }
385  for(i=1; i < incr; i++){
386  itest=0; ihld1=inford[0][i]; ihld2=inford[1][i]; ihld3=inford[2][i];
387  for(j=0; j<i; j++){
388  if(itest == 0 && inford[0][j] > ihld1){
389  itest=1;
390  for(k=i-1; k>=j; k--){
391  inford[0][k+1]=inford[0][k]; inford[1][k+1]=inford[1][k];
392  inford[2][k+1]=inford[2][k];
393  }
394  inford[0][j]=ihld1; inford[1][j]=ihld2; inford[2][j]=ihld3;
395  }
396  }
397  }
398 }
399 
400 //-------------------------------------------------
405 void StPmdClustering::refclust(StPmdDetector* m_pmd_det0,Int_t incr, Int_t supmod, Int_t idet,StPmdClusterCollection *pmdclus)
406 {
407  Int_t clno, i, j, k, i1, i2, id, icl, ncl[2000], iordR[2000], itest, ihld;
408  Int_t ig, nsupcl;
409  Double_t x1, y1, z1, x2, y2, z2, rr;
410  Double_t x[2000], y[2000], z[2000];
411  Double_t x_org[2000], y_org[2000];
412  Double_t xc[2000], yc[2000], zc[2000], d[96][72],rcl[2000],rcs[2000],cells[2000];
413  // Double_t xcl[2000], ycl[2000],clust_cell[100][2000];
414 
415 
416 
420 
421  clno=-1;
422  nsupcl=-1;
423 
424  for(i1=0; i1<96; i1++){
425  for(i2=0; i2<72; i2++){
426  d[i1][i2]=d1[i1][i2];
427  }
428  }
429  //subhasis 15/11/04 added to initialze arrays for centroidCalc refs
430  for(i=0; i<2000; i++){
431  x[i]=0.;
432  y[i]=0.;
433  z[i]=0.;
434  xc[i]=0.;
435  yc[i]=0.;
436  zc[i]=0.;
437  }
439 
440  for(i=0; i<2000; i++){
441  ncl[i]=-1; // ncl[i] --> initialization starts from '-1'
442  rcl[i] = 0.; // initialization of rcs and rcl
443  rcs[i] = 0.;
444  cells[i] = 0.;
445  }
446  for(i=0; i<incr; i++)
447  {
448  if(inford[0][i] != nsupcl)
449  {
450  nsupcl=nsupcl+1;
451  }
452  ncl[nsupcl]=ncl[nsupcl]+1;
453  }
454  id=-1;
455  icl=-1;
456 
457  for(i=0; i<=nsupcl; i++)
458  {
459  Int_t countr = 0;
460  if(ncl[i] == 0){
461  id=id+1; icl=icl+1;
462 
463  countr +=1;
465 
466  clno=clno+1; i1=inford[1][id]; i2=inford[2][id];
467 
468  clusters[0][clno]=coord[0][i1][i2];
469  clusters[1][clno]=coord[1][i1][i2];
470  clusters[2][clno]=d[i1][i2];
471  clusters[3][clno]=1.;
472  clusters[4][clno]=0.; // for single cell put sigma as 0.
473  clusters[5][clno]=0.; // for single cell put sigma as 0.
474 
475  //Create the StPmdCluster and add to the ClusterCollection
476  StPmdCluster *pclust = new StPmdCluster();
477  pmdclus->addCluster(pclust);
478 
479  printclust(supmod,clno,pclust);
480 
481  //Get StPmdHit* corresponding to the co-ordinate of single cell
482  StPmdHit* phit = GetHit(m_pmd_det0,supmod,crd_org[0][i1][i2],crd_org[1][i1][i2]);
483  if(phit)pclust->addHitCollection(phit);
484 
485  }
486  /*
487  // Set up this block to work even if the number of cells is>2 RR
488  else if((ncl[i] == 1 && mOptRefineCluster==kTRUE) || (ncl[i]>0 && mOptRefineCluster==kFALSE))
489  { // Super cluster having 'Two' cells when refined clustering is ON
490  // Super clusters having more than one cell when refined is OFF
491  //Create the StPmdCluster and add to the ClusterCollection
492  */
493  else if(ncl[i]==1)
494  {
495  StPmdCluster *pclust = new StPmdCluster();
496  pmdclus->addCluster(pclust);
497  id=id+1; icl=icl+1;
498  clno=clno+1; i1=inford[1][id]; i2=inford[2][id];
499  x1=coord[0][i1][i2]; y1=coord[1][i1][i2]; z1=d[i1][i2];
500 
501  countr +=1;
502  //Get StPmdHit* corresponding to the co-ordinate of first cell
503  StPmdHit* phit = GetHit(m_pmd_det0,supmod,crd_org[0][i1][i2],crd_org[1][i1][i2]);
504  if(phit)pclust->addHitCollection(phit);
505 
506  id=id+1; i1=inford[1][id]; i2=inford[2][id];
507  x2=coord[0][i1][i2]; y2=coord[1][i1][i2]; z2=d[i1][i2];
508 
509 
510  //Get StPmdHit* corresponding to the co-ordinate of second cell
511  phit = GetHit(m_pmd_det0,supmod,crd_org[0][i1][i2],crd_org[1][i1][i2]);
512  if(phit)pclust->addHitCollection(phit);
513 
514  /***** Adding for calculating Spread of the cluster *****/
515  Double_t xcell[2],ycell[2],zcell[2],xcl[2000],ycl[2000];
516  xcell[0] = x1;
517  ycell[0] = y1;
518  zcell[0] = z1;
519  xcell[1] = x2;
520  ycell[1] = y2;
521  zcell[1] = z2;
522 
523  xcl[i] = (x1*z1+x2*z2)/(z1+z2);
524  ycl[i] = (y1*z1+y2*z2)/(z1+z2);
525  Double_t sumxx,sumyy,sumxy;
526  sumxx = 0.; sumyy = 0.; sumxy = 0.;
527  for(j=0; j<2; j++)
528  {
529  sumxx = sumxx + zcell[j]*(xcell[j]-xcl[i])*(xcell[j]-xcl[i])/(z1+z2);
530  sumyy = sumyy + zcell[j]*(ycell[j]-ycl[i])*(ycell[j]-ycl[i])/(z1+z2);
531  sumxy = sumxy + zcell[j]*(xcell[j]-xcl[i])*(ycell[j]-ycl[i])/(z1+z2);
532  }
533  Double_t b1 = sumxx + sumyy;
534  Double_t c1 = sumxx*sumyy - sumxy*sumxy;
535  Double_t dis = b1*b1/4.-c1;
536  if (fabs(dis) < 1e-6) dis = 0.;
537  dis = sqrt(dis);
538  Double_t r1=b1/2.+dis;
539  Double_t r2=b1/2.-dis;
540 
541 
544  if(r1 < r2)
545  {
546  clusters[4][clno] = r2; //SigmaL
547  clusters[5][clno] = r1; //SigmaS
548  }
549  else
550  {
551  clusters[4][clno] = r1; //SigmaL
552  clusters[5][clno] = r2; //SigmaS
553  }
554 
555  clusters[0][clno]=(x1*z1+x2*z2)/(z1+z2);
556  clusters[1][clno]=(y1*z1+y2*z2)/(z1+z2);
557  clusters[2][clno]=z1+z2;
558  clusters[3][clno]=2.;
559 
560  printclust(supmod,clno,pclust);
561 
562  }
563  else
564  {
565  // cout<<" big supercluster ncell="<<ncl[i]+1<<endl;
566  id=id+1; iordR[0]=0;
567  /* super-cluster of more than two cells - broken up into smaller clusters gaussian centers computed. (peaks separated by > 1 cell). Start from top */
568 
569  i1=inford[1][id]; i2=inford[2][id];
570  x[0]=coord[0][i1][i2]; y[0]=coord[1][i1][i2]; z[0]=d[i1][i2];iordR[0]=0;
571  x_org[0]=crd_org[0][i1][i2]; y_org[0]=crd_org[1][i1][i2];
572  for(j=1;j<=ncl[i];j++)
573  {
574  id=id+1;
575  i1=inford[1][id]; i2=inford[2][id];iordR[j]=j;
576  x[j]=coord[0][i1][i2]; y[j]=coord[1][i1][i2]; z[j]=d[i1][i2];
577  x_org[j]=crd_org[0][i1][i2]; y_org[j]=crd_org[1][i1][i2];
578  }
580  for(j=1;j<=ncl[i];j++)
581  {
582  itest=0; ihld=iordR[j];
583  for(i1=0;i1<j;i1++){
584  if(itest == 0 && z[iordR[i1]] < z[ihld]){
585  itest=1;
586  for(i2=j-1;i2>=i1;i2--){
587  iordR[i2+1]=iordR[i2];
588  }
589  iordR[i1]=ihld;
590  }
591  }
592  }
594  ig=0;
595  xc[ig]=x[iordR[0]]; yc[ig]=y[iordR[0]]; zc[ig]=z[iordR[0]];
596 
597  //-----------------------------------------------------------
598  // If this block is executed then the superclusters larger than
599  // 2 cells are broken into smaller clusters if the right conditions are met.
600 
601  if(mOptRefineCluster==kTRUE){
602  //Find Local maxima
603  for(j=1;j<=ncl[i];j++)
604  {
605  itest=-1; x1=x[iordR[j]]; y1=y[iordR[j]];
606  for(k=0;k<=ig;k++)
607  {
608  x2=xc[k]; y2=yc[k]; rr=Dist(x1,y1,x2,y2);
609 
610  if( rr >= 1.1 && rr < 1.8 && z[iordR[j]] > zc[k]*0.30)itest=itest+1;
611  if( rr >= 1.8 && rr < 2.1 && z[iordR[j]] > zc[k]*0.15)itest=itest+1;
612  if( rr >= 2.1 && rr < 2.8 && z[iordR[j]] > zc[k]*0.05)itest=itest+1;
613  if( rr >= 2.8)itest=itest+1;
614  }
615  if(itest == ig)
616  {
617  ig=ig+1; xc[ig]=x1; yc[ig]=y1; zc[ig]=z[iordR[j]];
618  }
619  }
620  // End of finding l,ocal maxima
621  }
622 
623  //-----------------------------------------------------------------
624 
625  memset(cell_frac[0],0,sizeof(cell_frac));
626  //
627 
628  Int_t censtat=CentroidCal(ncl[i],ig,x[0],y[0],z[0],xc[0],yc[0],zc[0],rcl[0],rcs[0],cells[0]);
629  if(censtat==kStOK){
630 
631  icl=icl+ig+1;
632 
633  Float_t temp[2000];
634  Int_t take_cell[2000];
635  for(Int_t jk=0; jk<2000; jk++)
636  {
637  temp[jk]=0.;
638  take_cell[jk]=-999;
639  }
640 
641  //VP memset(temp,0,2000*sizeof(Float_t)); //it is already zeroed above
642 
643  for(Int_t pb=0;pb<=ig;pb++)
644  {
645  for(Int_t jk=0; jk<=ncl[i]; jk++)
646  {
647  if(cell_frac[pb][jk]>temp[jk])
648  {
649  take_cell[jk]=pb;
650  temp[jk]=cell_frac[pb][jk];
651  }
652  }
653  }
654 
655 
656 
657 
658 
660  for(k=0; k<=ig; k++)
661  {
662  countr +=1;
663  clno=clno+1;
664  clusters[0][clno]=xc[k];
665  clusters[1][clno]=yc[k];
666  clusters[2][clno]=zc[k];
667  clusters[3][clno]=cells[k];
668  clusters[4][clno]=rcl[k];
669  clusters[5][clno]=rcs[k];
670  if(clusters[3][clno]==1)
671  {
672  clusters[4][clno]=0.;
673  clusters[5][clno]=0.;
674  }
675 
676 
677  //subhasis , ig = no of gaussians.
678  //
679  // looping over clusters first and cells within to
680  // create StPmdCluster and attach cells to them
681 
682  StPmdCluster *pclust = new StPmdCluster();
683  pmdclus->addCluster(pclust);
684 
685  printclust(supmod,clno,pclust);
686 
687  for(Int_t jk=0; jk<=ncl[i]; jk++)
688  {// loop over all cells in supercluster
689 
690  // dist=Dist(x[jk], y[jk], xc[k], yc[k]);
691  // if(dist < 2.8){ //changed from 2.1 to 2.8 : dipak
692  // StPmdHit* phit = GetHit(m_pmd_det0,supmod,x_org[jk],y_org[jk]);
693  // if(phit)pclust->addHitCollection(phit);
694  // attach the hits
695  // } // if dist loop
696 
697  if(take_cell[jk]==k)
698  {
699  StPmdHit* phit = GetHit(m_pmd_det0,supmod,x_org[jk],y_org[jk]);
700  if(phit)pclust->addHitCollection(phit);
701  }
702 
703  } //for loop 'jk
704  } //for 'k' loop
705  }//censtat check
706  } // 'else' loop
707 
708  // if(idet == 2)fout << i <<"\t" << nsupcl << "\t" << clno << "\t" << ncl[i] << "\t" << countr <<"\t" << supmod << "\t" << idet << endl; // idet = 2 for cpv
709 
710  }// for loop 'i<nsupcl'
711  // cout<<"CLUSTER NUMBER IS "<<clno<<"supmod**"<<supmod<<endl;
712 
713 }
714 //-------------------------------------------------------
715 
716 
717 Int_t StPmdClustering::CentroidCal(Int_t ncell,Int_t nclust,Double_t &x,
718  Double_t &y,Double_t &z,Double_t &xc,
719  Double_t &yc,Double_t &zc,
720  Double_t &rcl,Double_t &rcs,Double_t &cells)
721 {
722  Int_t i1,i2;
723  Int_t cluster[2000][10];
724  Double_t sum,sumx,sum1,sumy,sumxy,sumxx,sumyy;
725  Double_t rr,x1,y1,x2,y2,b,c,r1,r2;
726  Double_t xx[2000],yy[2000],zz[2000];
727  Double_t xxc[2000],yyc[2000];
728  Double_t zzct[2000],cellsc[2000];
729  Double_t str[2000],str1[2000], cln[2000];
730  Double_t xcl[2000], ycl[2000],clust_cell[200][2000];
731  //Initialisation part starts
732  double dis;
733  if(nclust>=200 || ncell >=2000){
734  cout<<"Number of cluster of Ncell crosses limit "<<nclust<<" "<<ncell<<endl;
735  return kStWarn;
736  }
737 
738  for(int i=0;i<=nclust;i++)
739  {
740  xxc[i] = *(&xc+i);
741  yyc[i] = *(&yc+i);
742  cellsc[i] = 0.;
743  zzct[i] = 0.;
744  }
745  for(int i=0;i<=ncell;i++)
746  {
747  xx[i] = *(&x+i);
748  yy[i] = *(&y+i);
749  zz[i] = *(&z+i);
750  }
751 
752 memset(clust_cell[0],0,200*2000*sizeof(Double_t));
753 memset(cell_frac[0],0,200*2000*sizeof(Float_t));
754 memset(cluster[0],0,2000*10*sizeof(Int_t));
755 
756 
757  //If there is more than one local maxima
758  if(nclust>0)
759  {
760  for(int i=0;i<=ncell;i++)
761  {
762  x1 = xx[i];
763  y1 = yy[i];
764  cluster[i][0] = 0;
765  //Checking cells shared betn several clusters.
766  for(int j=0;j<=nclust;j++) // For checking the first layer neighbours
767  {
768  x2 = xxc[j];
769  y2 = yyc[j];
770  rr = Dist(x1,y1,x2,y2);
771  if(rr < 1.) //Check if cluster is one cell unit from the shared cell
772  {
773  cluster[i][0]++;
774  i1 = cluster[i][0];
775  cluster[i][i1] = j;
776  }
777  }
778 
779  if(cluster[i][0] ==0){ //For checking the second layer neighbours
780  for(int j=0;j<=nclust;j++){
781  x2=xxc[j];
782  y2=yyc[j];
783  rr=Dist(x1,y1,x2,y2);
784  // if(rr>=1.73&&rr<2.1)
785  if(rr<= sqrt(3.))
786  {
787  cluster[i][0]++;
788  i1=cluster[i][0];
789  cluster[i][i1] = j;
790  }
791  } //for loop 'nclust'
792  } //if loop 'cluster[i][0]
793  if(cluster[i][0] ==0){ //For checking the third layer neighbours
794  for(int j=0;j<=nclust;j++){
795  x2=xxc[j];
796  y2=yyc[j];
797  rr=Dist(x1,y1,x2,y2);
798  if(rr<= 2.)
799  {
800  cluster[i][0]++;
801  i1=cluster[i][0];
802  cluster[i][i1] = j;
803  }
804  }
805  }
806 
807  if(cluster[i][0] ==0){ //For checking the fourth layer neighbours
808  for(int j=0;j<=nclust;j++){
809  x2=xxc[j];
810  y2=yyc[j];
811  rr=Dist(x1,y1,x2,y2);
812  if(rr<= 2.7)
813  {
814  cluster[i][0]++;
815  i1=cluster[i][0];
816  cluster[i][i1] = j;
817  }
818  } // for loop 'j'
819  } // if loop 'cluster[i][0]'
820  } // for loop 'ncell'
822 
823 
824  //Compute the cluster strength.
825 memset(str,0,2000*sizeof(Double_t));
826 memset(str1,0,2000*sizeof(Double_t));
827 
828  for(int i=0;i<=ncell;i++){
829  if(cluster[i][0]!=0){
830  i1 = cluster[i][0];
831  for(int j=1;j<=i1;j++){
832  i2 = cluster[i][j];
833  str[i2] = str[i2] + zz[i]/i1;
834  }
835  }
836  }
837 
838  for(int k=0; k<5; k++){
839  for(int i=0; i<=ncell; i++){
840  if(cluster[i][0] != 0){
841  i1=cluster[i][0];
842  sum = 0;
843  for(int j=1;j<=i1;j++){
844  sum = sum + str[cluster[i][j]];
845  }
846  for(int j=1;j<=i1;j++){
847  i2 = cluster[i][j];
848  str1[i2] = str1[i2] + zz[i]*str[i2]/sum;
849  clust_cell[i2][i] = zz[i]*str[i2]/sum;
850  }
851  }
852  }
853  for(int j=0; j<=nclust; j++){
854  str[j]=str1[j];
855  str1[j]=0.;
856  }
857  }
858 
859 
860  for(int i=0;i<=nclust;i++){
861  sumx=0.;sumy=0.;sum=0.;sum1=0.;
862  for(int j=0;j<=ncell;j++){
863  if(clust_cell[i][j] !=0){
864 
865  sumx=sumx+clust_cell[i][j]*xx[j];
866  sumy=sumy+clust_cell[i][j]*yy[j];
867  sum=sum+clust_cell[i][j];
868  sum1=sum1+clust_cell[i][j]/zz[j];
869  //cell contribution
870  cell_frac[i][j]=clust_cell[i][j]/zz[j];
871 
872  //
873  //cout<<"icl,icell,zzi,clust_cell "<<i<<" "<<j<<" "<<zz[j]<<" "<<clust_cell[i][j]<<" "<<sum1<<endl;
874  // cout << xx[j] << " " << yy[j] << " " << clust_cell[i][j] << " " <<
875  // clust_cell[i][j]/zz[j] <<" NClust*** "<<nclust<<" "<<i<<" "<<j<<endl;
876  }
877  }
878  xcl[i]=sumx/sum; ycl[i]=sumy/sum; cln[i]=sum1;
879  // if(sum1>10)cout<<"many gaus, i,sum1, cln "<<i<<" "<<sum1<<" "<<cln[i]<<endl;
880  sumxx=0.; sumyy=0.; sumxy=0.;
881  for(int j=0; j<=ncell; j++){
882  sumxx=sumxx+clust_cell[i][j]*(xx[j]-xcl[i])*(xx[j]-xcl[i])/sum;
883  sumyy=sumyy+clust_cell[i][j]*(yy[j]-ycl[i])*(yy[j]-ycl[i])/sum;
884  sumxy=sumxy+clust_cell[i][j]*(xx[j]-xcl[i])*(yy[j]-ycl[i])/sum;
885  }
886  b=sumxx+sumyy;
887  c=sumxx*sumyy-sumxy*sumxy;
888  dis = b*b/4.-c;
889  if (fabs(dis) < 1e-6) dis = 0.;
890  dis = sqrt(dis);
891  r1=b/2.+dis;
892  r2=b/2.-dis;
893  if(r1 < r2){
894  *(&rcl+i) = r2;
895  *(&rcs+i) = r1;
896  }else{
897  *(&rcl+i) = r1;
898  *(&rcs+i) = r2;
899  }
900  //cout<< "1st Check########**** "<<neve<<" "<<xcl[i] << " " << ycl[i] << " " << str[i] <<" "<< xxc[i] << " " << yyc[i] <<" "<<r2<<" "<<r1<<" " << (sumxx) << " " << (sumyy) << " " << sumxy << " " << cln[i] << endl;
901 
902  // Final assignment
903  *(&xc + i) = xcl[i];
904  *(&yc + i) = ycl[i];
905  *(&zc + i) = str[i];
906  *(&cells + i) = cln[i];
907 
908  }
909  } //if loop 'nclust > 0'
910  else
911  {
912  // cout << " only one cluster " << endl;
913  sumx=0.; sumy=0.; sum=0.; sum1=0.;
914  int i=0;
915  for(int j=0; j<=ncell; j++)
916  {
917  sumx=sumx+zz[j]*xx[j];
918  sumy=sumy+zz[j]*yy[j];
919  sum=sum+zz[j];
920  sum1=sum1+1.;
921 // cell contribution to cluster
922  cell_frac[i][j]=1;
923  // cout << xx[j] << " " << yy[j] << " " << zz[j] << endl;
924  }
925  xcl[i]=sumx/sum; ycl[i]=sumy/sum; cln[i]=sum1; str[i]=sum;
926  // if(sum1>10)cout<<"one gaus, i,sum1, cln "<<i<<" "<<sum1<<" "<<cln[i]<<endl;
927  sumxx=0.; sumyy=0.; sumxy=0.;
928  for(int j=0; j<=ncell; j++)
929  {
930  sumxx=sumxx+zz[j]*(xx[j]-xcl[i])*(xx[j]-xcl[i])/sum;
931  sumyy=sumyy+zz[j]*(yy[j]-ycl[i])*(yy[j]-ycl[i])/sum;
932  sumxy=sumxy+zz[j]*(xx[j]-xcl[i])*(yy[j]-ycl[i])/sum;
933  }
934 
935  /**** Dipak: added for calculating 'rcl' & 'rcs' similar to above *****/
936  b=sumxx+sumyy;
937  c=sumxx*sumyy-sumxy*sumxy;
938  dis = b*b/4.-c;
939  if (fabs(dis) < 1e-6) dis = 0.;
940  dis = sqrt(dis);
941  r1=b/2.+dis;
942  r2=b/2.-dis;
943  if(r1 < r2){
944  *(&rcl+i) = r2;
945  *(&rcs+i) = r1;
946  }else{
947  *(&rcl+i) = r1;
948  *(&rcs+i) = r2;
949  }
950 
951  *(&xc + i) = xcl[i];
952  *(&yc + i) = ycl[i];
953  *(&zc + i) = str[i];
954  *(&cells + i) = cln[i];
955  // cout<<"cln,cells "<<cln[i]<<" "<<endl;
956 
957 
958  }
959  return kStOK;
960 
961 }
962 
963 //---------------------------------------------------
964 
966 Double_t StPmdClustering::Dist(Double_t x1, Double_t y1, Double_t x2, Double_t y2)
967 {
968  return sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
969 }
970 //----------------------------------------------
974 Int_t StPmdClustering::crclust(Double_t ave, Double_t cutoff, Int_t nmx1, Int_t idet)
975 {
976  Double_t xx, yy, d[96][72];
977  Int_t i,j,k,id1,id2,icl,clust[2][3000], numcell, cellcount;
978  Int_t jd1,jd2,icell;
979 
981  static Int_t neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
982 
983 
984  memcpy(d,d1,96*72*sizeof(Double_t));
985 
986  for (j=0; j < 96; j++){
987  for(k=0; k < 72; k++){
988  for (i=0; i < 2; i++){infocl[i][j][k] = 0;}
989  }
990  }
991  cellcount=0;
992  for(i=0; i < nmx; i++){
993  id1=iord[0][i]; id2=iord[1][i];
994  if(d[id1][id2] <= cutoff){infocl[0][id1][id2]=-1;}
995  }
998  icl=-1;
999  for(icell=0; icell <= nmx1; icell++){
1000 
1001  id1=iord[0][icell]; id2=iord[1][icell]; xx=id1+id2/2.; yy=sqrth*id2;
1002  if(infocl[0][id1][id2] == 0 ){
1003 
1008  icl=icl+1; numcell=0;
1009  for(i=0; i < 3000; i++){
1010  clust[0][i]=0; clust[1][i]=0;
1011  }
1012  clust[0][numcell]=id1; clust[1][numcell]=id2;
1013  infocl[0][id1][id2]=1; infocl[1][id1][id2]=icl;
1015  for(i=0; i<6; i++){
1016  jd1=id1+neibx[i]; jd2=id2+neiby[i];
1017  if( (jd1 >= 0 && jd1 < 96) && (jd2 >= 0 && jd2 < 72) &&
1018  d[jd1][jd2] > cutoff && infocl[0][jd1][jd2] == 0)
1019  {
1020  numcell=numcell+1; clust[0][numcell]=jd1; clust[1][numcell]=jd2;
1021  infocl[0][jd1][jd2]=2; infocl[1][jd1][jd2]=icl;
1022  xx=jd1+jd2/2.; yy=sqrth*jd2;
1023  }
1024  }
1025  for(i=1;i < 3000;i++){
1026  if(clust[0][i] != 0){
1027  id1=clust[0][i]; id2=clust[1][i];
1028  for(j=0; j<6 ; j++){
1029  jd1=id1+neibx[j]; jd2=id2+neiby[j];
1030  if( (jd1 >= 0 && jd1 < 96) && (jd2 >= 0 && jd2 < 72) &&
1031  d[jd1][jd2] > cutoff && infocl[0][jd1][jd2] == 0 ){
1032  infocl[0][jd1][jd2]=2; infocl[1][jd1][jd2]=icl;
1033  numcell=numcell+1;
1034  clust[0][numcell]=jd1; clust[1][numcell]=jd2;
1035  xx=jd1+jd2/2.; yy=sqrth*jd2;
1036  }
1037  }
1038  }
1039  }
1040  cellcount=cellcount+numcell+1;
1041  }
1042  }
1043  // cout<<idet<<" CRCLUST**, icl "<<icl<<" ********* "<<cellcount<<endl;
1044  return cellcount;
1045 }
1046 //------------------------------------------
1047 
1048 StPmdHit*
1049 StPmdClustering::GetHit(StPmdDetector* pdet,Int_t id,Double_t xc, Double_t yc)
1050 {
1051  Int_t xpad,ypad,super;
1052  Int_t nmh=pdet->module_hit(id);
1053  StPmdModule * pmod=pdet->module(id);
1054  TIter next(pmod->Hits());
1055  StPmdHit *spmcl;
1056  // Loop over hits for each SM
1057  for(Int_t im=0; im<nmh; im++)
1058  {
1059  spmcl = (StPmdHit*)next();
1060  if(spmcl){
1061  ypad=spmcl->Row();
1062  xpad=spmcl->Column();
1063  super = spmcl->Gsuper();
1064  if(int(yc+1)==ypad && int(xc+1)==xpad && id == super) return spmcl;
1065  }
1066  }
1067  return NULL;
1068 }
1069 //--------------------------------------------------------
1070 // This routine was copied from StPmdGeom with the difference that it uses
1071 // zreal ( PMDZ wrt vertex) instead of mzreal ( PMDZ wrt 0,0,0)
1072 void StPmdClustering::Cluster_Eta_Phi(Float_t xreal,Float_t yreal,Float_t zreal,Float_t& eta,Float_t& phi){
1073  Float_t rdist = (TMath::Sqrt(xreal*xreal + yreal*yreal))/TMath::Abs(zreal);
1074  Float_t theta = TMath::ATan(rdist);
1075  //Bedanga & pawan - added theta !=0.
1076  // The negative sign is avoided to take care of -ve PMDZ
1077  if(theta !=0.){ eta = TMath::Log(TMath::Tan(0.5*theta));}
1078  if( xreal==0) {
1079  if(yreal>0) {phi = 1.571428;}
1080  if(yreal<0) {phi = -1.571428;}
1081  }
1082  if(xreal != 0) {phi = atan2(yreal,xreal);}
1083 
1084 }
1085 
1086 //---------------------------------
1088 Bool_t StPmdClustering::findPmdClusters2(StPmdDetector *mdet)
1089 {
1090  cout<<"cutoff in findPmdClusters2="<<cutoff<<endl;
1091  //cout<<" Finding cluster for PMD detector"<<endl;
1092 
1093 
1094  // Declared in StpmdClusterin.h now
1095  // StPmdClusterCollection * pmdclus;
1096  pmdclus = new StPmdClusterCollection();
1097  mdet->setCluster(pmdclus);
1098 
1099  if(mdet){
1100 
1101 
1102  for(Int_t ism = 1; ism<=12; ism++){
1103 
1104  TClonesArray * mPmdSuperClusters = new TClonesArray("TList",1000);
1105 
1106  // Getting the pointer to the supermodule from the detector
1107  // notice that StPmdDetector->module(xx) takes xx from 1 to 12
1108  StPmdModule * pmd_mod = mdet->module(ism); //
1109  Int_t nhits = mdet->module_hit(ism);
1110  // cout<<" Number of entries in "<<ism<<" module is "<<nhits<<endl;
1111 
1112  // if no hits, go to the next supermodule
1113  if(nhits<=0) continue;
1114  // pass on the StPmdModule to make superclusters
1115  // All the superclusters are filled into the mPmdSuperClusters TClonesArray
1116  Int_t nSuperClusters = MakeSuperClusters(ism,pmd_mod,mPmdSuperClusters);
1117  // cout<<" number of superclusters on module "<<ism<<" is "<<nSuperClusters<<endl;
1118 
1119  // Going to form clusters now
1120  if (nSuperClusters==0|| nSuperClusters > nhits){
1121  cout<<"Warning!!! nent/nsuper "<<nhits<<"/"<<nSuperClusters<<endl;
1122  cout<<"The hits of this module will not be added to the outgoing TClonesArray of hits"<<endl;
1123  // return kStErr;
1124  }else{
1125  // cout<<"Going to MakeClusters()"<<endl;
1126  // Int_t gotclusters = MakeClusters(mPmdSuperClusters,pmdclus);
1127  Int_t gotclusters = MakeClusters(mPmdSuperClusters);
1128 
1129  if (gotclusters==0){
1130  cout<<" Did not get any clusters in this sm"<<endl;
1131  }
1132  // if (mOptDebug)
1133  }
1134  mPmdSuperClusters->Delete();
1135  delete mPmdSuperClusters;
1136  }
1137  }
1138  cout<<"Number of clusters in findPmdClusters2 = "<<pmdclus->Nclusters()<<endl;
1139 
1140  return kStOK;
1141 }
1142 
1143 //---------------------------------------------------------------------------------------
1144 
1145 Int_t StPmdClustering::MakeSuperClusters(Int_t modnum,StPmdModule * mpmdMod , TClonesArray* mPmdSuperClusters){
1146 
1147  Int_t Columnneigh[6] = {-1,-1,0,1,1,0};
1148  Int_t Rowneigh[6] = {0,1,1,0,-1,-1};
1149  // cout<<"In MakeSuperClusters: CutOff="<<cutoff<<endl;
1150 
1151  //Initializing the hitmap
1152 
1153  //Row and column numbers start from 0
1154  Int_t hitmap[NRowMax][NColumnMax];
1155  Int_t Flag[NRowMax][NColumnMax];
1156  for( Int_t ncellX = 0; ncellX < NRowMax; ncellX++ ){
1157  for( Int_t ncellY = 0; ncellY < NColumnMax; ncellY++ ){
1158  hitmap[ncellX][ncellY]=-1;
1159  Flag[ncellX][ncellY]=0;
1160  }
1161  }
1162  // Get pointer to the TObjArray of pmd hits in this module
1163  TObjArray * PmdHits = mpmdMod->Hits();
1164  // PmdHits->Sort();
1165  // Sorts in descending order
1166  PmdHits->Sort(0);
1167  // GetNumber of hits in the module
1168  Int_t Ncell_mod = PmdHits->GetEntries();
1169 
1170  // Just accountancy
1171  Int_t NWithin=0;
1172  Int_t NAccCut=0;
1173  // Filling hitmap will all the hits. We start making superclusters
1174  // with the hitmap.
1175  // cout<<" number of entries in module number "<<modnum<<" is "<<Ncell_mod<<endl;
1176  for (Int_t i = 0; i < Ncell_mod; i++){
1177 
1178  StPmdHit * mPmdHit = (StPmdHit*)PmdHits->At(i);
1179  // cout<<" mPmdHit deatils mod "<<mPmdHit->Gsuper()<<" detector="<<mPmdHit->SubDetector()<<endl;
1180 
1181  // Add those cells which have adc above cutoff to the hitmap.
1182  // The hitmap hold the index of the hit in the PmdHits TObjArray
1183  if((mPmdHit->Row()>NRowMax||mPmdHit->Row()<=0)||
1184  (mPmdHit->Column()>NColumnMax||mPmdHit->Column()<=0)){
1185  cout<<"hit out of range :"<<mPmdHit->Row()-1<<"/"<<mPmdHit->Column()-1<<endl;
1186  }else{
1187  if (mPmdHit->Edep() > cutoff ){
1188  hitmap[mPmdHit->Row()-1][mPmdHit->Column()-1]= i;
1189  NWithin++;
1190  }else {
1191  // cout<<"Adc less than cutoff Adc/cutoff "<<
1192  //mPmdHit->Edep()<<"/"<<cutoff<<endl;
1193  hitmap[mPmdHit->Row()-1][mPmdHit->Column()-1] = -1;
1194  Flag[mPmdHit->Row()-1][mPmdHit->Column()-1]=2;
1195  NAccCut++;
1196  }
1197  }
1198  }
1199  // hitmap ready here
1200 
1201  // Account statement
1202  if ((NWithin+NAccCut)!=Ncell_mod){
1203  cout<<" Some hits outside hitmap NWithin/NAccCut/NTotal "<<
1204  NWithin<<"/"<<NAccCut<<"/"<<Ncell_mod<<endl;
1205  }
1206  if (NWithin==0){
1207  return 0;
1208  }
1209  // cout<<" NWithin/NAccCut/NTotal "<<NWithin<<"/"<<NAccCut<<"/"<<Ncell_mod<<endl;
1210 
1211 
1212  //Make superclusters now
1213  //Initializing the number of Superclusters
1214  Int_t nPmdSuperCluster = -1;
1215 
1216  for(Int_t icell = 0; icell < Ncell_mod; icell++){
1217  StPmdHit *mPmdHit = (StPmdHit*)PmdHits->At(icell);
1218  // cout<<" flag of hit["<<icell<<"] is Flag["<<mPmdHit->Row()-1<<"]["<<mPmdHit->Column()-1<<"]="<<Flag[mPmdHit->Row()-1][mPmdHit->Column()-1]<<endl;
1219  if((mPmdHit->Row()>NRowMax||mPmdHit->Row()<=0)||
1220  (mPmdHit->Column()>NColumnMax||mPmdHit->Column()<=0)){
1221  cout<<"hit out of range :"<<mPmdHit->Row()-1<<"/"<<mPmdHit->Column()-1<<endl;
1222  continue;
1223  }
1224 
1225  if (Flag[mPmdHit->Row()-1][mPmdHit->Column()-1]==0){
1226  nPmdSuperCluster++;
1227  // cout<<"npmdsuperclusters="<<nPmdSuperCluster<<endl;
1228  TList * mPmdSuperCluster = (TList*)mPmdSuperClusters->New(nPmdSuperCluster);
1229  // cout<<" made a new supercluster with address "<<mPmdSuperCluster<<" in list with address "<<mPmdSuperClusters<<endl;
1230  TIter ClusterListIter(mPmdSuperCluster);
1231  //Add hit to the supercluster and put flag =1
1232  mPmdSuperCluster->Add(mPmdHit);
1233  // cout<<" added hit # "<<icell<<" to supercluster # "<<nPmdSuperCluster<<endl;
1234  Flag[mPmdHit->Row()-1][mPmdHit->Column()-1]=1;
1235  //make a list of its neighbours and their neighbours
1236 
1237  for(Int_t k=0;k<mPmdSuperCluster->GetSize();k++){
1238  StPmdHit* clusterhit = (StPmdHit*)mPmdSuperCluster->At(k);
1239  Int_t clusterhitRow = clusterhit->Row();
1240  Int_t clusterhitColumn = clusterhit->Column();
1241  if(Flag[clusterhitRow-1][clusterhitColumn-1]!=1) cout<<"???????"<<endl;
1242  //checking neighbours
1243  for( Int_t ineigh = 0;ineigh < 6; ineigh++ ){
1244  Int_t neighbourRow = clusterhitRow + Rowneigh[ineigh];
1245  Int_t neighbourColumn = clusterhitColumn + Columnneigh[ineigh];
1246  Int_t boundary = CheckBoundary(modnum,neighbourRow-1,neighbourColumn-1);
1247 
1248  if (boundary==0 && hitmap[neighbourRow-1][neighbourColumn-1]!=-1){
1249 
1250  // if (hitmap[neighbourRow-1][neighbourColumn-1]!=-1){
1251  // if(mOptDebug)cout<<"neigh row/col "<<neighbourRow<<"/"<<neighbourColumn<<endl;
1252  StPmdHit* neighbour = (StPmdHit*)PmdHits->At(hitmap[neighbourRow-1][neighbourColumn-1]);
1253  // cout<<" ineigh = "<<ineigh<<" Got hit number "<<hitmap[neighbourRow-1][neighbourColumn-1]<<" with adc"<<neighbour->Edep()<<" as connecting hit. k= "<<k<<" SpCl size is "<<mPmdSuperCluster->GetSize()<<endl;
1254  if (neighbour){
1255  if (neighbour->Edep()>cutoff
1256  && Flag[neighbourRow-1][neighbourColumn-1]==0){
1257  mPmdSuperCluster->Add(neighbour);
1258  // cout<<"The size of supercluster is "<<mPmdSuperCluster->GetSize()<<endl;
1259  // cout<<"added "<<hitmap[neighbourRow-1][neighbourColumn-1]<<" with adc"<<neighbour->Edep()<<" as connecting hit"<<endl;
1260  Flag[neighbourRow-1][neighbourColumn-1]=1;
1261  }//if neighbour is above cutoff
1262  }//if neighbour exists and is not already part of the supercluster
1263  }//if hitmap exists
1264  } //looping over 6 neighbours
1265  }//looping over superclusterlist (for neighbours)
1266  // if (mOptDebug && (mPmdSuperCluster->GetSize()>MaxSuperSize)){
1267  if (mPmdSuperCluster->GetSize()>MaxSuperSize){
1268  cout<<"The size of supercluster "<<mPmdSuperCluster->GetSize()<<" exceeds MaxSuperSize("<<MaxSuperSize<<")"<<endl;
1269  return 0;
1270  }
1271  mPmdSuperCluster->Sort(0);
1272 
1273  }//if flag=0
1274 
1275  } // looping over all hits TObjArray
1276 
1277  // cout<<"going out of MakerSuperClusters with ";
1278  Int_t Nentries = mPmdSuperClusters->GetEntries();
1279  // cout<<"Got "<<Nentries<<" Superclusters."<<endl;
1280 
1281  return Nentries;
1282 }
1283 
1284 //--------------------------------------------------------------------------------------
1285 
1286 //Int_t StPmdClustering::MakeClusters(TClonesArray * mPmdSuperClusters, StPmdClusterCollection * pmdclus)
1287 Int_t StPmdClustering::MakeClusters(TClonesArray * mPmdSuperClusters)
1288 {
1289 
1290  // cout<<" In MakeCluster"<<endl;
1291 
1292  Int_t nLocalPeak=0;
1293  Int_t sum_nclusters = 0;
1294  // cout<<"MakeCluster::# of superclusters is "<<mPmdSuperClusters->GetEntries()<<endl;
1295  // cout<<" address of mPmdSuperCluster list is "<<mPmdSuperClusters<<endl;
1296  for(Int_t i=0;i<mPmdSuperClusters->GetEntries();i++){
1297  // cout<<" Got supercluster number "<<i<<" / "<<mPmdSuperClusters->GetEntries()<<endl;
1298  TList *mPmdSuperCluster = (TList*)mPmdSuperClusters->At(i);
1299  Int_t ncell = mPmdSuperCluster->GetSize();
1300  // cout<<"# of cells in "<<i<<" superclusters is "<<ncell<<" address is "<<mPmdSuperCluster<<endl;
1301  //Initialzation of cluster parameters
1302  // Int_t nmod,det;
1303  if (ncell==0) {
1304  break;
1305  cout<<"zero cells in this supercluster!!!"<<endl;
1306  }else{
1307 
1308  if ((ncell<3 && mOptRefineCluster==kTRUE)||(mOptRefineCluster==kFALSE)){
1309 
1310  // cout<<" converting supercluster to cluster with ncell = "<<ncell<<endl;
1311 
1312  StPmdCluster *pcluster = new StPmdCluster();;
1313  pmdclus->addCluster(pcluster);
1314  Float_t hitfract[ncell];
1315 
1316  for(Int_t j = 0; j < ncell; j++ ){
1317  StPmdHit * phit = (StPmdHit*)mPmdSuperCluster->At(j);
1318  pcluster->addHitCollection(phit);
1319  hitfract[j]=1.0;
1320  }
1321  // Float_t nmem =
1322  BuildCluster(pcluster,hitfract);
1323  // cout<<" No refine and so I have a cluster with "<<nmem<<" cells"<<endl;
1324 
1325  sum_nclusters++;
1326  }
1327 
1328  // This option is used only if ncell >=3 and Refine Clustering is ON
1329  else{
1330  // cout<<" Got a big cluster and Refined clustering ON says I have to break it"<<endl;
1331  // jLocalPeak is array of the hits in a supercluster which are
1332  // Local peaks. Since the index of hits in a supercluster starts
1333  // with 0, the jLocalPeaks are initialised to -1 and not 0.
1334  Int_t jLocalPeak[MaxLocalPeaks];
1335  for (Int_t j=0;j<MaxLocalPeaks;j++){
1336  jLocalPeak[j] = -1; //
1337  }
1338  // Get local peaks for SuperCluster with ncell
1339  Int_t *pLocalPeak = &jLocalPeak[0];
1340  // cout<<"Choice="<<LocalPeakChoice<<BreakSuperClusterChoice<<SigmaChoice<<endl;
1341 
1342  nLocalPeak = GetLocalPeaks(pLocalPeak,mPmdSuperCluster);
1343 
1344  // if (mOptDebug)
1345  // cout<<i<<"\t"<<" pLocalPeak = "<< pLocalPeak <<"\t"<<"nLocalPeak = "<<nLocalPeak<<endl;
1346 
1347  // If there is only 1 local peak then we dont have to break the supercluster
1348  if(nLocalPeak==1){
1349 
1350  StPmdCluster *pcluster = new StPmdCluster();;
1351  pmdclus->addCluster(pcluster);
1352  Float_t hitfract[nLocalPeak];
1353 
1354  // put all the hits into this cluste collection and send it to builtcluster
1355  for(Int_t ihit = 0;ihit<mPmdSuperCluster->GetEntries();ihit++){
1356  hitfract[ihit]=1.0;
1357  StPmdHit * phit = (StPmdHit*)mPmdSuperCluster->At(ihit);
1358  pcluster->addHitCollection(phit);
1359  }
1360  Float_t nmem = BuildCluster(pcluster,hitfract);
1361 
1362  if(nmem==0) cout<<" 1 local peak only so I have a cluster with "<<nmem<<" cells"<<endl;
1363  }
1364 
1365  if (nLocalPeak>1 && nLocalPeak<MaxLocalPeaks){
1366  sum_nclusters+=nLocalPeak;
1367  /*
1368  cout<<" sum_clusters = "<<sum_nclusters<<endl;
1369  cout<<" number of hits in this SuperCluster ="<<mPmdSuperCluster->GetEntries()<<endl;
1370  cout<<" Number of clusters in ClusterCollection is "<<pmdclus->Nclusters()<<endl;
1371  cout<<" address of jLocalPeak = "<<pLocalPeak<<endl;
1372  cout<<" address of mPmdSuperCluster = "<<mPmdSuperCluster<<endl;
1373  cout<<" address of pmdclusCollection = "<<pmdclus<<endl;
1374  */
1375  // Bool_t broken = BreakSuperCluster(nLocalPeak,pLocalPeak,mPmdSuperCluster,pmdclus);
1376  Bool_t broken = BreakSuperCluster(nLocalPeak,pLocalPeak,mPmdSuperCluster);
1377  // if(mOptDebug)
1378 
1379  if(broken==1){
1380  cout<<" too many local peaks I think"<<endl;
1381  }
1382  }else{
1383  // the module was cut due to number of localpeaks>MaxLocalPeaks
1384  //if (mOptDebug)
1385  cout<<"Module cut due to nLocalPeaks>MaxLocalPeaks"<<endl;
1386  return 0;
1387  }
1388  }
1389  }
1390  }
1391  // cout<<"out of makeclusters"<<endl;
1392  return sum_nclusters;
1393 }
1394 
1395 
1396 //--------------------------------------------------------------------------------------
1397 //------------------------------------------------------------------
1398 //Returns cell xy in units of the cell size wrt to the cornermost cell.
1399 // The corner most cell has row=0;col=0;
1400 //CAUTION: Not absolute xy.
1401 
1402 void StPmdClustering::CellXY(Int_t row, Int_t column, Float_t& xx, Float_t& yy){
1403  // column = column -1;
1404  // row = row -1;
1405  xx = (Float_t)(column + 0.5*row);
1406  yy = 0.5*sqrt(3.)*row;
1407 }
1408 //---------------------------------------------------------------------
1409 
1410 Int_t StPmdClustering::CheckBoundary(Int_t nmod,Int_t numRow, Int_t numColumn){
1411  if ((numRow>=0 && numRow< NRowMax) &&
1412  (numColumn>=0 && numColumn<NColumnMax))
1413  { // cout<<"In checkboundary"<<endl;
1414  return 0;}
1415  else return 1;
1416 }
1417 //---------------------------------------------------------------------
1418 
1419 Int_t StPmdClustering::GetLocalPeaks(Int_t* jLocalPeak,TList *mPmdSuperCluster){
1420 
1421  Int_t nLocalPeak =-1;
1422  Int_t ncell = mPmdSuperCluster->GetSize();
1423  Float_t xLocalPeak[ncell];
1424  Float_t yLocalPeak[ncell];
1425  Float_t adcLocalPeak[ncell];
1426 
1427  //The highest peak of the supercluster is definately a peak so:
1428  nLocalPeak++;
1429  if(nLocalPeak>=MaxLocalPeaks) cout<<" Number of Local peaks is "<<nLocalPeak+1<<" which exceeds the array size MaxLocalPeaks ="<<MaxLocalPeaks<<endl;
1430  StPmdHit * phitj = (StPmdHit*)mPmdSuperCluster->At(0);
1431  Float_t hitx1 = 0; Float_t hity1 = 0;
1432  CellXY(phitj->Row()-1,phitj->Column()-1,hitx1,hity1);
1433  Int_t sm = phitj->Gsuper();
1434  jLocalPeak[nLocalPeak] = 0;
1435  xLocalPeak[nLocalPeak] = hitx1;
1436  yLocalPeak[nLocalPeak] = hity1;
1437  adcLocalPeak[nLocalPeak] = phitj->Edep();
1438  // cout<<"Out of ncell ="<<ncell<<" cell 0 has adc ="<<phitj->adc()<<endl;
1439  //cout<<"nLocalPeak = "<<nLocalPeak<<
1440  // " jLocalPeak ="<<jLocalPeak[nLocalPeak]<<
1441  // " adcLocalPeak = "<<adcLocalPeak[nLocalPeak]<<
1442  // " x,y = "<<xLocalPeak[nLocalPeak]<<","<<yLocalPeak[nLocalPeak]<<endl;
1443  // Check the next highest cell........
1444  for (Int_t j = 1; j < ncell; j++){
1445  StPmdHit * phitj = (StPmdHit*)mPmdSuperCluster->At(j);
1446  // cout<<"cell "<<j<<" has adc ="<<phitj->adc()<<endl;
1447  Float_t hitx = 0; Float_t hity = 0;
1448  CellXY(phitj->Row()-1,phitj->Column()-1,hitx,hity);
1449  //......against all the previous LocalPeaks
1450  Int_t nclear = -1;
1451  for (Int_t k = 0; k <= nLocalPeak; k++){
1452  Float_t gap = Dist(hitx,hity,xLocalPeak[k],yLocalPeak[k]);
1453  // Float_t fadc = 1.0*phitj->adc()/adcLocalPeak[k];
1454  //the distances are in units of cell sizes.
1455  if( gap >= 1.1 && gap < 1.8 && phitj->Edep() > adcLocalPeak[k]*0.30) nclear++;
1456  if( gap >= 1.8 && gap < 2.1 && phitj->Edep() > adcLocalPeak[k]*0.15) nclear++;
1457  if( gap >= 2.1 && gap < 2.8 && phitj->Edep() > adcLocalPeak[k]*0.05) nclear++;
1458  if( gap >= 2.8 ) nclear++;
1459  }
1460  if (nclear==nLocalPeak){//cleared against all previous local peaks
1461  // then it is another local peak so add to list
1462  nLocalPeak++;
1463  if (nLocalPeak<MaxLocalPeaks){
1464  jLocalPeak[nLocalPeak] = j;
1465  xLocalPeak[nLocalPeak] = hitx;
1466  yLocalPeak[nLocalPeak] = hity;
1467  adcLocalPeak[nLocalPeak] = phitj->Edep();
1468  }else{
1469  cout<<"number of local peaks is "<<nLocalPeak<<" which is more than "<<MaxLocalPeaks<<" for a supercluster of size ="<<ncell<<" in sm number ="<<sm<<endl;
1470  return 0;
1471  }
1472  }
1473  }
1474  return ++nLocalPeak;
1475 
1476 }
1477 //----------------------------------------------------------------------
1478 void StPmdClustering::Cell_eta_phi(Float_t xreal,Float_t yreal,Float_t& eta,Float_t& phi){
1479 
1480  // Converts xyz from STAR 0,0,0 to actual vertex
1481  xreal = xreal - mVertexPos.x();
1482  yreal = yreal - mVertexPos.y();
1483  Float_t zreal = -1.0*geom->GetPmdZ();
1484  zreal = zreal - mVertexPos.z();
1485  // cout<<"IN Clustering: xreal="<<xreal<<" yreal="<<yreal<<" zreal="<<zreal;
1486 
1487  Float_t rdist = (TMath::Sqrt(xreal*xreal + yreal*yreal))/zreal;
1488 
1489  Int_t flag = -1;
1490  if (rdist<0.) {
1491  flag = +1;
1492  rdist=TMath::Abs(rdist);
1493  }
1494 
1495  Float_t theta = TMath::ATan(rdist);
1496  if(theta !=0.)eta = flag*TMath::Log(TMath::Tan(0.5*theta));
1497  if( xreal==0) {
1498  if(yreal>0) phi = TMath::Pi()/2.0;
1499  if(yreal<0) phi = -1.0*TMath::Pi()/2.0;
1500  }
1501  if(xreal != 0) phi = atan2(yreal,xreal);
1502 
1503  // cout<<" theta="<<theta<<" eta="<<eta<<" phi="<<phi<<endl;
1504 }
1505 
1506 
1507 //-----------------------------------------------------------------------------
1508 Float_t StPmdClustering::BuildCluster( StPmdCluster *pcluster,Float_t* hitfract){
1509 
1510 
1511  // Get the hit collection from the pcluster and built the cluster
1512  TObjArray* hitCollection = pcluster->HitCollection();
1513  Int_t ncell = hitCollection->GetEntries();
1514 
1515  Float_t adchit[ncell],xhit[ncell],yhit[ncell];
1516  Float_t SigmaL = 0, SigmaS = 0;
1517  Float_t cluX=0, cluY=0, cluAdc=0;
1518  Float_t eta=0,phi=0;
1519 
1520 
1521  // cout<<" Got hit collection with size = "<<ncell<<endl;
1522  // Getting the module number of the cluster from the module and detector number of hit.
1523  // and setting the module and number of member hits of pcluster
1524  StPmdHit* phit = (StPmdHit*)hitCollection->At(0);
1525  Int_t nmod = phit->Gsuper();
1526  Int_t det = phit->SubDetector();//1 for pmd and 2 for cpv
1527  // cout<<" nmod = "<<nmod<<" det = "<<det<<endl;
1528  // nmod = 12*(2-det) + nmod-1;//to go to 0-23 format
1529  Int_t nmodclu = 12*(2-det) + nmod;//to go to 0-23 format
1530  // cout<<" New module number for clusters is "<<nmodclu<<endl;
1531  pcluster->setModule(nmodclu);
1532  Float_t nmem = Float_t(ncell*1.0);
1533  pcluster->setNumofMems(nmem);
1534 
1535 
1536  // Now to calculate the eta,phi and sigmaS,L and ADC of the cluster
1537 
1538  // For a single cell cluster
1539  if (ncell==1){
1540  SigmaL=0;
1541  SigmaS=0;
1542  StPmdHit * phit = (StPmdHit*)hitCollection->At(0);
1543 
1544  // cout<<" Building cluster for 1 cell with adc ="<<phit->Edep()<<endl;
1545  // Eta phi is useless here since it is wrt 0,0,0 and not wrt vertex
1546  geom->IntDetCell_xy(nmod,(phit->Row()),(phit->Column()),cluX,cluY,eta,phi);
1547  cluAdc = phit->Edep()*1.0;
1548  // cout<<" Cluster position and adc is is "<<cluX<<" "<<cluY<<" "<<cluAdc<<endl;
1549 
1550  }else{
1551  // If there are more than 1 hit then
1552  // loop over all the hits
1553  Float_t sumx=0,sumy=0;
1554  for(Int_t j = 0; j < ncell; j++ ){
1555  StPmdHit * phit = (StPmdHit*)hitCollection->At(j);
1556  if (!phit) cout<<" Did not get the phit I was hoping for "<<endl;
1557  if (!phit) break;
1558  // cout<<" Hi I am here with phit at address "<<phit<<endl;
1559 
1560 
1561  Float_t xreal=0,yreal=0;
1562  geom->IntDetCell_xy(nmod,(phit->Row()),(phit->Column()),xreal,yreal,eta,phi);
1563  adchit[j] = hitfract[j]*phit->Edep();
1564  xhit[j] = xreal;
1565  yhit[j] = yreal;
1566  // cout<<"Hit in MakeCluster:"<<nmodclu<<","<<det<<","<<phit->Row()<<","<<phit->Column()<<","<<phit->Edep()<<","<<xreal<<","<<yreal<<","<<eta<<","<<phi<<endl;
1567  // cout<<"Hit in BuildCluster:"<<nmodclu<<","<<det<<","<<phit->Row()<<","<<phit->Column()<<","<<hitfract[j]<<","<<phit->Edep()<<","<<xreal<<","<<yreal<<endl;
1568  sumx+=adchit[j]*xhit[j];
1569  sumy+=adchit[j]*yhit[j];
1570  cluAdc+=adchit[j];
1571 
1572  }
1573 
1574  cluX = sumx/cluAdc;
1575  cluY = sumy/cluAdc;
1576  // if(ncell>3) cout<<"Cluster: Xcen,Ycen,cluADC="<<cluX<<","<<cluY<<","<<cluAdc<<endl;
1577 
1578  Float_t sumxx=0,sumxy=0,sumyy=0;
1579  for(Int_t j=0;j<ncell;j++){
1580  sumxx+=adchit[j]*(xhit[j]-cluX)*(xhit[j]-cluX)/cluAdc;
1581  sumyy+=adchit[j]*(yhit[j]-cluY)*(yhit[j]-cluY)/cluAdc;
1582  sumxy+=adchit[j]*(yhit[j]-cluY)*(xhit[j]-cluX)/cluAdc;
1583  }
1584 
1585  Float_t b=sumxx+sumyy;
1586  Float_t c=sumxx*sumyy-sumxy*sumxy;
1587  // cout<<"sumxx,sumyy,sumxy,b,c="<<sumxx<<","<<sumyy<<","<<sumxy<<","<<b<<","<<c<<endl;
1588  Float_t r1=b/2.+sqrt(b*b/4.-c);
1589  Float_t r2=b/2.-sqrt(b*b/4.-c);
1590  // cout<<"r1,r2="<<r1<<","<<r2<<endl;
1591  if (r1<r2) {
1592  SigmaS = r1;
1593  SigmaL = r2;
1594  }else{
1595  SigmaS = r2;
1596  SigmaL = r1;
1597  }
1598  if(ncell==2) SigmaS=0;
1599  }
1600 
1601  // Get etaphi wrt vertex
1602  Cell_eta_phi(cluX,cluY,eta,phi);
1603 
1604  // I have all the cluster paramters here and now I will fill the cluster
1605 
1606  pcluster->setCluSigmaS(SigmaS);
1607  pcluster->setCluSigmaL(SigmaL);
1608 
1609  pcluster->setCluEta(eta);
1610  pcluster->setCluPhi(phi);
1611 
1612  Float_t edep = cluAdc*1.0;
1613  pcluster->setCluEdep(edep);
1614 
1615  pcluster->setCluX(cluX);
1616  pcluster->setCluY(cluY);
1617 
1618  Int_t pid = 0;
1619  Int_t edeppid = 0;
1620  Int_t mcpid = 0;
1621 
1622 
1623  pcluster->setMcCluPID(mcpid);
1624  pcluster->setCluPID(pid);
1625  pcluster->setCluEdepPID(edeppid);
1626  // cout<<"Cluster in BuildCluster :"<<nmod<<","<<ncell<<","<<cluX<<","<<cluY<<","<<eta<<","<<phi<<","<<edep<<endl;
1627 
1628  // cout<<" And have added it to the pmdclus (StPmdClusterCollection) which now has "<<pmdclus->Nclusters()<<endl;
1629 
1630  return nmem;
1631 
1632 }
1633 //-----------------------------------------------------------------------------------------
StPmdDetector::module_hit
Int_t module_hit(Int_t)
module number
Definition: StPmdDetector.cxx:83
StPmdClustering::findPmdClusters2
Bool_t findPmdClusters2(StPmdDetector *)
for Pmd clusters
Definition: StPmdClustering.cxx:1088
StPmdClusterCollection::Nclusters
Int_t Nclusters() const
destructor
Definition: StPmdClusterCollection.h:57
StPmdClustering::printclust
void printclust(Int_t, Int_t, StPmdCluster *)
distance between two clusters
Definition: StPmdClustering.cxx:281
StPmdDetector::setCluster
void setCluster(StPmdClusterCollection *)
number of clusters
Definition: StPmdDetector.cxx:128
StPmdClustering::GetHit
StPmdHit * GetHit(StPmdDetector *, Int_t, Double_t, Double_t)
for getting hits of each cluster
Definition: StPmdClustering.cxx:1049
StPmdHit::Row
Int_t Row() const
function for subdetector
Definition: StPmdHit.h:95
StPmdHit::Column
Int_t Column() const
function for row
Definition: StPmdHit.h:96
StPmdClustering::findPmdClusters
void findPmdClusters(StPmdDetector *)
for Pmd clusters
Definition: StPmdClustering.cxx:169
StPmdClustering::order
void order(Int_t)
for ordering
Definition: StPmdClustering.cxx:336
StPmdCluster::setModule
void setModule(Int_t)
hit collection
Definition: StPmdCluster.h:128
StPmdClustering::CentroidCal
Int_t CentroidCal(Int_t, Int_t, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &)
for Calculating cluster properties, those clusters having more then two cells
Definition: StPmdClustering.cxx:717
StPmdClustering::arrange
void arrange(Int_t)
ordering
Definition: StPmdClustering.cxx:371
StPmdClustering::refclust
void refclust(StPmdDetector *, Int_t, Int_t, Int_t, StPmdClusterCollection *)
refined clustering
Definition: StPmdClustering.cxx:405
StPmdClustering::~StPmdClustering
virtual ~StPmdClustering()
destructor
Definition: StPmdClustering.cxx:164
kStWarn
Definition: Stypes.h:42
StPmdHit::Gsuper
Int_t Gsuper() const
A destructor.
Definition: StPmdHit.h:91
kStOK
Definition: Stypes.h:40
StPmdClustering::crclust
Int_t crclust(Double_t, Double_t, Int_t, Int_t)
crude clustering
Definition: StPmdClustering.cxx:974
StPmdHit::Edep
Float_t Edep() const
function for col
Definition: StPmdHit.h:97
StPmdDetector::module
StPmdModule * module(unsigned int)
number of hits
Definition: StPmdDetector.cxx:104
StPmdHit::SubDetector
Int_t SubDetector() const
function for module
Definition: StPmdHit.h:94
StPmdClustering::Dist
Double_t Dist(Double_t, Double_t, Double_t, Double_t)
Definition: StPmdClustering.cxx:966
Generated by   doxygen 1.8.5