Index: trunk/Mars/fact/analysis/callisto.C =================================================================== --- trunk/Mars/fact/analysis/callisto.C (revision 17885) +++ trunk/Mars/fact/analysis/callisto.C (revision 17893) @@ -1,5 +1,5 @@ #include "MLogManip.h" -int callisto(const char *seqfile="seq/2012/01/23/20120123_023.seq", const char *outpath = "output") +int callisto9(const char *seqfile="seq/2012/01/23/20120123_023.seq", const char *outpath = "output", bool use_delays=false) { // ====================================================== @@ -12,182 +12,8 @@ // map file to use (get that from La Palma!) - const char *map = usemap ? "/home/fact/FACT++/FACTmap111030.txt" : NULL; -/* - Bool_t maximum = kTRUE; - - const char *lp_template = maximum ? - "template-lp-extractor-maximum.root" : - "template-lp-extractor-leading-edge.root"; - - const char *pulse_template = "template-pulse.root"; -*/ +// const char *map = usemap ? "/home/fact/FACT++/FACTmap111030.txt" : NULL; + const char *map = usemap ? "/scratch/fact/FACTmap111030.txt" : NULL; + // ------------------------------------------------------ - - //bool use_delays=false; - - int spike_removal=3; - - // The gain as extracted from our dark count spectra - double gain = 258; - - // ------------------------------------------------------ - - // Extraction range in slices. It will always(!) contain the full range - // of integration - const int first_slice = 25; // 10ns - const int last_slice = 225; // 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 int type = maximum ? (MExtralgoSpline::kAmplitudeRel) : (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 max0 = max; // Time marker - //Long_t max1 = max; // Light pulser - //Long_t max2 = 3000; // Calibration ratio - //Long_t max3 = max; // Pedestal Rndm - //Long_t max4 = max; // Pedestal Ext - Long_t max5 = max; // Data - - // ====================================================== - - if (map && gSystem->AccessPathName(map, kFileExists)) - { - gLog << err << "ERROR - Cannot access mapping file '" << map << "'" << endl; - return 1; - } - - // The sequence file which defines the files for the analysis - MSequence seq(seqfile); - if (!seq.IsValid()) - { - gLog << err << "ERROR - Sequence '" << seqfile << "' invalid!" << endl; - return 2; - } - - // -------------------------------------------------------------------------------- - - gLog.Separator("Callisto"); - gLog << all << "Calibrate data of sequence '" << seq.GetFileName() << "'" << endl; - gLog << endl; - - // ------------------------------------------------------ - - ostringstream drsname; - drsname << gSystem->DirName(seqfile) << "/"; - drsname << seq.GetNight().GetNightAsInt() << "_"; - drsname << Form("%03d", seq.GetDrsSequence()) << ".drs.seq"; - - MSequence drs(drsname.str().c_str()); - if (!drs.IsValid()) - { - gLog << err << "ERROR - DRS sequence invalid!" << endl; - return 3; - } - - gLog << all << "DRS sequence file: " << drsname.str() << '\n' << endl; - - TString drsfile = seq.GetFileName(0, MSequence::kRawDrs); - if (drsfile.IsNull()) - { - gLog << err << "No DRS file available in sequence." << endl; - return 4; - } - - TString timfile = drs.GetFileName(0, MSequence::kFitsDat); - TString drs1024 = drs.GetFileName(0, MSequence::kFitsDrs); - TString calfile = seq.GetFileName(0, MSequence::kFitsCal); - //TString pedfile = seq.GetFileName(0, MSequence::kFitsPed); - - gLog << all; - gLog << "DRS calib 300: " << drsfile << '\n'; - gLog << "DRS calib 1024: " << drs1024 << "\n\n"; - - MDrsCalibration drscalib300; - if (!drscalib300.ReadFits(drsfile.Data())) - return 5; - - MDrsCalibration drscalib1024; - if (!drscalib1024.ReadFits(drs1024.Data())) - return 6; - - gLog << all; - gLog << "Time calibration : " << timfile << '\n'; - gLog << "Light Pulser file: " << calfile << '\n' << endl; - //gLog << "Pedestal file: " << pedfile << '\n'; - - // ------------------------------------------------------ - - MDirIter iter; - if (seq.GetRuns(iter, MSequence::kFitsDat)<=0) - { - gLog << err << "ERROR - Sequence valid but without files." << endl; - return 7; - } - iter.Print(); - - // ====================================================== - - 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; @@ -201,5 +27,4 @@ badpixels[1208].SetUnsuitable(MBadPixelsPix::kUnsuitable); badpixels[1399].SetUnsuitable(MBadPixelsPix::kUnsuitable); - // Twin pixel // 113 @@ -210,5 +35,152 @@ // 1393 - MDrsCalibrationTime timecam; + // ------------------------------------------------------ + + // ------------------------------------------------------ + + // 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 = 25; // 10ns + const int last_slice = 225; // 125ns + + const double heighttm = 0.5; // IntegralAbs { 1.5pe * 9.6mV/pe } / IntegralRel { 0.5 } + + Long_t max = 0; // All + Long_t max0 = max; // Time marker + Long_t max1 = max; // Light pulser + //Long_t max2 = 3000; // Calibration ratio + Long_t max3 = max; // Pedestal Rndm + Long_t max4 = max; // Pedestal Ext + Long_t max5 = max; // Data + + // ========================= Result ================================== + + //double scale = 0.1; + double scale = 0.1024; + + // ====================================================== + + if (map && gSystem->AccessPathName(map, kFileExists)) + { + gLog << err << "ERROR - Cannot access mapping file '" << map << "'" << endl; + return 1; + } + + // The sequence file which defines the files for the analysis + MSequence seq(seqfile); + if (!seq.IsValid()) + { + gLog << err << "ERROR - Sequence '" << seqfile << "' invalid!" << endl; + return 2; + } + + // -------------------------------------------------------------------------------- + + gLog.Separator("Callisto"); + gLog << all << "Calibrate data of sequence '" << seq.GetFileName() << "'" << endl; + gLog << endl; + + // ------------------------------------------------------ + + ostringstream drsname; + drsname << gSystem->DirName(seqfile) << "/"; + drsname << seq.GetNight().GetNightAsInt() << "_"; + drsname << Form("%03d", seq.GetDrsSequence()) << ".drs.seq"; + + MSequence drs(drsname.str().c_str()); + if (!drs.IsValid()) + { + gLog << err << "ERROR - DRS sequence invalid!" << endl; + return 3; + } + + gLog << all << "DRS sequence file: " << drsname.str() << '\n' << endl; + + TString drsfile = seq.GetFileName(0, MSequence::kRawDrs); + if (drsfile.IsNull()) + { + gLog << err << "No DRS file available in sequence." << endl; + return 4; + } + + TString timfile = drs.GetFileName(0, MSequence::kFitsDat); + TString drs1024 = drs.GetFileName(0, MSequence::kFitsDrs); + TString pedfile = seq.GetFileName(0, MSequence::kFitsPed); + TString calfile = seq.GetFileName(0, MSequence::kFitsCal); + + gLog << all; + gLog << "DRS calib 300: " << drsfile << '\n'; + gLog << "DRS calib 1024: " << drs1024 << "\n\n"; + + MDrsCalibration drscalib300; + if (!drscalib300.ReadFits(drsfile.Data())) + return 5; + + MDrsCalibration drscalib1024; + if (!drscalib1024.ReadFits(drs1024.Data())) + return 6; + + gLog << all; + gLog << "Time calibration : " << timfile << '\n'; + gLog << "Pedestal file: " << pedfile << '\n'; + gLog << "Light Pulser file: " << calfile << '\n' << endl; + + // ------------------------------------------------------ + + MDirIter iter; + if (seq.GetRuns(iter, MSequence::kFitsDat)<=0) + { + gLog << err << "ERROR - Sequence valid but without files." << endl; + return 7; + } + iter.Print(); + + // ====================================================== + +/* + 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); +*/ + // ====================================================== // Plot the trigger pattern rates vs. run-number @@ -223,9 +195,5 @@ // hrate.DefaultLabelY("ERROR"); - Bool_t isinteg = - (type&MExtralgoSpline::kIntegral) || - (type&MExtralgoSpline::kFixedWidth) || - (type&MExtralgoSpline::kDynWidth) - ? kTRUE : kFALSE; + MDrsCalibrationTime timecam; gStyle->SetOptFit(kTRUE); @@ -241,5 +209,4 @@ plist0.AddToList(&tlist0); plist0.AddToList(&drscalib1024); - plist0.AddToList(&badpixels); plist0.AddToList(&timecam); @@ -257,5 +224,5 @@ MDrsCalibApply drsapply0; - drsapply0.SetRemoveSpikes(spike_removal); + //drsapply0.SetRemoveSpikes(4); MFillH fill0("MHDrsCalibrationTime"); @@ -281,5 +248,5 @@ // ====================================================== -/* + gLog << endl; gLog.Separator("Processing external light pulser run"); @@ -308,35 +275,13 @@ MDrsCalibApply drsapply1; - drsapply1.SetRemoveSpikes(spike_removal); + //drsapply1.SetRemoveSpikes(4); + + MFilterData filterdata1; // --- - MExtractTimeAndChargeSpline extractor1b("ExtractPulse"); + MExtractFACT extractor1b("ExtractPulse"); extractor1b.SetRange(first_slice, last_slice); - extractor1b.SetRiseTimeHiGain(rise_time_cal); - extractor1b.SetFallTimeHiGain(fall_time_cal); - extractor1b.SetHeightTm(heighttm); - extractor1b.SetChargeType(type); - extractor1b.SetSaturationLimit(600000); extractor1b.SetNoiseCalculation(kFALSE); - - MExtractTimeAndChargeSpline extractor1c("ExtractAmplitude"); - extractor1c.SetRange(first_slice, last_slice); - extractor1c.SetChargeType(MExtralgoSpline::kAmplitude); - extractor1c.SetSaturationLimit(600000); - extractor1c.SetNoiseCalculation(kFALSE); - extractor1c.SetNameSignalCam("Amplitude"); - extractor1c.SetNameTimeCam("AmplitudePos"); - - // --- - - MHCamEvent evt1a(5, "CalRatio", "Ratio per slice between integrated signal and amplitude;; r [1/n]"); - evt1a.SetNameSub("Amplitude", kTRUE); - MFillH fill1a(&evt1a, "MExtractedSignalCam", "FillRatio"); - fill1a.SetDrawOption("gaus"); - - MParameterD ratio1a; - ratio1a.SetVal(1./(fall_time_cal+rise_time_cal)); - fill1a.SetWeight(&ratio1a); // --- @@ -391,10 +336,6 @@ tlist1.AddToList(&drsapply1); tlist1.AddToList(&cont1); + tlist1.AddToList(&filterdata1); tlist1.AddToList(&extractor1b); - if (isinteg) - { - tlist1.AddToList(&extractor1c); - tlist1.AddToList(&fill1a); - } tlist1.AddToList(&calctm1a); tlist1.AddToList(&treat1); @@ -416,123 +357,6 @@ timecam.SetDelays(*evt1h.GetHist()); - // ========================= 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++) - { - Double_t g = gain->GetBinContent(i+1)>0.5 ? gain->GetBinContent(i+1) : 1; - if (evt1f.GetHist()->GetBinContent(i+1)>0 && !badpixels[i].IsUnsuitable()) - calib[i] = lpref->GetBinContent(i+1) / evt1f.GetHist()->GetBinContent(i+1) / g; - } - - 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(); -*/ - - TF1 f("f", "[0]*(1-1/(1+exp(x/[1])))*exp(-x/[2])", -15, 150); - f.SetParameter(0, gain*0.05143); - f.SetParameter(1, 1.075); - f.SetParameter(2, 19.3); - f.SetNpx(2*165); //2GHz - - // Convert to graph - TGraph g(&f); - - // Convert to float - MArrayF x(g.GetN()); - MArrayF y(g.GetN()); - for (int i=0; i