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