10 #include <arich/modules/arichBtest/arichBtestModule.h>
17 #include <framework/datastore/StoreObjPtr.h>
18 #include <framework/dataobjects/EventMetaData.h>
19 #include <framework/datastore/StoreArray.h>
20 #include <framework/logging/Logger.h>
21 #include <framework/gearbox/Unit.h>
24 #include <boost/algorithm/string.hpp>
25 #include <boost/foreach.hpp>
29 #include <arich/dataobjects/ARICHDigit.h>
30 #include <arich/dataobjects/ARICHAeroHit.h>
31 #include "arich/geometry/ARICHGeometryPar.h"
32 #include "arich/geometry/ARICHBtestGeometryPar.h"
34 #include "arich/modules/arichBtest/arichBtestData.h"
46 const char* record_strings[] = {
"Event",
"Begin_run",
"Pause",
"Resume",
"End_run" };
70 setDescription(
"Module for the ARICH Beamtest data analysis. It creates track form the MWPC hits and reads the HAPD hits");
73 addParam(
"outputFileName", m_outFile,
"Output Root Filename",
string(
"output.root"));
74 vector<string> defaultList;
75 addParam(
"runList", m_runList,
"Data Filenames.", defaultList);
76 vector<int> defaultMask;
77 addParam(
"mwpcTrackMask", m_MwpcTrackMask,
"Create track from MWPC layers", defaultMask);
79 addParam(
"beamMomentum", m_beamMomentum,
"Momentum of the beam [GeV]", 0.0);
81 for (
int i = 0; i < 32; i++) m_tdc[i] = 0;
97 void arichBtestModule::initialize()
99 B2INFO(
"arichBtestModule::initialize()");
102 aeroHits.registerInDataStore();
103 digits.registerInDataStore();
105 m_file =
new TFile(m_outFile.c_str(),
"RECREATE");
106 m_tuple =
new TNtuple(
"nt",
"Btest data hits",
"x:y:xrec:yrec:m:c:gx:gy");
108 for (
int i = 0; i < 6; i++) {
109 sprintf(hapdName,
"hapd%d", i);
110 hapd[i] =
new TH1F(hapdName, hapdName, 145, -0.5, 144.5);
114 for (
int i = 0; i < 4; i++) {
115 for (
int k = 0; k < 2; k++) {
116 sprintf(name,
"mwpc%d_a%d_", i, k);
117 mwpc_diff[i][k] =
new TH1F(strcat(name,
"diff"), name, 300, -150, 150);
118 sprintf(name,
"mwpc%d_a%d_", i, k);
119 mwpc_sum[i][k] =
new TH1F(strcat(name,
"sum"), name, 200, -0.5, 199.5);
120 sprintf(name,
"mwpc%d_a%d_", i, k);
121 mwpc_sum_cut[i][k] =
new TH1F(strcat(name,
"sum_cut"), name, 200, -0.5, 199.5);
122 sprintf(name,
"mwpc%d_a%d_", i, k);
123 mwpc_residuals[i][k] =
new TH1F(strcat(name,
"resd"), name, 200, -100, 100);
124 sprintf(name,
"mwpc%d_a%d_", i, k);
125 mwpc_residualsz[i][k] =
new TH2F(strcat(name,
"resd_z"), name, 200, -25, 25, 400, -1000, 1000);
127 for (
int k = 0; k < 5; k++) {
128 sprintf(name,
"mwpc%d_a%d_", i, k);
129 mwpc_tdc[i][k] =
new TH1F(strcat(name,
"tdc"), name, 1024, -1024, 1024);
131 sprintf(name,
"mwpc%d_", i);
132 mwpc_xy[i] =
new TH2F(strcat(name,
"xy"), name, 120, -30, 30, 120, -30, 30);
153 void arichBtestModule::beginRun()
156 B2INFO(
"arichBtestModule::beginRun()");
159 B2INFO(
"arichBtestModule::eventMetaDataPtr run:" << eventMetaDataPtr->getRun());
160 B2INFO(
"arichBtestModule::eventMetaDataPtr exp:" << eventMetaDataPtr->getExperiment());
163 m_mwpc = _arichbtgp->
getMwpc();
165 static int first = 1;
167 m_runCurrent = m_runList.begin();
177 if (m_runCurrent < m_runList.end()) {
179 m_eveList.push_back(m_events);
182 m_fp = gzopen((*m_runCurrent).c_str(),
"rb");
184 B2INFO(
"Cannot open " << *m_runCurrent);
187 B2INFO(
"File opened " << *m_runCurrent);
194 int arichBtestModule::skipdata(gzFile fp)
197 gzread(fp, &u,
sizeof(
unsigned int));
198 gzseek(fp, u *
sizeof(
unsigned int), SEEK_CUR);
202 void arichBtestModule::readmwpc(
unsigned int* dbuf,
unsigned int len,
int print1290)
205 const int unsigned MAXV1290 = 32;
212 for (
int i = 0; i < 32; i++) m_tdc[i] = 0XFFFF;
215 for (
unsigned int i = 0; i < len; i++) {
216 int rid = (dbuf[i] >> 27) & 0x1F;
219 if (print1290) printf(
"Filler 0x%08x %u.data\n", dbuf[i], i);
222 if (print1290) printf(
"Global Header 0x%08x %u.data\n", dbuf[i], i);
225 nhits = ((dbuf[i] >> 5) & 0xFFFF);
226 if (print1290) printf(
"Global Trailer 0x%08x %u.data STATUS=0x%03x nhits=%u\n", dbuf[i], i, (dbuf[i] >> 24) & 0x7, nhits);
229 if (print1290) printf(
"V1290 nhits!=len %u %u\n", nhits, len);
233 if (print1290) printf(
"Global Trigger TimeTag 0x%08x %u.data\n", dbuf[i], i);
236 if (print1290) printf(
"TDC header 0x%08x %u.data evid=%d wc=%u\n", dbuf[i], i, (dbuf[i] >> 12) & 0xFFF, dbuf[i] & 0xFFF);
240 if (print1290) printf(
"TDC trailer 0x%08x %u.data\n", dbuf[i], i);
244 if (print1290) printf(
"TDC Error 0x%08x %u.data ERR=0x%x\n", dbuf[i], i, dbuf[i] & 0x3FFF);
249 edge = (dbuf[i] >> 26) & 0x1;
250 tdcch = (dbuf[i] >> 21) & 0x1F;
251 tdcl = dbuf[i] & 0x1FFFFF;
252 if (tdcch < MAXV1290) {
256 printf(
"V1290 0x%08x %u. [ch=%2u] edge=%u data=%u \n", dbuf[i], i, tdcch, edge, tdcl);
257 m_tdc[tdcch] = tdcl / 40;
263 if (print1290) printf(
"Unknown 0x%08x %u.data\n", dbuf[i], i);
272 int arichBtestModule::readhapd(
unsigned int len,
unsigned int* data)
281 for (
int module = 0; module < 6; module++) {
282 int istart = 19 * module;
283 int istop = 19 * module + 18;
284 for (
int i = istart; i < istop; i++) {
285 for (
int j = 0; j < 8; j++) {
286 unsigned int kdata = data[i] & (bmask << (j * 4));
288 int channelID = j + 8 * (i - istart);
290 if (_arichgp->
isActive(module, channelID)) {
292 hapd[module]->Fill(channelID);
295 double globalTime = 0;
297 double rposx = 0, rposy = 0;
299 int channel = eposhapd.second;
301 if ((channel < 108 && channel > 71) || channel < 36) channel = 108 - (int(channel / 6) * 2 + 1) * 6 +
303 else channel = 144 - (int((channel - 36) / 6) * 2 + 1) * 6 + channel - 36;
305 if (abs(loc.X()) > 2.3 || abs(loc.Y()) > 2.3)
continue;
307 arichDigits.
appendNew(module + 1, channel, globalTime);
311 m_tuple ->Fill(-poshapd.first, poshapd.second, rechit.x(), rechit.y(), module, channelID, rposx, rposy);
322 int arichBtestModule::getTrack(
int mask, TVector3& r, TVector3& dir)
327 for (
int i = 0; i < 4; i++) {
330 for (
int k = 0; k < 4; k++)
mwpc_tdc[i][k]->Fill(m_tdc[w->tdc[k]] - t0);
331 mwpc_tdc[i][4]->Fill(m_tdc[w->atdc] - t0);
334 for (
int k = 0; k < 2; k++) {
336 w->diff[k] = m_tdc[w->tdc[2 * k]] - m_tdc[w->tdc[2 * k + 1]];
337 w->sum[k] = m_tdc[w->tdc[2 * k + 1]] + m_tdc[w->tdc[2 * k]] - 2 * m_tdc[w->atdc];
339 if (w->sum[k] < w->cutll[k] || w->sum[k] > w->cutul[k])
continue;
342 w->reco[k] = w->diff[k] * w->slp[k] + w->offset[k];
343 w->reco[k] += w->pos[k];
348 w->reco[2] = w->pos[2];
349 if (!w->status[0] && !w->status[1])
mwpc_xy[i]->Fill(w->reco[0], w->reco[1]);
351 if (mask & (1 << i)) {
352 if (!w->status[0] && !w->status[1]) {
361 int ii[4] = {0, 1, 0, 0};
363 for (
int i = 0; i < 4; i++) {
364 if (mask & (1 << i)) {
372 r.SetXYZ(m_mwpc[i0].reco[1], m_mwpc[i0].reco[0], m_mwpc[i0].reco[2]);
373 dir.SetXYZ(m_mwpc[i1].reco[1] - m_mwpc[i0].reco[1],
374 m_mwpc[i1].reco[0] - m_mwpc[i0].reco[0],
375 m_mwpc[i1].reco[2] - m_mwpc[i0].reco[2]);
381 for (
int i = 0; i < 4; i++) {
383 double l = (w->reco[2] - r.z()) / dir.z() ;
384 TVector3 rext = r + dir * l;
385 if (!w->status[0])
mwpc_residuals[i][0]->Fill(w->reco[0] - rext.y());
386 if (!w->status[1])
mwpc_residuals[i][1]->Fill(w->reco[1] - rext.x());
389 for (
int k = 0; k < axis->GetNbins(); k++) {
390 double ll = (w->reco[2] + axis->GetBinCenter(k + 1) - r.z()) / dir.z();
391 TVector3 rextt = r + dir * ll;
392 mwpc_residualsz[i][0]->Fill(w->reco[0] - rextt.y(), axis->GetBinCenter(k + 1));
393 mwpc_residualsz[i][1]->Fill(w->reco[1] - rextt.x(), axis->GetBinCenter(k + 1));
402 int arichBtestModule::readdata(gzFile fp,
int rec_id,
int)
405 unsigned int len, data[10000];
406 gzread(fp, &len,
sizeof(
unsigned int));
408 gzread(fp, data,
sizeof(
unsigned int)*len);
414 int retval = getTrack(*(m_MwpcTrackMask.begin()), r, dir);
422 dir *= m_beamMomentum * Unit::GeV;
435 TVector3 rrel = rc - rot * rc;
438 r.SetX(-r.X()); dir.SetX(-dir.X());
445 B2DEBUG(50,
"-----------> " << rc.x() <<
" " << rc.y() <<
" " << rc.z() <<
"::::" << rrel.x() <<
" " << rrel.y() <<
" " <<
446 rrel.z() <<
" ----> R " << r.x() <<
" " << r.y() <<
" " << r.z() <<
" ----> S " << dir.x() <<
" " << dir.y() <<
" " <<
450 arichAeroHits.
appendNew(particleId, r, dir);
459 for (
unsigned int j = 0; j < len; j++) {
470 void arichBtestModule::event()
482 static char msg[1024];
485 const int sint =
sizeof(
unsigned int);
487 if (gzread(m_fp, hdr, sint * 2) != 2 * sint && m_end) {
488 B2INFO(
"[" << m_events <<
"] End of file: " << *m_runCurrent);
490 if (m_runCurrent == m_runList.end()) {
492 eventMetaDataPtr->setEndOfData();
498 }
while (m_end && m_runCurrent != m_runList.end());
499 if (m_events % 1000 == 0) {
502 B2INFO(
"neve= [" << m_events <<
"] in " << (
double)(m_time - m_timestart) / 60. <<
" min (" <<
int(
503 m_time - m_timestart) <<
"s) from " << *m_runCurrent);
511 gzseek(m_fp, 2 *
sizeof(
unsigned int), SEEK_CUR);
513 case BEGIN_RECORD_TYPE: {
514 gzread(m_fp, &beginrec,
sizeof(beginrec));
515 time_t t = beginrec.
time;
516 sprintf(msg,
"BeginRec run %u time %s", beginrec.
runno, ctime(&t));
520 case END_RECORD_TYPE: {
521 gzread(m_fp, &endrec,
sizeof(endrec));
522 time_t t = endrec.
time;
523 sprintf(msg,
"EndRec run %u time %s", endrec.
runno, ctime(&t));
527 case EVENT_RECORD_TYPE: {
528 gzread(m_fp, &rec,
sizeof(rec));
529 int print = !(rec.evtno % 10000);
532 sprintf(msg,
"EventRec run %u evt %u mstime %u, time %s", rec.runno, rec.evtno, rec.mstime, ctime(&t));
536 int pos = gztell(m_fp);
538 int buflen = rec.len -
sizeof(rec) /
sizeof(
int);
543 while (i < buflen && j < 5) {
544 n[j] = readdata(m_fp, j, print);
549 if (gzseek(m_fp, pos + buflen *
sizeof(
unsigned int), SEEK_SET) < 0) B2ERROR(
"gzseek returns -1 ");
553 B2ERROR(
"IO error unknown record type " << type);
556 if (gzeof(m_fp)) m_end = 1;
559 void arichBtestModule::endRun()
561 B2INFO(
" arichBtestModule: End Run !!!");
568 m_eveList.push_back(m_events);
572 void arichBtestModule::terminate()
575 BOOST_FOREACH(
const string & fname, m_runList) {
576 B2INFO(m_eveList[j] <<
" events processed from file " << fname);
579 for (
int i = 0; i < 4; i++) {
581 B2INFO(i <<
" a1=" << m_mwpc[i].tdc[0] <<
" a2=" << m_mwpc[i].tdc[1] <<
" a3=" << m_mwpc[i].tdc[2] <<
" a2=" << m_mwpc[i].tdc[3]
582 <<
" A=" << m_mwpc[i].atdc);