Index: /trunk/MagicSoft/Mars/manalysis/MHillasCalc.cc =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MHillasCalc.cc (revision 1761) +++ /trunk/MagicSoft/Mars/manalysis/MHillasCalc.cc (revision 1762) @@ -84,5 +84,5 @@ } - fHillas = (MHillas*)pList->FindCreateObj(fHilName); + fHillas = (MHillas*)pList->FindCreateObj("MHillas",fHilName); if (!fHillas) return kFALSE; Index: /trunk/MagicSoft/Mars/manalysis/MHillasSrcCalc.cc =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MHillasSrcCalc.cc (revision 1761) +++ /trunk/MagicSoft/Mars/manalysis/MHillasSrcCalc.cc (revision 1762) @@ -56,6 +56,20 @@ // The default is "MHillasSrc" // +//MHillasSrcCalc::MHillasSrcCalc(const char *src, const char *hil, +// const char *name, const char *title) +// : fHillas(NULL), fSrcPos(NULL), fHillasSrc(NULL) +//{ +// fName = name ? name : gsDefName.Data(); +// fTitle = title ? title : gsDefTitle.Data(); + +// fSrcName = src; +// fHillasName = hil; +// fHillasInput = "MHillas"; +//} +// ------------------------------------------------------------------------- +// MHillasSrcCalc::MHillasSrcCalc(const char *src, const char *hil, - const char *name, const char *title) + const char *name, const char *title, + const char *hilinput) : fHillas(NULL), fSrcPos(NULL), fHillasSrc(NULL) { @@ -65,4 +79,5 @@ fSrcName = src; fHillasName = hil; + fHillasInput = hilinput; } @@ -71,5 +86,5 @@ Bool_t MHillasSrcCalc::PreProcess(MParList *pList) { - fHillas = (MHillas*)pList->FindObject("MHillas"); + fHillas = (MHillas*)pList->FindObject(fHillasInput); if (!fHillas) { Index: /trunk/MagicSoft/Mars/manalysis/MHillasSrcCalc.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MHillasSrcCalc.h (revision 1761) +++ /trunk/MagicSoft/Mars/manalysis/MHillasSrcCalc.h (revision 1762) @@ -19,4 +19,5 @@ TString fSrcName; TString fHillasName; + TString fHillasInput; Int_t fErrors; @@ -29,6 +30,10 @@ public: + // MHillasSrcCalc(const char *src="MSrcPosCam", const char *hil="MHillasSrc", + // const char *name=NULL, const char *title=NULL); + MHillasSrcCalc(const char *src="MSrcPosCam", const char *hil="MHillasSrc", - const char *name=NULL, const char *title=NULL); + const char *name=NULL, const char *title=NULL, + const char *hilinput="MHillas"); ClassDef(MHillasSrcCalc, 1) // task to calculate the source position depandant hillas parameters Index: /trunk/MagicSoft/Mars/manalysis/MSelBasic.cc =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MSelBasic.cc (revision 1762) +++ /trunk/MagicSoft/Mars/manalysis/MSelBasic.cc (revision 1762) @@ -0,0 +1,231 @@ +/* ======================================================================== *\ +! +! * +! * 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): Wolfgang Wittek 02/2003 +! +! Copyright: MAGIC Software Development, 2000-2002 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MSelBasic // +// // +// This is a task to evaluate basic cuts // +// // +// to be called after the calibration (when the number of photons is // +// available for all pixels) // +// // +///////////////////////////////////////////////////////////////////////////// + +#include "MSelBasic.h" + +#include "MParList.h" + +#include "MHillas.h" +#include "MCerPhotEvt.h" +#include "MMcEvt.hxx" +#include "MRawRunHeader.h" +#include "MGeomCam.h" +#include "MPedestalCam.h" +#include "MGeomPix.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MSelBasic); + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MSelBasic::MSelBasic(const char *name, const char *title) +{ + fName = name ? name : "MSelBasic"; + fTitle = title ? title : "Task to evaluate basic cuts"; +} + +// -------------------------------------------------------------------------- +// +// +// +// +// +Bool_t MSelBasic::PreProcess(MParList *pList) +{ + fRawRun = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); + if (!fRawRun) + { + *fLog << dbginf << "MRawRunHeader not found... aborting." << endl; + return kFALSE; + } + + fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); + if (!fMcEvt) + { + *fLog << dbginf << "MMcEvt not found... aborting." << endl; + return kFALSE; + } + + fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt"); + if (!fEvt) + { + *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl; + return kFALSE; + } + + + fCam = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!fCam) + { + *fLog << dbginf << "MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl; + return kFALSE; + } + + fPed = (MPedestalCam*)pList->FindObject("MPedestalCam"); + if (!fPed) + { + *fLog << dbginf << "MPedestalCam missing in Parameter List... aborting." << endl; + return kFALSE; + } + + memset(fErrors, 0, sizeof(fErrors)); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Evaluate basic cuts +// +// if cuts are fulfilled : return 0 +// if they are not fullfilled : skip remaining tasks for this event +// +Bool_t MSelBasic::Process() +{ + Int_t rc = 0; + + //if ( fRawRun->GetRunNumber() < 1025 ) + //{ + // rc = 1; + // return kCONTINUE; + //} + + Double_t Theta = kRad2Deg*fMcEvt->GetTelescopeTheta(); + if (Theta > 45.0 || !SwTrigger() ) + { + //*fLog << "MSelBasic::Process; Theta = " << Theta << endl; + rc = 1; + } + + fErrors[rc]++; + + return rc==0 ? kTRUE : kCONTINUE; +} +// -------------------------------------------------------------------------- +// +// Software trigger +// +// require 2 neighboring pixels (which are not in the outermost ring), +// each having at least 'minphotons' photons +// +// +Bool_t MSelBasic::SwTrigger() +{ + // minimum number of photons required + Double_t minphotons = 13.0; + + const Int_t entries = fEvt->GetNumPixels(); + + //$$$$$$$$$$$$$$$$$$ + //const Int_t nall = fPed->GetSize(); + //*fLog << "nall = " << nall << endl; + //for (Int_t id=0; id= minphotons) return kTRUE; + } + } + return kFALSE; +} + +// -------------------------------------------------------------------------- +// +// Prints some statistics about the Basic selections. +// +Bool_t MSelBasic::PostProcess() +{ + if (GetNumExecutions()==0) + return kTRUE; + + *fLog << inf << endl; + *fLog << GetDescriptor() << " execution statistics:" << endl; + *fLog << dec << setfill(' '); + *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) << (int)(fErrors[1]*100/GetNumExecutions()) << "%) Evts skipped due to: Basic selections are not fullfilled" << endl; + + *fLog << " " << fErrors[0] << " (" << (int)(fErrors[0]*100/GetNumExecutions()) << "%) Evts survived Basic selections!" << endl; + *fLog << endl; + + return kTRUE; +} Index: /trunk/MagicSoft/Mars/manalysis/MSelBasic.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MSelBasic.h (revision 1762) +++ /trunk/MagicSoft/Mars/manalysis/MSelBasic.h (revision 1762) @@ -0,0 +1,57 @@ +#ifndef MARS_MSelBasic +#define MARS_MSelBasic + +///////////////////////////////////////////////////////////////////////////// +// // +// MSelBasic // +// // +// Task to evaluate basic cuts // +// // +///////////////////////////////////////////////////////////////////////////// + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class MGeomCam; +class MPedestalCam; +class MCerPhotEvt; +class MHillas; +class MMcEvt; +class MRawRunHeader; + +class MSelBasic : public MTask +{ +private: + const MPedestalCam *fPed; // Pedestal information + const MGeomCam *fCam; // Camera Geometry + const MCerPhotEvt *fEvt; // Cerenkov Photon Event + const MMcEvt *fMcEvt; + const MRawRunHeader *fRawRun; + + Int_t fErrors[2]; + +public: + MSelBasic(const char *name=NULL, const char *title=NULL); + + Bool_t PreProcess(MParList *pList); + Bool_t Process(); + Bool_t PostProcess(); + + Bool_t SwTrigger(); + + ClassDef(MSelBasic, 0) // Task to evaluate basic cuts +}; + +#endif + + + + + + + + + + + Index: /trunk/MagicSoft/Mars/manalysis/MSelStandard.cc =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MSelStandard.cc (revision 1762) +++ /trunk/MagicSoft/Mars/manalysis/MSelStandard.cc (revision 1762) @@ -0,0 +1,157 @@ +/* ======================================================================== *\ +! +! * +! * 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): Wolfgang Wittek 02/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MSelStandard // +// // +// This is a task to evaluate the Standard Cuts // +// // +// to be called after the calculation of the image parameters // +// before the g/h separation // +// // +///////////////////////////////////////////////////////////////////////////// + +#include "MSelStandard.h" + +#include "MParList.h" + +#include "MHillas.h" +#include "MHillasSrc.h" +#include "MCerPhotEvt.h" +#include "MMcEvt.hxx" +#include "MGeomCam.h" +#include "MGeomPix.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MSelStandard); + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MSelStandard::MSelStandard(const MHillas *parhil, const MHillasSrc *parhilsrc, + const char *name, const char *title) +{ + fName = name ? name : "MSelStandard"; + fTitle = title ? title : "Task to evaluate the Standard Cuts"; + + fHil = parhil; + fHilsrc = parhilsrc; +} + +// -------------------------------------------------------------------------- +// +// +// +// +// +Bool_t MSelStandard::PreProcess(MParList *pList) +{ + fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); + if (!fMcEvt) + { + *fLog << dbginf << "MMcEvt not found... aborting." << endl; + return kFALSE; + } + + fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt"); + if (!fEvt) + { + *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl; + return kFALSE; + } + + + fCam = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!fCam) + { + *fLog << dbginf << "MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl; + return kFALSE; + } + fMm2Deg = fCam->GetConvMm2Deg(); + + *fLog << "fMm2Deg = " << fMm2Deg << endl; + + memset(fErrors, 0, sizeof(fErrors)); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Evaluate standard cuts +// +// if cuts are fulfilled : return 0 +// if they are not fullfilled : skip the remaining tasks for this event +// +// +Bool_t MSelStandard::Process() +{ + Int_t rc = 0; + + //Double_t fLength = fHil->GetLength() * fMm2Deg; + //Double_t fWidth = fHil->GetWidth() * fMm2Deg; + Double_t fDist = fHilsrc->GetDist()* fMm2Deg; + //Double_t fDelta = fHil->GetDelta() * kRad2Deg; + Double_t fSize = fHil->GetSize(); + Int_t fNumUsedPixels = fHil->GetNumUsedPixels(); + Int_t fNumCorePixels = fHil->GetNumCorePixels(); + + if ( fSize <= 60.0 || fDist< 0.4 || fDist > 1.1 + || fNumUsedPixels >= 92 || fNumCorePixels < 4) + { + //*fLog << "MSelStandard::Process; fSize, fDist, fNumUsedPixels, fNumCorePixels = " + // << fSize << ", " << fDist << ", " << fNumUsedPixels << ", " + // << fNumCorePixels << endl; + rc = 1; + } + + + fErrors[rc]++; + + return rc==0 ? kTRUE : kCONTINUE; +} + +// -------------------------------------------------------------------------- +// +// Prints some statistics about the Standard selections. +// +Bool_t MSelStandard::PostProcess() +{ + if (GetNumExecutions()==0) + return kTRUE; + + *fLog << inf << endl; + *fLog << GetDescriptor() << " execution statistics:" << endl; + *fLog << dec << setfill(' '); + *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) << (int)(fErrors[1]*100/GetNumExecutions()) << "%) Evts skipped due to: Standard selections are not fullfilled" << endl; + + *fLog << " " << fErrors[0] << " (" << (int)(fErrors[0]*100/GetNumExecutions()) << "%) Evts survived Standard selections!" << endl; + *fLog << endl; + + return kTRUE; +} Index: /trunk/MagicSoft/Mars/manalysis/MSelStandard.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MSelStandard.h (revision 1762) +++ /trunk/MagicSoft/Mars/manalysis/MSelStandard.h (revision 1762) @@ -0,0 +1,57 @@ +#ifndef MARS_MSelStandard +#define MARS_MSelStandard + +///////////////////////////////////////////////////////////////////////////// +// // +// MSelStandard // +// // +// Task to evaluate standard cuts // +// // +///////////////////////////////////////////////////////////////////////////// + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class MGeomCam; +class MCerPhotEvt; +class MHillas; +class MHillasSrc; +class MMcEvt; + +class MSelStandard : public MTask +{ +private: + const MGeomCam *fCam; // Camera Geometry + const MCerPhotEvt *fEvt; // Cerenkov Photon Event + const MMcEvt *fMcEvt; + const MHillas *fHil; + const MHillasSrc *fHilsrc; + + Double_t fMm2Deg; // conversion mm to degrees in camera + Int_t fErrors[2]; + +public: + MSelStandard(const MHillas *fHil, const MHillasSrc *fHilsrc, + const char *name=NULL, const char *title=NULL); + + Bool_t PreProcess(MParList *pList); + Bool_t Process(); + Bool_t PostProcess(); + + + ClassDef(MSelStandard, 0) // Task to evaluate standard cuts +}; + +#endif + + + + + + + + + + + Index: /trunk/MagicSoft/Mars/mfilter/MFEventSelector.cc =================================================================== --- /trunk/MagicSoft/Mars/mfilter/MFEventSelector.cc (revision 1761) +++ /trunk/MagicSoft/Mars/mfilter/MFEventSelector.cc (revision 1762) @@ -94,9 +94,12 @@ // the class description above. // -MFEventSelector::MFEventSelector(const char *name, const char *title) +MFEventSelector::MFEventSelector(const char *name, const char *title, + const char *read) : fNumTotalEvts(-1), fNumSelectEvts(-1), fSelRatio(-1), fNumSelectedEvts(0) { fName = name ? name : gsDefName.Data(); fTitle = title ? title : gsDefTitle.Data(); + + fRead = read; } @@ -127,5 +130,5 @@ return kFALSE; } - MRead *read = (MRead*)tlist->FindObject("MRead"); + MRead *read = (MRead*)tlist->FindObject(fRead); if (!read) { @@ -134,4 +137,7 @@ } fNumTotalEvts = read->GetEntries(); + + *fLog << "MFEventSelector::PreProcess; fNumTotalEvts = " + << fNumTotalEvts << endl; SetBit(kNumTotalFromFile); Index: /trunk/MagicSoft/Mars/mfilter/MFEventSelector.h =================================================================== --- /trunk/MagicSoft/Mars/mfilter/MFEventSelector.h (revision 1761) +++ /trunk/MagicSoft/Mars/mfilter/MFEventSelector.h (revision 1762) @@ -24,4 +24,5 @@ Bool_t fResult; + TString fRead; void StreamPrimitive(ofstream &out) const; @@ -38,5 +39,6 @@ public: // MFEventSelector(); - MFEventSelector(const char *name=NULL, const char *title=NULL); + MFEventSelector(const char *name=NULL, const char *title=NULL, + const char *read="MRead"); ~MFEventSelector(); Index: /trunk/MagicSoft/Mars/mhist/MHHillasExt.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHHillasExt.cc (revision 1761) +++ /trunk/MagicSoft/Mars/mhist/MHHillasExt.cc (revision 1762) @@ -57,5 +57,6 @@ // Setup four histograms for Width, Length // -MHHillasExt::MHHillasExt(const char *name, const char *title) +MHHillasExt::MHHillasExt(const char *name, const char *title, + const char *hil) : fMm2Deg(1), fUseMmScale(kTRUE) { @@ -65,4 +66,7 @@ fName = name ? name : "MHHillasExt"; fTitle = title ? title : "Container for extended Hillas histograms"; + + fHilName = hil; + //*fLog << "MHHillasExt : fHilName = " << fHilName << endl; // @@ -139,8 +143,9 @@ Bool_t MHHillasExt::SetupFill(const MParList *plist) { - TObject *obj = plist->FindObject("MHillas"); + TObject *obj = plist->FindObject(fHilName, "MHillas"); if (!obj) { - *fLog << err << dbginf << "Sorry 'MHillas' not found in parameter list... aborting." << endl; + *fLog << err << dbginf << "Sorry '" << fHilName + << "' not found in parameter list... aborting." << endl; return kFALSE; } Index: /trunk/MagicSoft/Mars/mhist/MHHillasExt.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHHillasExt.h (revision 1761) +++ /trunk/MagicSoft/Mars/mhist/MHHillasExt.h (revision 1762) @@ -25,6 +25,11 @@ Bool_t fUseMmScale; + TString fHilName; + public: - MHHillasExt(const char *name=NULL, const char *title=NULL); + //MHHillasExt(const char *name=NULL, const char *title=NULL); + MHHillasExt(const char *name=NULL, const char *title=NULL, + const char *hil="MHillas"); + ~MHHillasExt(); Index: /trunk/MagicSoft/Mars/mhist/MHMatrix.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHMatrix.cc (revision 1761) +++ /trunk/MagicSoft/Mars/mhist/MHMatrix.cc (revision 1762) @@ -16,7 +16,9 @@ ! ! -! Author(s): Thomas Bretz 2002 +! Author(s): Thomas Bretz 2002 +! Rudy Boeck 2003 ! -! Copyright: MAGIC Software Development, 2000-2002 +! Copyright: MAGIC Software Development, 2000-2003 ! ! @@ -52,4 +54,8 @@ #include +#include +#include +#include + #include "MLog.h" #include "MLogManip.h" @@ -134,4 +140,6 @@ } +// -------------------------------------------------------------------------- +// void MHMatrix::AddColumns(MDataArray *matrix) { @@ -341,4 +349,6 @@ } +// -------------------------------------------------------------------------- +// const TMatrix *MHMatrix::InvertPosDef() { @@ -541,4 +551,6 @@ } +// -------------------------------------------------------------------------- +// const TArrayI MHMatrix::GetIndexOfSortedColumn(Int_t ncol, Bool_t desc) const { @@ -558,4 +570,6 @@ } +// -------------------------------------------------------------------------- +// void MHMatrix::SortMatrixByColumn(Int_t ncol, Bool_t desc) { @@ -577,4 +591,6 @@ } +// -------------------------------------------------------------------------- +// Bool_t MHMatrix::Fill(MParList *plist, MTask *read) { @@ -617,6 +633,4 @@ // -------------------------------------------------------------------------- // -// Return a comma seperated list of all data members used in the matrix. -// This is mainly used in MTask::AddToBranchList // void MHMatrix::ReduceNumberOfRows(UInt_t numrows, const TString opt) @@ -674,4 +688,304 @@ } +// ------------------------------------------------------------------------ +// +// Define the reference matrix +// refcolumn number of the column containing the variable, for which a +// target distribution may be given; +// if refcolumn is negative the target distribution will be set +// equal to the real distribution; the events in the reference +// matrix will then be simply a random selection of the events +// in the original matrix. +// thsh histogram containing the target distribution of the variable +// nmaxevts maximum number of events in the reference matrix +// +Bool_t MHMatrix::DefRefMatrix(const Int_t refcolumn, const TH1F *thsh, + const Int_t nmaxevts) +{ + if (!fM.IsValid()) + { + *fLog << "MHMatrix::DefRefMatrix; matrix not initialized" << endl; + return kFALSE; + } + + const TH1F *fHthsh; // target distribution + TH1F *fHthd; // normalization factors + TH1F *fHth; // distribution before + TH1F *fHthaft; // distribution after renormalization + + Int_t frefcol; // number of the column (count from 0) containing + // the variable for which the target distribution is given + Int_t fnbins; // histogram parameters for thsh + Float_t ffrombin; // + Float_t ftobin; // + Float_t fdbin; // + + TArrayF fnormfac; // array of normalization factors + TMatrix fMnew; // matrix fM having the requested distribution + // and the requested number of rows; + // this is the matrix to be used in the g/h separation + + + //=================================================================== + //--------------------------------------------------------- + // Calculate normalization factors + // + + fHthsh = thsh; + + Int_t fNrows = fM.GetNrows(); + Int_t fNcols = fM.GetNcols(); + Int_t fColLwb = fM.GetColLwb(); + Int_t fRowLwb = fM.GetRowLwb(); + + // print characteristics + // + //*fLog << "MHMatrix::CalcNormFactors; Event matrix (input): " << endl; + //*fLog << " fNrows, fNcols, fColLwb, fRowLwb = " << fNrows << ", " + // << fNcols << ", " << fColLwb << ", " << fRowLwb << endl; + + + //--------------------------------------------------------- + // + // if refcol < 0 : select reference events randomly + // i.e. set the normaliztion factotrs equal to 1.0 + + if (refcolumn<0) + { + frefcol = -refcolumn; + } + else + { + frefcol = refcolumn; + } + + + //--------------------------------------------------------- + + if (fHthsh == NULL) + { + *fLog << "MHMatrix::DefRefMatrix; target distribution has to be defined" + << endl; + return kFALSE; + } + + fnbins = fHthsh->GetNbinsX(); + ffrombin = fHthsh->GetBinLowEdge(1); + ftobin = fHthsh->GetBinLowEdge(fnbins+1); + fdbin = (ftobin-ffrombin)/fnbins; + //*fLog << "CalcNormFactors; Reference column (count from 0) = " + // << frefcol << endl; + //*fLog << "CalcNormFactors; Histo params : " << fnbins << ", " << ffrombin + // <<", "<< ftobin << endl; + + + //--------------------------------------------------------- + // ===== set up the real histogram + // + fHth = new TH1F("th", "data at input", fnbins, ffrombin, ftobin); + for (Int_t j=1; j<=fNrows; j++) + { + fHth->Fill(fM(j-1,frefcol)); + //*fLog << "j, frefcol, fM(j-1,frefcol) = " << j << ", " << frefcol + // << ", " << fM(j-1,frefcol) << endl; + } + + // if refcolumn<0 set target distribution equal to the real distribution + // in order to obtain normalization factors = 1.0 + //if (refcolumn<0) + //{ + // for (Int_t j=1; j<=fnbins; j++) + // { + // float cont = fHth-> GetBinContent(j); + // fHthsh->SetBinContent(j, cont); + // } + //} + + //--------------------------------------------------------- + // ===== obtain correction factors + // + fHthd = new TH1F ("thd", "correction factors", fnbins, ffrombin, ftobin); + + Double_t cmax = 0.; + Float_t a; + Float_t b; + Float_t c; + for (Int_t j=1; j<=fnbins; j++) + { + a = fHth->GetBinContent(j); + + // if refcolumn < 0 set the correction factors equal to 1 + if ( refcolumn>=0.0 ) + b = fHthsh->GetBinContent(j); + else + b = a; + + if ( a<0.0 || b<0.0 ) + { + *fLog << "MHMatrix::DefRefMatrix; renormalization of input matrix is not possible" << endl; + *fLog << " a, b = " << a << ", " << b << endl; + return kFALSE; + } + + if (a == 0.0) + c = 0.0; + else + c = b/a; + + fHthd->SetBinContent(j, c); + if (cmax < c) cmax = c; + //*fLog << j << ", " << a << ", " << b << ", " << c << endl; + } + + //*fLog << "MHMatrix::CalcNormFactors; maximum ratio between histogram and target " + // << cmax << endl; + + if (cmax <= 0.0) + { + *fLog << "MHMatrix::CalcNormFactors; maximum ratio is LE zero" << endl; + return kFALSE; + } + + fnormfac.Set(fnbins); + Float_t minbin = 1.0; + for (Int_t j=1; j<=fnbins; j++) + { + float c = (fHthd->GetBinContent(j)) / cmax; + fnormfac[j-1] = c; + if (minbin > c && c>0) minbin = c; + fHthd->SetBinContent(j, c); + } + //*fLog << "MHMatrix::CalcNormFactors; maximum weight = 1, minimum non-zero weight = " + // << minbin <Rndm(i); + + TArrayI ind(fNrows); + TMath::Sort(fNrows, ranx.GetArray(), ind.GetArray(),kTRUE); + + //*fLog << "MHMatrix::DefRefMatrix; permuting array " <0.5) + { + cumul_weight[indref] -= 1.; + if (evtcount++ < nmaxevts) + { + fHthaft->Fill(fM(ind[ir],frefcol)); + for (Int_t ic=0; icGetBinContent(j); + if (a>0) *fLog << j << " "<< a << " "; + } + *fLog <Divide(2,2); + + th1->cd(1); + ((TH1F*)fHthsh)->Draw(); // target + + th1->cd(2); + ((TH1F*)fHth)->Draw(); // real histogram before + + th1->cd(3); + ((TH1F*)fHthd) -> Draw(); // correction factors + + th1->cd(4); + ((TH1F*)fHthaft) -> Draw(); // histogram after + + //--------------------------------------------------------- + // --- write onto output file + // + //TFile *outfile = new TFile(fileNameout, "RECREATE", ""); + //Mnew.Write(fileNameout); + //outfile->Close(); + + return kTRUE; +} + // -------------------------------------------------------------------------- // @@ -691,2 +1005,6 @@ + + + + Index: /trunk/MagicSoft/Mars/mhist/MHMatrix.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHMatrix.h (revision 1761) +++ /trunk/MagicSoft/Mars/mhist/MHMatrix.h (revision 1762) @@ -13,6 +13,10 @@ #endif +#include + class TArrayI; +class TArrayF; +class TH1F; class MTask; class MParList; @@ -32,4 +36,5 @@ MDataArray *fData; // List of data members (columns) + enum { @@ -94,4 +99,7 @@ Int_t Read(const char *name); + Bool_t DefRefMatrix (const Int_t refcolumn=-1, const TH1F *thsh=NULL, + const Int_t nmaxevts=0); + ClassDef(MHMatrix, 1) // Multidimensional Matrix to store events }; @@ -99,2 +107,3 @@ #endif +