Index: trunk/MagicSoft/Mars/mhist/MHCamera.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHCamera.cc (revision 3927) +++ 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 3927) +++ 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);