00001 00017 #include "StEEmcIUSmdCluster.h" 00018 #include <TMath.h> 00019 00020 #include "StEvent/StEmcCluster.h" 00021 00022 ClassImp(StEEmcIUSmdCluster); 00023 00024 // ---------------------------------------------------------------------------- 00025 StEEmcIUSmdCluster::StEEmcIUSmdCluster() 00026 { 00027 mEnergy=0.; 00028 mSumXW=0.; 00029 mSumX2W=0.; 00030 mEmcCluster = 0; 00031 } 00032 00033 void StEEmcIUSmdCluster::print() 00034 { 00035 // std::cout << "smd cluster plane=" << plane() << std::endl; 00036 std::cout << "key=" << key() << std::endl; 00037 std::cout << "energy=" << energy() << std::endl; 00038 std::cout << "mean=" << mean() << std::endl; 00039 } 00040 00041 StEEmcIUSmdCluster::StEEmcIUSmdCluster( const StEEmcIUSmdCluster &c ) 00042 { 00043 00044 mStrips=c.mStrips; 00045 mWeights=c.mWeights; 00046 mEnergy=c.mEnergy; 00047 mSumXW=c.mSumXW; 00048 mSumX2W=c.mSumX2W; 00049 mMean=c.mMean; 00050 mSigma=c.mSigma; 00051 mSize=c.mSize; 00052 mKey=c.mKey; 00053 mMatchedTowers=c.mMatchedTowers; 00054 mEmcCluster=c.mEmcCluster; 00055 mPlane=c.mPlane; 00056 mSector=c.mSector; 00057 } 00058 // ---------------------------------------------------------------------------- 00059 void StEEmcIUSmdCluster::add( StEEmcStrip strip, Float_t weight ) 00060 { 00061 00062 if ( mStrips.size() == 0 ) 00063 { 00064 mPlane=strip.plane(); 00065 mSector=strip.sector(); 00066 } 00067 00068 mWeights.push_back(weight); 00069 mStrips.push_back(strip); 00070 00071 Float_t energy = strip.energy() * weight; 00072 00073 mEnergy += energy; 00074 mSumXW += ( strip.index() + 0.5 ) * energy; 00075 mSumX2W += ( strip.index() + 0.5 ) * ( strip.index() + 0.5 ) * energy; 00076 00077 mMean = mSumXW / mEnergy; 00078 mSigma = mSumX2W / mEnergy - mMean*mMean; 00079 mSigma = TMath::Sqrt(mSigma); 00080 00081 mSize=(Int_t)mStrips.size(); 00082 00083 } 00084 00085 // ---------------------------------------------------------------------------- 00086 StEmcCluster *StEEmcIUSmdCluster::stemc() 00087 { 00088 00089 if ( mEmcCluster ) return mEmcCluster; 00090 00091 mEmcCluster=new StEmcCluster(); 00092 Int_t iuv=( mStrips[0].plane()==1 ); // 0=U, 1=V 00093 00094 if ( iuv==0 ) { 00095 mEmcCluster->setEta( mean() ); 00096 mEmcCluster->setSigmaEta( sigma() ); 00097 mEmcCluster->setPhi( -1. ); 00098 mEmcCluster->setSigmaPhi( -1. ); 00099 } 00100 else { 00101 mEmcCluster->setEta( -1. ); 00102 mEmcCluster->setSigmaEta( -1. ); 00103 mEmcCluster->setPhi( mean() ); 00104 mEmcCluster->setSigmaPhi( sigma() ); 00105 } 00106 mEmcCluster->setEnergy( energy() ); 00107 00108 for ( Int_t i=0; i<numberOfStrips(); i++ ) { 00109 StEmcRawHit *hit=mStrips[i].stemc(); 00110 assert(hit); 00111 mEmcCluster->addHit( hit ); 00112 } 00113 00114 return mEmcCluster; 00115 00116 } 00117 00118 // ---------------------------------------------------------------------------- 00119 Float_t StEEmcIUSmdCluster::energy(Int_t nmax, Option_t *opts) 00120 { 00121 00122 Float_t mymean; 00123 if ( TString(opts).Contains("seed") ) 00124 { 00125 mymean=mStrips[0].index()+0.5; 00126 } 00127 else 00128 { 00129 mymean=mean(); 00130 } 00131 Int_t min=(Int_t)(mymean-(Float_t)nmax); 00132 Int_t max=(Int_t)(mymean+(Float_t)nmax); 00133 00134 Float_t esum=0.; 00135 00136 for ( UInt_t i=0;i<mStrips.size();i++ ) 00137 { 00138 StEEmcStrip &s=mStrips[i]; 00139 if ( s.index() >= min && s.index() <= max ) esum+=s.energy(); 00140 } 00141 00142 return esum; 00143 00144 } 00145 00146 Float_t StEEmcIUSmdCluster::sigma(Int_t nmax, Option_t *opts) 00147 { 00148 00149 Float_t mymean; 00150 if ( TString(opts).Contains("seed") ) 00151 { 00152 mymean=mStrips[0].index()+0.5; 00153 } 00154 else 00155 { 00156 mymean=mean(); 00157 } 00158 Int_t min=(Int_t)(mymean-(Float_t)nmax); 00159 Int_t max=(Int_t)(mymean+(Float_t)nmax); 00160 00161 std::cout << "mymean=" << mymean << std::endl; 00162 std::cout << "min =" << min << std::endl; 00163 std::cout << "max =" << max << std::endl; 00164 00165 Float_t esum=0.; 00166 Float_t xsum=0.; 00167 Float_t x2sum=0.; 00168 00169 for ( UInt_t i=0;i<mStrips.size();i++ ) 00170 { 00171 StEEmcStrip &s=mStrips[i]; 00172 std::cout << i << " " << s.index() << std::endl; 00173 if ( s.index() >= min && s.index() <= max ) { 00174 esum+=s.energy(); 00175 xsum+=(s.index()+0.5)*s.energy(); 00176 x2sum+=(s.index()+0.5)*(s.index()+0.5)*s.energy(); 00177 } 00178 } 00179 Float_t sig=0.; 00180 if ( esum>0. ) 00181 { 00182 sig=TMath::Sqrt( x2sum/esum - (xsum/esum)*(xsum/esum) ); 00183 } 00184 00185 00186 return sig; 00187 00188 }
1.5.9