Changeset 3775 for trunk/MagicSoft/Mars


Ignore:
Timestamp:
04/17/04 04:45:15 (20 years ago)
Author:
gaug
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r3770 r3775  
    6565   * mcalib/MHPedestalPix.cc
    6666     - fixed a bug which made normalization to values per slice not happen
     67
     68   * macros/pedestalstudies.C
     69     - fixed and documented
     70
    6771
    6872 2004/04/15: Markus Gaug

  • trunk/MagicSoft/Mars/macros/pedestalstudies.C

    r3316 r3775  
    2222!
    2323\* ======================================================================== */
    24 
    25 //const TString pedfile = "/remote/home/pc2/operator/Crab20040214/20040215_16743_P_CrabOn_E.root";
    26 const TString pedfile = "../20040215_16770_P_OffCrab4_E.root";
    27 
    28 //const TString pedfile = "/mnt/users/mdoro/Mars/Data/20040201_14418_P_OffMrk421-1_E.root";
    29 //const TString pedfile = "/mnt/Data/rootdata/CrabNebula/2004_02_10/20040210_14607_P_CrabNebula_E.root";
    30 //const TString pedfile = "/mnt/Data/rootdata/CrabNebula/2004_01_26/20040125_10412_P_Crab-On_E.root";
    31 //const TString pedfile = "/mnt/Data/rootdata/Miscellaneous/2003_12_19/20031218_03522_P_Park_E.root";
    32 
    33 void pedestalstudies(TString pedname=pedfile)
     24//////////////////////////////////////////////////////////////////////////////
     25//
     26// pedestalstudies.C
     27//
     28// macro to observe the pedestals and pedestalRMS with the number of FADC
     29// slices summed up.
     30//
     31// In order to use this macro, you have to uncomment the following
     32// line in MPedCalcPedRun (line 214):
     33//
     34// fNumHiGainSamples = runheader->GetNumSamplesHiGain() & ~1;
     35//
     36//
     37/////////////////////////////////////////////////////////////////////////////////
     38const TString pedfile = "./20040303_20123_P_NewCalBoxTestLidOpen_E.root";
     39
     40void pedestalstudies(const TString pedname=pedfile)
    3441{
    3542
    36     gStyle->SetOptStat(1111);
    37     gStyle->SetOptFit();
    38 
    39     MStatusDisplay *display = new MStatusDisplay;
    40     display->SetUpdateTime(500);
    41     display->Resize(850,700);
    42 
    43     //
    44     // Create a empty Parameter List and an empty Task List
    45     // The tasklist is identified in the eventloop by its name
    46     //
    47     MParList  plist;
    48     MTaskList tlist;
    49     plist.AddToList(&tlist);
    50 
    51     //
    52     // Now setup the tasks and tasklist for the pedestals:
    53     // ---------------------------------------------------
    54     //
    55 
    56     MReadMarsFile read("Events", pedname);
    57     read.DisableAutoScheme();
    58 
    59     MGeomApply      geomapl;
    60     MExtractSignal  sigcalc;
    61 
    62     //
    63     // Set the extraction range higher:
    64     //         
    65     //sigcalc.SetRange(1,14,1,14);
    66 
    67     MPedCalcPedRun pedcalc;
    68 
    69     //
    70     // Additionally to calculating the pedestals,
    71     // you can fill histograms and look at them
    72     //
    73     MFillH fill("MHPedestalCam", "MExtractedSignalCam");
    74 
    75     tlist.AddToList(&read);
    76     tlist.AddToList(&geomapl);
    77     tlist.AddToList(&sigcalc);
    78     tlist.AddToList(&pedcalc);
    79     tlist.AddToList(&fill);
    80 
    81     MGeomCamMagic  geomcam;
    82     MPedestalCam   pedcam;
    83     MHPedestalCam  hpedcam;
    84     plist.AddToList(&geomcam);
    85     plist.AddToList(&pedcam);
    86     plist.AddToList(&hpedcam);
    87 
    88     //
    89     // Create and setup the eventloop
    90     //
    91     MEvtLoop evtloop;
    92 
    93     evtloop.SetParList(&plist);
    94     evtloop.SetDisplay(display);
    95  
    96     //
    97     // Execute first analysis
    98     //
    99     if (!evtloop.Eventloop())
     43  Int_t loops = 13;
     44  Int_t stepsize = 2;
     45
     46  gStyle->SetOptStat(1111);
     47  gStyle->SetOptFit();
     48 
     49  TArrayF *hmeandiffinn = new TArrayF(loops);
     50  TArrayF *hrmsdiffinn  = new TArrayF(loops);
     51  TArrayF *hmeandiffout = new TArrayF(loops);
     52  TArrayF *hrmsdiffout  = new TArrayF(loops);
     53  TArrayF *hmeaninn  = new TArrayF(loops);
     54  TArrayF *hmeanout  = new TArrayF(loops);
     55  TArrayF *hrmsinn   = new TArrayF(loops);
     56  TArrayF *hrmsout   = new TArrayF(loops);
     57  TArrayF *hmuinn    = new TArrayF(loops);
     58  TArrayF *hmuout    = new TArrayF(loops);
     59  TArrayF *hsigmainn = new TArrayF(loops);
     60  TArrayF *hsigmaout = new TArrayF(loops);
     61
     62  TArrayF *hmeandiffinnerr = new TArrayF(loops);
     63  TArrayF *hrmsdiffinnerr  = new TArrayF(loops);
     64  TArrayF *hmeandiffouterr = new TArrayF(loops);
     65  TArrayF *hrmsdiffouterr  = new TArrayF(loops);
     66  TArrayF *hmeaninnerr  = new TArrayF(loops);
     67  TArrayF *hmeanouterr  = new TArrayF(loops);
     68  TArrayF *hrmsinnerr   = new TArrayF(loops);
     69  TArrayF *hrmsouterr   = new TArrayF(loops);
     70  TArrayF *hmuinnerr    = new TArrayF(loops);
     71  TArrayF *hmuouterr    = new TArrayF(loops);
     72  TArrayF *hsigmainnerr = new TArrayF(loops);
     73  TArrayF *hsigmaouterr = new TArrayF(loops);
     74
     75
     76  MStatusDisplay *display = new MStatusDisplay;
     77  display->SetUpdateTime(500);
     78  display->Resize(850,700);
     79     
     80  //
     81  // Create a empty Parameter List and an empty Task List
     82  // The tasklist is identified in the eventloop by its name
     83  //
     84  MParList  plist;
     85  MTaskList tlist;
     86  plist.AddToList(&tlist);
     87
     88  for (Int_t samples=2;samples<stepsize*loops+1;samples=samples+stepsize)
     89    {
     90
     91      plist.Reset();
     92      tlist.Reset();
     93     
     94      //
     95      // Now setup the tasks and tasklist for the pedestals:
     96      // ---------------------------------------------------
     97      //
     98     
     99      MReadMarsFile read("Events", pedname);
     100      read.DisableAutoScheme();
     101     
     102      MGeomApply      geomapl;
     103      //
     104      // Set the extraction range higher:
     105      //               
     106      //sigcalc.SetRange(1,14,1,14);
     107      MPedCalcPedRun pedcalc;
     108      pedcalc.SetNumHiGainSamples((Byte_t)samples);
     109     
     110      MExtractSignal  sigcalc;
     111      sigcalc.SetRange(0,samples-1,0,samples-1);
     112 
     113      //
     114      // Additionally to calculating the pedestals,
     115      // you can fill histograms and look at them
     116      //
     117      MFillH fill("MHPedestalCam", "MExtractedSignalCam");
     118      fill.SetNameTab(Form("%s%2d","PedCam",samples));
     119
     120      tlist.AddToList(&read);
     121      tlist.AddToList(&geomapl);
     122      tlist.AddToList(&sigcalc);
     123      tlist.AddToList(&pedcalc);
     124      tlist.AddToList(&fill);
     125     
     126      MGeomCamMagic      geomcam;
     127      MPedestalCam       pedcam;
     128      MBadPixelsCam      badcam;
     129      badcam.AsciiRead("badpixels.dat"); 
     130     
     131      MHPedestalCam      hpedcam;
     132      MCalibrationPedCam cpedcam;
     133     
     134      plist.AddToList(&geomcam);
     135      plist.AddToList(&pedcam);
     136      plist.AddToList(&hpedcam);
     137      plist.AddToList(&cpedcam);
     138      plist.AddToList(&badcam);
     139     
     140      //
     141      // Create and setup the eventloop
     142      //
     143      MEvtLoop evtloop;
     144     
     145      evtloop.SetParList(&plist);
     146      evtloop.SetDisplay(display);
     147
     148      //
     149      // Execute first analysis
     150      //
     151      if (!evtloop.Eventloop())
    100152        return;
    101 
    102     tlist.PrintStatistics();
    103 
    104     //
    105     // Look at one specific pixel, after all the histogram manipulations:
    106     //
    107 //    hpedcam[9].DrawClone("fourierevents");
    108 
    109 
    110     MHCamera dispped0  (geomcam, "Ped;Pedestal",               "Mean per Slice");
    111     MHCamera dispped1  (geomcam, "Ped;PedestalErr",            "Mean Error per Slice");
    112     MHCamera dispped2  (geomcam, "Ped;PedestalRms",            "RMS per Slice");
    113     MHCamera dispped3  (geomcam, "Ped;PedestalRmsErr",         "RMS Error per Slice");
    114 
    115     MHCamera dispped4  (geomcam, "Ped;Mean",                   "Fitted Mean per Slice");
    116     MHCamera dispped5  (geomcam, "Ped;MeanErr",                "Fitted Error of Mean per Slice");
    117     MHCamera dispped6  (geomcam, "Ped;Sigma",                  "Fitted Sigma per Slice");
    118     MHCamera dispped7  (geomcam, "Ped;SigmaErr",               "Fitted Error of Sigma per Slice");
    119     MHCamera dispped8  (geomcam, "Ped;Prob",                   "Probability of Fit");
    120     MHCamera dispped9  (geomcam, "Ped;DeltaPedestalMean",      "Rel. Diff. Mean per Slice (Calc.-Fitte)");
    121     MHCamera dispped10 (geomcam, "Ped;DeltaPedestalMeanError", "Rel. Diff. Mean Error per Slice (Calc.-Fitted)");
    122     MHCamera dispped11  (geomcam, "Ped;DeltaRmsSigma",         "Rel. Diff. RMS per Slice (Calc.-Fitted)");
    123     MHCamera dispped12  (geomcam, "Ped;DeltaRmsSigmaError",    "Rel. Diff. RMS Error per Slice (Calc.-Fitted)");
    124     MHCamera dispped13  (geomcam, "Ped;FitOK",                 "Gaus Fit not OK");
    125     MHCamera dispped14  (geomcam, "Ped;FourierOK",             "Fourier Analysis not OK");
    126 
    127     dispped0.SetCamContent(  pedcam, 0);
    128     dispped0.SetCamError(    pedcam, 1);
    129     dispped1.SetCamContent(  pedcam, 1);
    130     dispped2.SetCamContent(  pedcam, 2);
    131     dispped2.SetCamError(    pedcam, 3);
    132     dispped3.SetCamContent(  pedcam, 3);
    133 
    134     dispped4.SetCamContent( hpedcam, 0);
    135     dispped4.SetCamError(   hpedcam, 1);
    136     dispped5.SetCamContent( hpedcam, 1);
    137     dispped6.SetCamContent( hpedcam, 2);
    138     dispped6.SetCamError(   hpedcam, 3);
    139     dispped7.SetCamContent( hpedcam, 3);
    140     dispped8.SetCamContent( hpedcam, 4);
    141     dispped9.SetCamContent( hpedcam, 5);
    142     dispped9.SetCamError(   hpedcam, 6);
    143     dispped10.SetCamContent(hpedcam, 7);
    144     dispped11.SetCamContent(hpedcam, 8);
    145     dispped11.SetCamError(  hpedcam, 9);
    146     dispped12.SetCamContent(hpedcam, 10);
    147     dispped13.SetCamContent(hpedcam, 11);
    148     dispped14.SetCamContent(hpedcam, 12);
    149 
    150     dispped0.SetYTitle("Calc. Pedestal per slice [FADC counts]");
    151     dispped1.SetYTitle("Calc. Pedestal Error per slice [FADC counts]");
    152     dispped2.SetYTitle("Calc. Pedestal RMS per slice [FADC counts]");
    153     dispped3.SetYTitle("Calc. Pedestal RMS Error per slice [FADC counts]");
    154     dispped4.SetYTitle("Fitted Mean per slice [FADC counts]");
    155     dispped5.SetYTitle("Error of Fitted Mean per slice [FADC counts]");
    156     dispped6.SetYTitle("Fitted Sigma per slice [FADC counts]");
    157     dispped7.SetYTitle("Error of Fitted Sigma per slice [FADC counts]");
    158     dispped8.SetYTitle("Fit Probability [1]");
    159     dispped9.SetYTitle("Rel. Diff. Pedestal Calc.-Fitted per slice [1]");
    160     dispped10.SetYTitle("Rel. Diff. Pedestal Error Calc.-Fitted per slice [1]");
    161     dispped11.SetYTitle("Rel. Diff. Pedestal RMS Calc.-Fitted per slice [1]");
    162     dispped12.SetYTitle("Rel. Diff. Pedestal RMS Error Calc.-Fitted per slice [1]");
    163     dispped13.SetYTitle("[1]");
    164     dispped14.SetYTitle("[1]");
    165    
    166     // Histogram values
    167     TCanvas &b1 = display->AddTab("Ped.Calc.");
    168     b1.Divide(4,3);
    169 
    170     CamDraw(b1,dispped0,pedcam,1,4,1);
    171     CamDraw(b1,dispped1,pedcam,2,4,2);
    172     CamDraw(b1,dispped2,pedcam,3,4,2);
    173     CamDraw(b1,dispped3,pedcam,4,4,2);
    174 
    175     // Fitted values
    176     TCanvas &b2 = display->AddTab("Ped.Fit");
    177     b2.Divide(4,3);
    178 
    179     CamDraw(b2,dispped4,hpedcam,1,4,1);
    180     CamDraw(b2,dispped5,hpedcam,2,4,2);
    181     CamDraw(b2,dispped6,hpedcam,3,4,2);
    182     CamDraw(b2,dispped7,hpedcam,4,4,2);
    183 
    184 
    185     // Fits Probability
    186     TCanvas &b3 = display->AddTab("Ped.Fit Prob.");
    187     b3.Divide(1,3);
    188 
    189     CamDraw(b3,dispped8,hpedcam,1,1,3);
    190 
    191     // Differences
    192     TCanvas &c4 = display->AddTab("Rel.Diff.Calc.-Fit");
    193     c4.Divide(4,3);
    194 
    195     CamDraw(c4,dispped9,hpedcam,1,4,1);
    196     CamDraw(c4,dispped10,hpedcam,2,4,1);
    197     CamDraw(c4,dispped11,hpedcam,3,4,1);
    198     CamDraw(c4,dispped12,hpedcam,4,4,1);
    199 
    200     // Defects
    201     TCanvas &c5 = display->AddTab("Defects");
    202     c5.Divide(2,2);
    203 
    204     CamDraw(c5,dispped13,hpedcam,1,2,0);
    205     CamDraw(c5,dispped14,hpedcam,2,2,0);
     153     
     154      //
     155      // Look at one specific pixel, after all the histogram manipulations:
     156      //
     157      /*
     158      MHGausEvents &hpix = hpedcam.GetAverageHiGainArea(0);
     159      hpix.DrawClone("fourierevents");
     160     
     161      MHGausEvents &lpix = hpedcam.GetAverageHiGainArea(1);
     162      lpix.DrawClone("fourierevents");
     163
     164      hpedcam[170].DrawClone("fourierevents");
     165 
     166      */
     167
     168      MHCamera dispped0  (geomcam, "Ped;Pedestal",       "Mean per Slice");
     169      MHCamera dispped2  (geomcam, "Ped;PedestalRms",    "RMS per Slice");
     170      MHCamera dispped4  (geomcam, "Ped;Mean",           "Fitted Mean per Slice");
     171      MHCamera dispped6  (geomcam, "Ped;Sigma",          "Fitted Sigma per Slice");
     172      MHCamera dispped9  (geomcam, "Ped;DeltaPedMean",   "Rel. Diff. Mean per Slice (Fit-Calc.)");
     173      MHCamera dispped11 (geomcam, "Ped;DeltaRmsSigma",  "Rel. Diff. RMS per Slice (Fit-Calc.)");
     174 
     175      dispped0.SetCamContent(  pedcam, 0);
     176      dispped0.SetCamError(    pedcam, 1);
     177      dispped2.SetCamContent(  pedcam, 2);
     178      dispped2.SetCamError(    pedcam, 3);
     179     
     180      dispped4.SetCamContent( hpedcam, 0);
     181      dispped4.SetCamError(   hpedcam, 1);
     182      dispped6.SetCamContent( hpedcam, 2);
     183      dispped6.SetCamError(   hpedcam, 3);
     184      dispped9.SetCamContent( hpedcam, 5);
     185      dispped9.SetCamError(   hpedcam, 6);
     186      dispped11.SetCamContent(hpedcam, 8);
     187      dispped11.SetCamError(  hpedcam, 9);
     188 
     189      dispped0.SetYTitle("Calc. Pedestal per slice [FADC counts]");
     190      dispped2.SetYTitle("Calc. Pedestal RMS per slice [FADC counts]");
     191      dispped4.SetYTitle("Fitted Mean per slice [FADC counts]");
     192      dispped6.SetYTitle("Fitted Sigma per slice [FADC counts]");
     193      dispped9.SetYTitle("Rel. Diff. Pedestal per slice Fit-Calc [1]");
     194      dispped11.SetYTitle("Rel. Diff. Pedestal RMS per slice Fit-Calc [1]");
     195
     196
     197      // Histogram values
     198      TCanvas &b1 = display->AddTab(Form("%s%d","MeanRMS",samples));
     199      b1.Divide(4,3);
     200
     201      CamDraw(b1,dispped0,1,4,*hmeaninn,*hmeanout,*hmeaninnerr,*hmeanouterr,samples,stepsize);
     202      CamDraw(b1,dispped2,2,4,*hrmsinn,*hrmsout,*hrmsinnerr,*hrmsouterr,samples,stepsize); 
     203      CamDraw(b1,dispped4,3,4,*hmuinn,*hmuout,*hmuinnerr,*hmuouterr,samples,stepsize);
     204      CamDraw(b1,dispped6,4,4,*hsigmainn,*hsigmaout,*hsigmainnerr,*hsigmaouterr,samples,stepsize);
     205     
     206      display->SaveAsGIF(3*((samples-1)/stepsize)+2,Form("%s%d","MeanRmsSamples",samples));
     207
     208      // Differences
     209      TCanvas &c4 = display->AddTab(Form("%s%d","RelDiff",samples));
     210      c4.Divide(2,3);
     211     
     212      CamDraw(c4,dispped9,1,2,*hmeandiffinn,*hmeandiffout,*hmeandiffinnerr,*hmeandiffouterr,samples,stepsize);
     213      CamDraw(c4,dispped11,2,2,*hrmsdiffinn,*hrmsdiffout,*hrmsdiffinnerr,*hrmsdiffouterr,samples,stepsize);
     214
     215      display->SaveAsGIF(3*((samples-1)/stepsize)+3,Form("%s%d","RelDiffSamples",samples));
     216
     217    }
     218
     219  TF1 *logg = new TF1("logg","[1]+TMath::Log(x-[0])",1.,30.,2);
     220  logg->SetParameters(1.,3.5);
     221  logg->SetParLimits(0,-1.,3.);
     222  logg->SetParLimits(1,-1.,7.);
     223  logg->SetLineColor(kRed);
     224
     225  TCanvas *canvas = new TCanvas("PedstudInner","Pedestal Studies Inner Pixels",600,900);
     226  canvas->Divide(2,3);
     227  canvas->cd(1);
     228
     229  TGraphErrors *gmeaninn = new TGraphErrors(hmeaninn->GetSize(),
     230                                            CreateXaxis(hmeaninn->GetSize(),stepsize),hmeaninn->GetArray(),
     231                                            CreateXaxisErr(hmeaninnerr->GetSize(),stepsize),hmeaninnerr->GetArray());
     232  gmeaninn->Draw("A*");
     233  gmeaninn->SetTitle("Calculated Mean per Slice Inner Pixels");
     234  gmeaninn->GetXaxis()->SetTitle("Nr. added FADC slices");
     235  gmeaninn->GetYaxis()->SetTitle("Calculated Mean per slice");
     236  gmeaninn->Fit("pol0");
     237  gmeaninn->GetFunction("pol0")->SetLineColor(kGreen);
     238  //  gmeaninn->Fit(logg);
     239
     240  canvas->cd(2);
     241
     242  TGraphErrors *gmuinn = new TGraphErrors(hmuinn->GetSize(),
     243                                          CreateXaxis(hmuinn->GetSize(),stepsize),hmuinn->GetArray(),
     244                                          CreateXaxisErr(hmuinnerr->GetSize(),stepsize),hmuinnerr->GetArray());
     245  gmuinn->Draw("A*");
     246  gmuinn->SetTitle("Fitted Mean per Slice Inner Pixels");
     247  gmuinn->GetXaxis()->SetTitle("Nr. added FADC slices");
     248  gmuinn->GetYaxis()->SetTitle("Fitted Mean per Slice");
     249  gmuinn->Fit("pol0");
     250  gmuinn->GetFunction("pol0")->SetLineColor(kGreen);
     251  //gmuinn->Fit(logg);
     252
     253
     254  canvas->cd(3);
     255
     256  TGraphErrors *grmsinn = new TGraphErrors(hrmsinn->GetSize(),
     257                                           CreateXaxis(hrmsinn->GetSize(),stepsize),hrmsinn->GetArray(),
     258                                           CreateXaxisErr(hrmsinnerr->GetSize(),stepsize),hrmsinnerr->GetArray());
     259  grmsinn->Draw("A*");
     260  grmsinn->SetTitle("Calculated Rms per Slice Inner Pixels");
     261  grmsinn->GetXaxis()->SetTitle("Nr. added FADC slices");
     262  grmsinn->GetYaxis()->SetTitle("Calculated Rms per Slice");
     263  //grmsinn->Fit("pol2");
     264  //grmsinn->GetFunction("pol2")->SetLineColor(kRed);
     265  grmsinn->Fit(logg);
     266
     267  canvas->cd(4);
     268
     269  TGraphErrors *gsigmainn = new TGraphErrors(hsigmainn->GetSize(),
     270                                             CreateXaxis(hsigmainn->GetSize(),stepsize),hsigmainn->GetArray(),
     271                                             CreateXaxisErr(hsigmainnerr->GetSize(),stepsize),hsigmainnerr->GetArray());
     272  gsigmainn->Draw("A*");
     273  gsigmainn->SetTitle("Fitted Sigma per Slice Inner Pixels");
     274  gsigmainn->GetXaxis()->SetTitle("Nr. added FADC slices");
     275  gsigmainn->GetYaxis()->SetTitle("Fitted Sigma per Slice");
     276  //  gsigmainn->Fit("pol2");
     277  //  gsigmainn->GetFunction("pol2")->SetLineColor(kRed);
     278  gsigmainn->Fit(logg);
     279
     280  canvas->cd(5);
     281
     282  TGraphErrors *gmeandiffinn = new TGraphErrors(hmeandiffinn->GetSize(),
     283                                                CreateXaxis(hmeandiffinn->GetSize(),stepsize),hmeandiffinn->GetArray(),
     284                                                CreateXaxisErr(hmeandiffinnerr->GetSize(),stepsize),hmeandiffinnerr->GetArray());
     285  gmeandiffinn->Draw("A*");
     286  gmeandiffinn->SetTitle("Rel. Difference  Mean per Slice Inner Pixels");
     287  gmeandiffinn->GetXaxis()->SetTitle("Nr. added FADC slices");
     288  gmeandiffinn->GetYaxis()->SetTitle("Rel. Difference Mean per Slice");
     289  //gmeandiffinn->Fit("pol2");
     290  //gmeandiffinn->GetFunction("pol2")->SetLineColor(kBlue);
     291  gmeandiffinn->Fit(logg);
     292
     293
     294  canvas->cd(6);
     295
     296  TGraphErrors *grmsdiffinn = new TGraphErrors(hrmsdiffinn->GetSize(),
     297                                               CreateXaxis(hrmsdiffinn->GetSize(),stepsize),hrmsdiffinn->GetArray(),
     298                                               CreateXaxisErr(hrmsdiffinnerr->GetSize(),stepsize),hrmsdiffinnerr->GetArray());
     299  grmsdiffinn->Draw("A*");
     300  grmsdiffinn->SetTitle("Rel. Difference Sigma per Slice-RMS Inner Pixels");
     301  grmsdiffinn->GetXaxis()->SetTitle("Nr. added FADC slices");
     302  grmsdiffinn->GetYaxis()->SetTitle("Rel. Difference Sigma per Slice-RMS");
     303  //grmsdiffinn->Fit("pol2");
     304  //grmsdiffinn->GetFunction("pol2")->SetLineColor(kBlue);
     305  grmsdiffinn->Fit(logg);
     306
     307
     308  TCanvas *canvas2 = new TCanvas("PedstudOut","Pedestal Studies Outer Pixels",600,900);
     309  canvas2->Divide(2,3);
     310  canvas2->cd(1);
     311
     312
     313  canvas2->cd(1);
     314
     315  TGraphErrors *gmeanout = new TGraphErrors(hmeanout->GetSize(),
     316                                            CreateXaxis(hmeanout->GetSize(),stepsize),hmeanout->GetArray(),
     317                                            CreateXaxisErr(hmeanouterr->GetSize(),stepsize),hmeanouterr->GetArray());
     318  gmeanout->Draw("A*");
     319  gmeanout->SetTitle("Calculated Mean per Slice Outer Pixels");
     320  gmeanout->GetXaxis()->SetTitle("Nr. added FADC slices");
     321  gmeanout->GetYaxis()->SetTitle("Calculated Mean per Slice");
     322  gmeanout->Fit("pol0");
     323  gmeanout->GetFunction("pol0")->SetLineColor(kGreen);
     324  //gmeanout->Fit(logg);
     325
     326  canvas2->cd(2);
     327
     328  TGraphErrors *gmuout = new TGraphErrors(hmuout->GetSize(),
     329                                          CreateXaxis(hmuout->GetSize(),stepsize),hmuout->GetArray(),
     330                                          CreateXaxisErr(hmuouterr->GetSize(),stepsize),hmuouterr->GetArray());
     331  gmuout->Draw("A*");
     332  gmuout->SetTitle("Fitted Mean per Slice Outer Pixels");
     333  gmuout->GetXaxis()->SetTitle("Nr. added FADC slices");
     334  gmuout->GetYaxis()->SetTitle("Fitted Mean per Slice");
     335  gmuout->Fit("pol0");
     336  gmuout->GetFunction("pol0")->SetLineColor(kGreen);
     337  //gmuout->Fit(logg);
     338
     339  canvas2->cd(3);
     340
     341  TGraphErrors *grmsout = new TGraphErrors(hrmsout->GetSize(),
     342                                           CreateXaxis(hrmsout->GetSize(),stepsize),hrmsout->GetArray(),
     343                                           CreateXaxisErr(hrmsouterr->GetSize(),stepsize),hrmsouterr->GetArray());
     344  grmsout->Draw("A*");
     345  grmsout->SetTitle("Calculated Rms per Slice Outer Pixels");
     346  grmsout->GetXaxis()->SetTitle("Nr. added FADC slices");
     347  grmsout->GetYaxis()->SetTitle("Calculated Rms per Slice");
     348  //grmsout->Fit("pol2");
     349  //grmsout->GetFunction("pol2")->SetLineColor(kRed);
     350  grmsout->Fit(logg);
     351
     352  canvas2->cd(4);
     353
     354  TGraphErrors *gsigmaout = new TGraphErrors(hsigmaout->GetSize(),
     355                                             CreateXaxis(hsigmaout->GetSize(),stepsize),hsigmaout->GetArray(),
     356                                             CreateXaxisErr(hsigmaouterr->GetSize(),stepsize),hsigmaouterr->GetArray());
     357  gsigmaout->Draw("A*");
     358  gsigmaout->SetTitle("Fitted Sigma per Slice Outer Pixels");
     359  gsigmaout->GetXaxis()->SetTitle("Nr. added FADC slices");
     360  gsigmaout->GetYaxis()->SetTitle("Fitted Sigma per Slice");
     361  //gsigmaout->Fit("pol2");
     362  //gsigmaout->GetFunction("pol2")->SetLineColor(kRed);
     363  gsigmaout->Fit(logg);
     364
     365
     366  canvas2->cd(5);
     367
     368  TGraphErrors *gmeandiffout = new TGraphErrors(hmeandiffout->GetSize(),
     369                                                CreateXaxis(hmeandiffout->GetSize(),stepsize),hmeandiffout->GetArray(),
     370                                                CreateXaxisErr(hmeandiffouterr->GetSize(),stepsize),hmeandiffouterr->GetArray());
     371  gmeandiffout->Draw("A*");
     372  gmeandiffout->SetTitle("Rel. Difference  Mean per Slice Outer Pixels");
     373  gmeandiffout->GetXaxis()->SetTitle("Nr. added FADC slices");
     374  gmeandiffout->GetYaxis()->SetTitle("Rel. Difference Mean per Slice");
     375  //gmeandiffout->Fit("pol2");
     376  //gmeandiffout->GetFunction("pol2")->SetLineColor(kBlue);
     377  gmeandiffout->Fit(logg);
     378
     379  canvas2->cd(6);
     380
     381  TGraphErrors *grmsdiffout = new TGraphErrors(hrmsdiffout->GetSize(),
     382                                               CreateXaxis(hrmsdiffout->GetSize(),stepsize),hrmsdiffout->GetArray(),
     383                                               CreateXaxisErr(hrmsdiffouterr->GetSize(),stepsize),hrmsdiffouterr->GetArray());
     384  grmsdiffout->Draw("A*");
     385  grmsdiffout->SetTitle("Rel. Difference Sigma per Slice-RMS Outer Pixels");
     386  grmsdiffout->GetXaxis()->SetTitle("Nr. added FADC slices");
     387  grmsdiffout->GetYaxis()->SetTitle("Rel. Difference Sigma per Slice-RMS");
     388  //grmsdiffout->Fit("pol2");
     389  //grmsdiffout->GetFunction("pol2")->SetLineColor(kBlue);
     390  grmsdiffout->Fit(logg);
     391
    206392
    207393}
    208394
    209 void CamDraw(TCanvas &c, MHCamera &cam, MCamEvent &evt, Int_t i, Int_t j, Int_t fit)
     395
     396void CamDraw(TCanvas &c, MHCamera &cam, Int_t i, Int_t j, TArrayF &a1, TArrayF &a2,
     397             TArrayF &a1err, TArrayF &a2err, Int_t samp, Int_t stepsize)
    210398{
    211399
    212400  c.cd(i);
    213   gPad->SetBorderMode(0);
    214401  MHCamera *obj1=(MHCamera*)cam.DrawCopy("hist");
    215   //  obj1->AddNotify(evt);
     402  obj1->SetDirectory(NULL);
    216403 
    217404  c.cd(i+j);
    218   gPad->SetBorderMode(0);
    219405  obj1->Draw();
    220406  ((MHCamera*)obj1)->SetPrettyPalette();
    221407
    222   if (fit != 0)
     408  c.cd(i+2*j);
     409  TH1D *obj2 = (TH1D*)obj1->Projection();
     410  obj2->SetDirectory(NULL);
     411 
     412  //      obj2->Sumw2();
     413  obj2->Draw();
     414  obj2->SetBit(kCanDelete);
     415 
     416  const Double_t min   = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
     417  const Double_t max   = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
     418  const Double_t integ = obj2->Integral("width")/2.5066283;
     419  const Double_t mean  = obj2->GetMean();
     420  const Double_t rms   = obj2->GetRMS();
     421  const Double_t width = max-min;
     422 
     423  if (rms == 0. || width == 0. )
     424    return;
     425 
     426  TArrayI s0(6);
     427  s0[0] = 6;
     428  s0[1] = 1;
     429  s0[2] = 2;
     430  s0[3] = 3;
     431  s0[4] = 4;
     432  s0[5] = 5;
     433 
     434  TArrayI inner(1);
     435  inner[0] = 0;
     436 
     437  TArrayI outer(1);
     438  outer[0] = 1;
     439     
     440  // Just to get the right (maximum) binning
     441  TH1D *half[2];
     442  half[0] = obj1->ProjectionS(s0, inner, "Inner");
     443  half[1] = obj1->ProjectionS(s0, outer, "Outer");
     444
     445  half[0]->SetDirectory(NULL);
     446  half[1]->SetDirectory(NULL);
     447 
     448  for (int i=0; i<2; i++)
    223449    {
    224       c.cd(i+2*j);
    225       gPad->SetBorderMode(0);
    226       TH1D *obj2 = (TH1D*)obj1->Projection();
    227      
    228 //      obj2->Sumw2();
    229       obj2->Draw();
    230       obj2->SetBit(kCanDelete);
    231 
    232       const Double_t min   = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
    233       const Double_t max   = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
    234       const Double_t integ = obj2->Integral("width")/2.5066283;
    235       const Double_t mean  = obj2->GetMean();
    236       const Double_t rms   = obj2->GetRMS();
    237       const Double_t width = max-min;
    238 
    239       if (rms == 0. || width == 0. )
    240         return;
    241      
    242       switch (fit)
     450      half[i]->SetLineColor(kRed+i);
     451      half[i]->SetDirectory(0);
     452      half[i]->SetBit(kCanDelete);
     453      half[i]->Draw("same");
     454      half[i]->Fit("gaus","Q+");
     455
     456      if (i==0)
    243457        {
    244         case 1:
    245           TF1 *sgaus = new TF1("sgaus","gaus(0)",min,max);
    246           sgaus->SetBit(kCanDelete);
    247           sgaus->SetParNames("Area","#mu","#sigma");
    248           sgaus->SetParameters(integ/rms,mean,rms);
    249           sgaus->SetParLimits(0,0.,integ);
    250           sgaus->SetParLimits(1,min,max);
    251           sgaus->SetParLimits(2,0,width/1.5);
    252           obj2->Fit("sgaus","QLR");
    253           obj2->GetFunction("sgaus")->SetLineColor(kYellow);
    254           break;
    255 
    256         case 2:
    257           TString dgausform = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
    258           dgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
    259           TF1 *dgaus = new TF1("dgaus",dgausform.Data(),min,max);
    260           dgaus->SetBit(kCanDelete);
    261           dgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}","A_{2}","#mu_{2}","#sigma_{2}");
    262           dgaus->SetParameters(integ,(min+mean)/2.,width/4.,
    263                                integ/width/2.,(max+mean)/2.,width/4.);
    264           // The left-sided Gauss
    265           dgaus->SetParLimits(0,integ-1.5,integ+1.5);
    266           dgaus->SetParLimits(1,min+(width/10.),mean);
    267           dgaus->SetParLimits(2,0,width/2.);
    268           // The right-sided Gauss
    269           dgaus->SetParLimits(3,0,integ);
    270           dgaus->SetParLimits(4,mean,max-(width/10.));
    271           dgaus->SetParLimits(5,0,width/2.);
    272           obj2->Fit("dgaus","QLRM");
    273           obj2->GetFunction("dgaus")->SetLineColor(kYellow);
    274           break;
    275          
    276         case 3:
    277           TString tgausform = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
    278           tgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
    279           tgausform += "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])";
    280           TF1 *tgaus = new TF1("tgaus",tgausform.Data(),min,max);
    281           tgaus->SetBit(kCanDelete);
    282           tgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}",
    283                              "A_{2}","#mu_{2}","#sigma_{2}",
    284                              "A_{3}","#mu_{3}","#sigma_{3}");
    285           tgaus->SetParameters(integ,(min+mean)/2,width/4.,
    286                                integ/width/3.,(max+mean)/2.,width/4.,
    287                                integ/width/3.,mean,width/2.);
    288           // The left-sided Gauss
    289           tgaus->SetParLimits(0,integ-1.5,integ+1.5);
    290           tgaus->SetParLimits(1,min+(width/10.),mean);
    291           tgaus->SetParLimits(2,width/15.,width/2.);
    292           // The right-sided Gauss
    293           tgaus->SetParLimits(3,0.,integ);
    294           tgaus->SetParLimits(4,mean,max-(width/10.));
    295           tgaus->SetParLimits(5,width/15.,width/2.);
    296           // The Gauss describing the outliers
    297           tgaus->SetParLimits(6,0.,integ);
    298           tgaus->SetParLimits(7,min,max);
    299           tgaus->SetParLimits(8,width/4.,width/1.5);
    300           obj2->Fit("tgaus","QLRM");
    301           obj2->GetFunction("tgaus")->SetLineColor(kYellow);
    302           break;
    303         case 4:
    304           obj2->Fit("pol0","Q");
    305           obj2->GetFunction("pol0")->SetLineColor(kYellow);
    306           break;
    307         case 9:
    308           break;
    309         default:
    310           obj2->Fit("gaus","Q");
    311           obj2->GetFunction("gaus")->SetLineColor(kYellow);
    312           break;
     458          a1[(samp-1)/stepsize] = half[i]->GetFunction("gaus")->GetParameter(1);
     459          a1err[(samp-1)/stepsize] = half[i]->GetFunction("gaus")->GetParError(1);
     460          if (a1err[(samp-1)/stepsize] > 3.)
     461            a1err[(samp-1)/stepsize] = 1.;
    313462        }
    314      
    315       gPad->Modified();
    316       gPad->Update();
    317      
     463     if (i==1)
     464       {
     465         a2[(samp-1)/stepsize] = half[i]->GetFunction("gaus")->GetParameter(1);
     466         a2err[(samp-1)/stepsize] = half[i]->GetFunction("gaus")->GetParError(1);
     467          if (a2err[(samp-1)/stepsize] > 3.)
     468            a2err[(samp-1)/stepsize] = 1.;
     469       }
    318470    }
     471 
     472 
    319473}
     474
     475// -----------------------------------------------------------------------------
     476//
     477// Create the x-axis for the event graph
     478//
     479Float_t *CreateXaxis(Int_t n, Int_t step)
     480{
     481
     482  Float_t *xaxis = new Float_t[n];
     483
     484  for (Int_t i=0;i<n;i++)
     485    xaxis[i] = 2. + step*i;
     486
     487  return xaxis;
     488                 
     489}
     490
     491// -----------------------------------------------------------------------------
     492//
     493// Create the x-axis for the event graph
     494//
     495Float_t *CreateXaxisErr(Int_t n, Int_t step)
     496{
     497
     498  Float_t *xaxis = new Float_t[n];
     499
     500  for (Int_t i=0;i<n;i++)
     501    xaxis[i] = step/2.;
     502
     503  return xaxis;
     504                 
     505}

  • trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.cc

    r3701 r3775  
    118118    AddToBranchList("fHiGainFadcSamples");
    119119
     120    fNumHiGainSamples = 0;
    120121    Clear();
    121122}
    … …  
    123124void MPedCalcPedRun::Clear(const Option_t *o)
    124125{
    125     fNumHiGainSamples = 0;
    126     fNumSamplesTot    = 0;
    127 
    128     fRawEvt    = NULL;
    129     fPedestals = NULL;
     126
     127  fNumSamplesTot    = 0;
     128
     129  fRawEvt    = NULL;
     130  fPedestals = NULL;
    130131}
    131132
    … …  
    239240      UInt_t sum = 0;
    240241      UInt_t sqr = 0;
    241        
     242
    242243      do
    243244        {
Note: See TracChangeset for help on using the changeset viewer.