00001
#include "StEEmcPointMaker.h"
00002
#include "StEEmcA2EMaker.h"
00003
#include "StEEmcClusterMaker.h"
00004
00005
#include "StEEmcUtil/EEmcGeom/EEmcGeomDefs.h"
00006
#include "StEEmcUtil/EEmcGeom/EEmcGeomSimple.h"
00007
#include "StEEmcUtil/StEEmcSmd/EEmcSmdGeom.h"
00008
00009
#include <iostream>
00010
#include <algorithm>
00011
#include <map>
00012
00013
#include "eeTowerFunction.h"
00014
00015
00016
#include "StEvent/StEvent.h"
00017
#include "StEvent/StEmcCollection.h"
00018
#include "StEvent/StEmcPoint.h"
00019
00020
00021
#include "TMatrixF.h"
00022
00023
00024
00025 ClassImp(
StEEmcPointMaker);
00026
00027
00028
StEEmcPointMaker::StEEmcPointMaker(
const Char_t *name):
StMaker(name)
00029 {
00030 std::cout <<
"StEEmcPointMaker("<<name<<
")" << std::endl << std::flush;
00031
00035
mEEtow=
new EEmcGeomSimple();
00036
mEEsmd=EEmcSmdGeom::instance();
00037
mEEmap=EEmcSmdMap::instance();
00038
00039
mTowerThreshold=0.;
00040
mFillStEvent=
false;
00041
mEnergyMode=1;
00042
mLimit=10;
00043
00044 }
00045
00046
00047 Int_t
StEEmcPointMaker::Init()
00048 {
00049
mEEanalysis=(
StEEmcA2EMaker*)GetMaker(mNameAnalysis);
00050
mEEclusters=(
StEEmcClusterMaker*)GetMaker(mNameClusters);
00051 assert(mEEanalysis);
00052 assert(mEEclusters);
00053
return StMaker::Init();
00054 }
00055
00056
00057 Int_t
StEEmcPointMaker::Make()
00058 {
00059
00066
00067
00068
for ( Int_t sector=0; sector<12; sector++ )
00069 {
00070
00072 StEEmcSmdClusterVec_t uclusters=
mEEclusters->
smdclusters(sector,0);
00073 StEEmcSmdClusterVec_t vclusters=
mEEclusters->
smdclusters(sector,1);
00074
00076 std::sort( uclusters.begin(), uclusters.end(), inner );
00077 std::sort( vclusters.begin(), vclusters.end(), inner );
00078
00079
findPoints( sector, uclusters, vclusters, mPoints );
00080
00081 }
00082
00083
00090
for ( UInt_t i=0;i<
mPoints.size();i++ )
00091 {
00092
mPoints[i].energy( mPoints[i].energy()/0.007/2. );
00093 }
00094
00095
00096
00097
00098 StEEmcPointVec_t orgpoints=
mPoints;
00099
00101
00102
00103
if (
mEnergyMode == 1 )
00104
shareEnergySimple();
00105
else
00106
shareEnergySmd();
00107
00108
00110
countRelatives();
00111
00112
00113
if (
mFillStEvent )
00114 {
00115
fillStEvent();
00116
verifyStEvent();
00117 }
00118
00119
00120
return kStOK;
00121 }
00122
00123
00124
00125 StEEmcPointVec_t
StEEmcPointMaker::buildSmdPoints( Int_t sector,
00126 StEEmcSmdClusterVec_t &u,
00127 StEEmcSmdClusterVec_t &v )
00128 {
00129
00130 StEEmcPointVec_t
points;
00131
00132
for ( UInt_t i=0; i<u.size(); i++ )
00133 {
00134
00135
StEEmcSmdCluster uc=u[i];
00136 Float_t xu=uc.
mean();
00137
00138
for ( UInt_t j=0;j<v.size(); j++ )
00139 {
00140
00141
00142
StEEmcSmdCluster vc=v[j];
00143 Float_t xv=vc.
mean();
00144
00146 TVector3 direct =
mEEsmd->
getIntersection( sector, xu, xv );
00147
00148 Int_t sec,sub,eta;
00149
if ( !
mEEtow->
getTower(direct,sec,sub,eta) )
00150 {
00151
continue;
00152 }
00153
else
00154 {
00156 }
00157
00160
if ( sector != sec )
00161
continue;
00162
00163
00164
00168 Bool_t good =
false;
00169
if (
mEEanalysis->
tower(sec,sub,eta).
energy() > 0. )
00170 good=
true;
00171
else if (
mEEanalysis->
tower(sec,sub,eta).
fail() )
00172 {
00173
for ( Int_t layer=1;layer<=3;layer++ )
00174 {
00175
if (
mEEanalysis->
tower(sec,sub,eta,layer).
energy() > 0. )
00176 good=
true;
00177 }
00178 }
00179
00180
if ( good ) {
00181
00182
StEEmcPoint p;
00183 p.
cluster( uc, 0 );
00184 p.
cluster( vc, 1 );
00185 p.
energy( uc.
energy() + vc.
energy() );
00186 p.
tower(
mEEanalysis->
tower(sec,sub,eta) );
00187 TVector3 position=
mEEsmd->
getIntersection(sector,uc.
mean(),vc.
mean());
00188 p.
position(position);
00189 points.push_back(p);
00190
00191 }
00192
00193 }
00194
00195 }
00196
00197
return points;
00198 }
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217 Bool_t
StEEmcPointMaker::findPoints( Int_t sector,
00218 StEEmcSmdClusterVec_t uclusters,
00219 StEEmcSmdClusterVec_t vclusters,
00220 StEEmcPointVec_t &points )
00221 {
00222
00223
00225 StEEmcPointVec_t mypoints;
00226
00228 std::sort( uclusters.begin(), uclusters.end(), inner );
00229 std::sort( vclusters.begin(), vclusters.end(), inner );
00230
00232 StEEmcPointVec_t smdpoints =
buildSmdPoints( sector, uclusters, vclusters );
00233
00234
00236
if ( smdpoints.size() < 1 )
return false;
00237
00239 std::sort( smdpoints.begin(), smdpoints.end(), chiSquare );
00240
00241
00245 std::map< Int_t, std::vector<Int_t> > u2p, v2p;
00246
for ( UInt_t i=0; i<smdpoints.size(); i++ )
00247 {
00248 u2p[ smdpoints[i].cluster(0).key() ].push_back( i );
00249 v2p[ smdpoints[i].cluster(1).key() ].push_back( i );
00250 }
00251
00256 Bool_t go =
false;
00257 StEEmcSmdClusterVec_t::iterator uiter=uclusters.begin();
00258 StEEmcSmdClusterVec_t::iterator viter=vclusters.begin();
00259
00260
00261
00262
00263
00264
00265
00266
while ( uiter<uclusters.end() || viter<vclusters.end() )
00267 {
00268
00270
StEEmcSmdCluster ucl;ucl.
key(-1);
00271
StEEmcSmdCluster vcl;vcl.
key(-1);
00272
if ( uiter<uclusters.end() ) ucl=(*uiter);
00273
if ( viter<vclusters.end() ) vcl=(*viter);
00274 Int_t iUV=-1;
00275
if ( ucl.
key()<0 )
00276 iUV=1;
00277
else if ( vcl.
key()<0 )
00278 iUV=0;
00279
else if ( (*uiter).mean() < (*viter).mean() )
00280 iUV=0;
00281
else
00282 iUV=1;
00283
00285
StEEmcSmdCluster &cluster=(iUV==0)?ucl:vcl;
00286
00288 std::vector<Int_t> matched=(iUV==0)?
00289 u2p[cluster.key()]:
00290 v2p[cluster.key()];
00291
00292
00295
if ( matched.size()==0 || matched.size() >1 )
00296 {
00297
if ( iUV==0 ) uiter++;
00298
if ( iUV==1 ) viter++;
00299
continue;
00300 }
00301
00306 go=
true;
00307
00310
StEEmcPoint p=smdpoints[matched.back()];
00311
00313 mypoints.push_back(p);
00314
00315
if ( iUV==0 ) uiter++;
00316
if ( iUV==1 ) viter++;
00317
00318 }
00319
00320
00323 Float_t chisq=9.0E9;
00324 Int_t imin=-1;
00325
for ( UInt_t i=0; i<mypoints.size(); i++ )
00326 {
00327 Float_t eu=mypoints[i].cluster(0).energy();
00328 Float_t ev=mypoints[i].cluster(1).energy();
00329 Float_t x2=(eu-ev)*(eu-ev);
00330
if ( x2 < chisq ) {
00331 imin=(Int_t)i;
00332 chisq=x2;
00333 }
00334 }
00335
00336
00343
if ( imin >= 0 ) {
00344
00345
StEEmcPoint p=mypoints[imin];
00346
removeCluster( uclusters, p.
cluster(0).key() );
00347
removeCluster( vclusters, p.
cluster(1).key() );
00348 points.push_back(p);
00349
findPoints(sector, uclusters, vclusters, points );
00350
return true;
00351
00352 }
00353
00354
00355
00356
00357
00359
00360
00367 uiter=uclusters.begin();
00368 viter=vclusters.begin();
00369
while ( uiter!=uclusters.end() || viter!=vclusters.end() )
00370 {
00371
00373
StEEmcSmdCluster ucl;ucl.
key(-1);
00374
StEEmcSmdCluster vcl;vcl.
key(-1);
00375
if ( uiter!=uclusters.end() ) ucl=(*uiter);
00376
if ( viter!=vclusters.end() ) vcl=(*viter);
00377 Int_t iUV=-1;
00378
if ( ucl.
key()<0 )
00379 iUV=1;
00380
else if ( vcl.
key()<0 )
00381 iUV=0;
00382
else if ( (*uiter).mean() < (*viter).mean() )
00383 iUV=0;
00384
else
00385 iUV=1;
00386
00388
StEEmcSmdCluster &cluster=(iUV==0)?ucl:vcl;
00389
00391 std::vector<Int_t> matched=(iUV==0)?
00392 u2p[cluster.key()]:
00393 v2p[cluster.key()];
00394
00397
if ( matched.size()==0 || matched.size()==1 )
00398 {
00399
if ( iUV==0 ) uiter++;
00400
if ( iUV==1 ) viter++;
00401
continue;
00402 }
00403
00406
StEEmcPoint p=smdpoints[matched.front()];
00407
00409 mypoints.push_back(p);
00410
00411
if ( iUV==0 ) uiter++;
00412
if ( iUV==1 ) viter++;
00413
00414 }
00415
00416
00419 chisq=9.0E9;
00420 imin=-1;
00421
for ( UInt_t i=0; i<mypoints.size(); i++ )
00422 {
00423 Float_t eu=mypoints[i].cluster(0).energy();
00424 Float_t ev=mypoints[i].cluster(1).energy();
00425 Float_t x2=(eu-ev)*(eu-ev);
00426
if ( x2 < chisq ) {
00427 imin=(Int_t)i;
00428 chisq=x2;
00429 }
00430 }
00431
00432
00433
00440
if ( imin >= 0 ) {
00441
00442
StEEmcPoint p=mypoints[imin];
00443
removeCluster( uclusters, p.
cluster(0).key() );
00444
removeCluster( vclusters, p.
cluster(1).key() );
00445 points.push_back(p);
00446
findPoints(sector, uclusters, vclusters, points );
00447
return true;
00448
00449 }
00450
00451
00452
return true;
00453
00454 }
00455
00456
00457
00458
00459
void StEEmcPointMaker::Clear( Option_t *opts )
00460 {
00461
00462
mEseen=0.;
00463
mPoints.clear();
00464
00465 }
00466
00467
00468
void StEEmcPointMaker::removeCluster( StEEmcSmdClusterVec_t &clusters, Int_t k )
00469 {
00470 StEEmcSmdClusterVec_t::iterator iter=clusters.begin();
00471
while ( iter != clusters.end() )
00472 {
00473
if ( (*iter).key() == k ) {
00474 clusters.erase(iter);
00475
return;
00476 }
00477 iter++;
00478 }
00479 }
00480
00481
00482
void StEEmcPointMaker::shareEnergy()
00483 {
00484
00490
00491
00492
00493 Int_t nrow=(Int_t)
mPoints.size();
00494 Int_t ncol=(Int_t)
mPoints.size();
00495
if ( nrow < 1 )
return;
00496
00497 std::vector<Float_t> Ef(nrow,0.);
00498
00499 TMatrixF fractions(nrow,ncol);
00500
00502
for ( Int_t k=0; k<
mEEanalysis->
numberOfHitTowers(0); k++ )
00503 {
00504
00506
StEEmcTower tower=
mEEanalysis->
hittower(k,0);
00507
00509
for ( UInt_t i=0; i<
mPoints.size(); i++ )
00510 {
00511
00514 Float_t fi=
fracp2t( mPoints[i], tower );
00515
if ( fi<=0. )
continue;
00516
00518 Ef[i] += fi * tower.
energy();
00519
00520
for ( UInt_t j=0; j<
mPoints.size(); j++ )
00521 {
00522
00525 Float_t fj=fracp2t( mPoints[j], tower );
00526
if (fi*fj<=0.)
continue;
00527
00528 fractions[i][j] += fi*fj;
00529
00530 }
00531
00532 }
00533
00534 }
00535
00536 fractions.Print();
00537
00539 Double_t det = 0.;
00540 TMatrixF invert= fractions;
00541 invert.Invert(&det);
00542
00543 invert.Print();
00544
00545 TMatrixF test= fractions * invert;
00546
00547 test.Print();
00548
00549
00551
00552 std::vector<Float_t> epoints(nrow,0.);
00553
00554
for ( Int_t i=0; i<nrow; i++ )
00555 {
00556
for ( Int_t j=0; j<ncol; j++ )
00557 {
00558 epoints[i] += invert[i][j] * Ef[j];
00559 }
00560
00561
00562 }
00563
00564
00565
00566
00567
00568
00569 }
00570
00571
00572
00573 Float_t
StEEmcPointMaker::fracp2t(
StEEmcPoint &p,
StEEmcTower &t )
00574 {
00575
00577
00578
00581
if ( !t.
isNeighbor( p.
tower(0) ) )
return 0.;
00582
00583
00584
00586 Float_t xeta=(Float_t)t.
etabin();
00587 Float_t xphi=(Float_t)t.
phibin();
00588 Double_t X[]={xphi,xeta};
00589
00591 Float_t xeta0=(Float_t)p.
tower(0).etabin();
00592 Float_t xphi0=(Float_t)p.
tower(0).phibin();
00593
00596 Int_t sec,sub,eta;
00597 Float_t dphi,deta;
00598
if ( !
mEEtow->
getTower(p.
position(), sec,sub,eta, dphi,deta ) )
return 0.;
00599 dphi/=2.0;
00600 deta/=2.0;
00601
00603 Float_t xetap=xeta0+deta;
00604 Float_t xphip=xphi0+dphi;
00605 Double_t P[]={xphip,xetap,1.0};
00606
00607
return eeTowerFunction( X, P );
00608
00609 }
00610
00611
00612
00613
void StEEmcPointMaker::shareEnergySimple()
00614 {
00615
00616 Int_t limit=
mLimit;
00617 Int_t count=0;
00618
00619
while ( count++ < limit )
00620 {
00621
00624 Float_t sumw[720];
00625
for (Int_t i=0;i<720;i++) sumw[i]=0.;
00626
00629
for ( UInt_t i=0;i<
mPoints.size();i++ )
00630 {
00631
00632
StEEmcPoint point=
mPoints[i];
00633
StEEmcTower tower=point.
tower(0);
00634
00635 sumw[ tower.
index() ] += point.
energy() *
fracp2t( point, tower );
00636
00638
for ( Int_t i=0; i<tower.
numberOfNeighbors(); i++ )
00639 {
00640
StEEmcTower neighbor=tower.
neighbor(i);
00641 sumw[ neighbor.
index() ] += point.
energy() *
fracp2t( point, neighbor );
00642 }
00643
00644 }
00645
00646
00650
00651
for ( UInt_t i=0;i<
mPoints.size();i++ )
00652 {
00653
00654
StEEmcPoint &
point=
mPoints[i];
00655
StEEmcTower tower=point.
tower(0);
00656 Float_t epoint=0.;
00657
00658 Float_t frac=0.;
00659
if ( !tower.
fail() && !tower.
stat() && sumw[tower.
index()]>0. )
00660 frac = point.
energy() *
fracp2t(point,tower) / sumw[ tower.
index() ];
00661
00662 epoint += tower.
energy() * frac;
00663
00666
for ( Int_t i=0; i<tower.
numberOfNeighbors(); i++ )
00667 {
00668
StEEmcTower neighbor=tower.
neighbor(i);
00669
if ( neighbor.
stat() || neighbor.
fail() || sumw[neighbor.
index()]<=0. )
continue;
00670 frac = point.
energy() *
fracp2t(point,neighbor) / sumw[ neighbor.
index() ];
00671 epoint += frac * neighbor.
energy();
00672 }
00673
00674
00676 point.
energy( epoint );
00677
00678 }
00679
00680 }
00681
00682
00683
00684
00685 std::vector<Bool_t> seen(720,
false);
00686
for ( UInt_t i=0;i<
mPoints.size();i++ )
00687 {
00688
00689
StEEmcTower tow=
mPoints[i].tower(0);
00690
if ( !seen[ tow.
index() ] )
mEseen+=tow.
energy();
00691 seen[ tow.
index() ] =
true;
00692
00693
for ( Int_t j=0;j<tow.
numberOfNeighbors();j++ )
00694 {
00695
StEEmcTower nei=tow.
neighbor(j);
00696
if ( !seen[ nei.
index() ] )
mEseen += nei.
energy();
00697 seen[ nei.
index() ] =
true;
00698 }
00699
00700 }
00701
00702
00703
00704
00705
00706
00707 }
00708
00709
00710
00711
00712
00713
00714
void StEEmcPointMaker::fillStEvent()
00715 {
00716
00717
00718
StEvent *stevent=(
StEvent*)GetInputDS(
"StEvent");
00719
if ( !stevent ) {
00720 Warning(
"fillStEvent",
"called, but no StEvent to be found");
00721
return;
00722 }
00723
00724
00726
for ( UInt_t i=0; i<
mPoints.size(); i++ )
00727 {
00728
00729
StEmcPoint *
point=
mPoints[i].stemc();
00730 stevent->
emcCollection()->
addEndcapPoint( point );
00731
00732
mEtoEE[ point ] = mPoints[i];
00733
00734
00735 }
00736
00737 }
00738
00739
00740
00741
00742
void StEEmcPointMaker::verifyStEvent()
00743 {
00744
00745 Float_t emc_sum_points = 0.;
00746 Float_t eemc_sum_points = 0.;
00747
00748
StEvent *stevent=(
StEvent*)GetInputDS(
"StEvent");
00749
if ( !stevent ) {
00750 Warning(
"verifyStEvent",
"called, but no StEvent to be found");
00751
return;
00752 }
00753
00754 StSPtrVecEmcPoint& emcpts = stevent->
emcCollection()->
endcapPoints();
00755
for ( UInt_t i=0;i<emcpts.size();i++ )
00756 {
00757
00758
StEmcPoint *p=emcpts[i];
00759 assert(p);
00760 emc_sum_points += p->
energy();
00761
00762 }
00763
00764
for ( UInt_t i=0;i<
mPoints.size();i++ )
00765 {
00766
00767
StEEmcPoint p=
mPoints[i];
00768 eemc_sum_points += p.
energy();
00769
00770 }
00771
00772 std::cout <<
"StEEmcPointMaker point checksum: ";
00773
if ( emc_sum_points == eemc_sum_points )
00774 {
00775 std::cout <<
"passed";
00776 }
00777
else
00778 std::cout <<
"FAILED";
00779 std::cout << std::endl;
00780
00781
00782 }
00783
00784
00785
00786
void StEEmcPointMaker::countRelatives()
00787 {
00788
00790 Int_t npoints[720];
00791
for ( Int_t i=0;i<720;i++ ) npoints[i]=0;
00792
00793
for ( UInt_t i=0;i<
mPoints.size();i++ )
00794 npoints[
mPoints[i].tower(0).index() ]++;
00795
00797
for ( UInt_t i=0;i<mPoints.size();i++ )
00798 {
00799
00800
StEEmcTower tower=mPoints[i].tower(0);
00801 Int_t nn=tower.
numberOfNeighbors();
00802
00803 Int_t nrel=npoints[ tower.
index() ] - 1;
00804 assert(nrel>=0);
00805
00806
for ( Int_t j=0;j<nn;j++ )
00807 {
00808
StEEmcTower t2=tower.
neighbor(j);
00809 nrel+=npoints[ t2.
index() ];
00810 }
00811
00812 mPoints[i].numberOfRelatives(nrel);
00813
00814 }
00815
00816
00817 }
00818
00819
00820
00821
00822
00823
void StEEmcPointMaker::shareEnergySmd()
00824 {
00825
00828 Float_t sumw[720];
for ( Int_t i=0;i<720;i++ )sumw[i]=0.;
00829
00830
for ( UInt_t ipoint=0;ipoint<
mPoints.size();ipoint++ )
00831 {
00832
00833
00834
StEEmcPoint point=
mPoints[ipoint];
00835
StEEmcTower tower=point.
tower(0);
00836 sumw[tower.
index()]+=point.
energy();
00837
00838
for ( Int_t itow=0;itow<tower.
numberOfNeighbors();itow++ )
00839 {
00840
StEEmcTower t=tower.
neighbor(itow);
00841 Int_t index=t.
index();
00842 sumw[index]+=point.
energy();
00843 }
00844
00845 }
00846
00849
00850
for ( UInt_t ipoint=0;ipoint<
mPoints.size();ipoint++ )
00851 {
00852
00853
StEEmcPoint point=
mPoints[ipoint];
00854
StEEmcTower tower = point.
tower(0);
00855 Float_t epoint = 0.;
00856 Int_t index = tower.
index();
00857 epoint += tower.
energy() * point.
energy() / sumw[index];
00858
00859
for ( Int_t itow=0;itow<tower.
numberOfNeighbors();itow++ )
00860 {
00861
StEEmcTower t=tower.
neighbor(itow);
00862 index = t.
index();
00863 epoint += t.
energy() * point.
energy() / sumw[index];
00864 }
00865
00866
00867 mPoints[ipoint].
energy(epoint);
00868
00869 }
00870
00871
00872 }
00873
00874