StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
EEmcMCData.cxx
1 // $Id: EEmcMCData.cxx,v 1.18 2010/08/26 22:48:54 ogrebeny Exp $
2 
3 #include "StEventTypes.h"
4 
5 //#include "emc_def.h" /* FIXME - move it to pams */
6 
7 
8 //#include "TError.h"
9 //
10 #include "tables/St_g2t_emc_hit_Table.h"
11 #include "tables/St_g2t_event_Table.h"
12 //
13 #include "StBFChain.h"
14 //
15 #include "EEmcMCData.h"
16 
17 #include "StEEmcUtil/EEmcGeom/EEmcGeomDefs.h"
18 #include "StEEmcUtil/EEevent/EEeventDst.h"
19 #include "StEEmcUtil/EEevent/EEsectorDst.h"
20 
21 #include "StMessMgr.h"
22 
23 ClassImp(EEmcMCData)
24 
25 
26 //-------------------------------------------------------------------------
27 //-------------------------------------------------------------------------
28 EEmcMCData::EEmcMCData()
29  : TObject()
30 {
31  mEventID = -1;
32  mLastHit = 0;
33  mEthr = kEEmcDefaultEnergyThreshold;
34  mSize = kEEmcDefaultMCHitSize;
35  mHit = new struct EEmcMCHit[mSize]; // for now a static array
36 }
37 
38 //-------------------------------------------------------------------------
39 //-------------------------------------------------------------------------
40 EEmcMCData::EEmcMCData(const EEmcMCData &e )
41  : TObject(e)
42 {
43  mEventID = e.getEventID();
44  mLastHit = e.getLastHit();
45  mEthr = e.getEnergyThreshold();
46  mSize = e.getSize();
47  mHit = new struct EEmcMCHit[mSize];
48  mSize = e.getHitArray(mHit,mSize);
49 }
50 
51 
52 //-------------------------------------------------------------------------
53 //-------------------------------------------------------------------------
54 EEmcMCData::~EEmcMCData() {
55  mSize = 0;
56  if(mHit) delete [] mHit;
57 }
58 
59 
60 //-------------------------------------------------------------------------
61 // decode Endcap Emc data
62 //-------------------------------------------------------------------------
63 Int_t
64 EEmcMCData::readEventFromChain(const StMaker *myMk)
65 {
66  St_g2t_event *g2t_event = NULL; // g2t event
67  St_g2t_emc_hit *emc_hit = NULL; // tower hits
68  St_g2t_emc_hit *smd_hit = NULL; // smd hits
69  g2t_event_st *event_head = NULL; // g2t event header
70 
71  // get info from event header
72  if( (g2t_event=(St_g2t_event *) myMk->GetDataSet("g2t_event")) == NULL ) {
73  LOG_ERROR << "missing MC event header" << endm;
74  return 0;
75  }
76  if( (event_head= g2t_event->GetTable()) == NULL ) {
77  LOG_ERROR << "missing MC event header table" << endm;
78  return 0;
79  }
80  mEventID=event_head->n_event;
81  emc_hit = (St_g2t_emc_hit *) myMk->GetDataSet("g2t_eem_hit");
82  smd_hit = (St_g2t_emc_hit *) myMk->GetDataSet("g2t_esm_hit");
83  unpackGeantHits(emc_hit, smd_hit);
84  return mLastHit;
85 
86 }
87 
88 
89 //-------------------------------------------------------------------------
90 // deal with GEANT volume ID mess
91 //-------------------------------------------------------------------------
92 void
93 EEmcMCData::unpackGeantHits(St_g2t_emc_hit* emc_hit, St_g2t_emc_hit* smd_hit ){
94  int err=0;
95  mLastHit = 0;
96  memset(mHit,0x0,sizeof(struct EEmcMCHit)*mSize);
97 
98  // get tower data
99  if( emc_hit !=NULL ) {
100  Int_t nhits = emc_hit->GetNRows();
101  if(nhits<=0) {
102  Warning("readEventFromChain","no tower hits (%d)",nhits);
103  goto skipTower;
104  }
105 
106  g2t_emc_hit_st *hit = emc_hit->GetTable();
107  for(Int_t ihit=0; ihit<nhits; ihit++,hit++) {
108 
109  Int_t ivid = hit->volume_id;
110  Short_t sec = 0;
111  Short_t ssec = 0;
112  // Short_t half = ivid/kEEmcTowerHalfId;
113  ivid %= kEEmcTowerHalfId;
114  Short_t phi = ivid/kEEmcTowerPhiId; ivid %= kEEmcTowerPhiId;
115  Short_t eta = ivid/kEEmcTowerEtaId; ivid %= kEEmcTowerEtaId;
116  Short_t depth = ivid/kEEmcTowerDepId; ivid %= kEEmcTowerDepId;
117  //printf("Tw hit->volume_id=%d hit->de=%g phi=%d\n",hit->volume_id,hit->de,phi);
118 
119  if(!ivid==0){err=1; goto crash;};
120 
121  ssec = (phi-1)%5 + 1; // sub-sector
122 
123  sec = (phi-1)/5 + 1;
124 
125  if(!( 0<sec && sec<=kEEmcNumSectors )){err=2; goto crash;};
126  if(!( 0<ssec && ssec<=kEEmcNumSubSectors)){err=3; goto crash;};
127  if(!( 0<eta && eta<=kEEmcNumEtas)){err=4; goto crash;};
128  if(!( 0<depth && depth<=kEEmcNumDepths)){err=5; goto crash;};
129 
130  switch(depth) {
131  case kPreShower1Depth:
132  mHit[mLastHit].detector = kEEmcMCPreShower1Id;
133  break;
134  case kPreShower2Depth:
135  mHit[mLastHit].detector = kEEmcMCPreShower2Id;
136  break;
137  case kTower1Depth:
138  case kTower2Depth:
139  mHit[mLastHit].detector = kEEmcMCTowerId;
140  break;
141  case kPostShowerDepth:
142  mHit[mLastHit].detector = kEEmcMCPostShowerId;
143  break;
144  default:
145  Warning("readEventFromChain","unknown depth %d",depth);
146  goto crash;
147  break;
148  }
149 
150  mHit[mLastHit].sector = sec;
151  mHit[mLastHit].tower.ssec = ssec;
152  mHit[mLastHit].tower.eta = eta;
153  mHit[mLastHit].de = hit->de;
154  mHit[mLastHit].track_p = hit->track_p;
155  mLastHit++;
156  // printf("depth=%d nH=%d\n",depth,mLastHit);
157 
158  if ((mLastHit >= mSize) && !expandMemory()) {
159  LOG_ERROR << "failed expandMemory() for tower tails" << endm;
160  goto crash;
161  }
162  } // end of tower hits
163  }
164 
165  skipTower:
166 
167  // get smd data
168  if( smd_hit != NULL ) {
169  Int_t nhits = smd_hit->GetNRows();
170 
171  if(nhits<=0) {
172  Warning("readEventFromChain","no smd hits (%d)",nhits);
173  goto done;
174  }
175  g2t_emc_hit_st *hit = smd_hit->GetTable();
176 
177  for(Int_t ihit=0; ihit<nhits ; ihit++,hit++) {
178  //printf("Smd hit->volume_id=%d hit->de=%g\n",hit->volume_id,hit->de);
179 
180  Int_t ivid = hit->volume_id;
181  Short_t det = 0;
182  Short_t sec = 0;
183  Short_t half = ivid/kEEmcSmdHalfId; ivid %= kEEmcSmdHalfId;
184  Short_t phi = ivid/kEEmcSmdPhiId; ivid %= kEEmcSmdPhiId;
185  Short_t plane = ivid/kEEmcSmdPlaneId;ivid %= kEEmcSmdPlaneId;
186  Short_t strip = ivid/kEEmcSmdStripId;ivid %= kEEmcSmdStripId;
187 
188  if(!ivid==0){err=10; goto crash;};
189 
190  switch(phi) { /* FIXME ONE DAY */
191  case 1:
192  case 4:
193  case 7:
194  case 10:
195  switch(plane) {
196  case 1: det = kEEmcMCSmdVStripId; break;
197  case 3: det = kEEmcMCSmdUStripId; break;
198  default: det = kUnknownId; break;
199  }
200  break;
201  case 2:
202  case 5:
203  case 8:
204  case 11:
205  switch(plane) {
206  case 2: det = kEEmcMCSmdVStripId; break;
207  case 1: det = kEEmcMCSmdUStripId; break;
208  default: det = kUnknownId; break;
209  }
210  break;
211  case 3:
212  case 6:
213  case 9:
214  case 12:
215  switch(plane) {
216  case 3: det = kEEmcMCSmdVStripId; break;
217  case 2: det = kEEmcMCSmdUStripId; break;
218  default: det = kUnknownId; break;
219  }
220  break;
221  default:
222  det = kUnknownId;
223  break;
224  }
225  if(det!=kEEmcMCSmdVStripId && det!=kEEmcMCSmdUStripId ) {
226  Warning("readEventFromChain","unknown smd layer %d %d-%d-%d-%d",det,half,phi,plane,strip);
227  goto crash;
228  }
229 
230  sec = phi;
231 
232  if(! ( 0<sec && sec <=kEEmcNumSectors )){err=12; goto crash;};
233  if(!( 0<strip && strip<=kEEmcNumStrips )){err=13; goto crash;};
234 
235  // printf("Smd hit->volume_id=%d hit->de=%g sec=%d U/V=%d, strip=%d\n",hit->volume_id,hit->de,sec, det,strip );
236  // fill in
237  mHit[mLastHit].detector = det;
238  mHit[mLastHit].sector = sec;
239  mHit[mLastHit].strip = strip;
240  mHit[mLastHit].de = hit->de;
241  mHit[mLastHit].track_p = hit->track_p;
242 
243  mLastHit++;
244  if ((mLastHit >= mSize) && !expandMemory()) {
245  LOG_ERROR << "failed expandMemory() for SMD strips" << endm;
246  goto crash;
247  }
248  }
249  }
250 
251  done:
252  return ;
253 
254  crash:
255  printf("\n\n==============================\nEEmcMCData::readEventFromChain() Fatal error while decoding .fzd hits for EEMC err=%d,\n all EEMC data erased \n=====================================\n",err);
256  mLastHit = 0;
257  memset(mHit,0x0,sizeof(struct EEmcMCHit)*mSize);
258  return;
259 
260 }
261 
262 
263 //-------------------------------------------------------------------------
264 //-------------------------------------------------------------------------
265 Int_t
266 EEmcMCData::getHitArray(EEmcMCHit *h, Int_t size) const
267 {
268  if(size<=0) {
269  LOG_ERROR << "invalid size: " << size << endm;
270  return 0;
271  }
272  if(size<mSize) {LOG_WARN << "truncating to " << size << " hits (out of " << mSize << ")" << endm;}
273  int n = size; if (n > mSize) n = mSize;
274  memcpy(h,mHit,size*sizeof(EEmcMCHit));
275  return size;
276 }
277 
278 Int_t
279 EEmcMCData::setHitArray(EEmcMCHit *h, Int_t size)
280 {
281  if(size<=0) {
282  LOG_ERROR << "invalid size: " << size << endm;
283  return 0;
284  }
285  if(size>mSize) {LOG_WARN << "truncating to " << mSize << " hits (out of " << size << ")" << endm;}
286  int n = size; if (n > mSize) n = mSize;
287  memcpy(h,mHit,n*sizeof(EEmcMCHit));
288  return mSize;
289 }
290 
291 //-------------------------------------------------------------------------
292 //-------------------------------------------------------------------------
293 void
294 EEmcMCData::print() const
295 {
296  TString detName;
297  Int_t detId;
298  EEmcMCHit *h = mHit;
299  printf("EndcapEmc MC event #%d\n",mEventID);
300  for(Int_t i=0; i<mLastHit; i++,h++) {
301  detId=h->detector;
302  switch(detId) {
303  case kEEmcMCTowerId: detName = "EndcapEmcTower "; break;
304  case kEEmcMCPreShower1Id: detName = "EndcapEmcPreShower1 "; break;
305  case kEEmcMCPreShower2Id: detName = "EndcapEmcPreShower2 "; break;
306  case kEEmcMCPostShowerId: detName = "EndcapEmcPostShower "; break;
307  case kEEmcMCSmdUStripId: detName = "EndcapEmcSmdUStrip "; break;
308  case kEEmcMCSmdVStripId: detName = "EndcapEmcSmdVStrip "; break;
309  default: detName = "EndcapEmcUnknown "; break;
310  }
311 
312 
313  switch(detId) {
314  case kEEmcMCTowerId:
315  case kEEmcMCPreShower1Id:
316  case kEEmcMCPreShower2Id:
317  case kEEmcMCPostShowerId:
318  printf("%s sec=%2d sub=%c eta=%d de=%g tr_p=%d\n",detName.Data(), h->sector,h->tower.ssec-1+'A',h->tower.eta,h->de, h->track_p);
319  break;
320  case kEEmcMCSmdUStripId:
321  case kEEmcMCSmdVStripId:
322  printf("%s sec=%2d strip=%d de=%g tr_p=%d\n",detName.Data(), h->sector,h->strip,h->de, h->track_p);
323  break;
324  default:
325  LOG_WARN << "detectorId=" << detId << " is unknown" << endm;
326  break;
327  }
328 
329  }
330 }
331 
332 //-------------------------------------------------------------------------
333 //
334 //-------------------------------------------------------------------------
335 Int_t
336 EEmcMCData::expandMemory()
337 {
338  Int_t newSize = mSize + kEEmcDefaultMCHitSize;
339  EEmcMCHit* newHit = new EEmcMCHit[newSize];
340 
341  if(mHit) {
342  if (newHit) memcpy(newHit,mHit,mSize*sizeof(EEmcMCHit));
343  delete [] mHit;
344  mHit = 0;
345  }
346  Info("expandMemory"," MCHit memory expanded from 0x%04x to 0x%04x",mSize,newSize);
347  mSize = newSize;
348  mHit = newHit;
349  return kTRUE;
350 }
351 
352 
353 //-------------------------------------------------------------------------
354 //
355 //-------------------------------------------------------------------------
356 Int_t
357 EEmcMCData::read (void *hit, int size)
358 {
359 
360  const int HitSize = sizeof(struct EEmcMCHit);
361  const int MaxSize = mSize * HitSize;
362  if(size>MaxSize) size=MaxSize;
363  memcpy(mHit,hit,size);
364  mLastHit =size/HitSize;
365  return size;
366 }
367 
368 //-------------------------------------------------------------------------
369 // writes binary hits
370 //-------------------------------------------------------------------------
371 Int_t
372 EEmcMCData::write (void *hit, int size)
373 {
374  Int_t nbytes = mLastHit * sizeof(struct EEmcMCHit);
375  if(size>nbytes) size=nbytes;
376  memcpy(hit,mHit,size);
377  return size;
378 }
379 
380 
381 //-------------------------------------------------------
382 //-------------------------------------------------------
383 Int_t EEmcMCData::write(EEeventDst *EEeve) {
384  //printf("Clear & Export EEeventDst to TTree\n");
385  EEeve->clear();
386 
387  EEeve->setType(EEeventDst::kRawMC);
388  EEeve->setID(mEventID);
389 
390  EEmcMCHit *h = mHit;
391  for(Int_t i=0; i<mLastHit; i++,h++) {
392  // Info("EEmcMCData().write(TTree)","depth=%2d sector=%2d sub=%2d eta=%2d energy=%f", h->depth,h->sec,h->ssec,h->eta,h->de);
393 
394  int secID=h->sector;
395  EEsectorDst *EEsec= (EEsectorDst *)EEeve->getSec(secID,1);
396 
397  // temporary projection, JB
398 
399  int subSec='A'+h->tower.ssec-1;
400  if (EEsec) switch( h->detector) {
401  case kEEmcMCTowerId:
402  EEsec->addTwHit(subSec,h->tower.eta,h->de); break;
403  case kEEmcMCPreShower1Id:
404  EEsec->addPre1Hit(subSec,h->tower.eta,h->de); break;
405  case kEEmcMCPreShower2Id:
406  EEsec->addPre2Hit(subSec,h->tower.eta,h->de); break;
407  case kEEmcMCPostShowerId:
408  EEsec->addPostHit(subSec,h->tower.eta,h->de); break;
409  case kEEmcMCSmdUStripId:
410  EEsec->addSmdUHit(h->strip,h->de); break;
411  case kEEmcMCSmdVStripId:
412  EEsec->addSmdVHit(h->strip,h->de); break;
413  default:
414  LOG_ERROR << "invalid MC depth" << endm;
415  break;
416  }
417 
418  } // end of loop over all hits
419 
420  return 0;
421 }
422 
423 // $Log: EEmcMCData.cxx,v $
424 // Revision 1.18 2010/08/26 22:48:54 ogrebeny
425 // Improved constness
426 //
427 // Revision 1.17 2009/04/29 22:22:51 ogrebeny
428 // Bug fixed - missing {}, thanks to Pibero.
429 //
430 // Revision 1.16 2009/02/11 20:37:37 ogrebeny
431 // *** empty log message ***
432 //
433 // Revision 1.15 2009/02/11 20:35:58 ogrebeny
434 // No asserts, no exceptions.
435 //
436 // Revision 1.14 2007/07/12 19:30:14 fisyak
437 // Add includes for ROOT 5.16
438 //
439 // Revision 1.13 2006/08/10 03:20:00 perev
440 // Assert==>assert
441 //
442 // Revision 1.12 2005/06/03 19:19:48 balewski
443 // for embedding, GEANT unpcker was split on 2 parts
444 //
445 // Revision 1.11 2004/04/08 21:34:55 perev
446 // Leak off
447 //
448 // Revision 1.10 2003/11/12 19:59:13 balewski
449 // *** empty log message ***
450 //
451 // Revision 1.9 2003/09/17 22:05:36 zolnie
452 // delete mumbo-jumbo
453 //
454 // Revision 1.8 2003/09/11 19:41:08 zolnie
455 // updates for gcc3.2
456 //
457 // Revision 1.7 2003/09/02 17:57:56 perev
458 // gcc 3.2 updates + WarnOff
459 //
460 // Revision 1.6 2003/03/26 21:16:59 balewski
461 // swap U & V on Wei-Ming request
462 //
463 // Revision 1.5 2003/03/22 19:37:33 balewski
464 // *** empty log message ***
465 //
466 // Revision 1.4 2003/02/21 15:31:38 balewski
467 // do not kill the chain (it is against my will, JB)
468 //
469 // Revision 1.3 2003/02/20 21:27:06 zolnie
470 // added simple geometry class
471 //
472 // Revision 1.2 2003/02/20 20:13:20 balewski
473 // fixxy
474 // xy
475 //
476 // Revision 1.1 2003/02/20 05:14:07 balewski
477 // reorganization
478 //
479 // Revision 1.2 2003/02/14 00:04:39 balewski
480 // remove few printouts
481 //
482 // Revision 1.1 2003/01/28 23:16:07 balewski
483 // start
484 //
485 // Revision 1.21 2002/12/17 19:43:31 balewski
486 // remove some dependecies, choose better name for maker
487 //
488 // Revision 1.20 2002/12/10 19:02:33 balewski
489 // fix sector ID according to SN #299
490 // tower: volID contains absolute subsectorID in the range 1-60, subSec #1 is the first subSec of sector 1, and advances clockwise
491 // SMD : volID contains absolute sector ID from 1-12,
492 // it does not matter what geometry has been used (i.e. 4 sectors, lower half or both halfs)
493 // The half ID info in volID is redundant now
494 //
495 // Revision 1.19 2002/12/05 14:21:58 balewski
496 // cleanup, time stamp corrected
497 //
498 // Revision 1.18 2002/12/04 15:03:36 balewski
499 // fix: read also just tower or just SMD data
500 //
501 // Revision 1.17 2002/11/13 21:53:51 zolnie
502 // restored "private" DetectorID in EEmcMCData
503 //
504 // Revision 1.16 2002/11/11 21:22:48 balewski
505 // EEMC added to StEvent
506 //
507 // Revision 1.15 2002/10/03 15:52:25 zolnie
508 // updates reflecting changes in *.g files
509 //
510 // Revision 1.14 2002/10/03 00:30:45 balewski
511 // tof taken away
512 //
513 // Revision 1.13 2002/10/01 14:41:54 balewski
514 // SMD added
515 //
516 // Revision 1.12 2002/10/01 06:03:15 balewski
517 // added smd & pre2 to TTree, tof removed
518 //
519 // Revision 1.11 2002/09/30 21:58:27 zolnie
520 // do we understand Oleg? (the depth problem)
521 //
522 // Revision 1.10 2002/09/30 21:04:02 zolnie
523 // SMD updates
524 //
525 // Revision 1.9 2002/09/30 20:15:55 zolnie
526 // Oleg's geometry updates
527 //
528 // Revision 1.8 2002/09/27 19:05:13 zolnie
529 // EEmcMCData updates
530 //
531 // Revision 1.5 2002/09/25 01:36:13 balewski
532 // fixed TOF in geant
533 //
534 // Revision 1.4 2002/09/24 22:47:35 zolnie
535 // major rewrite: SMD incorporated, use constants rather hard coded numbers
536 // introducing exceptions (rather assert)
537 //
538 // Revision 1.3 2002/09/20 21:58:13 balewski
539 // sum of MC hits over activ detectors
540 // produce total tower energy with weight 1 1 1 1
541 //
542 // Revision 1.2 2002/09/20 15:49:05 balewski
543 // add event ID
544 //
545 // Revision 1.1.1.1 2002/09/19 18:58:54 zolnie
546 // Imported sources
547 //
548 // Revision 1.1.1.1 2002/08/29 19:32:01 zolnie
549 // imported sources
550 // Revision 1.2 2002/08/28 01:42:59 zolnie
551 // version alpha ready
552 // Revision 1.1 2002/08/27 16:39:10 zolnie
553 // Initial revision
554 //
555 //
Generated by   doxygen 1.8.5