StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
EEsmdCal.cxx
1 // $Id: EEsmdCal.cxx,v 1.20 2009/02/04 20:33:22 ogrebeny Exp $
2 
3 #include <assert.h>
4 #include <stdlib.h>
5 #include <string.h>
6 
7 #include <TClonesArray.h>
8 #include <TObjArray.h>
9 #include <TH1.h>
10 #include <TH2.h>
11 #include <TFile.h>
12 
13 #include "EEsmdCal.h"
14 #include "EEsmdPlain.h"
15 #include "StEEmcUtil/EEmcGeom/EEmcGeomSimple.h"
16 #include "StEEmcUtil/StEEmcSmd/EEmcSmdGeom.h"
17 #include "StEEmcUtil/EEmcSmdMap/EEmcSmdMap.h"
18 #include "StEEmcUtil/database/cstructs/eemcConstDB.hh"
19 #include "StEEmcUtil/database/EEmcDbItem.h"
20 
21 #ifdef StRootFREE
22  #include "EEmcDb/EEmcDb.h"
23 #else
24  #include "StEEmcUtil/database/StEEmcDb.h"
25 #endif
26 
27 
28 ClassImp(EEsmdCal)
29 //--------------------------------------------------
30 //--------------------------------------------------
31 EEsmdCal::EEsmdCal(){
32  nInpEve=0;
33  HList=0;
34  eeDb=0;
35 
36  // clear all histo pointers just in case
37  memset(hT,0,sizeof(hT));
38  memset(hSs,0,sizeof(hSs));
39  // memset(hSp,0,sizeof(hSp));
40 
41  dbMapped=-1;
42  memset(dbT,0,sizeof(dbT));
43  memset(dbS,0,sizeof(dbS));
44 
45  // initialization
46  smdHitPl=new EEsmdPlain [MaxSmdPlains];
47  geoTw=new EEmcGeomSimple;
48  geoSmd= EEmcSmdGeom::instance();
49  mapSmd= EEmcSmdMap::instance();
50 
51  thrMipSmdE=-1; emptyStripCount=-2;
52  twMipRelEneLow=-3; twMipRelEneHigh=-4;
53  offCenter=0.7; thrMipPresAdc=-5;
54 
55  maxStripAdc=120; // suppress large jump in ped or sticky bits
56 
57  // chose which stat bits are fatal
58  killStat=EEMCSTAT_ONLPED | EEMCSTAT_HOTSTR ;
59 
60  printf("EEsmdCal() constructed\n");
61 }
62 
63 //--------------------------------------------------
64 //--------------------------------------------------
65 EEsmdCal::~EEsmdCal() {/* noop */}
66 
67 
68 //-------------------------------------------------
69 //-------------------------------------------------
70 void EEsmdCal::init( ){
71  printf("EEsmdCal() init , calibrate sector=%d\n",sectID);
72 
73  initAuxHisto();
74 
75  initTileHistoAdc('a',"inclusive ADC ",kBlack);
76  initTileHistoAdc('d',"ADC, tag: UxV 2*thr",kBlack);
77  initTileHistoAdc('e',"ADC, tag: mipT UxV 2*thr",kBlue);
78 
79  initTileHistoEne('f',"energy, tag: UxV 2*thr",kBlack);
80  initTileHistoEne('g',"energy, tag: mipT UxV 2*thr",kBlue);
81 
82  initSmdHist('a',"inclusive ADC");
83  initSmdHist('b',"ADC, tag: best MIP",kBlue);
84  initSmdEneHist('e',"Energy (K)+(K+1)*tgh(eta): best MIP",kBlack);
85 
86  //.................... initialize MIP finding algo for SMD
87  int i;
88  for(i=0;i<MaxSmdPlains;i++) {
89  smdHitPl[i].set(thrMipSmdE,emptyStripCount,i+'U');
90  }
91 
92  printf("use thrMipSmdE/MeV=%.2f emptyStripCount=%d twMipRelEne/high=%.2f/%.2f offCenter=%.2f maxStripAdc=%.1f thrMipPresAdc=%d\n", thrMipSmdE*1000.,emptyStripCount,twMipRelEneLow, twMipRelEneHigh,offCenter,maxStripAdc,thrMipPresAdc);
93  assert(sectID>0 && sectID<=MaxSectors);
94 
95  //....................... initilize MIP energy in towers
96  assert(twMipRelEneLow< twMipRelEneHigh);
97 
98  for(int i=0;i<MaxEtaBins;i++){
99  float etaValue=(geoTw->getEtaMean(i));
100  float tghEta=TMath::TanH(etaValue);
101  twTghEta[i]=tghEta;
102  towerMipE[i]= 1./(2.89*tghEta); // GeV EM
103  presMipE[i]=0.0009/tghEta; //GeV MIP, 5mm*1.8mm
104  }
105 
106  smdAvrMipE=0.0013; //GeV MIP; 7mm*1.8 MeV/cm
107 
108 #if 0 //smdMap - histos for finding mapping
109  int ij;
110  for( ij=0;ij<12;ij++) {
111  char t1[100];
112  sprintf(t1,"hM%c",'a'+ij);
113  hM[ij]=new TH1F(t1,t1,600,0.5,600.5);
114  HList->Add( hM[ij]);
115  }
116 #endif
117 
118 }
119 
120 
121 //-------------------------------------------------
122 //-------------------------------------------------
123 void EEsmdCal::initRun(int runID){
124  printf(" EEsmdCal::initRun(%d)\n",runID);
125 
126 #if 1 //smdMap verification
127  if(dbMapped>0) {
128  printf(" EEsmdCal::initRun(%d) N-th time, Ignore\n",runID);
129  return;
130  }
131 #endif
132 
133  assert(dbMapped<0); // at the moment DB reloading is not implemented/tested,JB
134  mapTileDb();
135  dbMapped=runID;
136 
137  // now init all what relays on DB inofrmation
138  addTwMipEbarsToHisto(kGreen,'g');
139 
140  addPresMipEbarsToHisto(kGreen,'P');
141  addPresMipEbarsToHisto(kGreen,'Q');
142  addPresMipEbarsToHisto(kGreen,'R');
143 
144  addSmdMipEbarsToHisto(kGreen,'U');
145  addSmdMipEbarsToHisto(kGreen,'V');
146  histoGains();
147 }
148 
149 //--------------------------------------------------
150 //--------------------------------------------------
151 void EEsmdCal:: mapTileDb(){
152  printf("EEsmdCal:: mapTileDb()\n");
153 
154  //.................... Tower like ....................
155  char cT[mxTile]={'T','P','Q','R'};
156  int iT=0;
157  for(iT=0;iT<mxTile;iT++) {
158  for(char iSub=0; iSub<MaxSubSec; iSub++){
159  for(int iEta=0; iEta<MaxEtaBins; iEta++){
160  int iPhi=iSect*MaxSubSec+iSub;
161  dbT[iT][iEta][iPhi]=eeDb->getTile(sectID,iSub+'A',iEta+1,cT[iT]);
162  //printf("%d %d %d '%s' \n",iT,iEta,iPhi,dbT[iT][iEta][iPhi]->name);
163  }
164  }
165  }
166 
167  //................ SMD strips
168  int iu;
169  for(iu=0;iu<MaxSmdPlains;iu++) {
170  int istr;
171  for(istr=0;istr<MaxSmdStrips;istr++) {
172  if(istr>=288) break; // ugly, but I gave up, JB
173  dbS[iu][istr]=eeDb->getByStrip(sectID,iu+'U',istr+1);
174  }
175  }
176 
177 }
178 
179 
180 //-------------------------------------------------
181 //-------------------------------------------------
182 void EEsmdCal::clear(){ // called for every event
183  memset(tileAdc,0,sizeof(tileAdc));
184  memset(tileEne,0,sizeof(tileEne));
185  memset(tileThr,0,sizeof(tileThr));
186  memset(smdAdc,0,sizeof(smdAdc));
187  memset(smdEne,0,sizeof(smdEne));
188 
189  memset(killT,true,sizeof(killT));// default is dead
190 
191  int i;
192  for(i=0;i<MaxSmdPlains;i++) {
193  smdHitPl[i].clear();
194  }
195 
196 }
197 //-------------------------------------------------
198 //-------------------------------------------------
199 void EEsmdCal::findSectorMip( ){
200  /*...................................
201  .... ....
202  .... main physics analysis ....
203  ................................... */
204 
205  hA[9]->Fill(1);
206  fillSmdHisto_a(); // inclusive SMD histos
207 
208  for(char iSub=0; iSub<MaxSubSec; iSub++){
209  for(int iEta=0; iEta<MaxEtaBins; iEta++){
210  int iPhi=iSect*MaxSubSec+iSub;
211  fillOneTailHisto('a', iEta,iPhi); // inclusive histos
212  }
213  }
214 
215  // searching for MIP in SMD
216  int ret=getUxVmip();
217 
218  if(ret>1) hA[9]->Fill(2);// counts multiple MIP's per both planes
219  if(ret==1) hA[9]->Fill(3);// counts multiple MIP's per both planes
220 
221  int nU,nV;
222  EEsmdPlain *plU=smdHitPl+0;
223  EEsmdPlain *plV=smdHitPl+1;
224 
225 
226 #if 0 //smdMap verification
227  // verify mapping for a subset of strips
228  // U-plan must have MIP
229  if(plU->nMatch>0) {
230  int ist;
231  for(ist=0;ist<12;ist++)
232  scanSpike(smdAdc[1][30+ist-1],hM[ist]);
233  }
234 #endif
235 
236  // calibrate P,Q,R with UxV
237  for(nU=0;nU<plU->nMatch;nU++){
238  for(nV=0;nV<plV->nMatch;nV++){
239  hA[9]->Fill(4);// any UxV pair
240  calibAllwithMip(plU->iStrip[nU],plV->iStrip[nV]);
241  }
242  }
243 
244  return;
245 }
246 
247 
248 //-------------------------------------------------
249 //-------------------------------------------------
250 void EEsmdCal::calibAllwithMip(int iStrU, int iStrV){
251 
252  // find MIP for a UxV pair of strips
253 
254  //printf("\n jj iSect=%d,iStrU=%d,iStrV=%d\n",iSect,iStrU,iStrV);
255  TVector3 r=geoSmd->getIntersection (iSect,iStrU,iStrV);
256 
257  //printf(" UxV = %f %f %f\n", r.x(),r.y(),r.z());
258 
259  int iSecX, iSubX, iEtaX;
260  Float_t dphi, deta;
261  int ret=geoTw->getTower(r, iSecX, iSubX, iEtaX,dphi, deta);
262  //printf("ret=%d, isecX=%d isubX=%d, ietaX=%d dphi=%f, deta=%f\n",ret,iSecX, iSubX, iEtaX, dphi, deta);
263 
264  if(ret==0 || iSecX!=iSect) return;
265  //................ UxV is in a tower boundary within selected sector
266 
267  hA[9]->Fill(5);
268  hA[10]->Fill(iStrU+1);
269  hA[11]->Fill(iStrV+1);
270 
271  // select central MIP in tower
272  int inCenter=(fabs(dphi)<offCenter && fabs(deta)<offCenter);
273  ((TH2F*) hA[21])->Fill( r.x(),r.y());
274 
275  if(!inCenter) return;
276  //...................... assure central UxV hit in a tower
277 
278  hA[9]->Fill(6);
279  int iPhiX=iSect*MaxSubSec+iSubX;
280  hA[24]->Fill(iPhiX+iEtaX*MaxPhiBins);
281 
282  //................ auxuiliary variables
283  // ped corrected adc
284  float adcT=tileAdc[kT][iEtaX][iPhiX];
285  float adcP=tileAdc[kP][iEtaX][iPhiX];
286  float adcQ=tileAdc[kQ][iEtaX][iPhiX];
287  float adcR=tileAdc[kR][iEtaX][iPhiX];
288 
289  // calibrated energy
290  float eneT=tileEne[kT][iEtaX][iPhiX]; // GeV
291  float eneP=tileEne[kP][iEtaX][iPhiX]*1000; // MeV
292  float eneQ=tileEne[kQ][iEtaX][iPhiX]*1000; // MeV
293  float eneR=tileEne[kR][iEtaX][iPhiX]*1000; // MeV
294 
295  // logical conditions:
296  // to recover few tiles mark all faild towers as with ADC>thres
297  bool thrP=false, thrQ=false, thrR=false, thrDum=false;
298  bool *thr_p[mxTile]={&thrDum,&thrP,&thrQ,&thrR};
299 
300  int iT;
301  for (iT=kP; iT<=kR;iT++) {// loop over pre/post
302  bool thr=tileThr[iT][iEtaX][iPhiX] ; // .... ADC must be above ped
303  thr= thr || killT[iT][iEtaX][iPhiX]; //.... account for masked pixels
304 #if 0 //.... after pre/post calibration is avaliable
305  float preMipRelEneLow=0.6, preMipRelEneHigh=3.0; // tmp
306  float ene2Mip=tileEne[iT][iEtaX][iPhiX]/presMipE[iEtaX];
307  thr= thr && ene2Mip>preMipRelEneLow && ene2Mip<preMipRelEneHigh;
308 #endif
309  (*thr_p[iT])=thr; // .... record final answer
310  }
311 
312  //printf("bb1 (adc) P=%.1f Q=%.1f R=%.1f iphiX=%d\n",tileAdc[1][iEtaX][iPhiX],tileAdc[2][iEtaX][iPhiX],tileAdc[3][iEtaX][iPhiX],iPhiX);
313  //printf("bb2 (bool) thr P=%d Q=%d R=%d\n",thrP,thrQ,thrR);
314 
315 
316  //printf("bb3 (bool) thr P=%d Q=%d R=%d\n",thrP,thrQ,thrR);
317 
318 
319  // check MIP upper/lower limits
320  float RelTwEne=tileEne[kT][iEtaX][iPhiX]/towerMipE[iEtaX];
321  bool mipT= RelTwEne>twMipRelEneLow && RelTwEne<twMipRelEneHigh;
322  mipT= mipT || killT[kT][iEtaX][iPhiX]; // recover dead tower
323  // printf("iphi=%d ieta=%d Tene=%f mipEne=%f mipT=%d adcT=%.1f\n",iPhiX,iEtaX,tileEne[kT][iEtaX][iPhiX],towerMipE[iEtaX],mipT,tileAdc[kT][iEtaX][iPhiX]);
324 
325 
326  if(thrR) hA[9]->Fill(7);
327 
328  if(mipT )hA[9]->Fill(8);
329 
330  //........... MIP ADC spectra .........................
331  int iCut='d'-'a';
332  if( thrP && thrQ && thrR ) hT[iCut][kT][iEtaX][iPhiX]->Fill(adcT);
333  if( mipT && thrQ && thrR ) hT[iCut][kP][iEtaX][iPhiX]->Fill(adcP);
334  if( mipT && thrP && thrR ) hT[iCut][kQ][iEtaX][iPhiX]->Fill(adcQ);
335  if( mipT && thrP && thrQ ) hT[iCut][kR][iEtaX][iPhiX]->Fill(adcR);
336 
337 
338  iCut='e'-'a';
339  if( mipT && thrP && thrQ && thrR ){
340  hT[iCut][kT][iEtaX][iPhiX]->Fill(adcT);
341  hT[iCut][kP][iEtaX][iPhiX]->Fill(adcP);
342  hT[iCut][kQ][iEtaX][iPhiX]->Fill(adcQ);
343  hT[iCut][kR][iEtaX][iPhiX]->Fill(adcR);
344  }
345 
346  // ................. calibrate (energy) MIP spectra ...............
347  iCut='f'-'a';
348  if( thrP && thrQ && thrR ) hT[iCut][kT][iEtaX][iPhiX]->Fill(eneT);
349  if( mipT && thrQ && thrR ) hT[iCut][kP][iEtaX][iPhiX]->Fill(eneP);
350  if( mipT && thrP && thrR ) hT[iCut][kQ][iEtaX][iPhiX]->Fill(eneQ);
351  if( mipT && thrP && thrQ ) hT[iCut][kR][iEtaX][iPhiX]->Fill(eneR);
352 
353  iCut='g'-'a';
354  if( mipT && thrP && thrQ && thrR ){
355  hT[iCut][kT][iEtaX][iPhiX]->Fill(eneT);
356  hT[iCut][kP][iEtaX][iPhiX]->Fill(eneP);
357  hT[iCut][kQ][iEtaX][iPhiX]->Fill(eneQ);
358  hT[iCut][kR][iEtaX][iPhiX]->Fill(eneR);
359  ((TH2F*) hA[22])->Fill( r.x(),r.y());
360  }
361 
362  //..................... calibration of SMD strips
363  if( mipT && thrP && thrQ && thrR ) {
364 
365  int iuv,i2;
366  const int mx=2; // # of strips from given plain
367  int iStr[MaxSmdPlains];
368  iStr[0]=iStrU; // working pointers
369  iStr[1]=iStrV;
370  float eUV=0;// sum from both plains
371  for(iuv=0;iuv<MaxSmdPlains;iuv++) {
372  float e12=0;
373  for(i2=0;i2<mx;i2++) {
374  int istrip=iStr[iuv]+i2;
375  float adc=smdAdc[iuv][istrip];
376  float ene=smdEne[iuv][istrip]*twTghEta[iEtaX]*1000; // MeV, at normal angle
377  if(adc<2) continue;// drop pedestal, tmp
378  e12+=ene; // sum pairs
379  // re-calibratiop of strips
380  hSs['b'-'a'][iuv][istrip]->Fill(adc);
381  }// end of loop over one plain
382  eUV+=e12;
383  // SMD energy for pairs
384  int istrip1=iStr[iuv];
385  ((TH2F*) hA[20])->Fill(istrip1+1,e12);
386  hSs['e'-'a'][iuv][istrip1]->Fill(e12);
387  hA[14+iuv]->Fill(istrip1+1);
388 
389  }// end of loop over UV plains
390  ((TH2F*) hA[23])->Fill(iEtaX+1,eUV);
391  }
392 
393 }
394 
395 #if 0
396  // SMD light attenuation for pair of strips
397  StructEEmcStrip* bar=geoSmd->getStripPtr(istrip,iuv,iSect);
398  TVector3 r1=bar->end1;
399  TVector3 r2=bar->end2;
400  float dx,dy;
401  if(r1.Eta()<r2.Eta()) {
402  dx=r1.x()-r.x();
403  dy=r1.y()-r.y();
404  } else {
405  dx=r2.x()-r.x();
406  dy=r2.y()-r.y();
407  }
408  float dist=sqrt(dx*dx+dy*dy);
409  int iG=istrip/stripGang;
410  assert(iG>=0 && iG <MaxAt);
411  hSc[iuv][iSubX][iG]->Fill(ene);
412  hSeta[iuv][iSubX][iG]->Fill(r.Eta());
413  hSdist[iuv][iSubX][iG]->Fill(dist);
414 #endif
415 
416 //-------------------------------------------------
417 //-------------------------------------------------
418 int EEsmdCal::getUxVmip(){
419  // printf("\n\n EEsmdCal::getUxVmip() eve=%d\n",nInpEve);
420 
421  int iuv;
422  for(iuv=0;iuv<MaxSmdPlains;iuv++) {
423  EEsmdPlain *pl=smdHitPl+iuv;
424  pl->scanAdc(smdEne[iuv], thrMipSmdE);// use if SMD gains are known
425  //pl->scanAdc(smdAdc[iuv], 15); // use if SMD gains are NOT known
426  pl->findMipPattern();
427  hA[12+iuv]->Fill(pl->nMatch);
428  // pl->print(1);
429  // printf("iuv=%d nM=%d\n",iuv, pl->nMatch);
430  }
431 
432  int ret=smdHitPl[0].nMatch*smdHitPl[1].nMatch;
433 
434  // printf("UxV ret=%d\n",ret);
435  return ret;
436 }
437 
438 
439 #include <TLine.h> // just for drawing
440 #include <TText.h> // just for drawing
441 #include <TCanvas.h> // just for drawing
442 
443 //-------------------------------------------------
444 //-------------------------------------------------
445 void EEsmdCal::finish(int k){
446  printf("\n EEsmdCal::finish(sec=%d) nInpEve=%d \n",sectID,nInpEve);
447 
448  if(k<=0) return;
449  // some test drawing:
450  printf("drawing ...\n");
451  // TFile* f=new TFile("outSec5b/mip05b-8zAB.hist.root");
452 
453  int iuv=1;
454  TString tt="sec"; tt+=sectID;
455  TCanvas *c=new TCanvas(tt,tt,400,400);
456  c->Divide(1,2);
457  // TH2F * h=(TH2F *)f->Get("xy05ct");
458  TH2F * h=(TH2F *)hA[22];
459 
460  for(iuv=0;iuv<2;iuv++) {
461  c->cd(iuv+1);
462  h->Draw("colz");
463  char uv='U'+iuv;
464  int i;
465  for(i=0;i<13;i++) {
466  int strip=i*20+21;
467  StructEEmcStrip* bar=geoSmd->getStripPtr(strip-1,iuv,iSect);
468  // printf("len=%f\n",bar->lenght);
469 
470  TVector3 r1=bar->end1;
471  TVector3 r2=bar->end2;
472 
473  TLine *ln=new TLine(r1.x(),r1.y(),r2.x(),r2.y());
474  ln->Draw();
475  ln->SetLineColor(kRed);
476 
477  if(strip==161 ||strip==201 ||strip==241 ) continue;
478  TString strT=uv; strT+=strip;
479  TText *tx;
480  if(uv=='U')
481  tx=new TText(r2.x()+1,r2.y()-5,strT);
482  else {
483  if(strip>140)
484  tx=new TText(r2.x()-10,r2.y()-8,strT);
485  else
486  tx=new TText(r2.x()-20,r2.y(),strT);
487  }
488  tx->Draw();
489  }
490  char sub;
491  for(sub='A'; sub<='E';sub++) {
492  int i=sub-'A';
493  float x=45 -i*8;
494  float y=-55 -i*5;
495  TString strT=sub;
496  TText *tx=new TText(x,y,strT);
497  tx->Draw();
498  }
499 
500  }
501 
502 
503 }
504 
505 
506 #if 0 //smdMap verification
507 //-------------------------------------------------
508 //-------------------------------------------------
509 void EEsmdCal:: scanSpike(float adc1, TH1F *h){
510  float adcTh=20;
511  if(adc1<adcTh) return;
512  int istrip;
513  for(istrip=0+3;istrip<MaxSmdStrips-3;istrip++) {
514  int iuv;
515  for(iuv=0;iuv<2;iuv++) {
516 #if 0
517  if(smdAdc[iuv][istrip-3]>=adcTh) continue;
518  if(smdAdc[iuv][istrip-2]>=adcTh) continue;
519  if(smdAdc[iuv][istrip+2]>=adcTh) continue;
520  if(smdAdc[iuv][istrip+3]>=adcTh) continue;
521 
522  if(smdAdc[iuv][istrip-1]>=adcTh) continue;
523  if(smdAdc[iuv][istrip+1]>=adcTh) continue;
524 #endif
525  if(smdAdc[iuv][istrip ]< adcTh) continue;
526  h->Fill(1+300*iuv+istrip);
527  }
528  }
529 }
530 #endif
531 
EEsmdCal::iSect
int iSect
calibrate only one sector
Definition: EEsmdCal.h:94
EEsmdCal::sectID
int sectID
no. of input events
Definition: EEsmdCal.h:93
EEsmdCal::smdAdc
float smdAdc[MaxSmdPlains][MaxSmdStrips]
30 deg (only for this sector)
Definition: EEsmdCal.h:105
EEsmdCal::tileAdc
float tileAdc[mxTile][MaxEtaBins][MaxPhiBins]
Definition: EEsmdCal.h:100
EEmcGeomSimple
EEMC simple geometry.
Definition: EEmcGeomSimple.h:22
EEmcSmdGeom::getIntersection
TVector3 getIntersection(Int_t iSec, Int_t iUStrip, Int_t iVStrip, const TVector3 &vertex) const
Definition: EEmcSmdGeom.cxx:628
EEmcSmdGeom::getStripPtr
StructEEmcStrip * getStripPtr(const Int_t iStrip, const Int_t iUV, const Int_t iSec)
return a strip pointer from indices
Definition: EEmcSmdGeom.cxx:375
EEmcGeomSimple::getTower
bool getTower(const TVector3 &r, int &sec, int &sub, int &etabin, Float_t &dphi, Float_t &deta) const
Definition: EEmcGeomSimple.cxx:159
EEmcGeomSimple::getEtaMean
Float_t getEtaMean(UInt_t eta) const
Definition: EEmcGeomSimple.h:85
Generated by   doxygen 1.8.5