Index: trunk/MagicSoft/Mars/macros/pedestalstudies.C =================================================================== --- trunk/MagicSoft/Mars/macros/pedestalstudies.C (revision 3316) +++ trunk/MagicSoft/Mars/macros/pedestalstudies.C (revision 3316) @@ -0,0 +1,319 @@ +/* ======================================================================== *\ +! +! * +! * 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, 11/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +//const TString pedfile = "/remote/home/pc2/operator/Crab20040214/20040215_16743_P_CrabOn_E.root"; +const TString pedfile = "../20040215_16770_P_OffCrab4_E.root"; + +//const TString pedfile = "/mnt/users/mdoro/Mars/Data/20040201_14418_P_OffMrk421-1_E.root"; +//const TString pedfile = "/mnt/Data/rootdata/CrabNebula/2004_02_10/20040210_14607_P_CrabNebula_E.root"; +//const TString pedfile = "/mnt/Data/rootdata/CrabNebula/2004_01_26/20040125_10412_P_Crab-On_E.root"; +//const TString pedfile = "/mnt/Data/rootdata/Miscellaneous/2003_12_19/20031218_03522_P_Park_E.root"; + +void pedestalstudies(TString pedname=pedfile) +{ + + gStyle->SetOptStat(1111); + gStyle->SetOptFit(); + + MStatusDisplay *display = new MStatusDisplay; + display->SetUpdateTime(500); + display->Resize(850,700); + + // + // 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); + + // + // Now setup the tasks and tasklist for the pedestals: + // --------------------------------------------------- + // + + MReadMarsFile read("Events", pedname); + read.DisableAutoScheme(); + + MGeomApply geomapl; + MExtractSignal sigcalc; + + // + // Set the extraction range higher: + // + //sigcalc.SetRange(1,14,1,14); + + MPedCalcPedRun pedcalc; + + // + // Additionally to calculating the pedestals, + // you can fill histograms and look at them + // + MFillH fill("MHPedestalCam", "MExtractedSignalCam"); + + tlist.AddToList(&read); + tlist.AddToList(&geomapl); + tlist.AddToList(&sigcalc); + tlist.AddToList(&pedcalc); + tlist.AddToList(&fill); + + MGeomCamMagic geomcam; + MPedestalCam pedcam; + MHPedestalCam hpedcam; + plist.AddToList(&geomcam); + plist.AddToList(&pedcam); + plist.AddToList(&hpedcam); + + // + // Create and setup the eventloop + // + MEvtLoop evtloop; + + evtloop.SetParList(&plist); + evtloop.SetDisplay(display); + + // + // Execute first analysis + // + if (!evtloop.Eventloop()) + return; + + tlist.PrintStatistics(); + + // + // Look at one specific pixel, after all the histogram manipulations: + // +// hpedcam[9].DrawClone("fourierevents"); + + + MHCamera dispped0 (geomcam, "Ped;Pedestal", "Mean per Slice"); + MHCamera dispped1 (geomcam, "Ped;PedestalErr", "Mean Error per Slice"); + MHCamera dispped2 (geomcam, "Ped;PedestalRms", "RMS per Slice"); + MHCamera dispped3 (geomcam, "Ped;PedestalRmsErr", "RMS Error per Slice"); + + MHCamera dispped4 (geomcam, "Ped;Mean", "Fitted Mean per Slice"); + MHCamera dispped5 (geomcam, "Ped;MeanErr", "Fitted Error of Mean per Slice"); + MHCamera dispped6 (geomcam, "Ped;Sigma", "Fitted Sigma per Slice"); + MHCamera dispped7 (geomcam, "Ped;SigmaErr", "Fitted Error of Sigma per Slice"); + MHCamera dispped8 (geomcam, "Ped;Prob", "Probability of Fit"); + MHCamera dispped9 (geomcam, "Ped;DeltaPedestalMean", "Rel. Diff. Mean per Slice (Calc.-Fitte)"); + MHCamera dispped10 (geomcam, "Ped;DeltaPedestalMeanError", "Rel. Diff. Mean Error per Slice (Calc.-Fitted)"); + MHCamera dispped11 (geomcam, "Ped;DeltaRmsSigma", "Rel. Diff. RMS per Slice (Calc.-Fitted)"); + MHCamera dispped12 (geomcam, "Ped;DeltaRmsSigmaError", "Rel. Diff. RMS Error per Slice (Calc.-Fitted)"); + MHCamera dispped13 (geomcam, "Ped;FitOK", "Gaus Fit not OK"); + MHCamera dispped14 (geomcam, "Ped;FourierOK", "Fourier Analysis not OK"); + + dispped0.SetCamContent( pedcam, 0); + dispped0.SetCamError( pedcam, 1); + dispped1.SetCamContent( pedcam, 1); + dispped2.SetCamContent( pedcam, 2); + dispped2.SetCamError( pedcam, 3); + dispped3.SetCamContent( pedcam, 3); + + dispped4.SetCamContent( hpedcam, 0); + dispped4.SetCamError( hpedcam, 1); + dispped5.SetCamContent( hpedcam, 1); + dispped6.SetCamContent( hpedcam, 2); + dispped6.SetCamError( hpedcam, 3); + dispped7.SetCamContent( hpedcam, 3); + dispped8.SetCamContent( hpedcam, 4); + dispped9.SetCamContent( hpedcam, 5); + dispped9.SetCamError( hpedcam, 6); + dispped10.SetCamContent(hpedcam, 7); + dispped11.SetCamContent(hpedcam, 8); + dispped11.SetCamError( hpedcam, 9); + dispped12.SetCamContent(hpedcam, 10); + dispped13.SetCamContent(hpedcam, 11); + dispped14.SetCamContent(hpedcam, 12); + + dispped0.SetYTitle("Calc. Pedestal per slice [FADC counts]"); + dispped1.SetYTitle("Calc. Pedestal Error per slice [FADC counts]"); + dispped2.SetYTitle("Calc. Pedestal RMS per slice [FADC counts]"); + dispped3.SetYTitle("Calc. Pedestal RMS Error per slice [FADC counts]"); + dispped4.SetYTitle("Fitted Mean per slice [FADC counts]"); + dispped5.SetYTitle("Error of Fitted Mean per slice [FADC counts]"); + dispped6.SetYTitle("Fitted Sigma per slice [FADC counts]"); + dispped7.SetYTitle("Error of Fitted Sigma per slice [FADC counts]"); + dispped8.SetYTitle("Fit Probability [1]"); + dispped9.SetYTitle("Rel. Diff. Pedestal Calc.-Fitted per slice [1]"); + dispped10.SetYTitle("Rel. Diff. Pedestal Error Calc.-Fitted per slice [1]"); + dispped11.SetYTitle("Rel. Diff. Pedestal RMS Calc.-Fitted per slice [1]"); + dispped12.SetYTitle("Rel. Diff. Pedestal RMS Error Calc.-Fitted per slice [1]"); + dispped13.SetYTitle("[1]"); + dispped14.SetYTitle("[1]"); + + // Histogram values + TCanvas &b1 = display->AddTab("Ped.Calc."); + b1.Divide(4,3); + + CamDraw(b1,dispped0,pedcam,1,4,1); + CamDraw(b1,dispped1,pedcam,2,4,2); + CamDraw(b1,dispped2,pedcam,3,4,2); + CamDraw(b1,dispped3,pedcam,4,4,2); + + // Fitted values + TCanvas &b2 = display->AddTab("Ped.Fit"); + b2.Divide(4,3); + + CamDraw(b2,dispped4,hpedcam,1,4,1); + CamDraw(b2,dispped5,hpedcam,2,4,2); + CamDraw(b2,dispped6,hpedcam,3,4,2); + CamDraw(b2,dispped7,hpedcam,4,4,2); + + + // Fits Probability + TCanvas &b3 = display->AddTab("Ped.Fit Prob."); + b3.Divide(1,3); + + CamDraw(b3,dispped8,hpedcam,1,1,3); + + // Differences + TCanvas &c4 = display->AddTab("Rel.Diff.Calc.-Fit"); + c4.Divide(4,3); + + CamDraw(c4,dispped9,hpedcam,1,4,1); + CamDraw(c4,dispped10,hpedcam,2,4,1); + CamDraw(c4,dispped11,hpedcam,3,4,1); + CamDraw(c4,dispped12,hpedcam,4,4,1); + + // Defects + TCanvas &c5 = display->AddTab("Defects"); + c5.Divide(2,2); + + CamDraw(c5,dispped13,hpedcam,1,2,0); + CamDraw(c5,dispped14,hpedcam,2,2,0); + +} + +void CamDraw(TCanvas &c, MHCamera &cam, MCamEvent &evt, Int_t i, Int_t j, Int_t fit) +{ + + c.cd(i); + gPad->SetBorderMode(0); + MHCamera *obj1=(MHCamera*)cam.DrawCopy("hist"); + // obj1->AddNotify(evt); + + c.cd(i+j); + gPad->SetBorderMode(0); + obj1->Draw(); + ((MHCamera*)obj1)->SetPrettyPalette(); + + if (fit != 0) + { + c.cd(i+2*j); + gPad->SetBorderMode(0); + TH1D *obj2 = (TH1D*)obj1->Projection(); + +// obj2->Sumw2(); + obj2->Draw(); + obj2->SetBit(kCanDelete); + + const Double_t min = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst()); + const Double_t max = obj2->GetBinCenter(obj2->GetXaxis()->GetLast()); + const Double_t integ = obj2->Integral("width")/2.5066283; + const Double_t mean = obj2->GetMean(); + const Double_t rms = obj2->GetRMS(); + const Double_t width = max-min; + + if (rms == 0. || width == 0. ) + return; + + switch (fit) + { + case 1: + TF1 *sgaus = new TF1("sgaus","gaus(0)",min,max); + sgaus->SetBit(kCanDelete); + sgaus->SetParNames("Area","#mu","#sigma"); + sgaus->SetParameters(integ/rms,mean,rms); + sgaus->SetParLimits(0,0.,integ); + sgaus->SetParLimits(1,min,max); + sgaus->SetParLimits(2,0,width/1.5); + obj2->Fit("sgaus","QLR"); + obj2->GetFunction("sgaus")->SetLineColor(kYellow); + break; + + case 2: + TString dgausform = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"; + dgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])"; + TF1 *dgaus = new TF1("dgaus",dgausform.Data(),min,max); + dgaus->SetBit(kCanDelete); + dgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}","A_{2}","#mu_{2}","#sigma_{2}"); + dgaus->SetParameters(integ,(min+mean)/2.,width/4., + integ/width/2.,(max+mean)/2.,width/4.); + // The left-sided Gauss + dgaus->SetParLimits(0,integ-1.5,integ+1.5); + dgaus->SetParLimits(1,min+(width/10.),mean); + dgaus->SetParLimits(2,0,width/2.); + // The right-sided Gauss + dgaus->SetParLimits(3,0,integ); + dgaus->SetParLimits(4,mean,max-(width/10.)); + dgaus->SetParLimits(5,0,width/2.); + obj2->Fit("dgaus","QLRM"); + obj2->GetFunction("dgaus")->SetLineColor(kYellow); + break; + + case 3: + TString tgausform = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"; + tgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])"; + tgausform += "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])"; + TF1 *tgaus = new TF1("tgaus",tgausform.Data(),min,max); + tgaus->SetBit(kCanDelete); + tgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}", + "A_{2}","#mu_{2}","#sigma_{2}", + "A_{3}","#mu_{3}","#sigma_{3}"); + tgaus->SetParameters(integ,(min+mean)/2,width/4., + integ/width/3.,(max+mean)/2.,width/4., + integ/width/3.,mean,width/2.); + // The left-sided Gauss + tgaus->SetParLimits(0,integ-1.5,integ+1.5); + tgaus->SetParLimits(1,min+(width/10.),mean); + tgaus->SetParLimits(2,width/15.,width/2.); + // The right-sided Gauss + tgaus->SetParLimits(3,0.,integ); + tgaus->SetParLimits(4,mean,max-(width/10.)); + tgaus->SetParLimits(5,width/15.,width/2.); + // The Gauss describing the outliers + tgaus->SetParLimits(6,0.,integ); + tgaus->SetParLimits(7,min,max); + tgaus->SetParLimits(8,width/4.,width/1.5); + obj2->Fit("tgaus","QLRM"); + obj2->GetFunction("tgaus")->SetLineColor(kYellow); + break; + case 4: + obj2->Fit("pol0","Q"); + obj2->GetFunction("pol0")->SetLineColor(kYellow); + break; + case 9: + break; + default: + obj2->Fit("gaus","Q"); + obj2->GetFunction("gaus")->SetLineColor(kYellow); + break; + } + + gPad->Modified(); + gPad->Update(); + + } +}