Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 2885) +++ /trunk/MagicSoft/Mars/Changelog (revision 2886) @@ -4,4 +4,30 @@ -*-*- END OF LINE -*-*- + + 2004/01/22: Abelardo Moralejo + + * manalysis/MMcCalibrationUpdate.[h,cc] + - Now ratio of high to low gain is taken from MCalibrationCam if + it existed previously in the parameter list, instead of being + read again from the MMcFadcHeader. Removed Setter function for + fADC2PhInner, no longer necessary. Fixed error regarding the + pedestal conversion to photons (did not read conversion factor + from preexisting MCalibrationCam object). + + * mcalib/MMcCalibrationCalc.cc + - Changed parameters of the histogram, and also the quantity being + histogrammed. Check that input data come from a noiseless camera + file before proceeding to do the calibration. Introduced lower + size in cut for calibration. Now rhe calibration constant is not + calculated from the mean of photons/ADC counts, but from the peak + of the histogram. + + * macros/starmc.C + - Introduced new scheme. Now there are two loops over two different + sets of files. First loop calculates the calibration constants, + second one does the analysis. Introduced comments. Now the + histogram used in the light calibration is written to the output + file. + 2004/01/22: Thomas Bretz Index: /trunk/MagicSoft/Mars/macros/starmc.C =================================================================== --- /trunk/MagicSoft/Mars/macros/starmc.C (revision 2885) +++ /trunk/MagicSoft/Mars/macros/starmc.C (revision 2886) @@ -1,26 +1,26 @@ /* ======================================================================== *\ -! -! * -! * 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): Abelardo Moralejo 1/2004 -! Thomas Bretz 5/2002 -! -! Copyright: MAGIC Software Development, 2000-2004 -! -! -\* ======================================================================== */ + ! + ! * + ! * 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): Abelardo Moralejo 1/2004 + ! Thomas Bretz 5/2002 + ! + ! Copyright: MAGIC Software Development, 2000-2004 + ! + ! + \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// @@ -37,115 +37,152 @@ void starmc() { - // - // This is a demonstration program which calculates the image - // parameters from Magic Monte Carlo files (output of camera). - - // - // Create a empty Parameter List and an empty Task List - // The tasklist is identified in the eventloop by its name - // - MParList plist; - - MTaskList tlist; - - plist.AddToList(&tlist); - - MSrcPosCam src; - src.SetReadyToSave(); - - plist.AddToList(&src); - - // - // Now setup the tasks and tasklist: - // --------------------------------- - // - MReadMarsFile read("Events"); - - // ------------- user change ----------------- - - read.AddFile("Gamma_zbin*.root"); - - read.DisableAutoScheme(); - - MGeomApply geom; // Reads in geometry from MC file and sets the right sizes for - // several parameter containers. - - MMcPedestalCopy pcopy; - // Copies pedestal data from the MC file run fadc header to the MPedestalCam container. - - MExtractSignal sigextract; - // Define ADC slices to be integrated in high and low gain: - sigextract.SetRange(0, 5, 0, 5); - - MMcCalibrationUpdate mccalibupdate; - - // - // Now introduce conversion factor from ADC counts to photons before camera for - // inner pixels. The corresponding value for outer pixels is then calculated - // automatically. Bear in mind that the conversion factor depend both on the - // parameters used in the camera simulation as well as on the number of - // integrated FADC slices. In the future it should be calculated in a previous - // event loop, either from the same MC file or from a MC calibration file - // written on purpose. - // (FIXME: the conversion must be calculated automatically from the analyzed file) - // - mccalibupdate.SetADC2PhInner(1.2586); + // + // This is a demonstration program which calculates the image + // parameters from Magic Monte Carlo files (output of camera). - MCalibrate calib; // Transforms signals from ADC counts into photons. + // ------------- user change ----------------- - // MBlindPixelCalc blind; - // blind.SetUseInterpolation(); + TString* CalibrationFilename; - MImgCleanStd clean(4.1,3.4); // Applies tail cuts to image. + // + // Comment out next line to disable calibration. In that case the units of the + // MHillas.fSize parameter will be ADC counts (rather, equivalent ADC counts in + // inner pixels, since we correct for the possible differences in gain of outer + // pixels) + // + CalibrationFilename = new TString("../../../nonoise/gammas/Gamma_zbin0_0*.root"); + // File to be used in the calibration (must be a camera file without added noise) + + Char_t* AnalysisFilename = "Proton_zbin*.root"; // File to be analyzed + Char_t* OutFilename = "star_mc.root"; // Output file name + + Float_t CleanLev[2] = {4., 3.}; // Tail cuts for image analysis + + Int_t BinsHigh[2] = {0, 5}; // First and last FADC bin of the range to be integrated, + Int_t BinsLow[2] = {0, 5}; // for high and low gain respectively. + + // ------------------------------------------- + + // + // Create a empty Parameter List and an empty Task List + // The tasklist is identified in the eventloop by its name + // + MParList plist; + + MTaskList tlist; + + plist.AddToList(&tlist); + + MSrcPosCam src; + src.SetReadyToSave(); + + plist.AddToList(&src); + + // + // Now setup the tasks and tasklist: + // --------------------------------- + // + MReadMarsFile read("Events"); + + if (CalibrationFilename) + read.AddFile(CalibrationFilename->Data()); + + read.DisableAutoScheme(); + + MGeomApply geom; // Reads in geometry from MC file and sets the right sizes for + // several parameter containers. + + MMcPedestalCopy pcopy; + // Copies pedestal data from the MC file run fadc header to the MPedestalCam container. + + MExtractSignal sigextract; + // Define ADC slices to be integrated in high and low gain: + sigextract.SetRange(BinsHigh[0], BinsHigh[1], BinsLow[0], BinsLow[1]); + + MMcCalibrationUpdate mccalibupdate; + + MCalibrate calib; // Transforms signals from ADC counts into photons. + + // MBlindPixelCalc blind; + // blind.SetUseInterpolation(); + + MImgCleanStd clean(CleanLev[0], CleanLev[1]); // Applies tail cuts to image. + + MHillasCalc hcalc; // Calculates Hillas parameters not dependent on source position. + MHillasSrcCalc scalc; // Calculates source-dependent Hillas parameters + + MMcCalibrationCalc mccalibcalc; + + tlist.AddToList(&read); + tlist.AddToList(&geom); + tlist.AddToList(&pcopy); + + tlist.AddToList(&sigextract); + tlist.AddToList(&mccalibupdate); + tlist.AddToList(&calib); + tlist.AddToList(&clean); + // tlist.AddToList(&blind); + tlist.AddToList(&hcalc); + + tlist.AddToList(&mccalibcalc); + + // + // Open output file: + // + MWriteRootFile write(OutFilename); // Writes output + write.AddContainer("MRawRunHeader", "RunHeaders"); + write.AddContainer("MMcRunHeader", "RunHeaders"); + write.AddContainer("MSrcPosCam", "RunHeaders"); + write.AddContainer("MMcEvt", "Events"); + write.AddContainer("MHillas", "Events"); + write.AddContainer("MHillasExt", "Events"); + write.AddContainer("MHillasSrc", "Events"); + write.AddContainer("MNewImagePar", "Events"); + + // + // First loop: Calibration loop + // + + MProgressBar bar; + bar.SetWindowName("Calibrating..."); + + MEvtLoop evtloop; + evtloop.SetProgressBar(&bar); + evtloop.SetParList(&plist); + + if (CalibrationFilename) + { + if (!evtloop.Eventloop()) + return; + mccalibcalc->GetHist()->Write(); + } + + // + // Second loop: analysis loop + // + + // + // Change the read task by another one which reads the file we want to analyze: + // + + MReadMarsFile read2("Events"); + read2.AddFile(AnalysisFilename); + read2.DisableAutoScheme(); + tlist.AddToListBefore(&read2, &read, "All"); + tlist.RemoveFromList(&read); + + bar.SetWindowName("Analyzing..."); + + tlist.RemoveFromList(&mccalibcalc); // Removes calibration task from list. + + tlist.AddToList(&scalc); // Calculates Source-dependent Hillas parameters + + tlist.AddToList(&write); // Add task to write output. + + if (!evtloop.Eventloop()) + return; - MHillasCalc hcalc; // Calculates Hillas parameters not dependent on source position. - MHillasSrcCalc scalc; // Calculates source-dependent Hillas parameters - // (!!Preliminary!! Will be removed later!) - - - // ------------- user change ----------------- - - MWriteRootFile write("star_mc.root"); // Writes output - - write.AddContainer("MRawRunHeader", "RunHeaders"); - write.AddContainer("MMcRunHeader", "RunHeaders"); - write.AddContainer("MSrcPosCam", "RunHeaders"); - write.AddContainer("MMcEvt", "Events"); - write.AddContainer("MHillas", "Events"); - write.AddContainer("MHillasExt", "Events"); - write.AddContainer("MHillasSrc", "Events"); - write.AddContainer("MNewImagePar", "Events"); - - - tlist.AddToList(&read); - tlist.AddToList(&geom); - tlist.AddToList(&pcopy); - - tlist.AddToList(&sigextract); - tlist.AddToList(&mccalibupdate); - tlist.AddToList(&calib); - tlist.AddToList(&clean); - // tlist.AddToList(&blind); - tlist.AddToList(&hcalc); - tlist.AddToList(&scalc); - tlist.AddToList(&write); - - // - // Create and set up the eventloop - // - MProgressBar bar; - - MEvtLoop evtloop; - evtloop.SetProgressBar(&bar); - evtloop.SetParList(&plist); - - // - // Execute your analysis - // - if (!evtloop.Eventloop()) - return; - - tlist.PrintStatistics(); + tlist.PrintStatistics(); } Index: /trunk/MagicSoft/Mars/manalysis/MMcCalibrationUpdate.cc =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MMcCalibrationUpdate.cc (revision 2885) +++ /trunk/MagicSoft/Mars/manalysis/MMcCalibrationUpdate.cc (revision 2886) @@ -37,8 +37,9 @@ // MMcFadcHeader // MRawRunHeader +// [MCalibrationCam] (if it existed previously) // // Output Containers: -// MCalibrationCam // MPedPhotCam +// [MCalibrationCam] (if it did not exist previously) // ///////////////////////////////////////////////////////////////////////////// @@ -95,16 +96,4 @@ } - -//--------------------------------------------------------------------------- -// -// Set ADC to photon conversion factors. -// -void MMcCalibrationUpdate::SetADC2PhInner(Float_t x) -{ - fADC2PhInner = x; - - return; -} - // -------------------------------------------------------------------------- // @@ -123,5 +112,5 @@ if ( !fCalCam ) { - *fLog << warn << dbginf << AddSerialNumber("MCalibrationCam") << " does not exist... Creating." << endl; + *fLog << inf << dbginf << AddSerialNumber("MCalibrationCam") << " does not exist... Creating." << endl; fCalCam = (MCalibrationCam*) pList->FindCreateObj(AddSerialNumber("MCalibrationCam")); @@ -135,5 +124,5 @@ { fFillCalibrationCam = kFALSE; - *fLog << warn << dbginf << AddSerialNumber("MCalibrationCam") << " already exists... " << endl; + *fLog << inf << AddSerialNumber("MCalibrationCam") << " already exists... " << endl; } @@ -189,15 +178,4 @@ // - // Set the ADC to photons conversion factor for outer pixels: - // - fADC2PhOuter = fADC2PhInner * - (fHeaderFadc->GetAmplitud() / fHeaderFadc->GetAmplitudOuter()); - - // - // Set the conversion factor between high and low gain: - // - fConversionHiLo = fHeaderFadc->GetLow2HighGain(); - - // // If MCalibrationCam already existed in the parameter list before // MMcCalibrationUpdate::PreProcess was executed (from a @@ -208,4 +186,14 @@ return kTRUE; + // + // Set the ADC to photons conversion factor for outer pixels: + // + fADC2PhOuter = fADC2PhInner * + (fHeaderFadc->GetAmplitud() / fHeaderFadc->GetAmplitudOuter()); + + // + // Set the conversion factor between high and low gain: + // + fConversionHiLo = fHeaderFadc->GetLow2HighGain(); const int num = fCalCam->GetSize(); @@ -264,5 +252,4 @@ // counts for the RMS per slice: // - const Float_t pedestrms = sigpix.IsLoGainUsed()? @@ -279,10 +266,11 @@ MPedPhotPix &pedpix = (*fPedPhotCam)[i]; - Float_t adc2phot = (fGeom->GetPixRatio(i) < fGeom->GetPixRatio(0))? - fADC2PhOuter : fADC2PhInner; + MCalibrationPix &calpix = (*fCalCam)[i]; + Float_t adc2phot = calpix.GetMeanConversionBlindPixelMethod(); + Float_t hi2lo = calpix.GetConversionHiLo(); if (sigpix.IsLoGainUsed()) - pedpix.Set(adc2phot*fConversionHiLo*pedestmean, - adc2phot*fConversionHiLo*pedestrms); + pedpix.Set(adc2phot*hi2lo*pedestmean, + adc2phot*hi2lo*pedestrms); else pedpix.Set(adc2phot*pedestmean, adc2phot*pedestrms); Index: /trunk/MagicSoft/Mars/manalysis/MMcCalibrationUpdate.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MMcCalibrationUpdate.h (revision 2885) +++ /trunk/MagicSoft/Mars/manalysis/MMcCalibrationUpdate.h (revision 2886) @@ -35,6 +35,4 @@ MMcCalibrationUpdate(const char *name=NULL, const char *title=NULL); - void SetADC2PhInner(Float_t x); - ClassDef(MMcCalibrationUpdate, 0) // Task which obtains, for MC files, the pedestal mean and rms, and the calibration factor from ADC counts to photons. }; Index: /trunk/MagicSoft/Mars/mcalib/MMcCalibrationCalc.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MMcCalibrationCalc.cc (revision 2885) +++ /trunk/MagicSoft/Mars/mcalib/MMcCalibrationCalc.cc (revision 2886) @@ -71,5 +71,6 @@ fEvents = 0; - fHistRatio = new TH1F("HistRatio", "Ratio fPassPhotCone/fSize", 1000., 0., 10.); + fHistRatio = new TH1F("HistRatio", "log10(fPassPhotCone/fSize)", 1500., -3., 3.); + fHistRatio->GetXaxis()->SetTitle("log_{10}(fPassPhotCone / fSize) (in photons/ADC count)"); } @@ -165,7 +166,17 @@ } + 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; + return kFALSE; + } + } + // // FIXME: Check that the relevant parameters in the FADC header are the - // same for all the analyzed data! + // same for all the analyzed data, both Calibration and Analyzed data! // @@ -181,11 +192,22 @@ { + // + // 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(); fEvents ++; - fHistRatio->Fill(fMcEvt->GetPassPhotCone()/fHillas->GetSize()); + fHistRatio->Fill(log10(fMcEvt->GetPassPhotCone()/fHillas->GetSize())); return kTRUE; @@ -207,4 +229,25 @@ } + // + // 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(); Index: /trunk/MagicSoft/Mars/mcalib/MMcCalibrationCalc.h =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MMcCalibrationCalc.h (revision 2885) +++ /trunk/MagicSoft/Mars/mcalib/MMcCalibrationCalc.h (revision 2886) @@ -28,5 +28,5 @@ Long_t fEvents; - TH1F* fHistRatio; + TH1F* fHistRatio; // Histogram for monitoring the calibration. Bool_t CheckRunType(MParList *pList) const; @@ -39,4 +39,6 @@ MMcCalibrationCalc(const char *name=NULL, const char *title=NULL); + TH1F* GetHist() { return fHistRatio; } + ClassDef(MMcCalibrationCalc, 0) // Task which obtains, for MC files, the calibration factor from ADC counts to photons. };