#include #include #include "MLog.h" #include "MLogManip.h" #if !defined(__CINT__) || defined(__MAKECINT__) #include "TH1F.h" #include "TFile.h" #include "TStyle.h" #include "TGraph.h" #include "TLine.h" #include "../mcore/DrsCalib.h" #include "MDrsCalibration.h" #include "MExtralgoSpline.h" #include "MSequence.h" #include "MStatusArray.h" #include "MHCamera.h" #include "MJob.h" #include "MWriteRootFile.h" #include "MHCamera.h" #include "MBadPixelsCam.h" #include "MBadPixelsPix.h" #include "MDirIter.h" #include "MTaskList.h" #include "MFDataPhrase.h" #include "MArrayF.h" #include "MBadPixelsTreat.h" #include "MCalibrateDrsTimes.h" #include "MHSectorVsTime.h" #include "MHCamEvent.h" #include "MExtractTimeAndChargeSpline.h" #include "MFillH.h" #include "MDrsCalibApply.h" #include "MGeomApply.h" #include "MContinue.h" #include "MRawFitsRead.h" #include "MReadMarsFile.h" #include "MEvtLoop.h" #include "MParList.h" #include "MStatusDisplay.h" #include "MDrsCalibrationTime.h" #include "MH3.h" #include "MGeomCamFACT.h" #include "MCalibrateFact.h" #include "MParameters.h" #include "MWriteAsciiFile.h" #endif using namespace std; /* Maybe you wanna use this macro like this: * * 0.) ---- call root ---- * root -b * * 1.) ---- compile the stuff ---- * .L fact/analysis/callisto_buildable_no_sequence_file.C++ * * * 2.) ---- you can call it then ---- * Therefore you need to specify all the paths ... see below. * * When you wanna call the stuff directly from the bash make sure to * escape the bracets and quotes correctly. * * your can do: * root -b -q callisto_buildable_no_sequence_file.C++'("path1","path2",...)' * or: * root -b -q callisto_buildable_no_sequence_file.C++(\"path1\",\"$HOME\",...) * using bash enviroment variables like $HOME is not possible in the upper variant. */ int callisto(const TString drsfile="test300samples.drs.fits", const TString pedfile="00000001.001_P_MonteCarlo000_Events.fits", const TString datfile="00000003.387_D_MonteCarlo010_Events.fits", TString outfile = "", TString displayfile = "", TString displaytitle = "") { // ====================================================== if (displaytitle.IsNull()) displaytitle = gSystem->BaseName(datfile); FileStat_t fstat; int rc = gSystem->GetPathInfo(outfile, fstat); bool isdir = !rc || R_ISDIR(fstat.fMode); const char *buf = gSystem->ConcatFileName(outfile, "callisto.root"); outfile = buf; delete [] buf; if (displayfile.IsNull()) { displayfile = outfile; displayfile.Insert(displayfile.Last('.'), "-display"); } else { if (isdir && gSystem->DirName(displayfile)==TString(".")) { buf = gSystem->ConcatFileName(outfile, displayfile); displayfile = buf; delete [] buf; } } // ====================================================== // true: Display correctly mapped pixels in the camera displays // but the value-vs-index plot is in software/spiral indices // false: Display pixels in hardware/linear indices, // but the order is the camera display is distorted bool usemap = true; // map file to use (get that from La Palma!) const char *pmap = usemap ? "TestForThomas/FACT/FACTmap111030.txt" : NULL; Bool_t maximum = kTRUE; //const char *lp_template = maximum ? // "/cm/shared/apps/fact/Mars_svn_LP/template-lp-extractor-maximum.root" : // "/cm/shared/apps/fact/Mars_svn_LP/template-lp-extractor-leading-edge.root"; const char *pulse_template = "TestForThomas/FACT/template-pulse.root"; // ------------------------------------------------------ // Calib: 51 / 90 / 197 (20% TH) // Data: 52 / 64 / 104 (20% TH) // Extraction range in slices. It will always(!) contain the full range // of integration const int first_slice = 20; // 10ns const int last_slice = 250; // 125ns // Note that rise and fall time mean different things whether you use IntegralFixed or IntegraRel: // // IntegralFixed: // * fRiseTime: Number of slices left from arrival time // * fFallTime: Number of slices right from arrival time // IntegralRel: // * fRiseTime: Number of slices left from maximum time // * fFallTime: Number of slices right from maximum time // const int rise_time_cal = maximum ? 40 : 10; // was 13; 5ns const int fall_time_cal = maximum ? 120 : 160; // was 23; 80ns const int rise_time_dat = maximum ? 10 : 2; // was 13; was 10; 1ns const int fall_time_dat = maximum ? 40 : 48; // was 23; was 40; 24ns // Extraction type: Extract integral and half leading edge const MExtralgoSpline::ExtractionType_t type = maximum ? (MExtralgoSpline::kIntegralRel) : (MExtralgoSpline::kIntegralFixed); //const int type = MExtralgoSpline::kIntegralFixed; const double heighttm = 0.5; // IntegralAbs { 1.5pe * 9.6mV/pe } / IntegralRel { 0.5 } Long_t max = 0; // All Long_t max3 = max; // Pedestal Rndm Long_t max4 = max; // Pedestal Ext // ====================================================== if (pmap && gSystem->AccessPathName(pmap, kFileExists)) { gLog << err << "ERROR - Cannot access mapping file '" << pmap << "'" << endl; return 1; } gLog.Separator("Callisto"); gLog << all; gLog << "Data File: " << datfile << '\n'; gLog << "DRS calib 300: " << drsfile << endl;; MDrsCalibration drscalib300; if (!drscalib300.ReadFits(drsfile.Data())) { gLog << err << "ERROR - Cannot access drscallib300 file '" << drsfile << "'" << endl; return 5; } gLog << all; gLog << "Pedestal file: " << pedfile << '\n'; gLog << "Output file: " << outfile << '\n'; gLog << "Display file: " << displayfile << '\n'; gLog << "Display title: " << displaytitle << endl; // ------------------------------------------------------ MStatusArray arrt, arrp; // TFile ft(lp_template); // if (arrt.Read()<=0) // { // gLog << err << "ERROR - Reading LP template from " << lp_template << endl; // return 100; // } // MHCamera *lpref = (MHCamera*)arrt.FindObjectInCanvas("ExtCalSig;avg", "MHCamera", "Cam"); // if (!lpref) // { // gLog << err << "ERROR - LP Template not found in " << lp_template << endl; // return 101; // } // lpref->SetDirectory(0); // MHCamera *gain = (MHCamera*)arrt.FindObjectInCanvas("gain", "MHCamera", "Gain"); // if (!gain) // { // gLog << err << "ERROR - Gain not found in " << lp_template << endl; // return 101; // } // gain->SetDirectory(0); TFile fp(pulse_template); if (arrp.Read()<=0) { gLog << err << "ERROR - Reading Pulse template from " << pulse_template << endl; return 102; } TH1F *hpulse = (TH1F*)arrp.FindObjectInCanvas("hPixelEdgeMean0_0", "TH1F", "cgpPixelPulses0"); if (!hpulse) { gLog << err << "ERROR - Pulse Template not found in " << pulse_template << endl; return 103; } hpulse->SetDirectory(0); // ====================================================== MStatusDisplay *d = new MStatusDisplay; MBadPixelsCam badpixels; badpixels.InitSize(1440); /* badpixels[ 424].SetUnsuitable(MBadPixelsPix::kUnsuitable); badpixels[ 583].SetUnsuitable(MBadPixelsPix::kUnsuitable); badpixels[ 830].SetUnsuitable(MBadPixelsPix::kUnsuitable); badpixels[ 923].SetUnsuitable(MBadPixelsPix::kUnsuitable); badpixels[1208].SetUnsuitable(MBadPixelsPix::kUnsuitable); badpixels[1399].SetUnsuitable(MBadPixelsPix::kUnsuitable); */ // Twin pixel // 113 // 115 // 354 // 423 // 1195 // 1393 MDrsCalibrationTime timecam; // Plot the trigger pattern rates vs. run-number MH3 hrate("MRawRunHeader.GetFileID", "MRawEvtHeader.GetTriggerID&0xff00"); hrate.SetWeight("1./TMath::Max(MRawRunHeader.GetRunLength,1)"); hrate.SetName("Rate"); hrate.SetTitle("Event rate [Hz];File Id;Trigger Type;"); hrate.InitLabels(MH3::kLabelsXY); hrate.DefineLabelY( 0, "Data"); // What if TriggerID==0 already??? hrate.DefineLabelY(0x100, "Cal"); hrate.DefineLabelY(0x400, "Ped"); // hrate.DefaultLabelY("ERROR"); gStyle->SetOptFit(kTRUE); // ========================= Result ================================== //~ Double_t avgS = evt1f.GetHist()->GetMean(); //~ Double_t medS = evt1f.GetHist()->GetMedian(); //~ Double_t rmsS = evt1f.GetHist()->GetRMS(); //~ Double_t maxS = evt1f.GetHist()->GetMaximum(); MArrayF der1(hpulse->GetNbinsX()); MArrayF der2(hpulse->GetNbinsX()); MExtralgoSpline spline(hpulse->GetArray()+1, hpulse->GetNbinsX(), der1.GetArray(), der2.GetArray()); spline.SetRiseFallTime(rise_time_dat, fall_time_dat); spline.SetExtractionType(type); spline.SetHeightTm(heighttm); spline.Extract(hpulse->GetMaximumBin()-1); // The pulser signal is most probably around 400mV/9.5mV // IntegraFixed 2/48 corresponds to roughly 215mV*50slices Double_t scale = 1./spline.GetSignal(); MArrayD calib(1440); for (int i=0; i<1440; i++) calib[i] =1.; gROOT->SetSelectedPad(0); d->AddTab("PulseTemp"); gPad->SetGrid(); hpulse->SetNameTitle("Pulse", "Single p.e. pulse template"); hpulse->SetDirectory(0); hpulse->SetLineColor(kBlack); hpulse->DrawCopy(); TAxis *ax = hpulse->GetXaxis(); Double_t w = hpulse->GetBinWidth(1); Double_t T = w*(spline.GetTime()+0.5) +ax->GetXmin(); //~ Double_t H = w*(hpulse->GetMaximumBin()+0.5)+ax->GetXmin(); TLine line; line.SetLineColor(kRed); line.DrawLine(T-rise_time_dat*w, spline.GetHeight(), T+fall_time_dat*w, spline.GetHeight()); line.DrawLine(T, spline.GetHeight()/4, T, 3*spline.GetHeight()/4); line.DrawLine(T-rise_time_dat*w, 0, T-rise_time_dat*w, spline.GetHeight()); line.DrawLine(T+fall_time_dat*w, 0, T+fall_time_dat*w, spline.GetHeight()); TGraph gg; for (int ix=1; ix<=hpulse->GetNbinsX(); ix++) for (int i=0; i<10; i++) { Double_t x = hpulse->GetBinLowEdge(ix)+i*hpulse->GetBinWidth(ix)/10.; gg.SetPoint(gg.GetN(), x+w/2, spline.EvalAt(ix-1+i/10.)); } gg.SetLineColor(kBlue); gg.SetMarkerColor(kBlue); gg.SetMarkerStyle(kFullDotMedium); gg.DrawClone("L"); gROOT->SetSelectedPad(0); d->AddTab("CalConst"); MGeomCamFACT fact; MHCamera hcalco(fact); hcalco.SetName("CalConst"); hcalco.SetTitle(Form("Relative calibration constant [%.0f/pe]", 1./scale)); hcalco.SetCamContent(calib); hcalco.SetAllUsed(); //hcalco.Scale(scale); hcalco.DrawCopy(); // ====================================================== gLog << endl; gLog.Separator("Extracting random pedestal"); MTaskList tlist3; MParList plist3; plist3.AddToList(&tlist3); plist3.AddToList(&drscalib300); plist3.AddToList(&badpixels); plist3.AddToList(&timecam); MEvtLoop loop3("DetermineRndmPed"); loop3.SetDisplay(d); loop3.SetParList(&plist3); // ------------------ Setup the tasks --------------- MRawFitsRead read3; read3.LoadMap(pmap); read3.AddFile(pedfile); MFillH fill3a(&hrate); MContinue cont3("(MRawEvtHeader.GetTriggerID&0xff00)!=0x400", "SelectPed"); MGeomApply apply3; MDrsCalibApply drsapply3; //--- MExtractTimeAndChargeSpline extractor3; extractor3.SetRange(first_slice, last_slice); extractor3.SetRiseTimeHiGain(rise_time_dat); extractor3.SetFallTimeHiGain(fall_time_dat); extractor3.SetHeightTm(heighttm); extractor3.SetChargeType(type); extractor3.SetSaturationLimit(600000); extractor3.SetNoiseCalculation(kTRUE); // MHCamEvent evt2a(0, "PedRdm", "Extracted Pedestal Signal;;S"); // MFillH fill2a(&evt2a, "MExtractedSignalCam", "FillPedRndm"); // Use this for data, but not for calibration events // evt2a.SetErrorSpread(kFALSE); /* MCalibrateData conv3; conv3.SetCalibrationMode(MCalibrateData::kNone); conv3.SetPedestalFlag(MCalibrateData::kNo); conv3.SetCalibConvMinLimit(0); conv3.SetCalibConvMaxLimit(10000); conv3.SetScaleFactor(scale); */ MCalibrateFact conv3; conv3.SetScale(scale); conv3.SetCalibConst(calib); MBadPixelsTreat treat3; treat3.SetProcessPedestalRun(kFALSE); treat3.SetProcessPedestalEvt(kFALSE); treat3.SetProcessTimes(kFALSE); MHCamEvent evt3b(0, "PedRdm","Interpolated random pedestal;;Signal [~phe]"); //evt2b.SetErrorSpread(kFALSE); MFillH fill3b(&evt3b, "MSignalCam", "FillPedRdm"); fill3b.SetDrawOption("gaus"); // ------------------ Setup eventloop and run analysis --------------- tlist3.AddToList(&read3); tlist3.AddToList(&apply3); tlist3.AddToList(&drsapply3); tlist3.AddToList(&cont3); tlist3.AddToList(&extractor3); // tlist3.AddToList(&fill3a); tlist3.AddToList(&conv3); tlist3.AddToList(&treat3); tlist3.AddToList(&fill3b); if (!loop3.Eventloop(max3)) return 14; if (!loop3.GetDisplay()) return 15; // ====================================================== gLog << endl; gLog.Separator("Extracting pedestal"); MTaskList tlist4; MParList plist4; plist4.AddToList(&tlist4); plist4.AddToList(&drscalib300); plist4.AddToList(&badpixels); plist4.AddToList(&timecam); MEvtLoop loop4("DetermineExtractedPed"); loop4.SetDisplay(d); loop4.SetParList(&plist4); // ------------------ Setup the tasks --------------- MRawFitsRead read4; read4.LoadMap(pmap); read4.AddFile(pedfile); MContinue cont4("(MRawEvtHeader.GetTriggerID&0xff00)!=0x400", "SelectPed"); MGeomApply apply4; MDrsCalibApply drsapply4; MExtractTimeAndChargeSpline extractor4; extractor4.SetRange(first_slice, last_slice); extractor4.SetRiseTimeHiGain(rise_time_dat); extractor4.SetFallTimeHiGain(fall_time_dat); extractor4.SetHeightTm(heighttm); extractor4.SetChargeType(type); extractor4.SetSaturationLimit(600000); extractor4.SetNoiseCalculation(kFALSE); // MHCamEvent evt3a(0, "PedExt", "Extracted Pedestal Signal;;S"); // MFillH fill3a(&evt3a, "MExtractedSignalCam", "FillPedExt"); // Use this for data, but not for calibration events // evt3a.SetErrorSpread(kFALSE); /* MCalibrateData conv4; conv4.SetCalibrationMode(MCalibrateData::kNone); conv4.SetPedestalFlag(MCalibrateData::kNo); conv4.SetCalibConvMinLimit(0); conv4.SetCalibConvMaxLimit(10000); conv4.SetScaleFactor(scale); */ MCalibrateFact conv4; conv4.SetScale(scale); conv4.SetCalibConst(calib); MBadPixelsTreat treat4; treat4.SetProcessPedestalRun(kFALSE); treat4.SetProcessPedestalEvt(kFALSE); MHCamEvent evt4b(0, "PedExt","Interpolated extracted pedestal;;Signal [~phe]"); //evt4b.SetErrorSpread(kFALSE); MFillH fill4b(&evt4b, "MSignalCam", "FillPedExt"); fill4b.SetDrawOption("gaus"); // ------------------ Setup eventloop and run analysis --------------- tlist4.AddToList(&read4); tlist4.AddToList(&apply4); tlist4.AddToList(&drsapply4); tlist4.AddToList(&cont4); tlist4.AddToList(&extractor4); // tlist4.AddToList(&fill4a); tlist4.AddToList(&conv4); tlist4.AddToList(&treat4); tlist4.AddToList(&fill4b); if (!loop4.Eventloop(max4)) return 15; if (!loop4.GetDisplay()) return 16; // =================================================================== gLog << endl; gLog.Separator("Extracting and calibration data"); MTaskList tlist5; MParList plist5; plist5.AddToList(&tlist5); plist5.AddToList(&drscalib300); plist5.AddToList(&badpixels); plist5.AddToList(&timecam); MEvtLoop loop5("CalibratingData"); loop5.SetDisplay(d); loop5.SetParList(&plist5); // ------------------ Setup the tasks --------------- MRawFitsRead read5a; MReadMarsFile read5b("Events"); read5a.LoadMap(pmap); read5a.AddFile(datfile); read5b.DisableAutoScheme(); read5b.AddFile(datfile); MRead &read5 = datfile.EndsWith(".root") ? static_cast(read5b) : static_cast(read5a); MFillH fill5a(&hrate); MGeomApply apply5; MDrsCalibApply drsapply5; MFDataPhrase filterdat("(MRawEvtHeader.GetTriggerID&0xff00)==0", "SelectDat"); MFDataPhrase filtercal("(MRawEvtHeader.GetTriggerID&0xff00)==0x100", "SelectCal"); MFDataPhrase filterped("(MRawEvtHeader.GetTriggerID&0xff00)==0x400", "SelectPed"); MFDataPhrase filterncl("(MRawEvtHeader.GetTriggerID&0xff00)!=0x100", "SelectNonCal"); //MContinue cont4("MRawEvtHeader.GetTriggerID!=4", "SelectData"); // --- MExtractTimeAndChargeSpline extractor5dat; extractor5dat.SetRange(first_slice, last_slice); extractor5dat.SetRiseTimeHiGain(rise_time_dat); extractor5dat.SetFallTimeHiGain(fall_time_dat); extractor5dat.SetHeightTm(heighttm); extractor5dat.SetChargeType(type); extractor5dat.SetSaturationLimit(600000); extractor5dat.SetNoiseCalculation(kFALSE); MExtractTimeAndChargeSpline extractor5cal; extractor5cal.SetRange(first_slice, last_slice); extractor5cal.SetRiseTimeHiGain(rise_time_cal); extractor5cal.SetFallTimeHiGain(fall_time_cal); extractor5cal.SetHeightTm(heighttm); extractor5cal.SetChargeType(type); extractor5cal.SetSaturationLimit(600000); extractor5cal.SetNoiseCalculation(kFALSE); MExtractTimeAndChargeSpline extractor5tm("ExtractTM"); extractor5tm.SetRange(last_slice, 294); extractor5tm.SetRiseTimeHiGain(1); extractor5tm.SetFallTimeHiGain(1); extractor5tm.SetHeightTm(0.5); extractor5tm.SetChargeType(MExtralgoSpline::kAmplitudeRel); extractor5tm.SetSaturationLimit(600000); extractor5tm.SetNoiseCalculation(kFALSE); extractor5tm.SetNameSignalCam("TimeMarkerAmplitude"); extractor5tm.SetNameTimeCam("TimeMarkerTime"); extractor5dat.SetFilter(&filterncl); extractor5cal.SetFilter(&filtercal); //extractor4tm.SetFilter(&filtercal); // --- /* MCalibrateData conv5; conv5.SetCalibrationMode(MCalibrateData::kNone); conv5.SetPedestalFlag(MCalibrateData::kNo); conv5.SetCalibConvMinLimit(0); conv5.SetCalibConvMaxLimit(10000); conv5.SetScaleFactor(scale); */ MCalibrateFact conv5; conv5.SetScale(scale); conv5.SetCalibConst(calib); MCalibrateDrsTimes calctm5; calctm5.SetNameUncalibrated("UncalibratedTimes"); MCalibrateDrsTimes calctm5tm("CalibrateTimeMarker"); calctm5tm.SetNameArrivalTime("TimeMarkerTime"); calctm5tm.SetNameUncalibrated("UncalTimeMarker"); calctm5tm.SetNameCalibrated("TimeMarker"); calctm5tm.SetTimeMarker(); //calctm4tm.SetFilter(&filtercal); MBadPixelsTreat treat5; treat5.SetProcessPedestalRun(kFALSE); treat5.SetProcessPedestalEvt(kFALSE); MHCamEvent evt5b(0, "ExtSig", "Extracted signal;;S [mV·sl]"); MHCamEvent evt5c(0, "CalSig", "Calibrated and interpolated signal;;S [~phe]"); MHCamEvent evt5d(4, "ExtSigTm", "Extracted time;;T [sl]"); MHCamEvent evt5e(6, "CalSigTm", "Calibrated and interpolated time;;T [ns]"); MFillH fill5b(&evt5b, "MExtractedSignalCam", "FillExtSig"); MFillH fill5c(&evt5c, "MSignalCam", "FillCalSig"); MFillH fill5d(&evt5d, "MArrivalTimeCam", "FillExtTm"); MFillH fill5e(&evt5e, "MSignalCam", "FillCalTm"); fill5c.SetDrawOption("gaus"); fill5d.SetDrawOption("gaus"); fill5e.SetDrawOption("gaus"); /* fill4b.SetFilter(&filterdat); fill4c.SetFilter(&filterdat); fill4d.SetFilter(&filterdat); fill4e.SetFilter(&filterdat); */ //MFSoftwareTrigger swtrig; //MContinue contsw(&swtrig, "FilterSwTrigger", "Software trigger"); //contsw.SetInverted(); // The second rule is for the case reading raw-files! MWriteRootFile write5(outfile, "RECREATE", "Calibrated Data", 2); write5.AddContainer("MRawRunHeader", "RunHeaders"); write5.AddContainer("MGeomCam", "RunHeaders"); write5.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE); write5.AddContainer("MCorsikaRunHeader", "RunHeaders", kFALSE); write5.AddContainer("MMcRunHeader", "RunHeaders", kFALSE); // Common events write5.AddContainer("MCorsikaEvtHeader", "Events", kFALSE); write5.AddContainer("MMcEvt", "Events", kFALSE); write5.AddContainer("IncidentAngle", "Events", kFALSE); write5.AddContainer("MPointingPos", "Events", kFALSE); write5.AddContainer("MSignalCam", "Events"); write5.AddContainer("MTime", "Events", kFALSE); write5.AddContainer("MRawEvtHeader", "Events"); //write.AddContainer("MTriggerPattern", "Events"); // ------------------ Setup histograms and fill tasks ---------------- MContinue test; test.SetFilter(&filterncl); MTaskList tlist5tm; tlist5tm.AddToList(&extractor5tm); tlist5tm.AddToList(&calctm5tm); tlist5tm.SetFilter(&filtercal); MTaskList tlist5dat; tlist5dat.AddToList(&fill5b); tlist5dat.AddToList(&fill5c); tlist5dat.AddToList(&fill5d); tlist5dat.AddToList(&fill5e); tlist5dat.SetFilter(&filterdat); tlist5.AddToList(&read5); tlist5.AddToList(&apply5); tlist5.AddToList(&drsapply5); tlist5.AddToList(&filterncl); //tlist5.AddToList(&test); tlist5.AddToList(&filterdat); tlist5.AddToList(&filtercal); tlist5.AddToList(&filterped); tlist5.AddToList(&fill5a); tlist5.AddToList(&extractor5dat); tlist5.AddToList(&extractor5cal); tlist5.AddToList(&calctm5); tlist5.AddToList(&tlist5tm); tlist5.AddToList(&conv5); tlist5.AddToList(&treat5); tlist5.AddToList(&tlist5dat); tlist5.AddToList(&write5); if (!loop5.Eventloop(max4)) return 18; if (!loop5.GetDisplay()) return 19; d->SetTitle(displaytitle, kFALSE); d->SaveAs(displayfile); return 0; }