StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StEmcEqualMaker.cxx
1 #include "StEmcEqualMaker.h"
2 #include "TFile.h"
3 #include "TROOT.h"
4 #include "TF1.h"
5 #include "StEventTypes.h"
6 #include "StEvent.h"
7 #include "StEmcUtil/geometry/StEmcGeom.h"
8 #include "StDbLib/StDbManager.hh"
9 #include "StDbLib/StDbTable.h"
10 #include "StDbLib/StDbConfigNode.hh"
11 #include "tables/St_emcCalib_Table.h"
12 
13 ClassImp(StEmcEqualMaker)
14 
15 //_____________________________________________________________________________
16 StEmcEqualMaker::StEmcEqualMaker(const char *name):StEmcCalibMaker(name)
17 {
18  setRange(700);
19 }
20 //_____________________________________________________________________________
21 StEmcEqualMaker::~StEmcEqualMaker()
22 {
23 }
24 //_____________________________________________________________________________
25 Int_t StEmcEqualMaker::Init()
26 {
27  float pi = 3.1415926;
28  mA = new TH1F("mA","",getNChannel(),0.5,getNChannel()+0.5);
29  mB = new TH1F("mB","",getNChannel(),0.5,getNChannel()+0.5);
30  mADistr = new TH1F("mADistr","",200,0,10);
31  mBDistr = new TH1F("mBDistr","",200,-10,10);
32  mStatus = new TH1F("mStatus","",getNChannel(),0.5,getNChannel()+0.5);
33  mRefSlopes = new TH1F("mRefSlopes","",getNChannel(),0.5,getNChannel()+0.5);
34  mRefAmps = new TH1F("mRefAmps","",getNChannel(),0.5,getNChannel()+0.5);
35  mSlopes = new TH1F("mSlopes","",getNChannel(),0.5,getNChannel()+0.5);
36  mAmps = new TH1F("mAmps","",getNChannel(),0.5,getNChannel()+0.5);
37  mSlopesTheta = new TH2F("mSlopesTheta","",40,45*pi/180,135*pi/180,100,0,100);
38 
39  mSpecName ="mSpecEqual";
40  mAcceptName = "mAcceptEqual";
41  return StEmcCalibMaker::Init();
42 }
43 //_____________________________________________________________________________
44 void StEmcEqualMaker::Clear(Option_t *option)
45 {
46 }
47 //_____________________________________________________________________________
48 Int_t StEmcEqualMaker::Make()
49 {
50  if(!accept()) return kStOk;
51 
52  for(int i=0;i<mNChannel;i++)
53  {
54  int id = i+1;
55  float rms = getCalib()->getPedRms(mDetector,id);
56  float adc = (float) getCalib()->getADCPedSub(mDetector,id);
57  if(adc!=0 && adc>3*rms) fill(id,adc);
58  }
59 
60  return kStOK;
61 }
62 //_____________________________________________________________________________
63 Int_t StEmcEqualMaker::Finish()
64 {
65  saveHist((char*)mFileName.Data());
66  return StMaker::Finish();
67 }
68 //_____________________________________________________________________________
69 void StEmcEqualMaker::equalize(int mode,int DeltaEta,bool Et)
70 {
71  mADistr->Reset();
72  mBDistr->Reset();
73  StEmcGeom *g = getGeom();
74  int Neta = g->NEta();
75  int Nsub = g->NSub();
76  int etabin = 0;
77  float avg[18000];
78  float MIN = 20;
79  float MAX = 500;
80  TF1 *f= new TF1("slopes","22*TMath::Sin(x)");
81 
82  for(int eta0 = -Neta; eta0<Neta; eta0+=DeltaEta)
83  {
84  etabin++;
85  float sum = 0;
86  float n = 0;
87  int m1 = 1;
88  int m2 = 60;
89  if(eta0<0) { m1 = 61; m2 = 120; }
90  if(DeltaEta == 2*Neta){ m1=1; m2=120;}
91  for(int j = 0;j<DeltaEta;j++)
92  {
93  int eta = eta0+j;
94  if(eta>=0) eta+=1;
95  cout <<"etabin = "<<etabin<<" eta = "<<eta<<endl;
96  eta = abs(eta);
97  //cout <<"Calculating averages ..."<<endl;
98  for(int m = m1;m<=m2;m++)
99  {
100  for(int sub = 1;sub<=Nsub;sub++)
101  {
102  int id;
103  g->getId(m,eta,sub,id);
104  TH1D * h = getSpec(id);
105  if(h->Integral()>0)
106  {
107  float m1,r1;
108  getMeanAndRms(h,MIN,MAX,&m1,&r1);
109  sum+=m1;
110  n++;
111  avg[id-1] = m1;
112  }
113  }
114  }
115  }
116  float AVG = 0;
117  if(n>0) AVG=sum/n;
118  if(AVG>0)
119  {
120  int ID = 0;
121  float DCA = 99999;
122  //cout <<"Looking for reference ..."<<endl;
123  for(int m = m1;m<=m2;m++) for(int j = 0;j<DeltaEta;j++) for(int sub = 1;sub<=Nsub;sub++)
124  {
125  int id;
126  int eta = eta0+j;
127  if(eta>=0) eta+=1;
128  eta = abs(eta);
129  g->getId(m,eta,sub,id);
130  if(avg[id-1]>0 && fabs(avg[id-1]-AVG)<DCA) { DCA = fabs(avg[id-1]-AVG); ID = id;}
131  }
132  if(ID>0)
133  {
134  cout <<"Reference found id = "<<ID<<" avg = "<<avg[ID-1]<<endl;
135  cout <<"Equalizing eta bin ..."<<endl;
136  for(int m = m1;m<=m2;m++) for(int j = 0;j<DeltaEta;j++) for(int sub = 1;sub<=Nsub;sub++)
137  {
138  int eta = eta0+j;
139  if(eta>=0) eta+=1;
140  eta = abs(eta);
141  int id;
142  g->getId(m,eta,sub,id);
143  if(mode!=5) equalizeRelative(ID,id,mode,Et);
144  else equalizeToFunction(id,f);
145 
146  }
147  }
148  }
149  }
150  delete f;
151  return;
152 }
153 //_____________________________________________________________________________
154 void StEmcEqualMaker::equalizeRelative(int ref, int channel,int mode, bool Et)
155 {
156  // mean and RMS modes
157  // 0 - mean and RMS with liear average
158  // 1 - mean with linear average
159  // 2 - mean and rms with log average
160  // 3 - mean with log average
161  // 4 - exponential fit
162 
163  TH1D *h1 = getSpec(ref,"ref");
164  TH1D *h2 = getSpec(channel,"channel");
165 
166  float integral1 = h1->Integral();
167  float integral2 = h2->Integral();
168 
169  mA->SetBinContent(channel,0);
170  mB->SetBinContent(channel,0);
171  mStatus->SetBinContent(channel,0);
172  mRefSlopes->SetBinContent(channel,0);
173  mRefAmps->SetBinContent(channel,0);
174  mSlopes->SetBinContent(channel,0);
175  mAmps->SetBinContent(channel,0);
176 
177  float MIN = 20;
178  float MAX = 500;
179 
180  //cout <<"ref = "<<ref<<" channel = "<<channel<<endl;
181  //cout <<"h1 = "<<h1<<" h2 = "<<h2<<endl;
182  //cout <<"integral 1 = "<<integral1<<" integral 2 = "<<integral2<<endl;
183 
184  if((integral1==0 || integral2==0) && h1!=h2)
185  {
186  delete h1;
187  delete h2;
188  return;
189  }
190  if((integral1==0 || integral2==0) && h1==h2)
191  {
192  delete h1;
193  return;
194  }
195 
196  bool EqDone=kFALSE;
197  float a=0,b=0;
198 
199  if(mode==-1)
200  {
201  a = 1;
202  b = 0;
203  EqDone=true;
204  }
205 
206  if(mode>=0 && mode <=3)
207  {
208  float m1,r1,m2,r2;
209  if(mode==0 || mode==1)
210  {
211  getMeanAndRms(h1,MIN,MAX,&m1,&r1);
212  getMeanAndRms(h2,MIN,MAX,&m2,&r2);
213  }
214  if(mode==2 || mode==3)
215  {
216  getLogMeanAndRms(h1,MIN,MAX,&m1,&r1);
217  getLogMeanAndRms(h2,MIN,MAX,&m2,&r2);
218  }
219  if(mode==0 || mode==2)
220  {
221  a=r1/r2;
222  b=m1-a*m2;
223  }
224  if(mode==1 || mode==3)
225  {
226  a=m1/m2;
227  b=0;
228  }
229  EqDone=kTRUE;
230  cout <<" id = "<<channel<<" ref = "<<ref<<" mean = "<<m1<<" , "<<m2
231  <<" rms = "<<r1<<" , "<<r2
232  <<" a = "<<a<<" b = "<<b<<endl;
233  }
234 
235  float m1=0,m2=0,A1=0,A2=0;
236  if(mode==4)
237  {
238  TF1 *f=new TF1("ff","[0]*exp(-x/[1])",MIN,MAX);
239  float I1 = h1->Integral(h1->FindBin(MIN),h1->FindBin(MAX));
240  f->SetParameter(1,10);
241  h1->Fit(f,"RQN");
242  m1 = f->GetParameter(1);
243  A1 = f->GetParameter(0);
244  mRefSlopes->SetBinContent(channel,m1);
245  mRefAmps->SetBinContent(channel,A1);
246 
247  float I2 = h2->Integral(h2->FindBin(MIN),h2->FindBin(MAX));
248  h2->Fit(f,"RQN");
249  m2 = f->GetParameter(1);
250  A2 = f->GetParameter(0);
251  mSlopes->SetBinContent(channel,m2);
252  mAmps->SetBinContent(channel,A2);
253 
254  a=m1/m2;
255  b=-::log((A2*I1)/(A1*I2));
256  EqDone=true;
257  if(!finite(a) || !finite(b) || a<=0 || b>1000) {EqDone = false; a = 0;}
258  b=0;
259  cout <<" id = "<<channel<<" ref = "<<ref<<" slopes = "<<m2<<" , "<<m1
260  <<" a = "<<a<<" b = "<<b<<" EQDONE = "<<(Int_t)EqDone<<endl;
261  delete f;
262  }
263 
264  if (EqDone)
265  {
266  mA->SetBinContent(channel,a);
267  mB->SetBinContent(channel,b);
268  mADistr->Fill(a);
269  mBDistr->Fill(b);
270  mStatus->SetBinContent(channel,1);
271 
272  float theta;
273  int mod,eta,sub;
274  getGeom()->getBin(channel,mod,eta,sub);
275  getGeom()->getTheta(mod,eta,theta);
276  mSlopesTheta->Fill(theta,m2);
277 
278  }
279  else mStatus->SetBinContent(channel,2);
280 
281  delete h1;
282  delete h2;
283  return;
284 }
285 //_____________________________________________________________________________
286 void StEmcEqualMaker::equalizeToFunction(int channel,TF1 *func)
287 {
288  TH1D *h2 = getSpec(channel,"channel");
289 
290  float integral2 = h2->Integral();
291 
292  mA->SetBinContent(channel,0);
293  mB->SetBinContent(channel,0);
294  mStatus->SetBinContent(channel,0);
295  mRefSlopes->SetBinContent(channel,0);
296  mRefAmps->SetBinContent(channel,0);
297  mSlopes->SetBinContent(channel,0);
298  mAmps->SetBinContent(channel,0);
299 
300  float MIN = 20;
301  float MAX = 500;
302 
303  if(integral2==0)
304  {
305  delete h2;
306  mA->SetBinContent(channel,0);
307  mB->SetBinContent(channel,0);
308  return;
309  }
310 
311  bool EqDone=kFALSE;
312  float a=0,b=0;
313 
314  float theta;
315  int mod,eta,sub;
316  getGeom()->getBin(channel,mod,eta,sub);
317  getGeom()->getTheta(mod,eta,theta);
318 
319  float m1=0,m2=0,A1=0,A2=0;
320  TF1 *f=new TF1("ff","[0]*exp(-x/[1])",MIN,MAX);
321  f->SetParameter(1,10);
322 
323  m1 = func->Eval(theta);
324  mRefSlopes->SetBinContent(channel,m1);
325  mRefAmps->SetBinContent(channel,A1);
326 
327  if(channel<=2400)
328  {
329  //float I2 = h2->Integral(h2->FindBin(MIN),h2->FindBin(MAX));
330  h2->Fit(f,"RQN");
331  m2 = f->GetParameter(1);
332  A2 = f->GetParameter(0);
333  mSlopes->SetBinContent(channel,m2);
334  mAmps->SetBinContent(channel,A2);
335 
336 
337  a=m1/m2;
338  } //else a =1; //protection because crates in the east side are not timed
339 
340  EqDone=true;
341  if(!finite(a) || !finite(b) || a<=0.01 || b>1000 || a > 10) { EqDone = false; a = 0;}
342  b=0;
343  cout <<" id = "<<channel<<" slopes = "<<m2<<" , "<<m1
344  <<" a = "<<a<<" b = "<<b<<" EQDONE = "<<(Int_t)EqDone<<endl;
345  delete f;
346 
347  if (EqDone)
348  {
349  mA->SetBinContent(channel,a);
350  mB->SetBinContent(channel,b);
351  mADistr->Fill(a);
352  mBDistr->Fill(b);
353  mStatus->SetBinContent(channel,1);
354 
355  mSlopesTheta->Fill(theta,m2);
356 
357  }
358  else mStatus->SetBinContent(channel,2);
359 
360  delete h2;
361  return;
362 }
363 //_____________________________________________________________________________
364 void StEmcEqualMaker::calcSlopes()
365 {
366 
367  float MIN = 20;
368  float MAX = 500;
369  mSlopes->Reset();
370  mAmps->Reset();
371  mSlopesTheta->Reset();
372  TF1 *f=new TF1("ff","[0]*exp(-x/[1])",MIN,MAX);
373  StEmcGeom *g = getGeom();
374 
375  for(int channel = 1;channel<=getNChannel();channel++)
376  {
377  TH1D *h2 = getSpec(channel,"channel");
378  float I2 = h2->Integral(h2->FindBin(MIN),h2->FindBin(MAX));
379  if(I2>0)
380  {
381  float theta;
382  int mod,eta,sub;
383  g->getBin(channel,mod,eta,sub);
384  g->getTheta(mod,eta,theta);
385 
386  float m2=0,A2=0;
387 
388  f->SetParameter(1,10);
389  h2->Fit(f,"RQN");
390  m2 = f->GetParameter(1);
391  A2 = f->GetParameter(0);
392  mSlopes->SetBinContent(channel,m2);
393  mAmps->SetBinContent(channel,A2);
394 
395  mSlopesTheta->Fill(theta,m2);
396  cout <<"Slope for channel "<<channel<<" is "<<m2<<endl;
397  }
398  delete h2;
399  }
400  delete f;
401 
402  return;
403 }
404 //_____________________________________________________________________________
405 void StEmcEqualMaker::saveEqual(int date,int time)
406 {
407  char ts[100];
408  TString n[] = {"bemcEqual","bprsEqual","bsmdeEqual","bsmdpEqual"};
409  sprintf(ts,"%s.%08d.%06d.root",n[getDetector()-1].Data(),date,time);
410  TFile *f = new TFile(ts,"RECREATE");
411  if(getSpec()) getSpec()->Write();
412  if(mA) mA->Write();
413  if(mB) mB->Write();
414  if(mADistr) mADistr->Write();
415  if(mBDistr) mBDistr->Write();
416  if(mStatus) mStatus->Write();
417  if(mRefSlopes) mRefSlopes->Write();
418  if(mRefAmps) mRefAmps->Write();
419  if(mSlopes) mSlopes->Write();
420  if(mAmps) mAmps->Write();
421  if(mSlopesTheta) mSlopesTheta->Write();
422  f->Close();
423  delete f;
424  return;
425 }
426 //_____________________________________________________________________________
427 void StEmcEqualMaker::loadEqual(char* file)
428 {
429  TFile *f = new TFile(file);
430  mA->Reset();
431  mB->Reset();
432  mADistr->Reset();
433  mBDistr->Reset();
434  mStatus->Reset();
435  mRefSlopes->Reset();
436  mRefAmps->Reset();
437  mSlopes->Reset();
438  mAmps->Reset();
439  mSlopesTheta->Reset();
440  getSpec()->Reset();
441 
442  if(f->Get("mSpec;1")) getSpec()->Add((TH2F*)f->Get("mSpec;1"),1);
443  if(f->Get("mA;1")) mA->Add((TH1F*)f->Get("mA;1"),1);
444  if(f->Get("mB;1")) mB->Add((TH1F*)f->Get("mB;1"),1);
445  if(f->Get("mADistr;1")) mADistr->Add((TH1F*)f->Get("mADistr;1"),1);
446  if(f->Get("mBDistr;1")) mBDistr->Add((TH1F*)f->Get("mBDistr;1"),1);
447  if(f->Get("mStatus;1")) mStatus->Add((TH1F*)f->Get("mStatus;1"),1);
448  if(f->Get("mRefSlopes;1")) mRefSlopes->Add((TH1F*)f->Get("mRefSlopes;1"),1);
449  if(f->Get("mRefAmps;1")) mRefAmps->Add((TH1F*)f->Get("mRefAmps;1"),1);
450  if(f->Get("mSlopes;1")) mSlopes->Add((TH1F*)f->Get("mSlopes;1"),1);
451  if(f->Get("mAmps;1")) mAmps->Add((TH1F*)f->Get("mAmps;1"),1);
452  if(f->Get("mSlopesTheta;1")) mSlopesTheta->Add((TH2F*)f->Get("mSlopesTheta;1"),1);
453  f->Close();
454  delete f;
455  return;
456 }
457 //_____________________________________________________________________________
458 TH1F* StEmcEqualMaker::getEtaBinSpec(int eta0, int eta1,TH2F* SPEC)
459 {
460  TH2F *spec = NULL;
461  if(!SPEC) spec = getSpec(); else spec = SPEC;
462  if(!spec) return NULL;
463  int mi = 1;
464  int mf = 60;
465  if(eta0<0 && eta1<0)
466  {
467  mi = 61;
468  mf = 120;
469  }
470  if(eta0*eta1<0)
471  {
472  mi =1;
473  mf = 120;
474  }
475  int ns = getGeom()->NSub();
476  int e0 = abs(eta0);
477  int e1 = abs(eta1);
478  TH1F *hist = (TH1F*)spec->ProjectionY("etabin",1,1);
479  hist->Reset();
480 
481  for(int m=mi;m<=mf;m++)
482  {
483  for(int e=e0;e<=e1;e++)
484  {
485  for(int s=1;s<=ns;s++)
486  {
487  int id;
488  getGeom()->getId(m,e,s,id);
489  char name[30];
490  sprintf(name,"Tower %d",id);
491  TH1F* tmp = rebin(id,name,spec);
492  if(tmp)
493  {
494  hist->Add(tmp,1);
495  delete tmp; tmp = NULL;
496  }
497  }
498  }
499  }
500  return hist;
501 }
502 //_____________________________________________________________________________
503 TH1F* StEmcEqualMaker::rebin(int id,const char *name, TH2F* SPEC)
504 {
505  TH2F *spec = NULL;
506  if(!SPEC) spec = getSpec(); else spec = SPEC;
507  if(!spec) return NULL;
508 
509  float seg = 50;
510  float a = mA->GetBinContent(id);
511  float b = mB->GetBinContent(id);
512  //cout <<"A = "<<a<<" B = "<<b<<endl;
513 
514  if(a<=0) return NULL;
515 
516  TH1F *h = (TH1F*)spec->ProjectionY(name,id,id);
517  if(a==1 && b==0) return h;
518 
519  float deltaBin=h->GetBinWidth(1)/seg;
520  int nbins=h->GetNbinsX();
521 
522  float Y[4096];
523  for(int i=0;i<4096;i++) Y[i] = 0;
524  for(int i=1;i<=nbins;i++) Y[i] = h->GetBinContent(i);
525 
526  h->Reset();
527 
528  for(int i=1;i<=nbins;i++)
529  {
530  float x = h->GetBinLowEdge(i);
531  float y = Y[i];
532  for(int j=0;j<(int)seg;j++)
533  {
534  float x1 = x+((float)j+0.5)*deltaBin;
535  float x2 = x1*a+b;
536  //if (a>1) cout <<position<<" "<<a<<" "<<x1<<" "<<x2<<endl;
537  h->Fill(x2,y/seg);
538  }
539  }
540  return h;
541 }
542 //_____________________________________________________________________________
543 // this method leaks memory because it draws a lot of histograms
544 // use it with care
545 // AASPSUAIDE
546 void StEmcEqualMaker::drawEtaBin(int eta0, int eta1,TH2F* SPEC)
547 {
548  TH2F *spec = NULL;
549  if(!SPEC) spec = getSpec(); else spec = SPEC;
550  if(!spec) return;
551  int mi = 1;
552  int mf = 60;
553  if(eta0<0 && eta1<0)
554  {
555  mi = 61;
556  mf = 120;
557  }
558  if(eta0*eta1<0)
559  {
560  mi =1;
561  mf = 120;
562  }
563  int ns = getGeom()->NSub();
564  int e0 = abs(eta0);
565  int e1 = abs(eta1);
566 
567  bool first = true;
568 
569  for(int m=mi;m<=mf;m++)
570  {
571  for(int e=e0;e<=e1;e++)
572  {
573  for(int s=1;s<=ns;s++)
574  {
575  int id;
576  getGeom()->getId(m,e,s,id);
577  char name[30];
578  sprintf(name,"Tower %d",id);
579  TH1F* tmp = rebin(id,name,spec);
580  if(tmp)
581  {
582  if(first) { tmp->Draw(); first = false;} else tmp->Draw("same");
583  }
584  }
585  }
586  }
587  return;
588 }
589 
590 
591 
592 
593 
StEmcEqualMaker::Make
virtual Int_t Make()
Definition: StEmcEqualMaker.cxx:48
kStOK
Definition: Stypes.h:40
StEmcEqualMaker::Clear
virtual void Clear(Option_t *option="")
User defined functions.
Definition: StEmcEqualMaker.cxx:44
StMaker::Finish
virtual Int_t Finish()
Definition: StMaker.cxx:776
StEmcGeom::getBin
Int_t getBin(const Float_t phi, const Float_t eta, Int_t &m, Int_t &e, Int_t &s) const
Definition: StEmcGeom.h:321
kStOk
Definition: Stypes.h:41
StEmcEqualMaker::Finish
virtual Int_t Finish()
Definition: StEmcEqualMaker.cxx:63
Generated by   doxygen 1.8.5