Changeset 6368

Timestamp:

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

Author:

rico

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mtemp/mifae

Files:

• library/MEffAreaAndCoeffCalc.cc (modified) (12 diffs)
• library/MEffAreaAndCoeffCalc.h (modified) (4 diffs)
• macros/computeCoeff.C (modified) (3 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;
 }
 

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

@@ -7 +7 @@
 
 class TF1;
-class TH1D;
+class TH1F;
+class TH2F;
 class MHillas;
 class MMcEvt;
@@ -18 +19 @@
 
   TF1* fSpec;        // function used to parametrized the spectrum
-  TH1D* fHorig;      // histogram with the original energy spectrum
+  TH1F* fHorig;      // histogram with the original energy spectrum
 
   Int_t fEbins;      // number of bins to build spectrum
@@ -25 +26 @@
   Double_t fEmax;    // Maximum energy in GeV
 
-  
   Double_t* fWeight; // array containing weights
-  Double_t* fCoeff;  // array containing coefficients
-  Double_t* fEffA;   // array containing effective areas
+  TH2F* fCoeff;      // histogram containing unfolding coefficients
+  TH2F* fEffA;       // histogram containing effective areas
 
-  TChain* fCini; // chain for initial MC events (even those not triggering)
-  TChain* fCcut; // chain for surviving MC events (after cuts)
+  TChain* fCini;     // chain for initial MC events (even those not triggering)
+  TChain* fCcut;     // chain for surviving MC events (after cuts)
 
-  MHillas* fHillas; // pointer to the MHillas Branch
+  MHillas* fHillas;  // pointer to the MHillas Branch
   MMcEvt*  fMcEvt;   // pointer to the MMcEvt Branch
 
  protected:
+
   void FillOriginalSpectrum();
   void ComputeCoefficients();
@@ -51 +52 @@
   void SetEbins(Int_t i)    {fEbins=i;}
   void SetEsubbins(Int_t i) {fEsubbins=i;}
-  void SetE(Float_t x)        {fEmin=x;}
-  void SetEbins(Float_t x)    {fEmax=x;}
+  void SetEmin(Float_t x)   {fEmin=x;}
+  void SetEmax(Float_t x)   {fEmax=x;}
 
   void AddFileToInitialMC(const Char_t* name) {fCini->Add(name);}
   void AddFileToFinalMC(const Char_t* name) {fCcut->Add(name);}
 
-  Double_t GetEffectiveArea(Int_t i) {return fEffA[i];}
-  Double_t GetCoefficient(Int_t i)   {return fCoeff[i];}
+  TH2F* GetEffectiveAreaHisto() {return fEffA;}
+  TH2F* GetCoefficientHisto()   {return fCoeff;}
 
   void ComputeAllFactors();

trunk/MagicSoft/Mars/mtemp/mifae/macros/computeCoeff.C

@@ -34 +34 @@
   const Int_t esubbins = 20;           // number of subbins per big bin
   const Double_t emin=10;
-  const Double_t emax=10000;
+  const Double_t emax=1000;
   const Double_t logemin = TMath::Log10(emin);
   const Double_t logemax = TMath::Log10(emax);
@@ -79 +79 @@
       UInt_t coefbin = UInt_t((TMath::Log10(estim)-logemin)/de);
       
-      Float_t effa  =  calc.GetEffectiveArea(effabin);
-      Float_t unfold = calc.GetCoefficient(coefbin);
+      Float_t effa  =  calc.GetEffectiveAreaHisto()->GetBinContent(effabin+1,1);
+      Float_t unfold = calc.GetCoefficientHisto()->GetBinContent(coefbin+1,1);
       
       if(effa)
@@ -89 +89 @@
   TFile file("prueba.root","RECREATE");
   hspec->Write();
-
+  calc.GetEffectiveAreaHisto()->Write();
+  calc.GetCoefficientHisto()->Write();
   return;
 }