StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
EEsmdCalHisto.cxx
1 // $Id: EEsmdCalHisto.cxx,v 1.18 2009/12/03 22:35:03 ogrebeny Exp $
2 
3 #include <assert.h>
4 #include <stdlib.h>
5 
6 #include <TClonesArray.h>
7 #include <TObjArray.h>
8 #include <TH1.h>
9 #include <TH2.h>
10 #include <TLine.h>
11 #include <TFile.h>
12 
13 #include "EEsmdCal.h"
14 #include "StEEmcUtil/database/EEmcDbItem.h"
15 
16 #ifdef StRootFREE
17  #include "EEmcDb/EEmcDb.h"
18 #else
19  #include "StEEmcUtil/database/StEEmcDb.h"
20 #endif
21 
22 //--------------------------------------------------
23 //--------------------------------------------------
24 void EEsmdCal::initTileHistoAdc(char cut, const char *title, int col) {
25  int iCut=cut -'a';
26  assert(iCut>=0 && iCut<kCut);
27  char tt1[100], tt2[500];
28 
29  const char *cTile[mxTile]={"Tower","Pres1","Pres2","Post"};
30  const char cT[mxTile]={'T','P','Q','R'};
31 
32  int iT=0;
33  for(iT=0;iT<mxTile;iT++) {
34  for(char iSub=0; iSub<MaxSubSec; iSub++){
35  for(int iEta=0; iEta<MaxEtaBins; iEta++){
36  char sub=iSub+'A';
37  int eta=iEta+1;
38  int iPhi=iSect*MaxSubSec+iSub;
39  char core[100];
40  sprintf(core,"%02d%c%c%02d",sectID,cT[iT],sub,eta);
41  sprintf(tt1,"%c%s",cut,core);
42  sprintf(tt2,"%s(%c) %s , %s; ADC-ped",cTile[iT],cut,core,title);
43  //printf("tt1=%s, tt2=%s\n",tt1,tt2);
44  TH1F *h;
45  if(cT[iT]=='T')
46  h=new TH1F(tt1,tt2,220,-20,200.);
47  else
48  h=new TH1F(tt1,tt2,400,-50,350.);
49  h->SetLineColor(col);
50  HList->Add(h);
51  hT[iCut][iT][iEta][iPhi]=h;
52  }
53  }
54  }
55 }
56 
57 
58 
59 
60 //--------------------------------------------------
61 //--------------------------------------------------
62 void EEsmdCal::initTileHistoEne(char cut, const char *title, int col) {
63  int iCut=cut -'a';
64  assert(iCut>=0 && iCut<kCut);
65  char tt1[100], tt2[500];
66 
67  const char *cTile[mxTile]={"Tower","Pres1","Pres2","Post"};
68  const char cT[mxTile]={'T','P','Q','R'};
69  const char *cUnits[mxTile]={"(GeV)_EM", "(MEV)_MIP", "(MEV)_MIP", "(MEV)_MIP"};
70 
71  int iT=0;
72  for(iT=0;iT<mxTile;iT++) {
73  for(char iSub=0; iSub<MaxSubSec; iSub++){
74  for(int iEta=0; iEta<MaxEtaBins; iEta++){
75  char sub=iSub+'A';
76  int eta=iEta+1;
77  int iPhi=iSect*MaxSubSec+iSub;
78  char core[100];
79  sprintf(core,"%02d%c%c%02d",sectID,cT[iT],sub,eta);
80  sprintf(tt1,"%c%s",cut,core);
81  sprintf(tt2,"%s(%c) %s , %s; ADC-ped/gain %s",cTile[iT],cut,core,title,cUnits[iT]);
82  //printf("tt1=%s, tt2=%s\n",tt1,tt2);
83  TH1F *h=new TH1F(tt1,tt2,200,-0.5,9.5);
84  h->SetLineColor(col);
85  HList->Add(h);
86  hT[iCut][iT][iEta][iPhi]=h;
87  }
88  }
89  }
90 }
91 
92 
93 
94 //--------------------------------------------------
95 //--------------------------------------------------
96 void EEsmdCal::addTwMipEbarsToHisto (int col, char mxC) {
97  // search all existing tower histo (with '05T' in name) and add
98  // bars for MIP limits
99 
100  assert(dbMapped>0);
101 
102  char core[100];
103  sprintf(core,"%02dT",sectID);
104  float yMax=1000;
105  TIterator *it=HList->MakeIterator();
106  TH1 *h;
107  while( (h=(TH1*) it->Next())) {
108  const char *name=h->GetName();
109  if(strstr(name,core)==0) continue;
110  if(name[0]>mxC) continue; // to skip energy-like histos declared later
111  //printf("%s\n",h->GetName());
112  int iSub=name[4]-'A';
113  int iEta=atoi(name+5)-1;
114  int iPhi=iSect*MaxSubSec+iSub;
115  const EEmcDbItem *x=dbT[kT][iEta][iPhi];
116  assert(x);
117  if(x->gain<=0) continue;
118  TList *L=h->GetListOfFunctions();
119 
120  float adcC=towerMipE[iEta]*x->gain;
121  if(name[0]>'e') adcC=towerMipE[iEta]; //ugly hack,JB, last two not ADC
122 
123  TLine *ln=new TLine(adcC,0,adcC,yMax);
124  ln->SetLineColor(kRed); ln->SetLineStyle(2);
125  L->Add(ln);
126 
127  float adc=adcC*twMipRelEneLow;
128  ln=new TLine(adc,0,adc,yMax);
129  ln->SetLineColor(col);
130  L->Add(ln);
131 
132  adc=adcC*twMipRelEneHigh;
133  ln=new TLine(adc,0,adc,yMax);
134  ln->SetLineColor(col);
135  L->Add(ln);
136  }
137 }
138 
139 //--------------------------------------------------
140 //--------------------------------------------------
141 void EEsmdCal::addPresMipEbarsToHisto (int col, char cT) {
142  // search all existing tower histo (with '05X' in name) and add
143  // bars for MIP limits
144 
145  int iT=1+cT-'P';
146  assert(iT>kT && iT<mxTile);
147  assert(dbMapped>0);
148 
149  char core[100];
150  sprintf(core,"%02d%c",sectID,cT);
151  float yMax=1000;
152  TIterator *it=HList->MakeIterator();
153  TH1 *h;
154  while( (h=(TH1*) it->Next())) {
155  const char *name=h->GetName();
156  if(strstr(name,core)==0) continue;
157  // printf("%s\n",h->GetName());
158  int iSub=name[4]-'A';
159  int iEta=atoi(name+5)-1;
160  int iPhi=iSect*MaxSubSec+iSub;
161  const EEmcDbItem *x=dbT[iT][iEta][iPhi];
162  assert(x);
163  if(x->gain<=0) continue;
164 
165  float adcC=presMipE[iEta]*x->gain;
166  if(name[0]>'e') adcC=presMipE[iEta]*1000; //ugly hack2,JB, last two not ADC
167 
168  TLine *ln=new TLine(adcC,0,adcC,yMax);
169  ln->SetLineColor(col); ln->SetLineStyle(2);
170  TList *L=h->GetListOfFunctions();
171  L->Add(ln);
172  }
173 }
174 
175 
176 //--------------------------------------------------
177 //--------------------------------------------------
178 void EEsmdCal::initSmdHist(char cut, const char *title, int col) {
179  int iCut=cut -'a';
180  assert(iCut>=0 && iCut<kCut);
181  char tt1[100], tt2[500];
182 
183  int iuv,istrip;
184  for(iuv=0;iuv<MaxSmdPlains;iuv++) {
185  for(istrip=0;istrip<MaxSmdStrips;istrip++) {
186  char core[100];
187  // single strip (ADC-ped) spectra
188  sprintf(core,"%02d%c%03d",sectID,iuv+'U',istrip+1);
189  sprintf(tt1,"%c%s",cut,core);
190  sprintf(tt2,"SMD(%c) %s , %s; ADC-ped",cut,core,title);
191  //printf("tt1=%s, tt2=%s\n",tt1,tt2);
192  TH1F *h=new TH1F(tt1,tt2,400,-50,350);
193  //TH1F *h=new TH1F(tt1,tt2,4400,-400,4000);//tmp
194  h->SetLineColor(col);
195 
196  HList->Add(h);
197  hSs[iCut][iuv][istrip]=h;
198  }
199  }
200 
201 }
202 
203 //--------------------------------------------------
204 //--------------------------------------------------
205 void EEsmdCal::initSmdEneHist(char cut, const char *title, int col) {
206  int iCut=cut -'a';
207  assert(iCut>=0 && iCut<kCut);
208  char tt1[100], tt2[500];
209 
210  int iuv,istrip;
211  for(iuv=0;iuv<MaxSmdPlains;iuv++) {
212  for(istrip=0;istrip<MaxSmdStrips;istrip++) {
213  char core[100];
214  sprintf(core,"%02d%c%03d",sectID,iuv+'U',istrip+1);
215  sprintf(tt1,"%c%s",cut,core);
216  sprintf(tt2,"SMD(%c) %s , %s; normal incident MIP energy (MeV)",cut,core,title);
217  //printf("tt1=%s, tt2=%s\n",tt1,tt2);
218  TH1F *h=new TH1F(tt1,tt2,200,-0.1,9.9);
219  h->SetLineColor(col);
220 
221  HList->Add(h);
222  hSs[iCut][iuv][istrip]=h;
223  }
224  }
225 
226 }
227 
228 //--------------------------------------------------
229 //--------------------------------------------------
230 void EEsmdCal::addSmdMipEbarsToHisto (int col, char cU) {
231  // search all existing tower histo (with 'x05U' in name) and add
232  // bars for MIP limits
233 
234  assert(cU=='U' || cU=='V');
235  assert(dbMapped>0);
236 
237  char core[100];
238  sprintf(core,"%02d%c",sectID,cU);
239  float yMax=1000;
240  TIterator *it=HList->MakeIterator();
241  TH1 *h;
242  while( (h=(TH1*) it->Next())) {
243  const char *name=h->GetName();
244  if(strstr(name,core)!=name+1) continue;
245  // printf("%s-%c\n",h->GetName(),h->GetName()[0]);
246  //printf("%s\n",h->GetTitle());
247  int iU=name[3]-'U';
248  int iStr=atoi(name+4)-1;
249  assert(iU>=0 && iU<MaxSmdPlains);
250  assert(iStr>=0 && iStr<MaxSmdStrips);
251  const EEmcDbItem *x=dbS[iU][iStr];
252  // printf("iU=%c iStr=%d\n",iU,iStr);
253  if(x==0) continue;
254 
255  if((h->GetName()[0])=='a') {// one more hack
256  char tt3[500];
257  sprintf(tt3,"%s, tube=%s",h->GetTitle(),x->tube);
258  h->SetTitle(tt3);
259  // printf("%s %p %s\n",h->GetTitle(),h,h->GetName());
260  }
261 
262  if(x->gain<=0) continue;
263  TList *L=h->GetListOfFunctions();
264  float adcC=smdAvrMipE*x->gain;
265  if((h->GetName()[0])>'b') adcC=smdAvrMipE*1000.; // now in MeV
266 
267  TLine *ln=new TLine(adcC,0,adcC,yMax);
268  ln->SetLineColor(col); ln->SetLineStyle(2);
269  L->Add(ln);
270 
271  // printf("%s iU=%c iStr=%d adcC=%f\n",name,iU+'U',iStr,adcC);
272 
273  }
274 }
275 
276 
277 //--------------------------------------------------
278 //--------------------------------------------------
279 void EEsmdCal::histoGains(){
280 
281  //.................... SMD ................
282  int iuv,istrip;
283  for(iuv=0;iuv<MaxSmdPlains;iuv++) {
284  for(istrip=0;istrip<MaxSmdStrips;istrip++) {
285  const EEmcDbItem *x=dbS[iuv][istrip];
286  if(x==0) continue;
287  //if(x->fail) continue; // use any non-zero gain
288  if(x->gain<=0) continue;
289  hA[16+iuv]->Fill(x->strip,x->gain);
290  // dig out MAPMT pixel
291  const char *tube=x->tube+2;
292  assert(tube[0]=='S');
293  int box=atoi(tube+1);
294  int pmt=atoi(tube+3);
295  int pix=atoi(tube+6);
296  int xx=(pmt-1)*16+pix;
297  hA[5+box-1]->Fill(xx,x->gain);
298  // printf("%d %d %d %d", box,pmt,pix,xx); x->print(); assert(1==3);
299  }
300  }
301 
302 
303  //.............. Pre/post
304  int iT,iEta,iPhi;
305  for(iT=kP;iT<=kR;iT++)
306  for(iEta=0;iEta<MaxEtaBins;iEta++)
307  for(iPhi=0;iPhi<MaxPhiBins;iPhi++){
308  const EEmcDbItem *x=dbT[iT][iEta][iPhi];
309  if(x==0) continue;
310  // if(x->fail) continue; // use any non-zero gain
311  if(x->gain<=0) continue;
312  // dig out the MAPMT pixel
313  assert(x->sec==sectID);
314  const char *tube=x->tube+2;
315  assert(tube[0]=='P');
316  int pmt=atoi(tube+3);
317  int pix=atoi(tube+6);
318  int xx=(pmt-1)*16+pix;
319  hA[5+3]->Fill(xx,x->gain);
320  }
321 
322 
323  //.............. Towers
324  iT=0;
325  for(iEta=0;iEta<MaxEtaBins;iEta++)
326  for(iPhi=0;iPhi<MaxPhiBins;iPhi++){
327  const EEmcDbItem *x=dbT[iT][iEta][iPhi];
328  if(x==0) continue;
329  if(x->gain<=0) continue;
330  assert(x->sec==sectID);
331  int ispir=(x->eta-1)*MaxSubSec + x->sub-'A';
332  hA[4]->Fill(ispir,x->gain);
333  }
334 
335 }
336 
337 
338 //--------------------------------------------------
339 //--------------------------------------------------
340 void EEsmdCal::initAuxHisto(){
341  int i;
342  // float Emax=2;
343  memset(hA,0,sizeof(hA));
344 
345  char tt1[100], tt2[500], tt0[100];
346  TH1F *h;
347  TH2F *h2;
348 
349  // tower gains
350  sprintf(tt1,"ug%02dT",sectID);
351  sprintf(tt2,"used tower gains sec=%d ; x=spiral=(eta-1)*5+subs-A; gain [ch/GeV]",sectID);
352  h=new TH1F(tt1,tt2, 60,-0.5,59.5);
353  hA[4]=h;
354 
355  for(i=0;i<4;i++) {
356  // gains vs. MAPMT box/pmt/ch, histo 5,6,7,8
357  if(i<3) {
358  sprintf(tt0,"%02d_S%d",sectID,i+1);
359  } else {
360  sprintf(tt0,"%02d_P1",sectID);
361  }
362  sprintf(tt1,"ug%s",tt0);
363  sprintf(tt2,"used gains MAPMT box %s ; chann=(pmt ID-1)*16+pix ; gain [ch/GeV]",tt0);
364 
365  h=new TH1F(tt1,tt2, MaxMapmtCrateCh,0.5, MaxMapmtCrateCh+0.5);
366  hA[5+i]=h;
367  }
368 
369  sprintf(tt1,"my%02dStat",sectID);
370  hA[9]=new TH1F (tt1,"type of events ",30,.5,30.5);
371 
372  for(i=0;i<MaxSmdPlains;i++) {
373  sprintf(tt1,"fr%02d%c",sectID,i+'U');
374  sprintf(tt2,"freq. of MIP, UxV only, plane %02d%c; strip ID",sectID,i+'U');
375  h=new TH1F(tt1,tt2,MaxSmdStrips,-0.5,MaxSmdStrips-0.5);
376  hA[10+i]=h;
377 
378  sprintf(tt1,"mm%02d%c",sectID,i+'U');
379  sprintf(tt2,"freq of 00xx00 pattern per plane %02d%c",sectID,i+'U');
380  h=new TH1F(tt1,tt2,20,-0.5,19.5);
381  hA[12+i]=h;
382 
383  sprintf(tt1,"fr%02d%cm",sectID,i+'U');
384  sprintf(tt2,"freq. of best MIP, plane %02d%c; strip ID",sectID,i+'U');
385  h=new TH1F(tt1,tt2,MaxSmdStrips,-0.5,MaxSmdStrips-0.5);
386  hA[14+i]=h;
387 
388  sprintf(tt1,"ug%02d%c",sectID,i+'U');
389  sprintf(tt2,"used gains for plane %02d%c; strip ID; gain [ch/GeV]",sectID,i+'U');
390  h=new TH1F(tt1,tt2,MaxSmdStrips,0.5,MaxSmdStrips+0.5);
391  hA[16+i]=h;
392  }
393 
394  //..................
395  sprintf(tt1,"ep%02dUorV",sectID);
396  sprintf(tt2,"#Sigma E of 2-strips, normal angle, best MIP , SMD %02dUorV; strip ID; MIP #Sigma E (MeV) ",sectID);
397  h2=new TH2F(tt1,tt2,30,0,300,100,-.1,10.);
398  hA[20]=(TH1F*)h2;
399 
400  //..................
401  sprintf(tt1,"xy%02d",sectID);
402  sprintf(tt2,"MIP position , UxV only, sect=%02d; X(cm); Y(cm) ",sectID);
403  h2=new TH2F(tt1,tt2,500,-250,250,500,-250,250);
404  hA[21]=(TH1F*)h2;
405 
406  sprintf(tt1,"xy%02dm",sectID);
407  sprintf(tt2,"MIP position, best MIP, sect=%02d; X(cm); Y(cm) ",sectID);
408  h2=new TH2F(tt1,tt2,250,-250,250,250,-250,250);
409  hA[22]=(TH1F*)h2;
410 
411  //..................
412  sprintf(tt1,"eq%02dUV",sectID);
413  sprintf(tt2,"#Sigma from 4-strips, best MIP, plane %02dU+V; eta bin; #DeltaE MeV",sectID);
414  h2=new TH2F(tt1,tt2,12,.5,12.5,50,-.1,7.5);
415  hA[23]=(TH1F*)h2;
416 
417  //................
418  sprintf(tt1,"ca%02d",sectID);
419  sprintf(tt2,"# UxV candidates per tower, sector=%d; x=spiral=iPhi+60*iEta",sectID);
420  int mxTw=MaxEtaBins*MaxPhiBins;
421  hA[24]=new TH1F(tt1,tt2,mxTw,-0.5,mxTw-0.5);
422 
423 
424  // add histos to the list (if provided)
425  if(HList) {
426  for(i=0;i<32;i++) {
427  if(hA[i]==0) continue;
428  HList->Add(hA[i]);
429  }
430  }
431 
432 }
433 
434 
435 //--------------------------------------------------
436 //--------------------------------------------------
437 void EEsmdCal::fillOneTailHisto(char cut, int iEta, int iPhi){
438  // Fill inclusive spectra for one sector
439  int iCut=cut-'a';
440  assert(iCut>=0 && iCut<kCut);
441 
442  int iT=0;
443  for(iT=0;iT<mxTile;iT++) {
444  TH1F *h=hT[iCut][iT][iEta][iPhi];
445  h->Fill(tileAdc[iT][iEta][iPhi]);
446  }
447 }
448 
449 //--------------------------------------------------
450 //--------------------------------------------------
451 void EEsmdCal::fillSmdHisto_a(){
452  // Fill inclusive spectra for one sector
453  int iCut='a'-'a';
454 
455  int iuv,istrip;
456  for(iuv=0;iuv<MaxSmdPlains;iuv++) {
457  float *adc=smdAdc[iuv];
458  TH1F **hs=hSs[iCut][iuv]; // single strip spectra
459  // note, to loop over N-2 strips to get right the sum of pairs
460  for(istrip=0;istrip<MaxSmdStrips;istrip++) {
461  hs[istrip]->Fill(adc[istrip]);
462  }
463  }
464 }
465 
466 
467 //-------------------------------------------------
468 //-------------------------------------------------
469 void EEsmdCal::saveHisto(TString fname){
470  TString outName=fname+".hist.root";
471  TFile f( outName,"recreate");
472  assert(f.IsOpen());
473  printf("%d histos are written to '%s' ...\n",HList->GetEntries(),outName.Data());
474  HList->Write();
475  f.Close();
476 }
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
EEmcDbItem::tube
char tube[StEEmcNameLen]
name of PMT or MAPMT pixel
Definition: EEmcDbItem.h:21
EEsmdCal::tileAdc
float tileAdc[mxTile][MaxEtaBins][MaxPhiBins]
Definition: EEsmdCal.h:100
Generated by   doxygen 1.8.5