Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 2820) +++ trunk/MagicSoft/Mars/Changelog (revision 2821) @@ -4,4 +4,16 @@ -*-*- END OF LINE -*-*- + 2004/01/15: Javier Rico + + * manalysis/MPedCalcPedRun.[h,cc] + - optimize the running time + - add (some) documentation + - correct treatment for the case of several input files + + * macros/pedvsevent.C + - added + - draw pedestal mean and rms vs event# for input pixel# and run + file, and compares them to the global pedestal mean and rms + 2004/01/15: Raquel de los Reyes Index: trunk/MagicSoft/Mars/macros/pedvsevent.C =================================================================== --- trunk/MagicSoft/Mars/macros/pedvsevent.C (revision 2821) +++ trunk/MagicSoft/Mars/macros/pedvsevent.C (revision 2821) @@ -0,0 +1,214 @@ + +///////////////////////////////////////////////////////////////////////////// +// // +// Author(s): S.C. Commichau, Javier Rico, 12/2003 // +// // +// Macro to generate pedestal vs time (event) plot for a single pixel // +// // +///////////////////////////////////////////////////////////////////////////// + + +const Int_t default_pixel = 1 ; + +void pedvsevent(Int_t pixel = default_pixel, + TString pname = "/disc02/Data/rootdata/Miscellaneous/2003_11_29/20031128_03118_P_Park-camera-closed_E.root") +{ + + TH1F* Pedestal = new TH1F("Pedestal","Pedestals",100,0,20); + TH1F* PedestalRMS = new TH1F("PedestalRMS","Pedestal RMS",100,0,10); + + //Create an empty parameter list and an empty task list + //the tasklist is identified in the eventloop by ist name + MParList plist; + MTaskList tlist; + + //Creates a MPedestalPix object for each pixel, i.e. it is a + //storage container for all Pedestal information in the camera + MPedestalCam cam; + + plist.AddToList(&cam); + + //MHCamEvent hist; + //hist.SetType(1); + //plist.AddToList(&hist); + + plist.AddToList(&tlist); + + //Setup task and tasklist for the pedestals + MReadMarsFile read("Events", pname); + read.DisableAutoScheme(); + + //Apply the geometry to geometry dependant containers. + //MGeomApply changes the size of the arrays in the containers to a size + //matching the number of pixels, eg: MPedestalCam, MBlindPixels + //Default geometry is MGeomCamMagic + MGeomApply geomapl; + + //Task to calculate pedestals + MPedCalcPedRun ped; + + tlist.AddToList(&read); + tlist.AddToList(&geomapl); + tlist.AddToList(&ped); + + //ped.SetPixel(pixel); + //ped.Draw(); + + //Create and setup the 1st Eventloop + MEvtLoop evtloop; + evtloop.SetParList(&plist); + + //Execute first analysis, pedestals... + if (!tlist.PreProcess(&plist)) + return; + + const Int_t nevents = read.GetEntries(); + + Float_t x[nevents], rms[nevents], mean[nevents]; + + + Int_t i = 0; + while (tlist.Process()) + { + mean[i] = cam[pixel].GetPedestal(); + rms[i] = cam[pixel].GetPedestalRms(); + x[i] = i; + i++; + } + + TGraphErrors* pedgraph = new TGraphErrors(nevents,x,mean,NULL,NULL); + TGraphErrors* rmsgraph = new TGraphErrors(nevents,x,rms,NULL,NULL); + + //plist.FindObject("MPedestalCam")->Print(); + + tlist.PostProcess(); + //if (!evtloop.Eventloop()) + // return; + + Float_t finalmean = cam[pixel].GetPedestal(); + Float_t finalrms = cam[pixel].GetPedestalRms(); + + TLine* pedline = new TLine(0,finalmean,nevents,finalmean); + TLine* rmsline = new TLine(0,finalrms,nevents,finalrms); + + //The draw area + + + + gROOT->Reset(); + gStyle->SetOptStat(0); + // Set statistics options, 1111111 + // ||||||| + // ||||||histogram name + // |||||number of entries + // ||||mean value + // |||RMS + // ||underflows + // |overflows + // sum of bins + + // Set gPad options + gStyle->SetFrameBorderMode(0); + gStyle->SetPalette(1); + // Delete Frames of subCanvas' does not work, hook it Mr. gPad! + gStyle->SetFrameBorderSize(0); + gStyle->SetCanvasColor(0); + gStyle->SetFrameFillColor(0); + gStyle->SetTitleFont(102); + gStyle->SetTitleFillColor(0); + gStyle->SetTitleBorderSize(0); + gStyle->SetStatColor(0); + gStyle->SetStatBorderSize(0); + + // Set Canvas options + TCanvas *MyC1 = new TCanvas("MyC","Pedestal vs Event", 0, 100, 900, 500); + MyC->SetBorderMode(0); // Delete the Canvas' border line + + MyC->cd(); + gPad->SetBorderMode(0); + + // TLine* pedline = new TLine(0,finalmean,nevents,finalmean); + //TLine* rmsline = new TLine(0,finalrms,nevents,finalrms); + + tlist.PrintStatistics(); + //plist.FindObject("MPedestalCam")->Print(); + Size_t pos = pname.Last('/')+10; + TString iRun = TString(pname(pos,5)); + + char str[64]; + + sprintf(str,"Run %s Pixel %d",iRun.Data(),pixel); + + pedgraph->SetMaximum(30); + pedgraph->SetMinimum(0); + pedgraph->SetMarkerStyle(24); + pedgraph->SetMarkerSize(.5); + pedgraph->GetXaxis()->SetTitleFont(102); + pedgraph->GetYaxis()->SetTitleFont(102); + pedgraph->GetXaxis()->SetLabelFont(102); + pedgraph->GetYaxis()->SetLabelFont(102); + pedgraph->SetTitle(str); +// pedgraph->SetTitleFont(102); + pedgraph->GetYaxis()->SetTitleFont(102); + pedgraph->GetXaxis()->SetTitle("Event"); + pedgraph->GetYaxis()->SetTitle("[FADC Counts]"); + pedgraph->GetXaxis()->SetLimits(0,nevents-1); + + rmsgraph->SetMarkerStyle(25); + rmsgraph->SetMarkerSize(.5); + rmsgraph->SetMarkerColor(8); + + pedline->SetLineColor(2); + rmsline->SetLineColor(4); + pedline->SetLineWidth(2); + rmsline->SetLineWidth(2); + + pedgraph->Draw("AP"); + rmsgraph->Draw("P"); + + pedline->Draw("same"); + rmsline->Draw("same"); + + TLegend* leg = new TLegend(.14,.68,.39,.88); + leg->SetHeader(""); + leg->AddEntry(pedgraph,"Event based Pedestal","P"); + leg->AddEntry(pedline,"Run based Pedestal","L"); + leg->AddEntry(rmsgraph,"Event based RMS","P"); + leg->AddEntry(rmsline,"Run based RMS","L"); + leg->SetFillColor(0); + leg->SetLineColor(1); + leg->SetBorderSize(1); + + + leg->Draw("same"); +} + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Index: trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.cc (revision 2820) +++ trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.cc (revision 2821) @@ -18,6 +18,8 @@ ! Author(s): Josep Flix 04/2001 ! Author(s): Thomas Bretz 05/2001 +! Author(s): Sebastian Commichau 12/2003 +! Author(s): Javier Rico 01/2004 ! -! Copyright: MAGIC Software Development, 2000-2001 +! Copyright: MAGIC Software Development, 2000-2004 ! ! @@ -28,4 +30,10 @@ // MPedCalcPedRun // // // +// This task takes a pedestal run file and fills MPedestalCam during // +// the Process() with the pedestal and rms computed in an event basis. // +// In the PostProcess() MPedestalCam is finally filled with the pedestal // +// mean and rms computed in a run basis. // +// More than one run (file) can be merged // +// // // Input Containers: // // MRawEvtData // @@ -33,4 +41,5 @@ // Output Containers: // // MPedestalCam // +// // // // ///////////////////////////////////////////////////////////////////////////// @@ -70,5 +79,5 @@ if (!fRawEvt) { - *fLog << dbginf << "MRawEvtData not found... aborting." << endl; + *fLog << err << "MRawEvtData not found... aborting." << endl; return kFALSE; } @@ -78,69 +87,75 @@ return kFALSE; - MGeomCamMagic magiccam; - - fSumx.Set(magiccam.GetNumPixels()); - fSumx2.Set(magiccam.GetNumPixels()); - - for(UInt_t i=0;iFindObject("MRawRunHeader"); + if (!runheader) + { + *fLog << warn << dbginf; + *fLog << "Warning - cannot check file type, MRawRunHeader not found." << endl; + } + else + if (runheader->GetRunType() == kRTMonteCarlo) + return kTRUE; - fRunheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); - if (!fRunheader) - { - *fLog << warn << dbginf << - "Warning - cannot check file type, MRawRunHeader not found." << endl; - } - else - if (fRunheader->GetRunType() == kRTMonteCarlo) - { - return kTRUE; - } + fNumPixels = fPedestals->GetSize(); - fNumHiGainSamples = fRunheader->GetNumSamplesHiGain(); + if(fSumx.GetSize()==0) + { + fSumx.Set(fNumPixels); + fSumx2.Set(fNumPixels); - fPedestals->InitSize(fRunheader->GetNumPixel()); + memset(fSumx.GetArray(), 0, sizeof(Float_t)*fNumPixels); + memset(fSumx2.GetArray(), 0, sizeof(Float_t)*fNumPixels); + } + + // Calculate an even number for the hi gain samples to avoid + // biases due to the fluctuation in pedestal from one slice to + // the other one + fNumHiGainSamples = runheader->GetNumSamplesHiGain() & ~1; return kTRUE; - } - Int_t MPedCalcPedRun::Process() { - MRawEvtPixelIter pixel(fRawEvt); while (pixel.Next()) { - Byte_t shift=(fNumHiGainSamples/2*2==fNumHiGainSamples) ? 0:1; Byte_t *ptr = pixel.GetHiGainSamples(); - const Byte_t *end = ptr + fRawEvt->GetNumHiGainSamples()-shift; + const Byte_t *end = ptr + fNumHiGainSamples; - const Float_t higainped = CalcHiGainMean(ptr, end); - const Float_t higainrms = CalcHiGainRms(ptr, end, higainped); - - const UInt_t pixid = pixel.GetPixelId(); - MPedestalPix &pix = (*fPedestals)[pixid]; + UInt_t sum = 0; + UInt_t sqr = 0; - // cumulate the sum of pedestals and of pedestal squares - fSumx.AddAt(higainped+fSumx.At(pixid),pixid); - fSumx2.AddAt(GetSumx2(ptr, end)+fSumx2.At(pixid),pixid); - - // set the value of the pedestal and rms computed from the processed event - pix.Set(higainped, higainrms); + do + { + sum += *ptr; + sqr += *ptr * *ptr; + } + while (++ptr != end); + + const Float_t msum = (Float_t)sum; + const Float_t msqr = (Float_t)sqr; + + const Float_t higainped = msum/fNumHiGainSamples; + const Float_t higainrms = TMath::Sqrt((msqr-msum*msum/fNumHiGainSamples)/(fNumHiGainSamples-1.)); + + const UInt_t idx = pixel.GetPixelId(); + (*fPedestals)[idx].Set(higainped, higainrms); + + fSumx[idx] += msum; + fSumx2[idx] += sqr; } - + fPedestals->SetReadyToSave(); - + fNumSamplesTot+=fNumHiGainSamples; + return kTRUE; } @@ -148,84 +163,24 @@ Int_t MPedCalcPedRun::PostProcess() { -// Compute pedestals and rms from the whole run - - MRawEvtPixelIter pixel(fRawEvt); - - while (pixel.Next()) + // Compute pedestals and rms from the whole run + const Int_t n = fNumSamplesTot; + + MRawEvtPixelIter pixel(fRawEvt); + + while (pixel.Next()) { - const UInt_t pixid = pixel.GetPixelId(); - MPedestalPix &pix = (*fPedestals)[pixid]; - - const Int_t N = GetNumExecutions(); - const Float_t sum = fSumx.At(pixid); - const Float_t sum2 = fSumx2.At(pixid); - const Float_t higainped = sum/N; - const Float_t higainrms = sqrt(1./(N-1.)*(sum2-sum*sum/N)); - pix.Set(higainped,higainrms); - + const UInt_t pixid = pixel.GetPixelId(); + + const Float_t sum = fSumx.At(pixid); + const Float_t sum2 = fSumx2.At(pixid); + + const Float_t higainped = sum/n; + const Float_t higainrms = TMath::Sqrt((sum2-sum*sum/n)/(n-1.)); + + (*fPedestals)[pixid].Set(higainped, higainrms); } - return kTRUE; - + + return kTRUE; } -Float_t MPedCalcPedRun::CalcHiGainMean(Byte_t *ptr, const Byte_t *end) const -{ - Int_t sum=0; - Byte_t EvenNumSamples=(fNumHiGainSamples/2*2==fNumHiGainSamples) - ? fNumHiGainSamples : fNumHiGainSamples-1; - - do sum += *ptr; - while (++ptr != end); - - return (Float_t)sum/EvenNumSamples; -} - - - -Float_t MPedCalcPedRun::GetSumx2(Byte_t *ptr, const Byte_t *end) const -{ - Float_t square = 0; - - // Take an even number of time slices to avoid biases due to A/B effect - Byte_t EvenNumSamples=(fNumHiGainSamples/2*2==fNumHiGainSamples) ? fNumHiGainSamples:fNumHiGainSamples-1; - - do - { - const Float_t val = (Float_t)(*ptr); - - square += val*val; - } while (++ptr != end); - - return square/EvenNumSamples; -} - - -Float_t MPedCalcPedRun::CalcHiGainRms(Byte_t *ptr, const Byte_t *end, Float_t higainped) const -{ - Float_t rms = 0; - Byte_t EvenNumSamples=(fNumHiGainSamples/2*2==fNumHiGainSamples) ? fNumHiGainSamples:fNumHiGainSamples-1; - - do - { - const Float_t diff = (Float_t)(*ptr)-higainped; - - rms += diff*diff; - } while (++ptr != end); - - return sqrt(rms/(EvenNumSamples-1)); -} - - - -/* -Float_t MPedCalcPedRun::CalcHiGainMeanErr(Float_t higainrms) const -{ - return higainrms/sqrt((float)fNumHiGainSamples); -} - -Float_t MPedCalcPedRun::CalcHiGainRmsErr(Float_t higainrms) const -{ - return higainrms/sqrt(2.*fNumHiGainSamples); -} -*/ Index: trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.h (revision 2820) +++ trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.h (revision 2821) @@ -16,5 +16,4 @@ #include -class MRawRunHeader; class MRawEvtData; class MPedestalCam; @@ -22,7 +21,8 @@ class MPedCalcPedRun : public MTask { - Byte_t fNumHiGainSamples; + Byte_t fNumHiGainSamples; + UShort_t fNumPixels; + ULong_t fNumSamplesTot; - MRawRunHeader *fRunheader; // raw event run header MRawEvtData *fRawEvt; // raw event data (time slices) MPedestalCam *fPedestals; // Pedestals of all pixels in the camera @@ -30,10 +30,4 @@ TArrayF fSumx; // sum of values TArrayF fSumx2; // sum of squared values - - Float_t CalcHiGainMean(Byte_t *ptr, const Byte_t *end) const; - Float_t CalcHiGainRms(Byte_t *ptr, const Byte_t *end, Float_t higainped) const; - Float_t GetSumx2(Byte_t* ptr, const Byte_t* end) const; - //Float_t CalcHiGainMeanErr(Float_t higainrms) const; - //Float_t CalcHiGainRmsErr(Float_t higainrms) const; Bool_t ReInit(MParList *pList); @@ -44,5 +38,4 @@ public: - MPedCalcPedRun(const char *name=NULL, const char *title=NULL);