source: trunk/Mars/mraw/MRawEvtData.cc@ 17050

Last change on this file since 17050 was 11870, checked in by tbretz, 13 years ago
When setting a pixel map in the event data subtract one to get from the id to the index.
File size: 20.1 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Thomas Bretz, 12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2008
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MRawEvtData
28//
29// Storage container to store the raw FADC values.
30//
31// MArrayS fHiGainPixId
32// ---------------------
33// Array of Pixel Numbers for their high voltage channel in the order the
34// FADC values are stored in fHiGainFadcSamples
35//
36// MArrayB fHiGainFadcSaples
37// -------------------------
38// FADC samples (hi gain) of all pixels
39//
40// MArrayS fLoGainPixId
41// --------------------
42// see fHiGainPixId
43//
44// MArrayB fLoGainFadcSamples
45// --------------------------
46// see fHiGainFadcSaples
47//
48//
49// Version 8
50// ------------------
51// + MArrayS *fStartCell
52// + Bool_t fIsSigned
53//
54// Version 7
55// ------------------
56// + UShort_t fNumBytesPerSample;
57//
58// Version 6
59// ------------------
60// - The data can now be stoe in a single array keeping he full
61// compatibility
62//
63// Version 5 (0.8.5):
64// ------------------
65// - Changed type of ABFlags from TArrayC to MArrayB
66//
67// Version 4 (0.8.4):
68// ------------------
69// - Changed derivement from MCamEvent to MParContainer and MCamEvent
70//
71// Version 3 (0.8.4):
72// ------------------
73// - Added fABFlags
74//
75// Version 2:
76// ----------
77// - Derives from MCamEvent now
78//
79// Version 1:
80// ----------
81// - First implementation
82//
83/////////////////////////////////////////////////////////////////////////////
84#include "MRawEvtData.h"
85
86#include <fstream>
87
88#include <TH1.h>
89#include <TGraph.h>
90#include <TArrayC.h>
91#include <TVirtualPad.h>
92
93#include "MLog.h"
94#include "MLogManip.h"
95
96#include "MArrayS.h"
97#include "MArrayB.h"
98#include "MArrayI.h"
99#include "MGeomCam.h"
100
101#include "MRawCrateArray.h"
102#include "MRawCrateData.h"
103
104#include "MRawRunHeader.h"
105#include "MRawEvtPixelIter.h"
106
107ClassImp(MRawEvtData);
108
109using namespace std;
110
111// --------------------------------------------------------------------------
112//
113// Default constructor. It initializes all arrays with zero size.
114//
115MRawEvtData::MRawEvtData(const char *name, const char *title)
116 : fRunHeader(0), fNumBytesPerSample(1), fIsSigned(false)
117{
118 fName = name ? name : "MRawEvtData";
119 fTitle = title ? title : "Raw Event Data Information";
120
121 InitArrays();
122}
123
124// --------------------------------------------------------------------------
125//
126// Destructor. Deletes all the arrays.
127//
128MRawEvtData::~MRawEvtData()
129{
130 DeleteArrays();
131}
132
133// --------------------------------------------------------------------------
134//
135// reset all arrays
136//
137void MRawEvtData::Clear(Option_t *)
138{
139 /*
140 FIXME:
141 Is Reset (set all entries to zero) what you want to do
142 or Set(0) (delete the array)
143 */
144 fHiGainPixId->Reset();
145 fLoGainPixId->Reset();
146 fHiGainFadcSamples->Reset();
147 fLoGainFadcSamples->Reset();
148 fABFlags->Reset();
149}
150
151// --------------------------------------------------------------------------
152//
153// return the number of hi gain samples per pixel
154//
155UShort_t MRawEvtData::GetNumHiGainSamples() const
156{
157 // If value<0 we are reading old MC files which don't have the
158 // value set so far. So we use the old methid to determin the
159 // numbers and calculate them from the length of the arrays.
160 return fHiGainPixId->GetSize() ? fHiGainFadcSamples->GetSize()/(fHiGainPixId->GetSize()*fNumBytesPerSample) : 0;
161}
162
163// --------------------------------------------------------------------------
164//
165// return the number of lo gain samples per pixel
166//
167UShort_t MRawEvtData::GetNumLoGainSamples() const
168{
169 // If value<0 we are reading old MC files which don't have the
170 // value set so far. So we use the old methid to determin the
171 // numbers and calculate them from the length of the arrays.
172 return fLoGainPixId->GetSize() ? fLoGainFadcSamples->GetSize()/(fLoGainPixId->GetSize()*fNumBytesPerSample) : 0;
173}
174
175// --------------------------------------------------------------------------
176//
177// return the number of stored pixel
178//
179UShort_t MRawEvtData::GetNumPixels() const
180{
181 return fHiGainPixId ? fHiGainPixId->GetSize() : 0;
182}
183
184// --------------------------------------------------------------------------
185//
186// Print out the onformation to *fLog.
187// Options:
188// "hex" Prints the time slices hexadecimal (default)
189// "dec" Prints the time slices decimal
190//
191void MRawEvtData::Print(Option_t *opt) const
192{
193 //
194 // print fadc inforation to screen
195 // Possible Options:
196 // - DEC: Print values decimal instead of hexadecimal (default)
197 //
198 const UShort_t nHiSamp = GetNumHiGainSamples();
199 const UShort_t nLoSamp = GetNumLoGainSamples();
200
201 const UShort_t bps = GetNumBytesPerSample();
202
203 fLog->unsetf(ios::showbase);
204
205 *fLog << dec << all;
206 *fLog << GetDescriptor() << ": " << endl;
207 *fLog << GetNumPixels() << " Pixels with " << (Int_t)nHiSamp << "/" << (Int_t)nLoSamp << " samples" << endl;
208
209 TString str(opt);
210 Int_t manip = str.Contains("dec", TString::kIgnoreCase);
211
212 MRawEvtPixelIter pixel(const_cast<MRawEvtData*>(this));
213 Int_t i = 0;
214 while (pixel.Next())
215 {
216 const UShort_t idx = pixel.GetPixelId();
217
218 *fLog << endl << dec << setfill(' ');
219 *fLog << " " << setw(3) << i++ << " - " << setw(3) << idx << " ";
220
221 if (pixel.IsABFlagValid())
222 *fLog << "<" << (pixel.HasABFlag()?"B":"A") << "> ";
223
224 *fLog << (manip?dec:hex) << setfill(manip?' ':'0');
225
226 const Byte_t *hi = (Byte_t*)pixel.GetHiGainSamples();
227 const Byte_t *lo = (Byte_t*)pixel.GetLoGainSamples();
228
229 for (int j=0; j<nHiSamp*bps; j++)
230 {
231 *fLog << setw(manip?3:2);
232 *fLog << (hi[j]&0xff);
233 if (manip)
234 *fLog << ' ';
235 }
236
237 for (int j=0; j<nLoSamp*bps; j++)
238 {
239 *fLog << setw(manip?3:2);
240 *fLog << ((UShort_t)lo[j]&0xff);
241 if (manip)
242 *fLog << ' ';
243 }
244 }
245 *fLog << dec << endl;
246}
247
248// --------------------------------------------------------------------------
249//
250// Draw a pixel. A Histogram or Graph is created and it's draw function is
251// called.
252// Options:
253// "GRAPH" A graph is drawn
254// "HIST" A histogram is drawn
255// <index> The pixel with the given index is drawn
256//
257void MRawEvtData::Draw(Option_t *opt)
258{
259 if (GetNumPixels()==0)
260 {
261 *fLog << warn << "Sorry, no pixel to draw!" << endl;
262 return;
263 }
264
265 TString str(opt);
266 str.ToLower();
267
268 UInt_t id = 0;
269
270 if (str.BeginsWith("graph"))
271 sscanf(str.Data()+5, "%d", &id);
272 if (str.BeginsWith("hist"))
273 sscanf(str.Data()+4, "%d", &id);
274
275 MRawEvtPixelIter pix(this);
276 if (!pix.Jump(id))
277 {
278 *fLog << warn << dec << "Pixel Idx #" << id << " doesn't exist!" << endl;
279 return;
280 }
281
282 const void *higains = pix.GetHiGainSamples();
283 const void *logains = pix.GetLoGainSamples();
284
285 const Int_t nh = GetNumHiGainSamples();
286 const Int_t nl = GetNumLoGainSamples();
287
288 TString name = "Pixel Idx.";
289 name += pix.GetPixelId();
290
291 fIsSigned = kTRUE;
292
293 Bool_t same = str.Contains("same");
294
295 cout << "MIN/MAX=" << fNumBytesPerSample << " " << GetMin() << " " << GetMax() << endl;
296
297 if (str.BeginsWith("graph"))
298 {
299 *fLog << inf << "Drawing Graph: Pixel Idx #" << dec << pix.GetPixelId();
300 *fLog << " of " << (int)GetNumPixels() << "Pixels" << endl;
301
302 TGraph *graphhi = new TGraph;
303
304 for (int i=0; i<nh; i++)
305 graphhi->SetPoint(graphhi->GetN(), i, GetSample(higains, i));
306 for (int i=0; i<nl; i++)
307 graphhi->SetPoint(graphhi->GetN(), i+nh, GetSample(logains, i));
308
309 graphhi->SetMaximum(GetMax()+0.5);
310 graphhi->SetMinimum(GetMin());
311
312 graphhi->SetBit(kCanDelete);
313 graphhi->Draw(same ? "C*" : "AC*");
314
315 TH1F *histhi = graphhi->GetHistogram();
316 histhi->SetMinimum(GetMin());
317 histhi->SetMaximum(GetMax()+0.5);
318
319 histhi->SetXTitle("Time/FADC Slices");
320 histhi->SetYTitle("Signal/FADC Units");
321
322 return;
323 }
324
325 if (str.BeginsWith("hist"))
326 {
327 // FIXME: Add Legend
328 *fLog << inf << "Drawing Histogram of Pixel with Idx #" << dec << pix.GetPixelId() << " to " << gPad << endl;
329
330 TH1F *histh = new TH1F(name, "FADC Samples", nh+nl, -0.5, nh+nl-.5);
331 histh->SetMinimum(GetMin());
332 histh->SetMaximum(GetMax()+0.5);
333 histh->SetXTitle("Time [FADC Slices]");
334 histh->SetYTitle("Signal [FADC Units]");
335 histh->SetDirectory(NULL);
336 for (int i=0; i<nh; i++)
337 histh->Fill(i, GetSample(higains, i));
338 for (int i=0; i<nl; i++)
339 histh->Fill(i+nl, GetSample(logains, i));
340 histh->SetBit(kCanDelete);
341 histh->Draw(same ? "same" : "");
342
343 return;
344 }
345
346 *fLog << warn << dbginf << "WARNING - You must specify either 'GRAPH' or 'HIST'" << endl;
347}
348
349// --------------------------------------------------------------------------
350//
351// Deletes all the arrays
352// The flag is for future usage.
353//
354void MRawEvtData::InitArrays(UShort_t numconnected, UShort_t maxid)
355{
356 // fRunHeader should not be set only in the constructor!
357 const Int_t numhi = fRunHeader ? fRunHeader->GetNumSamplesHiGain() : 0;
358 const Int_t numlo = fRunHeader ? fRunHeader->GetNumSamplesLoGain() : 0;
359
360 fNumBytesPerSample = fRunHeader ? fRunHeader->GetNumBytesPerSample() : 1;
361
362 fHiGainPixId = new MArrayS(numconnected);
363 fLoGainPixId = new MArrayS(0);
364 fHiGainFadcSamples = new MArrayB(numconnected*(numhi+numlo)*fNumBytesPerSample);
365 fLoGainFadcSamples = new MArrayB(0);
366
367 fABFlags = new MArrayB(maxid==0 ? 0 : maxid/8+1);
368 fStartCells = new MArrayS(0);
369
370 fConnectedPixels = 0;
371}
372
373void MRawEvtData::InitStartCells()
374{
375 delete fStartCells;
376 fStartCells = new MArrayS(fHiGainPixId->GetSize());
377}
378
379// --------------------------------------------------------------------------
380//
381// Deletes all the arrays
382//
383void MRawEvtData::DeleteArrays()
384{
385 delete fHiGainPixId;
386 delete fLoGainPixId;
387 delete fHiGainFadcSamples;
388 delete fLoGainFadcSamples;
389 delete fABFlags;
390 delete fStartCells;
391}
392
393// --------------------------------------------------------------------------
394//
395// Deletes all arrays describing the pixel Id and Samples in pixels.
396// The flag is for future usage.
397//
398void MRawEvtData::ResetPixels(UShort_t numconnected, UShort_t maxid)
399{
400 //const int npix = fRunHeader->GetNumCrates()*fRunHeader->GetNumPixInCrate();
401 if (fHiGainPixId && fHiGainPixId->GetSize()==numconnected && (UShort_t)fABFlags->GetSize()==(maxid/8+1))
402 {
403 fConnectedPixels = 0;
404 return;
405 }
406
407 DeleteArrays();
408 InitArrays(numconnected, maxid);
409}
410
411void MRawEvtData::ResetPixels()
412{
413 if (!fRunHeader)
414 return;
415
416 // FIXME: Better give NumNormalPixels if numconnected==0 ?
417
418 ResetPixels(fRunHeader->GetNumNormalPixels(), fRunHeader->GetNumNormalPixels()-1);
419}
420
421// --------------------------------------------------------------------------
422//
423// This is to fill the data of one pixel to the MRawEvtHeader Class.
424// The parameters are the pixelnumber and the FADC_SLICES values of ADCs
425//
426void MRawEvtData::AddPixel(UShort_t nOfPixel, const TArrayC &data)
427{
428 const Int_t n = fRunHeader->GetNumSamples();
429 if (data.GetSize()!=n)
430 {
431 *fLog << err << "RawEvtData::AddPixel: Error, number of samples in ";
432 *fLog << "TArrayC " << data.GetSize() << " doesn't match current number " << n << endl;
433 return;
434 }
435
436 //
437 // enhance pixel array by one
438 //
439 fHiGainPixId->Set(fHiGainPixId->GetSize()+1);
440
441 //
442 // add the number of the new pixel to the array as last entry
443 //
444 fHiGainPixId->AddAt(nOfPixel, fHiGainPixId->GetSize()-1);
445
446 //
447 // enhance the array by the number of new samples
448 //
449 fHiGainFadcSamples->Set(fHiGainFadcSamples->GetSize()+n);
450
451 //
452 // add the new slices as last entries to array
453 //
454 fHiGainFadcSamples->AddAt((Byte_t*)data.GetArray(), fHiGainFadcSamples->GetSize()-n, n);
455}
456
457// --------------------------------------------------------------------------
458//
459// Add the contents of the MArrayI to the fHiGainFadcSamples
460// One Integer is added to one sample in the array.
461//
462void MRawEvtData::Set(const MArrayI &data)
463{
464 fConnectedPixels = fHiGainPixId->GetSize();
465
466 UInt_t n = fConnectedPixels*fRunHeader->GetNumSamplesHiGain();
467 if (n!=data.GetSize())
468 {
469 *fLog << err << "MRawEvtData::Set: Size mismatch." << endl;
470 *fLog << " fConnectedPixels=" << fConnectedPixels << endl;
471 *fLog << " NumHiGainSamples=" << fRunHeader->GetNumSamplesHiGain() << endl;
472 *fLog << " data.GetSize()=" << data.GetSize() << endl;
473 return;
474 }
475
476 Byte_t *dest = fHiGainFadcSamples->GetArray();
477
478 UInt_t *src = reinterpret_cast<UInt_t*>(data.GetArray());
479 UInt_t *end = reinterpret_cast<UInt_t*>(data.GetArray()) + n;
480
481 switch (fNumBytesPerSample)
482 {
483 case 1:
484 {
485 Byte_t *ptr = reinterpret_cast<Byte_t*>(dest);
486 while (src<end)
487 *ptr++ = Byte_t(*src++);
488 }
489 return;
490
491 case 2:
492 {
493 UShort_t *ptr = reinterpret_cast<UShort_t*>(dest);
494 while (src<end)
495 *ptr++ = UShort_t(*src++);
496 }
497 return;
498
499 case 4:
500 memcpy(dest, data.GetArray(), n*4);
501 return;
502 }
503}
504
505// --------------------------------------------------------------------------
506//
507// Set indices according to the given array. The array must have the same
508// size as fHiGainPixId
509//
510void MRawEvtData::SetIndices(const MArrayS &idx)
511{
512 if (idx.GetSize()!=fHiGainPixId->GetSize())
513 return;
514
515 for (UInt_t i=0; i<fHiGainPixId->GetSize(); i++)
516 (*fHiGainPixId)[i] = idx[i]-1;
517}
518
519// --------------------------------------------------------------------------
520//
521// Set indices according to the length of the array fHiGainPixId
522// from 0 to n-1
523//
524void MRawEvtData::SetIndices()
525{
526 for (UInt_t i=0; i<fHiGainPixId->GetSize(); i++)
527 (*fHiGainPixId)[i] = i;
528}
529
530Byte_t *MRawEvtData::GetSamples() const { return fHiGainFadcSamples->GetArray(); }
531UShort_t *MRawEvtData::GetStartCells() const { return fStartCells->GetArray(); }
532UShort_t *MRawEvtData::GetPixelIds() const { return fHiGainPixId->GetArray(); }
533Bool_t MRawEvtData::HasStartCells() const { return fStartCells->GetSize()==fHiGainPixId->GetSize(); }
534
535/*
536void MRawEvtData::AddPixel(UShort_t nOfPixel, TArrayC *data, Bool_t lflag)
537{
538 MArrayS *arrpix = lflag ? fLoGainPixId : fHiGainPixId;
539 MArrayB *arrsam = lflag ? fLoGainFadcSamples : fHiGainFadcSamples;
540
541 //
542 // check whether we got the right number of new samples
543 // if there are no samples already stored: this is the new number of samples
544 //
545 const UShort_t ns = data->GetSize();
546 const UShort_t nSamp = lflag ? GetNumLoGainSamples() : GetNumHiGainSamples();
547 if (nSamp && ns!=nSamp)
548 {
549 *fLog << err << "RawEvtData::AddPixel: Error, number of samples in ";
550 *fLog << "TArrayC " << ns << " doesn't match current number " << nSamp << endl;
551 return;
552 }
553
554 //
555 // enhance pixel array by one
556 //
557 arrpix->Set(arrpix->GetSize()+1);
558
559 //
560 // add the number of the new pixel to the array as last entry
561 //
562 arrpix->AddAt(nOfPixel, arrpix->GetSize()-1);
563
564 //
565 // enhance the array by the number of new samples
566 //
567 arrsam->Set(arrsam->GetSize()+ns);
568
569 //
570 // add the new slices as last entries to array
571 //
572 arrsam->AddAt((Byte_t*)data->GetArray(), arrsam->GetSize()-ns, ns);
573}
574*/
575
576void MRawEvtData::ReadPixel(istream &fin, Int_t npix)
577{
578 // number of samples
579 const UInt_t ns = fRunHeader->GetNumSamples();
580
581 // bytes per sample
582 const Int_t bps = fRunHeader->GetNumBytesPerSample();
583
584 // Number of bytes
585 const Int_t nb = ns*bps;
586
587 // position in higain array
588 Byte_t *pos = fHiGainFadcSamples->GetArray() + fConnectedPixels*nb;
589
590 // Set pixel index
591 fHiGainPixId->AddAt(npix, fConnectedPixels++);
592
593 // Read data for one pixel
594 fin.read((char*)pos, nb);
595}
596
597void MRawEvtData::SetABFlag(Int_t npix, Bool_t ab)
598{
599 if (ab)
600 SETBIT((*fABFlags)[npix/8], npix%8);
601 else
602 CLRBIT((*fABFlags)[npix/8], npix%8);
603}
604
605// --------------------------------------------------------------------------
606//
607// Fills members with information from a magic binary file.
608// WARNING: you have to use Init() before you can do this
609//
610/*
611void MRawEvtData::ReadEvt(istream &fin, Int_t posinarray)
612{
613
614 const UShort_t npic = fRunHeader->GetNumPixInCrate();
615
616 const UShort_t npos = npic*posinarray;
617
618 //const Byte_t ab = fCrateArray->GetEntry(posinarray)->GetABFlags();
619
620 for (int i=npos; i<npic+npos; i++)
621 {
622 // calc the spiral hardware pixel number
623 const UShort_t ipos = i;
624
625 // Get Hardware Id
626 const Short_t hwid = fRunHeader->GetPixAssignment(ipos);
627
628 // Check whether the pixel is connected or not
629 if (hwid==0)
630 {
631 const UShort_t n = fRunHeader->GetNumSamplesLoGain()+fRunHeader->GetNumSamplesHiGain();
632 fin.seekg(n, ios::cur);
633 return;
634 }
635
636 // -1 converts the hardware pixel Id into the software pixel index
637 const Int_t npix = (Int_t)hwid-1;
638
639 const Byte_t ab = fCrateArray->GetEntry(posinarray)->GetABFlags();
640 AddPixel(fin, npix, TESTBIT(ab, i-npos));
641 }
642}
643*/
644
645// --------------------------------------------------------------------------
646//
647// Return the size in bytes of one event data block
648//
649Int_t MRawEvtData::GetNumBytes() const
650{
651 const UShort_t nlo = fRunHeader->GetNumSamplesLoGain();
652 const UShort_t nhi = fRunHeader->GetNumSamplesHiGain();
653 const UShort_t npic = fRunHeader->GetNumPixInCrate();
654 const UShort_t nbps = fRunHeader->GetNumBytesPerSample();
655
656 return (nhi+nlo)*npic*nbps;
657}
658
659// --------------------------------------------------------------------------
660//
661// Make sure, that you skip the whole event. This function only skips a part
662// of the event - see MRawRead::SkipEvent
663//
664void MRawEvtData::SkipEvt(istream &fin)
665{
666 fin.seekg(GetNumBytes(), ios::cur);
667}
668
669Bool_t MRawEvtData::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
670{
671 *fLog << warn << "WARNING - The use of MRawEvtData::GetPixelContent is deprecated!" << endl;
672
673 /*
674 MRawEvtPixelIter Next(const_cast<MRawEvtData*>(this));
675 if (!Next.Jump(idx))
676 return kFALSE;
677
678 switch (type)
679 {
680 case 0:
681 val = Next.GetSumHiGainSamples()-(float)GetNumHiGainSamples()*fHiGainFadcSamples->GetArray()[0];
682 val*= cam.GetPixRatio(idx);
683 break;
684 case 1:
685 val = Next.GetMaxHiGainSample();
686 break;
687 case 2:
688 val = Next.GetMaxLoGainSample();
689 break;
690 case 3:
691 val = Next.GetIdxMaxHiGainSample();
692 break;
693 case 4:
694 val = Next.GetIdxMaxLoGainSample();
695 break;
696 case 5:
697 val = Next.GetIdxMaxHiLoGainSample();
698 return val >= 0;
699 }
700*/
701 return kTRUE;
702}
703
704void MRawEvtData::Copy(TObject &named)
705#if ROOT_VERSION_CODE > ROOT_VERSION(3,04,01)
706const
707#endif
708{
709 MRawEvtData &evt = (MRawEvtData &)named;
710
711 *evt.fHiGainPixId = *fHiGainPixId;
712 *evt.fLoGainPixId = *fLoGainPixId;
713
714 *evt.fHiGainFadcSamples = *fHiGainFadcSamples;
715 *evt.fLoGainFadcSamples = *fLoGainFadcSamples;
716
717 *evt.fABFlags = *fABFlags;
718
719 evt.fConnectedPixels = fConnectedPixels;
720
721 evt.fNumBytesPerSample = fNumBytesPerSample;
722}
Note: See TracBrowser for help on using the repository browser.