Index: trunk/MagicSoft/Mars/mtemp/mifae/Changelog
===================================================================
--- trunk/MagicSoft/Mars/mtemp/mifae/Changelog	(revision 6209)
+++ trunk/MagicSoft/Mars/mtemp/mifae/Changelog	(revision 6359)
@@ -19,4 +19,12 @@
                                                  -*-*- END OF LINE -*-*-
 
+ 2005/02/10 Javier Rico
+    * library/MEffAreaAndCoeffCalc.[h,cc]
+      - added class to compute Effective areas and unfolding coefficients
+    * library/Makefile, IFAELinkDef.h
+      - add MEffAreaAndCoeffCalc
+    * macros/computeCoeff.C
+      - add macro example to use MEffAreaAndCoeffCalc	
+	
  2005/02/02 Oscar Blanch Bigas
     * library/MGainFluctuationPix.[h,cc],MGainFluctuationCam.[h,cc]
Index: trunk/MagicSoft/Mars/mtemp/mifae/library/IFAELinkDef.h
===================================================================
--- trunk/MagicSoft/Mars/mtemp/mifae/library/IFAELinkDef.h	(revision 6209)
+++ trunk/MagicSoft/Mars/mtemp/mifae/library/IFAELinkDef.h	(revision 6359)
@@ -32,4 +32,5 @@
 #pragma link C++ class MFDisp+;
 #pragma link C++ class MFindDisp+;
