Index: trunk/MagicSoft/Mars/mmain/MCameraDisplay.cc =================================================================== --- trunk/MagicSoft/Mars/mmain/MCameraDisplay.cc (revision 6855) +++ trunk/MagicSoft/Mars/mmain/MCameraDisplay.cc (revision 6856) @@ -57,4 +57,15 @@ fTop2->AddFrame(disp, laybut); +/* + const TGPicture *pic2 = fList->GetPicture("marslogo.xpm"); + if (pic2) + { + TGPictureButton *mars = new TGPictureButton(fTop2, pic2, -1); + fList->Add(mars); + TGLayoutHints *lay2 = new TGLayoutHints(kLHintsRight, 10, 10, 10, 10); + fList->Add(lay2); + fTop2->AddFrame(mars, lay2); + } + */ MapSubwindows(); Index: trunk/MagicSoft/Mars/mmain/MEventDisplay.cc =================================================================== --- trunk/MagicSoft/Mars/mmain/MEventDisplay.cc (revision 6855) +++ trunk/MagicSoft/Mars/mmain/MEventDisplay.cc (revision 6856) @@ -70,5 +70,5 @@ #include "MMcPedestalNSBAdd.h" // MMcPedestalNSBAdd -#include "MCerPhotAnal2.h" // MPedestalCalc +#include "MSignalCalc.h" // MPedestalCalc #include "MTaskEnv.h" // MTaskEnv #include "MExtractTimeAndChargeDigitalFilter.h" @@ -159,5 +159,5 @@ return -1; - if (t->FindBranch("MCerPhotEvt.")) + if (t->FindBranch("MSignalCam.")) return 1; @@ -253,8 +253,8 @@ taskenv1->SetOwner(); - MCerPhotAnal2 *nanal = new MCerPhotAnal2; - MFillH *fill01 = new MFillH(evt01, "MCerPhotEvt", "MFillH01"); + MSignalCalc *nanal = new MSignalCalc; + MFillH *fill01 = new MFillH(evt01, "MSignalCam", "MFillH01"); MImgCleanStd *clean = new MImgCleanStd; - MFillH *fill02 = new MFillH(evt02, "MCerPhotEvt", "MFillH02"); + MFillH *fill02 = new MFillH(evt02, "MSignalCam", "MFillH02"); MFillH *fill03 = new MFillH(evt03, type==1?"MPedPhotFundamental":"MPedestalCam", "MFillH03"); MFillH *fill04 = new MFillH(evt04, type==1?"MPedPhotFromExtractorRndm":"MPedestalCam", "MFillH04"); @@ -264,5 +264,5 @@ MMcTriggerLvl2Calc *trcal = new MMcTriggerLvl2Calc; MFillH *fill09 = new MFillH(evt09, "MMcTriggerLvl2", "MFillH09"); - MFillH *fill10 = new MFillH(evt10, "MCerPhotEvt", "MFillH10"); + MFillH *fill10 = new MFillH(evt10, "MSignalCam", "MFillH10"); MBadPixelsCalc *bcalc = new MBadPixelsCalc; @@ -306,5 +306,5 @@ tlist->AddToList(taskenv1); // Calculates MExtractedSignalCam, MArrivalTimeCam tlist->AddToList(mcupd); - tlist->AddToList(calib); // MExtractedSignalCam --> MCerPhotEvt + tlist->AddToList(calib); // MExtractedSignalCam --> MSignalCam tlist->AddToList(bcalc); // Produce MBadPixelsCam @@ -332,5 +332,5 @@ if (type==1) { - MFillH *fill08 = new MFillH(evt08, "MArrivalTime", "MFillH8"); + MFillH *fill08 = new MFillH(evt08, "MSignalCam", "MFillH8"); tlist->AddToList(fill08); } Index: trunk/MagicSoft/Mars/mmain/MainLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mmain/MainLinkDef.h (revision 6855) +++ trunk/MagicSoft/Mars/mmain/MainLinkDef.h (revision 6856) @@ -7,15 +7,6 @@ #pragma link C++ class MBrowser+; -#pragma link C++ class MMars+; -#pragma link C++ class MAnalysis+; -#pragma link C++ class MDataCheck+; -#pragma link C++ class MMonteCarlo+; - #pragma link C++ class MCameraDisplay+; #pragma link C++ class MEventDisplay+; -#pragma link C++ class MOnlineDisplay+; -#pragma link C++ class MOnlineDump+; - -#pragma link C++ class MGDisplayAdc+; #endif Index: trunk/MagicSoft/Mars/mmain/Makefile =================================================================== --- trunk/MagicSoft/Mars/mmain/Makefile (revision 6855) +++ trunk/MagicSoft/Mars/mmain/Makefile (revision 6856) @@ -26,12 +26,5 @@ SRCFILES = MBrowser.cc \ - MDataCheck.cc \ - MMars.cc \ - MAnalysis.cc \ - MMonteCarlo.cc \ - MGDisplayAdc.cc \ - MEventDisplay.cc \ - MOnlineDisplay.cc \ - MOnlineDump.cc \ + MEventDisplay.cc \ MCameraDisplay.cc Index: trunk/MagicSoft/Mars/msignal/MArrivalTimeCam.cc =================================================================== --- trunk/MagicSoft/Mars/msignal/MArrivalTimeCam.cc (revision 6855) +++ trunk/MagicSoft/Mars/msignal/MArrivalTimeCam.cc (revision 6856) @@ -153,4 +153,5 @@ break; case 4: + case 6: if ((*this)[idx].GetNumHiGainSaturated() > 0) val = (*this)[idx].GetArrivalTimeLoGain(); Index: trunk/MagicSoft/Mars/msignal/MSignalCalc.cc =================================================================== --- trunk/MagicSoft/Mars/msignal/MSignalCalc.cc (revision 6856) +++ trunk/MagicSoft/Mars/msignal/MSignalCalc.cc (revision 6856) @@ -0,0 +1,293 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Thomas Bretz 12/2000 +! +! Copyright: MAGIC Software Development, 2000-2001 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MSignalCalc +// +// This is a task which calculates the number of photons from the FADC +// time slices. At the moment it integrates simply the FADC values. +// +// Input Containers: +// MRawEvtData, MPedestalCam +// +// Output Containers: +// MSignalCam +// +////////////////////////////////////////////////////////////////////////////// + +#include "MSignalCalc.h" + +#include "MParList.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MRawRunHeader.h" +#include "MRawEvtData.h" // MRawEvtData::GetNumPixels +#include "MSignalCam.h" +#include "MPedestalPix.h" +#include "MPedestalCam.h" +#include "MRawEvtPixelIter.h" + +ClassImp(MSignalCalc); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor. b is the number of slices before the maximum slice, +// a the number of slices behind the maximum slice which is taken as signal. +// +MSignalCalc::MSignalCalc(Byte_t b, Byte_t a, const char *name, const char *title) + : fBefore(b), fAfter(a) +{ + fName = name ? name : "MSignalCalc"; + fTitle = title ? title : "Task to calculate Cerenkov photons from raw data"; + + AddToBranchList("MRawEvtData.fHiGainPixId"); + AddToBranchList("MRawEvtData.fLoGainPixId"); + AddToBranchList("MRawEvtData.fHiGainFadcSamples"); + AddToBranchList("MRawEvtData.fLoGainFadcSamples"); +} + +// -------------------------------------------------------------------------- +// +// The PreProcess searches for the following input containers: +// - MRawRunHeader +// - MRawEvtData +// - MPestalCam +// +// The following output containers are also searched and created if +// they were not found: +// - MSignalCam +// +Int_t MSignalCalc::PreProcess(MParList *pList) +{ + fSkip = 0; + + fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); + if (!fRunHeader) + { + *fLog << dbginf << "MRawRunHeader not found... aborting." << endl; + return kFALSE; + } + + fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData"); + if (!fRawEvt) + { + *fLog << dbginf << "MRawEvtData not found... aborting." << endl; + return kFALSE; + } + + fCerPhotEvt = (MSignalCam*)pList->FindCreateObj("MSignalCam"); + if (!fCerPhotEvt) + return kFALSE; + + return kTRUE; +} + +Bool_t MSignalCalc::ReInit(MParList *pList) +{ + fPedestals=NULL; + + // This must be done in ReInit because in PreProcess the + // headers are not available + const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); + if (!runheader) + *fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl; + else + { + if (runheader->IsMonteCarloRun()) + return kTRUE; + } + + fPedestals = (MPedestalCam*)pList->FindCreateObj("MPedestalCam"); + if (runheader && !fPedestals) + return kFALSE; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Calculate the integral of the FADC time slices and store them as a new +// pixel in the MSignalCam container. +// +Int_t MSignalCalc::Process() +{ + MRawEvtPixelIter pixel(fRawEvt); + + while (pixel.Next()) + { + Byte_t *ptr = pixel.GetHiGainSamples(); + Byte_t *max = ptr+pixel.GetIdxMaxHiGainSample(); + Byte_t *end = ptr+fRawEvt->GetNumHiGainSamples(); + Byte_t *first = max-fBefore; + Byte_t *last = max+fAfter; + + ULong_t sumb = 0; // sum background + ULong_t sqb = 0; // sum sqares background + ULong_t sumsb = 0; // sum signal+background + ULong_t sqsb = 0; // sum sqares signal+background + + Int_t sat = 0; // saturates? + Int_t ishi = 0; // has a high content? + Int_t nb = 0; + Int_t nsb = 0; + + if (*max==255) + sat++; + + if (*max>80) + ishi++; + + while (ptrlast) + { + sumb += *ptr; + sqb += *ptr* *ptr; + nb++; + } + else + { + sumsb += *ptr; + sqsb += *ptr* *ptr; + nsb++; + } + ptr++; + } + + if (sat==0 && ishi) + { + // Area: x9 + ptr = pixel.GetLoGainSamples(); + end = ptr+fRawEvt->GetNumLoGainSamples(); + + sumb = 0; // sum background + sqb = 0; // sum sqares background + nb = 0; + + while (ptr1 && !ishi) + { + // Area: x9 + ptr = pixel.GetLoGainSamples(); + max = ptr+pixel.GetIdxMaxLoGainSample(); + + if (*max>250) + { + fSkip++; + return kCONTINUE; + } + + end = ptr+fRawEvt->GetNumLoGainSamples(); + first = max-fBefore; + last = max+fAfter; + + sumsb = 0; // sum signal+background + sqsb = 0; // sum sqares signal+background + //sumb = 0; // sum background + //sqb = 0; // sum sqares background + + //nb = 0; + nsb = 0; + while (ptrlast) + { + /* + // Background already calced from hi-gains! + sumb += ptr[i]; + sqb += ptr[i]*ptr[i]; + nb++;*/ + } + else + { + sumsb += *ptr; + sqsb += *ptr* *ptr; + nsb++; + } + ptr++; + } + } + + Float_t b = (float)sumb/nb; // background + Float_t sb = (float)sumsb/nsb; // signal+background + + Float_t msb = (float)sqb/nb; // mean square background + //Float_t mssb = (float)sqsb/nsb; // mean square signal+background + + Float_t sigb = sqrt(msb-b*b); // sigma background + //Float_t sigsb = sqrt(mssb-sb*sb); // sigma signal+background + + Float_t s = sb-b; // signal + Float_t sqs = sqsb-nsb*b; // sum sqaures signal + + Float_t mss = (float)sqs/nsb; // mean quare signal + //Float_t sigs = sqrt(mss-s*s); // sigma signal + + if (sat>1) + s*=10; // tgb has measured 9, but Florian said it's 10. + + Int_t idx = pixel.GetPixelId(); +// fCerPhotEvt->AddPixel(idx, s, sigs); + + // Preliminary: Do not overwrite pedestals calculated by + // MMcPedestalCopy and MMcPedestalNSBAdd + if (fPedestals) + (*fPedestals)[idx].Set(b, sigb); + } + +// fCerPhotEvt->FixSize(); +// fCerPhotEvt->SetReadyToSave(); + + if (fPedestals) + fPedestals->SetReadyToSave(); + + return kTRUE; +} + +Int_t MSignalCalc::PostProcess() +{ + if (GetNumExecutions()==0 || fSkip==0) + return kTRUE; + + *fLog << inf << endl; + *fLog << GetDescriptor() << " execution statistics:" << endl; + *fLog << dec << setfill(' '); + *fLog << " " << setw(7) << fSkip << " (" << setw(3) << (int)(fSkip*100/GetNumExecutions()) << "%) Evts skipped due to: lo gain saturated." << endl; + + return kTRUE; +} Index: trunk/MagicSoft/Mars/msignal/MSignalCalc.h =================================================================== --- trunk/MagicSoft/Mars/msignal/MSignalCalc.h (revision 6856) +++ trunk/MagicSoft/Mars/msignal/MSignalCalc.h (revision 6856) @@ -0,0 +1,37 @@ +#ifndef MARS_MSignalCalc +#define MARS_MSignalCalc + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class MRawEvtData; +class MPedestalCam; +class MSignalCam; +class MRawRunHeader; + +class MSignalCalc : public MTask +{ +private: + MPedestalCam *fPedestals; // Pedestals of all pixels in the camera + MRawEvtData *fRawEvt; // raw event data (time slices) + MSignalCam *fCerPhotEvt; // Cerenkov Photon Event used for calculation + MRawRunHeader *fRunHeader; // RunHeader information + + Byte_t fBefore; + Byte_t fAfter; + + Int_t fSkip; //! + + Int_t PreProcess(MParList *pList); + Bool_t ReInit(MParList *pList); + Int_t Process(); + Int_t PostProcess(); + +public: + MSignalCalc(Byte_t b=2, Byte_t a=5, const char *name=NULL, const char *title=NULL); + + ClassDef(MSignalCalc, 0) // Task to calculate cerenkov photons from raw data +}; + +#endif Index: trunk/MagicSoft/Mars/msignal/MSignalCam.cc =================================================================== --- trunk/MagicSoft/Mars/msignal/MSignalCam.cc (revision 6856) +++ trunk/MagicSoft/Mars/msignal/MSignalCam.cc (revision 6856) @@ -0,0 +1,662 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Thomas Bretz, 03/2005 +! +! Copyright: MAGIC Software Development, 2000-2005 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MSignalCam +// +// Class Version 1: +// ---------------- +// - first version +// +///////////////////////////////////////////////////////////////////////////// +#include "MSignalCam.h" + +#include +#include +#include + +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MGeomCam.h" +#include "MGeomPix.h" + +ClassImp(MSignalCam); +ClassImp(MSignalCamIter); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Creates a MSignalPix object for each pixel in the event +// +MSignalCam::MSignalCam(const char *name, const char *title) /*: fNumPixels(0)*/ +{ + fName = name ? name : "MSignalCam"; + fTitle = title ? title : "(Number of Photon)-Event Information"; + + fPixels = new TClonesArray("MSignalPix", 0); + + Reset(); +} + +// -------------------------------------------------------------------------- +// +// This is not yet implemented like it should. +// +/* +void MSignalCam::Draw(Option_t* option) +{ + // + // FIXME!!! Here the Draw function of the CamDisplay + // should be called to add the CamDisplay to the Pad. + // The drawing should be done in MCamDisplay::Paint + // + + // MGeomCam *geom = fType ? new MGeomCamMagic : new MGeomCamCT1; + // MCamDisplay *disp = new MCamDisplay(geom); + // delete geom; + // disp->DrawPhotNum(this); +} +*/ + +// -------------------------------------------------------------------------- +// +// reset counter and delete netries in list. +// +void MSignalCam::Reset() +{ + //fNumPixels = 0; + fNumSinglePixels = 0; + fSizeSinglePixels = 0; + fSizeSubIslands = 0; + fSizeMainIsland = 0; + //fMaxIndex = -1; + fNumIslands = -1; + //fLut.Set(0); + + fNumPixelsSaturatedHiGain = -1; + fNumPixelsSaturatedLoGain = -1; + + fPixels->ForEach(TObject, Clear)(); +} +/* +void MSignalCam::FixSize() +{ + fLut.Set(fMaxIndex+1); + + if (fPixels->GetEntriesFast() == (Int_t)fNumPixels) + return; + + fPixels->ExpandCreateFast(fNumPixels); +}*/ + +// -------------------------------------------------------------------------- +// +// Dump the cerenkov photon event to *fLog +// +void MSignalCam::Print(Option_t *) const +{ + const Int_t entries = fPixels->GetEntries(); + + *fLog << GetDescriptor() << dec << endl; + *fLog << " Number of Pixels: " << GetNumPixels() << "(" << entries << ")" << endl; + + for (Int_t i=0; iIsPixelUsed() : kFALSE; +} + +// -------------------------------------------------------------------------- +// +// Checks if in the pixel list is an entry with pixel id +// +Bool_t MSignalCam::IsPixelCore(Int_t id) const +{ + const MSignalPix *pix = GetPixById(id); + return pix ? pix->IsPixelCore() : kFALSE; +} + +// -------------------------------------------------------------------------- +// +// get the minimum number of photons of all valid pixels in the list +// If you specify a geometry the number of photons is weighted with the +// area of the pixel +// +Float_t MSignalCam::GetNumPhotonsMin(const MGeomCam *geom) const +{ + if (GetNumPixels() <= 0) + return -5.; + + const UInt_t n = GetNumPixels(); + + Float_t minval = FLT_MAX; + + for (UInt_t i=0; i=n) + // continue; + + Float_t testval = pix.GetNumPhotons(); + + if (geom) + testval *= geom->GetPixRatio(i); + + if (testval < minval) + minval = testval; + } + + return minval; +} + +// -------------------------------------------------------------------------- +// +// get the maximum number of photons of all valid pixels in the list +// If you specify a geometry the number of photons is weighted with the +// area of the pixel +// +Float_t MSignalCam::GetNumPhotonsMax(const MGeomCam *geom) const +{ + if (GetNumPixels() <= 0) + return 50.; + + const UInt_t n = GetNumPixels(); + + Float_t maxval = -FLT_MAX; + + for (UInt_t i=0; i=n) + // continue; + + Float_t testval = pix.GetNumPhotons(); + if (geom) + testval *= geom->GetPixRatio(i); + + if (testval > maxval) + maxval = testval; + } + return maxval; +} + +// -------------------------------------------------------------------------- +// +// get the minimum ratio of photons/error +// +Float_t MSignalCam::GetRatioMin(const MGeomCam *geom) const +{ + if (GetNumPixels() <= 0) + return -5.; + + Float_t minval = FLT_MAX; + + for (UInt_t i=0; iGetPixRatio(i/*pix.GetPixId()*/)); + + if (testval < minval) + minval = testval; + } + + return minval; +} + +// -------------------------------------------------------------------------- +// +// get the maximum ratio of photons/error +// +Float_t MSignalCam::GetRatioMax(const MGeomCam *geom) const +{ + if (GetNumPixels() <= 0) + return -5.; + + Float_t maxval = -FLT_MAX; + + for (UInt_t i=0; iGetPixRatio(i/*pix.GetPixId()*/)); + + if (testval > maxval) + maxval = testval; + } + + return maxval; +} + +// -------------------------------------------------------------------------- +// +// get the minimum of error +// If you specify a geometry the number of photons is weighted with the +// area of the pixel +// +Float_t MSignalCam::GetErrorPhotMin(const MGeomCam *geom) const +{ + if (GetNumPixels() <= 0) + return 50.; + + Float_t minval = FLT_MAX; + + for (UInt_t i=0; iGetPixRatio(i/*pix.GetPixId()*/); + + if (testval < minval) + minval = testval; + } + return minval; +} + +// -------------------------------------------------------------------------- +// +// get the maximum ratio of photons/error +// If you specify a geometry the number of photons is weighted with the +// area of the pixel +// +Float_t MSignalCam::GetErrorPhotMax(const MGeomCam *geom) const +{ + if (GetNumPixels() <= 0) + return 50.; + + Float_t maxval = -FLT_MAX; + + for (UInt_t i=0; iGetPixRatio(i/*pix.GetPixId()*/); + + if (testval > maxval) + maxval = testval; + } + return maxval; +} +/* +void MSignalCam::RemoveUnusedPixels() +{ + // Create iterator + TIter Next(fPixels); + MSignalPix *pix = NULL; + + fMaxIndex = -1; + + // Remove all unused pixels from list, calculate new fMaxIndex + while ((pix=(MSignalPix*)Next())) + { + if (!pix->IsPixelUsed()) + fPixels->Remove(pix); + else + fMaxIndex = TMath::Max(fMaxIndex, pix->GetPixId()); + } + + // Crompress array + fPixels->Compress(); + + // Get new number of entries in array + fNumPixels=fPixels->GetEntriesFast(); + + // Rebuild lookup table + RebuildLut(); +} +*/ +// -------------------------------------------------------------------------- +// +// Return a pointer to the pixel with the requested idx. NULL if it doesn't +// exist. The Look-up-table fLut is used. If its size is zero (according +// to Rene this will happen if an old class object is loaded) we still +// try to search in the array. +// +MSignalPix *MSignalCam::GetPixById(Int_t idx) const +{ + return (MSignalPix*)(fPixels->UncheckedAt(idx)); + /* + if (idx<0) + return 0; + + if (fLut.GetSize()>0) + { + if (idx>=fLut.GetSize()) + return 0; + return fLut[idx]<0 ? 0 : (MSignalPix*)(fPixels->UncheckedAt(fLut[idx])); + } + + TIter Next(fPixels); + MSignalPix *pix = NULL; + + while ((pix=(MSignalPix*)Next())) + if (pix->GetPixId()==idx) + return pix; + + return NULL;*/ +} + +MSignalPix *MSignalCam::AddPixel(Int_t idx, Float_t nph, Float_t er) +{ + /* + // + // If this is too slow or takes to much space we might use + // MGeomApply and an InitSize member function instead. + // + if (idx>=fLut.GetSize()) + { + const Int_t n = fLut.GetSize(); + fLut.Set(idx*2+1); //idx+1 is slower than idx*2+1 + for (int i=n; ifMaxIndex) + fMaxIndex=idx; + + return new ((*fPixels)[fNumPixels++]) MSignalPix(idx, nph, er); + */ + //return new ((*fPixels)[idx]) MSignalPix(nph, er); + + MSignalPix *pix = static_cast((*fPixels)[idx]); + pix->Set(nph, er); + return pix; +} + +// -------------------------------------------------------------------------- +// +// This function recursively finds all pixels of one island and assigns +// the number num as island number to the pixel. +// +// 1) Check whether a pixel with the index idx exists, is unused +// and has not yet a island number assigned. +// 2) Assign the island number num to the pixel +// 3) Loop over all its neighbors taken from the geometry geom. For all +// neighbors recursively call this function (CalcIsland) +// 4) Sum the size of the pixel and all neighbors newly assigned +// (by CalcIsland) to this island +// +// Returns the sum of the pixel size. +// +Double_t MSignalCam::CalcIsland(const MGeomCam &geom, Int_t idx, Int_t num) +{ + // Try to get the pixel information of a pixel with this index + MSignalPix *pix = GetPixById(idx); + + // If a pixel with this index is not existing... do nothing. + if (!pix) + return 0; + + // If an island number was already assigned to this pixel... do nothing. + if (pix->GetIdxIsland()>=0) + return 0; + + // If the pixel is an unused pixel... do nothing. + if (!pix->IsPixelUsed()) + return 0; + + // Assign the new island number num to this used pixel + pix->SetIdxIsland(num); + + // Get the geometry information (neighbors) of this pixel + const MGeomPix &gpix = geom[idx]; + + // Get the size of this pixel + Double_t size = pix->GetNumPhotons(); + + // Now do the same with all its neighbors and sum the + // sizes which they correspond to + const Int_t n = gpix.GetNumNeighbors(); + for (int i=0; iGetIdxIsland. The total number of islands available +// can be accessed with MSignalCam->GetNumIslands. +// +// CalcIslands returns the number of islands found. If an error occurs, +// eg the geometry has less pixels than the highest index stored, -1 is +// returned. +// +Int_t MSignalCam::CalcIslands(const MGeomCam &geom) +{ + const UInt_t numpix = GetNumPixels(); + + if (/*fMaxIndex<0 ||*/ numpix==0) + { + *fLog << err << "ERROR - MSignalCam doesn't contain pixels!" << endl; + fNumIslands = 0; + return -1; + } +/* + if ((UInt_t)fMaxIndex>=geom.GetNumPixels()) + { + *fLog << err << "ERROR - MSignalCam::CalcIslands: Size mismatch - geometry too small!" << endl; + return -1; + } + */ + // Create a list to hold the sizes of the islands (The maximum + // number of islands possible is roughly fNumPixels/4) + TArrayD size(numpix/3); + + // Calculate Islands + Int_t n=0; + Float_t totsize = 0; + + for (UInt_t idx=0; idxGetPixId()*/, n); + size[n++] = sz; + totsize += sz; + } + } + + // Create an array holding the indices + TArrayI idxarr(n); + + // Sort the sizes descending + TMath::Sort(n, size.GetArray(), idxarr.GetArray(), kTRUE); + + // Replace island numbers by size indices -- After this + // islands indices are sorted by the island size + for (UInt_t idx=0; idx0 ? totsize-size[idxarr[0]] : 0; + fSizeMainIsland = n>0 ? size[idxarr[0]] : 0; + + // return number of island + return fNumIslands; +} + +// -------------------------------------------------------------------------- +// +// Returns, depending on the type flag: +// +// 0: Number of Photons*PixRatio +// 1: Error*sqrt(PixRatio) +// 2: Cleaning level = Num Photons*sqrt(PixRatio)/Error +// 3: Number of Photons +// 4: Error +// 5: Island index +// +Bool_t MSignalCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const +{ + MSignalPix *pix = GetPixById(idx); + + if ((!pix || !pix->IsPixelUsed()) && type<6) + return kFALSE; + + const Double_t ratio = cam.GetPixRatio(idx); + + switch (type) + { + case 1: + val = pix->GetErrorPhot()*TMath::Sqrt(ratio); + return kTRUE; + case 2: + if (pix->GetErrorPhot()<=0) + return kFALSE; + val = pix->GetNumPhotons()*TMath::Sqrt(ratio)/pix->GetErrorPhot(); + return kTRUE; + case 3: + val = pix->GetNumPhotons(); + break; + case 4: + val = pix->GetErrorPhot(); + break; + case 5: + val = pix->GetIdxIsland(); + break; + case 6: + val = pix->GetArrivalTime(); + break; + case 7: + if (pix->GetNumPhotons()*ratio<15) // maybe 15 as static data member? + return kFALSE; + val = pix->GetArrivalTime(); + break; + default: + val = pix->GetNumPhotons()*ratio; + return kTRUE; + } + return kTRUE; +} + +void MSignalCam::DrawPixelContent(Int_t num) const +{ + *fLog << warn << "MSignalCam::DrawPixelContent - not available." << endl; +} + +TObject *MSignalCamIter::Next() +{ + if (!fUsedOnly) + { + fIdx++; + return TObjArrayIter::Next(); + } + + MSignalPix *pix; + while ((pix = (MSignalPix*)TObjArrayIter::Next())) + { + fIdx++; + if (pix->IsPixelUsed()) + return pix; + } + return pix; +} Index: trunk/MagicSoft/Mars/msignal/MSignalCam.h =================================================================== --- trunk/MagicSoft/Mars/msignal/MSignalCam.h (revision 6856) +++ trunk/MagicSoft/Mars/msignal/MSignalCam.h (revision 6856) @@ -0,0 +1,119 @@ +#ifndef MARS_MSignalCam +#define MARS_MSignalCam + +#ifndef ROOT_TClonesArray +#include +#endif +#ifndef ROOT_TArrayI +#include +#endif +#ifndef MARS_MCamEvent +#include "MCamEvent.h" +#endif +#ifndef MARS_MSignalPix +#include "MSignalPix.h" +#endif + +class MGeomCam; +class MSignalPix; +class MSignalCamIter; + +class MSignalCam : public MParContainer, public MCamEvent +{ + friend class MSignalCamIter; +private: + Short_t fNumIslands; + Short_t fNumSinglePixels; + Float_t fSizeSinglePixels; + Float_t fSizeSubIslands; + Float_t fSizeMainIsland; + + Int_t fNumPixelsSaturatedHiGain; + Int_t fNumPixelsSaturatedLoGain; + TClonesArray *fPixels; //-> FIXME: Change TClonesArray away from a pointer? + + Double_t CalcIsland(const MGeomCam &geom, Int_t idx, Int_t num); + +public: + MSignalCam(const char *name=NULL, const char *title=NULL); + ~MSignalCam() { delete fPixels; } + + // Setter function to fill pixels + MSignalPix *AddPixel(Int_t idx, Float_t nph=0, Float_t er=0); + void InitSize(const UInt_t i) + { + fPixels->ExpandCreate(i); + } + + // Getter functions + UInt_t GetNumPixels() const { return fPixels->GetEntriesFast(); } + Short_t GetNumIslands() const { return fNumIslands; }; + Short_t GetNumSinglePixels() const { return fNumSinglePixels; } + Float_t GetSizeSinglePixels() const { return fSizeSinglePixels; } + Float_t GetSizeSubIslands() const { return fSizeSubIslands; } + Float_t GetSizeMainIsland() const { return fSizeMainIsland; } + Int_t GetNumPixelsSaturatedHiGain() const { return fNumPixelsSaturatedHiGain; } + Int_t GetNumPixelsSaturatedLoGain() const { return fNumPixelsSaturatedLoGain; } + + // Setter functions for use in image cleaning classes only + void SetSinglePixels(Short_t num, Float_t size) { fNumSinglePixels=num; fSizeSinglePixels=size; } + void SetNumPixelsSaturated(UInt_t hi, UInt_t lo) { fNumPixelsSaturatedHiGain=hi;fNumPixelsSaturatedLoGain=lo; } + + Bool_t IsPixelExisting(Int_t id) const; + Bool_t IsPixelUsed (Int_t id) const; + Bool_t IsPixelCore (Int_t id) const; + + Float_t GetNumPhotonsMin(const MGeomCam *geom=NULL) const; + Float_t GetNumPhotonsMax(const MGeomCam *geom=NULL) const; + + Float_t GetRatioMin(const MGeomCam *geom=NULL) const; + Float_t GetRatioMax(const MGeomCam *geom=NULL) const; + + Float_t GetErrorPhotMin(const MGeomCam *geom=NULL) const; + Float_t GetErrorPhotMax(const MGeomCam *geom=NULL) const; + + // Getter functions to access single pixels + MSignalPix &operator[](int i) { return *(MSignalPix*)(fPixels->UncheckedAt(i)); } + MSignalPix &operator[](int i) const { return *(MSignalPix*)(fPixels->UncheckedAt(i)); } + + MSignalPix *GetPixById(Int_t idx) const; + + // Functions to change the contained data + void Scale(Double_t f) { fPixels->ForEach(MSignalPix, Scale)(f); } + Int_t CalcIslands(const MGeomCam &geom); + + // class MParContainer + void Reset(); + + // class TObject + void Print(Option_t *opt=NULL) const; + void Clear(Option_t *opt=NULL) { Reset(); } + + // class MCamEvent + Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const; + void DrawPixelContent(Int_t num) const; + + // To build an iterator for this class + operator TIterator*() const; + + ClassDef(MSignalCam, 1) // class for an event containing cerenkov photons +}; + +class MSignalCamIter : public TObjArrayIter +{ +private: + Bool_t fUsedOnly; + Int_t fIdx; + +public: + MSignalCamIter(const MSignalCam *evt, Bool_t usedonly=kTRUE, Bool_t dir=kIterForward) : TObjArrayIter(evt->fPixels, dir), fUsedOnly(usedonly), fIdx(-1) { } + TObject *Next(); + void Reset() { fIdx=-1; TObjArrayIter::Reset(); } + TIterator &operator=(const TIterator &) { return *this; } + Int_t GetIdx() const { return fIdx; } + ClassDef(MSignalCamIter, 0) +}; + +inline MSignalCam::operator TIterator*() const { return new MSignalCamIter(this); } + +#endif Index: trunk/MagicSoft/Mars/msignal/MSignalPix.cc =================================================================== --- trunk/MagicSoft/Mars/msignal/MSignalPix.cc (revision 6856) +++ trunk/MagicSoft/Mars/msignal/MSignalPix.cc (revision 6856) @@ -0,0 +1,102 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Thomas Bretz 12/2000 +! Author(s): Harald Kornmayer 1/2001 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + + +///////////////////////////////////////////////////////////////////////////// +// +// MSignalPix +// +// Storage container for the signal in a pixel in number of photons. +// +// NOTE: This container is NOT ment for I/O. Write it to a file on your +// own risk! +// +// fIsSaturated: boolean variable set to kTRUE whenever one or more of +// the low gain FADC slices of the pixel is in saturation. +// +// Version 2: +// ---------- +// - added fIsSaturated +// +// Version 4: +// ---------- +// - added fIsHGSaturated +// +// Version 5: +// ---------- +// - added fIdxIsland +// +// Version 6: +// ---------- +// - put the '!' into the comment line for +// Bool_t fIsCore; //! the pixel is a Core pixel -> kTRUE +// Short_t fRing; //! NT: number of analyzed rings around the core pixels, fRing>0 means: used, fRing= 0 means: unused, fRing= -1 means: unmapped (no possible to use in the calculation of the image parameters) +// Short_t fIdxIsland; //! the pixel is a Core pixel -> kTRUE +// Bool_t fIsHGSaturated; //! the pixel's high gain is saturated +// This is a queick hack to gain storage space - the structure of +// the container for calibrated data will change soon. +// +//////////////////////////////////////////////////////////////////////////// +#include "MSignalPix.h" + +#include "MLog.h" + +ClassImp(MSignalPix); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor. The pixel is assumed as used and not a core pixel. +// NT 29/04/2003: A pixel is considered used when fRing > 0. +// +MSignalPix::MSignalPix(Float_t phot, Float_t errphot) : + fIsCore(kFALSE), fRing(1), fIdxIsland(-1), + fPhot(phot), fErrPhot(errphot) +{ + MMath::ReducePrecision(fPhot); + MMath::ReducePrecision(fErrPhot); +} + +void MSignalPix::Clear(Option_t *o) +{ + fIsCore = kFALSE; + fRing = 1; + fIdxIsland = -1; + fPhot = 0; + fErrPhot = 0; +} + +// -------------------------------------------------------------------------- +// +// Print information to gLog. +// +void MSignalPix::Print(Option_t *) const +{ + gLog << GetDescriptor();// << " Pixel: "<< fPixId; + gLog << (fRing>0?" Used ":" Unused "); + gLog << (fIsCore?" Core ":" "); + gLog << "Nphot= " << fPhot << " Error(Nphot)=" << fErrPhot << endl; +} Index: trunk/MagicSoft/Mars/msignal/MSignalPix.h =================================================================== --- trunk/MagicSoft/Mars/msignal/MSignalPix.h (revision 6856) +++ trunk/MagicSoft/Mars/msignal/MSignalPix.h (revision 6856) @@ -0,0 +1,61 @@ +#ifndef MARS_MSignalPix +#define MARS_MSignalPix + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif +#ifndef MARS_MMath +#include "MMath.h" +#endif + +class MSignalPix : public MParContainer +{ +private: + Bool_t fIsCore; //! the pixel is a Core pixel -> kTRUE + Short_t fRing; //! NT: number of analyzed rings around the core pixels, fRing>0 means: used, fRing= 0 means: unused, fRing= -1 means: unmapped (no possible to use in the calculation of the image parameters) + Short_t fIdxIsland; //! the pixel is a Core pixel -> kTRUE + + Float_t fPhot; // The number of Cerenkov photons + Float_t fErrPhot; // the error of fPhot + Float_t fArrivalTime; // Calibrated Arrival Time + // FIXME: arrival time t, and it's error sigma t + +public: + MSignalPix(/*Int_t pix=-1,*/ Float_t phot=0, Float_t errphot=0); + + void Clear(Option_t *o=0); + + Float_t GetNumPhotons() const { return fPhot; } + Float_t GetErrorPhot() const { return fErrPhot; } + Float_t GetArrivalTime() const { return fArrivalTime; } + + Bool_t IsPixelUsed() const { return fRing>0; } + Bool_t IsPixelUnmapped() const { return fRing==-1; } + void SetPixelUnused() { fRing=0; } + void SetPixelUsed() { fRing=1; } + void SetPixelUnmapped() { fRing=-1;} + void SetIdxIsland(Short_t num) { fIdxIsland=num; } + Short_t GetIdxIsland() const { return fIdxIsland; } + + void SetRing(UShort_t r) { fRing = r; } + Short_t GetRing() const { return fRing;} + + void SetPixelCore(Bool_t b=kTRUE){ fIsCore = b; } + Bool_t IsPixelCore() const { return fIsCore; } + + void SetNumPhotons(Float_t f) { MMath::ReducePrecision(f); fPhot = f; } + void SetErrorPhot(Float_t f) { MMath::ReducePrecision(f); fErrPhot = f; } + void Set(Float_t np, Float_t ep) { MMath::ReducePrecision(np); MMath::ReducePrecision(ep); fPhot = np; fErrPhot = ep; } + void SetArrivalTime(Float_t tm) { MMath::ReducePrecision(tm); fArrivalTime = tm; } + + void AddNumPhotons(Float_t f) { fPhot += f; } + + void Scale(Float_t f) { fPhot/=f; } + + void Print(Option_t *opt = NULL) const; + Bool_t IsSortable() const { return kTRUE; } + + ClassDef(MSignalPix, 6) // class containing information about the Cerenkov Photons in a pixel +}; + +#endif Index: trunk/MagicSoft/Mars/msignal/Makefile =================================================================== --- trunk/MagicSoft/Mars/msignal/Makefile (revision 6855) +++ trunk/MagicSoft/Mars/msignal/Makefile (revision 6856) @@ -53,4 +53,7 @@ MExtractTimeAndChargeDigitalFilter.cc \ MExtractTimeAndChargeDigitalFilterPeakSearch.cc \ + MSignalCam.cc \ + MSignalPix.cc \ + MSignalCalc.cc \ MArrivalTimeCam.cc \ MArrivalTimePix.cc Index: trunk/MagicSoft/Mars/msignal/SignalLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/msignal/SignalLinkDef.h (revision 6855) +++ trunk/MagicSoft/Mars/msignal/SignalLinkDef.h (revision 6856) @@ -35,7 +35,11 @@ #pragma link C++ class MExtractTimeAndChargeDigitalFilterPeakSearch+; - #pragma link C++ class MArrivalTimeCam+; #pragma link C++ class MArrivalTimePix+; +#pragma link C++ class MSignalPix+; +#pragma link C++ class MSignalCam+; +#pragma link C++ class MSignalCamIter+; +#pragma link C++ class MSignalCalc+; + #endif Index: trunk/MagicSoft/Mars/msql/SqlLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/msql/SqlLinkDef.h (revision 6855) +++ trunk/MagicSoft/Mars/msql/SqlLinkDef.h (revision 6856) @@ -5,6 +5,4 @@ #pragma link off all functions; -#pragma link C++ class MSqlInsertRun+; - #pragma link C++ class MSQLColumn+; #pragma link C++ class MSQLServer+;