Home » Groups » STAR Protected

Linear Fit of Time Dependency in BEMC Gains

Updated on Sun, 2009-10-18 16:10. Originally created by wwitzke on 2009-10-18 14:41.

To better characterize the time-dependence in the BEMC in 2006, we have fit each of the gains versus time data sets for each eta bin with a linear (non-constant) fit. These fits show marked improvement in the Chi^2 tests. The original analysis can be found here. A complete set of P1 fits, along with Chi^2 and slopes, can be found here.

Here are some samples from the raw P1 fits, along with the Chi^2 and slope graphs:

Figure 1: 2006 Eta Bin 9 mip Peak Positions Over Time, P1 Fit

Figure 2: 2006 P1 Fit Chi^2 versus Eta Bin

Figure 3: 2006 Eta Bin 9 mip Peak Positions Over Time Using Weighted Average Over 4 Values, P1 Fit

Figure 4: 2006 P1 Fit Chi^2 versus Eta Bin for Weighted Averages

For comparison, especially in the Chi^2 test, the original, constant fits are shown below:

Figure 5: 2006 Eta Bin 9 mip Peak Positions Over Time

Figure 6: 2006 Chi^2 versus Eta Bin

Figure 7: 2006 Eta Bin 9 mip Peak Positions Over Time Using Weighted Average Over 4 Values

Figure 8: 2006 Chi^2 versus Eta Bin for Weighted Averages

We also plotted the slopes of each P1 fit for each eta bin. The results for the non-weighted and weighted data sets are shown below:

Figure 9: 2006 P1 Slopes versus Eta Bin

Figure 10: 2006 P1 Slopes versus Eta Bin for Weighted Averages

Incidentally, the slopes for the rebinned data are significantly different from the slopes for the raw data only because the X-axis has been compressed, with each 4 X-values of the raw data graph becoming a single X-value in the graphs of the rebinned data. This causes the rebinned slopes to scale up by a factor of four.

For each of these plots, eta bin 1 corresponds to eta = -1, eta = 0 occurs between eta bins 20 and 21, and eta bin 40 corresponds to eta = 1.

For ease of use, a table of mappings is provided here:

2006 Fills

1 0
2 7603
3 7622
4 7627
5 7630
6 7632
7 7637
8 7639
9 7641
10 7642
11 7645
12 7646
13 7651
14 7652
15 7654
16 7655
17 7657
18 7658
19 7662
20 7671
21 7672
22 7673
23 7674
24 7681
25 7685
26 7688
27 7691
28 7697
29 7699
30 7718
31 7722
32 7724
33 7725
34 7729
35 7739
36 7740
37 7744
38 7745
39 7753
40 7756
41 7757
42 7780
43 7781
44 7785
45 7786
46 7788
47 7789
48 7790
49 7791
50 7792
51 7794
52 7795
53 7796
54 7797
55 7800
56 7803
57 7804
58 7805
59 7810
60 7811
61 7815
62 7817
63 7820
64 7823
65 7824
66 7825
67 7826
68 7827
69 7830
70 7831
71 7847
72 7850
73 7851
74 7852
75 7853
76 7855
77 7856
78 7858
79 7863
80 7864
81 7865
82 7871
83 7872
84 7883
85 7886
86 7887
87 7889
88 7890
89 7891
90 7892
91 7893
92 7896
93 7901
94 7908
95 7909
96 7911
97 7913
98 7915
99 7916
100 7918
101 7921
102 7922
103 7926
104 7944
105 7946
106 7949
107 7951
108 7952
109 7954
110 7957

We have also started looking at time dependency in the pedistals. This is very preliminary, but we have established that the method for calculating pedistals is robust. The code for generating pedistals, contained in StEmcPedistalMaker.cxx, is:

#include "StEmcPedestalMaker.h"

#include <fstream>
using namespace std;

#include <TFile.h>
#include <TROOT.h>
#include <TF1.h>

#include <StEventTypes.h>
#include <StEvent.h>
#include <StEmcUtil/geometry/StEmcGeom.h>
#include <StDbLib/StDbManager.hh>
#include <StDbLib/StDbTable.h>
#include <StDbLib/StDbConfigNode.hh>
#include <tables/St_emcPed_Table.h>
#include <tables/St_smdPed_Table.h>
#include <StDaqLib/EMC/StEmcDecoder.h>

