Changeset 2225


Ignore:
Timestamp:
06/24/03 10:15:42 (22 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r2224 r2225  
    11                                                 -*-*- END OF LINE -*-*-
     2
     3 2003/06/24: Thomas Bretz
     4
     5   * manalysis/MCT1SupercutsCalc.[h,cc]:
     6     - implemented Mapping for Supercuts
     7     - changed data member arrays to TArrayD
     8     
     9   * manalysis/MEnergyEstParam.h:
     10     - added a comment
     11
     12   * mhist/MHHadronness.[h,cc]:
     13     - implemented mapping
     14     - implemented calculating Acc_g/sqrt(Acc_h) for filtercuts
     15
     16
    217
    318 2003/06/23: Thomas Bretz

  • trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.cc

    r2206 r2225  
    3333#include <fstream>
    3434
    35 #include "MLog.h"
    36 #include "MLogManip.h"
    37 
    3835#include "MParList.h"
    3936#include "MHillasExt.h"
    … …  
    4340#include "MGeomCam.h"
    4441#include "MHadronness.h"
     42#include "MHMatrix.h"
     43
     44#include "MLog.h"
     45#include "MLogManip.h"
    4546
    4647ClassImp(MCT1SupercutsCalc);
    … …  
    5051void MCT1SupercutsCalc::InitParams()
    5152{
     53    fLengthUp.Set(8);
     54    fLengthLo.Set(8);
     55    fWidthUp.Set(8);
     56    fWidthLo.Set(8);
     57    fDistUp.Set(8);
     58    fDistLo.Set(8);
     59    fAsymUp.Set(8);
     60    fAsymLo.Set(8);
     61    fAlphaUp.Set(8);
     62
    5263    //---------------------------------
    5364    // cut parameters
    … …  
    151162Int_t MCT1SupercutsCalc::PreProcess(MParList *pList)
    152163{
     164    MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam");
     165    if (!cam)
     166    {
     167        *fLog << err << "MGeomCam (Camera Geometry) not found... aborting." << endl;
     168        return kFALSE;
     169    }
     170
     171    fMm2Deg = cam->GetConvMm2Deg();
     172
     173    fHadronness = (MHadronness*)pList->FindCreateObj("MHadronness", fHadronnessName);
     174    if (!fHadronness)
     175    {
     176        *fLog << err << fHadronnessName << " [MHadronness] not found... aborting." << endl;
     177        return kFALSE;
     178    }
     179
     180    if (fMatrix)
     181        return kTRUE;
     182
    153183    fHil = (MHillas*)pList->FindObject(fHilName, "MHillas");
    154184    if (!fHil)
    … …  
    172202    }
    173203
    174     MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam");
    175     if (!cam)
    176     {
    177         *fLog << err << "MGeomCam (Camera Geometry) not found... aborting." << endl;
    178         return kFALSE;
    179     }
    180 
    181     fMm2Deg = cam->GetConvMm2Deg();
    182 
    183     fHadronness = (MHadronness*)pList->FindCreateObj("MHadronness", fHadronnessName);
    184     if (!fHadronness)
    185     {
    186         *fLog << err << fHadronnessName << " [MHadronness] not found... aborting." << endl;
    187         return kFALSE;
    188     }
    189 
    190 
    191204    return kTRUE;
    192205}
    … …  
    196209// Calculation of upper and lower limits
    197210//
    198 Double_t MCT1SupercutsCalc::CtsMCut(Double_t *a,  Double_t ls, Double_t ct,
     211Double_t MCT1SupercutsCalc::CtsMCut(
     212#if ROOT_VERSION_CODE > ROOT_VERSION(3,05,00)
     213const
     214#endif
     215                                    TArrayD &a,  Double_t ls, Double_t ct,
    199216                                    Double_t ls2, Double_t dd2)
    200217{
    … …  
    209226    //     ct: Cos(ZA.) - dNOMCOSZA
    210227    //    dd2: DIST^2
    211 
    212228    const Double_t limit =
    213229        a[0] + a[1] * dd2 + a[2] * ct  +
    … …  
    225241
    226242    return limit;
     243}
     244
     245// --------------------------------------------------------------------------
     246//
     247// Returns the mapped value from the Matrix
     248//
     249Double_t MCT1SupercutsCalc::GetVal(Int_t i) const
     250{
     251    return (*fMatrix)[fMap[i]];
     252}
     253
     254// --------------------------------------------------------------------------
     255//
     256// You can use this function if you want to use a MHMatrix instead of the
     257// given containers. This function adds all necessary columns to the
     258// given matrix. Afterward you should fill the matrix with the corresponding
     259// data (eg from a file by using MHMatrix::Fill). If you now loop
     260// through the matrix (eg using MMatrixLoop) MEnergyEstParam::Process
     261// will take the values from the matrix instead of the containers.
     262//
     263void MCT1SupercutsCalc::InitMapping(MHMatrix *mat)
     264{
     265    if (fMatrix)
     266        return;
     267
     268    fMatrix = mat;
     269
     270    fMap[0] = fMatrix->AddColumn("MMcEvt.fTelescopeTheta");
     271    fMap[1] = fMatrix->AddColumn("MHillas.fWidth");
     272    fMap[2] = fMatrix->AddColumn("MHillas.fLength");
     273    fMap[3] = fMatrix->AddColumn("MHillas.fSize");
     274    fMap[4] = fMatrix->AddColumn("MHillas.fMeanX");
     275    fMap[5] = fMatrix->AddColumn("MHillas.fMeanY");
     276    fMap[6] = fMatrix->AddColumn("MHillasSrc.fDist");
    227277}
    228278
    … …  
    240290    const Double_t kNomCosZA   = 1.0;
    241291
    242     const Double_t newdist = fHilSrc->GetDist() * fMm2Deg;
    243 
    244     const Double_t xbar    = fHil->GetMeanX();
    245     const Double_t ybar    = fHil->GetMeanY();
    246     const Double_t dist    = sqrt(xbar*xbar + ybar*ybar) * fMm2Deg;
    247     const Double_t dd2     = dist * dist;
    248 
    249     const Double_t dmls    = log(fHil->GetSize()) - kNomLogSize;
     292    const Double_t theta   = fMatrix ? GetVal(0) : fMcEvt->GetTelescopeTheta();
     293    const Double_t width0  = fMatrix ? GetVal(1) : fHil->GetWidth();
     294    const Double_t length0 = fMatrix ? GetVal(2) : fHil->GetLength();
     295    const Double_t size    = fMatrix ? GetVal(3) : fHil->GetSize();
     296    const Double_t meanx   = fMatrix ? GetVal(4) : fHil->GetMeanX();
     297    const Double_t meany   = fMatrix ? GetVal(5) : fHil->GetMeanY();
     298    const Double_t dist0   = fMatrix ? GetVal(6) : fHilSrc->GetDist();
     299
     300    const Double_t newdist = dist0 * fMm2Deg;
     301
     302    const Double_t dist2   = meanx*meanx + meany*meany;
     303    const Double_t dd2     = dist2*fMm2Deg;
     304    const Double_t dist    = sqrt(dist2) * fMm2Deg;
     305
     306    const Double_t dmls    = log(size) - kNomLogSize;
    250307    const Double_t dmls2   = dmls * dmls;
    251308
    252     const Double_t dmcza   = cos(fMcEvt->GetTelescopeTheta()) - kNomCosZA;
    253 
    254     const Double_t length  = fHil->GetLength() * fMm2Deg;
    255     const Double_t width   = fHil->GetWidth()  * fMm2Deg;
    256 
    257     //*fLog << "MCT1SupercutsCalc::Process; dmls, dmcza, dmls2, dd2 = "
    258     //      << dmls << ",  " << dmcza << ",  " << dmls2 << ",  "
    259     //      << dd2 << endl;
    260 
    261     //*fLog << "MCT1SupercutsCalc::Process; newdist, dist, length, width = "
    262     //      << newdist << ",  " << dist << ",  " << length << ",  "
    263     //      << width << endl;
     309    const Double_t dmcza   = cos(theta) - kNomCosZA;
     310
     311    const Double_t length  = length0 * fMm2Deg;
     312    const Double_t width   = width0  * fMm2Deg;
    264313
    265314    if (newdist < 1.05                                         &&

  • trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.h

    r2206 r2225  
    22#define MARS_MCT1SupercutsCalc
    33
    4 /////////////////////////////////////////////////////////////////////////////
    5 //                                                                         //
    6 // MCT1SupercutsCalc                                                       //
    7 //                                                                         //
    8 /////////////////////////////////////////////////////////////////////////////
    9 
    104#ifndef MARS_MFilter
    115#include "MFilter.h"
     6#endif
     7
     8#ifndef ROOT_TArrayD
     9#include <TArrayD.h>
    1210#endif
    1311
    … …  
    1917class MGeomCam;
    2018class MHadronness;
     19class MHMatrix;
    2120
    2221
    … …  
    3534    Double_t    fMm2Deg;
    3635
     36    Int_t     fMap[7];
     37    MHMatrix *fMatrix;
     38
    3739    //---------------------------------
    3840    // cut parameters
    39     Double_t fLengthUp[8];
    40     Double_t fWidthUp[8];
    41     Double_t fDistUp[8];
    42     Double_t fLengthLo[8];
    43     Double_t fWidthLo[8];
    44     Double_t fDistLo[8];
    45     Double_t fAsymUp[8];
    46     Double_t fAsymLo[8];
    47     Double_t fAlphaUp[8];
     41    TArrayD fLengthUp;
     42    TArrayD fLengthLo;
     43    TArrayD fWidthUp;
     44    TArrayD fWidthLo;
     45    TArrayD fDistUp;
     46    TArrayD fDistLo;
     47    TArrayD fAsymUp;
     48    TArrayD fAsymLo;
     49    TArrayD fAlphaUp;
    4850    //---------------------------------
    4951
    … …  
    5355    Int_t Process();
    5456
     57    void Set(TArrayD &a, const TArrayD &b) { if (a.GetSize()==b.GetSize()) a=b; }
     58    Double_t GetVal(Int_t i) const;
     59
     60    Double_t CtsMCut(
     61#if ROOT_VERSION_CODE > ROOT_VERSION(3,05,00)
     62const
     63#endif
     64                     TArrayD &a, Double_t ls, Double_t ct,
     65                     Double_t ls2, Double_t dd2);
     66
    5567public:
    5668    MCT1SupercutsCalc(const char *hilname="MHillas",
    … …  
    5870                      const char *name=NULL, const char *title=NULL);
    5971
    60     Double_t CtsMCut(Double_t *a, Double_t ls, Double_t ct,
    61                      Double_t ls2, Double_t dd2);
    62 
    6372    void SetHadronnessName(const TString name) { fHadronnessName = name; }
    6473    TString GetHadronnessName() const { return fHadronnessName; }
     74
     75    void InitMapping(MHMatrix *mat);
     76    void StopMapping() { InitMapping(NULL); }
     77
     78    void SetLengthUp(const TArrayD &d) { Set(fLengthUp, d); }
     79    void SetLengthLo(const TArrayD &d) { Set(fLengthLo, d); }
     80    void SetWidthUp (const TArrayD &d) { Set(fWidthUp,  d); }
     81    void SetWidthLo (const TArrayD &d) { Set(fWidthLo,  d); }
     82    void SetDistUp  (const TArrayD &d) { Set(fDistUp,   d); }
     83    void SetDistLo  (const TArrayD &d) { Set(fDistLo,   d); }
     84    void SetAsymUp  (const TArrayD &d) { Set(fAsymUp,   d); }
     85    void SetAsymLo  (const TArrayD &d) { Set(fAsymLo,   d); }
     86    void SetAlphaUp (const TArrayD &d) { Set(fAlphaUp,  d); }
     87
     88    const TArrayD &GetLengthUp() const { return fLengthUp; }
     89    const TArrayD &GetLengthLo() const { return fLengthLo; }
     90    const TArrayD &GetWidthUp() const  { return fWidthUp; }
     91    const TArrayD &GetWithLo() const   { return fWidthLo; }
     92    const TArrayD &GetDistUp() const   { return fDistUp; }
     93    const TArrayD &GetDistLo() const   { return fDistLo; }
     94    const TArrayD &GetAsymUp() const   { return fAsymUp; }
     95    const TArrayD &GetAsymLo() const   { return fAsymLo; }
     96    const TArrayD &GetAlphaUp() const  { return fAlphaUp; }
    6597
    6698    ClassDef(MCT1SupercutsCalc, 0) // A class to evaluate the Supercuts

  • trunk/MagicSoft/Mars/manalysis/MEnergyEstParam.h

    r2206 r2225  
    1919{
    2020private:
    21     Int_t     fMap[100];
     21    Int_t     fMap[100]; // FIXME!
    2222
    2323    MHMatrix *fMatrix;

  • trunk/MagicSoft/Mars/mhist/MHHadronness.cc

    r2180 r2225  
    5656// MFillH.
    5757//
     58// If you are using filtercuts which gives you only two discrete values
     59// of the hadronness (0.25 and 0.75) you may want to use MHHadronness(2).
     60// If you request Q05() from such a MHHadronness instance you will get
     61// Acc_g/sqrt(Acc_h)
     62//
    5863////////////////////////////////////////////////////////////////////////////
    5964#include "MHHadronness.h"
    … …  
    6570#include <TMarker.h>
    6671
     72#include "MParList.h"
     73#include "MBinning.h"
     74#include "MHMatrix.h"
     75#include "MHadronness.h"
     76
    6777#include "MLog.h"
    6878#include "MLogManip.h"
    6979
    70 #include "MParList.h"
    71 #include "MBinning.h"
    72 #include "MHadronness.h"
    73 
    7480#include "MMcEvt.hxx"
    7581
    … …  
    8490//
    8591MHHadronness::MHHadronness(Int_t nbins, const char *name, const char *title)
     92    : fMatrix(NULL)
    8693{
    8794    //
    … …  
    95102    fGraph->SetMarkerStyle(kFullDotSmall);
    96103
    97     fGhness = new TH1D("Ghness", "H. Gammas",  nbins, 0, 1);
    98     fPhness = new TH1D("Phness", "H. Hadrons", nbins, 0, 1);
     104    fGhness = new TH1D("Ghness", "Acceptance vs. Hadronness (Gammas)",  nbins, 0, 1);
     105    fPhness = new TH1D("Phness", "Acceptance vs. Hadronness (Hadrons)", nbins, 0, 1);
    99106    fGhness->SetXTitle("Hadronness");
    100107    fPhness->SetXTitle("Hadronness");
    101     fGhness->SetYTitle("Counts");
    102     fPhness->SetYTitle("Counts");
     108    fGhness->SetYTitle("Acceptance");
     109    fPhness->SetYTitle("Acceptance");
    103110    fPhness->SetLineColor(kRed);
    104111    fGhness->SetDirectory(NULL);
    105112    fPhness->SetDirectory(NULL);
    106113   
    107     fIntGhness = new TH1D("AccGammas",  "Acceptance", nbins, 0, 1);
    108     fIntPhness = new TH1D("AccHadrons", "Acceptance", nbins, 0, 1);
     114    fIntGhness = new TH1D("AccGammas",  "Integral Acceptance vs. Hadronness (Gammas)", nbins, 0, 1);
     115    fIntPhness = new TH1D("AccHadrons", "Integral Acceptance vs. Hadronness (Hadrons)", nbins, 0, 1);
    109116    fIntGhness->SetXTitle("Hadronness");
    110117    fIntPhness->SetXTitle("Hadronness");
    … …  
    145152Bool_t MHHadronness::SetupFill(const MParList *plist)
    146153{
    147     fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt");
    148     if (!fMcEvt)
    149     {
    150         *fLog << err << dbginf << "MMcEvt not found... aborting." << endl;
    151         return kFALSE;
     154    if (!fMatrix)
     155    {
     156        fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt");
     157        if (!fMcEvt)
     158        {
     159            *fLog << err << dbginf << "MMcEvt not found... aborting." << endl;
     160            return kFALSE;
     161        }
    152162    }
    153163
    154164    fHadronness = (MHadronness*)plist->FindObject("MHadronness");
     165
     166    fGhness->Reset();
     167    fPhness->Reset();
    155168
    156169    /*
    … …  
    200213        return kCONTINUE;
    201214
    202     if (fMcEvt->GetPartId()==kGAMMA)
     215    const Int_t particleid = fMatrix ? (Int_t)(*fMatrix)[fMap] : fMcEvt->GetPartId();
     216
     217    if (particleid==kGAMMA)
    203218        fGhness->Fill(h, w);
    204219    else
    … …  
    208223}
    209224
     225// --------------------------------------------------------------------------
     226//
     227// Returns the quality factor at gamma acceptance 0.5.
     228//
     229// If the histogram containes only two bins we return the
     230// naive quality: ag/sqrt(ah)
     231// with ag the acceptance of gammas and ah the acceptance of hadrons.
     232//
     233// You can use this (nbins=2) in case of some kind of filter cuts giving
     234// only a result: gamma yes/no (means discrete values of hadronness 0.25
     235// or 0.75)
     236//
     237// FIXME: In the later case weights cannot be used!
     238//
    210239Float_t MHHadronness::GetQ05() const
    211240{
    212     Int_t n = fQfac->GetN();
     241    if (fGhness->GetNbinsX()==2)
     242    {
     243        // acceptance of all gamma-like gammas  (h<0.5)
     244        const Double_t ig = fGhness->GetBinContent(1);
     245
     246        // acceptance of all gamma-like hadrons (h<0.5)
     247        const Double_t ip = fPhness->GetBinContent(1);
     248
     249        if (ip==0)
     250            return 0; // FIXME!
     251
     252        // naive quality factor
     253        const Double_t q = ig / sqrt(ip);
     254
     255        *fLog << all << ip << "/" << ig << ": " << q << endl;
     256
     257        return q;
     258    }
     259
     260    const Int_t n = fQfac->GetN();
    213261
    214262    Double_t val1x=0;
    … …  
    479527    }
    480528}
     529
     530// --------------------------------------------------------------------------
     531//
     532// You can use this function if you want to use a MHMatrix instead of
     533// MMcEvt. This function adds all necessary columns to the
     534// given matrix. Afterward you should fill the matrix with the corresponding
     535// data (eg from a file by using MHMatrix::Fill). If you now loop
     536// through the matrix (eg using MMatrixLoop) MHHadronness::Fill
     537// will take the values from the matrix instead of the containers.
     538//
     539void MHHadronness::InitMapping(MHMatrix *mat)
     540{
     541    if (fMatrix)
     542        return;
     543
     544    fMatrix = mat;
     545    fMap = fMatrix->AddColumn("MMcEvt.fPartId");
     546}
     547
     548void MHHadronness::StopMapping()
     549{
     550    fMatrix = NULL;
     551}
     552

  • trunk/MagicSoft/Mars/mhist/MHHadronness.h

    r2043 r2225  
    1111class MMcEvt;
    1212class MHadronness;
     13class MHMatrix;
    1314
    1415class MHHadronness : public MH
    … …  
    1718    const MMcEvt *fMcEvt;           //!
    1819    const MHadronness *fHadronness; //!
     20    MHMatrix *fMatrix;        //!
     21    Int_t fMap;                     //!
    1922
    2023    TH1D* fPhness;    //-> Hadrons Hadronness
    … …  
    4851    Bool_t Finalize();
    4952
     53    void InitMapping(MHMatrix *mat);
     54    void StopMapping();
     55
    5056    void Print(Option_t *option="") const;
    5157
Note: See TracChangeset for help on using the changeset viewer.