Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.cc (revision 6363) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.cc (revision 6368) @@ -37,5 +37,5 @@ #include "MHillas.h" #include "MMcEvt.hxx" -#include "TH1D.h" +#include "TH2F.h" #include "MLog.h" @@ -46,4 +46,7 @@ using namespace std; +const Int_t fNTbins = 2; +const Double_t fTbin[fNTbins+1] = {0,10,20}; + // ------------------------------------------------------------------------- @@ -52,5 +55,6 @@ // 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,11 +68,10 @@ // 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,8 +95,8 @@ if(fCoeff) - delete [] fCoeff; + delete fCoeff; if(fEffA) - delete [] fEffA; + delete fEffA; if(fHorig) @@ -142,6 +145,6 @@ 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,11 +204,11 @@ // 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,4 +221,5 @@ 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,6 +229,6 @@ if(wbinFill(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,16 +237,18 @@ // 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; - } - } + for(Int_t j=0;jGetBinContent(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,4 +262,10 @@ 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,26 +277,31 @@ // 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;iGetBinContent(i+1)<=0) + for(Int_t j=0;jGetBinContent(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; } Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.h (revision 6363) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.h (revision 6368) @@ -7,5 +7,6 @@ class TF1; -class TH1D; +class TH1F; +class TH2F; class MHillas; class MMcEvt; @@ -18,5 +19,5 @@ 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,16 +26,16 @@ 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,12 +52,12 @@ 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();