source: trunk/Mars/manalysis/MSoftwareTriggerCalc.cc @ 19311

Last change on this file since 19311 was 19311, checked in by tbretz, 6 months ago
Added some comments
File size: 5.9 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, 2016 <mailto:tbretz@physik.rwth-aachen.de>
19!
20!   Copyright: MAGIC Software Development, 2016
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27//   MSoftwareTriggerCalc
28//
29//////////////////////////////////////////////////////////////////////////////
30#include "MSoftwareTriggerCalc.h"
31
32#include <algorithm>
33
34#include "MLog.h"
35#include "MLogManip.h"
36
37#include "MParList.h"
38
39#include "MRawEvtData.h"
40#include "MSoftwareTrigger.h"
41#include "MPedestalSubtractedEvt.h"
42
43ClassImp(MSoftwareTriggerCalc);
44
45using namespace std;
46
47// --------------------------------------------------------------------------
48//
49// Default constructor.
50//
51MSoftwareTriggerCalc::MSoftwareTriggerCalc(const char *name, const char *title)
52    : fRawEvt(0), fSignal(0), fTrigger(0)
53{
54    fName  = name  ? name  : "MSoftwareTriggerCalc";
55    fTitle = title ? title : "Calculate the FACT trigger in software";
56}
57
58// --------------------------------------------------------------------------
59//
60Int_t MSoftwareTriggerCalc::PreProcess(MParList *pList)
61{
62    fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");//, AddSerialNumber(fNamePedestalSubtractedEvt));
63    if (!fRawEvt)
64    {
65        *fLog << err << "MRawEvtData not found... aborting." << endl;
66        return kFALSE;
67    }
68
69    fSignal = (MPedestalSubtractedEvt*)pList->FindObject("MPedestalSubtractedEvt");//, AddSerialNumber(fNamePedestalSubtractedEvt));
70    if (!fSignal)
71    {
72        *fLog << err << "MPedestalSubtractedEvt not found... aborting." << endl;
73        return kFALSE;
74    }
75
76    fTrigger = (MSoftwareTrigger*)pList->FindCreateObj("MSoftwareTrigger");
77    if (!fTrigger)
78        return kFALSE;
79
80    return kTRUE;
81}
82
83// --------------------------------------------------------------------------
84//
85Bool_t MSoftwareTriggerCalc::ReInit(MParList *pList)
86{
87    // FIXME: Check number of samples!
88    return kTRUE;
89}
90
91// --------------------------------------------------------------------------
92//
93Int_t MSoftwareTriggerCalc::Process()
94{
95    const UInt_t nch  = fRawEvt->GetNumPixels(); // HW Channels in file [1440]
96    const UInt_t npix = fSignal->GetNumPixels(); // Pixel in camera [1440 or 64]
97
98    const UShort_t *ids = fRawEvt->GetPixelIds();
99
100    // FACT: Exclude time marker, FAMOUS take the whole range
101    const int beg =  10;
102    const int end = npix==1440 ? std::min(UInt_t(200), fSignal->GetNumSamples()) : fSignal->GetNumSamples()-10;
103    const int num = end-beg;
104
105    MArrayF buf(160*num);
106
107    double avg = 0;
108
109    for (int hw=0; hw<nch; hw++)
110    {
111        // The FAMOUS camera has less pixels than
112        // what is stored as channels in the file
113        if (ids[hw]>=npix)
114            continue;
115
116        if (npix==1440) // FACT[1440] else FAMOUS[64]
117        {
118            // dead pixels
119            if (hw==927 || hw==80 || hw==873)                   // 3
120                continue;
121            // crazy pixels
122            if (hw==863 || hw==297 || hw==868)                  // 3
123                continue;
124            // broken DRS board
125            if (hw>=720 && hw<=728)                             // 9
126                continue;
127            // broken/suspicious bias channels
128            if ((hw>=171 && hw<=174) || (hw>=184 && hw<=188))   // 9
129                continue;
130        }
131
132        const Float_t *raw = fSignal->GetSamples(ids[hw]);
133
134        Float_t *sum = buf.GetArray()+(hw/9)*num;
135
136        // Sum all nine pixel of one trigger patch
137        for (const Float_t *ptr=raw+beg; ptr<raw+end; ptr++, sum++)
138            *sum += *ptr;
139    }
140
141    // apply correction factor for patches
142    // 927/9, 80/9, 873/9, 863/9, 297/9, 868/9
143/*
144    const int excl[] = { 8, 33, 95, 96, 97, 103, -1};
145    for (const int *e=excl; *e>=0; e++)
146    {
147        Float_t *raw = buf.GetArray() + *e * num;
148        for (Float_t *ptr=raw; ptr<raw+num; ptr++)
149            *ptr *= 1.125; // 9/8
150    }
151    {
152        Float_t *raw = buf.GetArray() + 19 * num;   // 171-174
153        for (Float_t *ptr=raw; ptr<raw+num; ptr++)
154            *ptr *= 1.8;   // 9/5 // 5 channels left
155    }
156    {
157        Float_t *raw = buf.GetArray() + 20 * num;   // 184-188
158        for (Float_t *ptr=raw; ptr<raw+num; ptr++)
159            *ptr *= 2.25;  // 9/4 // 4 channels left
160    }
161*/
162
163    Float_t  max = -50000;
164    UShort_t pos = 0;
165    Short_t  patch = -1;
166
167    for (int i=0; i<160; i++)
168    {
169        int idx = 0;
170        Float_t v[4] = { 0, 0, 0, 0 };
171
172        Float_t *sum=buf.GetArray() + num*i;
173        for (Float_t *ptr=sum+15; ptr<sum+num; ptr++)
174        {
175            ptr[0] -= 0.5*ptr[-15];
176
177            avg += ptr[0];
178            v[idx++%4] = ptr[0];
179
180            const Float_t min = *std::min_element(v, v+4);
181            if (min<=max)
182                continue;
183
184            max = min;
185            pos = ptr-sum+beg;
186            patch = i;
187        }
188    }
189
190    avg /= num-15;
191    avg /= npix==1440 ? 1440-24 : npix;
192
193    fTrigger->SetData(patch, avg, pos, max);
194    fTrigger->SetReadyToSave();
195
196    return kTRUE;
197}
198
199// --------------------------------------------------------------------------
200//
201Int_t MSoftwareTriggerCalc::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
202{
203    return kTRUE;
204}
205
Note: See TracBrowser for help on using the repository browser.