/* ======================================================================== *\ ! ! * ! * 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 ! ! \* ======================================================================== */ ////////////////////////////////////////////////////////////////////////////// // // MCerPhotAnal2 // // 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: // MCerPhotEvt // ////////////////////////////////////////////////////////////////////////////// #include "MCerPhotAnal2.h" #include "MParList.h" #include "MLog.h" #include "MLogManip.h" #include "MRawRunHeader.h" #include "MRawEvtData.h" // MRawEvtData::GetNumPixels #include "MCerPhotEvt.h" #include "MPedestalPix.h" #include "MPedestalCam.h" #include "MRawEvtPixelIter.h" ClassImp(MCerPhotAnal2); 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. // MCerPhotAnal2::MCerPhotAnal2(Byte_t b, Byte_t a, const char *name, const char *title) : fBefore(b), fAfter(a) { fName = name ? name : "MCerPhotAnal2"; 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 // - MPedestalCam // // The following output containers are also searched and created if // they were not found: // - MCerPhotEvt // Int_t MCerPhotAnal2::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 = (MCerPhotEvt*)pList->FindCreateObj("MCerPhotEvt"); if (!fCerPhotEvt) return kFALSE; const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); if (!runheader) *fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl; else if (runheader->GetRunType() == kRTMonteCarlo) { fPedestals=NULL; 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 MCerPhotEvt container. // Int_t MCerPhotAnal2::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 nb = 0; Int_t nsb = 0; if (*max==255) sat++; while (ptrlast) { sumb += *ptr; sqb += *ptr* *ptr; nb++; } else { sumsb += *ptr; sqsb += *ptr* *ptr; nsb++; } ptr++; } if (sat>1) { // 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 MCerPhotAnal2::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; }