Index: trunk/MagicSoft/Mars/manalysis/MPad.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPad.cc (revision 4577) +++ trunk/MagicSoft/Mars/manalysis/MPad.cc (revision 4584) @@ -64,8 +64,9 @@ #include #include +#include #include "MBinning.h" #include "MSigmabar.h" -#include "MMcEvt.hxx" +#include "MPointingPos.h" #include "MLog.h" #include "MLogManip.h" @@ -75,4 +76,6 @@ #include "MCerPhotPix.h" #include "MCerPhotEvt.h" + +#include "MH.h" #include "MPedPhotCam.h" @@ -190,4 +193,6 @@ // from sigma_k to sigma_j + //*fLog << all << "MPad::Merge2Distributions(); Entry" << endl; + Double_t eps = 1.e-10; @@ -205,7 +210,16 @@ Double_t v, s, w, t, x, u, a, b, arg; + //*fLog << "Content of histcp : " << endl; + for (Int_t j=nbinssig; j >= 1; j--) { + //*fLog << "j, hista, histb , histap : " << j << ", " + // << hista->GetBinContent(j) << ", " + // << histb->GetBinContent(j) << ", " + // << histap->GetBinContent(j) << endl; + v = histcp->GetBinContent(j); + //*fLog << "cp : j, v = " << j << ", " << v << endl; + if ( fabs(v) < eps ) continue; if (v >= 0.0) @@ -255,6 +269,6 @@ } - *fLog << all << "k, j, arg = " << k << ", " << j - << ", " << arg << endl; + //*fLog << all << "k, j, arg = " << k << ", " << j + // << ", " << arg << endl; //...................................... @@ -302,4 +316,7 @@ // check results + + //*fLog << "Content of histap, histbp, histcp : " << endl; + for (Int_t j=1; j<=nbinssig; j++) { @@ -308,4 +325,7 @@ Double_t c = histcp->GetBinContent(j); + //*fLog << "j, a, b, c = " << j << ": " << a << ", " << b << ", " + // << c << endl; + if( fabs(a-b)>3.0*eps || fabs(c)>3.0*eps ) *fLog << err << "MPad::Merge2Distributions(); inconsistency in results; j, a, b, c = " @@ -314,22 +334,28 @@ //--------------------------------------------------------------- - TCanvas &c = *(MH::MakeDefCanvas("Merging", "", 600, 600)); - c.Divide(2, 2); + TCanvas *pad = MH::MakeDefCanvas("Merging", "", 600, 600); gROOT->SetSelectedPad(NULL); - c.cd(1); + pad->Divide(2, 2); + + pad->cd(1); + gPad->SetBorderMode(0); hista->SetDirectory(NULL); - hista->DrawCopy(); + hista->Draw(); hista->SetBit(kCanDelete); - c.cd(2); + pad->cd(2); + gPad->SetBorderMode(0); histb->SetDirectory(NULL); - histb->DrawCopy(); + histb->Draw(); histb->SetBit(kCanDelete); - c.cd(3); + pad->cd(3); + gPad->SetBorderMode(0); histap->SetDirectory(NULL); - histap->DrawCopy(); + histap->Draw(); histap->SetBit(kCanDelete); + + pad->DrawClone(); //-------------------------------------------------------------------- @@ -339,4 +365,6 @@ delete histcp; + //*fLog << all << "MPad::Merge2Distributions(); Exit" << endl; + return kTRUE; } @@ -357,69 +385,125 @@ Bool_t MPad::MergeONOFFMC( - TH2D *sigthmc, TH3D *diffpixthmc, TH3D *sigmapixthmc, - TH2D *blindidthmc, TH2D *blindnthmc, - TH2D *sigthon, TH3D *diffpixthon, TH3D *sigmapixthon, - TH2D *blindidthon, TH2D *blindnthon, - TH2D *sigthoff, TH3D *diffpixthoff, TH3D *sigmapixthoff, - TH2D *blindidthoff, TH2D *blindnthoff) + TH2D& sigthmc, TH3D& diffpixthmc, TH3D& sigmapixthmc, + TH2D& blindidthmc, TH2D& blindnthmc, + TH2D& sigthon, TH3D& diffpixthon, TH3D& sigmapixthon, + TH2D& blindidthon, TH2D& blindnthon, + TH2D& sigthoff, TH3D& diffpixthoff, TH3D& sigmapixthoff, + TH2D& blindidthoff, TH2D& blindnthoff) { + //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + // for testing + /* + TAxis *xa = sigthon.GetXaxis(); + Int_t nbitheta = xa->GetNbins(); + + TAxis *ya = sigthon.GetYaxis(); + Int_t nbisig = ya->GetNbins(); + for (Int_t i=1; i<=nbitheta; i++) + { + for (Int_t j=1; j<=nbisig; j++) + { + if (i>0) + { + sigthmc.SetBinContent( i, j, (Float_t) (625000.0+(nbisig*nbisig*nbisig-j*j*j)) ); + sigthon.SetBinContent( i, j, (Float_t)(1.0) ); + sigthoff.SetBinContent( i, j, (Float_t)( (0.5*nbisig-j)*(0.5*nbisig-j) ) ); + } + else + { + sigthmc.SetBinContent( i, j, 0.0); + sigthon.SetBinContent( i, j, 0.0); + sigthoff.SetBinContent(i, j, 0.0); + } + } + } + */ + //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + *fLog << all << "----------------------------------------------------------------------------------" << endl; *fLog << all << "MPad::MergeONOFFMC(); Merge the ON, OFF and MC histograms to obtain the target distributions for the padding" << endl; - fHSigmaTheta = new TH2D( (TH2D&)*sigthon ); + fHSigmaTheta = new TH2D( (const TH2D&) sigthon ); fHSigmaTheta->SetNameTitle("2D-ThetaSigmabar", "2D-ThetaSigmabar (target)"); - fHDiffPixTheta = new TH3D( (TH3D&)*diffpixthon ); + *fLog << "after booking fHSigmaTheta" << endl; + + fHDiffPixTheta = new TH3D( (const TH3D&) diffpixthon ); fHDiffPixTheta->SetNameTitle("3D-ThetaPixDiff", "3D-ThetaPixDiff (target)"); - fHSigmaPixTheta = new TH3D( (TH3D&)*sigmapixthon ); + *fLog << "after booking fHDiffPixTheta" << endl; + + fHSigmaPixTheta = new TH3D( (const TH3D&) sigmapixthon ); fHSigmaPixTheta->SetNameTitle("3D-ThetaPixSigma", "3D-ThetaPixSigma (target)"); - - fHBlindPixNTheta = new TH2D( (TH2D&)*blindnthon ); + *fLog << "after booking fHSigmaPixTheta" << endl; + + fHBlindPixNTheta = new TH2D( (const TH2D&) blindnthon ); fHBlindPixNTheta->SetNameTitle("2D-ThetaBlindN", "2D-ThetaBlindN (target)"); - fHBlindPixIdTheta = new TH2D( (TH2D&)*blindidthon ); + *fLog << "after booking fHBlindPixNTheta" << endl; + + fHBlindPixIdTheta = new TH2D( (const TH2D&) blindidthon ); fHBlindPixIdTheta->SetNameTitle("2D-ThetaBlindId", "2D-ThetaBlindId (target)"); + *fLog << "after booking fHBlindPixIdTheta" << endl; + *fLog << all << "Histograms for the merged padding plots were booked" + << endl; + //-------------------------- - fHSigmaThetaMC = sigthmc; + fHSigmaThetaMC = &sigthmc; fHSigmaThetaMC->SetNameTitle("2D-ThetaSigmabarMC", "2D-ThetaSigmabarMC (orig.)"); - fHSigmaThetaON = sigthon; + fHSigmaThetaON = &sigthon; fHSigmaThetaON->SetNameTitle("2D-ThetaSigmabarON", "2D-ThetaSigmabarON (orig.)"); - fHSigmaThetaOFF = sigthoff; + fHSigmaThetaOFF = &sigthoff; fHSigmaThetaOFF->SetNameTitle("2D-ThetaSigmabarOFF", "2D-ThetaSigmabarOFF (orig.)"); + //*fLog << all << "addresses of SigmaTheta padding plots were copied" + // << endl; + //-------------------------- - fHBlindPixNThetaMC = blindnthmc; + fHBlindPixNThetaMC = &blindnthmc; fHBlindPixNThetaMC->SetNameTitle("2D-ThetaBlindNMC", "2D-ThetaBlindNMC (orig.)"); - fHBlindPixNThetaON = blindnthon; + fHBlindPixNThetaON = &blindnthon; fHBlindPixNThetaON->SetNameTitle("2D-ThetaBlindNON", "2D-ThetaBlindNON (orig.)"); - fHBlindPixNThetaOFF = blindnthoff; + fHBlindPixNThetaOFF = &blindnthoff; fHBlindPixNThetaOFF->SetNameTitle("2D-ThetaBlindNOFF", "2D-ThetaBlindNOFF (orig.)"); + //*fLog << all << "addresses of BlindPixTheta padding plots were copied" + // << endl; + //-------------------------- - fHBlindPixIdThetaMC = blindidthmc; + fHBlindPixIdThetaMC = &blindidthmc; fHBlindPixIdThetaMC->SetNameTitle("2D-ThetaBlindIdMC", "2D-ThetaBlindIdMC (orig.)"); - fHBlindPixIdThetaON = blindidthon; + fHBlindPixIdThetaON = &blindidthon; fHBlindPixIdThetaON->SetNameTitle("2D-ThetaBlindIdON", "2D-ThetaBlindIdON (orig.)"); - fHBlindPixIdThetaOFF = blindidthoff; + fHBlindPixIdThetaOFF = &blindidthoff; fHBlindPixIdThetaOFF->SetNameTitle("2D-ThetaBlindIdOFF", "2D-ThetaBlindIdOFF (orig.)"); + //*fLog << all << "addresses of BlindPixIdTheta padding plots were copied" + // << endl; + //-------------------------- - fHDiffPixThetaMC = diffpixthmc; + fHDiffPixThetaMC = &diffpixthmc; fHDiffPixThetaMC->SetNameTitle("3D-ThetaPixDiffMC", "3D-ThetaPixDiffMC (orig.)"); - fHDiffPixThetaON = diffpixthon; + fHDiffPixThetaON = &diffpixthon; fHDiffPixThetaON->SetNameTitle("3D-ThetaPixDiffON", "3D-ThetaPixDiffON (orig.)"); - fHDiffPixThetaOFF = diffpixthoff; + fHDiffPixThetaOFF = &diffpixthoff; fHDiffPixThetaOFF->SetNameTitle("3D-ThetaPixDiffOFF", "3D-ThetaPixDiffOFF (orig.)"); + //*fLog << all << "addresses of DiffPixTheta padding plots were copied" + // << endl; + //-------------------------- - fHSigmaPixThetaMC = sigmapixthmc; + fHSigmaPixThetaMC = &sigmapixthmc; fHSigmaPixThetaMC->SetNameTitle("3D-ThetaPixSigmaMC", "3D-ThetaPixSigmaMC (orig.)"); - fHSigmaPixThetaON = sigmapixthon; + fHSigmaPixThetaON = &sigmapixthon; fHSigmaPixThetaON->SetNameTitle("3D-ThetaPixSigmaON", "3D-ThetaPixSigmaON (orig.)"); - fHSigmaPixThetaOFF = sigmapixthoff; + fHSigmaPixThetaOFF = &sigmapixthoff; fHSigmaPixThetaOFF->SetNameTitle("3D-ThetaPixSigmaOFF", "3D-ThetaPixSigmaOFF (orig.)"); + + //*fLog << all << "addresses of SigmaPixTheta padding plots were copied" + // << endl; + //-------------------------- @@ -427,5 +511,553 @@ // get binning for fHgON, fHgOFF and fHgMC from sigthon // binning in Theta - TAxis *ax = sigthon->GetXaxis(); + TAxis *ax = sigthon.GetXaxis(); + Int_t nbinstheta = ax->GetNbins(); + TArrayD edgesx; + edgesx.Set(nbinstheta+1); + for (Int_t i=0; i<=nbinstheta; i++) + { + edgesx[i] = ax->GetBinLowEdge(i+1); + //*fLog << all << "i, theta low edge = " << i << ", " << edgesx[i] + // << endl; + } + MBinning binth; + binth.SetEdges(edgesx); + + // binning in sigmabar + TAxis *ay = sigthon.GetYaxis(); + Int_t nbinssig = ay->GetNbins(); + TArrayD edgesy; + edgesy.Set(nbinssig+1); + for (Int_t i=0; i<=nbinssig; i++) + { + edgesy[i] = ay->GetBinLowEdge(i+1); + //*fLog << all << "i, sigmabar low edge = " << i << ", " << edgesy[i] + // << endl; + } + MBinning binsg; + binsg.SetEdges(edgesy); + + + fHgMC = new TH3D; + MH::SetBinning(fHgMC, &binth, &binsg, &binsg); + fHgMC->SetNameTitle("3D-PaddingMatrixMC", "3D-PadMatrixThetaMC"); + + fHgON = new TH3D; + MH::SetBinning(fHgON, &binth, &binsg, &binsg); + fHgON->SetNameTitle("3D-PaddingMatrixON", "3D-PadMatrixThetaON"); + + fHgOFF = new TH3D; + MH::SetBinning(fHgOFF, &binth, &binsg, &binsg); + fHgOFF->SetNameTitle("3D-PaddingMatrixOFF", "3D-PadMatrixThetaOFF"); + + //*fLog << all << "fHg histograms were booked" << endl; + + //............ loop over Theta bins ........................... + + + + //*fLog << all << "MPad::MergeONOFFMC(); bins of Theta, Sigmabarold, Sigmabarnew, fraction of events to be padded" << endl; + for (Int_t l=1; l<=nbinstheta; l++) + { + //------------------------------------------- + // merge ON and OFF distributions + // input : hista is the normalized 1D distr. of sigmabar for ON data + // histb is the normalized 1D distr. of sigmabar for OFF data + // output : histap will be the merged distribution (ON-OFF) + // fHga(k,j) will tell which fraction of ON events should be + // padded from bin k to bin j + // fHgb(k,j) will tell which fraction of OFF events should be + // padded from bin k to bin j + + TH2D * fHga = new TH2D; + MH::SetBinning(fHga, &binsg, &binsg); + TH2D * fHgb = new TH2D; + MH::SetBinning(fHgb, &binsg, &binsg); + + TH1D *hista = sigthon.ProjectionY("sigon_y", l, l, ""); + hista->SetNameTitle("ON-orig", "ON (orig)"); + Stat_t suma = hista->Integral(); + if (suma != 0.0) hista->Scale(1./suma); + + TH1D *histap = new TH1D( (const TH1D&)*hista ); + histap->SetNameTitle("ON-OFF-merged)", "ON-OFF (merged)"); + + TH1D *histb = sigthoff.ProjectionY("sigoff_y", l, l, ""); + histb->SetNameTitle("OFF-orig", "OFF (orig)"); + Stat_t sumb = histb->Integral(); + if (sumb != 0.0) histb->Scale(1./sumb); + + if (suma == 0.0 || sumb == 0.0) + { + delete hista; + delete histb; + delete fHga; + delete fHgb; + delete histap; + + *fLog << err + << "cannot call Merge2Distributions(ON, OFF) for theta bin l = " + << l << "; NevON, NevOFF = " << suma << ", " << sumb << endl; + + // go to next theta bin (l) + continue; + } + + //*fLog << "call Merge2Distributions(ON, OFF) for theta bin l = " + // << l << endl; + + Merge2Distributions(hista, histb, histap, fHga, fHgb, nbinssig); + + // fill fHgON and fHgOFF + for (Int_t k=1; k<=nbinssig; k++) + for (Int_t j=1; j<=nbinssig; j++) + { + Double_t a = fHga->GetBinContent(k,j); + fHgON->SetBinContent(l, k, j, a); + + Double_t b = fHgb->GetBinContent(k,j); + fHgOFF->SetBinContent(l, k, j, b); + } + + //------------------------------------------- + // now merge ON-OFF and MC distributions + // input : histe is the result of the merging of ON and OFF distributions + // histf is the normalized 1D distr. of sigmabar for MC data + // output : histep will be the merged distribution (target distribution) + // fHge(k,j) will tell which fraction of ON, OFF events should be + // padded from bin k to bin j + // fHgf(k,j) will tell which fraction of MC events should be + // padded from bin k to bin j + + TH2D * fHge = new TH2D; + MH::SetBinning(fHge, &binsg, &binsg); + TH2D * fHgf = new TH2D; + MH::SetBinning(fHgf, &binsg, &binsg); + + TH1D *histe = new TH1D( (const TH1D&)*histap); + histe->SetNameTitle("ON-OFF-merged", "ON-OFF (merged)"); + + TH1D *histep = new TH1D( (const TH1D&)*histe); + histep->SetNameTitle("ON-OFF-MC-merged)", "ON-OFF-MC (merged)"); + + TH1D *histf = sigthmc.ProjectionY("sigmc_y", l, l, ""); + histf->SetNameTitle("MC-orig", "MC (orig)"); + Double_t sumf = histf->Integral(); + if (sumf != 0.0) histf->Scale(1./sumf); + + if (sumf == 0.0) + { + delete hista; + delete histb; + delete fHga; + delete fHgb; + delete histap; + + delete histe; + delete histf; + delete fHge; + delete fHgf; + delete histep; + + *fLog << err + << "cannot call Merge2Distributions(ON-OFF,MC) for theta bin l = " + << l << "; NevMC = " << sumf << endl; + + // go to next theta bin (l) + continue; + } + + //*fLog << "call Merge2Distributions(ON-OFF, MC) for theta bin l = " + // << l << endl; + + Merge2Distributions(histe, histf, histep, fHge, fHgf, nbinssig); + + // update fHgON and fHgOFF + UpdateHg(fHga, histap, fHge, fHgON, nbinssig, l); + UpdateHg(fHgb, histap, fHge, fHgOFF, nbinssig, l); + + // fill fHgMC + for (Int_t k=1; k<=nbinssig; k++) + for (Int_t j=1; j<=nbinssig; j++) + { + Double_t f = fHgf->GetBinContent(k,j); + fHgMC->SetBinContent(l, k, j, f); + } + + // fill target distribution fHSigmaTheta + // (only for plotting, it will not be used in the padding) + for (Int_t j=1; j<=nbinssig; j++) + { + Double_t ep = histep->GetBinContent(j); + fHSigmaTheta->SetBinContent(l, j, ep); + } + //------------------------------------------- + + delete hista; + delete histb; + delete fHga; + delete fHgb; + delete histap; + + delete histe; + delete histf; + delete fHge; + delete fHgf; + delete histep; + } + //............ end of loop over Theta bins .................... + + + // Define the target distributions + // fHDiffPixTheta, fHSigmaPixTheta + // (they are calculated as + // averages of the ON and OFF distributions) + // fHBlindPixNTheta, fHBlindPixIdTheta + // (they are calculated as + // the sum of the ON and OFF distributions) + // (fHDiffPixTheta will be used in the padding of MC events) + + //............ start of new loop over Theta bins .................... + for (Int_t j=1; j<=nbinstheta; j++) + { + // fraction of ON/OFF/MC events to be padded : fracON, fracOFF, fracMC + + Double_t fracON = fHgON->Integral (j, j, 1, nbinssig, 1, nbinssig, ""); + Double_t fracOFF = fHgOFF->Integral(j, j, 1, nbinssig, 1, nbinssig, ""); + Double_t fracMC = fHgMC->Integral(j, j, 1, nbinssig, 1, nbinssig, ""); + + Double_t normON; + Double_t normOFF; + + // set ranges for 2D-projections of 3D-histograms + TAxis *ax = diffpixthon.GetXaxis(); + ax->SetRange(j, j); + ax = diffpixthoff.GetXaxis(); + ax->SetRange(j, j); + + TH2D *hist; + TH2D *histOFF; + + TH1D *hist1; + TH1D *hist1OFF; + + // weights for ON and OFF distrubtions when + // calculating the weighted average + Double_t facON = 0.5 * (1. - fracON + fracOFF); + Double_t facOFF = 0.5 * (1. + fracON - fracOFF); + + *fLog << all << "Theta bin j : fracON, fracOFF, fracMC, facON, facOFF = " + << j << ", " << fracON << ", " << fracOFF << ", " + << fracMC << ", " << facON << ", " << facOFF << endl; + + //------------------------------------------------------------------ + // define target distribution 'sigma^2-sigmavar^2 vs. pixel number' + ay = diffpixthon.GetYaxis(); + Int_t nbinspixel = ay->GetNbins(); + + TAxis *az = diffpixthon.GetZaxis(); + Int_t nbinsdiff = az->GetNbins(); + + hist = (TH2D*)diffpixthon.Project3D("zy"); + hist->SetName("dummy"); + histOFF = (TH2D*)diffpixthoff.Project3D("zy"); + + normON = hist->Integral(); + normOFF = histOFF->Integral(); + if (normON != 0.0) hist->Scale(1.0/normON); + if (normOFF != 0.0) histOFF->Scale(1.0/normOFF); + + // weighted average of ON and OFF distributions + hist->Scale(facON); + hist->Add(histOFF, facOFF); + + for (Int_t k=1; k<=nbinspixel; k++) + for (Int_t l=1; l<=nbinsdiff; l++) + { + Double_t cont = hist->GetBinContent(k,l); + fHDiffPixTheta->SetBinContent(j, k, l, cont); + } + + delete hist; + delete histOFF; + + //------------------------------------------------------------------ + // define target distribution 'sigma vs. pixel number' + ay = sigmapixthon.GetYaxis(); + nbinspixel = ay->GetNbins(); + + az = sigmapixthon.GetZaxis(); + Int_t nbinssigma = az->GetNbins(); + + hist = (TH2D*)sigmapixthon.Project3D("zy"); + hist->SetName("dummy"); + histOFF = (TH2D*)sigmapixthoff.Project3D("zy"); + + normON = hist->Integral(); + normOFF = histOFF->Integral(); + if (normON != 0.0) hist->Scale(1.0/normON); + if (normOFF != 0.0) histOFF->Scale(1.0/normOFF); + + // weighted average of ON and OFF distributions + hist->Scale(facON); + hist->Add(histOFF, facOFF); + + for (Int_t k=1; k<=nbinspixel; k++) + for (Int_t l=1; l<=nbinssigma; l++) + { + Double_t cont = hist->GetBinContent(k,l); + fHSigmaPixTheta->SetBinContent(j, k, l, cont); + } + + delete hist; + delete histOFF; + + + //------------------------------------------------------------------ + // define target distribution 'number of blind pixels per event' + ay = blindnthon.GetYaxis(); + Int_t nbinsn = ay->GetNbins(); + + hist1 = fHBlindPixNThetaON->ProjectionY("", j, j, ""); + hist1->SetName("dummy"); + hist1OFF = fHBlindPixNThetaOFF->ProjectionY("", j, j, ""); + + normON = hist1->Integral(); + normOFF = hist1OFF->Integral(); + if (normON != 0.0) hist1->Scale(1.0/normON); + if (normOFF != 0.0) hist1OFF->Scale(1.0/normOFF); + + // sum of ON and OFF distributions + hist1->Add(hist1OFF, 1.0); + Stat_t sum1 = hist1->Integral(); + if (sum1 != 0.0) hist1->Scale( 1.0/sum1 ); + + for (Int_t k=1; k<=nbinsn; k++) + { + Double_t cont = hist1->GetBinContent(k); + fHBlindPixNTheta->SetBinContent(j, k, cont); + } + + delete hist1; + delete hist1OFF; + + //------------------------------------------------------------------ + // define target distribution 'id of blind pixel' + ay = blindidthon.GetYaxis(); + Int_t nbinsid = ay->GetNbins(); + + hist1 = fHBlindPixIdThetaON->ProjectionY("", j, j, ""); + hist1->SetName("dummy"); + hist1OFF = fHBlindPixIdThetaOFF->ProjectionY("", j, j, ""); + + // divide by the number of events (from fHBlindPixNTheta) + if (normON != 0.0) hist1->Scale(1.0/normON); + if (normOFF != 0.0) hist1OFF->Scale(1.0/normOFF); + + // sum of ON and OFF distributions + hist1->Add(hist1OFF, 1.0); + + for (Int_t k=1; k<=nbinsid; k++) + { + Double_t cont = hist1->GetBinContent(k); + fHBlindPixIdTheta->SetBinContent(j, k, cont); + } + + delete hist1; + delete hist1OFF; + } + //............ end of new loop over Theta bins .................... + + *fLog << all + << "MPad::MergeONOFFMC(); The target distributions for the padding have been created" + << endl; + *fLog << all << "----------------------------------------------------------" + << endl; + //-------------------------------------------- + + fHSigmaThetaMC->SetDirectory(NULL); + fHSigmaThetaON->SetDirectory(NULL); + fHSigmaThetaOFF->SetDirectory(NULL); + + fHDiffPixThetaMC->SetDirectory(NULL); + fHDiffPixThetaON->SetDirectory(NULL); + fHDiffPixThetaOFF->SetDirectory(NULL); + + fHSigmaPixThetaMC->SetDirectory(NULL); + fHSigmaPixThetaON->SetDirectory(NULL); + fHSigmaPixThetaOFF->SetDirectory(NULL); + + fHBlindPixIdThetaMC->SetDirectory(NULL); + fHBlindPixIdThetaON->SetDirectory(NULL); + fHBlindPixIdThetaOFF->SetDirectory(NULL); + + fHBlindPixNThetaMC->SetDirectory(NULL); + fHBlindPixNThetaON->SetDirectory(NULL); + fHBlindPixNThetaOFF->SetDirectory(NULL); + + fHgMC->SetDirectory(NULL); + fHgON->SetDirectory(NULL); + fHgOFF->SetDirectory(NULL); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Update the matrix fHgA +// which contains the final padding prescriptions +// + +Bool_t MPad::UpdateHg(TH2D *fHga, TH1D *histap, TH2D *fHge, TH3D *fHgA, + Int_t nbinssig, Int_t l) +{ + // histap target distribution after ON-OFF merging + // fHga padding prescription for ON (or OFF) data to get to histap + // fHge padding prescription to get from histap to the target + // distribution after the ON-OFF-MC merging + // fHgA updated padding prescription for ON (or OFF) data to get + // from the original ON (or OFF) histogram to the target + // distribution after the ON-OFF-MC merging + // l Theta bin + + Double_t eps = 1.e-10; + + for (Int_t k=1; k<=nbinssig; k++) + { + Double_t na = fHga->Integral(1, nbinssig, k, k, " "); + Double_t ne = fHge->Integral(k, k, 1, nbinssig, " "); + Double_t nap = histap->GetBinContent(k); + + if (ne <= eps) + { + // go to next k + continue; + } + + Double_t r = nap - na; + + if (r >= ne-eps) + { + for (Int_t j=k+1; j<=nbinssig; j++) + { + Double_t e = fHge->GetBinContent(k,j); + Double_t A = fHgA->GetBinContent(l,k,j); + fHgA->SetBinContent(l, k, j, A + e); + } + // go to next k + continue; + } + + for (Int_t j=k+1; j<=nbinssig; j++) + { + Double_t e = fHge->GetBinContent(k,j); + Double_t A = fHgA->GetBinContent(l,k,j); + fHgA->SetBinContent(l, k, j, A + r*e/ne); + } + + if (na <= eps) + { + // go to next k + continue; + } + + for (Int_t i=1; iGetBinContent(i,k); + Double_t A = fHgA->GetBinContent(l,i,k); + fHgA->SetBinContent(l, i, k, A - (ne-r)*a/na); + } + + for (Int_t i=1; i<=nbinssig; i++) + { + if (i == k) continue; + for (Int_t j=i+1; j<=nbinssig; j++) + { + if (j == k) continue; + Double_t a = fHga->GetBinContent(i,k); + Double_t e = fHge->GetBinContent(k,j); + Double_t A = fHgA->GetBinContent(l,i,j); + fHgA->SetBinContent(l, i, j, A + (ne-r)*a*e/(na*ne) ); + } + } + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Merge the distributions +// fHSigmaTheta 2D-histogram (Theta, sigmabar) +// +// ===> of ON and MC data <============== +// +// and define the matrices fHgMC and fHgON +// +// These matrices tell which fraction of events should be padded +// from bin k of sigmabar to bin j +// + +Bool_t MPad::MergeONMC( + TH2D& sigthmc, TH3D& diffpixthmc, TH3D& sigmapixthmc, + TH2D& blindidthmc, TH2D& blindnthmc, + TH2D& sigthon, TH3D& diffpixthon, TH3D& sigmapixthon, + TH2D& blindidthon, TH2D& blindnthon) +{ + *fLog << all << "----------------------------------------------------------------------------------" << endl; + *fLog << all << "MPad::MergeONMC(); Merge the ON and MC histograms to obtain the target distributions for the padding" + << endl; + + fHSigmaTheta = new TH2D( (TH2D&) sigthon ); + fHSigmaTheta->SetNameTitle("2D-ThetaSigmabar", "2D-ThetaSigmabar (target)"); + + fHDiffPixTheta = new TH3D( (TH3D&) diffpixthon ); + fHDiffPixTheta->SetNameTitle("3D-ThetaPixDiff", "3D-ThetaPixDiff (target)"); + + fHSigmaPixTheta = new TH3D( (TH3D&) sigmapixthon ); + fHSigmaPixTheta->SetNameTitle("3D-ThetaPixSigma", "3D-ThetaPixSigma (target)"); + + fHBlindPixNTheta = new TH2D( (TH2D&) blindnthon ); + fHBlindPixNTheta->SetNameTitle("2D-ThetaBlindN", "2D-ThetaBlindN (target)"); + + fHBlindPixIdTheta = new TH2D( (TH2D&) blindidthon ); + fHBlindPixIdTheta->SetNameTitle("2D-ThetaBlindId", "2D-ThetaBlindId (target)"); + + //-------------------------- + fHSigmaThetaMC = &sigthmc; + fHSigmaThetaMC->SetNameTitle("2D-ThetaSigmabarMC", "2D-ThetaSigmabarMC (orig.)"); + fHSigmaThetaON = &sigthon; + fHSigmaThetaON->SetNameTitle("2D-ThetaSigmabarON", "2D-ThetaSigmabarON (orig.)"); + + //-------------------------- + fHBlindPixNThetaMC = &blindnthmc; + fHBlindPixNThetaMC->SetNameTitle("2D-ThetaBlindNMC", "2D-ThetaBlindNMC (orig.)"); + fHBlindPixNThetaON = &blindnthon; + fHBlindPixNThetaON->SetNameTitle("2D-ThetaBlindNON", "2D-ThetaBlindNON (orig.)"); + + //-------------------------- + fHBlindPixIdThetaMC = &blindidthmc; + fHBlindPixIdThetaMC->SetNameTitle("2D-ThetaBlindIdMC", "2D-ThetaBlindIdMC (orig.)"); + fHBlindPixIdThetaON = &blindidthon; + fHBlindPixIdThetaON->SetNameTitle("2D-ThetaBlindIdON", "2D-ThetaBlindIdON (orig.)"); + + //-------------------------- + fHDiffPixThetaMC = &diffpixthmc; + fHDiffPixThetaMC->SetNameTitle("3D-ThetaPixDiffMC", "3D-ThetaPixDiffMC (orig.)"); + fHDiffPixThetaON = &diffpixthon; + fHDiffPixThetaON->SetNameTitle("3D-ThetaPixDiffON", "3D-ThetaPixDiffON (orig.)"); + + //-------------------------- + fHSigmaPixThetaMC = &sigmapixthmc; + fHSigmaPixThetaMC->SetNameTitle("3D-ThetaPixSigmaMC", "3D-ThetaPixSigmaMC (orig.)"); + fHSigmaPixThetaON = &sigmapixthon; + fHSigmaPixThetaON->SetNameTitle("3D-ThetaPixSigmaON", "3D-ThetaPixSigmaON (orig.)"); + + //-------------------------- + + + // get binning for fHgON and fHgMC from sigthon + // binning in Theta + TAxis *ax = sigthon.GetXaxis(); Int_t nbinstheta = ax->GetNbins(); TArrayD edgesx; @@ -441,5 +1073,5 @@ // binning in sigmabar - TAxis *ay = sigthon->GetYaxis(); + TAxis *ay = sigthon.GetYaxis(); Int_t nbinssig = ay->GetNbins(); TArrayD edgesy; @@ -463,411 +1095,4 @@ fHgON->SetNameTitle("3D-PaddingMatrixON", "3D-PadMatrixThetaON"); - fHgOFF = new TH3D; - MH::SetBinning(fHgOFF, &binth, &binsg, &binsg); - fHgOFF->SetNameTitle("3D-PaddingMatrixOFF", "3D-PadMatrixThetaOFF"); - - //............ loop over Theta bins ........................... - - - - *fLog << all << "MPad::MergeONOFFMC(); bins of Theta, Sigmabarold, Sigmabarnew, fraction of events to be padded" << endl; - for (Int_t l=1; l<=nbinstheta; l++) - { - TH2D * fHga = new TH2D; - MH::SetBinning(fHga, &binsg, &binsg); - TH2D * fHgb = new TH2D; - MH::SetBinning(fHgb, &binsg, &binsg); - - //------------------------------------------- - // merge ON and OFF distributions - // input : hista is the normalized 1D distr. of sigmabar for ON data - // histb is the normalized 1D distr. of sigmabar for OFF data - // output : histap will be the merged distribution (ON-OFF) - // fHga(k,j) will tell which fraction of ON events should be - // padded from bin k to bin j - // fHgb(k,j) will tell which fraction of OFF events should be - // padded from bin k to bin j - - TH1D *hista = sigthon->ProjectionY("sigon_y", l, l, ""); - Stat_t suma = hista->Integral(); - hista->Scale(1./suma); - - TH1D *histap = new TH1D( (const TH1D&)*hista ); - - TH1D *histb = sigthoff->ProjectionY("sigoff_y", l, l, ""); - Stat_t sumb = histb->Integral(); - histb->Scale(1./sumb); - - Merge2Distributions(hista, histb, histap, fHga, fHgb, nbinssig); - delete hista; - delete histb; - - // fill fHgON and fHgOFF - for (Int_t k=1; k<=nbinssig; k++) - for (Int_t j=1; j<=nbinssig; j++) - { - Double_t a = fHga->GetBinContent(k,j); - fHga->SetBinContent (k, j, 0.0); - fHgON->SetBinContent(l, k, j, a); - - Double_t b = fHgb->GetBinContent(k,j); - fHgb->SetBinContent (k, j, 0.0); - fHgOFF->SetBinContent(l, k, j, b); - } - - //------------------------------------------- - // now merge ON-OFF and MC distributions - // input : hista is the result of the merging of ON and OFF distributions - // histb is the normalized 1D distr. of sigmabar for MC data - // output : histap will be the merged distribution (target distribution) - // fHga(k,j) will tell which fraction of ON, OFF events should be - // padded from bin k to bin j - // fHgb(k,j) will tell which fraction of MC events should be - // padded from bin k to bin j - - hista = new TH1D( (const TH1D&)*histap); - - histb = sigthmc->ProjectionY("sigmc_y", l, l, ""); - sumb = histb->Integral(); - histb->Scale(1./sumb); - - Merge2Distributions(hista, histb, histap, fHga, fHgb, nbinssig); - delete hista; - delete histb; - - // update fHgON and fHgOFF and fill fHgMC - for (Int_t k=1; k<=nbinssig; k++) - for (Int_t j=1; j<=nbinssig; j++) - { - Double_t a = fHga->GetBinContent (k,j); - - Double_t aON = fHgON->GetBinContent(l,k,j); - fHgON->SetBinContent(l, k, j, aON+a); - - Double_t aOFF = fHgOFF->GetBinContent(l,k,j); - fHgOFF->SetBinContent(l, k, j, aOFF+a); - - Double_t b = fHgb->GetBinContent(k,j); - fHgMC->SetBinContent(l, k, j, b); - } - - // fill target distribution fHSigmaTheta - // (only for plotting, it will not be used in the padding) - for (Int_t j=1; j<=nbinssig; j++) - { - Double_t a = histap->GetBinContent(j); - fHSigmaTheta->SetBinContent(l, j, a); - } - - delete fHga; - delete fHgb; - - delete histap; - } - //............ end of loop over Theta bins .................... - - - // Define the target distributions - // fHDiffPixTheta, fHSigmaPixTheta - // (they are calculated as - // averages of the ON and OFF distributions) - // fHBlindPixNTheta, fHBlindPixIdTheta - // (they are calculated as - // the sum of the ON and OFF distributions) - // (fHDiffPixTheta will be used in the padding of MC events) - - //............ start of new loop over Theta bins .................... - for (Int_t j=1; j<=nbinstheta; j++) - { - // fraction of ON/OFF/MC events to be padded : fracON, fracOFF, fracMC - - Double_t fracON = fHgON->Integral (j, j, 1, nbinssig, 1, nbinssig, ""); - Double_t fracOFF = fHgOFF->Integral(j, j, 1, nbinssig, 1, nbinssig, ""); - Double_t fracMC = fHgMC->Integral(j, j, 1, nbinssig, 1, nbinssig, ""); - - Double_t normON; - Double_t normOFF; - - // set ranges for 2D-projections of 3D-histograms - TAxis *ax = diffpixthon->GetXaxis(); - ax->SetRange(j, j); - ax = diffpixthoff->GetXaxis(); - ax->SetRange(j, j); - - TH2D *hist; - TH2D *histOFF; - - TH1D *hist1; - TH1D *hist1OFF; - - // weights for ON and OFF distrubtions when - // calculating the weighted average - Double_t facON = 0.5 * (1. - fracON + fracOFF); - Double_t facOFF = 0.5 * (1. + fracON - fracOFF); - - *fLog << all << "Theta bin j : fracON, fracOFF, facON, facOFF = " - << j << ", " << fracON << ", " << fracOFF << ", " - << fracMC << ", " << facON << ", " << facOFF << endl; - - //------------------------------------------------------------------ - // define target distribution 'sigma^2-sigmavar^2 vs. pixel number' - ay = diffpixthon->GetYaxis(); - Int_t nbinspixel = ay->GetNbins(); - - TAxis *az = diffpixthon->GetZaxis(); - Int_t nbinsdiff = az->GetNbins(); - - hist = (TH2D*)diffpixthon->Project3D("zy"); - hist->SetName("dummy"); - histOFF = (TH2D*)diffpixthoff->Project3D("zy"); - - normON = hist->Integral(); - normOFF = histOFF->Integral(); - if (normON != 0.0) hist->Scale(1.0/normON); - if (normOFF != 0.0) histOFF->Scale(1.0/normOFF); - - // weighted average of ON and OFF distributions - hist->Scale(facON); - hist->Add(histOFF, facOFF); - - for (Int_t k=1; k<=nbinspixel; k++) - for (Int_t l=1; l<=nbinsdiff; l++) - { - Double_t cont = hist->GetBinContent(k,l); - fHDiffPixTheta->SetBinContent(j, k, l, cont); - } - - delete hist; - delete histOFF; - - //------------------------------------------------------------------ - // define target distribution 'sigma vs. pixel number' - ay = sigmapixthon->GetYaxis(); - nbinspixel = ay->GetNbins(); - - az = sigmapixthon->GetZaxis(); - Int_t nbinssigma = az->GetNbins(); - - hist = (TH2D*)sigmapixthon->Project3D("zy"); - hist->SetName("dummy"); - histOFF = (TH2D*)sigmapixthoff->Project3D("zy"); - - normON = hist->Integral(); - normOFF = histOFF->Integral(); - if (normON != 0.0) hist->Scale(1.0/normON); - if (normOFF != 0.0) histOFF->Scale(1.0/normOFF); - - // weighted average of ON and OFF distributions - hist->Scale(facON); - hist->Add(histOFF, facOFF); - - for (Int_t k=1; k<=nbinspixel; k++) - for (Int_t l=1; l<=nbinssigma; l++) - { - Double_t cont = hist->GetBinContent(k,l); - fHSigmaPixTheta->SetBinContent(j, k, l, cont); - } - - delete hist; - delete histOFF; - - - //------------------------------------------------------------------ - // define target distribution 'number of blind pixels per event' - ay = blindnthon->GetYaxis(); - Int_t nbinsn = ay->GetNbins(); - - hist1 = fHBlindPixNThetaON->ProjectionY("", j, j, ""); - hist1->SetName("dummy"); - hist1OFF = fHBlindPixNThetaOFF->ProjectionY("", j, j, ""); - - normON = hist1->Integral(); - normOFF = hist1OFF->Integral(); - if (normON != 0.0) hist1->Scale(1.0/normON); - if (normOFF != 0.0) hist1OFF->Scale(1.0/normOFF); - - // sum of ON and OFF distributions - hist1->Add(hist1OFF, 1.0); - hist1->Scale( 1.0/hist1->Integral() ); - - for (Int_t k=1; k<=nbinsn; k++) - { - Double_t cont = hist1->GetBinContent(k); - fHBlindPixNTheta->SetBinContent(j, k, cont); - } - - delete hist1; - delete hist1OFF; - - //------------------------------------------------------------------ - // define target distribution 'id of blind pixel' - ay = blindidthon->GetYaxis(); - Int_t nbinsid = ay->GetNbins(); - - hist1 = fHBlindPixIdThetaON->ProjectionY("", j, j, ""); - hist1->SetName("dummy"); - hist1OFF = fHBlindPixIdThetaOFF->ProjectionY("", j, j, ""); - - // divide by the number of events (from fHBlindPixNTheta) - if (normON != 0.0) hist1->Scale(1.0/normON); - if (normOFF != 0.0) hist1OFF->Scale(1.0/normOFF); - - // sum of ON and OFF distributions - hist1->Add(hist1OFF, 1.0); - - for (Int_t k=1; k<=nbinsid; k++) - { - Double_t cont = hist1->GetBinContent(k); - fHBlindPixIdTheta->SetBinContent(j, k, cont); - } - - delete hist1; - delete hist1OFF; - } - //............ end of new loop over Theta bins .................... - - *fLog << all - << "MPad::MergeONOFFMC(); The target distributions for the padding have been created" - << endl; - *fLog << all << "----------------------------------------------------------" - << endl; - //-------------------------------------------- - - fHSigmaThetaMC->SetDirectory(NULL); - fHSigmaThetaON->SetDirectory(NULL); - fHSigmaThetaOFF->SetDirectory(NULL); - - fHDiffPixThetaMC->SetDirectory(NULL); - fHDiffPixThetaON->SetDirectory(NULL); - fHDiffPixThetaOFF->SetDirectory(NULL); - - fHSigmaPixThetaMC->SetDirectory(NULL); - fHSigmaPixThetaON->SetDirectory(NULL); - fHSigmaPixThetaOFF->SetDirectory(NULL); - - fHBlindPixIdThetaMC->SetDirectory(NULL); - fHBlindPixIdThetaON->SetDirectory(NULL); - fHBlindPixIdThetaOFF->SetDirectory(NULL); - - fHBlindPixNThetaMC->SetDirectory(NULL); - fHBlindPixNThetaON->SetDirectory(NULL); - fHBlindPixNThetaOFF->SetDirectory(NULL); - - fHgMC->SetDirectory(NULL); - fHgON->SetDirectory(NULL); - fHgOFF->SetDirectory(NULL); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Merge the distributions -// fHSigmaTheta 2D-histogram (Theta, sigmabar) -// -// ===> of ON and MC data <============== -// -// and define the matrices fHgMC and fHgON -// -// These matrices tell which fraction of events should be padded -// from bin k of sigmabar to bin j -// - -Bool_t MPad::MergeONMC( - TH2D *sigthmc, TH3D *diffpixthmc, TH3D *sigmapixthmc, - TH2D *blindidthmc, TH2D *blindnthmc, - TH2D *sigthon, TH3D *diffpixthon, TH3D *sigmapixthon, - TH2D *blindidthon, TH2D *blindnthon) -{ - *fLog << all << "----------------------------------------------------------------------------------" << endl; - *fLog << all << "MPad::MergeONMC(); Merge the ON and MC histograms to obtain the target distributions for the padding" - << endl; - - fHSigmaTheta = new TH2D( (TH2D&)*sigthon ); - fHSigmaTheta->SetNameTitle("2D-ThetaSigmabar", "2D-ThetaSigmabar (target)"); - - fHDiffPixTheta = new TH3D( (TH3D&)*diffpixthon ); - fHDiffPixTheta->SetNameTitle("3D-ThetaPixDiff", "3D-ThetaPixDiff (target)"); - - fHSigmaPixTheta = new TH3D( (TH3D&)*sigmapixthon ); - fHSigmaPixTheta->SetNameTitle("3D-ThetaPixSigma", "3D-ThetaPixSigma (target)"); - - fHBlindPixNTheta = new TH2D( (TH2D&)*blindnthon ); - fHBlindPixNTheta->SetNameTitle("2D-ThetaBlindN", "2D-ThetaBlindN (target)"); - - fHBlindPixIdTheta = new TH2D( (TH2D&)*blindidthon ); - fHBlindPixIdTheta->SetNameTitle("2D-ThetaBlindId", "2D-ThetaBlindId (target)"); - - //-------------------------- - fHSigmaThetaMC = sigthmc; - fHSigmaThetaMC->SetNameTitle("2D-ThetaSigmabarMC", "2D-ThetaSigmabarMC (orig.)"); - fHSigmaThetaON = sigthon; - fHSigmaThetaON->SetNameTitle("2D-ThetaSigmabarON", "2D-ThetaSigmabarON (orig.)"); - - //-------------------------- - fHBlindPixNThetaMC = blindnthmc; - fHBlindPixNThetaMC->SetNameTitle("2D-ThetaBlindNMC", "2D-ThetaBlindNMC (orig.)"); - fHBlindPixNThetaON = blindnthon; - fHBlindPixNThetaON->SetNameTitle("2D-ThetaBlindNON", "2D-ThetaBlindNON (orig.)"); - - //-------------------------- - fHBlindPixIdThetaMC = blindidthmc; - fHBlindPixIdThetaMC->SetNameTitle("2D-ThetaBlindIdMC", "2D-ThetaBlindIdMC (orig.)"); - fHBlindPixIdThetaON = blindidthon; - fHBlindPixIdThetaON->SetNameTitle("2D-ThetaBlindIdON", "2D-ThetaBlindIdON (orig.)"); - - //-------------------------- - fHDiffPixThetaMC = diffpixthmc; - fHDiffPixThetaMC->SetNameTitle("3D-ThetaPixDiffMC", "3D-ThetaPixDiffMC (orig.)"); - fHDiffPixThetaON = diffpixthon; - fHDiffPixThetaON->SetNameTitle("3D-ThetaPixDiffON", "3D-ThetaPixDiffON (orig.)"); - - //-------------------------- - fHSigmaPixThetaMC = sigmapixthmc; - fHSigmaPixThetaMC->SetNameTitle("3D-ThetaPixSigmaMC", "3D-ThetaPixSigmaMC (orig.)"); - fHSigmaPixThetaON = sigmapixthon; - fHSigmaPixThetaON->SetNameTitle("3D-ThetaPixSigmaON", "3D-ThetaPixSigmaON (orig.)"); - - //-------------------------- - - - // get binning for fHgON and fHgMC from sigthon - // binning in Theta - TAxis *ax = sigthon->GetXaxis(); - Int_t nbinstheta = ax->GetNbins(); - TArrayD edgesx; - edgesx.Set(nbinstheta+1); - for (Int_t i=0; i<=nbinstheta; i++) - { - edgesx[i] = ax->GetBinLowEdge(i+1); - *fLog << all << "i, theta low edge = " << i << ", " << edgesx[i] - << endl; - } - MBinning binth; - binth.SetEdges(edgesx); - - // binning in sigmabar - TAxis *ay = sigthon->GetYaxis(); - Int_t nbinssig = ay->GetNbins(); - TArrayD edgesy; - edgesy.Set(nbinssig+1); - for (Int_t i=0; i<=nbinssig; i++) - { - edgesy[i] = ay->GetBinLowEdge(i+1); - *fLog << all << "i, sigmabar low edge = " << i << ", " << edgesy[i] - << endl; - } - MBinning binsg; - binsg.SetEdges(edgesy); - - - fHgMC = new TH3D; - MH::SetBinning(fHgMC, &binth, &binsg, &binsg); - fHgMC->SetNameTitle("3D-PaddingMatrixMC", "3D-PadMatrixThetaMC"); - - fHgON = new TH3D; - MH::SetBinning(fHgON, &binth, &binsg, &binsg); - fHgON->SetNameTitle("3D-PaddingMatrixON", "3D-PadMatrixThetaON"); - //............ loop over Theta bins ........................... @@ -893,13 +1118,13 @@ // padded from bin k to bin j - TH1D *hista = sigthon->ProjectionY("sigon_y", l, l, ""); + TH1D *hista = sigthon.ProjectionY("sigon_y", l, l, ""); Stat_t suma = hista->Integral(); - hista->Scale(1./suma); + if (suma != 0.0) hista->Scale(1./suma); TH1D *histap = new TH1D( (const TH1D&)*hista ); - TH1D *histb = sigthmc->ProjectionY("sigmc_y", l, l, ""); + TH1D *histb = sigthmc.ProjectionY("sigmc_y", l, l, ""); Stat_t sumb = histb->Integral(); - histb->Scale(1./sumb); + if (sumb != 0.0) histb->Scale(1./sumb); Merge2Distributions(hista, histb, histap, fHga, fHgb, nbinssig); @@ -958,7 +1183,7 @@ // set ranges for 2D-projections of 3D-histograms - TAxis *ax = diffpixthon->GetXaxis(); + TAxis *ax = diffpixthon.GetXaxis(); ax->SetRange(j, j); - ax = diffpixthmc->GetXaxis(); + ax = diffpixthmc.GetXaxis(); ax->SetRange(j, j); @@ -980,13 +1205,13 @@ //------------------------------------------------------------------ // define target distribution 'sigma^2-sigmavar^2 vs. pixel number' - ay = diffpixthon->GetYaxis(); + ay = diffpixthon.GetYaxis(); Int_t nbinspixel = ay->GetNbins(); - TAxis *az = diffpixthon->GetZaxis(); + TAxis *az = diffpixthon.GetZaxis(); Int_t nbinsdiff = az->GetNbins(); - hist = (TH2D*)diffpixthon->Project3D("zy"); + hist = (TH2D*)diffpixthon.Project3D("zy"); hist->SetName("dummy"); - histMC = (TH2D*)diffpixthmc->Project3D("zy"); + histMC = (TH2D*)diffpixthmc.Project3D("zy"); normON = hist->Integral(); @@ -1011,18 +1236,18 @@ //------------------------------------------------------------------ // define target distribution 'sigma vs. pixel number' - ay = sigmapixthon->GetYaxis(); + ay = sigmapixthon.GetYaxis(); nbinspixel = ay->GetNbins(); - az = sigmapixthon->GetZaxis(); + az = sigmapixthon.GetZaxis(); Int_t nbinssigma = az->GetNbins(); - hist = (TH2D*)sigmapixthon->Project3D("zy"); + hist = (TH2D*)sigmapixthon.Project3D("zy"); hist->SetName("dummy"); - histMC = (TH2D*)sigmapixthmc->Project3D("zy"); + histMC = (TH2D*)sigmapixthmc.Project3D("zy"); normON = hist->Integral(); normMC = histMC->Integral(); if (normON != 0.0) hist->Scale(1.0/normON); - if (normMC != 0.0) histMC->Scale(1.0/normMC); + if (normMC != 0.0) histMC->Scale(1.0/normMC); // weighted average of ON and MC distributions @@ -1043,5 +1268,5 @@ //------------------------------------------------------------------ // define target distribution 'number of blind pixels per event' - ay = blindnthon->GetYaxis(); + ay = blindnthon.GetYaxis(); Int_t nbinsn = ay->GetNbins(); @@ -1053,9 +1278,10 @@ normMC = hist1MC->Integral(); if (normON != 0.0) hist1->Scale(1.0/normON); - if (normMC != 0.0) hist1MC->Scale(1.0/normMC); + if (normMC != 0.0) hist1MC->Scale(1.0/normMC); // sum of ON and MC distributions hist1->Add(hist1MC, 1.0); - hist1->Scale( 1.0/hist1->Integral() ); + Stat_t sum1 = hist1->Integral(); + if (sum1 != 0.0) hist1->Scale( 1.0/sum1 ); for (Int_t k=1; k<=nbinsn; k++) @@ -1070,5 +1296,5 @@ //------------------------------------------------------------------ // define target distribution 'id of blind pixel' - ay = blindidthon->GetYaxis(); + ay = blindidthon.GetYaxis(); Int_t nbinsid = ay->GetNbins(); @@ -1079,5 +1305,5 @@ // divide by the number of events (from fHBlindPixNTheta) if (normON != 0.0) hist1->Scale(1.0/normON); - if (normMC != 0.0) hist1MC->Scale(1.0/normMC); + if (normMC != 0.0) hist1MC->Scale(1.0/normMC); // sum of ON and MC distributions @@ -1138,4 +1364,6 @@ //------------------------------------ + + /* fHSigmaThetaMC = (TH2D*) gROOT->FindObject("2D-ThetaSigmabarMC"); @@ -1147,7 +1375,11 @@ return kFALSE; } + */ + fHSigmaThetaMC = new TH2D; + fHSigmaThetaMC->Read("2D-ThetaSigmabarMC"); *fLog << all << "MPad : Object '2D-ThetaSigmabarMC' was read in" << endl; + /* fHSigmaThetaON = (TH2D*) gROOT->FindObject("2D-ThetaSigmabarON"); @@ -1159,7 +1391,11 @@ return kFALSE; } + */ + fHSigmaThetaON = new TH2D; + fHSigmaThetaON->Read("2D-ThetaSigmabarON"); *fLog << all << "MPad : Object '2D-ThetaSigmabarON' was read in" << endl; + /* fHSigmaThetaOFF = (TH2D*) gROOT->FindObject("2D-ThetaSigmabarOFF"); @@ -1171,8 +1407,12 @@ return kFALSE; } + */ + fHSigmaThetaOFF = new TH2D; + fHSigmaThetaOFF->Read("2D-ThetaSigmabarOFF"); *fLog << all << "MPad:Object '2D-ThetaSigmabarOFF' was read in" << endl; //------------------------------------ + /* fHDiffPixTheta = (TH3D*) gROOT->FindObject("3D-ThetaPixDiff"); @@ -1184,7 +1424,11 @@ return kFALSE; } + */ + fHDiffPixTheta = new TH3D; + fHDiffPixTheta->Read("3D-ThetaPixDiff"); *fLog << all << "MPad : Object '3D-ThetaPixDiff' was read in" << endl; + /* fHDiffPixThetaMC = (TH3D*) gROOT->FindObject("3D-ThetaPixDiffMC"); @@ -1196,7 +1440,11 @@ return kFALSE; } + */ + fHDiffPixThetaMC = new TH3D; + fHDiffPixThetaMC->Read("3D-ThetaPixDiffMC"); *fLog << all << "MPad : Object '3D-ThetaPixDiffMC' was read in" << endl; + /* fHDiffPixThetaON = (TH3D*) gROOT->FindObject("3D-ThetaPixDiffON"); @@ -1208,7 +1456,11 @@ return kFALSE; } + */ + fHDiffPixThetaON = new TH3D; + fHDiffPixThetaON->Read("3D-ThetaPixDiffON"); *fLog << all << "MPad : Object '3D-ThetaPixDiffON' was read in" << endl; + /* fHDiffPixThetaOFF = (TH3D*) gROOT->FindObject("3D-ThetaPixDiffOFF"); @@ -1220,4 +1472,7 @@ return kFALSE; } + */ + fHDiffPixThetaOFF = new TH3D; + fHDiffPixThetaOFF->Read("3D-ThetaPixDiffOFF"); *fLog << all << "MPad : Object '3D-ThetaPixDiffOFF' was read in" << endl; @@ -1225,4 +1480,5 @@ //------------------------------------ + /* fHSigmaPixTheta = (TH3D*) gROOT->FindObject("3D-ThetaPixSigma"); @@ -1234,7 +1490,11 @@ return kFALSE; } + */ + fHSigmaPixTheta = new TH3D; + fHSigmaPixTheta->Read("3D-ThetaPixSigma"); *fLog << all << "MPad : Object '3D-ThetaPixSigma' was read in" << endl; + /* fHSigmaPixThetaMC = (TH3D*) gROOT->FindObject("3D-ThetaPixSigmaMC"); @@ -1246,7 +1506,11 @@ return kFALSE; } + */ + fHSigmaPixThetaMC = new TH3D; + fHSigmaPixThetaMC->Read("3D-ThetaPixSigmaMC"); *fLog << all << "MPad : Object '3D-ThetaPixSigmaMC' was read in" << endl; + /* fHSigmaPixThetaON = (TH3D*) gROOT->FindObject("3D-ThetaPixSigmaON"); @@ -1258,7 +1522,11 @@ return kFALSE; } + */ + fHSigmaPixThetaON = new TH3D; + fHSigmaPixThetaON->Read("3D-ThetaPixSigmaON"); *fLog << all << "MPad : Object '3D-ThetaPixSigmaON' was read in" << endl; + /* fHSigmaPixThetaOFF = (TH3D*) gROOT->FindObject("3D-ThetaPixSigmaOFF"); @@ -1270,8 +1538,12 @@ return kFALSE; } + */ + fHSigmaPixThetaOFF = new TH3D; + fHSigmaPixThetaOFF->Read("3D-ThetaPixSigmaOFF"); *fLog << all << "MPad : Object '3D-ThetaPixSigmaOFF' was read in" << endl; //------------------------------------ + /* fHgMC = (TH3D*) gROOT->FindObject("3D-PaddingMatrixMC"); @@ -1283,7 +1555,11 @@ return kFALSE; } + */ + fHgMC = new TH3D; + fHgMC->Read("3D-PaddingMatrixMC"); *fLog << all << "MPad : Object '3D-PaddingMatrixMC' was read in" << endl; + /* fHgON = (TH3D*) gROOT->FindObject("3D-PaddingMatrixON"); @@ -1295,7 +1571,11 @@ return kFALSE; } + */ + fHgON = new TH3D; + fHgON->Read("3D-PaddingMatrixON"); *fLog << all << "MPad : Object '3D-PaddingMatrixON' was read in" << endl; + /* fHgOFF = (TH3D*) gROOT->FindObject("3D-PaddingMatrixOFF"); @@ -1307,8 +1587,28 @@ return kFALSE; } + */ + fHgOFF = new TH3D; + fHgOFF->Read("3D-PaddingMatrixOFF"); *fLog << all << "MPad : Object '3D-PaddingMatrixOFF' was read in" << endl; //------------------------------------ + /* + fHBlindPixIdTheta = + (TH2D*) gROOT->FindObject("2D-ThetaBlindId"); + if (!fHBlindPixIdTheta) + { + *fLog << all + << "MPad : Object '2D-ThetaBlindId' not found on root file" + << endl; + return kFALSE; + } + */ + fHBlindPixIdTheta = new TH2D; + fHBlindPixIdTheta->Read("2D-ThetaBlindId"); + *fLog << all + << "MPad : Object '2D-ThetaBlindId' was read in" << endl; + + /* fHBlindPixIdThetaMC = (TH2D*) gROOT->FindObject("2D-ThetaBlindIdMC"); @@ -1320,7 +1620,11 @@ return kFALSE; } + */ + fHBlindPixIdThetaMC = new TH2D; + fHBlindPixIdThetaMC->Read("2D-ThetaBlindIdMC"); *fLog << all << "MPad : Object '2D-ThetaBlindIdMC' was read in" << endl; + /* fHBlindPixIdThetaON = (TH2D*) gROOT->FindObject("2D-ThetaBlindIdON"); @@ -1332,7 +1636,11 @@ return kFALSE; } + */ + fHBlindPixIdThetaON = new TH2D; + fHBlindPixIdThetaON->Read("2D-ThetaBlindIdON"); *fLog << all << "MPad : Object '2D-ThetaBlindIdON' was read in" << endl; + /* fHBlindPixIdThetaOFF = (TH2D*) gROOT->FindObject("2D-ThetaBlindIdOFF"); @@ -1344,8 +1652,28 @@ return kFALSE; } + */ + fHBlindPixIdThetaOFF = new TH2D; + fHBlindPixIdThetaOFF->Read("2D-ThetaBlindIdOFF"); *fLog << all - << "MPad : Object '2D-ThetaBlindIdMC' was read in" << endl; + << "MPad : Object '2D-ThetaBlindIdOFF' was read in" << endl; //------------------------------------ + /* + fHBlindPixNTheta = + (TH2D*) gROOT->FindObject("2D-ThetaBlindN"); + if (!fHBlindPixNTheta) + { + *fLog << all + << "MPad : Object '2D-ThetaBlindN' not found on root file" + << endl; + return kFALSE; + } + */ + fHBlindPixNTheta = new TH2D; + fHBlindPixNTheta->Read("2D-ThetaBlindN"); + *fLog << all + << "MPad : Object '2D-ThetaBlindN' was read in" << endl; + + /* fHBlindPixNThetaMC = (TH2D*) gROOT->FindObject("2D-ThetaBlindNMC"); @@ -1357,7 +1685,11 @@ return kFALSE; } + */ + fHBlindPixNThetaMC = new TH2D; + fHBlindPixNThetaMC->Read("2D-ThetaBlindNMC"); *fLog << all << "MPad : Object '2D-ThetaBlindNMC' was read in" << endl; + /* fHBlindPixNThetaON = (TH2D*) gROOT->FindObject("2D-ThetaBlindNON"); @@ -1369,7 +1701,11 @@ return kFALSE; } + */ + fHBlindPixNThetaON = new TH2D; + fHBlindPixNThetaON->Read("2D-ThetaBlindNON"); *fLog << all << "MPad : Object '2D-ThetaBlindNON' was read in" << endl; + /* fHBlindPixNThetaOFF = (TH2D*) gROOT->FindObject("2D-ThetaBlindNOFF"); @@ -1381,4 +1717,7 @@ return kFALSE; } + */ + fHBlindPixNThetaOFF = new TH2D; + fHBlindPixNThetaOFF->Read("2D-ThetaBlindNOFF"); *fLog << all << "MPad : Object '2D-ThetaBlindNOFF' was read in" << endl; @@ -1402,8 +1741,10 @@ TFile outfile(namefileout, "RECREATE"); + fHBlindPixNTheta->Write(); fHBlindPixNThetaMC->Write(); fHBlindPixNThetaON->Write(); fHBlindPixNThetaOFF->Write(); + fHBlindPixIdTheta->Write(); fHBlindPixIdThetaMC->Write(); fHBlindPixIdThetaON->Write(); @@ -1470,11 +1811,11 @@ } - fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); - if (!fMcEvt) + fPointPos = (MPointingPos*)pList->FindObject("MPointingPos"); + if (!fPointPos) { - *fLog << err << dbginf << "MPad : MMcEvt not found... aborting." + *fLog << err << dbginf << "MPad : MPointingPos not found... aborting." << endl; return kFALSE; - } + } fPed = (MPedPhotCam*)pList->FindObject("MPedPhotCam"); @@ -1510,6 +1851,6 @@ } - fBlinds = (MBlindPixels*)pList->FindCreateObj("MBlindPixels"); - if (!fBlinds) + fBlindPix = (MBlindPixels*)pList->FindCreateObj("MBlindPixels"); + if (!fBlindPix) { *fLog << err << "MPad : MBlindPixels not found... aborting." @@ -1532,5 +1873,5 @@ // plot pedestal sigmas fHSigmaPedestal = new TH2D("SigPed","Sigma: after vs. before padding", - 100, 0.0, 5.0, 100, 0.0, 5.0); + 100, 0.0, 75.0, 100, 0.0, 75.0); fHSigmaPedestal->SetXTitle("Pedestal sigma before padding"); fHSigmaPedestal->SetYTitle("Pedestal sigma after padding"); @@ -1538,10 +1879,10 @@ // plot no.of photons (before vs. after padding) fHPhotons = new TH2D("Photons","Photons: after vs.before padding", - 100, -10.0, 90.0, 100, -10, 90); + 100, -100.0, 300.0, 100, -100, 300); fHPhotons->SetXTitle("no.of photons before padding"); fHPhotons->SetYTitle("no.of photons after padding"); // plot of added NSB - fHNSB = new TH1D("NSB","Distribution of added NSB", 100, -10.0, 10.0); + fHNSB = new TH1D("NSB","Distribution of added NSB", 100, -100.0, 200.0); fHNSB->SetXTitle("no.of added NSB photons"); fHNSB->SetYTitle("no.of pixels"); @@ -1550,5 +1891,5 @@ //-------------------------------------------------------------------- - fIter = 20; + fIter = 10; memset(fErrors, 0, sizeof(fErrors)); @@ -1570,5 +1911,5 @@ Int_t MPad::Process() { - *fLog << all << "Entry MPad::Process();" << endl; + //*fLog << all << "Entry MPad::Process();" << endl; // this is the conversion factor from the number of photons @@ -1576,10 +1917,16 @@ // later to be taken from MCalibration // fPEperPhoton = PW * LG * QE * 1D - Double_t fPEperPhoton = 0.95 * 0.90 * 0.13 * 0.9; + Double_t fPEperPhoton = 0.92 * 0.94 * 0.23 * 0.9; //------------------------------------------------ - // add pixels to MCerPhotEvt which are not yet in; - // set their number of photons equal to zero - + // Add pixels to MCerPhotEvt which are not yet in; + // set their number of photons equal to zero. + // Purpose : by the padding the number of photons is changed + // so that cleaning and cuts may be affected. + // However, no big effect is expectec unless the padding is strong + // (strong increase of the noise level) + // At present comment out this code + + /* UInt_t imaxnumpix = fCam->GetNumPixels(); @@ -1603,5 +1950,5 @@ fEvt->AddPixel(i, 0.0, (*fPed)[i].GetRms()); } - + */ @@ -1615,12 +1962,22 @@ // Calculate average sigma of the event // - Double_t sigbarold_phot = fSigmabar->Calc(*fCam, *fPed, *fEvt); - *fLog << all << "MPad::Process(); before padding : " << endl; - fSigmabar->Print(""); - Double_t sigbarold = sigbarold_phot * fPEperPhoton; - Double_t sigbarold2 = sigbarold*sigbarold; - - const Double_t theta = kRad2Deg*fMcEvt->GetTelescopeTheta(); - *fLog << all << "theta = " << theta << endl; + + fSigmabar->Calc(*fCam, *fPed, *fEvt); + //*fLog << all << "MPad::Process(); before padding : " << endl; + //fSigmabar->Print(""); + + // inner sigmabar/sqrt(area) + Double_t sigbarInnerold_phot = fSigmabar->GetSigmabarInner(); + Double_t sigbarInnerold = sigbarInnerold_phot * fPEperPhoton; + Double_t sigbarInnerold2 = sigbarInnerold*sigbarInnerold; + + // outer sigmabar/sqrt(area) + Double_t sigbarOuterold_phot = fSigmabar->GetSigmabarOuter(); + Double_t sigbarOuterold = sigbarOuterold_phot * fPEperPhoton; + Double_t sigbarOuterold2 = sigbarOuterold*sigbarOuterold; + + const Double_t theta = fPointPos->GetZd(); + + //*fLog << all << "theta = " << theta << endl; Int_t binTheta = fHBlindPixNTheta->GetXaxis()->FindBin(theta); @@ -1629,6 +1986,5 @@ *fLog << warn << "MPad::Process(); binNumber out of range : theta, binTheta = " - << theta << ", " << binTheta << "; Skip event " << endl; - // event cannot be padded; skip event + << theta << ", " << binTheta << "; aborting " << endl; rc = 2; @@ -1637,7 +1993,10 @@ } + //*fLog << all << "binTheta = " << binTheta << endl; + + //------------------------------------------- // get number of events in this theta bin - // and number of events in this sigbarold_phot bin + // and number of events in this sigbarInnerold_phot bin Double_t nTheta; @@ -1657,16 +2016,24 @@ TH1D *hn; hn = st->ProjectionY("", binTheta, binTheta, ""); + + //*fLog << "Title of histogram : " << st->GetTitle() << endl; + //for (Int_t i=1; i<=hn->GetNbinsX(); i++) + //{ + // *fLog << "bin " << i << " : no.of entries = " << hn->GetBinContent(i) + // << endl; + //} + nTheta = hn->Integral(); - Int_t binSigma = hn->FindBin(sigbarold_phot); + Int_t binSigma = hn->FindBin(sigbarInnerold_phot); nSigma = hn->GetBinContent(binSigma); - *fLog << all - << "Theta, sigbarold_phot, binTheta, binSigma, nTheta, nSigma = " - << theta << ", " << sigbarold_phot << ", " << binTheta << ", " - << binSigma << ", " << nTheta << ", " << nSigma << endl; + //*fLog << all + // << "Theta, sigbarold_phot, binTheta, binSigma, nTheta, nSigma = " + // << theta << ", " << sigbarInnerold_phot << ", " << binTheta + // << ", " + // << binSigma << ", " << nTheta << ", " << nSigma << endl; delete hn; - - + //------------------------------------------- // for the current theta : @@ -1688,6 +2055,7 @@ if ( nblind->Integral() == 0.0 ) { - *fLog << warn << "MPad::Process(); projection for Theta bin " - << binTheta << " has no entries; Skip event " << endl; + *fLog << warn << "MPad::Process(); projection of '" + << fHBlindPixNThetaOFF->GetName() << "' for Theta bin " + << binTheta << " has no entries; aborting " << endl; // event cannot be padded; skip event delete nblind; @@ -1732,8 +2100,8 @@ nlist++; - fBlinds->SetPixelBlind(idBlind); - *fLog << all << "idBlind = " << idBlind << endl; + fBlindPix->SetPixelBlind(idBlind); + //*fLog << all << "idBlind = " << idBlind << endl; } - fBlinds->SetReadyToSave(); + fBlindPix->SetReadyToSave(); delete hblind; @@ -1752,7 +2120,7 @@ if ( nblind->Integral() == 0.0 ) { - *fLog << warn << "MPad::Process(); projection for Theta bin " + *fLog << warn << "MPad::Process(); projection of '" + << fHBlindPixNThetaON->GetName() << "' for Theta bin " << binTheta << " has no entries; Skip event " << endl; - // event cannot be padded; skip event delete nblind; @@ -1796,8 +2164,8 @@ nlist++; - fBlinds->SetPixelBlind(idBlind); - *fLog << all << "idBlind = " << idBlind << endl; + fBlindPix->SetPixelBlind(idBlind); + //*fLog << all << "idBlind = " << idBlind << endl; } - fBlinds->SetReadyToSave(); + fBlindPix->SetReadyToSave(); delete hblind; @@ -1810,7 +2178,7 @@ // - Int_t binSig = fHgON->GetYaxis()->FindBin(sigbarold_phot); - *fLog << all << "binSig, sigbarold_phot = " << binSig << ", " - << sigbarold_phot << endl; + Int_t binSig = fHgON->GetYaxis()->FindBin(sigbarInnerold_phot); + //*fLog << all << "binSig, sigbarInnerold_phot = " << binSig << ", " + // << sigbarInnerold_phot << endl; Double_t prob; @@ -1843,6 +2211,6 @@ prob = 0.0; - *fLog << all << "nTheta, nSigma, prob = " << nTheta << ", " << nSigma - << ", " << prob << endl; + //*fLog << all << "nTheta, nSigma, prob = " << nTheta << ", " << nSigma + // << ", " << prob << endl; if ( prob <= 0.0 || gRandom->Uniform() > prob ) @@ -1850,5 +2218,5 @@ delete hpad; // event should not be padded - *fLog << all << "event is not padded" << endl; + //*fLog << all << "event is not padded" << endl; rc = 8; @@ -1857,5 +2225,5 @@ } // event should be padded - *fLog << all << "event will be padded" << endl; + //*fLog << all << "event will be padded" << endl; @@ -1864,6 +2232,6 @@ // for MC/ON/OFF data according to the matrix fHgMC/ON/OFF // - Double_t sigmabar_phot = 0; - Double_t sigmabar = 0; + Double_t sigmabarInner_phot = 0; + Double_t sigmabarInner = 0; @@ -1876,23 +2244,25 @@ //} - sigmabar_phot = hpad->GetRandom(); - sigmabar = sigmabar_phot * fPEperPhoton; - - *fLog << all << "sigmabar_phot = " << sigmabar_phot << endl; + sigmabarInner_phot = hpad->GetRandom(); + sigmabarInner = sigmabarInner_phot * fPEperPhoton; + + //*fLog << all << "sigmabarInner_phot = " << sigmabarInner_phot << endl; delete hpad; - const Double_t sigmabar2 = sigmabar*sigmabar; + // new inner sigmabar2/area + const Double_t sigmabarInner2 = sigmabarInner*sigmabarInner; //------------------------------------------- - *fLog << all << "MPad::Process(); sigbarold, sigmabar = " - << sigbarold << ", "<< sigmabar << endl; - - // Skip event if target sigmabar is <= sigbarold - if (sigmabar <= sigbarold) + //*fLog << all << "MPad::Process(); sigbarInnerold, sigmabarInner = " + // << sigbarInnerold << ", "<< sigmabarInner << endl; + + // Stop if target sigmabar is <= sigbarold + if (sigmabarInner <= sigbarInnerold) { - *fLog << all << "MPad::Process(); target sigmabar is less than sigbarold : " - << sigmabar << ", " << sigbarold << ", Skip event" << endl; + *fLog << all << "MPad::Process(); target sigmabarInner is less than sigbarInnerold : " + << sigmabarInner << ", " << sigbarInnerold << ", aborting" + << endl; rc = 4; @@ -1923,16 +2293,21 @@ } - // Get RMS of (Sigma^2-sigmabar^2) in this Theta bin. - // The average electronic noise (to be added) has to be well above this RMS, - // otherwise the electronic noise of an individual pixel (elNoise2Pix) - // may become negative + // In this Theta bin, get RMS of (Sigma^2-sigmabar^2) for inner pixels + // The average electronic noise (to be added) has to be in the order of + // this RMS, otherwise the electronic noise of an individual pixel + // (elNoise2Pix) may become negative + + // Attention : maximum ID of inner pixels hard coded !!! + Int_t idmaxpixinner = 396; + Int_t binmaxpixinner = + fHDiffPixTheta->GetYaxis()->FindBin( (Double_t)idmaxpixinner ); TH1D *hnoise; if (fType == "MC") - hnoise = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, 0, 9999, ""); + hnoise = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, 0, binmaxpixinner, ""); else if (fType == "ON") - hnoise = fHDiffPixThetaOFF->ProjectionZ("", binTheta, binTheta, 0, 9999, ""); + hnoise = fHDiffPixThetaOFF->ProjectionZ("", binTheta, binTheta, 0, binmaxpixinner, ""); else if (fType == "OFF") - hnoise = fHDiffPixThetaON->ProjectionZ("", binTheta, binTheta, 0, 9999, ""); + hnoise = fHDiffPixThetaON->ProjectionZ("", binTheta, binTheta, 0, binmaxpixinner, ""); else { @@ -1946,8 +2321,8 @@ delete hnoise; - elNoise2 = TMath::Min(RMS, sigmabar2 - sigbarold2); - *fLog << all << "elNoise2 = " << elNoise2 << endl; - - lambdabar = (sigmabar2 - sigbarold2 - elNoise2) / F2excess; + elNoise2 = TMath::Min(2.0*RMS, sigmabarInner2 - sigbarInnerold2); + //*fLog << all << "elNoise2 = " << elNoise2 << endl; + + lambdabar = (sigmabarInner2 - sigbarInnerold2 - elNoise2) / F2excess; // This value of lambdabar is the same for all pixels; @@ -1957,5 +2332,5 @@ // do the padding for each pixel // - // pad only pixels - which are used (before image cleaning) + // pad only pixels - which are used and not blind (before image cleaning) // @@ -1971,6 +2346,29 @@ + // throw the Sigma for the pixels from the distribution fHDiffPixTheta + // MC : from fHDiffPixTheta + // ON : from fHDiffPixThetaOFF + // OFF : from fHDiffPixThetaON + + TH3D *sp = NULL; + + if (fType == "MC") sp = fHDiffPixTheta; + else if (fType == "ON") sp = fHDiffPixThetaOFF; + else if (fType == "OFF") sp = fHDiffPixThetaON; + else + { + *fLog << err << "MPad::Process(); illegal data type... aborting" + << endl; + return kERROR; + } + + Double_t sigbarold2; + Double_t sigmabar2; + Double_t sigmabarOuter2; + + for (UInt_t i=0; iGetPixRatio(j); + if (ratioArea < 1.5) + { + sigbarold2 = sigbarInnerold2; + sigmabar2 = sigmabarInner2; + } + else + { + sigbarold2 = sigbarOuterold2; + //sigbarOuter2 = sigmabarInner2 * sigbarOuterold2 / sigbarInnerold2; + sigmabarOuter2 = sigmabarInner2 + (sigbarOuterold2 - sigbarInnerold2); + sigmabar2 = sigmabarOuter2; + } + MPedPhotPix &ppix = (*fPed)[j]; + + if ( fBlindPix != NULL && fBlindPix->IsBlind(j) ) + { + // this should never occur, because blind pixels should have + // been set unused by MBlindPixelsCalc2::UnMap() + //*fLog << all << "MPad::Process; blind pixel found which is used, j = " + // << j << "... go to next pixel." << endl; + continue; + } + + // count number of pixels treated + fWarnings[0]++; + + Double_t oldsigma_phot = ppix.GetRms(); Double_t oldsigma = oldsigma_phot * fPEperPhoton; - Double_t oldsigma2 = oldsigma*oldsigma; + Double_t oldsigma2 = oldsigma*oldsigma / ratioArea; //--------------------------------- @@ -2005,69 +2430,64 @@ TH1D *hdiff; - // throw the Sigma for this pixel from the distribution fHDiffPixTheta - // MC : from fHDiffPixTheta - // ON : from fHDiffPixThetaOFF - // OFF : from fHDiffPixThetaON - - TH3D *sp = NULL; - - if (fType == "MC") sp = fHDiffPixTheta; - else if (fType == "ON") sp = fHDiffPixThetaOFF; - else if (fType == "OFF") sp = fHDiffPixThetaON; - else - { - *fLog << err << "MPad::Process(); illegal data type... aborting" - << endl; - return kERROR; - } - - hdiff = sp->ProjectionZ("", binTheta, binTheta, - binPixel, binPixel, ""); - - if ( hdiff->GetEntries() == 0 ) - { - *fLog << warn << "MPad::Process(); projection for Theta bin " - << binTheta << " and pixel bin " << binPixel - << " has no entries; aborting " << endl; - delete hdiff; - - rc = 5; - fErrors[rc]++; - return kCONTINUE; - } - - count = 0; - ok = kFALSE; - for (Int_t m=0; mGetRandom(); - diff = diff_phot * fPEperPhoton * fPEperPhoton; + + hdiff = sp->ProjectionZ("", binTheta, binTheta, + binPixel, binPixel, ""); + Double_t integral = hdiff->Integral(); + // if there are no entries in hdiff, diff cannot be thrown + // in this case diff will be set to zero + if ( integral == 0 ) + { + //*fLog << warn << "MPad::Process(); fType = " << fType + // << ", projection of '" + // << sp->GetName() << "' for Theta bin " + // << binTheta << " and pixel " << j + // << " has no entries; set diff equal to the old difference " + // << endl; + + diff = TMath::Max(oldsigma2 - sigbarold2, + lambdabar*F2excess + oldsigma2 - sigmabar2); + + rc = 2; + fWarnings[rc]++; + } + + // start of else ------------------- + else + { + count = 0; + ok = kFALSE; + for (Int_t m=0; mGetRandom(); + + //*fLog << "after GetRandom : j, m, diff_phot = " << j << " : " + // << m << ", " << diff_phot << endl; + + diff = diff_phot * fPEperPhoton * fPEperPhoton; - // the following condition ensures that elNoise2Pix > 0.0 - if ( (diff + sigmabar2 - oldsigma2/ratioArea - - lambdabar*F2excess) > 0.0 ) - { - ok = kTRUE; - break; - } - } - - if (!ok) - { - *fLog << all << "theta, j, count, sigmabar, diff = " << theta - << ", " - << j << ", " << count << ", " << sigmabar << ", " - << diff << endl; - diff = lambdabar*F2excess + oldsigma2/ratioArea - sigmabar2; - - rw = 1; - fWarnings[rw]++; - } - else - { - rw = 0; - fWarnings[rw]++; - } + // the following condition ensures that elNoise2Pix > 0.0 + if ( (diff + sigmabar2 - oldsigma2 + - lambdabar*F2excess) > 0.0 ) + { + ok = kTRUE; + break; + } + } + + if (!ok) + { + //*fLog << all << "theta, j, count, sigmabar2, diff, oldsigma2, ratioArea, lambdabar = " + // << theta << ", " + // << j << ", " << count << ", " << sigmabar2 << ", " + // << diff << ", " << oldsigma2 << ", " << ratioArea << ", " + // << lambdabar << endl; + diff = lambdabar*F2excess + oldsigma2 - sigmabar2; + + rw = 1; + fWarnings[rw]++; + } + } + // end of else -------------------- delete hdiff; @@ -2078,5 +2498,5 @@ // get the additional sigma^2 for this pixel (due to the padding) - addSig2 = sigma2*ratioArea - oldsigma2; + addSig2 = (sigma2 - oldsigma2) * ratioArea; addSig2_phot = addSig2 / (fPEperPhoton * fPEperPhoton); @@ -2139,5 +2559,5 @@ - Double_t newsigma = sqrt( oldsigma2 + addSig2 ); + Double_t newsigma = sqrt( sigma2 * ratioArea ); Double_t newsigma_phot = newsigma / fPEperPhoton; ppix.SetRms( newsigma_phot ); @@ -2150,8 +2570,8 @@ fSigmabar->Calc(*fCam, *fPed, *fEvt); - *fLog << all << "MPad::Process(); after padding : " << endl; - fSigmabar->Print(""); - - *fLog << all << "Exit MPad::Process();" << endl; + //*fLog << all << "MPad::Process(); after padding : " << endl; + //fSigmabar->Print(""); + + //*fLog << all << "Exit MPad::Process();" << endl; rc = 0; @@ -2173,18 +2593,15 @@ *fLog << dec << setfill(' '); - int temp; - if (fWarnings[0] != 0.0) - temp = (int)(fWarnings[1]*100/fWarnings[0]); - else - temp = 100; + if (fWarnings[0] == 0 ) fWarnings[0] = 1; *fLog << " " << setw(7) << fWarnings[1] << " (" << setw(3) - << temp - << "%) Warning : iteration to find acceptable sigma failed" + << (int)(fWarnings[1]*100/fWarnings[0]) + << "%) Pixel: iteration to find acceptable sigma failed" << ", fIter = " << fIter << endl; - *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) - << (int)(fErrors[1]*100/GetNumExecutions()) - << "%) Evts skipped due to: Sigmabar_old > 0" << endl; + *fLog << " " << setw(7) << fWarnings[2] << " (" << setw(3) + << (int)(fWarnings[2]*100/fWarnings[0]) + << "%) Pixel: No data for generating Sigma^2-Sigmabar^2" << endl; + *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3) @@ -2192,19 +2609,7 @@ << "%) Evts skipped due to: Zenith angle out of range" << endl; - *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3) - << (int)(fErrors[3]*100/GetNumExecutions()) - << "%) Evts skipped due to: No data for generating Sigmabar" << endl; - *fLog << " " << setw(7) << fErrors[4] << " (" << setw(3) << (int)(fErrors[4]*100/GetNumExecutions()) << "%) Evts skipped due to: Target sigma <= Sigmabar_old" << endl; - - *fLog << " " << setw(7) << fErrors[5] << " (" << setw(3) - << (int)(fErrors[5]*100/GetNumExecutions()) - << "%) Evts skipped due to: No data for generating Sigma^2-Sigmabar^2" << endl; - - *fLog << " " << setw(7) << fErrors[6] << " (" << setw(3) - << (int)(fErrors[6]*100/GetNumExecutions()) - << "%) Evts skipped due to: No data for generating Sigma" << endl; *fLog << " " << setw(7) << fErrors[7] << " (" << setw(3) @@ -2224,6 +2629,6 @@ //--------------------------------------------------------------- - TCanvas &c = *(MH::MakeDefCanvas("Pad", "", 900, 1500)); - c.Divide(3, 5); + TCanvas &c = *(MH::MakeDefCanvas("Pad", "", 900, 600)); + c.Divide(3, 2); gROOT->SetSelectedPad(NULL); @@ -2245,5 +2650,5 @@ //-------------------------------------------------------------------- - + /* c.cd(4); fHSigmaTheta->SetDirectory(NULL); @@ -2251,9 +2656,9 @@ fHSigmaTheta->DrawCopy(); fHSigmaTheta->SetBit(kCanDelete); - + */ //-------------------------------------------------------------------- - + /* c.cd(7); fHBlindPixNTheta->SetDirectory(NULL); @@ -2261,8 +2666,9 @@ fHBlindPixNTheta->DrawCopy(); fHBlindPixNTheta->SetBit(kCanDelete); + */ //-------------------------------------------------------------------- - + /* c.cd(10); fHBlindPixIdTheta->SetDirectory(NULL); @@ -2270,9 +2676,10 @@ fHBlindPixIdTheta->DrawCopy(); fHBlindPixIdTheta->SetBit(kCanDelete); - + */ //-------------------------------------------------------------------- // draw the 3D histogram (target): Theta, pixel, Sigma^2-Sigmabar^2 + /* c.cd(5); TH2D *l1 = NULL; @@ -2296,8 +2703,10 @@ l2->DrawCopy("box"); l2->SetBit(kCanDelete);; + */ //-------------------------------------------------------------------- // draw the 3D histogram (target): Theta, pixel, Sigma + /* c.cd(6); TH2D *k1 = NULL; @@ -2321,13 +2730,14 @@ k2->DrawCopy("box"); k2->SetBit(kCanDelete);; - + */ //-------------------------------------------------------------------- - c.cd(13); + + c.cd(4); TH2D *m1; m1 = (TH2D*) ((TH3*)fHgON)->Project3D("zy"); m1->SetDirectory(NULL); - m1->SetTitle("(Target) Sigmabar-new vs. Sigmabar-old (ON, all \\Theta)"); + m1->SetTitle("(fHgON) Sigmabar-new vs. Sigmabar-old (ON, all \\Theta)"); m1->SetXTitle("Sigmabar-old"); m1->SetYTitle("Sigmabar-new"); @@ -2336,9 +2746,9 @@ m1->SetBit(kCanDelete);; - c.cd(14); + c.cd(5); TH2D *m2; m2 = (TH2D*) ((TH3*)fHgOFF)->Project3D("zy"); m2->SetDirectory(NULL); - m2->SetTitle("(Target) Sigmabar-new vs. Sigmabar-old (OFF, all \\Theta)"); + m2->SetTitle("(fHgOFF) Sigmabar-new vs. Sigmabar-old (OFF, all \\Theta)"); m2->SetXTitle("Sigmabar-old"); m2->SetYTitle("Sigmabar-new"); @@ -2347,9 +2757,9 @@ m2->SetBit(kCanDelete);; - c.cd(15); + c.cd(6); TH2D *m3; m3 = (TH2D*) ((TH3*)fHgMC)->Project3D("zy"); m3->SetDirectory(NULL); - m3->SetTitle("(Target) Sigmabar-new vs. Sigmabar-old (MC, all \\Theta)"); + m3->SetTitle("(fHgMC) Sigmabar-new vs. Sigmabar-old (MC, all \\Theta)"); m3->SetXTitle("Sigmabar-old"); m3->SetYTitle("Sigmabar-new"); @@ -2369,2 +2779,6 @@ + + + + Index: trunk/MagicSoft/Mars/manalysis/MPad.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPad.h (revision 4577) +++ trunk/MagicSoft/Mars/manalysis/MPad.h (revision 4584) @@ -17,5 +17,5 @@ class MCerPhotEvt; class MPedPhotCam; -class MMcEvt; +class MPointingPos; class MSigmabar; class MParList; @@ -31,7 +31,7 @@ MCerPhotEvt *fEvt; MSigmabar *fSigmabar; - MMcEvt *fMcEvt; + MPointingPos *fPointPos; MPedPhotCam *fPed; - MBlindPixels *fBlinds; + MBlindPixels *fBlindPix; TString fType; // type of data to be padded @@ -41,5 +41,5 @@ Int_t fErrors[9]; - Int_t fWarnings[2]; + Int_t fWarnings[3]; //---------------------------------- @@ -89,7 +89,9 @@ - Bool_t Merge2Distributions(TH1D *hista, TH1D * histb, TH1D * histap, - TH2D *fHga, TH2D * fHgb, Int_t nbinssig); + Bool_t Merge2Distributions(TH1D *hista, TH1D *histb, TH1D *histap, + TH2D *fHga, TH2D *fHgb, Int_t nbinssig); + Bool_t UpdateHg(TH2D *fHga, TH1D *histap, TH2D *fHge, TH3D *fHgA, + Int_t nbinssig, Int_t l); public: @@ -98,16 +100,16 @@ Bool_t MergeONOFFMC( - TH2D *sigthmc, TH3D *diffpixthmc, TH3D *sigmapixthmc, - TH2D *blindidthmc, TH2D *blindnthmc, - TH2D *sigthon, TH3D *diffpixthon, TH3D *sigmapixthon, - TH2D *blindidthon, TH2D *blindnthon, - TH2D *sigthoff=NULL, TH3D *diffpixthoff=NULL,TH3D *sigmapixthoff=NULL, - TH2D *blindidthoff=NULL, TH2D *blindnthoff=NULL); + TH2D& sigthmc, TH3D& diffpixthmc, TH3D& sigmapixthmc, + TH2D& blindidthmc, TH2D& blindnthmc, + TH2D& sigthon, TH3D& diffpixthon, TH3D& sigmapixthon, + TH2D& blindidthon, TH2D& blindnthon, + TH2D& sigthoff, TH3D& diffpixthoff,TH3D& sigmapixthoff, + TH2D& blindidthoff, TH2D& blindnthoff); Bool_t MergeONMC( - TH2D *sigthmc, TH3D *diffpixthmc, TH3D *sigmapixthmc, - TH2D *blindidthmc, TH2D *blindnthmc, - TH2D *sigthon, TH3D *diffpixthon, TH3D *sigmapixthon, - TH2D *blindidthon, TH2D *blindnthon); + TH2D& sigthmc, TH3D& diffpixthmc, TH3D& sigmapixthmc, + TH2D& blindidthmc, TH2D& blindnthmc, + TH2D& sigthon, TH3D&diffpixthon, TH3D& sigmapixthon, + TH2D& blindidthon, TH2D& blindnthon); Bool_t ReadPaddingDist(const char *filein); Index: trunk/MagicSoft/Mars/manalysis/MSigmabar.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MSigmabar.cc (revision 4577) +++ trunk/MagicSoft/Mars/manalysis/MSigmabar.cc (revision 4584) @@ -30,7 +30,10 @@ // // // This is the storage container to hold information about // -// "average" pedestal sigmas // -// // -// In calculating averages all sigmas are normalized to the area of pixel 0// +// "average" pedestal RMS // +// // +// The "average" pedestal RMS is calculated as // +// = sqrt( sum_i( (pedRMS_i)^2/area_i ) / no.of pixels ) // +// // +// which is the sqrt of the average (pedRMS^2 per area) // // // ///////////////////////////////////////////////////////////////////////////// @@ -109,5 +112,5 @@ // - // sum up sigma/sqrt(area) for each sector, + // sum up sigma^2/area for each sector, // separately for inner and outer region; // @@ -147,4 +150,6 @@ const Double_t ratio = geom.GetPixRatio(idx); + //*fLog << "pixel id, ratio = " << idx << ", " << ratio << endl; + const MGeomPix &gpix = geom[idx]; @@ -171,10 +176,10 @@ { innerPixels[sector]++; - innerSum[sector]+= sigma * sqrt(ratio); + innerSum[sector]+= sigma*sigma * ratio; } else { outerPixels[sector]++; - outerSum[sector]+= sigma * sqrt(ratio); + outerSum[sector]+= sigma*sigma * ratio; } } @@ -192,12 +197,12 @@ } - if (fInnerPixels > 0) fSigmabarInner = fSumInner / fInnerPixels; - if (fOuterPixels > 0) fSigmabarOuter = fSumOuter / fOuterPixels; - - // - // this is the average sigma/sqrt(area) (over all pixels) + if (fInnerPixels > 0) fSigmabarInner = sqrt(fSumInner / fInnerPixels); + if (fOuterPixels > 0) fSigmabarOuter = sqrt(fSumOuter / fOuterPixels); + + // + // this is the sqrt of the average sigma^2/area // fSigmabar = (fInnerPixels+fOuterPixels)<=0 ? 0: - (fSumInner+fSumOuter)/(fInnerPixels+fOuterPixels); + sqrt( (fSumInner+fSumOuter)/(fInnerPixels+fOuterPixels) ); for (UInt_t i=0; i<6; i++) @@ -209,7 +214,7 @@ const Double_t sum = ip + op; - fSigmabarInnerSector[i] = ip <=0 ? 0 : iss/ip; - fSigmabarOuterSector[i] = op <=0 ? 0 : oss/op; - fSigmabarSector[i] = sum<=0 ? 0 : (iss+oss)/sum; + fSigmabarInnerSector[i] = ip <=0 ? 0 : sqrt(iss/ip); + fSigmabarOuterSector[i] = op <=0 ? 0 : sqrt(oss/op); + fSigmabarSector[i] = sum<=0 ? 0 : sqrt((iss+oss)/sum); } @@ -217,5 +222,5 @@ //Print(opt); - return fSigmabar; + return fSigmabarInner; } Index: trunk/MagicSoft/Mars/manalysis/MSigmabar.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MSigmabar.h (revision 4577) +++ trunk/MagicSoft/Mars/manalysis/MSigmabar.h (revision 4584) @@ -14,5 +14,5 @@ { private: - Float_t fSigmabar; // Sigmabar ("average" RMS) of pedestal + Float_t fSigmabar; // Sigmabar ( sqrt(average pedestalRMS^2) ) Float_t fSigmabarInner; // --only for inner pixels Float_t fSigmabarOuter; // --only for outer pixels Index: trunk/MagicSoft/Mars/manalysis/MSigmabarCalc.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MSigmabarCalc.cc (revision 4577) +++ trunk/MagicSoft/Mars/manalysis/MSigmabarCalc.cc (revision 4584) @@ -37,5 +37,5 @@ // MPedPhotCam // MRawRunHeader -// MMcEvt (FIXME: Must be replaced by a 'real-data' container) +// MPointingPos // MCerPhotEvt // @@ -58,5 +58,5 @@ #include "MSigmabarParam.h" -#include "MMcEvt.hxx" +#include "MPointingPos.h" ClassImp(MSigmabarCalc); @@ -105,8 +105,8 @@ // This is needed for determining min/max Theta - fMcEvt = (MMcEvt*)pList->FindObject(AddSerialNumber("MMcEvt")); - if (!fMcEvt) + fPointPos = (MPointingPos*)pList->FindObject(AddSerialNumber("MPointingPos")); + if (!fPointPos) { - *fLog << warn << "MMcEvt not found... aborting." << endl; + *fLog << warn << "MPointingPos not found... aborting." << endl; fThetaMin = 0; fThetaMax = 90; @@ -148,11 +148,11 @@ if (rcGetTelescopeTheta()*kRad2Deg; + const Double_t theta = fPointPos->GetZd(); - if (theta>fSigmabarMax) fSigmabarMax=theta; - if (thetafThetaMax) fThetaMax=theta; + if (theta