/* ======================================================================== *\ ! ! * ! * 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 ! ! \* ======================================================================== */ //const TString pedfile = "/mnt/users/mdoro/Mars/Data/20040201_14418_P_OffMrk421-1_E.root"; //const TString calfile = "/mnt/users/mdoro/Mars/Data/20040201_1441*_C_OffMrk421-1_E.root"; const TString pedfile = "/mnt/Data/rootdata/CrabNebula/2004_02_10/20040210_14607_P_CrabNebula_E.root"; const TString calfile = "/mnt/Data/rootdata/CrabNebula/2004_02_10/20040210_14608_C_CrabNebula_E.root"; //const TString pedfile = "/mnt/Data/rootdata/CrabNebula/2004_01_26/20040125_10412_P_Crab-On_E.root"; //const TString calfile = "/mnt/Data/rootdata/CrabNebula/2004_01_26/20040125_1041*_C_Crab-On_E.root"; //const TString pedfile = "/mnt/Data/rootdata/Miscellaneous/2003_12_19/20031218_03522_P_Park_E.root"; //const TString calfile = "/mnt/Data/rootdata/Miscellaneous/2003_12_19/20031218_03527_C_Park_E.root"; void calibration(TString pedname=pedfile, TString calname=calfile) { gStyle->SetOptStat(1111); gStyle->SetOptFit(); MStatusDisplay *display = new MStatusDisplay; display->SetUpdateTime(3000); display->Resize(850,700); // // Create a empty Parameter List and an empty Task List // The tasklist is identified in the eventloop by its name // MParList plist; MTaskList tlist; plist.AddToList(&tlist); // // Now setup the tasks and tasklist for the pedestals: // --------------------------------------------------- // MReadMarsFile read("Events", pedname); read.DisableAutoScheme(); MGeomApply geomapl; MExtractSignal sigcalc; // // Set the extraction range higher: // //sigcalc.SetRange(1,14,1,14); MPedCalcPedRun pedcalc; // // Additionally to calculating the pedestals, // you can fill histograms and look at them // MFillH fill("MHPedestalCam", "MExtractedSignalCam"); tlist.AddToList(&read); tlist.AddToList(&geomapl); tlist.AddToList(&sigcalc); tlist.AddToList(&pedcalc); tlist.AddToList(&fill); MGeomCamMagic geomcam; MPedestalCam pedcam; MHPedestalCam hpedcam; plist.AddToList(&geomcam); plist.AddToList(&pedcam); plist.AddToList(&hpedcam); // // Create and setup the eventloop // MEvtLoop evtloop; evtloop.SetParList(&plist); evtloop.SetDisplay(display); // // Execute first analysis // if (!evtloop.Eventloop()) return; tlist.PrintStatistics(); // // Look at one specific pixel, after all the histogram manipulations: // hpedcam[9].DrawClone("fourier"); MHCamera dispped0 (geomcam, "Ped;Pedestal", "Mean per Slice"); MHCamera dispped1 (geomcam, "Ped;PedestalErr", "Mean Error per Slice"); MHCamera dispped2 (geomcam, "Ped;PedestalRms", "RMS per Slice"); MHCamera dispped3 (geomcam, "Ped;PedestalRmsErr", "RMS Error per Slice"); MHCamera dispped4 (geomcam, "Ped;Mean", "Fitted Mean per Slice"); MHCamera dispped5 (geomcam, "Ped;MeanErr", "Fitted Error of Mean per Slice"); MHCamera dispped6 (geomcam, "Ped;Sigma", "Fitted Sigma per Slice"); MHCamera dispped7 (geomcam, "Ped;SigmaErr", "Fitted Error of Sigma per Slice"); MHCamera dispped8 (geomcam, "Ped;Prob", "Probability of Fit"); MHCamera dispped9 (geomcam, "Ped;DeltaPedestalMean", "Rel. Diff. Mean per Slice (Calc.-Fitte)"); MHCamera dispped10 (geomcam, "Ped;DeltaPedestalMeanError", "Rel. Diff. Mean Error per Slice (Calc.-Fitted)"); MHCamera dispped11 (geomcam, "Ped;DeltaRmsSigma", "Rel. Diff. RMS per Slice (Calc.-Fitted)"); MHCamera dispped12 (geomcam, "Ped;DeltaRmsSigmaError", "Rel. Diff. RMS Error per Slice (Calc.-Fitted)"); MHCamera dispped13 (geomcam, "Ped;FitOK", "Gaus Fit not OK"); MHCamera dispped14 (geomcam, "Ped;FourierOK", "Fourier Analysis not OK"); dispped0.SetCamContent( pedcam, 0); dispped0.SetCamError( pedcam, 1); dispped1.SetCamContent( pedcam, 1); dispped2.SetCamContent( pedcam, 2); dispped2.SetCamError( pedcam, 3); dispped3.SetCamContent( pedcam, 3); dispped4.SetCamContent( hpedcam, 0); dispped4.SetCamError( hpedcam, 1); dispped5.SetCamContent( hpedcam, 1); dispped6.SetCamContent( hpedcam, 2); dispped6.SetCamError( hpedcam, 3); dispped7.SetCamContent( hpedcam, 3); dispped8.SetCamContent( hpedcam, 4); dispped9.SetCamContent( hpedcam, 5); dispped9.SetCamError( hpedcam, 6); dispped10.SetCamContent(hpedcam, 7); dispped11.SetCamContent(hpedcam, 8); dispped11.SetCamError( hpedcam, 9); dispped12.SetCamContent(hpedcam, 10); dispped13.SetCamContent(hpedcam, 11); dispped14.SetCamContent(hpedcam, 12); dispped0.SetYTitle("Calc. Pedestal per slice [FADC counts]"); dispped1.SetYTitle("Calc. Pedestal Error per slice [FADC counts]"); dispped2.SetYTitle("Calc. Pedestal RMS per slice [FADC counts]"); dispped3.SetYTitle("Calc. Pedestal RMS Error per slice [FADC counts]"); dispped4.SetYTitle("Fitted Mean per slice [FADC counts]"); dispped5.SetYTitle("Error of Fitted Mean per slice [FADC counts]"); dispped6.SetYTitle("Fitted Sigma per slice [FADC counts]"); dispped7.SetYTitle("Error of Fitted Sigma per slice [FADC counts]"); dispped8.SetYTitle("Fit Probability [1]"); dispped9.SetYTitle("Rel. Diff. Pedestal Calc.-Fitted per slice [1]"); dispped10.SetYTitle("Rel. Diff. Pedestal Error Calc.-Fitted per slice [1]"); dispped11.SetYTitle("Rel. Diff. Pedestal RMS Calc.-Fitted per slice [1]"); dispped12.SetYTitle("Rel. Diff. Pedestal RMS Error Calc.-Fitted per slice [1]"); dispped13.SetYTitle("[1]"); dispped14.SetYTitle("[1]"); // Histogram values TCanvas &b1 = display->AddTab("Ped.Calc."); b1.Divide(4,3); CamDraw(b1,dispped0,pedcam,1,4,1); CamDraw(b1,dispped1,pedcam,2,4,2); CamDraw(b1,dispped2,pedcam,3,4,2); CamDraw(b1,dispped3,pedcam,4,4,2); // Fitted values TCanvas &b2 = display->AddTab("Ped.Fit"); b2.Divide(4,3); CamDraw(b2,dispped4,hpedcam,1,4,1); CamDraw(b2,dispped5,hpedcam,2,4,2); CamDraw(b2,dispped6,hpedcam,3,4,2); CamDraw(b2,dispped7,hpedcam,4,4,2); // Fits Probability TCanvas &b3 = display->AddTab("Ped.Fit Prob."); b3.Divide(1,3); CamDraw(b3,dispped8,hpedcam,1,1,3); // Differences TCanvas &c4 = display->AddTab("Rel.Diff.Calc.-Fit"); c4.Divide(4,3); CamDraw(c4,dispped9,hpedcam,1,4,1); CamDraw(c4,dispped10,hpedcam,2,4,1); CamDraw(c4,dispped11,hpedcam,3,4,1); CamDraw(c4,dispped12,hpedcam,4,4,1); // Defects TCanvas &c5 = display->AddTab("Defects"); c5.Divide(2,2); CamDraw(c5,dispped13,hpedcam,1,2,0); CamDraw(c5,dispped14,hpedcam,2,2,0); // // Create a empty Parameter List and an empty Task List // MParList plist2; MTaskList tlist2; plist2.AddToList(&tlist2); MExtractedSignalCam sigcam; MCalibrationCam calcam; MHCalibrationRelTimeCam timecam; // // Get the previously created MPedestalCam into the new Parameter List // plist2.AddToList(&geomcam); plist2.AddToList(&pedcam); plist2.AddToList(&sigcam); plist2.AddToList(&calcam); plist2.AddToList(&timecam); // // Get the MAGIC geometry // tlist2.AddToList(&geomapl); // // Now setup the new tasks and tasklist for the calibration // --------------------------------------------------- // MReadMarsFile read2("Events", calname); read2.DisableAutoScheme(); // // We saw that the signal jumps between slices, // thus take the sliding window // MExtractPINDiode pincalc; MExtractBlindPixel blindcalc; MExtractSignal2 sigcalc2; MArrivalTimeCalc2 timecalc; MCalibrationCalc calcalc; MFillH filltime("MHCalibrationRelTimeCam", "MArrivalTime"); MFillH fillpin("MHCalibrationChargePINDiode", "MExtractedSignalPINDiode"); // // Set the range (other than default) // of FADC slices for the blind pixel // // calcalc.SetBlindPixelRange(10,25); // // Set the cut upon which a superposition of the blind pixel // FADC slices will be filled into the SinglePHE histogram // // calcalc.SetBlindPixelSinglePheCut(500); // // Skip the HiGain vs. LoGain calibration // calcalc.SkipHiLoGainCalibration(); // // As long, as we don't have digital modules, // we have to set the color of the pulser LED by hand // calcalc.SetPulserColor(MCalibrationCalc::kECT1); // // In case, we want to exclude a pre-defined list of bad pixels: // (This is a preliminary feature) // // calcalc.ExcludePixelsFromAsciiFile("badpixels.dat"); // // In case, you want to skip the quality checks // (NOT RECOMMENDED!!!) // // calcalc.SkipQualityChecks(); // // In case, we want to apply another fit function to the // blind pixel // // MCalibrationBlindPix *bp = calcam.GetBlindPixel(); // bp->ChangeFitFunc(MHCalibrationBlindPixel::kEPoisson5); // // Apply a filter against cosmics // (was directly in MCalibrationCalc in earlier versions) // MFCosmics cosmics; MContinue cont(&cosmics); tlist2.AddToList(&read2); tlist2.AddToList(&blindcalc); tlist2.AddToList(&pincalc); tlist2.AddToList(&sigcalc2); // // 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 // tlist2.AddToList(&timecalc); tlist2.AddToList(&filltime); tlist2.AddToList(&fillpin); // tlist2.AddToList(&calcalc); // // Create and setup the eventloop // MEvtLoop evtloop2; evtloop2.SetParList(&plist2); evtloop2.SetDisplay(display); // // Execute second analysis // if (!evtloop2.Eventloop()) return; tlist2.PrintStatistics(); // // print the most important results of all pixels // // calcam.Print(); // // just one example how to get the plots of individual pixels // // calcam[563].DrawClone(); // calcam[564].DrawClone(); // Create histograms to display MHCamera disp1 (geomcam, "Cal;Charge", "Fitted Mean Charges"); MHCamera disp2 (geomcam, "Cal;SigmaCharge", "Sigma of Fitted Charges"); MHCamera disp3 (geomcam, "Cal;FitProb", "Probability of Fit"); MHCamera disp4 (geomcam, "Cal;RSigma", "Reduced Sigmas"); MHCamera disp5 (geomcam, "Cal;RSigma/Charge", "Reduced Sigma per Charge"); MHCamera disp6 (geomcam, "Cal;FFactorPhe", "Nr. of Phe's (F-Factor Method)"); MHCamera disp7 (geomcam, "Cal;FFactorConv", "Conversion Factor (F-Factor Method)"); MHCamera disp8 (geomcam, "Cal;FFactorFFactor", "Total F-Factor (F-Factor Method)"); MHCamera disp9 (geomcam, "Cal;BlindPixPh", "Photon flux inside plexiglass (Blind Pixel Method)"); MHCamera disp10 (geomcam, "Cal;BlindPixConv", "Conversion Factor (Blind Pixel Method)"); MHCamera disp11 (geomcam, "Cal;BlindPixFFactor","Total F-Factor (Blind Pixel Method)"); MHCamera disp12 (geomcam, "Cal;PINDiodePh", "Photon flux outside plexiglass (PIN Diode Method)"); MHCamera disp13 (geomcam, "Cal;PINDiodeConv", "Conversion Factor (PIN Diode Method)"); MHCamera disp14 (geomcam, "Cal;PINDiodeFFactor","Total F-Factor (PIN Diode Method)"); MHCamera disp15 (geomcam, "Cal;Excluded", "Pixels previously excluded"); MHCamera disp16 (geomcam, "Cal;NotFitted", "Pixels that could not be fitted"); MHCamera disp17 (geomcam, "Cal;NotFitValid", "Pixels with not valid fit results"); MHCamera disp18 (geomcam, "Cal;Oscillating", "Oscillating Pixels"); MHCamera disp19 (geomcam, "Cal;Saturation", "Pixels with saturated Hi Gain"); MHCamera disp20 (geomcam, "Cal;Ped", "Pedestals"); MHCamera disp21 (geomcam, "Cal;PedRms", "Pedestal RMS"); MHCamera disp22 (geomcam, "time;Time", "Rel. Arrival Times"); MHCamera disp23 (geomcam, "time;SigmaTime", "Sigma of Rel. Arrival Times"); MHCamera disp24 (geomcam, "time;TimeProb", "Probability of Time Fit"); MHCamera disp25 (geomcam, "time;NotFitValid", "Pixels with not valid fit results"); MHCamera disp26 (geomcam, "time;Oscillating", "Oscillating Pixels"); MHCamera disp27 (geomcam, "Cal;AbsTimeMean", "Abs. Arrival Times"); MHCamera disp28 (geomcam, "Cal;AbsTimeRms", "RMS of Arrival Times"); // Fitted charge means and sigmas disp1.SetCamContent(calcam, 0); disp1.SetCamError( calcam, 1); disp2.SetCamContent(calcam, 2); disp2.SetCamError( calcam, 3); // Fit probabilities disp3.SetCamContent(calcam, 4); // Reduced Sigmas and reduced sigmas per charge disp4.SetCamContent(calcam, 5); disp4.SetCamError( calcam, 6); disp5.SetCamContent(calcam, 7); disp5.SetCamError( calcam, 8); // F-Factor Method disp6.SetCamContent(calcam, 9); disp6.SetCamError( calcam, 10); disp7.SetCamContent(calcam, 11); disp7.SetCamError( calcam, 12); disp8.SetCamContent(calcam, 13); disp8.SetCamError( calcam, 14); /// Blind Pixel Method disp9.SetCamContent(calcam, 15); disp9.SetCamError( calcam, 16); disp10.SetCamContent(calcam,17); disp10.SetCamError( calcam,18); disp11.SetCamContent(calcam,19); disp11.SetCamError( calcam,20); // PIN Diode Method disp12.SetCamContent(calcam,21); disp12.SetCamError( calcam,22); disp13.SetCamContent(calcam,23); disp13.SetCamError( calcam,24); disp14.SetCamContent(calcam,25); disp14.SetCamError( calcam,26); // Pixels with defects disp15.SetCamContent(calcam,27); disp16.SetCamContent(calcam,28); disp17.SetCamContent(calcam,29); disp18.SetCamContent(calcam,30); // Lo Gain calibration disp19.SetCamContent(calcam,31); // Pedestals disp20.SetCamContent(calcam,35); disp20.SetCamError( calcam,36); disp21.SetCamContent(calcam,37); disp21.SetCamError( calcam,38); // Relative Times disp22.SetCamContent(timecam,0); disp22.SetCamError( timecam,1); disp23.SetCamContent(timecam,2); disp23.SetCamError( timecam,3); disp24.SetCamContent(timecam,4); disp25.SetCamContent(timecam,5); disp26.SetCamContent(timecam,6); // Absolute Times disp27.SetCamContent(calcam,39); disp27.SetCamError( calcam,40); disp28.SetCamContent(calcam,41); disp28.SetCamError( calcam,42); disp1.SetYTitle("Charge [FADC units]"); disp2.SetYTitle("\\sigma_{Charge} [FADC units]"); disp3.SetYTitle("P_{Charge} [1]"); disp4.SetYTitle("\\sqrt{\\sigma^{2}_{Charge} - RMS^{2}_{Ped}} [FADC Counts]"); disp5.SetYTitle("Reduced Sigma / Mean Charge [1]"); disp6.SetYTitle("Nr. Photo-Electrons [1]"); disp7.SetYTitle("Conversion Factor [PhE/FADC Count]"); disp8.SetYTitle("\\sqrt{N_{PhE}}*\\sigma_{Charge}/\\mu_{Charge} [1]"); disp9.SetYTitle("Photon flux [ph/mm^2]"); disp10.SetYTitle("Conversion Factor [Phot/FADC Count]"); disp11.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]"); disp12.SetYTitle("Photon flux [ph/mm^2]"); disp13.SetYTitle("Conversion Factor [Phot/FADC Count]"); disp14.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]"); disp15.SetYTitle("[1]"); disp16.SetYTitle("[1]"); disp17.SetYTitle("[1]"); disp18.SetYTitle("[1]"); disp19.SetYTitle("[1]"); disp20.SetYTitle("Ped [FADC Counts ]"); disp21.SetYTitle("RMS_{Ped} [FADC Counts ]"); disp22.SetYTitle("Time Offset [ns]"); disp23.SetYTitle("Timing resolution [ns]"); disp24.SetYTitle("P_{Time} [1]"); disp25.SetYTitle("[1]"); disp26.SetYTitle("[1]"); disp27.SetYTitle("Mean Abs. Time [FADC slice]"); disp28.SetYTitle("RMS Abs. Time [FADC slices]"); gStyle->SetOptStat(1111); gStyle->SetOptFit(); // Charges TCanvas &c1 = display->AddTab("Fit.Charge"); c1.Divide(2, 3); CamDraw(c1, disp1,calcam,1, 2 , 2); CamDraw(c1, disp2,calcam,2, 2 , 2); // Fit Probability TCanvas &c2 = display->AddTab("Fit.Prob"); c2.Divide(1,3); CamDraw(c2, disp3,calcam,1, 1 , 4); // Reduced Sigmas TCanvas &c3 = display->AddTab("Red.Sigma"); c3.Divide(2,3); CamDraw(c3, disp4,calcam,1, 2 , 2); CamDraw(c3, disp5,calcam,2, 2 , 2); // F-Factor Method TCanvas &c4 = display->AddTab("F-Factor"); c4.Divide(3,3); CamDraw(c4, disp6,calcam,1, 3 , 2); CamDraw(c4, disp7,calcam,2, 3 , 2); CamDraw(c4, disp8,calcam,3, 3 , 2); // Blind Pixel Method TCanvas &c5 = display->AddTab("BlindPix"); c5.Divide(3, 3); CamDraw(c5, disp9,calcam,1, 3 , 9); CamDraw(c5, disp10,calcam,2, 3 , 2); CamDraw(c5, disp11,calcam,3, 3 , 2); // PIN Diode Method TCanvas &c6 = display->AddTab("PINDiode"); c6.Divide(3,3); CamDraw(c6, disp12,calcam,1, 3 , 9); CamDraw(c6, disp13,calcam,2, 3 , 2); CamDraw(c6, disp14,calcam,3, 3 , 2); // Defects TCanvas &c7 = display->AddTab("Defects"); c7.Divide(4,2); CamDraw(c7, disp15,calcam,1,4, 0); CamDraw(c7, disp16,calcam,2,4, 0); CamDraw(c7, disp17,calcam,3,4, 0); CamDraw(c7, disp18,calcam,4,4,0); // Lo Gain Calibration TCanvas &c8 = display->AddTab("LowGain"); c8.Divide(1,3); CamDraw(c8, disp19,calcam,1,4,0); // Pedestals TCanvas &c9 = display->AddTab("Pedestals"); c9.Divide(2,3); CamDraw(c9,disp20,calcam,1,2,1); CamDraw(c9,disp21,calcam,2,2,2); // Rel. Times TCanvas &c10 = display->AddTab("Fitted Rel. Times"); c10.Divide(3,3); CamDraw(c10,disp22,calcam,1,3,2); CamDraw(c10,disp23,calcam,2,3,2); CamDraw(c10,disp24,calcam,3,3,4); // Time Defects TCanvas &c11 = display->AddTab("Time Def."); c11.Divide(2,2); CamDraw(c11, disp25,calcam,1,2, 0); CamDraw(c11, disp26,calcam,2,2, 0); // Abs. Times TCanvas &c12 = display->AddTab("Abs. Times"); c12.Divide(2,3); CamDraw(c12,disp27,calcam,1,2,2); CamDraw(c12,disp28,calcam,2,2,2); } void CamDraw(TCanvas &c, MHCamera &cam, MCamEvent &evt, Int_t i, Int_t j, Int_t fit) { c.cd(i); gPad->SetBorderMode(0); MHCamera *obj1=(MHCamera*)cam.DrawCopy("hist"); // obj1->AddNotify(evt); c.cd(i+j); gPad->SetBorderMode(0); obj1->Draw(); ((MHCamera*)obj1)->SetPrettyPalette(); if (fit != 0) { c.cd(i+2*j); gPad->SetBorderMode(0); TH1D *obj2 = (TH1D*)obj1->Projection(); obj2->Draw(); obj2->SetBit(kCanDelete); const Double_t min = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst()); const Double_t max = obj2->GetBinCenter(obj2->GetXaxis()->GetLast()); const Double_t integ = obj2->Integral("width")/2.5066283; const Double_t mean = obj2->GetMean(); const Double_t rms = obj2->GetRMS(); const Double_t width = max-min; if (rms == 0. || width == 0. ) return; switch (fit) { case 1: TF1 *sgaus = new TF1("sgaus","gaus(0)",min,max); sgaus->SetBit(kCanDelete); sgaus->SetParNames("Area","#mu","#sigma"); sgaus->SetParameters(integ/rms,mean,rms); sgaus->SetParLimits(0,0.,integ); sgaus->SetParLimits(1,min,max); sgaus->SetParLimits(2,0,width/1.5); obj2->Fit("sgaus","QLR"); obj2->GetFunction("sgaus")->SetLineColor(kYellow); break; case 2: TString dgausform = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"; dgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])"; TF1 *dgaus = new TF1("dgaus",dgausform.Data(),min,max); dgaus->SetBit(kCanDelete); dgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}","A_{2}","#mu_{2}","#sigma_{2}"); dgaus->SetParameters(integ,(min+mean)/2.,width/4., integ/width/2.,(max+mean)/2.,width/4.); // The left-sided Gauss dgaus->SetParLimits(0,integ-1.5,integ+1.5); dgaus->SetParLimits(1,min+(width/10.),mean); dgaus->SetParLimits(2,0,width/2.); // The right-sided Gauss dgaus->SetParLimits(3,0,integ); dgaus->SetParLimits(4,mean,max-(width/10.)); dgaus->SetParLimits(5,0,width/2.); obj2->Fit("dgaus","QLRM"); obj2->GetFunction("dgaus")->SetLineColor(kYellow); break; case 3: TString tgausform = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"; tgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])"; tgausform += "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])"; TF1 *tgaus = new TF1("tgaus",tgausform.Data(),min,max); tgaus->SetBit(kCanDelete); tgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}", "A_{2}","#mu_{2}","#sigma_{2}", "A_{3}","#mu_{3}","#sigma_{3}"); tgaus->SetParameters(integ,(min+mean)/2,width/4., integ/width/3.,(max+mean)/2.,width/4., integ/width/3.,mean,width/2.); // The left-sided Gauss tgaus->SetParLimits(0,integ-1.5,integ+1.5); tgaus->SetParLimits(1,min+(width/10.),mean); tgaus->SetParLimits(2,width/15.,width/2.); // The right-sided Gauss tgaus->SetParLimits(3,0.,integ); tgaus->SetParLimits(4,mean,max-(width/10.)); tgaus->SetParLimits(5,width/15.,width/2.); // The Gauss describing the outliers tgaus->SetParLimits(6,0.,integ); tgaus->SetParLimits(7,min,max); tgaus->SetParLimits(8,width/4.,width/1.5); obj2->Fit("tgaus","QLRM"); obj2->GetFunction("tgaus")->SetLineColor(kYellow); break; case 4: obj2->Fit("pol0","Q"); obj2->GetFunction("pol0")->SetLineColor(kYellow); break; case 9: break; default: obj2->Fit("gaus","Q"); obj2->GetFunction("gaus")->SetLineColor(kYellow); break; } gPad->Modified(); gPad->Update(); } }