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