StRoot: StEEmcPool/StEEmcIUPi0/StEEmcIUClusterMaker.cxx Source File

00001 
00050 #include "StEEmcIUClusterMaker.h"
00051 
00052 #include <algorithm>
00053 #include <iostream>
00054 
00055 #include "TMath.h"
00056 
00057 /* StEvent stuff */
00058 #include "StEvent/StEvent.h"
00059 #include "StEvent/StEmcCollection.h"
00060 #include "StEvent/StEmcDetector.h"
00061 #include "StEvent/StEmcModule.h"
00062 #include "StEvent/StEmcClusterCollection.h"
00063 #include "StEvent/StEmcCluster.h"
00064 
00065 #define DEBUG 0
00066 //#define MonteCarloS
00067 
00068 ClassImp(StEEmcIUClusterMaker);
00069 
00070 // ----------------------------------------------------------------------------
00071 StEEmcIUClusterMaker::StEEmcIUClusterMaker(const Char_t *name):StMaker(name)
00072 {
00073 
00075   mFillStEvent = 0;
00076 
00078   mSuppress = 0;
00079 
00080   mLoose=false;
00081   mSkip=true;
00082   //Meanclust=0;
00085   mClusterId = 0;
00086 
00088   const Float_t eseeds[] = { 0.6, 1.0/1000, 1.0/1000, 1.0/1000., 0.3/1000., 0.3/1000. };
00089   for ( Int_t i = 0; i < 6; i++ ) seedEnergy( eseeds[i], i );
00090 
00091   mMaxExtent = 3; 
00092   mSeedFloor = 2.0; 
00093 
00106 
00107   StEEmcIUClusterVec_t t;
00108   std::vector< StEEmcIUClusterVec_t > layers;
00109   for ( Int_t i = 0; i < 4; i++ ) layers.push_back(t);
00110   for ( Int_t i = 0; i < 12; i++ ) mTowerClusters.push_back(layers);
00111 
00112   StEEmcIUSmdClusterVec_t s;
00113   std::vector< StEEmcIUSmdClusterVec_t > planes;
00114   planes.push_back(s);
00115   planes.push_back(s);
00116   for ( Int_t i = 0; i < 12; i++ ) mSmdClusters.push_back(planes);
00117   for ( Int_t i = 0; i < 12; i++ ) TmSmdClusters.push_back(planes);
00118 
00119   StEEmcTowerVec_t tow;
00120   std::vector< StEEmcTowerVec_t > lay;
00121   for ( Int_t i = 0; i < 4; i++ ) lay.push_back(tow);
00122   for ( Int_t i = 0; i < 12; i++ ) mSeedTowers.push_back(lay);
00123 
00124   mEEtow=new EEmcGeomSimple();
00125   mEEsmd=EEmcSmdGeom::instance();
00126   mEEmap=EEmcSmdMap::instance();
00127 
00128   mMinStrips=3;
00129   
00130 }
00131 
00132 // ----------------------------------------------------------------------------
00133 Int_t StEEmcIUClusterMaker::Init()
00134 {
00135   mEEanalysis=(StEEmcA2EMaker*)GetMaker(mAnalysisName);
00136   assert(mEEanalysis);
00137   clusize=new TH1D("clustersize", "smd cluster size ",10,0,10);
00138   tclusize=new TH1D("tclustersize", "smd cluster size ",10,0,10);
00139   cludis=new TH1D("cludis", "cluster distance ",20,0,20);
00140   return StMaker::Init();
00141 }
00142 
00143 // ----------------------------------------------------------------------------
00144 Int_t StEEmcIUClusterMaker::Make()
00145 {
00146 
00148   if ( !buildTowerClusters() ) return kStWarn;
00149 
00151   if ( !buildSmdClusters() ) return kStWarn; 
00152   //fillStEvent();
00153 
00155   if ( mFillStEvent ) fillStEvent();
00156 
00158   if ( mFillStEvent ) 
00159     if ( !verifyStEvent() ) Warning("Make","StEvent not properly copied");
00160   //print();  
00161 
00162 
00163   return kStOK;
00164 }
00165 
00166 // ----------------------------------------------------------------------------
00167 void StEEmcIUClusterMaker::Clear( Option_t *opts )
00168 {
00169 
00171   for ( Int_t sector=0; sector<12; sector++ ) {
00172     for ( Int_t layer=0; layer<4; layer++ )
00173       mTowerClusters[sector][layer].clear();    
00174     for ( Int_t plane=0; plane<2; plane++ )
00175       {
00176         mSmdClusters[sector][plane].clear();
00177         TmSmdClusters[sector][plane].clear();
00178       }
00179   }
00180 
00181   for ( Int_t i=0;i<6;i++ ) 
00182     {
00183       mNumberOfClusters[i]=0;
00184       TmNumberOfClusters[i]=0;
00185     }
00186   mClusterId = 0;
00187 
00188   return;
00189 }
00190 
00191 
00192 // ----------------------------------------------------------------------------
00193 Bool_t StEEmcIUClusterMaker::buildTowerClusters()
00194 {
00195 
00196   // This part is not important in the reconstruction of pi0. It only provides some fundamental checks and basic idea about the tower clusters.
00199   //static const Int_t layer=0;
00200 
00201 
00203 
00204   eeen=0;
00205   ep1=0;
00206   ep2=0;
00207   ep3=0;
00208   for ( Int_t layer=0;layer<4;layer++)
00209   {
00210 
00213     Float_t weights[720]; for ( Int_t i=0;i<720;i++ ) weights[i]=0.;
00214 
00216     StEEmcIUClusterVec_t myClusters;
00217 
00218 
00220     StEEmcTowerVec_t towers = mEEanalysis -> towers( layer );
00221 
00222     
00223 
00225     std::sort(towers.begin(),towers.end());
00227     std::reverse(towers.begin(),towers.end());
00228    
00229 
00232     StEEmcTowerVec_t::iterator last = towers.begin();
00233     while ( last != towers.end() ) {
00234       if(layer==0){
00235         eeen+=(*last).energy();
00236       }
00237       if(layer==1){
00238         ep1+=(*last).energy();
00239       }
00240       if(layer==2){
00241         ep2+=(*last).energy();
00242       }
00243       if(layer==3){
00244         ep3+=(*last).energy();
00245       }
00246       if ( (*last).energy() < mSeedEnergy[layer] ) break;
00251       //$$$      mSeedTowers[layer][ (*last).sector() ].push_back((*last));
00252 #if DEBUG
00253       std::cout << "-- Seed tower ----------------------" << std::endl;
00254       (*last).print();
00255 #endif
00256       last++;
00257     }
00258     //std::cout<<"tower energy total = "<< eeen<< std::endl;
00259     
00260     StEEmcTowerVec_t::iterator iter = towers.begin();
00261     while ( iter != towers.end() ) {
00262 
00264       if ( iter==last ) break;
00265       
00266 
00273       if ( weights[ (*iter).index() ] > 0. ) {
00274         iter++;
00275         continue;
00276       }
00277 
00280       for ( Int_t in=0; in < (*iter).numberOfNeighbors(); in++ ) {
00281         StEEmcTower t=(*iter).neighbor(in);
00282         weights[ t.index() ] += (*iter).energy();
00283       }
00284 
00285      
00286       iter++;
00287 
00288     }// loop over towers to init weights
00289 
00290 
00292     iter=towers.begin();
00293     while ( iter != towers.end() ) {
00294 
00296       if ( iter==last ) break;
00297 
00301       if ( weights[ (*iter).index() ] > 0. ) {
00302         iter++;
00303         continue;
00304       }
00305 
00306       StEEmcIUCluster cluster;
00307       StEEmcTower seed=(*iter);
00308 #if DEBUG
00309       std::cout << "--- Clustering ----------------" << std::endl;
00310       seed.print();
00311 #endif
00312 
00313       TVector3 momentum;
00314       cluster.add(seed,1.0);
00315       UInt_t sec,sub,eta;
00316       sec=(UInt_t)seed.sector(); 
00317       sub=(UInt_t)seed.subsector(); 
00318       eta=(UInt_t)seed.etabin();
00319       TVector3 d=mEEtow->getTowerCenter(sec,sub,eta).Unit();
00320       momentum += ( seed.energy() * d );
00321 
00322       for ( Int_t in=0; in<seed.numberOfNeighbors(); in++ ) {
00323         StEEmcTower t=seed.neighbor(in);
00324         sec=(UInt_t)t.sector(); 
00325         sub=(UInt_t)t.subsector(); 
00326         eta=(UInt_t)t.etabin();
00327         d=mEEtow->getTowerCenter(sec,sub,eta).Unit();
00328         Float_t weight = seed.energy() / weights[ t.index() ];
00329         momentum += ( t.energy() * d * weight );
00330         cluster.add(t, weight);
00331 #if DEBUG
00332         std::cout << "adding " << t.name() << " E=" << t.energy() << " W=" << weight << std::endl;
00333 #endif
00334       }
00335 
00336 
00338       cluster.momentum( momentum );
00339       cluster.key( mClusterId++ );
00340 #if DEBUG
00341       cluster.print();
00342 #endif
00343 
00344       mTowerClusters[ seed.sector() ][ layer ].push_back( cluster );
00345       mNumberOfClusters[layer]++;
00346      
00347       iter++;
00348 
00349     }
00350           
00351   }// end of eventual loop over layers
00352   
00353 
00354   return true;
00355 
00356 }
00357 
00358 // ----------------------------------------------------------------------------
00359 Bool_t StEEmcIUClusterMaker::buildSmdClusters()
00360 {
00361 
00362   //this part is essential in the reconstruction of pi0 in the EEMC.
00363 
00364   Int_t max_extent = mMaxExtent;
00365   countUseed=0;
00366   countVseed=0;
00367   uswidth=0;
00368   vswidth=0;
00369   uamp=0;
00370   vamp=0;
00371   esmdu=0;
00372   esmdv=0;
00373 
00374 
00376   for ( Int_t sector=0; sector<12; sector++ ) {
00378 
00379     for ( Int_t plane=0; plane<2; plane++ ) {
00380       float ampe=0;
00382         Float_t floor[288]; for ( Int_t i=0; i<288; i++ ) floor[i]=0.;
00383 
00385         Float_t energy[288]; for ( Int_t i=0; i<288; i++ ) energy[i]=0.;
00386         //averaging energies
00387         //Float_t newenergy[288]; for ( Int_t i=0; i<288; i++ ) newenergy[i]=0.;
00389         StEEmcStripVec_t strips=mEEanalysis->strips(sector,plane);
00390         StEEmcStripVec_t seeds;
00392         std::sort(strips.begin(),strips.end());
00394         std::reverse(strips.begin(),strips.end());
00395 
00397         StEEmcStripVec_t::iterator istrip=strips.begin();
00398 
00399         while ( istrip != strips.end() ) {
00400           if ( (*istrip).stat()||(*istrip).fail() ) {
00401             istrip++;
00402             continue;
00403           }
00404           ampe=(*(strips.begin())).energy();
00405           energy[ (*istrip).index() ] = (*istrip).energy();
00406           //cout<<"energy="<<(*istrip).energy()<<endl;
00407           
00408           if((*istrip).energy()>=0)
00409             {
00410               //fprintf(fout2,"data%02d%d %d  %f\n",sector,plane,(*istrip).index(),(*istrip).energy());
00411               if(plane==0)
00412               {
00413                 esmdu+=(*istrip).energy();
00414                 uswidth++;
00415                 uamp=ampe;
00416               }
00417               if(plane==1)
00418               {
00419                 esmdv+=(*istrip).energy();
00420                 vswidth++;
00421                 vamp=ampe;
00422               }
00423             }
00424           istrip++;
00425         }
00426 
00427         //averaging algo for strip energies. We decide not using the averaging algo for 2006 long pp data after analyses. But the algo is still kept here for future user reference.
00428         //istrip=strips.begin();
00429         //for(int i=2;i<=283;i++){
00430         //newenergy[i]=energy[i-1]*0.25+energy[i]*0.5+energy[i+1]*0.25;
00431         //}
00432         //while ( istrip != strips.end() ) {
00433         //if ( (*istrip).stat()||(*istrip).fail() ) {
00434         //istrip++;
00435         //continue;
00436         //}
00437         //(*istrip).energy(newenergy[(*istrip).index()]);
00438           
00439         //istrip++;
00440         //}
00441         //end of averaging algo
00442 
00444         std::vector<Bool_t> seed_flags( strips.size(), false );
00445 
00446 
00448         Int_t nstrip=0;
00449         Int_t nseeds=0;
00450         istrip=strips.begin();
00451         while ( istrip!=strips.end() ) {
00452 
00453           Int_t index=(*istrip).index();
00454           Float_t eseed=(*istrip).energy();
00455           
00457           if ( index <= 3 || index >= 283 ) {
00458             istrip++;
00459             nstrip++;
00460             continue;
00461           }
00462 
00464           if ( mSkip )
00465             if ( (*istrip).fail() ) {
00466               istrip++;
00467               nstrip++;
00468               continue;
00469             }
00470 
00471 
00474           if ( eseed < (mSeedFloor*floor[ index ] + mSeedEnergy[4+plane]) ) {
00475             istrip++;
00476             nstrip++;
00477             continue;
00478           }
00479           //if(energy[index+1]<floor[index+1] && energy[index-1]<floor[index-1]){
00480           //istrip++;
00481           //nstrip++;
00482           //continue;
00483           //}
00485           //cout<<"succeed!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1"<<endl;
00486           seeds.push_back( (*istrip) );
00487           if(plane==0){countUseed++;}
00488           if(plane==1){countVseed++;}
00489           nstrip++;
00490           nseeds++;
00491 
00492           // Now we want to set up special floor values for the SMD shape based on the real data fluctuation. The floor setting, which was basically applied to minimize fluctuation around a used seed strip, is now re-set for real data. Based on the fluctuatin feature we found from real data analyses, the second plane of SMD layers always suppress twice fluctuations than the first plane by using two parameters floor_para1 and floor_para2 properly. This part works specifically for real data and pythia sample. For Monte-Carlo sample, we can use the same floor setting for both SMD planes because we don't find fluctuation from MC sample.
00493           float floor_para1=0.0;
00494           float floor_para2=0.0;
00495           if(sector==0||sector==3||sector==6||sector==9)
00496             {
00497               if(plane==0)
00498                 {
00499                   floor_para1=0.4;
00500                   floor_para2=0.2;
00501                 }
00502               else{
00503                 floor_para1=0.2;
00504                 floor_para2=0.1;
00505               }
00506             }
00507           if(sector==2 || sector ==1 || sector ==5 || sector==4 || sector==8 || sector==7 || sector==11 || sector==10)
00508             {
00509               if(plane==1)
00510                 {
00511                   floor_para1=0.4;
00512                   floor_para2=0.2;
00513                 }
00514               else{
00515                 floor_para1=0.2;
00516                 floor_para2=0.1;
00517               }
00518             }
00519 #ifdef MonteCarloS
00520           floor_para1=0.2;
00521           floor_para2=0.1;
00522 #endif
00528 #ifndef LOOSE_CUTS        
00530           for ( Int_t i=0; i < 288; i++ ) 
00531             {
00532               Int_t dx=TMath::Abs(index-i);
00534               if(dx<3) {floor[i]+=1.0*eseed;}
00536               if(dx>=3&&dx<5){floor[i]+=floor_para1*eseed;}
00538               if(dx>=5&&dx<11){floor[i]+=floor_para2*eseed;}
00540               if(dx>=11&&dx<21){floor[i]+=0.05*eseed;}
00541             }
00542 #else
00544           for ( Int_t i=0; i < 288; i++ ) {
00545 
00546             Int_t dx=TMath::Abs(index-i);
00547 
00549             if ( dx<7) 
00550               if ( 0.05 * eseed > floor[i] ) floor[i]+= 0.05 * eseed;
00551 
00552           }
00553 #endif
00554           istrip++;
00555 
00556         }// found all seed strips in that plane of that sector
00557 
00563         Bool_t owned[288]; for (Int_t i=0;i<288;i++) owned[i]=false;
00564         Bool_t xseed[288]; for (Int_t i=0;i<288;i++) xseed[i]=false;
00565 
00566 
00567         StEEmcStripVec_t::iterator iseed=seeds.begin();
00568         while ( iseed != seeds.end() ) {
00569 
00570           int flcount=0;
00572           Int_t index=(*iseed).index();
00573           if ( owned[index] || (mSuppress&&xseed[index]) ) {
00574             iseed++;
00575             continue;
00576           }
00577 
00579           owned[index]=true;
00580           xseed[index]=true;
00582           StEEmcIUSmdCluster cluster;
00584           cluster.add( (*iseed) );
00585           flcount++;
00586 
00587           Int_t ind_max = TMath::Min(287,index+max_extent); 
00588           Int_t ind_min = TMath::Max(0,  index-max_extent); 
00589 
00592           for ( Int_t i=index+1; i<=ind_max; i++ ) {
00594             StEEmcStrip strip=mEEanalysis->strip(sector,plane,i);
00595 
00596             if ( strip.energy()<=0. ) strip.energy(0.0);
00597 
00599             owned[ strip.index() ] = true;
00600             xseed[ strip.index() ] = true;
00602             cluster.add(strip);
00603             if(strip.energy()>0) flcount++;     
00604             
00605           }
00606           for ( Int_t i=index-1; i>=ind_min; i-- ) {
00608             StEEmcStrip strip=mEEanalysis->strip(sector,plane,i);
00609             //strip.energy(newenergy[i]);
00610 
00611             if (strip.energy()<=0.) strip.energy(0.0);
00612 
00614             owned[ strip.index() ] = true;
00615             xseed[ strip.index() ] = true;
00617             cluster.add(strip);
00618             if(strip.energy()>0) flcount++;       
00619             
00620           }
00621          
00623           
00624           if ( cluster.size() >= mMinStrips && flcount>=3 ) {
00625             tclusize->Fill(flcount);
00626            
00627             TmSmdClusters[ sector ][ plane ].push_back(cluster);
00628             TmNumberOfClusters[4+plane]++;
00629            
00630             // disallow strips on either side of the cluster
00631             // from forming seeds
00632             Int_t ns=cluster.numberOfStrips();
00633             Int_t left=999,right=-999;
00634             for ( Int_t is=0;is<ns;is++ )
00635             {
00636                 StEEmcStrip s=cluster.strip(is);
00637                 //printf("strip.energy=%f\n",s.energy()); 
00638                 if ( s.index()<left ) left=s.index();
00639                 if ( s.index()>right ) right=s.index();
00640             }
00641 
00642             for ( Int_t ii=0;ii<=mSuppress;ii++ )
00643               {
00644                 if ( left-ii>=0   ) xseed[left-ii] = true;
00645                 if ( right+ii<288 ) xseed[right+ii]= true;
00646               }
00647 
00648           }
00649           
00650 
00651           iseed++;
00652         }//loop over seeds
00653 
00654 
00655     }}// loop over planes/sectors
00656 
00657   //overlapped SMD cluster area divided here, so that we get better resolution.
00658  
00659   for ( Int_t sector=0; sector<12; sector++ ) 
00660     {
00661       for(Int_t plane=0;plane<=1;plane++)
00662         {
00663 
00664           Int_t para1[288];for (Int_t i=0;i<288;i++) para1[i]=0;
00665           for ( UInt_t fclust=0; fclust<TmSmdClusters[sector][plane].size(); fclust++ ) 
00666             {
00667               StEEmcIUSmdCluster cllu = (TmSmdClusters[sector][plane])[fclust];
00668               Int_t fnums=cllu.numberOfStrips();
00669               for ( Int_t fis=0;fis<fnums;fis++ )
00670                 {
00671                   StEEmcStrip stri1=cllu.strip(fis);
00672                   
00673                   int fstindex=stri1.index();
00674                   para1[fstindex]++;
00675                   
00676                 }
00677             }
00678           int flag2;
00679           float EI,EIII;
00680           Int_t Narr=TmSmdClusters[sector][plane].size();
00681           int flag3[288];for (Int_t i=0;i<Narr;i++) flag3[i]=0;
00682           Int_t flag1[288];for (Int_t i=0;i<288;i++) flag1[i]=0;                
00683           for ( UInt_t iclust=0; iclust<TmSmdClusters[sector][plane].size(); iclust++ ) 
00684             {
00685               flag2=0;
00686               
00687               
00688               flag3[iclust]++;
00689               StEEmcIUSmdCluster cl1 = (TmSmdClusters[sector][plane])[iclust];
00690               Int_t nums1=cl1.numberOfStrips();
00691               UInt_t oindex1=cl1.strip(0).index();
00692 
00693               for ( UInt_t jclust=0; jclust<TmSmdClusters[sector][plane].size(); jclust++ ) 
00694                 {
00695                   StEEmcIUSmdCluster cl2 = (TmSmdClusters[sector][plane])[jclust];
00696                   Int_t nums2=cl2.numberOfStrips();
00697                   UInt_t oindex2=cl2.strip(0).index();
00698                   if(oindex1==oindex2)
00699                     {continue;}
00700                   int seedds=abs(int(oindex1-oindex2));
00701                   cludis->Fill(seedds);
00702                   if(seedds>2*max_extent)
00703                     {continue;}
00704                   flag2++;
00705                   flag3[jclust]++;
00706                   //printf("flag2=%d\n",flag2);
00707                   if(flag3[iclust]>=2 && flag3[jclust]>=2)
00708                     {continue;}
00709                   EI=0;
00710                   EIII=0;
00711                   if(oindex1<oindex2)
00712                     {
00713                       for ( Int_t is=0;is<nums1;is++ )
00714                         {
00715                           StEEmcStrip stri=cl1.strip(is); 
00716                           UInt_t sindex=stri.index();
00717                           if(sindex<=oindex1)
00718                             {
00719                               EI=EI+stri.energy();
00720                             }
00721                         }
00722                       for ( Int_t ist=0;ist<nums2;ist++ )
00723                         {
00724                           StEEmcStrip strit=cl2.strip(ist); 
00725                           UInt_t tindex=strit.index();
00726                           if(tindex>=oindex2)
00727                             {
00728                               EIII=EIII+strit.energy();
00729                             }
00730                         }
00731                     }
00732                   else
00733                     {
00734                       for ( Int_t is=0;is<nums1;is++ )
00735                         {
00736                           StEEmcStrip stri=cl1.strip(is); 
00737                           UInt_t sindex=stri.index();
00738                           if(sindex>=oindex1)
00739                             {
00740                               EIII=EIII+stri.energy();
00741                             }
00742                         }
00743                       for ( Int_t ist=0;ist<nums2;ist++ )
00744                         {
00745                           StEEmcStrip strit=cl2.strip(ist); 
00746                           UInt_t tindex=strit.index();
00747                           if(tindex<=oindex2)
00748                             {
00749                               EI=EI+strit.energy();
00750                             }
00751                         }
00752                    
00753                     }
00754                   
00755                   StEEmcIUSmdCluster tclust1;
00756                   StEEmcIUSmdCluster tclust2;
00757 
00758 
00759                   if(oindex1<oindex2)
00760                     {
00761                       for ( Int_t is=0;is<nums1;is++ )
00762                         {
00763                           StEEmcStrip stri=cl1.strip(is); 
00764                           UInt_t sindex=stri.index();
00765                           if(para1[sindex]==2 )
00766                             {
00767                               float tpe1=stri.energy()*EI/(EI+EIII);
00768                               stri.energy(tpe1);
00769                               tclust1.add(stri);
00770                               flag1[sindex]++;
00771                              
00772                             }
00773                           if(para1[sindex]==1 )
00774                             {
00775                               tclust1.add(stri);
00776                               flag1[sindex]++;
00777                             
00778                             }
00779                         }
00780                       for ( Int_t ist=0;ist<nums2;ist++ )
00781                         {
00782                           StEEmcStrip strit=cl2.strip(ist); 
00783                           UInt_t tindex=strit.index();
00784                           if(para1[tindex]==2 )
00785                             {
00786                               float tpe2=strit.energy()*EIII/(EI+EIII);
00787                               strit.energy(tpe2);
00788                               tclust2.add(strit);
00789                               flag1[tindex]++;
00790                               
00791                             }
00792                           if(para1[tindex]==1 )
00793                             {
00794                               tclust2.add(strit);
00795                               flag1[tindex]++;
00796                               
00797                             }
00798                         }
00799                     } //end of if(oindex1<oindex2)
00800 
00801                   if(oindex1>oindex2)
00802                     {
00803                       for ( Int_t is=0;is<nums1;is++ )
00804                         {
00805                           StEEmcStrip stri=cl1.strip(is); 
00806                           UInt_t sindex=stri.index();
00807                           if(para1[sindex]==2 )
00808                             {
00809                               float tpe2=stri.energy()*EIII/(EI+EIII);
00810                               stri.energy(tpe2);
00811                               tclust2.add(stri);
00812                               flag1[sindex]++;
00813                               
00814                             }
00815                           if(para1[sindex]==1 )
00816                             {
00817                               tclust2.add(stri);
00818                               flag1[sindex]++;
00819                               
00820                             }
00821                         }
00822                       for ( Int_t ist=0;ist<nums2;ist++ )
00823                         {
00824                           StEEmcStrip strit=cl2.strip(ist); 
00825                           UInt_t tindex=strit.index();
00826                           if(para1[tindex]==2 )
00827                             {
00828                               float tpe1=strit.energy()*EI/(EI+EIII);
00829                               strit.energy(tpe1);
00830                               tclust1.add(strit);
00831                               flag1[tindex]++;
00832                               
00833                             }
00834                           if(para1[tindex]==1 )
00835                             {
00836                               tclust1.add(strit);
00837                               flag1[tindex]++;
00838                               
00839                             }
00840                         }
00841                     } //end of if(oindex1>oindex2)
00842 
00843                   // read the new clusters(originally overlapped) into storage
00844                   int seedind1=tclust1.strip(0).index();
00845                   int seedind2=tclust2.strip(0).index();
00846                   if ( flag1[seedind1]<2) 
00847                     {
00848                       tclust1.key( mClusterId++ );
00849                       mSmdClusters[ sector ][ plane ].push_back(tclust1); 
00850                       clusize->Fill(tclust1.size());
00851                       mNumberOfClusters[4+plane]++;
00852                     }
00853                   if ( flag1[seedind2]<2 ) 
00854                     {
00855                       tclust2.key( mClusterId++ );
00856                       mSmdClusters[ sector ][ plane ].push_back(tclust2); 
00857                       clusize->Fill(tclust2.size());
00858                       mNumberOfClusters[4+plane]++;
00859                     }
00860 
00861                 }//end of first round looping clusters
00862 
00863               if(flag2==0)
00864                 {         
00865                       
00866                       cl1.key( mClusterId++ );
00867                       mSmdClusters[ sector ][ plane ].push_back(cl1); 
00868                       clusize->Fill(cl1.size());              
00869                       mNumberOfClusters[4+plane]++;
00870                     
00871                 }
00872             }//endo of second round looping clusters
00873 
00874 
00875 
00876 
00877         }//end of plane
00878       
00879     }//end of sector
00880 
00881 
00882   return true;
00883 
00884 }
00885 
00886 
00887 
00888 // ----------------------------------------------------------------------------
00889 
00890 void StEEmcIUClusterMaker::fillStEvent()
00891 {
00892 
00893   StEvent *stevent=(StEvent*)GetInputDS("StEvent");
00894   if ( !stevent ) {
00895     Warning("fillStEvent","called, but no StEvent to be found");
00896     return;
00897   }
00898 
00899   std::cout << "Adding tower clusters to StEvent at " << stevent << std::endl << std::flush;
00900 
00904   StEmcDetector *detector=stevent->emcCollection()->detector(kEndcapEmcTowerId);
00905   if ( !detector )
00906     {
00907       Warning("fillStEvent","detector == NULL, MAJOR StEvent problem, continuing");
00908       return;
00909     }
00910   
00911 
00919   
00920   if ( mNumberOfClusters[0] > 0 )
00921     {
00922 
00923 
00928       StEmcClusterCollection *collect = detector -> cluster();
00929       if ( !collect ) 
00930         {
00931           //Warning("fillStEvent","StEmcClusterCollection (towers) was NULL, so I'm creating one.");
00932           collect = new StEmcClusterCollection();
00933           detector->setCluster( collect );
00934         }
00935 
00936       assert(collect);
00937       collect->setDetector( kEndcapEmcTowerId );
00938       collect->setClusterFinderId( 123 );
00939       collect->setClusterFinderParamVersion( 123 );//123 or 321?
00940   
00942       for ( Int_t isector=0; isector<12; isector++ )
00943         {
00944 
00946           for ( UInt_t iclust=0; iclust<mTowerClusters[isector][0].size(); iclust++ ) 
00947             {
00948 
00949               StEEmcIUCluster cl=(mTowerClusters[isector][0])[iclust]; 
00950               
00954               StEmcCluster *emccluster = new StEmcCluster();
00955               emccluster->setEta( cl.momentum().Eta() );
00956               emccluster->setPhi( cl.momentum().Phi() );
00957               emccluster->setSigmaEta(-1.);
00958               emccluster->setSigmaPhi(-1.);
00959               emccluster->setEnergy( cl.energy() );
00960               emccluster->SetUniqueID( cl.key() );
00961 #if 1
00962               for ( Int_t i=0; i< cl.numberOfTowers(); i++ ) 
00963                 {
00964                   StEmcRawHit *hit=cl.tower(i).stemc();
00965                   assert( hit );         
00966                   emccluster->addHit( hit );
00967                 }
00968 #endif
00969             
00970               collect->addCluster( emccluster );
00971 
00972               
00974               mEtoEE[ emccluster ] = cl;
00975               cl.stemc( emccluster );
00976 
00977             }
00978           
00979         }
00980 
00981     }
00982 
00983   else // if ( mNumberOfClusters[0] == 0 )
00984     {
00985 
00986       detector->setCluster( NULL );
00987 
00988     }
00989 
00990 
00991 
00992 
00994   //  detector->setCluster( collect );
00995 
00999   detector=stevent->emcCollection()->detector(kEndcapEmcPreShowerId);
01000   if ( !detector )
01001     {
01002       Warning("fillStEvent","detector == NULL for pre/post, no clusters for you");
01003     }
01004   else if ( mNumberOfClusters[1] > 0 ||
01005             mNumberOfClusters[2] > 0 || 
01006             mNumberOfClusters[3] > 0 )
01007     {
01008 
01009       StEmcClusterCollection *pqr = detector -> cluster();
01010       if ( !pqr )
01011         {
01012           //Warning("fillStEvent","StEmcClusterCollection (pre/post) was NULL, so I'm creating one.");
01013           pqr = new StEmcClusterCollection();
01014           detector->setCluster( pqr );
01015         }
01016       assert(pqr); 
01017       pqr -> setDetector( kEndcapEmcPreShowerId );
01018       pqr -> setClusterFinderId( 123 );
01019       pqr -> setClusterFinderParamVersion( 321 );
01020   
01022       for ( Int_t isector=0; isector<12; isector++ )
01023     
01025       for ( Int_t ilayer=1; ilayer<4; ilayer++ ) 
01026     
01028       for ( UInt_t iclust=0; iclust<mTowerClusters[isector][ilayer].size(); iclust++ ) 
01029         {
01030 
01031           StEEmcIUCluster cl=(mTowerClusters[isector][ilayer])[iclust];
01032 
01036           StEmcCluster *emccluster = new StEmcCluster();
01037           emccluster->setEta( cl.momentum().Eta() );
01038           emccluster->setPhi( cl.momentum().Phi() );
01039           emccluster->setSigmaEta(-1.);
01040           emccluster->setSigmaPhi(-1.);
01041           emccluster->setEnergy( cl.energy() );
01042           emccluster->SetUniqueID( cl.key() );
01043 #if 1
01044           for ( Int_t i=0; i< cl.numberOfTowers(); i++ ) 
01045             {
01046               StEmcRawHit *hit=cl.tower(i).stemc();
01047               assert( hit );         
01048               emccluster->addHit( hit );
01049             }
01050 #endif
01051 
01052 
01053           pqr->addCluster( emccluster );
01054 
01055           mEtoEE[ emccluster ] = cl;
01056           cl.stemc( emccluster );
01057 
01058         }
01059 
01060     }
01061   else // if ( mNumberOfClusters[1, 2 OR 3] == 0 )
01062     {
01063 
01064       detector->setCluster( NULL );
01065 
01066     }
01067 
01068   
01072   StDetectorId ids[]={ kEndcapSmdUStripId, kEndcapSmdVStripId };
01073 
01074 
01075   for ( Int_t iplane=0; iplane<2; iplane++ ) 
01076     {
01077 
01078       detector=stevent->emcCollection()->detector(ids[iplane]);
01079       if ( !detector ) 
01080         {
01081           Warning("fillStEvent","detector == NULL for smd plane, no clusters for you");
01082         }
01083       else if ( mNumberOfClusters[4+iplane] > 0 )
01084         {
01085           
01086           StEmcClusterCollection *smdc = detector -> cluster();
01087           if ( !smdc ) 
01088             {
01089               //Warning("fillStEvent","StEmcClusterCollection (smd) was NULL, so I'm creating one.");
01090               smdc = new StEmcClusterCollection();
01091               detector->setCluster( smdc );
01092             }
01093           
01094           smdc->setDetector( ids[iplane] );
01095           smdc->setClusterFinderId( 123 );
01096           smdc->setClusterFinderParamVersion( 321 );
01097 
01098           
01099           for ( Int_t isector=0; isector<12; isector++ ) 
01100             {
01101               
01102               for ( UInt_t iclust=0; iclust<mSmdClusters[isector][iplane].size(); iclust++ ) 
01103                 {
01104                   
01105                   StEEmcIUSmdCluster cl = (mSmdClusters[isector][iplane])[iclust];
01106                   
01107 
01108                   StEmcCluster *emccluster = new StEmcCluster();
01109                   emccluster->setEta( cl.mean() );
01110                   emccluster->setPhi( (Float_t)iplane );
01111                   emccluster->setSigmaEta(-1.);
01112                   emccluster->setSigmaPhi(-1.);
01113                   emccluster->setEnergy( cl.energy() );
01114                   emccluster->SetUniqueID( cl.key() );
01115                   printf("fjaoshjaojgnfaj----------------------------");
01116                   printf("uniq id=%d\n",cl.key());
01117                   for ( Int_t i=0; i< cl.numberOfStrips(); i++ ) 
01118                     {
01119                       StEmcRawHit *hit=cl.strip(i).stemc();
01120                       assert( hit );         
01121                       emccluster->addHit( hit );
01122                     }             
01123                   smdc->addCluster( emccluster );
01124 
01125                   mEtoEEsmd[ emccluster ] = cl;
01126                   cl.stemc( emccluster );
01127                   
01128                 }
01129           
01130               
01131             }
01132           
01133         }
01134 
01135       else
01136         {
01137           detector -> setCluster( NULL );
01138         }
01139 
01140 
01141 
01142     }
01143 
01144 
01145 }
01146 
01147 
01148 
01149 // ----------------------------------------------------------------------------
01150 Bool_t StEEmcIUClusterMaker::verifyStEvent()
01151   {
01152 
01154   StEvent *stevent=(StEvent*)GetInputDS("StEvent");
01155   if ( !stevent ) {
01156     Warning("verifyStEvent","No StEvent found.");
01157     return true;
01158   }
01159   //StEmcCollection *emccollection=stevent->emcCollection();
01160 
01161   Bool_t go = true;
01162  
01163   StEmcDetector *detector=stevent->emcCollection()->detector(kEndcapEmcTowerId);
01164   if ( !detector )
01165     {
01166       Warning("verifyStEvent","detector == NULL for towers");
01167       return false;
01168     }
01169 
01170   
01171   StEmcClusterCollection *cc=detector->cluster();
01172   if ( cc )
01173     {
01174 
01175 
01180       
01181       Float_t emc_sum_towers = 0.;
01182       Float_t eemc_sum_towers = 0.;
01183   
01184       StSPtrVecEmcCluster &emcClusters=cc->clusters();
01185       for ( UInt_t i=0; i<emcClusters.size(); i++ ) 
01186         {
01187           StEmcCluster *cl=emcClusters[i];
01188           assert(cl);
01189           emc_sum_towers += cl->energy();
01190         }
01191       
01192       for ( Int_t sec=0; sec<12; sec++ ) 
01193         {
01194           
01195           for ( Int_t i=0; i<numberOfClusters(sec,0); i++ )
01196             eemc_sum_towers += cluster(sec,0,i).energy();
01197         }
01198       
01199       std::cout << "StEEmcIUClusterMaker tower checksum: ";
01200       if ( emc_sum_towers == eemc_sum_towers ) {
01201         std::cout << "passed";
01202       }
01203       else {
01204         std::cout << "FAILED"; go=false;
01205       }
01206       std::cout << std::endl;
01207 
01208     }
01209   else {
01210 
01211     std::cout << "StEEmcIUClusterMaker tower checksum: NULL collection, nclust=" << mNumberOfClusters[0] << std::endl;
01212     go &= (mNumberOfClusters[0]==0);
01213 
01214   }
01215 
01220   Float_t emc_sum_smdu = 0.;
01221   Float_t eemc_sum_smdu = 0.;
01222 
01223   detector=stevent->emcCollection()->detector(kEndcapSmdUStripId);
01224   if ( !detector )
01225     {
01226       Warning("verifyStEvent","detector == NULL for smdu");
01227       return false;
01228     }
01229 
01230   cc=detector->cluster();
01231   if ( cc ) 
01232     {
01233 
01234       StSPtrVecEmcCluster &smduClusters=cc->clusters();
01235 
01236       for ( UInt_t i=0; i<smduClusters.size(); i++ ) 
01237         {
01238           StEmcCluster *cl=smduClusters[i];
01239           assert(cl);
01240           emc_sum_smdu += cl->energy();
01241         }
01242       
01243   
01244       for ( Int_t sec=0; sec<12; sec++ ) 
01245         {
01246 
01247           for ( Int_t i=0; i<numberOfSmdClusters(sec,0); i++ )
01248         
01249             {
01250               eemc_sum_smdu += smdcluster(sec,0,i).energy();
01251 
01252             }
01253 
01254         }
01255 
01256       std::cout << "StEEmcIUClusterMaker smdu checksum: ";
01257       if ( emc_sum_smdu == eemc_sum_smdu ) {
01258         std::cout << "passed";
01259       }
01260       else {
01261         std::cout << "FAILED"; go=false;
01262       }
01263       std::cout << std::endl;
01264 
01265     }  
01266   else 
01267     {
01268       std::cout << "StEEmcIUClusterMaker smdu checksum: NULL collection, nclust=" << mNumberOfClusters[4] << std::endl;
01269       go &= (mNumberOfClusters[4]==0);
01270     }
01271 
01272 
01274 
01275   Float_t emc_sum_smdv = 0.;
01276   Float_t eemc_sum_smdv = 0.;
01277 
01278   detector=stevent->emcCollection()->detector(kEndcapSmdVStripId);
01279   if (!detector)
01280     {
01281       Warning("verifyStEvent","detector == NULL for smdv");
01282       return false;
01283     }
01284 
01285   cc=detector->cluster();
01286 
01287   if ( cc )
01288     {
01289 
01290       StSPtrVecEmcCluster &smdvClusters=cc->clusters();
01291       
01292       for ( UInt_t i=0; i<smdvClusters.size(); i++ ) 
01293         {
01294           StEmcCluster *cl=smdvClusters[i];
01295           assert(cl);
01296           emc_sum_smdv += cl->energy();
01297         }
01298       
01299       
01300       for ( Int_t sec=0; sec<12; sec++ ) 
01301         {
01302           
01303           for ( Int_t i=0; i<numberOfSmdClusters(sec,1); i++ )
01304             
01305             {
01306               eemc_sum_smdv += smdcluster(sec,1,i).energy();
01307               
01308             }
01309           
01310         }
01311       
01312       std::cout << "StEEmcIUClusterMaker smdv checksum: ";
01313       if ( emc_sum_smdv == eemc_sum_smdv ) {
01314         std::cout << "passed";
01315       }
01316       else {
01317         std::cout << "FAILED"; go=false;
01318       }
01319       std::cout << std::endl;
01320 
01321     }
01322   else 
01323     {
01324       std::cout << "StEEmcIUClusterMaker smdv checksum: NULL collection, nclust=" << mNumberOfClusters[5] << std::endl;
01325       go &= (mNumberOfClusters[5]==0);
01326     }
01327   
01328 
01329   
01330   return go;
01331 }
01332 
01333 
01334 
01335 // ----------------------------------------------------------------------------
01336 void StEEmcIUClusterMaker::print()
01337 {
01338 
01339   std::cout << "StEEmcIUClusterMaker::print()" << std::endl;
01340   const Char_t *names[] = { "tower", "pre1", "pre2", "post", "smdu", "smdv" };
01341   for ( Int_t i=0;i<6;i++ )
01342     {
01343 
01344       std::cout << "Number of " << names[i] 
01345                 << " clusters = " << mNumberOfClusters[i] 
01346                 << std::endl;
01347 
01348     }
01349 
01350   std::cout << "printout of tower clusters follows:" << std::endl;
01351   for ( Int_t sec=0;sec<12;sec++){
01352   for ( Int_t i=0;i<numberOfSmdClusters(sec,0);i++ )
01353     {
01354       StEEmcIUSmdCluster clust=(mSmdClusters[sec][0])[i];
01355       clust.print();
01356       //std::cout << "cluster.key()=" << clust.key() << std::endl;
01357       //      std::cout << "cluster.eta()=" << clust.momentum().Eta() << std::endl;
01358       //      std::cout << "cluster.phi()=" << clust.momentum().Phi() << std::endl;      
01359       //      std::cout << "cluster.energy()=" << clust.energy() << std::endl;
01360     }
01361  for ( Int_t i=0;i<numberOfSmdClusters(sec,1);i++ )
01362     {
01363       StEEmcIUSmdCluster clust=(mSmdClusters[sec][1])[i];
01364       clust.print();
01365     }
01366   }
01367 
01368 }