Changeset 6368 for trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.cc

Timestamp:

02/10/05 19:54:14 (20 years ago)

Author:

rico

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.cc (modified) (12 diffs)

trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.cc

@@ -37 +37 @@
 #include "MHillas.h"
 #include "MMcEvt.hxx"
-#include "TH1D.h"
+#include "TH2F.h"
 
 #include "MLog.h"
@@ -46 +46 @@
 using namespace std;
 
+const Int_t fNTbins = 2;
+const Double_t fTbin[fNTbins+1] = {0,10,20};
+
 
 // -------------------------------------------------------------------------
@@ -52 +55 @@
 //
 MEffAreaAndCoeffCalc::MEffAreaAndCoeffCalc()
-  : fSpec(NULL), fHorig(NULL), fEbins(10), fEsubbins(20), fEmin(10.), fEmax(10000.)
+  : fSpec(NULL), fHorig(NULL), fEbins(10), fEsubbins(20), fEmin(10.), fEmax(10000.),
+    fCoeff(NULL), fEffA(NULL)
 {
   // set the default function
@@ -64 +68 @@
 
   // define some usefull aliases
-  fCini->SetAlias("energy","MMcEvt.fEnergy");
   fCini->SetAlias("logenergy","log10(MMcEvt.fEnergy)");
-  fCini->SetAlias("hastriggered","MMcTrig.fNumFirstLevel");
+  fCini->SetAlias("theta","MMcEvt.fTheta*180./3.14159");
 
   // OJO!! the estimated energy must be changed for a real one
-  fCcut->SetAlias("logestenergy","log10(MHillas.fSize/15.)");
   fCcut->SetAlias("logenergy","log10(MMcEvt.fEnergy)");
+  fCcut->SetAlias("theta","MMcEvt.fTheta*180./3.14159");
 
   fCcut->SetBranchStatus("*",0);
@@ -92 +95 @@
 
   if(fCoeff)
-    delete [] fCoeff;
+    delete fCoeff;
 
   if(fEffA)
-    delete [] fEffA;
+    delete fEffA;
 
   if(fHorig) 
@@ -142 +145 @@
   
   if(fHorig) delete fHorig;
-  fHorig   = new TH1D("horig","Original energy spectrum",esbins,logemin,logemax);
-  fCini->Draw("logenergy>>horig","","goff");
+  fHorig   = new TH1F("fhorig","Original energy spectrum",esbins,logemin,logemax);
+  fCini->Draw("logenergy>>fhorig","","goff");
 }
 
@@ -201 +204 @@
  
   // declare needed histos
-  TH1D* hest = new TH1D("hest","Estimated energy",fEbins,logemin,logemax);
-  TH1D* hmc  = new TH1D("hmc","MC energy",fEbins,logemin,logemax);
+  TH2F* hest = new TH2F("hest","Estimated energy",fEbins,logemin,logemax,fNTbins,fTbin);
+  TH2F* hmc  = new TH2F("hmc","MC energy",fEbins,logemin,logemax,fNTbins,fTbin);
 
   // reset the coefficients
   if(fCoeff)
-    delete [] fCoeff;
-  fCoeff = new Double_t[fEbins];
+    delete fCoeff;
+  fCoeff = new TH2F("fcoeff","Coefficients for unfolding",fEbins,logemin,logemax,fNTbins,fTbin);
 
   // fill the histograms of estimated and measured energy for weighted events
@@ -218 +221 @@
       Float_t emc   = fMcEvt->GetEnergy();
       Float_t estim = fHillas->GetSize()/15.;
+      Float_t theta = fMcEvt->GetTheta()*180./3.14159; // zenith angle (in degrees)
 
       // compute the bin where the weight is taken from
@@ -225 +229 @@
       if(wbin<esbins)
         {
-          hest->Fill(TMath::Log10(estim),fWeight[wbin]); 
-          hmc->Fill(TMath::Log10(emc),fWeight[wbin]);
+          hest->Fill(TMath::Log10(estim),theta,fWeight[wbin]); 
+          hmc->Fill(TMath::Log10(emc),theta,fWeight[wbin]);
         }
       else
