L rich newdecodeRICH
From Hall A Wiki
Revision as of 12:19, 27 September 2009 by Cisbani (Talk | contribs) (New page: newdecodeRICH.C (replay pedestal raw data) <pre> // // Run analyzer // > .L decodeRICH.C // > decodeRICH("/home/adaq/ped159",1000) // // #include <string.h> #include <stdio.h> #includ...)
newdecodeRICH.C (replay pedestal raw data)
//
// Run analyzer
// > .L decodeRICH.C
// > decodeRICH("/home/adaq/ped159",1000)
//
//
#include <string.h>
#include <stdio.h>
#include <TROOT.h>
#include <TMath.h>
#include <TStyle.h>
#include <TH1D.h>
#include <TH1I.h>
#include <TH2I.h>
#include <TTree.h>
#include <TCanvas.h>
#include <TArray.h>
#include <TProfile.h>
#include <TGraph.h>
#include "THaCodaDecoder.h"
#include "THaCodaFile.h"
#define EVT_HEADER_TAG 0x73530000
#define EVT_HEADER_OFFSET 2
#define ADC_HEADER_TAG 0xcd000000
#define ADC_HEADER_NUL 0xca000000
#define ADC_HEADER_OFFSET 3
#define EVT_COUNT_MASK 0xffff
#define EVT_COUNT_OFFSET 1
#define ADC_VALUE_MASK 0xfff
#define ADC_CH_MASK 0x7ff000
#define ADC_VALID_BIT 0x40000000
#define ADC_OVER_BIT 0x80000000
#define CH_COUNT_MASK 0x7FF
#define NUMBER_OF_ROC 2
#define NUMBER_ADC 20
#define CHANNEL_X_CARD 240
#define CARD_X_ADC 2
#define CHANNEL_X_ADC (CHANNEL_X_CARD*CARD_X_ADC)
#define ADC_X_PANEL 8
#define ACCEPTED_NSIGMA 3
// level of noise above pedestal
#define XPAD_X_ADC 6
#define YPAD_X_ADC 80
#define ADC_ALONG_X NUMBER_ADC
#define ADC_ALONG_Y 1
//
// Return the (absolute row) that serve the given channel
Int_t ch2arow(Int_t ch) {
return ch/CHANNEL_X_ADC;
}
//
// Interlieved ADC in NUMBER_OF_ROC crates
Int_t ch2ADC(Int_t ch) {
Int_t ach;
ach = ch/CHANNEL_X_ADC;
ach = ach/NUMBER_OF_ROC + NUMBER_ADC*(ach % NUMBER_OF_ROC);
return ach;
}
//
// Return the channel relative to the corresponding ADC
Int_t ch2rch(Int_t ch) {
return ch % CHANNEL_X_ADC;
}
//
// Return the panel that contains the given channel
Int_t ch2panel(Int_t ch) {
return ch/(CHANNEL_X_ADC*ADC_X_PANEL);
}
//
// x = long side
// y = short side
Int_t xy2ch(Int_t x, Int_t y) {
// Int_t xch[6]={4,2,0,1,3,5};
Int_t xch[6]={5,3,1,0,2,4};
return x/XPAD_X_ADC*CHANNEL_X_ADC + y*XPAD_X_ADC + xch[x % XPAD_X_ADC];
}
//
// x = long side
Int_t ch2x(Int_t ch) {
Int_t chx[6]={2,3,1,4,0,5}; //do not know which one is correct!
//Int_t chx[6]={3,2,4,1,5,0};
Int_t adc = ch2arow(ch);
return (adc*XPAD_X_ADC + chx[ch % XPAD_X_ADC]);
}
//
// y = short side
Int_t ch2y(Int_t ch) {
return ((ch / XPAD_X_ADC) % YPAD_X_ADC);
}
//
// Return the row relative to the panel containing the channel
Int_t ch2prow(Int_t ch) {
return (ch - ch2panel(ch)*CHANNEL_X_ADC*ADC_X_PANEL)/CHANNEL_X_ADC;
}
//
// Return the gassiplex card 0 or 1 of the corresponding ADC
Int_t ch2rcard(Int_t ch) {
return (ch-ch2arow(ch)*CHANNEL_X_ADC)/CHANNEL_X_CARD;
}
//
// Return the ROC number of the given channel
Int_t ch2roc(Int_t ch) {
Int_t ach;
ach = ch/CHANNEL_X_ADC % NUMBER_OF_ROC;
return ach;
}
//=====================================================
//
// event_preset : number of events to process (minimum)
// channel_offset : first channel in adc (since we can acquire half the detector at present)
// mode = 0 pedestal mode, mode = 1 acq mode
Int_t newdecodeRICH(const Int_t nrun=0, Int_t event_preset=1000, Int_t mode=0, Double_t nsigma=2.)
{
gStyle->SetPalette(1);
THaCodaFile * coda = NULL;
Int_t * buffer;
THaCodaDecoder * ev = NULL;
TTree * tdata = NULL;
TProfile *pdata = NULL;
TH2I *mdata = NULL;
TGraph *grrms, *grmean, *grscat;
Int_t channel_offset, roc;
Int_t iroc, ie;
Int_t evt_count=0;
Int_t evt_count_valid[NUMBER_OF_ROC], evt_count_bad[NUMBER_OF_ROC];
Short_t adc_count;
Short_t ch_count;
Int_t adc_ch;
Short_t adc_ch_value;
Short_t adc_ch_flag;
Int_t eff_min_ch, eff_max_ch; // min and max channel analyzed
// tree variables
Int_t evt; // event number
Int_t nch; // number of channel above thr in event
Int_t tcharge; // total charge in event
Short_t charge[NUMBER_OF_ROC*NUMBER_ADC*CHANNEL_X_ADC]; // value of adc for each ch.
Short_t ch[NUMBER_OF_ROC*NUMBER_ADC*CHANNEL_X_ADC]; // value of adc for each ch.
Int_t total_channel;
total_channel = NUMBER_OF_ROC*NUMBER_ADC*CHANNEL_X_ADC;
eff_max_ch = 0;
eff_min_ch = total_channel+1;
int found = 0;
// const char* path = 0;
FILE *fd;
char filename[100];
while( found==0 ) {
if( nrun<0 ) break;
int disk = 1, disk_max = 5; // 5 = /cache/mss
while (disk <= disk_max) {
sprintf(filename,"/adaql2/data%d/e06010_%d.dat.0",disk,nrun);
if (disk == 5) {
sprintf(filename,"/cache/mss/halla/e06010/raw/e06010_%d.dat.0",nrun);
}
fd = fopen(filename,"r");
if (fd != NULL) {
found = 1;
cout << "Will analyze file "<<filename<<endl;
fclose(fd);
break;
}
disk++;
}
if ( disk >= disk_max && !found ) {
cout << "File not found, try to retrieve from cache using:" << endl;
cout << Form(" jcache /mss/halla/e06010/raw/e06010_%d.dat.0",nrun) << endl;
cout << " check the status of jcache at: http://scicomp.jlab.org/scicomp/index.jsp" << endl;
return -1;
}
}
coda = new THaCodaFile(filename);
if (coda != NULL) {
ev = new THaCodaDecoder;
//
tdata = new TTree("tdata", Form("%s decoded pedestal data",filename));
tdata->Branch("evt",&evt,"evt/I");
tdata->Branch("nch",&nch,"nch/I");
tdata->Branch("tcharge",&tcharge,"tcharge/I");
tdata->Branch("ch",ch,"ch[nch]/S");
tdata->Branch("charge",charge,"charge[nch]/S");
//
pdata = new TProfile("pdata","Pedestal / Noise",total_channel,-0.5,total_channel-0.5,"s");
//
mdata = new TH2I("mdata","RICH Intensity Map",
(XPAD_X_ADC*ADC_ALONG_X),-0.5,(XPAD_X_ADC*ADC_ALONG_X)-0.5, (YPAD_X_ADC*ADC_ALONG_Y),-0.5,(YPAD_X_ADC*ADC_ALONG_Y)-0.5);
//
for (iroc=0; iroc<NUMBER_OF_ROC;iroc++) {
evt_count_valid[iroc] = 0;
evt_count_bad[iroc] = 0;
}
Int_t idum;
while (coda->codaRead()==0 && ((evt_count_valid[0]+evt_count_bad[0]) < event_preset)) {
buffer = coda->getEvBuffer();
idum = ev->THaEvData::LoadEvent(buffer);
//if (ev->getEvType()==1)
// int nevt=ev->GetEvLength();
// for (int i = 0;i<nevt;i++) {
for (iroc = 0; iroc< NUMBER_OF_ROC; iroc++) {
if (iroc == 0) {
roc = 14;
channel_offset = 0;
} else {
roc = 15;
channel_offset = CHANNEL_X_ADC;
}
if ((ev->GetRawData(roc,EVT_HEADER_OFFSET) & EVT_HEADER_TAG) == EVT_HEADER_TAG) {
// memset(ch,0,total_channel*sizeof(Short_t)); // clear chv
evt_count_valid[iroc]++;
// evt = evt_count_valid[iroc];
evt_count = ev->GetRawData(roc, EVT_COUNT_OFFSET) & EVT_COUNT_MASK;
evt = evt_count;
// cout << "EVENT COUNT = " << evt_count << " " << evt << endl;
Int_t idx;
idx = ADC_HEADER_OFFSET;
unsigned int adctag;
adc_count = 0;
nch = 0;
tcharge = 0;
Int_t cch;
cch = 0;
for (Int_t iadc = 0; iadc <NUMBER_ADC;iadc++) { // loop on adc
adctag = ev->GetRawData(roc,idx);
if ((adctag & ADC_HEADER_TAG) == ADC_HEADER_TAG) {
adc_count++;
ch_count = adctag & CH_COUNT_MASK;
nch += ch_count;
idx++;
for (Int_t js = 0; js<ch_count;js++) { // loop on channels
Int_t duv = ev->GetRawData(roc,idx);
if ((duv & ADC_VALID_BIT) ) {
adc_ch_value = (Short_t) (duv & ADC_VALUE_MASK);
adc_ch = (Int_t) (duv & ADC_CH_MASK) >> 12;
adc_ch_flag = (Short_t) (duv & (ADC_VALID_BIT | ADC_OVER_BIT)) >> 30;
Int_t ich =adc_ch+iadc*CHANNEL_X_ADC*2+channel_offset; // assume interleaved ROC14/ROC15 ADC
charge[cch]=(Int_t) adc_ch_value;
ch[cch]=(Int_t) ich;
tcharge += (Int_t) adc_ch_value;
cch++;
pdata->Fill(ich,adc_ch_value);
if (mode != 0) {
mdata->Fill(ch2x(ich),ch2y(ich),adc_ch_value);
//mdata->Fill(ch2x(ich),ch2y(ich),1);
}
eff_max_ch = (ich > eff_max_ch) ? ich : eff_max_ch;
eff_min_ch = (ich < eff_min_ch) ? ich : eff_min_ch;
idx++;
} else {
cout << adc_count << " " << js << " not valid/overrange " << hex << duv << endl;
}
}
} else {
if ((adctag & ADC_HEADER_NUL) == ADC_HEADER_NUL) {
// cout << "No event/adc " << evt_count << " / " << iadc
// << ", HEADER= 0x"<< hex << adctag << endl;
idx += 2;
}
}
}
if (mode != 0) tdata->Fill();
if (((evt_count_valid[iroc]+evt_count_bad[iroc]) % 100) == 0) {
cout << "Roc " << roc << " : events so far (valid/bad): " << dec << evt_count_valid[iroc] << " / " << evt_count_bad[iroc] << endl;
}
} else {
cout << "WARNING: event " << (evt_count+evt_count_bad[iroc]) << " with different header = " << hex << ev->GetRawData(roc,EVT_HEADER_OFFSET) << endl;
evt_count_bad[iroc]++;
}
} // roc loop
} // while on buffer
} // coda
//
cout << " Total channels : " << nch << endl;
for (iroc=0;iroc<NUMBER_OF_ROC;iroc++) {
cout << iroc << " Total valid/bad events processed: " << dec << evt_count_valid[iroc] << " / " << evt_count_bad[iroc] << endl;
if (eff_min_ch > eff_max_ch) {
cout << " WARNING: No RICH data in events !!! " << endl;
return;
}
if (evt_count_valid[iroc] == 0) {
return;
}
}
cout << " Effective Channel range from " << dec << eff_min_ch << " to " << eff_max_ch << endl;
//
//
// analysis
//
//
Int_t iea;
Int_t eff_range_ch = eff_max_ch - eff_min_ch + 1;
Int_t bad_ch = 0;
Int_t pbad_ch = 0;
if (eff_range_ch < total_channel) {
cout << "WARNING: Effective detected channels : " << eff_range_ch << " expected " << total_channel << endl;
}
if (mode == 0) { // pedestals
Double_t amean = pdata->GetMean(2);
Double_t arms = pdata->GetRMS(2);
TArrayD mean = TArrayD(eff_range_ch);
TArrayD rms = TArrayD(eff_range_ch);
TArrayD tadx = TArrayD(eff_range_ch);
cout << " Average mean and RMS = " << amean << " " << arms << endl;
for (ie = 0; ie<eff_range_ch; ie++) {
iea = ie + eff_min_ch;
tadx.AddAt((Double_t) iea, ie);
rms.SetAt( pdata->GetBinError(iea + 1), ie);
mean.SetAt( pdata->GetBinContent(iea + 1), ie);
if (mode == 0) { // pedestal mode
if ((mean.GetAt(ie) == 0) || (rms.GetAt(ie) == 0)) {
cout << " Ch " << iea << " in ADC " << ch2ADC(iea)
<< " (panel " << ch2panel(iea) << " row " << ch2prow(iea) << " card " << ch2rcard(iea)
<< " ) looks bad: mean= "
<< mean.GetAt(ie) << " , rms= " << rms[ie] << endl;
bad_ch++;
}
if (TMath::Abs(mean.GetAt(ie)-amean) > (ACCEPTED_NSIGMA * arms)) {
cout << " Ch " << iea << " in ADC " << ch2ADC(iea)
<< " (panel " << ch2panel(iea) << " row " << ch2prow(iea) << " card " << ch2rcard(iea)
<< " ) maybe bad: mean= "
<< mean.GetAt(ie) << " , rms= " << rms.GetAt(ie) << endl;
pbad_ch++;
}
}
}
cout << "Total Bad Channels : " << bad_ch << " ( " << ((Double_t) bad_ch)/eff_range_ch << " % ) " << endl;
cout << "Total maybe Bad Chs: " << pbad_ch << " ( " << ((Double_t) pbad_ch)/eff_range_ch << " % ) " << endl;
grmean = new TGraph(mean.GetSize(), tadx.GetArray(), mean.GetArray() );
grrms = new TGraph(rms.GetSize(), tadx.GetArray(), rms.GetArray() );
grscat = new TGraph(rms.GetSize(), mean.GetArray(), rms.GetArray() );
//
// Plot data
//
TCanvas *cdm = new TCanvas("cdm","Mean Distribution per ADC");
cdm->Divide(ADC_X_PANEL,NUMBER_OF_ROC*NUMBER_ADC/ADC_X_PANEL);
cdm->Update();
TCanvas *cdr = new TCanvas("cdr","RMS Distribution per ADC");
cdr->Divide(ADC_X_PANEL,NUMBER_OF_ROC*NUMBER_ADC/ADC_X_PANEL);
cdr->Update();
TCanvas *cdp = new TCanvas("cdp","Profile per ADC");
cdp->Divide(ADC_X_PANEL,NUMBER_OF_ROC*NUMBER_ADC/ADC_X_PANEL);
cdp->Update();
TH1D *hdm;
TH1D *hdr;
Double_t hmin, hmax;
hmin = amean - 1.2*ACCEPTED_NSIGMA*arms;
hmax = amean + 1.2*ACCEPTED_NSIGMA*arms;
hdm = new TH1D("hdm","Mean of ADC",100,hmin,hmax);
hdr = new TH1D("hdr","RMS of ADC",100,0,arms*1.2); // TBC
Int_t je_off=0;
// je_off = (roc == 14) ? 0 : 8;
for (Int_t ie = 0; ie<eff_range_ch; ie+=CHANNEL_X_ADC) {
iea = ie + eff_min_ch;
Int_t dmax = eff_range_ch > (ie + CHANNEL_X_ADC) ? (ie + CHANNEL_X_ADC) : eff_range_ch;
hdm->Reset();
for (Int_t je = ie; je < dmax + je_off; je++) {
hdm->Fill(mean.GetAt(je));
hdr->Fill(rms.GetAt(je));
}
hdr->SetLineColor(ch2roc(iea)+1);
hdm->SetLineColor(ch2roc(iea)+1);
pdata->SetLineColor(ch2roc(iea)+1);
pdata->SetMarkerColor(ch2roc(iea)+1);
cdm->cd(ch2arow(iea)+1);
hdm->SetTitle(Form("ADC%2d p%1d r%1d c%1d",
ch2ADC(iea), ch2panel(iea), ch2prow(iea), ch2rcard(iea)));
hdm->DrawCopy();
cdm->Update();
cdr->cd(ch2arow(iea)+1);
hdr->SetTitle(Form("ADC%2d p%1d r%1d c%1d",
ch2ADC(iea), ch2panel(iea), ch2prow(iea), ch2rcard(iea)));
hdr->DrawCopy();
cdr->Update();
cdp->cd(ch2arow(iea)+1);
pdata->SetTitle(Form("ADC%2d p%1d r%1d c%1d",
ch2ADC(iea), ch2panel(iea), ch2prow(iea), ch2rcard(iea)));
pdata->SetAxisRange(iea,iea+CHANNEL_X_ADC-1);
pdata->DrawCopy();
cdp->Update();
}
//cdm->SaveAs(Form("%s_mean.jpg",filename));
//cdr->SaveAs(Form("%s_rms.jpg",filename));
//cdp->SaveAs(Form("%s_profile.jpg",filename));
TCanvas *cgr = new TCanvas("cgr","Pedestal and RMS");
cgr->Divide(2,2);
cgr->Update();
cgr->cd(1);
grmean->Draw("PA");
grmean->SetTitle("Pedestal");
cgr->Update();
cgr->cd(2);
grrms->Draw("PA");
grrms->SetTitle("Noise");
cgr->Update();
cgr->cd(3);
grscat->Draw("PA");
grscat->SetTitle("Noise vs Pedestal");
cgr->Update();
pdata->SetMarkerColor(1);
pdata->SetLineColor(1);
cgr->cd(4);
pdata->SetTitle("Pedestal with Noise");
pdata->SetAxisRange(eff_min_ch-0.5,eff_max_ch-0.5);
pdata->Draw();
cgr->Update();
/*
TCanvas *cpedt = new TCanvas("cpedt","Pedestal / Noise");
cpedt->Divide(1,1);
cpedt->Update();
pdata->SetTitle("Pedestal / Noise");
pdata->SetAxisRange(eff_min_ch-0.5,eff_max_ch-0.5);
cpedt->cd(1);
pdata->Draw();
cpedt->Update();
*/
//cpedt->SaveAs(Form("%s_total.pdf",filename));
//cpedt->SaveAs(Form("%s_total.jpg",filename));
}
//
// Output on file
FILE *fout[NUMBER_OF_ROC];
fout[0] = fopen(Form("rich_ped_bottom.dat",filename),"w");
fout[1] = fopen(Form("rich_ped_top.dat",filename),"w");
for (iroc=0;iroc<NUMBER_OF_ROC;iroc++) {
fprintf(fout[iroc],"# ADC ch pedestal rms (ADC unit, ROC%d , run %d , sigma %f )\n",iroc+14, nrun, nsigma);
}
for (ie = 0; ie<eff_range_ch; ie++) {
iea = ie + eff_min_ch;
iroc = ch2roc(iea);
fprintf(fout[iroc],"%3d %4d %4d %4d\n",ch2ADC(iea), ch2rch(iea), (Int_t) mean.GetAt(ie), (Int_t) (mean.GetAt(ie)+nsigma*rms.GetAt(ie)));
}
fclose(fout[iroc]);
// plotting
// Whole pedestals and noise
if (mode == 1) { // standard data
// map
TCanvas *cmap = new TCanvas("cmap","RICH HIT MAP");
cmap->Divide(1,1);
cmap->Update();
cmap->cd(1);
mdata->Draw("zcol");
cmap->Update();
//cmap->SaveAs(Form("%s_map.pdf",filename));
//cmap->SaveAs(Form("%s_map.jpg",filename));
TCanvas *cmap1 = new TCanvas("cmap1","RICH HIT MAP Projections");
cmap1->Divide(1,2);
cmap1->Update();
cmap1->cd(1);
mdata->ProjectionX()->Draw();
cmap1->cd(2);
mdata->ProjectionY()->Draw();
//cmap1->SaveAs(Form("%s_map_proj.pdf",filename));
//cmap1->SaveAs(Form("%s_map_proj.jpg",filename));
TCanvas *cpedt = new TCanvas("cpedt","Data");
cpedt->Divide(2,2);
cpedt->Update();
(cpedt->cd(1))->SetLogy();
tdata->Draw("nch>>hnch","nch<100");
TH1 *hnch = (TH1*) gROOT->FindObject("hnch");
hnch->Fit("gaus");
cpedt->Update();
cpedt->cd(2);
tdata->Draw("nch:evt");
cpedt->Update();
(cpedt->cd(3))->SetLogy();
tdata->Draw("charge>>hcharge");
TH1 *hcharge = (TH1*) gROOT->FindObject("hcharge");
hcharge->Fit("landau");
cpedt->Update();
(cpedt->cd(4))->SetLogy();
tdata->Draw("tcharge>>htcharge","tcharge<6000");
TH1 *htcharge = (TH1*) gROOT->FindObject("htcharge");
htcharge->Fit("landau");
cpedt->Update();
// cpedt->SaveAs(Form("%s_stat1.pdf",filename));
// cpedt->SaveAs(Form("%s_stat1.jpg",filename));
TCanvas *ctree = new TCanvas("ctree","Data Selected");
ctree->Divide(2,2);
ctree->Update();
ctree->cd(1);
tdata->Draw("charge:ch");
ctree->Update();
ctree->cd(2);
tdata->Draw("charge:nch");
ctree->Update();
ctree->cd(3);
tdata->Draw("charge:evt");
ctree->Update();
ctree->cd(4);
tdata->Draw("tcharge:nch");
ctree->Update();
// ctree->SaveAs(Form("%s_stat2.pdf",filename));
// ctree->SaveAs(Form("%s_stat2.jpg",filename));
}
//
}
