StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StEEmcIUPointMaker.cxx
1 
26 #include "StEEmcIUPointMaker.h"
27 #include "StEEmcPool/StEEmcA2EMaker/StEEmcA2EMaker.h"
28 #include "StEEmcIUClusterMaker.h"
29 #include "StEEmcIUCluster.h"
30 #include "StEEmcIUSmdCluster.h"
31 #include "StEEmcUtil/EEmcGeom/EEmcGeomDefs.h"
32 #include "StEEmcUtil/EEmcGeom/EEmcGeomSimple.h"
33 #include "StEEmcUtil/StEEmcSmd/EEmcSmdGeom.h"
34 #include "TH1F.h"
35 #include "TH2F.h"
36 #include <iostream>
37 #include <algorithm>
38 #include <map>
39 
40 #include "StEEmcPool/StEEmcPointMaker/eeTowerFunction.h"
41 
42 /* StEvent stuff */
43 #include "StEvent/StEvent.h"
44 #include "StEvent/StEmcCollection.h"
45 #include "StEvent/StEmcPoint.h"
46 
47 /* Root's linear algebra package */
48 #include "TMatrixF.h"
49 
50 //#define DEBUG_buildPoints
51 
52 ClassImp(StEEmcIUPointMaker);
53 
54 // ----------------------------------------------------------------------------
55 StEEmcIUPointMaker::StEEmcIUPointMaker(const Char_t *name):StMaker(name)
56 {
57  std::cout << "StEEmcIUPointMaker("<<name<<")" << std::endl;
58 
62  mEEtow=new EEmcGeomSimple();
63  mEEsmd=EEmcSmdGeom::instance();
64  mEEmap=EEmcSmdMap::instance();
65 
66  mTowerThreshold=0.;
67  mFillStEvent=false;
68  mEnergyMode=1;
69  mLimit=10;
70  mSmdMatch = 0.;
71  ukey=0;
72  vkey=0;
73 
74 }
75 
76 // ----------------------------------------------------------------------------
77 Int_t StEEmcIUPointMaker::Init()
78 {
79  mEEanalysis=(StEEmcA2EMaker*)GetMaker(mNameAnalysis);
80  mEEclusters=(StEEmcIUClusterMaker*)GetMaker(mNameClusters);
81  assert(mEEanalysis);
82  assert(mEEclusters);
83  hZEratio=new TH1F("smdezratio","energy Zratio between smd clusters",10,0.,1.);
84  pointet=new TH1F("pointet","transverse energy of every point candidate",40,0,10);
85  pointgeo=new TH2F("pointgeo", "Reconstructed Point Geometry; phi / deg; Eta",360,-180.,180,20,1.,2.);
86  return StMaker::Init();
87 }
88 
89 // ----------------------------------------------------------------------------
90 Int_t StEEmcIUPointMaker::Make()
91 {
92 
99 
100  // Loop over all 12 EEMC sectors
101  for ( Int_t sector=0; sector<12; sector++ )
102  {
103 
105  StEEmcIUSmdClusterVec_t uclusters=mEEclusters->smdclusters(sector,0);
106  StEEmcIUSmdClusterVec_t vclusters=mEEclusters->smdclusters(sector,1);
107 
109  std::sort( uclusters.begin(), uclusters.end(), inner );
110  std::sort( vclusters.begin(), vclusters.end(), inner );
111 
112  findPoints( sector, uclusters, vclusters, mPoints );
113 
114  }
115 
116 
123  for ( UInt_t i=0;i<mPoints.size();i++ )
124  {
125 
126  mPoints[i].energy( (Float_t)(mPoints[i].energy()/0.007/2.) );
127 
128  TVector3 pointpos=mPoints[i].position();
129  pointpos=pointpos.Unit();
130  Float_t pet = (mPoints[i].energy()*pointpos).Perp();
131  pointet->Fill(pet);
132  pointgeo->Fill(mPoints[i].position().Phi()/3.1416*180,mPoints[i].position().Eta());
133  }
134 
135 
136 
137 
138  StEEmcIUPointVec_t orgpoints=mPoints;
139 
141  // shareEnergy(); // <<<<<<< leads to negative point energies... rethink
142 
143  if ( mEnergyMode == 1 )
144  shareEnergySimple();
145  else
146  shareEnergySmd();
147 
148 
150  countRelatives();
151 
152 
153  if ( mFillStEvent )
154  {
155  fillStEvent();
156  verifyStEvent();
157  }
158 
159  /*
160  for ( UInt_t i=0;i<mPoints.size();i++ )
161  {
162  mPoints[i].print();
163  }
164  for ( Int_t i=0;i<10;i++ ) std::cout << std::endl;
165 */
166 
167  return kStOK;
168 }
169 
170 // ----------------------------------------------------------------------------
171 
172 StEEmcIUPointVec_t StEEmcIUPointMaker::buildSmdPoints( Int_t sector,
173  StEEmcIUSmdClusterVec_t &u,
174  StEEmcIUSmdClusterVec_t &v )
175 {
176 
177  StEEmcIUPointVec_t points;
178  //printf("u.size=%d v.size=%d\n",u.size(),v.size());
179  for ( UInt_t i=0; i<u.size(); i++ )
180  {
181 
182  StEEmcIUSmdCluster uc=u[i];
183 
184  Float_t xu=uc.mean();
185 
186 
187  for ( UInt_t j=0;j<v.size(); j++ )
188  {
189 
190 
191  StEEmcIUSmdCluster vc=v[j];
192  Float_t xv=vc.mean();
193 
195  TVector3 direct = mEEsmd->getIntersection( sector, xu, xv );
196 
197  Int_t sec,sub,eta;
198  if ( !mEEtow->getTower(direct,sec,sub,eta) )
199  {
200  continue;
201  }
202  else
203  {
205  }
206 
209  if ( sector != sec )
210  continue;
211 
212 
213 
217  Bool_t good = false;
218  if ( mEEanalysis->tower(sec,sub,eta).energy() > mTowerThreshold )
219  good=true;
220  else if ( mEEanalysis->tower(sec,sub,eta).fail() )
221  {
222  for ( Int_t layer=1;layer<=3;layer++ )
223  {
224  if ( mEEanalysis->tower(sec,sub,eta,layer).energy() > 0. )
225  good=true;
226  }
227  }
228 
231  /*
232  Float_t uvdiff = TMath::Abs( uc.energy() - vc.energy() );
233  Float_t uvavg = 0.5 * ( uc.energy() + vc.energy() );
234  if ( uvdiff > mSmdMatch * uvavg ) good = false;
235  */
236  Float_t Eu=uc.energy();
237  Float_t Ev=vc.energy();
238  if ( Eu < mSmdMatch * Ev || Ev < mSmdMatch * Eu ) good=false;
239 
240  if ( good ) {
241 
242  StEEmcIUPoint p;
243  p.cluster( uc, 0 );
244  p.cluster( vc, 1 );
245  p.energy( uc.energy() + vc.energy() );
246  p.tower( mEEanalysis->tower(sec,sub,eta) );
247  TVector3 position=mEEsmd->getIntersection(sector,uc.mean(),vc.mean());
248  p.position(position);
249  points.push_back(p);
250 
252  mSmdPoints.push_back(p);
253 
254 
255  }
256 
257  }
258 
259  }
260 
261  return points;
262 
263 }
264 
265 
266 
267 
268 
269 
270 
271 
272 
273 
274 
275 
276 
277 
278 
279 
280 
281 // ----------------------------------------------------------------------------
282 Bool_t StEEmcIUPointMaker::findPoints( Int_t sector,
283  StEEmcIUSmdClusterVec_t uclusters,
284  StEEmcIUSmdClusterVec_t vclusters,
285  StEEmcIUPointVec_t &points )
286 {
287 
288 
290  StEEmcIUPointVec_t Tmypoints;
291  StEEmcIUPointVec_t mypoints;
293  std::sort( uclusters.begin(), uclusters.end(), inner );
294  std::sort( vclusters.begin(), vclusters.end(), inner );
295 
297  StEEmcIUPointVec_t smdpoints = buildSmdPoints( sector, uclusters, vclusters );
298 
299 
301  if ( smdpoints.size() < 1 ) return false;
302 
304  std::sort( smdpoints.begin(), smdpoints.end(), chiSquare );
305 
306 
310  std::map< Int_t, std::vector<Int_t> > u2p, v2p;
311  for ( UInt_t i=0; i<smdpoints.size(); i++ )
312  {
313  u2p[ smdpoints[i].cluster(0).key() ].push_back( i );
314  v2p[ smdpoints[i].cluster(1).key() ].push_back( i );
315 
316 
317  }
318 
320  // Bool_t go = false;
321  StEEmcIUSmdClusterVec_t::iterator uiter=uclusters.begin();
322  StEEmcIUSmdClusterVec_t::iterator viter=vclusters.begin();
323 
324 
329 
330  // ----------<<<<<<<<< stage one >>>>>>>>>-------------
331 
332  int iii=0,jjj=0;
333 
334  while ( uiter<uclusters.end() || viter<vclusters.end() )
335  {
336 
338  StEEmcIUSmdCluster ucl;ucl.key(-1);
339  StEEmcIUSmdCluster vcl;vcl.key(-1);
340  if ( uiter<uclusters.end() ) ucl=(*uiter);
341  if ( viter<vclusters.end() ) vcl=(*viter);
342  Int_t iUV=-1;
343  if ( ucl.key()<0 )
344  iUV=1;
345  else if ( vcl.key()<0 )
346  iUV=0;
347  else if ( (*uiter).mean() < (*viter).mean() )
348  iUV=0;
349  else
350  iUV=1;
351 
353  StEEmcIUSmdCluster &cluster=(iUV==0)?ucl:vcl;
354 
356  std::vector<Int_t> matched=(iUV==0)?
357  u2p[cluster.key()]:
358  v2p[cluster.key()];
359 
360 
363 
364  if ( matched.size()==0 || matched.size() >1 )
365  {
366  if ( iUV==0 ) uiter++;
367  if ( iUV==1 ) viter++;
368  continue;
369  }
370 
376 
379  StEEmcIUPoint p=smdpoints[matched.back()];
380 
382  Tmypoints.push_back(p);
383 
384  if ( iUV==0 ) {uiter++;iii++;}
385  if ( iUV==1 ) {viter++;jjj++;}
386 
387  }//end of while
388 
389 
390  //splitting algo starts. When two points ramp into one lane, it only shows us one Smd cluster information in one of the Smd planes, thus can only reconstruct one point. We now split the overlapped U/V Smd cluster based on the other two Smd cluster in V/U accordingly.
391  std::sort( Tmypoints.begin(), Tmypoints.end(), chiSquare );
392  //>=2 to turn on splitting algo, <0 to off splitting algo.
393  if(Tmypoints.size()>=2)
394  {
395  int aa=0,bb=0;
396  for ( UInt_t i=0; i<Tmypoints.size(); i++ )
397  {
398  StEEmcIUSmdCluster uone=Tmypoints[i].cluster(0);
399 
400  StEEmcIUSmdCluster vone=Tmypoints[i].cluster(1);
401 
402 
403  for ( UInt_t j=i+1; j<Tmypoints.size(); j++ )
404  {
405  StEEmcIUSmdCluster utwo=Tmypoints[j].cluster(0);
406  StEEmcIUSmdCluster vtwo=Tmypoints[j].cluster(1);
407 
408  if((uone.mean()==utwo.mean())||(vone.mean()==vtwo.mean()))
409  {
410  //case 1
411  if(uone.mean()==utwo.mean()&&(aa==0)&&(vone.mean()!=vtwo.mean()))
412  {
413  float uzr=fabs(uone.energy()-vone.energy()-vtwo.energy())/(uone.energy()+vone.energy()+vtwo.energy());
414 
415  hZEratio->Fill(uzr);
416 
417  if(uzr<=0.2)
418  {
419  float vratio1=vone.energy()/(vone.energy()+vtwo.energy());
420  float vratio2=vtwo.energy()/(vone.energy()+vtwo.energy());
421 
422  int ukey=uone.key();
423  int ukey1=ukey+100;
424  int ukey2=ukey+200;
425  float u1energy=uone.energy()*vratio1;
426  float u2energy=uone.energy()*vratio2;
427  uone.key(ukey1);
428  utwo.key(ukey2);
429  uone.energy(u1energy);
430  utwo.energy(u2energy);
431  float p1energy=uone.energy()+vone.energy();
432  Tmypoints[i].energy(p1energy);
433  float p2energy=utwo.energy()+vtwo.energy();
434  Tmypoints[j].energy(p2energy);
435  Tmypoints[i].cluster(uone,0);
436  Tmypoints[j].cluster(utwo,0);
437  StEEmcIUPoint p1=Tmypoints[i];
438 
439  //mypoints.push_back(p1);
440  StEEmcIUPoint p2=Tmypoints[j];
441  //mypoints.push_back(p2);
442 
443  removeCluster( uclusters, p1.cluster(0).key() );
444  removeCluster( vclusters, p1.cluster(1).key() );
445  points.push_back(p1);
446  removeCluster( uclusters, p2.cluster(0).key() );
447  removeCluster( vclusters, p2.cluster(1).key() );
448  points.push_back(p2);
449  removeCluster( uclusters, ukey );
450 
451  findPoints(sector, uclusters, vclusters, points );
452  return true;
453 
454  aa++;
455  }
456  } //end of case 1
457  //case 2
458  if(vone.mean()==vtwo.mean()&&(bb==0)&&(uone.mean()!=utwo.mean()))
459  {
460  float vzr=fabs(vone.energy()-uone.energy()-utwo.energy())/(vone.energy()+uone.energy()+utwo.energy());
461  hZEratio->Fill(vzr);
462  if(vzr<=0.2)
463  {
464  float uratio1=uone.energy()/(uone.energy()+utwo.energy());
465  float uratio2=utwo.energy()/(uone.energy()+utwo.energy());
466 
467 
468  int vkey=vone.key();
469  int vkey1=vkey+100;
470  int vkey2=vkey+200;
471  float v1energy=vone.energy()*uratio1;
472  float v2energy=vone.energy()*uratio2;
473  vone.key(vkey1);
474  vtwo.key(vkey2);
475  vone.energy(v1energy);
476  vtwo.energy(v2energy);
477  float p1energy=vone.energy()+uone.energy();
478  Tmypoints[i].energy(p1energy);
479  float p2energy=vtwo.energy()+utwo.energy();
480  Tmypoints[j].energy(p2energy);
481  Tmypoints[i].cluster(vone,1);
482  Tmypoints[j].cluster(vtwo,1);
483  StEEmcIUPoint p1=Tmypoints[i];
484 
485  //mypoints.push_back(p1);
486  StEEmcIUPoint p2=Tmypoints[j];
487  //mypoints.push_back(p2);
488 
489  removeCluster( uclusters, p1.cluster(0).key() );
490  removeCluster( vclusters, p1.cluster(1).key() );
491  points.push_back(p1);
492  removeCluster( uclusters, p2.cluster(0).key() );
493  removeCluster( vclusters, p2.cluster(1).key() );
494  points.push_back(p2);
495  removeCluster( vclusters, vkey );
496 
497  findPoints(sector, uclusters, vclusters, points );
498  return true;
499 
500 
501  bb++;
502  }
503  }//end of case 2
504  //case 3
505  if(uone.mean()==utwo.mean()&&(vone.mean()==vtwo.mean()))
506  {
507  int cc=0;
508  for ( UInt_t i=0; i<Tmypoints.size(); i++ )
509  {
510  StEEmcIUPoint p=Tmypoints[i];
511  mypoints.push_back(p);
512  cc++;
513  }
514 
515  //printf("caase 3 points.size=%d\n",mypoints.size());
516  }//end of case 3
517  }
518  }
519  }
520 
521  //case 4 catch those isolated points
522  for ( UInt_t i=0; i<Tmypoints.size(); i++ )
523  {
524  StEEmcIUPoint pp=Tmypoints[i];
525  mypoints.push_back(pp);
526  }//end of case 4
527  }
528  else
529  {
530  for ( UInt_t i=0; i<Tmypoints.size(); i++ )
531  {
532  StEEmcIUPoint p=Tmypoints[i];
533  mypoints.push_back(p);
534  }
535  }
538  Float_t chisq=9.0E9;
539  Int_t imin=-1;
540  for ( UInt_t i=0; i<mypoints.size(); i++ )
541  {
542  Float_t eu=mypoints[i].cluster(0).energy();
543  Float_t ev=mypoints[i].cluster(1).energy();
544  Float_t x2=fabs((eu-ev)/(eu+ev));
545  if ( x2 < chisq ) {
546  imin=(Int_t)i;
547  chisq=x2;
548  }
549  }
550 
551 
558  if ( imin >= 0 ) {
559 
560  StEEmcIUPoint p=mypoints[imin];
561 
562  removeCluster( uclusters, p.cluster(0).key() );
563  removeCluster( vclusters, p.cluster(1).key() );
564  points.push_back(p);
565  findPoints(sector, uclusters, vclusters, points );
566  return true;
567 
568  }
569 
570 
571 
572 
573 
575 
576 
583  uiter=uclusters.begin();
584  viter=vclusters.begin();
585  while ( uiter!=uclusters.end() || viter!=vclusters.end() )
586  {
587 
589  StEEmcIUSmdCluster ucl;ucl.key(-1);
590  StEEmcIUSmdCluster vcl;vcl.key(-1);
591  if ( uiter!=uclusters.end() ) ucl=(*uiter);
592  if ( viter!=vclusters.end() ) vcl=(*viter);
593  Int_t iUV=-1;
594  if ( ucl.key()<0 )
595  iUV=1;
596  else if ( vcl.key()<0 )
597  iUV=0;
598  else if ( (*uiter).mean() < (*viter).mean() )
599  iUV=0;
600  else
601  iUV=1;
602 
604  StEEmcIUSmdCluster &cluster=(iUV==0)?ucl:vcl;
605 
607  std::vector<Int_t> matched=(iUV==0)?
608  u2p[cluster.key()]:
609  v2p[cluster.key()];
610 
613  if ( matched.size()==0 || matched.size()==1 )
614  {
615  if ( iUV==0 ) uiter++;
616  if ( iUV==1 ) viter++;
617  continue;
618  }
619 
622  StEEmcIUPoint p=smdpoints[matched.front()];
623  //StEEmcIUPoint p2=smdpoints[matched.front()++];
625  mypoints.push_back(p);
626  //mypoints.push_back(p2);
627 
628  if ( iUV==0 ) uiter++;
629  if ( iUV==1 ) viter++;
630 
631  }
632 
633 
636  chisq=9.0E9;
637  imin=-1;
638  for ( UInt_t i=0; i<mypoints.size(); i++ )
639  {
640 
641  Float_t eu=mypoints[i].cluster(0).energy();
642  Float_t ev=mypoints[i].cluster(1).energy();
643  Float_t x2=fabs((eu-ev)/(eu+ev));
644  if ( x2 < chisq ) {
645  imin=(Int_t)i;
646  chisq=x2;
647  }
648  }
649  float chisq2=9.0E9;
650  int inn=-1;
651 
652  for ( UInt_t i=0; i<mypoints.size(); i++ )
653  {
654  inn=i;
655  if(inn!=imin)
656  {
657 
658  Float_t eu=mypoints[i].cluster(0).energy();
659  Float_t ev=mypoints[i].cluster(1).energy();
660  Float_t x22=(eu-ev)*(eu-ev);
661  if ( x22 < chisq2 ) {
662  chisq2=x22;
663  }
664  }
665  }
666 
667 
668 
675  if ( imin >= 0 ) {
676 
677  StEEmcIUPoint p1=mypoints[imin];
678 
679  removeCluster( uclusters, p1.cluster(0).key() );
680  removeCluster( vclusters, p1.cluster(1).key() );
681 
682  points.push_back(p1);
683 
684  findPoints(sector, uclusters, vclusters, points );
685  return true;
686 
687  }
688 
689 
690  return true;
691 
692 }
693 
694 
695 
696 // ----------------------------------------------------------------------------
697 void StEEmcIUPointMaker::Clear( Option_t *opts )
698 {
699 
700  mEseen=0.;
701  mPoints.clear();
702  mSmdPoints.clear();
703 
704 }
705 
706 // ----------------------------------------------------------------------------
707 void StEEmcIUPointMaker::removeCluster( StEEmcIUSmdClusterVec_t &clusters, Int_t k )
708 {
709  StEEmcIUSmdClusterVec_t::iterator iter=clusters.begin();
710  while ( iter != clusters.end() )
711  {
712  if ( (*iter).key() == k ) {
713  clusters.erase(iter);
714  return;
715  }
716  iter++;
717  }
718 }
719 
720 // ----------------------------------------------------------------------------
721 void StEEmcIUPointMaker::shareEnergy()
722 {
723 
729 
730 
731 
732  Int_t nrow=(Int_t)mPoints.size();
733  Int_t ncol=(Int_t)mPoints.size();
734  if ( nrow < 1 ) return;
735 
736  std::vector<Float_t> Ef(nrow,0.);
737 
738  TMatrixF fractions(nrow,ncol);
739 
741  for ( Int_t k=0; k<mEEanalysis->numberOfHitTowers(0); k++ )
742  {
743 
745  StEEmcTower tower=mEEanalysis->hittower(k,0);
746 
748  for ( UInt_t i=0; i< mPoints.size(); i++ )
749  {
750 
753  Float_t fi=fracp2t( mPoints[i], tower );
754  if ( fi<=0. ) continue;
755 
757  Ef[i] += fi * tower.energy();
758 
759  for ( UInt_t j=0; j<mPoints.size(); j++ )
760  {
761 
764  Float_t fj=fracp2t( mPoints[j], tower );
765  if (fi*fj<=0.) continue;
766 
767  fractions[i][j] += fi*fj;
768 
769  }
770 
771  }
772 
773  }
774 
775  fractions.Print();
776 
778  Double_t det = 0.;
779  TMatrixF invert= fractions;
780  invert.Invert(&det);
781 
782  invert.Print();
783 
784  TMatrixF test= fractions * invert;
785 
786  test.Print();
787 
788 
790 
791  std::vector<Float_t> epoints(nrow,0.);
792 
793  for ( Int_t i=0; i<nrow; i++ )
794  {
795  for ( Int_t j=0; j<ncol; j++ )
796  {
797  epoints[i] += invert[i][j] * Ef[j];
798  }
799 
800 
801  }
802 
803 
804 
805 
806 
807 
808 }
809 
810 
811 // ----------------------------------------------------------------------------
812 Float_t StEEmcIUPointMaker::fracp2t( StEEmcIUPoint &p, StEEmcTower &t )
813 {
814 
816 
817 
820  if ( !t.isNeighbor( p.tower(0) ) ) return 0.;
821  //if ( !(t.index()==p.tower(0).index()) ) return 0.;
822 
823 
825  Float_t xeta=(Float_t)t.etabin();
826  Float_t xphi=(Float_t)t.phibin();
827  Double_t X[]={xphi,xeta};
828 
830  Float_t xeta0=(Float_t)p.tower(0).etabin();
831  Float_t xphi0=(Float_t)p.tower(0).phibin();
832 
835  Int_t sec,sub,eta;
836  Float_t dphi,deta;
837  if ( !mEEtow->getTower(p.position(), sec,sub,eta, dphi,deta ) ) return 0.;
838  dphi/=2.0;
839  deta/=2.0;
840 
842  Float_t xetap=xeta0+deta;
843  Float_t xphip=xphi0+dphi;
844  Double_t P[]={xphip,xetap,1.0};
845 
846  return eeTowerFunction( X, P );
847 
848 }
849 
850 
851 // ----------------------------------------------------------------------------
852 void StEEmcIUPointMaker::shareEnergySimple()
853 {
854 
855  Int_t limit=mLimit; // algo quickly converges on energy
856  Int_t count=0;
857 
858  while ( count++ < limit )
859  {
860 
863  Float_t sumw[720];
864  for (Int_t i=0;i<720;i++) sumw[i]=0.;
865 
868  for ( UInt_t i=0;i<mPoints.size();i++ )
869  {
870 
871  StEEmcIUPoint point=mPoints[i];
872  StEEmcTower tower=point.tower(0);
873 
874  sumw[ tower.index() ] += point.energy() * fracp2t( point, tower );
875 
877  for ( Int_t i=0; i<tower.numberOfNeighbors(); i++ )
878  {
879  StEEmcTower neighbor=tower.neighbor(i);
880  sumw[ neighbor.index() ] += point.energy() * fracp2t( point, neighbor );
881  }
882 
883  }
884 
885 
889 
890  for ( UInt_t i=0;i<mPoints.size();i++ )
891  {
892 
893  StEEmcIUPoint &point=mPoints[i]; // note the reference
894  StEEmcTower tower=point.tower(0);
895  Float_t epoint=0.;
896 
897  Float_t frac=0.;
898  if ( !tower.fail() && !tower.stat() && sumw[tower.index()]>0. )
899  frac = point.energy() * fracp2t(point,tower) / sumw[ tower.index() ];
900 
901  epoint += tower.energy() * frac;
902 
905  for ( Int_t i=0; i<tower.numberOfNeighbors(); i++ )
906  {
907  StEEmcTower neighbor=tower.neighbor(i);
908  if ( neighbor.stat() || neighbor.fail() || sumw[neighbor.index()]<=0. ) continue;
909  frac = point.energy() * fracp2t(point,neighbor) / sumw[ neighbor.index() ];
910  epoint += frac * neighbor.energy();
911  }
912 
913 
915  point.energy( epoint );
916 
917  }
918 
919  }
920 
921 
922 
923 
924  std::vector<Bool_t> seen(720,false);
925  for ( UInt_t i=0;i<mPoints.size();i++ )
926  {
927 
928  StEEmcTower tow=mPoints[i].tower(0);
929  if ( !seen[ tow.index() ] ) mEseen+=tow.energy();
930  seen[ tow.index() ] = true;
931 
932  for ( Int_t j=0;j<tow.numberOfNeighbors();j++ )
933  {
934  StEEmcTower nei=tow.neighbor(j);
935  if ( !seen[ nei.index() ] ) mEseen += nei.energy();
936  seen[ nei.index() ] = true;
937  }
938 
939  }
940 
941 
942 
943 
944 
945 
946 }
947 
948 
949 
950 
951 
952 // ----------------------------------------------------------------------------
953 void StEEmcIUPointMaker::fillStEvent()
954 {
955 
956 
957  StEvent *stevent=(StEvent*)GetInputDS("StEvent");
958  if ( !stevent ) {
959  Warning("fillStEvent","called, but no StEvent to be found");
960  return;
961  }
962 
963 
965  for ( UInt_t i=0; i<mPoints.size(); i++ )
966  {
967 
968  StEmcPoint *point=mPoints[i].stemc();
969  stevent->emcCollection()->addEndcapPoint( point );
970 
971  mEtoEE[ point ] = mPoints[i];
972 
973 
974  }
975 
976 }
977 
978 
979 
980 // ----------------------------------------------------------------------------
981 void StEEmcIUPointMaker::verifyStEvent()
982 {
983 
984  Float_t emc_sum_points = 0.;
985  Float_t eemc_sum_points = 0.;
986 
987  StEvent *stevent=(StEvent*)GetInputDS("StEvent");
988  if ( !stevent ) {
989  Warning("verifyStEvent","called, but no StEvent to be found");
990  return;
991  }
992 
993  StSPtrVecEmcPoint& emcpts = stevent->emcCollection()->endcapPoints();
994  for ( UInt_t i=0;i<emcpts.size();i++ )
995  {
996 
997  StEmcPoint *p=emcpts[i];
998  assert(p);
999  emc_sum_points += p->energy();
1000 
1001  }
1002 
1003  for ( UInt_t i=0;i<mPoints.size();i++ )
1004  {
1005 
1006  StEEmcIUPoint p=mPoints[i];
1007  eemc_sum_points += p.energy();
1008 
1009  }
1010 
1011  std::cout << "StEEmcIUPointMaker point checksum: ";
1012  if ( emc_sum_points == eemc_sum_points )
1013  {
1014  std::cout << "passed";
1015  }
1016  else
1017  std::cout << "FAILED";
1018  std::cout << std::endl;
1019 
1020 
1021 }
1022 
1023 
1024 // ----------------------------------------------------------------------------
1025 void StEEmcIUPointMaker::countRelatives()
1026 {
1027 
1029  Int_t npoints[720];
1030  for ( Int_t i=0;i<720;i++ ) npoints[i]=0;
1031 
1032  for ( UInt_t i=0;i<mPoints.size();i++ )
1033  npoints[ mPoints[i].tower(0).index() ]++;
1034 
1036  for ( UInt_t i=0;i<mPoints.size();i++ )
1037  {
1038 
1039  StEEmcTower tower=mPoints[i].tower(0);
1040  Int_t nn=tower.numberOfNeighbors();
1041 
1042  Int_t nrel=npoints[ tower.index() ] - 1; // don't count self
1043  assert(nrel>=0); // pbck
1044 
1045  for ( Int_t j=0;j<nn;j++ )
1046  {
1047  StEEmcTower t2=tower.neighbor(j);
1048  nrel+=npoints[ t2.index() ];
1049  }
1050 
1051  mPoints[i].numberOfRelatives(nrel);
1052 
1053  }
1054 
1055 
1056 }
1057 
1058 
1059 
1060 
1061 // ----------------------------------------------------------------------------
1062 void StEEmcIUPointMaker::shareEnergySmd()
1063 {
1064 
1067  Float_t sumw[720];for ( Int_t i=0;i<720;i++ )sumw[i]=0.;
1068  Float_t sumw1[720];for ( Int_t i=0;i<720;i++ )sumw1[i]=0.;
1069 
1070  for ( UInt_t ipoint=0;ipoint<mPoints.size();ipoint++ )
1071  {
1072 
1073  StEEmcIUPoint point=mPoints[ipoint];
1074  StEEmcTower tower=point.tower(0);
1075  sumw[tower.index()]+=point.energy();
1076  sumw1[tower.index()]+=point.energy();
1077  //printf("en==%f index=%d sumw=%f\n",point.energy(),tower.index(),sumw[tower.index()]);
1078  for ( Int_t itow=0;itow<tower.numberOfNeighbors();itow++ )
1079  {
1080  StEEmcTower t=tower.neighbor(itow);
1081  Int_t index=t.index();
1082 
1083  sumw[index]+=point.energy();
1084  //printf("index=%d sumw=%f\n",index,sumw[index]);
1085  }
1086 
1087  }
1088  //printf("end of first loop.........................\n");
1091 
1092  for ( UInt_t ipoint=0;ipoint<mPoints.size();ipoint++ )
1093  {
1094 
1095  StEEmcIUPoint point=mPoints[ipoint]; // note reference
1096  StEEmcTower tower = point.tower(0);
1097  StEEmcTower pre1 = mEEanalysis->tower(tower.index(),1);
1098  StEEmcTower pre2 = mEEanalysis->tower(tower.index(),2);
1099  StEEmcTower post = mEEanalysis->tower(tower.index(),3);
1100  //hPostE->Fill(post.energy()*1000.0);
1101  Float_t epoint = 0.;
1102  Float_t epre1 = 0.;
1103  Float_t epre2 = 0.;
1104  Float_t epost = 0.;
1105 
1106  Int_t index = tower.index();
1107  epoint += tower.energy() * point.energy() / sumw[index];
1108  //printf("start of the calculation of a point\n");
1109  //printf("tower.energy=%f epoint=%f smde=%f sumw=%f index=%d\n",tower.energy(),epoint,point.energy(),sumw[index],index);
1110  epre1 += pre1.energy() * point.energy() / sumw1[index];
1111  epre2 += pre2.energy() * point.energy() / sumw1[index];
1112  epost += post.energy() * point.energy() / sumw[index];
1113  //printf("start of the neighbor calculation of a point\n");
1114  for ( Int_t itow=0;itow<tower.numberOfNeighbors();itow++ )
1115  {
1116  StEEmcTower t=tower.neighbor(itow);
1117 // StEEmcTower p=pre1.neighbor(itow);
1118 // StEEmcTower q=pre2.neighbor(itow);
1119  StEEmcTower r=post.neighbor(itow);
1120 
1121  index = t.index();
1122 
1123  epoint += t.energy() * point.energy() / sumw[index];
1124 
1125  //printf("neighbortower.energy=%f epoint=%f smde=%f sumw=%f index=%d\n",t.energy(),epoint,point.energy(),sumw[index],index);
1126 // epre1 += p.energy() * point.energy() / sumw[index];
1127 // epre2 += q.energy() * point.energy() / sumw[index];
1128 // epost += r.energy() * point.energy() / sumw[index];
1129  }
1130 
1131 
1132  mPoints[ipoint].energy(epoint);
1133  mPoints[ipoint].energy(epre1,1);
1134  mPoints[ipoint].energy(epre2,2);
1135  mPoints[ipoint].energy(epost,3);
1136 
1137 
1138  }
1139 
1140 
1141 }
1142 
1143 
StEEmcIUPointMaker::findPoints
Bool_t findPoints(Int_t sector, StEEmcIUSmdClusterVec_t u, StEEmcIUSmdClusterVec_t v, StEEmcIUPointVec_t &p)
find points in the endcap
Definition: StEEmcIUPointMaker.cxx:282
StEEmcTower::numberOfNeighbors
Int_t numberOfNeighbors() const
get the number of neighboring towers
Definition: StEEmcTower.h:54
StEEmcA2EMaker
EEmc ADC –&gt; energy maker.
Definition: StEEmcA2EMaker.h:19
StEEmcIUSmdCluster::key
Int_t key()
Return a unique key assigned by the cluster maker.
Definition: StEEmcIUSmdCluster.h:26
StEEmcElement::stat
void stat(unsigned s)
Set a status bit for this element.
Definition: StEEmcElement.h:23
StEEmcTower::isNeighbor
Bool_t isNeighbor(const StEEmcTower &t) const
Definition: StEEmcTower.cxx:98
StEEmcIUPointMaker::shareEnergySimple
void shareEnergySimple()
Divide energy of eemc towers between identified smd points (doesn&#39;t work as well as smd algo) ...
Definition: StEEmcIUPointMaker.cxx:852
StEEmcA2EMaker::hittower
StEEmcTower & hittower(Int_t hit, Int_t layer)
Definition: StEEmcA2EMaker.h:63
StEEmcIUPointMaker::Clear
void Clear(Option_t *opts="")
Clear old points.
Definition: StEEmcIUPointMaker.cxx:697
StEEmcIUPointMaker::mSmdPoints
StEEmcIUPointVec_t mSmdPoints
SMD only points.
Definition: StEEmcIUPointMaker.h:123
StEEmcIUPointMaker::mEseen
Float_t mEseen
Energy seen by the algorithm.
Definition: StEEmcIUPointMaker.h:163
StEEmcTower::neighbor
void neighbor(StEEmcTower *n)
add a tower to list of neighbors
Definition: StEEmcTower.h:52
StEEmcIUPoint::energy
void energy(Float_t e, Int_t layer=0)
Set the energy of this point.
Definition: StEEmcIUPoint.h:25
StEEmcIUClusterMaker::smdclusters
StEEmcIUSmdClusterVec_t smdclusters(Int_t sec, Int_t plane)
Return a vector of smd clusters.
Definition: StEEmcIUClusterMaker.h:220
StEEmcIUPointMaker::mFillStEvent
Bool_t mFillStEvent
Option to fill StEvent.
Definition: StEEmcIUPointMaker.h:165
StEEmcIUPoint::position
void position(TVector3 p)
Set the position of this point at the SMD plane.
Definition: StEEmcIUPoint.h:23
StEEmcIUPointMaker::points
StEEmcIUPointVec_t points()
Return a unique key assigned by the cluster maker.
Definition: StEEmcIUPointMaker.h:182
StEEmcIUPointMaker::mEEmap
EEmcSmdMap * mEEmap
Tower to smd map.
Definition: StEEmcIUPointMaker.h:110
StEEmcIUPointMaker::mEtoEE
std::map< StEmcPoint *, StEEmcIUPoint > mEtoEE
Map connecting StEEmcIUPoint to StEmcPoint.
Definition: StEEmcIUPointMaker.h:172
StEEmcIUPointMaker::Init
Int_t Init()
Initialize.
Definition: StEEmcIUPointMaker.cxx:77
StEEmcIUPointMaker::buildSmdPoints
StEEmcIUPointVec_t buildSmdPoints(Int_t sector, StEEmcIUSmdClusterVec_t &u, StEEmcIUSmdClusterVec_t &v)
build smd points and associations between smd points and clusters
Definition: StEEmcIUPointMaker.cxx:172
StEEmcIUPointMaker
Class for building points from smd clusters.
Definition: StEEmcIUPointMaker.h:21
StEEmcIUPointMaker::fillStEvent
void fillStEvent()
Fills the StEmcPoint collection.
Definition: StEEmcIUPointMaker.cxx:953
StEEmcIUClusterMaker
A cluster maker for the EEMC.
Definition: StEEmcIUClusterMaker.h:17
StEEmcIUPointMaker::removeCluster
void removeCluster(StEEmcIUSmdClusterVec_t &clusters, Int_t key)
Remove a cluster from the list of clusters.
Definition: StEEmcIUPointMaker.cxx:707
StEEmcA2EMaker::numberOfHitTowers
Int_t numberOfHitTowers(Int_t layer) const
Definition: StEEmcA2EMaker.h:58
StEEmcTower::etabin
Int_t etabin() const
Returns the etabin of this tower, pre- or postshower element.
Definition: StEEmcTower.h:45
StEEmcElement::fail
void fail(unsigned f)
Set a fail bit for this element.
Definition: StEEmcElement.h:25
StEEmcIUPointMaker::verifyStEvent
void verifyStEvent()
Checks that StEvent is properly saved.
Definition: StEEmcIUPointMaker.cxx:981
StEEmcIUPointMaker::Make
Int_t Make()
Build points for this event.
Definition: StEEmcIUPointMaker.cxx:90
StEEmcIUPoint::cluster
void cluster(StEEmcIUSmdCluster c, Int_t plane)
Add an smd cluster to this point.
Definition: StEEmcIUPoint.h:31
StEEmcIUPointMaker::mEnergyMode
Int_t mEnergyMode
Option for dividing energy.
Definition: StEEmcIUPointMaker.h:167
StEEmcTower::index
void index(Int_t i)
Definition: StEEmcTower.cxx:76
StEEmcTower
Base class for representing tower, preshower and postshower elements.
Definition: StEEmcTower.h:11
StEEmcIUPointMaker::shareEnergySmd
void shareEnergySmd()
Divide energy of eemc towers between identified smd points in proportion to the smd energy...
Definition: StEEmcIUPointMaker.cxx:1062
StEEmcIUPointMaker::mEEanalysis
StEEmcA2EMaker * mEEanalysis
ADC2E.
Definition: StEEmcIUPointMaker.h:101
StEEmcIUPointMaker::StEEmcIUPointMaker
StEEmcIUPointMaker(const Char_t *name)
Definition: StEEmcIUPointMaker.cxx:55
StEEmcIUSmdCluster::mean
Float_t mean()
Return the mean strip number of this cluster.
Definition: StEEmcIUSmdCluster.h:153
StEEmcIUPointMaker::point
StEEmcIUPoint point(Int_t ipoint)
Return a specified point.
Definition: StEEmcIUPointMaker.h:184
StEEmcIUPointMaker::mLimit
Int_t mLimit
How many iterations for the tower energy sharing mode.
Definition: StEEmcIUPointMaker.h:169
EEmcGeomSimple
EEMC simple geometry.
Definition: EEmcGeomSimple.h:22
kStOK
Definition: Stypes.h:40
StEEmcIUSmdCluster
A base class for representing clusters of EEMC smd strips.
Definition: StEEmcIUSmdCluster.h:11
StEEmcIUPoint
Base class for representing EEMC points.
Definition: StEEmcIUPoint.h:12
StEEmcIUPointMaker::mEEsmd
EEmcSmdGeom * mEEsmd
Smd geometry.
Definition: StEEmcIUPointMaker.h:108
StEEmcTower::phibin
Int_t phibin() const
Returns the phibin of this tower.
Definition: StEEmcTower.h:47
EEmcSmdGeom::getIntersection
TVector3 getIntersection(Int_t iSec, Int_t iUStrip, Int_t iVStrip, const TVector3 &vertex) const
Definition: EEmcSmdGeom.cxx:628
StEEmcIUPointMaker::mEEtow
EEmcGeomSimple * mEEtow
Tower geometry.
Definition: StEEmcIUPointMaker.h:106
StEEmcIUPoint::tower
void tower(StEEmcTower t, Float_t w=1.)
Add a tower with specified weight to the point.
Definition: StEEmcIUPoint.h:29
StEEmcA2EMaker::tower
StEEmcTower & tower(Int_t index, Int_t layer=0)
Definition: StEEmcA2EMaker.h:69
StEEmcIUPointMaker::mPoints
StEEmcIUPointVec_t mPoints
All fully reconstructed points.
Definition: StEEmcIUPointMaker.h:120
StEEmcIUPointMaker::fracp2t
Float_t fracp2t(StEEmcIUPoint &p, StEEmcTower &t)
Definition: StEEmcIUPointMaker.cxx:812
EEmcGeomSimple::getTower
bool getTower(const TVector3 &r, int &sec, int &sub, int &etabin, Float_t &dphi, Float_t &deta) const
Definition: EEmcGeomSimple.cxx:159
StEEmcElement::energy
void energy(Float_t e)
Set the energy (adc-ped+0.5)/gain for this element.
Definition: StEEmcElement.h:21
StEEmcIUPointMaker::shareEnergy
void shareEnergy()
Divide energy of eemc towers between identified smd points using fit (doesn&#39;t work) ...
Definition: StEEmcIUPointMaker.cxx:721
StEEmcIUPointMaker::mEEclusters
StEEmcIUClusterMaker * mEEclusters
Clusters.
Definition: StEEmcIUPointMaker.h:103
StEEmcIUSmdCluster::energy
Float_t energy()
Return the energy of this cluster.
Definition: StEEmcIUSmdCluster.h:152
StEEmcIUPointMaker::countRelatives
void countRelatives()
Determine the number of points which share tower energy with another point.
Definition: StEEmcIUPointMaker.cxx:1025
Generated by   doxygen 1.8.5