Index: trunk/MagicSoft/Mars/mcalib/MCalibrate.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrate.cc (revision 3182) +++ trunk/MagicSoft/Mars/mcalib/MCalibrate.cc (revision 3183) @@ -17,5 +17,5 @@ ! ! Author(s): Javier Lopez 12/2003 -! Modified by: Javier Rico 01/2004 +! Author(s): Javier Rico 01/2004 ! ! Copyright: MAGIC Software Development, 2000-2004 Index: trunk/MagicSoft/Mars/mcalib/MCalibrateData.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrateData.cc (revision 3182) +++ trunk/MagicSoft/Mars/mcalib/MCalibrateData.cc (revision 3183) @@ -17,6 +17,6 @@ ! ! Author(s): Javier Lopez 12/2003 -! Modified by: Javier Rico 01/2004 -! Wolfgang Wittek 02/2004 +! Author(s): Javier Rico 01/2004 +! Author(s): Wolfgang Wittek 02/2004 ! ! Copyright: MAGIC Software Development, 2000-2004 Index: trunk/MagicSoft/Mars/mcalib/MMcCalibrationCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MMcCalibrationCalc.cc (revision 3182) +++ trunk/MagicSoft/Mars/mcalib/MMcCalibrationCalc.cc (revision 3183) @@ -18,5 +18,5 @@ ! Author(s): Abelardo Moralejo, 12/2003 ! -! Copyright: MAGIC Software Development, 2000-2003 +! Copyright: MAGIC Software Development, 2000-2004 ! ! @@ -65,12 +65,9 @@ MMcCalibrationCalc::MMcCalibrationCalc(const char *name, const char *title) { - fName = name ? name : "MMcCalibrationCalc"; - fTitle = title ? title : "Calculate and write conversion factors into MCalibrationCam Container"; - - fADC2Phot = 0.; - fEvents = 0; - - fHistRatio = new TH1F(AddSerialNumber("HistRatio"), "log10(fPassPhotCone/fSize)", 1500, -3., 3.); - fHistRatio->GetXaxis()->SetTitle("log_{10}(fPassPhotCone / fSize) (in photons/ADC count)"); + fName = name ? name : "MMcCalibrationCalc"; + fTitle = title ? title : "Calculate and write conversion factors into MCalibrationCam Container"; + + fHistRatio = new TH1F(AddSerialNumber("HistRatio"), "log10(fPassPhotCone/fSize)", 1500, -3., 3.); + fHistRatio->SetXTitle("log_{10}(fPassPhotCone / fSize) [phot/ADC count]"); } @@ -83,12 +80,12 @@ Bool_t MMcCalibrationCalc::CheckRunType(MParList *pList) const { - const MRawRunHeader *run = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); - if (!run) - { - *fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl; - return kTRUE; - } - - return run->GetRunType() == kRTMonteCarlo; + const MRawRunHeader *run = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); + if (!run) + { + *fLog << warn << "Warning - cannot check file type, MRawRunHeader not found." << endl; + return kTRUE; + } + + return run->GetRunType() == kRTMonteCarlo; } @@ -99,36 +96,36 @@ Int_t MMcCalibrationCalc::PreProcess(MParList *pList) { - - fCalCam = (MCalibrationCam*) pList->FindObject(AddSerialNumber("MCalibrationCam")); - - if ( !fCalCam ) - { - *fLog << err << dbginf << "Cannot find " << AddSerialNumber("MCalibrationCam") << "... aborting." << endl; - return kFALSE; - } - - fHillas = (MHillas*) pList->FindCreateObj(AddSerialNumber("MHillas")); - if ( !fHillas) - { - *fLog << err << dbginf << "Cannot find " << AddSerialNumber("MHillas") << "... aborting." << endl; - return kFALSE; - } - - fNew = (MNewImagePar*) pList->FindCreateObj(AddSerialNumber("MNewImagePar")); - if ( !fNew) - { - *fLog << err << dbginf << "Cannot find " << AddSerialNumber("MNewImagePar") << "... aborting." << endl; - return kFALSE; - } - - fMcEvt = (MMcEvt*) pList->FindCreateObj(AddSerialNumber("MMcEvt")); - if ( !fMcEvt) - { - *fLog << err << dbginf << "Cannot find " << AddSerialNumber("MMcEvt") << "... aborting." << endl; - return kFALSE; - } - - return kTRUE; - + fHistRatio->Reset(); + fADC2Phot = 0; + + fCalCam = (MCalibrationCam*) pList->FindObject(AddSerialNumber("MCalibrationCam")); + if (!fCalCam) + { + *fLog << err << AddSerialNumber("MCalibrationCam") << "not found... aborting." << endl; + return kFALSE; + } + + fHillas = (MHillas*) pList->FindObject(AddSerialNumber("MHillas")); + if ( !fHillas) + { + *fLog << err << AddSerialNumber("MHillas") << "not found... aborting." << endl; + return kFALSE; + } + + fNew = (MNewImagePar*)pList->FindObject(AddSerialNumber("MNewImagePar")); + if (!fNew) + { + *fLog << err << AddSerialNumber("MNewImagePar") << "not found... aborting." << endl; + return kFALSE; + } + + fMcEvt = (MMcEvt*) pList->FindObject(AddSerialNumber("MMcEvt")); + if (!fMcEvt) + { + *fLog << err << AddSerialNumber("MMcEvt") << "not found... aborting." << endl; + return kFALSE; + } + + return kTRUE; } @@ -144,36 +141,36 @@ // if (!CheckRunType(pList)) - { - *fLog << err << dbginf << "This is no MC file... aborting." << endl; + { + *fLog << err << "MMcCalibrationCalc can only used with MC files... aborting." << endl; return kFALSE; - } + } // // Now check the existence of all necessary containers. // - fGeom = (MGeomCam*) pList->FindObject(AddSerialNumber("MGeomCam")); - if ( ! fGeom ) - { - *fLog << err << dbginf << "Cannot find " << AddSerialNumber("MGeomCam") << "... aborting." << endl; + if (!fGeom) + { + *fLog << err << AddSerialNumber("MGeomCam") << " mot found... aborting." << endl; return kFALSE; - } + } fHeaderFadc = (MMcFadcHeader*)pList->FindObject(AddSerialNumber("MMcFadcHeader")); if (!fHeaderFadc) - { - *fLog << err << dbginf << AddSerialNumber("MMcFadcHeader") << " not found... aborting." << endl; + { + *fLog << err << AddSerialNumber("MMcFadcHeader") << " not found... aborting." << endl; return kFALSE; - } + } for (UInt_t ipix = 0; ipix < fGeom->GetNumPixels(); ipix++) - { + { if (fHeaderFadc->GetPedestalRmsHigh(ipix) > 0 || - fHeaderFadc->GetPedestalRmsLow(ipix) > 0 ) - { - *fLog << err << endl << endl << dbginf << "You are trying to calibrate the data using a Camera file produced with added noise. Please use a noiseless file for calibration. Aborting..." << endl << endl; + fHeaderFadc->GetPedestalRmsLow(ipix) > 0 ) + { + *fLog << err << "Trying to calibrate the data using a Camera file produced with added noise." << endl; + *fLog << "Please use a noiseless file for calibration... aborting." << endl << endl; return kFALSE; - } - } + } + } return kTRUE; @@ -187,76 +184,73 @@ Int_t MMcCalibrationCalc::Process() { - - // - // Exclude events with some saturated pixel - // - if ( fNew->GetNumSaturatedPixels() > 0 ) + // + // Exclude events with some saturated pixel + // + if (fNew->GetNumSaturatedPixels()>0) + return kTRUE; + + // + // Exclude events with low Size (larger fluctuations) + // FIXME? The present cut (1000 "inner-pixel-counts") is somehow + // arbitrary. Might it be optimized? + // + if (fHillas->GetSize()<1000) + return kTRUE; + + fADC2Phot += fMcEvt->GetPassPhotCone()/fHillas->GetSize(); + + fHistRatio->Fill(TMath::Log10(fMcEvt->GetPassPhotCone()/fHillas->GetSize())); + return kTRUE; - - // - // Exclude events with low Size (larger fluctuations) - // FIXME? The present cut (1000 "inner-pixel-counts") is somehow arbitrary. - // Might it be optimized? - // - if ( fHillas->GetSize() < 1000 ) +} + +// -------------------------------------------------------------------------- +// +// Fill the MCalibrationCam object +// +Int_t MMcCalibrationCalc::PostProcess() +{ + const Stat_t n = fHistRatio->GetEntries(); + if (n<1) + { + *fLog << err << "No events read... aborting." << endl; + return kFALSE; + } + + fADC2Phot /= n; + + // + // For the calibration we no longer use the mean, + // but the peak of the distribution: + // + const Int_t reach = 2; + + Stat_t summax = 0; + Int_t mode = 0; + + // FIXME: Is this necessary? We could use GetMaximumBin instead.. + for (Int_t ibin = 1+reach; ibin <= fHistRatio->GetNbinsX()-reach; ibin++) + { + const Stat_t sum = fHistRatio->Integral(ibin-reach, ibin+reach); + + if (sum <= summax) + continue; + + summax = sum; + mode = ibin; + } + + fADC2Phot = TMath::Power(10, fHistRatio->GetBinCenter(mode)); + + const Int_t num = fCalCam->GetSize(); + for (int i=0; iGetPassPhotCone() / fHillas->GetSize(); - fEvents ++; - - fHistRatio->Fill(log10(fMcEvt->GetPassPhotCone()/fHillas->GetSize())); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Fill the MCalibrationCam object -// -Int_t MMcCalibrationCalc::PostProcess() -{ - - if (fEvents > 0) - fADC2Phot /= fEvents; - else - { - *fLog << err << dbginf << "No events were read! Aborting." << endl; - return kFALSE; - } - - // - // For the calibration we no longer use the mean, but thepeak of the distribution: - // - - Float_t summax = 0.; - Int_t mode = 0; - Int_t reach = 2; - for (Int_t ibin = 1+reach; ibin <= fHistRatio->GetNbinsX()-reach; ibin++) - { - Float_t sum = 0; - for(Int_t k = ibin-reach; k <= ibin+reach; k++) - sum += fHistRatio->GetBinContent(k); - if (sum > summax) - { - summax = sum; - mode = ibin; - } - } - - fADC2Phot = pow(10., fHistRatio->GetXaxis()->GetBinCenter(mode)); - - const int num = fCalCam->GetSize(); - - for (int i=0; i