Changeset 7094


Ignore:
Timestamp:
05/27/05 09:46:22 (20 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft
Files:
2 added
6 deleted
16 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r7091 r7094  
    2929     - implement usage of extractor resolution (only if UseExtractorRes
    3030       is set).
     31
     32
     33 2005/05/27 Thomas Bretz
     34
     35   * Makefile:
     36     - removed mmontecarlo directory
     37
     38   * mmontecarlo/MMcEnergyEst.[h,cc],
     39     mmontecarlo/MMcTimeGenerate.[h,cc],
     40     mmontecarlo/MMcWeightEnergySpecCalc.[h,cc]:
     41     - removed
     42
     43   * sponde.rc:
     44     - added new line for weighted spectral index
     45
     46   * mbadpixels/MBadPixelsCalc.[h,cc]:
     47     - added an option to perform the checks also in PostProcess
     48
     49   * mhbase/MFillH.h:
     50     - added default argument to SetWeight
     51
     52   * mhbase/MH3.h:
     53     - added Sumw2() member function
     54
     55   * mhflux/MHCollectionArea.[h,cc]:
     56     - added TLatex output to plots
     57     - added some Getter
     58
     59   * mjobs/MJSpectrum.[h,cc]:
     60     - implemented the possibility to weight the monte carlo spectrum
     61       to a new index or function. More details can be found
     62       in MMcSpectrumWeight
     63     - slightly changed the plot comparing the size distributions
     64     - scale the comparsison plots by the resulting spectrum
     65
     66   * mjobs/MJob.[h,cc]:
     67     - added a member function to check ReadEnv of a single
     68       container
     69
     70   * mjobs/Makefile:
     71     - added -I../mmc
     72
     73   * mmc/MMcEvt.[hxx,cxx], mmc/McEvtBasic.[hxx,cxx]:
     74     - changed the inheritance: MMcEvt now derives from MMcEvtBasic
     75       so that both classes are interchangable
     76     - increased both class versions
     77     - chaged the default partictle in MMcEvtBasic from
     78       kGAMMA to kUNDEFINED
     79     - added new particle type: kUNDEFINED
     80
     81   * mhflux/MMcSpectrumWeight.[h,cc]:
     82     - added
    3183
    3284

  • trunk/MagicSoft/Mars/Makefile

    r6979 r7094  
    5656          msql \
    5757          mimage \
    58           mmontecarlo \
    5958          mhft \
    6059          mmc \

  • trunk/MagicSoft/Mars/NEWS

    r7082 r7094  
    22
    33 *** Version <cvs>
     4
     5   - general: MMcEvt now derived from MMcEvtBasic which should
     6     have no influence on compatibility with older camera files
     7
     8   - sponde: The input MC spectrum can now be weighted to fake a
     9     different spectrum. This is done via MMcSpectrumWeight. For
     10     more details see the class description and sponde.rc
     11
     12   - sponde: The paremetr comparsion plots are not scaled by
     13     their entries anymore. Instead the MC plot is scaled by using
     14     the result spectrum of the analysis. If the input MC spectrum
     15     and the result spectrum has different slopes the absolut
     16     normalization is normally wrong.
    417
    518

  • trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCalc.cc

    r7035 r7094  
    8686//
    8787MBadPixelsCalc::MBadPixelsCalc(const char *name, const char *title)
    88     : fPedestalLevel(3), fPedestalLevelVariance(5), fNamePedPhotCam("MPedPhotCam")
     88    : fPedestalLevel(3), fPedestalLevelVariance(5), fNamePedPhotCam("MPedPhotCam"),
     89    fCheckInProcess(kTRUE), fCheckInPostProcess(kFALSE)
    8990{
    9091    fName  = name  ? name  : gsDefName.Data();
    … …  
    129130// --------------------------------------------------------------------------
    130131//
    131 // Check the pedestal RMS of every pixel with respect to the mean pedestal RMS
    132 // of the camera (Sigmabar).
    133 //
    134 // The pixels can be set as blind if the pedestalRMS is too big or 0.
    135 //
    136 // If you don't want to use this option set the PedestalLevel<=0;
    137 //
    138 //     MBadPixelsCalc calc;
    139 //     calc.SetPedestalLevel(-1);
    140 /*
    141 void MBadPixelsCalc::CheckPedestalRMS() const
    142 {
    143     const Int_t entries = fPedPhotCam->GetSize();
    144 
    145     const Float_t meanPedRMS = fSigmabar->GetSigmabar();
    146 
    147     for (Int_t i=0; i<entries; i++)
    148     {
    149         //
    150         // get pixel as entry from list
    151         //
    152         const Double_t nratio    = fGeomCam->GetPixRatio(i);
    153         const Double_t pixPedRms = (*fPedPhotCam)[i].GetRms();
    154 
    155         if (pixPedRms*nratio > fPedestalLevel * meanPedRMS || pixPedRms == 0)
    156             (*fBadPixels)[i].SetUnsuitable(MBadPixelsPix::kUnsuitableEvt);
    157     }
    158 }
    159 */
    160 
    161 // --------------------------------------------------------------------------
    162 //
    163132// Check the pedestal Rms of the pixels: compute with 2 iterations the mean
    164133// for inner and outer pixels. Set as blind the pixels with too small or
    165134// too high pedestal Rms with respect to the mean.
    166135//
    167 Bool_t MBadPixelsCalc::CheckPedestalRms() const
    168 {
     136Bool_t MBadPixelsCalc::CheckPedestalRms(MBadPixelsPix::UnsuitableType_t type) const
     137{
     138    if (fPedestalLevel<=0 && fPedestalLevelVariance<=0)
     139        return kTRUE;
     140
    169141    const Int_t entries = fPedPhotCam->GetSize();
    170142
    … …  
    266238            continue;
    267239
    268         (*fBadPixels)[i].SetUnsuitable(MBadPixelsPix::kUnsuitableEvt);
     240        (*fBadPixels)[i].SetUnsuitable(type);
     241        if (type==MBadPixelsPix::kUnsuitableRun)
     242            (*fBadPixels)[i].SetUncalibrated(MBadPixelsPix::kDeadPedestalRms);
     243
    269244        bads++;
    270245    }
    … …  
    281256}
    282257
    283 
    284258// --------------------------------------------------------------------------
    285259//
    … …  
    287261Int_t MBadPixelsCalc::Process()
    288262{
    289     if (fPedestalLevel>0 || fPedestalLevelVariance>0)
    290         CheckPedestalRms();
    291 
    292     return kTRUE;
     263    return fCheckInProcess ? CheckPedestalRms(MBadPixelsPix::kUnsuitableEvt) : kTRUE;
     264}
     265
     266// --------------------------------------------------------------------------
     267//
     268//
     269Int_t MBadPixelsCalc::PostProcess()
     270{
     271    return fCheckInPostProcess ? CheckPedestalRms(MBadPixelsPix::kUnsuitableRun) : kTRUE;
    293272}
    294273

  • trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCalc.h

    r5777 r7094  
    44#ifndef MARS_MTask
    55#include "MTask.h"
     6#endif
     7#ifndef MARS_MBadPixelsPix
     8#include "MBadPixelsPix.h"
    69#endif
    710
    … …  
    2124
    2225    TString fNamePedPhotCam; // name of the 'MPedPhotCam' container
    23    
    24     //    void CheckPedestalRMS() const;
    25     Bool_t CheckPedestalRms() const;
    2626
     27    Bool_t fCheckInProcess;
     28    Bool_t fCheckInPostProcess;
     29
     30    // MBadPixelsCalc
     31    Bool_t CheckPedestalRms(MBadPixelsPix::UnsuitableType_t t) const;
     32
     33    // MTask
    2734    Int_t PreProcess(MParList *pList);
    2835    Int_t Process();
     36    Int_t PostProcess();
     37
     38    // MParContainer
    2939    Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print);
    3040
    … …  
    3242    MBadPixelsCalc(const char *name=NULL, const char *title=NULL);
    3343
     44    // Setter
    3445    void SetPedestalLevel(Float_t f)         { fPedestalLevel=f; }
    3546    void SetPedestalLevelVariance(Float_t f) { fPedestalLevelVariance=f; }
    3647    void SetNamePedPhotCam(const char *name) { fNamePedPhotCam = name; }
     48
     49    void EnableCheckInProcess(Bool_t b=kTRUE)     { fCheckInProcess = b; }
     50    void EnableCheckInPostProcess(Bool_t b=kTRUE) { fCheckInPostProcess = b; }
    3751
    3852    ClassDef(MBadPixelsCalc, 1) // Task to find bad pixels (star, broken pixels, etc)

  • trunk/MagicSoft/Mars/mhbase/MFillH.h

    r6915 r7094  
    6363
    6464    void SetWeight(MParameterD *w)   { fWeight = w; }
    65     void SetWeight(const char *name) { fWeightName = name; }
     65    void SetWeight(const char *name="MWeight") { fWeightName = name; }
    6666
    6767    void      SetDrawOption(Option_t *option="");

  • trunk/MagicSoft/Mars/mhbase/MH3.h

    r6978 r7094  
    2626    Double_t    fScale[3];       // Scale for the three axis (eg unit)
    2727
    28     // TString     fNameProfX;      //! This should make sure, that gROOT doen't confuse the profile with something else
    29     // TString     fNameProfY;      //! This should make sure, that gROOT doen't confuse the profile with something else
    30 
    3128    void StreamPrimitive(ofstream &out) const;
    3229
    … …  
    3633        kIsLogz = BIT(19)
    3734    };
    38 
    39 //    Int_t DistancetoPrimitive(Int_t px, Int_t py);
    40 //    void  ExecuteEvent(Int_t event, Int_t px, Int_t py);
    4135
    4236public:
    … …  
    4842    ~MH3();
    4943
     44    // Setter
    5045    void SetScaleX(Double_t scale) { fScale[0] = scale; }
    5146    void SetScaleY(Double_t scale) { fScale[1] = scale; }
    … …  
    5651    void SetLogz(Bool_t b=kTRUE) { b ? fHist->SetBit(kIsLogz) : fHist->ResetBit(kIsLogz); }
    5752
     53    void Sumw2() const { if (fHist) fHist->Sumw2(); }
     54
     55    // Getter
    5856    Int_t GetDimension() const { return fDimension; }
    5957    Int_t GetNbins() const;
    6058    Int_t FindFixBin(Double_t x, Double_t y=0, Double_t z=0) const;
    6159
    62     void SetName(const char *name);
    63     void SetTitle(const char *title);
    64 
    65     Bool_t SetupFill(const MParList *pList);
    66     Bool_t Fill(const MParContainer *par, const Stat_t w=1);
     60    TH1 &GetHist() { return *fHist; }
     61    const TH1 &GetHist() const { return *fHist; }
    6762
    6863    TString GetDataMember() const;
    6964    TString GetRule(const Char_t axis='x') const;
    7065
    71     TH1 &GetHist() { return *fHist; }
    72     const TH1 &GetHist() const { return *fHist; }
     66    // MH
     67    Bool_t SetupFill(const MParList *pList);
     68    Bool_t Fill(const MParContainer *par, const Stat_t w=1);
    7369
    7470    TH1 *GetHistByName(const TString name="") const { return fHist; }
    … …  
    7975    }
    8076
     77    // MParContainer
     78    MParContainer *New() const;
     79
     80    // TObject
     81    void SetName(const char *name);
     82    void SetTitle(const char *title);
     83
    8184    void SetColors() const;
    8285    void Draw(Option_t *opt=NULL);
    8386    void Paint(Option_t *opt="");
    84 
    85     MParContainer *New() const;
    8687
    8788    ClassDef(MH3, 1) // Generalized 1/2/3D-histogram for Mars variables

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

    r6976 r7094  
    3131#include "MHCollectionArea.h"
    3232
     33#include <TLatex.h>
    3334#include <TCanvas.h>
    3435#include <TPaveStats.h>
    … …  
    8182    fHistAll.SetTitle("Number of events produced");
    8283    fHistAll.SetXTitle("\\Theta [deg]");
    83     fHistAll.SetYTitle("E [GeV]");
     84    fHistAll.SetYTitle("E_{mc} [GeV]");
    8485    fHistAll.SetDirectory(NULL);
    8586    fHistAll.UseCurrentStyle();
    8687
    8788    fHEnergy.SetName("CollEnergy");
    88     fHEnergy.SetTitle("Collection Area vs Energy");
     89    fHEnergy.SetTitle("Collection Area");
    8990    fHEnergy.SetXTitle("E [GeV]");
    90     fHEnergy.SetYTitle("A [m^{2}]");
     91    fHEnergy.SetYTitle("A_{eff} [m^{2}]");
    9192    fHEnergy.SetDirectory(NULL);
    9293    fHEnergy.UseCurrentStyle();
    … …  
    119120    // read from run header, but it is not yet in!!
    120121    //
    121     const Float_t r1 = 0;
    122     const Float_t r2 = fMcAreaRadius;
    123 
    124     //*fLog << warn << endl << dbginf << "WARNING! I will assume a maximum impact parameter of " << fMCAreaRadius << " meters for the MC events. Check that this is the true one!" <<endl<<endl;
    125     const Float_t total_area = TMath::Pi() * (r2*r2 - r1*r1);
     122    const Float_t totalarea = GetCollectionAreaAbs();//TMath::Pi() * (r2*r2 - r1*r1);
    126123
    127124    for (Int_t ix=1; ix<=nbinx; ix++)
    … …  
    141138            continue;
    142139
    143         const Double_t eff = Ns/Na;
    144 
    145         const Double_t efferr = TMath::Sqrt((1.-eff)*Ns)/Na;
     140        const Double_t eff         = Ns/Na;
     141        const Double_t efferr      = TMath::Sqrt((1.-eff)*Ns)/Na;
    146142       
    147         const Float_t col_area  =  eff * total_area;
    148         const Float_t col_area_error =  efferr * total_area;
    149 
    150         fHEnergy.SetBinContent(ix, col_area);
    151         fHEnergy.SetBinError(ix, col_area_error);
     143        const Float_t colarea      =  eff    * totalarea;
     144        const Float_t colareaerror =  efferr * totalarea;
     145
     146        fHEnergy.SetBinContent(ix, colarea);
     147        fHEnergy.SetBinError(ix,   colareaerror);
    152148    }
    153149
    … …  
    275271void MHCollectionArea::Paint(Option_t *option)
    276272{
     273    if (TString(option)=="paint3")
     274    {
     275        /*
     276        TH1 *h = dynamic_cast<TH1*>(gPad->FindObject("Efficiency"));
     277        if (h)
     278        {
     279            const TString txt = Form("N/N_{0}=%.2f",
     280                                 GetCollectionAreaEff(),
     281                                 GetCollectionAreaAbs(), fMcAreaRadius);
     282
     283        TLatex text(0.31, 0.95, txt);
     284        text.SetBit(TLatex::kTextNDC);
     285        text.SetTextSize(0.04);
     286        text.Paint();*/
     287        return;
     288    }
     289    if (TString(option)=="paint4")
     290    {
     291        const TString txt = Form("A_{eff}=%.0fm^{2}  A_{abs}=%.0fm^{2}  r=%.0fm",
     292                                 GetCollectionAreaEff(),
     293                                 GetCollectionAreaAbs(), fMcAreaRadius);
     294
     295        TLatex text(0.31, 0.95, txt);
     296        text.SetBit(TLatex::kTextNDC);
     297        text.SetTextSize(0.04);
     298        text.Paint();
     299        return;
     300    }
     301
    277302    TVirtualPad *pad;
    278303
    … …  
    405430        h->SetNameTitle("Efficiency", "Combined cut and trigger efficiency");
    406431        h->SetDirectory(NULL);
     432        AppendPad("paint3");
    407433    }
    408434    else
    … …  
    416442        gPad->SetGridy();
    417443        fHEnergy.Draw();
     444        AppendPad("paint4");
    418445    }
    419446    else
    … …  
    428455    //*fLog << energy << " " << theta << endl;
    429456
    430     fHistSel.Fill(theta, energy);
     457    fHistSel.Fill(theta, energy, weight);
    431458
    432459    return kTRUE;

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

    r6976 r7094  
    5656    const TH1D &GetHEnergy() const { return fHEnergy; }
    5757
     58    Double_t GetEntries() const { return fHistAll.Integral(); }
     59    Double_t GetCollectionAreaEff() const { return fHEnergy.Integral(); }
     60    Double_t GetCollectionAreaAbs() const { return TMath::Pi()*fMcAreaRadius*fMcAreaRadius; }
     61
    5862/*
    5963    void DrawAll(Option_t *option="");

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

    r7071 r7094  
    5252#include "MBinning.h"
    5353#include "MDataSet.h"
     54#include "MMcCorsikaRunHeader.h"
    5455
    5556// Spectrum
    … …  
    5859#include "../mhflux/MHCollectionArea.h"
    5960#include "../mhflux/MHEnergyEst.h"
     61#include "../mhflux/MMcSpectrumWeight.h"
    6062
    6163// Eventloop
    … …  
    155157}
    156158
     159// --------------------------------------------------------------------------
     160//
     161// Read the first MMcCorsikaRunHeader from the RunHeaders tree in
     162// the dataset.
     163// The simulated energy range and spectral slope is initialized from
     164// there.
     165// In the following eventloops the forced check in MMcSpectrumWeight
     166// ensures, that the spectral slope and energy range doesn't change.
     167//
     168Bool_t  MJSpectrum::InitWeighting(const MDataSet &set, MMcSpectrumWeight &w) const
     169{
     170    fLog->Separator("Initialize energy weighting");
     171
     172    MParList l;
     173    l.AddToList(&w);
     174    if (!CheckEnv(l))
     175    {
     176        *fLog << err << "ERROR - Reading resources for MMcSpectrumWeight failed." << endl;
     177        return kFALSE;
     178    }
     179
     180    TChain chain("RunHeaders");
     181    set.AddFilesOn(chain);
     182
     183    MMcCorsikaRunHeader *h=0;
     184    chain.SetBranchAddress("MMcCorsikaRunHeader.", &h);
     185    chain.GetEntry(1);
     186
     187    if (!h)
     188    {
     189        *fLog << err << "ERROR - Couldn't read MMcCorsikaRunHeader from DataSet." << endl;
     190        return kFALSE;
     191    }
     192
     193    if (!w.Set(*h))
     194    {
     195        *fLog << err << "ERROR - Initializing MMcSpectrumWeight failed." << endl;
     196        return kFALSE;
     197    }
     198
     199    w.Print();
     200    return kTRUE;
     201}
     202
    157203Float_t MJSpectrum::ReadInput(MParList &plist, TH1D &h1, TH1D &h2)
    158204{
    … …  
    211257}
    212258
    213 Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h) const
     259Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h, MMcSpectrumWeight &weight) const
    214260{
    215261    // Some debug output
    216262    fLog->Separator("Compiling original MC distribution");
    217263
    218     *fLog << inf << "Please stand by, this may take a while..." << flush;
     264    weight.SetNameMcEvt("MMcEvtBasic");
     265    const TString w(weight.GetFormulaWeights());
     266    weight.SetNameMcEvt();
     267
     268    *fLog << inf << "Using weights: " << w << endl;
     269    *fLog << "Please stand by, this may take a while..." << flush;
    219270
    220271    if (fDisplay)
    … …  
    236287        h.SetYTitle("E [GeV]");
    237288        h.SetZTitle("Counts");
    238         chain.Draw("MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC", "", "goff");
     289        chain.Draw("MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC", w, "goff");
    239290    }
    240291    else
    … …  
    243294        h.SetXTitle("\\Theta [\\circ]");
    244295        h.SetYTitle("Counts");
    245         chain.Draw("MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC", "", "goff");
     296        chain.Draw("MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC", w, "goff");
    246297    }
    247298    h.SetDirectory(0);
    … …  
    413464}
    414465
    415 Bool_t MJSpectrum::IntermediateLoop(MParList &plist, MH3 &mh1, TH1D &temp1, const MDataSet &set) const
     466Bool_t MJSpectrum::IntermediateLoop(MParList &plist, MH3 &mh1, TH1D &temp1, const MDataSet &set, MMcSpectrumWeight &weight) const
    416467{
    417468    MTaskList tlist1;
    … …  
    439490    }
    440491    tlist1.AddToList(&readmc);
     492    tlist1.AddToList(&weight);
    441493
    442494    temp1.SetXTitle("MMcEvtBasic.fTelescopeTheta*kRad2Deg");
    … …  
    477529        bins3->SetName("BinningTheta");
    478530    }
     531
    479532    return kTRUE;
    480533}
    481534
    482 void MJSpectrum::DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const
     535TArrayD MJSpectrum::DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const
    483536{
    484537    TH1D collarea(area.GetHEnergy());
    … …  
    546599    f.SetParameter(1, 1.9e-6);
    547600    f.SetLineColor(kGreen);
    548     spectrum.Fit(&f, "NI", "", 55, 2e4);
     601    spectrum.Fit(&f, "NIM", "", 55, 2e4);
    549602    f.DrawCopy("same");
    550603
    … …  
    558611     tex.DrawLatex(2e2, 7e-5, str);
    559612     */
    560 }
    561 
    562 Bool_t MJSpectrum::PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot) const
    563 {
    564     cout << name << endl;
     613
     614    TArrayD res(2);
     615    res[0] = f.GetParameter(0);
     616    res[1] = f.GetParameter(1);
     617    return res;
     618}
     619
     620Bool_t MJSpectrum::PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot, Double_t scale) const
     621{
    565622    TString same(name);
    566623    same += "Same";
    … …  
    587644
    588645    const MAlphaFitter *fit = (MAlphaFitter*)plist.FindObject("MAlphaFitter");
    589     const Double_t scale = fit ? fit->GetScaleFactor() : 1;
     646    const Double_t ascale = fit ? fit->GetScaleFactor() : 1;
    590647
    591648    gPad->SetBorderMode(0);
    592649    h2->SetLineColor(kBlack);
    593650    h3->SetLineColor(kBlue);
    594     h2->Add(h1, -scale);
    595 
    596     h2->Scale(1./h2->Integral());
    597     h3->Scale(1./h3->Integral());
     651    h2->Add(h1, -ascale);
     652
     653    //h2->Scale(1./ontime);   //h2->Integral());
     654    h3->Scale(scale);    //h3->Integral());
    598655
    599656    h2->SetMaximum(1.05*TMath::Max(h2->GetMaximum(), h3->GetMaximum()));
    … …  
    609666}
    610667
    611 Bool_t MJSpectrum::DisplaySize(MParList &plist) const
     668Bool_t MJSpectrum::DisplaySize(MParList &plist, Double_t scale) const
    612669{
    613670    *fLog << inf << "Reading from file: " << fPathIn << endl;
    … …  
    644701
    645702    excess->SetTitle("Number of excess events vs Size (data, mc/blue)");
    646     excess->Scale(1./excess->Integral());
    647     excess = excess->DrawCopy();
     703    excess = excess->DrawCopy("E2");
    648704    // Don't do this on the original object!
    649705    excess->SetStats(kFALSE);
     706    excess->SetMarkerStyle(kFullDotMedium);
     707    excess->SetFillColor(kBlack);
     708    excess->SetFillStyle(0);
     709    excess->SetName("Excess  ");
     710    excess->SetDirectory(0);
    650711
    651712    TObject *o=0;
    652     if ((o=plist.FindObject("ExcessSize")))
     713    if ((o=plist.FindObject("ExcessMC")))
    653714    {
    654715        TH1 *histsel = (TH1F*)o->FindObject("");
    655716        if (histsel)
    656717        {
    657             histsel->Scale(1./histsel->Integral());
     718            if (scale<0)
     719                scale = excess->Integral()/histsel->Integral();
     720
     721            histsel->Scale(scale);
    658722            histsel->SetLineColor(kBlue);
    659723            histsel->SetBit(kCanDelete);
    660             histsel = histsel->DrawCopy("same");
     724            histsel = histsel->DrawCopy("E1 same");
    661725            // Don't do this on the original object!
    662726            histsel->SetStats(kFALSE);
     727
     728            fLog->Separator("Kolmogorov Test");
     729            histsel->KolmogorovTest(excess, "DX");
     730            fLog->Separator("Chi^2 Test");
     731            histsel->Chi2Test(excess, "P");
    663732        }
    664733    }
    … …  
    666735    // -------------- Comparison of Image Parameters --------------
    667736    c.cd(2);
    668     PlotSame(arr, plist, "Dist",   "HilSrc",  "MHHilSrcMCPost");
     737    PlotSame(arr, plist, "Dist",   "HilSrc",  "MHHilSrcMCPost", scale);
    669738
    670739    c.cd(3);
    671     PlotSame(arr, plist, "Length", "PostCut", "MHHillasMCPost");
     740    PlotSame(arr, plist, "Length", "PostCut", "MHHillasMCPost", scale);
    672741
    673742    c.cd(4);
    674     PlotSame(arr, plist, "M3l",    "HilExt",  "MHHilExtMCPost");
     743    PlotSame(arr, plist, "M3l",    "HilExt",  "MHHilExtMCPost", scale);
    675744
    676745    c.cd(5);
    677     PlotSame(arr, plist, "Conc1",  "NewPar",  "MHNewParMCPost");
     746    PlotSame(arr, plist, "Conc1",  "NewPar",  "MHNewParMCPost", scale);
    678747
    679748    c.cd(6);
    680     PlotSame(arr, plist, "Width",  "PostCut", "MHHillasMCPost");
     749    PlotSame(arr, plist, "Width",  "PostCut", "MHHillasMCPost", scale);
    681750
    682751    return kTRUE;
    … …  
    734803    }
    735804
     805    MMcSpectrumWeight weight;
     806    if (!InitWeighting(set, weight))
     807        return kFALSE;
     808
    736809    PrintSetup(fit);
    737810    bins3.SetEdges(temp1, 'x');
    738811
    739812    TH1D temp2(temp1);
    740     if (!ReadOrigMCDistribution(set, temp2))
     813    if (!ReadOrigMCDistribution(set, temp2, weight))
    741814        return kFALSE;
    742815
    … …  
    754827        hist.UseCurrentStyle();
    755828        MH::SetBinning(&hist, &bins3/*temp1.GetXaxis()*/, &bins2/*excess.GetXaxis()*/);
    756         if (!ReadOrigMCDistribution(set, hist))
     829        if (!ReadOrigMCDistribution(set, hist, weight))
    757830            return kFALSE;
    758831
    … …  
    762835                for (int x=0; x<hist.GetNbinsX(); x++)
    763836                    hist.SetBinContent(x, y, hist.GetBinContent(x, y)*temp1.GetBinContent(x));
     837            //hist.SetEntries(hist.Integral());
    764838        }
    765839    }
    766840    else
    767         if (!IntermediateLoop(plist, mh1, temp1, set))
     841    {
     842        weight.SetNameMcEvt("MMcEvtBasic");
     843        if (!IntermediateLoop(plist, mh1, temp1, set, weight))
    768844            return kFALSE;
     845        weight.SetNameMcEvt();
     846    }
    769847
    770848    DisplayResult(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
    … …  
    813891    MFillH fill3(&area, "", "FillCollectionArea");
    814892    MFillH fill4(&hest, "", "FillEnergyEst");
     893    fill3.SetWeight();
     894    fill4.SetWeight();
    815895
    816896    MH3 hsize("MHillas.fSize");
    817     //MH3 henergy("MEnergyEst.fVal");
    818     hsize.SetName("ExcessSize");
    819     //henergy.SetName("EnergyEst");
    820     MBinning bins(size, "BinningExcessSize");
     897    hsize.SetName("ExcessMC");
     898    hsize.Sumw2();
     899
     900    MBinning bins(size, "BinningExcessMC");
    821901    plist.AddToList(&hsize);
    822     //plist.AddToList(&henergy);
    823902    plist.AddToList(&bins);
    824903
    … …  
    830909    MFillH fill6a("MHImgParMCPost [MHImagePar]",    "MImagePar",    "FillImgParPost");
    831910    MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
    832     MFillH fill8a("ExcessSize     [MH3]",           "",             "FillExcessSize");
    833     //MFillH fill9a("EnergyEst      [MH3]",           "",             "FillExcessEEst");
     911    MFillH fill8a("ExcessMC       [MH3]",           "",             "FillExcessMC");
    834912    fill1a.SetNameTab("PreCut");
    835913    fill2a.SetNameTab("PostCut");
    … …  
    840918    fill7a.SetNameTab("NewPar");
    841919    fill8a.SetBit(MFillH::kDoNotDisplay);
    842     //fill9a.SetBit(MFillH::kDoNotDisplay);
     920    fill1a.SetWeight();
     921    fill2a.SetWeight();
     922    fill3a.SetWeight();
     923    fill4a.SetWeight();
     924    fill5a.SetWeight();
     925    fill6a.SetWeight();
     926    fill7a.SetWeight();
     927    fill8a.SetWeight();
    843928
    844929    MEnergyEstimate est;
    … …  
    852937    tlist2.AddToList(&hcalc1);
    853938    tlist2.AddToList(&hcalc2);
     939    tlist2.AddToList(&weight);
    854940    tlist2.AddToList(&fill1a);
    855941    tlist2.AddToList(fCut0);
    … …  
    891977    // -------------------------- Spectrum ----------------------------
    892978
    893     DisplaySpectrum(area, excess, hest, ontime);
    894     DisplaySize(plist);
    895 
     979    // Calculate and display spectrum (N/TeVsm^2 at 1TeV)
     980    TArrayD res(DisplaySpectrum(area, excess, hest, ontime));
     981
     982    // Spectrum fitted (convert res[1] from TeV to GeV)
     983    TF1 flx("flx", Form("%e*pow(x/1000, %f)", res[1]/1000, res[0]));
     984
     985    // Number of events this spectrum would produce per s and m^2
     986    Double_t n = flx.Integral(weight.GetEnergyMin(), weight.GetEnergyMax());
     987
     988    // scale with effective collection area to get the event rate (N/s)
     989    // scale with the effective observation time to absolute observed events
     990    n *= area.GetCollectionAreaAbs()*ontime; // N
     991
     992    // Now calculate the scale factor from the number of events
     993    // produced and the number of events which should have been
     994    // observed with our telescope in the time ontime
     995    const Double_t scale = n/area.GetEntries();
     996
     997    // Print normalization constant
     998    cout << "MC normalization factor:  " << scale << endl;
     999
     1000    // Overlay normalized plots
     1001    DisplaySize(plist, scale);
     1002
     1003    // check if output should be written
    8961004    if (!fPathOut.IsNull())
    8971005        fDisplay->SaveAsRoot(fPathOut);

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

    r6978 r7094  
    1818class MStatusArray;
    1919class MHCollectionArea;
     20class MMcSpectrumWeight;
    2021
    2122class MJSpectrum : public MJob
    … …  
    3536    Bool_t  ReadTask(MTask* &task, const char *name) const;
    3637    Float_t ReadInput(MParList &plist, TH1D &h1, TH1D &size);
    37     Bool_t  ReadOrigMCDistribution(const MDataSet &set, TH1 &h) const;
     38    Bool_t  ReadOrigMCDistribution(const MDataSet &set, TH1 &h, MMcSpectrumWeight &w) const;
    3839    Bool_t  GetThetaDistribution(TH1D &temp1, TH1D &temp2) const;
    3940    Bool_t  Refill(MParList &plist, TH1D &h) const;
     41    Bool_t  InitWeighting(const MDataSet &set, MMcSpectrumWeight &w) const;
    4042
    4143    // Display Output
    4244    void    PrintSetup(const MAlphaFitter &fit) const;
    4345    void    DisplayResult(const TH2D &mh1) const;
    44     Bool_t  IntermediateLoop(MParList &plist, MH3 &h1, TH1D &temp1, const MDataSet &set) const;
    45     void    DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const;
    46     Bool_t  DisplaySize(MParList &plist) const;
    47     Bool_t  PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot) const;
     46    Bool_t  IntermediateLoop(MParList &plist, MH3 &h1, TH1D &temp1, const MDataSet &set, MMcSpectrumWeight &w) const;
     47    TArrayD DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const;
     48    Bool_t  DisplaySize(MParList &plist, Double_t scale) const;
     49    Bool_t  PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot, Double_t scale) const;
    4850
    4951public:

  • trunk/MagicSoft/Mars/sponde.rc

    r6977 r7094  
    11EstimateEnergy.Rule: (0.380075+(MPointingPos.fZd*MPointingPos.fZd*0.00109028))*pow(MHillas.fSize,0.892462)
     2
     3#MMcSpectrumWeight.NewSlope: -2.26

  • trunk/MagicSoft/include-Classes/MMcFormat/MMcEvt.cxx

    r5321 r7094  
    1 #include "MMcEvt.hxx"
    2 
    3 #include "MLog.h"
    4 #include "MLogManip.h"
    5 
    6 //==========
    7 // MMcEvt
    8 //   
     1/* ======================================================================== *\
     2!
     3! *
     4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
     5! * Software. It is distributed to you in the hope that it can be a useful
     6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
     7! * It is distributed WITHOUT ANY WARRANTY.
     8! *
     9! * Permission to use, copy, modify and distribute this software and its
     10! * documentation for any purpose is hereby granted without fee,
     11! * provided that the above copyright notice appear in all copies and
     12! * that both that copyright notice and this permission notice appear
     13! * in supporting documentation. It is provided "as is" without express
     14! * or implied warranty.
     15! *
     16!
     17!
     18!   Author(s):
     19!
     20!   Copyright: MAGIC Software Development, 2000-2005
     21!
     22!
     23\* ======================================================================== */
     24
     25/////////////////////////////////////////////////////////////////////////////
     26//
     27//  MMcEvt
     28//
    929// This class handles and contains the MonteCarlo information
    1030// with which the events have been generated
    … …  
    1939// So, fTelescopeTheta=90, fTelescopePhi = 0 means the telescope is
    2040// pointing horizontally towards South. For an explanation, see also
    21 // TDAS 02-11.
    22 //
     41// TDAS 02-11.
     42//
     43// Version 4:
     44//   - Added member fFadcTimeJitter
     45//
     46// Version 5:
     47//   - removed fPartId, fEnergy, fImpact, fTelescopeTheta, fTelescopePhi
     48//   - derives now from MMcEvtBasic which contains all these values
     49//   - moved ParticleId_t to base class MMcEvtBasic
     50//
     51/////////////////////////////////////////////////////////////////////////////
     52#include "MMcEvt.hxx"
     53
     54#include "MLog.h"
     55#include "MLogManip.h"
     56
    2357ClassImp(MMcEvt);
    2458
    … …  
    2660
    2761
     62// --------------------------------------------------------------------------
     63//
     64// Default constructor. Calls Clear()
     65//
    2866MMcEvt::MMcEvt()
    2967{
    30     //
    31     //  default constructor
    32     //  set all values to zero
    3368    fName  = "MMcEvt";
    3469    fTitle = "Event info from Monte Carlo";
    … …  
    3772}
    3873
    39 MMcEvt::MMcEvt( UInt_t  fEvtNum,
    40                 UShort_t usPId,
    41                 Float_t  fEner,
    42                 Float_t  fThi0,
    43                 Float_t  fFirTar,
    44                 Float_t  fzFirInt,
    45                 Float_t  fThet,
    46                 Float_t  fPhii,
    47                 Float_t  fCorD,
    48                 Float_t  fCorX,
    49                 Float_t  fCorY,
    50                 Float_t  fImpa,
    51                 Float_t  fTPhii,
    52                 Float_t  fTThet,
    53                 Float_t  fTFirst,
    54                 Float_t  fTLast,
    55                 Float_t  fL_Nmax,
    56                 Float_t  fL_t0,
    57                 Float_t  fL_tmax,
    58                 Float_t  fL_a,
    59                 Float_t  fL_b,
    60                 Float_t  fL_c,
    61                 Float_t  fL_chi2,
    62                 UInt_t   uiPin,
    63                 UInt_t   uiPat, 
    64                 UInt_t   uiPre,
    65                 UInt_t   uiPco, 
    66                 UInt_t   uiPelS,
    67                 UInt_t   uiPelC,
    68                 Float_t  elec,
    69                 Float_t  muon,
    70                 Float_t  other,
    71                 Float_t  fadc_jitter) {
    72 
     74// --------------------------------------------------------------------------
     75//
     76// Constructor. Use this to set all data members
     77//
     78// THIS FUNCTION IS FOR THE SIMULATION OLNY.
     79// DON'T USE THIS MEMBERFUNCTION IN THE ANALYSIS.
     80//
     81MMcEvt::MMcEvt(UInt_t  fEvtNum,    ParticleId_t usPId, Float_t  fEner,
     82               Float_t  fThi0,     Float_t  fFirTar,   Float_t  fzFirInt,
     83               Float_t  fThet,     Float_t  fPhii,     Float_t  fCorD,
     84               Float_t  fCorX,     Float_t  fCorY,     Float_t  fImpa,
     85               Float_t  fTPhii,    Float_t  fTThet,    Float_t  fTFirst,
     86               Float_t  fTLast,    Float_t  fL_Nmax,   Float_t  fL_t0,
     87               Float_t  fL_tmax,   Float_t  fL_a,      Float_t  fL_b,
     88               Float_t  fL_c,      Float_t  fL_chi2,   UInt_t   uiPin,
     89               UInt_t   uiPat,     UInt_t   uiPre,     UInt_t   uiPco,
     90               UInt_t   uiPelS,    UInt_t   uiPelC,    Float_t  elec,
     91               Float_t  muon,      Float_t  other,     Float_t  fadc_jitter)
     92{
    7393    fName  = "MMcEvt";
    7494    fTitle = "Event info from Monte Carlo";
    75   //
    76   //  constuctor II
    77   //
    78   //  All datamembers are parameters.
    79   //
    80   //  Don't use this memberfunction in analysis
    81   // 
    82 
    83   fEvtNumber = fEvtNum;
    84   fPartId = usPId  ;
    85   fEnergy  = fEner  ;
    86   fThick0 = fThi0;
    87   fFirstTarget = fFirTar;
    88   fZFirstInteraction = fzFirInt;
    89 
    90   fTheta   = fThet ;
    91   fPhi     = fPhii ;
    92 
    93   fCoreD   = fCorD ;
    94   fCoreX   = fCorX ;
    95   fCoreY   = fCorY ;
    96   fImpact  = fImpa ;
    97 
    98   fTelescopePhi = fTPhii;
    99   fTelescopeTheta = fTThet;
    100   fTimeFirst = fTFirst;
    101   fTimeLast = fTLast;
    102   fLongiNmax = fL_Nmax;
    103   fLongit0 = fL_t0;
    104   fLongitmax = fL_tmax;
    105   fLongia = fL_a;
    106   fLongib = fL_b;
    107   fLongic = fL_c;
    108   fLongichi2 = fL_chi2;
    109 
    110 
    111   fPhotIni      = uiPin ;
    112   fPassPhotAtm  = uiPat ;
    113   fPassPhotRef  = uiPre ;
    114   fPassPhotCone = uiPco ;
    115   fPhotElfromShower = uiPelS ;
    116   fPhotElinCamera   = uiPelC ;
    117 
    118   fElecCphFraction=elec;
    119   fMuonCphFraction=muon;
    120   fOtherCphFraction=other;
    121 
    122   fFadcTimeJitter = fadc_jitter;
    123 }
    124 
    125 
    126 
    127 MMcEvt::~MMcEvt() {
    128   //
    129   //  default destructor
    130   //
    131 }
    132 
    133 
    134 
    135 
     95
     96    Fill(fEvtNum, usPId, fEner, fThi0, fFirTar, fzFirInt, fThet,
     97         fPhii, fCorD, fCorX, fCorY, fImpa, fTPhii, fTThet, fTFirst,
     98         fTLast, fL_Nmax, fL_t0, fL_tmax, fL_a, fL_b, fL_c, fL_chi2,
     99         uiPin, uiPat, uiPre, uiPco, uiPelS, uiPelC, elec, muon, other,
     100         fadc_jitter);
     101}
     102
     103// --------------------------------------------------------------------------
     104//
     105//  reset all values to values as nonsense as possible
     106//
    136107void MMcEvt::Clear(Option_t *opt)
    137108{
    138     //
    139     //  reset all values to values as nonsense as possible
    140     //
    141     fPartId = 0;
     109    fPartId = kUNDEFINED;
    142110    fEnergy = -1;
    143111
    … …  
    162130}
    163131
    164 void MMcEvt::Fill( UInt_t  fEvtNum,
    165                    UShort_t usPId,
    166                    Float_t  fEner,
    167                    Float_t  fThi0,
    168                    Float_t  fFirTar,
    169                    Float_t  fzFirInt,
    170                    Float_t  fThet,
    171                    Float_t  fPhii,
    172                    Float_t  fCorD,
    173                    Float_t  fCorX,
    174                    Float_t  fCorY,
    175                    Float_t  fImpa,
    176                    Float_t  fTPhii,
    177                    Float_t  fTThet,
    178                    Float_t  fTFirst,
    179                    Float_t  fTLast,
    180                    Float_t  fL_Nmax,
    181                    Float_t  fL_t0,
    182                    Float_t  fL_tmax,
    183                    Float_t  fL_a,
    184                    Float_t  fL_b,
    185                    Float_t  fL_c,
    186                    Float_t  fL_chi2,
    187                    UInt_t   uiPin,
    188                    UInt_t   uiPat, 
    189                    UInt_t   uiPre,
    190                    UInt_t   uiPco, 
    191                    UInt_t   uiPelS, 
    192                    UInt_t   uiPelC,
    193                    Float_t  elec,
    194                    Float_t  muon,
    195                    Float_t  other,
    196                    Float_t  fadc_jitter) {
    197   //
    198   //  All datamembers are filled with the correspondin parameters.
    199   //
    200   //  Don't use this memberfunction in analysis
    201   // 
    202 
    203   fEvtNumber = fEvtNum;
    204   fPartId = usPId  ;
    205   fEnergy = fEner  ;
    206   fThick0 = fThi0;
    207   fFirstTarget = fFirTar;
    208   fZFirstInteraction = fzFirInt;
    209 
    210   fTheta  = fThet ;
    211   fPhi    = fPhii ;
    212 
    213   fCoreD  = fCorD ;
    214   fCoreX  = fCorX ;
    215   fCoreY  = fCorY ;
    216   fImpact = fImpa ;
    217 
    218   fTelescopePhi = fTPhii;
    219   fTelescopeTheta = fTThet;
    220   fTimeFirst = fTFirst;
    221   fTimeLast = fTLast;
    222   fLongiNmax = fL_Nmax;
    223   fLongit0 = fL_t0;
    224   fLongitmax = fL_tmax;
    225   fLongia = fL_a;
    226   fLongib = fL_b;
    227   fLongic = fL_c;
    228   fLongichi2 = fL_chi2;
    229 
    230   fPhotIni      = uiPin ;
    231   fPassPhotAtm  = fPhotIni-uiPat ;
    232   fPassPhotRef  = fPassPhotAtm-uiPre ;
    233   fPassPhotCone = uiPco ;
    234   fPhotElfromShower = uiPelS ;
    235   fPhotElinCamera = uiPelC ;
    236 
    237   fElecCphFraction=elec;
    238   fMuonCphFraction=muon;
    239   fOtherCphFraction=other;
    240 
    241   fFadcTimeJitter = fadc_jitter;
    242 }
    243 
    244 /*
    245 void MMcEvt::AsciiWrite(ofstream &fout) const
    246 {
    247     fout << fEnergy << " ";
    248     fout << fTheta ;
    249 }
    250 */
     132// --------------------------------------------------------------------------
     133//
     134// Use this to set all data members
     135//
     136// THIS FUNCTION IS FOR THE SIMULATION OLNY.
     137// DON'T USE THIS MEMBERFUNCTION IN THE ANALYSIS.
     138//
     139void MMcEvt::Fill( UInt_t  fEvtNum,    ParticleId_t usPId, Float_t  fEner,
     140                   Float_t  fThi0,     Float_t  fFirTar,   Float_t  fzFirInt,
     141                   Float_t  fThet,     Float_t  fPhii,     Float_t  fCorD,
     142                   Float_t  fCorX,     Float_t  fCorY,     Float_t  fImpa,
     143                   Float_t  fTPhii,    Float_t  fTThet,    Float_t  fTFirst,
     144                   Float_t  fTLast,    Float_t  fL_Nmax,   Float_t  fL_t0,
     145                   Float_t  fL_tmax,   Float_t  fL_a,      Float_t  fL_b,
     146                   Float_t  fL_c,      Float_t  fL_chi2,   UInt_t   uiPin,
     147                   UInt_t   uiPat,     UInt_t   uiPre,     UInt_t   uiPco, 
     148                   UInt_t   uiPelS,    UInt_t   uiPelC,    Float_t  elec,
     149                   Float_t  muon,      Float_t  other,     Float_t  fadc_jitter)
     150{
     151    //
     152    //  All datamembers are filled with the correspondin parameters.
     153    //
     154    //  Don't use this memberfunction in analysis
     155    //
     156    fEvtNumber = fEvtNum;
     157    fPartId = usPId  ;
     158    fEnergy = fEner  ;
     159    fThick0 = fThi0;
     160    fFirstTarget = fFirTar;
     161    fZFirstInteraction = fzFirInt;
     162
     163    fTheta  = fThet ;
     164    fPhi    = fPhii ;
     165
     166    fCoreD  = fCorD ;
     167    fCoreX  = fCorX ;
     168    fCoreY  = fCorY ;
     169    fImpact = fImpa ;
     170
     171    fTelescopePhi = fTPhii;
     172    fTelescopeTheta = fTThet;
     173    fTimeFirst = fTFirst;
     174    fTimeLast = fTLast;
     175    fLongiNmax = fL_Nmax;
     176    fLongit0 = fL_t0;
     177    fLongitmax = fL_tmax;
     178    fLongia = fL_a;
     179    fLongib = fL_b;
     180    fLongic = fL_c;
     181    fLongichi2 = fL_chi2;
     182
     183    fPhotIni      = uiPin ;
     184    fPassPhotAtm  = fPhotIni-uiPat ;
     185    fPassPhotRef  = fPassPhotAtm-uiPre ;
     186    fPassPhotCone = uiPco ;
     187    fPhotElfromShower = uiPelS ;
     188    fPhotElinCamera = uiPelC ;
     189
     190    fElecCphFraction=elec;
     191    fMuonCphFraction=muon;
     192    fOtherCphFraction=other;
     193
     194    fFadcTimeJitter = fadc_jitter;
     195}
    251196
    252197// --------------------------------------------------------------------------
    … …  
    260205void MMcEvt::Print(Option_t *opt) const
    261206{
    262     //
    263     //  print out the data member on screen
    264     //
     207    MMcEvtBasic::Print(opt);
     208
    265209    TString str(opt);
    266210    if (str.IsNull())
    267211    {
    268         *fLog << all << endl;
    269         *fLog << "Monte Carlo output:" << endl;
    270         *fLog << " Particle Id:    ";
    271         switch(fPartId)
    272         {
    273         case kGAMMA:
    274             *fLog << "Gamma" << endl;
    275             break;
    276         case kPROTON:
    277             *fLog << "Proton" << endl;
    278             break;
    279         case kHELIUM:
    280             *fLog << "Helium" << endl;
    281             break;
    282         }
    283         *fLog << " Energy:         " << fEnergy << "GeV" << endl;
    284         *fLog << " Impactpar.:     " << fImpact/100 << "m" << endl;
    285212        *fLog << " Photoelectrons: " << fPhotElfromShower << endl;
    286         *fLog << endl;
    287213        return;
    288214    }
    289     if (str.Contains("id", TString::kIgnoreCase))
    290         switch(fPartId)
    291         {
    292         case kGAMMA:
    293             *fLog << "Particle: Gamma" << endl;
    294             break;
    295         case kPROTON:
    296             *fLog << "Particle: Proton" << endl;
    297             break;
    298         case kHELIUM:
    299             *fLog << "Particle: Helium" << endl;
    300             break;
    301         }
    302     if (str.Contains("energy", TString::kIgnoreCase))
    303         *fLog << "Energy: " << fEnergy << "GeV" << endl;
    304     if (str.Contains("impact", TString::kIgnoreCase))
    305         *fLog << "Impact: " << fImpact << "cm" << endl;
    306 }
     215}

  • trunk/MagicSoft/include-Classes/MMcFormat/MMcEvt.hxx

    r5423 r7094  
    1 #ifndef __MMcEvt__
    2 #define __MMcEvt__
     1#ifndef MARS_MMcEvt
     2#define MARS_MMcEvt
    33
    4 #ifndef MARS_MParContainer
    5 #include "MParContainer.h"
     4#ifndef MARS_MMcEvtBasic
     5#include "MMcEvtBasic.h"
    66#endif
    77
    8 //
    9 // Version 4:
    10 //   Added member fFadcTimeJitter
    11 //
    128
    13 class MMcEvt : public MParContainer
     9class MMcEvt : public MMcEvtBasic
    1410{
     11private:
     12    UInt_t  fEvtNumber;
     13    Float_t fThick0;             // [g/cm2]
     14    Float_t fFirstTarget;        // []
     15    Float_t fZFirstInteraction;  // [cm]
     16
     17    Float_t fTheta;              // [rad] Theta angle of event
     18    Float_t fPhi;                // [rad] Phi angle of event (see class description)
     19
     20    Float_t fCoreD;              // [cm] Core d pos
     21    Float_t fCoreX;              // [cm] Core x pos
     22    Float_t fCoreY;              // [cm] Core y pos
     23
     24    // Up to here, the info from the CORSIKA event header.
     25
     26    // Time of first and last photon:
     27    Float_t fTimeFirst;          // [ns]
     28    Float_t fTimeLast;           // [ns]
     29
     30    // 6 parameters and chi2 of the NKG fit to the longitudinal
     31    // particle distribution. See CORSIKA manual for explanation,
     32    // section 4.42 "Longitudinal shower development":
     33    //
     34    Float_t fLongiNmax;          // [particles]
     35    Float_t fLongit0;            // [g/cm2]
     36    Float_t fLongitmax;          // [g/cm2]
     37    Float_t fLongia;             // [g/cm2]
     38    Float_t fLongib;             // []
     39    Float_t fLongic;             // [cm2/g]
     40    Float_t fLongichi2;
     41
     42    UInt_t fPhotIni;             // [ph] Initial number of photons
     43    UInt_t fPassPhotAtm;         // [ph] Passed atmosphere
     44    UInt_t fPassPhotRef;         // [ph] Passed reflector(reflectivity + effective area)
     45    UInt_t fPassPhotCone;        // [ph]  Within any valid pixel, before plexiglas
     46    UInt_t fPhotElfromShower;    // [phe] Passed qe, coming from the shower
     47    UInt_t fPhotElinCamera;      // [phe] usPhotElfromShower + mean of phe from NSB
     48
     49    // Now follow the fraction of photons reaching the camera produced by
     50    // electrons, muons and other particles respectively:
     51
     52    Float_t  fElecCphFraction;
     53    Float_t  fMuonCphFraction;
     54    Float_t  fOtherCphFraction;
     55
     56    Float_t  fFadcTimeJitter;
     57
    1558public:
    16     //
    17     //     ParticleId for Monte Carlo simulation
    18     //
    19     enum ParticleId_t
    20     {
    21         kGAMMA    =  1,
    22         kPOSITRON =  2,
    23         kELECTRON =  3,
    24         kANTIMUON =  5,
    25         kMUON     =  6,
    26         kPI0      =  7,
    27         kNEUTRON  = 13,
    28         kPROTON =   14,
    29         kHELIUM =  402,
    30         kOXYGEN = 1608,
    31         kIRON   = 5626
    32     };
     59    MMcEvt();
     60    MMcEvt(UInt_t, ParticleId_t, Float_t, Float_t, Float_t,
     61           Float_t, Float_t, Float_t, Float_t, Float_t, Float_t,
     62           Float_t, Float_t, Float_t, Float_t, Float_t, Float_t,
     63           Float_t, Float_t, Float_t, Float_t, Float_t, Float_t,
     64           UInt_t, UInt_t, UInt_t, UInt_t, UInt_t, UInt_t,
     65           Float_t, Float_t, Float_t, Float_t) ;
    3366
    34 private:
    35   UInt_t      fEvtNumber;
    36   UShort_t     fPartId;             // Type of particle
    37   Float_t      fEnergy;             // [GeV] Energy
    38   Float_t      fThick0;             // [g/cm2]
    39   Float_t      fFirstTarget;        // []
    40   Float_t      fZFirstInteraction;  // [cm]
     67    // Getter
     68    UInt_t  GetEvtNumber() const { return fEvtNumber; }  //Get Event Number
     69    Float_t GetTheta() const { return fTheta; }          //Get Theta angle
     70    Float_t GetPhi() const { return fPhi ;  }            //Get Phi angle
    4171
    42   Float_t fTheta;           // [rad] Theta angle of event
    43   Float_t fPhi;             // [rad] Phi angle of event (see class description)
     72    Float_t GetCoreX() const { return fCoreX; }          //Get Core x pos
     73    Float_t GetCoreY() const { return fCoreY; }          //Get Core y pos
    4474
    45   Float_t fCoreD;           // [cm] Core d pos
    46   Float_t fCoreX;           // [cm] Core x pos
    47   Float_t fCoreY;           // [cm] Core y pos
    48   Float_t fImpact;          // [cm] impact parameter
     75    UInt_t  GetPhotIni() const { return fPhotIni; }           //Get Initial photons
     76    UInt_t  GetPassPhotAtm() const { return fPassPhotAtm;}    //Get Passed atmosphere
     77    UInt_t  GetPassPhotRef() const { return fPassPhotRef; }   //Get Passed reflector
     78    UInt_t  GetPassPhotCone() const { return fPassPhotCone; } //Get Passed glas
     79    UInt_t  GetPhotElfromShower() const { return fPhotElfromShower; }   //Get Passed qe from shower
     80    UInt_t  GetPhotElinCamera() const { return fPhotElinCamera; }       //Get Passed qe total
     81    Float_t GetZFirstInteraction() const { return fZFirstInteraction; }
    4982
    50   // Up to here, the info from the CORSIKA event header.
     83    Float_t GetOtherCphFraction() const { return fOtherCphFraction; }
    5184
    52   // Telescope orientation:
    53   Float_t       fTelescopePhi;    // [rad] (see class description)
    54   Float_t       fTelescopeTheta;  // [rad]
     85    Float_t GetLongiNmax() const { return fLongiNmax; }
     86    Float_t GetLongia()    const { return fLongia; }
     87    Float_t GetLongib()    const { return fLongib; }
     88    Float_t GetLongic()    const { return fLongic; }
     89    Float_t GetLongichi2() const { return fLongichi2; }
     90    Float_t GetLongit0()   const { return fLongit0; }
     91    Float_t GetLongitmax() const { return fLongitmax; }
    5592
    56   // Time of first and last photon:
    57   Float_t      fTimeFirst;   // [ns]
    58   Float_t      fTimeLast;    // [ns]
     93    Float_t GetFadcTimeJitter() const { return fFadcTimeJitter; }
    5994
    60   // 6 parameters and chi2 of the NKG fit to the longitudinal
    61   // particle distribution. See CORSIKA manual for explanation,
    62   // section 4.42 "Longitudinal shower development":
    63   //
    64   Float_t       fLongiNmax;   // [particles]
    65   Float_t       fLongit0;     // [g/cm2]
    66   Float_t       fLongitmax;   // [g/cm2]
    67   Float_t       fLongia;      // [g/cm2]
    68   Float_t       fLongib;      // []
    69   Float_t       fLongic;      // [cm2/g]
    70   Float_t       fLongichi2;
     95    Float_t GetMuonCphFraction() const { return fMuonCphFraction; }
    7196
    72   UInt_t fPhotIni;          // [ph] Initial number of photons
    73   UInt_t fPassPhotAtm;      // [ph] Passed atmosphere
    74   UInt_t fPassPhotRef;      // [ph] Passed reflector(reflectivity + effective area)
    75   UInt_t fPassPhotCone;     // [ph]  Within any valid pixel, before plexiglas
    76   UInt_t fPhotElfromShower; // [phe] Passed qe, coming from the shower
    77   UInt_t fPhotElinCamera;   // [phe] usPhotElfromShower + mean of phe
    78                             // from NSB
     97    // Setter
     98    void SetTheta(Float_t Theta) { fTheta=Theta; }                //Set Theta angle
     99    void SetPhi(Float_t Phi)     { fPhi=Phi;  }                   //Set Phi angle
     100    void SetCoreD(Float_t CoreD) { fCoreD=CoreD; }                //Set Core d pos
     101    void SetCoreX(Float_t CoreX) { fCoreX=CoreX; }                //Set Core x pos
     102    void SetCoreY(Float_t CoreY) { fCoreY=CoreY; }                //Set Core y pos
    79103
    80   // Now follow the fraction of photons reaching the camera produced by
    81   // electrons, muons and other particles respectively:
     104    void Fill( UInt_t, ParticleId_t, Float_t, Float_t, Float_t,
     105               Float_t, Float_t, Float_t, Float_t, Float_t, Float_t,
     106               Float_t, Float_t, Float_t, Float_t, Float_t, Float_t,
     107               Float_t, Float_t, Float_t, Float_t, Float_t, Float_t,
     108               UInt_t, UInt_t, UInt_t, UInt_t, UInt_t, UInt_t,
     109               Float_t, Float_t, Float_t, Float_t) ;
    82110
    83   Float_t  fElecCphFraction;
    84   Float_t  fMuonCphFraction;
    85   Float_t  fOtherCphFraction;
    86  
    87   Float_t  fFadcTimeJitter;
     111    // TObject
     112    void Print(Option_t *opt=NULL) const;
     113    void Clear(Option_t *opt=NULL);
    88114
    89  public:
    90   MMcEvt() ;
    91  
    92   MMcEvt( UInt_t, UShort_t, Float_t, Float_t, Float_t,
    93           Float_t, Float_t, Float_t, Float_t, Float_t, Float_t,
    94           Float_t, Float_t, Float_t, Float_t, Float_t, Float_t,
    95           Float_t, Float_t, Float_t, Float_t, Float_t, Float_t,
    96           UInt_t, UInt_t, UInt_t, UInt_t, UInt_t, UInt_t,
    97           Float_t, Float_t, Float_t, Float_t) ;
    98  
    99   ~MMcEvt();
    100 
    101   void Clear(Option_t *opt=NULL);
    102 
    103   void Fill( UInt_t, UShort_t, Float_t, Float_t, Float_t,
    104              Float_t, Float_t, Float_t, Float_t, Float_t, Float_t,
    105              Float_t, Float_t, Float_t, Float_t, Float_t, Float_t,
    106              Float_t, Float_t, Float_t, Float_t, Float_t, Float_t,
    107              UInt_t, UInt_t, UInt_t, UInt_t, UInt_t, UInt_t,
    108              Float_t, Float_t, Float_t, Float_t) ;
    109 
    110   //virtual void AsciiWrite(ofstream &fout) const;
    111 
    112   void Print(Option_t *opt=NULL) const;
    113 
    114   UInt_t GetEvtNumber() const { return fEvtNumber; }  //Get Event Number
    115   Short_t GetPartId() const { return fPartId; }       //Get Type of particle
    116   Float_t GetEnergy() const { return fEnergy; }        //Get Energy
    117 
    118   Float_t GetTheta() const { return fTheta; }          //Get Theta angle
    119   Float_t GetPhi() const { return fPhi ;  }            //Get Phi angle
    120 
    121 /*    Float_t GetCoreD() { return fCoreD; }          //Get Core d pos */
    122   Float_t GetCoreX() { return fCoreX; }          //Get Core x pos
    123   Float_t GetCoreY() { return fCoreY; }          //Get Core y pos
    124   Float_t GetImpact() const { return fImpact;}         //Get impact parameter
    125 
    126   UInt_t GetPhotIni() { return fPhotIni; }           //Get Initial photons
    127   UInt_t GetPassPhotAtm() { return fPassPhotAtm;}    //Get Passed atmosphere
    128   UInt_t GetPassPhotRef() { return fPassPhotRef; }   //Get Passed reflector
    129   UInt_t GetPassPhotCone() { return fPassPhotCone; } //Get Passed glas
    130   UInt_t GetPhotElfromShower() { return fPhotElfromShower; }   //Get Passed qe from shower
    131   UInt_t GetPhotElinCamera() { return fPhotElinCamera; }   //Get Passed qe total
    132   Float_t GetZFirstInteraction() const { return fZFirstInteraction; }
    133 
    134   Float_t GetTelescopePhi() const { return fTelescopePhi; }
    135   Float_t GetTelescopeTheta() const { return fTelescopeTheta; }
    136   Float_t GetOtherCphFraction() const { return fOtherCphFraction; }
    137 
    138   Float_t GetLongiNmax() const { return fLongiNmax; }
    139   Float_t GetLongia()    const { return fLongia; }
    140   Float_t GetLongib()    const { return fLongib; }
    141   Float_t GetLongic()    const { return fLongic; }
    142   Float_t GetLongichi2() const { return fLongichi2; }
    143   Float_t GetLongit0()   const { return fLongit0; }
    144   Float_t GetLongitmax() const { return fLongitmax; }
    145 
    146   Float_t GetFadcTimeJitter() const { return fFadcTimeJitter; }
    147 
    148   Float_t GetMuonCphFraction() const { return fMuonCphFraction; }
    149 
    150   void SetPartId(Short_t PartId)
    151   {fPartId=PartId;}             //Set Type of particle
    152 
    153   TString GetParticleName() const
    154   {
    155       switch (fPartId)
    156       {
    157       case kGAMMA:    return "Gamma";
    158       case kPOSITRON: return "Positron";
    159       case kELECTRON: return "Electron";
    160       case kANTIMUON: return "Anti-Muon";
    161       case kMUON:     return "Muon";
    162       case kPI0:      return "Pi-0";
    163       case kNEUTRON:  return "Neutron";
    164       case kPROTON:   return "Proton";
    165       case kHELIUM:   return "Helium";
    166       case kOXYGEN:   return "Oxygen";
    167       case kIRON:     return "Iron";
    168       }
    169 
    170       return Form("Id:%d", fPartId);
    171   }
    172   TString GetParticleSymbol() const
    173   {
    174       switch (fPartId)
    175       {
    176       case kGAMMA:    return "\\gamma";
    177       case kPOSITRON: return "e^{+}";
    178       case kELECTRON: return "e^{-}";
    179       case kANTIMUON: return "\\mu^{+}";
    180       case kMUON:     return "\\mu^{-}";
    181       case kPI0:      return "\\pi^{0}";
    182       case kNEUTRON:  return "n";
    183       case kPROTON:   return "p";
    184       case kHELIUM:   return "He";
    185       case kOXYGEN:   return "O";
    186       case kIRON:     return "Fe";
    187       }
    188 
    189       return Form("Id:%d", fPartId);
    190   }
    191   TString GetEnergyStr() const
    192   {
    193       if (fEnergy>1000)
    194           return Form("%.1fTeV", fEnergy/1000);
    195 
    196       if (fEnergy>10)
    197           return Form("%dGeV", (Int_t)(fEnergy+.5));
    198 
    199       if (fEnergy>1)
    200           return Form("%.1fGeV", fEnergy);
    201 
    202       return Form("%dMeV", (Int_t)(fEnergy*1000+.5));
    203   }
    204 
    205   void SetEnergy(Float_t Energy)
    206   { fEnergy=Energy; }              //Set Energy
    207  
    208   void SetTheta(Float_t Theta)
    209   { fTheta=Theta; }                //Set Theta angle
    210 
    211   void SetPhi(Float_t Phi)
    212   { fPhi=Phi;  }                   //Set Phi angle
    213  
    214   void SetCoreD(Float_t CoreD)
    215   { fCoreD=CoreD; }                //Set Core d pos
    216 
    217   void SetCoreX(Float_t CoreX)
    218   { fCoreX=CoreX; }                //Set Core x pos
    219 
    220   void SetCoreY(Float_t CoreY )
    221   { fCoreY=CoreY; }                //Set Core y pos
    222 
    223   void SetImpact(Float_t Impact)
    224   { fImpact=Impact;}               //Set impact parameter
    225 
    226   // DO NOT USE THIS IS A WORKAROUND!
    227   void SetTelescopeTheta(Float_t Theta) { fTelescopeTheta=Theta; }
    228   void SetTelescopePhi(Float_t Phi)     { fTelescopePhi=Phi; }
    229 
    230  
    231 /*    void SetPhotIni(Short_t PhotIni)  */
    232 /*      { fPhotIni=PhotIni; }                 //Set Initial photons */
    233 /*    void SetPassPhotAtm(Short_t PassPhotAtm)  */
    234 /*      { fPassPhotAtm=PassPhotAtm;}         //Set Passed atmosphere */
    235 /*    void SetPassPhotRef(Short_t PassPhotRef)  */
    236 /*      { fPassPhotRef=PassPhotRef ; }       //Set Passed reflector */
    237 /*    void SetPassPhotCone(Short_t PhotCon)  */
    238 /*      { fPassPhotCone=PhotCon; }           //Set Passed glas */
    239 
    240 
    241   ClassDef(MMcEvt, 4)  //Stores Montecarlo Information of one event (eg. the energy)
    242 
     115    ClassDef(MMcEvt, 5)  //Stores Montecarlo Information of one event (eg. the energy)
    243116};
    244117

  • trunk/MagicSoft/include-Classes/MMcFormat/MMcEvtBasic.cc

    r6375 r7094  
    4040// South.
    4141//
     42//
     43// Version 1:
     44//  New container to keep the very basic informations on the
     45//  original MC events produced by Corsika
     46//
     47// Version 2:
     48//  - added typedef for ParticleId_t from MMcEvt
     49//  - replaced MMcEvt::ParticleId_t by ParticleId_t
     50//
    4251/////////////////////////////////////////////////////////////////////////////
    43 
    4452#include "MMcEvtBasic.h"
    4553
    … …  
    5260
    5361// --------------------------------------------------------------------------
    54 //  Default constructor set all values to zero
     62//
     63// Default constructor. Calls Clear()
    5564//
    5665MMcEvtBasic::MMcEvtBasic()
    57 {   
    58   fName  = "MMcEvtBasic";
    59   fTitle = "Basic event info from Monte Carlo";
     66{
     67    fName  = "MMcEvtBasic";
     68    fTitle = "Basic event info from Monte Carlo";
    6069
    61   Clear();
     70    Clear();
    6271}
    6372
    6473// --------------------------------------------------------------------------
    65 //  constuctor II
    6674//
    67 //  All datamembers are parameters.
     75// Constructor. Use this to set all data members
    6876//
    69 MMcEvtBasic::MMcEvtBasic(MMcEvt::ParticleId_t usPId,
    70                          Float_t              fEner,
    71                          Float_t              fImpa,
    72                          Float_t              fTPhii,
    73                          Float_t             fTThet)
     77// THIS FUNCTION IS FOR THE SIMULATION OLNY.
     78// DON'T USE THIS MEMBERFUNCTION IN THE ANALYSIS.
     79//
     80MMcEvtBasic::MMcEvtBasic(ParticleId_t usPId, Float_t fEner,
     81                         Float_t fImpa, Float_t fTPhii, Float_t fTThet)
    7482{
    75   fName  = "MMcEvtBasic";
    76   fTitle = "Basic event info from Monte Carlo";
     83    fName  = "MMcEvtBasic";
     84    fTitle = "Basic event info from Monte Carlo";
    7785
    78   fPartId         = usPId;
    79   fEnergy         = fEner;
    80   fImpact         = fImpa;
    81   fTelescopePhi   = fTPhii;
    82   fTelescopeTheta = fTThet;
     86    Fill(usPId, fEner, fImpa, fTPhii, fTThet);
    8387}
    84 
    85 
    86 // --------------------------------------------------------------------------
    87 //  default destructor
    88 //
    89 MMcEvtBasic::~MMcEvtBasic() {
    90 }
    91 
    9288
    9389// --------------------------------------------------------------------------
    9490//
    95 //  Reset all values. We have no "error" value for fPartId,
    96 //  we just set kGAMMA
     91//  Reset all values: Fill(kUNDEFINED, -1, -1, 0, 0)
    9792//
    9893void MMcEvtBasic::Clear(Option_t *opt)
    9994{
    100   fPartId         = MMcEvt::kGAMMA;
    101   fEnergy         = -1.;
    102   fImpact         = -1.;
    103   fTelescopeTheta = -1.;
    104   fTelescopePhi   = -1.;
     95    Fill(kUNDEFINED, -1, -1, 0, 0);
    10596}
    10697
    … …  
    109100// Fill all data members
    110101//
    111 void MMcEvtBasic::Fill(MMcEvt::ParticleId_t usPId,
    112                        Float_t  fEner,
    113                        Float_t  fImpa,
    114                        Float_t  fTPhii,
    115                        Float_t  fTThet)
     102void MMcEvtBasic::Fill(ParticleId_t usPId, Float_t fEner,
     103                       Float_t fImpa, Float_t fTPhii, Float_t fTThet)
    116104{
    117   fPartId = usPId;
    118   fEnergy = fEner;
    119   fImpact = fImpa;
    120   fTelescopePhi = fTPhii;
    121   fTelescopeTheta = fTThet;
     105    fPartId         = usPId;
     106
     107    fEnergy         = fEner;
     108    fImpact         = fImpa;
     109
     110    fTelescopePhi   = fTPhii;
     111    fTelescopeTheta = fTThet;
    122112}
    123113
    … …  
    131121void MMcEvtBasic::Print(Option_t *opt) const
    132122{
    133   //
    134   //  print out the data member on screen
    135   //
    136   TString str(opt);
    137   if (str.IsNull())
     123    //
     124    //  print out the data member on screen
     125    //
     126    TString str(opt);
     127    if (str.IsNull())
    138128    {
    139       *fLog << all << endl;
    140       *fLog << "Monte Carlo output:" << endl;
    141       *fLog << " Particle Id:    ";
    142       switch(fPartId)
    143         {
    144         case MMcEvt::kGAMMA:
    145           *fLog << "Gamma" << endl;
    146           break;
    147         case MMcEvt::kPROTON:
    148           *fLog << "Proton" << endl;
    149           break;
    150         case MMcEvt::kHELIUM:
    151           *fLog << "Helium" << endl;
    152           break;
    153         default:
    154           break;
    155         }
    156       *fLog << " Energy:         " << fEnergy << "GeV" << endl;
    157       *fLog << " Impactpar.:     " << fImpact/100 << "m" << endl;
    158       *fLog << endl;
    159       return;
     129        *fLog << all << endl;
     130        *fLog << "Monte Carlo output:" << endl;
     131        *fLog << " Particle Id:    " << GetParticleName() << endl;
     132        *fLog << " Energy:         " << fEnergy << "GeV" << endl;
     133        *fLog << " Impactparam.:   " << fImpact/100 << "m" << endl;
     134        *fLog << endl;
     135        return;
    160136    }
    161   if (str.Contains("id", TString::kIgnoreCase))
    162     switch(fPartId)
    163       {
    164       case MMcEvt::kGAMMA:
    165         *fLog << "Particle: Gamma" << endl;
    166         break;
    167       case MMcEvt::kPROTON:
    168         *fLog << "Particle: Proton" << endl;
    169         break;
    170       case MMcEvt::kHELIUM:
    171         *fLog << "Particle: Helium" << endl;
    172         break;
    173       default:
    174         break;
    175       }
    176   if (str.Contains("energy", TString::kIgnoreCase))
    177     *fLog << "Energy: " << fEnergy << "GeV" << endl;
    178   if (str.Contains("impact", TString::kIgnoreCase))
    179     *fLog << "Impact: " << fImpact << "cm" << endl;
     137    if (str.Contains("id", TString::kIgnoreCase))
     138        *fLog << "Particle: " << GetParticleName() << endl;
     139    if (str.Contains("energy", TString::kIgnoreCase))
     140        *fLog << "Energy: " << fEnergy << "GeV" << endl;
     141    if (str.Contains("impact", TString::kIgnoreCase))
     142        *fLog << "Impact: " << fImpact << "cm" << endl;
    180143}

  • trunk/MagicSoft/include-Classes/MMcFormat/MMcEvtBasic.h

    r6375 r7094  
    1 #ifndef __MMcEvtBasic__
    2 #define __MMcEvtBasic__
     1#ifndef MARS_MMcEvtBasic
     2#define MARS_MMcEvtBasic
    33
    44#ifndef MARS_MParContainer
    … …  
    66#endif
    77
    8 #include "MMcEvt.hxx"
    9 
    10 //
    11 // Version 1:
    12 // New container to keep the very basic informations on the
    13 // original MC events produced by Corsika
    14 //
    158
    169class MMcEvtBasic : public MParContainer
    1710{
     11public:
     12    enum ParticleId_t
     13    {
     14        kUNDEFINED = -1,
     15        kGAMMA     =  1,
     16        kPOSITRON  =  2,
     17        kELECTRON  =  3,
     18        kANTIMUON  =  5,
     19        kMUON      =  6,
     20        kPI0       =  7,
     21        kNEUTRON   = 13,
     22        kPROTON  =   14,
     23        kHELIUM  =  402,
     24        kOXYGEN  = 1608,
     25        kIRON    = 5626
     26    };
    1827
    1928private:
    20   MMcEvt::ParticleId_t fPartId;  // Type of particle
    21   Float_t              fEnergy;  // [GeV] Energy
    22   Float_t              fImpact;  // [cm] impact parameter
     29  ParticleId_t fPartId;  // Type of particle
     30  Float_t      fEnergy;  // [GeV] Energy
     31  Float_t      fImpact;  // [cm] impact parameter
    2332
    2433  // Telescope orientation (see TDAS 02-11 regarding the
    2534  // precise meaning of these angles):
    26   Float_t              fTelescopePhi;    // [rad]
    27   Float_t              fTelescopeTheta;  // [rad]
    28 
     35  Float_t      fTelescopePhi;    // [rad]
     36  Float_t      fTelescopeTheta;  // [rad]
    2937 
    3038 public:
    3139  MMcEvtBasic();
    32  
    33   MMcEvtBasic(MMcEvt::ParticleId_t, Float_t, Float_t, Float_t, Float_t);
    34   ~MMcEvtBasic();
     40  MMcEvtBasic(ParticleId_t, Float_t, Float_t, Float_t, Float_t);
    3541
    36   void Clear(Option_t *opt=NULL);
     42  // Getter
     43  ParticleId_t GetPartId() const { return fPartId; }
    3744
    38   void Fill(MMcEvt::ParticleId_t, Float_t, Float_t, Float_t, Float_t);
    39 
    40   void Print(Option_t *opt=NULL) const;
    41 
    42   MMcEvt::ParticleId_t GetPartId() const { return fPartId; }
    4345  Float_t GetEnergy()  const { return fEnergy; }
    4446  Float_t GetImpact()  const { return fImpact; }
     47
    4548  Float_t GetTelescopePhi() const { return fTelescopePhi; }
    4649  Float_t GetTelescopeTheta() const { return fTelescopeTheta; }
    … …  
    5053      switch (fPartId)
    5154      {
    52       case MMcEvt::kGAMMA:    return "Gamma";
    53       case MMcEvt::kPOSITRON: return "Positron";
    54       case MMcEvt::kELECTRON: return "Electron";
    55       case MMcEvt::kANTIMUON: return "Anti-Muon";
    56       case MMcEvt::kMUON:     return "Muon";
    57       case MMcEvt::kPI0:      return "Pi-0";
    58       case MMcEvt::kNEUTRON:  return "Neutron";
    59       case MMcEvt::kPROTON:   return "Proton";
    60       case MMcEvt::kHELIUM:   return "Helium";
    61       case MMcEvt::kOXYGEN:   return "Oxygen";
    62       case MMcEvt::kIRON:     return "Iron";
     55      case kUNDEFINED:return "Undefined";
     56      case kGAMMA:    return "Gamma";
     57      case kPOSITRON: return "Positron";
     58      case kELECTRON: return "Electron";
     59      case kANTIMUON: return "Anti-Muon";
     60      case kMUON:     return "Muon";
     61      case kPI0:      return "Pi-0";
     62      case kNEUTRON:  return "Neutron";
     63      case kPROTON:   return "Proton";
     64      case kHELIUM:   return "Helium";
     65      case kOXYGEN:   return "Oxygen";
     66      case kIRON:     return "Iron";
    6367      }
    6468
    … …  
    7074      switch (fPartId)
    7175      {
    72       case MMcEvt::kGAMMA:    return "\\gamma";
    73       case MMcEvt::kPOSITRON: return "e^{+}";
    74       case MMcEvt::kELECTRON: return "e^{-}";
    75       case MMcEvt::kANTIMUON: return "\\mu^{+}";
    76       case MMcEvt::kMUON:     return "\\mu^{-}";
    77       case MMcEvt::kPI0:      return "\\pi^{0}";
    78       case MMcEvt::kNEUTRON:  return "n";
    79       case MMcEvt::kPROTON:   return "p";
    80       case MMcEvt::kHELIUM:   return "He";
    81       case MMcEvt::kOXYGEN:   return "O";
    82       case MMcEvt::kIRON:     return "Fe";
     76      case kUNDEFINED:return "N/A";
     77      case kGAMMA:    return "\\gamma";
     78      case kPOSITRON: return "e^{+}";
     79      case kELECTRON: return "e^{-}";
     80      case kANTIMUON: return "\\mu^{+}";
     81      case kMUON:     return "\\mu^{-}";
     82      case kPI0:      return "\\pi^{0}";
     83      case kNEUTRON:  return "n";
     84      case kPROTON:   return "p";
     85      case kHELIUM:   return "He";
     86      case kOXYGEN:   return "O";
     87      case kIRON:     return "Fe";
    8388      }
    8489
    … …  
    101106
    102107
    103   void SetPartId(MMcEvt::ParticleId_t id)
    104     { fPartId = id; }
    105 
    106   void SetEnergy(Float_t Energy)
    107   { fEnergy=Energy; }              //Set Energy
    108  
    109   void SetImpact(Float_t Impact)
    110   { fImpact=Impact;}               //Set impact parameter
     108  // Setter
     109  void SetPartId(ParticleId_t id) { fPartId = id; }
     110  void SetEnergy(Float_t Energy)  { fEnergy=Energy; }              //Set Energy
     111  void SetImpact(Float_t Impact)  { fImpact=Impact;}               //Set impact parameter
    111112
    112113  void SetTelescopeTheta(Float_t Theta) { fTelescopeTheta=Theta; }
    113 
    114114  void SetTelescopePhi  (Float_t Phi)   { fTelescopePhi=Phi; }
    115115
    116   ClassDef(MMcEvtBasic, 1) //Stores Basic Montecarlo Information of one event
     116  void Fill(ParticleId_t, Float_t, Float_t, Float_t, Float_t);
     117
     118  // TObject
     119  void Clear(Option_t *opt=NULL);
     120  void Print(Option_t *opt=NULL) const;
     121
     122  ClassDef(MMcEvtBasic, 2) //Stores Basic Montecarlo Information of one event
    117123
    118124};
Note: See TracChangeset for help on using the changeset viewer.