Index: trunk/MagicSoft/Mars/mtemp/mifae/Changelog =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/Changelog (revision 6355) +++ 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 6355) +++ 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 +! +! 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 +#include + +#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 || emaxDraw("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;iGetBinContent(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;iGetEntry(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(wbinFill(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;iGetBinContent(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;iGetBinContent(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 6355) +++ 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;iGetEntry(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; +}