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

Last change on this file since 18549 was 18549, checked in by tbretz, 3 years ago
The wavelength is now kept as a negative value so that it can be checked on the photon level if CEFFIC was turned on. It is however returned without sign for convenience. The loop to simulate the wavelength is only executed if the first photon has a zero wavelength -- it is assumed they are all consistent.
File size: 19.3 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(Float_t * data, Int_t numEvents, Int_t arrayIdx)
616{
617   Int_t n = 0;
618
619   for (Int_t event = 0; event < numEvents; event++)
620      {
621      const Int_t rc = Add(n).FillCorsika(data + 7 * event, arrayIdx);
622
623      switch (rc)
624      {
625      case kCONTINUE:  continue;        // No data in this bunch... skip it.
626      case kERROR:     return kERROR;   // Error occured
627      }
628
629      // This is a photon we would like to keep later.
630      // Increase the counter by one
631      n++;
632      }
633
634   Resize(n);
635   fData.UnSort();
636
637   SetReadyToSave();
638
639   //*fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl;
640   return kTRUE;
641
642}
643
644// --------------------------------------------------------------------------
645//
646// Print the array
647//
648void MPhotonEvent::Print(Option_t *) const
649{
650    // This is much faster than looping over all entries and discarding
651    // the empty ones
652    const UInt_t n = GetNumPhotons();
653    for (UInt_t i=0; i<n; i++)
654        operator[](i).Print();
655}
656
657// ------------------------------------------------------------------------
658//
659// You can call Draw() to add the photons to the current pad.
660// The photons are painted each tim ethe pad is updated.
661// Make sure that you use the right (world) coordinate system,
662// like created, eg. by the MHCamera histogram.
663//
664void MPhotonEvent::Paint(Option_t *)
665{
666    MPhotonData *ph=NULL;
667
668    TMarker m;
669    m.SetMarkerStyle(kFullDotMedium); // Gtypes.h
670
671    TIter Next(&fData);
672    while ((ph=(MPhotonData*)Next()))
673    {
674        m.SetX(ph->GetPosY()*10);  // north
675        m.SetY(ph->GetPosX()*10);  // east
676        m.Paint();
677    }
678}
Note: See TracBrowser for help on using the repository browser.