ClassImp(StEmcPedestalMaker)

//_____________________________________________________________________________
StEmcPedestalMaker::StEmcPedestalMaker(const Char_t *name)
: StEmcCalibMaker(name)
{
setRange(300);
setMaxTracks(100);
mLastPedDate = 2000;
mLastPedTime = 0;

setNPedEvents(2000);
setSaveTables(false);
setPedDiffSaveDB(2.0);
setPedDiffSaveNum(3);
setPedDiffSaveMinTime(60 * 60 * 24);
setCompareLastTableDB(false);
setPedCrateFilenameFormat("pedestal_crate0x%02x.dat");
setBemcStatusFilename("bemcStatus.txt");
setUseBemcStatus(true);
}
//_____________________________________________________________________________
StEmcPedestalMaker::~StEmcPedestalMaker()
{
}
//_____________________________________________________________________________
Int_t StEmcPedestalMaker::Init()
{
mPedestal = new TH1F("mPed","",getNChannel(),0.5,getNChannel()+0.5);
mRms = new TH1F("mRms","",getNChannel(),0.5,getNChannel()+0.5);
mChi = new TH1F("mChi","",getNChannel(),0.5,getNChannel()+0.5);
mStatus = new TH1F("mStatus","",getNChannel(),0.5,getNChannel()+0.5);

mSpecName = "mSpecPed";
mAcceptName = "mAcceptPed";

mStarted = false;

return StEmcCalibMaker::Init();
}
//_____________________________________________________________________________
void StEmcPedestalMaker::Clear(Option_t *option)
{
}
//_____________________________________________________________________________
Int_t StEmcPedestalMaker::Make()
{
LOG_DEBUG << "starting make" << endm;
if(!accept()) return kStOk;

if(!mStarted)
{
if(getTimeInterval(mLastPedDate,mLastPedTime)>mPedInterval)
{
mLastPedDate = getDate();
mLastPedTime = getTime();
mStarted = true;
reset();
}
else
{
if(isDebug()) cout <<"Time remaining for a new pedestal run is "
<<mPedInterval-getTimeInterval(mLastPedDate,mLastPedTime)
<<" hours for detector number "<<getDetector()<<endl;
}
}

if(!mStarted) return kStOk;
//cout <<"Started\n";
for(int i=0;i<mNChannel;i++)
{
int id = i+1;
float adc = (float) getCalib()->getADC(mDetector,id);
LOG_DEBUG <<"Detector = "<<mDetector<<" id = "<<id<<" adc = "<<adc<<endm;
if(adc!=0) fill(id,adc);
}

if(getNEvents()>getNPedEvents())
{
calcPedestals();
savePedestals(mLastPedDate,mLastPedTime,isAutoSaveDB());
if(isAutoSaveDB() || getSaveTables()) saveToDb(mLastPedDate,mLastPedTime);
mStarted = false;
}

return kStOK;
}
//_____________________________________________________________________________
Int_t StEmcPedestalMaker::Finish()
{
return StMaker::Finish();
}
//_____________________________________________________________________________
void StEmcPedestalMaker::calcPedestals()
{
cout <<"***** Calculating pedestals for detector "<<getDetector()<<" ...\n";
mPedestal->Reset();
mRms->Reset();
mChi->Reset();
mStatus->Reset();
float left = 3;
float right= 2;

int ngood=0,nped=0,nrms=0,nchi=0,nbad=0,nodata=0;

for(int id = 1;id<=getNChannel();id++) {
TH1D *h = getSpec(id);
int ibin = h->GetMaximumBin();
float avg = (float)h->GetBinCenter(ibin);
float max = (float)h->GetMaximum();
float rms = 1.5;
if(getDetector()>2) rms = h->GetRMS();
float integral = (float)h->Integral();
if((max!=0) && (integral > (getNPedEvents() / 2.0))) {
TF1 func("ped","gaus(0)");
float rmsInit = rms;
func.SetParameter(0,max);
func.SetParameter(1,avg);
func.SetParameter(2,rms);
func.SetParLimits(2,0,100000);
float seed = avg;
float fitleft = avg-left*rms;
if(fitleft<0) fitleft = 0;
float fitright = avg+right*rms;
func.SetRange(fitleft,fitright);

int npt = (Int_t)((left+right+1.0)*rms);
int ndg = (Int_t)((float)npt-3.0);

h->Fit(&func,"RQN"); // pre fit
max = func.GetParameter(0);
avg = func.GetParameter(1);
rms = func.GetParameter(2);
int status = 1; // data present
float chi = func.GetChisquare()/(float)ndg;
float res = avg-seed;

if(avg<0) {status+= 2; nped++; avg = 0;}// negative pedestal
if(rms<0 || rms >7*rmsInit) {status+= 4; nrms++;}// bad rms
if(fabs(res)>1.5*rms) {status+= 8; nchi++;}// large distance to seed
if(status==1) ngood++; else nbad++;
mPedestal->Fill((float)id,avg);
mRms->Fill((float)id,rms);
mChi->Fill((float)id,chi);
mStatus->Fill((float)id,(float)status);
if(status>1) cout <<"det = "<<getDetector()<<" id = "<<id <<" max = "<<seed<<" initRms = "<<rmsInit
<<" peakY = "<<max
<<" ped = "<<avg <<" res = " <<res<<" rms = "<<rms
<<" chi = "<<chi<<" status = "<<status<<endl;
} else {
mPedestal->Fill((float)id,0);
mRms->Fill((float)id,0);
mChi->Fill((float)id,0);
mStatus->Fill((float)id,0);
nodata++;
nbad++;
cout << "det = " << getDetector() << " id = " << id << " peakY = "<< max
<< " ped = " << avg << " rms = " << rms
<< " int = " << integral << endl;
}
if(h) {delete h; h = 0;}
}
cout <<"nGood = "<<ngood<<"; nBad = "<<nbad<<": neg Ped = "<<nped<<", bad rms = "<<nrms<<", large res = "<<nchi << ", no data = " << nodata <<endl;
}
//_____________________________________________________________________________
void StEmcPedestalMaker::saveToDb(const Char_t *timeStamp, const Char_t *tableFilename) const {
cout << "=================================================" << endl;
cout << "Saving pedestal table for detector " << getDetector() << endl;
cout << "TimeStamp = " << timeStamp << endl;

TString n[] = {"bemcPed", "bprsPed", "bsmdePed", "bsmdpPed"};

Bool_t saveThisTable = true;
TString lastTableFilename = getLastTablePath();
lastTableFilename += "/";
lastTableFilename += n[getDetector() - 1];
lastTableFilename += ".last";
TString lastTableTimestampFilename = lastTableFilename + ".timestamp";
lastTableFilename += ".root";

if (getCompareLastTableDB()) {
cout << "Comparing to the last table file " << lastTableFilename << endl;
TFile *lastTableFile = new TFile(lastTableFilename);
if (lastTableFile && lastTableFile->IsOpen()) {
cout << "Opened last table file " << lastTableFilename << endl;
saveThisTable = false;
ifstream ifstr_timestamp(lastTableTimestampFilename);
if (ifstr_timestamp.good()) {
Char_t lastTableTimestampStr[2048];
ifstr_timestamp.getline(lastTableTimestampStr, 2048);
ifstr_timestamp.close();
TDatime tsLast(&lastTableTimestampStr[0]);
TDatime tsCurrent(timeStamp);
UInt_t timeLast = tsLast.Convert();
UInt_t timeCurrent = tsCurrent.Convert();
Double_t diffTime = difftime(timeCurrent, timeLast);
cout << "Timestamps: Last: " << lastTableTimestampStr << ", " << tsLast.AsSQLString() << endl;
cout << "Timestamps: Current: " << timeStamp << ", " << tsCurrent.AsSQLString() << "; Diff: " << diffTime << endl;
if (diffTime > this->getPedDiffSaveMinTime()) {
saveThisTable = true;
cout << "Saving this table because min time passed" << endl;
}
} else {
saveThisTable = true;
cout << "Saving this table because cannot open last timestamp file " << lastTableTimestampFilename << endl;
}
const emcPed_st *ped_t_st = 0;
const smdPed_st *ped_s_st = 0;
if (getDetector() < 3) {
const St_emcPed *pedT = (const St_emcPed *)lastTableFile->Get(n[getDetector() - 1]);
ped_t_st = pedT ? pedT->GetTable() : 0;
} else {
const St_smdPed *pedT = (const St_smdPed *)lastTableFile->Get(n[getDetector() - 1]);
ped_s_st = pedT ? pedT->GetTable() : 0;
}
Float_t maxPedDiff = 0;
Float_t numPedDiff = 0;
Int_t *bemcStatus = 0;
if ((getDetector() == 1) && this->getUseBemcStatus() && this->getBemcStatusFilename()) {
TString bemcStatusFilename = this->getBemcStatusFilename();
ifstream ifstr(bemcStatusFilename);
if (ifstr.good()) {
cout << "Reading BEMC trigger status file " << bemcStatusFilename << endl;
} else {
cout << "Cannot open BEMC trigger status file! " << bemcStatusFilename << endl;
}
if (ifstr.good()) {
bemcStatus = new Int_t[4800];
if (bemcStatus) {
for (Int_t i = 4800;i;bemcStatus[--i] = 1);
while (ifstr.good()) {
string token;
do {
if (token == "#") {
char dummy[4096];
ifstr.getline(dummy, sizeof(dummy));
}
ifstr >> token;
} while (ifstr.good() && (token != "SoftId"));
if (ifstr.good()) {
if (token == "SoftId") {
int softId, crate, crateSeq, unmaskTower, unmaskHT, unmaskPA;
float ped;
ifstr >> softId >> crate >> crateSeq >> unmaskTower >> unmaskHT >> unmaskPA >> ped;
if ((softId >= 1) && (softId <= 4800)) {
bemcStatus[softId - 1] = (unmaskTower && (unmaskHT || unmaskPA)) ? 1 : 0;
if (bemcStatus[softId - 1] != 1) {
//cout << "tower " << softId << " masked out: " << unmaskTower << " " << unmaskHT << " " << unmaskPA << endl;
}
}
}
}
}
}
ifstr.close();
}
}
for (Int_t i = 0;i < getNChannel();i++) {
Int_t id = i + 1;
Float_t pedNew = getPedestal(id);
Float_t rmsNew = getRms(id);
Float_t chiNew = getChi(id);
Int_t statusNew = (int)getStatus(id);
Float_t pedLast = 0;
Float_t rmsLast = 0;
Float_t chiLast = 0;
Int_t statusLast = 0;
if (getDetector() < 3) {
if (ped_t_st) {
pedLast = (short)ped_t_st->AdcPedestal[i] / 100.0;
rmsLast = (short)ped_t_st->AdcPedestalRMS[i] / 100.0;
chiLast = ped_t_st->ChiSquare[i];
statusLast = ped_t_st->Status[i];
}
} else {
if (ped_s_st) {
pedLast = (short)ped_s_st->AdcPedestal[i][0] / 100.0;
rmsLast = (short)ped_s_st->AdcPedestalRMS[i][0] / 100.0;
statusLast = ped_s_st->Status[i];
}
}
Float_t pedDiff = pedNew - pedLast;
if (maxPedDiff < TMath::Abs(pedDiff)) maxPedDiff = TMath::Abs(pedDiff);
if (/*(statusLast == 1) || */(statusNew == 1)) {
if (TMath::Abs(pedDiff) >= TMath::Abs(getPedDiffSaveDB() * rmsNew)) {
TString statusLabel = "";
if (!bemcStatus || ((getDetector() == 1) && bemcStatus && (id >= 1) && (id <= 4800) && (bemcStatus[id-1] == 1))) {
numPedDiff++;
} else {
statusLabel = " (bad in trigger)";
}
cout << "det " << getDetector() << ", id " << id << statusLabel << ": ped diff " << pedDiff << ", last " << pedLast << " " << rmsLast << " " << (Int_t)statusLast << ", new " << pedNew << " " << rmsNew << " " << (Int_t)statusNew << endl;
}
}
}
if (bemcStatus) delete [] bemcStatus; bemcStatus = 0;
cout << "max ped diff " << maxPedDiff << endl;
cout << "num ped diff " << numPedDiff << endl;
saveThisTable |= (numPedDiff >= this->getPedDiffSaveNum());
lastTableFile->Close();
if (!saveThisTable) cout << "This table does not need to be saved" << endl;
}
if (lastTableFile) delete lastTableFile; lastTableFile = 0;
}

St_emcPed *pedT_t = new St_emcPed(n[getDetector() - 1].Data(), 1);
emcPed_st *tnew = pedT_t ? pedT_t->GetTable() : 0;

St_smdPed *pedS_t = new St_smdPed(n[getDetector() - 1].Data(), 1);
smdPed_st *snew = pedS_t ? pedS_t->GetTable() : 0;

for (int i = 0;i < getNChannel();i++) {
int id = i + 1;
float ped = getPedestal(id);
float rms = getRms(id);
float chi = getChi(id);
int status = (int)getStatus(id);
if (getDetector() < 3) {
if (tnew) {
tnew->Status[i] = (char)status;
tnew->AdcPedestal[i] = (short)(ped * 100.0);
tnew->AdcPedestalRMS[i] = (short)(rms * 100.0);
tnew->ChiSquare[i] = chi;
//cout << "tnew i = " << i << ": Status " << (Int_t)tnew->Status[i] << ", ped " << tnew->AdcPedestal[i] << ", rms " << tnew->AdcPedestalRMS[i] << ", chi " << tnew->ChiSquare[i] << endl;
}
} else {
if (snew) {
snew->Status[i] = (char)status;
snew->AdcPedestal[i][0] = (short)(ped * 100.0);
snew->AdcPedestalRMS[i][0] = (short)(rms * 100.0);
snew->AdcPedestal[i][1] = 0;
snew->AdcPedestalRMS[i][1] = 0;
snew->AdcPedestal[i][2] = 0;
snew->AdcPedestalRMS[i][2] = 0;
//cout << "snew i = " << i << ": Status " << (Int_t)snew->Status[i] << ", ped " << snew->AdcPedestal[i][0] << ", rms " << snew->AdcPedestalRMS[i][0] << endl;
}
}
}
if (isAutoSaveDB() && (!getCompareLastTableDB() || saveThisTable)) {
StDbManager* mgr = StDbManager::Instance();
cout << "mgr = " << mgr << endl;
StDbConfigNode* node = mgr ? mgr->initConfig(dbCalibrations, dbEmc) : 0;
cout << "node = " << node << endl;
StDbTable* table = node ? node->addDbTable(n[getDetector() - 1].Data()) : 0;
cout << "table = " << table << endl;
if (table) {
table->setFlavor("ofl");
if (getDetector() < 3) {
table->SetTable((char*)tnew, 1);
} else {
table->SetTable((char*)snew, 1);
}
}
cout << "table set" << endl;
if (mgr && table) {
cout << "setStoreTime " << timeStamp << endl;
mgr->setStoreTime(timeStamp);
cout << "Storing " << n[getDetector() - 1] << " " << timeStamp << endl;
mgr->storeDbTable(table);
cout << "Stored." << endl;
}
}
if (getSaveTables() && (!getCompareLastTableDB() || saveThisTable) && tableFilename) {
cout << "Saving DB table into " << tableFilename << endl;
TFile *f = new TFile(tableFilename, "RECREATE");
if (f) {
if (getDetector() < 3) {
pedT_t->AddAt(tnew, 0);
pedT_t->Write();
} else {
pedS_t->AddAt(snew, 0);
pedS_t->Write();
}
f->Close();
delete f; f = 0;
}
}
if (getCompareLastTableDB() && saveThisTable) {
cout << "Saving last DB table into " << lastTableFilename << endl;
TFile *f = new TFile(lastTableFilename, "RECREATE");
if (f) {
if (getDetector() < 3) {
pedT_t->AddAt(tnew, 0);
pedT_t->Write();
} else {
pedS_t->AddAt(snew, 0);
pedS_t->Write();
}
f->Close();
delete f; f = 0;
ofstream ofstr_timestamp(lastTableTimestampFilename);
ofstr_timestamp << timeStamp << endl;
}
if (getDetector() == 1) {
StEmcDecoder *d = new StEmcDecoder();
for (Int_t crate = 1;crate <= 30;crate++) {
TString pedCrateFilename = getLastTablePath();
pedCrateFilename += "/";
pedCrateFilename += Form(getPedCrateFilenameFormat(), crate);
TString pedCrateTimestampFilename = pedCrateFilename + ".timestamp";
cout << "Saving pedestals for crate " << crate << ": " << pedCrateFilename << endl;
ofstream ofstr(pedCrateFilename);
for (Int_t ch = 0;ch < 160;ch++) {
Int_t softId = -1;
if (d) d->GetTowerIdFromCrate(crate, ch, softId);
if ((softId >= 1) && (softId <= 4800)) {
TString line = Form("%i %.2f %.2f", ch, getPedestal(softId), getRms(softId));
ofstr << line.Data() << endl;
}
}
ofstream ofstr_timestamp(pedCrateTimestampFilename);
ofstr_timestamp << timeStamp << endl;
}
if (d) delete d; d = 0;
}
}
}
//_____________________________________________________________________________
void StEmcPedestalMaker::saveToDb(Int_t date, Int_t time) const
{
Int_t year = (Int_t)(date/10000);
Int_t month = (Int_t)(date-year*10000)/100;
Int_t day = (Int_t)(date-year*10000-month*100);
Int_t hour = (Int_t)(time/10000);
Int_t minute= (Int_t)(time-hour*10000)/100;
Int_t second= (Int_t)(time-hour*10000-minute*100);
Char_t ts[1024]; ts[0] = 0;
Char_t tf[1024]; tf[0] = 0;
TString d[] = {"bemc","bprs","bsmde","bsmdp"};
TString n[] = {"bemcPed","bprsPed","bsmdePed","bsmdpPed"};
sprintf(ts,"%04d-%02d-%02d %02d:%02d:%02d",year,month,day,hour,minute,second);
//sprintf(ts,"%04d%02d%02d%02d%02d%02d",year,month,day,hour,minute,second);
sprintf(tf, "%s/y3%s/%s.%08d.%06d.root", getTablesPath(), d[getDetector() - 1].Data(), n[getDetector() - 1].Data(), date, time);
saveToDb(ts, tf);
}
//_____________________________________________________________________________
void StEmcPedestalMaker::savePedestals(Int_t date, Int_t time, Bool_t DB) const
{
Char_t ts[1024]; ts[0] = 0;
TString n[] = {"bemcPed","bprsPed","bsmdePed","bsmdpPed"};
if (DB) sprintf(ts,"%s/%s.%08d.%06d.root", getSavePath(), n[getDetector()-1].Data(), date, time);
else sprintf(ts,"%s/%s.%08d.%06d.NO_DB.root", getSavePath(), n[getDetector()-1].Data(), date, time);
cout << "Saving pedestal histograms to " << ts << endl;
TFile *f = new TFile(ts,"RECREATE");
if(getSpec()) getSpec()->Write();
if(mPedestal) mPedestal->Write();
if(mRms) mRms->Write();
if(mChi) mChi->Write();
if(mStatus) mStatus->Write();
f->Close();
delete f;
}
//_____________________________________________________________________________
void StEmcPedestalMaker::loadPedestals(const Char_t* file)
{
TFile *f = new TFile(file);
if(getSpec()) getSpec()->Reset();
if(mPedestal) mPedestal->Reset();
if(mRms) mRms->Reset();
if(mChi) mChi->Reset();
if(mStatus) mStatus->Reset();
TH2F* h = (TH2F*)f->Get("mSpec;1");
if(h && getSpec()) getSpec()->Add(h,1);
TH1F *g=(TH1F*)f->Get("mPed;1");
if(g && mPedestal) mPedestal->Add(g,1);
g=(TH1F*)f->Get("mRms;1");
if(g && mRms) mRms->Add(g,1);
g=(TH1F*)f->Get("mChi;1");
if(g && mChi) mChi->Add(g,1);
g=(TH1F*)f->Get("mStatus;1");
if(g && mStatus) mStatus->Add(g,1);
f->Close();
delete f;
return;
}
Groups: