Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 6548) +++ trunk/MagicSoft/Mars/Changelog (revision 6549) @@ -21,9 +21,26 @@ -*-*- END OF LINE -*-*- + 2005/02/16 Javier Rico + + * mhistmc/MHMcUnfoldCoeff.[cc,h], mhistmc/Makefile, + mhistmc/HistMcLinkDef.h + - Added container class holding the histograms of the coefficients + for the (iterative) unfolding. + + * mmontecarlo/MMcUnfoldCoeffCalc.[cc,h], mmontecarlo/Makefile, + mmontecarlo/MonteCarloLinkDef.h + - Added task class to compute the coefficients for the (iterative) + unfolding. + + * macros/unfoldCoeff.C + - added macro that computes the effective areas and coefficients + from a MC input file (with OriginalMC tree and MEnergyEst + branch containing the energy estimator). This may be used + as layout for spectrum computation program. + 2005/02/16 Abelardo Moralejo * mtemp/mpadova/macros/resize.C - Changed weights - 2005/02/16 Thomas Hengstebeck Index: trunk/MagicSoft/Mars/macros/unfoldCoeff.C =================================================================== --- trunk/MagicSoft/Mars/macros/unfoldCoeff.C (revision 6549) +++ trunk/MagicSoft/Mars/macros/unfoldCoeff.C (revision 6549) @@ -0,0 +1,200 @@ +/* ======================================================================== *\ +! +! * +! * 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, 12/2005 +! +! Copyright: MAGIC Software Development, 2000-2005 +! +! +\* ======================================================================== */ + +// +// Example macro on how to calculate coefficients for unfolding and +// effective areas. +//.The input is a file of the MC events surviving a certain kind +// of analysis. It must contain the branch MEnergyEst with energy +// estimator. +// It must also contain a tree called "OriginalMC" with a branch +// "MMcEvtBasic" and one entry per original simulated shower used to produce +// the events in that input file. +// +// The macro performs 2 loops. +// In the first one the histograms of the original MC events are filled, +// needed by the coefficient (MMcUnfoldCoeffCalc) and effective area +// (MMcCollectionAreaCalc) computation +// In the second one we deal with the survivors and compute effective areas +// and coefficients + +void unfoldCoeff(TString filename="/data/19990101_10002_Q_MCGamTestLZA_E_10_5.root", TString outname="coeff.root") +{ + + /////////////////////////////////////////////////////////////////// + // First loop: over all events processed by corsika + /////////////////////////////////////////////////////////////////// + MParList parlist; + MTaskList tasklist; + parlist.AddToList(&tasklist); + + // theta binning + Int_t zbins = 2; + TArrayD edges(zbins+1); + edges[0] = 0; + edges[1] = 10; + edges[2] = 20; + + MBinning binstheta("binsTheta"); + binstheta.SetEdges(edges); // Theta [deg] + + // energy binning + const Int_t ebins=10; + const Float_t emin=10; + const Float_t emax=1000; + + MBinning binsenergy("binsEnergy"); + binsenergy.SetEdgesLog(ebins,emin,emax); // Energy [GeV] + + parlist.AddToList(&binsenergy); + parlist.AddToList(&binstheta); + + // Tentative spectrum + TF1 spectrum("func","pow(x,-1.6)",2.,20000.); + + // Setup out tasks: + MReadMarsFile reader("OriginalMC", filename); + reader.DisableAutoScheme(); + + // Create unfolding coefficients object and add it to the parameter list + MHMcCollectionArea aeff; + MHMcUnfoldCoeff coeff; + parlist.AddToList(&coeff); + parlist.AddToList(&aeff); + + // Task which fills the necessary histograms in MHMcCollectionArea + MMcUnfoldCoeffCalc coeffcalc; + MMcCollectionAreaCalc aeffcalc; + + coeffcalc.SetSpectrum(&spectrum); + aeffcalc.SetSpectrum(&spectrum); + + tasklist.AddToList(&reader); + tasklist.AddToList(&coeffcalc); + tasklist.AddToList(&aeffcalc); + + // set up the loop for the processing + MEvtLoop magic; + magic.SetParList(&parlist); + + if (!magic.Eventloop()) + return; + + tasklist.PrintStatistics(); + + /////////////////////////////////////////////////////////////////// + // Second loop: now fill the histograms of the survivors + /////////////////////////////////////////////////////////////////// + + MParList parlist2; + MTaskList tasklist2; + parlist2.AddToList(&tasklist2); + parlist2.AddToList(&coeff); + parlist2.AddToList(&aeff); + + MReadMarsFile reader2("Events", filename); + reader2.DisableAutoScheme(); + + tasklist2.AddToList(&reader2); + tasklist2.AddToList(&aeffcalc); + tasklist2.AddToList(&coeffcalc); + + // + // set up the loop for the processing + // + MEvtLoop magic2; + magic2.SetParList(&parlist2); + + if (!magic2.Eventloop()) + return; + + tasklist2.PrintStatistics(); + + + + // Dummy cross-check, see if we recover the original spectrum + const UInt_t ncutfiles=1; + Char_t* cutName[ncutfiles]={"/data/19990101_10002_Q_MCGamTestLZA_E_10_5.root"}; + + TChain* ccut = new TChain("Events"); + for(Int_t i = 0; i < ncutfiles; i++) + ccut->Add(cutName[i]); + // ccut->SetAlias("logestenergy","log10(MHillas.fSize/0.18/15.)"); + ccut->SetAlias("logestenergy","log10(MEnergyEst.fEnergy)"); + ccut->SetAlias("theta","MMcEvt.fTelescopeTheta*180./3.14159"); + + const Double_t logemin = TMath::Log10(emin); + const Double_t logemax = TMath::Log10(emax); + TH1D* hspec = new TH1D("hspec","Spectrum",ebins,logemin,logemax); + hspec->Sumw2(); + ccut->Draw("logestenergy>>hspec","theta<10","goff"); + + for(Int_t i=0;iGetBinContent(i+1); + const Float_t effa = ((TH2D*) aeff.GetHistCoarse())->GetBinContent(i+1,1); + const Float_t unfold = ((TH2D*)coeff.GetHistCoeff())->GetBinContent(i+1,1); + + const Float_t euncorrval = hspec->GetBinError(i+1); + const Float_t eeffa = ((TH2D*) aeff.GetHistCoarse())->GetBinError(i+1,1); + const Float_t eunfold = ((TH2D*)coeff.GetHistCoeff())->GetBinError(i+1,1); + + Float_t corrval,ecorrval; + if(effa>0) + { + corrval = uncorrval*unfold/effa; + if(uncorrval>0 && effa>0 & unfold>0) + ecorrval = corrval*TMath::Sqrt(euncorrval/uncorrval*euncorrval/uncorrval+ + eeffa/effa*eeffa/effa+ + eunfold/unfold*eunfold/unfold); + else + ecorrval = 0; + } + else + { + corrval = 0; + ecorrval = 0; + } + hspec->SetBinContent(i+1,corrval); + hspec->SetBinError(i+1,ecorrval); + } + + // Write histogram to a file in case an output + // filename has been supplie + if (outname.IsNull()) + return; + + TFile f(outname, "recreate"); + if (!f) + return; + + coeff.GetHistAll()->Write(); + coeff.GetHistWeight()->Write(); + coeff.GetHistMcE()->Write(); + coeff.GetHistEstE()->Write(); + coeff.GetHistCoeff()->Write(); + aeff.GetHistCoarse()->Write(); + aeff.Write(); + hspec->Write(); +} Index: trunk/MagicSoft/Mars/mhistmc/HistMcLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mhistmc/HistMcLinkDef.h (revision 6548) +++ trunk/MagicSoft/Mars/mhistmc/HistMcLinkDef.h (revision 6549) @@ -17,4 +17,5 @@ #pragma link C++ class MHMcCT1CollectionArea+; #pragma link C++ class MHMcTriggerLvl2+; +#pragma link C++ class MHMcUnfoldCoeff+; #endif Index: trunk/MagicSoft/Mars/mhistmc/MHMcUnfoldCoeff.cc =================================================================== --- trunk/MagicSoft/Mars/mhistmc/MHMcUnfoldCoeff.cc (revision 6549) +++ trunk/MagicSoft/Mars/mhistmc/MHMcUnfoldCoeff.cc (revision 6549) @@ -0,0 +1,352 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MHMcUnfoldCoeff +// Container that stores the coefficients to convert from estimated to real +// (unfolded) energy. +// It is filled by MMcUnfoldCoeffCalc +// Events are weighted from any initial spectrum to a tentative spectrum +// that must be set through MMcUnfoldCoeffCalc::SetSpectrum method. +// +////////////////////////////////////////////////////////////////////////////// + +#include +#include + +#include "TH1D.h" +#include "TH2D.h" +#include "TF1.h" +#include "TCanvas.h" +#include "TArrayD.h" + +#include "MHMcUnfoldCoeff.h" + +#include "MBinning.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MHMcUnfoldCoeff); + +using namespace std; + + +// ------------------------------------------------------------------------- +// +// Constructor: Define and configure the different histograms to be filled +// +MHMcUnfoldCoeff::MHMcUnfoldCoeff(const char *name, const char *title) : + fNsubbins(20), fEmin(1.), fEmax(100000.), fNumMin(10), fFineBinning(NULL) +{ + fName = name ? name : "MHMcUnfoldCoeff"; + fTitle = title ? title : "Unfolding Coefficients vs. Theta vs. Energy"; + + // define histos + fHistAll = new TH1D(); + fHistWeight = new TH1D(); + fHistMcE = new TH2D(); + fHistEstE = new TH2D(); + fHistCoeff = new TH2D(); + + // set histo names + fHistAll->SetName("AllEvents"); + fHistWeight->SetName("Weights"); + fHistMcE->SetName("MCEnergy"); + fHistEstE->SetName("EstimatedEnergy"); + fHistCoeff->SetName("UnfoldCoeff"); + + // set histo titles + fHistAll->SetTitle("MC energy for all generated events"); + fHistWeight->SetTitle("Weights"); + fHistMcE->SetTitle("MC energy for survivors"); + fHistEstE->SetTitle("Estimate energy for survivors"); + fHistCoeff->SetTitle(fTitle); + + // histos directory + fHistAll->SetDirectory(NULL); + fHistWeight->SetDirectory(NULL); + fHistMcE->SetDirectory(NULL); + fHistEstE->SetDirectory(NULL); + fHistCoeff->SetDirectory(NULL); + + // histos style + fHistAll->UseCurrentStyle(); + fHistWeight->UseCurrentStyle(); + fHistMcE->UseCurrentStyle(); + fHistEstE->UseCurrentStyle(); + fHistCoeff->UseCurrentStyle(); + + // labels + fHistAll->SetXTitle("E [GeV]"); + fHistWeight->SetXTitle("E [GeV]"); + fHistMcE->SetXTitle("E_{MC} [GeV]"); + fHistEstE->SetXTitle("E_{est} [GeV]"); + fHistCoeff->SetXTitle("E_{est} [GeV]"); + + fHistWeight->SetYTitle("weight"); + fHistMcE->SetYTitle("\\theta [\\circ]"); + fHistEstE->SetYTitle("\\theta [\\circ]"); + fHistCoeff->SetYTitle("\\theta [\\circ]"); + + fHistMcE->SetZTitle("weighted entries"); + fHistEstE->SetZTitle("weighted entries"); + fHistCoeff->SetZTitle("Coefficient"); +} + +// ------------------------------------------------------------------------- +// +// Destructor: Delete all the histograms +// +MHMcUnfoldCoeff::~MHMcUnfoldCoeff() +{ + delete fHistAll; + delete fHistWeight; + delete fHistMcE; + delete fHistEstE; + delete fHistCoeff; + + if(fFineBinning) + delete fFineBinning; +} + +// -------------------------------------------------------------------------- +// +// Set the histograms binning. The coarse binning has to be provided +// from outside (in the parameter list) whereas the fine binning for +// energy is computed from the former. +// The subbinning is set logarithmicaly (linearly) for a logarithmic +// (linear) coarse binning. +// It is extended down to fEmin and up to fEmax +// The number of subbins per coarse bins may be changed using +// SetNsubbins() method +// The maximum and minimum energies of the fine binning histograms +// may be changed with the SetEmin() and SetEmax() methods +// +void MHMcUnfoldCoeff::SetCoarseBinnings(const MBinning &binsEnergy, + const MBinning &binsTheta) +{ + // set histo binning (for coarse bin histos) + MH::SetBinning(fHistMcE,&binsEnergy,&binsTheta); + MH::SetBinning(fHistEstE,&binsEnergy,&binsTheta); + MH::SetBinning(fHistCoeff,&binsEnergy,&binsTheta); + + // good errors + fHistMcE->Sumw2(); + fHistEstE->Sumw2(); + + // define fine binning from coarse one + const Int_t nbins = binsEnergy.GetNumBins(); + const TArrayD edges = binsEnergy.GetEdgesD(); + const Bool_t isLinear = binsEnergy.IsLinear(); + + // compute bins needed to extend down to fEmin and up to fEmax + Double_t dedown, deup; + Int_t binsdown,binsup; + MBinning dbin; + MBinning ubin; + if(isLinear) + { + dedown = (edges[1]-edges[0])/fNsubbins; + deup = (edges[nbins]-edges[nbins-1])/fNsubbins; + binsdown = (fEminedges[nbins]) ? Int_t((fEmax-edges[nbins])/deup) : 0; + if(binsdown) dbin.SetEdges(binsdown,fEmin,edges[0]); + if(binsup) ubin.SetEdges(binsup,edges[nbins],fEmax); + } + else + { + dedown = (TMath::Log10(edges[1])-TMath::Log10(edges[0]))/fNsubbins; + deup = (TMath::Log10(edges[nbins])-TMath::Log10(edges[nbins-1]))/fNsubbins; + binsdown = (fEminedges[nbins])? Int_t((TMath::Log10(fEmax)-TMath::Log10(edges[nbins]))/deup) : 0; + if(binsdown) dbin.SetEdgesLog(binsdown,fEmin,edges[0]); + if(binsup) ubin.SetEdgesLog(binsup,edges[nbins],fEmax); + } + + // fill the subins' edges + TArrayD fineedges(binsdown+nbins*fNsubbins+binsup+1); + + for(Int_t i=0;iSetEdges(fineedges); + + // apply fine binning to histograms that need it + fFineBinning->Apply(*fHistAll); + fFineBinning->Apply(*fHistWeight); +} + +// ------------------------------------------------------------------------- +// +// Compute the weights for a particular tentative spectrum. +// For each energy (fine) bin we compute it as the value of the input function +// divided by the number of entries in the actual energy histogram. +// First and last bin are not filled since they could be biased +// +void MHMcUnfoldCoeff::ComputeWeights(TF1* spectrum) +{ + Bool_t first=kTRUE; + Int_t lastb=0; + for(Int_t i=0;iGetNbinsX();i++) + { + const Float_t denom = fHistAll->GetBinContent(i+1); // number of events in initial spectrum + const Float_t ew = fHistAll->GetBinCenter(i+1); // energy associated to the bin + const Float_t numer = spectrum->Eval(ew); // number of events for the required spectrum + + if(denom>fNumMin) + { + // do not fill it if it is the first one + if(first) + { + fHistWeight->SetBinContent(i+1,-1); + first=kFALSE; + } + else + { + fHistWeight->SetBinContent(i+1,numer/denom); + lastb=i+1; + } + } + else + fHistWeight->SetBinContent(i+1,-1); + } + + //remove last filled bin + if(lastb) + fHistWeight->SetBinContent(lastb,-1); +} + +// -------------------------------------------------------------- +// +// Compute the coefficients used for the (iterative) unfolding +// The coefficients are the ratio between the contents of the +// mc energy and estimate energy histograms (filled with weighted +// events +// Errors are computed as if estimated and MC energy histos where +// uncorrelated +// +void MHMcUnfoldCoeff::ComputeCoefficients() +{ + for(Int_t j=0;jGetNbinsY();j++) + for(Int_t i=0;iGetNbinsX();i++) + { + const Float_t num = fHistMcE->GetBinContent(i+1,j+1); + const Float_t den = fHistEstE->GetBinContent(i+1,j+1); + + if(den>0) + { + const Float_t cf = num/den; + fHistCoeff->SetBinContent(i+1,j+1,cf); + + if(num>0) + { + const Float_t ernum = fHistMcE->GetBinError(i+1,j+1); + const Float_t erden = fHistEstE->GetBinError(i+1,j+1); + const Float_t ercf = cf*(TMath::Sqrt(ernum/num*ernum/num+erden/den*erden/den)); + fHistCoeff->SetBinError(i+1,j+1,ercf); + } + else + fHistCoeff->SetBinError(i+1,j+1,0); + } + else + { + *fLog << warn << "MHMcUnfoldCoeff::ComputeCoefficients Warning: energy bin " << i << ", thetabin " << j << " has no survivors, setting coefficient to 0" << endl; + fHistCoeff->SetBinContent(i+1,j+1,0.); + fHistCoeff->SetBinError(i+1,j+1,0.); + } + } +} + +// -------------------------------------------------------------------------- +// +// Fill data into the histogram which contains all showers +// +void MHMcUnfoldCoeff::FillAll(Double_t energy) +{ + fHistAll->Fill(energy); +} + +// -------------------------------------------------------------------------- +// +// Fill data into the histograms which contain survivors +// each event is introduced with a weight depending on its (MC) energy +// +void MHMcUnfoldCoeff::FillSel(Double_t mcenergy,Double_t estenergy,Double_t theta) +{ + // bin of the corresponding weight + const Int_t wbin = fFineBinning->FindHiEdge(mcenergy); + + if(wbin>0) + { + // weight + const Double_t weight = fHistWeight->GetBinContent(wbin); + + if(weight>0) + { + fHistMcE->Fill(mcenergy,theta,weight); + fHistEstE->Fill(estenergy,theta,weight); + } + // else + // *fLog << warn << "MHMcUnfoldCoeff::FillSel Warning: event with energy " << mcenergy << " has no computed weight (lack of MC statistics), skipping" << endl; + } + // the event has an energy out of the considered range + else + *fLog << warn << "MHMcUnfoldCoeff::FillSel Warning: event with energy " << mcenergy << " has energy out of bounds, skipping" << endl; +} + +// -------------------------------------------------------------------------- +// +// Draw +// +void MHMcUnfoldCoeff::Draw(Option_t* option) +{ + TCanvas *c1 = new TCanvas(); + c1->SetLogx(); + c1->SetLogy(); + c1->SetGridx(); + c1->SetGridy(); + + fHistCoeff->Draw(""); +} + Index: trunk/MagicSoft/Mars/mhistmc/MHMcUnfoldCoeff.h =================================================================== --- trunk/MagicSoft/Mars/mhistmc/MHMcUnfoldCoeff.h (revision 6549) +++ trunk/MagicSoft/Mars/mhistmc/MHMcUnfoldCoeff.h (revision 6549) @@ -0,0 +1,59 @@ +#ifndef MARS_MHMcUnfoldCoeff +#define MARS_MHMcUnfoldCoeff + +#ifndef MARS_MH +#include "MH.h" +#endif + +class TF1; +class TH1D; +class TH2D; +class MBinning; + +class MHMcUnfoldCoeff : public MH +{ + private: + + TH1D* fHistAll; // histogram of all events (fine bins) + TH1D* fHistWeight; // histogram of weights (fine bins) + + TH2D* fHistMcE; // histogram of MC energy for survivors (coarse bins) + TH2D* fHistEstE; // histogram of estimated energy for survivors (coarse bins) + TH2D* fHistCoeff; // histogram of coefficients for iterative unfolding + + Int_t fNsubbins; // number of subbins per coarse bin + Double_t fEmin; // [GeV] absolute minimum energy where weights are computed + Double_t fEmax; // [GeV] absolute maximum energy where weights are computed + Int_t fNumMin; // minimum number of events in a fine bin required to compute the weight + + MBinning* fFineBinning; // fine binning (used for weights) + + public: + + MHMcUnfoldCoeff(const char *name=NULL, const char *title=NULL); + + virtual ~MHMcUnfoldCoeff(); + + void FillAll(Double_t energy); + void FillSel(Double_t mcenergy,Double_t estenergy,Double_t theta); + void Draw(Option_t* option=""); + void ComputeWeights(TF1* spectrum); + void ComputeCoefficients(); + + void SetCoarseBinnings(const MBinning& binsEnergy,const MBinning& binsTheta); + void SetNsubbins(Int_t n) {fNsubbins=n;} + void SetEmax(Double_t emax) {fEmax=emax;} + void SetEmin(Double_t emin) {fEmin=emin;} + void SetNumMin(Int_t n) {fNumMin=n;} + + const TH1D* GetHistAll() const {return fHistAll;} + const TH1D* GetHistWeight() const {return fHistWeight;} + const TH2D* GetHistMcE() const {return fHistMcE;} + const TH2D* GetHistEstE() const {return fHistEstE;} + const TH2D* GetHistCoeff() const {return fHistCoeff;} + + ClassDef(MHMcUnfoldCoeff, 1) // Data Container to store Unfolding Coefficients +}; + +#endif + Index: trunk/MagicSoft/Mars/mhistmc/Makefile =================================================================== --- trunk/MagicSoft/Mars/mhistmc/Makefile (revision 6548) +++ trunk/MagicSoft/Mars/mhistmc/Makefile (revision 6549) @@ -33,5 +33,6 @@ MHMcEnergyMigration.cc \ MHMcCT1CollectionArea.cc \ - MHMcTriggerLvl2.cc + MHMcTriggerLvl2.cc \ + MHMcUnfoldCoeff.cc ############################################################ Index: trunk/MagicSoft/Mars/mmontecarlo/MMcUnfoldCoeffCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcUnfoldCoeffCalc.cc (revision 6549) +++ trunk/MagicSoft/Mars/mmontecarlo/MMcUnfoldCoeffCalc.cc (revision 6549) @@ -0,0 +1,206 @@ +/* ======================================================================== *\ +! +! * +! * 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 2/2005 +! +! Copyright: MAGIC Software Development, 2000-2005 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MMcUnfoldCoeffCalc +// Task to compute the coefficients for energy unfolding (to be used +// typically in an iterative process). +// +// Input containers: +// MMcEvt +// MMcEvtBasic +// MEnergyEst +// binsTheta [MBinning] +// binsEnergy [MBinning] +// +// Output containers: +// MHMcUnfodCoeff +// +// The binning for energy and theta are looked in the parameter list +// The tentative spectrum MUST be provided using SetSpectrum, otherwise +// it is used for default pow(energy,-1.6) [~Crab at 1 TeV] +// +// This task is supposed to be in a two-looped macro/program: +// In the first one MMcEvtBasic is read in order to compute the +// weights to convert from original to tentative spectrum. +// In the second one MMcEvt and MEnergyEst are read and fill the +// MC and estimated energy that are divided in the postprocess to +// compute the coefficients +// +////////////////////////////////////////////////////////////////////////////// + +#include "TF1.h" + +#include "MMcUnfoldCoeffCalc.h" +#include "MHMcUnfoldCoeff.h" + +#include "MParList.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MMcEvt.hxx" +#include "MMcEvtBasic.h" +#include "MEnergyEst.h" + +ClassImp(MMcUnfoldCoeffCalc); + +using namespace std; + +//////////////////////////////////////////////////////////////////////////////// +// +// Constructor +// +MMcUnfoldCoeffCalc::MMcUnfoldCoeffCalc(const char* name, const char* title) : + fSpectrum(NULL) +{ + fName = name ? name : "MMcUnfoldCoeffCalc"; + fTitle = title ? title : "Task to calculate the unfolding coefficients"; +} + +//////////////////////////////////////////////////////////////////////////////// +// +// PreProcess. +// Create MHMcUnfoldCoeff object if necessary. Search for other necessary +// containers: if MMcEvt is found, it means we are in the loop over the Events +// tree, and so we must fill the histogram MHMcUnfoldCoeff::fHistMcE and +// MHMcUnfoldCoeff::fHistEstE (using MHMcUnfoldCoeff::FillSel()). +// If MMcEvt is not found, it means that we are in the loop over the +// "OriginalMC" tree, containing all the original Corsika events +// produced, and hence we must fill the histogram fHistAll through +// MHMcUnfoldCoeff::FillAll() +// +Int_t MMcUnfoldCoeffCalc::PreProcess (MParList *pList) +{ + fUnfoldCoeff = (MHMcUnfoldCoeff*)pList->FindCreateObj("MHMcUnfoldCoeff"); + + if (!fBinsTheta) + { + fBinsTheta = (MBinning*) pList->FindObject("binsTheta"); + if (!fBinsTheta) + { + *fLog << err << "Coarse Theta binning not found... Aborting." + << endl; + return kFALSE; + } + } + + if (!fBinsEnergy) + { + fBinsEnergy = (MBinning*) pList->FindObject("binsEnergy"); + if (!fBinsEnergy) + { + *fLog << err << "Coarse Energy binning not found... Aborting." + << endl; + return kFALSE; + } + } + + fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); + if (fMcEvt) + { + *fLog << inf << "MMcEvt found... Filling MHMcUnfoldCoeff.fHistSel..." << endl; + + fEnergyEst = (MEnergyEst*)pList->FindObject("MEnergyEst"); + if(!fEnergyEst) + { + *fLog << err << "MEnergyEst not found... aborting..." << endl; + return kFALSE; + } + + return kTRUE; + } + + *fLog << inf << "MMcEvt not found... looking for MMcEvtBasic..." << endl; + + fMcEvtBasic = (MMcEvtBasic*) pList->FindObject("MMcEvtBasic"); + + if (fMcEvtBasic) + { + fUnfoldCoeff->SetCoarseBinnings(*fBinsEnergy,*fBinsTheta); + + *fLog << inf << "MMcEvtBasic found... Filling MHMcUnfoldCoeff.fHistAll..." << endl; + return kTRUE; + } + + *fLog << err << "MMcEvtBasic not found. Aborting..." << endl; + + return kFALSE; +} + +//////////////////////////////////////////////////////////////////////////////// +// +// Depending on whether fMcEvt or fMcEvtBasic are available, fill +// MHMcUnfoldCoeff::fHistAll, else, if fMcEvt is available, fill +// MHMcUnfoldCoeff::fHistMcE and MHMcUnfoldCoeff::fHistEstE +// +Int_t MMcUnfoldCoeffCalc::Process() +{ + Double_t mcenergy = fMcEvt ? fMcEvt->GetEnergy() : fMcEvtBasic->GetEnergy(); + Double_t estenergy = fEnergyEst ? fEnergyEst->GetEnergy() : 0; + Double_t theta = fMcEvt? fMcEvt->GetTelescopeTheta()*TMath::RadToDeg() : 0; // in deg + + if (fMcEvt) + fUnfoldCoeff->FillSel(mcenergy, estenergy, theta); + else + fUnfoldCoeff->FillAll(mcenergy); + + return kTRUE; +} + +/////////////////////////////////////////////////////////////////////////////// +// +// Postprocess. +// If both fHistMcE and fHistEstE are already filled, calculate +// the coefficients. +// Else it means we still have to run one more loop, and then the +// weights are computed +// +Int_t MMcUnfoldCoeffCalc::PostProcess() +{ + if(((TH2D*)fUnfoldCoeff->GetHistMcE())->GetEntries() > 0 && + ((TH2D*)fUnfoldCoeff->GetHistEstE())->GetEntries() > 0) + { + fUnfoldCoeff->ComputeCoefficients(); + return kTRUE; + } + + if(((TH1D*)fUnfoldCoeff->GetHistAll())->GetEntries() <= 0) + { + *fLog << err << "Postprocess reached with no histogram computed. Aborting" << endl; + return kFALSE; + } + + + if(fSpectrum) + fUnfoldCoeff->ComputeWeights(fSpectrum); + else + { + // default spectrum in case no other one is provided + TF1 spectrum("func","pow(x,-1.6)",2.,20000.); + fUnfoldCoeff->ComputeWeights(&spectrum); + } + + return kTRUE; +} Index: trunk/MagicSoft/Mars/mmontecarlo/MMcUnfoldCoeffCalc.h =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcUnfoldCoeffCalc.h (revision 6549) +++ trunk/MagicSoft/Mars/mmontecarlo/MMcUnfoldCoeffCalc.h (revision 6549) @@ -0,0 +1,46 @@ +#ifndef MARS_MMcUnfoldCoeffCalc +#define MARS_MMcUnfoldCoeffCalc + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#include + +class MParList; +class MMcEvt; +class MMcEvtBasic; +class MEnergyEst; +class MBinning; +class MHMcUnfoldCoeff; + +class MMcUnfoldCoeffCalc : public MTask +{ + private: + const MMcEvt* fMcEvt; + const MMcEvtBasic* fMcEvtBasic; + const MEnergyEst* fEnergyEst; + + MBinning* fBinsTheta; // coarse zenith angle binning + MBinning* fBinsEnergy; // coarse energy binning + + TF1* fSpectrum; // Tentative energy spectrum. + + MHMcUnfoldCoeff* fUnfoldCoeff; // container holding coefficients histogram + + TString fObjName; + + Int_t PreProcess(MParList *pList); + Int_t Process(); + Int_t PostProcess(); + + public: + MMcUnfoldCoeffCalc(const char *name = NULL, const char *title = NULL); + + void SetSpectrum(TF1 *f) { fSpectrum = f; } + + ClassDef(MMcUnfoldCoeffCalc, 0) // Task to calculate the coefficients for unfolding +}; + +#endif + Index: trunk/MagicSoft/Mars/mmontecarlo/Makefile =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/Makefile (revision 6548) +++ trunk/MagicSoft/Mars/mmontecarlo/Makefile (revision 6549) @@ -21,5 +21,5 @@ INCLUDES = -I. -I../mbase -I../mmc -I../mhbase -I../mhist -I../mhistmc \ -I../mgeom -I../manalysis -I../mtools -I../mfileio \ - -I../mfilter -I../mdata -I../mfbase + -I../mfilter -I../mdata -I../mfbase SRCFILES = MMcCollectionAreaCalc.cc \ @@ -28,5 +28,6 @@ MMcTriggerRateCalc.cc \ MMcEnergyEst.cc \ - MMcWeightEnergySpecCalc.cc + MMcWeightEnergySpecCalc.cc \ + MMcUnfoldCoeffCalc.cc ############################################################ Index: trunk/MagicSoft/Mars/mmontecarlo/MonteCarloLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MonteCarloLinkDef.h (revision 6548) +++ trunk/MagicSoft/Mars/mmontecarlo/MonteCarloLinkDef.h (revision 6549) @@ -11,4 +11,5 @@ #pragma link C++ class MMcEnergyEst+; #pragma link C++ class MMcWeightEnergySpecCalc+; +#pragma link C++ class MMcUnfoldCoeffCalc+; #endif