00001 #include "StEmcPedestalMaker.h"
00002
00003 #include <fstream>
00004 using namespace std;
00005
00006 #include <TFile.h>
00007 #include <TROOT.h>
00008 #include <TF1.h>
00009
00010 #include <StEventTypes.h>
00011 #include <StEvent.h>
00012 #include <StEmcUtil/geometry/StEmcGeom.h>
00013 #include <StDbLib/StDbManager.hh>
00014 #include <StDbLib/StDbTable.h>
00015 #include <StDbLib/StDbConfigNode.hh>
00016 #include <tables/St_emcPed_Table.h>
00017 #include <tables/St_smdPed_Table.h>
00018 #include <StDaqLib/EMC/StEmcDecoder.h>
00019
00020 ClassImp(StEmcPedestalMaker)
00021
00022
00023 StEmcPedestalMaker::StEmcPedestalMaker(const Char_t *name)
00024 : StEmcCalibMaker(name)
00025 {
00026 setRange(300);
00027 setMaxTracks(100);
00028 mLastPedDate = 2000;
00029 mLastPedTime = 0;
00030
00031 setNPedEvents(2000);
00032 setSaveTables(false);
00033 setPedDiffSaveDB(2.0);
00034 setPedDiffSaveNum(3);
00035 setPedDiffSaveMinTime(60 * 60 * 24);
00036 setCompareLastTableDB(false);
00037 setPedCrateFilenameFormat("pedestal_crate0x%02x.dat");
00038 setBemcStatusFilename("bemcStatus.txt");
00039 setUseBemcStatus(true);
00040 }
00041
00042 StEmcPedestalMaker::~StEmcPedestalMaker()
00043 {
00044 }
00045
00046 Int_t StEmcPedestalMaker::Init()
00047 {
00048 mPedestal = new TH1F("mPed","",getNChannel(),0.5,getNChannel()+0.5);
00049 mRms = new TH1F("mRms","",getNChannel(),0.5,getNChannel()+0.5);
00050 mChi = new TH1F("mChi","",getNChannel(),0.5,getNChannel()+0.5);
00051 mStatus = new TH1F("mStatus","",getNChannel(),0.5,getNChannel()+0.5);
00052
00053 mSpecName = "mSpecPed";
00054 mAcceptName = "mAcceptPed";
00055
00056 mStarted = false;
00057
00058 return StEmcCalibMaker::Init();
00059 }
00060
00061 void StEmcPedestalMaker::Clear(Option_t *option)
00062 {
00063 }
00064
00065 Int_t StEmcPedestalMaker::Make()
00066 {
00067 if(!accept()) return kStOk;
00068
00069 if(!mStarted)
00070 {
00071 if(getTimeInterval(mLastPedDate,mLastPedTime)>mPedInterval)
00072 {
00073 mLastPedDate = getDate();
00074 mLastPedTime = getTime();
00075 mStarted = true;
00076 reset();
00077 }
00078 else
00079 {
00080 if(isDebug()) cout <<"Time remaining for a new pedestal run is "
00081 <<mPedInterval-getTimeInterval(mLastPedDate,mLastPedTime)
00082 <<" hours for detector number "<<getDetector()<<endl;
00083 }
00084 }
00085
00086 if(!mStarted) return kStOk;
00087
00088 for(int i=0;i<mNChannel;i++)
00089 {
00090 int id = i+1;
00091 float adc = (float) getCalib()->getADC(mDetector,id);
00092
00093 if(adc!=0) fill(id,adc);
00094 }
00095
00096 if(getNEvents()>getNPedEvents())
00097 {
00098 calcPedestals();
00099 savePedestals(mLastPedDate,mLastPedTime,isAutoSaveDB());
00100 if(isAutoSaveDB() || getSaveTables()) saveToDb(mLastPedDate,mLastPedTime);
00101 mStarted = false;
00102 }
00103
00104 return kStOK;
00105 }
00106
00107 Int_t StEmcPedestalMaker::Finish()
00108 {
00109 return StMaker::Finish();
00110 }
00111
00112 void StEmcPedestalMaker::calcPedestals()
00113 {
00114 cout <<"***** Calculating pedestals for detector "<<getDetector()<<" ...\n";
00115 mPedestal->Reset();
00116 mRms->Reset();
00117 mChi->Reset();
00118 mStatus->Reset();
00119 float left = 3;
00120 float right= 2;
00121
00122 int ngood=0,nped=0,nrms=0,nchi=0,nbad=0,nodata=0;
00123
00124 for(int id = 1;id<=getNChannel();id++) {
00125 TH1D *h = getSpec(id);
00126 int ibin = h->GetMaximumBin();
00127 float avg = (float)h->GetBinCenter(ibin);
00128 float max = (float)h->GetMaximum();
00129 float rms = 1.5;
00130 if(getDetector()>2) rms = h->GetRMS();
00131 float integral = (float)h->Integral();
00132 if((max!=0) && (integral > (getNPedEvents() / 2.0))) {
00133 TF1 func("ped","gaus(0)");
00134 float rmsInit = rms;
00135 func.SetParameter(0,max);
00136 func.SetParameter(1,avg);
00137 func.SetParameter(2,rms);
00138 func.SetParLimits(2,0,100000);
00139 float seed = avg;
00140 float fitleft = avg-left*rms;
00141 if(fitleft<0) fitleft = 0;
00142 float fitright = avg+right*rms;
00143 func.SetRange(fitleft,fitright);
00144
00145 int npt = (Int_t)((left+right+1.0)*rms);
00146 int ndg = (Int_t)((float)npt-3.0);
00147
00148 h->Fit(&func,"RQN");
00149 max = func.GetParameter(0);
00150 avg = func.GetParameter(1);
00151 rms = func.GetParameter(2);
00152 int status = 1;
00153 float chi = func.GetChisquare()/(float)ndg;
00154 float res = avg-seed;
00155
00156 if(avg<0) {status+= 2; nped++; avg = 0;}
00157 if(rms<0 || rms >7*rmsInit) {status+= 4; nrms++;}
00158 if(fabs(res)>1.5*rms) {status+= 8; nchi++;}
00159 if(status==1) ngood++; else nbad++;
00160 mPedestal->Fill((float)id,avg);
00161 mRms->Fill((float)id,rms);
00162 mChi->Fill((float)id,chi);
00163 mStatus->Fill((float)id,(float)status);
00164 if(status>1) cout <<"det = "<<getDetector()<<" id = "<<id <<" max = "<<seed<<" initRms = "<<rmsInit
00165 <<" peakY = "<<max
00166 <<" ped = "<<avg <<" res = " <<res<<" rms = "<<rms
00167 <<" chi = "<<chi<<" status = "<<status<<endl;
00168 } else {
00169 mPedestal->Fill((float)id,0);
00170 mRms->Fill((float)id,0);
00171 mChi->Fill((float)id,0);
00172 mStatus->Fill((float)id,0);
00173 nodata++;
00174 nbad++;
00175 cout << "det = " << getDetector() << " id = " << id << " peakY = "<< max
00176 << " ped = " << avg << " rms = " << rms
00177 << " int = " << integral << endl;
00178 }
00179 if(h) {delete h; h = 0;}
00180 }
00181 cout <<"nGood = "<<ngood<<"; nBad = "<<nbad<<": neg Ped = "<<nped<<", bad rms = "<<nrms<<", large res = "<<nchi << ", no data = " << nodata <<endl;
00182 }
00183
00184 void StEmcPedestalMaker::saveToDb(const Char_t *timeStamp, const Char_t *tableFilename) const {
00185 cout << "=================================================" << endl;
00186 cout << "Saving pedestal table for detector " << getDetector() << endl;
00187 cout << "TimeStamp = " << timeStamp << endl;
00188
00189 TString n[] = {"bemcPed", "bprsPed", "bsmdePed", "bsmdpPed"};
00190
00191 Bool_t saveThisTable = true;
00192 TString lastTableFilename = getLastTablePath();
00193 lastTableFilename += "/";
00194 lastTableFilename += n[getDetector() - 1];
00195 lastTableFilename += ".last";
00196 TString lastTableTimestampFilename = lastTableFilename + ".timestamp";
00197 lastTableFilename += ".root";
00198
00199 if (getCompareLastTableDB()) {
00200 cout << "Comparing to the last table file " << lastTableFilename << endl;
00201 TFile *lastTableFile = new TFile(lastTableFilename);
00202 if (lastTableFile && lastTableFile->IsOpen()) {
00203 cout << "Opened last table file " << lastTableFilename << endl;
00204 saveThisTable = false;
00205 ifstream ifstr_timestamp(lastTableTimestampFilename);
00206 if (ifstr_timestamp.good()) {
00207 Char_t lastTableTimestampStr[2048];
00208 ifstr_timestamp.getline(lastTableTimestampStr, 2048);
00209 ifstr_timestamp.close();
00210 TDatime tsLast(&lastTableTimestampStr[0]);
00211 TDatime tsCurrent(timeStamp);
00212 UInt_t timeLast = tsLast.Convert();
00213 UInt_t timeCurrent = tsCurrent.Convert();
00214 Double_t diffTime = difftime(timeCurrent, timeLast);
00215 cout << "Timestamps: Last: " << lastTableTimestampStr << ", " << tsLast.AsSQLString() << endl;
00216 cout << "Timestamps: Current: " << timeStamp << ", " << tsCurrent.AsSQLString() << "; Diff: " << diffTime << endl;
00217 if (diffTime > this->getPedDiffSaveMinTime()) {
00218 saveThisTable = true;
00219 cout << "Saving this table because min time passed" << endl;
00220 }
00221 } else {
00222 saveThisTable = true;
00223 cout << "Saving this table because cannot open last timestamp file " << lastTableTimestampFilename << endl;
00224 }
00225 const emcPed_st *ped_t_st = 0;
00226 const smdPed_st *ped_s_st = 0;
00227 if (getDetector() < 3) {
00228 const St_emcPed *pedT = (const St_emcPed *)lastTableFile->Get(n[getDetector() - 1]);
00229 ped_t_st = pedT ? pedT->GetTable() : 0;
00230 } else {
00231 const St_smdPed *pedT = (const St_smdPed *)lastTableFile->Get(n[getDetector() - 1]);
00232 ped_s_st = pedT ? pedT->GetTable() : 0;
00233 }
00234 Float_t maxPedDiff = 0;
00235 Float_t numPedDiff = 0;
00236 Int_t *bemcStatus = 0;
00237 if ((getDetector() == 1) && this->getUseBemcStatus() && this->getBemcStatusFilename()) {
00238 TString bemcStatusFilename = this->getBemcStatusFilename();
00239 ifstream ifstr(bemcStatusFilename);
00240 if (ifstr.good()) {
00241 cout << "Reading BEMC trigger status file " << bemcStatusFilename << endl;
00242 } else {
00243 cout << "Cannot open BEMC trigger status file! " << bemcStatusFilename << endl;
00244 }
00245 if (ifstr.good()) {
00246 bemcStatus = new Int_t[4800];
00247 if (bemcStatus) {
00248 for (Int_t i = 4800;i;bemcStatus[--i] = 1);
00249 while (ifstr.good()) {
00250 string token;
00251 do {
00252 if (token == "#") {
00253 char dummy[4096];
00254 ifstr.getline(dummy, sizeof(dummy));
00255 }
00256 ifstr >> token;
00257 } while (ifstr.good() && (token != "SoftId"));
00258 if (ifstr.good()) {
00259 if (token == "SoftId") {
00260 int softId, crate, crateSeq, unmaskTower, unmaskHT, unmaskPA;
00261 float ped;
00262 ifstr >> softId >> crate >> crateSeq >> unmaskTower >> unmaskHT >> unmaskPA >> ped;
00263 if ((softId >= 1) && (softId <= 4800)) {
00264 bemcStatus[softId - 1] = (unmaskTower && (unmaskHT || unmaskPA)) ? 1 : 0;
00265 if (bemcStatus[softId - 1] != 1) {
00266
00267 }
00268 }
00269 }
00270 }
00271 }
00272 }
00273 ifstr.close();
00274 }
00275 }
00276 for (Int_t i = 0;i < getNChannel();i++) {
00277 Int_t id = i + 1;
00278 Float_t pedNew = getPedestal(id);
00279 Float_t rmsNew = getRms(id);
00280
00281 Int_t statusNew = (int)getStatus(id);
00282 Float_t pedLast = 0;
00283 Float_t rmsLast = 0;
00284 Float_t chiLast = 0;
00285 Int_t statusLast = 0;
00286 if (getDetector() < 3) {
00287 if (ped_t_st) {
00288 pedLast = (short)ped_t_st->AdcPedestal[i] / 100.0;
00289 rmsLast = (short)ped_t_st->AdcPedestalRMS[i] / 100.0;
00290 chiLast = ped_t_st->ChiSquare[i];
00291 statusLast = ped_t_st->Status[i];
00292 }
00293 } else {
00294 if (ped_s_st) {
00295 pedLast = (short)ped_s_st->AdcPedestal[i][0] / 100.0;
00296 rmsLast = (short)ped_s_st->AdcPedestalRMS[i][0] / 100.0;
00297 statusLast = ped_s_st->Status[i];
00298 }
00299 }
00300 Float_t pedDiff = pedNew - pedLast;
00301 if (maxPedDiff < TMath::Abs(pedDiff)) maxPedDiff = TMath::Abs(pedDiff);
00302 if ((statusNew == 1)) {
00303 if (TMath::Abs(pedDiff) >= TMath::Abs(getPedDiffSaveDB() * rmsNew)) {
00304 TString statusLabel = "";
00305 if (!bemcStatus || ((getDetector() == 1) && bemcStatus && (id >= 1) && (id <= 4800) && (bemcStatus[id-1] == 1))) {
00306 numPedDiff++;
00307 } else {
00308 statusLabel = " (bad in trigger)";
00309 }
00310 cout << "det " << getDetector() << ", id " << id << statusLabel << ": ped diff " << pedDiff << ", last " << pedLast << " " << rmsLast << " " << (Int_t)statusLast << ", new " << pedNew << " " << rmsNew << " " << (Int_t)statusNew << endl;
00311 }
00312 }
00313 }
00314 if (bemcStatus) delete [] bemcStatus; bemcStatus = 0;
00315 cout << "max ped diff " << maxPedDiff << endl;
00316 cout << "num ped diff " << numPedDiff << endl;
00317 saveThisTable |= (numPedDiff >= this->getPedDiffSaveNum());
00318 lastTableFile->Close();
00319 if (!saveThisTable) cout << "This table does not need to be saved" << endl;
00320 }
00321 if (lastTableFile) delete lastTableFile; lastTableFile = 0;
00322 }
00323
00324 St_emcPed *pedT_t = new St_emcPed(n[getDetector() - 1].Data(), 1);
00325 emcPed_st *tnew = pedT_t ? pedT_t->GetTable() : 0;
00326
00327 St_smdPed *pedS_t = new St_smdPed(n[getDetector() - 1].Data(), 1);
00328 smdPed_st *snew = pedS_t ? pedS_t->GetTable() : 0;
00329
00330 for (int i = 0;i < getNChannel();i++) {
00331 int id = i + 1;
00332 float ped = getPedestal(id);
00333 float rms = getRms(id);
00334 float chi = getChi(id);
00335 int status = (int)getStatus(id);
00336 if (getDetector() < 3) {
00337 if (tnew) {
00338 tnew->Status[i] = (char)status;
00339 tnew->AdcPedestal[i] = (short)(ped * 100.0);
00340 tnew->AdcPedestalRMS[i] = (short)(rms * 100.0);
00341 tnew->ChiSquare[i] = chi;
00342
00343 }
00344 } else {
00345 if (snew) {
00346 snew->Status[i] = (char)status;
00347 snew->AdcPedestal[i][0] = (short)(ped * 100.0);
00348 snew->AdcPedestalRMS[i][0] = (short)(rms * 100.0);
00349 snew->AdcPedestal[i][1] = 0;
00350 snew->AdcPedestalRMS[i][1] = 0;
00351 snew->AdcPedestal[i][2] = 0;
00352 snew->AdcPedestalRMS[i][2] = 0;
00353
00354 }
00355 }
00356 }
00357 if (isAutoSaveDB() && (!getCompareLastTableDB() || saveThisTable)) {
00358 StDbManager* mgr = StDbManager::Instance();
00359 cout << "mgr = " << mgr << endl;
00360 StDbConfigNode* node = mgr ? mgr->initConfig(dbCalibrations, dbEmc) : 0;
00361 cout << "node = " << node << endl;
00362 StDbTable* table = node ? node->addDbTable(n[getDetector() - 1].Data()) : 0;
00363 cout << "table = " << table << endl;
00364 if (table) {
00365 table->setFlavor("ofl");
00366 if (getDetector() < 3) {
00367 table->SetTable((char*)tnew, 1);
00368 } else {
00369 table->SetTable((char*)snew, 1);
00370 }
00371 }
00372 cout << "table set" << endl;
00373 if (mgr && table) {
00374 cout << "setStoreTime " << timeStamp << endl;
00375 mgr->setStoreTime(timeStamp);
00376 cout << "Storing " << n[getDetector() - 1] << " " << timeStamp << endl;
00377 mgr->storeDbTable(table);
00378 cout << "Stored." << endl;
00379 }
00380 }
00381 if (getSaveTables() && (!getCompareLastTableDB() || saveThisTable) && tableFilename) {
00382 cout << "Saving DB table into " << tableFilename << endl;
00383 TFile *f = new TFile(tableFilename, "RECREATE");
00384 if (f) {
00385 if (getDetector() < 3) {
00386 pedT_t->AddAt(tnew, 0);
00387 pedT_t->Write();
00388 } else {
00389 pedS_t->AddAt(snew, 0);
00390 pedS_t->Write();
00391 }
00392 f->Close();
00393 delete f; f = 0;
00394 }
00395 }
00396 if (getCompareLastTableDB() && saveThisTable) {
00397 cout << "Saving last DB table into " << lastTableFilename << endl;
00398 TFile *f = new TFile(lastTableFilename, "RECREATE");
00399 if (f) {
00400 if (getDetector() < 3) {
00401 pedT_t->AddAt(tnew, 0);
00402 pedT_t->Write();
00403 } else {
00404 pedS_t->AddAt(snew, 0);
00405 pedS_t->Write();
00406 }
00407 f->Close();
00408 delete f; f = 0;
00409 ofstream ofstr_timestamp(lastTableTimestampFilename);
00410 ofstr_timestamp << timeStamp << endl;
00411 }
00412 if (getDetector() == 1) {
00413 StEmcDecoder *d = new StEmcDecoder();
00414 for (Int_t crate = 1;crate <= 30;crate++) {
00415 TString pedCrateFilename = getLastTablePath();
00416 pedCrateFilename += "/";
00417 pedCrateFilename += Form(getPedCrateFilenameFormat(), crate);
00418 TString pedCrateTimestampFilename = pedCrateFilename + ".timestamp";
00419 cout << "Saving pedestals for crate " << crate << ": " << pedCrateFilename << endl;
00420 ofstream ofstr(pedCrateFilename);
00421 for (Int_t ch = 0;ch < 160;ch++) {
00422 Int_t softId = -1;
00423 if (d) d->GetTowerIdFromCrate(crate, ch, softId);
00424 if ((softId >= 1) && (softId <= 4800)) {
00425 TString line = Form("%i %.2f %.2f", ch, getPedestal(softId), getRms(softId));
00426 ofstr << line.Data() << endl;
00427 }
00428 }
00429 ofstream ofstr_timestamp(pedCrateTimestampFilename);
00430 ofstr_timestamp << timeStamp << endl;
00431 }
00432 if (d) delete d; d = 0;
00433 }
00434 }
00435 }
00436
00437 void StEmcPedestalMaker::saveToDb(Int_t date, Int_t time) const
00438 {
00439 Int_t year = (Int_t)(date/10000);
00440 Int_t month = (Int_t)(date-year*10000)/100;
00441 Int_t day = (Int_t)(date-year*10000-month*100);
00442 Int_t hour = (Int_t)(time/10000);
00443 Int_t minute= (Int_t)(time-hour*10000)/100;
00444 Int_t second= (Int_t)(time-hour*10000-minute*100);
00445 Char_t ts[1024]; ts[0] = 0;
00446 Char_t tf[1024]; tf[0] = 0;
00447 TString d[] = {"bemc","bprs","bsmde","bsmdp"};
00448 TString n[] = {"bemcPed","bprsPed","bsmdePed","bsmdpPed"};
00449 sprintf(ts,"%04d-%02d-%02d %02d:%02d:%02d",year,month,day,hour,minute,second);
00450
00451 sprintf(tf, "%s/y3%s/%s.%08d.%06d.root", getTablesPath(), d[getDetector() - 1].Data(), n[getDetector() - 1].Data(), date, time);
00452 saveToDb(ts, tf);
00453 }
00454
00455 void StEmcPedestalMaker::savePedestals(Int_t date, Int_t time, Bool_t DB) const
00456 {
00457 Char_t ts[1024]; ts[0] = 0;
00458 TString n[] = {"bemcPed","bprsPed","bsmdePed","bsmdpPed"};
00459 if (DB) sprintf(ts,"%s/%s.%08d.%06d.root", getSavePath(), n[getDetector()-1].Data(), date, time);
00460 else sprintf(ts,"%s/%s.%08d.%06d.NO_DB.root", getSavePath(), n[getDetector()-1].Data(), date, time);
00461 cout << "Saving pedestal histograms to " << ts << endl;
00462 TFile *f = new TFile(ts,"RECREATE");
00463 if(getSpec()) getSpec()->Write();
00464 if(mPedestal) mPedestal->Write();
00465 if(mRms) mRms->Write();
00466 if(mChi) mChi->Write();
00467 if(mStatus) mStatus->Write();
00468 f->Close();
00469 delete f;
00470 }
00471
00472 void StEmcPedestalMaker::loadPedestals(const Char_t* file)
00473 {
00474 TFile *f = new TFile(file);
00475 if(getSpec()) getSpec()->Reset();
00476 if(mPedestal) mPedestal->Reset();
00477 if(mRms) mRms->Reset();
00478 if(mChi) mChi->Reset();
00479 if(mStatus) mStatus->Reset();
00480 TH2F* h = (TH2F*)f->Get("mSpec;1");
00481 if(h && getSpec()) getSpec()->Add(h,1);
00482 TH1F *g=(TH1F*)f->Get("mPed;1");
00483 if(g && mPedestal) mPedestal->Add(g,1);
00484 g=(TH1F*)f->Get("mRms;1");
00485 if(g && mRms) mRms->Add(g,1);
00486 g=(TH1F*)f->Get("mChi;1");
00487 if(g && mChi) mChi->Add(g,1);
00488 g=(TH1F*)f->Get("mStatus;1");
00489 if(g && mStatus) mStatus->Add(g,1);
00490 f->Close();
00491 delete f;
00492 return;
00493 }
00494
