Index: trunk/MagicSoft/Mars/mtemp/Changelog =================================================================== --- trunk/MagicSoft/Mars/mtemp/Changelog (revision 5126) +++ trunk/MagicSoft/Mars/mtemp/Changelog (revision 5127) @@ -18,4 +18,14 @@ -*-*- END OF LINE -*-*- + +2004/09/24: Hendrik Bartko + * mtemp/mmpi/macros/calibration_shape.C + - fine tuning of the plots + * mtemp/mmpi/macros/data_shape.C + - added new macro to compute the signal shape for cosmics runs + * mtemp/mmpi/macros/calculate_of_weights.C + - added new macro to compute the weights for the optimal filter + + 2004/09/19: Hendrik Bartko Index: trunk/MagicSoft/Mars/mtemp/mmpi/macros/calculate_of_weights.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mmpi/macros/calculate_of_weights.C (revision 5127) +++ trunk/MagicSoft/Mars/mtemp/mmpi/macros/calculate_of_weights.C (revision 5127) @@ -0,0 +1,152 @@ +const int no_samples=6; + + +void calculate_of_weights(TH1F *shape, TString noise_file = "/mnt/home/pcmagic17/hbartko/mars/results/04sep12/noise_autocorr_AB_36038.root", Int_t t_offset = -10) +{ + + + TH1F * derivative = new TH1F("derivative","derivative",161,-1.05,15.05); + + + for (int i = 1; i<162;i++){ + derivative->SetBinContent(i,((p->GetBinContent(i+1)-p->GetBinContent(i-1))/0.2)); + derivative->SetBinError(i,(sqrt(p->GetBinError(i+1)*p->GetBinError(i+1)+p->GetBinError(i-1)*p->GetBinError(i-1))/0.2)); + } + + // normalize the shape, such that the integral for 6 slices is one! + + float sum = 0; + for (int i=40; i<101; i++){sum+=shape->GetBinContent(i);} + sum /= 10; + + shape->Scale(1./sum); + derivative->Scale(1./sum); + + + TCanvas * c1 = new TCanvas("c1","c1",600,400); + //c1= canvas(); + shape->Draw(); + + TCanvas *c2 = new TCanvas("c2","c2",600,400); + //c2=canvas(); + derivative->Draw(); + + + // book the histograms for the weights + + TH1F * hw_amp = new TH1F("hw_amp","hw_amp",no_samples*10,-0.5,no_samples-0.5); + TH1F * hw_time = new TH1F("hw_time","hw_time",no_samples*10,-0.5,no_samples-0.5); + TH1F * hshape = new TH1F("hshape","hshape",no_samples*10,-0.5,no_samples-0.5); + TH1F * hderivative = new TH1F("hderivative","hderivative",no_samples*10,-0.5,no_samples-0.5); + + + // read in the noise auto-correlation function: + + TMatrix B(no_samples,no_samples); + + f = new TFile(noise_file); + TH2F * noise_corr = (TH2F*)c_corr->FindObject("hcorr"); + for (int i=0; iGetBinContent(i+1,j+1); + } + } + f->Close(); + + B.Invert(); + + // now the loop over t_{rel} in [-0.4;0.6[ : + + for (int i=-4; i<6; i++){ + + + TMatrix g(no_samples,1); + TMatrix gT(1,no_samples); + TMatrix d(no_samples,1); + TMatrix dT(1,no_samples); + + + for (int count=0; count < no_samples; count++){ + + g[count][0]=shape->GetBinContent(55+t_offset-int((10*no_samples-50)/2.)+10*count+i); + gT[0][count]=shape->GetBinContent(55+t_offset-int((10*no_samples-50)/2.)+10*count+i); + d[count][0]=derivative->GetBinContent(55+t_offset-int((10*no_samples-50)/2.)+10*count+i); + dT[0][count]=derivative->GetBinContent(55+t_offset-int((10*no_samples-50)/2.)+10*count+i); + + hshape->SetBinContent(i+5+10*count,shape->GetBinContent(55+t_offset-int((10*no_samples-50)/2.)+10*count+i)); + hderivative->SetBinContent(i+5+10*count,derivative->GetBinContent(55+t_offset-int((10*no_samples-50)/2.)+10*count+i)); + } + + + TMatrix m_denom = (gT*(B*g))*(dT*(B*d)) - (dT*(B*g))*(dT*(B*g)); + float denom = m_denom[0][0]; // m_denom is a real number + + TMatrix m_first = dT*(B*d); // m_first is a real number + float first = m_first[0][0]/denom; + + TMatrix m_last = gT*(B*d); // m_last is a real number + float last = m_last[0][0]/denom; + + TMatrix m1 = gT*B; + + m1 *=first; + + TMatrix m2 = dT*B; + + m2 *=last; + + TMatrix w_amp= m1 - m2; + + + TMatrix m_first1 = gT*(B*g); + float first1 = m_first1[0][0]/denom; + + TMatrix m_last1 = gT*(B*d); + float last1 = m_last1[0][0]/denom; + + TMatrix m11 = dT*B; + + m11 *=first1; + + TMatrix m21 = gT*B; + + m21 *=last1; + + + TMatrix w_time= m11 - m21; + + + float amp = 0; + float amp_time = 0; + + for (int count=0; count < no_samples; count++){ + hw_amp->SetBinContent(i+5+10*count,w_amp[0][count]); + hw_time->SetBinContent(i+5+10*count,w_time[0][count]); + } + + } // end loop over t_rel + + + TCanvas * c3 = new TCanvas("c3","c3",600,400); + hw_amp->Draw(); + + TCanvas * c4 = new TCanvas("c4","c4",600,400); + hw_time->Draw(); + + TCanvas * c5 = new TCanvas("c5","c5",600,400); + hshape->Draw(); + + TCanvas * c6 = new TCanvas("c6","c6",600,400); + hderivative->Draw(); + + /* + f_out = new TFile("/home/pcmagic17/hbartko/mars/results/04sep15/calibration_weights_corr_iterated.root","RECREATE"); //"../of_weights_bin4.root" + hw_amp->Write(); + hw_time->Write(); + hshape->Write(); + hderivative->Write(); + shape->Write(); + derivative->Write(); + f_out->Close(); + */ +} Index: trunk/MagicSoft/Mars/mtemp/mmpi/macros/calibration_shape.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mmpi/macros/calibration_shape.C (revision 5126) +++ trunk/MagicSoft/Mars/mtemp/mmpi/macros/calibration_shape.C (revision 5127) @@ -30,10 +30,10 @@ Int_t pedr = 36038; -Int_t calr = 36040; +Int_t calr = 36042; //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) +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=316) { @@ -112,8 +112,5 @@ plist2.AddToList(&badcam); - gLog << endl;; - gLog << "Calculate MCalibrationCam from Runs " << cruns.GetRunsAsString() << endl; - gLog << endl; - + MReadMarsFile read2("Events"); read2.DisableAutoScheme(); @@ -149,4 +146,5 @@ MExtractTimeHighestIntegral timeext; timeext.SetRange(0, 14, 0, 14); + //timeext.SetWindowSize(10,10); timeext.SetWindowSize(WindowSize,WindowSize); @@ -195,5 +193,5 @@ int count = 1; - TH2F * shape = new TH2F("shape","shape",30,0,30,1000,-20,250); + TH2F * shape = new TH2F("shape","shape",300,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); @@ -225,4 +223,7 @@ sigcam.GetPixelContent(sig, ipix, geomcam, 0); + Byte_t sat = (sigcam)[ipix].GetNumHiGainSaturated(); + + htime->Fill(time); hsig->Fill(sig); @@ -243,4 +244,5 @@ const Float_t ABoffs = pedcam[ipix].GetPedestalABoffset(); + Float_t PedMean[2]; PedMean[0] = ped_mean + ABoffs; @@ -254,7 +256,9 @@ 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); + shape_corr->Fill(sample+15.5+5.-time,logains[sample]-PedMean[(sample+nh+ABFlag)&0x1]); + if (sat==0){ + 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+nh+ABFlag)&0x1])/sig*250); + } } @@ -350,4 +354,24 @@ htime->Draw(); + + TCanvas * c7 = new TCanvas("c7","c7",600,400); + c7->cd(); + c7->SetGridx(); + c7->SetGridy(); + c7->SetFillColor(0); + // c7->SetBordermode(0); + + TProfile * shape_corr_all_pfx = (TProfile*)shape_corr_all->ProfileX(); + + shape_corr_all_pfx->SetStats(0); + TString title6 = "Average Calibration Signal Shape, Amplitude + Arrival Time Corrected, Pixel"; + title6+=ipix; + shape_corr_all_pfx->SetTitle(title6); + shape_corr_all_pfx->SetXTitle("FADC sample no. + ( - t_{arrival}) / T_{ADC}"); + shape_corr_all_pfx->SetYTitle("signal [a.u.]"); + shape_corr_all_pfx->SetMarkerStyle(20); + shape_corr_all_pfx->SetMarkerColor(4); + shape_corr_all_pfx->SetMarkerSize(0.5); + shape_corr_all_pfx->Draw(); } Index: trunk/MagicSoft/Mars/mtemp/mmpi/macros/data_shape.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mmpi/macros/data_shape.C (revision 5127) +++ trunk/MagicSoft/Mars/mtemp/mmpi/macros/data_shape.C (revision 5127) @@ -0,0 +1,377 @@ +/* ======================================================================== *\ +! +! * +! * 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 = 36058; +Int_t calr = 36058; +//Int_t datar = 12517; + + + +void data_shape(const TString inpath = "/.magic/magicserv01/scratch/hbartko/rootdata/2004_09_06/", Int_t pedruns = pedr, Int_t calruns = calr, Int_t ipix=316) + +{ + + + 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); + + + 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,-50,250); + TH2F * shape_corr = new TH2F("shape_corr","shape_corr",300,0,30,1000,-50,250); + TH2F * shape_corr_all = new TH2F("shape_corr_all","shape_corr_all",300,0,30,1000,-50,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); + + Byte_t sat = (sigcam)[ipix].GetNumHiGainSaturated(); + + hsig->Fill(sig); + + if (sig>60 && sat==0){ + + htime->Fill(time); + + 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+6.-time,higains[sample]-PedMean[(sample+ABFlag)&0x1]); + shape_corr->Fill(sample+15.5+6.-time,logains[sample]-PedMean[(sample+nh+ABFlag)&0x1]); + shape_corr_all->Fill(sample+0.5+6.-time,(higains[sample]-PedMean[(sample+ABFlag)&0x1])/sig*250); + shape_corr_all->Fill(sample+15.5+6.-time,(logains[sample]-PedMean[(sample+nh+ABFlag)&0x1])/sig*250); + } + }// end if (sig > 50) + + } // 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 Data 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 = "Data 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 = "Data 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 = "Data 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 = "Data 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(); + + TCanvas * c7 = new TCanvas("c7","c7",600,400); + c7->cd(); + c7->SetGridx(); + c7->SetGridy(); + c7->SetFillColor(0); + // c7->SetBordermode(0); + + TProfile * shape_corr_all_pfx = (TProfile*)shape_corr_all->ProfileX(); + + shape_corr_all_pfx->SetStats(0); + TString title6 = "Average Data Signal Shape, Amplitude + Arrival Time Corrected, Pixel"; + title6+=ipix; + shape_corr_all_pfx->SetTitle(title6); + shape_corr_all_pfx->SetXTitle("FADC sample no. + ( - t_{arrival}) / T_{ADC}"); + shape_corr_all_pfx->SetYTitle("signal [a.u.]"); + shape_corr_all_pfx->SetMarkerStyle(20); + shape_corr_all_pfx->SetMarkerColor(4); + shape_corr_all_pfx->SetMarkerSize(0.5); + shape_corr_all_pfx->Draw(); + +} + + + +