+#pragma link C++ class MEffAreaAndCoeffCalc+;
 #pragma link C++ class MGainFluctuationPix+;
 #pragma link C++ class MGainFluctuationCam+;
Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.cc
===================================================================
--- trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.cc	(revision 6359)
+++ trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.cc	(revision 6359)
@@ -0,0 +1,300 @@
+/* ======================================================================== *\
+!
+! *
+! * 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): Javier Rico     02/2005 <mailto:jrico@ifae.es>
+!
+!   Copyright: MAGIC Software Development, 2000-2005
+!
+!
+\* ======================================================================== */
+
+//////////////////////////////////////////////////////////////////////////////
+//
+// Computes the Effective areas and coefficients for unfolding for a given
+// spectrum that can be parametrized by a function
+//
+//////////////////////////////////////////////////////////////////////////////
+
+#include <fstream>
+#include <math.h>
+
+#include "MEffAreaAndCoeffCalc.h"
+
+#include "TF1.h"
+#include "MHillas.h"
+#include "MMcEvt.hxx"
+#include "TH1D.h"
+
+#include "MLog.h"
+#include "MLogManip.h"
+
+ClassImp(MEffAreaAndCoeffCalc);
+
+using namespace std;
+
+
+// -------------------------------------------------------------------------
+//
+// Constructor
+//
+MEffAreaAndCoeffCalc::MEffAreaAndCoeffCalc()
+  : fSpec(NULL), fHorig(NULL), fEbins(10), fEsubbins(20), fEmin(10.), fEmax(10000.)
+{
+  // set the default function
+  // OJO!! this should be changed
+  SetFunction("4.e9*pow(x,-2.6+1)");  
+
+  // create the TChains 
+  // OJO!! make sure they read the appropriate branch
+  fCini = new TChain("Events");
+  fCcut = new TChain("Events");
+
+  // define some usefull aliases
+  fCini->SetAlias("energy","MMcEvt.fEnergy");
+  fCini->SetAlias("logenergy","log10(MMcEvt.fEnergy)");
+  fCini->SetAlias("hastriggered","MMcTrig.fNumFirstLevel");
+
+  // 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->SetBranchStatus("*",0);
+  fCcut->SetBranchStatus("MHillas.*",1);
+  fCcut->SetBranchStatus("MMcEvt.*",1);
+  fCcut->SetBranchAddress("MHillas.",&fHillas);
+  fCcut->SetBranchAddress("MMcEvt.",&fMcEvt);
+}
+
+// -------------------------------------------------------------------------
+//
+// Destructor
+//
+MEffAreaAndCoeffCalc::~MEffAreaAndCoeffCalc()
+{
+  if(fSpec)
+    delete fSpec;
+
+  if(fWeight)
+    delete [] fWeight;
+
+  if(fCoeff)
+    delete [] fCoeff;
+
+  if(fEffA)
+    delete [] fEffA;
+
+  if(fHorig) 
+    delete fHorig;
+
+  delete fCini;
+  delete fCcut;
+}
+
+// -------------------------------------------------------------------------
+//
+// Set the function by expression and minimum and maximum energies
+//
+void MEffAreaAndCoeffCalc::SetFunction(const Char_t* chfunc, Float_t emin, Float_t emax)
+{
+  if(fSpec)
+    delete fSpec;
+  if(emin<=0 || emax<=0 || emax<emin)
+    {
+      emin = fEmin;
+      emax = fEmax;
+    }
+  fSpec = new TF1("fspec",chfunc,emin,emax);
+}
+
+// -------------------------------------------------------------------------
+//
+// Set the function by function pointer
+//
+void MEffAreaAndCoeffCalc::SetFunction(TF1* newf)
+{
+  if(fSpec)
+    delete fSpec;
+  fSpec = newf;
+}
+
+// -------------------------------------------------------------------------
+//
+// fill the histogram containing the original sample energy spectrum
+//
+void MEffAreaAndCoeffCalc::FillOriginalSpectrum()
+{  
+  const Double_t logemin = TMath::Log10(fEmin);
+  const Double_t logemax = TMath::Log10(fEmax);
+  const Int_t esbins = fEbins*fEsubbins; // total number of subbins
+  
+  if(fHorig) delete fHorig;
+  fHorig   = new TH1D("horig","Original energy spectrum",esbins,logemin,logemax);
+  fCini->Draw("logenergy>>horig","","goff");
+}
+
+// -------------------------------------------------------------------------
+//
+// compute the weights for a particular input spectrum
+//
+void MEffAreaAndCoeffCalc::ComputeWeights()
+{  
+  const Double_t logemin = TMath::Log10(fEmin);
+  const Double_t logemax = TMath::Log10(fEmax);
+  const Double_t de = (logemax-logemin)/fEbins; // bin size (in log)
+  const Double_t desub = de/fEsubbins; // subbin size (in log)
+  const Int_t esbins = fEbins*fEsubbins; // total number of subbins
+
+  // reset the weights array
+  if(fWeight)
+    delete [] fWeight;
+  fWeight = new Double_t[esbins];
+    
+  if(!fHorig)
+    FillOriginalSpectrum();
+  
+  for(Int_t i=0;i<esbins;i++)
+    {
+      const Double_t denom = fHorig->GetBinContent(i+1);              // number of events in initial spectrum
+      const Double_t ew    = TMath::Power(10,logemin+(i+0.5)*desub); // real energy
+      const Double_t numer = fSpec->Eval(ew);                        // number of events for the required spectrum
+      if(denom)
+	fWeight[i]=numer/denom;
+      else
+	{
+	  cout << "MEffAreaAndCoeffCalc::ComputeWeights Warning: no statistic to compute weight for energy " << ew << ", setting it to -1 "  << endl;
+	  fWeight[i]=-1;
+	}
+    }
+}
+
+// --------------------------------------------------------------
+//
+// compute the coefficients used for the (iterative) unfolding 
+//
+void MEffAreaAndCoeffCalc::ComputeCoefficients()
+{ 
+  if(!fWeight)
+    {
+      cout << "MEffAreaAndCoeffCalc::ComputeCoefficients Warning: No weights computed! nothing done" << endl;
+      return;
+    }
+
+  const Double_t logemin = TMath::Log10(fEmin);
+  const Double_t logemax = TMath::Log10(fEmax);
+  const Double_t de = (logemax-logemin)/fEbins; // bin size (in log)
+  const Double_t desub = de/fEsubbins; // subbin size (in log)
+  const Int_t esbins = fEbins*fEsubbins; // total number of subbins
+  const Int_t nentries = Int_t(fCcut->GetEntries());
+ 
+  // declare needed histos
+  TH1D* hest = new TH1D("hest","Estimated energy",fEbins,logemin,logemax);
+  TH1D* hmc  = new TH1D("hmc","MC energy",fEbins,logemin,logemax);
+
+  // reset the coefficients
+  if(fCoeff)
+    delete [] fCoeff;
+  fCoeff = new Double_t[fEbins];
+
+  // fill the histograms of estimated and measured energy for weighted events
+  for(Int_t i=0;i<nentries;i++)
+    {
+      fCcut->GetEntry(i);  
+
+      // mc and estimated energy
+      // OJO!! Estimated energy will be taken directly from some input container
+      Float_t emc   = fMcEvt->GetEnergy();
+      Float_t estim = fHillas->GetSize()/15.;
+
+      // compute the bin where the weight is taken from
+      Int_t wbin = Int_t((TMath::Log10(emc)-logemin)/desub);
+
+      // fill the histograms
+      if(wbin<esbins)
+	{
+	  hest->Fill(TMath::Log10(estim),fWeight[wbin]); 
+	  hmc->Fill(TMath::Log10(emc),fWeight[wbin]);
+	}
+      else
+	cout << "MEffAreaAndCoeffCalc::ComputeCoefficients Warning: event " << i << " with energy " << emc << " has energy out of bounds, skpping" << endl;
+    }
+
+  // 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;
+	}
+    }
+
+  delete hmc;
+  delete hest;
+}
+
+// --------------------------------------------------------------
+//
+// compute the coefficients used for the (iterative) unfolding 
+//
+void MEffAreaAndCoeffCalc::ComputeEffectiveAreas()
+{
+  //OJO!! radius should be read from somewhere!
+  const Double_t radius = 30000.; // generation radius (in cm)
+  const Double_t Atot = 3.14159*radius*radius; //total area (in cm^2)
+  const Double_t logemin = TMath::Log10(fEmin);
+  const Double_t logemax = TMath::Log10(fEmax);
+  const Int_t esbins = fEbins*fEsubbins; // total number of subbins
+
+  // reset the effective areas
+  if(fEffA)
+    delete [] fEffA;
+  fEffA = new Double_t[esbins];
+
+  // check if the original spectrum is filled
+  if(!fHorig)
+    FillOriginalSpectrum();
+
+  // fill the spectrum of the survivors
+  TH1F* hpass= new TH1F("hpass","Survivors",esbins,logemin,logemax);
+  fCcut->Draw("logenergy>>hpass","","goff");
+
+  // compute the effective areas
+  for(Int_t i=0;i<esbins;i++)
+    if(fHorig->GetBinContent(i+1)<=0)
+      {
+	cout << "MEffAreaAndCoeffCalc::ComputeEffectiveAreas Warning: energy bin " << i << " contains no events setting effective area to 0" << endl;
+	fEffA[i]=0;
+      }
+    else
+      fEffA[i]=Atot*hpass->GetBinContent(i+1)/fHorig->GetBinContent(i+1);
+
+  delete hpass;
+}
+
+// --------------------------------------------------------------
+//
+// Call the internal functions to compute all the factors
+//
+void MEffAreaAndCoeffCalc::ComputeAllFactors()
+{
+  ComputeWeights();
+  ComputeCoefficients();
+  ComputeEffectiveAreas();
+}
Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.h
===================================================================
--- trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.h	(revision 6359)
+++ trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.h	(revision 6359)
@@ -0,0 +1,68 @@
+#ifndef MARS_MEffAreaAndCoeffCalc
+#define MARS_MEffAreaAndCoeffCalc
+
+#ifndef MARS_MTask
+#include "MTask.h"
+#endif
+
+class TF1;
+class TH1D;
+class MHillas;
+class MMcEvt;
+
+#include "TChain.h"
+
+class MEffAreaAndCoeffCalc
+{
+ private:
+
+  TF1* fSpec;        // function used to parametrized the spectrum
+  TH1D* fHorig;      // histogram with the original energy spectrum
+
+  Int_t fEbins;      // number of bins to build spectrum
+  Int_t fEsubbins;   // number of subbins per big bin
+  Double_t fEmin;    // Minimum energy in GeV
+  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
+
+  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
+  MMcEvt*  fMcEvt;   // pointer to the MMcEvt Branch
+
+ protected:
+  void FillOriginalSpectrum();
+  void ComputeCoefficients();
+  void ComputeWeights();
+  void ComputeEffectiveAreas();
+
+ public:
+  MEffAreaAndCoeffCalc();
+
+  virtual ~MEffAreaAndCoeffCalc();
+
+  void SetFunction(const Char_t* chfunc, Float_t emin=0., Float_t emax=0.);
+  void SetFunction(TF1*);
+  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 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];}
+
+  void ComputeAllFactors();
+
+  ClassDef(MEffAreaAndCoeffCalc, 0) // task to compute the Effective areas and Coefficients for the unfolding 
+};
+
+#endif
+
Index: trunk/MagicSoft/Mars/mtemp/mifae/library/Makefile
===================================================================
--- trunk/MagicSoft/Mars/mtemp/mifae/library/Makefile	(revision 6209)
+++ trunk/MagicSoft/Mars/mtemp/mifae/library/Makefile	(revision 6359)
@@ -79,4 +79,5 @@
 	MFDisp.cc \
 	MFindDisp.cc \
