Ignore:
Timestamp:
05/03/04 10:11:15 (21 years ago)
Author:
gaug
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mjobs
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/mjobs/MGCamDisplays.cc

    r3889 r3929  
    3939#include "MGeomCam.h"
    4040#include "TVirtualPad.h"
    41 #include "TProfile.h"
    42 #include "TF1.h"
    4341
    4442#include "MLog.h"
     
    9997        gPad->SetBorderMode(0);
    10098        gPad->SetTicks();
    101         DrawRadialProfile(obj1);
     99        obj1->DrawRadialProfile();
    102100      }
    103101   
     
    109107    gPad->SetBorderMode(0);
    110108    gPad->SetTicks();
    111     DrawProjection(obj1, fit);
     109    obj1->DrawProjection(fit);
    112110}
    113111
    114 // --------------------------------------------------------------------------
    115 //
    116 // Draw a projection of MHCamera onto the y-axis values. Depending on the
    117 // variable fit, the following fits are performed:
    118 //
    119 // 1: Single Gauss (for distributions flat-fielded over the whole camera)
    120 // 2: Double Gauss (for distributions different for inner and outer pixels)
    121 // 3: Triple Gauss (for distributions with inner, outer pixels and outliers)
    122 // 4: flat         (for the probability distributions)
    123 // (1-4:) Moreover, sectors 6,1 and 2 of the camera and sectors 3,4 and 5 are
    124 //        drawn separately, for inner and outer pixels.
    125 // 5: Fit Inner and Outer pixels separately by a single Gaussian
    126 //                 (only for MAGIC cameras)
    127 // 6: Fit Inner and Outer pixels separately by a single Gaussian and display
    128 //                 additionally the two camera halfs separately (for MAGIC camera)
    129 //
    130 //
    131 void MGCamDisplays::DrawProjection(MHCamera *obj, Int_t fit) const
    132 {
    133  
    134   TArrayI inner(1);
    135   inner[0] = 0;
    136  
    137   TArrayI outer(1);
    138   outer[0] = 1;
    139          
    140   if (fit==5 || fit==6)
    141     {
    142      
    143       if (obj->GetGeomCam().InheritsFrom("MGeomCamMagic"))
    144         {
    145           TArrayI s0(6);
    146           s0[0] = 6;
    147           s0[1] = 1;
    148           s0[2] = 2;
    149           s0[3] = 3;
    150           s0[4] = 4;
    151           s0[5] = 5;
    152112
    153           TArrayI s1(3);
    154           s1[0] = 6;
    155           s1[1] = 1;
    156           s1[2] = 2;
    157          
    158           TArrayI s2(3);
    159           s2[0] = 3;
    160           s2[1] = 4;
    161           s2[2] = 5;
    162 
    163           gPad->Clear();
    164           TVirtualPad *pad = gPad;
    165           pad->Divide(2,1);
    166          
    167           TH1D *inout[2];
    168           inout[0] = obj->ProjectionS(s0, inner, "Inner");
    169           inout[1] = obj->ProjectionS(s0, outer, "Outer");
    170          
    171           inout[0]->SetDirectory(NULL);
    172           inout[1]->SetDirectory(NULL);
    173          
    174           for (int i=0; i<2; i++)
    175             {
    176               pad->cd(i+1);
    177               inout[i]->SetLineColor(kRed+i);
    178               inout[i]->SetBit(kCanDelete);
    179               inout[i]->Draw();
    180               inout[i]->Fit("gaus","Q");
    181 
    182               if (fit == 6)
    183                 {
    184                   TH1D *half[2];
    185                   half[0] = obj->ProjectionS(s1, i==0 ? inner : outer , "Sector 6-1-2");
    186                   half[1] = obj->ProjectionS(s2, i==0 ? inner : outer , "Sector 3-4-5");
    187                  
    188                   for (int j=0; j<2; j++)
    189                     {
    190                       half[j]->SetLineColor(kRed+i+j);
    191                       half[j]->SetDirectory(0);
    192                       half[j]->SetBit(kCanDelete);
    193                       half[j]->Draw("same");
    194                     }
    195                 }
    196              
    197             }
    198          
    199           gLog << all << obj->GetName()
    200                << Form("%s%5.3f%s%3.2f"," Mean: Inner Pixels: ",
    201                        inout[0]->GetFunction("gaus")->GetParameter(1),"+-",
    202                        inout[0]->GetFunction("gaus")->GetParError(1));
    203           gLog << Form("%s%5.3f%s%3.2f","  Outer Pixels: ",
    204                        inout[1]->GetFunction("gaus")->GetParameter(1),"+-",
    205                        inout[1]->GetFunction("gaus")->GetParError(1)) << endl;
    206           gLog << all << obj->GetName()
    207                << Form("%s%5.3f%s%3.2f"," Sigma: Inner Pixels: ",
    208                        inout[0]->GetFunction("gaus")->GetParameter(2),"+-",
    209                        inout[0]->GetFunction("gaus")->GetParError(2));
    210           gLog << Form("%s%5.3f%s%3.2f","  Outer Pixels: ",
    211                        inout[1]->GetFunction("gaus")->GetParameter(2),"+-",
    212                        inout[1]->GetFunction("gaus")->GetParError(2)) << endl;
    213 
    214         }
    215       return;
    216     }
    217  
    218   TH1D *obj2 = (TH1D*)obj->Projection(obj->GetName());
    219   obj2->SetDirectory(0);
    220   obj2->Draw();
    221   obj2->SetBit(kCanDelete);
    222  
    223  
    224   if (obj->GetGeomCam().InheritsFrom("MGeomCamMagic"))
    225     {
    226       TArrayI s0(3);
    227       s0[0] = 6;
    228       s0[1] = 1;
    229       s0[2] = 2;
    230      
    231       TArrayI s1(3);
    232       s1[0] = 3;
    233       s1[1] = 4;
    234       s1[2] = 5;
    235      
    236      
    237       TH1D *halfInOut[4];
    238      
    239       // Just to get the right (maximum) binning
    240       halfInOut[0] = obj->ProjectionS(s0, inner, "Sector 6-1-2 Inner");
    241       halfInOut[1] = obj->ProjectionS(s1, inner, "Sector 3-4-5 Inner");
    242       halfInOut[2] = obj->ProjectionS(s0, outer, "Sector 6-1-2 Outer");
    243       halfInOut[3] = obj->ProjectionS(s1, outer, "Sector 3-4-5 Outer");
    244      
    245       for (int i=0; i<4; i++)
    246         {
    247           halfInOut[i]->SetLineColor(kRed+i);
    248           halfInOut[i]->SetDirectory(0);
    249           halfInOut[i]->SetBit(kCanDelete);
    250           halfInOut[i]->Draw("same");
    251         }
    252     }
    253  
    254   const Double_t min   = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
    255   const Double_t max   = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
    256   const Double_t integ = obj2->Integral("width")/2.5;
    257   const Double_t mean  = obj2->GetMean();
    258   const Double_t rms   = obj2->GetRMS();
    259   const Double_t width = max-min;
    260  
    261   const TString dgausformula = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"
    262     "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
    263  
    264   const TString tgausformula = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"
    265     "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])"
    266     "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])";
    267   TF1 *f=0;
    268   switch (fit)
    269     {
    270     case 1:
    271       f = new TF1("sgaus", "gaus(0)", min, max);
    272       f->SetLineColor(kYellow);
    273       f->SetBit(kCanDelete);
    274       f->SetParNames("Area", "#mu", "#sigma");
    275       f->SetParameters(integ/rms, mean, rms);
    276       f->SetParLimits(0, 0,   integ);
    277       f->SetParLimits(1, min, max);
    278       f->SetParLimits(2, 0,   width/1.5);
    279      
    280       obj2->Fit(f, "QLR");
    281       break;
    282      
    283     case 2:
    284       f = new TF1("dgaus",dgausformula.Data(),min,max);
    285       f->SetLineColor(kYellow);
    286       f->SetBit(kCanDelete);
    287       f->SetParNames("A_{tot}", "#mu1", "#sigma1", "A2", "#mu2", "#sigma2");
    288       f->SetParameters(integ,(min+mean)/2.,width/4.,
    289                        integ/width/2.,(max+mean)/2.,width/4.);
    290       // The left-sided Gauss
    291       f->SetParLimits(0,integ-1.5      , integ+1.5);
    292       f->SetParLimits(1,min+(width/10.), mean);
    293       f->SetParLimits(2,0              , width/2.);
    294       // The right-sided Gauss
    295       f->SetParLimits(3,0   , integ);
    296       f->SetParLimits(4,mean, max-(width/10.));
    297       f->SetParLimits(5,0   , width/2.);
    298       obj2->Fit(f,"QLRM");
    299       break;
    300      
    301     case 3:
    302       f = new TF1("tgaus",tgausformula.Data(),min,max);
    303       f->SetLineColor(kYellow);
    304       f->SetBit(kCanDelete);
    305       f->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}",
    306                      "A_{2}","#mu_{2}","#sigma_{2}",
    307                      "A_{3}","#mu_{3}","#sigma_{3}");
    308       f->SetParameters(integ,(min+mean)/2,width/4.,
    309                        integ/width/3.,(max+mean)/2.,width/4.,
    310                        integ/width/3.,mean,width/2.);
    311       // The left-sided Gauss
    312       f->SetParLimits(0,integ-1.5,integ+1.5);
    313       f->SetParLimits(1,min+(width/10.),mean);
    314       f->SetParLimits(2,width/15.,width/2.);
    315       // The right-sided Gauss
    316       f->SetParLimits(3,0.,integ);
    317       f->SetParLimits(4,mean,max-(width/10.));
    318       f->SetParLimits(5,width/15.,width/2.);
    319       // The Gauss describing the outliers
    320       f->SetParLimits(6,0.,integ);
    321       f->SetParLimits(7,min,max);
    322       f->SetParLimits(8,width/4.,width/1.5);
    323       obj2->Fit(f,"QLRM");
    324       break;
    325      
    326     case 4:
    327       obj2->Fit("pol0", "Q");
    328       obj2->GetFunction("pol0")->SetLineColor(kYellow);
    329       break;
    330      
    331     case 9:
    332       break;
    333      
    334     default:
    335       obj2->Fit("gaus", "Q");
    336       obj2->GetFunction("gaus")->SetLineColor(kYellow);
    337       break;
    338     }
    339 }
    340 
    341 // --------------------------------------------------------------------------
    342 //
    343 // Draw a projection of MHCamera vs. the radius from the central pixel.
    344 //
    345 // The inner and outer pixels are drawn separately, both fitted by a polynomial
    346 // of grade 1.
    347 //
    348 void MGCamDisplays::DrawRadialProfile(MHCamera *obj) const
    349 {
    350  
    351   TProfile *obj2 = (TProfile*)obj->RadialProfile(obj->GetName());
    352   obj2->SetDirectory(0);
    353   obj2->Draw();
    354   obj2->SetBit(kCanDelete);
    355  
    356   if (obj->GetGeomCam().InheritsFrom("MGeomCamMagic"))
    357     {
    358 
    359       TArrayI s0(6);
    360       s0[0] = 1;
    361       s0[1] = 2;
    362       s0[2] = 3;
    363       s0[3] = 4;
    364       s0[4] = 5;
    365       s0[5] = 6;
    366 
    367       TArrayI inner(1);
    368       inner[0] = 0;
    369      
    370       TArrayI outer(1);
    371       outer[0] = 1;
    372      
    373       // Just to get the right (maximum) binning
    374       TProfile *half[2];
    375       half[0] = obj->RadialProfileS(s0, inner,Form("%s%s",obj->GetName(),"Inner"));
    376       half[1] = obj->RadialProfileS(s0, outer,Form("%s%s",obj->GetName(),"Outer"));
    377      
    378       for (Int_t i=0; i<2; i++)
    379         {
    380           Double_t min = obj->GetGeomCam().GetMinRadius(i);
    381           Double_t max = obj->GetGeomCam().GetMaxRadius(i);
    382 
    383           half[i]->SetLineColor(kRed+i);
    384           half[i]->SetDirectory(0);
    385           half[i]->SetBit(kCanDelete);
    386           half[i]->Draw("same");
    387           half[i]->Fit("pol1","Q","",min,max);
    388           half[i]->GetFunction("pol1")->SetLineColor(kRed+i);
    389           half[i]->GetFunction("pol1")->SetLineWidth(1);
    390         }
    391     }
    392 }
    393 
  • trunk/MagicSoft/Mars/mjobs/MGCamDisplays.h

    r3875 r3929  
    88class TCanvas;
    99class MHCamera;
    10 class MCamEvent;
    1110class MGCamDisplays
    1211{
    1312protected:
    1413
    15   void DrawProjection   ( MHCamera *obj, Int_t fit) const;     // Draw projection of pixels values
    16   void DrawRadialProfile( MHCamera *obj)            const;      // Draw projection of pixels values onto camera radius
    1714  void CamDraw(TCanvas &c, const Int_t x, const Int_t y, const MHCamera &cam1,
    1815               const Int_t fit, const Int_t rad=0,
Note: See TracChangeset for help on using the changeset viewer.