/* ======================================================================== *\ ! ! * ! * This file is part of MARS, the MAGIC Analysis and Reconstruction ! * Software. It is distributed to you in the hope that it can be a useful ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. ! * It is distributed WITHOUT ANY WARRANTY. ! * ! * Permission to use, copy, modify and distribute this software and its ! * documentation for any purpose is hereby granted without fee, ! * provided that the above copyright notice appear in all copies and ! * that both that copyright notice and this permission notice appear ! * in supporting documentation. It is provided "as is" without express ! * or implied warranty. ! * ! ! ! Author(s): Wolfgang Wittek 5/2002 ! ! Copyright: MAGIC Software Development, 2000-2002 ! ! \* ======================================================================== */ ////////////////////////////////////////////////////////////////////////////// // // // MHFlux // // // // calculates absolute photon fluxes // // from the distributions of the estimated energy // // for the different bins in some variable 'Var' // // (Var = Theta or time) // // // ////////////////////////////////////////////////////////////////////////////// #include "MHFlux.h" #include #include #include #include #include #include "MTime.h" #include "MBinning.h" #include "MParList.h" #include "MLog.h" #include "MLogManip.h" ClassImp(MHFlux); // -------------------------------------------------------------------------- // // Default Constructor. It sets the variable name (Theta or time) // and the units for the variable // MHFlux::MHFlux(const TH2D &h2d, const Bool_t Draw, // const char *varname, const char *unit) const TString varname, const TString unit) : fHOrig(), fHUnfold(), fHFlux(), fVarname(), fUnit() { // if (varname == NULL || unit == NULL) if (varname == "" || unit == "") { *fLog << "MHFlux : varname or unit not defined" << endl; } fVarname = varname; fUnit = unit; TString strg(varname); strg += unit; // char txt[100]; // original distribution of E-est for different bins // of the variable (Theta or time) // sprintf(txt, "gammas vs. E-est and %s",varname); TString strg1 = "no.of gammas vs. E-est and "; strg1 += varname; ((TH2D&)h2d).Copy(fHOrig); fHOrig.SetName("E-est"); fHOrig.SetTitle(strg1); fHOrig.SetDirectory(NULL); fHOrig.SetXTitle("E-est [GeV] "); fHOrig.SetYTitle(strg); fHOrig.Sumw2(); SetBinning((TH2*)&fHOrig, (TH2*)&h2d); // unfolded distribution of E-unfold for different bins // of the variable (Theta or time) // sprintf(txt, "gammas vs. E-unfold and %s",varname); TString strg2 = "no.of gammas vs. E-unfold and "; strg2 += varname; fHUnfold.SetName("E-unfolded"); fHUnfold.SetTitle(strg2); fHUnfold.SetDirectory(NULL); fHUnfold.SetXTitle("E-unfold [GeV] "); fHUnfold.SetYTitle(strg); fHUnfold.Sumw2(); SetBinning((TH2*)&fHUnfold, (TH2*)&fHOrig); // absolute photon flux vs. E-unfold // for different bins of the variable (Theta or time) // // sprintf(txt, "gamma flux [1/(s m2 GeV) vs. E-unfold and %s",varname); TString strg3 = "gamma flux [1/(s m2 GeV) vs. E-unfold and "; strg3 += varname; fHFlux.SetName("photon flux"); fHFlux.SetTitle(strg3); fHFlux.SetDirectory(NULL); fHFlux.SetXTitle("E-unfold [GeV] "); fHFlux.SetYTitle(strg); fHFlux.Sumw2(); SetBinning((TH2*)&fHFlux, (TH2*)&fHUnfold); // copy fHOrig into fHUnfold in case no unfolding is done const Int_t nEunf = fHUnfold.GetNbinsX(); const Int_t nVar = fHUnfold.GetNbinsY(); for (int m=1; m<=nEunf; m++) { for (int n=1; n<=nVar; n++) { Double_t cont = fHOrig.GetBinContent(m,n); Double_t dcont = fHOrig.GetBinError(m,n); fHUnfold.SetBinContent(m,n,cont); fHUnfold.SetBinError(m,n,dcont); } } //.............................................. // draw the No.of photons vs. E-est // for the individual bins of the variable Var if (Draw == kTRUE) { const Int_t nVar = fHOrig.GetNbinsY(); for (int n=1; n<=nVar; n++) { TString strg0("Orig-"); strg0 += fVarname; TH1D &h = *fHOrig.ProjectionX(strg0, n, n, "E"); char txt0[100]; strg0 = fVarname; strg0 += "-bin %d"; sprintf(txt0, strg0, n); TString strg1("No.of photons vs. E-est for "); strg1 += txt0; new TCanvas(txt0,strg1); // TCanvas &c = *MakeDefCanvas(txt0, strg1); // gROOT->SetSelectedPad(NULL); gPad->SetLogx(); h.SetName(txt0); h.SetTitle(strg1); h.SetXTitle("E-est [GeV] "); h.SetYTitle("No.of photons"); h.DrawCopy(); // c.Modified(); // c.Update(); } } //........................ } // ------------------------------------------------------------------------- // // Dummy Fill (has to be included because in base class MH Fill is set to 0); // without the dummy Fill one gets the error message : // // Error: Can't call MHFlux::MHFlux(evttime,"time","[s]") in current scope // FILE:macros/flux.C LINE:465 // Possible candidates are... // filename line:size busy function type and name (in MHFlux) // filename line:size busy function type and name (in MH) // filename line:size busy function type and name (in MParContainer) // filename line:size busy function type and name (in TObject) // Bool_t MHFlux::Fill(const MParContainer *par) { return kTRUE; } // ------------------------------------------------------------------------- // // Unfold the distribution in E-est // void MHFlux::Unfold(const Bool_t Draw) { //.............................................. // draw the No.of photons vs. E-unfold // for the individual bins of the variable Var if (Draw == kTRUE) { const Int_t nVar = fHUnfold.GetNbinsY(); for (int n=1; n<=nVar; n++) { TString strg0("Unfold-"); strg0 += fVarname; TH1D &h = *fHUnfold.ProjectionX(strg0, n, n, "E"); char txt0[100]; strg0 = fVarname; strg0 += "-bin %d"; sprintf(txt0, strg0, n); TString strg1("No.of photons vs. E-unfold for "); strg1 += txt0; new TCanvas(txt0,strg1); // TCanvas &c = *MakeDefCanvas(txt0, strg1); // gROOT->SetSelectedPad(NULL); gPad->SetLogx(); h.SetName(txt0); h.SetTitle(strg1); h.SetXTitle("E-unfold [GeV] "); h.SetYTitle("No.of photons"); h.DrawCopy(); // c.Modified(); // c.Update(); } } //........................ } // ------------------------------------------------------------------------- // // Calculate photon flux by dividing the distribution in Eunf (fHUnfold) by // the width of the energy interval (deltaE) // the effective ontime (*teff) // and the effective collection area (*aeff) // void MHFlux::CalcFlux(const TH1D *teff, const TProfile *thetabar, const TH2D *aeff, const Bool_t Draw) { // Note that fHUnfold has bins in Eunf and Var // *teff has bins in Var (the same bins in Var as fHUnfold) // *thetabar has bins in Var (the same bins in Var as fHUnfold) // *aeff has bins in Etru and Theta // (where in general the binning in Etru is different // from the binning in Eunf) // The variable Var may be 'time' or 'Theta' // Draw = kTRUE means the differential flux vs E-unf should be drawn // for the individual bins of the variable Var //.................................... // define dummy histogram *aeff ((TH1*)aeff)->Sumw2(); MBinning binsetru("BinningEtru"); binsetru.SetEdgesLog(10, 10, 1e3); MBinning binsthetatru("BinningThetatru"); binsthetatru.SetEdges(7, -2.5, 32.5); //SetBinning((TH1*)aeff, &binsetru, &binsthetatru); SetBinning((TH2*)aeff, &binsetru, &binsthetatru); const Int_t netru = aeff->GetNbinsX(); const Int_t ntheta = aeff->GetNbinsY(); for (int j=1; j<=netru; j++) { for (int k=1; k<=ntheta; k++) { Double_t cont = 10000.0;; ((TH1*)aeff)->SetBinContent(j,k,cont); Double_t dcont = 100.0; ((TH1*)aeff)->SetBinError(j,k,dcont); } } // *fLog << "Dummy aeff : netru =" << netru << ", ntheta = " << ntheta << endl; //.................................... // number of Eunf and Var bins (histograms : fHUnfold, fHFlux) const Int_t nEunf = fHFlux.GetNbinsX(); const Int_t nVar = fHFlux.GetNbinsY(); // number of Etru and Theta bins (histogram *aeff of collection area) const Int_t nEtru = aeff->GetNbinsX(); const Int_t nTheta = aeff->GetNbinsY(); //*fLog << "nEunf =" << nEunf << ", nVar = " << nVar << endl; //*fLog << "nEtru =" << nEtru << ", nTheta = " << nTheta << endl; //................................................................... // calculate effective collection area // for the Eunf and Var bins of the histogram fHUnfold // from the histogram *aeff, which has bins in Etru and Theta // the result is the histogram fHAeff // TH2D fHAeff; fHAeff.Sumw2(); SetBinning((TH2*)&fHAeff, (TH2*)&fHUnfold); Int_t errflag; Double_t c0, c1, c2; Double_t t1, t2, t3; Double_t a1, a2, a3; Double_t *aeffbar; aeffbar = new Double_t[nEtru]; Double_t *daeffbar; daeffbar = new Double_t[nEtru]; Double_t aeffEunfVar; Double_t daeffEunfVar; //------ start n loop ------ for (int n=1; n<=nVar; n++) { Double_t Thetabar = thetabar->GetBinContent(n); Double_t cosThetabar = cos(Thetabar); // determine Theta bins (k1, k2, k3) for interpolation in Theta // k0 denotes the Theta bin from whicvh the error is copied Int_t k0=0, k1=0, k2=0, k3=0; for (int k=3; k<=nTheta; k++) { Double_t Thetalow = ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(k); if (Thetabar < Thetalow) { k1 = k-2; k2 = k-1; k3 = k; k0 = k2; break; } } if (k3 == 0) { k1 = nTheta-2; k2 = nTheta-1; k3 = nTheta; k0 = k2; } if (Thetabar < ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(2)) k0 = 1; else if (Thetabar > ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(nTheta)) k0 = nTheta; Double_t Thetamin = ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(1); Double_t Thetamax = ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(nTheta+1); if (Thetabar < Thetamin || Thetabar > Thetamax) { *fLog << "MHFlux.cc : extrapolation in Theta; Thetabar = " << Thetabar << ", Thetamin =" << Thetamin << ", Thetamax =" << Thetamax << endl; } //*fLog << "Var bin " << n << ":" << endl; //*fLog << "Thetabar= " << Thetabar // << ", k0= " << k0 // << ", k1= " << k1 // << ", k2= " << k2 // << ", k3= " << k3 << endl; // calculate effective collection area at Theta = Thetabar // by quadratic interpolation in cos(Theta); // do this for each bin of Etru // for (int j=1; j<=nEtru; j++) { c0 = 0.0; c1 = 0.0; c2 = 0.0; t1 = cos( ((TH1*)aeff)->GetYaxis()->GetBinCenter (k1) ); t2 = cos( ((TH1*)aeff)->GetYaxis()->GetBinCenter (k2) ); t3 = cos( ((TH1*)aeff)->GetYaxis()->GetBinCenter (k3) ); a1 = aeff->GetBinContent(j,k1); a2 = aeff->GetBinContent(j,k2); a3 = aeff->GetBinContent(j,k3); Parab(t1, t2, t3, a1, a2, a3, &c0, &c1, &c2, &errflag); aeffbar[j] = c0 + c1*cosThetabar + c2*cosThetabar*cosThetabar; daeffbar[j] = aeff->GetBinError(j,k0); //*fLog << "Etru bin " << j << ": tbar= " << Thetabar // << ", abar= " << aeffbar[j] // << ", dabar= " << daeffbar[j] << endl; } //--- start m loop --- // calculate effective collection area at (E = Ebar, Theta = Thetabar) // by quadratic interpolation in log10(Etru) // do this for each bin of Eunf // for (int m=1; m<=nEunf; m++) { Double_t log10Ebar = 0.5 * ( log10( fHUnfold.GetXaxis()->GetBinLowEdge(m) ) +log10( fHUnfold.GetXaxis()->GetBinLowEdge(m+1)) ); Double_t Ebar = pow(10.0, log10Ebar); // determine Etru bins (j1, j2, j3) for interpolation in E // j0 denotes the Etru bin from which the error is copied Int_t j0=0, j1=0, j2=0, j3=0; for (int j=3; j<=nEtru; j++) { Double_t Elow = ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(j); if (Ebar < Elow) { j1 = j-2; j2 = j-1; j3 = j; j0 = j2; break; } } if (j3 == 0) { j1 = nEtru-2; j2 = nEtru-1; j3 = nEtru; j0 = j2; } if (Ebar < ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(2)) j0 = 1; else if (Ebar > ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(nEtru)) j0 = nEtru; Double_t Etrumin = ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(1); Double_t Etrumax = ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(nEtru+1); if (Ebar < Etrumin || Ebar > Etrumax) { *fLog << "MHFlux.cc : extrapolation in Energy; Ebar = " << Ebar << ", Etrumin =" << Etrumin << ", Etrumax =" << Etrumax << endl; } //*fLog << "Var bin " << n << ":" << endl; //*fLog << "Ebar= " << Ebar // << ", j1= " << j1 // << ", j2= " << j2 // << ", j3= " << j3 << endl; c0=0.0; c1=0.0; c2=0.0; t1 = 0.5 * ( log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j1) ) +log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j1+1)) ); t2 = 0.5 * ( log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j2) ) +log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j2+1)) ); t3 = 0.5 * ( log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j3) ) +log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j3+1)) ); a1 = aeffbar[j1]; a2 = aeffbar[j2]; a3 = aeffbar[j3]; Parab(t1, t2, t3, a1, a2, a3, &c0, &c1, &c2, &errflag); aeffEunfVar = c0 + c1*log10(Ebar) + c2*log10(Ebar)*log10(Ebar); daeffEunfVar = daeffbar[j0]; //*fLog << "Eunf bin " << m << ": Ebar= " << Ebar // << ", aeffEunfVar = " << aeffEunfVar // << ", daeffEunfVar = " << daeffEunfVar << endl; fHAeff.SetBinContent(m,n,aeffEunfVar); fHAeff.SetBinError(m,n,daeffEunfVar); } //--- end m loop --- } //------ end n loop ------ delete aeffbar; //................................................................... // now calculate the absolute gamma flux // for (int m=1; m<=nEunf; m++) { Double_t DeltaE = fHFlux.GetXaxis()->GetBinWidth(m); for (int n=1; n<=nVar; n++) { Double_t Ngam = fHUnfold.GetBinContent(m,n); Double_t dNgam = fHUnfold.GetBinError(m,n); Double_t Aeff = fHAeff.GetBinContent(m,n); Double_t dAeff = fHAeff.GetBinError(m,n); Double_t Effon = teff->GetBinContent(n); Double_t dEffon = teff->GetBinError(n); Double_t Cont, dCont; if (Ngam > 0.0 && DeltaE > 0.0 && Effon > 0.0 && Aeff > 0.0) { Cont = Ngam / (DeltaE * Effon * Aeff); dCont = Cont * sqrt( dNgam*dNgam / (Ngam*Ngam) + dEffon*dEffon / (Effon*Effon) + dAeff*dAeff / (Aeff*Aeff) ); } else { Cont = 1.e-20; dCont = 1.e-20; } fHFlux.SetBinContent(m,n,Cont); fHFlux.SetBinError(m,n,dCont); //*fLog << "Eunf bin " << m << ", Var bin " << n // << ": Ngam = " << Ngam << ", Flux = " // << Cont << ", dFlux = " << dCont << endl; //*fLog << ", DeltaE = " << DeltaE << ", Effon = " << Effon // << ", Aeff = " << Aeff << endl; } } //.............................................. // draw the differential photon flux vs. E-unf // for the individual bins of the variable Var if (Draw == kTRUE) { for (int n=1; n<=nVar; n++) { TString strg0("Flux-"); strg0 += fVarname; TH1D &h = *fHFlux.ProjectionX(strg0, n, n, "E"); char txt[100]; TString strg1("Photon flux vs. E-unfold for "); TString strg2 = fVarname; strg2 += "-bin %d"; sprintf(txt, strg2, n); TString strg3 = strg1 + txt; new TCanvas(txt, strg3); // TCanvas &c = *MakeDefCanvas(txt, txt); // gROOT->SetSelectedPad(NULL); gPad->SetLogx(); h.SetName(txt); h.SetTitle(strg3); h.SetXTitle("E-unfold [GeV] "); h.SetYTitle("photons / (s m2 GeV)"); h.DrawCopy(); // c.Modified(); // c.Update(); } } //........................ } // ------------------------------------------------------------------------- // // Draw copies of the histograms // TObject *MHFlux::DrawClone(Option_t *opt) const { TCanvas &c = *MakeDefCanvas("flux", "Orig - Unfold - Flux plots"); c.Divide(2, 2); gROOT->SetSelectedPad(NULL); c.cd(1); ((TH2*)&fHOrig)->DrawCopy(""); gPad->SetLogx(); c.cd(2); ((TH2*)&fHUnfold)->DrawCopy(""); gPad->SetLogx(); c.cd(3); ((TH2*)&fHFlux)->DrawCopy(""); gPad->SetLogx(); c.Modified(); c.Update(); return &c; } // ------------------------------------------------------------------------- // // Draw the histograms // void MHFlux::Draw(Option_t *opt) { if (!gPad) MakeDefCanvas("flux", "orig-unfold-flux plots"); gPad->Divide(2,2); gPad->cd(1); fHOrig.Draw(opt); gPad->cd(2); fHUnfold.Draw(opt); gPad->cd(3); fHFlux.Draw(opt); gPad->Modified(); gPad->Update(); } // ------------------------------------------------------------------------- // // Quadratic interpolation // // *** calculate the parameters of a parabula // y = a + b*x + c*x**2 = F(x) // such that yi = F(xi) for (i=1,3) // void MHFlux::Parab(Double_t x1, Double_t x2, Double_t x3, Double_t y1, Double_t y2, Double_t y3, Double_t *a, Double_t *b, Double_t *c, Int_t *errflag) { double ai11,ai12,ai13,ai21,ai22,ai23,ai31,ai32,ai33; double det,det1; //double yt1,yt2,yt3; det = x2*x3*x3 + x1*x2*x2 + x3*x1*x1 - x2*x1*x1 - x3*x2*x2 - x1*x3*x3; if (det != 0.0) { *errflag = 0; det1 = 1.0/det; ai11 = x2*x3*x3 - x3*x2*x2; ai12 = x3*x1*x1 - x1*x3*x3; ai13 = x1*x2*x2 - x2*x1*x1; ai21 = x2*x2 - x3*x3; ai22 = x3*x3 - x1*x1; ai23 = x1*x1 - x2*x2; ai31 = x3 - x2; ai32 = x1 - x3; ai33 = x2 - x1; *a = (ai11*y1 + ai12*y2 + ai13*y3) * det1; *b = (ai21*y1 + ai22*y2 + ai23*y3) * det1; *c = (ai31*y1 + ai32*y2 + ai33*y3) * det1; //yt1 = *a + *b * x1 + *c * x1*x1; //yt2 = *a + *b * x2 + *c * x2*x2; //yt3 = *a + *b * x3 + *c * x3*x3; //*fLog << "x1 = " << x1 << ", x2 = " << x2 << ", x3 = " << x3 << endl; //*fLog << "y1 = " << y1 << ", y2 = " << y2 << ", y3 = " << y3 << endl; //*fLog << "yt1 = " << yt1 << ", yt2 = " << yt2 // << ", yt3 = " << yt3 << endl; //*fLog << "*a = " << *a << ", *b = " << *b << ", *c= " << *c // << ", *errflag = " << *errflag << endl; return; } *errflag = 1; *a = 0.0; *b = 0.0; *c = 0.0; return; }