/* ======================================================================== *\ ! ! * ! * 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): Hendrik Bartko 09/2004 ! ! Copyright: MAGIC Software Development, 2000-2001 ! ! \* ======================================================================== */ const TString defname = "/.magic/magicserv01/scratch/hbartko/rootdata/2004_09_03/20040903_35698_P_multiplexFADC_E.root"; const TString defpedname = "/.magic/magicserv01/scratch/hbartko/rootdata/2004_09_03/20040903_35698_P_multiplexFADC_E.root"; // works with one pedestal file to calculate the noise autocorrelation void noise_autocorrelation_AB(const TString fname = defpedname, const TString pedname = defpedname, Int_t ipix = 67) { MParList plist_ped; MTaskList tlist_ped; plist_ped.AddToList(&tlist_ped); MPedestalCam pedcam; plist_ped.AddToList(&pedcam); MReadMarsFile read("Events", pedname); read.DisableAutoScheme(); tlist_ped.AddToList(&read); MGeomApply geomapl_ped; MGeomCamMagic geomcam; tlist_ped.AddToList(&geomapl_ped); MPedCalcFromLoGain pedcalc_ped; //pedcalc_ped.SetMaxHiGainVar(20); //pedcalc_ped.SetRange(0, 11, 1, 14); //pedcalc_ped.SetWindowSize(12,14); pedcalc_ped.SetPedestalUpdate(kFALSE); tlist_ped.AddToList(&pedcalc_ped); MEvtLoop evtloop_ped; evtloop_ped.SetParList(&plist_ped); if (!evtloop_ped.Eventloop()) return; tlist_ped.PrintStatistics(); // now the event loop for the signal reconstruction with pedestals subtracted MParList plist; MTaskList tlist; plist.AddToList(&tlist); plist.AddToList(&pedcam); MRawRunHeader runheader; plist.AddToList(&runheader); MRawEvtData evtdata; plist.AddToList(&evtdata); MExtractedSignalCam sigcam; plist.AddToList(&sigcam); MArrivalTimeCam timecam; plist.AddToList(&timecam); MReadMarsFile read("Events", fname); read.DisableAutoScheme(); 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); tlist.AddToList(&read); tlist.AddToList(&geomapl); tlist.AddToList(&extractor); tlist.AddToList(&timeext); MEvtLoop evtloop; evtloop.SetParList(&plist); if (!evtloop.PreProcess()) return; TString title = "Noise Autocorrelation B_{ij} = - *, Pixel "; title += ipix; Int_t First = 1; TH2F * hcorr = new TH2F("hcorr","hcorr",15,0,15,15,0,15); TH2F * shape = new TH2F("shape","shape",311,-1.05,30.05,1000,-20,250); TH2F * shape_corr = new TH2F("shape_corr","shape_corr",311,-1.05,30.05,1000,-20,250); TH1F * htime = new TH1F("htime","htime",150,0,15); TH1F * hsig = new TH1F("hsig","hsig",300,-50,550); Float_t B[15][15]; Float_t bi[15]; Float_t bj[15]; for (int i=0; i<15; i++){ bi[i]=0; bj[i]=0; for (int j=0; j<15; j++){ B[i][j]=0; } } const Float_t ped_mean = pedcam[ipix].GetPedestal(); const Float_t pedRMS = pedcam[ipix].GetPedestalRms(); cout << " ped_mean " << ped_mean << " pedRMS " << pedRMS << endl; int count = 0; while (tlist.Process()) { if (First) { First = 0; const Int_t nh = runheader.GetNumSamplesHiGain(); const Int_t nl = runheader.GetNumSamplesLoGain(); const Int_t nt = nh+nl; } htime->Fill(timecam[ipix].GetArrivalTimeHiGain()); hsig->Fill(sigcam[ipix].GetExtractedSignalHiGain()); MRawEvtPixelIter pixel(&evtdata); pixel.Jump(ipix); Bool_t ABFlag = pixel.HasABFlag(); count++; // couter for the number of events 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; // calculate the noise autocorrelation matrix for (int slice=0; sliceFill(slice,higains[slice]-PedMean[(slice+ABFlag)&0x1]); shape->Fill(slice+nh,logains[slice]-PedMean[(slice+ABFlag)&0x1]); if (trigger==0 && (higains[slice]-PedMean[(slice+ABFlag)&0x1]) >= 4.0) trigger = slice; } for (int slice=0; slice=1 && slice<=trigger+2){ charge += higains[slice]-PedMean[(slice+ABFlag)&0x1]; charge_time += slice*(higains[slice]-PedMean[(slice+ABFlag)&0x1]); } if (charge>0) charge_time /= charge; for (int slice=0; slice=1) shape_corr->Fill(slice+3-charge_time,(higains[slice]-PedMean[(slice+ABFlag)&0x1])); } // end eventloop.Process() evtloop.PostProcess(); for (int i=0; i<15; i++){ for (int j=0; j<15; j++){ hcorr->Fill(i,j,B[i][j]/count-bi[i]/count*bj[j]/count); // cout << "i " << i << " j " << j << " B[i][j] " << B[i][j] << " bi[i] " << bi[i] << " bj[j] " << bj[j] << endl; } } cout << " pedmean " << ped_mean << " pedRMS " << pedRMS << endl; TCanvas * c_corr = new TCanvas("c_corr","c_corr",600,400); c_corr->SetFillColor(0); c_corr->SetBorderMode(0); c_corr->SetGridx(); c_corr->SetGridy(); gStyle->SetPalette(1); hcorr->SetStats(0); hcorr->SetTitle(title); hcorr->SetXTitle("slice i"); hcorr->SetYTitle("slice j"); hcorr->Draw("colz"); TCanvas * c_corr_px = new TCanvas("c_corr_px","c_corr_px",600,400); c_corr_px->SetFillColor(0); c_corr_px->SetBorderMode(0); c_corr_px->SetGridx(); c_corr_px->SetGridy(); TH1F * hcorr_px = (TH1F*)hcorr->ProjectionX("hcorr_px",6,6); hcorr_px->SetLineColor(4); hcorr_px->SetLineWidth(2); hcorr_px->SetStats(0); TString title1 = "Noise Autocorrelation, FADC sample 6, Pixel"; title1+=ipix; hcorr_px->SetTitle(title); hcorr_px->SetXTitle("slice i"); hcorr_px->SetYTitle("B_{ij} = - * [FADC counts]"); hcorr_px->Draw(); 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 = "Noise Signal Shape, Signal #geq 3 FADC counts, 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 = "Noise Signal Shape, Signal #geq 3 FADC counts, 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_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 = "Noise 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 = "Noise 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(); }