StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
MuEzSmdCalMaker.cxx
1 // *-- Author : Jan Balewski
2 //
3 // $Id: MuEzSmdCalMaker.cxx,v 1.7 2009/02/04 20:33:22 ogrebeny Exp $
4 
5 #include <TFile.h>
6 #include <TH1.h>
7 #include <TH2.h>
8 #include <StMessMgr.h>
9 
10 #include "MuEzSmdCalMaker.h"
11 
12 #include "StMuDSTMaker/COMMON/StMuEvent.h"
13 #include "StMuDSTMaker/COMMON/StMuDst.h"
14 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
15 
16 #include "StMuDSTMaker/EZTREE/EztEventHeader.h"
17 #include "StMuDSTMaker/EZTREE/EztTrigBlob.h"
18 #include "StMuDSTMaker/EZTREE/EztEmcRawData.h"
19 #include "StMuDSTMaker/EZTREE/StTriggerDataMother.h"
20 //tmp
21 #include "StTriggerData2005.h" // tmp
22 
23 #include "StEEmcUtil/database/StEEmcDb.h"
24 #include "StEEmcUtil/database/EEmcDbItem.h"
25 
26 ClassImp(MuEzSmdCalMaker)
27 
28 //________________________________________________
29 //________________________________________________
30 MuEzSmdCalMaker::MuEzSmdCalMaker( const char* self ,const char* muDstMakerName) : StMaker(self){
31  mMuDstMaker = (StMuDstMaker*)GetMaker(muDstMakerName);
32  assert(mMuDstMaker);
33 
34  trgAkio=0;
35  nAcceptEve=nTrigEve=nCorrEve=0;
36  setHList(0);
37  setTrigIdFilter(0);
38  setMaxCtbSum(0);
39  setEZtree();
40 }
41 
42 //________________________________________________
43 //________________________________________________
44 void MuEzSmdCalMaker::setSector(int sec){
45  EEsmdCal::setSector(sec);
46  TString name=GetName();
47  name+="-";
48  name+=sec;
49  // printf("change name to %s sec=%d\n", name.Data(),sec);
50  SetName(name);
51 }
52 
53 //___________________ _____________________________
54 //________________________________________________
55 MuEzSmdCalMaker::~MuEzSmdCalMaker(){
56  delete trgAkio;
57 }
58 
59 //___________________ _____________________________
60 //________________________________________________
61 void
62 MuEzSmdCalMaker::saveHisto(TString fname){
63  TString outName=fname+".hist.root";
64  TFile f( outName,"recreate");
65  assert(f.IsOpen());
66  printf("%d histos are written to '%s' ...\n",HList->GetEntries(),outName.Data());
67 
68  HList->Write();
69  f.Close();
70 
71 }
72 
73 //________________________________________________
74 //________________________________________________
75 Int_t
76 MuEzSmdCalMaker::Init(){
77 
78  assert(HList);
79  eeDb = (StEEmcDb*)this->GetDataSet("StEEmcDb");
80  assert(eeDb);
81  EEsmdCal::init();
82 
83  gMessMgr->Message("","I") <<GetName()<<"::Init() filter trigID="<<trigID<<" maxCtbSum="<<maxCtbSum<<endm;
84  return StMaker::Init();
85 }
86 
87 //________________________________________________
88 //________________________________________________
89 Int_t
90 MuEzSmdCalMaker::InitRun(int runNo){
91  if(runNo==0) {
92  gMessMgr->Message("","W")<<GetName()<<"::InitRun("<<runNo<<") ??? changed to 555, it s OK for M-C - perhaps, JB"<<endm;
93  runNo=555;
94  }
95 
96  initRun(runNo);
97 
98  return kStOK;
99 }
100 
101 //________________________________________________
102 //________________________________________________
103 Int_t
104 MuEzSmdCalMaker::Finish(){
105  finish(1);// do not draw
106  gMessMgr->Message("","I") <<GetName()<<"::Finish()\n inputEve="<<nInpEve<<" trigFilterEve="<<nTrigEve<<" nCorrEve="<<nCorrEve<<" nAcceptEve="<<nAcceptEve<<endm;
107  return kStOK;
108 }
109 
110 //________________________________________________
111 //________________________________________________
112 void
113 MuEzSmdCalMaker::Clear(const Option_t*){
114  eHead=0;
115  eETow=0;
116  eESmd=0;
117  eTrig=0;
118  // delete trgAkio; //JAN: perhaps is a memory leak? But crashes
119 }
120 
121 //________________________________________________
122 //________________________________________________
123 Int_t
124 MuEzSmdCalMaker::Make(){
125  clear();
126  nInpEve++;
127 // gMessMgr->Message("","D") <<GetName()<<"::Make() is called ,
128 // useEZtree="<<useEZtree<<endm;
129  if(useEZtree) return MakeEZtree();
130 
131  return MakeRegular();
132 }
133 
134 //________________________________________________
135 //________________________________________________
136 Int_t
137 MuEzSmdCalMaker::MakeEZtree(){
138  //.......... acquire EztHeader
139  eHead= mMuDstMaker->muDst()->eztHeader();
140  if(eHead==0) {
141  gMessMgr->Message("","E") <<GetName()<<"::Make() no EztEventHeader, skip event "<<endm; return kStOK;
142  }
143 
144  if(nInpEve==1) eHead->print();
145 
146  if(trigID ) {// filter by triggerID on demand
147  if (! mMuDstMaker->muDst()->event()->triggerIdCollection().nominal().isTrigger(trigID)) return kStOK;
148  }
149  nTrigEve++;
150  //.... get data .....
151  eETow=mMuDstMaker->muDst()->eztETow();
152  eESmd=mMuDstMaker->muDst()->eztESmd();
153  eTrig=mMuDstMaker->muDst()->eztTrig();
154  // printf("pp %p %p %p\n",eETow, eESmd, eTrig);
155  if(!eETow || !eESmd || !eTrig) return kStOK;
156 
157  // trgAkio=new StTriggerDataMother(eTrig);
158  // trgAkio->dump();
159  //eETow->print(0);
160  // eESmd->print(0);
161  // eHead->print();
162 
163  StMuEvent *muEve = mMuDstMaker -> muDst() -> event();
164  assert(muEve);
165  StEventInfo &info=muEve->eventInfo();
166  int runId=info.runId();
167 
168 
169  // .... process adata ......
170  void *blob=eTrig->trgd->GetArray();
171 
172  StTriggerData2005 trgAkio5( (const TrgDataType2005 *)blob,runId);
173  if(eETow->doTowerHeadCorruptionTest(trgAkio5.token())||
174  eESmd->doMapmtHeadCorruptionTest(trgAkio5.token())
175  ) {
176  nCorrEve++;
177  return kStOK;
178  }
179 
180  int ctbSum=trgAkio5.ctbSum();
181  if(maxCtbSum>0 && (ctbSum>maxCtbSum || ctbSum<maxCtbSum/2.)) return kStOK;
182  nAcceptEve++;
183 
184  // printf("\n ctbSum=%d \n",ctbSum);
185  //x hA[7]->Fill(ctbSum);
186 
187 
188  unpackMuEzt(eETow);
189  unpackMuEzt(eESmd);
190 
191  findSectorMip();// do real analysis
192 
193  return kStOK;
194 }
195 
196 //________________________________________________
197 //________________________________________________
198 Int_t
199 MuEzSmdCalMaker::MakeRegular(){
200  gMessMgr->Message("","D") <<GetName()<<"::MakeRegular() is called , useEZtree="<<useEZtree<<endm;
201 
202 #if 0
203  vector<unsigned int> trgL=mMuDstMaker->muDst()->event()->triggerIdCollection().nominal().triggerIds();
204  printf("trigL len=%d\n",trgL.size());
205  int ii;
206  for(ii=0;ii<trgL.size();ii++) printf("ii=%d trigID=%d\n",ii,trgL[ii]);
207 #endif
208 
209  if(trigID ) {// filter by triggerID on demand
210  if (! mMuDstMaker->muDst()->event()->triggerIdCollection().nominal().isTrigger(trigID)) return kStOK;
211  }
212  nTrigEve++;
213 
214  // M-C case, when muDst data are zero supressed
215  memset(killT,false,sizeof(killT));// default is working if zero suppressed
216 
217 
218  unpackMuTails();
219  unpackMuSmd();
220  nAcceptEve++;
221 
222  // printf("UU (bool) killThr T=%d P=%d Q=%d R=%d\n",killT[0][3][14],killT[1][3][14],killT[2][3][14],killT[3][3][14]);
223 
224  findSectorMip();// do real analysis
225 
226  return kStOK;
227 }
228 
229 
230 //________________________________________________
231 //________________________________________________
232 void
233 MuEzSmdCalMaker::tileReMap( int &iT,int &sec , char &sub , int &eta){
234  assert(1==2); //use only by expert - disabled, JB
235  if(sec==1 ) {
236  if (iT==99 ) { ;
237  } else if( iT==2 && sub=='A' && eta==11 ) { // QA11<==>QB2
238  iT=2; sub='B'; eta=2;
239  } else if( iT==2 && sub=='B' && eta==2 ) {
240  iT=2; sub='A'; eta=11;
241  } else if( iT==0 && sub=='A' && eta==4 ) { //TA4 <==> TA5
242  iT=0; sub='A'; eta=5;
243  } else if( iT==0 && sub=='A' && eta==5 ) {
244  iT=0; sub='A'; eta=4;
245  }
246  }
247  return;
248 }
249 
250 //________________________________________________
251 //________________________________________________
252 int
253 MuEzSmdCalMaker:: stripReMap(const EEmcDbItem *x){
254  assert(1==2); //use only by expert - disabled, JB
255  int str=x->strip;
256  if(x->sec==8 && x->plane=='V' ) {
257  switch( x->strip) {
258  // connector & pair swap
259  case 209: str=216; break;
260  case 210: str=215; break;
261  case 211: str=214; break;
262  case 212: str=213; break;
263  case 213: str=212; break;
264  case 214: str=211; break;
265  case 215: str=210; break;
266  case 216: str=280; break;
267  case 280: str=209; break;
268 
269  // another connector
270  case 265: str=272; break;
271  case 266: str=271; break;
272  case 267: str=270; break;
273  case 268: str=269; break;
274  case 269: str=268; break;
275  case 270: str=267; break;
276  case 271: str=266; break;
277  case 272: str=265; break;
278  }
279  }
280 
281  return str;
282 }
283 
284 //________________________________________________
285 //________________________________________________
286 void
287 MuEzSmdCalMaker::unpackMuEzt(EztEmcRawData *eRaw){
288  int n1=0,n2=0,n3=0;
289 
290  if(eRaw==0) return ; // no data block
291  int icr;
292  for(icr=0;icr<eRaw->getNBlocks();icr++) {
293  if(eRaw->isCrateVoid(icr)) continue;
294  assert(!eRaw ->getCorruption(icr)); // zero-tolerance
295 
296  int crateID=eRaw->getCrateID(icr);
297  int chan;
298  const UShort_t* data=eRaw->data(icr);
299  for(chan=0;chan<eRaw->sizeData(icr);chan++) {
300  const EEmcDbItem *x=eeDb->getByCrate(crateID,chan);
301  if(x==0) continue;
302  if(x->sec!=sectID && crateID>6 ) break;// assumes crates do not cross sectors for non-towers, faster code
303  if(x->fail ) continue; // drop broken channels
304  if(x->stat & killStat) continue; // drop masked channels
305  // accept this hit
306  float rawAdc=data[chan];
307  float adc=rawAdc-x->ped;
308 
309  if(x->isSMD()) {
310  //........................ SMD U or V .......
311  if(rawAdc>x->thr) n3++;
312  int iuv=x->plane-'U';
313  int istr=x->strip -1;
314 
315  // istr=stripReMap(x)-1; //<<==== S W A P S, not use it
316 
317  assert(iuv>=0 && iuv<MaxSmdPlains);
318  assert(istr>=0 && istr<MaxSmdStrips);
319  smdAdc[iuv][istr]=adc;
320  if(x->gain<=0)continue; // drop channels w/o gains
321 
322  smdEne[iuv][istr]=adc/x->gain;
323  // if(rawAdc>x->thr) printf("%s %f %f \n",x->name,smdEne[iuv][istr],x->gain);
324  } else {
325  //............................... tower/pre/post crates
326  int iT=-1;// store T,P,Q,R depending on 'iT'
327  if(x->name[2]=='T'){
328  iT=0;
329  } else{
330  iT=x->name[2]-'P'+1;
331  }
332  assert(iT>=0 && iT<mxTile);
333  bool aboveThr=rawAdc>x->thr;
334  if(iT==1 || iT==2) {
335  if( adc<=thrMipPresAdc || adc>(thrMipPresAdc+100) ) aboveThr=false;
336  } else if (iT==3) {
337  if( adc<=(thrMipPresAdc/2.) || adc>(thrMipPresAdc/2.+100) ) aboveThr=false;
338  }
339 
340  // if(iT==1 || iT==3) continue; // mask alomst all
341  int sec=x->sec;
342  char sub=x->sub;
343  int eta=x->eta;
344 
345  // tileReMap( iT,sec,sub,eta); //<<==== S W A P S , not use it
346 
347  int iphi=(sec-1)*MaxSubSec+(sub-'A');
348  int ieta=eta-1;
349  assert(iphi>=0 && iphi<MaxPhiBins);
350  assert(ieta>=0 && ieta<MaxEtaBins);
351  tileAdc[iT][ieta][iphi]=adc;
352  tileThr[iT][ieta][iphi]=aboveThr;
353  killT[iT][ieta][iphi]=false; // it is alive
354  if(aboveThr) {
355  if(iT==0)
356  n1++;
357  else
358  n2++;
359  }
360 
361  if(x->gain<=0) continue;// drop channels w/o gains
362  tileEne[iT][ieta][iphi]=adc/x->gain;
363  }
364  } // end of loop over one data block
365  }// end of loop over blocks
366  // printf("%s-->nTow=%d nPQR=%d nSMD=%d\n",GetName(),n1,n2,n3);
367 }
368 
369 //________________________________________________
370 //________________________________________________
371 void
372 MuEzSmdCalMaker::unpackMuTails(){
373 
374  // Access to muDst .......................
375  StMuEmcCollection* emc = mMuDstMaker->muDst()->muEmcCollection();
376  if (!emc) {
377  gMessMgr->Warning() <<"No EMC data for this event"<<endm; return;
378  }
379 
380  int i ;
381  //printf("aaa %d %d \n",eeDb->mfirstSecID,eeDb->mlastSecID);
382  //......................... T O W E R S .....................
383  for (i=0; i< emc->getNEndcapTowerADC(); i++) {
384  int sec,eta,sub,rawAdc; //muDst ranges:sec:1-12, sub:1-5, eta:1-12
385  emc->getEndcapTowerADC(i,rawAdc,sec,sub,eta);
386  if(sec!=sectID) continue;
387  //Db ranges: sec=1-12,sub=A-E,eta=1-12,type=T,P-R ; slow method
388  const EEmcDbItem *x=eeDb->getTile(sec,'A'+sub-1,eta,'T');
389  assert(x); // it should never happened for muDst
390  if(x->fail || x->stat & killStat) {// drop not working channels
391  killTail(x,0);
392  continue;
393  }
394  float adc=rawAdc-x->ped; // ped subtracted ADC
395 
396  //W A R N !! HARDCODDED SF correction
397  // adc/=0.8; //replace SF of 5% used in M_C generation by 4% believed to be true as of October 2005
398 
399  bool aboveThr=rawAdc>x->thr;
400 
401  int iT=0; // towers
402  recordTail( x,adc,aboveThr,iT);
403 
404  }// end of loop over towers
405 
406  //......................... P R E - P O S T .....................
407  int pNh= emc->getNEndcapPrsHits();
408  for (i=0; i < pNh; i++) {
409  int pre, sec,eta,sub;
410  //muDst ranges: sec:1-12, sub:1-5, eta:1-12 ,pre:1-3==>pre1/pre2/post
411  StMuEmcHit *hit=emc->getEndcapPrsHit(i,sec,sub,eta,pre);
412  float rawAdc=hit->getAdc();
413  //Db ranges: sec=1-12,sub=A-E,eta=1-12,type=T,P-R ; slow method
414  const EEmcDbItem *x=eeDb-> getTile(sec,sub-1+'A', eta, pre-1+'P');
415  assert(x); // it should never happened for muDst
416  if(sec!=sectID) continue;
417  int iT=pre;// store P,Q,R depending on 'iT'
418  assert(iT>=0 && iT<mxTile);
419  if(x->fail || x->stat & killStat) {// drop not working channels
420  killTail(x,iT);
421  continue;
422  }
423 
424  float adc=rawAdc-x->ped; // ped subtracted ADC
425 
426  bool aboveThr=rawAdc > x->thr;
427  //pre1,2
428  if( adc<=thrMipPresAdc || adc>(thrMipPresAdc+100) ) aboveThr=false;
429  if(iT==33){ //post, tmp for M-C ADC is 2x larger
430  if( adc<=(thrMipPresAdc/2.) || adc>(thrMipPresAdc/2.+100) ) aboveThr=false;
431  }
432 
433  recordTail( x,adc,aboveThr,iT);
434  }
435 
436  return ;
437 }
438 
439 //________________________________________________
440 //________________________________________________
441 void
442 MuEzSmdCalMaker::unpackMuSmd(){
443 
444  // Access to muDst .......................
445  StMuEmcCollection* emc = mMuDstMaker->muDst()->muEmcCollection();
446  if (!emc) {
447  gMessMgr->Warning() <<"No EMC data for this event"<<endm; return;
448  }
449 
450  //....................... S M D ................................
451  char uv='U';
452  for(uv='U'; uv <='V'; uv++) {
453  int sec,strip;
454  int nh= emc->getNEndcapSmdHits(uv);
455  int i;
456  for (i=0; i < nh; i++) {
457  StMuEmcHit *hit=emc->getEndcapSmdHit(uv,i,sec,strip);
458  float rawAdc=hit->getAdc();
459  const EEmcDbItem *x=eeDb->getByStrip(sec,uv,strip);
460  assert(x); // it should never happened for muDst
461  if(sec!=sectID) continue;
462 
463  if(x->fail ) continue; // drop broken channels
464  if(x->stat & killStat) continue; // drop not working channels
465  float adc=rawAdc-x->ped; // ped subtracted ADC
466  // x->print(); printf("adc=%f\n",adc);
467  int iuv=x->plane-'U';
468  int istr=x->strip -1;
469 
470  assert(iuv>=0 && iuv<MaxSmdPlains);
471  assert(istr>=0 && istr<MaxSmdStrips);
472  smdAdc[iuv][istr]=adc;
473  if(x->gain<=0)continue; // drop channels w/o gains
474 
475  smdEne[iuv][istr]=adc/x->gain;
476  }
477  }
478 }
479 
480 //________________________________________________
481 //________________________________________________
482 void
483 MuEzSmdCalMaker::recordTail( const EEmcDbItem *x, float adc, bool aboveThr, int iT) {
484  int sec=x->sec;
485  char sub=x->sub;
486  int eta=x->eta;
487 
488 #if 0
489  if(strstr(x->name,"C10")) {
490  printf("\nadc=%f iT=%d abTthr=%d \n",adc,iT,aboveThr); x->print();
491  }
492 #endif
493 
494  // tileReMap( iT,sec,sub,eta); //<<==== S W A P S , not use it
495 
496  int iphi=(sec-1)*MaxSubSec+(sub-'A');
497  int ieta=eta-1;
498  assert(iphi>=0 && iphi<MaxPhiBins);
499  assert(ieta>=0 && ieta<MaxEtaBins);
500  tileAdc[iT][ieta][iphi]=adc;
501  tileThr[iT][ieta][iphi]=aboveThr;
502  killT[iT][ieta][iphi]=false; // it is alive
503 
504  if(x->gain<=0) return;// drop channels w/o gains
505  tileEne[iT][ieta][iphi]=adc/x->gain;
506 }
507 
508 
509 //________________________________________________
510 //________________________________________________
511 void
512 MuEzSmdCalMaker::killTail( const EEmcDbItem *x, int iT) {
513  int sec=x->sec;
514  char sub=x->sub;
515  int eta=x->eta;
516 
517  int iphi=(sec-1)*MaxSubSec+(sub-'A');
518  int ieta=eta-1;
519  assert(iphi>=0 && iphi<MaxPhiBins);
520  assert(ieta>=0 && ieta<MaxEtaBins);
521  killT[iT][ieta][iphi]=true; // is dead
522 }
523 
524 //---------------------------------------------------
525 // $Log: MuEzSmdCalMaker.cxx,v $
526 // Revision 1.7 2009/02/04 20:33:22 ogrebeny
527 // Moved the EEMC database functionality from StEEmcDbMaker to StEEmcUtil/database. See ticket http://www.star.bnl.gov/rt2/Ticket/Display.html?id=1388
528 //
529 // Revision 1.6 2007/08/21 13:10:04 balewski
530 // final, used in 2006 offline calibration by soScott
531 //
532 //
533 // VS: ----------------------------------------------------------------------
534 //
535 // Revision 1.5 2006/09/15 01:45:34 balewski
536 // add run# to trg-data unpaker
537 //
538 // Revision 1.4 2005/09/29 13:57:57 balewski
539 // after SMD gains were rescaled
540 //
541 // Revision 1.3 2005/08/09 18:46:31 balewski
542 // after smd calib in 2005
543 //
544 // Revision 1.2 2005/05/04 17:00:32 balewski
545 // tuned for MIP detection in CuCu200
546 //
547 // Revision 1.1 2005/03/11 15:44:25 balewski
548 // works with muEzt, cucu200
549 //
550 
StMuDstMaker::muDst
StMuDst * muDst()
Definition: StMuDstMaker.h:425
EEsmdCal::sectID
int sectID
no. of input events
Definition: EEsmdCal.h:93
EEmcDbItem::name
char name[StEEmcNameLen]
ASCII name of the channel, see Readme.
Definition: EEmcDbItem.h:20
StMuDst::eztESmd
static EztEmcRawData * eztESmd()
returns pointer to eztESmd +pre/post
Definition: StMuDst.h:459
MuEzSmdCalMaker::Finish
virtual Int_t Finish()
Definition: MuEzSmdCalMaker.cxx:104
MuEzSmdCalMaker::Make
virtual Int_t Make()
Definition: MuEzSmdCalMaker.cxx:124
StMuDst::eztTrig
static EztTrigBlob * eztTrig()
returns pointer to eztTrig
Definition: StMuDst.h:448
EEsmdCal::smdAdc
float smdAdc[MaxSmdPlains][MaxSmdStrips]
30 deg (only for this sector)
Definition: EEsmdCal.h:105
StMuEmcHit::getAdc
int getAdc() const
Return ADC value.
Definition: StMuEmcHit.h:19
StMuDst::eztHeader
static EztEventHeader * eztHeader()
returns pointer to eztHeader
Definition: StMuDst.h:442
EEsmdCal::killStat
unsigned short killStat
DB access point.
Definition: EEsmdCal.h:113
StMuDst::eztETow
static EztEmcRawData * eztETow()
returns pointer to ETOW
Definition: StMuDst.h:456
EEsmdCal::setSector
void setSector(int x)
the same info, counted from 0
Definition: EEsmdCal.h:95
StMuDst::muEmcCollection
static StMuEmcCollection * muEmcCollection()
returns pointer to current StMuEmcCollection
Definition: StMuDst.h:389
StMuDst::event
static StMuEvent * event()
returns pointer to current StMuEvent (class holding the event wise information, e.g. event number, run number)
Definition: StMuDst.h:320
StMaker::GetName
virtual const char * GetName() const
special overload
Definition: StMaker.cxx:237
EEsmdCal::tileAdc
float tileAdc[mxTile][MaxEtaBins][MaxPhiBins]
Definition: EEsmdCal.h:100
kStOK
Definition: Stypes.h:40
Generated by   doxygen 1.8.5