Newer
Older
TB_Chris / TbUT / src / alibava / .svn / text-base / TbAsciiRoot.cpp.svn-base
  1. #include <iostream>
  2. #include <fstream>
  3. #include <iomanip>
  4. #include <ctime>
  5. #include <cmath>
  6. #include <sstream>
  7. #include <csignal>
  8. #include <vector>
  9. #include <cstdlib>
  10. #include <cerrno>
  11. #include <climits>
  12. #include <TROOT.h>
  13. #include <TCanvas.h>
  14. #include <TProfile2D.h>
  15. #include <TF1.h>
  16. #include <TH1.h>
  17. #include <TH2.h>
  18. #include "TbAsciiRoot.h"
  19. #include "utils.h"
  20. #include "Tracer.h"
  21. #include "TbAlibavaHit.h"
  22. #include "ChanList.h"
  23. /*
  24.  
  25. #ifdef __APPLE__
  26. #define sighandler_t sig_t
  27. #endif
  28.  
  29. bool _A_do_run = true;
  30. void _A_got_intr(int)
  31. {
  32. _A_do_run = false;
  33. }
  34. */
  35. // decodes the header and returns a vector with the integers found
  36. std::vector<int> decode_header(const std::string &h, TbAsciiRoot::XtraValues &xtra)
  37. {
  38. std::vector<int> vout;
  39. std::istringstream istr(h);
  40. char *endptr;
  41. char buf[256];
  42. long val;
  43.  
  44. xtra.clear();
  45.  
  46. while (istr)
  47. {
  48. istr.getline(buf, sizeof(buf), ';');
  49. if (!istr)
  50. break;
  51.  
  52. errno = 0;
  53. val = strtol(buf, &endptr, 0);
  54.  
  55. if ((errno == ERANGE && (val == LONG_MAX || val == LONG_MIN)) || (errno != 0 && val == 0))
  56. {
  57. std::string sval(buf), sout;
  58. sout = trim_str(sval);
  59. if (!sout.empty())
  60. xtra.push_back( sout );
  61. }
  62. else if ( endptr == buf || *endptr != '\0' )
  63. {
  64. std::string sval(buf), sout;
  65. sout = trim_str(sval);
  66. if (!sout.empty())
  67. xtra.push_back( sout );
  68. }
  69. else
  70. {
  71. vout.push_back(atoi(buf) );
  72. }
  73. }
  74. return vout;
  75. }
  76.  
  77. TbAsciiRoot::TbAsciiRoot(const char * /*nam*/, const char * /*pedfile*/, const char * /*gainfile*/) :
  78. _nchan(max_nchan),_seedcut(10.), _neighcut(5.), _average_gain(1.), _version(2), _polarity(1)
  79. {
  80. int i;
  81. for (i=0;i<max_nchan;i++)
  82. {
  83. _ped[i] = 0.;
  84. _gain[i] = 1.;
  85. _noise[i] = 1.;
  86. _data.data[i] = 0;
  87. _mask[i] = false;
  88. }
  89. ifile=NULL;
  90. }
  91.  
  92. TbAsciiRoot::~TbAsciiRoot()
  93. {
  94. if (ifile)
  95. {
  96. ifile->close();
  97. }
  98. }
  99.  
  100. void TbAsciiRoot::open(const char *name)
  101. {
  102. ifile = new std::ifstream(name);
  103. if (!(*ifile))
  104. {
  105. std::cout << "Could not open data file: " << name << std::endl;
  106. delete ifile;
  107. ifile = 0;
  108. return;
  109. }
  110. std::string header;
  111. unsigned int ic, lheader;
  112. char c;
  113. ifile->read((char *)&_t0, sizeof(time_t));
  114. //std::cout << "64b: " << ctime(&_t0) << std::endl;
  115. ifile->read((char *)&_type, sizeof(int));
  116. //std::cout << "type_ " << _type << std::endl;
  117. if ( _type > 15 )
  118. {
  119. ifile->seekg(0, std::ios::beg);
  120. ifile->read((char *)&_t0, sizeof(int));
  121. ifile->read((char *)&_type, sizeof(int));
  122. //std::cout << "32b: " << ctime(&_t0) << std::endl;
  123. }
  124.  
  125.  
  126. ifile->read((char *)&lheader, sizeof(unsigned int));
  127. for (ic=0; ic<80; ic++)
  128. {
  129. ifile->read(&c, sizeof(char));
  130. header.append(1, c);
  131. }
  132. header = trim_str(header);
  133.  
  134. if (header[0]!='V' && header[0]!='v')
  135. {
  136. _version = 0;
  137. }
  138. else
  139. {
  140. _version = int(header[1]-'0');
  141. header = header.substr(5);
  142. }
  143.  
  144. std::cout << "type: " << _type << " header: " << header << std::endl;
  145. std::vector<int> param = decode_header(header, _xtra);
  146. ifile->read((char *)_ped, max_nchan*sizeof(double));
  147. ifile->read((char *)_noise, max_nchan*sizeof(double));
  148. switch (_type)
  149. {
  150. case 1: // calibration
  151. case 2: // laser sync
  152. _npoints = param[0];
  153. _from = param[1];
  154. _to = param[2];
  155. _step = param[3];
  156. break;
  157. case 3: // laser run
  158. case 4: // source run
  159. case 5: // pedestal run
  160. if (param.empty())
  161. _nevts = 100000;
  162. else
  163. _nevts = param[0];
  164. _npoints = _from = _to = _step = 0;
  165. break;
  166. }
  167. data_start = ifile->tellg();
  168. }
  169.  
  170. void TbAsciiRoot::rewind()
  171. {
  172. if (ifile)
  173. {
  174. ifile->clear();
  175. ifile->seekg(data_start, std::ios::beg);
  176. }
  177. }
  178.  
  179. void TbAsciiRoot::close()
  180. {
  181. if (ifile)
  182. {
  183. ifile->close();
  184. delete ifile;
  185. ifile = 0;
  186. }
  187. }
  188.  
  189. void TbAsciiRoot::reset_data()
  190. {
  191. memset(&_data, 0, sizeof(_data));
  192. }
  193.  
  194. int TbAsciiRoot::read_event(std::string & error_code)
  195. {
  196. if (ifile)
  197. {
  198. unsigned int header, size, user=0, code=0;
  199. char *block_data=0;
  200. if (_version)
  201. {
  202. do
  203. {
  204. do
  205. {
  206. ifile->read((char *)&header, sizeof(unsigned int));
  207. if (ifile->bad() || ifile->eof()){
  208. error_code="read Header";
  209. return -1;
  210. }
  211.  
  212. code = (header>>16) & 0xFFFF;
  213. } while ( code != 0xcafe );
  214.  
  215. code = header & 0x0fff;
  216. user = header & 0x1000;
  217. switch (code)
  218. {
  219. case NewFile:
  220. ifile->read((char *)&size, sizeof(unsigned int));
  221. block_data = new char[size];
  222. ifile->read(block_data, size);
  223. new_file(size, block_data);
  224. break;
  225. case StartOfRun:
  226. ifile->read((char *)&size, sizeof(unsigned int));
  227. block_data = new char[size];
  228. ifile->read(block_data, size);
  229. start_of_run(size, block_data);
  230. break;
  231. case DataBlock:
  232. ifile->read((char *)&size, sizeof(unsigned int));
  233. if (user)
  234. {
  235. reset_data();
  236. block_data = new char[size];
  237. ifile->read(block_data, size);
  238. new_data_block(size, block_data);
  239. }
  240. else
  241. {
  242. if ( _version == 1 )
  243. {
  244. ifile->read((char *)&_data, sizeof(EventData));
  245. for (int ii=0; ii<2; ii++)
  246. memset(_header[ii], 0, 16*sizeof(unsigned short));
  247.  
  248. }
  249. else
  250. {
  251. ifile->read((char *)&_data.value, sizeof(double));
  252. ifile->read((char *)&_data.time, sizeof(unsigned int));
  253. ifile->read((char *)&_data.temp, sizeof(unsigned short));
  254. for (int ii=0; ii<2; ii++)
  255. {
  256. ifile->read((char *)_header[ii], 16*sizeof(unsigned short));
  257. ifile->read((char *)&_data.data[ii*128], 128*sizeof(unsigned short));
  258. }
  259. }
  260. }
  261.  
  262. break;
  263. case CheckPoint:
  264. ifile->read((char *)&size, sizeof(unsigned int));
  265. block_data = new char[size];
  266. ifile->read(block_data, size);
  267. check_point(size, block_data);
  268. break;
  269. case EndOfRun:
  270. ifile->read((char *)&size, sizeof(unsigned int));
  271. block_data = new char[size];
  272. ifile->read(block_data, size);
  273. end_of_run(size, block_data);
  274. break;
  275. default:
  276. std::cout << "Unknown block data type: " << std::hex << header << " - " << code << std::dec << std::endl;
  277. }
  278. if (block_data)
  279. {
  280. delete [] block_data;
  281. block_data = 0;
  282. }
  283.  
  284. } while ( code != DataBlock && !(ifile->bad() || ifile->eof()) );
  285. }
  286. else
  287. {
  288. ifile->read((char *)&_data, sizeof(EventData));
  289. for (int ii=0; ii<2; ii++)
  290. memset(_header[ii], 0, 16*sizeof(unsigned short));
  291. }
  292.  
  293. if (ifile->eof())
  294. {
  295. std::cout << "End of file" << std::endl;
  296. error_code="End of file";
  297. return -1;
  298. }
  299. else if (ifile->bad())
  300. {
  301. std::cout << "Problems with data file" << std::endl;
  302. error_code="Problems with data file";
  303. return -1;
  304. }
  305. else
  306. {
  307. error_code="No problem";
  308. //process_event();
  309. return 0;
  310. }
  311. }
  312. else{
  313. error_code="big file problem";
  314. return -1;
  315. }
  316. }
  317.  
  318. void TbAsciiRoot::set_data(int nchan, const unsigned short int *data)
  319. {
  320. int i;
  321. _nchan = nchan;
  322. for (i=0;i<_nchan;i++)
  323. _data.data[i] = data[i];
  324. }
  325.  
  326.  
  327. double TbAsciiRoot::time() const
  328. {
  329. unsigned short fpart = _data.time & 0xffff;
  330. short ipart = (_data.time & 0xffff0000)>>16;
  331. if (ipart<0)
  332. fpart *= -1;
  333. //double tt = 100.*(1. -(ipart + (fpart/65535.)));
  334. double tt = 100.0*(ipart + (fpart/65535.));
  335. return tt;
  336. }
  337.  
  338. double TbAsciiRoot::temp() const
  339. {
  340. return 0.12*_data.temp - 39.8;
  341. }
  342.  
  343.  
  344. /*
  345. TH1 *TbAsciiRoot::show_pedestals()
  346. {
  347. int ic;
  348. TH1 *hst = create_h1("hPed","Pedestals",nchan(),-0.5, nchan()-0.5);
  349. hst->SetYTitle("ADCs");
  350. hst->SetXTitle("Channel no.");
  351. for (ic=0; ic<nchan(); ic++)
  352. hst->SetBinContent(ic+1, _ped[ic]);
  353.  
  354. return hst;
  355. }*/
  356.  
  357. /*
  358. TH1 *TbAsciiRoot::show_noise()
  359. {
  360. int ic;
  361. TH1 *hst = create_h1("hNoise","Noise",nchan(),-0.5, nchan()-0.5);
  362. if (gain()==1)
  363. {
  364. hst->SetYTitle("ADCs");
  365. }
  366. else
  367. {
  368. hst->SetYTitle("e^{-} ENC");
  369. }
  370. hst->SetXTitle("Channel no.");
  371. for (ic=0; ic<nchan(); ic++)
  372. hst->SetBinContent(ic+1, noise(ic));
  373.  
  374. return hst;
  375. }
  376. */
  377.  
  378. void TbAsciiRoot::compute_pedestals_fast(int mxevts, double wped, double wnoise)
  379. {
  380. if (!ifile)
  381. return;
  382.  
  383. if (mxevts<0)
  384. mxevts = 100000000;
  385.  
  386. for (int i=0;i<max_nchan;i++)
  387. _ped[i] = _noise[i] = 0.;
  388.  
  389. // std::ifstream::pos_type here = ifile->tellg();
  390. std::cout << "Computing fast pedestals..." << std::endl;
  391. std::string error="";
  392. for (int ievt=0; read_event(error)==0 && ievt<mxevts; ievt++)
  393. {
  394. if (!(ievt%100))
  395. {
  396. std::cout << "\revent " << std::setw(10) << ievt << std::flush;
  397. }
  398. common_mode();
  399. for (int i=0; i<nchan(); i++)
  400. {
  401. // TODO: figure out how to determine the chip number when
  402. // Plugin::filter_event has been called
  403. int ichip = i/128;
  404. // IF noise is 0, set it arbitrarily to 1.
  405. if (_noise[i]==0.)
  406. _noise[i] = 1.;
  407.  
  408. if (_ped[i]==0.)
  409. {
  410. // If pedestal is not yet computed we assume the current
  411. // channel value should not be too far
  412. _ped[i] = _data.data[i];
  413. }
  414. else
  415. {
  416. // Do the pedestal and noise correction
  417. double corr;
  418. double xs;
  419.  
  420. _signal[i] = _data.data[i] - _ped[i];
  421. corr = _signal[i] * wped;
  422.  
  423. xs = (_signal[i]-_cmmd[ichip])/_noise[i];
  424. if (corr > 1.)
  425. corr = 1.;
  426.  
  427. if (corr < -1)
  428. corr = -1.;
  429.  
  430. _ped[i] += corr;
  431.  
  432. if (fabs(xs) < 3.)
  433. {
  434. _noise[i] = _noise[i]*(1.0-wnoise) + xs*xs*wnoise;
  435. }
  436. }
  437. }
  438. }
  439. std::cout << "\nDone" << std::endl;
  440. rewind();
  441. }
  442.  
  443. /*
  444. TH2 *TbAsciiRoot::compute_pedestals(int mxevts, bool do_cmmd)
  445. {
  446. if (!ifile)
  447. return 0;
  448.  
  449. if (mxevts<0)
  450. mxevts = 100000000;
  451.  
  452. int ievt, ichan;
  453. TH2 *hst = create_h2("hRaw","Raw data",nchan(), -0.5,nchan()-0.5, 256, -0.5,1023.5);
  454. TH2 *hsts = create_h2("hSig","Signal",nchan(), -0.5,nchan()-0.5,256, -127.5,127.5);
  455.  
  456.  
  457. std::ifstream::pos_type here = ifile->tellg();
  458. std::cout << "Computing pedestas..." << std::endl;
  459. for (ievt=0; read_event()==0 && ievt<mxevts; ievt++)
  460. {
  461. process_event(do_cmmd);
  462. for (ichan=0; ichan<nchan(); ichan++)
  463. // TODO: get right chip number in all situations (after calling set_data)
  464. hst->Fill(ichan, data(ichan)-get_cmmd(ichan/128));
  465.  
  466. if (!(ievt%100))
  467. {
  468. std::cout << "\revent " << std::setw(10) << ievt << std::flush;
  469. }
  470. }
  471. std::cout << "\nDone" << std::endl;
  472. rewind();
  473.  
  474. // TODO: _nchan can be updated in an event by event basis
  475. // while here we are assuming that it is the same
  476. // for all the events
  477. for (ichan=0; ichan<nchan(); ichan++)
  478. {
  479. TF1 *g = new TF1("g1", "gaus");
  480. TH1 *h1 = hst->ProjectionY("__hx__", ichan+1, ichan+1);
  481. g->SetParameters(h1->GetSumOfWeights(), h1->GetMean(), h1->GetRMS());
  482. g->SetRange(h1->GetMean()-2.5*h1->GetRMS(), h1->GetMean()+2.5*h1->GetRMS());
  483. h1->Fit("g1", "q0wr");
  484. _ped[ichan] = h1->GetFunction("g1")->GetParameter(1);
  485. _noise[ichan] = h1->GetFunction("g1")->GetParameter(2);
  486. delete h1;
  487. delete g;
  488. }
  489.  
  490. rewind();
  491. for (ievt=0; read_event()==0 && ievt<mxevts; ievt++)
  492. {
  493. process_event(do_cmmd);
  494. for (ichan=0; ichan<nchan(); ichan++)
  495. hsts->Fill(ichan, signal(ichan));
  496.  
  497. if (!(ievt%100))
  498. {
  499. std::cout << "\revent " << std::setw(10) << ievt << std::flush;
  500. }
  501. }
  502. std::cout << "\nDone" << std::endl;
  503. rewind();
  504.  
  505. return hst;
  506. }
  507. */
  508.  
  509.  
  510.  
  511. void TbAsciiRoot::find_clusters(int ichip)
  512. {
  513. int chan0=0;
  514. int chan1=255;
  515. if (ichip>=0 && ichip<2)
  516. {
  517. chan0 = ichip*128;
  518. chan1 = (ichip+1)*128 -1;
  519. }
  520.  
  521. std::ostringstream ostr;
  522. ostr << chan0 << '-' << chan1;
  523. ChanList C(ostr.str().c_str());
  524.  
  525. clear();
  526. find_clusters(C);
  527. _hits = C.hit_list();
  528. }
  529.  
  530. void TbAsciiRoot::find_clusters(ChanList &C)
  531. {
  532. // TODO: figure out how to determine the chip number in
  533. // all the situations
  534. int i, j, imax=-1, left, right;
  535. double mxsig=-1.e20, sg, val;
  536. std::vector<bool> used(C.Nch());
  537.  
  538. for (i=0;i<C.Nch();i++)
  539. {
  540. used[i]= _mask[C[i]] ? true : false;
  541. }
  542.  
  543.  
  544.  
  545. while (true)
  546. {
  547. /*
  548. * Find the highest
  549. */
  550. imax = -1;
  551. for (j=0; j<C.Nch(); j++)
  552. {
  553. i = C[j];
  554. if (used[j] || _signal[i]*polarity()<0.)
  555. continue;
  556.  
  557. if ( polarity()*sn(i) > _seedcut)
  558. {
  559. val = fabs(signal(i));
  560. if (mxsig<val)
  561. {
  562. mxsig = val;
  563. imax = j;
  564. }
  565. }
  566. }
  567.  
  568. if (imax<0 || imax >= C.Nch() )
  569. break;
  570.  
  571. sg = signal(C[imax]);
  572. used[imax]=true;
  573. // Now look at neighbors
  574. // first to the left
  575. left = imax;
  576. for (j=imax-1;j>=0;j--)
  577. {
  578. i = C[j];
  579. if ( used[j] || _signal[i]*polarity()<0.)
  580. break;
  581.  
  582. if ( fabs(sn(i)) > _neighcut )
  583. {
  584. used[j] = true;
  585. sg += signal(i);
  586. left = j;
  587. }
  588. else
  589. // TODO: this needs to be removed
  590. // The idea is to merge to clusters that have only one strip in between
  591. // In the laser runs this is a consequences of reflections...
  592. {
  593. int jx = j-1;
  594. if (jx>=0 )
  595. {
  596. if ( fabs(sn(C[jx])) > _neighcut )
  597. continue;
  598. }
  599. break;
  600. }
  601. }
  602.  
  603. // now to the right
  604. right = imax;
  605. for (j=imax+1;j<C.Nch();j++)
  606. {
  607. i = C[j];
  608. if ( used[j] || _signal[i]*polarity()<0.)
  609. break;
  610. if ( fabs(sn(i))>_neighcut )
  611. {
  612. used[j] = true;
  613. sg += signal(i);
  614. right = j;
  615. }
  616. else
  617. // TODO: this needs to be removed
  618. // The idea is to merge to clusters that hanve only one strip in between
  619. // In the laser runs this is a consequences of reflections...
  620. {
  621. int jx = i+1;
  622. if (jx<C.Nch())
  623. {
  624. if ( fabs(sn(C[jx])) > _neighcut )
  625. continue;
  626. }
  627. break;
  628. }
  629. }
  630. C.add_hit(TbAlibavaHit(imax, left, right, sg));
  631. }
  632. }
  633.  
  634. void TbAsciiRoot::save_pedestals(const char *fnam)
  635. {
  636. std::ofstream ofile(fnam);
  637. if (!ofile)
  638. {
  639. std::cout << "Could not open " << fnam << " to save pedestals." << std::endl;
  640. return;
  641. }
  642.  
  643. // TODO: _nchan can be updated in an event by event basis
  644. // while here we are assuming that it is the same
  645. // for all the events
  646. int i;
  647. for (i=0; i<nchan(); i++)
  648. {
  649. ofile << _ped[i] << "\t" << _noise[i] << "\n";
  650. }
  651. ofile.close();
  652. }
  653.  
  654. void TbAsciiRoot::load_pedestals(const char *fnam)
  655. {
  656. std::ifstream ifile(fnam);
  657. if (!ifile)
  658. {
  659. std::cout << "Could not open " << fnam << " to load pedestals." << std::endl;
  660. return;
  661. }
  662. int i;
  663. for (i=0; i<max_nchan; i++)
  664. {
  665. if (ifile.eof())
  666. break;
  667.  
  668. ifile >> _ped[i] >> std::ws >> _noise[i] >> std::ws;
  669. _mask[i] = (_noise[i]>20. || _noise[i]<=0.);
  670. }
  671. ifile.close();
  672. /*
  673. TCanvas *pedcanvas = create_canvas("Pedcanvas", "Pedestal Values", 600, 400);
  674. TH1 *pedestalhisto = create_h1("pedestalhisto", "Pedestal Values", 256, -0.5, 255.5);
  675. for (i=0; i<256; i++)
  676. {
  677. pedestalhisto->Fill(i, _ped[i]);
  678. }
  679. pedcanvas->cd(1);
  680. pedestalhisto->Draw();
  681. */
  682. }
  683.  
  684. void TbAsciiRoot::load_masking(const char *fnam)
  685. {
  686. std::ifstream ifile(fnam);
  687. if (!ifile)
  688. {
  689. std::cout << "Could not open masked.txt. " << std::endl;
  690. return;
  691. }
  692. int val;
  693. for (int i=0; i<500; i++)
  694. {
  695. ifile >> val >> std::ws;
  696. if (ifile.eof())
  697. break;
  698. if (val>255)
  699. {
  700. std::cout << "A value is greater than 255, causing an overflow crash. Please check the text file again. It has been set to 1 for continuation purposes. " << std::endl;
  701. val = 1;
  702. }
  703. _mask[val] = true;
  704. }
  705. }
  706.  
  707. void TbAsciiRoot::load_gain(const char *fnam)
  708. {
  709. std::ifstream ifile(fnam);
  710. if (!ifile)
  711. {
  712. std::cout << "Could not open " << fnam << " to load the gain." << std::endl;
  713. return;
  714. }
  715. int i;
  716. int ichan;
  717. double val, xn, xm;
  718. xn=xm=0.;
  719. for (i=0; i<max_nchan; i++)
  720. {
  721. ifile >> ichan >> std::ws;
  722. if (ifile.eof())
  723. break;
  724. ifile >> val;
  725. if (ifile.eof())
  726. break;
  727.  
  728. xn++;
  729.  
  730. xm += val;
  731. _gain[ichan] = val;
  732.  
  733. ifile >> std::ws;
  734. if (ifile.eof())
  735. break;
  736. }
  737. if (xn>0)
  738. {
  739. _average_gain = xm/xn;
  740. }
  741. ifile.close();
  742. }
  743.  
  744. void TbAsciiRoot::process_event(bool do_cmmd)
  745. {
  746. int i;
  747. for (i=0; i<nchan(); i++)
  748. {
  749. _signal[i] = _data.data[i]-_ped[i];
  750. _sn[i] = _noise[i]>1. && !_mask[i] ? _signal[i]/_noise[i] : 0.;
  751. }
  752. if (do_cmmd)
  753. {
  754. int ichip=-1;
  755. common_mode();
  756.  
  757. for (i=0; i<nchan(); i++)
  758. {
  759. // TODO: figure out the right chip number
  760. if (!(i%128))
  761. ichip ++;
  762.  
  763. _signal[i] = _data.data[i]-_ped[i] - _cmmd[ichip];
  764. _sn[i] = (_noise[i] >1. && !_mask[i] ? _signal[i]/_noise[i] : 0.);
  765. }
  766. }
  767. }
  768.  
  769. void TbAsciiRoot::add_channel_list(const ChanList &C)
  770. {
  771. chan_list.push_back(C);
  772. }
  773.  
  774.  
  775. void TbAsciiRoot::common_mode()
  776. {
  777. ChanList C("0-127");
  778. common_mode(C);
  779.  
  780. _cmmd[0] = C.CommonMode();
  781. _cnoise[0] = C.Noise();
  782.  
  783.  
  784. C.Set("128-255");
  785. common_mode(C);
  786.  
  787. _cmmd[1] = C.CommonMode();
  788. _cnoise[1] = C.Noise();
  789. }
  790.  
  791. void TbAsciiRoot::common_mode(ChanList &C, bool correct)
  792. {
  793. int ip, i, j;
  794.  
  795. double mean, sm, xn, xx, xs, xm, tmp;
  796. bool use_it;
  797. mean = sm = 0.;
  798. for (ip=0;ip<3;ip++)
  799. {
  800. xn = xs = xm = 0.;
  801. for (j=0; j<C.Nch(); j++)
  802. {
  803. i = C[j];
  804. if (_mask[i])
  805. continue;
  806.  
  807. use_it = true;
  808. xx = data(i) - _ped[i];
  809. if (ip)
  810. {
  811. tmp = fabs((xx-mean)/sm);
  812. use_it = (tmp<2.5);
  813. }
  814. if (use_it)
  815. {
  816. xn++;
  817. xm += xx;
  818. xs += xx * xx;
  819. }
  820. }
  821. if (xn>0.)
  822. {
  823. mean = xm / xn;
  824. sm = sqrt( xs/xn - mean*mean);
  825. }
  826. // std::cout << "...iter " << ip << ": xm " << mean << " xs: " << sm << std::endl;
  827. }
  828. C.CommonMode(mean);
  829. C.Noise(sm);
  830.  
  831. if (correct)
  832. {
  833. for ( j=0; j<C.Nch(); j++ )
  834. {
  835. i = C[j];
  836. _signal[i] = _data.data[i]-_ped[i] - C.CommonMode();
  837. _sn[i] = (_noise[i] >1. && !_mask[i] ? _signal[i]/_noise[i] : 0.);
  838. }
  839. }
  840. }
  841.  
  842.  
  843.  
  844.  
  845. /*
  846. void TbAsciiRoot::spy_data(bool with_signal, int nevt)
  847. {
  848. TVirtualPad *pad;
  849. if (!ifile)
  850. return;
  851.  
  852. sighandler_t old_handler = ::signal(SIGINT, _A_got_intr);
  853. _A_do_run = true;
  854.  
  855. TCanvas *cnvs = (TCanvas *)gROOT->FindObject("cnvs");
  856. if (cnvs)
  857. {
  858. cnvs->Clear();
  859. }
  860. else
  861. cnvs = new TCanvas("cnvs","cnvs", 700, 800);
  862.  
  863. cnvs->Divide(2,3);
  864.  
  865.  
  866. TH1 *hsignal = create_h1("hsignal","signal (ADC)",256, -0.5, 255.0);
  867. hsignal->SetXTitle("Channel");
  868. hsignal->SetYTitle("ADC");
  869. hsignal->SetMinimum(-300);
  870. hsignal->SetMaximum(300);
  871.  
  872. TH1 *helec = create_h1("helec","signal (elec)", 256, -0.5, 255.5);
  873. helec->SetXTitle("Channel");
  874. helec->SetYTitle("electrons");
  875. helec->SetMinimum(-300/gain());
  876. helec->SetMaximum(300/gain());
  877.  
  878. TH1 *hraw = create_h1("hraw","Raw Data (around 512.)",256, 0., 256.);
  879. hraw->SetXTitle("Channel");
  880. hraw->SetYTitle("ADC");
  881. hraw->SetMinimum(-300);
  882. hraw->SetMaximum(+300);
  883.  
  884. TH1 *hrawc = create_h1("hrawc","Raw Data (no commd)",256, 0., 256.);
  885. hrawc->SetXTitle("Channel");
  886. hrawc->SetYTitle("ADC");
  887. hrawc->SetMinimum(-300);
  888. hrawc->SetMaximum(+300);
  889.  
  890.  
  891. TH1 *hcmmd[2];
  892. hcmmd[0] = create_h1("hcmmd0","Common mode (Chip 0)",50,-100.,100.);
  893. hcmmd[0]->SetXTitle("Common mode");
  894. hcmmd[1] = create_h1("hcmmd1","Common mode (Chip 1)",50,-100.,100.);
  895. hcmmd[1]->SetXTitle("Common mode");
  896.  
  897. int ievt,jevt;
  898. for (ievt=jevt=0; read_event()==0 && _A_do_run && ievt<nevt;jevt++)
  899. {
  900. process_event();
  901. find_clusters();
  902. if ( with_signal && empty())
  903. continue;
  904.  
  905. int i,ichip=-1;
  906. for (i=0; i<nchan(); i++)
  907. {
  908. // TODO: figure out chip number
  909. if (!(i%128))
  910. ichip++;
  911.  
  912. hsignal->SetBinContent(i+1, _signal[i]);
  913. helec->SetBinContent(i+1, signal(i));
  914. hraw->SetBinContent(i+1,data(i)-512.);
  915. hrawc->SetBinContent(i+1, data(i)-_ped[i]);
  916. // TODO: why we draw the signal + common mode ?
  917. // May be cause signal should be ~0...
  918. hcmmd[ichip]->Fill(_signal[i]+get_cmmd(ichip));
  919. }
  920. pad = cnvs->cd(1);
  921. pad->SetGrid(1,1);
  922. hsignal->Draw();
  923. pad = cnvs->cd(2);
  924. pad->SetGrid(1,1);
  925. helec->Draw();
  926.  
  927. pad = cnvs->cd(3);
  928. pad->SetGrid(1,1);
  929. hraw->Draw();
  930.  
  931. pad = cnvs->cd(4);
  932. pad->SetGrid(1,1);
  933. hrawc->Draw();
  934.  
  935. pad = cnvs->cd(5);
  936. pad->SetGrid(1,1);
  937. hcmmd[0]->Draw();
  938.  
  939. pad = cnvs->cd(6);
  940. pad->SetGrid(1,1);
  941. hcmmd[1]->Draw();
  942.  
  943. std::cout << std::setiosflags(std::ios::fixed);
  944. std::cout << "*** Event " << jevt << " *****" << std::endl;
  945. std::cout << "Common Mode:" << std::endl
  946. << " Chip 0 " << std::setw(6) << std::setprecision(1) << get_cmmd(0) << " noise: " << get_cnoise(0)
  947. << std::endl
  948. << " Chip 1 " << std::setw(6) << std::setprecision(1) << get_cmmd(1) << " noise: " << get_cnoise(1)
  949. << std::endl;
  950.  
  951. std::cout << "Time: " << time() << " ns" << std::endl;
  952. std::cout << "Signal chan(0) << " << signal(0) << " chan(1) " << signal(1) << std::endl;
  953. std::cout << "Clusters: " << std::endl;
  954.  
  955. TbAlibavaHitList::iterator ip;
  956. for (ip=begin(); ip!=end(); ++ip)
  957. {
  958. std::cout << " chan: " << ip->center()
  959. << " sig: "
  960. << std::setw(6) << std::setprecision(1) << ip->signal()
  961. << " left: " << ip->left() << " right: " << ip->right()
  962. << std::endl;
  963. std::cout << '\t' << "channels: " << std::endl;
  964. int j;
  965. for (j=ip->left();j<=ip->right();j++)
  966. std::cout << "\t " << j << " sn: " << _sn[j] << " signal: " << _signal[j] << " noise: " << _noise[j] << '\n';
  967. std::cout << std::endl;
  968. }
  969.  
  970. cnvs->Update();
  971. ievt++;
  972. }
  973. std::cout << std::endl;
  974. _A_do_run= true;
  975. ::signal(SIGINT, old_handler);
  976. }
  977. */
  978.  
  979. bool is_text(const char *fnam)
  980. {
  981. int nc;
  982. char buffer[1024];
  983. std::ifstream ifile(fnam);
  984. if (!fnam)
  985. return false;
  986.  
  987. ifile.read(buffer, sizeof(buffer));
  988. nc = ifile.gcount();
  989. ifile.close();
  990. if (!nc) // empty files are text
  991. {
  992. return true;
  993. }
  994.  
  995. std::string ss(buffer, nc);
  996. ifile.close();
  997.  
  998. if ( ss.find('\0') != ss.npos )
  999. return false;
  1000.  
  1001. double nontext = 0.;
  1002. double ntotal = 0.;
  1003. std::string::iterator ip;
  1004. for (ip=ss.begin(); ip!=ss.end(); ++ip)
  1005. {
  1006. ntotal++;
  1007. char c = *ip;
  1008. if ( (c<' ' || c >'~') && !strchr("\n\t\r\b", c) )
  1009. nontext++;
  1010. }
  1011. if ( nontext/ntotal > 0.3 )
  1012. return false;
  1013.  
  1014. return true;
  1015. }
  1016.