source: trunk/Mars/msim/MPhotonEvent.cc @ 19338

Last change on this file since 19338 was 19338, checked in by tbretz, 8 months ago
This allows to distinguish between the thinning file format (8 floats per bunch) and the normal one (7 floats per bunch)
File size: 19.7 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of CheObs, the Modular 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 appears 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!   Author(s): Qi Zhe,       06/2007 <mailto:qizhe@astro.uni-wuerzburg.de>
20!
21!   Copyright: CheObs Software Development, 2000-2009
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27//
28//  MPhotonEvent
29//
30// Storage container to store photon collections
31//
32// The class is designed to be extremely fast which is important taking into
33// account the extremely high number of photons. This has some impacts on
34// its handling.
35//
36// The list has to be kept consistent, i.e. without holes.
37//
38// There are two ways to achieve this:
39//
40//   a) Use RemoveAt to remove an entry somewhere
41//   b) Compress() the TClonesArray afterwards
42//
43// Compress is not the fastes, so there is an easier way.
44//
45//   a) When you loop over the list and want to remove an entry copy all
46//      following entry backward in the list, so that the hole will
47//      be created at its end.
48//   b) Call Shrink(n) with n the number of valid entries in your list.
49//
50// To loop over the TClonesArray you can use a TIter which has some
51// unnecessary overhead and therefore is slower than necessary.
52//
53// Since the list is kept consistent you can use a simple loop saving
54// a lot of CPU time taking into account the high number of calls to
55// TObjArrayIter::Next which you would create.
56//
57// Here is an example (how to remove every second entry)
58//
59//  ---------------------------------------------------------------------
60//
61//    Int_t cnt = 0;
62//
63//    const Int_t num = event->GetNumPhotons();
64//    for (Int_t idx=0; idx<num; idx++)
65//    {
66//        if (idx%2==0)
67//           continue;
68//
69//        MPhotonData *dat = (*event)[idx];
70//
71//        (*event)[cnt++] = *dat;
72//     }
73//
74//     event->Shrink(cnt);
75//
76//  ---------------------------------- or -------------------------------
77//
78//    TClonesArray &arr = MPhotonEvent->GetArray();
79//
80//    Int_t cnt = 0;
81//
82//    const Int_t num = arr.GetEntriesFast();
83//    for (Int_t idx=0; idx<num; idx++)
84//    {
85//        if (idx%2==0)
86//           continue;
87//
88//        MPhotonData *dat = static_cast<MPhotonData*>(arr.UncheckedAt(idx));
89//
90//        *static_cast<MPhotonData*>(arr.UncheckedAt(cnt++)) = *dat;
91//     }
92//
93//     MPhotonEvent->Shrink(cnt);
94//
95//  ---------------------------------------------------------------------
96//
97// The flag for a sorted array is for speed reasons not in all conditions
98// maintained automatically. Especially Add() doesn't reset it.
99//
100// So be sure that if you want to sort your array it is really sorted.
101//
102//
103//   Version 1:
104//   ----------
105//    - First implementation
106//
107/////////////////////////////////////////////////////////////////////////////
108#include "MPhotonEvent.h"
109
110#include <fstream>
111#include <iostream>
112
113#include <TMarker.h>
114
115#include <MMath.h>
116
117#include "MArrayF.h"
118
119#include "MLog.h"
120#include "MLogManip.h"
121
122#include "MPhotonData.h"
123#include "MCorsikaFormat.h"
124
125ClassImp(MPhotonEvent);
126
127using namespace std;
128
129// ==========================================================================
130
131class MyClonesArray : public TClonesArray
132{
133public:
134    TObject **FirstRef() { return fCont; }
135
136    // --------------------------------------------------------------------------
137    //
138    // This is an extremly optimized version of ExpandCreateFast. It only resets
139    // the marker for the last element (fLast) to n-1 and doen't change anything
140    // else. This implicitly assumes that the stored objects don't allocate
141    // memory. It does not necessarily mean that the slots after fLast
142    // are empty (set to 0). This is what is assumed in the TClonesArray.
143    // We also don't call Changed() because it would reset Sorted. If the
144    // array was sorted before it is also sorted now. You MUST make sure
145    // that you only set n in a range for which valid entries have been
146    // created before (e.g. by ExpandCreateFast).
147    //
148    void FastShrink(Int_t n)
149    {
150        fLast = n - 1;
151    }
152
153    // --------------------------------------------------------------------------
154    //
155    // This is a optimized (faster) version of Delete which deletes consequtive
156    // entries from index idx1 to idx2 (both included) and calls their
157    // destructor. Note that there is no range checking done!
158    //
159    void FastRemove(Int_t idx1, Int_t idx2)
160    {
161        // Remove object at index idx.
162
163        //if (!BoundsOk("RemoveAt", idx1)) return 0;
164        //if (!BoundsOk("RemoveAt", idx2)) return 0;
165
166        Long_t dtoronly = TObject::GetDtorOnly();
167
168        idx1 -= fLowerBound;
169        idx2 -= fLowerBound;
170
171        for (TObject **obj=fCont+idx1; obj<=fCont+idx2; obj++)
172        {
173            if (!*obj)
174                continue;
175
176            if ((*obj)->TestBit(kNotDeleted)) {
177                // Tell custom operator delete() not to delete space when
178                // object fCont[i] is deleted. Only destructors are called
179                // for this object.
180                TObject::SetDtorOnly(*obj);
181                delete *obj;
182            }
183
184            *obj = 0;
185            // recalculate array size
186        }
187        TObject::SetDtorOnly((void*)dtoronly);
188
189        if (idx1<=fLast && fLast<=idx2)
190        {
191            do {
192                fLast--;
193            } while (fLast >= 0 && fCont[fLast] == 0);
194        }
195
196        Changed();
197    }
198
199
200    //void SetSorted() { fSorted = kTRUE; }
201
202    // --------------------------------------------------------------------------
203    //
204    // This is an optimized version of Sort which doesn't check the
205    // IsSortable flag before. It only sorts the entries from 0
206    // to GetEntriesFast().
207    //
208    void UncheckedSort()
209    {
210        if (fSorted)
211            return;
212
213        const Int_t nentries = GetEntriesFast();
214        if (nentries <= 0)
215            return;
216
217        QSort(GetObjectRef(First()), fKeep->GetObjectRef(fKeep->First()),
218              0, TMath::Min(nentries, kMaxInt-fLowerBound));
219
220        fSorted = kTRUE;
221    }
222};
223
224// ==========================================================================
225
226// --------------------------------------------------------------------------
227//
228// Default constructor. It initializes all arrays with zero size.
229//
230MPhotonEvent::MPhotonEvent(const char *name, const char *title)
231    : fData("MPhotonData", 1)
232{
233    fName  = name  ? name  : "MPhotonEvent";
234    fTitle = title ? title : "Corsika Event Data Information";
235
236    fData.SetBit(TClonesArray::kForgetBits);
237    fData.BypassStreamer(kFALSE);
238}
239
240// --------------------------------------------------------------------------
241//
242// This is an extremly optimized version of ExpandCreateFast. It only resets
243// the marker for the last element (fLast) to n-1 and doen't change anything
244// else. This has the advantage that the allocated memory is kept but only
245// valid entries are written to a file.
246//
247// If the array was sorted before it is also sorted now. You MUST make sure
248// that you only set n in a range for which valid entries have been
249// created before (e.g. by ExpandCreateFast).
250//
251Int_t MPhotonEvent::Shrink(Int_t n)
252{
253    /*
254    // The number of object written by the streamer is defined by
255    // GetEntriesFast(), i.e. this number must be shrinked to
256    // the real array size. We use RemoveAt instead of ExpandCreate
257    // because RemoveAt doesn't free memory. Thus in the next
258    // iteration it can be reused and doesn't need to be reallocated.
259    // Do not change this. It is optimized for execution speed
260    //        for (int i=fData.GetSize()-1; i>=n; i--)
261    //          fData.RemoveAt(i);
262    const Bool_t sorted = fData.IsSorted();
263
264    MyClonesArray &loc = static_cast<MyClonesArray&>(fData);
265
266    loc.FastRemove(n, fData.GetSize()-1);
267
268    // If it was sorted before it is also sorted now.
269    if (sorted)
270        loc.SetSorted();
271    */
272
273    // fData.ExpandCreateFast(n);  // Just set fLast = n -1
274
275    // Just set fLast = n -1
276    static_cast<MyClonesArray&>(fData).FastShrink(n);
277    return fData.GetEntriesFast();
278}
279
280// --------------------------------------------------------------------------
281//
282// The resized the array. If it has to be increased in size it is done
283// with ExpandCreateFast. If it should be shrinked it is done with
284// ExpandCreateFast if n>fData.GetSize()/100 or n==0. This keeps the allocated
285// memory and just sets the marker for the last element in the array (fLast).
286//
287// If the allocated memory grew to huge we reset the allocated memory
288// by calling ExpandCreate(n) (frees the allocated storage for the
289// objects) and Expand(n) (frees the allocated memory for the list
290// of pointers to the objects)
291//
292// 100 hundred is an arbitrary number taking into account that todays
293// computers have a lot of memory and we don't want to free and allocate
294// new memory too often.
295//
296// In priciple there might be more clever methods to handle the memory.
297//
298void MPhotonEvent::Resize(Int_t n)
299{
300    if (n==0 || n*100>fData.GetSize())
301        fData.ExpandCreateFast(n);  // Just set fLast = n -1
302    else
303    {
304        fData.ExpandCreate(n);      // Free memory of allocated MPhotonData
305        fData.Expand(n);            // Free memory of allocated pointers
306    }
307}
308
309// Overload the AsciiWrite method to store the informations of the photons onto disc.
310// For each corsika/pedestal... event there is a simple header written. After this
311// the informations for each photon are written in one line.
312Bool_t MPhotonEvent::AsciiWrite(ostream &out) const
313{
314    // Write a simple header. Be aware there will be written one header for each
315    // corsika/pedestal... event. Each photon of the event will be written in one
316    // line.
317    out << "#";
318    out << "photonID" << ",";
319    out << "tag" << ",";
320    out << "mirrorTag" << ",";
321    out << "posX" << ",";
322    out << "posY" << ",";
323    out << "cosU" << ",";
324    out << "cosV" << ",";
325    out << "time" << ",";
326    out << "weight" << ",";
327    out << "wavelength" << ",";
328    out << "productionHeight" << ",";
329    out << "primaryID" << endl;
330
331    // Get number of photons
332    const Int_t num = GetNumPhotons();
333
334    // Loop over all photons
335    for (Int_t i=0; i<num; i++)
336    {
337        // Get i-th photon
338        const MPhotonData &ph = operator[](i);
339
340        out << i << "," ;
341        out << ph.GetTag() << "," ;
342        out << ph.GetMirrorTag() << "," ;
343        out << ph.GetPosX() << "," ;
344        out << ph.GetPosY() << "," ;
345        out << ph.GetCosU() << "," ;
346        out << ph.GetCosV() << "," ;
347        out << ph.GetTime() << "," ;
348        out << ph.GetWeight() << "," ;
349        out << ph.GetWavelength() << "," ;
350        out << ph.GetProductionHeight() << "," ;
351        out << ph.GetPrimary() << endl ;
352    }
353    return kTRUE;
354}
355
356// --------------------------------------------------------------------------
357//
358// If n is smaller than the current allocated array size a reference to
359// the n-th entry is returned, otherwise an entry at n is created
360// calling the default constructor. Note, that there is no range check
361// but it is not recommended to call this function with
362// n>fData.GetSize()
363//
364MPhotonData &MPhotonEvent::Add(Int_t n)
365{
366    // Do not modify this. It is optimized for execution
367    // speed and flexibility!
368    TObject *o = 0;
369    if (n<fData.GetSize() && fData.UncheckedAt(n))
370    {
371        o=fData.UncheckedAt(n);
372        static_cast<MyClonesArray&>(fData).FastShrink(n+1);
373    }
374    else
375    {
376        o=fData.New(n);
377    }
378    return static_cast<MPhotonData&>(*o);
379}
380
381// --------------------------------------------------------------------------
382//
383// Add a new photon (MPhtonData) at the end of the array.
384// In this case the default constructor of MPhotonData is called.
385//
386// A reference to the new object is returned.
387//
388MPhotonData &MPhotonEvent::Add()
389{
390    return Add(GetNumPhotons());
391}
392
393void MPhotonEvent::Sort(Bool_t force)
394{
395    if (force)
396        fData.UnSort();
397
398    static_cast<MyClonesArray&>(fData).UncheckedSort(); /*Sort(GetEntriesFast())*/
399}
400
401// --------------------------------------------------------------------------
402//
403// Get the i-th photon from the array. Not, for speed reasons there is no
404// range check so you are responsible that you do not excess the number
405// of photons (GetNumPhotons)
406//
407MPhotonData &MPhotonEvent::operator[](UInt_t idx)
408{
409    return *static_cast<MPhotonData*>(fData.UncheckedAt(idx));
410}
411
412// --------------------------------------------------------------------------
413//
414// Get the i-th photon from the array. Not, for speed reasons there is no
415// range check so you are responsible that you do not excess the number
416// of photons (GetNumPhotons)
417//
418const MPhotonData &MPhotonEvent::operator[](UInt_t idx) const
419{
420    return *static_cast<MPhotonData*>(fData.UncheckedAt(idx));
421}
422
423// --------------------------------------------------------------------------
424//
425// Return a pointer to the first photon if available.
426//
427MPhotonData *MPhotonEvent::GetFirst() const
428{
429    return fData.GetEntriesFast()==0 ? 0 : static_cast<MPhotonData*>(fData.First());
430}
431
432// --------------------------------------------------------------------------
433//
434// Return a pointer to the last photon if available.
435//
436MPhotonData *MPhotonEvent::GetLast() const
437{
438    return fData.GetEntriesFast()==0 ? 0 : static_cast<MPhotonData*>(fData.Last());
439}
440
441// --------------------------------------------------------------------------
442//
443// Return the number of "external" photons, i.e. which are not NightSky
444//
445Int_t MPhotonEvent::GetNumExternal() const
446{
447    const Int_t n=GetNumPhotons();
448
449    Int_t rc = 0;
450    for (int i=0; i<n; i++)
451        if ((*this)[i].GetPrimary()!=MMcEvtBasic::kNightSky && (*this)[i].GetPrimary()!=MMcEvtBasic::kArtificial)
452            rc++;
453
454    return rc;
455}
456
457// --------------------------------------------------------------------------
458//
459// Return time of first photon, 0 if none in array.
460// Note: If you want this to be the earliest make sure that the array
461// is properly sorted.
462//
463Float_t MPhotonEvent::GetTimeFirst() const
464{
465    const MPhotonData *dat=GetFirst();
466    return dat ? dat->GetTime() : 0;
467}
468
469// --------------------------------------------------------------------------
470//
471// Return time of first photon, 0 if none in array.
472// Note: If you want this to be the latest make sure that the array
473// is properly sorted.
474//
475Float_t MPhotonEvent::GetTimeLast() const
476{
477    const MPhotonData *dat=GetLast();
478    return dat ? dat->GetTime() : 0;
479}
480
481// --------------------------------------------------------------------------
482//
483// Return the median devian from the median of all arrival times.
484// The median deviation is calculated using MMath::MedianDev.
485// It is the half width in which one sigma (~68%) of all times are
486// contained around the median.
487//
488Double_t MPhotonEvent::GetTimeMedianDev() const
489{
490    const UInt_t n = GetNumPhotons();
491
492    MArrayF arr(n);
493    for (UInt_t i=0; i<n; i++)
494        arr[i] = operator[](i).GetTime();
495
496    return MMath::MedianDev(n, arr.GetArray()/*, Double_t &med*/);
497}
498
499Double_t MPhotonEvent::GetMeanX() const
500{
501    const UInt_t n = GetNumPhotons();
502    if (n==0)
503        return 0;
504
505    Double_t mean = 0;
506    for (UInt_t i=0; i<n; i++)
507        mean += operator[](i).GetPosX();
508
509    return mean / n;
510}
511
512Double_t MPhotonEvent::GetMeanY() const
513{
514    const UInt_t n = GetNumPhotons();
515    if (n==0)
516        return 0;
517
518    Double_t mean = 0;
519    for (UInt_t i=0; i<n; i++)
520        mean += operator[](i).GetPosY();
521
522    return mean / n;
523}
524
525Double_t MPhotonEvent::GetMeanT() const
526{
527    const UInt_t n = GetNumPhotons();
528    if (n==0)
529        return 0;
530
531    Double_t mean = 0;
532    for (UInt_t i=0; i<n; i++)
533        mean += operator[](i).GetTime();
534
535    return mean / n;
536}
537
538void MPhotonEvent::AddXY(Double_t x, Double_t y)
539{
540    const UInt_t n = GetNumPhotons();
541
542    for (UInt_t i=0; i<n; i++)
543    {
544        MPhotonData &p = operator[](i);
545        p.SetPosition(p.GetPosX()+x, p.GetPosY()+y);
546    }
547}
548
549void MPhotonEvent::SimWavelength(Float_t wmin, Float_t wmax)
550{
551    const UInt_t n = GetNumPhotons();
552    if (n==0)
553        return;
554
555    // Corsika has already produced and written the wavelength
556    if (operator[](0).GetWavelength()>0)
557        return;
558
559    for (UInt_t i=0; i<n; i++)
560        operator[](i).SimWavelength(wmin, wmax);
561}
562
563Int_t MPhotonEvent::ReadEventIoEvtCompact(MCorsikaFormat *fInFormat)
564{
565   Int_t bunchHeader[3];
566   fInFormat->Read(bunchHeader, 3 * sizeof(Int_t));
567
568   Int_t n = 0;
569
570   for (int bunch = 0; bunch < bunchHeader[2]; bunch++)
571   {
572       Short_t buffer[8];
573       fInFormat->Read(buffer, 8 * sizeof(Short_t));
574
575       if (Add(n).FillEventIO(buffer))
576           n++;
577   }
578
579   Resize(n);
580   fData.UnSort();
581
582   SetReadyToSave();
583
584   //*fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl;
585   return kTRUE;
586
587}
588
589Int_t MPhotonEvent::ReadEventIoEvt(MCorsikaFormat *fInFormat)
590{
591   Int_t  bunchHeader[3];
592   fInFormat->Read(bunchHeader, 3 * sizeof(Int_t));
593
594   Int_t n = 0;
595
596   for (int bunch = 0; bunch < bunchHeader[2]; bunch++)
597      {
598      Float_t buffer[8];
599      fInFormat->Read(buffer, 8 * sizeof(Float_t));
600
601      if (Add(n).FillEventIO(buffer))
602         n++;
603      }
604
605   Resize(n);
606   fData.UnSort();
607
608   SetReadyToSave();
609
610   //*fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl;
611   return kTRUE;
612
613}
614
615Int_t MPhotonEvent::ReadCorsikaEvt(const vector<Float_t> &data, const uint32_t &blockLength, const Int_t &arrayIdx)
616{
617    if (data.size()%blockLength)
618    {
619        gLog << err << "ERROR - Size mismatch in photon bunch data." << endl;
620        return kERROR;   // Error occured
621    }
622
623    const uint32_t N = data.size()/blockLength;
624
625   Int_t n = 0;
626
627   for (const Float_t *ptr=data.data(); ptr<data.data()+data.size(); ptr+=blockLength)
628   {
629       const Int_t rc = blockLength==7 ?
630           Add(n).FillCorsika(ptr, arrayIdx) :
631           Add(n).FillCorsikaThin(ptr, arrayIdx);
632
633      switch (rc)
634      {
635      case kCONTINUE:  continue;        // No data in this bunch... skip it.
636      case kERROR:     return kERROR;   // Error occured
637      }
638
639      // This is a photon we would like to keep later.
640      // Increase the counter by one
641      n++;
642      }
643
644   Resize(n);
645   fData.UnSort();
646
647   SetReadyToSave();
648
649   //*fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl;
650   return kTRUE;
651
652}
653
654// --------------------------------------------------------------------------
655//
656// Print the array
657//
658void MPhotonEvent::Print(Option_t *) const
659{
660    // This is much faster than looping over all entries and discarding
661    // the empty ones
662    const UInt_t n = GetNumPhotons();
663    for (UInt_t i=0; i<n; i++)
664        operator[](i).Print();
665}
666
667// ------------------------------------------------------------------------
668//
669// You can call Draw() to add the photons to the current pad.
670// The photons are painted each tim ethe pad is updated.
671// Make sure that you use the right (world) coordinate system,
672// like created, eg. by the MHCamera histogram.
673//
674void MPhotonEvent::Paint(Option_t *)
675{
676    MPhotonData *ph=NULL;
677
678    TMarker m;
679    m.SetMarkerStyle(kFullDotMedium); // Gtypes.h
680
681    TIter Next(&fData);
682    while ((ph=(MPhotonData*)Next()))
683    {
684        m.SetX(ph->GetPosY()*10);  // north
685        m.SetY(ph->GetPosX()*10);  // east
686        m.Paint();
687    }
688}
Note: See TracBrowser for help on using the repository browser.