StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StEEmcIUClusterMaker.cxx
1 
50 #include "StEEmcIUClusterMaker.h"
51 
52 #include <algorithm>
53 #include <iostream>
54 
55 #include "TMath.h"
56 
57 /* StEvent stuff */
58 #include "StEvent/StEvent.h"
59 #include "StEvent/StEmcCollection.h"
60 #include "StEvent/StEmcDetector.h"
61 #include "StEvent/StEmcModule.h"
62 #include "StEvent/StEmcClusterCollection.h"
63 #include "StEvent/StEmcCluster.h"
64 
65 #define DEBUG 0
66 //#define MonteCarloS
67 
68 ClassImp(StEEmcIUClusterMaker);
69 
70 // ----------------------------------------------------------------------------
72 {
73 
75  mFillStEvent = 0;
76 
78  mSuppress = 0;
79 
80  mLoose=false;
81  mSkip=true;
82  //Meanclust=0;
85  mClusterId = 0;
86 
88  const Float_t eseeds[] = { 0.6, 1.0/1000, 1.0/1000, 1.0/1000., 0.3/1000., 0.3/1000. };
89  for ( Int_t i = 0; i < 6; i++ ) seedEnergy( eseeds[i], i );
90 
91  mMaxExtent = 3;
92  mSeedFloor = 2.0;
93 
106 
107  StEEmcIUClusterVec_t t;
108  std::vector< StEEmcIUClusterVec_t > layers;
109  for ( Int_t i = 0; i < 4; i++ ) layers.push_back(t);
110  for ( Int_t i = 0; i < 12; i++ ) mTowerClusters.push_back(layers);
111 
112  StEEmcIUSmdClusterVec_t s;
113  std::vector< StEEmcIUSmdClusterVec_t > planes;
114  planes.push_back(s);
115  planes.push_back(s);
116  for ( Int_t i = 0; i < 12; i++ ) mSmdClusters.push_back(planes);
117  for ( Int_t i = 0; i < 12; i++ ) TmSmdClusters.push_back(planes);
118 
119  StEEmcTowerVec_t tow;
120  std::vector< StEEmcTowerVec_t > lay;
121  for ( Int_t i = 0; i < 4; i++ ) lay.push_back(tow);
122  for ( Int_t i = 0; i < 12; i++ ) mSeedTowers.push_back(lay);
123 
124  mEEtow=new EEmcGeomSimple();
125  mEEsmd=EEmcSmdGeom::instance();
126  mEEmap=EEmcSmdMap::instance();
127 
128  mMinStrips=3;
129 
130 }
131 
132 // ----------------------------------------------------------------------------
134 {
136  assert(mEEanalysis);
137  clusize=new TH1D("clustersize", "smd cluster size ",10,0,10);
138  tclusize=new TH1D("tclustersize", "smd cluster size ",10,0,10);
139  cludis=new TH1D("cludis", "cluster distance ",20,0,20);
140  return StMaker::Init();
141 }
142 
143 // ----------------------------------------------------------------------------
145 {
146 
148  if ( !buildTowerClusters() ) return kStWarn;
149 
151  if ( !buildSmdClusters() ) return kStWarn;
152  //fillStEvent();
153 
155  if ( mFillStEvent ) fillStEvent();
156 
158  if ( mFillStEvent )
159  if ( !verifyStEvent() ) Warning("Make","StEvent not properly copied");
160  //print();
161 
162 
163  return kStOK;
164 }
165 
166 // ----------------------------------------------------------------------------
167 void StEEmcIUClusterMaker::Clear( Option_t *opts )
168 {
169 
171  for ( Int_t sector=0; sector<12; sector++ ) {
172  for ( Int_t layer=0; layer<4; layer++ )
173  mTowerClusters[sector][layer].clear();
174  for ( Int_t plane=0; plane<2; plane++ )
175  {
176  mSmdClusters[sector][plane].clear();
177  TmSmdClusters[sector][plane].clear();
178  }
179  }
180 
181  for ( Int_t i=0;i<6;i++ )
182  {
183  mNumberOfClusters[i]=0;
184  TmNumberOfClusters[i]=0;
185  }
186  mClusterId = 0;
187 
188  return;
189 }
190 
191 
192 // ----------------------------------------------------------------------------
194 {
195 
196  // This part is not important in the reconstruction of pi0. It only provides some fundamental checks and basic idea about the tower clusters.
199  //static const Int_t layer=0;
200 
201 
203 
204  eeen=0;
205  ep1=0;
206  ep2=0;
207  ep3=0;
208  for ( Int_t layer=0;layer<4;layer++)
209  {
210 
213  Float_t weights[720]; for ( Int_t i=0;i<720;i++ ) weights[i]=0.;
214 
216  StEEmcIUClusterVec_t myClusters;
217 
218 
220  StEEmcTowerVec_t towers = mEEanalysis -> towers( layer );
221 
222 
223 
225  std::sort(towers.begin(),towers.end());
227  std::reverse(towers.begin(),towers.end());
228 
229 
232  StEEmcTowerVec_t::iterator last = towers.begin();
233  while ( last != towers.end() ) {
234  if(layer==0){
235  eeen+=(*last).energy();
236  }
237  if(layer==1){
238  ep1+=(*last).energy();
239  }
240  if(layer==2){
241  ep2+=(*last).energy();
242  }
243  if(layer==3){
244  ep3+=(*last).energy();
245  }
246  if ( (*last).energy() < mSeedEnergy[layer] ) break;
251  //$$$ mSeedTowers[layer][ (*last).sector() ].push_back((*last));
252 #if DEBUG
253  std::cout << "-- Seed tower ----------------------" << std::endl;
254  (*last).print();
255 #endif
256  last++;
257  }
258  //std::cout<<"tower energy total = "<< eeen<< std::endl;
259 
260  StEEmcTowerVec_t::iterator iter = towers.begin();
261  while ( iter != towers.end() ) {
262 
264  if ( iter==last ) break;
265 
266 
273  if ( weights[ (*iter).index() ] > 0. ) {
274  iter++;
275  continue;
276  }
277 
280  for ( Int_t in=0; in < (*iter).numberOfNeighbors(); in++ ) {
281  StEEmcTower t=(*iter).neighbor(in);
282  weights[ t.index() ] += (*iter).energy();
283  }
284 
285 
286  iter++;
287 
288  }// loop over towers to init weights
289 
290 
292  iter=towers.begin();
293  while ( iter != towers.end() ) {
294 
296  if ( iter==last ) break;
297 
301  if ( weights[ (*iter).index() ] > 0. ) {
302  iter++;
303  continue;
304  }
305 
307  StEEmcTower seed=(*iter);
308 #if DEBUG
309  std::cout << "--- Clustering ----------------" << std::endl;
310  seed.print();
311 #endif
312 
313  TVector3 momentum;
314  cluster.add(seed,1.0);
315  UInt_t sec,sub,eta;
316  sec=(UInt_t)seed.sector();
317  sub=(UInt_t)seed.subsector();
318  eta=(UInt_t)seed.etabin();
319  TVector3 d=mEEtow->getTowerCenter(sec,sub,eta).Unit();
320  momentum += ( seed.energy() * d );
321 
322  for ( Int_t in=0; in<seed.numberOfNeighbors(); in++ ) {
323  StEEmcTower t=seed.neighbor(in);
324  sec=(UInt_t)t.sector();
325  sub=(UInt_t)t.subsector();
326  eta=(UInt_t)t.etabin();
327  d=mEEtow->getTowerCenter(sec,sub,eta).Unit();
328  Float_t weight = seed.energy() / weights[ t.index() ];
329  momentum += ( t.energy() * d * weight );
330  cluster.add(t, weight);
331 #if DEBUG
332  std::cout << "adding " << t.name() << " E=" << t.energy() << " W=" << weight << std::endl;
333 #endif
334  }
335 
336 
338  cluster.momentum( momentum );
339  cluster.key( mClusterId++ );
340 #if DEBUG
341  cluster.print();
342 #endif
343 
344  mTowerClusters[ seed.sector() ][ layer ].push_back( cluster );
345  mNumberOfClusters[layer]++;
346 
347  iter++;
348 
349  }
350 
351  }// end of eventual loop over layers
352 
353 
354  return true;
355 
356 }
357 
358 // ----------------------------------------------------------------------------
360 {
361 
362  //this part is essential in the reconstruction of pi0 in the EEMC.
363 
364  Int_t max_extent = mMaxExtent;
365  countUseed=0;
366  countVseed=0;
367  uswidth=0;
368  vswidth=0;
369  uamp=0;
370  vamp=0;
371  esmdu=0;
372  esmdv=0;
373 
374 
376  for ( Int_t sector=0; sector<12; sector++ ) {
378 
379  for ( Int_t plane=0; plane<2; plane++ ) {
380  float ampe=0;
382  Float_t floor[288]; for ( Int_t i=0; i<288; i++ ) floor[i]=0.;
383 
385  //Float_t energy[288]; for ( Int_t i=0; i<288; i++ ) energy[i]=0.;
386  //averaging energies
387  //Float_t newenergy[288]; for ( Int_t i=0; i<288; i++ ) newenergy[i]=0.;
389  StEEmcStripVec_t strips=mEEanalysis->strips(sector,plane);
390  StEEmcStripVec_t seeds;
392  std::sort(strips.begin(),strips.end());
394  std::reverse(strips.begin(),strips.end());
395 
397  StEEmcStripVec_t::iterator istrip=strips.begin();
398 
399  while ( istrip != strips.end() ) {
400  if ( (*istrip).stat()||(*istrip).fail() ) {
401  istrip++;
402  continue;
403  }
404  ampe=(*(strips.begin())).energy();
405  //energy[ (*istrip).index() ] = (*istrip).energy();
406  //cout<<"energy="<<(*istrip).energy()<<endl;
407 
408  if((*istrip).energy()>=0)
409  {
410  //fprintf(fout2,"data%02d%d %d %f\n",sector,plane,(*istrip).index(),(*istrip).energy());
411  if(plane==0)
412  {
413  esmdu+=(*istrip).energy();
414  uswidth++;
415  uamp=ampe;
416  }
417  if(plane==1)
418  {
419  esmdv+=(*istrip).energy();
420  vswidth++;
421  vamp=ampe;
422  }
423  }
424  istrip++;
425  }
426 
427  //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.
428  //istrip=strips.begin();
429  //for(int i=2;i<=283;i++){
430  //newenergy[i]=energy[i-1]*0.25+energy[i]*0.5+energy[i+1]*0.25;
431  //}
432  //while ( istrip != strips.end() ) {
433  //if ( (*istrip).stat()||(*istrip).fail() ) {
434  //istrip++;
435  //continue;
436  //}
437  //(*istrip).energy(newenergy[(*istrip).index()]);
438 
439  //istrip++;
440  //}
441  //end of averaging algo
442 
444  std::vector<Bool_t> seed_flags( strips.size(), false );
445 
446 
448  Int_t nstrip=0;
449  Int_t nseeds=0;
450  istrip=strips.begin();
451  while ( istrip!=strips.end() ) {
452 
453  Int_t index=(*istrip).index();
454  Float_t eseed=(*istrip).energy();
455 
457  if ( index <= 3 || index >= 283 ) {
458  istrip++;
459  nstrip++;
460  continue;
461  }
462 
464  if ( mSkip )
465  if ( (*istrip).fail() ) {
466  istrip++;
467  nstrip++;
468  continue;
469  }
470 
471 
474  if ( eseed < (mSeedFloor*floor[ index ] + mSeedEnergy[4+plane]) ) {
475  istrip++;
476  nstrip++;
477  continue;
478  }
479  //if(energy[index+1]<floor[index+1] && energy[index-1]<floor[index-1]){
480  //istrip++;
481  //nstrip++;
482  //continue;
483  //}
485  //cout<<"succeed!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1"<<endl;
486  seeds.push_back( (*istrip) );
487  if(plane==0){countUseed++;}
488  if(plane==1){countVseed++;}
489  nstrip++;
490  nseeds++;
491 
492  // 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.
493  float floor_para1=0.0;
494  float floor_para2=0.0;
495  if(sector==0||sector==3||sector==6||sector==9)
496  {
497  if(plane==0)
498  {
499  floor_para1=0.4;
500  floor_para2=0.2;
501  }
502  else{
503  floor_para1=0.2;
504  floor_para2=0.1;
505  }
506  }
507  if(sector==2 || sector ==1 || sector ==5 || sector==4 || sector==8 || sector==7 || sector==11 || sector==10)
508  {
509  if(plane==1)
510  {
511  floor_para1=0.4;
512  floor_para2=0.2;
513  }
514  else{
515  floor_para1=0.2;
516  floor_para2=0.1;
517  }
518  }
519 #ifdef MonteCarloS
520  floor_para1=0.2;
521  floor_para2=0.1;
522 #endif
523 #ifndef LOOSE_CUTS
529  for ( Int_t i=0; i < 288; i++ )
531  {
532  Int_t dx=TMath::Abs(index-i);
534  if(dx<3) {floor[i]+=1.0*eseed;}
536  if(dx>=3&&dx<5){floor[i]+=floor_para1*eseed;}
538  if(dx>=5&&dx<11){floor[i]+=floor_para2*eseed;}
540  if(dx>=11&&dx<21){floor[i]+=0.05*eseed;}
541  }
542 #else
543  for ( Int_t i=0; i < 288; i++ ) {
545 
546  Int_t dx=TMath::Abs(index-i);
547 
549  if ( dx<7)
550  if ( 0.05 * eseed > floor[i] ) floor[i]+= 0.05 * eseed;
551 
552  }
553 #endif
554  istrip++;
555 
556  }// found all seed strips in that plane of that sector
557 
563  Bool_t owned[288]; for (Int_t i=0;i<288;i++) owned[i]=false;
564  Bool_t xseed[288]; for (Int_t i=0;i<288;i++) xseed[i]=false;
565 
566 
567  StEEmcStripVec_t::iterator iseed=seeds.begin();
568  while ( iseed != seeds.end() ) {
569 
570  int flcount=0;
572  Int_t index=(*iseed).index();
573  if ( owned[index] || (mSuppress&&xseed[index]) ) {
574  iseed++;
575  continue;
576  }
577 
579  owned[index]=true;
580  xseed[index]=true;
584  cluster.add( (*iseed) );
585  flcount++;
586 
587  Int_t ind_max = TMath::Min(287,index+max_extent);
588  Int_t ind_min = TMath::Max(0, index-max_extent);
589 
592  for ( Int_t i=index+1; i<=ind_max; i++ ) {
594  StEEmcStrip strip=mEEanalysis->strip(sector,plane,i);
595 
596  if ( strip.energy()<=0. ) strip.energy(0.0);
597 
599  owned[ strip.index() ] = true;
600  xseed[ strip.index() ] = true;
602  cluster.add(strip);
603  if(strip.energy()>0) flcount++;
604 
605  }
606  for ( Int_t i=index-1; i>=ind_min; i-- ) {
608  StEEmcStrip strip=mEEanalysis->strip(sector,plane,i);
609  //strip.energy(newenergy[i]);
610 
611  if (strip.energy()<=0.) strip.energy(0.0);
612 
614  owned[ strip.index() ] = true;
615  xseed[ strip.index() ] = true;
617  cluster.add(strip);
618  if(strip.energy()>0) flcount++;
619 
620  }
621 
623 
624  if ( cluster.size() >= mMinStrips && flcount>=3 ) {
625  tclusize->Fill(flcount);
626 
627  TmSmdClusters[ sector ][ plane ].push_back(cluster);
628  TmNumberOfClusters[4+plane]++;
629 
630  // disallow strips on either side of the cluster
631  // from forming seeds
632  Int_t ns=cluster.numberOfStrips();
633  Int_t left=999,right=-999;
634  for ( Int_t is=0;is<ns;is++ )
635  {
636  StEEmcStrip s=cluster.strip(is);
637  //printf("strip.energy=%f\n",s.energy());
638  if ( s.index()<left ) left=s.index();
639  if ( s.index()>right ) right=s.index();
640  }
641 
642  for ( Int_t ii=0;ii<=mSuppress;ii++ )
643  {
644  if ( left-ii>=0 ) xseed[left-ii] = true;
645  if ( right+ii<288 ) xseed[right+ii]= true;
646  }
647 
648  }
649 
650 
651  iseed++;
652  }//loop over seeds
653 
654 
655  }}// loop over planes/sectors
656 
657  //overlapped SMD cluster area divided here, so that we get better resolution.
658 
659  for ( Int_t sector=0; sector<12; sector++ )
660  {
661  for(Int_t plane=0;plane<=1;plane++)
662  {
663 
664  Int_t para1[288];for (Int_t i=0;i<288;i++) para1[i]=0;
665  for ( UInt_t fclust=0; fclust<TmSmdClusters[sector][plane].size(); fclust++ )
666  {
667  StEEmcIUSmdCluster cllu = (TmSmdClusters[sector][plane])[fclust];
668  Int_t fnums=cllu.numberOfStrips();
669  for ( Int_t fis=0;fis<fnums;fis++ )
670  {
671  StEEmcStrip stri1=cllu.strip(fis);
672 
673  int fstindex=stri1.index();
674  para1[fstindex]++;
675 
676  }
677  }
678  int flag2;
679  float EI,EIII;
680  Int_t Narr=TmSmdClusters[sector][plane].size();
681  int flag3[288];for (Int_t i=0;i<Narr;i++) flag3[i]=0;
682  Int_t flag1[288];for (Int_t i=0;i<288;i++) flag1[i]=0;
683  for ( UInt_t iclust=0; iclust<TmSmdClusters[sector][plane].size(); iclust++ )
684  {
685  flag2=0;
686 
687 
688  flag3[iclust]++;
689  StEEmcIUSmdCluster cl1 = (TmSmdClusters[sector][plane])[iclust];
690  Int_t nums1=cl1.numberOfStrips();
691  UInt_t oindex1=cl1.strip(0).index();
692 
693  for ( UInt_t jclust=0; jclust<TmSmdClusters[sector][plane].size(); jclust++ )
694  {
695  StEEmcIUSmdCluster cl2 = (TmSmdClusters[sector][plane])[jclust];
696  Int_t nums2=cl2.numberOfStrips();
697  UInt_t oindex2=cl2.strip(0).index();
698  if(oindex1==oindex2)
699  {continue;}
700  int seedds=abs(int(oindex1-oindex2));
701  cludis->Fill(seedds);
702  if(seedds>2*max_extent)
703  {continue;}
704  flag2++;
705  flag3[jclust]++;
706  //printf("flag2=%d\n",flag2);
707  if(flag3[iclust]>=2 && flag3[jclust]>=2)
708  {continue;}
709  EI=0;
710  EIII=0;
711  if(oindex1<oindex2)
712  {
713  for ( Int_t is=0;is<nums1;is++ )
714  {
715  StEEmcStrip stri=cl1.strip(is);
716  UInt_t sindex=stri.index();
717  if(sindex<=oindex1)
718  {
719  EI=EI+stri.energy();
720  }
721  }
722  for ( Int_t ist=0;ist<nums2;ist++ )
723  {
724  StEEmcStrip strit=cl2.strip(ist);
725  UInt_t tindex=strit.index();
726  if(tindex>=oindex2)
727  {
728  EIII=EIII+strit.energy();
729  }
730  }
731  }
732  else
733  {
734  for ( Int_t is=0;is<nums1;is++ )
735  {
736  StEEmcStrip stri=cl1.strip(is);
737  UInt_t sindex=stri.index();
738  if(sindex>=oindex1)
739  {
740  EIII=EIII+stri.energy();
741  }
742  }
743  for ( Int_t ist=0;ist<nums2;ist++ )
744  {
745  StEEmcStrip strit=cl2.strip(ist);
746  UInt_t tindex=strit.index();
747  if(tindex<=oindex2)
748  {
749  EI=EI+strit.energy();
750  }
751  }
752 
753  }
754 
755  StEEmcIUSmdCluster tclust1;
756  StEEmcIUSmdCluster tclust2;
757 
758 
759  if(oindex1<oindex2)
760  {
761  for ( Int_t is=0;is<nums1;is++ )
762  {
763  StEEmcStrip stri=cl1.strip(is);
764  UInt_t sindex=stri.index();
765  if(para1[sindex]==2 )
766  {
767  float tpe1=stri.energy()*EI/(EI+EIII);
768  stri.energy(tpe1);
769  tclust1.add(stri);
770  flag1[sindex]++;
771 
772  }
773  if(para1[sindex]==1 )
774  {
775  tclust1.add(stri);
776  flag1[sindex]++;
777 
778  }
779  }
780  for ( Int_t ist=0;ist<nums2;ist++ )
781  {
782  StEEmcStrip strit=cl2.strip(ist);
783  UInt_t tindex=strit.index();
784  if(para1[tindex]==2 )
785  {
786  float tpe2=strit.energy()*EIII/(EI+EIII);
787  strit.energy(tpe2);
788  tclust2.add(strit);
789  flag1[tindex]++;
790 
791  }
792  if(para1[tindex]==1 )
793  {
794  tclust2.add(strit);
795  flag1[tindex]++;
796 
797  }
798  }
799  } //end of if(oindex1<oindex2)
800 
801  if(oindex1>oindex2)
802  {
803  for ( Int_t is=0;is<nums1;is++ )
804  {
805  StEEmcStrip stri=cl1.strip(is);
806  UInt_t sindex=stri.index();
807  if(para1[sindex]==2 )
808  {
809  float tpe2=stri.energy()*EIII/(EI+EIII);
810  stri.energy(tpe2);
811  tclust2.add(stri);
812  flag1[sindex]++;
813 
814  }
815  if(para1[sindex]==1 )
816  {
817  tclust2.add(stri);
818  flag1[sindex]++;
819 
820  }
821  }
822  for ( Int_t ist=0;ist<nums2;ist++ )
823  {
824  StEEmcStrip strit=cl2.strip(ist);
825  UInt_t tindex=strit.index();
826  if(para1[tindex]==2 )
827  {
828  float tpe1=strit.energy()*EI/(EI+EIII);
829  strit.energy(tpe1);
830  tclust1.add(strit);
831  flag1[tindex]++;
832 
833  }
834  if(para1[tindex]==1 )
835  {
836  tclust1.add(strit);
837  flag1[tindex]++;
838 
839  }
840  }
841  } //end of if(oindex1>oindex2)
842 
843  // read the new clusters(originally overlapped) into storage
844  int seedind1=tclust1.strip(0).index();
845  int seedind2=tclust2.strip(0).index();
846  if ( flag1[seedind1]<2)
847  {
848  tclust1.key( mClusterId++ );
849  mSmdClusters[ sector ][ plane ].push_back(tclust1);
850  clusize->Fill(tclust1.size());
851  mNumberOfClusters[4+plane]++;
852  }
853  if ( flag1[seedind2]<2 )
854  {
855  tclust2.key( mClusterId++ );
856  mSmdClusters[ sector ][ plane ].push_back(tclust2);
857  clusize->Fill(tclust2.size());
858  mNumberOfClusters[4+plane]++;
859  }
860 
861  }//end of first round looping clusters
862 
863  if(flag2==0)
864  {
865 
866  cl1.key( mClusterId++ );
867  mSmdClusters[ sector ][ plane ].push_back(cl1);
868  clusize->Fill(cl1.size());
869  mNumberOfClusters[4+plane]++;
870 
871  }
872  }//endo of second round looping clusters
873 
874 
875 
876 
877  }//end of plane
878 
879  }//end of sector
880 
881 
882  return true;
883 
884 }
885 
886 
887 
888 // ----------------------------------------------------------------------------
889 
891 {
892 
893  StEvent *stevent=(StEvent*)GetInputDS("StEvent");
894  if ( !stevent ) {
895  Warning("fillStEvent","called, but no StEvent to be found");
896  return;
897  }
898 
899  std::cout << "Adding tower clusters to StEvent at " << stevent << std::endl;
900 
904  StEmcDetector *detector=stevent->emcCollection()->detector(kEndcapEmcTowerId);
905  if ( !detector )
906  {
907  Warning("fillStEvent","detector == NULL, MAJOR StEvent problem, continuing");
908  return;
909  }
910 
911 
919 
920  if ( mNumberOfClusters[0] > 0 )
921  {
922 
923 
928  StEmcClusterCollection *collect = detector -> cluster();
929  if ( !collect )
930  {
931  //Warning("fillStEvent","StEmcClusterCollection (towers) was NULL, so I'm creating one.");
932  collect = new StEmcClusterCollection();
933  detector->setCluster( collect );
934  }
935 
936  assert(collect);
937  collect->setDetector( kEndcapEmcTowerId );
938  collect->setClusterFinderId( 123 );
939  collect->setClusterFinderParamVersion( 123 );//123 or 321?
940 
942  for ( Int_t isector=0; isector<12; isector++ )
943  {
944 
946  for ( UInt_t iclust=0; iclust<mTowerClusters[isector][0].size(); iclust++ )
947  {
948 
949  StEEmcIUCluster cl=(mTowerClusters[isector][0])[iclust];
950 
954  StEmcCluster *emccluster = new StEmcCluster();
955  emccluster->setEta( cl.momentum().Eta() );
956  emccluster->setPhi( cl.momentum().Phi() );
957  emccluster->setSigmaEta(-1.);
958  emccluster->setSigmaPhi(-1.);
959  emccluster->setEnergy( cl.energy() );
960  emccluster->SetUniqueID( cl.key() );
961 #if 1
962  for ( Int_t i=0; i< cl.numberOfTowers(); i++ )
963  {
964  StEmcRawHit *hit=cl.tower(i).stemc();
965  assert( hit );
966  emccluster->addHit( hit );
967  }
968 #endif
969 
970  collect->addCluster( emccluster );
971 
972 
974  mEtoEE[ emccluster ] = cl;
975  cl.stemc( emccluster );
976 
977  }
978 
979  }
980 
981  }
982 
983  else // if ( mNumberOfClusters[0] == 0 )
984  {
985 
986  detector->setCluster( NULL );
987 
988  }
989 
990 
991 
992 
994  // detector->setCluster( collect );
995 
999  detector=stevent->emcCollection()->detector(kEndcapEmcPreShowerId);
1000  if ( !detector )
1001  {
1002  Warning("fillStEvent","detector == NULL for pre/post, no clusters for you");
1003  }
1004  else if ( mNumberOfClusters[1] > 0 ||
1005  mNumberOfClusters[2] > 0 ||
1006  mNumberOfClusters[3] > 0 )
1007  {
1008 
1009  StEmcClusterCollection *pqr = detector -> cluster();
1010  if ( !pqr )
1011  {
1012  //Warning("fillStEvent","StEmcClusterCollection (pre/post) was NULL, so I'm creating one.");
1013  pqr = new StEmcClusterCollection();
1014  detector->setCluster( pqr );
1015  }
1016  assert(pqr);
1017  pqr -> setDetector( kEndcapEmcPreShowerId );
1018  pqr -> setClusterFinderId( 123 );
1019  pqr -> setClusterFinderParamVersion( 321 );
1020 
1022  for ( Int_t isector=0; isector<12; isector++ )
1023 
1025  for ( Int_t ilayer=1; ilayer<4; ilayer++ )
1026 
1028  for ( UInt_t iclust=0; iclust<mTowerClusters[isector][ilayer].size(); iclust++ )
1029  {
1030 
1031  StEEmcIUCluster cl=(mTowerClusters[isector][ilayer])[iclust];
1032 
1036  StEmcCluster *emccluster = new StEmcCluster();
1037  emccluster->setEta( cl.momentum().Eta() );
1038  emccluster->setPhi( cl.momentum().Phi() );
1039  emccluster->setSigmaEta(-1.);
1040  emccluster->setSigmaPhi(-1.);
1041  emccluster->setEnergy( cl.energy() );
1042  emccluster->SetUniqueID( cl.key() );
1043 #if 1
1044  for ( Int_t i=0; i< cl.numberOfTowers(); i++ )
1045  {
1046  StEmcRawHit *hit=cl.tower(i).stemc();
1047  assert( hit );
1048  emccluster->addHit( hit );
1049  }
1050 #endif
1051 
1052 
1053  pqr->addCluster( emccluster );
1054 
1055  mEtoEE[ emccluster ] = cl;
1056  cl.stemc( emccluster );
1057 
1058  }
1059 
1060  }
1061  else // if ( mNumberOfClusters[1, 2 OR 3] == 0 )
1062  {
1063 
1064  detector->setCluster( NULL );
1065 
1066  }
1067 
1068 
1072  StDetectorId ids[]={ kEndcapSmdUStripId, kEndcapSmdVStripId };
1073 
1074 
1075  for ( Int_t iplane=0; iplane<2; iplane++ )
1076  {
1077 
1078  detector=stevent->emcCollection()->detector(ids[iplane]);
1079  if ( !detector )
1080  {
1081  Warning("fillStEvent","detector == NULL for smd plane, no clusters for you");
1082  }
1083  else if ( mNumberOfClusters[4+iplane] > 0 )
1084  {
1085 
1086  StEmcClusterCollection *smdc = detector -> cluster();
1087  if ( !smdc )
1088  {
1089  //Warning("fillStEvent","StEmcClusterCollection (smd) was NULL, so I'm creating one.");
1090  smdc = new StEmcClusterCollection();
1091  detector->setCluster( smdc );
1092  }
1093 
1094  smdc->setDetector( ids[iplane] );
1095  smdc->setClusterFinderId( 123 );
1096  smdc->setClusterFinderParamVersion( 321 );
1097 
1098 
1099  for ( Int_t isector=0; isector<12; isector++ )
1100  {
1101 
1102  for ( UInt_t iclust=0; iclust<mSmdClusters[isector][iplane].size(); iclust++ )
1103  {
1104 
1105  StEEmcIUSmdCluster cl = (mSmdClusters[isector][iplane])[iclust];
1106 
1107 
1108  StEmcCluster *emccluster = new StEmcCluster();
1109  emccluster->setEta( cl.mean() );
1110  emccluster->setPhi( (Float_t)iplane );
1111  emccluster->setSigmaEta(-1.);
1112  emccluster->setSigmaPhi(-1.);
1113  emccluster->setEnergy( cl.energy() );
1114  emccluster->SetUniqueID( cl.key() );
1115  printf("fjaoshjaojgnfaj----------------------------");
1116  printf("uniq id=%d\n",cl.key());
1117  for ( Int_t i=0; i< cl.numberOfStrips(); i++ )
1118  {
1119  StEmcRawHit *hit=cl.strip(i).stemc();
1120  assert( hit );
1121  emccluster->addHit( hit );
1122  }
1123  smdc->addCluster( emccluster );
1124 
1125  mEtoEEsmd[ emccluster ] = cl;
1126  cl.stemc( emccluster );
1127 
1128  }
1129 
1130 
1131  }
1132 
1133  }
1134 
1135  else
1136  {
1137  detector -> setCluster( NULL );
1138  }
1139 
1140 
1141 
1142  }
1143 
1144 
1145 }
1146 
1147 
1148 
1149 // ----------------------------------------------------------------------------
1151  {
1152 
1154  StEvent *stevent=(StEvent*)GetInputDS("StEvent");
1155  if ( !stevent ) {
1156  Warning("verifyStEvent","No StEvent found.");
1157  return true;
1158  }
1159  //StEmcCollection *emccollection=stevent->emcCollection();
1160 
1161  Bool_t go = true;
1162 
1163  StEmcDetector *detector=stevent->emcCollection()->detector(kEndcapEmcTowerId);
1164  if ( !detector )
1165  {
1166  Warning("verifyStEvent","detector == NULL for towers");
1167  return false;
1168  }
1169 
1170 
1171  StEmcClusterCollection *cc=detector->cluster();
1172  if ( cc )
1173  {
1174 
1175 
1180 
1181  Float_t emc_sum_towers = 0.;
1182  Float_t eemc_sum_towers = 0.;
1183 
1184  StSPtrVecEmcCluster &emcClusters=cc->clusters();
1185  for ( UInt_t i=0; i<emcClusters.size(); i++ )
1186  {
1187  StEmcCluster *cl=emcClusters[i];
1188  assert(cl);
1189  emc_sum_towers += cl->energy();
1190  }
1191 
1192  for ( Int_t sec=0; sec<12; sec++ )
1193  {
1194 
1195  for ( Int_t i=0; i<numberOfClusters(sec,0); i++ )
1196  eemc_sum_towers += cluster(sec,0,i).energy();
1197  }
1198 
1199  std::cout << "StEEmcIUClusterMaker tower checksum: ";
1200  if ( emc_sum_towers == eemc_sum_towers ) {
1201  std::cout << "passed";
1202  }
1203  else {
1204  std::cout << "FAILED"; go=false;
1205  }
1206  std::cout << std::endl;
1207 
1208  }
1209  else {
1210 
1211  std::cout << "StEEmcIUClusterMaker tower checksum: NULL collection, nclust=" << mNumberOfClusters[0] << std::endl;
1212  go &= (mNumberOfClusters[0]==0);
1213 
1214  }
1215 
1220  Float_t emc_sum_smdu = 0.;
1221  Float_t eemc_sum_smdu = 0.;
1222 
1223  detector=stevent->emcCollection()->detector(kEndcapSmdUStripId);
1224  if ( !detector )
1225  {
1226  Warning("verifyStEvent","detector == NULL for smdu");
1227  return false;
1228  }
1229 
1230  cc=detector->cluster();
1231  if ( cc )
1232  {
1233 
1234  StSPtrVecEmcCluster &smduClusters=cc->clusters();
1235 
1236  for ( UInt_t i=0; i<smduClusters.size(); i++ )
1237  {
1238  StEmcCluster *cl=smduClusters[i];
1239  assert(cl);
1240  emc_sum_smdu += cl->energy();
1241  }
1242 
1243 
1244  for ( Int_t sec=0; sec<12; sec++ )
1245  {
1246 
1247  for ( Int_t i=0; i<numberOfSmdClusters(sec,0); i++ )
1248 
1249  {
1250  eemc_sum_smdu += smdcluster(sec,0,i).energy();
1251 
1252  }
1253 
1254  }
1255 
1256  std::cout << "StEEmcIUClusterMaker smdu checksum: ";
1257  if ( emc_sum_smdu == eemc_sum_smdu ) {
1258  std::cout << "passed";
1259  }
1260  else {
1261  std::cout << "FAILED"; go=false;
1262  }
1263  std::cout << std::endl;
1264 
1265  }
1266  else
1267  {
1268  std::cout << "StEEmcIUClusterMaker smdu checksum: NULL collection, nclust=" << mNumberOfClusters[4] << std::endl;
1269  go &= (mNumberOfClusters[4]==0);
1270  }
1271 
1272 
1274 
1275  Float_t emc_sum_smdv = 0.;
1276  Float_t eemc_sum_smdv = 0.;
1277 
1278  detector=stevent->emcCollection()->detector(kEndcapSmdVStripId);
1279  if (!detector)
1280  {
1281  Warning("verifyStEvent","detector == NULL for smdv");
1282  return false;
1283  }
1284 
1285  cc=detector->cluster();
1286 
1287  if ( cc )
1288  {
1289 
1290  StSPtrVecEmcCluster &smdvClusters=cc->clusters();
1291 
1292  for ( UInt_t i=0; i<smdvClusters.size(); i++ )
1293  {
1294  StEmcCluster *cl=smdvClusters[i];
1295  assert(cl);
1296  emc_sum_smdv += cl->energy();
1297  }
1298 
1299 
1300  for ( Int_t sec=0; sec<12; sec++ )
1301  {
1302 
1303  for ( Int_t i=0; i<numberOfSmdClusters(sec,1); i++ )
1304 
1305  {
1306  eemc_sum_smdv += smdcluster(sec,1,i).energy();
1307 
1308  }
1309 
1310  }
1311 
1312  std::cout << "StEEmcIUClusterMaker smdv checksum: ";
1313  if ( emc_sum_smdv == eemc_sum_smdv ) {
1314  std::cout << "passed";
1315  }
1316  else {
1317  std::cout << "FAILED"; go=false;
1318  }
1319  std::cout << std::endl;
1320 
1321  }
1322  else
1323  {
1324  std::cout << "StEEmcIUClusterMaker smdv checksum: NULL collection, nclust=" << mNumberOfClusters[5] << std::endl;
1325  go &= (mNumberOfClusters[5]==0);
1326  }
1327 
1328 
1329 
1330  return go;
1331 }
1332 
1333 
1334 
1335 // ----------------------------------------------------------------------------
1337 {
1338 
1339  std::cout << "StEEmcIUClusterMaker::print()" << std::endl;
1340  const Char_t *names[] = { "tower", "pre1", "pre2", "post", "smdu", "smdv" };
1341  for ( Int_t i=0;i<6;i++ )
1342  {
1343 
1344  std::cout << "Number of " << names[i]
1345  << " clusters = " << mNumberOfClusters[i]
1346  << std::endl;
1347 
1348  }
1349 
1350  std::cout << "printout of tower clusters follows:" << std::endl;
1351  for ( Int_t sec=0;sec<12;sec++){
1352  for ( Int_t i=0;i<numberOfSmdClusters(sec,0);i++ )
1353  {
1354  StEEmcIUSmdCluster clust=(mSmdClusters[sec][0])[i];
1355  clust.print();
1356  //std::cout << "cluster.key()=" << clust.key() << std::endl;
1357  // std::cout << "cluster.eta()=" << clust.momentum().Eta() << std::endl;
1358  // std::cout << "cluster.phi()=" << clust.momentum().Phi() << std::endl;
1359  // std::cout << "cluster.energy()=" << clust.energy() << std::endl;
1360  }
1361  for ( Int_t i=0;i<numberOfSmdClusters(sec,1);i++ )
1362  {
1363  StEEmcIUSmdCluster clust=(mSmdClusters[sec][1])[i];
1364  clust.print();
1365  }
1366  }
1367 
1368 }
void print()
Prints cluster data.
TVector3 getTowerCenter(const UInt_t sec, const UInt_t sub, const UInt_t etabin) const
Int_t mClusterId
Keep track of clusters.
Int_t Init()
Initialize.
Int_t numberOfNeighbors() const
get the number of neighboring towers
Definition: StEEmcTower.h:54
EEmc ADC –&gt; energy maker.
Int_t key()
Return a unique key assigned by the cluster maker.
std::map< StEmcCluster *, StEEmcIUCluster > mEtoEE
Map StEEmcIUClusters to StEmcClusters.
StEmcCluster * stemc()
void print()
Event summary.
TVector3 momentum()
Get the momentum of this cluster.
Int_t mSuppress
Supress seeds adjacent to clusters.
Int_t numberOfStrips()
Returns the number of SMD strips in the cluster.
Int_t mMinStrips
Minimum number of smd strips to form seed.
void fillStEvent()
Fills StEvent cluster collections if the option is selected.
void neighbor(StEEmcTower *n)
add a tower to list of neighbors
Definition: StEEmcTower.h:52
StEEmcStrip & strip(Int_t sector, Int_t plane, Int_t strip)
std::map< StEmcCluster *, StEEmcIUSmdCluster > mEtoEEsmd
... and for smd clusters
void add(StEEmcTower, Float_t weight=1.0)
void name(const Char_t *n)
Set the name for this element.
Definition: StEEmcElement.h:27
Bool_t buildSmdClusters()
Constructs smd clusters.
Int_t numberOfSmdClusters(Int_t sec, Int_t plane)
Return number of smd clusters for a given sector, plane.
StEEmcIUSmdCluster smdcluster(Int_t sec, Int_t plane, Int_t index)
return a specific cluster from a given sector, plane
std::vector< std::vector< StEEmcTowerVec_t > > mSeedTowers
A base class for describing clusters of EEMC towers.
A cluster maker for the EEMC.
Int_t etabin() const
Returns the etabin of this tower, pre- or postshower element.
Definition: StEEmcTower.h:45
Bool_t buildTowerClusters()
Constructs tower clusters.
Int_t subsector() const
Returns subsector of this tower, pre- or postshower element.
Definition: StEEmcTower.h:43
StEEmcStrip strip(Int_t s)
Returns the specified smd strip w/in the cluster.
StEEmcStripVec_t & strips(Int_t sec, Int_t pln)
Returns a vector of hit strips, given the sector and plane.
void print() const
Print a summary of this tower.
Definition: StEEmcTower.cxx:58
Bool_t mLoose
Loose cuts option.
void index(Int_t i)
Sets the index for this SMD strip, 0..287.
Definition: StEEmcStrip.cxx:32
void index(Int_t i)
Definition: StEEmcTower.cxx:76
Int_t key()
Returns unique id of the cluster.
Base class for representing tower, preshower and postshower elements.
Definition: StEEmcTower.h:11
StEEmcA2EMaker * mEEanalysis
ADC–&gt;E maker.
TString mAnalysisName
ADC–&gt;E maker name.
Int_t mMaxExtent
Maximum distance from SMD seed strips.
StEEmcIUClusterMaker(const Char_t *name="mEEclusters")
Int_t numberOfClusters(Int_t sec, Int_t layer)
Return number of clusters for a given sector, layer.
Int_t size()
Return the size (number of strips) of the cluster.
Int_t sector() const
Returns sector of this tower, pre- or postshower element.
Definition: StEEmcTower.h:41
Float_t mean()
Return the mean strip number of this cluster.
std::vector< std::vector< StEEmcIUClusterVec_t > > mTowerClusters
std::vector< std::vector< StEEmcIUSmdClusterVec_t > > mSmdClusters
StEEmcTower tower(Int_t t)
Get the specified tower.
Int_t mNumberOfClusters[6]
Counts clusters for full eemc, 0=T, 1=P, 2=Q, 3=R, 4=U, 5=V.
Definition: Stypes.h:42
Int_t numberOfTowers()
Get the number of towers in cluster.
EEMC simple geometry.
Definition: Stypes.h:40
Int_t Make()
Make clusters for this event.
A base class for representing clusters of EEMC smd strips.
Float_t energy()
Get energy of this cluster.
void Clear(Option_t *opts="")
Clear clusters for next event.
void seedEnergy(Float_t energy, Int_t layer=0)
void stemc(StEmcRawHit *h)
Sets pointer to the StEmcRawHit when processing an StEvent file.
Definition: StEEmcElement.h:43
Bool_t mFillStEvent
Option to fill StEvent.
Float_t mSeedEnergy[6]
Seed energy for 0=T, 1=P, 2=Q, 3=R, 4=U, 5=V.
Bool_t mSkip
Skip strips if failbit set.
StEEmcIUCluster cluster(Int_t sec, Int_t layer, Int_t index)
Return a specific cluster from a given sector, layer.
void energy(Float_t e)
Set the energy (adc-ped+0.5)/gain for this element.
Definition: StEEmcElement.h:21
Base class for describing an endcap SMD strip.
Definition: StEEmcStrip.h:8
Float_t energy()
Return the energy of this cluster.
StEmcCluster * stemc()
void add(StEEmcStrip strip, Float_t weight=1.0)