Ignore:
Timestamp:
02/08/04 20:45:43 (21 years ago)
Author:
gaug
Message:
*** empty log message ***
File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/mjobs/MJCalibration.cc

    r2992 r3055  
    5252#include "MGeomApply.h"
    5353#include "MExtractSignal.h"
     54#include "MExtractSignal2.h"
    5455#include "MCalibrationCalc.h"
    5556
    … …  
    5859
    5960ClassImp(MJCalibration);
    60 
    6161using namespace std;
    6262
    … …  
    6969void MJCalibration::DrawProjection(MHCamera *obj1, Int_t fit) const
    7070{
    71     TH1D *obj2 = (TH1D*)obj1->Projection();
    72     obj2->Draw();
    73     obj2->SetBit(kCanDelete);
    74 
    75     const Double_t min   = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
    76     const Double_t max   = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
    77     const Double_t integ = obj2->Integral("width")/2.5;
    78     const Double_t mean  = obj2->GetMean();
    79     const Double_t rms   = obj2->GetRMS();
    80     const Double_t width = max-min;
    81 
    82     TF1 *f=0;
    83     switch (fit)
    84     {
    85     case 0:
    86         f = new TF1("sgaus", "gaus(0)", min, max);
    87         f->SetLineColor(kYellow);
    88         f->SetBit(kCanDelete);
    89         f->SetParNames("Area", "#mu", "#sigma");
    90         f->SetParameters(integ/rms, mean, rms);
    91         f->SetParLimits(0, 0,   integ);
    92         f->SetParLimits(1, min, max);
    93         f->SetParLimits(2, 0,   width/1.5);
    94 
    95         obj2->Fit(f, "QLR");
    96         break;
    97 
     71
     72  TH1D *obj2 = (TH1D*)obj1->Projection();
     73  obj2->Draw();
     74  obj2->SetBit(kCanDelete);
     75 
     76  const Double_t min   = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
     77  const Double_t max   = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
     78  const Double_t integ = obj2->Integral("width")/2.5;
     79  const Double_t mean  = obj2->GetMean();
     80  const Double_t rms   = obj2->GetRMS();
     81  const Double_t width = max-min;
     82 
     83  const TString dgausformula = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"
     84                               "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
     85
     86  const TString tgausformula = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"
     87                               "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])"
     88                               "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])";
     89  TF1 *f=0;
     90  switch (fit)
     91    {
    9892    case 1:
    99         f = new TF1("dgaus", "gaus(0)+gaus(3)", min, max);
    100         f->SetLineColor(kYellow);
    101         f->SetBit(kCanDelete);
    102         f->SetParNames("A1", "#mu1", "#sigma1", "A2", "#mu2", "#sigma2");
    103         f->SetParameters(integ/width, max-width/6, width/4, integ/width, min+width/6, width/4);
    104         f->SetParLimits(0, 0,   integ);
    105         f->SetParLimits(1, min, max);
    106         f->SetParLimits(2, 0,   width/3);
    107         f->SetParLimits(3, 0,   integ);
    108         f->SetParLimits(4, min, max);
    109         f->SetParLimits(5, 0,   width/3);
    110 
    111         obj2->Fit(f, "QLR");
    112         break;
    113 
     93      f = new TF1("sgaus", "gaus(0)", min, max);
     94      f->SetLineColor(kYellow);
     95      f->SetBit(kCanDelete);
     96      f->SetParNames("Area", "#mu", "#sigma");
     97      f->SetParameters(integ/rms, mean, rms);
     98      f->SetParLimits(0, 0,   integ);
     99      f->SetParLimits(1, min, max);
     100      f->SetParLimits(2, 0,   width/1.5);
     101     
     102      obj2->Fit(f, "QLR");
     103      break;
     104     
    114105    case 2:
    115         f = new TF1("tgaus", "gaus(0)+gaus(3)+gaus(6)", min, max);
    116         f->SetLineColor(kYellow);
    117         f->SetBit(kCanDelete);
    118         f->SetParNames("A1", "#mu1", "#sigma1", "A2", "#mu2", "#sigma2", "A3", "#mu3", "#sigma3");
    119         f->SetParameters(integ/width, max-width/6, width/4, integ/width, min+width/6, width/4, integ/width,min+width/6, width/2);
    120         f->SetParLimits(0, 0,   integ);
    121         f->SetParLimits(1, min, max);
    122         f->SetParLimits(2, 0,   width/4);
    123         f->SetParLimits(3, 0,   integ);
    124         f->SetParLimits(4, min, max);
    125         f->SetParLimits(5, 0,   width/4);
    126         f->SetParLimits(6, 0,   integ);
    127         f->SetParLimits(7, min, max);
    128         f->SetParLimits(8, 0,   width/2);
    129 
    130         obj2->Fit(f, "QLR");
    131         break;
    132 
     106      f = new TF1("dgaus",dgausformula.Data(),min,max);
     107      f->SetLineColor(kYellow);
     108      f->SetBit(kCanDelete);
     109      f->SetParNames("A_{tot}", "#mu1", "#sigma1", "A2", "#mu2", "#sigma2");
     110      f->SetParameters(integ,(min+mean)/2.,width/4.,
     111                       integ/width/2.,(max+mean)/2.,width/4.);
     112      // The left-sided Gauss
     113      f->SetParLimits(0,integ-1.5      , integ+1.5);
     114      f->SetParLimits(1,min+(width/10.), mean);
     115      f->SetParLimits(2,0              , width/2.);
     116      // The right-sided Gauss
     117      f->SetParLimits(3,0   , integ);
     118      f->SetParLimits(4,mean, max-(width/10.));
     119      f->SetParLimits(5,0   , width/2.);
     120      obj2->Fit(f,"QLRM");
     121      break;
     122     
    133123    case 3:
    134         obj2->Fit("pol0", "Q");
    135         obj2->GetFunction("pol0")->SetLineColor(kYellow);
    136         break;
    137 
     124      f = new TF1("tgaus",tgausformula.Data(),min,max);
     125      f->SetLineColor(kYellow);
     126      f->SetBit(kCanDelete);
     127      f->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}",
     128                     "A_{2}","#mu_{2}","#sigma_{2}",
     129                     "A_{3}","#mu_{3}","#sigma_{3}");
     130      f->SetParameters(integ,(min+mean)/2,width/4.,
     131                       integ/width/3.,(max+mean)/2.,width/4.,
     132                       integ/width/3.,mean,width/2.);
     133      // The left-sided Gauss
     134      f->SetParLimits(0,integ-1.5,integ+1.5);
     135      f->SetParLimits(1,min+(width/10.),mean);
     136      f->SetParLimits(2,width/15.,width/2.);
     137      // The right-sided Gauss
     138      f->SetParLimits(3,0.,integ);
     139      f->SetParLimits(4,mean,max-(width/10.));
     140      f->SetParLimits(5,width/15.,width/2.);
     141      // The Gauss describing the outliers
     142      f->SetParLimits(6,0.,integ);
     143      f->SetParLimits(7,min,max);
     144      f->SetParLimits(8,width/4.,width/1.5);
     145      obj2->Fit(f,"QLRM");
     146      break;
     147
     148    case 4:
     149      obj2->Fit("pol0", "Q");
     150      obj2->GetFunction("pol0")->SetLineColor(kYellow);
     151      break;
     152     
    138153    case 9:
    139         break;
    140 
     154      break;
     155       
    141156    default:
    142         obj2->Fit("gaus", "Q");
    143         obj2->GetFunction("gaus")->SetLineColor(kYellow);
    144         break;
    145     }
    146 }
    147 
    148 void MJCalibration::CamDraw(TCanvas &c, Int_t x, Int_t y, MHCamera &cam1, Int_t fit)
    149 {
     157      obj2->Fit("gaus", "Q");
     158      obj2->GetFunction("gaus")->SetLineColor(kYellow);
     159      break;
     160    }
     161}
     162
     163void MJCalibration::CamDraw(TCanvas &c, const Int_t x, const Int_t y, const MHCamera &cam1, const Int_t fit)
     164{
     165
     166
    150167    c.cd(x);
    151168    gPad->SetBorderMode(0);
    … …  
    157174    obj1->Draw();
    158175
    159     c.cd(x+2*y);
    160     gPad->SetBorderMode(0);
    161     DrawProjection(obj1, fit);
    162 
    163     // tbretz: Using gStyle in this context is totally useless!
    164     //const Float_t he = gStyle->GetStatH();
    165     //const Float_t wi = gStyle->GetStatH();
    166     //gStyle->SetStatH(0.4);
    167     //gStyle->SetStatW(0.25);
    168     //gStyle->SetStatH(he);
    169     //gStyle->SetStatW(wi);
    170 }
     176    if (fit)
     177      {
     178        c.cd(x+2*y);
     179        gPad->SetBorderMode(0);
     180        DrawProjection(obj1, fit);
     181      }
     182}
     183
    171184
    172185void MJCalibration::DisplayResult(MParList &plist)
    … …  
    190203
    191204    // Create histograms to display
    192     MHCamera disp1 (geomcam, "Cal;Charge",        "Fitted Mean Charges");
    193     MHCamera disp3 (geomcam, "Cal;SigmaCharge",   "Sigma of Fitted Charges");
    194     MHCamera disp5 (geomcam, "Cal;FitProb",       "Probability of Fit");
    195     /*
    196      MHCamera disp6 (geomcam, "Cal;Time",          "Arrival Times");
    197      MHCamera disp7 (geomcam, "Cal;SigmaTime",     "Sigma of Arrival Times");
    198      MHCamera disp8 (geomcam, "Cal;TimeChiSquare", "Chi Square of Time Fit");
    199      */
    200 //    MHCamera disp9 (geomcam, "Cal;Ped",           "Pedestals");
    201 //    MHCamera disp10(geomcam, "Cal;PedRms",        "Pedestal RMS");
    202 //    MHCamera disp11(geomcam, "Cal;RSigma",         "Reduced Sigmas");
    203     MHCamera disp12(geomcam, "Cal;FFactorPhe",     "Nr. of Phe's (F-Factor Method)");
    204     MHCamera disp13(geomcam, "Cal;FFactorConv",    "Conversion Factor (F-Factor Method)");
    205     MHCamera disp14(geomcam, "Cal;BlindPixPhe",    "Nr. of Photons inside plexiglass (Blind Pixel Method)");
    206     MHCamera disp15(geomcam, "Cal;BlindPixConv",   "Conversion Factor (Blind Pixel Method)");
    207 //    MHCamera disp16(geomcam, "Cal;RSigma/Charge",  "Reduced Sigma per Charge");
    208 
    209     disp1.SetCamContent(fCalibrationCam, 0);
    210     disp1.SetCamError(fCalibrationCam, 1);
    211 
    212     disp3.SetCamContent(fCalibrationCam, 2);
    213     disp3.SetCamError(fCalibrationCam, 3);
    214 
    215     disp5.SetCamContent(fCalibrationCam, 4);
    216 
    217     /*
    218     disp6.SetCamContent(fCalibrationCam, 5);
    219     disp6.SetCamError(fCalibrationCam, 6);
    220     disp7.SetCamContent(fCalibrationCam, 6);
    221     disp8.SetCamContent(fCalibrationCam, 7);
    222     */
    223 /*
    224     disp9.SetCamContent(calcam, 8);
    225     disp9.SetCamError(calcam, 9);
    226     disp10.SetCamContent(calcam, 9);
    227     disp11.SetCamContent(fCalibrationCam, 10);
    228  */
    229     disp12.SetCamContent(fCalibrationCam, 11);
    230     disp12.SetCamError(fCalibrationCam, 12);
    231     disp13.SetCamContent(fCalibrationCam, 13);
    232     disp13.SetCamError(fCalibrationCam, 14);
    233 
    234     disp14.SetCamContent(fCalibrationCam, 15);
    235     disp15.SetCamContent(fCalibrationCam, 16);
    236  //   disp16.SetCamContent(fCalibrationCam, 17);
     205    MHCamera disp1  (geomcam, "Cal;Charge",         "Fitted Mean Charges");
     206    MHCamera disp2  (geomcam, "Cal;SigmaCharge",    "Sigma of Fitted Charges");
     207    MHCamera disp3  (geomcam, "Cal;FitProb",        "Probability of Fit");
     208    MHCamera disp4  (geomcam, "Cal;RSigma",         "Reduced Sigmas");
     209    MHCamera disp5  (geomcam, "Cal;RSigma/Charge",  "Reduced Sigma per Charge");
     210    MHCamera disp6  (geomcam, "Cal;FFactorPhe",     "Nr. of Phe's (F-Factor Method)");
     211    MHCamera disp7  (geomcam, "Cal;FFactorConv",    "Conversion Factor (F-Factor Method)");
     212    MHCamera disp8  (geomcam, "Cal;FFactorFFactor", "Total F-Factor (F-Factor Method)");
     213    MHCamera disp9  (geomcam, "Cal;BlindPixPh",     "Nr. of Photons inside plexiglass (Blind Pixel Method)");
     214    MHCamera disp10 (geomcam, "Cal;BlindPixConv",   "Conversion Factor (Blind Pixel Method)");
     215    MHCamera disp11 (geomcam, "Cal;BlindPixFFactor","Total F-Factor (Blind Pixel Method)");
     216    MHCamera disp12 (geomcam, "Cal;PINDiodePh",     "Nr. of Photons outside plexiglass (PIN Diode Method)");
     217    MHCamera disp13 (geomcam, "Cal;PINDiodeConv",   "Conversion Factor (PIN Diode Method)");
     218    MHCamera disp14 (geomcam, "Cal;PINDiodeFFactor","Total F-Factor (PIN Diode Method)");
     219    MHCamera disp15 (geomcam, "Cal;Excluded",       "Pixels previously excluded");
     220    MHCamera disp16 (geomcam, "Cal;NotFitted",      "Pixels that could not be fitted");
     221    MHCamera disp17 (geomcam, "Cal;NotFitValid",    "Pixels with not valid fit results");
     222    MHCamera disp18 (geomcam, "Cal;Oscillating",    "Oscillating Pixels");
     223    MHCamera disp19 (geomcam, "Cal;Saturation",     "Pixels with saturated Hi Gain");
     224
     225
     226    // Fitted charge means and sigmas
     227    disp1.SetCamContent(fCalibrationCam,  0);
     228    disp1.SetCamError(  fCalibrationCam,  1);
     229    disp2.SetCamContent(fCalibrationCam,  2);
     230    disp2.SetCamError(  fCalibrationCam,  3);
     231    // Fit probabilities
     232    disp3.SetCamContent(fCalibrationCam,  4);
     233
     234    // Reduced Sigmas and reduced sigmas per charge
     235    disp4.SetCamContent(fCalibrationCam,  5);
     236    disp4.SetCamError(  fCalibrationCam,  6);
     237    disp5.SetCamContent(fCalibrationCam,  7);
     238    disp5.SetCamError(  fCalibrationCam,  8);
     239
     240    // F-Factor Method
     241    disp6.SetCamContent(fCalibrationCam,  9);
     242    disp6.SetCamError(  fCalibrationCam, 10);
     243    disp7.SetCamContent(fCalibrationCam, 11);
     244    disp7.SetCamError(  fCalibrationCam, 12);
     245    disp8.SetCamContent(fCalibrationCam, 13);
     246    disp8.SetCamError(  fCalibrationCam, 14);
     247
     248    /// Blind Pixel Method
     249    disp9.SetCamContent(fCalibrationCam, 15);
     250    disp9.SetCamError(  fCalibrationCam, 16);
     251    disp10.SetCamContent(fCalibrationCam,17);
     252    disp10.SetCamError(  fCalibrationCam,18);
     253    disp11.SetCamContent(fCalibrationCam,19);
     254    disp11.SetCamError(  fCalibrationCam,20);
     255
     256    // PIN Diode Method
     257    disp12.SetCamContent(fCalibrationCam,21);
     258    disp12.SetCamError(  fCalibrationCam,22);
     259    disp13.SetCamContent(fCalibrationCam,23);
     260    disp13.SetCamError(  fCalibrationCam,24);
     261    disp14.SetCamContent(fCalibrationCam,25);
     262    disp14.SetCamError(  fCalibrationCam,26);
     263
     264    // Pixels with defects
     265    disp15.SetCamContent(fCalibrationCam,27);
     266    disp16.SetCamContent(fCalibrationCam,28);
     267    disp17.SetCamContent(fCalibrationCam,29);
     268    disp18.SetCamContent(fCalibrationCam,30);
     269
     270    // Lo Gain calibration
     271    disp19.SetCamContent(fCalibrationCam,31);
     272
    237273
    238274    disp1.SetYTitle("Charge [FADC units]");
    239     disp3.SetYTitle("\\sigma_{Charge} [FADC units]");
    240     disp5.SetYTitle("P_{Charge} [1]");
    241     /*
    242      disp6.SetYTitle("Arr. Time [FADC slice Nr.]");
    243      disp7.SetYTitle("\\sigma_{Time} [FADC slices]");
    244      disp8.SetYTitle("\\chi^{2}_{Time} [1]");
    245     disp9.SetYTitle("Ped [FADC Counts ]");
    246     disp10.SetYTitle("RMS_{Ped} [FADC Counts ]");
    247     disp11.SetYTitle("\\sqrt{\\sigma^{2}_{Charge} - RMS^{2}_{Ped}} [FADC units]");
    248      */
    249     disp12.SetYTitle("Nr. Photo-Electrons [1]");
    250     disp13.SetYTitle("Conversion Factor [PhE/FADC Count]");
    251     disp14.SetYTitle("Nr. Photons [1]");
    252     disp15.SetYTitle("Conversion Factor [Phot/FADC Count]");
    253 //    disp16.SetYTitle("Reduced Sigma / Charge [1]");
     275    disp2.SetYTitle("\\sigma_{Charge} [FADC units]");
     276    disp3.SetYTitle("P_{Charge} [1]");
     277
     278    disp4.SetYTitle("\\sqrt{\\sigma^{2}_{Charge} - RMS^{2}_{Ped}} [FADC Counts]");
     279    disp5.SetYTitle("Reduced Sigma / Mean Charge [1]");
     280
     281    disp6.SetYTitle("Nr. Photo-Electrons [1]");
     282    disp7.SetYTitle("Conversion Factor [PhE/FADC Count]");
     283    disp8.SetYTitle("\\sqrt{N_{PhE}}*\\sigma_{Charge}/\\mu_{Charge} [1]");
     284
     285    disp9.SetYTitle("Nr. Photons [1]");
     286    disp10.SetYTitle("Conversion Factor [Phot/FADC Count]");
     287    disp11.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]");
     288
     289    disp12.SetYTitle("Nr. Photons [1]");
     290    disp13.SetYTitle("Conversion Factor [Phot/FADC Count]");
     291    disp14.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]");
     292
     293    disp15.SetYTitle("[1]");
     294    disp16.SetYTitle("[1]");
     295    disp17.SetYTitle("[1]");
     296    disp18.SetYTitle("[1]");
    254297
    255298    gStyle->SetOptStat(1111);
    256299    gStyle->SetOptFit();
    257300
    258     TCanvas &c1 = fDisplay->AddTab("FitCharge");
     301    // Charges
     302    TCanvas &c1 = fDisplay->AddTab("Fit.Charge");
    259303    c1.Divide(2, 3);
    260304
    261     CamDraw(c1, 1, 2, disp1, 1);
    262     CamDraw(c1, 2, 2, disp3, 2);
     305    CamDraw(c1, 1, 2, disp1, 2);
     306    CamDraw(c1, 2, 2, disp2, 2);
    263307
    264308    // Fit Probability
    265     TCanvas &c2 = fDisplay->AddTab("FitProb");
     309    TCanvas &c2 = fDisplay->AddTab("Fit.Prob");
    266310    c2.Divide(1,3);
    267311
    268     CamDraw(c2, 1, 1, disp5, 3);
    269 /*
    270     // Times
    271     TCanvas &c3 = fDisplay->AddTab("FitTimes");
    272     c3.Divide(3,3);
    273 
    274     CamDraw(c3, 1, 3, disp6, 1);
    275     CamDraw(c3, 2, 3, disp7, 0);
    276     CamDraw(c3, 3, 3, disp8, 0);
    277 
    278     // Pedestals
    279     TCanvas &c4 = d3->AddTab("Pedestals");
    280     c4.Divide(2,3);
    281 
    282     CamDraw(c4, 1, 2, disp9,  0);
    283     CamDraw(c4, 2, 2, disp10, 1);
     312    CamDraw(c2, 1, 1, disp3, 4);
    284313
    285314    // Reduced Sigmas
    286     TCanvas &c5 = fDisplay->AddTab("RedSigma");
    287     c5.Divide(2,3);
    288 
    289     CamDraw(c5, 1, 2, disp11, 2);
    290     CamDraw(c5, 2, 2, disp16, 2);
    291     */
     315    TCanvas &c3 = fDisplay->AddTab("Red.Sigma");
     316    c3.Divide(2,3);
     317
     318    CamDraw(c3, 1, 2, disp4, 2);
     319    CamDraw(c3, 2, 2, disp5, 2);
    292320
    293321    // F-Factor Method
    294     TCanvas &c6 = fDisplay->AddTab("F-Factor");
    295     c6.Divide(2,3);
    296 
    297     CamDraw(c6, 1, 2, disp12, 1);
    298     CamDraw(c6, 2, 2, disp13, 1);
     322    TCanvas &c4 = fDisplay->AddTab("F-Factor");
     323    c4.Divide(3,3);
     324
     325    CamDraw(c4, 1, 3, disp6, 2);
     326    CamDraw(c4, 2, 3, disp7, 2);
     327    CamDraw(c4, 3, 3, disp8, 2);
    299328
    300329    // Blind Pixel Method
    301     TCanvas &c7 = fDisplay->AddTab("BlindPix");
    302     c7.Divide(2, 3);
    303 
    304     CamDraw(c7, 1, 2, disp14, 9);
    305     CamDraw(c7, 2, 2, disp15, 1);
     330    TCanvas &c5 = fDisplay->AddTab("BlindPix");
     331    c5.Divide(3, 3);
     332
     333    CamDraw(c5, 1, 3, disp9,  9);
     334    CamDraw(c5, 2, 3, disp10, 2);
     335    CamDraw(c5, 3, 3, disp11, 2);
     336
     337    // PIN Diode Method
     338    TCanvas &c6 = fDisplay->AddTab("PINDiode");
     339    c6.Divide(3,3);
     340
     341    CamDraw(c6, 1, 3, disp12, 9);
     342    CamDraw(c6, 2, 3, disp13, 2);
     343    CamDraw(c6, 3, 3, disp14, 2);
     344
     345    // Defects
     346    TCanvas &c7 = fDisplay->AddTab("Defects");
     347    c7.Divide(4,2);
     348
     349    CamDraw(c7, 1, 4, disp15, 0);
     350    CamDraw(c7, 2, 4, disp16, 0);
     351    CamDraw(c7, 3, 4, disp17, 0);
     352    CamDraw(c7, 4, 4, disp18, 0);
     353
     354    // Lo Gain Calibration
     355    TCanvas &c8 = fDisplay->AddTab("LowGain");
     356    c8.Divide(1,3);
     357
     358    CamDraw(c8, 1, 4, disp19, 0);
    306359
    307360}
    … …  
    413466
    414467    MGeomApply       apply;
    415     MExtractSignal   extract;
     468    MExtractSignal2  extract;
    416469    MCalibrationCalc calcalc;
    417470
Note: See TracChangeset for help on using the changeset viewer.