Ignore:
Timestamp:
03/13/03 22:25:00 (22 years ago)
Author:
moralejo
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mhist
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/mhist/MHMcCT1CollectionArea.cc

    r1821 r1823  
    5252    //   we set the energy range from 10 Gev to 30000 GeV (in log 4.5 orders
    5353    //   of magnitude) and for each order we take 10 subdivision --> 45 xbins
    54     //
    55     //   we set the radius range from 0 m to 500 m with 10 m bin --> 50 ybins
     54    //   we set the theta range from 12.5 to 48 deg, with 6 bins.
    5655    //
    5756    fName  = name  ? name  : "MHMcCT1CollectionArea";
     
    6059    MBinning binsx;
    6160    MBinning binsy;
     61
    6262    binsx.SetEdgesLog(45, 10, 30000);
    63     binsy.SetEdges   (50,  0,   500);
     63    const Double_t yedge[7] = {12.5, 17.5, 23.5, 29.5, 35.5, 42., 48.};
     64    const TArrayD yed(7,yedge);
     65    binsy.SetEdges(yed);
    6466
    6567    fHistAll = new TH2D;
    6668    fHistSel = new TH2D;
    67     fHistCol = new TH1D;
    68  
     69    fHistCol = new TH2D;
     70
    6971    MH::SetBinning(fHistAll, &binsx, &binsy);
    7072    MH::SetBinning(fHistSel, &binsx, &binsy);
     73    MH::SetBinning(fHistCol, &binsx, &binsy);
     74
    7175
    7276    fHistCol->SetName(fName);
     
    7579
    7680    fHistCol->SetTitle(fTitle);
    77     fHistAll->SetTitle("All showers - Radius vs Energy distribution");
    78     fHistSel->SetTitle("Selected showers - Radius vs Energy distribution");
     81    fHistAll->SetTitle("All showers - Theta vs Energy distribution");
     82    fHistSel->SetTitle("Selected showers - Theta vs Energy distribution");
    7983
    8084    fHistAll->SetDirectory(NULL);
     
    8387
    8488    fHistAll->SetXTitle("E [GeV]");
    85     fHistAll->SetYTitle("r [m]");
     89    fHistAll->SetYTitle("theta [deg]");
    8690    fHistAll->SetZTitle("N");
    8791
    8892    fHistSel->SetXTitle("E [GeV]");
    89     fHistSel->SetYTitle("r [m]");
     93    fHistSel->SetYTitle("theta [deg]");
    9094    fHistSel->SetYTitle("N");
    9195
    9296    fHistCol->SetXTitle("E [GeV]");
    93     fHistCol->SetYTitle("A [m^{2}]");
     97    fHistCol->SetYTitle("theta [deg]");
     98    fHistCol->SetZTitle("A [m^{2}]");
    9499}
    95100
     
    109114// Fill data into the histogram which contains the selected showers
    110115//
    111 void MHMcCT1CollectionArea::FillSel(Double_t energy, Double_t radius)
    112 {
    113     fHistSel->Fill(energy, radius);
     116void MHMcCT1CollectionArea::FillSel(Double_t energy, Double_t theta)
     117{
     118    fHistSel->Fill(energy, theta);
    114119}
    115120
     
    159164
    160165    //
    161     // This is necessary to get the expected bahviour of DrawClone
     166    // This is necessary to get the expected behaviour of DrawClone
    162167    //
    163168    gROOT->SetSelectedPad(NULL);
     
    190195//  and set the 'ReadyToSave' flag
    191196//
    192 void MHMcCT1CollectionArea::CalcEfficiency(Float_t theta)
     197void MHMcCT1CollectionArea::CalcEfficiency()
    193198{
    194199  //
     
    205210  //
    206211
    207     TH1D &histsel = *fHistSel->ProjectionX();
    208 
    209     TH1D histall;
    210 
    211     TAxis &xaxis = *histsel.GetXaxis();
    212 
    213     MH::SetBinning(&histall, &xaxis);
    214     MH::SetBinning(fHistCol, &xaxis);
    215 
    216     Float_t emin1, emax1, emin2, emax2;
    217     Float_t index, expo, k1, k2;
    218     Float_t numevts1, numevts2;
    219     Float_t r1, r2;        // Impact parameter range (on ground).
    220 
    221     emin1 = 0;  emax1 = 0; emin2 = 0;  emax2 = 0;
    222     expo = 0.; k1 = 0.; k2 = 0.; r1 = 0.; r2 = 0.;
    223     numevts1 = 0; numevts2 = 0;
    224 
    225     if (theta > 0.26 && theta < 0.27)  // 15 degrees
    226       {
    227         r1 = 0.;
    228         r2 = 250.;     //meters
    229         emin1 = 300.;
    230         emax1 = 400.;  // Energies in GeV.
    231         emin2 = 400.;
    232         emax2 = 30000.;
    233         numevts1 = 4000.;
    234         numevts2 = 25740.;
    235       }
    236     else if (theta > 0.34 && theta < 0.35)  // 20 degrees
    237       {
    238         r1 = 0.;
    239         r2 = 263.;     //meters
    240         emin1 = 300.;
    241         emax1 = 400.;  // Energies in GeV.
    242         emin2 = 400.;
    243         emax2 = 30000.;
    244         numevts1 = 6611.;
    245         numevts2 = 24448.;
    246       }
    247     else if (theta > 0.47 && theta < 0.48)  // 27 degrees
    248       {
    249         r1 = 0.;
    250         r2 = 290.;     //meters
    251         emin1 = 300.;
    252         emax1 = 400.;  // Energies in GeV.
    253         emin2 = 400.;
    254         emax2 = 30000.;
    255         numevts1 = 4000.;
    256         numevts2 = 26316.;
    257       }
    258     else if (theta > 0.55 && theta < 0.56)  // 32 degrees
    259       {
    260         r1 = 0.;
    261         r2 = 350.;     //meters
    262         emin1 = 300.;
    263         emax1 = 30000.;  // Energies in GeV.
    264         numevts1 = 33646.;
    265       }
    266     else if (theta > 0.68 && theta < 0.69)  // 39 degrees
    267       {
    268         r1 = 0.;
    269         r2 = 380.;     //meters
    270         emin1 = 300.;
    271         emax1 = 30000.;  // Energies in GeV.
    272         numevts1 = 38415.;
    273       }
    274     else if (theta > 0.78 && theta < 0.79)  // 45 degrees
    275       {
    276         r1 = 0.;
    277         r2 = 565.;     //meters
    278         emin1 = 300.;
    279         emax1 = 50000.;  // Energies in GeV.
    280         numevts1 = 30197.;
    281       }
    282 
    283     index = 1.5;     // Differential spectral Index.
    284     expo = 1.-index;
    285     k1 = numevts1 / (pow(emax1,expo) - pow(emin1,expo));
    286     k2 = numevts2 / (pow(emax2,expo) - pow(emin2,expo));
    287 
    288     const Int_t nbinx = xaxis.GetNbins();
    289 
    290     for (Int_t i=1; i<=nbinx; i++)
     212  for (Int_t thetabin = 1; thetabin <= fHistAll->GetNbinsY(); thetabin++)
    291213    {
    292         const Float_t e1 = histall.GetBinLowEdge(i);
    293         const Float_t e2 = histall.GetBinLowEdge(i+1);
    294 
    295         if (e1 < emin1 || e2 > emax2)
     214      // This theta is not exactly the one of the MC events, just about
     215      // the same:
     216      Float_t theta = fHistAll->GetYaxis()->GetBinCenter(thetabin);
     217
     218      Float_t emin1, emax1, emin2, emax2;
     219      Float_t index, expo, k1, k2;
     220      Float_t numevts1, numevts2;
     221      Float_t r1, r2;        // Impact parameter range (on ground).
     222
     223      emin1 = 0;  emax1 = 0; emin2 = 0;  emax2 = 0;
     224      expo = 0.; k1 = 0.; k2 = 0.; r1 = 0.; r2 = 0.;
     225      numevts1 = 0; numevts2 = 0;
     226
     227      if (theta > 14 && theta < 16)   // 15 deg
     228        {
     229          r1 = 0.;
     230          r2 = 250.;     //meters
     231          emin1 = 300.;
     232          emax1 = 400.;  // Energies in GeV.
     233          emin2 = 400.;
     234          emax2 = 30000.;
     235          numevts1 = 4000.;
     236          numevts2 = 25740.;
     237        }
     238      else if (theta > 20 && theta < 21)  // 20.5 deg
     239        {
     240          r1 = 0.;
     241          r2 = 263.;     //meters
     242          emin1 = 300.;
     243          emax1 = 400.;  // Energies in GeV.
     244          emin2 = 400.;
     245          emax2 = 30000.;
     246          numevts1 = 6611.;
     247          numevts2 = 24448.;
     248        }
     249      else if (theta > 26 && theta < 27)  // 26.5 degrees
     250        {
     251          r1 = 0.;
     252          r2 = 290.;     //meters
     253          emin1 = 300.;
     254          emax1 = 400.;  // Energies in GeV.
     255          emax2 = emax1;          emin2 = 400.;
     256          emax2 = 30000.;
     257          numevts1 = 4000.;
     258          numevts2 = 26316.;
     259        }
     260      else if (theta > 32 && theta < 33)  // 32.5 degrees
     261        {
     262          r1 = 0.;
     263          r2 = 350.;     //meters
     264          emin1 = 300.;
     265          emax1 = 30000.;  // Energies in GeV.
     266          emax2 = emax1;
     267          numevts1 = 33646.;
     268        }
     269      else if (theta > 38 && theta < 39)  // 38.75 degrees
     270        {
     271          r1 = 0.;
     272          r2 = 380.;     //meters
     273          emin1 = 300.;
     274          emax1 = 30000.;  // Energies in GeV.
     275          emax2 = emax1;
     276          numevts1 = 38415.;
     277        }
     278      else if (theta > 44 && theta < 46)  // 45 degrees
     279        {
     280          r1 = 0.;
     281          r2 = 565.;     //meters
     282          emin1 = 300.;
     283          emax1 = 50000.;  // Energies in GeV.
     284          emax2 = emax1;
     285          numevts1 = 30197.;
     286        }
     287
     288      index = 1.5;     // Differential spectral Index.
     289      expo = 1.-index;
     290      k1 = numevts1 / (pow(emax1,expo) - pow(emin1,expo));
     291      k2 = numevts2 / (pow(emax2,expo) - pow(emin2,expo));
     292
     293      for (Int_t i=1; i <= fHistAll->GetNbinsX(); i++)
     294        {
     295          const Float_t e1 = fHistAll->GetXaxis()->GetBinLowEdge(i);
     296          const Float_t e2 = fHistAll->GetXaxis()->GetBinLowEdge(i+1);
     297
     298          if (e1 < emin1 || e2 > emax2)
    296299            continue;
    297300
    298         Float_t events;
    299 
    300         if (e2<emax1)
    301           events = k1 * (pow(e2, expo) - pow(e1, expo));
    302         else if (e1>emin2)
    303           events = k2 * (pow(e2, expo) - pow(e1, expo));
    304         else
    305           events =
    306             k1 * (pow(emax1, expo) - pow(e1, expo))+
    307             k2 * (pow(e2, expo) - pow(emin2, expo));
    308 
    309         histall.SetBinContent(i, events);
    310         histall.SetBinError(i, sqrt(events));
     301          Float_t events;
     302
     303          if (e2 <= emax1)
     304            events = k1 * (pow(e2, expo) - pow(e1, expo));
     305          else if (e1 >= emin2)
     306            events = k2 * (pow(e2, expo) - pow(e1, expo));
     307          else
     308            events =
     309              k1 * (pow(emax1, expo) - pow(e1, expo))+
     310              k2 * (pow(e2, expo) - pow(emin2, expo));
     311
     312          fHistAll->SetBinContent(i, thetabin, events);
     313          fHistAll->SetBinError(i, thetabin, sqrt(events));
     314        }
     315
     316      // -----------------------------------------------------------
     317
     318      const Float_t dr = TMath::Pi() * (r2*r2 - r1*r1);
     319
     320      for (Int_t ix = 1; ix <= fHistAll->GetNbinsX(); ix++)
     321        {
     322          const Float_t Na = fHistAll->GetBinContent(ix,thetabin);
     323
     324          if (Na <= 0)
     325            continue;
     326
     327          const Float_t Ns = fHistSel->GetBinContent(ix,thetabin);
     328
     329          // Since Na is an estimate of the total number of showers generated
     330          // in the energy bin, it may happen that Ns (triggered showers) is
     331          // larger than Na. In that case, the bin is skipped:
     332
     333          if (Na < Ns)
     334            continue;
     335
     336          const Double_t eff = Ns/Na;
     337          const Double_t err = sqrt((1.-eff)*Ns)/Na;
     338
     339          const Float_t area = dr * cos(theta*TMath::Pi()/180.);
     340
     341          fHistCol->SetBinContent(ix, thetabin, eff*area);
     342          fHistCol->SetBinError(ix, thetabin, err*area);
     343        }
    311344    }
    312345
    313     // -----------------------------------------------------------
    314 
    315     const Float_t dr = TMath::Pi() * (r2*r2 - r1*r1);
    316 
    317     for (Int_t ix=1; ix<=nbinx; ix++)
    318       {
    319         const Float_t Na = histall.GetBinContent(ix);
    320 
    321         if (Na <= 0)
    322           continue;
    323 
    324         const Float_t Ns = histsel.GetBinContent(ix);
    325 
    326         // Since Na is an estimate of the total number of showers generated
    327         // in the energy bin, it may happen that Ns (triggered showers) is
    328         // larger than Na. In that case, the bin is skipped:
    329 
    330         if (Na < Ns)
    331           continue;
    332 
    333         const Double_t eff = Ns/Na;
    334 
    335         const Double_t err = sqrt((1.-eff)*Ns)/Na;
    336 
    337         const Float_t area = dr * cos(theta);
    338 
    339         fHistCol->SetBinContent(ix, eff*area);
    340         fHistCol->SetBinError(ix, err*area);
    341       }
    342 
    343     delete &histsel;
    344 
    345346    SetReadyToSave();
    346347}
     348
  • trunk/MagicSoft/Mars/mhist/MHMcCT1CollectionArea.h

    r1821 r1823  
    1515    TH2D *fHistAll; //! all simulated showers
    1616    TH2D *fHistSel; //! the selected showers
    17 
    18     TH1D *fHistCol; //  the collection area
     17    TH2D *fHistCol; //  the collection area
    1918
    2019public:
     
    2726    void DrawSel(Option_t *option="");
    2827
    29     const TH1D *GetHist() const { return fHistCol; }
     28    const TH2D *GetHist() const { return fHistCol; }
    3029
    3130    void Draw(Option_t *option="");
    3231    TObject *DrawClone(Option_t *option="") const;
    3332
    34     void CalcEfficiency(Float_t theta);
     33    void CalcEfficiency();
    3534
    3635    ClassDef(MHMcCT1CollectionArea, 1)  // Data Container to calculate Collection Area
     
    3837
    3938#endif
     39
     40
Note: See TracChangeset for help on using the changeset viewer.