50 #include "StEEmcIUClusterMaker.h"
58 #include "StEvent/StEvent.h"
59 #include "StEvent/StEmcCollection.h"
60 #include "StEvent/StEmcDetector.h"
61 #include "StEvent/StEmcModule.h"
62 #include "StEvent/StEmcClusterCollection.h"
63 #include "StEvent/StEmcCluster.h"
88 const Float_t eseeds[] = { 0.6, 1.0/1000, 1.0/1000, 1.0/1000., 0.3/1000., 0.3/1000. };
89 for ( Int_t i = 0; i < 6; i++ )
seedEnergy( eseeds[i], i );
107 StEEmcIUClusterVec_t t;
108 std::vector< StEEmcIUClusterVec_t > layers;
109 for ( Int_t i = 0; i < 4; i++ ) layers.push_back(t);
110 for ( Int_t i = 0; i < 12; i++ )
mTowerClusters.push_back(layers);
112 StEEmcIUSmdClusterVec_t s;
113 std::vector< StEEmcIUSmdClusterVec_t > planes;
116 for ( Int_t i = 0; i < 12; i++ )
mSmdClusters.push_back(planes);
117 for ( Int_t i = 0; i < 12; i++ ) TmSmdClusters.push_back(planes);
119 StEEmcTowerVec_t tow;
120 std::vector< StEEmcTowerVec_t > lay;
121 for ( Int_t i = 0; i < 4; i++ ) lay.push_back(tow);
122 for ( Int_t i = 0; i < 12; i++ )
mSeedTowers.push_back(lay);
125 mEEsmd=EEmcSmdGeom::instance();
126 mEEmap=EEmcSmdMap::instance();
137 clusize=
new TH1D(
"clustersize",
"smd cluster size ",10,0,10);
138 tclusize=
new TH1D(
"tclustersize",
"smd cluster size ",10,0,10);
139 cludis=
new TH1D(
"cludis",
"cluster distance ",20,0,20);
140 return StMaker::Init();
159 if ( !
verifyStEvent() ) Warning(
"Make",
"StEvent not properly copied");
171 for ( Int_t sector=0; sector<12; sector++ ) {
172 for ( Int_t layer=0; layer<4; layer++ )
174 for ( Int_t plane=0; plane<2; plane++ )
177 TmSmdClusters[sector][plane].clear();
181 for ( Int_t i=0;i<6;i++ )
184 TmNumberOfClusters[i]=0;
208 for ( Int_t layer=0;layer<4;layer++)
213 Float_t weights[720];
for ( Int_t i=0;i<720;i++ ) weights[i]=0.;
216 StEEmcIUClusterVec_t myClusters;
220 StEEmcTowerVec_t towers =
mEEanalysis -> towers( layer );
225 std::sort(towers.begin(),towers.end());
227 std::reverse(towers.begin(),towers.end());
232 StEEmcTowerVec_t::iterator last = towers.begin();
233 while ( last != towers.end() ) {
235 eeen+=(*last).energy();
238 ep1+=(*last).energy();
241 ep2+=(*last).energy();
244 ep3+=(*last).energy();
246 if ( (*last).energy() <
mSeedEnergy[layer] )
break;
253 std::cout <<
"-- Seed tower ----------------------" << std::endl;
260 StEEmcTowerVec_t::iterator iter = towers.begin();
261 while ( iter != towers.end() ) {
264 if ( iter==last )
break;
273 if ( weights[ (*iter).index() ] > 0. ) {
280 for ( Int_t in=0; in < (*iter).numberOfNeighbors(); in++ ) {
282 weights[ t.
index() ] += (*iter).energy();
293 while ( iter != towers.end() ) {
296 if ( iter==last )
break;
301 if ( weights[ (*iter).index() ] > 0. ) {
309 std::cout <<
"--- Clustering ----------------" << std::endl;
314 cluster.
add(seed,1.0);
316 sec=(UInt_t)seed.
sector();
318 eta=(UInt_t)seed.
etabin();
320 momentum += ( seed.
energy() * d );
328 Float_t weight = seed.
energy() / weights[ t.
index() ];
329 momentum += ( t.
energy() * d * weight );
330 cluster.
add(t, weight);
332 std::cout <<
"adding " << t.
name() <<
" E=" << t.
energy() <<
" W=" << weight << std::endl;
376 for ( Int_t sector=0; sector<12; sector++ ) {
379 for ( Int_t plane=0; plane<2; plane++ ) {
382 Float_t floor[288];
for ( Int_t i=0; i<288; i++ ) floor[i]=0.;
390 StEEmcStripVec_t seeds;
392 std::sort(strips.begin(),strips.end());
394 std::reverse(strips.begin(),strips.end());
397 StEEmcStripVec_t::iterator istrip=strips.begin();
399 while ( istrip != strips.end() ) {
400 if ( (*istrip).stat()||(*istrip).fail() ) {
404 ampe=(*(strips.begin())).energy();
408 if((*istrip).energy()>=0)
413 esmdu+=(*istrip).energy();
419 esmdv+=(*istrip).energy();
444 std::vector<Bool_t> seed_flags( strips.size(), false );
450 istrip=strips.begin();
451 while ( istrip!=strips.end() ) {
453 Int_t index=(*istrip).index();
454 Float_t eseed=(*istrip).energy();
457 if ( index <= 3 || index >= 283 ) {
465 if ( (*istrip).fail() ) {
486 seeds.push_back( (*istrip) );
487 if(plane==0){countUseed++;}
488 if(plane==1){countVseed++;}
493 float floor_para1=0.0;
494 float floor_para2=0.0;
495 if(sector==0||sector==3||sector==6||sector==9)
507 if(sector==2 || sector ==1 || sector ==5 || sector==4 || sector==8 || sector==7 || sector==11 || sector==10)
529 for ( Int_t i=0; i < 288; i++ )
532 Int_t dx=TMath::Abs(index-i);
534 if(dx<3) {floor[i]+=1.0*eseed;}
536 if(dx>=3&&dx<5){floor[i]+=floor_para1*eseed;}
538 if(dx>=5&&dx<11){floor[i]+=floor_para2*eseed;}
540 if(dx>=11&&dx<21){floor[i]+=0.05*eseed;}
543 for ( Int_t i=0; i < 288; i++ ) {
546 Int_t dx=TMath::Abs(index-i);
550 if ( 0.05 * eseed > floor[i] ) floor[i]+= 0.05 * eseed;
563 Bool_t owned[288];
for (Int_t i=0;i<288;i++) owned[i]=
false;
564 Bool_t xseed[288];
for (Int_t i=0;i<288;i++) xseed[i]=
false;
567 StEEmcStripVec_t::iterator iseed=seeds.begin();
568 while ( iseed != seeds.end() ) {
572 Int_t index=(*iseed).index();
573 if ( owned[index] || (
mSuppress&&xseed[index]) ) {
584 cluster.
add( (*iseed) );
587 Int_t ind_max = TMath::Min(287,index+max_extent);
588 Int_t ind_min = TMath::Max(0, index-max_extent);
592 for ( Int_t i=index+1; i<=ind_max; i++ ) {
599 owned[ strip.
index() ] =
true;
600 xseed[ strip.
index() ] =
true;
603 if(strip.
energy()>0) flcount++;
606 for ( Int_t i=index-1; i>=ind_min; i-- ) {
614 owned[ strip.
index() ] =
true;
615 xseed[ strip.
index() ] =
true;
618 if(strip.
energy()>0) flcount++;
625 tclusize->Fill(flcount);
627 TmSmdClusters[ sector ][ plane ].push_back(cluster);
628 TmNumberOfClusters[4+plane]++;
633 Int_t left=999,right=-999;
634 for ( Int_t is=0;is<ns;is++ )
644 if ( left-ii>=0 ) xseed[left-ii] =
true;
645 if ( right+ii<288 ) xseed[right+ii]=
true;
659 for ( Int_t sector=0; sector<12; sector++ )
661 for(Int_t plane=0;plane<=1;plane++)
664 Int_t para1[288];
for (Int_t i=0;i<288;i++) para1[i]=0;
665 for ( UInt_t fclust=0; fclust<TmSmdClusters[sector][plane].size(); fclust++ )
669 for ( Int_t fis=0;fis<fnums;fis++ )
673 int fstindex=stri1.
index();
680 Int_t Narr=TmSmdClusters[sector][plane].size();
681 int flag3[288];
for (Int_t i=0;i<Narr;i++) flag3[i]=0;
682 Int_t flag1[288];
for (Int_t i=0;i<288;i++) flag1[i]=0;
683 for ( UInt_t iclust=0; iclust<TmSmdClusters[sector][plane].size(); iclust++ )
693 for ( UInt_t jclust=0; jclust<TmSmdClusters[sector][plane].size(); jclust++ )
700 int seedds=abs(
int(oindex1-oindex2));
701 cludis->Fill(seedds);
702 if(seedds>2*max_extent)
707 if(flag3[iclust]>=2 && flag3[jclust]>=2)
713 for ( Int_t is=0;is<nums1;is++ )
716 UInt_t sindex=stri.
index();
722 for ( Int_t ist=0;ist<nums2;ist++ )
725 UInt_t tindex=strit.
index();
734 for ( Int_t is=0;is<nums1;is++ )
737 UInt_t sindex=stri.
index();
743 for ( Int_t ist=0;ist<nums2;ist++ )
746 UInt_t tindex=strit.
index();
761 for ( Int_t is=0;is<nums1;is++ )
764 UInt_t sindex=stri.
index();
765 if(para1[sindex]==2 )
767 float tpe1=stri.
energy()*EI/(EI+EIII);
773 if(para1[sindex]==1 )
780 for ( Int_t ist=0;ist<nums2;ist++ )
783 UInt_t tindex=strit.
index();
784 if(para1[tindex]==2 )
786 float tpe2=strit.
energy()*EIII/(EI+EIII);
792 if(para1[tindex]==1 )
803 for ( Int_t is=0;is<nums1;is++ )
806 UInt_t sindex=stri.
index();
807 if(para1[sindex]==2 )
809 float tpe2=stri.
energy()*EIII/(EI+EIII);
815 if(para1[sindex]==1 )
822 for ( Int_t ist=0;ist<nums2;ist++ )
825 UInt_t tindex=strit.
index();
826 if(para1[tindex]==2 )
828 float tpe1=strit.
energy()*EI/(EI+EIII);
834 if(para1[tindex]==1 )
846 if ( flag1[seedind1]<2)
850 clusize->Fill(tclust1.
size());
853 if ( flag1[seedind2]<2 )
857 clusize->Fill(tclust2.
size());
868 clusize->Fill(cl1.
size());
895 Warning(
"fillStEvent",
"called, but no StEvent to be found");
899 std::cout <<
"Adding tower clusters to StEvent at " << stevent << std::endl;
904 StEmcDetector *detector=stevent->emcCollection()->detector(kEndcapEmcTowerId);
907 Warning(
"fillStEvent",
"detector == NULL, MAJOR StEvent problem, continuing");
933 detector->setCluster( collect );
937 collect->setDetector( kEndcapEmcTowerId );
938 collect->setClusterFinderId( 123 );
939 collect->setClusterFinderParamVersion( 123 );
942 for ( Int_t isector=0; isector<12; isector++ )
946 for ( UInt_t iclust=0; iclust<
mTowerClusters[isector][0].size(); iclust++ )
955 emccluster->setEta( cl.
momentum().Eta() );
956 emccluster->setPhi( cl.
momentum().Phi() );
957 emccluster->setSigmaEta(-1.);
958 emccluster->setSigmaPhi(-1.);
959 emccluster->setEnergy( cl.
energy() );
960 emccluster->SetUniqueID( cl.
key() );
966 emccluster->addHit( hit );
970 collect->addCluster( emccluster );
974 mEtoEE[ emccluster ] = cl;
975 cl.
stemc( emccluster );
986 detector->setCluster( NULL );
999 detector=stevent->emcCollection()->detector(kEndcapEmcPreShowerId);
1002 Warning(
"fillStEvent",
"detector == NULL for pre/post, no clusters for you");
1014 detector->setCluster( pqr );
1017 pqr -> setDetector( kEndcapEmcPreShowerId );
1018 pqr -> setClusterFinderId( 123 );
1019 pqr -> setClusterFinderParamVersion( 321 );
1022 for ( Int_t isector=0; isector<12; isector++ )
1025 for ( Int_t ilayer=1; ilayer<4; ilayer++ )
1028 for ( UInt_t iclust=0; iclust<
mTowerClusters[isector][ilayer].size(); iclust++ )
1037 emccluster->setEta( cl.
momentum().Eta() );
1038 emccluster->setPhi( cl.
momentum().Phi() );
1039 emccluster->setSigmaEta(-1.);
1040 emccluster->setSigmaPhi(-1.);
1041 emccluster->setEnergy( cl.
energy() );
1042 emccluster->SetUniqueID( cl.
key() );
1048 emccluster->addHit( hit );
1053 pqr->addCluster( emccluster );
1055 mEtoEE[ emccluster ] = cl;
1056 cl.
stemc( emccluster );
1064 detector->setCluster( NULL );
1072 StDetectorId ids[]={ kEndcapSmdUStripId, kEndcapSmdVStripId };
1075 for ( Int_t iplane=0; iplane<2; iplane++ )
1078 detector=stevent->emcCollection()->detector(ids[iplane]);
1081 Warning(
"fillStEvent",
"detector == NULL for smd plane, no clusters for you");
1091 detector->setCluster( smdc );
1094 smdc->setDetector( ids[iplane] );
1095 smdc->setClusterFinderId( 123 );
1096 smdc->setClusterFinderParamVersion( 321 );
1099 for ( Int_t isector=0; isector<12; isector++ )
1102 for ( UInt_t iclust=0; iclust<
mSmdClusters[isector][iplane].size(); iclust++ )
1109 emccluster->setEta( cl.
mean() );
1110 emccluster->setPhi( (Float_t)iplane );
1111 emccluster->setSigmaEta(-1.);
1112 emccluster->setSigmaPhi(-1.);
1113 emccluster->setEnergy( cl.
energy() );
1114 emccluster->SetUniqueID( cl.
key() );
1115 printf(
"fjaoshjaojgnfaj----------------------------");
1116 printf(
"uniq id=%d\n",cl.
key());
1121 emccluster->addHit( hit );
1123 smdc->addCluster( emccluster );
1126 cl.
stemc( emccluster );
1137 detector -> setCluster( NULL );
1156 Warning(
"verifyStEvent",
"No StEvent found.");
1163 StEmcDetector *detector=stevent->emcCollection()->detector(kEndcapEmcTowerId);
1166 Warning(
"verifyStEvent",
"detector == NULL for towers");
1181 Float_t emc_sum_towers = 0.;
1182 Float_t eemc_sum_towers = 0.;
1184 StSPtrVecEmcCluster &emcClusters=cc->clusters();
1185 for ( UInt_t i=0; i<emcClusters.size(); i++ )
1189 emc_sum_towers += cl->energy();
1192 for ( Int_t sec=0; sec<12; sec++ )
1199 std::cout <<
"StEEmcIUClusterMaker tower checksum: ";
1200 if ( emc_sum_towers == eemc_sum_towers ) {
1201 std::cout <<
"passed";
1204 std::cout <<
"FAILED"; go=
false;
1206 std::cout << std::endl;
1211 std::cout <<
"StEEmcIUClusterMaker tower checksum: NULL collection, nclust=" <<
mNumberOfClusters[0] << std::endl;
1220 Float_t emc_sum_smdu = 0.;
1221 Float_t eemc_sum_smdu = 0.;
1223 detector=stevent->emcCollection()->detector(kEndcapSmdUStripId);
1226 Warning(
"verifyStEvent",
"detector == NULL for smdu");
1230 cc=detector->cluster();
1234 StSPtrVecEmcCluster &smduClusters=cc->clusters();
1236 for ( UInt_t i=0; i<smduClusters.size(); i++ )
1240 emc_sum_smdu += cl->energy();
1244 for ( Int_t sec=0; sec<12; sec++ )
1256 std::cout <<
"StEEmcIUClusterMaker smdu checksum: ";
1257 if ( emc_sum_smdu == eemc_sum_smdu ) {
1258 std::cout <<
"passed";
1261 std::cout <<
"FAILED"; go=
false;
1263 std::cout << std::endl;
1268 std::cout <<
"StEEmcIUClusterMaker smdu checksum: NULL collection, nclust=" <<
mNumberOfClusters[4] << std::endl;
1275 Float_t emc_sum_smdv = 0.;
1276 Float_t eemc_sum_smdv = 0.;
1278 detector=stevent->emcCollection()->detector(kEndcapSmdVStripId);
1281 Warning(
"verifyStEvent",
"detector == NULL for smdv");
1285 cc=detector->cluster();
1290 StSPtrVecEmcCluster &smdvClusters=cc->clusters();
1292 for ( UInt_t i=0; i<smdvClusters.size(); i++ )
1296 emc_sum_smdv += cl->energy();
1300 for ( Int_t sec=0; sec<12; sec++ )
1312 std::cout <<
"StEEmcIUClusterMaker smdv checksum: ";
1313 if ( emc_sum_smdv == eemc_sum_smdv ) {
1314 std::cout <<
"passed";
1317 std::cout <<
"FAILED"; go=
false;
1319 std::cout << std::endl;
1324 std::cout <<
"StEEmcIUClusterMaker smdv checksum: NULL collection, nclust=" <<
mNumberOfClusters[5] << std::endl;
1339 std::cout <<
"StEEmcIUClusterMaker::print()" << std::endl;
1340 const Char_t *names[] = {
"tower",
"pre1",
"pre2",
"post",
"smdu",
"smdv" };
1341 for ( Int_t i=0;i<6;i++ )
1344 std::cout <<
"Number of " << names[i]
1350 std::cout <<
"printout of tower clusters follows:" << std::endl;
1351 for ( Int_t sec=0;sec<12;sec++){
void print()
Prints cluster data.
TVector3 getTowerCenter(const UInt_t sec, const UInt_t sub, const UInt_t etabin) const
Int_t mClusterId
Keep track of clusters.
Int_t numberOfNeighbors() const
get the number of neighboring towers
EEmc ADC –> energy maker.
Int_t key()
Return a unique key assigned by the cluster maker.
std::map< StEmcCluster *, StEEmcIUCluster > mEtoEE
Map StEEmcIUClusters to StEmcClusters.
void print()
Event summary.
TVector3 momentum()
Get the momentum of this cluster.
Int_t mSuppress
Supress seeds adjacent to clusters.
Int_t numberOfStrips()
Returns the number of SMD strips in the cluster.
Int_t mMinStrips
Minimum number of smd strips to form seed.
void fillStEvent()
Fills StEvent cluster collections if the option is selected.
void neighbor(StEEmcTower *n)
add a tower to list of neighbors
StEEmcStrip & strip(Int_t sector, Int_t plane, Int_t strip)
std::map< StEmcCluster *, StEEmcIUSmdCluster > mEtoEEsmd
... and for smd clusters
void add(StEEmcTower, Float_t weight=1.0)
void name(const Char_t *n)
Set the name for this element.
Bool_t buildSmdClusters()
Constructs smd clusters.
Int_t numberOfSmdClusters(Int_t sec, Int_t plane)
Return number of smd clusters for a given sector, plane.
StEEmcIUSmdCluster smdcluster(Int_t sec, Int_t plane, Int_t index)
return a specific cluster from a given sector, plane
std::vector< std::vector< StEEmcTowerVec_t > > mSeedTowers
A base class for describing clusters of EEMC towers.
A cluster maker for the EEMC.
Int_t etabin() const
Returns the etabin of this tower, pre- or postshower element.
Bool_t buildTowerClusters()
Constructs tower clusters.
Int_t subsector() const
Returns subsector of this tower, pre- or postshower element.
StEEmcStrip strip(Int_t s)
Returns the specified smd strip w/in the cluster.
StEEmcStripVec_t & strips(Int_t sec, Int_t pln)
Returns a vector of hit strips, given the sector and plane.
void print() const
Print a summary of this tower.
Bool_t mLoose
Loose cuts option.
void index(Int_t i)
Sets the index for this SMD strip, 0..287.
Int_t key()
Returns unique id of the cluster.
Base class for representing tower, preshower and postshower elements.
StEEmcA2EMaker * mEEanalysis
ADC–>E maker.
TString mAnalysisName
ADC–>E maker name.
Int_t mMaxExtent
Maximum distance from SMD seed strips.
StEEmcIUClusterMaker(const Char_t *name="mEEclusters")
Int_t numberOfClusters(Int_t sec, Int_t layer)
Return number of clusters for a given sector, layer.
Int_t size()
Return the size (number of strips) of the cluster.
Int_t sector() const
Returns sector of this tower, pre- or postshower element.
Float_t mean()
Return the mean strip number of this cluster.
std::vector< std::vector< StEEmcIUClusterVec_t > > mTowerClusters
std::vector< std::vector< StEEmcIUSmdClusterVec_t > > mSmdClusters
StEEmcTower tower(Int_t t)
Get the specified tower.
Int_t mNumberOfClusters[6]
Counts clusters for full eemc, 0=T, 1=P, 2=Q, 3=R, 4=U, 5=V.
Int_t numberOfTowers()
Get the number of towers in cluster.
Int_t Make()
Make clusters for this event.
A base class for representing clusters of EEMC smd strips.
Float_t energy()
Get energy of this cluster.
void Clear(Option_t *opts="")
Clear clusters for next event.
void seedEnergy(Float_t energy, Int_t layer=0)
void stemc(StEmcRawHit *h)
Sets pointer to the StEmcRawHit when processing an StEvent file.
Bool_t mFillStEvent
Option to fill StEvent.
Float_t mSeedEnergy[6]
Seed energy for 0=T, 1=P, 2=Q, 3=R, 4=U, 5=V.
Bool_t mSkip
Skip strips if failbit set.
StEEmcIUCluster cluster(Int_t sec, Int_t layer, Int_t index)
Return a specific cluster from a given sector, layer.
void energy(Float_t e)
Set the energy (adc-ped+0.5)/gain for this element.
Float_t mSeedFloor
blah...
Base class for describing an endcap SMD strip.
Float_t energy()
Return the energy of this cluster.
void add(StEEmcStrip strip, Float_t weight=1.0)