Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 3759) +++ trunk/MagicSoft/Mars/Changelog (revision 3760) @@ -26,4 +26,5 @@ * macros/bootcampstandardanalysis.C + * macros/pedphotcalc.C - updated the macro to the cvs-standard used at the Udine bootcamp Index: trunk/MagicSoft/Mars/macros/pedphotcalc.C =================================================================== --- trunk/MagicSoft/Mars/macros/pedphotcalc.C (revision 3759) +++ trunk/MagicSoft/Mars/macros/pedphotcalc.C (revision 3760) @@ -26,31 +26,53 @@ ! \* ======================================================================== */ - +///////////////////////////////////////////////////////////////////////////// +// +// pedphotcalc.C +// +// This macro is an example of the use of MPedPhotCalc. It computes +// the pedestal mean and rms from pedestal files undergoing the +// signal extraction + calibration chain, in units of photons. It +// produces plots of the relevant computed quantities. +// +// Needs as arguments the run number of a pedestal file ("*_P_*.root"), +// one of a calibration file ("*_C_*.root"). +// performs the pedestal calculation, the calibration +/// constants calculation and the calibration of the pedestals. +// +// The TString inpath has to be set correctly. +// +// The macro searches for the pulser colour which corresponds to the calibration +// run number. If the run number is smaller than 20000, pulser colour "CT1" +// is assumed, otherwise, it searches for the strings "green", "blue", "uv" or +// "ct1" in the filenames. If no colour or multiple colours are found, the +// execution is aborted. +// +//////////////////////////////////////////////////////////////////////////////////// #include "MParList.h" #include "MTaskList.h" +#include "MJPedestal.h" +#include "MJCalibration.h" #include "MPedestalCam.h" +#include "MPedestalPix.h" #include "MReadMarsFile.h" #include "MGeomApply.h" -#include "MPedCalcPedRun.h" #include "MGeomCamMagic.h" #include "MEvtLoop.h" +#include "MCalibrationCam.h" #include "MCalibrationChargeCam.h" -#include "MHCalibrationChargeCam.h" -#include "MCalibrationChargePINDiode.h" -#include "MHCalibrationChargePINDiode.h" -#include "MCalibrationChargeBlindPix.h" -#include "MHCalibrationChargeBlindPix.h" -#include "MCalibrationChargeCalc.h" +#include "MCalibrationChargePix.h" +#include "MCalibrationQECam.h" +#include "MCalibrationQEPix.h" #include "MExtractedSignalCam.h" #include "MExtractSignal.h" -#include "MExtractedSignalPINDiode.h" -#include "MExtractPINDiode.h" -#include "MExtractedSignalBlindPixel.h" -#include "MExtractBlindPixel.h" #include "MCerPhotEvt.h" #include "MCalibrate.h" #include "MPedPhotCalc.h" #include "MPedPhotCam.h" -#include "MBadPixelsCam.h" +#include "MPedPhotPix.h" +#include "MHCamera.h" +#include "MRunIter.h" +#include "MDirIter.h" +#include "MStatusDisplay.h" #include "MHCamera.h" @@ -58,56 +80,279 @@ #include "TString.h" #include "TCanvas.h" - -#include - - -const TString pedfile="/disc02/Data/rootdata/Miscellaneous/2003_12_19/20031218_03522_P_Park_E.root"; -const TString calfile="/disc02/Data/rootdata/Miscellaneous/2003_12_19/20031218_03527_C_Park_E.root"; - -void pedphotcalc(TString pedname=pedfile,TString calname=calfile) +#include "TStyle.h" +#include "TF1.h" +#include "TLegend.h" + +#include +#include "Getline.h" + +const TString inpath = "/mnt/Data/rootdata/CrabNebula/2004_02_10/"; +const Int_t dpedrun = 14607; +const Int_t dcalrun1 = 14608; +const Int_t dcalrun2 = 0; +const Bool_t usedisplay = kTRUE; + +static MCalibrationCam::PulserColor_t FindColor(MDirIter* run) { - /* - This macro is an example of the use of MPedPhotCalc. It computes - the pedestal mean and rms from pedestal files undergoing the - signal extraction + calibration chain, in units of photons. It - produces plots of the relevant computed quantities. - */ - + MCalibrationCam::PulserColor_t col = MCalibrationCam::kNONE; + + TString filenames; + + while (!(filenames=run->Next()).IsNull()) + { + + filenames.ToLower(); + + if (filenames.Contains("green")) + if (col == MCalibrationCam::kNONE) + { + cout << "Found colour: Green in " << filenames << endl; + col = MCalibrationCam::kGREEN; + } + else if (col != MCalibrationCam::kGREEN) + { + cout << "Different colour found in " << filenames << "... abort" << endl; + return MCalibrationCam::kNONE; + } + + if (filenames.Contains("blue")) + if (col == MCalibrationCam::kNONE) + { + cout << "Found colour: Blue in " << filenames << endl; + col = MCalibrationCam::kBLUE; + } + else if (col != MCalibrationCam::kBLUE) + { + cout << "Different colour found in " << filenames << "... abort" << endl; + return MCalibrationCam::kNONE; + } + + if (filenames.Contains("uv")) + if (col == MCalibrationCam::kNONE) + { + cout << "Found colour: Uv in " << filenames << endl; + col = MCalibrationCam::kUV; + } + else if (col != MCalibrationCam::kUV) + { + cout << "Different colour found in " << filenames << "... abort" << endl; + return MCalibrationCam::kNONE; + } + + if (filenames.Contains("ct1")) + if (col == MCalibrationCam::kNONE) + { + cout << "Found colour: Ct1 in " << filenames << endl; + col = MCalibrationCam::kCT1; + } + else if (col != MCalibrationCam::kCT1) + { + cout << "Different colour found in " << filenames << "... abort" << endl; + return MCalibrationCam::kNONE; + } + + } + + + + if (col == MCalibrationCam::kNONE) + cout << "No colour found in filenames of runs: " << ((MRunIter*)run)->GetRunsAsString() + << "... abort" << endl; + + return col; +} + + + +void DrawProjection(MHCamera *obj1, Int_t fit) +{ + TH1D *obj2 = (TH1D*)obj1->Projection(obj1->GetName()); + obj2->SetDirectory(0); + obj2->Draw(); + obj2->SetBit(kCanDelete); + gPad->SetLogy(); + + if (obj1->GetGeomCam().InheritsFrom("MGeomCamMagic")) + { + TArrayI s0(3); + s0[0] = 6; + s0[1] = 1; + s0[2] = 2; + + TArrayI s1(3); + s1[0] = 3; + s1[1] = 4; + s1[2] = 5; + + TArrayI inner(1); + inner[0] = 0; + + TArrayI outer(1); + outer[0] = 1; + + // Just to get the right (maximum) binning + TH1D *half[4]; + half[0] = obj1->ProjectionS(s0, inner, "Sector 6-1-2 Inner"); + half[1] = obj1->ProjectionS(s1, inner, "Sector 3-4-5 Inner"); + half[2] = obj1->ProjectionS(s0, outer, "Sector 6-1-2 Outer"); + half[3] = obj1->ProjectionS(s1, outer, "Sector 3-4-5 Outer"); + + for (int i=0; i<4; i++) + { + half[i]->SetLineColor(kRed+i); + half[i]->SetDirectory(0); + half[i]->SetBit(kCanDelete); + half[i]->Draw("same"); + } + } + + 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; + + const TString dgausformula("([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])" + "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])"); + + TF1 *f=0; + switch (fit) + { + // FIXME: MAYBE add function to TH1? + case 0: + f = new TF1("sgaus", "gaus(0)", min, max); + f->SetLineColor(kYellow); + f->SetBit(kCanDelete); + f->SetParNames("Area", "#mu", "#sigma"); + f->SetParameters(integ/rms, mean, rms); + f->SetParLimits(0, 0, integ); + f->SetParLimits(1, min, max); + f->SetParLimits(2, 0, width/1.5); + obj2->Fit(f, "QLRM"); + break; + + case 1: + f = new TF1("dgaus", dgausformula, min, max); + f->SetLineColor(kYellow); + f->SetBit(kCanDelete); + f->SetParNames("A_{tot}", "#mu_{1}", "#sigma_{1}", "A_{2}", "#mu_{2}", "#sigma_{2}"); + f->SetParameters(integ, (min+mean)/2, width/4, + integ/width/2, (max+mean)/2, width/4); + // The left-sided Gauss + f->SetParLimits(0, integ-1.5, integ+1.5); + f->SetParLimits(1, min+(width/10), mean); + f->SetParLimits(2, 0, width/2); + // The right-sided Gauss + f->SetParLimits(3, 0, integ); + f->SetParLimits(4, mean, max-(width/10)); + f->SetParLimits(5, 0, width/2); + obj2->Fit(f, "QLRM"); + break; + + default: + obj2->Fit("gaus", "Q"); + obj2->GetFunction("gaus")->SetLineColor(kYellow); + break; + } +} + +void CamDraw(TCanvas &c, Int_t x, Int_t y, MHCamera &cam1, Int_t fit) +{ + c.cd(x); + gPad->SetBorderMode(0); + MHCamera *obj1=(MHCamera*)cam1.DrawCopy("hist"); + + c.cd(x+y); + gPad->SetBorderMode(0); + obj1->Draw(); + + c.cd(x+2*y); + gPad->SetBorderMode(0); + DrawProjection(obj1, fit); +} + +void pedphotcalc(const Int_t prun=dpedrun, // pedestal file + const Int_t crun1=dcalrun1, const Int_t crun2=dcalrun2 // calibration file(s) + ) +{ + + MRunIter pruns; + MRunIter cruns; + + pruns.AddRun(prun,inpath); + + if (crun2==0) + cruns.AddRun(crun1,inpath); + else + cruns.AddRuns(crun1,crun2,inpath); + + MCalibrationCam::PulserColor_t color; + + if (crun1 < 20000) + color = MCalibrationCam::kCT1; + else + color = FindColor((MDirIter*)&cruns); + + if (color == MCalibrationCam::kNONE) + { + TTimer timer("gSystem->ProcessEvents();", 50, kFALSE); + + while (1) + { + timer.TurnOn(); + TString input = Getline("Could not find the correct colour: Type 'q' to exit, " + "green, blue, uv or ct1 to go on: "); + timer.TurnOff(); + + if (input=="q\n") + return ; + + if (input=="green") + color = MCalibrationCam::kGREEN; + if (input=="blue") + color = MCalibrationCam::kBLUE; + if (input=="uv") + color = MCalibrationCam::kUV; + if (input=="ct1") + color = MCalibrationCam::kCT1; + } + } + + // + // Now setup the tasks and tasklist for the pedestals: + // --------------------------------------------------- + // + MGeomCamMagic geomcam; + MGeomApply geomapl; + MStatusDisplay *display = NULL; + /************************************/ /* FIRST LOOP: PEDESTAL COMPUTATION */ /************************************/ - MParList plist; - MTaskList tlist; - plist.AddToList(&tlist); - - // containers - MPedestalCam pedcam; - plist.AddToList(&pedcam); - - //tasks - MReadMarsFile read("Events", pedname); - MGeomApply geomapl; - MPedCalcPedRun pedcalc; - MGeomCamMagic geomcam; - - read.DisableAutoScheme(); - - tlist.AddToList(&read); - tlist.AddToList(&geomapl); - tlist.AddToList(&pedcalc); - - // Create and execute the event looper - MEvtLoop evtloop; - evtloop.SetParList(&plist); + MJPedestal pedloop; + pedloop.SetInput(&pruns); + if (usedisplay) + { + display = new MStatusDisplay; + display->SetUpdateTime(3000); + display->Resize(850,700); + display->SetBit(kCanDelete); + pedloop.SetDisplay(display); + } cout << "*************************" << endl; cout << "** COMPUTING PEDESTALS **" << endl; cout << "*************************" << endl; - if (!evtloop.Eventloop()) + + if (!pedloop.Process()) return; - - tlist.PrintStatistics(); + + MPedestalCam pedcam = pedloop.GetPedestalCam(); /****************************/ @@ -115,76 +360,34 @@ /****************************/ - MParList plist2; - MTaskList tlist2; - plist2.AddToList(&tlist2); - - // containers - MCalibrationChargeCam calcam; - MExtractedSignalCam sigcam; - MBadPixelsCam badcam; - - MCalibrationChargeCam chargecam; - MCalibrationChargePINDiode chargepin; - MCalibrationChargeBlindPixel chargeblind; - chargeblind.SetColor(MCalibrationChargeBlindPix::kECT1); - chargepin.SetColor( MCalibrationChargePINDiode::kECT1); - // - plist2.AddToList(&geomcam); - plist2.AddToList(&pedcam); - plist2.AddToList(&calcam); - plist2.AddToList(&sigcam); - plist2.AddToList(&badcam); - plist2.AddToList(&chargecam); - plist2.AddToList(&chargepin); - plist2.AddToList(&chargeblind); - - //tasks - MReadMarsFile read2("Events", calname); - read2.DisableAutoScheme(); - - MExtractPINDiode pincalc; - MExtractBlindPixel blindcalc; - MExtractSignal sigcalc; - MCalibrationChargeCalc calcalc; - - MFillH fillpin("MHCalibrationChargePINDiode", "MExtractedSignalPINDiode"); - MFillH fillblind("MHCalibrationChargeBlindPixel", "MExtractedSignalBlindPixel"); - MFillH fillcam("MHCalibrationChargeCam" , "MExtractedSignalCam"); - // - // Apply a filter against cosmics - // (was directly in MCalibrationChargeCalc in earlier versions) - // - MFCosmics cosmics; - MContinue cont(&cosmics); - - tlist2.AddToList(&read2); - tlist2.AddToList(&geomapl); - tlist2.AddToList(&blindcalc); - tlist2.AddToList(&pincalc); - tlist2.AddToList(&sigcalc); - // - // In case, you want to skip the cosmics rejection, - // uncomment the next line - // - tlist2.AddToList(&cont); - // - tlist2.AddToList(&fillpin); - tlist2.AddToList(&fillblind); - tlist2.AddToList(&fillcam); - tlist2.AddToList(&calcalc); - - // Execute second loop (calibration) - MEvtLoop evtloop2; - evtloop2.SetParList(&plist2); + // + // Now setup the new tasks for the calibration: + // --------------------------------------------------- + // + MJCalibration calloop; + calloop.SetColor(color); + calloop.SetInput(&cruns); + // + // Use as signal extractor MExtractSignal: + // + calloop.SetExtractorLevel(1); + // + // The next two commands are for the display: + // + if (usedisplay) + { + calloop.SetDisplay(display); + calloop.SetDataCheck(); + } cout << "***********************************" << endl; cout << "** COMPUTING CALIBRATION FACTORS **" << endl; cout << "***********************************" << endl; - if (!evtloop2.Eventloop()) + + if (!calloop.Process(pedcam)) return; - tlist2.PrintStatistics(); - - + MCalibrationChargeCam &calcam = calloop.GetCalibrationCam(); + MCalibrationQECam &qecam = calloop.GetQECam(); + /************************************************************************/ /* THIRD LOOP: PEDESTAL COMPUTATION USING EXTRACTED SIGNAL (IN PHOTONS) */ @@ -202,6 +405,7 @@ plist3.AddToList(&geomcam); - plist3.AddToList(&pedcam); - plist3.AddToList(&calcam); + plist3.AddToList(&pedcam ); + plist3.AddToList(&calcam ); + plist3.AddToList(&qecam ); plist3.AddToList(&photcam); plist3.AddToList(&extedsig); @@ -209,10 +413,13 @@ //tasks - MReadMarsFile read3("Events", pedname); + MReadMarsFile read3("Events"); + read3.DisableAutoScheme(); + static_cast(read3).AddFiles(pruns); + + MExtractSignal sigcalc; MCalibrate photcalc; + photcalc.SetCalibrationMode(MCalibrate::kFfactor); MPedPhotCalc pedphotcalc; - read3.DisableAutoScheme(); - tlist3.AddToList(&read3); tlist3.AddToList(&geomapl); @@ -233,9 +440,35 @@ tlist3.PrintStatistics(); - /**********************/ /* PRODUCE NICE PLOTS */ /**********************/ + if (usedisplay) + { + + MHCamera disp0(geomcam, "MPedPhotCam;ped", "Pedestals"); + MHCamera disp1(geomcam, "MPedPhotCam;rms", "Sigma Pedestal"); + + disp0.SetCamContent(pedphotcam, 0); + disp0.SetCamError (pedphotcam, 1); + + disp1.SetCamContent(pedphotcam, 1); + + disp0.SetYTitle("Pedestal signal [photons]"); + disp1.SetYTitle("Pedestal signal RMS [photons]"); + + // + // Display data + // + TCanvas &c1 = display->AddTab("PedPhotCam"); + c1.Divide(2,3); + + CamDraw(c1, 1, 2, disp0, 0); + CamDraw(c1, 2, 2, disp1, 1); + + } + + +#if 0 const UShort_t npix = 577; @@ -250,4 +483,6 @@ TH1F* calhist = new TH1F("calhist","Calibration factors",npix,0,npix); TH1F* caldist = new TH1F("caldist","Calibration factors",100,0,1); + TH1F* qehist = new TH1F("qehist", "Quantrum efficiencies",npix,0,npix); + TH1F* qedist = new TH1F("qedist", "Quantrum efficiencies",100,0,1); // pedestal signals @@ -260,8 +495,11 @@ for(int i=0;iFill(i,cal); caldist->Fill(cal); + qehist->Fill(i,qe); + qedist->Fill(qe); pedphothist->Fill(i,pedphot); @@ -287,7 +527,6 @@ canvas->SetBorderMode(0); // Delete the Canvas' border line canvas->cd(); - canvas->SetBorderMode(0); - - canvas->Divide(2,4); + + canvas->Divide(2,5); // draw pedestal histo @@ -322,5 +561,4 @@ pedrmsdist->Draw("same"); - TLegend* leg2 = new TLegend(.14,.68,.39,.88); leg2->SetHeader(""); @@ -339,5 +577,5 @@ calhist->SetMaximum(1); calhist->SetMinimum(0); - calhist->GetYaxis()->SetTitle("Calibration factor (photon/ADC count)"); + calhist->GetYaxis()->SetTitle("Calibration factor (phe/ADC count)"); calhist->SetStats(kFALSE); calhist->Draw(); @@ -347,9 +585,27 @@ gPad->cd(); gPad->SetBorderMode(0); - caldist->GetXaxis()->SetTitle("Calibration factor (photon/ADC count)"); + caldist->GetXaxis()->SetTitle("Calibration factor (phe/ADC count)"); caldist->Draw(); + // draw qe histo + canvas->cd(5); + gPad->cd(); + gPad->SetBorderMode(0); + qehist->GetXaxis()->SetTitle("Pixel SW Id"); + qehist->SetMaximum(1); + qehist->SetMinimum(0); + qehist->GetYaxis()->SetTitle("Quantum efficiency for cascades"); + qehist->SetStats(kFALSE); + qehist->Draw(); + + // draw qe distribution + canvas->cd(6); + gPad->cd(); + gPad->SetBorderMode(0); + qedist->GetXaxis()->SetTitle("Quantum efficiency for cascades"); + qedist->Draw(); + // draw pedestal signal histo - canvas->cd(5); + canvas->cd(7); gPad->cd(); gPad->SetBorderMode(0); @@ -360,5 +616,5 @@ // draw pedestal signal distribution - canvas->cd(6); + canvas->cd(8); gPad->cd(); gPad->SetBorderMode(0); @@ -367,5 +623,5 @@ // draw pedestal signal rms histo - canvas->cd(7); + canvas->cd(9); gPad->cd(); gPad->SetBorderMode(0); @@ -376,5 +632,5 @@ // draw pedestal signal rms distribution - canvas->cd(8); + canvas->cd(10); gPad->cd(); gPad->SetBorderMode(0); @@ -382,5 +638,6 @@ pedphotrmsdist->Draw(); - return; + canvas->SaveAs("pedphotcalc.root"); + +#endif } -