@@ -233 +237 @@
 
   // divide the previous histos to get the coefficients for unfolding
-  for(Int_t i=0;i<fEbins;i++)
-    {
-      const Float_t num = hmc->GetBinContent(i+1);
-      const Float_t den = hest->GetBinContent(i+1);
-      if(den)
-        fCoeff[i]=num/den;
-      else
-        {
-          cout << "MEffAreaAndCoeffCalc::ComputeCoefficients Warning: energy bin " << i << " has no survivors, setting coefficient to 0" << endl;
-          fCoeff[i]=0;
-        }
-    }
+  for(Int_t j=0;j<fNTbins;j++)
+    for(Int_t i=0;i<fEbins;i++)
+      {
+        const Float_t num = hmc->GetBinContent(i+1,j+1);
+        const Float_t den = hest->GetBinContent(i+1,j+1);
+        if(den)
+          fCoeff->SetBinContent(i+1,j+1,num/den);
+        else
+          {
+            cout << "MEffAreaAndCoeffCalc::ComputeCoefficients Warning: energy bin " << i << ", thetabin " << j << " has no survivors, setting coefficient to 0" << endl;
+            fCoeff->SetBinContent(i+1,j+1,0.);
+          }
+      }
+
 
   delete hmc;
@@ -256 +262 @@
 void MEffAreaAndCoeffCalc::ComputeEffectiveAreas()
 {
+  if(!fWeight)
+    {
+      cout << "MEffAreaAndCoeffCalc::ComputeEffectiveAreas Warning: No weights computed! nothing done" << endl;
+      return;
+    }
+
   //OJO!! radius should be read from somewhere!
   const Double_t radius = 30000.; // generation radius (in cm)
@@ -265 +277 @@
   // reset the effective areas
   if(fEffA)
-    delete [] fEffA;
-  fEffA = new Double_t[esbins];
-
-  // check if the original spectrum is filled
-  if(!fHorig)
-    FillOriginalSpectrum();
+    delete fEffA;
+  fEffA = new TH2F("feffa","Effective area",esbins,logemin,logemax,fNTbins,fTbin);
 
   // fill the spectrum of the survivors
-  TH1F* hpass= new TH1F("hpass","Survivors",esbins,logemin,logemax);
-  fCcut->Draw("logenergy>>hpass","","goff");
-
+  TH2F* hpass= new TH2F("hpass","Survivors",esbins,logemin,logemax,fNTbins,fTbin);
+  TH2F* horig= new TH2F("horig","Original events",esbins,logemin,logemax,fNTbins,fTbin);
+
+  fCcut->Draw("theta:logenergy>>hpass","","goff");
+  fCini->Draw("theta:logenergy>>horig","","goff");
+  
   // compute the effective areas
-  for(Int_t i=0;i<esbins;i++)
-    if(fHorig->GetBinContent(i+1)<=0)
+  for(Int_t j=0;j<fNTbins;j++)
+    for(Int_t i=0;i<esbins;i++)
       {
-        cout << "MEffAreaAndCoeffCalc::ComputeEffectiveAreas Warning: energy bin " << i << " contains no events setting effective area to 0" << endl;
-        fEffA[i]=0;
+        if(horig->GetBinContent(i+1,j+1)<=0)
+          {
+            cout << "MEffAreaAndCoeffCalc::ComputeEffectiveAreas Warning: energy bin " << i <<", theta bin " << j << " contains no events setting effective area to 0" << endl;
+            fEffA->SetBinContent(i+1,j+1,0);
+          }
+        else
+          fEffA->SetBinContent(i+1,j+1,Atot*hpass->GetBinContent(i+1,j+1)/horig->GetBinContent(i+1,j+1));
       }
-    else
-      fEffA[i]=Atot*hpass->GetBinContent(i+1)/fHorig->GetBinContent(i+1);
-
+  
+
+  
   delete hpass;
+  delete horig;
 }