00001
00027 #include "StEEmcCluster.h"
00028 #include "StEvent/StEmcCluster.h"
00029
00030 #include <iostream>
00031
00032 #include "StEEmcUtil/EEmcGeom/EEmcGeomDefs.h"
00033 #include "StEEmcUtil/EEmcGeom/EEmcGeomSimple.h"
00034
00035 ClassImp(StEEmcCluster);
00036
00037
00038 StEEmcCluster::StEEmcCluster()
00039 : StEEmcBaseCluster()
00040 {
00041 mEmcCluster=0;
00042 mEnergy=0.;
00043 mfPhibin=0.0;
00044 mfEtabin=0.0;
00045 mSumEta2W=0.;
00046 mSumEtaW=0.;
00047 mSumPhi2W=0.;
00048 mSumPhiW=0.;
00049 mMomentum=TVector3(0,0,0);
00050 }
00051
00052
00053 StEEmcCluster::StEEmcCluster( const StEEmcCluster &other )
00054 : StEEmcBaseCluster()
00055 {
00056 mEmcCluster=other.mEmcCluster;
00057 mKey=other.mKey;
00058
00059 mEnergy=other.mEnergy;
00060 mfPhibin=other.mfPhibin;
00061 mfEtabin=other.mfEtabin;
00062 mSumEta2W=other.mSumEta2W;
00063 mSumEtaW=other.mSumEtaW;
00064 mSumPhi2W=other.mSumPhi2W;
00065 mSumPhiW=other.mSumPhi2W;
00066
00067 mWeights=other.mWeights;
00068
00069 mEnergy=other.mEnergy;
00070 mMomentum=other.mMomentum;
00071 mPosition=other.mPosition;
00072 mMatched=other.mMatched;
00073
00074 for ( UInt_t ii=0;ii<other.mTowers.size();ii++ )
00075 mTowers.push_back( other.mTowers[ii] );
00076 }
00077
00078
00079 StEEmcCluster::~StEEmcCluster(){
00081
00082 mTowers.clear();
00083 mWeights.clear();
00084 }
00085
00086
00087 void StEEmcCluster::add(const StEEmcTower &tower, Float_t weight)
00088 {
00089
00090 if ( weight * tower.energy() <= 0. ) return;
00091
00092 mWeights.push_back( weight );
00093 mTowers.push_back( tower );
00094
00095 Float_t energy = weight * tower.energy();
00096 mEnergy+=energy;
00097
00098 mfEtabin+=energy*weight*tower.etabin();
00099 int myphi = (tower.phibin() - mTowers[0].phibin());
00100 mfPhibin+=energy*weight*myphi;
00101
00102 mSumEta2W = (tower.etabin()*tower.etabin())*tower.energy()*weight;
00103 mSumPhi2W = (myphi*myphi)*tower.energy()*weight;
00104 mSumEtaW = (tower.etabin())*tower.energy()*weight;
00105 mSumPhiW = (myphi)*tower.energy()*weight;
00106
00107
00108
00109 EEmcGeomSimple *geom=new EEmcGeomSimple();
00110 mMomentum += energy * ( geom->getTowerCenter( (UInt_t)tower.sector(), (UInt_t)tower.subsector(), (UInt_t)tower.etabin() ).Unit() );
00111 delete geom;
00112
00113 mPosition = mMomentum.Unit();
00114 mPosition *= ( kEEmcZSMD / mMomentum.CosTheta() );
00115
00116 if ( TMath::Abs(mPosition.Eta())>100.0 ) {
00117 Warning("add","Cluster eta is crazy!");
00118 print();
00119 }
00120
00121 mNumberOfElements=mTowers.size();
00122
00123 }
00124
00125
00126 Float_t StEEmcCluster::fracEtabin() const
00127 {
00128 return mfEtabin/mEnergy;
00129 }
00130
00131
00132 Float_t StEEmcCluster::fracPhibin() const
00133 {
00134 return mfPhibin/mEnergy+mTowers[0].phibin();
00135 }
00136
00137
00138 StEmcCluster *StEEmcCluster::stemc()
00139 {
00140
00141 if ( mEmcCluster ) return mEmcCluster;
00142 mEmcCluster=new StEmcCluster();
00143
00144 mEmcCluster->setEta( momentum().Eta() );
00145 mEmcCluster->setPhi( momentum().Phi() );
00146 mEmcCluster->setSigmaEta(-1.);
00147 mEmcCluster->setSigmaPhi(-1.);
00148 mEmcCluster->setEnergy( energy() );
00149 mEmcCluster->SetUniqueID( mKey );
00150 #if 1
00151 for ( Int_t i=0; i< numberOfTowers(); i++ )
00152 {
00153 StEmcRawHit *hit=mTowers[i].stemc();
00154 if (hit) mEmcCluster->addHit( hit );
00155 }
00156 #endif
00157
00158 return mEmcCluster;
00159 }
00160
00161
00162 void StEEmcCluster::print() const
00163 {
00164 std::cout << "cluster key: " << mKey << std::endl;
00165 std::cout << "seed tower: " << mTowers[0].name() << std::endl;
00166
00167
00168 std::cout << "fphi: " << fracPhibin() << std::endl;
00169 std::cout << "energy: " << mEnergy << std::endl;
00170 std::cout << "pt: " << mMomentum.Perp() << std::endl;
00171 std::cout << "eta: " << mMomentum.Eta() << std::endl;
00172 std::cout << "phi: " << mMomentum.Phi() << std::endl;
00173
00174 for ( UInt_t i=0;i<mTowers.size();i++ )
00175 {
00176 mTowers[i].printLine(); std::cout << " W=" << mWeights[i] << std::endl;
00177 }
00178 }
00179
00180
00181 void StEEmcCluster::printLine(Bool_t Endl) const
00182 {
00183 std::cout << "key="<<mKey<<" seed="<<mTowers[0].name()<<" ntow="<<mTowers.size()<<" pt=" <<mMomentum.Perp()<<" eta="<<mMomentum.Eta()<<" phi="<<mMomentum.Phi();
00184 if ( Endl ) std::cout << std::endl;
00185 }
00186
00187
00188 Bool_t StEEmcCluster::isNeighbor(const StEEmcTower &tower) const
00189 {
00190
00191 for ( UInt_t i=0;i<mTowers.size();i++ )
00192 {
00193 StEEmcTower myTower=mTowers[i];
00194 if ( myTower.isNeighbor(tower) ) return true;
00195 }
00196
00197 return false;
00198
00199 }
00200
00201
00202 Bool_t StEEmcCluster::hasTower(const StEEmcTower &tower) const
00203 {
00204
00205 for ( UInt_t i=0;i<mTowers.size();i++ )
00206 {
00207 if ( tower.index() == mTowers[i].index() ) return true;
00208 }
00209 return false;
00210 }
00211
00212
00213 Float_t StEEmcCluster::sigmaE() const
00214 {
00215 Float_t sumE=0.;
00216 Float_t sumE2=0.;
00217 Float_t sumw=0.;
00218 for ( UInt_t i=0;i<mTowers.size();i++ )
00219 {
00220 StEEmcTower t=mTowers[i];
00221 Float_t w=mWeights[i];
00222 sumE += t.energy()*w;
00223 sumE2 += t.energy()*t.energy()*w;
00224 sumw+=w;
00225 }
00226 Float_t mean=sumE/sumw;
00227 Float_t var=sumE2/sumw-mean*mean;
00228 return TMath::Sqrt(var);
00229 }
00230
00231
00232 Int_t StEEmcCluster::numberOfEtabins() const
00233 {
00234 Int_t etabins[12];for ( Int_t ii=0;ii<12;ii++ ) etabins[ii]=0;
00235 for ( UInt_t ii=0;ii<mTowers.size();ii++ ) etabins[ mTowers[ii].etabin() ]++;
00236 Int_t count=0;
00237 for ( Int_t ii=0;ii<12;ii++ ) if ( etabins[ii] ) count++;
00238 return count;
00239 }
00240
00241
00242 Int_t StEEmcCluster::numberOfPhibins() const
00243 {
00244 Int_t phibins[60];for ( Int_t ii=0;ii<12;ii++ ) phibins[ii]=0;
00245 for ( UInt_t ii=0;ii<mTowers.size();ii++ ) phibins[ mTowers[ii].phibin() ]++;
00246 Int_t count=0;
00247 for ( Int_t ii=0;ii<12;ii++ ) if ( phibins[ii] ) count++;
00248 return count;
00249 }
