26 #include "StEEmcPointMaker.h"
27 #include "StEEmcPool/StEEmcA2EMaker/StEEmcA2EMaker.h"
28 #include "StEEmcPool/StEEmcClusterMaker/StEEmcClusterMaker.h"
30 #include "StEEmcUtil/EEmcGeom/EEmcGeomDefs.h"
31 #include "StEEmcUtil/EEmcGeom/EEmcGeomSimple.h"
32 #include "StEEmcUtil/StEEmcSmd/EEmcSmdGeom.h"
38 #include "eeTowerFunction.h"
41 #include "StEvent/StEvent.h"
42 #include "StEvent/StEmcCollection.h"
43 #include "StEvent/StEmcPoint.h"
55 std::cout <<
"StEEmcPointMaker("<<name<<
")" << std::endl;
61 mEEsmd=EEmcSmdGeom::instance();
62 mEEmap=EEmcSmdMap::instance();
79 return StMaker::Init();
94 for ( Int_t sector=0; sector<12; sector++ )
102 std::sort( uclusters.begin(), uclusters.end(), inner );
103 std::sort( vclusters.begin(), vclusters.end(), inner );
116 for ( UInt_t i=0;i<
mPoints.size();i++ )
124 StEEmcPointVec_t orgpoints=
mPoints;
159 StEEmcSmdClusterVec_t &u,
160 StEEmcSmdClusterVec_t &v )
165 for ( UInt_t i=0; i<u.size(); i++ )
169 Float_t xu=uc.mean();
171 for ( UInt_t j=0;j<v.size(); j++ )
176 Float_t xv=vc.mean();
206 for ( Int_t layer=1;layer<=3;layer++ )
220 Float_t Eu=uc.energy();
221 Float_t Ev=vc.energy();
229 p.
energy( uc.energy() + vc.energy() );
266 StEEmcSmdClusterVec_t uclusters,
267 StEEmcSmdClusterVec_t vclusters,
268 StEEmcPointVec_t &points )
273 StEEmcPointVec_t mypoints;
276 std::sort( uclusters.begin(), uclusters.end(), inner );
277 std::sort( vclusters.begin(), vclusters.end(), inner );
280 StEEmcPointVec_t smdpoints =
buildSmdPoints( sector, uclusters, vclusters );
284 if ( smdpoints.size() < 1 )
return false;
287 std::sort( smdpoints.begin(), smdpoints.end(), chiSquare );
293 std::map< Int_t, std::vector<Int_t> > u2p, v2p;
294 for ( UInt_t i=0; i<smdpoints.size(); i++ )
296 u2p[ smdpoints[i].cluster(0).key() ].push_back( i );
297 v2p[ smdpoints[i].cluster(1).key() ].push_back( i );
305 StEEmcSmdClusterVec_t::iterator uiter=uclusters.begin();
306 StEEmcSmdClusterVec_t::iterator viter=vclusters.begin();
314 while ( uiter<uclusters.end() || viter<vclusters.end() )
320 if ( uiter<uclusters.end() ) ucl=(*uiter);
321 if ( viter<vclusters.end() ) vcl=(*viter);
325 else if ( vcl.
key()<0 )
327 else if ( (*uiter).mean() < (*viter).mean() )
336 std::vector<Int_t> matched=(iUV==0)?
343 if ( matched.size()==0 || matched.size() >1 )
345 if ( iUV==0 ) uiter++;
346 if ( iUV==1 ) viter++;
361 mypoints.push_back(p);
363 if ( iUV==0 ) uiter++;
364 if ( iUV==1 ) viter++;
373 for ( UInt_t i=0; i<mypoints.size(); i++ )
375 Float_t eu=mypoints[i].cluster(0).energy();
376 Float_t ev=mypoints[i].cluster(1).energy();
377 Float_t x2=(eu-ev)*(eu-ev);
397 findPoints(sector, uclusters, vclusters, points );
415 uiter=uclusters.begin();
416 viter=vclusters.begin();
417 while ( uiter!=uclusters.end() || viter!=vclusters.end() )
423 if ( uiter!=uclusters.end() ) ucl=(*uiter);
424 if ( viter!=vclusters.end() ) vcl=(*viter);
428 else if ( vcl.
key()<0 )
430 else if ( (*uiter).mean() < (*viter).mean() )
439 std::vector<Int_t> matched=(iUV==0)?
445 if ( matched.size()==0 || matched.size()==1 )
447 if ( iUV==0 ) uiter++;
448 if ( iUV==1 ) viter++;
457 mypoints.push_back(p);
459 if ( iUV==0 ) uiter++;
460 if ( iUV==1 ) viter++;
469 for ( UInt_t i=0; i<mypoints.size(); i++ )
471 Float_t eu=mypoints[i].cluster(0).energy();
472 Float_t ev=mypoints[i].cluster(1).energy();
473 Float_t x2=(eu-ev)*(eu-ev);
494 findPoints(sector, uclusters, vclusters, points );
519 StEEmcSmdClusterVec_t::iterator iter=clusters.begin();
520 while ( iter != clusters.end() )
522 if ( (*iter).key() == k ) {
523 clusters.erase(iter);
542 Int_t nrow=(Int_t)
mPoints.size();
543 Int_t ncol=(Int_t)
mPoints.size();
544 if ( nrow < 1 )
return;
546 std::vector<Float_t> Ef(nrow,0.);
548 TMatrixF fractions(nrow,ncol);
558 for ( UInt_t i=0; i<
mPoints.size(); i++ )
564 if ( fi<=0. )
continue;
567 Ef[i] += fi * tower.
energy();
569 for ( UInt_t j=0; j<
mPoints.size(); j++ )
575 if (fi*fj<=0.)
continue;
577 fractions[i][j] += fi*fj;
589 TMatrixF invert= fractions;
594 TMatrixF test= fractions * invert;
601 std::vector<Float_t> epoints(nrow,0.);
603 for ( Int_t i=0; i<nrow; i++ )
605 for ( Int_t j=0; j<ncol; j++ )
607 epoints[i] += invert[i][j] * Ef[j];
635 Float_t xeta=(Float_t)t.
etabin();
636 Float_t xphi=(Float_t)t.
phibin();
637 Double_t X[]={xphi,xeta};
640 Float_t xeta0=(Float_t)p.
tower(0).etabin();
641 Float_t xphi0=(Float_t)p.
tower(0).phibin();
652 Float_t xetap=xeta0+deta;
653 Float_t xphip=xphi0+dphi;
654 Double_t P[]={xphip,xetap,1.0};
656 return eeTowerFunction( X, P );
668 while ( count++ < limit )
674 for (Int_t i=0;i<720;i++) sumw[i]=0.;
678 for ( UInt_t i=0;i<
mPoints.size();i++ )
700 for ( UInt_t i=0;i<
mPoints.size();i++ )
708 if ( !tower.
fail() && !tower.
stat() && sumw[tower.
index()]>0. )
711 epoint += tower.
energy() * frac;
718 if ( neighbor.
stat() || neighbor.
fail() || sumw[neighbor.
index()]<=0. )
continue;
720 epoint += frac * neighbor.
energy();
734 std::vector<Bool_t> seen(720,
false);
735 for ( UInt_t i=0;i<
mPoints.size();i++ )
740 seen[ tow.
index() ] =
true;
746 seen[ nei.
index() ] =
true;
769 Warning(
"fillStEvent",
"called, but no StEvent to be found");
775 for ( UInt_t i=0; i<
mPoints.size(); i++ )
779 stevent->emcCollection()->addEndcapPoint( point );
794 Float_t emc_sum_points = 0.;
795 Float_t eemc_sum_points = 0.;
799 Warning(
"verifyStEvent",
"called, but no StEvent to be found");
803 StSPtrVecEmcPoint& emcpts = stevent->emcCollection()->endcapPoints();
804 for ( UInt_t i=0;i<emcpts.size();i++ )
809 emc_sum_points += p->energy();
813 for ( UInt_t i=0;i<
mPoints.size();i++ )
817 eemc_sum_points += p.
energy();
821 std::cout <<
"StEEmcPointMaker point checksum: ";
822 if ( emc_sum_points == eemc_sum_points )
824 std::cout <<
"passed";
827 std::cout <<
"FAILED";
828 std::cout << std::endl;
840 for ( Int_t i=0;i<720;i++ ) npoints[i]=0;
842 for ( UInt_t i=0;i<
mPoints.size();i++ )
843 npoints[
mPoints[i].tower(0).index() ]++;
846 for ( UInt_t i=0;i<
mPoints.size();i++ )
852 Int_t nrel=npoints[ tower.
index() ] - 1;
855 for ( Int_t j=0;j<nn;j++ )
858 nrel+=npoints[ t2.
index() ];
861 mPoints[i].numberOfRelatives(nrel);
877 Float_t sumw[720];
for ( Int_t i=0;i<720;i++ )sumw[i]=0.;
878 Float_t sumw1[720];
for ( Int_t i=0;i<720;i++ )sumw1[i]=0.;
880 for ( UInt_t ipoint=0;ipoint<
mPoints.size();ipoint++ )
892 Int_t index=t.
index();
893 sumw[index]+=point.
energy();
901 for ( UInt_t ipoint=0;ipoint<
mPoints.size();ipoint++ )
915 Int_t index = tower.
index();
916 epoint += tower.
energy() * point.
energy() / sumw[index];
936 mPoints[ipoint].energy(epoint);
937 mPoints[ipoint].energy(epre1,1);
938 mPoints[ipoint].energy(epre2,2);
939 mPoints[ipoint].energy(epost,3);
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
StEEmcPointVec_t mSmdPoints
SMD only points.
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.