Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 3928) +++ trunk/MagicSoft/Mars/Changelog (revision 3929) @@ -20,4 +20,8 @@ 2004/05/03: Markus Gaug + + * mjobs/MGCamDisplays.[h,cc] + * mhist/MHCamera.[h,cc] + - put DrawProjection and DrawRadialProfile in MHCamera * macros/dohtml.C Index: trunk/MagicSoft/Mars/mhist/MHCamera.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHCamera.cc (revision 3928) +++ trunk/MagicSoft/Mars/mhist/MHCamera.cc (revision 3929) @@ -57,4 +57,6 @@ #include #include +#include +#include #include "MLog.h" @@ -68,5 +70,4 @@ #include "MCamEvent.h" - #define kItemsLegend 48 // see SetPalette(1,0) @@ -1528,2 +1529,288 @@ return TH1D::DrawCopy(fName+";cpy"); } + + +// -------------------------------------------------------------------------- +// +// Draw a projection of MHCamera onto the y-axis values. Depending on the +// variable fit, the following fits are performed: +// +// 0: No fit, simply draw the projection +// 1: Single Gauss (for distributions flat-fielded over the whole camera) +// 2: Double Gauss (for distributions different for inner and outer pixels) +// 3: Triple Gauss (for distributions with inner, outer pixels and outliers) +// 4: flat (for the probability distributions) +// (1-4:) Moreover, sectors 6,1 and 2 of the camera and sectors 3,4 and 5 are +// drawn separately, for inner and outer pixels. +// 5: Fit Inner and Outer pixels separately by a single Gaussian +// (only for MAGIC cameras) +// 6: Fit Inner and Outer pixels separately by a single Gaussian and display +// additionally the two camera halfs separately (for MAGIC camera) +// +// +void MHCamera::DrawProjection(Int_t fit) const +{ + + TArrayI inner(1); + inner[0] = 0; + + TArrayI outer(1); + outer[0] = 1; + + if (fit==5 || fit==6) + { + + if (GetGeomCam().InheritsFrom("MGeomCamMagic")) + { + TArrayI s0(6); + s0[0] = 6; + s0[1] = 1; + s0[2] = 2; + s0[3] = 3; + s0[4] = 4; + s0[5] = 5; + + TArrayI s1(3); + s1[0] = 6; + s1[1] = 1; + s1[2] = 2; + + TArrayI s2(3); + s2[0] = 3; + s2[1] = 4; + s2[2] = 5; + + gPad->Clear(); + TVirtualPad *pad = gPad; + pad->Divide(2,1); + + TH1D *inout[2]; + inout[0] = ProjectionS(s0, inner, "Inner"); + inout[1] = ProjectionS(s0, outer, "Outer"); + + inout[0]->SetDirectory(NULL); + inout[1]->SetDirectory(NULL); + + for (int i=0; i<2; i++) + { + pad->cd(i+1); + inout[i]->SetLineColor(kRed+i); + inout[i]->SetBit(kCanDelete); + inout[i]->Draw(); + inout[i]->Fit("gaus","Q"); + + if (fit == 6) + { + TH1D *half[2]; + half[0] = ProjectionS(s1, i==0 ? inner : outer , "Sector 6-1-2"); + half[1] = ProjectionS(s2, i==0 ? inner : outer , "Sector 3-4-5"); + + for (int j=0; j<2; j++) + { + half[j]->SetLineColor(kRed+i+j); + half[j]->SetDirectory(0); + half[j]->SetBit(kCanDelete); + half[j]->Draw("same"); + } + } + + } + + gLog << all << GetName() + << Form("%s%5.3f%s%3.2f"," Mean: Inner Pixels: ", + inout[0]->GetFunction("gaus")->GetParameter(1),"+-", + inout[0]->GetFunction("gaus")->GetParError(1)); + gLog << Form("%s%5.3f%s%3.2f"," Outer Pixels: ", + inout[1]->GetFunction("gaus")->GetParameter(1),"+-", + inout[1]->GetFunction("gaus")->GetParError(1)) << endl; + gLog << all << GetName() + << Form("%s%5.3f%s%3.2f"," Sigma: Inner Pixels: ", + inout[0]->GetFunction("gaus")->GetParameter(2),"+-", + inout[0]->GetFunction("gaus")->GetParError(2)); + gLog << Form("%s%5.3f%s%3.2f"," Outer Pixels: ", + inout[1]->GetFunction("gaus")->GetParameter(2),"+-", + inout[1]->GetFunction("gaus")->GetParError(2)) << endl; + + } + return; + } + + TH1D *obj2 = (TH1D*)Projection(GetName()); + obj2->SetDirectory(0); + obj2->Draw(); + obj2->SetBit(kCanDelete); + + if (fit == 0) + return; + + if (GetGeomCam().InheritsFrom("MGeomCamMagic")) + { + TArrayI s0(3); + s0[0] = 6; + s0[1] = 1; + s0[2] = 2; + + TArrayI s1(3); + s1[0] = 3; + s1[1] = 4; + s1[2] = 5; + + + TH1D *halfInOut[4]; + + // Just to get the right (maximum) binning + halfInOut[0] = ProjectionS(s0, inner, "Sector 6-1-2 Inner"); + halfInOut[1] = ProjectionS(s1, inner, "Sector 3-4-5 Inner"); + halfInOut[2] = ProjectionS(s0, outer, "Sector 6-1-2 Outer"); + halfInOut[3] = ProjectionS(s1, outer, "Sector 3-4-5 Outer"); + + for (int i=0; i<4; i++) + { + halfInOut[i]->SetLineColor(kRed+i); + halfInOut[i]->SetDirectory(0); + halfInOut[i]->SetBit(kCanDelete); + halfInOut[i]->Draw("same"); + } + } + + 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.5; + const Double_t mean = obj2->GetMean(); + const Double_t rms = obj2->GetRMS(); + const Double_t width = max-min; + + const TString dgausformula = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])" + "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])"; + + const TString tgausformula = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])" + "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])" + "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])"; + TF1 *f=0; + switch (fit) + { + case 1: + f = new TF1("sgaus", "gaus(0)", min, max); + f->SetLineColor(kYellow); + f->SetBit(kCanDelete); + f->SetParNames("Area", "#mu", "#sigma"); + f->SetParameters(integ/rms, mean, rms); + f->SetParLimits(0, 0, integ); + f->SetParLimits(1, min, max); + f->SetParLimits(2, 0, width/1.5); + + obj2->Fit(f, "QLR"); + break; + + case 2: + f = new TF1("dgaus",dgausformula.Data(),min,max); + f->SetLineColor(kYellow); + f->SetBit(kCanDelete); + f->SetParNames("A_{tot}", "#mu1", "#sigma1", "A2", "#mu2", "#sigma2"); + f->SetParameters(integ,(min+mean)/2.,width/4., + integ/width/2.,(max+mean)/2.,width/4.); + // The left-sided Gauss + f->SetParLimits(0,integ-1.5 , integ+1.5); + f->SetParLimits(1,min+(width/10.), mean); + f->SetParLimits(2,0 , width/2.); + // The right-sided Gauss + f->SetParLimits(3,0 , integ); + f->SetParLimits(4,mean, max-(width/10.)); + f->SetParLimits(5,0 , width/2.); + obj2->Fit(f,"QLRM"); + break; + + case 3: + f = new TF1("tgaus",tgausformula.Data(),min,max); + f->SetLineColor(kYellow); + f->SetBit(kCanDelete); + f->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}", + "A_{2}","#mu_{2}","#sigma_{2}", + "A_{3}","#mu_{3}","#sigma_{3}"); + f->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 + f->SetParLimits(0,integ-1.5,integ+1.5); + f->SetParLimits(1,min+(width/10.),mean); + f->SetParLimits(2,width/15.,width/2.); + // The right-sided Gauss + f->SetParLimits(3,0.,integ); + f->SetParLimits(4,mean,max-(width/10.)); + f->SetParLimits(5,width/15.,width/2.); + // The Gauss describing the outliers + f->SetParLimits(6,0.,integ); + f->SetParLimits(7,min,max); + f->SetParLimits(8,width/4.,width/1.5); + obj2->Fit(f,"QLRM"); + 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; + } +} + + +// -------------------------------------------------------------------------- +// +// Draw a projection of MHCamera vs. the radius from the central pixel. +// +// The inner and outer pixels are drawn separately, both fitted by a polynomial +// of grade 1. +// +void MHCamera::DrawRadialProfile() const +{ + + TProfile *obj2 = (TProfile*)RadialProfile(GetName()); + obj2->SetDirectory(0); + obj2->Draw(); + obj2->SetBit(kCanDelete); + + if (GetGeomCam().InheritsFrom("MGeomCamMagic")) + { + + TArrayI s0(6); + s0[0] = 1; + s0[1] = 2; + s0[2] = 3; + s0[3] = 4; + s0[4] = 5; + s0[5] = 6; + + TArrayI inner(1); + inner[0] = 0; + + TArrayI outer(1); + outer[0] = 1; + + // Just to get the right (maximum) binning + TProfile *half[2]; + half[0] = RadialProfileS(s0, inner,Form("%s%s",GetName(),"Inner")); + half[1] = RadialProfileS(s0, outer,Form("%s%s",GetName(),"Outer")); + + for (Int_t i=0; i<2; i++) + { + Double_t min = GetGeomCam().GetMinRadius(i); + Double_t max = GetGeomCam().GetMaxRadius(i); + + half[i]->SetLineColor(kRed+i); + half[i]->SetDirectory(0); + half[i]->SetBit(kCanDelete); + half[i]->Draw("same"); + half[i]->Fit("pol1","Q","",min,max); + half[i]->GetFunction("pol1")->SetLineColor(kRed+i); + half[i]->GetFunction("pol1")->SetLineWidth(1); + } + } +} + Index: trunk/MagicSoft/Mars/mhist/MHCamera.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHCamera.h (revision 3928) +++ trunk/MagicSoft/Mars/mhist/MHCamera.h (revision 3929) @@ -188,4 +188,7 @@ void Draw(Option_t *option=""); TObject *DrawClone(Option_t *option="") const; + void DrawProjection (Int_t fit=0) const; + void DrawRadialProfile() const; + void SavePrimitive(ofstream &out, Option_t *); Int_t DistancetoPrimitive(Int_t px, Int_t py); Index: trunk/MagicSoft/Mars/mjobs/MGCamDisplays.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MGCamDisplays.cc (revision 3928) +++ trunk/MagicSoft/Mars/mjobs/MGCamDisplays.cc (revision 3929) @@ -39,6 +39,4 @@ #include "MGeomCam.h" #include "TVirtualPad.h" -#include "TProfile.h" -#include "TF1.h" #include "MLog.h" @@ -99,5 +97,5 @@ gPad->SetBorderMode(0); gPad->SetTicks(); - DrawRadialProfile(obj1); + obj1->DrawRadialProfile(); } @@ -109,285 +107,6 @@ gPad->SetBorderMode(0); gPad->SetTicks(); - DrawProjection(obj1, fit); + obj1->DrawProjection(fit); } -// -------------------------------------------------------------------------- -// -// Draw a projection of MHCamera onto the y-axis values. Depending on the -// variable fit, the following fits are performed: -// -// 1: Single Gauss (for distributions flat-fielded over the whole camera) -// 2: Double Gauss (for distributions different for inner and outer pixels) -// 3: Triple Gauss (for distributions with inner, outer pixels and outliers) -// 4: flat (for the probability distributions) -// (1-4:) Moreover, sectors 6,1 and 2 of the camera and sectors 3,4 and 5 are -// drawn separately, for inner and outer pixels. -// 5: Fit Inner and Outer pixels separately by a single Gaussian -// (only for MAGIC cameras) -// 6: Fit Inner and Outer pixels separately by a single Gaussian and display -// additionally the two camera halfs separately (for MAGIC camera) -// -// -void MGCamDisplays::DrawProjection(MHCamera *obj, Int_t fit) const -{ - - TArrayI inner(1); - inner[0] = 0; - - TArrayI outer(1); - outer[0] = 1; - - if (fit==5 || fit==6) - { - - if (obj->GetGeomCam().InheritsFrom("MGeomCamMagic")) - { - TArrayI s0(6); - s0[0] = 6; - s0[1] = 1; - s0[2] = 2; - s0[3] = 3; - s0[4] = 4; - s0[5] = 5; - TArrayI s1(3); - s1[0] = 6; - s1[1] = 1; - s1[2] = 2; - - TArrayI s2(3); - s2[0] = 3; - s2[1] = 4; - s2[2] = 5; - - gPad->Clear(); - TVirtualPad *pad = gPad; - pad->Divide(2,1); - - TH1D *inout[2]; - inout[0] = obj->ProjectionS(s0, inner, "Inner"); - inout[1] = obj->ProjectionS(s0, outer, "Outer"); - - inout[0]->SetDirectory(NULL); - inout[1]->SetDirectory(NULL); - - for (int i=0; i<2; i++) - { - pad->cd(i+1); - inout[i]->SetLineColor(kRed+i); - inout[i]->SetBit(kCanDelete); - inout[i]->Draw(); - inout[i]->Fit("gaus","Q"); - - if (fit == 6) - { - TH1D *half[2]; - half[0] = obj->ProjectionS(s1, i==0 ? inner : outer , "Sector 6-1-2"); - half[1] = obj->ProjectionS(s2, i==0 ? inner : outer , "Sector 3-4-5"); - - for (int j=0; j<2; j++) - { - half[j]->SetLineColor(kRed+i+j); - half[j]->SetDirectory(0); - half[j]->SetBit(kCanDelete); - half[j]->Draw("same"); - } - } - - } - - gLog << all << obj->GetName() - << Form("%s%5.3f%s%3.2f"," Mean: Inner Pixels: ", - inout[0]->GetFunction("gaus")->GetParameter(1),"+-", - inout[0]->GetFunction("gaus")->GetParError(1)); - gLog << Form("%s%5.3f%s%3.2f"," Outer Pixels: ", - inout[1]->GetFunction("gaus")->GetParameter(1),"+-", - inout[1]->GetFunction("gaus")->GetParError(1)) << endl; - gLog << all << obj->GetName() - << Form("%s%5.3f%s%3.2f"," Sigma: Inner Pixels: ", - inout[0]->GetFunction("gaus")->GetParameter(2),"+-", - inout[0]->GetFunction("gaus")->GetParError(2)); - gLog << Form("%s%5.3f%s%3.2f"," Outer Pixels: ", - inout[1]->GetFunction("gaus")->GetParameter(2),"+-", - inout[1]->GetFunction("gaus")->GetParError(2)) << endl; - - } - return; - } - - TH1D *obj2 = (TH1D*)obj->Projection(obj->GetName()); - obj2->SetDirectory(0); - obj2->Draw(); - obj2->SetBit(kCanDelete); - - - if (obj->GetGeomCam().InheritsFrom("MGeomCamMagic")) - { - TArrayI s0(3); - s0[0] = 6; - s0[1] = 1; - s0[2] = 2; - - TArrayI s1(3); - s1[0] = 3; - s1[1] = 4; - s1[2] = 5; - - - TH1D *halfInOut[4]; - - // Just to get the right (maximum) binning - halfInOut[0] = obj->ProjectionS(s0, inner, "Sector 6-1-2 Inner"); - halfInOut[1] = obj->ProjectionS(s1, inner, "Sector 3-4-5 Inner"); - halfInOut[2] = obj->ProjectionS(s0, outer, "Sector 6-1-2 Outer"); - halfInOut[3] = obj->ProjectionS(s1, outer, "Sector 3-4-5 Outer"); - - for (int i=0; i<4; i++) - { - halfInOut[i]->SetLineColor(kRed+i); - halfInOut[i]->SetDirectory(0); - halfInOut[i]->SetBit(kCanDelete); - halfInOut[i]->Draw("same"); - } - } - - 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.5; - const Double_t mean = obj2->GetMean(); - const Double_t rms = obj2->GetRMS(); - const Double_t width = max-min; - - const TString dgausformula = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])" - "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])"; - - const TString tgausformula = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])" - "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])" - "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])"; - TF1 *f=0; - switch (fit) - { - case 1: - f = new TF1("sgaus", "gaus(0)", min, max); - f->SetLineColor(kYellow); - f->SetBit(kCanDelete); - f->SetParNames("Area", "#mu", "#sigma"); - f->SetParameters(integ/rms, mean, rms); - f->SetParLimits(0, 0, integ); - f->SetParLimits(1, min, max); - f->SetParLimits(2, 0, width/1.5); - - obj2->Fit(f, "QLR"); - break; - - case 2: - f = new TF1("dgaus",dgausformula.Data(),min,max); - f->SetLineColor(kYellow); - f->SetBit(kCanDelete); - f->SetParNames("A_{tot}", "#mu1", "#sigma1", "A2", "#mu2", "#sigma2"); - f->SetParameters(integ,(min+mean)/2.,width/4., - integ/width/2.,(max+mean)/2.,width/4.); - // The left-sided Gauss - f->SetParLimits(0,integ-1.5 , integ+1.5); - f->SetParLimits(1,min+(width/10.), mean); - f->SetParLimits(2,0 , width/2.); - // The right-sided Gauss - f->SetParLimits(3,0 , integ); - f->SetParLimits(4,mean, max-(width/10.)); - f->SetParLimits(5,0 , width/2.); - obj2->Fit(f,"QLRM"); - break; - - case 3: - f = new TF1("tgaus",tgausformula.Data(),min,max); - f->SetLineColor(kYellow); - f->SetBit(kCanDelete); - f->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}", - "A_{2}","#mu_{2}","#sigma_{2}", - "A_{3}","#mu_{3}","#sigma_{3}"); - f->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 - f->SetParLimits(0,integ-1.5,integ+1.5); - f->SetParLimits(1,min+(width/10.),mean); - f->SetParLimits(2,width/15.,width/2.); - // The right-sided Gauss - f->SetParLimits(3,0.,integ); - f->SetParLimits(4,mean,max-(width/10.)); - f->SetParLimits(5,width/15.,width/2.); - // The Gauss describing the outliers - f->SetParLimits(6,0.,integ); - f->SetParLimits(7,min,max); - f->SetParLimits(8,width/4.,width/1.5); - obj2->Fit(f,"QLRM"); - 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; - } -} - -// -------------------------------------------------------------------------- -// -// Draw a projection of MHCamera vs. the radius from the central pixel. -// -// The inner and outer pixels are drawn separately, both fitted by a polynomial -// of grade 1. -// -void MGCamDisplays::DrawRadialProfile(MHCamera *obj) const -{ - - TProfile *obj2 = (TProfile*)obj->RadialProfile(obj->GetName()); - obj2->SetDirectory(0); - obj2->Draw(); - obj2->SetBit(kCanDelete); - - if (obj->GetGeomCam().InheritsFrom("MGeomCamMagic")) - { - - TArrayI s0(6); - s0[0] = 1; - s0[1] = 2; - s0[2] = 3; - s0[3] = 4; - s0[4] = 5; - s0[5] = 6; - - TArrayI inner(1); - inner[0] = 0; - - TArrayI outer(1); - outer[0] = 1; - - // Just to get the right (maximum) binning - TProfile *half[2]; - half[0] = obj->RadialProfileS(s0, inner,Form("%s%s",obj->GetName(),"Inner")); - half[1] = obj->RadialProfileS(s0, outer,Form("%s%s",obj->GetName(),"Outer")); - - for (Int_t i=0; i<2; i++) - { - Double_t min = obj->GetGeomCam().GetMinRadius(i); - Double_t max = obj->GetGeomCam().GetMaxRadius(i); - - half[i]->SetLineColor(kRed+i); - half[i]->SetDirectory(0); - half[i]->SetBit(kCanDelete); - half[i]->Draw("same"); - half[i]->Fit("pol1","Q","",min,max); - half[i]->GetFunction("pol1")->SetLineColor(kRed+i); - half[i]->GetFunction("pol1")->SetLineWidth(1); - } - } -} - Index: trunk/MagicSoft/Mars/mjobs/MGCamDisplays.h =================================================================== --- trunk/MagicSoft/Mars/mjobs/MGCamDisplays.h (revision 3928) +++ trunk/MagicSoft/Mars/mjobs/MGCamDisplays.h (revision 3929) @@ -8,11 +8,8 @@ class TCanvas; class MHCamera; -class MCamEvent; class MGCamDisplays { protected: - void DrawProjection ( MHCamera *obj, Int_t fit) const; // Draw projection of pixels values - void DrawRadialProfile( MHCamera *obj) const; // Draw projection of pixels values onto camera radius void CamDraw(TCanvas &c, const Int_t x, const Int_t y, const MHCamera &cam1, const Int_t fit, const Int_t rad=0,