00001 #include "StEmcMipMaker.h"
00002 #include "TFile.h"
00003 #include "TROOT.h"
00004 #include "TF1.h"
00005 #include "StEventTypes.h"
00006 #include "StEvent.h"
00007 #include "StEmcUtil/geometry/StEmcGeom.h"
00008 #include "StDbLib/StDbManager.hh"
00009 #include "StDbLib/StDbTable.h"
00010 #include "StDbLib/StDbConfigNode.hh"
00011 #include "tables/St_emcCalib_Table.h"
00012 #include "StEmcEqualMaker.h"
00013 #include "TCanvas.h"
00014
00015 ClassImp(StEmcMipMaker)
00016
00017
00018 StEmcMipMaker::StEmcMipMaker(const char *name):StEmcCalibMaker(name)
00019 {
00020 mNPoints = 10;
00021 mPmin = 1.0;
00022 setRange(200);
00023 for(int i=0;i<10;i++) mNFailed[i]=0;
00024 mMipPos = 0;
00025 mMipPosErr = 0;
00026 mMipWid = 0;
00027 mMipWidErr = 0;
00028 mGain = 0;
00029 mGainDistr = 0;
00030 mChi2 = 0;
00031 mIntegral = 0;
00032
00033 funcFit = 0;
00034 funcFitPeak = 0;
00035 funcFitBack = 0;
00036
00037 }
00038
00039 StEmcMipMaker::~StEmcMipMaker()
00040 {
00041 delete mMipPos;
00042 delete mMipPosErr;
00043 delete mMipWid;
00044 delete mMipWidErr;
00045 delete mGain;
00046 delete mGainDistr;
00047 delete mChi2;
00048 delete mIntegral;
00049
00050 delete funcFit;
00051 delete funcFitPeak;
00052 delete funcFitBack;
00053
00054 }
00055
00056 Int_t StEmcMipMaker::Init()
00057 {
00058 mMipPos = new TH1F("mMipPos","",getNChannel(),0.5,getNChannel()+0.5);
00059 mMipPosErr = new TH1F("mMipPosErr","",getNChannel(),0.5,getNChannel()+0.5);
00060 mMipWid = new TH1F("mMipWid","",getNChannel(),0.5,getNChannel()+0.5);
00061 mMipWidErr = new TH1F("mMipWidErr","",getNChannel(),0.5,getNChannel()+0.5);
00062 mGain = new TH1F("mGain","",getNChannel(),0.5,getNChannel()+0.5);
00063 mGainDistr = new TH1F("mGainDistr","",100,0,0.1);
00064 mChi2 = new TH1F("mChi2","",getNChannel(),0.5,getNChannel()+0.5);
00065 mIntegral = new TH1F("mIntegral","",getNChannel(),0.5,getNChannel()+0.5);
00066
00067 mSpecName="mSpecMip";
00068 mAcceptName = "mAcceptMip";
00069
00070 return StEmcCalibMaker::Init();
00071 }
00072
00073 void StEmcMipMaker::Clear(Option_t *option)
00074 {
00075 }
00076
00077 Int_t StEmcMipMaker::Make()
00078 {
00079 if(!accept()) return kStOk;
00080
00081 cout <<"Number of tracks = "<<getCalib()->getNTracks()<<endl;
00082 for(int i=0;i<getCalib()->getNTracks();i++)
00083 {
00084 int id = getCalib()->getTrackTower(i);
00085 int id2 = getCalib()->getTrackTowerExit(i);
00086 float p = getCalib()->getTrackP(i);
00087
00088 float adc = (float) getCalib()->getADCPedSub(1,id);
00089
00090 if(p>mPmin && id!=0 && id2==id )
00091 {
00092 int nt = getCalib()->getNTracksInTower(id);
00093 float rms = getCalib()->getPedRms(1,id);
00094 if(adc!=0 && adc>2*rms && nt==1)
00095 {
00096 cout <<"Found MIP. p = "<<p<<" Tower id = "<<id<<" ADC = "<<adc<<endl;
00097 fill(id,adc);
00098 }
00099 }
00100 }
00101
00102 return kStOK;
00103 }
00104
00105 Int_t StEmcMipMaker::Finish()
00106 {
00107 saveHist((char*)mFileName.Data());
00108 return StMaker::Finish();
00109 }
00110
00111 void StEmcMipMaker::fit(TH1F* h)
00112 {
00113 float WIDTH = 10;
00114 float max = h->GetBinCenter(h->GetMaximumBin());
00115 max = 40;
00116 if(max<15) WIDTH = 3;
00117 if(max<5) max = 15;
00118 float maxy = h->GetBinContent(h->GetMaximumBin());
00119 float maxback = h->GetBinContent((int)(h->GetMaximumBin()+6*WIDTH));
00120
00121 int NPars = 6;
00122 float xmin = 7;
00123 float xmax = (max*3>120) ? max*3 : 120;
00124
00125 float p[] = {maxy-maxback,max,WIDTH,maxback,0,50};
00126
00127 delete funcFitPeak;
00128 funcFitPeak = new TF1("peak","gaus(0)",0,1000);
00129 funcFitPeak->SetLineColor(2);
00130 funcFitPeak->SetLineWidth(1);
00131 delete funcFitBack;
00132 funcFitBack = new TF1("back","gaus(0)",0,1000);
00133 funcFitBack->SetLineWidth(1);
00134 funcFitBack->SetLineColor(3);
00135 delete funcFit;
00136 funcFit = new TF1("func","gaus(0)+gaus(3)",0,1000);
00137 funcFit->SetLineWidth(1);
00138 funcFit->SetLineColor(4);
00139
00140 TF1 *func = funcFit;
00141 func->SetRange(xmin,xmax);
00142
00143 for(int i=0;i<NPars;i++) func->SetParameter(i,p[i]);
00144 func->SetParLimits(0,0,10000);
00145 func->SetParLimits(3,0,10000);
00146 func->SetParLimits(5,30,100);
00147
00148 h->Fit(func,"RNQ");
00149
00150 int nb = h->FindBin(xmax);
00151 for(int i=h->FindBin(xmin);i<nb;i++)
00152 {
00153
00154 float y = h->GetBinContent(i);
00155 float stat = sqrt(y);
00156
00157
00158 float sigma = sqrt(stat*stat);
00159 h->SetBinError(i,sigma);
00160 }
00161
00162 h->Fit(func,"RNQ");
00163
00164 return;
00165 }
00166
00167 TH1F* StEmcMipMaker::findMip(int eta0,int eta1,StEmcEqualMaker* equal)
00168 {
00169 if(!equal) return NULL;
00170 TH1F* h = equal->getEtaBinSpec(eta0,eta1,getSpec());
00171 char name[60];
00172 sprintf(name,"Mip.%d.%d",eta0,eta1);
00173 h->SetName(name);
00174 sprintf(name,"Mip for %d <= eta <= %d",eta0,eta1);
00175 h->SetTitle(name);
00176 if(h) fit(h);
00177
00178 TF1* func = funcFit;
00179
00180 int NPars = 6;
00181
00182 for(int i=0;i<NPars;i++)
00183 {
00184 float a = func->GetParameter(i);
00185 if(i<3) funcFitPeak->SetParameter(i,a); else funcFitBack->SetParameter(i-3,a);
00186 }
00187
00188 float chi=func->GetChisquare()/(func->GetNDF());
00189 if(chi<0 || chi>1000) chi = 1000;
00190 float peak = func->GetParameter(1);
00191 if(peak<5) peak = 0 ;
00192 float epeak = func->GetParError(1);
00193 float w = func->GetParameter(2);
00194 float ew = func->GetParError(2);
00195
00196
00197
00198 int m1 = 1;
00199 int m2 = 60;
00200 if(eta0<0 && eta1<0) {m1 = 61; m2 = 120;}
00201 if(eta0*eta1<0) {m1 = 1; m2 = 120;}
00202 eta0 = abs(eta0);
00203 eta1 = abs(eta1);
00204
00205 StEmcGeom *geo = getGeom();
00206 int ns = geo->NSub();
00207
00208 for(int m = m1;m<=m2;m++)
00209 for(int e = eta0; e<=eta1; e++)
00210 for(int s = 1; s<=ns; s++)
00211 {
00212 int id;
00213 geo->getId(m,e,s,id);
00214 float a = equal->getA()->GetBinContent(id);
00215 float b = equal->getB()->GetBinContent(id);
00216
00217 if(a>10) {a = 0; b = 0;}
00218 if(a<0.1) {a = 0; b = 0;}
00219
00220 mMipPos->SetBinContent(id,peak*a+b);
00221 mMipPosErr->SetBinContent(id,epeak*a);
00222 mMipWid->SetBinContent(id,w*a);
00223 mMipWidErr->SetBinContent(id,ew*a);
00224 mChi2->SetBinContent(id,chi);
00225 mIntegral->SetBinContent(id,0);
00226 cout <<"MIP peak = "<<peak<<"+-"<<epeak<<" width = "<<w<<"+-"<<ew<<endl;
00227 cout <<"Final chi2 = "<<chi<<endl;
00228 findGain(id,true);
00229 }
00230 return h;
00231 }
00232
00233 TH1F* StEmcMipMaker::findMip(int id, int rebin, bool print)
00234 {
00235 float mip = 0;
00236 float width =0;
00237 float chi = 0;
00238 delete funcFitPeak; funcFitPeak = NULL;
00239 delete funcFitBack; funcFitBack = NULL;
00240 delete funcFit ; funcFit = NULL;
00241 if(!mSpec) return NULL;
00242 if(gROOT->FindObject("id")) delete (TH1D*)gROOT->FindObject("id");
00243 TH1F *h = (TH1F*)mSpec->ProjectionY("id",id,id);
00244 if(!h) { mNFailed[0]++; return NULL;}
00245 h->Rebin(rebin);
00246
00247 float integral= h->Integral();
00248 if(integral==0) { mNFailed[0]++; mNFailed[1]++; return h;}
00249
00250 if(print) cout <<"Fitting MIP for id = "<<id<<endl;
00251 if(print) cout <<"Integral = "<<integral<<endl;
00252
00253 fit(h);
00254
00255 TF1 *func = funcFit;
00256
00257 chi=func->GetChisquare()/(func->GetNDF());
00258 if(chi<0 || chi>1000) chi = 1000;
00259
00260 int NPars = 6;
00261
00262 for(int i=0;i<NPars;i++)
00263 {
00264 float a = func->GetParameter(i);
00265 if(i<3) funcFitPeak->SetParameter(i,a); else funcFitBack->SetParameter(i-3,a);
00266 }
00267
00268 float peak = func->GetParameter(1);
00269 if(peak<5) peak = 0 ;
00270 float epeak = func->GetParError(1);
00271 float w = func->GetParameter(2);
00272 float ew = func->GetParError(2);
00273 float XMAX = func->GetMaximumX(7,200);
00274 if(fabs(XMAX-peak)>2)
00275 {
00276 mNFailed[0]++;
00277 mNFailed[3]++;
00278 if(print) cout << "****************** FAILED fabs(XMAX-peak) *****************\n";
00279 if(print) cout << "fabs(XMAX-peak) = "<<fabs(XMAX-peak)<<endl;
00280 chi = 0;
00281 }
00282 mMipPos->Fill(id,peak);
00283 mMipPosErr->Fill(id,epeak);
00284 mMipWid->Fill(id,w);
00285 mMipWidErr->Fill(id,ew);
00286 mChi2->Fill(id,chi);
00287 mIntegral->Fill(id,integral);
00288 if(print) cout <<"MIP peak = "<<peak<<"+-"<<epeak<<" width = "<<w<<"+-"<<ew<<endl;
00289 if(print) cout <<"Final chi2 = "<<chi<<endl;
00290 mip = peak;
00291 width = w;
00292 return h;
00293 }
00294
00295 float StEmcMipMaker::findGain(int id,bool print)
00296 {
00297 float mip = getMipPosition(id);
00298 if(mip==0) return 0;
00299 StEmcGeom *geom = StEmcGeom::instance("bemc");
00300 float eta,phi;
00301 geom->getEtaPhi(id,eta,phi);
00302 float theta=2.*atan(exp(-eta));
00303 float EoverMIP = 0.261;
00304 float MipE = EoverMIP*(1.+0.056*eta*eta)/sin(theta);
00305 float gain = MipE/mip;
00306 float chi = getChi2(id);
00307 if(chi<0.001 || chi>30) gain = 0;
00308 if(gain <0 || gain > 10000) gain = 0;
00309 if(gain==0) mNFailed[4]++;
00310 mGain->Fill(id,gain);
00311 mGainDistr->Fill(gain);
00312 if(print) cout <<"gain for id "<<id<<" = "<<gain<<endl;
00313 if(print) cout <<"-------------------------------------------------------\n";
00314 return gain;
00315 }
00316
00317 void StEmcMipMaker::mipCalib()
00318 {
00319 for(int id = 0;id<4800;id++)
00320 {
00321 if((id-1)%100==0) cout <<"Doing MIP calibration for id "<<id<<endl;
00322 delete findMip(id,3,false);
00323 delete funcFit;
00324 delete funcFitPeak;
00325 delete funcFitBack;
00326 findGain(id,false);
00327 }
00328 }
00329
00330 void StEmcMipMaker::mipCalib(int eta0, int eta1, int etabin, StEmcEqualMaker* equal, bool draw)
00331 {
00332 TCanvas *c1 = NULL;
00333
00334 int pad = 1;
00335 for(int e = eta0; e<=eta1;e++)
00336 {
00337 TH1F *h = findMip(e,e+etabin-1,equal);
00338 if(draw)
00339 {
00340 if(pad==1) {c1 = new TCanvas(); c1->Divide(2,2);}
00341 c1->cd(pad);
00342 h->Draw();
00343 funcFitBack->Draw("same");
00344 funcFitPeak->Draw("same");
00345 funcFit->Draw("same");
00346 pad++;
00347 if(pad>4) pad=1;
00348 }
00349 else
00350 {
00351 delete h;
00352 delete funcFit;
00353 delete funcFitPeak;
00354 delete funcFitBack;
00355 }
00356 }
00357 }
00358
00359 TH1F* StEmcMipMaker::compareToDb(char* timeStamp,int mode)
00360 {
00361
00362 StDbManager* mgr=StDbManager::Instance();
00363 StDbConfigNode* node=mgr->initConfig(dbCalibrations,dbEmc);
00364 StDbTable* table = node->addDbTable("bemcCalib");
00365 mgr->setRequestTime(timeStamp);
00366 mgr->fetchDbTable(table);
00367 emcCalib_st *calib = (emcCalib_st*) table->GetTable();
00368 if(!calib) return NULL;
00369
00370 char line[100];
00371 sprintf(line,"Comparison with DB - %s",timeStamp);
00372 TH1F *h = NULL;
00373 if(mode == 0) h = new TH1F("h",line,4800,0.5,4800.5);
00374 else h = new TH1F("h",line,100,0,4);
00375 for(int i=0;i<4800;i++)
00376 {
00377 int id = i+1;
00378 float gain = calib[0].AdcToE[i][1];
00379 float gain2 = getGain(id);
00380 float ratio = 0;
00381 if(gain!=0) ratio = gain2/gain;
00382 if(mode==0) h->Fill(id,ratio); else h->Fill(ratio);
00383 }
00384 return h;
00385 }
00386
00387 void StEmcMipMaker::saveToDb(char* timeStamp)
00388 {
00389 emcCalib_st tnew;
00390 for(int i=0;i<4800;i++)
00391 {
00392 int id = i+1;
00393 float gain = getGain(id);
00394 if(gain==0) tnew.Status[i] = 0; else tnew.Status[i] = 1;
00395 tnew.AdcToE[i][0] = 0;
00396 tnew.AdcToE[i][1] = gain;
00397 tnew.AdcToE[i][2] = 0;
00398 tnew.AdcToE[i][3] = 0;
00399 tnew.AdcToE[i][4] = 0;
00400 }
00401 StDbManager* mgr=StDbManager::Instance();
00402 StDbConfigNode* node=mgr->initConfig(dbCalibrations,dbEmc);
00403 StDbTable* table = node->addDbTable("bemcCalib");
00404 table->SetTable((char*)&tnew,1);
00405 mgr->setStoreTime(timeStamp);
00406 mgr->storeDbTable(table);
00407 return;
00408 }
00409
00410
00411
00412
00413
