StRoot: StEEmcPool/StEEmcIUPi0/StEEmcIUPointMaker.cxx Source File

00001 
00026 #include "StEEmcIUPointMaker.h" 
00027 #include "StEEmcPool/StEEmcA2EMaker/StEEmcA2EMaker.h"
00028 #include "StEEmcIUClusterMaker.h"
00029 #include "StEEmcIUCluster.h"
00030 #include "StEEmcIUSmdCluster.h"
00031 #include "StEEmcUtil/EEmcGeom/EEmcGeomDefs.h"
00032 #include "StEEmcUtil/EEmcGeom/EEmcGeomSimple.h"
00033 #include "StEEmcUtil/StEEmcSmd/EEmcSmdGeom.h" 
00034 #include "TH1F.h"
00035 #include "TH2F.h"
00036 #include <iostream> 
00037 #include <algorithm>
00038 #include <map>
00039 
00040 #include "StEEmcPool/StEEmcPointMaker/eeTowerFunction.h"
00041 
00042 /* StEvent stuff */
00043 #include "StEvent/StEvent.h"
00044 #include "StEvent/StEmcCollection.h"
00045 #include "StEvent/StEmcPoint.h"
00046 
00047 /* Root's linear algebra package */
00048 #include "TMatrixF.h"
00049 
00050 //#define DEBUG_buildPoints
00051 
00052 ClassImp(StEEmcIUPointMaker);
00053 
00054 // ----------------------------------------------------------------------------
00055 StEEmcIUPointMaker::StEEmcIUPointMaker(const Char_t *name):StMaker(name)
00056 {
00057     std::cout << "StEEmcIUPointMaker("<<name<<")" << std::endl << std::flush; 
00058 
00062     mEEtow=new EEmcGeomSimple();
00063     mEEsmd=EEmcSmdGeom::instance();
00064     mEEmap=EEmcSmdMap::instance(); 
00065 
00066     mTowerThreshold=0.;
00067     mFillStEvent=false;
00068     mEnergyMode=1;
00069     mLimit=10;
00070     mSmdMatch = 0.;
00071     ukey=0;
00072     vkey=0;
00073 
00074 }
00075 
00076 // ----------------------------------------------------------------------------
00077 Int_t StEEmcIUPointMaker::Init()
00078 {
00079     mEEanalysis=(StEEmcA2EMaker*)GetMaker(mNameAnalysis);
00080     mEEclusters=(StEEmcIUClusterMaker*)GetMaker(mNameClusters);
00081     assert(mEEanalysis);
00082     assert(mEEclusters);
00083     hZEratio=new TH1F("smdezratio","energy Zratio between smd clusters",10,0.,1.);
00084     pointet=new TH1F("pointet","transverse energy of every point candidate",40,0,10);
00085     pointgeo=new TH2F("pointgeo", "Reconstructed Point Geometry; phi / deg; Eta",360,-180.,180,20,1.,2.);
00086     return StMaker::Init();
00087 }
00088 
00089 // ----------------------------------------------------------------------------
00090 Int_t StEEmcIUPointMaker::Make()
00091 {
00092  
00099   
00100   // Loop over all 12 EEMC sectors
00101   for ( Int_t sector=0; sector<12; sector++ )
00102     {
00103 
00105       StEEmcIUSmdClusterVec_t uclusters=mEEclusters->smdclusters(sector,0);
00106       StEEmcIUSmdClusterVec_t vclusters=mEEclusters->smdclusters(sector,1);
00107 
00109       std::sort( uclusters.begin(), uclusters.end(), inner );
00110       std::sort( vclusters.begin(), vclusters.end(), inner );
00111 
00112       findPoints( sector, uclusters, vclusters, mPoints );
00113       
00114     }
00115 
00116 
00123   for ( UInt_t i=0;i<mPoints.size();i++ )
00124     {
00125       
00126       mPoints[i].energy( (Float_t)(mPoints[i].energy()/0.007/2.) );
00127      
00128       TVector3 pointpos=mPoints[i].position();
00129       pointpos=pointpos.Unit();
00130       Float_t pet = (mPoints[i].energy()*pointpos).Perp();
00131       pointet->Fill(pet);
00132       pointgeo->Fill(mPoints[i].position().Phi()/3.1416*180,mPoints[i].position().Eta());
00133     }
00134   
00135 
00136 
00137 
00138   StEEmcIUPointVec_t orgpoints=mPoints;
00139 
00141   //  shareEnergy(); //  <<<<<<< leads to negative point energies... rethink
00142 
00143   if ( mEnergyMode == 1 )
00144     shareEnergySimple();
00145   else
00146     shareEnergySmd();
00147 
00148 
00150   countRelatives();
00151   
00152 
00153   if ( mFillStEvent )
00154     {
00155       fillStEvent();
00156       verifyStEvent();
00157     }
00158 
00159   /*
00160   for ( UInt_t i=0;i<mPoints.size();i++ )
00161     {
00162       mPoints[i].print();
00163     }
00164   for ( Int_t i=0;i<10;i++ ) std::cout << std::endl;
00165 */
00166     
00167   return kStOK;
00168 }
00169 
00170 // ----------------------------------------------------------------------------
00171 
00172 StEEmcIUPointVec_t StEEmcIUPointMaker::buildSmdPoints( Int_t sector, 
00173                                                    StEEmcIUSmdClusterVec_t &u, 
00174                                                    StEEmcIUSmdClusterVec_t &v )
00175 {
00176 
00177   StEEmcIUPointVec_t points;
00178   //printf("u.size=%d v.size=%d\n",u.size(),v.size());
00179   for ( UInt_t i=0; i<u.size(); i++ ) 
00180     {
00181 
00182       StEEmcIUSmdCluster uc=u[i];
00183       
00184       Float_t xu=uc.mean();
00185 
00186 
00187       for ( UInt_t j=0;j<v.size(); j++ ) 
00188         {
00189 
00190           
00191           StEEmcIUSmdCluster vc=v[j];
00192           Float_t xv=vc.mean();
00193 
00195           TVector3 direct = mEEsmd->getIntersection( sector, xu, xv );
00196           
00197           Int_t sec,sub,eta;
00198           if ( !mEEtow->getTower(direct,sec,sub,eta) )
00199             {
00200               continue; 
00201             }
00202           else
00203             {
00205             }
00206 
00209           if ( sector != sec )
00210             continue;
00211 
00212 
00213 
00217           Bool_t good = false;
00218           if ( mEEanalysis->tower(sec,sub,eta).energy() > mTowerThreshold ) 
00219             good=true;
00220           else if ( mEEanalysis->tower(sec,sub,eta).fail() ) 
00221             {
00222               for ( Int_t layer=1;layer<=3;layer++ ) 
00223                 {
00224                   if ( mEEanalysis->tower(sec,sub,eta,layer).energy() > 0. )
00225                     good=true;          
00226                 }
00227             }
00228 
00231           /*
00232           Float_t uvdiff = TMath::Abs( uc.energy() - vc.energy() );
00233           Float_t uvavg  = 0.5 * ( uc.energy() + vc.energy() );
00234           if ( uvdiff > mSmdMatch * uvavg ) good = false;
00235           */
00236           Float_t Eu=uc.energy();
00237           Float_t Ev=vc.energy();
00238           if ( Eu < mSmdMatch * Ev || Ev < mSmdMatch * Eu ) good=false;
00239           
00240           if ( good ) {
00241 
00242             StEEmcIUPoint p;
00243             p.cluster( uc, 0 );
00244             p.cluster( vc, 1 );
00245             p.energy( uc.energy() + vc.energy() );
00246             p.tower( mEEanalysis->tower(sec,sub,eta) );
00247             TVector3 position=mEEsmd->getIntersection(sector,uc.mean(),vc.mean());
00248             p.position(position);
00249             points.push_back(p);
00250 
00252             mSmdPoints.push_back(p);
00253 
00254 
00255           }
00256                         
00257         }
00258 
00259     }
00260   
00261   return points;
00262 
00263 }
00264 
00265 
00266 
00267 
00268 
00269 
00270 
00271 
00272 
00273 
00274 
00275 
00276 
00277 
00278 
00279 
00280 
00281 // ----------------------------------------------------------------------------
00282 Bool_t StEEmcIUPointMaker::findPoints( Int_t sector, 
00283                                      StEEmcIUSmdClusterVec_t uclusters, 
00284                                      StEEmcIUSmdClusterVec_t vclusters,
00285                                      StEEmcIUPointVec_t &points )
00286 {
00287 
00288 
00290   StEEmcIUPointVec_t Tmypoints;
00291   StEEmcIUPointVec_t mypoints;
00293   std::sort( uclusters.begin(), uclusters.end(), inner );
00294   std::sort( vclusters.begin(), vclusters.end(), inner );
00295 
00297   StEEmcIUPointVec_t smdpoints = buildSmdPoints( sector, uclusters, vclusters );
00298 
00299  
00301   if ( smdpoints.size() < 1 ) return false;
00302   
00304   std::sort( smdpoints.begin(), smdpoints.end(), chiSquare );
00305   
00306 
00310   std::map< Int_t, std::vector<Int_t> > u2p, v2p;            
00311   for ( UInt_t i=0; i<smdpoints.size(); i++ ) 
00312     {
00313       u2p[ smdpoints[i].cluster(0).key() ].push_back( i );
00314       v2p[ smdpoints[i].cluster(1).key() ].push_back( i );
00315       
00316           
00317     }
00318 
00320   Bool_t go = false;
00321   StEEmcIUSmdClusterVec_t::iterator uiter=uclusters.begin();
00322   StEEmcIUSmdClusterVec_t::iterator viter=vclusters.begin();
00323 
00324 
00329   
00330   //    ----------<<<<<<<<< stage one >>>>>>>>>-------------
00331  
00332   int iii=0,jjj=0;
00333 
00334   while ( uiter<uclusters.end() || viter<vclusters.end() )
00335     {
00336 
00338       StEEmcIUSmdCluster ucl;ucl.key(-1);
00339       StEEmcIUSmdCluster vcl;vcl.key(-1);
00340       if ( uiter<uclusters.end() ) ucl=(*uiter);
00341       if ( viter<vclusters.end() ) vcl=(*viter);
00342       Int_t iUV=-1;
00343       if ( ucl.key()<0 )
00344         iUV=1;
00345       else if ( vcl.key()<0 )
00346         iUV=0;
00347       else if ( (*uiter).mean() < (*viter).mean() )
00348         iUV=0;
00349       else
00350         iUV=1;
00351 
00353       StEEmcIUSmdCluster &cluster=(iUV==0)?ucl:vcl;
00354 
00356       std::vector<Int_t> matched=(iUV==0)?
00357         u2p[cluster.key()]:
00358         v2p[cluster.key()];
00359 
00360 
00363       
00364       if (  matched.size()==0 || matched.size() >1 )
00365         {
00366           if ( iUV==0 ) uiter++;
00367           if ( iUV==1 ) viter++;
00368           continue;
00369         }
00370 
00375       go=true;
00376         
00379       StEEmcIUPoint p=smdpoints[matched.back()];
00380       
00382       Tmypoints.push_back(p);
00383       
00384       if ( iUV==0 ) {uiter++;iii++;}
00385       if ( iUV==1 ) {viter++;jjj++;}
00386       
00387     }//end of while
00388  
00389 
00390   //splitting algo starts. When two points ramp into one lane, it only shows us one Smd cluster information in one of the Smd planes, thus can only reconstruct one point. We now split the overlapped U/V Smd cluster based on the other two Smd cluster in V/U accordingly. 
00391   std::sort( Tmypoints.begin(), Tmypoints.end(), chiSquare );
00392   //>=2 to turn on splitting algo, <0 to off splitting algo.
00393   if(Tmypoints.size()>=2)
00394     { 
00395       int aa=0,bb=0;
00396       for ( UInt_t i=0; i<Tmypoints.size(); i++ )
00397         {
00398           StEEmcIUSmdCluster uone=Tmypoints[i].cluster(0);
00399 
00400           StEEmcIUSmdCluster vone=Tmypoints[i].cluster(1);
00401 
00402           
00403           for ( UInt_t j=i+1; j<Tmypoints.size(); j++ )
00404             {     
00405               StEEmcIUSmdCluster utwo=Tmypoints[j].cluster(0);
00406               StEEmcIUSmdCluster vtwo=Tmypoints[j].cluster(1);
00407               
00408               if((uone.mean()==utwo.mean())||(vone.mean()==vtwo.mean()))
00409                 {
00410                   //case 1
00411                   if(uone.mean()==utwo.mean()&&(aa==0)&&(vone.mean()!=vtwo.mean()))
00412                     {
00413                       float uzr=fabs(uone.energy()-vone.energy()-vtwo.energy())/(uone.energy()+vone.energy()+vtwo.energy());
00414                      
00415                       hZEratio->Fill(uzr);
00416                      
00417                       if(uzr<=0.2)
00418                         {
00419                           float vratio1=vone.energy()/(vone.energy()+vtwo.energy());
00420                           float vratio2=vtwo.energy()/(vone.energy()+vtwo.energy());
00421                            
00422                           int ukey=uone.key();
00423                           int ukey1=ukey+100;
00424                           int ukey2=ukey+200;
00425                           float u1energy=uone.energy()*vratio1;
00426                           float u2energy=uone.energy()*vratio2;
00427                           uone.key(ukey1);
00428                           utwo.key(ukey2);
00429                           uone.energy(u1energy);
00430                           utwo.energy(u2energy);
00431                           float p1energy=uone.energy()+vone.energy();
00432                           Tmypoints[i].energy(p1energy);
00433                           float p2energy=utwo.energy()+vtwo.energy();
00434                           Tmypoints[j].energy(p2energy);
00435                           Tmypoints[i].cluster(uone,0);
00436                           Tmypoints[j].cluster(utwo,0);
00437                           StEEmcIUPoint p1=Tmypoints[i];
00438                           
00439                           //mypoints.push_back(p1);
00440                           StEEmcIUPoint p2=Tmypoints[j];
00441                           //mypoints.push_back(p2);
00442                           
00443                           removeCluster( uclusters, p1.cluster(0).key() );
00444                           removeCluster( vclusters, p1.cluster(1).key() );
00445                           points.push_back(p1);
00446                           removeCluster( uclusters, p2.cluster(0).key() );
00447                           removeCluster( vclusters, p2.cluster(1).key() );
00448                           points.push_back(p2);
00449                           removeCluster( uclusters, ukey );
00450                           
00451                           findPoints(sector, uclusters, vclusters, points );
00452                           return true;
00453                           
00454                           aa++;
00455                         }
00456                     } //end of case 1
00457                   //case 2
00458                   if(vone.mean()==vtwo.mean()&&(bb==0)&&(uone.mean()!=utwo.mean()))
00459                     {
00460                       float vzr=fabs(vone.energy()-uone.energy()-utwo.energy())/(vone.energy()+uone.energy()+utwo.energy());
00461                       hZEratio->Fill(vzr);
00462                       if(vzr<=0.2)
00463                         {
00464                           float uratio1=uone.energy()/(uone.energy()+utwo.energy());
00465                           float uratio2=utwo.energy()/(uone.energy()+utwo.energy());
00466                           
00467                           
00468                           int vkey=vone.key();
00469                           int vkey1=vkey+100;
00470                           int vkey2=vkey+200;
00471                           float v1energy=vone.energy()*uratio1;
00472                           float v2energy=vone.energy()*uratio2;
00473                           vone.key(vkey1);
00474                           vtwo.key(vkey2);
00475                           vone.energy(v1energy);
00476                           vtwo.energy(v2energy);
00477                           float p1energy=vone.energy()+uone.energy();
00478                           Tmypoints[i].energy(p1energy);
00479                           float p2energy=vtwo.energy()+utwo.energy();
00480                           Tmypoints[j].energy(p2energy);
00481                           Tmypoints[i].cluster(vone,1);
00482                           Tmypoints[j].cluster(vtwo,1);
00483                           StEEmcIUPoint p1=Tmypoints[i];
00484                           
00485                           //mypoints.push_back(p1);
00486                           StEEmcIUPoint p2=Tmypoints[j];
00487                           //mypoints.push_back(p2);
00488                           
00489                           removeCluster( uclusters, p1.cluster(0).key() );
00490                           removeCluster( vclusters, p1.cluster(1).key() );
00491                           points.push_back(p1);
00492                           removeCluster( uclusters, p2.cluster(0).key() );
00493                           removeCluster( vclusters, p2.cluster(1).key() );
00494                           points.push_back(p2);
00495                           removeCluster( vclusters, vkey );
00496                           
00497                           findPoints(sector, uclusters, vclusters, points );
00498                           return true;
00499                           
00500                           
00501                           bb++;
00502                         }
00503                     }//end of case 2
00504                   //case 3
00505                   if(uone.mean()==utwo.mean()&&(vone.mean()==vtwo.mean()))
00506                     {
00507                       int cc=0;
00508                       for ( UInt_t i=0; i<Tmypoints.size(); i++ )
00509                         {
00510                           StEEmcIUPoint p=Tmypoints[i];
00511                           mypoints.push_back(p);
00512                           cc++;
00513                         }
00514                 
00515                       //printf("caase 3 points.size=%d\n",mypoints.size());
00516                     }//end of case 3
00517                 }
00518             }
00519         }
00520 
00521       //case 4 catch those isolated points
00522       for ( UInt_t i=0; i<Tmypoints.size(); i++ )
00523         {
00524           StEEmcIUPoint pp=Tmypoints[i];
00525           mypoints.push_back(pp);
00526         }//end of case 4
00527     }
00528   else
00529     {
00530       for ( UInt_t i=0; i<Tmypoints.size(); i++ )
00531         {
00532           StEEmcIUPoint p=Tmypoints[i];
00533           mypoints.push_back(p);
00534         }
00535     }
00538   Float_t chisq=9.0E9;
00539   Int_t  imin=-1;
00540   for ( UInt_t i=0; i<mypoints.size(); i++ )
00541     {
00542       Float_t eu=mypoints[i].cluster(0).energy();
00543       Float_t ev=mypoints[i].cluster(1).energy();
00544       Float_t x2=fabs((eu-ev)/(eu+ev));
00545       if ( x2 < chisq ) {
00546         imin=(Int_t)i;
00547         chisq=x2;
00548       }
00549     }
00550 
00551 
00558   if ( imin >= 0 ) {
00559 
00560     StEEmcIUPoint p=mypoints[imin];
00561  
00562     removeCluster( uclusters, p.cluster(0).key() );
00563     removeCluster( vclusters, p.cluster(1).key() );
00564     points.push_back(p);
00565     findPoints(sector, uclusters, vclusters, points );
00566     return true;
00567 
00568   }
00569 
00570 
00571 
00572 
00573 
00575 
00576 
00583   uiter=uclusters.begin();
00584   viter=vclusters.begin();
00585   while ( uiter!=uclusters.end() || viter!=vclusters.end() )
00586     {
00587 
00589       StEEmcIUSmdCluster ucl;ucl.key(-1);
00590       StEEmcIUSmdCluster vcl;vcl.key(-1);
00591       if ( uiter!=uclusters.end() ) ucl=(*uiter);
00592       if ( viter!=vclusters.end() ) vcl=(*viter);
00593       Int_t iUV=-1;
00594       if ( ucl.key()<0 )
00595         iUV=1;
00596       else if ( vcl.key()<0 )
00597         iUV=0;
00598       else if ( (*uiter).mean() < (*viter).mean() )
00599         iUV=0;
00600       else
00601         iUV=1;
00602 
00604       StEEmcIUSmdCluster &cluster=(iUV==0)?ucl:vcl;
00605 
00607       std::vector<Int_t> matched=(iUV==0)?
00608         u2p[cluster.key()]:
00609         v2p[cluster.key()];
00610 
00613       if ( matched.size()==0 || matched.size()==1 )
00614         {
00615           if ( iUV==0 ) uiter++;
00616           if ( iUV==1 ) viter++;
00617           continue;
00618         }
00619       
00622       StEEmcIUPoint p=smdpoints[matched.front()];
00623       //StEEmcIUPoint p2=smdpoints[matched.front()++];
00625       mypoints.push_back(p);
00626       //mypoints.push_back(p2);
00627     
00628       if ( iUV==0 ) uiter++;
00629       if ( iUV==1 ) viter++;
00630 
00631     }
00632 
00633 
00636   chisq=9.0E9;
00637   imin=-1;
00638   for ( UInt_t i=0; i<mypoints.size(); i++ )
00639     {
00640  
00641       Float_t eu=mypoints[i].cluster(0).energy();
00642       Float_t ev=mypoints[i].cluster(1).energy();
00643       Float_t x2=fabs((eu-ev)/(eu+ev));
00644       if ( x2 < chisq ) {
00645         imin=(Int_t)i;
00646         chisq=x2;
00647       }
00648     }
00649   float chisq2=9.0E9;
00650   int imin2=-1,inn=-1;
00651   
00652   for ( UInt_t i=0; i<mypoints.size(); i++ )
00653     {
00654       inn=i;
00655       if(inn!=imin)
00656         {
00657 
00658           Float_t eu=mypoints[i].cluster(0).energy();
00659           Float_t ev=mypoints[i].cluster(1).energy();
00660           Float_t x22=(eu-ev)*(eu-ev);
00661           if ( x22 < chisq2 ) {
00662             imin2=i;
00663             chisq2=x22;
00664           }
00665         }
00666     }
00667 
00668 
00669 
00676   if ( imin >= 0 ) {
00677 
00678     StEEmcIUPoint p1=mypoints[imin];
00679  
00680     removeCluster( uclusters, p1.cluster(0).key() );
00681     removeCluster( vclusters, p1.cluster(1).key() );
00682  
00683     points.push_back(p1);
00684  
00685     findPoints(sector, uclusters, vclusters, points );
00686     return true;
00687 
00688   }
00689 
00690 
00691   return true;
00692 
00693 }
00694 
00695 
00696 
00697 // ----------------------------------------------------------------------------
00698 void  StEEmcIUPointMaker::Clear( Option_t *opts )
00699 {
00700 
00701   mEseen=0.;
00702   mPoints.clear(); 
00703   mSmdPoints.clear(); 
00704 
00705 }
00706 
00707 // ----------------------------------------------------------------------------
00708 void StEEmcIUPointMaker::removeCluster( StEEmcIUSmdClusterVec_t &clusters, Int_t k )
00709 {
00710   StEEmcIUSmdClusterVec_t::iterator iter=clusters.begin();
00711   while ( iter != clusters.end() )
00712     {
00713       if ( (*iter).key() == k ) {
00714         clusters.erase(iter);
00715         return;
00716       }
00717       iter++;
00718     }
00719 }
00720 
00721 // ----------------------------------------------------------------------------
00722 void StEEmcIUPointMaker::shareEnergy()
00723 {
00724 
00730 
00731 
00732 
00733   Int_t nrow=(Int_t)mPoints.size();
00734   Int_t ncol=(Int_t)mPoints.size();
00735   if ( nrow < 1 ) return; 
00736 
00737   std::vector<Float_t> Ef(nrow,0.);
00738 
00739   TMatrixF fractions(nrow,ncol);
00740 
00742   for ( Int_t k=0; k<mEEanalysis->numberOfHitTowers(0); k++ ) 
00743     {
00744 
00746       StEEmcTower tower=mEEanalysis->hittower(k,0);
00747 
00749       for ( UInt_t i=0; i< mPoints.size(); i++ ) 
00750         {
00751 
00754           Float_t fi=fracp2t( mPoints[i], tower );
00755           if ( fi<=0. ) continue;
00756 
00758           Ef[i] += fi * tower.energy();
00759 
00760           for ( UInt_t j=0; j<mPoints.size(); j++ )
00761             {
00762 
00765               Float_t fj=fracp2t( mPoints[j], tower );
00766               if (fi*fj<=0.) continue;
00767 
00768               fractions[i][j] += fi*fj;
00769               
00770             }
00771           
00772         }
00773 
00774     }
00775 
00776   fractions.Print();
00777 
00779   Double_t det = 0.;
00780   TMatrixF invert= fractions;
00781   invert.Invert(&det);
00782 
00783   invert.Print();
00784 
00785   TMatrixF test= fractions * invert;
00786 
00787   test.Print();
00788 
00789 
00791   
00792   std::vector<Float_t> epoints(nrow,0.);
00793 
00794   for ( Int_t i=0; i<nrow; i++ ) 
00795     { 
00796       for ( Int_t j=0; j<ncol; j++ )
00797         {
00798           epoints[i] += invert[i][j] * Ef[j];
00799         }      
00800 
00801 
00802     }
00803 
00804  
00805 
00806   
00807 
00808 
00809 }
00810 
00811 
00812 // ----------------------------------------------------------------------------
00813 Float_t StEEmcIUPointMaker::fracp2t( StEEmcIUPoint &p, StEEmcTower &t )
00814 {
00815 
00817 
00818 
00821   if ( !t.isNeighbor( p.tower(0) ) ) return 0.;
00822   //if ( !(t.index()==p.tower(0).index()) ) return 0.;
00823 
00824 
00826   Float_t xeta=(Float_t)t.etabin();
00827   Float_t xphi=(Float_t)t.phibin();
00828   Double_t X[]={xphi,xeta}; 
00829 
00831   Float_t xeta0=(Float_t)p.tower(0).etabin();
00832   Float_t xphi0=(Float_t)p.tower(0).phibin();
00833 
00836   Int_t sec,sub,eta;
00837   Float_t dphi,deta;
00838   if ( !mEEtow->getTower(p.position(), sec,sub,eta, dphi,deta ) ) return 0.;
00839   dphi/=2.0;
00840   deta/=2.0;
00841 
00843   Float_t xetap=xeta0+deta;
00844   Float_t xphip=xphi0+dphi;
00845   Double_t P[]={xphip,xetap,1.0};
00846 
00847   return eeTowerFunction( X, P );
00848 
00849 }
00850 
00851 
00852 // ----------------------------------------------------------------------------
00853 void StEEmcIUPointMaker::shareEnergySimple()
00854 {
00855 
00856   Int_t limit=mLimit; // algo quickly converges on energy
00857   Int_t count=0;
00858 
00859   while ( count++ < limit ) 
00860     {
00861 
00864       Float_t sumw[720]; 
00865       for (Int_t i=0;i<720;i++) sumw[i]=0.;
00866 
00869       for ( UInt_t i=0;i<mPoints.size();i++ )
00870         {
00871 
00872           StEEmcIUPoint point=mPoints[i];
00873           StEEmcTower tower=point.tower(0);
00874 
00875           sumw[ tower.index() ] += point.energy() * fracp2t( point, tower );
00876 
00878           for ( Int_t i=0; i<tower.numberOfNeighbors(); i++ )
00879             {
00880               StEEmcTower neighbor=tower.neighbor(i);
00881               sumw[ neighbor.index() ] += point.energy() * fracp2t( point, neighbor );
00882             }
00883 
00884         }
00885 
00886 
00890       
00891       for ( UInt_t i=0;i<mPoints.size();i++ )
00892         {
00893 
00894           StEEmcIUPoint &point=mPoints[i]; // note the reference
00895           StEEmcTower tower=point.tower(0);
00896           Float_t epoint=0.;
00897 
00898           Float_t frac=0.;
00899           if ( !tower.fail() && !tower.stat() && sumw[tower.index()]>0. )
00900           frac = point.energy() * fracp2t(point,tower) / sumw[ tower.index() ];
00901 
00902           epoint += tower.energy() * frac;
00903           
00906           for ( Int_t i=0; i<tower.numberOfNeighbors(); i++ )
00907             {
00908               StEEmcTower neighbor=tower.neighbor(i);
00909               if ( neighbor.stat() || neighbor.fail() || sumw[neighbor.index()]<=0. ) continue;
00910               frac = point.energy() * fracp2t(point,neighbor) / sumw[ neighbor.index() ];
00911               epoint += frac * neighbor.energy();
00912             }
00913           
00914           
00916           point.energy( epoint );
00917 
00918         }
00919 
00920     }
00921 
00922 
00923 
00924 
00925   std::vector<Bool_t> seen(720,false);
00926   for ( UInt_t i=0;i<mPoints.size();i++ )
00927     {
00928 
00929       StEEmcTower tow=mPoints[i].tower(0);
00930       if ( !seen[ tow.index() ] ) mEseen+=tow.energy();
00931       seen[ tow.index() ] = true;
00932 
00933       for ( Int_t j=0;j<tow.numberOfNeighbors();j++ )
00934         {
00935           StEEmcTower nei=tow.neighbor(j);
00936           if ( !seen[ nei.index() ] ) mEseen += nei.energy();
00937           seen[ nei.index() ] = true;
00938         }
00939 
00940     }
00941   
00942 
00943 
00944 
00945 
00946 
00947 }
00948 
00949 
00950 
00951 
00952 
00953 // ----------------------------------------------------------------------------
00954 void StEEmcIUPointMaker::fillStEvent()
00955 {
00956 
00957 
00958   StEvent *stevent=(StEvent*)GetInputDS("StEvent");
00959   if ( !stevent ) {
00960     Warning("fillStEvent","called, but no StEvent to be found");
00961     return;
00962   }
00963 
00964 
00966   for ( UInt_t i=0; i<mPoints.size(); i++ )
00967     {
00968 
00969       StEmcPoint *point=mPoints[i].stemc();
00970       stevent->emcCollection()->addEndcapPoint( point );
00971 
00972       mEtoEE[ point ] = mPoints[i];
00973 
00974 
00975     }
00976 
00977 }
00978 
00979 
00980 
00981 // ----------------------------------------------------------------------------
00982 void StEEmcIUPointMaker::verifyStEvent()
00983 {
00984 
00985   Float_t emc_sum_points = 0.;
00986   Float_t eemc_sum_points = 0.;
00987 
00988   StEvent *stevent=(StEvent*)GetInputDS("StEvent");
00989   if ( !stevent ) {
00990     Warning("verifyStEvent","called, but no StEvent to be found");
00991     return;
00992   }
00993 
00994   StSPtrVecEmcPoint& emcpts = stevent->emcCollection()->endcapPoints();
00995   for ( UInt_t i=0;i<emcpts.size();i++ )
00996     {
00997 
00998       StEmcPoint *p=emcpts[i];
00999       assert(p);
01000       emc_sum_points += p->energy();
01001 
01002     }
01003 
01004   for ( UInt_t i=0;i<mPoints.size();i++ )
01005     {
01006 
01007       StEEmcIUPoint p=mPoints[i];
01008       eemc_sum_points += p.energy();
01009 
01010     }
01011 
01012   std::cout << "StEEmcIUPointMaker point checksum: ";
01013   if ( emc_sum_points == eemc_sum_points )
01014     {
01015       std::cout << "passed";
01016     }
01017   else
01018     std::cout << "FAILED";
01019   std::cout << std::endl;
01020     
01021 
01022 }
01023 
01024 
01025 // ----------------------------------------------------------------------------
01026 void StEEmcIUPointMaker::countRelatives()
01027 {
01028 
01030   Int_t npoints[720];
01031   for ( Int_t i=0;i<720;i++ ) npoints[i]=0;
01032 
01033   for ( UInt_t i=0;i<mPoints.size();i++ )    
01034       npoints[ mPoints[i].tower(0).index() ]++;
01035 
01037   for ( UInt_t i=0;i<mPoints.size();i++ )
01038     {
01039 
01040       StEEmcTower tower=mPoints[i].tower(0);
01041       Int_t nn=tower.numberOfNeighbors();
01042 
01043       Int_t nrel=npoints[ tower.index() ] - 1; // don't count self
01044       assert(nrel>=0); // pbck
01045 
01046       for ( Int_t j=0;j<nn;j++ )
01047         {
01048           StEEmcTower t2=tower.neighbor(j);
01049           nrel+=npoints[ t2.index() ];
01050         }
01051 
01052       mPoints[i].numberOfRelatives(nrel);
01053 
01054     }
01055     
01056 
01057 }
01058 
01059 
01060 
01061 
01062 // ----------------------------------------------------------------------------
01063 void StEEmcIUPointMaker::shareEnergySmd()
01064 {
01065 
01068   Float_t sumw[720];for ( Int_t i=0;i<720;i++ )sumw[i]=0.;
01069   Float_t sumw1[720];for ( Int_t i=0;i<720;i++ )sumw1[i]=0.;
01070 
01071   for ( UInt_t ipoint=0;ipoint<mPoints.size();ipoint++ )
01072     {
01073 
01074       StEEmcIUPoint point=mPoints[ipoint];
01075       StEEmcTower tower=point.tower(0);
01076       sumw[tower.index()]+=point.energy();
01077       sumw1[tower.index()]+=point.energy(); 
01078       //printf("en==%f index=%d sumw=%f\n",point.energy(),tower.index(),sumw[tower.index()]);
01079       for ( Int_t itow=0;itow<tower.numberOfNeighbors();itow++ )
01080         {
01081           StEEmcTower t=tower.neighbor(itow);
01082           Int_t index=t.index();
01083          
01084           sumw[index]+=point.energy();
01085           //printf("index=%d sumw=%f\n",index,sumw[index]); 
01086         }
01087      
01088     }
01089   //printf("end of first loop.........................\n");
01092 
01093   for ( UInt_t ipoint=0;ipoint<mPoints.size();ipoint++ )
01094     {
01095       
01096       StEEmcIUPoint point=mPoints[ipoint]; // note reference
01097       StEEmcTower tower = point.tower(0);
01098       StEEmcTower pre1  = mEEanalysis->tower(tower.index(),1); 
01099       StEEmcTower pre2  = mEEanalysis->tower(tower.index(),2); 
01100       StEEmcTower post  = mEEanalysis->tower(tower.index(),3); 
01101       //hPostE->Fill(post.energy()*1000.0);
01102       Float_t epoint = 0.;
01103       Float_t epre1  = 0.; 
01104       Float_t epre2  = 0.; 
01105       Float_t epost  = 0.;  
01106 
01107       Int_t index = tower.index();
01108       epoint += tower.energy() * point.energy() / sumw[index];
01109       //printf("start of the calculation of a point\n");
01110       //printf("tower.energy=%f epoint=%f smde=%f sumw=%f index=%d\n",tower.energy(),epoint,point.energy(),sumw[index],index);
01111       epre1  += pre1.energy()  * point.energy() / sumw1[index]; 
01112       epre2  += pre2.energy()  * point.energy() / sumw1[index]; 
01113       epost  += post.energy()  * point.energy() / sumw[index]; 
01114       //printf("start of the neighbor calculation of a point\n");
01115       for ( Int_t itow=0;itow<tower.numberOfNeighbors();itow++ )
01116         {
01117           StEEmcTower t=tower.neighbor(itow);
01118 //        StEEmcTower p=pre1.neighbor(itow);
01119 //        StEEmcTower q=pre2.neighbor(itow);
01120           StEEmcTower r=post.neighbor(itow); 
01121 
01122           index = t.index();
01123           
01124           epoint += t.energy() * point.energy() / sumw[index];
01125           
01126           //printf("neighbortower.energy=%f epoint=%f smde=%f sumw=%f index=%d\n",t.energy(),epoint,point.energy(),sumw[index],index);
01127 //        epre1  += p.energy() * point.energy() / sumw[index];
01128 //        epre2  += q.energy() * point.energy() / sumw[index];
01129 //        epost  += r.energy() * point.energy() / sumw[index];
01130         }
01131 
01132 
01133       mPoints[ipoint].energy(epoint);
01134       mPoints[ipoint].energy(epre1,1);
01135       mPoints[ipoint].energy(epre2,2);
01136       mPoints[ipoint].energy(epost,3);
01137  
01138       
01139     }
01140   
01141   
01142 }
01143 
01144