Index: /trunk/MagicSoft/Mars/mtemp/Changelog =================================================================== --- /trunk/MagicSoft/Mars/mtemp/Changelog (revision 5118) +++ /trunk/MagicSoft/Mars/mtemp/Changelog (revision 5119) @@ -23,4 +23,6 @@ - added new macro to compute the noise auto correlation matrix from a pedestal run + * mtemp/mmpi/macros/calibration_shape.C + - added new macro to compute the signal shape for calibration runs Index: /trunk/MagicSoft/Mars/mtemp/mmpi/macros/calibration_shape.C =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mmpi/macros/calibration_shape.C (revision 5119) +++ /trunk/MagicSoft/Mars/mtemp/mmpi/macros/calibration_shape.C (revision 5119) @@ -0,0 +1,356 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ +#include "MAGIC.h" + + +const TString pedfile="/.magic/magicserv01/MAGIC/rootdata/2004_01_27/20040126_12149_P_Cab-On_E.root"; +const TString calfile="/.magic/magicserv01/MAGIC/rootdata/2004_01_27/20040126_12526_C_Cab-On_E.root"; +const TString datafile="/.magic/magicserv01/MAGIC/rootdata/2004_01_27/20040126_12517_D_Cab-On_E.root"; + +Int_t pedr = 36038; +Int_t calr = 36040; +//Int_t datar = 12517; + + + +void calibration_shape(const TString inpath = "/.magic/magicserv01/scratch/hbartko/rootdata/2004_09_06/", Int_t pedruns = pedr, Int_t calruns = calr, Int_t ipix=297) + +{ + + + MRunIter pruns; + MRunIter cruns; + + pruns.AddRun(pedruns,inpath); + + cruns.AddRun(calruns,inpath); + + + gStyle->SetOptStat(1111); + gStyle->SetOptFit(); + + /************************************/ + /* FIRST LOOP: PEDESTAL COMPUTATION */ + /************************************/ + + MParList plist1; + MTaskList tlist1; + plist1.AddToList(&tlist1); + + // containers + MPedestalCam pedcam; + MBadPixelsCam badcam; + + + plist1.AddToList(&pedcam); + plist1.AddToList(&badcam); + + //tasks + MReadMarsFile read("Events"); + read.DisableAutoScheme(); + static_cast(read).AddFiles(pruns); + + MGeomApply geomapl; + // MPedCalcPedRun pedcalc; + MPedCalcFromLoGain pedcalc_ped; + pedcalc_ped.SetMaxHiGainVar(20); + pedcalc_ped.SetRange(0, 11, 1, 14); + pedcalc_ped.SetWindowSize(12,14); + pedcalc_ped.SetPedestalUpdate(kFALSE); + + + + MGeomCamMagic geomcam; + + tlist1.AddToList(&read); + tlist1.AddToList(&geomapl); + tlist1.AddToList(&pedcalc_ped); + + // Create and execute the event looper + MEvtLoop pedloop; + pedloop.SetParList(&plist1); + + cout << "*************************" << endl; + cout << "** COMPUTING PEDESTALS **" << endl; + cout << "*************************" << endl; + + if (!pedloop.Eventloop()) + return; + + tlist1.PrintStatistics(); + + + + // + // Create a empty Parameter List and an empty Task List + // + MParList plist2; + MTaskList tlist2; + plist2.AddToList(&tlist2); + plist2.AddToList(&pedcam); + plist2.AddToList(&badcam); + + gLog << endl;; + gLog << "Calculate MCalibrationCam from Runs " << cruns.GetRunsAsString() << endl; + gLog << endl; + + MReadMarsFile read2("Events"); + read2.DisableAutoScheme(); + static_cast(read2).AddFiles(cruns); + + MGeomCamMagic geomcam; + MExtractedSignalCam sigcam; + MArrivalTimeCam timecam; + MRawEvtData evtdata; + + // + // Get the previously created MPedestalCam into the new Parameter List + // + plist2.AddToList(&geomcam); + plist2.AddToList(&sigcam); + plist2.AddToList(&timecam); + plist2.AddToList(&evtdata); + + // + // We saw that the signal jumps between slices, + // thus take the sliding window + // + + MGeomApply geomapl; + + + Int_t WindowSize = 6; // size of the chosen integration window for signal and time extraction + + MExtractFixedWindowPeakSearch extractor; + extractor.SetRange(0, 14, 0, 14); + extractor.SetWindows( WindowSize, WindowSize, 4); + + MExtractTimeHighestIntegral timeext; + timeext.SetRange(0, 14, 0, 14); + timeext.SetWindowSize(WindowSize,WindowSize); + + + // timecalc->SetRange(8,14,8,14,6); + + // filter against cosmics + MFCosmics cosmics; + MContinue cont(&cosmics); + + tlist2.AddToList(&read2); + tlist2.AddToList(&geomapl); + tlist2.AddToList(&extractor); + tlist2.AddToList(&timeext); + + // + // In case, you want to skip the cosmics rejection, + // uncomment the next line + // + tlist2.AddToList(&cont); + // + // In case, you want to skip the somewhat lengthy calculation + // of the arrival times using a spline, uncomment the next two lines + // + + + // + // Create and setup the eventloop + // + MEvtLoop evtloop; + evtloop.SetParList(&plist2); + // evtloop.SetDisplay(display); + + cout << "***************************" << endl; + cout << "** COMPUTING CALIBRATION **" << endl; + cout << "***************************" << endl; + + // + // Execute second analysis + // + + + if (!evtloop.PreProcess()) + return; + + int count = 1; + + TH2F * shape = new TH2F("shape","shape",30,0,30,1000,-20,250); + TH2F * shape_corr = new TH2F("shape_corr","shape_corr",300,0,30,1000,-20,250); + TH2F * shape_corr_all = new TH2F("shape_corr_all","shape_corr_all",300,0,30,1000,-20,250); + TH1F * htime = new TH1F("htime","htime",150,0,15); + TH1F * hsig = new TH1F("hsig","hsig",300,-50,550); + + // cout << " hello " << endl; + + const Float_t ped_mean = pedcam[ipix].GetPedestal(); + const Float_t pedRMS = pedcam[ipix].GetPedestalRms(); + + while (tlist2.Process()){ //loop over the events + + count++; + + if (count%1000==0) cout << "Event #" << count << endl; + + MRawEvtPixelIter pixel(&evtdata); + + + double value = 0; + pedcam.GetPixelContent(value, ipix, geomcam, 0); + + + double time = 0; + timecam.GetPixelContent(time, ipix, geomcam, 0); // was 0) + + double sig = 0; + sigcam.GetPixelContent(sig, ipix, geomcam, 0); + + htime->Fill(time); + hsig->Fill(sig); + + + pixel.Jump(ipix); + + const Byte_t *higains = pixel.GetHiGainSamples(); + const Byte_t *logains = pixel.GetLoGainSamples(); + const Int_t nh = evtdata.GetNumHiGainSamples(); + const Int_t nl = evtdata.GetNumLoGainSamples(); + + Bool_t ABFlag = pixel.HasABFlag(); + + const Byte_t *higains = pixel.GetHiGainSamples(); + const Byte_t *logains = pixel.GetLoGainSamples(); + + const Float_t ABoffs = pedcam[ipix].GetPedestalABoffset(); + + Float_t PedMean[2]; + PedMean[0] = ped_mean + ABoffs; + PedMean[1] = ped_mean - ABoffs; + + + + + for (int sample=0; sampleFill(sample+0.5,higains[sample]-PedMean[(sample+ABFlag)&0x1]); + shape->Fill(sample+15.5,logains[sample]-PedMean[(sample+ABFlag)&0x1]); + shape_corr->Fill(sample+0.5+5.-time,higains[sample]-PedMean[(sample+ABFlag)&0x1]); + shape_corr->Fill(sample+15.5+5.-time,logains[sample]-PedMean[(sample+ABFlag)&0x1]); + shape_corr_all->Fill(sample+0.5+5.-time,(higains[sample]-PedMean[(sample+ABFlag)&0x1])/sig*250); + shape_corr_all->Fill(sample+15.5+5.-time,(logains[sample]-PedMean[(sample+ABFlag)&0x1])/sig*250); + } + + + } // end loop over the entries + + tlist2.PrintStatistics(); + + + TCanvas * c_shape = new TCanvas("c_shape","c_shape",600,400); + c_shape->SetFillColor(0); + c_shape->SetBorderMode(0); + c_shape->SetGridx(); + c_shape->SetGridy(); + shape->SetStats(0); + TString title2 = "Raw Calibration Signal Shape, Pixel "; + title2+=ipix; + shape->SetTitle(title2); + shape->SetXTitle("FADC sample no."); + shape->SetYTitle("signal [FADC counts]"); + shape->SetMarkerStyle(20); + shape->SetMarkerColor(4); + shape->SetMarkerSize(0.3); + shape->Draw(); + + + TCanvas * c_shape_corr = new TCanvas("c_shape_corr","c_shape_corr",600,400); + c_shape_corr->SetFillColor(0); + c_shape_corr->SetBorderMode(0); + c_shape_corr->SetGridx(); + c_shape_corr->SetGridy(); + shape_corr->SetStats(0); + TString title3 = "Calibration Signal Shape, Arrival Time Corrected, Pixel "; + title3+=ipix; + shape_corr->SetTitle(title3); + shape_corr->SetXTitle("FADC sample no. + ( - t_{arrival}) / T_{ADC}"); + shape_corr->SetYTitle("signal [FADC counts]"); + shape_corr->SetMarkerStyle(20); + shape_corr->SetMarkerColor(4); + shape_corr->SetMarkerSize(0.3); + shape_corr->Draw(); + + + TCanvas * c_shape_corr_all = new TCanvas("c_shape_corr_all","c_shape_corr_all",600,400); + c_shape_corr_all->SetFillColor(0); + c_shape_corr_all->SetBorderMode(0); + c_shape_corr_all->SetGridx(); + c_shape_corr_all->SetGridy(); + shape_corr_all->SetStats(0); + TString title3 = "Calibration Signal Shape, Arrival Time + Amplitude Corrected, Pixel "; + title3+=ipix; + shape_corr_all->SetTitle(title3); + shape_corr_all->SetXTitle("FADC sample no. + ( - t_{arrival}) / T_{ADC}"); + shape_corr_all->SetYTitle("signal [a.u.]"); + shape_corr_all->SetMarkerStyle(20); + shape_corr_all->SetMarkerColor(4); + shape_corr_all->SetMarkerSize(0.3); + shape_corr_all->Draw(); + + + + + TCanvas * c_sig = new TCanvas("c_sig","c_sig",600,400); + c_sig->SetFillColor(0); + c_sig->SetBorderMode(0); + c_sig->SetGridx(); + c_sig->SetGridy(); + + hsig->SetStats(0); + TString title4 = "Calibration Extracted Charge, MExtractFixedWindowPeakSearch, Pixel "; + title4+=ipix; + hsig->SetTitle(title4); + hsig->SetXTitle("reco charge [FADC counts]"); + hsig->SetYTitle("events / 2 FADC counts"); + hsig->SetLineColor(4); + hsig->SetLineWidth(2); + hsig->Draw(); + + TCanvas * c_time = new TCanvas("c_time","c_time",600,400); + c_time->SetFillColor(0); + c_time->SetBorderMode(0); + c_time->SetGridx(); + c_time->SetGridy(); + + htime->SetStats(0); + TString title5 = "Calibration Arrival Time, MExtractTimeHighestIntegral, Pixel "; + title5+=ipix; + htime->SetTitle(title5); + htime->SetXTitle("arrival time [T_{ADC}]"); + htime->SetYTitle("events / 0.1 T_{ADC}"); + htime->SetLineColor(4); + htime->SetLineWidth(2); + htime->Draw(); + +} + + + +