+	MEffAreaAndCoeffCalc.cc
 	MGainFluctuationPix.cc \
 	MGainFluctuationCam.cc \
Index: trunk/MagicSoft/Mars/mtemp/mifae/macros/computeCoeff.C
===================================================================
--- trunk/MagicSoft/Mars/mtemp/mifae/macros/computeCoeff.C	(revision 6359)
+++ trunk/MagicSoft/Mars/mtemp/mifae/macros/computeCoeff.C	(revision 6359)
@@ -0,0 +1,93 @@
+
+
+void computeCoeff()
+{    
+  MEffAreaAndCoeffCalc calc;
+
+  // initial MC files
+  calc.AddFileToInitialMC("/data/camera_NCL/Gamma_zbin0_0_7_1000to1009_w0.root");
+  calc.AddFileToInitialMC("/data/camera_NCL/Gamma_zbin0_0_7_1130to1139_w0.root");
+  calc.AddFileToInitialMC("/data/camera_NCL/Gamma_zbin0_0_7_1520to1529_w0.root");
+  calc.AddFileToInitialMC("/data/camera_NCL/Gamma_zbin0_90_7_1260to1269_w0.root");
+  calc.AddFileToInitialMC("/data/camera_NCL/Gamma_zbin0_90_7_1390to1399_w0.root");
+  calc.AddFileToInitialMC("/data/camera_NCL/Gamma_zbin0_90_7_1650to1659_w0.root");
+  
+
+  // files with remaining events after cuts (or where to apply cuts) see warning above
+  calc.AddFileToFinalMC("/data/camera_CL/Gamma_zbin0_0_7_1000to1009_w0_cleaned4035.root");
+  calc.AddFileToFinalMC("/data/camera_CL/Gamma_zbin0_0_7_1130to1139_w0_cleaned4035.root");
+  calc.AddFileToFinalMC("/data/camera_CL/Gamma_zbin0_0_7_1520to1529_w0_cleaned4035.root");
+  calc.AddFileToFinalMC("/data/camera_CL/Gamma_zbin0_90_7_1260to1269_w0_cleaned4035.root");
+  calc.AddFileToFinalMC("/data/camera_CL/Gamma_zbin0_90_7_1390to1399_w0_cleaned4035.root");
+  calc.AddFileToFinalMC("/data/camera_CL/Gamma_zbin0_90_7_1650to1659_w0_cleaned4035.root");
+
+
+  // define the funtion of the desired spectrum
+  calc.SetFunction("4.e9*pow(x,-2.6+1)",10.,10000.);
+  calc.ComputeAllFactors();
+  
+  /************************************************/
+  /* Build spectrum                               */
+  /* now it's just a cross check with the same MC */
+  /************************************************/
+  const Int_t ebins = 10;              // number of bins to build spectrum
+  const Int_t esubbins = 20;           // number of subbins per big bin
+  const Double_t emin=10;
+  const Double_t emax=10000;
+  const Double_t logemin = TMath::Log10(emin);
+  const Double_t logemax = TMath::Log10(emax);
+  const Double_t de = (logemax-logemin)/ebins; // bin size (in log)
+  const Double_t desub = de/esubbins; // subbin size (in log)
+  const Int_t esbins = ebins*esubbins; // total number of subbins
+
+  // remaining events after cuts
+  const UInt_t ncutfiles=6;
+  Char_t* cutName[ncutfiles]={
+    "/data/camera_CL/Gamma_zbin0_0_7_1000to1009_w0_cleaned4035.root",
+    "/data/camera_CL/Gamma_zbin0_0_7_1130to1139_w0_cleaned4035.root",
+    "/data/camera_CL/Gamma_zbin0_0_7_1520to1529_w0_cleaned4035.root",
+    "/data/camera_CL/Gamma_zbin0_90_7_1260to1269_w0_cleaned4035.root",
+    "/data/camera_CL/Gamma_zbin0_90_7_1390to1399_w0_cleaned4035.root",
+    "/data/camera_CL/Gamma_zbin0_90_7_1650to1659_w0_cleaned4035.root"
+  };
+
+  TChain* ccut = new TChain("Events");
+  for(Int_t i = 0; i < ncutfiles; i++)
+    ccut->Add(cutName[i]);
+  ccut->SetAlias("logestenergy","log10(MHillas.fSize/15.)");
+  const Int_t nentries = Int_t(ccut->GetEntries());
+
+  MHillas*  hillas;
+  MMcEvt*   mcevt;
+  ccut->SetBranchStatus("*",0);
+  ccut->SetBranchStatus("MHillas.*",1);
+  ccut->SetBranchStatus("MMcEvt.*",1);
+  ccut->SetBranchAddress("MHillas.",&hillas);
+  ccut->SetBranchAddress("MMcEvt.",&mcevt);
+
+  TH1D* hspec = new TH1D("hspec","Spectrum",ebins,logemin,logemax);
+  hspec->Sumw2();
+
+  for(Int_t i=0;i<nentries;i++)
+    {
+      ccut->GetEntry(i);  
+      
+      // OJO!! Estimated energy will be taken directly from some input container
+      Float_t estim  = hillas->GetSize()/15.;
+
+      UInt_t effabin = UInt_t((TMath::Log10(estim)-logemin)/desub);
+      UInt_t coefbin = UInt_t((TMath::Log10(estim)-logemin)/de);
+      
+      Float_t effa  =  calc.GetEffectiveArea(effabin);
+      Float_t unfold = calc.GetCoefficient(coefbin);
+      
+      if(effa)
+	hspec->Fill(TMath::Log10(estim),unfold/effa*1e9);
+    }      
+
+  // SAVE RESULTS
+  TFile file("prueba.root","RECREATE");
+  hspec->Write();
+
+  return;
+}
