1 #include "StEEmcPointMaker.h"
2 #include "StEEmcA2EMaker.h"
3 #include "StEEmcClusterMaker.h"
5 #include "StEEmcUtil/EEmcGeom/EEmcGeomDefs.h"
6 #include "StEEmcUtil/EEmcGeom/EEmcGeomSimple.h"
7 #include "StEEmcUtil/StEEmcSmd/EEmcSmdGeom.h"
13 #include "eeTowerFunction.h"
16 #include "StEvent/StEvent.h"
17 #include "StEvent/StEmcCollection.h"
18 #include "StEvent/StEmcPoint.h"
30 std::cout <<
"StEEmcPointMaker("<<name<<
")" << std::endl;
36 mEEsmd=EEmcSmdGeom::instance();
37 mEEmap=EEmcSmdMap::instance();
53 return StMaker::Init();
68 for ( Int_t sector=0; sector<12; sector++ )
76 std::sort( uclusters.begin(), uclusters.end(), inner );
77 std::sort( vclusters.begin(), vclusters.end(), inner );
90 for ( UInt_t i=0;i<
mPoints.size();i++ )
98 StEEmcPointVec_t orgpoints=
mPoints;
126 StEEmcSmdClusterVec_t &u,
127 StEEmcSmdClusterVec_t &v )
132 for ( UInt_t i=0; i<u.size(); i++ )
136 Float_t xu=uc.mean();
138 for ( UInt_t j=0;j<v.size(); j++ )
143 Float_t xv=vc.mean();
173 for ( Int_t layer=1;layer<=3;layer++ )
185 p.
energy( uc.energy() + vc.energy() );
218 StEEmcSmdClusterVec_t uclusters,
219 StEEmcSmdClusterVec_t vclusters,
220 StEEmcPointVec_t &points )
225 StEEmcPointVec_t mypoints;
228 std::sort( uclusters.begin(), uclusters.end(), inner );
229 std::sort( vclusters.begin(), vclusters.end(), inner );
232 StEEmcPointVec_t smdpoints =
buildSmdPoints( sector, uclusters, vclusters );
236 if ( smdpoints.size() < 1 )
return false;
239 std::sort( smdpoints.begin(), smdpoints.end(), chiSquare );
245 std::map< Int_t, std::vector<Int_t> > u2p, v2p;
246 for ( UInt_t i=0; i<smdpoints.size(); i++ )
248 u2p[ smdpoints[i].cluster(0).key() ].push_back( i );
249 v2p[ smdpoints[i].cluster(1).key() ].push_back( i );
257 StEEmcSmdClusterVec_t::iterator uiter=uclusters.begin();
258 StEEmcSmdClusterVec_t::iterator viter=vclusters.begin();
266 while ( uiter<uclusters.end() || viter<vclusters.end() )
272 if ( uiter<uclusters.end() ) ucl=(*uiter);
273 if ( viter<vclusters.end() ) vcl=(*viter);
277 else if ( vcl.
key()<0 )
279 else if ( (*uiter).mean() < (*viter).mean() )
288 std::vector<Int_t> matched=(iUV==0)?
295 if ( matched.size()==0 || matched.size() >1 )
297 if ( iUV==0 ) uiter++;
298 if ( iUV==1 ) viter++;
313 mypoints.push_back(p);
315 if ( iUV==0 ) uiter++;
316 if ( iUV==1 ) viter++;
325 for ( UInt_t i=0; i<mypoints.size(); i++ )
327 Float_t eu=mypoints[i].
cluster(0).energy();
328 Float_t ev=mypoints[i].cluster(1).energy();
329 Float_t x2=(eu-ev)*(eu-ev);
349 findPoints(sector, uclusters, vclusters, points );
367 uiter=uclusters.begin();
368 viter=vclusters.begin();
369 while ( uiter!=uclusters.end() || viter!=vclusters.end() )
375 if ( uiter!=uclusters.end() ) ucl=(*uiter);
376 if ( viter!=vclusters.end() ) vcl=(*viter);
380 else if ( vcl.
key()<0 )
382 else if ( (*uiter).mean() < (*viter).mean() )
391 std::vector<Int_t> matched=(iUV==0)?
397 if ( matched.size()==0 || matched.size()==1 )
399 if ( iUV==0 ) uiter++;
400 if ( iUV==1 ) viter++;
409 mypoints.push_back(p);
411 if ( iUV==0 ) uiter++;
412 if ( iUV==1 ) viter++;
421 for ( UInt_t i=0; i<mypoints.size(); i++ )
423 Float_t eu=mypoints[i].
cluster(0).energy();
424 Float_t ev=mypoints[i].cluster(1).energy();
425 Float_t x2=(eu-ev)*(eu-ev);
446 findPoints(sector, uclusters, vclusters, points );
470 StEEmcSmdClusterVec_t::iterator iter=clusters.begin();
471 while ( iter != clusters.end() )
473 if ( (*iter).key() == k ) {
474 clusters.erase(iter);
493 Int_t nrow=(Int_t)
mPoints.size();
494 Int_t ncol=(Int_t)
mPoints.size();
495 if ( nrow < 1 )
return;
497 std::vector<Float_t> Ef(nrow,0.);
499 TMatrixF fractions(nrow,ncol);
509 for ( UInt_t i=0; i<
mPoints.size(); i++ )
515 if ( fi<=0. )
continue;
518 Ef[i] += fi * tower.
energy();
520 for ( UInt_t j=0; j<
mPoints.size(); j++ )
526 if (fi*fj<=0.)
continue;
528 fractions[i][j] += fi*fj;
540 TMatrixF invert= fractions;
545 TMatrixF test= fractions * invert;
552 std::vector<Float_t> epoints(nrow,0.);
554 for ( Int_t i=0; i<nrow; i++ )
556 for ( Int_t j=0; j<ncol; j++ )
558 epoints[i] += invert[i][j] * Ef[j];
586 Float_t xeta=(Float_t)t.
etabin();
587 Float_t xphi=(Float_t)t.
phibin();
588 Double_t X[]={xphi,xeta};
591 Float_t xeta0=(Float_t)p.
tower(0).etabin();
592 Float_t xphi0=(Float_t)p.
tower(0).phibin();
603 Float_t xetap=xeta0+deta;
604 Float_t xphip=xphi0+dphi;
605 Double_t P[]={xphip,xetap,1.0};
607 return eeTowerFunction( X, P );
619 while ( count++ < limit )
625 for (Int_t i=0;i<720;i++) sumw[i]=0.;
629 for ( UInt_t i=0;i<
mPoints.size();i++ )
651 for ( UInt_t i=0;i<
mPoints.size();i++ )
659 if ( !tower.
fail() && !tower.
stat() && sumw[tower.
index()]>0. )
662 epoint += tower.
energy() * frac;
669 if ( neighbor.
stat() || neighbor.
fail() || sumw[neighbor.
index()]<=0. )
continue;
671 epoint += frac * neighbor.
energy();
685 std::vector<Bool_t> seen(720,
false);
686 for ( UInt_t i=0;i<
mPoints.size();i++ )
691 seen[ tow.
index() ] =
true;
697 seen[ nei.
index() ] =
true;
720 Warning(
"fillStEvent",
"called, but no StEvent to be found");
726 for ( UInt_t i=0; i<
mPoints.size(); i++ )
730 stevent->emcCollection()->addEndcapPoint( point );
745 Float_t emc_sum_points = 0.;
746 Float_t eemc_sum_points = 0.;
750 Warning(
"verifyStEvent",
"called, but no StEvent to be found");
754 StSPtrVecEmcPoint& emcpts = stevent->emcCollection()->endcapPoints();
755 for ( UInt_t i=0;i<emcpts.size();i++ )
760 emc_sum_points += p->energy();
764 for ( UInt_t i=0;i<
mPoints.size();i++ )
768 eemc_sum_points += p.
energy();
772 std::cout <<
"StEEmcPointMaker point checksum: ";
773 if ( emc_sum_points == eemc_sum_points )
775 std::cout <<
"passed";
778 std::cout <<
"FAILED";
779 std::cout << std::endl;
791 for ( Int_t i=0;i<720;i++ ) npoints[i]=0;
793 for ( UInt_t i=0;i<
mPoints.size();i++ )
797 for ( UInt_t i=0;i<
mPoints.size();i++ )
803 Int_t nrel=npoints[ tower.
index() ] - 1;
806 for ( Int_t j=0;j<nn;j++ )
809 nrel+=npoints[ t2.
index() ];
812 mPoints[i].numberOfRelatives(nrel);
828 Float_t sumw[720];
for ( Int_t i=0;i<720;i++ )sumw[i]=0.;
830 for ( UInt_t ipoint=0;ipoint<
mPoints.size();ipoint++ )
841 Int_t index=t.
index();
842 sumw[index]+=point.
energy();
850 for ( UInt_t ipoint=0;ipoint<
mPoints.size();ipoint++ )
856 Int_t index = tower.
index();
857 epoint += tower.
energy() * point.
energy() / sumw[index];
867 mPoints[ipoint].energy(epoint);
StEEmcA2EMaker * mEEanalysis
ADC2E.
void cluster(const StEEmcSmdCluster &c, Int_t plane)
Add an smd cluster to this point.
void shareEnergySmd()
Divide energy of eemc towers between identified smd points in proportion to the smd energy...
Int_t numberOfNeighbors() const
get the number of neighboring towers
EEmc ADC –> energy maker.
Base class for representing EEMC points.
void stat(unsigned s)
Set a status bit for this element.
Bool_t isNeighbor(const StEEmcTower &t) const
Float_t mEseen
Energy seen by the algorithm.
StEEmcTower & hittower(Int_t hit, Int_t layer)
void energy(Float_t e, Int_t layer=0)
Set the energy of this point.
StEEmcPointVec_t mPoints
All fully reconstructed points.
Bool_t mFillStEvent
Option to fill StEvent.
void shareEnergySimple()
Divide energy of eemc towers between identified smd points (doesn't work as well as smd algo) ...
void neighbor(StEEmcTower *n)
add a tower to list of neighbors
void verifyStEvent()
Checks that StEvent is properly saved.
StEEmcPoint point(Int_t ipoint)
Return a specified point.
EEmcSmdMap * mEEmap
Tower to smd map.
StEEmcPointVec_t points()
Return vector of all points found in endcap.
void shareEnergy()
Divide energy of eemc towers between identified smd points using fit (doesn't work) ...
Class for building points from smd clusters.
Int_t numberOfHitTowers(Int_t layer) const
Int_t etabin() const
Returns the etabin of this tower, pre- or postshower element.
void fail(unsigned f)
Set a fail bit for this element.
StEEmcPointMaker(const Char_t *name)
std::map< StEmcPoint *, StEEmcPoint > mEtoEE
Map connecting StEEmcPoint to StEmcPoint.
Base class for representing tower, preshower and postshower elements.
StEEmcPointVec_t buildSmdPoints(Int_t sector, StEEmcSmdClusterVec_t &u, StEEmcSmdClusterVec_t &v)
build smd points and associations between smd points and clusters
Bool_t findPoints(Int_t sector, StEEmcSmdClusterVec_t u, StEEmcSmdClusterVec_t v, StEEmcPointVec_t &p)
find points in the endcap
void countRelatives()
Determine the number of points which share tower energy with another point.
Int_t key()
Return a unique key assigned by the cluster maker.
A base class for representing clusters of EEMC smd strips.
void removeCluster(StEEmcSmdClusterVec_t &clusters, Int_t key)
Remove a cluster from the list of clusters.
void tower(const StEEmcTower &t, Float_t w=1.)
Add a tower with specified weight to the point.
Int_t phibin() const
Returns the phibin of this tower.
TVector3 getIntersection(Int_t iSec, Int_t iUStrip, Int_t iVStrip, const TVector3 &vertex) const
void position(const TVector3 &p)
Set the position of this point at the SMD plane.
EEmcGeomSimple * mEEtow
Tower geometry.
StEEmcTower & tower(Int_t index, Int_t layer=0)
bool getTower(const TVector3 &r, int &sec, int &sub, int &etabin, Float_t &dphi, Float_t &deta) const
void Clear(Option_t *opts="")
Clear old points.
EEmcSmdGeom * mEEsmd
Smd geometry.
Int_t mLimit
How many iterations for the tower energy sharing mode.
void energy(Float_t e)
Set the energy (adc-ped+0.5)/gain for this element.
Float_t fracp2t(StEEmcPoint &p, StEEmcTower &t)
void fillStEvent()
Fills the StEmcPoint collection.
StEEmcClusterMaker * mEEclusters
Clusters.
Int_t Make()
Build points for this event.
Int_t mEnergyMode
Option for dividing energy.
StEEmcSmdClusterVec_t & smdclusters(Int_t sec, Int_t plane)
Return a vector of smd clusters.
A cluster maker for the EEMC.