Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 2626) +++ trunk/MagicSoft/Mars/Changelog (revision 2627) @@ -4,4 +4,33 @@ -*-*- END OF LINE -*-*- + 2003/12/09: Markus Gaug + + * mhist/MHCamera: + - added feature MHCamError to set errors in histograms of MStatusDisplay + + * mhist/MHCalibration*: + - all hists doubles to separate High Gain and Low Gain + + * manalysis/MCalibration* + - separate Hi Gain and Lo Gain + - use MExtractSignals to extract charges + + * manalysis/MExtractSignals: + - new class to extract signals + + * manalysis/MExtractedSignalCam: + - new class to hold extracted signals + + * manalysis/MExtractedSignalPix: + - new class to hold extracted signals for pixel + + * manalysis/Makefile + * manalysis/AnalysisLinkDef: + - contain MExtractedSignalCam, MExtractedSignalPix, MExtractedSignal + + * macros/calibration.C + - EventDisplay which allows to get plot by clicking on pixel + + 2003/12/08: Thomas Bretz Index: trunk/MagicSoft/Mars/macros/calibration.C =================================================================== --- trunk/MagicSoft/Mars/macros/calibration.C (revision 2626) +++ trunk/MagicSoft/Mars/macros/calibration.C (revision 2627) @@ -23,6 +23,6 @@ \* ======================================================================== */ -void calibration(TString pedname="../../Mars-0.8.2/20031102_02399_P_Unavailable_E.root", - TString calname="../../Mars-0.8.2/20031102_02400_D_Flip500Hz_E.root") +void calibration(TString pedname="/data/MAGIC/rootdata/2003_12_01/20031130_03340_P_CrabNebula_E.root", + TString calname="/data/MAGIC/rootdata/2003_12_01/20031130_03341_C_CrabNebula_E.root") { @@ -62,8 +62,8 @@ tlist.AddToList(&fill); - // MStatusDisplay *d1 = new MStatusDisplay; + MStatusDisplay *d1 = new MStatusDisplay; // Set update time to 3s - // d1->SetUpdateTime(3000); + d1->SetUpdateTime(3000); // @@ -72,5 +72,5 @@ MEvtLoop evtloop; evtloop.SetParList(&plist); - // evtloop.SetDisplay(d1); + evtloop.SetDisplay(d1); // @@ -114,5 +114,5 @@ // we have to set the color of the pulser LED by hand // - calcalc.SetPulserColor(MCalibrationCalc::kEBlue); + calcalc.SetPulserColor(MCalibrationCalc::kECT1); tlist2.AddToList(&read2); @@ -145,22 +145,4 @@ cam.Print(); - Int_t pixnr; - - while (1) - { - - cout << "Which pixel number do you want to display? (Press 0 to exit)" << endl; - cin >> pixnr; - if (pixnr == 0) - break; - - if (pixnr >= 577) - break; - - MCalibrationPix *pix = cam->GetCalibrationPix(pixnr); - pix->Draw(); - - } - // // Here we are confronted to a serious bug in ROOT: @@ -174,7 +156,7 @@ MHCamera &disp0 = *h.GetHistByName(); MHCamera disp1 (geomcam, "MCalibrationCam;q", "Fitted Mean Charges"); - MHCamera disp2 (geomcam, "MCalibrationCam;errq", "Error of Fitted Mean Charges"); + // MHCamera disp2 (geomcam, "MCalibrationCam;errq", "Error of Fitted Mean Charges"); MHCamera disp3 (geomcam, "MCalibrationCam;sigmaq", "Sigma of Fitted Mean Charges"); - MHCamera disp4 (geomcam, "MCalibrationCam;errsigmaq", "Error of Sigma of Fitted Mean Charges"); + // MHCamera disp4 (geomcam, "MCalibrationCam;errsigmaq", "Error of Sigma of Fitted Mean Charges"); MHCamera disp5 (geomcam, "MCalibrationCam;probq", "Probability of Fit"); MHCamera disp6 (geomcam, "MCalibrationCam;t", "Arrival Times"); @@ -183,33 +165,41 @@ MHCamera disp9 (geomcam, "MCalibrationCam;ped", "Pedestals"); MHCamera disp10 (geomcam, "MCalibrationCam;pedrms", "Pedestal RMS"); - MHCamera disp11 (geomcam, "MCalibrationCam;rq", "Reduced Charges"); - MHCamera disp12 (geomcam, "MCalibrationCam;rsigma", "Reduced Sigmas"); - MHCamera disp13 (geomcam, "MCalibrationCam;phe", "Nr. of Phe's (F-Factor Method)"); - MHCamera disp14 (geomcam, "MCalibrationCam;convphe", "Conversion Factor (F-Factor Method)"); - MHCamera disp15 (geomcam, "MCalibrationCam;photons", "Nr. of Photons (Blind Pixel Method)"); - MHCamera disp16 (geomcam, "MCalibrationCam;convphot", "Conversion Factor (Blind Pixel Method)"); + MHCamera disp11 (geomcam, "MCalibrationCam;rsigma", "Reduced Sigmas"); + MHCamera disp12 (geomcam, "MCalibrationCam;phe", "Nr. of Phe's (F-Factor Method)"); + MHCamera disp13 (geomcam, "MCalibrationCam;convphe", "Conversion Factor (F-Factor Method)"); + MHCamera disp14 (geomcam, "MCalibrationCam;photons", "Nr. of Photons (Blind Pixel Method)"); + MHCamera disp15 (geomcam, "MCalibrationCam;convphot", "Conversion Factor (Blind Pixel Method)"); disp1.SetCamContent(*cam, 0); - disp2.SetCamContent(*cam, 1); + disp1.SetCamError(*cam,1); + // disp2.SetCamContent(*cam, 1); + disp3.SetCamContent(*cam, 2); - disp4.SetCamContent(*cam, 3); + disp3.SetCamError(*cam,3); + // disp4.SetCamContent(*cam, 3); disp5.SetCamContent(*cam, 4); + disp6.SetCamContent(*cam, 5); + disp6.SetCamError(*cam,6); disp7.SetCamContent(*cam, 6); disp8.SetCamContent(*cam, 7); + disp9.SetCamContent(*cam, 8); + disp9.SetCamError(*cam,9); disp10.SetCamContent(*cam, 9); + disp11.SetCamContent(*cam, 10); + disp12.SetCamContent(*cam, 11); disp13.SetCamContent(*cam, 12); disp14.SetCamContent(*cam, 13); disp15.SetCamContent(*cam, 14); - disp16.SetCamContent(*cam, 15); + // disp16.SetCamError(*cam, 16); disp1.SetYTitle("Q [FADC counts]"); - disp2.SetYTitle("\\Delta Q [FADC counts]"); + // disp2.SetYTitle("\\Delta Q [FADC counts]"); disp3.SetYTitle("\\sigma_{Q} [FADC counts]"); - disp4.SetYTitle("\\Delta {\\sigma_{Q}} [FADC counts]"); + // disp4.SetYTitle("\\Delta {\\sigma_{Q}} [FADC counts]"); disp5.SetYTitle("P [au]"); disp6.SetYTitle("T [FADC slices]"); @@ -218,10 +208,9 @@ disp9.SetYTitle("P [FADC counts/ slice ]"); disp10.SetYTitle("RMS_{P} [FADC counts / slice ]"); - disp11.SetYTitle("Q [FADC counts]"); - disp12.SetYTitle("\\sigma^2_{Q} - RMS^2_{P} [FADC counts^2]"); - disp13.SetYTitle("Nr Phe's"); - disp14.SetYTitle("Conversion Factor [Phe/FADC count]"); - disp15.SetYTitle("Nr Photons"); - disp16.SetYTitle("Conversion Factor [Ph/FADC count]"); + disp11.SetYTitle("\\sigma^2_{Q} - RMS^2_{P} [FADC counts^2]"); + disp12.SetYTitle("Nr Phe's"); + disp13.SetYTitle("Conversion Factor [Phe/FADC count]"); + disp14.SetYTitle("Nr Photons"); + disp15.SetYTitle("Conversion Factor [Ph/FADC count]"); MStatusDisplay *d2 = new MStatusDisplay; @@ -230,5 +219,5 @@ d2->SetUpdateTime(1000); - TCanvas *c1 = &d2->AddTab("Charges Mean"); + TCanvas *c1 = &d2->AddTab("Fitted Charges"); c1->Divide(2, 2); @@ -238,4 +227,5 @@ gStyle->SetOptStat(1111); obj=disp1.DrawCopy("hist"); + ((MHCamera*)obj)->AddNotify(*cam); c1->cd(3); @@ -245,5 +235,6 @@ c1->cd(2); gStyle->SetOptStat(1101); - obj=disp2.DrawCopy("hist"); + obj=disp3.DrawCopy("hist"); + ((MHCamera*)obj)->AddNotify(*cam); c1->cd(4); @@ -251,23 +242,4 @@ obj->Draw(); - TCanvas *c11 = &d2->AddTab("Charges Sigma"); - c11->Divide(2, 2); - - c11->cd(1); - gStyle->SetOptStat(1101); - obj=disp3.DrawCopy("hist"); - - c11->cd(3); - gPad->SetBorderMode(0); - obj->Draw(); - - c11->cd(2); - gStyle->SetOptStat(1101); - obj=disp4.DrawCopy("hist"); - - c11->cd(4); - gPad->SetBorderMode(0); - obj->Draw(); - TCanvas *c12 = &d2->AddTab("Fit Prob."); @@ -277,4 +249,5 @@ gStyle->SetOptStat(1101); obj=disp5.DrawCopy("hist"); + ((MHCamera*)obj)->AddNotify(*cam); c12->cd(2); @@ -288,4 +261,5 @@ gStyle->SetOptStat(1111); obj=disp6.DrawCopy("hist"); + ((MHCamera*)obj)->AddNotify(*cam); c2->cd(4); @@ -295,4 +269,5 @@ gStyle->SetOptStat(1101); obj=disp7.DrawCopy("hist"); + ((MHCamera*)obj)->AddNotify(*cam); c2->cd(5); @@ -302,4 +277,5 @@ gStyle->SetOptStat(1101); obj=disp8.DrawCopy("hist"); + ((MHCamera*)obj)->AddNotify(*cam); c2->cd(6); @@ -312,4 +288,5 @@ gStyle->SetOptStat(1111); obj=disp9.DrawCopy("hist"); + ((MHCamera*)obj)->AddNotify(*cam); c3->cd(3); @@ -319,4 +296,5 @@ gStyle->SetOptStat(1111); obj=disp10.DrawCopy("hist"); + ((MHCamera*)obj)->AddNotify(*cam); c3->cd(4); @@ -324,18 +302,12 @@ TCanvas *c4 = &d2->AddTab("Reduced Charges"); - c4->Divide(2, 2); + c4->Divide(2,1); c4->cd(1); gStyle->SetOptStat(1111); obj=disp11.DrawCopy("hist"); + ((MHCamera*)obj)->AddNotify(*cam); c4->cd(3); - obj->Draw(); - - c4->cd(2); - gStyle->SetOptStat(1101); - obj=disp12.DrawCopy("hist"); - - c4->cd(4); obj->Draw(); @@ -345,12 +317,14 @@ c5->cd(1); gStyle->SetOptStat(1111); + obj=disp12.DrawCopy("hist"); + ((MHCamera*)obj)->AddNotify(*cam); + + c5->cd(3); + obj->Draw(); + + c5->cd(2); + gStyle->SetOptStat(1101); obj=disp13.DrawCopy("hist"); - - c5->cd(3); - obj->Draw(); - - c5->cd(2); - gStyle->SetOptStat(1101); - obj=disp14.DrawCopy("hist"); + ((MHCamera*)obj)->AddNotify(*cam); c5->cd(4); @@ -362,19 +336,18 @@ c6->cd(1); gStyle->SetOptStat(1111); + obj=disp14.DrawCopy("hist"); + ((MHCamera*)obj)->AddNotify(*cam); + + c6->cd(3); + obj->Draw(); + + c6->cd(2); + gStyle->SetOptStat(1101); obj=disp15.DrawCopy("hist"); - - c6->cd(3); - obj->Draw(); - - c6->cd(2); - gStyle->SetOptStat(1101); - obj=disp16.DrawCopy("hist"); + ((MHCamera*)obj)->AddNotify(*cam); c6->cd(4); obj->Draw(); - -#endif - } Index: trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 2626) +++ trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 2627) @@ -69,4 +69,8 @@ #pragma link C++ class MCalibrationCalc+; +#pragma link C++ class MExtractedSignalCam+; +#pragma link C++ class MExtractedSignalPix+; +#pragma link C++ class MExtractSignal+; + #endif Index: trunk/MagicSoft/Mars/manalysis/MCalibrationBlindPix.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCalibrationBlindPix.h (revision 2626) +++ trunk/MagicSoft/Mars/manalysis/MCalibrationBlindPix.h (revision 2627) @@ -48,12 +48,12 @@ Float_t GetErrSigma1() const { return fErrSigma1; } - Float_t GetTime() const { return fTime; } - Float_t GetErrTime() const { return fErrTime; } + Float_t GetTime() const { return fTime; } + Float_t GetErrTime() const { return fErrTime; } - Bool_t FillCharge(Int_t q) { return fHist->FillBlindPixelCharge(q); } - Bool_t FillTime(Int_t t) { return fHist->FillBlindPixelTime(t); } + Bool_t FillCharge(Float_t q) { return fHist->FillBlindPixelCharge(q); } + Bool_t FillTime(Int_t t) { return fHist->FillBlindPixelTime(t); } Bool_t FillRChargevsTime(Float_t rq, Int_t t) { return fHist->FillBlindPixelChargevsN(rq,t); } - Bool_t IsValid() { return fLambda > 0. || fErrLambda > 0.; } + Bool_t IsValid() { return fLambda > 0. || fErrLambda > 0.; } Bool_t FitCharge(); Index: trunk/MagicSoft/Mars/manalysis/MCalibrationCalc.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCalibrationCalc.cc (revision 2626) +++ trunk/MagicSoft/Mars/manalysis/MCalibrationCalc.cc (revision 2627) @@ -53,20 +53,20 @@ // Every MCalibrationPix holds a histogram class, // // MHCalibrationPixel which itself hold histograms of type: // -// HCharge(npix) (distribution of summed FADC time slice entries) // -// HTime(npix) (distribution of position of maximum) // -// HChargevsN(npix) (distribution of charges vs. event number. // -// // -// PostProcess: All histograms HCharge(npix) are fitted to a Gaussian // -// All histograms HTime(npix) are fitted to a Gaussian // -// The histogram HBlindPixelCharge (blind pixel) is fitted to a single // -// PhE fit // -// The histogram HBlindPixelTime (blind pixel) is fitted to a Gaussian // -// The histograms of the PIN Diode are fitted to Gaussians // -// // -// Fits can be excluded via the commands: // -// MalibrationCam::SetSkipTimeFits() (skip all time fits) // -// MalibrationCam::SetSkipBlindPixelFits() (skip all blind pixel fits) // -// MalibrationCam::SetSkipPinDiodeFits() (skip all PIN Diode fits) // -// // +// HCharge(npix) (distribution of summed FADC time slice entries) +// HTime(npix) (distribution of position of maximum) +// HChargevsN(npix) (distribution of charges vs. event number. +// +// PostProcess: All histograms HCharge(npix) are fitted to a Gaussian +// All histograms HTime(npix) are fitted to a Gaussian +// The histogram HBlindPixelCharge (blind pixel) is fitted to a single +// PhE fit +// The histogram HBlindPixelTime (blind pixel) is fitted to a Gaussian +// The histograms of the PIN Diode are fitted to Gaussians +// +// Fits can be excluded via the commands: +// MalibrationCam::SetSkipTimeFits() (skip all time fits) +// MalibrationCam::SetSkipBlindPixelFits() (skip all blind pixel fits) +// MalibrationCam::SetSkipPinDiodeFits() (skip all PIN Diode fits) +// ////////////////////////////////////////////////////////////////////////////// @@ -169,13 +169,17 @@ case kEBlue: fCalibrations->SetColor(MCalibrationCam::kECBlue); - break; + break; case kEGreen: fCalibrations->SetColor(MCalibrationCam::kECGreen); - break; + break; case kEUV: fCalibrations->SetColor(MCalibrationCam::kECUV); + break; + case kECT1: + fCalibrations->SetColor(MCalibrationCam::kECCT1); + break; + default: + fCalibrations->SetColor(MCalibrationCam::kECCT1); } - - fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam"); @@ -185,5 +189,4 @@ return kFALSE; } - return kTRUE; @@ -245,4 +248,5 @@ MRawEvtPixelIter pixel(fRawEvt); + // // Create a first loop to sort out the cosmics ... // @@ -266,5 +270,5 @@ Float_t pedrms = ped.GetPedestalRms(); - if ((float)sum < (pedes*fNumHiGainSamples)+(1.5*pedrms) ) + if ((float)sum < ((pedes*fNumHiGainSamples)+(30.*pedrms)) ) cosmicpix++; } @@ -274,64 +278,60 @@ { fCosmics++; - return kTRUE; + return kCONTINUE; } pixel.Reset(); + // + // Create a second loop to do fill the calibration histograms + // + while (pixel.Next()) { UShort_t sat = 0; - UShort_t lowgainoverflow = 0; const Int_t pixid = pixel.GetPixelId(); - Byte_t *ptr = pixel.GetHiGainSamples(); Byte_t mid = pixel.GetIdxMaxHiGainSample(); UInt_t max = pixel.GetMaxHiGainSample(); - Int_t sum = (max > gkSaturationLimit // overflow ? - ? sat++, pixel.GetSumLoGainSamples() // take Low Gain - : pixel.GetSumHiGainSamples()); // no overflow - - if (sat) - { - - ptr = pixel.GetLoGainSamples(); - max = pixel.GetMaxLoGainSample(); - mid = pixel.GetIdxMaxLoGainSample(); - - // - // FIXME: It seems the conversion HiGain LoGain is already - // performed in the data?!? - // - sum = (max > gkSaturationLimit // overflow of LoGain ??? -> GimmeABreak!!! - ? lowgainoverflow++, gkLoGainOverFlow // OUCH (Florian was maybe right) - : sum*gkConversionHiLo ); // OUFF (Florian was wrong) !! - // : sum ); - - if (lowgainoverflow) - { - *fLog << err << dbginf << "Warning: Saturation of LoGain reached in pixel: " << pixid << " " - << " sum = " << sum << endl; - fHistOverFlow++; - } - } - - MPedestalPix &ped = (*fPedestals)[pixid]; - MCalibrationPix &pix = (*fCalibrations)[pixid]; - + MPedestalPix &ped = (*fPedestals)[pixid]; + MCalibrationPix &pix = (*fCalibrations)[pixid]; + Float_t pedes = ped.GetPedestal(); - Float_t pedrms = ped.GetPedestalRms(); - + Float_t sum; + // // FIXME: This is preliminary, we will change to pedestals per slice!!! // Assume pedestals per time slice ==> multiply with number of slices // - pedes *= (sat ? fNumLoGainSamples : fNumHiGainSamples ); - pedrms *= (sat ? fNumLoGainSamples : fNumHiGainSamples ); - - Float_t rsum = (float)sum - pedes; - + + if (max < gkSaturationLimit) // take Hi Gain + { + sum = (float)pixel.GetSumHiGainSamples() - pedes*fNumHiGainSamples; + } + else // Lo Gain + { + + sat++; + pix.SetHiGainSaturation(); + + sum = (float)pixel.GetSumLoGainSamples() - pedes*fNumLoGainSamples ; + sum *= gkConversionHiLo; + + max = pixel.GetMaxLoGainSample(); + mid = pixel.GetIdxMaxLoGainSample(); + + if (max > gkSaturationLimit) + { + *fLog << err << dbginf + << "Warning: Saturation of Lo Gain reached in pixel: " + << pixid << " " << " sum = " << sum << endl; + fHistOverFlow++; + } + + } + switch(pixid) { @@ -339,8 +339,4 @@ case gkCalibrationBlindPixelId: - // - // FIXME: This works only when the blind pixel ID is much larger than - // the rest of the pixels (which is the case right now) - // if (!blindpixel.FillCharge(sum)) *fLog << warn << @@ -351,5 +347,5 @@ "Overflow or Underflow occurred filling Blind Pixel time = " << (int)mid << endl; - if (!blindpixel.FillRChargevsTime(rsum,fEvents)) + if (!blindpixel.FillRChargevsTime(sum,fEvents)) *fLog << warn << "Overflow or Underflow occurred filling Blind Pixel eventnr = " << fEvents << endl; @@ -362,5 +358,5 @@ *fLog << warn << "Overflow or Underflow occurred filling HTime: time = " << (int)mid << endl; - if (!pindiode.FillRChargevsTime(rsum,fEvents)) + if (!pindiode.FillRChargevsTime(sum,fEvents)) *fLog << warn << "Overflow or Underflow occurred filling HChargevsN: eventnr = " << fEvents << endl; @@ -368,20 +364,38 @@ default: - if (!pix.FillCharge(sum)) - *fLog << warn << "Could not fill Charge of pixel: " << pixid - << " signal = " << sum << endl; - - // - // Fill the reduced charge into the control histo for better visibility - // - if (!pix.FillRChargevsTime(rsum,fEvents)) - *fLog << warn << "Could not fill red. Charge vs. EvtNr of pixel: " << pixid - << " signal = " << rsum << " event Nr: " << fEvents << endl; - - - if (!pix.FillTime((int)mid)) - *fLog << warn << "Could not fill Time of pixel: " << pixid << " time = " << (int)mid << endl; - - + if (sat) + { + + if (!pix.FillChargeLoGain(sum)) + *fLog << warn << "Could not fill Lo Gain Charge of pixel: " << pixid + << " signal = " << sum << endl; + + if (!pix.FillTimeLoGain((int)mid)) + *fLog << warn << "Could not fill Lo Gain Time of pixel: " + << pixid << " time = " << (int)mid << endl; + + // + // Fill the reduced charge into the control histo for better visibility + // + if (!pix.FillRChargevsTimeLoGain(sum,fEvents)) + *fLog << warn << "Could not fill Lo Gain Charge vs. EvtNr of pixel: " + << pixid << " signal = " << sum << " event Nr: " << fEvents << endl; + + } + else + { + if (!pix.FillChargeHiGain(sum)) + *fLog << warn << "Could not fill Hi Gain Charge of pixel: " << pixid + << " signal = " << sum << endl; + + if (!pix.FillTimeHiGain((int)mid)) + *fLog << warn << "Could not fill Hi Gain Time of pixel: " + << pixid << " time = " << (int)mid << endl; + + if (!pix.FillRChargevsTimeHiGain(sum,fEvents)) + *fLog << warn << "Could not fill Hi Gain Charge vs. EvtNr of pixel: " + << pixid << " signal = " << sum << " event Nr: " << fEvents << endl; + } + } /* switch(pixid) */ @@ -396,4 +410,5 @@ *fLog << inf << endl; *fLog << GetDescriptor() << " Cut Histogram Edges" << endl; + // // Cut edges to make fits and viewing of the hists easier @@ -414,9 +429,13 @@ { if (blindpixel.FitCharge()) - if (!fCalibrations->CalcNrPhotInnerPixel()) - *fLog << err << dbginf << "Could not calculate Number of photons from the blind pixel " << endl; + { + if (!fCalibrations->CalcNrPhotInnerPixel()) + *fLog << err << dbginf << "Could not calculate Number of photons from the blind pixel " << endl; + } else - *fLog << err << dbginf << "Could not fit the blind pixel " << endl; - + { + *fLog << err << dbginf << "Could not fit the blind pixel " << endl; + } + if (!blindpixel.FitTime()) *fLog << warn << "Could not the Times of the blind pixel " << endl; Index: trunk/MagicSoft/Mars/manalysis/MCalibrationCalc.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCalibrationCalc.h (revision 2626) +++ trunk/MagicSoft/Mars/manalysis/MCalibrationCalc.h (revision 2627) @@ -54,5 +54,5 @@ public: - enum PulserColor_t { kEGreen, kEBlue, kEUV }; + enum PulserColor_t { kEGreen, kEBlue, kEUV, kECT1 }; private: Index: trunk/MagicSoft/Mars/manalysis/MCalibrationCam.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCalibrationCam.cc (revision 2626) +++ trunk/MagicSoft/Mars/manalysis/MCalibrationCam.cc (revision 2627) @@ -57,5 +57,6 @@ MCalibrationCam::MCalibrationCam(const char *name, const char *title) : fMeanNrPhotAvailable(kFALSE), - fMeanNrPhotInnerPix(-1.) + fMeanNrPhotInnerPix(-1.), + fMeanNrPhotInnerPixErr(-1.) { fName = name ? name : "MCalibrationCam"; @@ -153,5 +154,5 @@ Bool_t MCalibrationCam::IsPixelFitted(Int_t idx) const { - return ((*this)[idx].GetRCharge() > 0. && (*this)[idx].GetErrRCharge() > 0.); + return ((*this)[idx].GetCharge() > 0. && (*this)[idx].GetErrCharge() > 0.); } @@ -359,4 +360,7 @@ int id = 0; + *fLog << "Succesfully calibrated pixels:" << endl; + *fLog << endl; + TIter Next(fPixels); MCalibrationPix *pix; @@ -364,10 +368,34 @@ { - *fLog << id << " Pedestals: " << pix->GetPed() << " +- " << pix->GetPedRms() << " Charge: " - << pix->GetCharge() << " Reduced Charge: " << pix->GetRCharge() << " +- " - << pix->GetSigmaCharge() << " Reduced Sigma: " << pix->GetRSigma() << endl; - - id++; - } + if (pix->GetCharge() >= 0.) + { + *fLog << pix->GetPixId() << " Pedestals: " << pix->GetPed() << " +- " << pix->GetPedRms() + << " Reduced Charge: " << pix->GetCharge() << " +- " + << pix->GetSigmaCharge() << " Reduced Sigma: " << pix->GetRSigma() << endl; + id++; + } + } + + *fLog << id << " succesful pixels :-))" << endl; + id = 0; + + *fLog << endl; + *fLog << "Pixels with errors:" << endl; + *fLog << endl; + + TIter Next2(fPixels); + while ((pix=(MCalibrationPix*)Next2())) + { + + if (pix->GetCharge() == -1.) + { + *fLog << pix->GetPixId() << " Pedestals: " << pix->GetPed() << " +- " << pix->GetPedRms() + << " Reduced Charge: " << pix->GetCharge() << " +- " + << pix->GetSigmaCharge() << " Reduced Sigma: " << pix->GetRSigma() << endl; + id++; + } + } + *fLog << id << " pixels with errors :-((" << endl; + } @@ -379,53 +407,60 @@ { case 0: - val = (*this)[idx].GetCharge(); - break; + val = (*this)[idx].GetCharge(); + break; case 1: - val = (*this)[idx].GetErrCharge(); - break; + val = (*this)[idx].GetErrCharge(); + break; case 2: - val = (*this)[idx].GetSigmaCharge(); - break; + val = (*this)[idx].GetSigmaCharge(); + break; case 3: - val = (*this)[idx].GetErrSigmaCharge(); - break; + val = (*this)[idx].GetErrSigmaCharge(); + break; case 4: - val = (*this)[idx].GetChargeProb(); - break; + val = (*this)[idx].GetChargeProb(); + break; case 5: - val = (*this)[idx].GetTime(); - break; + val = (*this)[idx].GetTime(); + break; case 6: - val = (*this)[idx].GetSigmaTime(); - break; + val = (*this)[idx].GetSigmaTime(); + break; case 7: - val = (*this)[idx].GetTimeProb(); - break; + val = (*this)[idx].GetTimeProb(); + break; case 8: - val = (*this)[idx].GetPed(); - break; + val = (*this)[idx].GetPed(); + break; case 9: - val = (*this)[idx].GetPedRms(); - break; + val = (*this)[idx].GetPedRms(); + break; case 10: - val = (*this)[idx].GetRCharge(); + val = (*this)[idx].GetRSigma(); break; case 11: - val = (*this)[idx].GetRSigma(); + val = (*this)[idx].GetPheFFactorMethod(); break; case 12: - val = (*this)[idx].GetPheFFactorMethod(); + val = (*this)[idx].GetConversionFFactorMethod(); break; case 13: - val = (*this)[idx].GetConversionFFactorMethod(); + if (idx < 397) + val = (double)fMeanNrPhotInnerPix; + else + val = (double)fMeanNrPhotInnerPix*gkCalibrationOuterPixelArea; break; case 14: - val = (double)fMeanNrPhotInnerPix; - break; - case 15: - if ((fMeanNrPhotInnerPix > 0. ) && ((*this)[idx].GetRCharge() > 100.)) - val = fMeanNrPhotInnerPix / (*this)[idx].GetRCharge(); - else - val = -1.; + if ((fMeanNrPhotInnerPix > 0. ) && ((*this)[idx].GetCharge() != -1.)) + { + if (idx < 397) + val = fMeanNrPhotInnerPix / (*this)[idx].GetCharge(); + else + val = fMeanNrPhotInnerPix*gkCalibrationOuterPixelArea / (*this)[idx].GetCharge(); + } + else + { + val = -1.; + } break; default: @@ -467,4 +502,10 @@ / gkCalibrationBlindPixelArea; break; + case kECCT1: + default: + fMeanNrPhotInnerPix = (mean / gkCalibrationBlindPixelQECT1 ) + *TMath::Power(10,gkCalibrationBlindPixelAttCT1) + / gkCalibrationBlindPixelArea; + break; } @@ -474,6 +515,5 @@ - -Bool_t MCalibrationCam::GetConversionFactor(Int_t ipx, Float_t &mean, Float_t &err) +Bool_t MCalibrationCam::GetConversionFactorBlindPixel(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma) { @@ -485,11 +525,35 @@ return kFALSE; - mean = fMeanNrPhotInnerPix / (*this)[ipx].GetRCharge(); + mean = fMeanNrPhotInnerPix / (*this)[ipx].GetCharge(); // // Not yet ready , sorry // - err = 100000000000.; + err = -1.; + sigma = -1.; return kTRUE; } + + +Bool_t MCalibrationCam::GetConversionFactorFFactor(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma) +{ + + if (ipx < 0 || !IsPixelFitted(ipx)) + return kFALSE; + + Float_t conv = (*this)[ipx].GetConversionFFactorMethod(); + + if (conv < 0.) + return kFALSE; + + mean = conv; + + // + // Not yet ready , sorry + // + err = -1.; + sigma = -1.; + + return kTRUE; +} Index: trunk/MagicSoft/Mars/manalysis/MCalibrationCam.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCalibrationCam.h (revision 2626) +++ trunk/MagicSoft/Mars/manalysis/MCalibrationCam.h (revision 2627) @@ -38,11 +38,11 @@ MCalibrationPINDiode *fPINDiode; //! containing PIN Diode data with fit results + Bool_t fMeanNrPhotAvailable; Float_t fMeanNrPhotInnerPix; // The mean number of photons in an inner pixel Float_t fMeanNrPhotInnerPixErr; // The uncertainty about the number of photons in an inner pixel - Bool_t fMeanNrPhotAvailable; public: - enum CalibrationColor_t { kECGreen, kECBlue, kECUV }; + enum CalibrationColor_t { kECGreen, kECBlue, kECUV, kECCT1 }; private: @@ -89,5 +89,6 @@ void SetColor(CalibrationColor_t color) { fColor = color; } - Bool_t GetConversionFactor(Int_t ipx, Float_t &mean, Float_t &err); + Bool_t GetConversionFactorFFactor(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma); + Bool_t GetConversionFactorBlindPixel(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma); ClassDef(MCalibrationCam, 1) // Storage Container for all calibration information of the camera Index: trunk/MagicSoft/Mars/manalysis/MCalibrationConfig.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCalibrationConfig.h (revision 2626) +++ trunk/MagicSoft/Mars/manalysis/MCalibrationConfig.h (revision 2627) @@ -15,5 +15,5 @@ // The conversion factor between High Gain and Low Gain -const UShort_t gkConversionHiLo = 10.; +const Float_t gkConversionHiLo = 10.; // The penalty constant to produce overflow in the histogram @@ -29,4 +29,5 @@ const Float_t gkCalibrationBlindPixelQEBlue = 0.226; const Float_t gkCalibrationBlindPixelQEUV = 0.247; +const Float_t gkCalibrationBlindPixelQECT1 = 0.247; // Attenuation factor Blind Pixel (three colours) @@ -34,7 +35,11 @@ const Float_t gkCalibrationBlindPixelAttBlue = 1.96; const Float_t gkCalibrationBlindPixelAttUV = 1.95; +const Float_t gkCalibrationBlindPixelAttCT1 = 1.95; // Area of Blind Pixel w.r.t. Inner Pixel const Float_t gkCalibrationBlindPixelArea = 0.25; + +// Area of Outer Pixel w.r.t. Inner Pixel +const Float_t gkCalibrationOuterPixelArea = 4.00; // ----- PIN DIODE ------------------------// Index: trunk/MagicSoft/Mars/manalysis/MCalibrationPINDiode.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCalibrationPINDiode.cc (revision 2626) +++ trunk/MagicSoft/Mars/manalysis/MCalibrationPINDiode.cc (revision 2627) @@ -79,5 +79,5 @@ Bool_t MCalibrationPINDiode::FitCharge() { - if(!fHist->FitCharge()) + if(!fHist->FitChargeHiGain()) return kFALSE; @@ -94,5 +94,5 @@ { - if(!fHist->FitTime()) + if(!fHist->FitTimeHiGain()) return kFALSE; Index: trunk/MagicSoft/Mars/manalysis/MCalibrationPINDiode.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCalibrationPINDiode.h (revision 2626) +++ trunk/MagicSoft/Mars/manalysis/MCalibrationPINDiode.h (revision 2627) @@ -33,12 +33,12 @@ void Clear(Option_t *o=""); - void SetPed(Float_t ped) { fPed = ped; } + void SetPed(Float_t ped) { fPed = ped; } void SetPedRms(Float_t pedrms) { fPedRms = pedrms; } Bool_t IsValid() const { return fRCharge >=0 || fErrRCharge >= 0; } - Bool_t FillCharge(Int_t q) { return fHist->FillCharge(q); } - Bool_t FillTime(Int_t t) { return fHist->FillTime(t); } - Bool_t FillRChargevsTime(Float_t rq, Int_t t) { return fHist->FillChargevsN(rq,t); } + Bool_t FillCharge(Float_t q) { return fHist->FillChargeHiGain(q); } + Bool_t FillTime(Int_t t) { return fHist->FillTimeHiGain(t); } + Bool_t FillRChargevsTime(Float_t rq, Int_t t) { return fHist->FillChargevsNHiGain(rq,t); } Bool_t FitCharge(); Index: trunk/MagicSoft/Mars/manalysis/MCalibrationPix.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCalibrationPix.cc (revision 2626) +++ trunk/MagicSoft/Mars/manalysis/MCalibrationPix.cc (revision 2627) @@ -50,4 +50,5 @@ fSigmaCharge(-1.), fErrSigmaCharge(-1.), + fRSigma(-1.), fChargeProb(-1.), fPed(-1.), @@ -56,10 +57,8 @@ fSigmaTime(-1.), fTimeProb(-1.), - fRCharge(-1.), - fErrRCharge(-1.), - fRSigma(-1.), fFactor(1.3), fPheFFactorMethod(-1.), - fConversionFFactorMethod(-1.) + fConversionFFactorMethod(-1.), + fHiGainSaturation(kFALSE) { @@ -102,15 +101,30 @@ if (fPed && fPedRms) - fHist->SetLowerFitRange(fPed + 1.5*fPedRms); + fHist->SetLowerFitRange(3.5*fPedRms); + // fHist->SetLowerFitRange(fPed + 3.5*fPedRms); else *fLog << warn << "Cannot set lower fit range to suppress cosmics: Pedestals not available" << endl; - if(!fHist->FitCharge()) + + if (fHiGainSaturation) { - *fLog << warn << "Could not fit charges of pixel " << fPixId << endl; - fHist->PrintChargeFitResult(); - return kFALSE; + if(!fHist->FitChargeLoGain()) + { + *fLog << warn << "Could not fit Lo Gain charges of pixel " << fPixId << endl; + fHist->PrintChargeFitResult(); + return kFALSE; + } + } + else + { + if(!fHist->FitChargeHiGain()) + { + *fLog << warn << "Could not fit Hi Gain charges of pixel " << fPixId << endl; + fHist->PrintChargeFitResult(); + return kFALSE; + } } + fCharge = fHist->GetChargeMean(); fErrCharge = fHist->GetChargeMeanErr(); @@ -122,13 +136,10 @@ { - fRCharge = fCharge - fPed; - fErrRCharge = TMath::Sqrt(fErrCharge*fErrCharge + fPedRms*fPedRms); - fRSigma = (fSigmaCharge*fSigmaCharge) - (fPedRms*fPedRms); if (fRSigma > 0. ) { - fPheFFactorMethod = fFactor * fRCharge*fRCharge / fRSigma; - fConversionFFactorMethod = fPheFFactorMethod / fRCharge ; + fPheFFactorMethod = fFactor * fCharge*fCharge / fRSigma; + fConversionFFactorMethod = fPheFFactorMethod / fCharge ; } else @@ -144,4 +155,5 @@ } + void MCalibrationPix::SetPedestal(Float_t ped, Float_t pedrms) { @@ -150,9 +162,4 @@ fPedRms = pedrms; - if ((fRCharge == -1.) && (fCharge > 0.)) - fRCharge = fCharge - fPed; - if ((fErrRCharge == -1.) && (fErrCharge > 0.)) - fErrRCharge = TMath::Sqrt(fErrCharge*fErrCharge + fPedRms*fPedRms); - } @@ -160,11 +167,23 @@ { - if(!fHist->FitTime()) + if (fHiGainSaturation) { - *fLog << warn << "Could not fit times of pixel " << fPixId << endl; - fHist->PrintTimeFitResult(); - return kFALSE; + if(!fHist->FitTimeLoGain()) + { + *fLog << warn << "Could not fit Lo Gain times of pixel " << fPixId << endl; + fHist->PrintTimeFitResult(); + return kFALSE; + } } - + else + { + if(!fHist->FitTimeHiGain()) + { + *fLog << warn << "Could not fit Hi Gain times of pixel " << fPixId << endl; + fHist->PrintTimeFitResult(); + return kFALSE; + } + } + fTime = fHist->GetTimeMean(); fSigmaTime = fHist->GetTimeSigma(); Index: trunk/MagicSoft/Mars/manalysis/MCalibrationPix.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCalibrationPix.h (revision 2626) +++ trunk/MagicSoft/Mars/manalysis/MCalibrationPix.h (revision 2627) @@ -14,8 +14,9 @@ Int_t fPixId; // the pixel Id - Float_t fCharge; // The mean charge after the fit - Float_t fErrCharge; // The error of mean charge after the fit + Float_t fCharge; // The mean reduced charge after the fit + Float_t fErrCharge; // The error of reduced mean charge after the fit Float_t fSigmaCharge; // The sigma of the mean charge after the fit Float_t fErrSigmaCharge; // The error of the sigma of the mean charge after the fit + Float_t fRSigma; // The reduced squares of sigmas after the fit Float_t fChargeProb; // The probability of the fit function @@ -27,11 +28,9 @@ Float_t fTimeProb; // The probability of the fit function - Float_t fRCharge; // The reduced mean charge after the fit - Float_t fErrRCharge; // The error of the reduced mean charge after the fit - Float_t fRSigma; // The reduced squares of sigmas after the fit - Float_t fFactor; // The laboratory F-factor Float_t fPheFFactorMethod; // The number of Phe's calculated after the F-factor method Float_t fConversionFFactorMethod; // The conversion factor to Phe's calculated after the F-factor method + + Bool_t fHiGainSaturation; // Is Lo-Gain used at all? MHCalibrationPixel *fHist; //! Pointer to the histograms performing the fits, etc. @@ -45,28 +44,31 @@ Float_t GetCharge() const { return fCharge; } - Float_t GetRCharge() const { return fRCharge; } Float_t GetRSigma() const { return fRSigma; } Float_t GetErrCharge() const { return fErrCharge; } - Float_t GetErrRCharge() const { return fErrRCharge; } Float_t GetChargeProb() const { return fChargeProb; } Float_t GetSigmaCharge() const { return fSigmaCharge; } Float_t GetErrSigmaCharge() const { return fErrSigmaCharge; } - Float_t GetTime() const { return fTime; } - Float_t GetSigmaTime() const { return fSigmaTime; } - Float_t GetTimeProb() const { return fTimeProb; } + Float_t GetTime() const { return fTime; } + Float_t GetSigmaTime() const { return fSigmaTime; } + Float_t GetTimeProb() const { return fTimeProb; } - Float_t GetPed() const { return fPed; } - Float_t GetPedRms() const { return fPedRms; } + Float_t GetPed() const { return fPed; } + Float_t GetPedRms() const { return fPedRms; } void SetPedestal(Float_t ped, Float_t pedrms); + void SetHiGainSaturation() { fHiGainSaturation = kTRUE; fHist->SetUseLoGain(); } - Bool_t FillCharge(Int_t q) { return fHist->FillCharge(q); } - Bool_t FillTime(Int_t t) { return fHist->FillTime(t); } - Bool_t FillRChargevsTime(Float_t rq, Int_t t) { return fHist->FillChargevsN(rq,t); } + Bool_t FillChargeHiGain(Float_t q) { return fHist->FillChargeHiGain(q); } + Bool_t FillTimeHiGain(Int_t t) { return fHist->FillTimeHiGain(t); } + Bool_t FillRChargevsTimeHiGain(Float_t rq, Int_t t) { return fHist->FillChargevsNHiGain(rq,t); } + + Bool_t FillChargeLoGain(Float_t q) { return fHist->FillChargeLoGain(q); } + Bool_t FillTimeLoGain(Int_t t) { return fHist->FillTimeLoGain(t); } + Bool_t FillRChargevsTimeLoGain(Float_t rq, Int_t t) { return fHist->FillChargevsNLoGain(rq,t); } - Bool_t IsValid() const { return fRCharge >=0 || fErrRCharge >= 0; } - Int_t GetPixId() const { return fPixId; } + Bool_t IsValid() const { return fCharge >=0 || fErrCharge >= 0; } + Int_t GetPixId() const { return fPixId; } void DefinePixId(Int_t i); @@ -74,6 +76,6 @@ Bool_t FitTime(); - MHCalibrationPixel *GetHist() const { return fHist; } - virtual void Draw(Option_t *opt="") { fHist->Draw(opt); } + MHCalibrationPixel *GetHist() const { return fHist; } + void Draw(Option_t *opt="") { fHist->Draw(opt); } Float_t GetPheFFactorMethod() const { return fPheFFactorMethod; } Index: trunk/MagicSoft/Mars/manalysis/MExtractSignal.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MExtractSignal.cc (revision 2627) +++ trunk/MagicSoft/Mars/manalysis/MExtractSignal.cc (revision 2627) @@ -0,0 +1,217 @@ +/* ======================================================================== *\ +! +! * +! * 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): Markus Gaug 09/2003 +! +! Copyright: MAGIC Software Development, 2000-2001 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// MExtractSignal // +// // +////////////////////////////////////////////////////////////////////////////// + +#include "MExtractSignal.h" + +#include "MExtractedSignalCam.h" +#include "MExtractedSignalPix.h" + +#include "MPedestalCam.h" +#include "MPedestalPix.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" +#include "MH.h" + +#include "MRawRunHeader.h" +#include "MRawEvtData.h" // MRawEvtData::GetNumPixels +#include "MRawEvtPixelIter.h" + +#include "TMath.h" + +ClassImp(MExtractSignal); + +using namespace std; +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MExtractSignal::MExtractSignal(const char *name, const char *title) + : fSaturationLimit(250), + fConversionHiLo(10.) +{ + + fName = name ? name : "MExtractSignal"; + fTitle = title ? title : "Task to extract the signal from the FADC slices"; + + AddToBranchList("MRawEvtData.*"); +} + +// -------------------------------------------------------------------------- +// +// The PreProcess searches for the following input containers: +// - MRawEvtData +// - MPedestalCam +// +// The following output containers are also searched and created if +// they were not found: +// +// - MExtractedSignalCam +// +Int_t MExtractSignal::PreProcess(MParList *pList) +{ + + fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData"); + if (!fRawEvt) + { + *fLog << dbginf << "MRawEvtData not found... aborting." << endl; + return kFALSE; + } + + const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); + if (!runheader) + *fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl; + + + fSignals = (MExtractedSignalCam*)pList->FindCreateObj("MExtractedSignalCam"); + if (!fSignals) + return kFALSE; + + fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam"); + if (!fPedestals) + { + *fLog << err << dbginf << "Cannot find MPedestalCam ... aborting" << endl; + return kFALSE; + } + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// The ReInit searches for the following input containers: +// - MRawRunHeader +// +Bool_t MExtractSignal::ReInit(MParList *pList ) +{ + + fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); + if (!fRunHeader) + { + *fLog << dbginf << "MRawRunHeader not found... aborting." << endl; + return kFALSE; + } + + fNumHiGainSamples = fRunHeader->GetNumSamplesHiGain(); + fNumLoGainSamples = fRunHeader->GetNumSamplesLoGain(); + fSqrtHiGainSamples = TMath::Sqrt((float)fNumHiGainSamples); + fSqrtLoGainSamples = TMath::Sqrt((float)fNumLoGainSamples); + + // + // FIXME: The next statement simply does not work: + // fRawEvt->GetNumPixels() returns always 0 + // + + fSignals->InitSize(577); + // fExtractedSignals->InitSize(fRawEvt->GetNumPixels()); + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Calculate the integral of the FADC time slices and store them as a new +// pixel in the MExtractedSignalCam container. +// +Int_t MExtractSignal::Process() +{ + + MRawEvtPixelIter pixel(fRawEvt); + fSignals->Clear(); + + while (pixel.Next()) + { + + UShort_t satHi = 0; + UShort_t satLo = 0; + Float_t sum = 0.; + Float_t sumerr = 0.; + + const Int_t pixid = pixel.GetPixelId(); + const MPedestalPix &ped = (*fPedestals)[pixid]; + MExtractedSignalPix &pix = (*fSignals)[pixid]; + + const Float_t pedes = ped.GetPedestal(); + const Float_t pedrms = ped.GetPedestalRms(); + + // + // FIXME: This is preliminary, we will change to pedestals per slice!!! + // Assume pedestals per time slice ==> multiply with number of slices + // + + UInt_t max = pixel.GetMaxHiGainSample(); + + if (max < fSaturationLimit) // take Hi Gain, no saturation + { + sum = (float)pixel.GetSumHiGainSamples() - pedes*fNumHiGainSamples; + sumerr = pedrms*fSqrtHiGainSamples; + } + else // Lo Gain + { + + satHi++; + + sum = (float)pixel.GetSumLoGainSamples() - pedes*fNumLoGainSamples ; + sum *= fConversionHiLo; + + sumerr = pedrms*fSqrtLoGainSamples; + + max = pixel.GetMaxLoGainSample(); + + if (max > fSaturationLimit) + { + *fLog << err << dbginf + << "Warning: Saturation of Lo Gain reached in pixel: " + << pixid << " " << " sum = " << sum << endl; + satLo++; + } + + } + + pix.SetGainSaturation((satHi),(satHi),(satLo)); + pix.SetExtractedSignal(sum,sumerr); + + } /* while (pixel.Next()) */ + + fSignals->SetReadyToSave(); + + return kTRUE; +} + +Int_t MExtractSignal::PostProcess() +{ + + return kTRUE; + +} Index: trunk/MagicSoft/Mars/manalysis/MExtractSignal.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MExtractSignal.h (revision 2627) +++ trunk/MagicSoft/Mars/manalysis/MExtractSignal.h (revision 2627) @@ -0,0 +1,57 @@ +#ifndef MARS_MExtractSignal +#define MARS_MExtractSignal + +///////////////////////////////////////////////////////////////////////////// +// // +// MExtractSignal // +// // +// Integrates the time slices of the all pixels of a calibration event // +// and substract the pedestal value // +// // +///////////////////////////////////////////////////////////////////////////// + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class MRawEvtData; +class MRawRunHeader; + +class MPedestalCam; +class MExtractedSignalCam; + +class MExtractSignal : public MTask +{ +private: + + MPedestalCam *fPedestals; // Pedestals of all pixels in the camera + MExtractedSignalCam *fSignals; // Extracted signal of all pixels in the camera + + MRawEvtData *fRawEvt; // raw event data (time slices) + MRawRunHeader *fRunHeader; // RunHeader information + + Byte_t fNumHiGainSamples; + Byte_t fNumLoGainSamples; + + Float_t fSqrtHiGainSamples; + Float_t fSqrtLoGainSamples; + + Byte_t fSaturationLimit; + Float_t fConversionHiLo; + + Bool_t ReInit(MParList *pList); + Int_t PreProcess(MParList *pList); + Int_t Process(); + Int_t PostProcess(); + +public: + + MExtractSignal(const char *name=NULL, const char *title=NULL); + + void SetSaturationLimit(Byte_t lim) { fSaturationLimit = lim; } + void SetConversionHiLo(Float_t con) { fConversionHiLo = con; } + + ClassDef(MExtractSignal, 0) // Task to fill the Extracted Signal Containers from raw data +}; + +#endif Index: trunk/MagicSoft/Mars/manalysis/MExtractedSignalCam.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MExtractedSignalCam.cc (revision 2627) +++ trunk/MagicSoft/Mars/manalysis/MExtractedSignalCam.cc (revision 2627) @@ -0,0 +1,149 @@ +/* ======================================================================== *\ +! +! * +! * 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): Markus Gaug 12/2003 +! Author(s): Thomas Bretz 12/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MExtractedSignalCam // +// // +// Hold the Extracted Signal information for all pixels in the camera // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MExtractedSignalCam.h" +#include "MExtractedSignalPix.h" + +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MGeomCam.h" + +ClassImp(MExtractedSignalCam); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor. Creates a MExtractedSignalPix object for each pixel +// +MExtractedSignalCam::MExtractedSignalCam(const char *name, const char *title) +{ + fName = name ? name : "MExtractedSignalCam"; + fTitle = title ? title : "Storage container for all Extracted Signal Information in the camera"; + + fArray = new TClonesArray("MExtractedSignalPix", 1); + +} + +// -------------------------------------------------------------------------- +// +// Delete the array conatining the pixel pedest information +// +MExtractedSignalCam::~MExtractedSignalCam() +{ + delete fArray; +} + +// -------------------------------------------------------------------------- +// +// Set the size of the camera +// +void MExtractedSignalCam::InitSize(const UInt_t i) +{ + fArray->ExpandCreate(i); +} + +// -------------------------------------------------------------------------- +// +// Get the size of the MExtractedSignalCam +// +Int_t MExtractedSignalCam::GetSize() const +{ + return fArray->GetEntriesFast(); +} + +// -------------------------------------------------------------------------- +// +// Get i-th pixel (pixel number) +// +MExtractedSignalPix &MExtractedSignalCam::operator[](Int_t i) +{ + return *static_cast(fArray->UncheckedAt(i)); +} + +// -------------------------------------------------------------------------- +// +// Get i-th pixel (pixel number) +// +MExtractedSignalPix &MExtractedSignalCam::operator[](Int_t i) const +{ + return *static_cast(fArray->UncheckedAt(i)); +} + + +void MExtractedSignalCam::Clear(Option_t *o) +{ + fArray->ForEach(TObject, Clear)(); +} + +void MExtractedSignalCam::Print(Option_t *o) const +{ + *fLog << all << GetDescriptor() << ":" << endl; + int id = 0; + + TIter Next(fArray); + MExtractedSignalPix *pix; + while ((pix=(MExtractedSignalPix*)Next())) + { + id++; + + if (!pix->IsValid()) + continue; + + *fLog << id-1 << ": " << endl; + pix->PrintOut(); + } +} + +Bool_t MExtractedSignalCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const +{ + switch (type) + { + case 0: + val = (*this)[idx].GetExtractedSignal(); + break; + case 1: + val = (*this)[idx].GetExtractedSignalError(); + break; + default: + return kFALSE; + } + return val>=0; +} + +void MExtractedSignalCam::DrawPixelContent(Int_t num) const +{ + *fLog << warn << "MExtractedSignaCam::DrawPixelContent - not available." << endl; +} Index: trunk/MagicSoft/Mars/manalysis/MExtractedSignalCam.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MExtractedSignalCam.h (revision 2627) +++ trunk/MagicSoft/Mars/manalysis/MExtractedSignalCam.h (revision 2627) @@ -0,0 +1,40 @@ +#ifndef MARS_MExtractedSignalCam +#define MARS_MExtractedSignalCam + +#ifndef MARS_MCamEvent +#include "MCamEvent.h" +#endif + +class TClonesArray; + +class MGeomCam; +class MExtractedSignalPix; + +class MExtractedSignalCam : public MCamEvent +{ +private: + TClonesArray *fArray; // FIXME: Change TClonesArray away from a pointer? + +public: + + MExtractedSignalCam(const char *name=NULL, const char *title=NULL); + ~MExtractedSignalCam(); + + void Clear(Option_t *o=""); + + void InitSize(const UInt_t i); + Int_t GetSize() const; + + MExtractedSignalPix &operator[](Int_t i); + MExtractedSignalPix &operator[](Int_t i) const; + + void Print(Option_t *o="") const; + + Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const; + void DrawPixelContent(Int_t num) const; + + ClassDef(MExtractedSignalCam, 0) // Storage Container for extracted signals in the camera +}; + +#endif + Index: trunk/MagicSoft/Mars/manalysis/MExtractedSignalPix.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MExtractedSignalPix.cc (revision 2627) +++ trunk/MagicSoft/Mars/manalysis/MExtractedSignalPix.cc (revision 2627) @@ -0,0 +1,85 @@ +/* ======================================================================== *\ +! +! * +! * 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): Markus Gaug 12/2003 +! Author(s): Thomas Bretz 12/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MExtractedSignalPix // +// // +// This is the storage container to hold informations about the pedestal // +// (offset) value of one Pixel (PMT). // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MExtractedSignalPix.h" +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MExtractedSignalPix); + +using namespace std; + +MExtractedSignalPix::MExtractedSignalPix(const char* name, const char* title) + : fExtractedSignal(-1.), + fExtractedSignalError(-1.), + fIsLoGainUsed(kFALSE), + fNumHiGainSaturated(0), + fNumLoGainSaturated(0) +{ + + fName = name ? name : "MExtractedSignalPix"; + fTitle = title ? title : "Container of the Extracted Signals"; + +} + + + +// ------------------------------------------------------------------------ +// +// Invalidate values +// +void MExtractedSignalPix::Clear(Option_t *o) +{ + + fExtractedSignal = -1.; + fExtractedSignalError = -1.; + fIsLoGainUsed = kFALSE; + fNumHiGainSaturated = 0; + fNumLoGainSaturated = 0; + +} + +void MExtractedSignalPix::PrintOut() +{ + + *fLog << all << GetDescriptor() << ":" << endl; + + *fLog << " Signal: " << fExtractedSignal + << " +- " << fExtractedSignalError + << " LoGain? " << fIsLoGainUsed + << " Nr. Sat. Hi Gain: " << fNumHiGainSaturated + << " Nr. Sat. Lo Gain: " << fNumLoGainSaturated + << endl; + +} + Index: trunk/MagicSoft/Mars/manalysis/MExtractedSignalPix.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MExtractedSignalPix.h (revision 2627) +++ trunk/MagicSoft/Mars/manalysis/MExtractedSignalPix.h (revision 2627) @@ -0,0 +1,43 @@ +#ifndef MARS_MExtractedSignalPix +#define MARS_MExtractedSignalPix + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +class MExtractedSignalPix : public MParContainer +{ +private: + + Float_t fExtractedSignal; // mean value of the extracted signal + Float_t fExtractedSignalError; // error of the mean value of the extracted signal + Bool_t fIsLoGainUsed; + Byte_t fNumHiGainSaturated; + Byte_t fNumLoGainSaturated; + +public: + + MExtractedSignalPix(const char* name=NULL, const char* title=NULL); + + void Clear(Option_t *o=""); + + Float_t GetExtractedSignal() const { return fExtractedSignal; } + Float_t GetExtractedSignalError() const { return fExtractedSignalError; } + + void SetExtractedSignal(Float_t sig, Float_t sigerr) { fExtractedSignal = sig; fExtractedSignalError = sigerr; } + void SetGainSaturation(Bool_t sat, Byte_t higain, Byte_t logain) + { fIsLoGainUsed = sat; + fNumHiGainSaturated = higain; + fNumLoGainSaturated = logain; } + + Bool_t IsLoGainUsed() { return fIsLoGainUsed; } + + Bool_t IsValid() const { return fExtractedSignal >= 0 || fExtractedSignalError >= 0; } + + void PrintOut(); + + ClassDef(MExtractedSignalPix, 0) // Storage Container for Extracted Signal information of one pixel +}; + +#endif + Index: trunk/MagicSoft/Mars/manalysis/Makefile =================================================================== --- trunk/MagicSoft/Mars/manalysis/Makefile (revision 2626) +++ trunk/MagicSoft/Mars/manalysis/Makefile (revision 2627) @@ -75,5 +75,8 @@ MCalibrationBlindPix.cc \ MCalibrationPINDiode.cc \ - MCalibrationCam.cc + MCalibrationCam.cc \ + MExtractedSignalCam.cc \ + MExtractedSignalPix.cc \ + MExtractSignal.cc SRCS = $(SRCFILES) Index: trunk/MagicSoft/Mars/mhist/MHCalibrationBlindPixel.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHCalibrationBlindPixel.cc (revision 2626) +++ trunk/MagicSoft/Mars/mhist/MHCalibrationBlindPixel.cc (revision 2627) @@ -69,9 +69,10 @@ // Create a large number of bins, later we will rebin - fBlindPixelChargefirst = 0; + fBlindPixelChargefirst = -1000.; fBlindPixelChargelast = gkStartBlindPixelBinNr; - fBlindPixelChargenbins = gkStartBlindPixelBinNr; - - fHBlindPixelCharge = new TH1I("HBlindPixelCharge","Distribution of Summed FADC Slices",fBlindPixelChargenbins,fBlindPixelChargefirst,fBlindPixelChargelast); + fBlindPixelChargenbins = gkStartBlindPixelBinNr+(int)fBlindPixelChargefirst; + + fHBlindPixelCharge = new TH1F("HBlindPixelCharge","Distribution of Summed FADC Slices", + fBlindPixelChargenbins,fBlindPixelChargefirst,fBlindPixelChargelast); fHBlindPixelCharge->SetXTitle("Sum FADC Slices"); fHBlindPixelCharge->SetYTitle("Nr. of events"); @@ -313,5 +314,5 @@ const Double_t lambda_guess = 0.2; const Double_t mu_0_guess = fHBlindPixelCharge->GetBinCenter(fHBlindPixelCharge->GetMaximumBin()); - const Double_t si_0_guess = mu_0_guess/10.; + const Double_t si_0_guess = mu_0_guess/5.; const Double_t mu_1_guess = mu_0_guess + 50.; const Double_t si_1_guess = si_0_guess + si_0_guess; Index: trunk/MagicSoft/Mars/mhist/MHCalibrationBlindPixel.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHCalibrationBlindPixel.h (revision 2626) +++ trunk/MagicSoft/Mars/mhist/MHCalibrationBlindPixel.h (revision 2627) @@ -32,5 +32,5 @@ private: - TH1I* fHBlindPixelCharge; //-> Histogram with the single Phe spectrum + TH1F* fHBlindPixelCharge; //-> Histogram with the single Phe spectrum TH1F* fHBlindPixelErrCharge; //-> Variance of summed FADC slices TH1I* fHBlindPixelTime; //-> Variance of summed FADC slices @@ -83,7 +83,7 @@ ~MHCalibrationBlindPixel(); - Bool_t FillBlindPixelCharge(Int_t q) { return fHBlindPixelCharge->Fill(q) > -1; } - Bool_t FillErrBlindPixelCharge(Float_t errq) { return fHBlindPixelErrCharge->Fill(errq) > -1; } - Bool_t FillBlindPixelTime(Int_t t) { return fHBlindPixelTime->Fill(t) > -1; } + Bool_t FillBlindPixelCharge(Float_t q) { return fHBlindPixelCharge->Fill(q) > -1; } + Bool_t FillErrBlindPixelCharge(Float_t errq) { return fHBlindPixelErrCharge->Fill(errq) > -1; } + Bool_t FillBlindPixelTime(Int_t t) { return fHBlindPixelTime->Fill(t) > -1; } Bool_t FillBlindPixelChargevsN(Stat_t rq, Int_t t) { return fHBlindPixelChargevsN->Fill(t,rq) > -1; } @@ -123,9 +123,8 @@ void ChangeFitFunc(BlindPixelFitFunc fitfunc, Int_t par=5); - void CutAllEdges(); void Draw(Option_t *option=""); - ClassDef(MHCalibrationBlindPixel, 0) + ClassDef(MHCalibrationBlindPixel, 1) }; Index: trunk/MagicSoft/Mars/mhist/MHCalibrationConfig.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHCalibrationConfig.h (revision 2626) +++ trunk/MagicSoft/Mars/mhist/MHCalibrationConfig.h (revision 2627) @@ -7,5 +7,5 @@ // // // Contains all configuration data of the Calibration // -// // +// i // ///////////////////////////////////////////////////////////////////////////// Index: trunk/MagicSoft/Mars/mhist/MHCalibrationPINDiode.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHCalibrationPINDiode.cc (revision 2626) +++ trunk/MagicSoft/Mars/mhist/MHCalibrationPINDiode.cc (revision 2627) @@ -67,9 +67,10 @@ // Create a large number of bins, later we will rebin - fChargeFirst = 0; - fChargeLast = gkStartPINDiodeBinNr; - fChargeNbins = gkStartPINDiodeBinNr; + fChargeFirstHiGain = -1000.; + fChargeLastHiGain = gkStartPINDiodeBinNr; + fChargeNbinsHiGain = gkStartPINDiodeBinNr; - fHPCharge = new TH1I("HPCharge","Distribution of Summed FADC Slices",fChargeNbins,fChargeFirst,fChargeLast); + fHPCharge = new TH1I("HPCharge","Distribution of Summed FADC Slices", + fChargeNbinsHiGain,fChargeFirstHiGain,fChargeLastHiGain); fHPCharge->SetXTitle("Sum FADC Slices"); fHPCharge->SetYTitle("Nr. of events"); Index: trunk/MagicSoft/Mars/mhist/MHCalibrationPixel.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHCalibrationPixel.cc (revision 2626) +++ trunk/MagicSoft/Mars/mhist/MHCalibrationPixel.cc (revision 2627) @@ -65,5 +65,11 @@ fTimeGausFit(NULL), fFitLegend(NULL), - fLowerFitRange(0.), + fLowerFitRange(-2000.), + fChargeFirstHiGain(-2000.5), + fChargeLastHiGain(9999.5), + fChargeNbinsHiGain(12000), + fChargeFirstLoGain(-2000.5), + fChargeLastLoGain(9999.5), + fChargeNbinsLoGain(1200), fFitOK(kFALSE), fChargeChisquare(-1.), @@ -72,5 +78,6 @@ fTimeChisquare(-1.), fTimeProb(-1.), - fTimeNdf(-1) + fTimeNdf(-1), + fUseLoGain(kFALSE) { @@ -79,15 +86,19 @@ // Create a large number of bins, later we will rebin - fChargeFirst = -0.5; - fChargeLast = 10000. - 0.5; - fChargeNbins = 20000; - - fHCharge = new TH1I("HCharge","Distribution of Summed FADC Slices Pixel ", - fChargeNbins,fChargeFirst,fChargeLast); - fHCharge->SetXTitle("Sum FADC Slices"); - fHCharge->SetYTitle("Nr. of events"); - fHCharge->Sumw2(); - - fHCharge->SetDirectory(NULL); + fHChargeHiGain = new TH1F("HChargeHiGain","Distribution of Summed FADC Hi Gain Slices Pixel ", + fChargeNbinsHiGain,fChargeFirstHiGain,fChargeLastHiGain); + fHChargeHiGain->SetXTitle("Sum FADC Slices (Hi Gain)"); + fHChargeHiGain->SetYTitle("Nr. of events"); + fHChargeHiGain->Sumw2(); + + fHChargeHiGain->SetDirectory(NULL); + + fHChargeLoGain = new TH1F("HChargeLoGain","Distribution of Summed FADC Lo Gain Slices Pixel ", + fChargeNbinsLoGain,fChargeFirstLoGain,fChargeLastLoGain); + fHChargeLoGain->SetXTitle("Sum FADC Slices (Lo Gain)"); + fHChargeLoGain->SetYTitle("Nr. of events"); + fHChargeLoGain->Sumw2(); + + fHChargeLoGain->SetDirectory(NULL); Axis_t tfirst = -0.5; @@ -95,13 +106,19 @@ Int_t ntbins = 16; - fHTime = new TH1I("HTime","Distribution of Mean Arrival Times Pixel ", - ntbins,tfirst,tlast); - fHTime->SetXTitle("Mean Arrival Times [FADC slice nr]"); - fHTime->SetYTitle("Nr. of events"); - fHTime->Sumw2(); - - fHTime->SetDirectory(NULL); - - TString qvsntitle = "Sum of Charges vs. Event Number Pixel "; + fHTimeHiGain = new TH1I("HTimeHiGain","Distribution of Mean Arrival Hi Gain Times Pixel ", + ntbins,tfirst,tlast); + fHTimeHiGain->SetXTitle("Mean Arrival Times [Hi Gain FADC slice nr]"); + fHTimeHiGain->SetYTitle("Nr. of events"); + fHTimeHiGain->Sumw2(); + + fHTimeHiGain->SetDirectory(NULL); + + fHTimeLoGain = new TH1I("HTimeLoGain","Distribution of Mean Arrival Lo Gain Times Pixel ", + ntbins,tfirst,tlast); + fHTimeLoGain->SetXTitle("Mean Arrival Times [Lo Gain FADC slice nr]"); + fHTimeLoGain->SetYTitle("Nr. of events"); + fHTimeLoGain->Sumw2(); + + fHTimeLoGain->SetDirectory(NULL); // We define a reasonable number and later enlarge it if necessary @@ -110,10 +127,17 @@ Axis_t nlast = (Axis_t)nqbins - 0.5; - fHChargevsN = new TH1I("HChargevsN",qvsntitle.Data(), - nqbins,nfirst,nlast); - fHChargevsN->SetXTitle("Event Nr."); - fHChargevsN->SetYTitle("Sum of FADC slices"); - - fHChargevsN->SetDirectory(NULL); + fHChargevsNHiGain = new TH1I("HChargevsNHiGain","Sum of Hi Gain Charges vs. Event Number Pixel ", + nqbins,nfirst,nlast); + fHChargevsNHiGain->SetXTitle("Event Nr."); + fHChargevsNHiGain->SetYTitle("Sum of Hi Gain FADC slices"); + + fHChargevsNHiGain->SetDirectory(NULL); + + fHChargevsNLoGain = new TH1I("HChargevsNLoGain","Sum of Lo Gain Charges vs. Event Number Pixel ", + nqbins,nfirst,nlast); + fHChargevsNLoGain->SetXTitle("Event Nr."); + fHChargevsNLoGain->SetYTitle("Sum of Lo Gain FADC slices"); + + fHChargevsNLoGain->SetDirectory(NULL); } @@ -122,7 +146,11 @@ { - delete fHCharge; - delete fHTime; - delete fHChargevsN; + delete fHChargeHiGain; + delete fHTimeHiGain; + delete fHChargevsNHiGain; + + delete fHChargeLoGain; + delete fHTimeLoGain; + delete fHChargevsNLoGain; if (fChargeGausFit) @@ -141,27 +169,52 @@ fPixId = id; - TString nameQ = TString(fHCharge->GetName()); - nameQ += id; - fHCharge->SetName(nameQ.Data()); - - TString nameT = TString(fHTime->GetName()); - nameT += id; - fHTime->SetName(nameT.Data()); - - TString nameQvsN = TString(fHChargevsN->GetName()); - nameQvsN += id; - fHChargevsN->SetName(nameQvsN.Data()); - - TString titleQ = TString(fHCharge->GetTitle()); - titleQ += id; - fHCharge->SetTitle(titleQ.Data()); - - TString titleT = TString(fHTime->GetTitle()); - titleT += id; - fHTime->SetTitle(titleT.Data()); - - TString titleQvsN = TString(fHChargevsN->GetTitle()); - titleQvsN += id; - fHChargevsN->SetTitle(titleQvsN.Data()); + TString nameQHiGain = TString(fHChargeHiGain->GetName()); + nameQHiGain += id; + fHChargeHiGain->SetName(nameQHiGain.Data()); + + TString nameTHiGain = TString(fHTimeHiGain->GetName()); + nameTHiGain += id; + fHTimeHiGain->SetName(nameTHiGain.Data()); + + TString nameQvsNHiGain = TString(fHChargevsNHiGain->GetName()); + nameQvsNHiGain += id; + fHChargevsNHiGain->SetName(nameQvsNHiGain.Data()); + + TString titleQHiGain = TString(fHChargeHiGain->GetTitle()); + titleQHiGain += id; + fHChargeHiGain->SetTitle(titleQHiGain.Data()); + + TString titleTHiGain = TString(fHTimeHiGain->GetTitle()); + titleTHiGain += id; + fHTimeHiGain->SetTitle(titleTHiGain.Data()); + + TString titleQvsNHiGain = TString(fHChargevsNHiGain->GetTitle()); + titleQvsNHiGain += id; + fHChargevsNHiGain->SetTitle(titleQvsNHiGain.Data()); + + + TString nameQLoGain = TString(fHChargeLoGain->GetName()); + nameQLoGain += id; + fHChargeLoGain->SetName(nameQLoGain.Data()); + + TString nameTLoGain = TString(fHTimeLoGain->GetName()); + nameTLoGain += id; + fHTimeLoGain->SetName(nameTLoGain.Data()); + + TString nameQvsNLoGain = TString(fHChargevsNLoGain->GetName()); + nameQvsNLoGain += id; + fHChargevsNLoGain->SetName(nameQvsNLoGain.Data()); + + TString titleQLoGain = TString(fHChargeLoGain->GetTitle()); + titleQLoGain += id; + fHChargeLoGain->SetTitle(titleQLoGain.Data()); + + TString titleTLoGain = TString(fHTimeLoGain->GetTitle()); + titleTLoGain += id; + fHTimeLoGain->SetTitle(titleTLoGain.Data()); + + TString titleQvsNLoGain = TString(fHChargevsNLoGain->GetTitle()); + titleQvsNLoGain += id; + fHChargevsNLoGain->SetTitle(titleQvsNLoGain.Data()); } @@ -170,14 +223,25 @@ { - for (Int_t i = fHCharge->FindBin(fChargeFirst); - i <= fHCharge->FindBin(fChargeLast); i++) - fHCharge->SetBinContent(i, 1.e-20); + for (Int_t i = fHChargeHiGain->FindBin(fChargeFirstHiGain); + i <= fHChargeHiGain->FindBin(fChargeLastHiGain); i++) + fHChargeHiGain->SetBinContent(i, 1.e-20); for (Int_t i = 0; i < 16; i++) - fHTime->SetBinContent(i, 1.e-20); - - fChargeLast = 9999.5; - - fHCharge->GetXaxis()->SetRangeUser(0.,fChargeLast); + fHTimeHiGain->SetBinContent(i, 1.e-20); + + fChargeLastHiGain = 9999.5; + + fHChargeLoGain->GetXaxis()->SetRangeUser(0.,fChargeLastLoGain); + + for (Int_t i = fHChargeLoGain->FindBin(fChargeFirstLoGain); + i <= fHChargeLoGain->FindBin(fChargeLastLoGain); i++) + fHChargeLoGain->SetBinContent(i, 1.e-20); + + for (Int_t i = 0; i < 16; i++) + fHTimeLoGain->SetBinContent(i, 1.e-20); + + fChargeLastLoGain = 9999.5; + + fHChargeLoGain->GetXaxis()->SetRangeUser(0.,fChargeLastLoGain); return; @@ -192,6 +256,8 @@ { - fHCharge->Reset(); - fHTime->Reset(); + fHChargeHiGain->Reset(); + fHTimeHiGain->Reset(); + fHChargeLoGain->Reset(); + fHTimeLoGain->Reset(); return kTRUE; @@ -240,5 +306,4 @@ fFitLegend->AddText(line8); - if (fFitOK) fFitLegend->AddText("Result of the Fit: OK"); @@ -266,5 +331,5 @@ gROOT->SetSelectedPad(NULL); - c->Divide(2,2); + c->Divide(2,4); c->cd(1); @@ -273,51 +338,115 @@ gPad->SetTicks(); - fHCharge->Draw(opt); - - if (fChargeGausFit) - { - if (fFitOK) - fChargeGausFit->SetLineColor(kGreen); - else - fChargeGausFit->SetLineColor(kRed); - - fChargeGausFit->Draw("same"); - c->Modified(); - c->Update(); - } - - - c->cd(2); - DrawLegend(); - c->Update(); - - c->cd(3); - gStyle->SetOptStat(1111111); - - gPad->SetLogy(1); - fHTime->Draw(opt); - - if (fTimeGausFit) - { - if (fTimeChisquare > 1.) - fTimeGausFit->SetLineColor(kRed); - else - fTimeGausFit->SetLineColor(kGreen); - - fTimeGausFit->Draw("same"); - c->Modified(); - c->Update(); - } - + fHChargeHiGain->DrawCopy(opt); + c->Modified(); c->Update(); - c->cd(4); - fHChargevsN->Draw(opt); -} - - - -Bool_t MHCalibrationPixel::FitTime(Axis_t rmin, Axis_t rmax, Option_t *option) + if (fUseLoGain) + { + c->cd(2); + gPad->SetLogy(1); + gPad->SetTicks(); + fHChargeLoGain->DrawCopy(opt); + + if (fChargeGausFit) + { + if (fFitOK) + fChargeGausFit->SetLineColor(kGreen); + else + fChargeGausFit->SetLineColor(kRed); + + fChargeGausFit->DrawCopy("same"); + } + + c->cd(4); + DrawLegend(); + + } + else + { + if (fChargeGausFit) + { + if (fFitOK) + fChargeGausFit->SetLineColor(kGreen); + else + fChargeGausFit->SetLineColor(kRed); + + fChargeGausFit->DrawCopy("same"); + } + c->cd(2); + gPad->SetLogy(1); + gPad->SetTicks(); + fHChargeLoGain->DrawCopy(opt); + c->cd(3); + DrawLegend(); + } + + c->Modified(); + c->Update(); + + c->cd(5); + gStyle->SetOptStat(1111111); + + gPad->SetLogy(1); + fHTimeHiGain->DrawCopy(opt); + c->Modified(); + c->Update(); + + if (fUseLoGain) + { + + c->cd(6); + gPad->SetLogy(1); + fHTimeLoGain->DrawCopy(opt); + c->Modified(); + c->Update(); + + if (fTimeGausFit) + { + if (fTimeChisquare > 1.) + fTimeGausFit->SetLineColor(kRed); + else + fTimeGausFit->SetLineColor(kGreen); + + fTimeGausFit->DrawCopy("same"); + } + } + else + { + if (fTimeGausFit) + { + if (fTimeChisquare > 1.) + fTimeGausFit->SetLineColor(kRed); + else + fTimeGausFit->SetLineColor(kGreen); + + fTimeGausFit->DrawCopy("same"); + c->Modified(); + c->Update(); + } + + c->cd(6); + gPad->SetLogy(1); + fHTimeLoGain->DrawCopy(opt); + } + c->Modified(); + c->Update(); + + c->cd(7); + fHChargevsNHiGain->DrawCopy(opt); + c->Modified(); + c->Update(); + + c->cd(8); + fHChargevsNLoGain->DrawCopy(opt); + c->Modified(); + c->Update(); + +} + + + +Bool_t MHCalibrationPixel::FitTimeHiGain(Axis_t rmin, Axis_t rmax, Option_t *option) { @@ -328,6 +457,6 @@ rmax = (rmax != 0.) ? rmax : 9.; - const Stat_t entries = fHTime->GetEntries(); - const Double_t mu_guess = fHTime->GetBinCenter(fHTime->GetMaximumBin()); + const Stat_t entries = fHTimeHiGain->GetEntries(); + const Double_t mu_guess = fHTimeHiGain->GetBinCenter(fHTimeHiGain->GetMaximumBin()); const Double_t sigma_guess = (rmax - rmin)/2.; const Double_t area_guess = entries/gkSq2Pi; @@ -350,5 +479,5 @@ fTimeGausFit->SetRange(rmin,rmax); - fHTime->Fit(fTimeGausFit,option); + fHTimeHiGain->Fit(fTimeGausFit,option); rmin = fTimeGausFit->GetParameter(1) - 3.*fTimeGausFit->GetParameter(2); @@ -356,5 +485,5 @@ fTimeGausFit->SetRange(rmin,rmax); - fHTime->Fit(fTimeGausFit,option); + fHTimeHiGain->Fit(fTimeGausFit,option); fTimeChisquare = fTimeGausFit->GetChisquare(); @@ -375,5 +504,61 @@ } -Bool_t MHCalibrationPixel::FitCharge(Option_t *option) +Bool_t MHCalibrationPixel::FitTimeLoGain(Axis_t rmin, Axis_t rmax, Option_t *option) +{ + + if (fTimeGausFit) + return kFALSE; + + rmin = (rmin != 0.) ? rmin : 4.; + rmax = (rmax != 0.) ? rmax : 9.; + + const Stat_t entries = fHTimeLoGain->GetEntries(); + const Double_t mu_guess = fHTimeLoGain->GetBinCenter(fHTimeLoGain->GetMaximumBin()); + const Double_t sigma_guess = (rmax - rmin)/2.; + const Double_t area_guess = entries/gkSq2Pi; + + TString name = TString("GausTime"); + name += fPixId; + fTimeGausFit = new TF1(name.Data(),"gaus",rmin,rmax); + + if (!fTimeGausFit) + { + *fLog << err << dbginf << "Could not create fit function for Gauss fit" << endl; + return kFALSE; + } + + fTimeGausFit->SetParameters(area_guess,mu_guess,sigma_guess); + fTimeGausFit->SetParNames("Area","#mu","#sigma"); + fTimeGausFit->SetParLimits(0,0.,entries); + fTimeGausFit->SetParLimits(1,rmin,rmax); + fTimeGausFit->SetParLimits(2,0.,(rmax-rmin)); + fTimeGausFit->SetRange(rmin,rmax); + + fHTimeLoGain->Fit(fTimeGausFit,option); + + rmin = fTimeGausFit->GetParameter(1) - 3.*fTimeGausFit->GetParameter(2); + rmax = fTimeGausFit->GetParameter(1) + 3.*fTimeGausFit->GetParameter(2); + fTimeGausFit->SetRange(rmin,rmax); + + fHTimeLoGain->Fit(fTimeGausFit,option); + + fTimeChisquare = fTimeGausFit->GetChisquare(); + fTimeNdf = fTimeGausFit->GetNDF(); + fTimeProb = fTimeGausFit->GetProb(); + + fTimeMean = fTimeGausFit->GetParameter(1); + fTimeSigma = fTimeGausFit->GetParameter(2); + + if (fTimeChisquare > 1.) + { + *fLog << warn << "Fit of the Arrival times failed ! " << endl; + return kFALSE; + } + + return kTRUE; + +} + +Bool_t MHCalibrationPixel::FitChargeHiGain(Option_t *option) { @@ -384,5 +569,5 @@ // Get the fitting ranges // - Axis_t rmin = (fLowerFitRange != 0.) ? fLowerFitRange : fChargeFirst; + Axis_t rmin = (fLowerFitRange != 0.) ? fLowerFitRange : fChargeFirstHiGain; Axis_t rmax = 0.; @@ -391,7 +576,7 @@ // otherwise the fit goes gaga because of high number of dimensions ... // - const Stat_t entries = fHCharge->GetEntries(); + const Stat_t entries = fHChargeHiGain->GetEntries(); const Double_t area_guess = entries/gkSq2Pi; - const Double_t mu_guess = fHCharge->GetBinCenter(fHCharge->GetMaximumBin()); + const Double_t mu_guess = fHChargeHiGain->GetBinCenter(fHChargeHiGain->GetMaximumBin()); const Double_t sigma_guess = mu_guess/15.; @@ -399,5 +584,5 @@ name += fPixId; - fChargeGausFit = new TF1(name.Data(),"gaus",rmin,fChargeLast); + fChargeGausFit = new TF1(name.Data(),"gaus",rmin,fChargeLastHiGain); if (!fChargeGausFit) @@ -410,9 +595,9 @@ fChargeGausFit->SetParNames("Area","#mu","#sigma"); fChargeGausFit->SetParLimits(0,0.,entries); - fChargeGausFit->SetParLimits(1,rmin,fChargeLast); - fChargeGausFit->SetParLimits(2,0.,fChargeLast-rmin); - fChargeGausFit->SetRange(rmin,fChargeLast); - - fHCharge->Fit(fChargeGausFit,option); + fChargeGausFit->SetParLimits(1,rmin,fChargeLastHiGain); + fChargeGausFit->SetParLimits(2,0.,fChargeLastHiGain-rmin); + fChargeGausFit->SetRange(rmin,fChargeLastHiGain); + + fHChargeHiGain->Fit(fChargeGausFit,option); Axis_t rtry = fChargeGausFit->GetParameter(1) - 2.5*fChargeGausFit->GetParameter(2); @@ -422,5 +607,5 @@ fChargeGausFit->SetRange(rmin,rmax); - fHCharge->Fit(fChargeGausFit,option); + fHChargeHiGain->Fit(fChargeGausFit,option); // rmin = fChargeGausFit->GetParameter(1) - 2.5*fChargeGausFit->GetParameter(2); @@ -428,5 +613,5 @@ // fChargeGausFit->SetRange(rmin,rmax); - // fHCharge->Fit(fChargeGausFit,option); + // fHChargeHiGain->Fit(fChargeGausFit,option); fChargeChisquare = fChargeGausFit->GetChisquare(); @@ -455,17 +640,107 @@ } + +Bool_t MHCalibrationPixel::FitChargeLoGain(Option_t *option) +{ + + if (fChargeGausFit) + return kFALSE; + + // + // Get the fitting ranges + // + Axis_t rmin = (fLowerFitRange != 0.) ? fLowerFitRange : fChargeFirstLoGain; + Axis_t rmax = 0.; + + // + // First guesses for the fit (should be as close to reality as possible, + // otherwise the fit goes gaga because of high number of dimensions ... + // + const Stat_t entries = fHChargeLoGain->GetEntries(); + const Double_t area_guess = entries/gkSq2Pi; + const Double_t mu_guess = fHChargeLoGain->GetBinCenter(fHChargeLoGain->GetMaximumBin()); + const Double_t sigma_guess = mu_guess/15.; + + TString name = TString("ChargeGausFit"); + name += fPixId; + + fChargeGausFit = new TF1(name.Data(),"gaus",rmin,fChargeLastLoGain); + + if (!fChargeGausFit) + { + *fLog << err << dbginf << "Could not create fit function for Gauss fit" << endl; + return kFALSE; + } + + fChargeGausFit->SetParameters(area_guess,mu_guess,sigma_guess); + fChargeGausFit->SetParNames("Area","#mu","#sigma"); + fChargeGausFit->SetParLimits(0,0.,entries); + fChargeGausFit->SetParLimits(1,rmin,fChargeLastLoGain); + fChargeGausFit->SetParLimits(2,0.,fChargeLastLoGain-rmin); + fChargeGausFit->SetRange(rmin,fChargeLastLoGain); + + fHChargeLoGain->Fit(fChargeGausFit,option); + + Axis_t rtry = fChargeGausFit->GetParameter(1) - 2.5*fChargeGausFit->GetParameter(2); + + rmin = (rtry < rmin ? rmin : rtry); + rmax = fChargeGausFit->GetParameter(1) + 2.5*fChargeGausFit->GetParameter(2); + fChargeGausFit->SetRange(rmin,rmax); + + fHChargeLoGain->Fit(fChargeGausFit,option); + + // rmin = fChargeGausFit->GetParameter(1) - 2.5*fChargeGausFit->GetParameter(2); + // rmax = fChargeGausFit->GetParameter(1) + 2.5*fChargeGausFit->GetParameter(2); + // fChargeGausFit->SetRange(rmin,rmax); + + // fHChargeLoGain->Fit(fChargeGausFit,option); + + fChargeChisquare = fChargeGausFit->GetChisquare(); + fChargeNdf = fChargeGausFit->GetNDF(); + fChargeProb = fChargeGausFit->GetProb(); + fChargeMean = fChargeGausFit->GetParameter(1); + fChargeMeanErr = fChargeGausFit->GetParError(1); + fChargeSigma = fChargeGausFit->GetParameter(2); + fChargeSigmaErr = fChargeGausFit->GetParError(2); + + // + // The fit result is accepted under condition + // The Probability is greater than gkProbLimit (default 0.01 == 99%) + // + if (fChargeProb < gkProbLimit) + { + *fLog << warn << "Prob: " << fChargeProb << " is smaller than the allowed value: " << gkProbLimit << endl; + fFitOK = kFALSE; + return kFALSE; + } + + + fFitOK = kTRUE; + + return kTRUE; +} + void MHCalibrationPixel::CutAllEdges() { - Int_t nbins = 50; - - CutEdges(fHCharge,nbins); - - fChargeFirst = fHCharge->GetBinLowEdge(fHCharge->GetXaxis()->GetFirst()); - fChargeLast = fHCharge->GetBinLowEdge(fHCharge->GetXaxis()->GetLast())+fHCharge->GetBinWidth(0); - fChargeNbins = nbins; - - CutEdges(fHChargevsN,0); + Int_t nbins = 30; + + CutEdges(fHChargeHiGain,nbins); + + fChargeFirstHiGain = fHChargeHiGain->GetBinLowEdge(fHChargeHiGain->GetXaxis()->GetFirst()); + fChargeLastHiGain = fHChargeHiGain->GetBinLowEdge(fHChargeHiGain->GetXaxis()->GetLast()) + +fHChargeHiGain->GetBinWidth(0); + fChargeNbinsHiGain = nbins; + + CutEdges(fHChargeLoGain,nbins); + + fChargeFirstLoGain = fHChargeLoGain->GetBinLowEdge(fHChargeLoGain->GetXaxis()->GetFirst()); + fChargeLastLoGain = fHChargeLoGain->GetBinLowEdge(fHChargeLoGain->GetXaxis()->GetLast()) + +fHChargeLoGain->GetBinWidth(0); + fChargeNbinsLoGain = nbins; + + CutEdges(fHChargevsNHiGain,0); + CutEdges(fHChargevsNLoGain,0); } Index: trunk/MagicSoft/Mars/mhist/MHCalibrationPixel.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHCalibrationPixel.h (revision 2626) +++ trunk/MagicSoft/Mars/mhist/MHCalibrationPixel.h (revision 2627) @@ -31,7 +31,11 @@ protected: - TH1I* fHCharge; //-> Summed FADC slices - TH1I* fHTime; //-> Mean arrival time in number of FADC sice - TH1I* fHChargevsN; //-> Summed Charge vs. Event Nr. + TH1F* fHChargeHiGain; //-> Summed FADC slices + TH1I* fHTimeHiGain; //-> Mean arrival time in number of FADC sice + TH1I* fHChargevsNHiGain; //-> Summed Charge vs. Event Nr. + + TH1F* fHChargeLoGain; //-> Summed FADC slices + TH1I* fHTimeLoGain; //-> Mean arrival time in number of FADC sice + TH1I* fHChargevsNLoGain; //-> Summed Charge vs. Event Nr. TF1* fChargeGausFit; @@ -41,7 +45,11 @@ Axis_t fLowerFitRange; - Axis_t fChargeFirst; - Axis_t fChargeLast; - Int_t fChargeNbins; + Axis_t fChargeFirstHiGain; + Axis_t fChargeLastHiGain; + Int_t fChargeNbinsHiGain; + + Axis_t fChargeFirstLoGain; + Axis_t fChargeLastLoGain; + Int_t fChargeNbinsLoGain; Bool_t fFitOK; @@ -62,4 +70,6 @@ Double_t fTimeMean; Double_t fTimeSigma; + + Bool_t fUseLoGain; virtual void DrawLegend(); @@ -75,10 +85,16 @@ Bool_t Fill(const MParContainer *, const Stat_t w=1) { return kTRUE; } - Bool_t FillCharge(Int_t q) { return fHCharge->Fill(q) > -1; } - Bool_t FillTime(Int_t t) { return fHTime->Fill(t) > -1; } - Bool_t FillChargevsN(Float_t q, Int_t n) { return fHChargevsN->Fill(n,q) > -1; } + Bool_t FillChargeLoGain(Float_t q) { return fHChargeLoGain->Fill(q) > -1; } + Bool_t FillTimeLoGain(Int_t t) { return fHTimeLoGain->Fill(t) > -1; } + Bool_t FillChargevsNLoGain(Float_t q, Int_t n) { return fHChargevsNLoGain->Fill(n,q) > -1; } - const TH1I *GetHCharge() { return fHCharge; } - const TH1I *GetHCharge() const { return fHCharge; } + Bool_t FillChargeHiGain(Float_t q) { return fHChargeHiGain->Fill(q) > -1; } + Bool_t FillTimeHiGain(Int_t t) { return fHTimeHiGain->Fill(t) > -1; } + Bool_t FillChargevsNHiGain(Float_t q, Int_t n) { return fHChargevsNHiGain->Fill(n,q) > -1; } + + void SetUseLoGain() { fUseLoGain = kTRUE; } + + const TH1F *GetHCharge() { return fHChargeHiGain; } + const TH1F *GetHCharge() const { return fHChargeHiGain; } const Double_t GetChargeMean() const { return fChargeMean; } @@ -100,12 +116,15 @@ const Int_t GetTimeNdf() const { return fTimeNdf; } - const TH1I *GetHTime() { return fHTime; } - const TH1I *GetHTime() const { return fHTime; } + const TH1I *GetHTime() { return fHTimeHiGain; } + const TH1I *GetHTime() const { return fHTimeHiGain; } - const TH1I *GetHChargevsN() { return fHChargevsN; } - const TH1I *GetHChargevsN() const { return fHChargevsN; } + const TH1I *GetHChargevsN() { return fHChargevsNHiGain; } + const TH1I *GetHChargevsN() const { return fHChargevsNHiGain; } - Bool_t FitCharge(Option_t *option="RQ0"); - Bool_t FitTime(Axis_t rmin=0, Axis_t rmax=0, Option_t *option="RQ0"); + Bool_t FitChargeHiGain(Option_t *option="RQ0"); + Bool_t FitTimeHiGain(Axis_t rmin=0, Axis_t rmax=0, Option_t *option="RQ0"); + + Bool_t FitChargeLoGain(Option_t *option="RQ0"); + Bool_t FitTimeLoGain(Axis_t rmin=0, Axis_t rmax=0, Option_t *option="RQ0"); virtual void Draw(Option_t *option=""); Index: trunk/MagicSoft/Mars/mhist/MHCamera.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHCamera.cc (revision 2626) +++ trunk/MagicSoft/Mars/mhist/MHCamera.cc (revision 2627) @@ -739,4 +739,27 @@ // ------------------------------------------------------------------------ // +// Call this function to add a MCamEvent on top of the present contents. +// +void MHCamera::SetCamError(const MCamEvent &evt, Int_t type) +{ + + if (fNcells<=1 || IsFreezed()) + return; + + // FIXME: Security check missing! + for (Int_t idx=0; idx