StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StEmcMipMaker.cxx
1 #include "StEmcMipMaker.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 #include "StEmcEqualMaker.h"
13 #include "TCanvas.h"
14 
15 ClassImp(StEmcMipMaker)
16 
17 //_____________________________________________________________________________
18 StEmcMipMaker::StEmcMipMaker(const char *name):StEmcCalibMaker(name)
19 {
20  mNPoints = 10;
21  mPmin = 1.0;
22  setRange(200);
23  for(int i=0;i<10;i++) mNFailed[i]=0;
24  mMipPos = 0;
25  mMipPosErr = 0;
26  mMipWid = 0;
27  mMipWidErr = 0;
28  mGain = 0;
29  mGainDistr = 0;
30  mChi2 = 0;
31  mIntegral = 0;
32 
33  funcFit = 0;
34  funcFitPeak = 0;
35  funcFitBack = 0;
36 
37 }
38 //_____________________________________________________________________________
39 StEmcMipMaker::~StEmcMipMaker()
40 {
41  delete mMipPos;
42  delete mMipPosErr;
43  delete mMipWid;
44  delete mMipWidErr;
45  delete mGain;
46  delete mGainDistr;
47  delete mChi2;
48  delete mIntegral;
49 
50  delete funcFit;
51  delete funcFitPeak;
52  delete funcFitBack;
53 
54 }
55 //_____________________________________________________________________________
56 Int_t StEmcMipMaker::Init()
57 {
58  mMipPos = new TH1F("mMipPos","",getNChannel(),0.5,getNChannel()+0.5);
59  mMipPosErr = new TH1F("mMipPosErr","",getNChannel(),0.5,getNChannel()+0.5);
60  mMipWid = new TH1F("mMipWid","",getNChannel(),0.5,getNChannel()+0.5);
61  mMipWidErr = new TH1F("mMipWidErr","",getNChannel(),0.5,getNChannel()+0.5);
62  mGain = new TH1F("mGain","",getNChannel(),0.5,getNChannel()+0.5);
63  mGainDistr = new TH1F("mGainDistr","",100,0,0.1);
64  mChi2 = new TH1F("mChi2","",getNChannel(),0.5,getNChannel()+0.5);
65  mIntegral = new TH1F("mIntegral","",getNChannel(),0.5,getNChannel()+0.5);
66 
67  mSpecName="mSpecMip";
68  mAcceptName = "mAcceptMip";
69 
70  return StEmcCalibMaker::Init();
71 }
72 //_____________________________________________________________________________
73 void StEmcMipMaker::Clear(Option_t *option)
74 {
75 }
76 //_____________________________________________________________________________
77 Int_t StEmcMipMaker::Make()
78 {
79  if(!accept()) return kStOk;
80 
81  cout <<"Number of tracks = "<<getCalib()->getNTracks()<<endl;
82  for(int i=0;i<getCalib()->getNTracks();i++)
83  {
84  int id = getCalib()->getTrackTower(i);
85  int id2 = getCalib()->getTrackTowerExit(i);
86  float p = getCalib()->getTrackP(i);
87  //int np = mCalib->getTrackNPoints(i);
88  float adc = (float) getCalib()->getADCPedSub(1,id);
89  //cout <<"candidate. p = "<<p<<" Tower id = "<<id<<" ADC = "<<adc<<endl;
90  if(p>mPmin && id!=0 && id2==id /*&& mCalib->isIsolatedTower(id)*/)
91  {
92  int nt = getCalib()->getNTracksInTower(id);
93  float rms = getCalib()->getPedRms(1,id);
94  if(adc!=0 && adc>2*rms && nt==1)
95  {
96  cout <<"Found MIP. p = "<<p<<" Tower id = "<<id<<" ADC = "<<adc<<endl;
97  fill(id,adc);
98  }
99  }
100  }
101 
102  return kStOK;
103 }
104 //_____________________________________________________________________________
105 Int_t StEmcMipMaker::Finish()
106 {
107  saveHist((char*)mFileName.Data());
108  return StMaker::Finish();
109 }
110 //_____________________________________________________________________________
111 void StEmcMipMaker::fit(TH1F* h)
112 {
113  float WIDTH = 10;
114  float max = h->GetBinCenter(h->GetMaximumBin());
115  max = 40;
116  if(max<15) WIDTH = 3;
117  if(max<5) max = 15;
118  float maxy = h->GetBinContent(h->GetMaximumBin());
119  float maxback = h->GetBinContent((int)(h->GetMaximumBin()+6*WIDTH));
120 
121  int NPars = 6;
122  float xmin = 7;
123  float xmax = (max*3>120) ? max*3 : 120;
124 
125  float p[] = {maxy-maxback,max,WIDTH,maxback,0,50};
126 
127  delete funcFitPeak;
128  funcFitPeak = new TF1("peak","gaus(0)",0,1000);
129  funcFitPeak->SetLineColor(2);
130  funcFitPeak->SetLineWidth(1);
131  delete funcFitBack;
132  funcFitBack = new TF1("back","gaus(0)",0,1000);
133  funcFitBack->SetLineWidth(1);
134  funcFitBack->SetLineColor(3);
135  delete funcFit;
136  funcFit = new TF1("func","gaus(0)+gaus(3)",0,1000);
137  funcFit->SetLineWidth(1);
138  funcFit->SetLineColor(4);
139 
140  TF1 *func = funcFit;
141  func->SetRange(xmin,xmax);
142 
143  for(int i=0;i<NPars;i++) func->SetParameter(i,p[i]);
144  func->SetParLimits(0,0,10000);
145  func->SetParLimits(3,0,10000);
146  func->SetParLimits(5,30,100);
147 
148  h->Fit(func,"RNQ");
149 
150  int nb = h->FindBin(xmax);
151  for(int i=h->FindBin(xmin);i<nb;i++)
152  {
153  //float x = h->GetBinCenter(i);
154  float y = h->GetBinContent(i);
155  float stat = sqrt(y);
156  //float binwidth = h->GetBinWidth(i);
157  //float deriv = (float) func->Derivative(x);
158  float sigma = sqrt(stat*stat);//+(deriv*binwidth/2)*(deriv*binwidth/2));
159  h->SetBinError(i,sigma);
160  }
161 
162  h->Fit(func,"RNQ");
163 
164  return;
165 }
166 //_____________________________________________________________________________
167 TH1F* StEmcMipMaker::findMip(int eta0,int eta1,StEmcEqualMaker* equal)
168 {
169  if(!equal) return NULL;
170  TH1F* h = equal->getEtaBinSpec(eta0,eta1,getSpec());
171  char name[60];
172  sprintf(name,"Mip.%d.%d",eta0,eta1);
173  h->SetName(name);
174  sprintf(name,"Mip for %d <= eta <= %d",eta0,eta1);
175  h->SetTitle(name);
176  if(h) fit(h);
177 
178  TF1* func = funcFit;
179 
180  int NPars = 6;
181 
182  for(int i=0;i<NPars;i++)
183  {
184  float a = func->GetParameter(i);
185  if(i<3) funcFitPeak->SetParameter(i,a); else funcFitBack->SetParameter(i-3,a);
186  }
187 
188  float chi=func->GetChisquare()/(func->GetNDF());
189  if(chi<0 || chi>1000) chi = 1000;
190  float peak = func->GetParameter(1);
191  if(peak<5) peak = 0 ;
192  float epeak = func->GetParError(1);
193  float w = func->GetParameter(2);
194  float ew = func->GetParError(2);
195  //float XMAX = func->GetMaximumX(7,200);
196  //if(fabs(XMAX-peak)>2) chi = 0;
197 
198  int m1 = 1;
199  int m2 = 60;
200  if(eta0<0 && eta1<0) {m1 = 61; m2 = 120;}
201  if(eta0*eta1<0) {m1 = 1; m2 = 120;}
202  eta0 = abs(eta0);
203  eta1 = abs(eta1);
204 
205  StEmcGeom *geo = getGeom();
206  int ns = geo->NSub();
207 
208  for(int m = m1;m<=m2;m++)
209  for(int e = eta0; e<=eta1; e++)
210  for(int s = 1; s<=ns; s++)
211  {
212  int id;
213  geo->getId(m,e,s,id);
214  float a = equal->getA()->GetBinContent(id);
215  float b = equal->getB()->GetBinContent(id);
216 
217  if(a>10) {a = 0; b = 0;}
218  if(a<0.1) {a = 0; b = 0;}
219 
220  mMipPos->SetBinContent(id,peak*a+b);
221  mMipPosErr->SetBinContent(id,epeak*a);
222  mMipWid->SetBinContent(id,w*a);
223  mMipWidErr->SetBinContent(id,ew*a);
224  mChi2->SetBinContent(id,chi);
225  mIntegral->SetBinContent(id,0);
226  cout <<"MIP peak = "<<peak<<"+-"<<epeak<<" width = "<<w<<"+-"<<ew<<endl;
227  cout <<"Final chi2 = "<<chi<<endl;
228  findGain(id,true);
229  }
230  return h;
231 }
232 //_____________________________________________________________________________
233 TH1F* StEmcMipMaker::findMip(int id, int rebin, bool print)
234 {
235  float mip = 0;
236  float width =0;
237  float chi = 0;
238  delete funcFitPeak; funcFitPeak = NULL;
239  delete funcFitBack; funcFitBack = NULL;
240  delete funcFit ; funcFit = NULL;
241  if(!mSpec) return NULL;
242  if(gROOT->FindObject("id")) delete (TH1D*)gROOT->FindObject("id");
243  TH1F *h = (TH1F*)mSpec->ProjectionY("id",id,id);
244  if(!h) { mNFailed[0]++; return NULL;}
245  h->Rebin(rebin);
246 
247  float integral= h->Integral();
248  if(integral==0) { mNFailed[0]++; mNFailed[1]++; return h;}
249 
250  if(print) cout <<"Fitting MIP for id = "<<id<<endl;
251  if(print) cout <<"Integral = "<<integral<<endl;
252 
253  fit(h);
254 
255  TF1 *func = funcFit;
256 
257  chi=func->GetChisquare()/(func->GetNDF());
258  if(chi<0 || chi>1000) chi = 1000;
259 
260  int NPars = 6;
261 
262  for(int i=0;i<NPars;i++)
263  {
264  float a = func->GetParameter(i);
265  if(i<3) funcFitPeak->SetParameter(i,a); else funcFitBack->SetParameter(i-3,a);
266  }
267 
268  float peak = func->GetParameter(1);
269  if(peak<5) peak = 0 ;
270  float epeak = func->GetParError(1);
271  float w = func->GetParameter(2);
272  float ew = func->GetParError(2);
273  float XMAX = func->GetMaximumX(7,200);
274  if(fabs(XMAX-peak)>2)
275  {
276  mNFailed[0]++;
277  mNFailed[3]++;
278  if(print) cout << "****************** FAILED fabs(XMAX-peak) *****************\n";
279  if(print) cout << "fabs(XMAX-peak) = "<<fabs(XMAX-peak)<<endl;
280  chi = 0;
281  }
282  mMipPos->Fill(id,peak);
283  mMipPosErr->Fill(id,epeak);
284  mMipWid->Fill(id,w);
285  mMipWidErr->Fill(id,ew);
286  mChi2->Fill(id,chi);
287  mIntegral->Fill(id,integral);
288  if(print) cout <<"MIP peak = "<<peak<<"+-"<<epeak<<" width = "<<w<<"+-"<<ew<<endl;
289  if(print) cout <<"Final chi2 = "<<chi<<endl;
290  mip = peak;
291  width = w;
292  return h;
293 }
294 //_____________________________________________________________________________
295 float StEmcMipMaker::findGain(int id,bool print)
296 {
297  float mip = getMipPosition(id);
298  if(mip==0) return 0;
299  StEmcGeom *geom = StEmcGeom::instance("bemc");
300  float eta,phi;
301  geom->getEtaPhi(id,eta,phi);
302  float theta=2.*atan(exp(-eta));
303  float EoverMIP = 0.261;
304  float MipE = EoverMIP*(1.+0.056*eta*eta)/sin(theta);
305  float gain = MipE/mip;
306  float chi = getChi2(id);
307  if(chi<0.001 || chi>30) gain = 0;
308  if(gain <0 || gain > 10000) gain = 0;
309  if(gain==0) mNFailed[4]++;
310  mGain->Fill(id,gain);
311  mGainDistr->Fill(gain);
312  if(print) cout <<"gain for id "<<id<<" = "<<gain<<endl;
313  if(print) cout <<"-------------------------------------------------------\n";
314  return gain;
315 }
316 //_____________________________________________________________________________
317 void StEmcMipMaker::mipCalib()
318 {
319  for(int id = 0;id<4800;id++)
320  {
321  if((id-1)%100==0) cout <<"Doing MIP calibration for id "<<id<<endl;
322  delete findMip(id,3,false);
323  delete funcFit;
324  delete funcFitPeak;
325  delete funcFitBack;
326  findGain(id,false);
327  }
328 }
329 //_____________________________________________________________________________
330 void StEmcMipMaker::mipCalib(int eta0, int eta1, int etabin, StEmcEqualMaker* equal, bool draw)
331 {
332  TCanvas *c1 = NULL;
333  //int ne = (eta1-eta0+1)/etabin;
334  int pad = 1;
335  for(int e = eta0; e<=eta1;e++)
336  {
337  TH1F *h = findMip(e,e+etabin-1,equal);
338  if(draw)
339  {
340  if(pad==1) {c1 = new TCanvas(); c1->Divide(2,2);}
341  c1->cd(pad);
342  h->Draw();
343  funcFitBack->Draw("same");
344  funcFitPeak->Draw("same");
345  funcFit->Draw("same");
346  pad++;
347  if(pad>4) pad=1;
348  }
349  else
350  {
351  delete h;
352  delete funcFit;
353  delete funcFitPeak;
354  delete funcFitBack;
355  }
356  }
357 }
358 //_____________________________________________________________________________
359 TH1F* StEmcMipMaker::compareToDb(char* timeStamp,int mode)
360 {
361 
362  StDbManager* mgr=StDbManager::Instance();
363  StDbConfigNode* node=mgr->initConfig(dbCalibrations,dbEmc);
364  StDbTable* table = node->addDbTable("bemcCalib");
365  mgr->setRequestTime(timeStamp);
366  mgr->fetchDbTable(table);
367  emcCalib_st *calib = (emcCalib_st*) table->GetTable();
368  if(!calib) return NULL;
369 
370  char line[100];
371  sprintf(line,"Comparison with DB - %s",timeStamp);
372  TH1F *h = NULL;
373  if(mode == 0) h = new TH1F("h",line,4800,0.5,4800.5);
374  else h = new TH1F("h",line,100,0,4);
375  for(int i=0;i<4800;i++)
376  {
377  int id = i+1;
378  float gain = calib[0].AdcToE[i][1];
379  float gain2 = getGain(id);
380  float ratio = 0;
381  if(gain!=0) ratio = gain2/gain;
382  if(mode==0) h->Fill(id,ratio); else h->Fill(ratio);
383  }
384  return h;
385 }
386 //_____________________________________________________________________________
387 void StEmcMipMaker::saveToDb(char* timeStamp)
388 {
389  emcCalib_st tnew;
390  for(int i=0;i<4800;i++)
391  {
392  int id = i+1;
393  float gain = getGain(id);
394  if(gain==0) tnew.Status[i] = 0; else tnew.Status[i] = 1;
395  tnew.AdcToE[i][0] = 0;
396  tnew.AdcToE[i][1] = gain;
397  tnew.AdcToE[i][2] = 0;
398  tnew.AdcToE[i][3] = 0;
399  tnew.AdcToE[i][4] = 0;
400  }
401  StDbManager* mgr=StDbManager::Instance();
402  StDbConfigNode* node=mgr->initConfig(dbCalibrations,dbEmc);
403  StDbTable* table = node->addDbTable("bemcCalib");
404  table->SetTable((char*)&tnew,1);
405  mgr->setStoreTime(timeStamp);
406  mgr->storeDbTable(table);
407  return;
408 }
409 
410 
411 
412 
413 
StEmcMipMaker::Clear
virtual void Clear(Option_t *option="")
User defined functions.
Definition: StEmcMipMaker.cxx:73
StDbTable::SetTable
virtual void SetTable(char *data, int nrows, int *idList=0)
calloc&#39;d version of data for StRoot
Definition: StDbTable.cc:550
kStOK
Definition: Stypes.h:40
StEmcMipMaker::Make
virtual Int_t Make()
Definition: StEmcMipMaker.cxx:77
StDbManager::Instance
static StDbManager * Instance()
strdup(..) is not ANSI
Definition: StDbManager.cc:155
StMaker::Finish
virtual Int_t Finish()
Definition: StMaker.cxx:776
StEmcMipMaker::Finish
virtual Int_t Finish()
Definition: StEmcMipMaker.cxx:105
kStOk
Definition: Stypes.h:41
Generated by   doxygen 1.8.5