Index: /trunk/MagicSoft/Mars/macros/calibrate_data.C =================================================================== --- /trunk/MagicSoft/Mars/macros/calibrate_data.C (revision 3489) +++ /trunk/MagicSoft/Mars/macros/calibrate_data.C (revision 3489) @@ -0,0 +1,737 @@ +/* ======================================================================== *\ +! +! * +! * 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, 03/2004 +! Markus Gaug, 03/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ +#include "MAGIC.h" + +const TString defpath = "/mnt/Data/rootdata/CrabNebula/2004_01_27/"; +const TString defrout = "output_test.root"; + +const Int_t defpedr [] = {12461}; +const Int_t defcalr [] = {12526,12527,12528,12529}; +const Int_t defdatar[] = {12517,12518}; + +void calibrate_data(const TString inpath=defpath, + const Int_t psize=1, const Int_t pedruns[]=defpedr, + const Int_t csize=4, const Int_t calruns[]=defcalr, + const Int_t dsize=2, const Int_t dataruns[]=defdatar, + const TString resname=defrout) + +{ + + MRunIter pruns; + MRunIter cruns; + MRunIter druns; + + for (Int_t i=0;iSetUpdateTime(3000); + display->Resize(850,700); + + gStyle->SetOptStat(1111); + gStyle->SetOptFit(); + + /************************************/ + /* FIRST LOOP: PEDESTAL COMPUTATION */ + /************************************/ + + MParList plist1; + MTaskList tlist1; + plist1.AddToList(&tlist1); + + // containers + MPedestalCam pedcam; + MBadPixelsCam badcam; + + // + // for excluding pixels from the beginning: + // + // badcam1.AsciiRead("badpixels.dat"); + + plist1.AddToList(&pedcam); + plist1.AddToList(&badcam); + + //tasks + MReadMarsFile read("Events"); + read.DisableAutoScheme(); + static_cast(read).AddFiles(pruns); + + MGeomApply geomapl; + MPedCalcPedRun pedcalc; + MGeomCamMagic geomcam; + + tlist1.AddToList(&read); + tlist1.AddToList(&geomapl); + tlist1.AddToList(&pedcalc); + + // Create and execute the event looper + MEvtLoop pedloop; + pedloop.SetParList(&plist1); + pedloop.SetDisplay(display); + + cout << "*************************" << endl; + cout << "** COMPUTING PEDESTALS **" << endl; + cout << "*************************" << endl; + + if (!pedloop.Eventloop()) + return; + + tlist1.PrintStatistics(); + + // + // Now the short version: + // +/* + MJCalibration calloop; + calloop.SetInput(&cruns); + calloop.SetDisplay(display); + if (!calloop.Process(pedloop.GetPedestalCam())) + return; +#if 0 +*/ + // + // The longer version: + // + + // + // Create a empty Parameter List and an empty Task List + // + MParList plist2; + MTaskList tlist2; + plist2.AddToList(&tlist2); + plist2.AddToList(&pedcam); + plist2.AddToList(&badcam); + + gLog << endl;; + gLog << "Calculate MCalibrationCam from Runs " << cruns.GetRunsAsString() << endl; + gLog << endl; + + MReadMarsFile read2("Events"); + read2.DisableAutoScheme(); + static_cast(read2).AddFiles(cruns); + + MGeomCamMagic geomcam; + MExtractedSignalCam sigcam; + MArrivalTimeCam timecam; + MCalibrationChargeCam calcam; + // MCalibrationChargePINDiode pindiode; + // MCalibrationChargeBlindPix blindpix; + + MHCalibrationRelTimeCam histtime; + MHCalibrationChargeCam histcharge; + // MHCalibrationChargePINDiode histpin; + // MHCalibrationChargeBlindPix histblind; + // + // + // Get the previously created MPedestalCam into the new Parameter List + // + plist2.AddToList(&geomcam); + plist2.AddToList(&sigcam); + plist2.AddToList(&timecam); + plist2.AddToList(&calcam); + plist2.AddToList(&histtime); + plist2.AddToList(&histcharge); +// plist2.AddToList(&histpin); +// plist2.AddToList(&histblind); + + // + // We saw that the signal jumps between slices, + // thus take the sliding window + // + MExtractSignal2 sigcalc; + MExtractPINDiode pincalc; + MExtractBlindPixel blindcalc; + + MArrivalTimeCalc2 timecalc; + MCalibrationChargeCalc calcalc; + MGeomApply geomapl; + + MFillH filltime( "MHCalibrationRelTimeCam" , "MArrivalTimeCam"); +// MFillH fillpin ("MHCalibrationChargePINDiode", "MExtractedSignalPINDiode"); +// MFillH fillblind("MHCalibrationChargeBlindPix", "MExtractedSignalBlindPixel"); + MFillH fillcam ("MHCalibrationChargeCam" , "MExtractedSignalCam"); + + // + // Skip the HiGain vs. LoGain calibration + // + calcalc.SkipHiLoGainCalibration(); + + // + // Apply a filter against cosmics + // (was directly in MCalibrationCalc in earlier versions) + // + MFCosmics cosmics; + MContinue cont(&cosmics); + + tlist2.AddToList(&read2); + tlist2.AddToList(&geomapl); + tlist2.AddToList(&sigcalc); + // tlist2.AddToList(&blindcalc); + // tlist2.AddToList(&pincalc); + // + // 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(&fillblind); + tlist2.AddToList(&fillcam); + // + tlist2.AddToList(&calcalc); + // + // 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.Eventloop()) + return; + + tlist2.PrintStatistics(); + + // + // print the most important results of all pixels to a file + // + // MLog gauglog; + // gauglog.SetOutputFile(Form("%s%s",calcam.GetName(),".txt"),1); + // calcam.SetLogStream(&gauglog); + // calcam.Print(); + // calcam.SetLogStream(&gLog); + // + // just one example how to get the plots of individual pixels + // + // histblind.DrawClone("all"); + // histcharge[5].DrawClone("all"); + // histcharge(5).DrawClone("all"); + // histtime[5].DrawClone("fourierevents"); + + // 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;FFactorPh", "Nr. of Photo-electrons (F-Factor Method)"); + MHCamera disp7 (geomcam, "Cal;FFactorConv", "Conversion Factor to photons (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 to photons (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 tp photons (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;HiGainOscillating", "Oscillating Pixels HI Gain"); + MHCamera disp19 (geomcam, "Cal;LoGainOscillating", "Oscillating Pixels LO Gain"); + MHCamera disp20 (geomcam, "Cal;HiGainPickup", "Number Pickup events Hi Gain"); + MHCamera disp21 (geomcam, "Cal;LoGainPickup", "Number Pickup events Lo Gain"); + MHCamera disp22 (geomcam, "Cal;Saturation", "Pixels with saturated Hi Gain"); + MHCamera disp23 (geomcam, "Cal;FFactorValid", "Pixels with valid F-Factor calibration"); + MHCamera disp24 (geomcam, "Cal;BlindPixelValid", "Pixels with valid BlindPixel calibration"); + MHCamera disp25 (geomcam, "Cal;PINdiodeFFactorValid", "Pixels with valid PINDiode calibration"); + + MHCamera disp26 (geomcam, "Cal;Ped", "Pedestals"); + MHCamera disp27 (geomcam, "Cal;PedRms", "Pedestal RMS"); + + MHCamera disp28 (geomcam, "time;Time", "Rel. Arrival Times"); + MHCamera disp29 (geomcam, "time;SigmaTime", "Sigma of Rel. Arrival Times"); + MHCamera disp30 (geomcam, "time;TimeProb", "Probability of Time Fit"); + MHCamera disp31 (geomcam, "time;NotFitValid", "Pixels with not valid fit results"); + MHCamera disp32 (geomcam, "time;Oscillating", "Oscillating Pixels"); + + MHCamera disp33 (geomcam, "Cal;AbsTimeMean", "Abs. Arrival Times"); + MHCamera disp34 (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(badcam,9); + disp18.SetCamContent(badcam,16); + disp19.SetCamContent(badcam,15); + disp20.SetCamContent(calcam,29); + disp21.SetCamContent(calcam,30); + + // Lo Gain calibration + disp22.SetCamContent(calcam,31); + + // Valid flags + disp23.SetCamContent(calcam,32); + disp24.SetCamContent(calcam,33); + disp25.SetCamContent(calcam,34); + + // Pedestals + disp26.SetCamContent(calcam,35); + disp26.SetCamError( calcam,36); + disp27.SetCamContent(calcam,37); + disp27.SetCamError( calcam,38); + + + // Relative Times + disp28.SetCamContent(histtime,0); + disp28.SetCamError( histtime,1); + disp29.SetCamContent(histtime,2); + disp29.SetCamError( histtime,3); + disp30.SetCamContent(histtime,4); + disp31.SetCamContent(histtime,5); + disp32.SetCamContent(histtime,6); + + // Absolute Times + disp33.SetCamContent(calcam,39); + disp33.SetCamError( calcam,40); + disp34.SetCamContent(calcam,41); + + 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 [Ph/FADC Count]"); + disp8.SetYTitle("\\sqrt{N_{Ph}}*\\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("[1]"); + disp21.SetYTitle("[1]"); + disp22.SetYTitle("[1]"); + disp23.SetYTitle("[1]"); + disp24.SetYTitle("[1]"); + disp25.SetYTitle("[1]"); + + disp26.SetYTitle("Ped [FADC Counts ]"); + disp27.SetYTitle("RMS_{Ped} [FADC Counts ]"); + + disp28.SetYTitle("Time Offset [ns]"); + disp29.SetYTitle("Timing resolution [ns]"); + disp30.SetYTitle("P_{Time} [1]"); + + disp31.SetYTitle("[1]"); + disp32.SetYTitle("[1]"); + + disp33.SetYTitle("Mean Abs. Time [FADC slice]"); + disp34.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, disp20,calcam,3,4, 0); + CamDraw(c7, disp21,calcam,4,4, 0); + + // BadCam + TCanvas &c8 = display->AddTab("Defects"); + c8.Divide(3,2); + + CamDraw(c8, disp17,badcam,1,3, 0); + CamDraw(c8, disp18,badcam,2,3, 0); + CamDraw(c8, disp19,badcam,3,3, 0); + + + // Valid flags + TCanvas &c9 = display->AddTab("Validity"); + c9.Divide(4,2); + + CamDraw(c9, disp22,calcam,1,4,0); + CamDraw(c9, disp23,calcam,2,4,0); + CamDraw(c9, disp24,calcam,3,4,0); + CamDraw(c9, disp25,calcam,4,4,0); + + + // Pedestals + TCanvas &c10 = display->AddTab("Pedestals"); + c10.Divide(2,3); + + CamDraw(c10,disp26,calcam,1,2,1); + CamDraw(c10,disp27,calcam,2,2,2); + + // Rel. Times + TCanvas &c11 = display->AddTab("Fitted Rel. Times"); + c11.Divide(3,3); + + CamDraw(c11,disp28,calcam,1,3,2); + CamDraw(c11,disp29,calcam,2,3,2); + CamDraw(c11,disp30,calcam,3,3,4); + + // Time Defects + TCanvas &c12 = display->AddTab("Time Def."); + c12.Divide(2,2); + + CamDraw(c12, disp31,calcam,1,2, 0); + CamDraw(c12, disp32,calcam,2,2, 0); + + // Abs. Times + TCanvas &c13 = display->AddTab("Abs. Times"); + c13.Divide(2,3); + + CamDraw(c13,disp33,calcam,1,2,2); + CamDraw(c13,disp34,calcam,2,2,2); + + /************************************************************************/ + /* THIRD LOOP: DATA CALIBRATION INTO PHOTONS */ + /************************************************************************/ + + // Create an empty Parameter List and an empty Task List + MParList plist3; + MTaskList tlist3; + plist3.AddToList(&tlist3); + + // containers + MCerPhotEvt photevt; + MPedPhotCam pedphotcam; + MSrcPosCam srccam; + MRawRunHeader runhead; + + plist3.AddToList(&geomcam ); + plist3.AddToList(&pedcam ); + plist3.AddToList(&calcam ); + plist3.AddToList(&badcam ); + plist3.AddToList(&timecam ); + plist3.AddToList(&sigcam ); + plist3.AddToList(&histtime); + plist3.AddToList(&photevt); + plist3.AddToList(&pedphotcam); + plist3.AddToList(&srccam); + plist3.AddToList(&runhead); + + //tasks + MReadMarsFile read3("Events"); + read3.DisableAutoScheme(); + static_cast(read3).AddFiles(druns); + + MCalibrateData photcalc; + photcalc.SetCalibrationMode(MCalibrateData::kFfactor); // !!! was only MCalibrate + // MPedPhotCalc pedphotcalc; // already done by MCalibrate Data + // MCerPhotCalc cerphotcalc; // already done by MCalibrate Data + + tlist3.AddToList(&read3); + tlist3.AddToList(&geomapl); + tlist3.AddToList(&sigcalc); + tlist3.AddToList(&timecalc); + // tlist3.AddToList(&cerphotcalc); // already done by MCalibrate Data + tlist3.AddToList(&photcalc); + // tlist3.AddToList(&pedphotcalc); // already done by MCalibrate Data + + MWriteRootFile write(resname); + + write.AddContainer("MGeomCam" , "RunHeaders"); + write.AddContainer("MRawRunHeader" , "RunHeaders"); + write.AddContainer("MSrcPosCam" , "RunHeaders"); + write.AddContainer("MCalibrationChargeCam" , "RunHeaders"); + // write.AddContainer("MPedPhotCam","RunHeaders"); // Attention, was in Events - Tree!! + write.AddContainer("MPedestalCam" , "RunHeaders"); + write.AddContainer("MHCalibrationRelTimeCam","RunHeaders"); + + write.AddContainer("MCerPhotEvt" , "Events"); + write.AddContainer("MRawEvtHeader" , "Events"); + write.AddContainer("MBadPixelsCam" , "Events"); + write.AddContainer("MPedPhotCam" , "Events"); + + tlist3.AddToList(&write); + + // Create and execute eventloop + MEvtLoop evtloop3; + evtloop3.SetParList(&plist3); + + cout << "*************************************************************" << endl; + cout << "*** COMPUTING DATA USING EXTRACTED SIGNAL (IN PHOTONS) ***" << endl; + cout << "*************************************************************" << endl; + + if (!evtloop3.Eventloop()) + return; + tlist3.PrintStatistics(); + +} + +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(obj1.GetName()); + +// obj2->Sumw2(); + 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; + } + + TArrayI s0(3); + s0[0] = 6; + s0[1] = 1; + s0[2] = 2; + + TArrayI s1(3); + s1[0] = 3; + s1[1] = 4; + s1[2] = 5; + + TArrayI inner(1); + inner[0] = 0; + + TArrayI outer(1); + outer[0] = 1; + + // Just to get the right (maximum) binning + TH1D *half[4]; + half[0] = obj1->ProjectionS(s0, inner, "Sector 6-1-2 Inner"); + half[1] = obj1->ProjectionS(s1, inner, "Sector 3-4-5 Inner"); + half[2] = obj1->ProjectionS(s0, outer, "Sector 6-1-2 Outer"); + half[3] = obj1->ProjectionS(s1, outer, "Sector 3-4-5 Outer"); + + for (int i=0; i<4; i++) + { + half[i]->SetLineColor(kRed+i); + half[i]->SetDirectory(0); + half[i]->SetBit(kCanDelete); + half[i]->Draw("same"); + } + + gPad->Modified(); + gPad->Update(); + + } +} + +