StRoot: StPreEclMaker/StEmcPreCluster.cxx Source File

00001 #include "StEmcPreCluster.h"
00002 #include "StEventTypes.h"
00003 ClassImp(StEmcPreCluster)
00004 
00005 StEmcPreCluster::StEmcPreCluster(Int_t detector)
00006 {
00007     mEnergy = 0.0;
00008     mEta = 0.0;
00009     mPhi = 0.0;
00010     mSigmaEta =0.0;
00011     mSigmaPhi=0.0;
00012     mDetector=detector;
00013     mMatchingId = 0;
00014     mGeom = StEmcGeom::instance(detector);
00015 }
00016 StEmcPreCluster::StEmcPreCluster(StEmcPreCluster& cluster)
00017 {
00018     Int_t nh = cluster.nHits();
00019     for(Int_t i = 0;i<nh;i++)
00020         addHit(cluster.getHit(i));
00021     mDetector = cluster.detector();
00022     mMatchingId = cluster.matchingId();
00023     mGeom = StEmcGeom::instance(mDetector);
00024     update();
00025 }
00026 StEmcPreCluster::StEmcPreCluster(StEmcCluster* cluster)
00027 {
00028     StPtrVecEmcRawHit& hits=cluster->hit();
00029     for(Int_t i=0;i<(Int_t)hits.size();i++)
00030         addHit(hits[i]);
00031     mDetector = (Int_t)(hits[0]->detector()-kBarrelEmcTowerId+1);
00032     mMatchingId = cluster->GetUniqueID();
00033     mGeom = StEmcGeom::instance(mDetector);
00034     update();
00035 }
00036 StEmcPreCluster::~StEmcPreCluster()
00037 {
00038     mHits.Clear("nodelete"); // StEvent hits are not owned by this cluster
00039 }
00040 void StEmcPreCluster::addHit(StEmcRawHit* hit)
00041 {
00042     if(!mHits.FindObject(hit))
00043         mHits.Add(hit);
00044 }
00045 void StEmcPreCluster::removeHit(StEmcRawHit* hit)
00046 {
00047     if(mHits.FindObject(hit))
00048         mHits.Remove(hit);
00049 }
00050 void StEmcPreCluster::removeHit(Int_t hitId)
00051 {
00052     if(mHits.At(hitId))
00053         mHits.Remove(mHits.At(hitId));
00054 }
00055 void StEmcPreCluster::update()
00056 {
00057     Int_t nH = nHits();
00058     mEnergy = 0.0;
00059     mEta = 0.0;
00060     mPhi = 0.0;
00061     mSigmaEta =0.0;
00062     mSigmaPhi=0.0;
00063     if(nH==0)
00064         return;
00065 
00066     Int_t m,e,s;
00067     Float_t E,P,energy,phi0;
00068     phi0 = -999; //this is always initialized below if there are any hits at all
00069     Bool_t firstHit = true;
00070     for(Int_t i = 0;i<nH;i++)
00071     {
00072         StEmcRawHit* hit = getHit(i);
00073         if(hit)
00074         {
00075             m=(Int_t)hit->module();
00076             e=(Int_t)hit->eta();
00077             s=abs(hit->sub());
00078             energy=hit->energy();
00079 
00080             mGeom->getEta(m,e, E);
00081             mGeom->getPhi(m,s, P);
00082             // Rotate to the system of first hit
00083             if(firstHit)
00084             {
00085                 phi0 =  P;
00086                 P = 0.0;
00087                 firstHit = false;
00088             }
00089             else
00090                 P-=phi0;
00091             P = StEmcMath::getPhiPlusMinusPi(P);
00092 
00093             mEta+=E*energy;
00094             mPhi+=P*energy;
00095             mSigmaEta+=E*E*energy;
00096             mSigmaPhi+=P*P*energy;
00097             mEnergy+=energy;
00098         }
00099     }
00100     mEta /= mEnergy;
00101     mSigmaEta = mSigmaEta/mEnergy - mEta*mEta;
00102     if(mSigmaEta <= 0.0)
00103         mSigmaEta = 0.0;
00104     else
00105         mSigmaEta = sqrt(mSigmaEta);
00106 
00107     mPhi /= mEnergy;
00108     mSigmaPhi = mSigmaPhi/mEnergy - mPhi*mPhi;
00109     if(mSigmaPhi <= 0.0)
00110         mSigmaPhi = 0.0;
00111     else
00112         mSigmaPhi = sqrt(mSigmaPhi);
00113 
00114     // Rotate back to STAR system
00115     mPhi += phi0;
00116     mPhi  = StEmcMath::getPhiPlusMinusPi(mPhi);
00117 
00118     return;
00119 }
00120 void StEmcPreCluster::addCluster(StEmcPreCluster* cluster)
00121 {
00122     if(!cluster)
00123         return;
00124     if(mDetector!=cluster->detector())
00125         return;
00126     Int_t nH = cluster->nHits();
00127     for(Int_t i = 0;i<nH;i++)
00128         addHit(cluster->getHit(i));
00129     update();
00130 }
00131 void StEmcPreCluster::addCluster(StEmcCluster* cluster)
00132 {
00133     if(!cluster)
00134         return;
00135     StPtrVecEmcRawHit& hits=cluster->hit();
00136     if(mDetector!=(Int_t)(hits[0]->detector()-kBarrelEmcTowerId+1))
00137         return;
00138     for(Int_t i=0;i<(Int_t)hits.size();i++)
00139         addHit(hits[i]);
00140     update();
00141 }
00142 void StEmcPreCluster::subtractCluster(StEmcPreCluster* cluster)
00143 {
00144     if(!cluster)
00145         return;
00146     if(mDetector!=cluster->detector())
00147         return;
00148     Int_t nH = cluster->nHits();
00149     for(Int_t i = 0;i<nH;i++)
00150         removeHit(cluster->getHit(i));
00151     update();
00152 }
00153 void StEmcPreCluster::subtractCluster(StEmcCluster* cluster)
00154 {
00155     if(!cluster)
00156         return;
00157     StPtrVecEmcRawHit& hits=cluster->hit();
00158     if(mDetector!=(Int_t)(hits[0]->detector()-kBarrelEmcTowerId+1))
00159         return;
00160     for(Int_t i=0;i<(Int_t)hits.size();i++)
00161         removeHit(hits[i]);
00162     update();
00163 }
00164 StEmcPreCluster* StEmcPreCluster::splitInEta(Float_t eta)
00165 {
00166     Int_t nH = nHits();
00167     if(nH<2)
00168         return NULL; // no way to split only with 1 hit
00169     StEmcPreCluster *cluster = NULL;
00170     TList *above = new TList();
00171     Float_t E;
00172     Int_t m,e;
00173     Bool_t hasBelow = kFALSE;
00174     for(Int_t i=0;i<nH;i++)
00175     {
00176         StEmcRawHit* hit = getHit(i);
00177         if(hit)
00178         {
00179             m=(Int_t)hit->module();
00180             e=(Int_t)hit->eta();
00181             mGeom->getEta(m,e, E);
00182             if(E<=eta)
00183                 hasBelow=kTRUE;
00184             else
00185                 above->Add(hit);
00186         }
00187     }
00188     if(hasBelow && above->GetSize()>0)
00189     {
00190         cluster = new StEmcPreCluster(detector());
00191         Int_t na = above->GetSize();
00192         for(Int_t i = 0;i<na;i++)
00193         {
00194             StEmcRawHit* hit = (StEmcRawHit*)above->At(i);
00195             removeHit(hit);
00196             cluster->addHit(hit);
00197         }
00198         update();
00199         cluster->update();
00200     }
00201     delete above;
00202     return cluster;
00203 }
00204 StEmcPreCluster* StEmcPreCluster::splitInPhi(Float_t phi)
00205 {
00206     Int_t nH = nHits();
00207     if(nH<2)
00208         return NULL; // no way to split only with 1 hit
00209     StEmcPreCluster *cluster = NULL;
00210     TList *above = new TList();
00211     Float_t P;
00212     Int_t m,s;
00213     Bool_t hasBelow = kFALSE;
00214     for(Int_t i=0;i<nH;i++)
00215     {
00216         StEmcRawHit* hit = getHit(i);
00217         if(hit)
00218         {
00219             m=(Int_t)hit->module();
00220             s=abs((Int_t)hit->sub());
00221             mGeom->getPhi(m,s, P);
00222             P-= phi;
00223             P = StEmcMath::getPhiPlusMinusPi(P);
00224             if(P<=0)
00225                 hasBelow=kTRUE;
00226             else
00227                 above->Add(hit);
00228         }
00229     }
00230     if(hasBelow && above->GetSize()>0)
00231     {
00232         cluster = new StEmcPreCluster(detector());
00233         Int_t na = above->GetSize();
00234         for(Int_t i = 0;i<na;i++)
00235         {
00236             StEmcRawHit* hit = (StEmcRawHit*)above->At(i);
00237             removeHit(hit);
00238             cluster->addHit(hit);
00239         }
00240         update();
00241         cluster->update();
00242     }
00243     delete above;
00244     return cluster;
00245 }
00246 StEmcCluster* StEmcPreCluster::makeStCluster()
00247 {
00248     Int_t nH = nHits();
00249     if(nH==0)
00250         return NULL;
00251     StEmcCluster *cluster = new StEmcCluster();
00252     for(Int_t i = 0;i<nH;i++)
00253         cluster->addHit(getHit(i));
00254     update();
00255     cluster->setEta(mEta);
00256     cluster->setPhi(mPhi);
00257     cluster->setSigmaEta(mSigmaEta);
00258     cluster->setSigmaPhi(mSigmaPhi);
00259     cluster->setEnergy(mEnergy);
00260     cluster->SetUniqueID(mMatchingId);
00261     return cluster;
00262 }