StRoot: StGammaMaker/StGammaCandidate.cxx Source File

00001 
00002 //                                                        //
00003 //  Note that, when calculating the isolation sums,       //
00004 //  the candidate momentum, tower eta and track eta       //
00005 //  have been corrected for vertex so that the            //
00006 //  dR comparisons are in fact consistent                 //
00007 //                                                        //
00009 
00010 #include "StEEmcUtil/EEmcGeom/EEmcGeomSimple.h"
00011 #include "StEEmcUtil/StEEmcSmd/EEmcSmdGeom.h"
00012 #include "StGammaEEmcLeakage.h"
00013 
00014 #include "StGammaCandidate.h"
00015 
00016 #include "TVector3.h"
00017 
00018 ClassImp(StGammaCandidate);
00019 
00021 //                 Constructor                  //
00023 StGammaCandidate::StGammaCandidate()
00024 {
00025 
00026     SetId(0);
00027     SetTowerId(0);
00028     SetTowerClusterId(0);
00029     SetSmduClusterId(0);
00030     SetSmdvClusterId(0);
00031     SetDetectorId(0);
00032     
00033     SetMomentum( TVector3(0.,0.,0.) );
00034     SetPosition( TVector3(0.,0.,0.) );
00035     SetEnergy(0.);
00036     SetSeedEnergy(0.);
00037     SetPre1Energy(0.);
00038     SetPre2Energy(0.);
00039     SetPostEnergy(0.);
00040     SetSmduEnergy(0.);
00041     SetSmdvEnergy(0.);
00042   
00043 }
00044 
00046 //                  Destructor                  //
00048 StGammaCandidate::~StGammaCandidate() 
00049 {}
00050 
00052 //    Return the total track pT and tower eT    //
00053 //   within a cone of the given radius around   //
00054 //                 the candidate                //
00056 
00057 Float_t StGammaCandidate::sumPt(Float_t radius, Float_t threshold, thresholdCut cut)
00058 {
00059     Float_t sumTracks = sumTrackPt(radius, threshold, cut);
00060     Float_t sumTowers = sumTowerPt(radius, threshold, cut);
00061     return sumTracks + sumTowers;
00062 }
00063 
00065 //    Return the total track pT within a cone   //
00066 //   of the given radius around the candidate,  //
00067 //  subject to a threshold requirement (in GeV) //
00069 
00070 Float_t StGammaCandidate::sumTrackPt(Float_t radius, Float_t threshold, thresholdCut cut)
00071 {
00072 
00073     Float_t sum = 0.0;
00074     Float_t myEta = mMomentum.Eta();
00075     Float_t myPhi = mMomentum.Phi();
00076   
00077     for(Int_t i = 0; i < numberOfTracks(); i++)
00078     {
00079     
00080         StGammaTrack *t = track(i);
00081         
00082         Float_t dEta = myEta - t->eta();
00083         Float_t dPhi = TVector2::Phi_mpi_pi(myPhi - t->phi());
00084         Float_t r = TMath::Sqrt(dEta * dEta + dPhi * dPhi);
00085         
00086         if(cut == kMagnitude)  if(r <= radius && t->momentum.Mag() > threshold) sum += t->pt();
00087         if(cut == kTransverse) if(r <= radius && t->pt() > threshold)           sum += t->pt();
00088     
00089     }
00090 
00091     return sum;
00092   
00093 }
00094 
00096 //    Return the total tower Et within a cone   //
00097 //   of the given radius around the candidate   //
00098 //  subject to a threshold requirement (in GeV) //
00100 
00101 Float_t StGammaCandidate::sumTowerPt(Float_t radius, Float_t threshold, thresholdCut cut)
00102 {
00103 
00104     Float_t sum = 0.0;
00105     Float_t myEta = mMomentum.Eta();
00106     Float_t myPhi = mMomentum.Phi();
00107   
00108     for(Int_t i = 0; i < numberOfTowers(); i++)
00109     {
00110     
00111         StGammaTower *t = tower(i);
00112         if(t->fail) continue;
00113         
00114         Float_t dEta = myEta - t->eta;
00115         Float_t dPhi = TVector2::Phi_mpi_pi(myPhi - t->phi);
00116         Float_t r = TMath::Sqrt(dEta * dEta + dPhi * dPhi);
00117         
00118         if(cut == kMagnitude)  if(r <= radius && t->energy > threshold) sum += t->pt();
00119         if(cut == kTransverse) if(r <= radius && t->pt() > threshold)   sum += t->pt();
00120     
00121     }
00122     
00123     return sum;
00124   
00125 }
00126 
00128 //  Return the total preshower Et within a cone //
00129 //    of the given radius around the candidate  //
00130 //  subject to a threshold requirement (in GeV) //
00132 
00133 Float_t StGammaCandidate::sumPreshower1(Float_t radius, Float_t threshold)
00134 {
00135 
00136     Float_t sum = 0.;
00137     Float_t myEta = mMomentum.Eta();
00138     Float_t myPhi = mMomentum.Phi();
00139     
00140     for(Int_t i = 0; i < numberOfPreshower1(); i++)
00141     {
00142     
00143         StGammaTower *t = preshower1(i);
00144         if(t->fail) continue;
00145         
00146         Float_t dEta = myEta - t->eta;
00147         Float_t dPhi = TVector2::Phi_mpi_pi(myPhi - t->phi);
00148         Float_t r = TMath::Sqrt(dEta * dEta + dPhi * dPhi);
00149         
00150         if(r <= radius && t->energy > threshold) sum += t->pt();
00151         
00152     }
00153     
00154     return sum;
00155   
00156 }
00157 
00159 //  Return the total preshower Et within a cone //
00160 //    of the given radius around the candidate  //
00161 //  subject to a threshold requirement (in GeV) //
00163 
00164 Float_t StGammaCandidate::sumPreshower2(Float_t radius, Float_t threshold)
00165 {
00166 
00167   Float_t sum = 0.0;
00168   Float_t myEta = mMomentum.Eta();
00169   Float_t myPhi = mMomentum.Phi();
00170   
00171   for(Int_t i = 0; i < numberOfPreshower2(); i++)
00172     {
00173     
00174       StGammaTower *t = preshower2(i);
00175       if(t->fail) continue;
00176       
00177       Float_t dEta = myEta - t->eta;
00178       Float_t dPhi = TVector2::Phi_mpi_pi(myPhi - t->phi);
00179       Float_t r = TMath::Sqrt(dEta * dEta + dPhi * dPhi);
00180       
00181       if(r <= radius && t->energy > threshold) sum += t->pt();
00182         
00183     }
00184     
00185   return sum;
00186   
00187 }
00188 
00190 // Return the total postshower Et within a cone //
00191 //    of the given radius around the candidate  //
00192 //  subject to a threshold requirement (in GeV) //
00194 
00195 Float_t StGammaCandidate::sumPostshower(Float_t radius, Float_t threshold)
00196 {
00197 
00198     Float_t sum = 0.0;
00199     Float_t myEta = mMomentum.Eta();
00200     Float_t myPhi = mMomentum.Phi();
00201     
00202     for(Int_t i = 0; i < numberOfPostshower(); i++)
00203     {
00204     
00205         StGammaTower *t = postshower(i);
00206         if (t->fail) continue;
00207         
00208         Float_t dEta = myEta - t->eta;
00209         Float_t dPhi = TVector2::Phi_mpi_pi(myPhi - t->phi);
00210         Float_t r = TMath::Sqrt(dEta * dEta + dPhi * dPhi);
00211         
00212         if (r <= radius && t->energy > threshold) sum += t->pt();
00213       
00214     }
00215     
00216     return sum;
00217   
00218 }
00219 
00221 //   Return the number of tracks within a cone  //
00222 //    of the given radius around the candidate  //
00223 //  subject to a threshold requirement (in GeV) //
00225 
00226 Int_t StGammaCandidate::numberOfTracks(Float_t radius, Float_t threshold, thresholdCut cut)
00227 {
00228 
00229     Int_t count = 0;
00230     Float_t myEta = mMomentum.Eta();
00231     Float_t myPhi = mMomentum.Phi();
00232     
00233     for (Int_t i = 0; i < numberOfTracks(); i++)
00234     {
00235     
00236       StGammaTrack *t = track(i);
00237       
00238       Float_t dEta = myEta - t->eta();
00239       Float_t dPhi = TVector2::Phi_mpi_pi(myPhi - t->phi());
00240       Float_t r = TMath::Sqrt(dEta * dEta + dPhi * dPhi);
00241       
00242       if(cut == kMagnitude)  if(r <= radius && t->momentum.Mag() > threshold) ++count;
00243       if(cut == kTransverse) if(r <= radius && t->pt() > threshold)           ++count;
00244       
00245     }
00246     
00247     return count;
00248   
00249 }
00250 
00252 //   Return the number of towers within a cone  //
00253 //    of the given radius around the candidate  //
00254 //  subject to a threshold requirement (in GeV) //
00256 
00257 Int_t StGammaCandidate::numberOfTowers(Float_t radius, Float_t threshold, thresholdCut cut)
00258 {
00259 
00260     Int_t count = 0;
00261     Float_t myEta = mMomentum.Eta();
00262     Float_t myPhi = mMomentum.Phi();
00263     
00264     for(Int_t i = 0; i < numberOfTowers(); i++)
00265     {
00266     
00267         StGammaTower *t = tower(i);
00268         if(t->fail) continue;
00269         
00270         Float_t dEta = myEta - t->eta;
00271         Float_t dPhi = TVector2::Phi_mpi_pi(myPhi - t->phi);
00272         Float_t r = TMath::Sqrt(dEta * dEta + dPhi * dPhi);
00273         
00274         if(cut == kMagnitude)  if(r <= radius && t->energy > threshold) ++count;
00275         if(cut == kTransverse) if(r <= radius && t->pt() > threshold)   ++count;
00276         
00277     }
00278     
00279     return count;
00280   
00281 }
00282 
00284 // Return the number of preshowers within a cone//
00285 //    of the given radius around the candidate  //
00286 //  subject to a threshold requirement (in GeV) //
00288 
00289 Int_t StGammaCandidate::numberOfPreshower1(Float_t radius, Float_t threshold)
00290 {
00291 
00292     Int_t count = 0;
00293     Float_t myEta = mMomentum.Eta();
00294     Float_t myPhi = mMomentum.Phi();
00295   
00296     for(Int_t i = 0; i < numberOfPreshower1(); i++)
00297     {
00298     
00299         StGammaTower *t = preshower1(i);
00300       
00301         Float_t dEta = myEta - t->eta;
00302         Float_t dPhi = TVector2::Phi_mpi_pi(myPhi - t->phi);
00303         Float_t r = TMath::Sqrt(dEta * dEta + dPhi * dPhi);
00304         
00305         if(r <= radius && t->energy > threshold) ++count;
00306       
00307     }
00308     
00309   return count;
00310   
00311 }
00312 
00314 // Return the number of preshowers within a cone//
00315 //    of the given radius around the candidate  //
00316 //  subject to a threshold requirement (in GeV) //
00318 
00319 Int_t StGammaCandidate::numberOfPreshower2(Float_t radius, Float_t threshold)
00320 {
00321 
00322     Int_t count = 0;
00323     Float_t myEta = mMomentum.Eta();
00324     Float_t myPhi = mMomentum.Phi();
00325     
00326     for(Int_t i = 0; i < numberOfPreshower2(); i++)
00327     {
00328     
00329         StGammaTower *t = preshower2(i);
00330       
00331         Float_t dEta = myEta - t->eta;
00332         Float_t dPhi = TVector2::Phi_mpi_pi(myPhi - t->phi);
00333         Float_t r = TMath::Sqrt(dEta * dEta + dPhi * dPhi);
00334       
00335         if(r <= radius && t->energy > threshold) ++count;
00336         
00337     }
00338     
00339     return count;
00340     
00341 }
00342 
00344 //  Return the number of postshowers within a   //
00345 // cone of the given radius around the candidate//
00346 //  subject to a threshold requirement (in GeV) //
00348 
00349 Int_t StGammaCandidate::numberOfPostshower(Float_t radius, Float_t threshold)
00350 {
00351 
00352     Int_t count = 0;
00353     Float_t myEta = mMomentum.Eta();
00354     Float_t myPhi = mMomentum.Phi();
00355     
00356     for(Int_t i = 0; i < numberOfPostshower(); i++)
00357     {
00358     
00359         StGammaTower *t = postshower(i);
00360       
00361         Float_t dEta = myEta - t->eta;
00362         Float_t dPhi = TVector2::Phi_mpi_pi(myPhi - t->phi);
00363         Float_t r = TMath::Sqrt(dEta * dEta + dPhi * dPhi);
00364         
00365         if(r <= radius && t->energy > threshold) ++count;
00366     }
00367     
00368     return count;
00369   
00370 }
00371 
00373 //       Return cluster preshower 1 energy      //
00374 // subject to a threshold requirement (in GeV)  //
00376 
00377 Float_t StGammaCandidate::pre1Energy(Float_t threshold)
00378 {
00379 
00380     mPre1Energy = 0;
00381 
00382     for(Int_t i = 0; i < numberOfMyPreshower1(); i++)
00383     {
00384     
00385         StGammaTower *t = mypreshower1(i);
00386         if(t->fail) continue;
00387         
00388         if(t->energy > threshold) mPre1Energy += t->energy;
00389         
00390     }
00391     
00392     return mPre1Energy;
00393   
00394 }
00395 
00397 //       Return cluster preshower 2 energy      //
00398 // subject to a threshold requirement (in GeV)  //
00400 
00401 Float_t StGammaCandidate::pre2Energy(Float_t threshold)
00402 {
00403 
00404     mPre2Energy = 0;
00405 
00406     for(Int_t i = 0; i < numberOfMyPreshower2(); i++)
00407     {
00408     
00409         StGammaTower *t = mypreshower2(i);
00410         if(t->fail) continue;
00411         
00412         if(t->energy > threshold) mPre2Energy += t->energy;
00413         
00414     }
00415     
00416     return mPre2Energy;
00417   
00418 }
00419 
00421 //       Return cluster postshower energy       //
00422 // subject to a threshold requirement (in GeV)  //
00424 
00425 Float_t StGammaCandidate::postEnergy(Float_t threshold)
00426 {
00427 
00428     mPostEnergy = 0;
00429 
00430     for(Int_t i = 0; i < numberOfMyPostshower(); i++)
00431     {
00432     
00433         StGammaTower *t = mypostshower(i);
00434         if(t->fail) continue;
00435         
00436         if(t->energy > threshold) mPostEnergy += t->energy;
00437         
00438     }
00439     
00440     return mPostEnergy;
00441   
00442 }
00443 
00445 //       Return cluster Smd(U/Eta) energy       //
00446 // subject to a threshold requirement (in GeV)  //
00448 
00449 Float_t StGammaCandidate::smduEnergy(Float_t threshold)
00450 {
00451 
00452     mSmduEnergy = 0;
00453 
00454     for(Int_t i = 0; i < numberOfSmdu(); i++)
00455     {
00456     
00457         StGammaStrip *t = smdu(i);
00458         if(t->fail) continue;
00459         
00460         if(t->energy > threshold) mSmduEnergy += t->energy;
00461         
00462     }
00463     
00464     return mSmduEnergy;
00465   
00466 }
00467 
00469 //       Return cluster Smd(V/Phi) energy       //
00470 // subject to a threshold requirement (in GeV)  //
00472 
00473 Float_t StGammaCandidate::smdvEnergy(Float_t threshold)
00474 {
00475 
00476     mSmdvEnergy = 0;
00477 
00478     for(Int_t i = 0; i < numberOfSmdv(); i++)
00479     {
00480     
00481         StGammaStrip *t = smdv(i);
00482         if(t->fail) continue;
00483         
00484         if(t->energy > threshold) mSmdvEnergy += t->energy;
00485         
00486     }
00487     
00488     return mSmdvEnergy;
00489   
00490 }
00491 
00493 //       Auxiliary struct and function for      //
00494 //           sorting towers by distance         //
00496 struct Tower 
00497 {
00498     StGammaTower *tower;
00499     Float_t dR;
00500 };
00501 
00502 Bool_t SortDistance(const Tower& t1, const Tower& t2)
00503 {
00504     return (t1.dR < t2.dR);
00505 }
00506 
00508 //     Recluster, computing cluster positon,    //
00509 //  momentum, and energy after exluding towers  //
00510 //      below a given threshold (in GeV)        //
00512 
00513 void StGammaCandidate::recluster(TVector3 vertex, Float_t threshold, thresholdCut cut)
00514 {
00515 
00516     mEnergy = 0.0;
00517     mPosition.SetPtEtaPhi(0, 0, 0);
00518 
00519     for(Int_t i = 0; i < numberOfMyTowers(); ++i)
00520     {
00521 
00522         StGammaTower *t = mytower(i);
00523  
00524         if(t->fail)continue;
00525 
00526         if(cut == kMagnitude)  if(t->energy < threshold) continue;
00527         if(cut == kTransverse) if(t->pt() < threshold)   continue;
00528         
00529         TVector3 tower;
00530         tower.SetPtEtaPhi(t->pt(), t->eta, t->phi);
00531 
00532         mPosition += tower * t->energy;
00533         mEnergy += t->energy;
00534 
00535     }
00536 
00537     mPosition *= 1.0 / mEnergy;
00538     
00539     mMomentum = mPosition - vertex;
00540     mMomentum.SetMag(mEnergy);
00541 
00542     return;
00543 
00544 }
00545 
00547 //   Return momentum using the seed tower only  //
00549 TVector3 StGammaCandidate::momentum1x1()
00550 {
00551     return mMomentum.Unit() * mSeedEnergy;
00552 }
00553 
00555 //     Return momentum using the seed tower     //
00556 //      with leakage correction (EEMC only)     //
00558 TVector3 StGammaCandidate::momentum1x1c()
00559 {
00560 
00561     // Perform light leakage correction w/in the EEmc
00562     if(detectorId() == kEEmc) 
00563     {
00564     
00565         static StGammaEEmcLeakage *shape = StGammaEEmcLeakage::instance();
00566         static EEmcGeomSimple    &geom  = EEmcGeomSimple::Instance();
00567         
00568         Int_t sec,sub,eta;
00569         if(!geom.getTower( position(), sec, sub, eta )) return TVector3(0.0, 0.0, -999.0);
00570         
00571         Int_t mysec,mysub,myeta;
00572         myeta = mTowerId % 12;
00573         mysec = (mTowerId / 12) / 5;
00574         mysub = (mTowerId / 12) % 5;
00575         
00576         TVector3 tower = geom.getTowerCenter( (UInt_t)sec, (UInt_t)sub, (UInt_t)eta );
00577         
00578         Float_t frac = shape->expectation( position() );
00579         TVector3 p=momentum1x1();
00580         p *= 1.0 / frac;
00581         
00582         return p;
00583     
00584     }
00585     else
00586     {
00587         // to be implemented
00588     }
00589   
00590     return momentum1x1();
00591       
00592 }
00593 
00595 //     Return momentum using the seed tower     //
00596 //     with leakage correction plus the most    //
00597 //    energetic neighboring tower(EEMC only)    //
00599 TVector3 StGammaCandidate::momentum2x1()
00600 {
00601 
00602     Float_t eta = mMomentum.Eta();
00603     Float_t phi = mMomentum.Phi();
00604     
00605     // Create Tower structs for each cluster tower
00606     std::vector<Tower> listOfTowers;
00607     for(Int_t i = 0; i < numberOfMyTowers(); i++) 
00608     {
00609     
00610         StGammaTower *t = mytower(i);
00611         Float_t deta = eta - t->eta;
00612         Float_t dphi = TVector2::Phi_mpi_pi(phi - t->phi);
00613         Float_t dR   = TMath::Sqrt(deta * deta + dphi * dphi);
00614         Tower T      = { t, dR };
00615         listOfTowers.push_back(T);
00616         
00617     }
00618     
00619     // Sort Tower structs by energy
00620     std::sort(listOfTowers.begin(), listOfTowers.end(), SortDistance);
00621     std::vector<Tower>::iterator iter;
00622     
00623     Float_t energy_sum = 0.0;
00624     Int_t count = 0;
00625         
00626     // Sum the energy from the two highest towers          
00627     for(iter=listOfTowers.begin(); iter != listOfTowers.end(); iter++)
00628     {
00629         if(count < 2) energy_sum += (*iter).tower->energy;
00630         count++;
00631     }
00632     
00633     TVector3 p = mMomentum;
00634     p.SetMag(energy_sum);
00635     
00636     return p;
00637   
00638 }
00639