StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StEEmcA2EMaker.cxx
1 #include "StEEmcA2EMaker.h"
2 
3 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
4 #include "StMuDSTMaker/COMMON/StMuEmcCollection.h"
5 #include "StMuDSTMaker/COMMON/StMuDst.h"
6 
7 #include "StEventMaker/StEventMaker.h"
8 
9 #include "StEEmcUtil/database/StEEmcDb.h"
10 #include "StEEmcUtil/database/EEmcDbItem.h"
11 
12 
14 #include "StEvent/StEvent.h"
15 #include "StEvent/StEmcDetector.h"
16 #include "StEvent/StEmcModule.h"
17 #include "StEvent/StEmcRawHit.h"
18 #include "StEvent/StEmcCollection.h"
19 
20 #include <iostream>
21 
22 ClassImp(StEEmcA2EMaker);
23 
24 // ----------------------------------------------------------------------------
25 StEEmcA2EMaker::StEEmcA2EMaker(const Char_t *name) : StMaker(name)
26 {
27 
28  scale(1.0);
29 
30  mEEgeom=new EEmcGeomSimple();
31 
32  //$$$ std::cout << "StEEmcA2EMaker(" << name << ")" << std::endl;
34  for ( Int_t tower=0; tower < 720; tower++ ) {
35  for ( Int_t layer=0; layer < 4; layer++ ) {
36  mTowers[tower][layer].Clear("");
37  mTowers[tower][layer].layer(layer);
38  mTowers[tower][layer].index(tower);
39  }
40  }
41 
42 
43 
44  for ( Int_t i=0;i<12;i++ )
45  for ( Int_t j=0;j<6;j++ )
46  {
47  mEnergy[i][j] = 0.;
48  mHits[i][j] = 0;
49  }
50 
51 
52 
54  for ( Int_t tower=0; tower < 720; tower++ )
55  {
56  for ( Int_t layer=0; layer < 4; layer++ )
57  {
58 
59  Int_t phibin=mTowers[tower][layer].phibin();
60  Int_t etabin=mTowers[tower][layer].etabin();
61 
62  //$$$std::cout << "tower " << mTowers[tower][layer].name() << ":";
63 
64  for ( Int_t phi=phibin-1;phi<=phibin+1;phi++ )
65  for ( Int_t eta=etabin-1;eta<=etabin+1;eta++ )
66  {
67 
69  if ( phi==phibin && eta==etabin ) continue;
70 
72  if ( eta<0 || eta>11) continue;
73 
75  Int_t kphi=phi;
76  if ( kphi<0 ) kphi+=60;
77  if ( kphi>59 ) kphi-=60;
78 
79  Int_t neighbor=12*kphi+eta;
80 
81  assert(neighbor>=0 && neighbor<720 );
82 
84  mTowers[tower][layer].neighbor( &mTowers[neighbor][layer] );
85 
86  }
87 
88 
89  }
90  }
91 
92 
93 
95  mHighTower[0]=&mTowers[0][0];
96  mHighTower[1]=&mTowers[0][1];
97  mHighTower[2]=&mTowers[0][2];
98  mHighTower[3]=&mTowers[0][3];
99 
101  for ( Int_t sec=0; sec<12; sec++ ) {
102  for ( Int_t plane=0; plane<2; plane++ ) {
103  for ( Int_t strip=0; strip<288; strip++ ) {
104  mStrips[sec][plane][strip].Clear("");
105  mStrips[sec][plane][strip].sector(sec);
106  mStrips[sec][plane][strip].plane(plane);
107  mStrips[sec][plane][strip].index(strip);
108  }
109  }
110  }
111 
113  mDbMaker = 0;
114  mMuDstMaker = 0;
115  mEventMaker = 0;
116 
117 
118 
120  threshold(3.0,0);
121  threshold(3.0,1);
122  threshold(3.0,2);
123  threshold(3.0,3);
124  threshold(3.0,4);
125  threshold(3.0,5);
126 
129 
130 
131 
133  StEEmcTowerVec_t t;
134  for ( Int_t i = 0; i < 4; i++ ) mHitTowers.push_back(t);
135 
137  StEEmcStripVec_t s;
138  std::vector< StEEmcStripVec_t > plane;
139  plane.push_back(s); // U plane
140  plane.push_back(s); // V plane
141  for ( Int_t i=0; i < 12; i++ ) mHitStrips.push_back(plane);
142 
143 
144 
145 }
146 
147 // ----------------------------------------------------------------------------
148 Int_t StEEmcA2EMaker::Init()
149 {
150  mDbMaker = (StEEmcDb*)this->GetDataSet("StEEmcDb");
151  assert(mDbMaker);
152 
154  if ( mInputType == 1 ) {
155  mMuDstMaker = (StMuDstMaker *)GetMaker(mInputName);
156  }
157  if ( mInputType == 2 ) {
158  mEventMaker = (StEventMaker *)GetMaker(mInputName);
159  }
160  if ( mInputType == 3 ) {
161  mMuDstMaker = (StMuDstMaker *)GetMaker(mInputName);
162  }
163 
164  return StMaker::Init();
165 }
166 
167 // ----------------------------------------------------------------------------
168 Int_t StEEmcA2EMaker::Make()
169 {
170 
172  if ( !readData() ) return kStWarn;
173 
174  return kStOK;
175 }
176 
177 // ----------------------------------------------------------------------------
179 {
180 
182  if ( mInputType == 1 ) {
183  if ( mMuDstMaker ) return readMuDst();
184  }
185  if ( mInputType == 2 ) {
186  // if ( mEventMaker )
187  return readStEvent();
188  }
189  if ( mInputType == 3 ) {
190  if ( mMuDstMaker ) return readEzt();
191  }
192 
193  return false;
194 }
195 
196 // ----------------------------------------------------------------------------
198 {
199  assert(mMuDstMaker);
200 
202  StMuEmcCollection *emc = mMuDstMaker->muDst()->muEmcCollection();
203  if ( !emc ) return false;
204 
205  fillFromMuDst(emc);
206 
207  return true;
208 
209 }
210 
212 {
213 
214  assert(emc);
215 
217  for ( Int_t ihit = 0; ihit < emc -> getNEndcapTowerADC(); ihit++ ) {
218 
219  Int_t adc, sec, sub, eta;
220  emc -> getEndcapTowerADC ( ihit, adc, sec, sub, eta );
221 
222  sec--; // Counting from 1 is insane in c++! It ends here.
223  sub--; // Eveerything in my classes assume indexes from 0.
224  eta--; //
225 
226  assert( sec >= 0 && sec < 12 ); // Indexing errors detected
227  assert( sub >= 0 && sub < 5 ); // Indexing errors detected
228  assert( eta >= 0 && eta < 12 ); // Indexing errors detected
229 
230  addTowerHit(sec,sub,eta,adc,0);
231 
232  }
233 
235  for ( Int_t ihit = 0; ihit < emc -> getNEndcapPrsHits(); ihit++ ) {
236 
237  Int_t adc, sec, sub, eta, det;
238  StMuEmcHit *hit = emc -> getEndcapPrsHit(ihit, sec, sub, eta, det);
239  adc = hit -> getAdc();
240 
241  sec--; // Counting from 1 is insane in c++! It ends here.
242  sub--; // Eveerything in my classes assume indexes from 0.
243  eta--; //
244 
245  assert( sec >= 0 && sec < 12 ); // Indexing errors detected
246  assert( sub >= 0 && sub < 5 ); // Indexing errors detected
247  assert( eta >= 0 && eta < 12 ); // Indexing errors detected
248  assert( det >= 1 && det < 4 ); // Indexing errors detected
249 
250  addTowerHit(sec,sub,eta,(Float_t)adc,det);
251 
252  }
253 
255  Char_t cpl[] = { 'U','V' };
256  for ( Int_t plane = 0; plane < 2; plane++ )
257  for ( Int_t ihit = 0; ihit < emc -> getNEndcapSmdHits(cpl[plane]); ihit++ ) {
258 
259  Int_t sec, strip, adc;
260  StMuEmcHit *hit = emc -> getEndcapSmdHit( cpl[plane], ihit, sec, strip );
261  adc = hit -> getAdc();
262 
263  sec--; // Counting from 1 is insane in c++! It ends here.
264  strip--; // Eveerything in my classes assume indexes from 0.
265 
266  assert( sec >= 0 && sec < 12 ); // Indexing errors detected
267  assert( strip >= 0 && strip < 288 ); // Indexing errors detected
268 
269  addSmdHit(sec,plane,strip,(Float_t)adc);
270 
271  }
272 
273  return true;
274 }
275 
276 // ---
277 Bool_t StEEmcA2EMaker::readEzt()
278 {
279  assert(0); // Please implement me
280  return true;
281 }
282 
283 // ----------------------------------------------------------------------------
284 void StEEmcA2EMaker::addTowerHit( Int_t sec, Int_t sub, Int_t eta, Float_t adc, Int_t layer )
285 {
286  assert(sec>=0 && sec<12);
287  assert(sub>=0 && sub<5);
288  assert(eta>=0 && eta<12);
289  assert(layer>=0 && layer < 4);
290 
291  Int_t index=60*sec+12*sub+eta;
292  assert(0<=index && index<720); // Just can't stay in bounds today
293 
294  static Char_t subsectors[] = { 'A','B','C','D','E' };
295  static Char_t detectors[] = { 'T', 'P', 'Q', 'R' };
296 
300 
301  const EEmcDbItem *dbitem = mDbMaker -> getTile( sec+1,subsectors[sub],eta+1, detectors[layer] );
302  assert(dbitem);
303 
304 
305 
307  mTowers[index][layer].fail( dbitem -> fail );
308  mTowers[index][layer].stat( dbitem -> stat );
309 
311 #if 1
312  if ( dbitem -> fail ) return;
313  //if ( dbitem -> stat ) return;
314 #endif
315 
316  Float_t ped = dbitem -> ped;
317  Float_t gain = dbitem -> gain;
318  Float_t threshold = ped + dbitem -> sigPed * mSigmaPed[layer];
319 
321  if ( adc < threshold ) return;
322  mHits[sec][layer]++;
323 
325  mTowers[index][layer].raw(adc);
326  mTowers[index][layer].adc(adc-ped);
327 
329  if ( gain <= 0. ) return;
331  Float_t energy = ( adc - ped + 0.5 ) / gain;
332 
334  mTowers[index][layer].energy( energy * mScale );
335 
337  UInt_t s=(UInt_t)mTowers[index][layer].sector();
338  UInt_t ss=(UInt_t)mTowers[index][layer].subsector();
339  UInt_t eb=(UInt_t)mTowers[index][layer].etabin();
340  TVector3 momentum=mEEgeom -> getTowerCenter( s,ss,eb );
341  momentum=momentum.Unit();
342  momentum*=energy;
343  mTowers[index][layer].et( (Float_t)momentum.Perp() );
344 
345  if ( !mTowers[index][layer].fail() )
346  if ( adc - ped > mHighTower[layer]->adc() ) {
347  mHighTower[layer] = &mTowers[index][layer];
348  }
349 
351 
352  mHitTowers[layer].push_back( mTowers[index][layer] );
353  mEnergy[sec][layer] += energy;
354 
355 #if 0
356  mTowers[index][layer].print();
357 #endif
358 
359 }
360 //---
361 void StEEmcA2EMaker::addSmdHit( Int_t sec, Int_t plane, Int_t strip, Float_t adc )
362 {
363  assert(0<=sec && sec<12);
364  assert(0<=plane && plane<2);
365  assert(0<=strip && strip<288);
366 
368  static Char_t planes[] = { 'U', 'V' };
369  const EEmcDbItem *dbitem = mDbMaker -> getByStrip ( sec+1, planes[plane], strip+1 );
370 
372  if ( dbitem == 0 ) return;
373 
375  mStrips[sec][plane][strip].fail( dbitem -> fail );
376  mStrips[sec][plane][strip].stat( dbitem -> stat );
377 
379  if ( dbitem -> fail ) return;
380  //if ( dbitem -> stat ) return;
381 
382  Float_t ped = dbitem -> ped;
383  Float_t gain = dbitem -> gain;
384  Float_t threshold = ped + dbitem -> sigPed * mSigmaPed[4+plane];
385 
387  if ( adc < threshold ) return;
388  mHits[sec][plane+4]++;
389 
391  mStrips[sec][plane][strip].raw(adc);
392  mStrips[sec][plane][strip].adc(adc-ped);
393 
395  if ( gain <= 0. ) return;
396 
398  Float_t energy = ( adc - ped + 0.5 ) / gain;
399 
400  mStrips[sec][plane][strip].energy(energy);
401 
402  mHitStrips[sec][plane].push_back( mStrips[sec][plane][strip] );
403 
404  mEnergy[sec][plane+4] += energy;
405 
406 }
407 
408 
409 // ----------------------------------------------------------------------------
410 void StEEmcA2EMaker::Clear(Option_t *opts)
411 {
412 
413 
415  for ( Int_t layer = 0; layer < 4; layer++ )
416  mHitTowers[layer].clear();
417 
418  for ( Int_t sector=0; sector<12; sector++ )
419  for ( Int_t plane = 0; plane < 2; plane++ )
420  mHitStrips[sector][plane].clear();
421 
422 
426 
427  for ( Int_t i = 0; i < 720; i++ )
428  for ( Int_t j = 0; j < 4; j++ )
429  mTowers[i][j].Clear("");
430 
431  for ( Int_t i = 0; i < 12; i++ )
432  for ( Int_t j = 0; j < 2; j++ )
433  for ( Int_t k = 0; k < 288; k++ ){
434  mStrips[i][j][k].Clear("");
435  }
436 
438  for ( Int_t i = 0; i < 12; i++ )
439  for ( Int_t j = 0; j < 6; j++ ) {
440  mEnergy[i][j] = 0.;
441  mHits[i][j]=0;
442  }
443 
444 
445 
446 }
447 
448 
449 // ---
450 Bool_t StEEmcA2EMaker::readStEvent()
451 {
452 
453  StEvent *event=(StEvent*)GetInputDS("StEvent");
454  assert(event);
455  StEmcCollection *emc=event->emcCollection();
456  assert(emc);
457 
458  fillFromSt( emc );
459 
460  return true;
461 }
462 
464 {
465 
466 
470  StEmcDetector *detector=emc->detector(kEndcapEmcTowerId);
471  if ( !detector )
472  {
473  Warning("fillFromSt","\n**********\nemc->detector() NULL for eemc towers. MAJOR PROBLEM, trying to procede.\n**********\n\n");
474  return false;
475  }
476 
477  for ( UInt_t sec = 0; sec < detector -> numberOfModules(); sec++ )
478  {
479 
482  StEmcModule *sector = detector -> module( sec+1 );
483  assert(sector);
484  StSPtrVecEmcRawHit &hits = sector->hits();
485 
487  for ( UInt_t ihit=0; ihit<hits.size(); ihit++ )
488  {
489 
490  assert(hits[ihit]);
491 
494  Int_t isector = hits[ihit]->module() - 1;
495  Int_t isubsector = hits[ihit]->sub() - 1;
496  Int_t ietabin = hits[ihit]->eta() - 1;
497  Int_t adc = hits[ihit]->adc();
498  assert(isector>=0 && isector<12);
499  assert(isubsector>=0 && isubsector<5);
500  assert(ietabin>=0 && ietabin<12);
501 
502  Int_t iphibin = 5 * isector + isubsector;
503  assert(iphibin>=0 && iphibin<60);
504  Int_t index = 12 * iphibin + ietabin;
505  assert(index>=0 && index<720);
506 
507 #if 0
508  std::cout << "isector=" << isector
509  << " isubsector=" << isubsector
510  << " ietabin=" << ietabin
511  << " " << hits[ihit]->module()
512  << " " << hits[ihit]->sub()
513  << " " << hits[ihit]->eta()
514  << " " << hits[ihit]->adc()
515  << std::endl;
516 #endif
517 
519  mTowers[index][0].stemc( hits[ihit] );
520 
521 
527  addTowerHit(isector,isubsector,ietabin,(Float_t)adc,0);
528 
529  }
530 
531  }
532 
533 
534 #if 0
535  for ( UInt_t ii=0; ii < mHitTowers.size(); ii++ )
537  {
538  std::cout << "mHitTowers[0][" << ii << "] "
539  << mHitTowers[0][ii].sector() << " "
540  << mHitTowers[0][ii].subsector() << " "
541  << mHitTowers[0][ii].etabin() << " "
542  << mHitTowers[0][ii].adc() << " "
543  <<std::endl << " ====> "
544  << mHitTowers[0][ii].stemc()
545  << std::endl;
546  }
548 #endif
549 
550 
551 
555 
556  detector=emc->detector(kEndcapEmcPreShowerId);
557 
558  if ( !detector )
559  {
560  Warning("fillFromSt","\n**********\nemc->detector() NULL for eemc pre/post. MAJOR PROBLEM, trying to procede.\n**********\n\n");
561  }
562  else
563  for ( UInt_t sec = 0; sec < detector -> numberOfModules(); sec++ )
564  {
565 
568  StEmcModule *sector = detector -> module( sec+1 );
569  StSPtrVecEmcRawHit &hits = sector->hits();
570 
572  for ( UInt_t ihit=0; ihit<hits.size(); ihit++ )
573  {
574 
575  assert(hits[ihit]);
576 
579  Int_t isector = hits[ihit]->module() - 1;
580  Int_t isubsector = (hits[ihit]->sub() - 1) % 5;
581  Int_t ietabin = hits[ihit]->eta() - 1;
582  Int_t adc = hits[ihit]->adc();
583 
584  Int_t ilayer = (hits[ihit]->sub() - 1) / 5 + 1;
585  assert(ilayer >= 1 && ilayer < 4 );
586 
587  Int_t iphibin = 5*isector+isubsector;
588  Int_t index = 12*iphibin + ietabin;
589 
590 
591  mTowers[index][ilayer].stemc( hits[ihit] );
592 
595  addTowerHit( isector, isubsector, ietabin, (Float_t)adc, ilayer );
596 
597  }
598 
599  }
600 
601 
605  StDetectorId ids[] = { kEndcapSmdUStripId, kEndcapSmdVStripId };
607  for ( Int_t iplane=0; iplane<2; iplane++ )
608  {
609 
611  detector=emc->detector( ids[iplane] );
612  if ( !detector )
613  {
614  Warning("fillFromSt","\n**********\nemc->detector() NULL for esmd plane. MAJOR PROBLEM, trying to procede.\n**********\n\n");
615  }
616  else
617  for ( UInt_t sec = 0; sec < detector -> numberOfModules(); sec++ )
618  {
619 
622  StEmcModule *sector = detector -> module( sec+1 );
623  StSPtrVecEmcRawHit &hits = sector->hits();
624 
626  for ( UInt_t ihit=0; ihit<hits.size(); ihit++ )
627  {
628  Int_t isector=hits[ihit]->module()-1;
629  Int_t istrip=hits[ihit]->eta()-1;
630  Int_t adc =hits[ihit]->adc();
631  assert(isector>=0 && isector<12);
632  assert(istrip>=0 && istrip < 288);
633 
634  mStrips[isector][iplane][istrip].stemc( hits[ihit] );
635 
637  addSmdHit( isector, iplane, istrip, (Float_t)adc);
638 
639  }
640 
641  }
642 
643  }
644 
645 
646 
647  return true;
648 }
649 
650 
651 
652 Float_t StEEmcA2EMaker::energy(Int_t layer)
653 {
654  Float_t sum=0.;
655  for ( Int_t sector=0;sector<12;sector++ ) sum+= energy(sector,layer);
656  return sum;
657 }
658 
659 
660 // ----------------------------------------------------------------------------
661 StEEmcTower *StEEmcA2EMaker::tower( TVector3 &r, Int_t layer )
662 {
663  Int_t sec=-1,sub=-1,eta=-1;
664  if ( !mEEgeom->getTower(r,sec,sub,eta) ) return NULL;
665  return &mTowers[ index(sec,sub,eta) ][layer];
666 }
StEEmcA2EMaker(const Char_t *name="mEEanalysis")
Float_t mEnergy[kEEmcNumSectors][6]
EEmc ADC –&gt; energy maker.
void stat(unsigned s)
Set a status bit for this element.
Definition: StEEmcElement.h:23
Bool_t fillFromMuDst(const StMuEmcCollection *emc)
virtual void Clear(Option_t *opts="")
Clears the element.
StMuDst * muDst()
Definition: StMuDstMaker.h:425
Float_t addTowerHit(Int_t sec, Int_t sub, Int_t eta, Float_t adc, Int_t layer)
void layer(Int_t l)
Sets the layer, [0-4]=[TPQR].
Definition: StEEmcTower.h:31
void scale(Float_t s)
StEEmcTower * mHighTower[4]
std::vector< StEEmcTowerVec_t > mHitTowers
Float_t addSmdHit(Int_t sec, Int_t plane, Int_t str, Float_t adc)
Int_t mHits[kEEmcNumSectors][6]
Number of hits in layer.
StEEmcStrip mStrips[kEEmcNumSectors][kEEmcNumSmdUVs][kEEmcNumStrips]
Array of 12x2x288 smd strip objects.
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)
virtual void Clear(Option_t *opts="")
Clears the tower.
Definition: StEEmcTower.h:72
const EEmcGeomSimple * mEEgeom
Int_t index(Int_t sector, Int_t subsector, Int_t etabin) const
Given tower sector, subsector, etabin, translate into index.
void raw(Float_t r)
Set the raw ADC for this element.
Definition: StEEmcElement.h:17
Bool_t fillFromSt(const StEmcCollection *emc)
If StEvent is used, we will fill additional parts of StEEmcElement.
Int_t etabin() const
Returns the etabin of this tower, pre- or postshower element.
Definition: StEEmcTower.h:45
void fail(unsigned f)
Set a fail bit for this element.
Definition: StEEmcElement.h:25
Int_t phibin(Int_t sector, Int_t subsector) const
Given tower sector, subsector, translate to phibin.
void threshold(Float_t cut, Int_t layer)
const StEEmcDb * mDbMaker
std::vector< std::vector< StEEmcStripVec_t > > mHitStrips
Float_t energy(Int_t sector, Int_t layer) const
void print() const
Print a summary of this tower.
Definition: StEEmcTower.cxx:58
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
Base class for representing tower, preshower and postshower elements.
Definition: StEEmcTower.h:11
virtual Int_t Init()
Initialize.
static StMuEmcCollection * muEmcCollection()
returns pointer to current StMuEmcCollection
Definition: StMuDst.h:389
Definition: Stypes.h:42
EEMC simple geometry.
Definition: Stypes.h:40
Int_t phibin() const
Returns the phibin of this tower.
Definition: StEEmcTower.h:47
void et(Float_t e)
Definition: StEEmcTower.h:34
void adc(Float_t a)
Set the pedestal-subtracted ADC for this element.
Definition: StEEmcElement.h:19
virtual void Clear(Option_t *opts="")
Clear the maker for next event.
StEEmcTower & tower(Int_t index, Int_t layer=0)
Float_t mSigmaPed[6]
void stemc(StEmcRawHit *h)
Sets pointer to the StEmcRawHit when processing an StEvent file.
Definition: StEEmcElement.h:43
StEEmcTower mTowers[kEEmcNumSectors *kEEmcNumSubSectors *kEEmcNumEtas][4]
Array of 720 x 4 tower objects.
void plane(Int_t p)
Sets the plane for this SMD strip, 0=U, 1=V.
Definition: StEEmcStrip.h:19
bool getTower(const TVector3 &r, int &sec, int &sub, int &etabin, Float_t &dphi, Float_t &deta) const
void sector(Int_t s)
Sets the sector for this SMD strip.
Definition: StEEmcStrip.h:17
void energy(Float_t e)
Set the energy (adc-ped+0.5)/gain for this element.
Definition: StEEmcElement.h:21
virtual Int_t Make()
Read and process one event.