Changeset 6988


Ignore:
Timestamp:
05/02/05 10:17:41 (20 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
1 added
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r6987 r6988  
    2121
    2222                                                 -*-*- END OF LINE -*-*-
     23
     24 2005/05/02 Thomas Bretz
     25
     26   * ganymed.rc:
     27     - updated to the latest changes
     28
     29   * mhflux/MHAlpha.[h,cc]:
     30     - made a base class with the necessary interface to derived classes
     31       to support other variables like dca or theta
     32
     33   * mjobs/MJCut.[h,cc]:
     34     - added support for MHAlpha derived classes
     35
     36
    2337
    2438 2005/04/29 Thomas Bretz

  • trunk/MagicSoft/Mars/NEWS

    r6985 r6988  
    4646   - MEventRateCalc handles the calculation of the event rate more
    4747     accurate now in case of the start of a new run inside a sequence
     48
     49   - ganymed got support for using other variables than Alpha, eg. Theta.
     50     Therefor you need a class deriving from MHAlpha which supports
     51     this variable (one is already existing: MHTheta). It is setup
     52     through ganymed.rc
    4853
    4954

  • trunk/MagicSoft/Mars/ganymed.rc

    r6356 r6988  
    4747#MAlphaFitter.ScaleMax:             80
    4848#MAlphaFitter.PolynomOrder:         2
    49 #/*MISSING*/ ScaleMode
     49#MAlphaFitter.ScaleMode:            kSignificance
     50
     51# -------------------------------------------------------------------------
     52# Define here which histogram class to use to determin the signal.
     53# Currently availyble: MHAlpha, MHTheta <default>
     54# -------------------------------------------------------------------------
     55#MJCut.NameHist: MHAlpha
     56
    5057
    5158# -------------------------------------------------------------------------
    … …  
    6673# Get rid of triangular events...
    6774Cut0.Condition: {0}
    68 Cut0.0:  0.33 * log10(MHillas.fSize) < 1.05 - log10(MNewImagePar.fConc)
     75Cut0.0: log10(MNewImagePar.fConc1) < (-0.467)*log10(MHillas.fSize) +0.75
    6976
    7077# If you could setup MFEventSelector by resource file you could use it here
    … …  
    9299# Setup cuts in false source plot
    93100# -------------------------------------------------------------------------
    94 MHFalseSource.DistMin: 0.55
     101#MHFalseSource.DistMin: 0.55
    95102#MHFalseSource.DistMax: 0.55
    96103#MHFalseSource.DWMin: 0.55
    97104#MHFalseSource.DWMax: 0.55
    98 
    99 # -------------------------------------------------------------------------
    100 # Energy Estimation
    101 # -------------------------------------------------------------------------
    102 
    103 #EstimateEnergy: MEnergyEstimate
    104 EstimateEnergy.Rule: {0} + {1} + {2} + {3} + {4} + {5}
    105 EstimateEnergy.0:  6.3
    106 EstimateEnergy.1:  0.04*MHillasSrc.fDist*MGeomCam.fConvMm2Deg
    107 EstimateEnergy.2: -0.13*MHillas.fLength*MGeomCam.fConvMm2Deg
    108 EstimateEnergy.3:  0.15*MHillas.fSize
    109 EstimateEnergy.4:  0.0000519*MHillas.fSize*MHillasSrc.fDist*MGeomCam.fConvMm2Deg
    110 EstimateEnergy.5:  0.0000519*MHillas.fLength*MGeomCam.fConvMm2Deg*MHillasSrc.fDist*MGeomCam.fConvMm2Deg
    111 
    112105
    113106# -------------------------------------------------------------------------

  • trunk/MagicSoft/Mars/mhflux/MHAlpha.cc

    r6977 r6988  
    1818!   Author(s): Thomas Bretz, 3/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
    1919!
    20 !   Copyright: MAGIC Software Development, 2000-2004
     20!   Copyright: MAGIC Software Development, 2000-2005
    2121!
    2222!
    … …  
    8181//
    8282MHAlpha::MHAlpha(const char *name, const char *title)
    83     : fOffData(0), fResult(0), /*fExcess(0),*/ fEnergy(0), fHillas(0),
     83    : fNameParameter("MHillasSrc"), fParameter(0),
     84    fOffData(0), fResult(0), /*fExcess(0),*/ fEnergy(0), fHillas(0),
    8485    fPointPos(0), fTimeEffOn(0), fTime(0),
    8586    fSkipHistTime(kFALSE), fSkipHistTheta(kFALSE), fSkipHistEnergy(kFALSE),
    … …  
    9394    fTitle = title ? title : "Alpha plot";
    9495
    95     fHAlpha.SetName("Alpha");
    96     fHAlpha.SetTitle("Alpha");
    97     fHAlpha.SetXTitle("\\Theta [deg]");
    98     //fHAlpha.SetYTitle("E_{est} [GeV]");
    99     fHAlpha.SetZTitle("|\\alpha| [\\circ]");
    100     fHAlpha.SetDirectory(NULL);
    101     fHAlpha.UseCurrentStyle();
     96    fHist.SetName("Alpha");
     97    fHist.SetTitle("Alpha");
     98    fHist.SetXTitle("\\Theta [deg]");
     99    //fHist.SetYTitle("E_{est} [GeV]");
     100    fHist.SetZTitle("|\\alpha| [\\circ]");
     101    fHist.SetDirectory(NULL);
     102    fHist.UseCurrentStyle();
    102103
    103104    // Main histogram
    104     fHAlphaTime.SetName("Alpha");
    105     fHAlphaTime.SetXTitle("|\\alpha| [\\circ]");
    106     fHAlphaTime.SetYTitle("Counts");
    107     fHAlphaTime.UseCurrentStyle();
    108     fHAlphaTime.SetDirectory(NULL);
     105    fHistTime.SetName("Alpha");
     106    fHistTime.SetXTitle("|\\alpha| [\\circ]");
     107    fHistTime.SetYTitle("Counts");
     108    fHistTime.UseCurrentStyle();
     109    fHistTime.SetDirectory(NULL);
    109110
    110111
    … …  
    142143    binse.Apply(fHEnergy);
    143144    binst.Apply(fHTheta);
    144     binsa.Apply(fHAlphaTime);
    145 
    146     MH::SetBinning(&fHAlpha, &binst, &binse, &binsa);
    147 }
    148 
    149 void MHAlpha::FitEnergySpec(Bool_t paint)
    150 {
    151     /*
    152     TF1 f("Spectrum", "pow(x, [0])*exp(-x/[1])*[2]*x");//
    153     f.SetParameter(0,  -2.1);
    154     f.SetParameter(1,  1400); // 50
    155     f.SetParameter(2,  fHEnergy.Integral()); // 50
    156     f.SetLineWidth(2);
    157 
    158     fHEnergy.Fit(&f,  "0QI", "", 100, 2400); // Integral? 90, 2800
    159 
    160     const Double_t chi2 = f.GetChisquare();
    161     const Int_t    NDF  = f.GetNDF();
    162 
    163     // was fit successful ?
    164     const Bool_t ok = NDF>0 && chi2<2.5*NDF;
    165     f.SetLineColor(ok ? kGreen : kRed);
    166 
    167     if (paint)
    168     {
    169         f.Paint("same");
    170 
    171         if (ok)
    172         {
    173         TString s;
    174         s += Form("\\alpha=%.1f ",   f.GetParameter(0));
    175         s += Form("E_{c}=%.1f TeV ", f.GetParameter(1)/1000);
    176         s += Form("N=%.1e",          f.GetParameter(2));
    177 
    178         TLatex text(0.25, 0.94, s.Data());
    179         text.SetBit(TLatex::kTextNDC);
    180         text.SetTextSize(0.04);
    181         text.Paint();
    182         }
    183     }*/
     145    binsa.Apply(fHistTime);
     146
     147    MH::SetBinning(&fHist, &binst, &binse, &binsa);
     148}
     149
     150Double_t MHAlpha::GetVal() const
     151{
     152    return static_cast<const MHillasSrc*>(fParameter)->GetAlpha();
    184153}
    185154
    … …  
    189158    MAlphaFitter fit(fFit);
    190159
    191     const Int_t n = fHAlpha.GetNbinsY();
     160    const Int_t n = fHist.GetNbinsY();
    192161
    193162    fHEnergy.SetEntries(0);
    … …  
    197166    for (int i=1; i<=n; i++)
    198167    {
    199         if (fit.FitEnergy(fHAlpha, fOffData, i))
     168        if (fit.FitEnergy(fHist, fOffData, i))
    200169        {
    201170            fHEnergy.SetBinContent(i, fit.GetEventsExcess());
    … …  
    220189    MAlphaFitter fit(fFit);
    221190
    222     const Int_t n = fHAlpha.GetNbinsX();
     191    const Int_t n = fHist.GetNbinsX();
    223192
    224193    fHTheta.SetEntries(0);
    … …  
    226195    for (int i=1; i<=n; i++)
    227196    {
    228         if (fit.FitTheta(fHAlpha, fOffData, i))
     197        if (fit.FitTheta(fHist, fOffData, i))
    229198        {
    230199            fHTheta.SetBinContent(i, fit.GetEventsExcess());
    … …  
    237206}
    238207
     208Bool_t MHAlpha::GetParameter(const MParList &pl)
     209{
     210    fParameter = (MParContainer*)pl.FindObject(fNameParameter, "MHillasSrc");
     211    if (fParameter)
     212        return kTRUE;
     213
     214    *fLog << err << fNameParameter << " [MHillasSrc] not found... abort." << endl;
     215    return kFALSE;
     216}
     217
    239218Bool_t MHAlpha::SetupFill(const MParList *pl)
    240219{
    241     fHAlpha.Reset();
     220    fHist.Reset();
    242221    fHEnergy.Reset();
    243222    fHTheta.Reset();
    244223    fHTime.Reset();
    245224
    246     if (fName!=(TString)"MHAlphaOff" && fOffData==NULL)
    247     {
    248         MHAlpha *hoff = (MHAlpha*)pl->FindObject("MHAlphaOff", "MHAlpha");
     225    const TString off(Form("%sOff", ClassName()));
     226    if (fName!=off && fOffData==NULL)
     227    {
     228        const TString desc(Form("%s [%s] found... using ", off.Data(), ClassName()));
     229        MHAlpha *hoff = (MHAlpha*)pl->FindObject(off, ClassName());
    249230        if (!hoff)
    250             *fLog << inf << "No MHAlphaOff [MHAlpha] found... using current data only!" << endl;
     231            *fLog << inf << "No " << desc << "current data only!" << endl;
    251232        else
    252233        {
    253             *fLog << inf << "MHAlphaOff [MHAlpha] found... using on-off mode!" << endl;
     234            *fLog << inf << desc << "on-off mode!" << endl;
    254235            SetOffData(*hoff);
    255236        }
    256237    }
     238
     239    if (!GetParameter(*pl))
     240        return kFALSE;
    257241
    258242    fHillas = 0;
    … …  
    271255        fHEnergy.SetTitle(" N_{exc} vs. Size ");
    272256        fHEnergy.SetXTitle("Size [phe]");
    273         fHAlpha.SetYTitle("Size [phe]");
     257        fHist.SetYTitle("Size [phe]");
    274258    }
    275259    else
    … …  
    277261        fHEnergy.SetTitle(" N_{exc} vs. E_{est} ");
    278262        fHEnergy.SetXTitle("E_{est} [GeV]");
    279         fHAlpha.SetYTitle("E_{est} [GeV]");
     263        fHist.SetYTitle("E_{est} [GeV]");
    280264    }
    281265
    … …  
    306290
    307291    MBinning binst, binse, binsa;
    308     binst.SetEdges(fOffData ? *fOffData : fHAlpha, 'x');
    309     binse.SetEdges(fOffData ? *fOffData : fHAlpha, 'y');
    310     binsa.SetEdges(fOffData ? *fOffData : fHAlpha, 'z');
     292    binst.SetEdges(fOffData ? *fOffData : fHist, 'x');
     293    binse.SetEdges(fOffData ? *fOffData : fHist, 'y');
     294    binsa.SetEdges(fOffData ? *fOffData : fHist, 'z');
    311295    if (!fOffData)
    312296    {
    … …  
    324308                binse.SetEdges(*pl, "BinningSize");
    325309
    326         binsa.SetEdges(*pl, "BinningAlpha");
     310        binsa.SetEdges(*pl, Form("Binning%s", ClassName()+2));
    327311    }
    328312    else
    … …  
    330314        fHEnergy.SetTitle(fOffData->GetTitle());
    331315        fHEnergy.SetXTitle(fOffData->GetYaxis()->GetTitle());
    332         fHAlpha.SetYTitle(fOffData->GetYaxis()->GetTitle());
     316        fHist.SetYTitle(fOffData->GetYaxis()->GetTitle());
    333317    }
    334318
    335319    binse.Apply(fHEnergy);
    336320    binst.Apply(fHTheta);
    337     binsa.Apply(fHAlphaTime);
    338     MH::SetBinning(&fHAlpha, &binst, &binse, &binsa);
     321    binsa.Apply(fHistTime);
     322    MH::SetBinning(&fHist, &binst, &binse, &binsa);
    339323
    340324    MAlphaFitter *fit = (MAlphaFitter*)pl->FindObject("MAlphaFitter");
    … …  
    402386
    403387    TH1D *h = fOffData ? fOffData->ProjectionZ("ProjTimeTemp", -1, 9999, -1, 9999, "E") : 0;
    404     const Bool_t rc = fit.ScaleAndFit(fHAlphaTime, h);
     388    const Bool_t rc = fit.ScaleAndFit(fHistTime, h);
    405389    if (h)
    406390        delete h;
    … …  
    410394
    411395    // Reset Histogram
    412     fHAlphaTime.Reset();
    413     fHAlphaTime.SetEntries(0);
     396    fHistTime.Reset();
     397    fHistTime.SetEntries(0);
    414398
    415399    //
    … …  
    469453    else
    470454    {
    471         const MHillasSrc *hil = dynamic_cast<const MHillasSrc*>(par);
    472         if (!par)
    473         {
    474             *fLog << err << dbginf << "MHillasSrc not found... abort." << endl;
    475             return kFALSE;
    476         }
    477 
    478         alpha  = hil->GetAlpha();
     455        alpha  = GetVal();
    479456
    480457        if (fHillas)
    … …  
    491468
    492469    // Fill histograms
    493     fHAlpha.Fill(theta, energy, TMath::Abs(alpha), w);
     470    fHist.Fill(theta, energy, TMath::Abs(alpha), w);
    494471
    495472    // Check cuts
    … …  
    503480
    504481    if (!fSkipHistTime)
    505         fHAlphaTime.Fill(TMath::Abs(alpha), w);
     482        fHistTime.Fill(TMath::Abs(alpha), w);
    506483
    507484    return kTRUE;
    … …  
    557534
    558535        padsave->cd(1);
    559         TH1D *hon = fHAlpha.ProjectionZ("ProjAlpha");
     536        TH1D *hon = fHist.ProjectionZ("Proj");
    560537        if (fOffData)
    561538        {
    562             TH1D *hoff = fOffData->ProjectionZ("ProjAlphaOff");
     539            TH1D *hoff = fOffData->ProjectionZ("ProjOff");
    563540            const Double_t alpha = fFit.Scale(*hoff, *hon);
    564541
    … …  
    569546            hoff->SetMaximum(alpha);
    570547
    571             if ((h0=(TH1D*)gPad->FindObject("ProjAlphaOnOff")))
     548            if ((h0=(TH1D*)gPad->FindObject("ProjOnOff")))
    572549            {
    573550                h0->Reset();
    … …  
    583560    }
    584561
    585     if (o==(TString)"alpha")
    586         if ((h0 = (TH1D*)gPad->FindObject("ProjAlpha")))
     562    if (o==(TString)"variable")
     563        if ((h0 = (TH1D*)gPad->FindObject("Proj")))
    587564        {
    588565            // Check whether we are dealing with on-off analysis
    589             TH1D *hoff = (TH1D*)gPad->FindObject("ProjAlphaOff");
     566            TH1D *hoff = (TH1D*)gPad->FindObject("ProjOff");
    590567
    591568            // BE CARFEULL: This is a really weird workaround!
    … …  
    603580    if (o==(TString)"theta")
    604581    {
    605         TH1 *h = (TH1*)gPad->FindObject("Alpha_x");
     582        TH1 *h = (TH1*)gPad->FindObject(Form("%s_x", fHist.GetName()));
    606583        if (h)
    607584        {
    608             TH1D *h2 = (TH1D*)fHAlpha.Project3D("dum_x");
     585            TH1D *h2 = (TH1D*)fHist.Project3D("dum_x");
    609586            h2->SetDirectory(0);
    610587            h2->Scale(fHTheta.Integral()/h2->Integral());
    … …  
    618595    if (o==(TString)"energy")
    619596    {
    620         TH1 *h = (TH1*)gPad->FindObject("Alpha_y");
     597        TH1 *h = (TH1*)gPad->FindObject(Form("%s_y", fHist.GetName()));
    621598        if (h)
    622599        {
    623             TH1D *h2 = (TH1D*)fHAlpha.Project3D("dum_y");
     600            TH1D *h2 = (TH1D*)fHist.Project3D("dum_y");
    624601            h2->SetDirectory(0);
    625602            h2->Scale(fHEnergy.Integral()/h2->Integral());
    … …  
    634611            gPad->SetLogy();
    635612        }
    636         FitEnergySpec(kTRUE);
    637613    }
    638614
    … …  
    660636    gPad->SetBorderMode(0);
    661637
    662     h = fHAlpha.ProjectionZ("ProjAlpha", -1, 9999, -1, 9999, "E");
     638    h = fHist.ProjectionZ("Proj", -1, 9999, -1, 9999, "E");
    663639    h->SetBit(TH1::kNoTitle);
    664640    h->SetStats(kTRUE);
    665     h->SetXTitle("\\alpha [\\circ]");
     641    h->SetXTitle(fHist.GetZaxis()->GetTitle());
    666642    h->SetYTitle("Counts");
    667643    h->SetDirectory(NULL);
    … …  
    674650        h->SetMarkerColor(kGreen);
    675651
    676         h = fOffData->ProjectionZ("ProjAlphaOff", -1, 9999, -1, 9999, "E");
     652        h = fOffData->ProjectionZ("ProjOff", -1, 9999, -1, 9999, "E");
    677653        h->SetBit(TH1::kNoTitle);
    678         h->SetXTitle("\\alpha [\\circ]");
     654        h->SetXTitle(fHist.GetZaxis()->GetTitle());
    679655        h->SetYTitle("Counts");
    680656        h->SetDirectory(NULL);
    … …  
    684660        h->Draw("same");
    685661
    686         h = (TH1D*)h->Clone("ProjAlphaOnOff");
     662        h = (TH1D*)h->Clone("ProjOnOff");
    687663        h->SetBit(TH1::kNoTitle);
    688         h->SetXTitle("\\alpha [\\circ]");
     664        h->SetXTitle(fHist.GetZaxis()->GetTitle());
    689665        h->SetYTitle("Counts");
    690666        h->SetDirectory(NULL);
    … …  
    696672        TLine lin;
    697673        lin.SetLineStyle(kDashed);
    698         lin.DrawLine(0, 0, 90, 0);
     674        lin.DrawLine(h->GetXaxis()->GetXmin(), 0, h->GetXaxis()->GetXmax(), 0);
    699675    }
    700676
    701677    // After the Alpha-projection has been drawn. Fit the histogram
    702678    // and paint the result into this pad
    703     AppendPad("alpha");
     679    AppendPad("variable");
    704680
    705681    if (fHEnergy.GetNbinsX()>1)
    … …  
    713689        AppendPad("energy");
    714690
    715         h = (TH1D*)fHAlpha.Project3D("y");
     691        h = (TH1D*)fHist.Project3D("y");
    716692        h->SetBit(TH1::kNoTitle|TH1::kNoStats);
    717693        h->SetXTitle("E [GeV]");
    … …  
    749725        AppendPad("theta");
    750726
    751         h = (TH1D*)fHAlpha.Project3D("x");
     727        h = (TH1D*)fHist.Project3D("x");
    752728        h->SetBit(TH1::kNoTitle|TH1::kNoStats);
    753729        h->SetXTitle("\\theta [\\circ]");
    … …  
    768744    // FIXME: Do in Paint if special option given!
    769745    TCanvas *c = new TCanvas;
    770     Int_t n = fHAlpha.GetNbinsY();
     746    Int_t n = fHist.GetNbinsY();
    771747    Int_t nc = (Int_t)(TMath::Sqrt((Float_t)n-1)+1);
    772748    c->Divide(nc, nc, 1e-10, 1e-10);
    … …  
    775751    MAlphaFitter fit(fFit);
    776752
    777     for (int i=1; i<=fHAlpha.GetNbinsY(); i++)
     753    for (int i=1; i<=fHist.GetNbinsY(); i++)
    778754    {
    779755        c->cd(i);
    780756
    781         TH1D *hon = fHAlpha.ProjectionZ("ProjAlpha", -1, 9999, i, i, "E");
     757        TH1D *hon = fHist.ProjectionZ("Proj", -1, 9999, i, i, "E");
    782758        hon->SetBit(TH1::kNoTitle);
    783759        hon->SetStats(kTRUE);
    784         hon->SetXTitle("\\alpha [\\circ]");
     760        hon->SetXTitle(fHist.GetZaxis()->GetTitle());
    785761        hon->SetYTitle("Counts");
    786762        hon->SetDirectory(NULL);
    … …  
    796772            hon->SetMarkerColor(kGreen);
    797773
    798             hof = fOffData->ProjectionZ("ProjAlphaOff", -1, 9999, i, i, "E");
     774            hof = fOffData->ProjectionZ("ProjOff", -1, 9999, i, i, "E");
    799775            hof->SetBit(TH1::kNoTitle|TH1::kNoStats);
    800             hof->SetXTitle("\\alpha [\\circ]");
     776            hof->SetXTitle(fHist.GetZaxis()->GetTitle());
    801777            hof->SetYTitle("Counts");
    802778            hof->SetDirectory(NULL);
    … …  
    814790            diff->Add(hof, -1);
    815791            diff->SetBit(TH1::kNoTitle);
    816             diff->SetXTitle("\\alpha [\\circ]");
     792            diff->SetXTitle(fHist.GetZaxis()->GetTitle());
    817793            diff->SetYTitle("Counts");
    818794            diff->SetDirectory(NULL);
    … …  
    824800            TLine lin;
    825801            lin.SetLineStyle(kDashed);
    826             lin.DrawLine(0, 0, 90, 0);
     802            lin.DrawLine(diff->GetXaxis()->GetXmin(), 0, diff->GetXaxis()->GetXmax(), 0);
    827803
    828804            const Float_t min = diff->GetMinimum()*1.05;
    … …  
    833809            *fLog << dbg << "Bin " << i << ": sigma=" << fit.GetSignificance() << " omega=" << fit.GetGausSigma() << " events=" << fit.GetEventsExcess() << endl;
    834810        /*
    835         if (fit.FitEnergy(fHAlpha, fOffData, i, kTRUE))
     811        if (fit.FitEnergy(fHist, fOffData, i, kTRUE))
    836812        {
    837813            fHEnergy.SetBinContent(i, fit.GetEventsExcess());
    … …  
    845821Bool_t MHAlpha::Finalize()
    846822{
    847     //TH1D *h = fHAlpha.ProjectionZ("AlphaExc_px", -1, 9999, -1, 9999, "E");
     823    //TH1D *h = fHist.ProjectionZ("AlphaExc_px", -1, 9999, -1, 9999, "E");
    848824    //h->SetDirectory(0);
    849825    //Bool_t rc = fFit.Fit(*h);
    850826    //delete h;
    851827
    852     if (!fFit.FitAlpha(fHAlpha, fOffData))
     828    if (!fFit.FitAlpha(fHist, fOffData))
    853829    {
    854830        *fLog << warn << "MAlphaFitter - Fit failed..." << endl;
    … …  
    900876    fMatrix = mat;
    901877
    902     fMap[0] = fMatrix->AddColumn("MHillasSrc.fAlpha");
     878    fMap[0] = fMatrix->AddColumn(GetParameterRule());
    903879    fMap[1] = -1;
    904880    fMap[2] = -1;

  • trunk/MagicSoft/Mars/mhflux/MHAlpha.h

    r6977 r6988  
    2424class MPointingPos;
    2525
    26 /*
    27 class MAlphaFitter : public MParContainer
    28 {
    29 private:
    30     Float_t fSigInt;
    31     Float_t fSigMax;
    32     Float_t fBgMin;
    33     Float_t fBgMax;
    34     Int_t   fPolynom;
    35 
    36     Double_t fSignificance;
    37     Double_t fExcessEvents;
    38     Double_t fChiSqSignal;
    39     Double_t fChiSqBg;
    40     Double_t fSigmaGaus;
    41     Double_t fIntegralMax;
    42 
    43 public:
    44     MAlphaFitter() : fSigInt(10), fSigMax(75), fBgMin(45), fBgMax(85), fPolynom(2)
    45     {
    46     }
    47 
    48     MAlphaFitter(const MAlphaFitter &f)
    49     {
    50         f.Copy(*this);
    51     }
    52 
    53     void Copy(TObject &o) const
    54     {
    55         MAlphaFitter &f = static_cast<MAlphaFitter&>(o);
    56 
    57         f.fSigInt  = fSigInt;
    58         f.fSigMax  = fSigMax;
    59         f.fBgMin   = fBgMin;
    60         f.fBgMax   = fBgMax;
    61         f.fPolynom = fPolynom;
    62     }
    63 
    64     void SetSignalIntegralMax(Float_t s) { fSigInt = s; }
    65     void SetSignalFitMax(Float_t s)      { fSigMax = s; }
    66     void SetBackgroundFitMin(Float_t s)  { fBgMin = s; }
    67     void SetBackgroundFitMax(Float_t s)  { fBgMax = s; }
    68     void SetNumPolynom(Int_t s)          { fPolynom = s; }
    69 
    70     Double_t GetExcessEvents() const { return fExcessEvents; }
    71     Double_t GetSignificance() const { return fSignificance; }
    72     Double_t GetChiSqSignal() const  { return fChiSqSignal; }
    73     Double_t GetChiSqBg() const      { return fChiSqBg; }
    74     Double_t GetSigmaGaus() const    { return fSigmaGaus; }
    75 
    76     void PaintResult(Float_t x=0.04, Float_t y=0.94, Float_t size=0.035) const;
    77 
    78     Bool_t Fit(TH1D &h, Bool_t paint=kFALSE);
    79     ClassDef(MAlphaFitter, 1)
    80 };
    81 */
    82 
    8326class MHAlpha : public MH
    8427{
     28protected:
     29    TH3D    fHist;              // Alpha vs. theta and energy
     30    TH1D    fHistTime;          //! temporary histogram to get alpha vs. time
     31
     32    TString fNameParameter;
     33
     34    MParContainer *fParameter;  //!
     35
    8536private:
    8637    const TH3D *fOffData;
    … …  
    8839    MAlphaFitter fFit;          // SEEMS THAT STREAMER HAS SOME PROBLEMS... MAYBE IF FUNC IS USED AT THE SAME TIME FOR FITS (PAINT)
    8940
    90     TH3D fHAlpha;               // Alpha vs. theta and energy
    9141    TH1D fHEnergy;              // excess events vs energy
    9242    TH1D fHTheta;               // excess events vs theta
    9343    TH1D fHTime;                // excess events vs time;
    94     TH1D fHAlphaTime;           //! temporary histogram to get alpha vs. time
    9544
    96     MParameterD  *fResult;      //!
    97     MParameterD  *fEnergy;      //!
    98     MHillas      *fHillas;      //!
    99     MPointingPos *fPointPos;    //!
     45    MParameterD   *fResult;     //!
     46    MParameterD   *fEnergy;     //!
     47    MHillas       *fHillas;     //!
     48    MPointingPos  *fPointPos;   //!
    10049
    10150    MTime   *fTimeEffOn;        //! Time to steer filling of fHTime
    … …  
    11665    Int_t fMap[5];              //!
    11766
    118     void FitEnergySpec(Bool_t paint=kFALSE);
    11967    Float_t FitEnergyBins(Bool_t paint=kFALSE);
    12068    void FitThetaBins(Bool_t paint=kFALSE);
    … …  
    12876    Int_t DistancetoPrimitive(Int_t px, Int_t py) { return 0; }
    12977
     78    virtual Bool_t      GetParameter(const MParList &pl);
     79    virtual Double_t    GetVal() const;
     80    virtual const char *GetParameterRule() const
     81    {
     82        return "MHillasSrc.fAlpha";
     83    }
     84
    13085public:
    13186    MHAlpha(const char *name=NULL, const char *title=NULL);
    13287
     88    // MH
    13389    Bool_t SetupFill(const MParList *pl);
    13490    Bool_t Fill(const MParContainer *par, const Stat_t w=1);
    13591    Bool_t Finalize();
    13692
    137     const TH3D &GetHist() const { return fHAlpha; }
     93    // MHAlpha
     94    const TH3D &GetHist() const { return fHist; }
    13895    const MAlphaFitter &GetAlphaFitter() const { return fFit; }
    13996
    14097    const TH1D &GetHEnergy() const { return fHEnergy; }
    14198
     99    // Setter
     100    void SetNameParameter(const char *name) { fNameParameter=name; }
    142101    void SetOffData(const MHAlpha &h)
    143102    {
    144         fOffData        = &h.fHAlpha;
     103        fOffData        = &h.fHist;
    145104        fForceUsingSize =  h.fForceUsingSize;
    146105        fNumTimeBins    =  h.fNumTimeBins;

  • trunk/MagicSoft/Mars/mjobs/MJCut.cc

    r6979 r6988  
    7777MJCut::MJCut(const char *name, const char *title)
    7878    : fStoreSummary(kFALSE), fStoreResult(kTRUE), fWriteOnly(kFALSE),
    79     fIsWobble(kFALSE), fIsMonteCarlo(kFALSE),  fFullDisplay(kFALSE), /*fSubstraction(kFALSE),*/
    80     /*fEstimateEnergy(0),*/ fCalcHadronness(0)
     79    fIsWobble(kFALSE), fIsMonteCarlo(kFALSE),  fFullDisplay(kFALSE),
     80    fNameHist("MHTheta"), fCalcHadronness(0)
    8181{
    8282    fName  = name  ? name  : "MJCut";
    … …  
    258258
    259259    // Save also the result, not only the setup
    260     const MHAlpha *halpha = (MHAlpha*)plist.FindObject("MHAlpha");
     260    const MHAlpha *halpha = (MHAlpha*)plist.FindObject(fNameHist, "MHAlpha");
    261261    if (halpha)
    262262        arr.Add((TObject*)(&halpha->GetAlphaFitter()));
    … …  
    285285//   MJCut.WriteResult:  yes, no
    286286//   MJCut.ResultFile:   filename
     287//   MJCut.HistName:     MHAlpha
    287288//
    288289Bool_t MJCut::CheckEnvLocal()
    … …  
    301302    EnableFullDisplay(GetEnv("FullDisplay", fFullDisplay));
    302303    //EnableSubstraction(GetEnv("HistogramSubstraction", fSubstraction));
     304
     305    SetNameHist(GetEnv("NameHist", fNameHist));
    303306
    304307    return kTRUE;
    … …  
    341344}
    342345
    343 void MJCut::DisplayResult(const MParList &plist) const
    344 {
    345     /*
    346      TObject *h1 = plist.FindObject("MHHillasOffPre",  "MHHillas");
    347      TObject *h2 = plist.FindObject("MHHillasOffPost", "MHHillas");
    348      TObject *h3 = plist.FindObject("MHVsSizeOffPost", "MHVsSize");
    349      TObject *h4 = plist.FindObject("MHHilExtOffPost", "MHHillasExt");
    350      TObject *h5 = plist.FindObject("MHHilSrcOffPost", "MHHillasSrc");
    351      TObject *h6 = plist.FindObject("MHImgParOffPost", "MHImagePar");
    352      TObject *h7 = plist.FindObject("MHNewParOffPost", "MHNewImagePar");
    353      */
     346// --------------------------------------------------------------------------
     347//
     348// Create a new instance of an object with name name of class
     349// type fNameHist in parlist. It must derive from MHAlpha.
     350// Call ForceUsingSize for it and return its pointer.
     351// If something fails NULL is returned.
     352//
     353MHAlpha *MJCut::CreateNewHist(MParList &plist, const char *name) const
     354{
     355    TClass *cls = gROOT->GetClass(fNameHist);
     356    if (!cls)
     357    {
     358        *fLog << err << "Class " << fNameHist << " not found in dictionary... abort." << endl;
     359        return NULL;
     360    }
     361    if (!cls->InheritsFrom(MHAlpha::Class()))
     362    {
     363        *fLog << err << "Class " << fNameHist << " doesn't inherit from MHAlpha... abort." << endl;
     364        return NULL;
     365    }
     366
     367    const TString objname(Form("%s%s", fNameHist.Data(), name));
     368    MHAlpha *h = (MHAlpha*)plist.FindCreateObj(fNameHist, objname);
     369    if (!h)
     370        return NULL;
     371
     372    h->ForceUsingSize();
     373
     374    return h;
    354375}
    355376
    … …  
    424445        fit.SetScaleMode(MAlphaFitter::kNone);
    425446
    426     MHAlpha halphaoff("MHAlphaOff");
    427     halphaoff.ForceUsingSize();
    428     plist.AddToList(&halphaoff);
    429 
    430     MFillH falpha("MHAlphaOff [MHAlpha]", "MHillasSrc", "FillAlpha");
     447    MHAlpha *halphaoff = CreateNewHist(plist, "Off");
     448    MFillH falpha(halphaoff, "MHillasSrc", "FillAlpha");
    431449    MFillH ffs("MHFalseSourceOff [MHFalseSource]", "MHillas", "FillFS");
    432450
    … …  
    655673    }
    656674    */
    657     MHAlpha halphaon("MHAlpha");
    658     halphaon.ForceUsingSize();
    659     plist.AddToList(&halphaon);
    660 
    661     MFillH falpha2("MHAlpha", "MHillasSrc", "FillAlpha");
     675    MHAlpha *halphaon=CreateNewHist(plist);
     676    MFillH falpha2(halphaon, "MHillasSrc", "FillAlpha");
    662677    MFillH ffs2("MHFalseSource", "MHillas", "FillFS");
    663678
    … …  
    706721    tlist.PrintStatistics();
    707722
    708     DisplayResult(plist);
    709 
    710723    // FIXME: Perform fit and plot energy dependant alpha plots
    711724    // and fit result to new tabs!

  • trunk/MagicSoft/Mars/mjobs/MJCut.h

    r6979 r6988  
    99class MDataSet;
    1010class MParList;
     11class MHAlpha;
    1112class MWriteRootFile;
    1213
    … …  
    2627    TString fNameOutput;
    2728
     29    TString fNameHist;
     30
    2831    //MTask *fEstimateEnergy;
    2932    MTask *fCalcHadronness;
    3033
    31     TString GetOutputFile(UInt_t num) const;
    32     Bool_t  CheckEnvLocal();
    33     void    SetupWriter(MWriteRootFile *write, const char *name) const;
    34     Bool_t  WriteTasks(UInt_t num, TObjArray &cont) const;
    35     Bool_t  WriteResult(const MParList &plist, UInt_t num) const;
    36     void    DisplayResult(const MParList &plist) const;
     34    TString  GetOutputFile(UInt_t num) const;
     35    Bool_t   CheckEnvLocal();
     36    void     SetupWriter(MWriteRootFile *write, const char *name) const;
     37    Bool_t   WriteTasks(UInt_t num, TObjArray &cont) const;
     38    Bool_t   WriteResult(const MParList &plist, UInt_t num) const;
     39    MHAlpha *CreateNewHist(MParList &plist, const char *name="") const;
    3740
    38     Bool_t  CanStoreSummary() const { return !fPathOut.IsNull() && fStoreSummary; }
    39     Bool_t  CanStoreResult() const  { return !fPathOut.IsNull() && fStoreResult;  }
     41    Bool_t   CanStoreSummary() const { return !fPathOut.IsNull() && fStoreSummary; }
     42    Bool_t   CanStoreResult() const  { return !fPathOut.IsNull() && fStoreResult;  }
    4043
    4144public:
    … …  
    5760    void SetNameOutFile(const char *name="")      { fNameOutput=name; }
    5861
     62    void SetNameHist(const char *name) { fNameHist=name; }
     63
    5964    //void SetEnergyEstimator(const MTask *task=0);
    6065    void SetHadronnessCalculator(const MTask *task=0);
Note: See TracChangeset for help on using the changeset viewer.