StRoot: StGammaMaker/StGammaCandidateMaker.cxx Source File

00001 #include "TVector2.h"
00002 
00003 #include "StEEmcPool/StEEmcA2EMaker/StEEmcA2EMaker.h"
00004 #include "StEEmcPool/StEEmcClusterMaker/StEEmcGenericClusterMaker.h"
00005 
00006 #include "SystemOfUnits.h"
00007 #include "StEmcUtil/geometry/StEmcGeom.h"
00008 #include "StEmcUtil/projection/StEmcPosition.h"
00009 
00010 #include "StBarrelEmcCluster.h"
00011 #include "StBarrelEmcClusterMaker.h"
00012 #include "StGammaFitter.h"
00013 #include "StGammaEvent.h"
00014 #include "StGammaEventMaker.h"
00015 #include "StGammaRawMaker.h"
00016 
00017 #include "StGammaCandidateMaker.h"
00018 
00019 ClassImp(StGammaCandidate);
00020 
00022 //                 Constructor                  //
00024 StGammaCandidateMaker::StGammaCandidateMaker(const char *name): StMaker(name)
00025 {
00026 
00027     mUseBemc = false;
00028     mUseEemc = false;
00029     
00030     mStrictBemcStatus = true;
00031     
00032     mCompressLevel = 2;
00033 
00034     mMinimumEt = 0.0;
00035     mRadius = 0.7; 
00036     mBsmdRange = 0.05237;
00037     mEsmdRange = 20.0;   
00038     
00039 }
00040 
00042 //                  Destructor                  //
00044 StGammaCandidateMaker::~StGammaCandidateMaker()
00045 {}
00046 
00048 //                 Maker Init                   //
00050 Int_t StGammaCandidateMaker::Init()
00051 {
00052     return StMaker::Init();
00053 }
00054 
00056 //                 Maker Clear                  //
00058 void StGammaCandidateMaker::Clear( Option_t *opts )
00059 {
00060     mId = 0;
00061 }
00062 
00064 //                  Maker Make                  //
00066 Int_t StGammaCandidateMaker::Make()
00067 { 
00068 
00069     if(mUseBemc) MakeBarrel();
00070     if(mUseEemc) MakeEndcap();
00071     
00072     Compress();
00073     
00074     return kStOK;
00075   
00076 }
00077 
00078 
00080 //     Return the vertex corrected momentum     //
00081 //         for the given EEMC cluster           //
00083 TVector3 momentum(StEEmcCluster cluster, TVector3 vertex)
00084 {
00085 
00086     TVector3 d = cluster.momentum().Unit();
00087     
00088     d.SetMag( kEEmcZSMD / d.CosTheta() );
00089     d = d - vertex;
00090     d = d.Unit();
00091     d *= cluster.energy();
00092     
00093     return d;
00094   
00095 }
00096 
00098 //            Compile BEMC Clusters             //
00100 Int_t StGammaCandidateMaker::MakeBarrel()
00101 {
00102 
00103     // Retrieve all those gamma makers from the chain
00104     StGammaEventMaker *mGammaEventMaker = dynamic_cast<StGammaEventMaker*>(GetMakerInheritsFrom("StGammaEventMaker"));
00105     if(!mGammaEventMaker)
00106     {
00107         LOG_WARN << "MakeBarrel() - No StGammaEventMaker found!" << endm;
00108         return kStWarn;
00109     }
00110 
00111     StGammaEvent *mGammaEvent = mGammaEventMaker->event();
00112     if(!mGammaEvent) 
00113     {
00114         LOG_WARN << "MakeBarrel() - No StGammaEvent found!" << endm;
00115         return kStWarn;
00116     }
00117 
00118     StGammaRawMaker *mGammaRawMaker = dynamic_cast<StGammaRawMaker*>(GetMakerInheritsFrom("StGammaRawMaker"));
00119     if(!mGammaRawMaker)
00120     {
00121         LOG_WARN << "MakeBarrel() - No StGammaRawMaker found!" << endm;
00122         return kStWarn;
00123     }
00124     
00125     StBarrelEmcClusterMaker *mBemcClusterMaker = dynamic_cast<StBarrelEmcClusterMaker*>(GetMakerInheritsFrom("StBarrelEmcClusterMaker"));
00126     if(!mBemcClusterMaker)
00127     {
00128         LOG_WARN << "MakeBarrel() - No StBarrelEmcClusterMaker found!" << endm;
00129         return kStWarn;
00130     }
00131 
00132     // Create a StEmcGeom instance for track projections
00133     StEmcGeom* geom = StEmcGeom::instance("bemc");
00134 
00135     StEmcGeom *smdEtaGeom = new StEmcGeom("bsmde");
00136     StEmcGeom *smdPhiGeom = new StEmcGeom("bsmdp");
00137 
00138     StEmcPosition emcPosition;
00139 
00140     // Loop over BEMC clusters
00141     for(unsigned int i = 0; i < mBemcClusterMaker->clusters().size(); ++i) 
00142     {
00143     
00144         StBarrelEmcCluster* cluster = mBemcClusterMaker->clusters()[i];
00145         
00146         // Require minimum Et
00147         if(cluster->momentum().Pt() < mMinimumEt) continue;
00148         
00149         // Require that seed towers have non fail status
00150         if(cluster->seed()->fail) continue;
00151         
00152         // If strict status is observed, require that all
00153         // neighboring towers also have a non fail status
00154         if(mStrictBemcStatus)
00155         {
00156         
00157             bool fail = false;
00158         
00159             for(int deta = -1; deta <= 1; ++deta)
00160             {
00161                 for(int dphi = -1; dphi <= 1; ++dphi)
00162                 {
00163                     if(StGammaTower *tower = cluster->tower(deta, dphi))
00164                     {
00165                         if(tower->fail) fail = true;
00166                     }
00167                 }
00168             }
00169             
00170             if(fail) continue;
00171         
00172         }
00173         
00174         // Create gamma candidate
00175         StGammaCandidate *candidate = mGammaEvent->newCandidate();
00176         candidate->SetTowerClusterId(i);
00177         candidate->SetId(nextId());
00178         candidate->SetDetectorId(StGammaCandidate::kBEmc);
00179         candidate->SetEnergy(cluster->energy());
00180         candidate->SetPosition(cluster->position());
00181         candidate->SetMomentum(cluster->momentum());
00182         
00183         // Start with the seed tower
00184         StGammaTower *seed = cluster->seed();
00185         candidate->SetSeedEnergy(seed->energy);
00186         candidate->SetTowerId(seed->id);
00187 
00188         // Loop over neighboring towers in the cluster
00189         float energy = 0;
00190         
00191         for(int deta = -1; deta <= 1; ++deta) 
00192         {
00193         
00194             for(int dphi = -1; dphi <= 1; ++dphi) 
00195             {
00196             
00197                 if(StGammaTower *tower = cluster->tower(deta, dphi)) 
00198                 {
00199                 
00200                     // Skip any failing towers
00201                     if(tower->fail) continue;
00202                     
00203                     // Associate tower with the candidate
00204                     candidate->addMyTower(tower);
00205                     tower->candidates.Add(candidate);
00206                     
00207                     // Associate preshower with the candidate
00208                     if(StGammaTower *preshower = mGammaRawMaker->tower(tower->id, kBEmcPres)) 
00209                     {
00210                         candidate->addMyPreshower1(preshower);
00211                         preshower->candidates.Add(candidate);
00212                         energy += preshower->energy;
00213                     }
00214                     
00215                     // Associate tracks with the candidate
00216                     for(int k = 0; k < mGammaEvent->numberOfTracks(); ++k) 
00217                     {
00218                     
00219                         StGammaTrack* track = mGammaEvent->track(k);
00220                         if(!track) continue;
00221                     
00222                         TVector3 position = track->positionAtRadius(geom->Radius());
00223                         if(position != TVector3(0, 0, 0))
00224                         {
00225                             
00226                             int id;
00227                             if(geom->getId(position.Phi(), position.Eta(), id) == 0 && id == static_cast<int>(tower->id)) 
00228                             {
00229                                 candidate->addMyTrack(track);
00230                                 track->candidates.Add(candidate);
00231                             }
00232                     
00233                         }
00234                     
00235                     }
00236 
00237                 } // if(tower)
00238             
00239             } // dphi
00240 
00241         } // deta
00242 
00243 
00244         // Set candidate preshower sum
00245         candidate->SetBPrsEnergy(energy);
00246         
00247         // Associate tracks falling within a cone of mRadius around the candidate
00248         for (int k = 0; k < mGammaEvent->numberOfTracks(); ++k) 
00249         {
00250     
00251             if (StGammaTrack* track = mGammaEvent->track(k)) 
00252             {
00253             
00254                 float deta = cluster->momentum().Eta() - track->eta();
00255                 float dphi = cluster->momentum().Phi() - track->phi();
00256                 dphi = TVector2::Phi_mpi_pi(dphi);
00257                 float r = hypot(deta, dphi);
00258                 
00259                 if(r <= mRadius) 
00260                 {
00261                     candidate->addTrack(track);
00262                     if(!track->candidates.FindObject(candidate)) track->candidates.Add(candidate);
00263                 }
00264                 
00265             }
00266             
00267         }
00268 
00269         // Associate towers falling within a cone of radius mRadius around the candidate
00270         for(int k = 0; k < mGammaEvent->numberOfTowers(); ++k) 
00271         {
00272         
00273             StGammaTower* tower = mGammaEvent->tower(k);
00274             if(tower && tower->layer == kBEmcTower) 
00275             {
00276             
00277                 float deta = cluster->momentum().Eta() - tower->eta;
00278                 float dphi = cluster->momentum().Phi() - tower->phi;
00279                 dphi = TVector2::Phi_mpi_pi(dphi);
00280                 float r = hypot(deta, dphi);
00281                 
00282                 if(r <= mRadius) 
00283                 {
00284                     candidate->addTower(tower);
00285                     if(!tower->candidates.FindObject(candidate)) tower->candidates.Add(candidate);
00286                     
00287                     // Associate tracks projecting to the tower with the tower
00288                     for(int t = 0; t < mGammaEvent->numberOfTracks(); ++t) 
00289                     {
00290                     
00291                         StGammaTrack* track = mGammaEvent->track(t);
00292                         if(!track) continue;
00293                         
00294                         TVector3 position = track->positionAtRadius(geom->Radius());
00295                         if(position != TVector3(0, 0, 0))
00296                         {
00297                             
00298                             int id;
00299                             if(geom->getId(position.Phi(), position.Eta(), id) == 0 && id == static_cast<int>(tower->id)) 
00300                             {
00301                                 if(!tower->tracks.FindObject(track)) tower->tracks.Add(track);
00302                             }
00303                     
00304                         }
00305                     
00306                     }
00307                     
00308                 }
00309                 
00310             }
00311             
00312         }
00313 
00314         // Associate preshower hits falling within a cone of radius mRadius around the candidate
00315         for(int k = 0; k < mGammaEvent->numberOfPreshower1(); ++k) 
00316         {
00317         
00318             StGammaTower* preshower = mGammaEvent->preshower1(k);
00319             if(preshower && preshower->layer == kBEmcPres) 
00320             {
00321             
00322                 float deta = cluster->momentum().Eta() - preshower->eta;
00323                 float dphi = cluster->momentum().Phi() - preshower->phi;
00324                 dphi = TVector2::Phi_mpi_pi(dphi);
00325                 float r = hypot(deta, dphi);
00326                 
00327                 if (r <= mRadius) 
00328                 {
00329                     candidate->addPreshower1(preshower);
00330                     if(!preshower->candidates.FindObject(candidate)) preshower->candidates.Add(candidate);
00331                 }
00332                 
00333             }
00334             
00335         }
00336 
00337         // Associate BSMDE and BSMDP strips falling within mSmdRange of the candidate
00338         float smdEtaEnergy = 0;
00339         float smdPhiEnergy = 0;
00340         
00341         int centralId = 0;
00342         int module = 0;
00343         int eta = 0;
00344         int sub = 0;
00345 
00346         smdEtaGeom->getId(cluster->position().Phi(), cluster->position().Eta(), centralId);         
00347 
00348         int sector = 0; // Dummy variable in the BEMC
00349 
00350         for(int dPhiStrip = -1; dPhiStrip <= 1; ++dPhiStrip)
00351         {
00352 
00353             for(int dEtaStrip = -20; dEtaStrip <= 20; ++dEtaStrip)
00354             {
00355 
00356                 int id = emcPosition.getNextId(3, centralId, dEtaStrip, dPhiStrip);
00357                 if(id == 0) continue;        
00358 
00359                 smdEtaGeom->getBin(id, module, eta, sub);
00360 
00361                 // Find vertex corrected position of the strip
00362                 float x, y, z;
00363                 smdEtaGeom->getXYZ(id, x, y, z);
00364                 TVector3 vEta(x, y, z);
00365 
00366                 float deta = cluster->position().Eta() - vEta.Eta();
00367                 float dphi = cluster->position().Phi() - vEta.Phi();
00368                 dphi = TVector2::Phi_mpi_pi(dphi);
00369                 float r = hypot(deta, dphi);
00370             
00371                 // If the strip is within range then associate the strip
00372                 if(r <= mBsmdRange)
00373                 {                      
00374                 
00375                     // Use strip is already created in StGammaRawMaker
00376                     if(StGammaStrip* strip = mGammaRawMaker->strip(sector, kBEmcSmdEta, id)) 
00377                     {
00378                 
00379                         candidate->addSmdEta(strip);
00380                         strip->candidates.Add(candidate);
00381                         smdEtaEnergy += strip->energy;
00382                 
00383                     }
00384                     // Otherwise create a new, empty strip
00385                     else
00386                     {
00387 
00388                         float eta;
00389 
00390                         smdEtaGeom->getEta(id, eta);
00391              
00392                         StGammaStrip *strip = mGammaEvent->newStrip();
00393          
00394                         strip->index = id;
00395                         strip->sector = module;
00396                         strip->plane  = kBEmcSmdEta;
00397                         // strip->stat Filled in StGammaRawMaker::AddEtaStrip()
00398                         // strip->fail Filled in StGammaRawMaker::AddEtaStrip()
00399                         strip->energy = 0;
00400                         strip->position = 230.705 * sinh(eta);
00401          
00402                         double offset = fabs(eta) < 0.5 ? 0.73 : 0.94;
00403               
00404                         strip->left = strip->position - offset;
00405                         strip->right = strip->position + offset;
00406 
00407                         mGammaRawMaker->AddEtaStrip(strip);
00408             
00409                         candidate->addSmdEta(strip);
00410                         strip->candidates.Add(candidate);
00411 
00412                     } // if exist
00413 
00414                 } // if range
00415 
00416             } // dEtaStrips
00417 
00418         } // dPhiStrips
00419 
00420         centralId = 0;
00421         smdPhiGeom->getId(cluster->position().Phi(), cluster->position().Eta(), centralId);
00422                         
00423         for(int dEtaStrip = -1; dEtaStrip <= 1; ++dEtaStrip)
00424         {
00425                         
00426             for(int dPhiStrip = -20; dPhiStrip <= 20; ++dPhiStrip)
00427             {
00428             
00429                 int id = emcPosition.getNextId(4, centralId, dEtaStrip, dPhiStrip);
00430                 if(id == 0) continue;
00431 
00432                 smdEtaGeom->getBin(id, module, eta, sub);
00433 
00434                 // Find vertex corrected position of the strip
00435                 float x, y, z;
00436                 smdPhiGeom->getXYZ(id, x, y, z);
00437                 TVector3 vPhi(x, y, z);
00438 
00439                 float deta = cluster->position().Eta() - vPhi.Eta();
00440                 float dphi = cluster->position().Phi() - vPhi.Phi();
00441                 dphi = TVector2::Phi_mpi_pi(dphi);
00442                 float r = hypot(deta, dphi);
00443 
00444                 // If the strip is within range then associate the strip
00445                 if(r < mBsmdRange)
00446                 {
00447 
00448                     // Use strip is already created in StGammaRawMaker
00449                     if(StGammaStrip* strip = mGammaRawMaker->strip(sector, kBEmcSmdPhi, id)) 
00450                     {
00451             
00452                         candidate->addSmdPhi(strip);
00453                         strip->candidates.Add(candidate);
00454                         smdPhiEnergy += strip->energy;
00455                   
00456                     }
00457                     // Otherwise create a new, empty strip
00458                     else
00459                     {
00460 
00461                         float phi;
00462 
00463                         smdPhiGeom->getPhi(id, phi);
00464 
00465                         StGammaStrip *strip = mGammaEvent->newStrip();
00466 
00467                         strip->index = id;
00468                         strip->sector = module;
00469                         strip->plane  = kBEmcSmdPhi;
00470                         //strip->stat Filled in StGammaRawMaker::AddPhiStrip()
00471                         //strip->fail Filled in StGammaRawMaker::AddPhiStrip()
00472                         strip->energy = 0;
00473                         strip->position = phi;
00474 
00475                         double offset = 0.00293;
00476 
00477                         strip->left = phi - offset;
00478                         strip->right = phi + offset;
00479 
00480                         mGammaRawMaker->AddPhiStrip(strip);
00481 
00482                         candidate->addSmdPhi(strip);
00483                         strip->candidates.Add(candidate);
00484 
00485                     } // if exist
00486                   
00487                 } // if range
00488 
00489             } // dEtaStrips
00490 
00491         } // dPhiStrips
00492 
00493         candidate->SetSmdEtaEnergy(smdEtaEnergy);
00494         candidate->SetSmdPhiEnergy(smdPhiEnergy);
00495         
00496     } // Clusters
00497     
00498     return kStOK;
00499     
00500 }
00501 
00503 //            Compile EEMC Clusters             //
00505 Int_t StGammaCandidateMaker::MakeEndcap()
00506 {
00507 
00508     // Retrieve all those gamma makers from the chain
00509     StGammaEventMaker *mGammaEventMaker = dynamic_cast<StGammaEventMaker*>(GetMakerInheritsFrom("StGammaEventMaker"));
00510     if(!mGammaEventMaker)
00511     {
00512         LOG_WARN << "MakeEndcap() - No StGammaEventMaker found!" << endm;
00513         return kStWarn;
00514     }
00515 
00516     StGammaEvent *mGammaEvent = mGammaEventMaker->event();
00517     if(!mGammaEvent) 
00518     {
00519         LOG_WARN << "MakeEndcap() - No StGammaEvent found!" << endm;
00520         return kStWarn;
00521     }
00522 
00523     StGammaRawMaker *mGammaRawMaker = dynamic_cast<StGammaRawMaker*>(GetMakerInheritsFrom("StGammaRawMaker"));
00524     if(!mGammaRawMaker)
00525     {
00526         LOG_WARN << "MakeEndcap() - No StGammaRawMaker found!" << endm;
00527         return kStWarn;
00528     }
00529     
00530     StEEmcGenericClusterMaker *mEEclusters = dynamic_cast<StEEmcGenericClusterMaker*>(GetMakerInheritsFrom("StEEmcGenericClusterMaker"));
00531     if(!mEEclusters)
00532     {
00533         LOG_DEBUG << "MakeEndcap() - No StEEmcGenericClusterMaker found!" << endm;
00534         return kStWarn;
00535     }
00536 
00537     StEEmcA2EMaker *mEEanalysis = dynamic_cast<StEEmcA2EMaker*>(GetMakerInheritsFrom("StEEmcA2EMaker"));
00538     if(!mEEanalysis)
00539     {
00540         LOG_DEBUG << "MakeEndcap() - No StEEmcA2EMaker found!" << endm;
00541         return kStWarn;
00542     }
00543 
00544 
00545     // Loop over each sector, 
00546     for(Int_t sector = 0; sector < kEEmcNumSectors; sector++)
00547     {
00548 
00549         // then the clusters in each sector
00550         StEEmcClusterVec_t clusters = mEEclusters->clusters(sector, 0);
00551         for(UInt_t i = 0; i < clusters.size(); i++)
00552         {
00553 
00554             // Correct the cluster momentum for vertex
00555             StEEmcCluster cluster = clusters[i];
00556             TVector3 position = getEEmcClusterPosition(cluster);
00557             if(position == TVector3(0,0,0)) position = cluster.position();
00558             TVector3 pcluster = position - mGammaEvent->vertex();
00559             pcluster.SetMag(cluster.energy());
00560             
00561             // Require minimum transverse energy
00562             Float_t ET = pcluster.Perp();
00563             if(ET < mMinimumEt) continue;
00564             
00565             // Creat a new candidate
00566             StGammaCandidate *mCandidate = mGammaEvent->newCandidate();
00567             
00568             mCandidate->SetId( nextId() );
00569             mCandidate->SetTowerClusterId( cluster.key() );
00570             mCandidate->SetDetectorId( StGammaCandidate::kEEmc );
00571             
00572             // Set candidate kinematics
00573             mCandidate->SetMomentum( pcluster );
00574             mCandidate->SetPosition( cluster.position() );
00575             mCandidate->SetEnergy( cluster.energy() );
00576             
00577             // Set seed information
00578             StEEmcTower seed = cluster.tower(0);
00579             mCandidate->SetSeedEnergy( seed.energy() );
00580             mCandidate->SetTowerId( seed.index() );
00581           
00582 
00583             // Associate towers, preshowers, and postshowers with the candidate
00584             // while summing up the total energy response in each
00585             Float_t epre1 = 0.;
00586             Float_t epre2 = 0.;
00587             Float_t epost = 0.;
00588             
00589             // Loop over neighboring towers
00590             for(Int_t j = 0; j < cluster.numberOfTowers(); j++)
00591             {
00592             
00593                 // Grab the jth neighboring tower
00594                 StEEmcTower tower = cluster.tower(j);
00595                 Int_t index = tower.index();
00596                 
00597                 // and its pre/postshowers
00598                 StEEmcTower pre1 = mEEanalysis->tower(index, 1);
00599                 StEEmcTower pre2 = mEEanalysis->tower(index, 2);
00600                 StEEmcTower post = mEEanalysis->tower(index, 3);
00601                 
00602                 // Increment the energy sums
00603                 if( !pre1.fail() ) epre1 += pre1.energy();
00604                 if( !pre2.fail() ) epre2 += pre2.energy();
00605                 if( !post.fail() ) epost += post.energy();
00606 
00607                 // Remember, if the respective detector element did not have
00608                 // sufficient energy deposited in it then these pointers
00609                 // will be null
00610                 StGammaTower *gtower = mGammaRawMaker->tower(index, kEEmcTower);
00611                 StGammaTower *gpre1  = mGammaRawMaker->tower(index, kEEmcPre1);
00612                 StGammaTower *gpre2  = mGammaRawMaker->tower(index, kEEmcPre2);
00613                 StGammaTower *gpost  = mGammaRawMaker->tower(index, kEEmcPost);
00614                 
00615                 
00616                 // Associate clustered towers with the candidate
00617                 if(gtower)
00618                 {
00619                     mCandidate->addMyTower(gtower);
00620                     gtower->candidates.Add(mCandidate);
00621                 }
00622             
00623                 // Associate clustered preshower1 tiles with the candidate
00624                 if(gpre1 && !pre1.fail() && pre1.energy() > 0.0) 
00625                 {
00626                     mCandidate->addMyPreshower1(gpre1);
00627                     gpre1->candidates.Add(mCandidate);
00628                 }
00629                 
00630                 // Associate clustered preshower2 tiles with the candidate
00631                 if(gpre2 && !pre2.fail() && pre2.energy() > 0.0) 
00632                 {
00633                     mCandidate->addMyPreshower2(gpre2);
00634                     gpre2->candidates.Add(mCandidate);
00635                 }
00636                 
00637                 // Associate clustered postshower tiles with the candidate
00638                 if(gpost && !post.fail() && post.energy() > 0.0) 
00639                 {
00640                     mCandidate->addMyPostshower(gpost);
00641                     gpost->candidates.Add(mCandidate);
00642                 }
00643 
00644                 // Associate tracks extrapolating to the tower to the cluster
00645                 if(gtower) 
00646                 {
00647                 
00648                     // Loop over all tracks in the event
00649                     for(int k = 0; k < mGammaEvent->numberOfTracks(); ++k) 
00650                     {
00651                     
00652                         StGammaTrack* track = mGammaEvent->track(k);
00653                         
00654                         // Skip misbehaving or null tracks
00655                         if (!track || track->pz() < 0) continue;
00656                         
00657                         // Extrapolate track to the SMD distance
00658                         EEmcGeomSimple& geom = EEmcGeomSimple::Instance();
00659                         TVector3 position = track->positionAtZ(geom.getZ1());
00660                     
00661                         if(position != TVector3(0, 0, 0))
00662                         {
00663                     
00664                             // If the tower to which the track extrapolates is equal
00665                             // to the current tower then associate track with cluster
00666                             int sector, subsector, etabin;
00667                             
00668                             if (geom.getTower(position, sector, subsector, etabin))
00669                             {
00670                             
00671                                 bool match = gtower->sector() == sector;
00672                                 match &= gtower->subsector() == subsector;
00673                                 match &= gtower->etabin() == etabin;
00674                                 
00675                                 if(match)
00676                                 {
00677                                     mCandidate->addMyTrack(track);
00678                                     track->candidates.Add(mCandidate);
00679                                 }
00680                     
00681                             }
00682                             
00683                         }
00684                     
00685                     
00686                     } // Tracks
00687                     
00688                 } // if(gtower)
00689                 
00690             } // Cluster towers
00691 
00692 
00693             // Set preshower and postshower energy sums
00694             mCandidate->SetPre1Energy(epre1);
00695             mCandidate->SetPre2Energy(epre2);
00696             mCandidate->SetPostEnergy(epost);
00697 
00698             // Now associate tracks, towers, etc, lying within a cone
00699             // of mRadium around the candidate to the candidate
00700             // Only these objects will be stored in the StGammaEvent
00701 
00702             Float_t candidateEta = pcluster.Eta();
00703             Float_t candidatePhi = pcluster.Phi();
00704             
00705             // Check all tracks
00706             for(Int_t i = 0; i < mGammaEvent->numberOfTracks(); i++)
00707             {
00708             
00709                 StGammaTrack *track = mGammaEvent->track(i);
00710                 if (!track) continue;
00711                 
00712                 Float_t trackEta = track->eta();
00713                 Float_t trackPhi = track->phi();
00714                 Float_t dEta = trackEta - candidateEta;
00715                 Float_t dPhi = TVector2::Phi_mpi_pi(trackPhi - candidatePhi);
00716                 
00717                 Float_t r = TMath::Sqrt(dPhi * dPhi + dEta * dEta);
00718                 
00719                 if(r <= mRadius) 
00720                 {
00721                     mCandidate->addTrack(track);
00722                     if(!track->candidates.FindObject(mCandidate)) track->candidates.Add(mCandidate);
00723                 }
00724             
00725             }
00726 
00727             // Check all towers
00728             for(Int_t i = 0; i < mGammaEvent->numberOfTowers(); i++)
00729             {
00730             
00731                 StGammaTower *tower = mGammaEvent->tower(i);
00732                 if(!tower) continue;
00733                 
00734                 Float_t towerEta = tower->eta;
00735                 Float_t towerPhi = tower->phi;
00736                 Float_t dEta = towerEta - candidateEta;
00737                 Float_t dPhi = TVector2::Phi_mpi_pi(towerPhi - candidatePhi);
00738                 
00739                 Float_t r = TMath::Sqrt(dPhi * dPhi + dEta * dEta);
00740                 
00741                 if(r <= mRadius) 
00742                 {
00743                     mCandidate->addTower(tower);
00744                     if(!tower->candidates.FindObject(mCandidate)) tower->candidates.Add(mCandidate);
00745                 }
00746             
00747             }
00748 
00749             // Check all preshowers
00750             for(Int_t i = 0; i < mGammaEvent->numberOfPreshower1(); i++)
00751             {
00752             
00753                 StGammaTower *tower = mGammaEvent->preshower1(i);
00754                 if(!tower) continue;
00755                 
00756                 Float_t towerEta = tower->eta;
00757                 Float_t towerPhi = tower->phi;
00758                 Float_t dEta = towerEta - candidateEta;
00759                 Float_t dPhi = TVector2::Phi_mpi_pi(towerPhi - candidatePhi);
00760                 
00761                 Float_t r = TMath::Sqrt(dPhi * dPhi + dEta * dEta);
00762                 
00763                 if(r <= mRadius) 
00764                 {
00765                     mCandidate->addPreshower1(tower);
00766                     if(!tower->candidates.FindObject(mCandidate)) tower->candidates.Add(mCandidate);
00767                 }
00768             
00769             }
00770 
00771 
00772             for(Int_t i = 0; i < mGammaEvent->numberOfPreshower2(); i++)
00773             {
00774             
00775                 StGammaTower *tower = mGammaEvent->preshower2(i);
00776                 if(!tower) continue;
00777             
00778                 Float_t towerEta = tower->eta;
00779                 Float_t towerPhi = tower->phi;
00780                 Float_t dEta = towerEta - candidateEta;
00781                 Float_t dPhi = TVector2::Phi_mpi_pi(towerPhi - candidatePhi);
00782                 
00783                 Float_t r = TMath::Sqrt(dPhi * dPhi + dEta * dEta);
00784             
00785                 if(r <= mRadius) 
00786                 {
00787                     mCandidate->addPreshower2(tower);
00788                     if(!tower->candidates.FindObject(mCandidate)) tower->candidates.Add(mCandidate);
00789                 }
00790             
00791             }
00792     
00793             // Check all postshowers
00794             for(Int_t i = 0; i < mGammaEvent->numberOfPostshower(); i++)
00795             {
00796             
00797                 StGammaTower *tower = mGammaEvent->postshower(i);
00798                 if (!tower) continue;
00799                 
00800                 Float_t towerEta = tower->eta;
00801                 Float_t towerPhi = tower->phi;
00802                 Float_t dEta = towerEta - candidateEta;
00803                 Float_t dPhi = TVector2::Phi_mpi_pi(towerPhi - candidatePhi);
00804                 
00805                 Float_t r = TMath::Sqrt(dPhi * dPhi + dEta * dEta);
00806                 
00807                 if(r <= mRadius) 
00808                     {
00809                     mCandidate->addPostshower(tower);
00810                     if(!tower->candidates.FindObject(mCandidate)) tower->candidates.Add(mCandidate);
00811                 }
00812             
00813             }
00814 
00815             // Now associate SMD strips within mSMDRange of the cluster, being
00816             // care in how one defines the center from which mSMDRange is defined
00817 
00818             Float_t umin[kEEmcNumSectors], vmin[kEEmcNumSectors];
00819             Float_t umax[kEEmcNumSectors], vmax[kEEmcNumSectors];
00820             Float_t umid[kEEmcNumSectors], vmid[kEEmcNumSectors];
00821             Int_t ntow[kEEmcNumSectors];
00822 
00823             // First, find the geometric center of the tower cluster in each sector
00824             for(Int_t ii = 0; ii < 12; ii++)
00825             { 
00826                 umin[ii] = 0.0; 
00827                 vmin[ii] = 0.0; 
00828                 umax[ii] = 0.0; 
00829                 vmax[ii] = 0.0; 
00830                 umid[ii] = -1.0; 
00831                 vmid[ii] = -1.0; 
00832                 ntow[ii] = 0; 
00833             }
00834 
00835             EEmcSmdMap *eemap = EEmcSmdMap::instance();
00836             for(Int_t itow = 0; itow < cluster.numberOfTowers(); itow++)
00837             {
00838             
00839                 StEEmcTower tower = cluster.tower(itow);
00840                 Int_t U, V;
00841                 eemap->getMiddleU( tower.sector(), tower.subsector(), tower.etabin(), U );
00842                 eemap->getMiddleV( tower.sector(), tower.subsector(), tower.etabin(), V );
00843                 ntow[ tower.sector() ]++;
00844                 umid[ tower.sector() ]+=0.5+(Float_t)U;// from middle of strip
00845                 vmid[ tower.sector() ]+=0.5+(Float_t)V;
00846             
00847             }
00848 
00849 
00850             for(Int_t isec = 0; isec < 12; isec++)
00851             {
00852             
00853                 if(ntow[isec]) 
00854                 { 
00855                     umid[isec] /= ntow[isec]; 
00856                     vmid[isec] /= ntow[isec]; 
00857                     umin[isec] = TMath::Max(umid[isec] - mEsmdRange * 2.0, 0. );
00858                     vmin[isec] = TMath::Max(vmid[isec] - mEsmdRange * 2.0, 0. );
00859                     umax[isec] = TMath::Min(umid[isec] + mEsmdRange * 2.0, 287. );
00860                     vmax[isec] = TMath::Min(vmid[isec] + mEsmdRange * 2.0, 287. );
00861                 }
00862                 
00863             }
00864 
00865             // Now associate SMD strips within mSMDRange of the center
00866             for(Int_t isec = 0; isec < 12; isec++)
00867             {
00868             
00869                 if(!ntow[isec]) continue;
00870                 
00871                 // Start with the U plane
00872                 for(Int_t i = (Int_t)umin[isec]; i < (Int_t)umax[isec]; i++)
00873                 {
00874                 
00875                     StGammaStrip *strip = mGammaRawMaker->strip(isec, 0, i);
00876                     if(strip) 
00877                     {
00878                         strip->position = i;
00879                         mCandidate->addSmdu(strip);
00880                         strip->candidates.Add(mCandidate);
00881                     }
00882                     
00883                 }
00884                 
00885                 // Finish with the V plane
00886                 for(Int_t i = (Int_t)vmin[isec]; i < (Int_t)vmax[isec]; i++)
00887                 {
00888                 
00889                     StGammaStrip *strip = mGammaRawMaker->strip(isec,1,i);
00890                     if(strip) 
00891                     {
00892                     strip->position = i;
00893                     mCandidate->addSmdv(strip);
00894                     strip->candidates.Add(mCandidate);
00895                     }
00896                 
00897                 }
00898                 
00899             } // Sectors (SMD)
00900 
00901             // Sum SMDU energy
00902             float smduEnergy = 0;
00903             
00904             for(int i = 0; i < mCandidate->numberOfSmdu(); ++i)
00905             {
00906                 StGammaStrip* strip = mCandidate->smdu(i);
00907                 smduEnergy += strip->energy;
00908             }
00909             
00910             mCandidate->SetSmduEnergy(smduEnergy);
00911             
00912             // Sum SMDV energy
00913             float smdvEnergy = 0;
00914             for (int i = 0; i < mCandidate->numberOfSmdv(); ++i)
00915             {
00916                 StGammaStrip* strip = mCandidate->smdv(i);
00917                 smdvEnergy += strip->energy;
00918             }
00919             
00920             mCandidate->SetSmdvEnergy(smdvEnergy);
00921 
00922         } // Clusters
00923 
00924     } // Sectors (Clusters)
00925 
00926     // Run the gamma fitter
00927     for(int i = 0; i < mGammaEvent->numberOfCandidates(); ++i) 
00928     {
00929     
00930         StGammaCandidate* candidate = mGammaEvent->candidate(i);
00931         if (candidate->detectorId() == StGammaCandidate::kEEmc) 
00932         {
00933             StGammaFitterResult ufit;
00934             StGammaFitterResult vfit;
00935             int ustatus = StGammaFitter::instance()->fit(candidate, &ufit, 0);
00936             if(ustatus == 0) candidate->SetSmdFit(ufit,0);
00937             int vstatus = StGammaFitter::instance()->fit(candidate, &vfit, 1);
00938             if(vstatus == 0) candidate->SetSmdFit(vfit,1);
00939         }
00940     
00941     }
00942     
00943     return kStOK;
00944     
00945 }
00946 
00948 //   Return the position of an EEMC Cluster     //
00950 TVector3 StGammaCandidateMaker::getEEmcClusterPosition(const StEEmcCluster& cluster)
00951 {
00952 
00953     // Get gamma raw maker
00954     StGammaRawMaker *mGammaRawMaker = dynamic_cast<StGammaRawMaker*>(GetMakerInheritsFrom("StGammaRawMaker"));
00955     if(!mGammaRawMaker)
00956     {
00957         LOG_WARN << "MakeEndcap() - No StGammaRawMaker found!" << endm;
00958         return kStWarn;
00959     }
00960     
00961     // Get cluster seed tower
00962     StEEmcTower seed = cluster.tower(0);
00963     
00964     // Get ranges of ESMD strips under seed tower
00965     int sector    = seed.sector();
00966     int subsector = seed.subsector();
00967     int etabin    = seed.etabin();
00968     int xmin[2], xmax[2];
00969     
00970     EEmcSmdMap::instance()->getRangeU(sector, subsector, etabin, xmin[0], xmax[0]);
00971     EEmcSmdMap::instance()->getRangeV(sector, subsector, etabin, xmin[1], xmax[1]);
00972     
00973     // Find max strips within ranges
00974     StGammaStrip* maxStrips[2] = {0, 0};
00975     for(int plane = 0; plane < 2; ++plane) 
00976     {
00977     
00978         for(int i = xmin[plane]; i <= xmax[plane]; ++i) 
00979         {
00980         
00981             StGammaStrip* strip = mGammaRawMaker->strip(sector, plane, i);
00982             if(strip) 
00983             {
00984                 if(!maxStrips[plane] || strip->energy > maxStrips[plane]->energy) 
00985                 {
00986                     maxStrips[plane] = strip;
00987                 }
00988             }
00989             
00990         }
00991         
00992     }
00993 
00994     // Get intersection of max strips and check that it lies within
00995     // the fiducial volume of the seed tower.
00996     if(maxStrips[0] && maxStrips[1]) 
00997     {
00998     
00999         TVector3 position = EEmcSmdGeom::instance()->getIntersection(sector, maxStrips[0]->index, maxStrips[1]->index);
01000         int sec, sub, eta;
01001         
01002         bool success = position.z() != -999;
01003         success &= EEmcGeomSimple::Instance().getTower(position, sec, sub, eta);
01004         success &= sector == sec;
01005         success &= subsector == sub;
01006         success &= etabin == eta;
01007         
01008         if(success) return position;
01009         
01010     }
01011     
01012     // Failed to calculate cluster position using ESMD strips
01013     return TVector3(0, 0, 0);
01014   
01015 }
01016 
01018 // Clean up the StGammaEvent, dropping detector //
01019 //    objects not associated with a candidate   //
01021 Int_t StGammaCandidateMaker::Compress()
01022 {
01023 
01024     // No compression
01025     if(mCompressLevel == 0) return kStOk;
01026 
01027     // Retrieve the StGammaEvent via StGammaEventMaker
01028      StGammaEventMaker *mGammaEventMaker = dynamic_cast<StGammaEventMaker*>(GetMakerInheritsFrom("StGammaEventMaker"));
01029     if(!mGammaEventMaker)
01030     {
01031         LOG_WARN << "Compress() - No StGammaEventMaker found!" << endm;
01032         return kStWarn;
01033     }
01034 
01035     StGammaEvent *mGammaEvent = mGammaEventMaker->event();
01036     if(!mGammaEvent) 
01037     {
01038         LOG_WARN << "Compress() - No StGammaEvent found!" << endm;
01039         return kStWarn;
01040     }
01041 
01042     Compress<StGammaStrip>(mGammaEvent->mStrips);
01043     
01044     // Compress SMD strips only
01045     if(mCompressLevel == 1) return kStOk;  
01046     
01047     // Drop strips, tracks, towers without candidates
01048     Compress<StGammaTrack>(mGammaEvent->mTracks);
01049     Compress<StGammaTower>(mGammaEvent->mTowers);
01050     Compress<StGammaTower>(mGammaEvent->mPreshower1);
01051     Compress<StGammaTower>(mGammaEvent->mPreshower2);
01052     Compress<StGammaTower>(mGammaEvent->mPostshower);
01053     
01054     // Compress all
01055     if(mCompressLevel == 2) return kStOk;
01056     
01057     LOG_WARN << Form("Unknown compression level (%d)", mCompressLevel) << endm;
01058     
01059     return kStWarn;
01060     
01061 }
01062 
01064 //  Template function for dropping objects in   //
01065 //   the TClonesArray without an association    //
01066 //             to a StGammaCandidate            //
01068 template<class T>
01069 void StGammaCandidateMaker::Compress(TClonesArray* clones)
01070 {
01071     TIter next(clones);
01072     while (T* x = (T*)next()) if (x->candidates.IsEmpty()) clones->Remove(x);
01073     clones->Compress();
01074 }