Changeset 1415 for trunk/MagicSoft/Mars


Ignore:
Timestamp:
07/16/02 15:35:15 (22 years ago)
Author:
wittek
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mhist
Files:
2 added
5 edited

Legend:

Unmodified
Added
Removed

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

    r1336 r1415  
    1818#pragma link C++ class MHEnergyTime+;
    1919#pragma link C++ class MHEnergyTheta+;
    20 #pragma link C++ class MHAlphaEnergyTime;
    21 #pragma link C++ class MHAlphaEnergyTheta;
    22 #pragma link C++ class MHEffOnTimeTime+;
    23 #pragma link C++ class MHEffOnTimeTheta+;
     20#pragma link C++ class MHAlphaEnergyTime+;
     21#pragma link C++ class MHAlphaEnergyTheta+;
     22#pragma link C++ class MHGamma+;
     23#pragma link C++ class MHEffOnTime+;
     24
    2425#pragma link C++ class MHTimeDiffTime+;
    2526#pragma link C++ class MHTimeDiffTheta+;
    … …  
    4041#pragma link C++ class MHMcCollectionArea+;
    4142#pragma link C++ class MHMcEnergyMigration+;
     43#pragma link C++ class MHFlux;
    4244
    43 #pragma link C++ class MHCompProb+;
    44 #pragma link C++ class MHHadroness+;
    4545
    4646#endif

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

    r1330 r1415  
    6666    fHist.SetTitle("3D-plot of alpha, E-est, Theta");
    6767    fHist.SetXTitle("\\alpha [\\circ]");
    68     fHist.SetYTitle("E_{est} [GeV]");
     68    fHist.SetYTitle("E-est [GeV]            ");
    6969    fHist.SetZTitle("\\Theta [\\circ]");
    7070}
    … …  
    7979   if (!fEnergy)
    8080   {
    81        *fLog << err << dbginf << "MEnergyEst not found... aborting." << endl;
     81       *fLog << err << dbginf << "MHAlphaEnergyTheta : MEnergyEst not found... aborting." << endl;
    8282       return kFALSE;
    8383   }
    … …  
    8686   if (!fMcEvt)
    8787   {
    88        *fLog << err << dbginf << "MMcEvt not found... aborting." << endl;
     88       *fLog << err << dbginf << "MHAlphaEnergyTheta : MMcEvt not found... aborting." << endl;
    8989       return kFALSE;
    9090   }
    9191
    9292   MBinning* binsenergy = (MBinning*)plist->FindObject("BinningE");
    93    MBinning* binsalpha  = (MBinning*)plist->FindObject("BinningAlpha");
     93   MBinning* binsalphaflux  = (MBinning*)plist->FindObject("BinningAlphaFlux");
    9494   MBinning* binstheta  = (MBinning*)plist->FindObject("BinningTheta");
    95    if (!binsenergy || !binsalpha || !binstheta)
     95   if (!binsenergy || !binsalphaflux || !binstheta)
    9696   {
    97        *fLog << err << dbginf << "At least one MBinning not found... aborting." << endl;
     97       *fLog << err << dbginf << "MHAlphaEnergyTheta : At least one MBinning not found... aborting." << endl;
    9898       return kFALSE;     
    9999   }
    100100
    101    SetBinning(&fHist, binsalpha, binsenergy, binstheta);
     101   SetBinning(&fHist, binsalphaflux, binsenergy, binstheta);
    102102
    103103   fHist.Sumw2();
    104104
     105   fHist.Sumw2();
     106
    105107   return kTRUE;
    106108}
    … …  
    114116    MHillasSrc &hil = *(MHillasSrc*)par;
    115117
    116     fHist.Fill(hil.GetAlpha(), fEnergy->GetEnergy(), fMcEvt->GetTheta()*kRad2Deg);
     118    fHist.Fill(fabs(hil.GetAlpha()), fEnergy->GetEnergy(), fMcEvt->GetTheta()*kRad2Deg);
    117119
    118120    return kTRUE;
    … …  
    146148
    147149    h->SetTitle("Distribution of E-est [GeV]");
    148     h->SetXTitle("E_{est} [GeV]");
     150    h->SetXTitle("E-est [GeV]            ");
    149151    h->SetYTitle("Counts");
    150152
    … …  
    223225}
    224226
    225 // --------------------------------------------------------------------------
    226 //
    227 // Calculate the histogram as the difference of two histograms :
    228 //          fHist(gamma) = h1(source) - h2(antisource)
    229 //
    230 void MHAlphaEnergyTheta::Subtract(const TH3D *h1, const TH3D *h2)
    231 {
    232     // MH::SetBinning(&fHist, (TH1*)h1);
    233 
    234     //    fHist.Sumw2();
    235     fHist.Add((TH1*)h1, (TH1*)h2, 1, -1); // Root: FIXME
    236 }
    237 
    238 
    239 // --------------------------------------------------------------------------
    240 //
    241 // Integrate fHist(gamma) in the alpha range (lo, up)
    242 //
    243 TH2D *MHAlphaEnergyTheta::GetAlphaProjection(Axis_t lo, Axis_t up)
    244 {
    245     if (up < lo)
    246     {
    247         *fLog << err << fName << ": Alpha projection not possible: lo=" << lo << " up=" << up << endl;
    248         return NULL;
    249     }
    250 
    251     TAxis &axe = *fHist.GetXaxis();
    252 
    253     Int_t ilo = axe.FindFixBin(lo);
    254     Int_t iup = axe.FindFixBin(up);
    255 
    256     const Double_t epslo1 = lo-axe.GetBinLowEdge(ilo);
    257     const Double_t epslo2 = axe.GetBinUpEdge(ilo)-lo;
    258 
    259     const Double_t epsup1 = up-axe.GetBinLowEdge(iup);
    260     const Double_t epsup2 = axe.GetBinUpEdge(iup)-up;
    261 
    262     const Double_t epslo = epslo1<epslo2 ? epslo1 : epslo2;
    263     const Double_t epsup = epsup1<epsup2 ? epsup1 : epsup2;
    264 
    265     if (epslo1>epslo2)
    266         ilo++;
    267 
    268     if (epsup1<epsup2)
    269         iup--;
    270 
    271     if (epslo>0.01*axe.GetBinWidth(ilo) || epsup>0.01*axe.GetBinWidth(iup))
    272     {
    273         *fLog << err << fName << ": binning is not adequate for the requested projection:" << endl;
    274         *fLog << "Please specify a lower or upper limit which is not more than 1% away from a bin edge" << endl;
    275         *fLog << " epslo = " << epslo << endl;
    276         *fLog << " epsup = " << epsup << endl;
    277         *fLog << " dwl   = " << axe.GetBinWidth(ilo) << endl;
    278         *fLog << " dwu   = " << axe.GetBinWidth(iup) << endl;
    279         return NULL;
    280     }
    281 
    282     axe.SetRange(ilo, iup);
    283 
    284     TH2D &h2D = *(TH2D *)fHist.Project3D("ezypro");
    285 
    286     h2D.SetTitle("2D-plot  of Theta vs. E-est");
    287     h2D.SetXTitle("E-est [GeV]            ");
    288     h2D.SetYTitle("\\Theta [\\circ]");
    289 
    290     return &h2D;
    291 }
    292 
    293 // --------------------------------------------------------------------------
    294 //
    295 // Draw the integrated histogram
    296 //
    297 TH2D *MHAlphaEnergyTheta::DrawAlphaProjection(Axis_t lo, Axis_t up, Option_t *opt)
    298 {
    299     TH2D *h2D = GetAlphaProjection(lo, up);
    300 
    301     if (!h2D)
    302         return NULL;
    303 
    304     char txt[100];
    305     sprintf(txt, "No.of Gammas vs. E-est and Theta (%.1f < alpha < %.1f deg)", lo, up);
    306 
    307     //    TCanvas *c = MakeDefCanvas("AlphaEnergyTheta", "2D histogram of gamma signal in energy and theta");
    308     TCanvas &c = *MakeDefCanvas("AlphaEnergyTheta", txt);
    309 
    310     c.Divide(2, 2);
    311 
    312     gROOT->SetSelectedPad(NULL);
    313 
    314     TH1 *h;
    315 
    316     c.cd(1);
    317     h = h2D->ProjectionX("Eest", -1, 9999, "E");
    318     h->SetTitle("Distribution of E-est [GeV]");
    319     h->SetXTitle("E-est [GeV]            ");
    320     h->SetYTitle("Counts");
    321 
    322     h->Draw(opt);
    323     h->SetBit(kCanDelete);
    324     gPad->SetLogx();
    325 
    326     c.cd(2);
    327     h = h2D->ProjectionY("theta", -1, 9999, "E");
    328     h->SetTitle("Distribution of \\Theta [\\circ]");
    329     h->SetXTitle("\\Theta [\\circ]");
    330     h->SetYTitle("Counts");
    331 
    332     h->Draw(opt);
    333     h->SetBit(kCanDelete);
    334 
    335     c.cd(3);
    336 
    337     h2D->DrawCopy(opt);
    338     gPad->SetLogx();
    339 
    340     c.Modified();
    341     c.Update();
    342 
    343     return h2D;
    344 }
    345 
    346 
    347 
    348 
    349 
    350 
    351 
    352 
    353 
    354 
    355 
    356 

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

    r1292 r1415  
    4242    TObject *DrawClone(Option_t *option="") const;
    4343
    44     void Subtract(const TH3D *h1, const TH3D *h2);
    45     void Subtract(const MHAlphaEnergyTheta *h1, const MHAlphaEnergyTheta *h2)
    46     {
    47         Subtract(h1->GetHist(), h2->GetHist());
    48     }
    49 
    50     TH2D *DrawAlphaProjection(Axis_t lo, Axis_t up, Option_t *opt="");
    51     TH2D *GetAlphaProjection(Axis_t lo, Axis_t up);
    52 
    5344    ClassDef(MHAlphaEnergyTheta, 1) //3D-histogram in alpha, Energy and theta
    5445};

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

    r1330 r1415  
    6868    fHist.SetTitle("3D-plot of alpha, E-est, time");
    6969    fHist.SetXTitle("\\alpha [\\circ]");
    70     fHist.SetYTitle("E_{est} [GeV]");
     70    fHist.SetYTitle("E-est [GeV]            ");
    7171    fHist.SetZTitle("time [s]");
    7272}
    … …  
    8181   if (!fEnergy)
    8282   {
    83        *fLog << err << dbginf << "MEnergyEst not found... aborting." << endl;
     83       *fLog << err << dbginf << "MHAlphaEnergyTime : MEnergyEst not found... aborting." << endl;
    8484       return kFALSE;
    8585   }
    … …  
    8888   if (!fTime)
    8989   {
    90        *fLog << err << dbginf << "MTime not found... aborting." << endl;
     90       *fLog << err << dbginf << "MHAlphaEnergyTime : MTime not found... aborting." << endl;
    9191       return kFALSE;
    9292   }
    9393
    9494   MBinning* binsenergy = (MBinning*)plist->FindObject("BinningE");
    95    MBinning* binsalpha  = (MBinning*)plist->FindObject("BinningAlpha");
     95   MBinning* binsalphaflux  = (MBinning*)plist->FindObject("BinningAlphaFlux");
    9696   MBinning* binstime   = (MBinning*)plist->FindObject("BinningTime");
    97    if (!binsenergy || !binsalpha || !binstime)
     97   if (!binsenergy || !binsalphaflux || !binstime)
    9898   {
    99        *fLog << err << dbginf << "At least one MBinning not found... aborting." << endl;
     99       *fLog << err << dbginf << "MHAlphaEnergyTime : At least one MBinning not found... aborting." << endl;
    100100       return kFALSE;     
    101101   }
    102102
    103    SetBinning(&fHist, binsalpha, binsenergy, binstime);
     103   SetBinning(&fHist, binsalphaflux, binsenergy, binstime);
    104104
    105105   fHist.Sumw2();
    … …  
    116116    MHillasSrc &hil = *(MHillasSrc*)par;
    117117
    118     fHist.Fill(hil.GetAlpha(), fEnergy->GetEnergy(), 0.0001*fTime->GetTimeLo());
     118    fHist.Fill(fabs(hil.GetAlpha()), fEnergy->GetEnergy(), 0.0001*fTime->GetTimeLo());
    119119    return kTRUE;
    120120}
    … …  
    147147
    148148    h->SetTitle("Distribution of E-est [GeV]");
    149     h->SetXTitle("E_{est} [GeV]");
     149    h->SetXTitle("E-est [GeV]            ");
    150150    h->SetYTitle("Counts");
    151151
    … …  
    226226}
    227227
    228 // --------------------------------------------------------------------------
    229 //
    230 // Calculate the histogram as the difference of two histograms :
    231 //          fHist(gamma) = h1(source) - h2(antisource)
    232 //
    233 void MHAlphaEnergyTime::Subtract(const TH3D *h1, const TH3D *h2)
    234 {
    235   //    MH::SetBinning(&fHist, (TH1*)h1);
    236 
    237   //    fHist.Sumw2();
    238     fHist.Add((TH1*)h1, (TH1*)h2, 1, -1); // ROOT: FIXME!
    239 }
    240 
    241 // --------------------------------------------------------------------------
    242 //
    243 // Integrate fHist(gamma) in the alpha range (lo, up)
    244 //
    245 TH2D *MHAlphaEnergyTime::GetAlphaProjection(Axis_t lo, Axis_t up)
    246 {
    247     if (up < lo)
     228
     229// --------------------------------------------------------------------------
     230//
     231// Integrate fHist     (Alpha,E-est,Time) over the Time to get
     232//           fAlphaEest(Alpha,E-est)
     233//
     234TH2D *MHAlphaEnergyTime::IntegrateTime(const char *title, Bool_t draw)
     235{
     236    Int_t nzbins = fHist.GetNbinsZ();
     237    TAxis &axez  = *fHist.GetZaxis();
     238    axez.SetRange(1,nzbins);
     239
     240    TH2D &fAlphaEest = *(TH2D *)fHist.Project3D("exy");
     241
     242    fAlphaEest.SetTitle(title);
     243    fAlphaEest.SetXTitle("E-est [GeV]            ");
     244    fAlphaEest.SetYTitle("\\alpha  [  \\circ]");
     245
     246    if (draw == kTRUE)
    248247    {
    249         *fLog << err << fName << ": Alpha projection not possible: lo=" << lo << " up=" << up << endl;
    250         return NULL;
     248      TCanvas &c = *MakeDefCanvas(title, title);
     249
     250      gROOT->SetSelectedPad(NULL);
     251
     252      fAlphaEest.DrawCopy();
     253      gPad->SetLogx();
     254
     255      c.Modified();
     256      c.Update();
    251257    }
    252258
    253     TAxis &axe = *fHist.GetXaxis();
    254 
    255     Int_t ilo = axe.FindFixBin(lo);
    256     Int_t iup = axe.FindFixBin(up);
    257 
    258     const Double_t epslo1 = lo-axe.GetBinLowEdge(ilo);
    259     const Double_t epslo2 = axe.GetBinUpEdge(ilo)-lo;
    260 
    261     const Double_t epsup1 = up-axe.GetBinLowEdge(iup);
    262     const Double_t epsup2 = axe.GetBinUpEdge(iup)-up;
    263 
    264     const Double_t epslo = epslo1<epslo2 ? epslo1 : epslo2;
    265     const Double_t epsup = epsup1<epsup2 ? epsup1 : epsup2;
    266 
    267     if (epslo1>epslo2)
    268         ilo++;
    269 
    270     if (epsup1<epsup2)
    271         iup--;
    272 
    273     if (epslo>0.01*axe.GetBinWidth(ilo) || epsup>0.01*axe.GetBinWidth(iup))
     259    return &fAlphaEest;
     260}
     261
     262// --------------------------------------------------------------------------
     263//
     264// Integrate fHist     (Alpha,E-est,Time) over E-est to get
     265//           fAlphaTime(Alpha,Time)
     266//
     267TH2D *MHAlphaEnergyTime::IntegrateEest(const char *title, Bool_t draw)
     268{
     269    Int_t nybins = fHist.GetNbinsY();
     270    TAxis &axey  = *fHist.GetYaxis();
     271    axey.SetRange(1,nybins);
     272
     273    TH2D &fAlphaTime = *(TH2D *)fHist.Project3D("exz");
     274
     275    fAlphaTime.SetTitle(title);
     276    fAlphaTime.SetXTitle("Time [s]");
     277    fAlphaTime.SetYTitle("\\alpha  [  \\circ]");
     278
     279    if (draw == kTRUE)
    274280    {
    275         *fLog << err << fName << ": binning is not adequate for the requested projection:" << endl;
    276         *fLog << "Please specify a lower or upper limit which is not more than 1% away from a bin edge" << endl;
    277         *fLog << " epslo = " << epslo << endl;
    278         *fLog << " epsup = " << epsup << endl;
    279         *fLog << " dwl   = " << axe.GetBinWidth(ilo) << endl;
    280         *fLog << " dwu   = " << axe.GetBinWidth(iup) << endl;
    281         return NULL;
     281      TCanvas &c = *MakeDefCanvas(title, title);
     282
     283      gROOT->SetSelectedPad(NULL);
     284
     285      fAlphaTime.DrawCopy();
     286
     287      c.Modified();
     288      c.Update();
    282289    }
    283290
    284     axe.SetRange(ilo, iup);
    285 
    286     TH2D &h2D = *(TH2D *)fHist.Project3D("ezy");
    287 
    288     h2D.SetTitle("2D-plot  of time vs. E-est");
    289     h2D.SetXTitle("E-est [GeV]            ");
    290     h2D.SetYTitle("time [s]");
    291 
    292     return &h2D;
    293 }
    294 
    295 //---------------------------------------------------------
    296 //
    297 // Draw the projected histogram
    298 //
    299 TH2D *MHAlphaEnergyTime::DrawAlphaProjection(Axis_t lo, Axis_t up, Option_t *opt)
    300 {
    301     TH2D *h2D = GetAlphaProjection(lo, up);
    302 
    303     if (!h2D)
    304         return NULL;
    305 
    306     char txt[100];
    307     sprintf(txt, "No.of Gammas vs. E-est and Time (%.1f < alpha < %.1f deg)", lo, up);
    308 
    309     //    TCanvas *c = MakeDefCanvas("AlphaEnergyTime", "2D histogram of gamma signal in energy and time");
    310     TCanvas &c = *MakeDefCanvas("AlphaEnergyTime", txt);
    311 
    312     c.Divide(2, 2);
    313 
    314     gROOT->SetSelectedPad(NULL);
    315 
    316     TH1 *h;
    317 
    318     c.cd(1);
    319     h = h2D->ProjectionX("Eest", -1, 9999, "E");
    320     h->SetTitle("Distribution of E-est [GeV]");
    321     h->SetXTitle("E-est [GeV]            ");
    322     h->SetYTitle("Counts");
    323 
    324     h->Draw(opt);
    325     h->SetBit(kCanDelete);
    326     gPad->SetLogx();
    327 
    328     c.cd(2);
    329     h = h2D->ProjectionY("time", -1, 9999, "E");
    330     h->SetTitle("Distribution of time [s]");
    331     h->SetXTitle("time [s]");
    332     h->SetYTitle("Counts");
    333 
    334     h->Draw(opt);
    335     h->SetBit(kCanDelete);
    336 
    337     c.cd(3);
    338 
    339     h2D->DrawCopy(opt);
    340     gPad->SetLogx();
    341 
    342     c.Modified();
    343     c.Update();
    344 
    345     return h2D;
    346 }
    347 
     291    return &fAlphaTime;
     292}
     293
     294// --------------------------------------------------------------------------
     295//
     296// Integrate fHist (Alpha,E-est,Time) over Eest and Time to get
     297//           fAlpha(Alpha)
     298//
     299TH1D *MHAlphaEnergyTime::IntegrateEestTime(const char *title, Bool_t draw)
     300{
     301    Int_t nybins = fHist.GetNbinsY();
     302    TAxis &axey  = *fHist.GetYaxis();
     303    axey.SetRange(1,nybins);
     304
     305    Int_t nzbins = fHist.GetNbinsZ();
     306    TAxis &axez  = *fHist.GetZaxis();
     307    axez.SetRange(1,nzbins);
     308
     309    TH1D &fAlpha = *(TH1D *)fHist.Project3D("ex");
     310
     311    fAlpha.SetTitle(title);
     312    fAlpha.SetXTitle("\\alpha  [  \\circ]");
     313    fAlpha.SetYTitle("Counts");
     314
     315    if (draw == kTRUE)
     316    {
     317      TCanvas &c = *MakeDefCanvas(title, title);
     318
     319      gROOT->SetSelectedPad(NULL);
     320
     321      fAlpha.DrawCopy();
     322
     323      c.Modified();
     324      c.Update();
     325    }
     326
     327    return &fAlpha;
     328}
     329
     330
     331
     332
     333
     334
     335

  • trunk/MagicSoft/Mars/mhist/Makefile

    r1336 r1415  
    4444           MHAlphaEnergyTime.cc \
    4545           MHAlphaEnergyTheta.cc \
    46            MHEffOnTimeTime.cc \
    47            MHEffOnTimeTheta.cc \
     46           MHEffOnTime.cc \
    4847           MHTimeDiffTime.cc \
    4948           MHTimeDiffTheta.cc \
    … …  
    5554           MHMcEfficiencyEnergy.cc \
    5655           MHMcEnergyImpact.cc \
     56           MHMcRate.cc \
     57           MHThetabarTime.cc \
     58           MHThetabarTheta.cc \
    5759           MHMcEnergyMigration.cc \
    58            MHThetabarTime.cc \
    59            MHMcRate.cc \
    60            MHMcIntRate.cc \
    61            MHMcDifRate.cc
    62 
     60           MHGamma.cc \
     61           MHFlux.cc
     62 
    6363
    6464SRCS    = $(SRCFILES)
    … …  
    7777
    7878# @endcode
     79
     80
     81
     82
     83
     84
     85
     86
     87
     88
Note: See TracChangeset for help on using the changeset viewer.