Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 1822) +++ trunk/MagicSoft/Mars/Changelog (revision 1823) @@ -4,8 +4,10 @@ * mhist/MHMcCT1CollectionArea.[h,cc] - - Added for calculations of collection area for CT1. + - Added for calculations of collection area for CT1.Contains three 2-d + histograms with axis energy vs theta angle: one histogram for all + events, one for analyzed events, one for the collection area. * mmontecarlo/MMcCT1CollectionAreaCalc.[h,cc] - - Added for the same reason. + - Added for the same reason. * macros/CT1collarea.C Index: trunk/MagicSoft/Mars/macros/CT1collarea.C =================================================================== --- trunk/MagicSoft/Mars/macros/CT1collarea.C (revision 1823) +++ trunk/MagicSoft/Mars/macros/CT1collarea.C (revision 1823) @@ -0,0 +1,82 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Thomas Bretz 12/2000 (tbretz@uni-sw.gwdg.de) +! +! Copyright: MAGIC Software Development, 2000-2001 +! +! +\* ======================================================================== */ + + +void CT1collarea(TString filename="MC_ON1.root", TString outname="") +{ + // + // first we have to create our empty lists + // + MParList parlist; + MTaskList tasklist; + + parlist.AddToList(&tasklist); + + // + // Setup out tasks: + // - First we have to read the events + // - Then we can fill the efficiency histograms + // + MReadMarsFile reader("Events", filename); + MMcCT1CollectionAreaCalc effi; + + tasklist.AddToList(&reader); + + tasklist.AddToList(&effi); + + // + // set up the loop for the processing + // + MEvtLoop magic; + magic.SetParList(&parlist); + + // + // Start to loop over all events + // + MProgressBar bar; + magic.SetProgressBar(&bar); + if (!magic.Eventloop()) + return; + + tasklist.PrintStatistics(); + + // + // Now the histogram we wanted to get out of the data is + // filled and can be displayed + // + parlist.FindObject("MHMcCT1CollectionArea")->DrawClone(); + + // + // Write histogram to a file in case an output + // filename has been supplied + // + if (outname.IsNull()) + return; + + TFile f(outname, "recreate"); + if (!f) + return; + parlist.FindObject("MHMcCT1CollectionArea")->Write(); +} + Index: trunk/MagicSoft/Mars/mhist/MHMcCT1CollectionArea.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHMcCT1CollectionArea.cc (revision 1822) +++ trunk/MagicSoft/Mars/mhist/MHMcCT1CollectionArea.cc (revision 1823) @@ -52,6 +52,5 @@ // we set the energy range from 10 Gev to 30000 GeV (in log 4.5 orders // of magnitude) and for each order we take 10 subdivision --> 45 xbins - // - // we set the radius range from 0 m to 500 m with 10 m bin --> 50 ybins + // we set the theta range from 12.5 to 48 deg, with 6 bins. // fName = name ? name : "MHMcCT1CollectionArea"; @@ -60,13 +59,18 @@ MBinning binsx; MBinning binsy; + binsx.SetEdgesLog(45, 10, 30000); - binsy.SetEdges (50, 0, 500); + const Double_t yedge[7] = {12.5, 17.5, 23.5, 29.5, 35.5, 42., 48.}; + const TArrayD yed(7,yedge); + binsy.SetEdges(yed); fHistAll = new TH2D; fHistSel = new TH2D; - fHistCol = new TH1D; - + fHistCol = new TH2D; + MH::SetBinning(fHistAll, &binsx, &binsy); MH::SetBinning(fHistSel, &binsx, &binsy); + MH::SetBinning(fHistCol, &binsx, &binsy); + fHistCol->SetName(fName); @@ -75,6 +79,6 @@ fHistCol->SetTitle(fTitle); - fHistAll->SetTitle("All showers - Radius vs Energy distribution"); - fHistSel->SetTitle("Selected showers - Radius vs Energy distribution"); + fHistAll->SetTitle("All showers - Theta vs Energy distribution"); + fHistSel->SetTitle("Selected showers - Theta vs Energy distribution"); fHistAll->SetDirectory(NULL); @@ -83,13 +87,14 @@ fHistAll->SetXTitle("E [GeV]"); - fHistAll->SetYTitle("r [m]"); + fHistAll->SetYTitle("theta [deg]"); fHistAll->SetZTitle("N"); fHistSel->SetXTitle("E [GeV]"); - fHistSel->SetYTitle("r [m]"); + fHistSel->SetYTitle("theta [deg]"); fHistSel->SetYTitle("N"); fHistCol->SetXTitle("E [GeV]"); - fHistCol->SetYTitle("A [m^{2}]"); + fHistCol->SetYTitle("theta [deg]"); + fHistCol->SetZTitle("A [m^{2}]"); } @@ -109,7 +114,7 @@ // Fill data into the histogram which contains the selected showers // -void MHMcCT1CollectionArea::FillSel(Double_t energy, Double_t radius) -{ - fHistSel->Fill(energy, radius); +void MHMcCT1CollectionArea::FillSel(Double_t energy, Double_t theta) +{ + fHistSel->Fill(energy, theta); } @@ -159,5 +164,5 @@ // - // This is necessary to get the expected bahviour of DrawClone + // This is necessary to get the expected behaviour of DrawClone // gROOT->SetSelectedPad(NULL); @@ -190,5 +195,5 @@ // and set the 'ReadyToSave' flag // -void MHMcCT1CollectionArea::CalcEfficiency(Float_t theta) +void MHMcCT1CollectionArea::CalcEfficiency() { // @@ -205,142 +210,139 @@ // - TH1D &histsel = *fHistSel->ProjectionX(); - - TH1D histall; - - TAxis &xaxis = *histsel.GetXaxis(); - - MH::SetBinning(&histall, &xaxis); - MH::SetBinning(fHistCol, &xaxis); - - Float_t emin1, emax1, emin2, emax2; - Float_t index, expo, k1, k2; - Float_t numevts1, numevts2; - Float_t r1, r2; // Impact parameter range (on ground). - - emin1 = 0; emax1 = 0; emin2 = 0; emax2 = 0; - expo = 0.; k1 = 0.; k2 = 0.; r1 = 0.; r2 = 0.; - numevts1 = 0; numevts2 = 0; - - if (theta > 0.26 && theta < 0.27) // 15 degrees - { - r1 = 0.; - r2 = 250.; //meters - emin1 = 300.; - emax1 = 400.; // Energies in GeV. - emin2 = 400.; - emax2 = 30000.; - numevts1 = 4000.; - numevts2 = 25740.; - } - else if (theta > 0.34 && theta < 0.35) // 20 degrees - { - r1 = 0.; - r2 = 263.; //meters - emin1 = 300.; - emax1 = 400.; // Energies in GeV. - emin2 = 400.; - emax2 = 30000.; - numevts1 = 6611.; - numevts2 = 24448.; - } - else if (theta > 0.47 && theta < 0.48) // 27 degrees - { - r1 = 0.; - r2 = 290.; //meters - emin1 = 300.; - emax1 = 400.; // Energies in GeV. - emin2 = 400.; - emax2 = 30000.; - numevts1 = 4000.; - numevts2 = 26316.; - } - else if (theta > 0.55 && theta < 0.56) // 32 degrees - { - r1 = 0.; - r2 = 350.; //meters - emin1 = 300.; - emax1 = 30000.; // Energies in GeV. - numevts1 = 33646.; - } - else if (theta > 0.68 && theta < 0.69) // 39 degrees - { - r1 = 0.; - r2 = 380.; //meters - emin1 = 300.; - emax1 = 30000.; // Energies in GeV. - numevts1 = 38415.; - } - else if (theta > 0.78 && theta < 0.79) // 45 degrees - { - r1 = 0.; - r2 = 565.; //meters - emin1 = 300.; - emax1 = 50000.; // Energies in GeV. - numevts1 = 30197.; - } - - index = 1.5; // Differential spectral Index. - expo = 1.-index; - k1 = numevts1 / (pow(emax1,expo) - pow(emin1,expo)); - k2 = numevts2 / (pow(emax2,expo) - pow(emin2,expo)); - - const Int_t nbinx = xaxis.GetNbins(); - - for (Int_t i=1; i<=nbinx; i++) + for (Int_t thetabin = 1; thetabin <= fHistAll->GetNbinsY(); thetabin++) { - const Float_t e1 = histall.GetBinLowEdge(i); - const Float_t e2 = histall.GetBinLowEdge(i+1); - - if (e1 < emin1 || e2 > emax2) + // This theta is not exactly the one of the MC events, just about + // the same: + Float_t theta = fHistAll->GetYaxis()->GetBinCenter(thetabin); + + Float_t emin1, emax1, emin2, emax2; + Float_t index, expo, k1, k2; + Float_t numevts1, numevts2; + Float_t r1, r2; // Impact parameter range (on ground). + + emin1 = 0; emax1 = 0; emin2 = 0; emax2 = 0; + expo = 0.; k1 = 0.; k2 = 0.; r1 = 0.; r2 = 0.; + numevts1 = 0; numevts2 = 0; + + if (theta > 14 && theta < 16) // 15 deg + { + r1 = 0.; + r2 = 250.; //meters + emin1 = 300.; + emax1 = 400.; // Energies in GeV. + emin2 = 400.; + emax2 = 30000.; + numevts1 = 4000.; + numevts2 = 25740.; + } + else if (theta > 20 && theta < 21) // 20.5 deg + { + r1 = 0.; + r2 = 263.; //meters + emin1 = 300.; + emax1 = 400.; // Energies in GeV. + emin2 = 400.; + emax2 = 30000.; + numevts1 = 6611.; + numevts2 = 24448.; + } + else if (theta > 26 && theta < 27) // 26.5 degrees + { + r1 = 0.; + r2 = 290.; //meters + emin1 = 300.; + emax1 = 400.; // Energies in GeV. + emax2 = emax1; emin2 = 400.; + emax2 = 30000.; + numevts1 = 4000.; + numevts2 = 26316.; + } + else if (theta > 32 && theta < 33) // 32.5 degrees + { + r1 = 0.; + r2 = 350.; //meters + emin1 = 300.; + emax1 = 30000.; // Energies in GeV. + emax2 = emax1; + numevts1 = 33646.; + } + else if (theta > 38 && theta < 39) // 38.75 degrees + { + r1 = 0.; + r2 = 380.; //meters + emin1 = 300.; + emax1 = 30000.; // Energies in GeV. + emax2 = emax1; + numevts1 = 38415.; + } + else if (theta > 44 && theta < 46) // 45 degrees + { + r1 = 0.; + r2 = 565.; //meters + emin1 = 300.; + emax1 = 50000.; // Energies in GeV. + emax2 = emax1; + numevts1 = 30197.; + } + + index = 1.5; // Differential spectral Index. + expo = 1.-index; + k1 = numevts1 / (pow(emax1,expo) - pow(emin1,expo)); + k2 = numevts2 / (pow(emax2,expo) - pow(emin2,expo)); + + for (Int_t i=1; i <= fHistAll->GetNbinsX(); i++) + { + const Float_t e1 = fHistAll->GetXaxis()->GetBinLowEdge(i); + const Float_t e2 = fHistAll->GetXaxis()->GetBinLowEdge(i+1); + + if (e1 < emin1 || e2 > emax2) continue; - Float_t events; - - if (e2emin2) - events = k2 * (pow(e2, expo) - pow(e1, expo)); - else - events = - k1 * (pow(emax1, expo) - pow(e1, expo))+ - k2 * (pow(e2, expo) - pow(emin2, expo)); - - histall.SetBinContent(i, events); - histall.SetBinError(i, sqrt(events)); + Float_t events; + + if (e2 <= emax1) + events = k1 * (pow(e2, expo) - pow(e1, expo)); + else if (e1 >= emin2) + events = k2 * (pow(e2, expo) - pow(e1, expo)); + else + events = + k1 * (pow(emax1, expo) - pow(e1, expo))+ + k2 * (pow(e2, expo) - pow(emin2, expo)); + + fHistAll->SetBinContent(i, thetabin, events); + fHistAll->SetBinError(i, thetabin, sqrt(events)); + } + + // ----------------------------------------------------------- + + const Float_t dr = TMath::Pi() * (r2*r2 - r1*r1); + + for (Int_t ix = 1; ix <= fHistAll->GetNbinsX(); ix++) + { + const Float_t Na = fHistAll->GetBinContent(ix,thetabin); + + if (Na <= 0) + continue; + + const Float_t Ns = fHistSel->GetBinContent(ix,thetabin); + + // Since Na is an estimate of the total number of showers generated + // in the energy bin, it may happen that Ns (triggered showers) is + // larger than Na. In that case, the bin is skipped: + + if (Na < Ns) + continue; + + const Double_t eff = Ns/Na; + const Double_t err = sqrt((1.-eff)*Ns)/Na; + + const Float_t area = dr * cos(theta*TMath::Pi()/180.); + + fHistCol->SetBinContent(ix, thetabin, eff*area); + fHistCol->SetBinError(ix, thetabin, err*area); + } } - // ----------------------------------------------------------- - - const Float_t dr = TMath::Pi() * (r2*r2 - r1*r1); - - for (Int_t ix=1; ix<=nbinx; ix++) - { - const Float_t Na = histall.GetBinContent(ix); - - if (Na <= 0) - continue; - - const Float_t Ns = histsel.GetBinContent(ix); - - // Since Na is an estimate of the total number of showers generated - // in the energy bin, it may happen that Ns (triggered showers) is - // larger than Na. In that case, the bin is skipped: - - if (Na < Ns) - continue; - - const Double_t eff = Ns/Na; - - const Double_t err = sqrt((1.-eff)*Ns)/Na; - - const Float_t area = dr * cos(theta); - - fHistCol->SetBinContent(ix, eff*area); - fHistCol->SetBinError(ix, err*area); - } - - delete &histsel; - SetReadyToSave(); } + Index: trunk/MagicSoft/Mars/mhist/MHMcCT1CollectionArea.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHMcCT1CollectionArea.h (revision 1822) +++ trunk/MagicSoft/Mars/mhist/MHMcCT1CollectionArea.h (revision 1823) @@ -15,6 +15,5 @@ TH2D *fHistAll; //! all simulated showers TH2D *fHistSel; //! the selected showers - - TH1D *fHistCol; // the collection area + TH2D *fHistCol; // the collection area public: @@ -27,10 +26,10 @@ void DrawSel(Option_t *option=""); - const TH1D *GetHist() const { return fHistCol; } + const TH2D *GetHist() const { return fHistCol; } void Draw(Option_t *option=""); TObject *DrawClone(Option_t *option="") const; - void CalcEfficiency(Float_t theta); + void CalcEfficiency(); ClassDef(MHMcCT1CollectionArea, 1) // Data Container to calculate Collection Area @@ -38,2 +37,4 @@ #endif + + Index: trunk/MagicSoft/Mars/mmontecarlo/MMcCT1CollectionAreaCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcCT1CollectionAreaCalc.cc (revision 1822) +++ trunk/MagicSoft/Mars/mmontecarlo/MMcCT1CollectionAreaCalc.cc (revision 1823) @@ -51,5 +51,4 @@ AddToBranchList("MMcEvt.fEnergy"); - AddToBranchList("MMcEvt.fImpact"); AddToBranchList("MMcEvt.fTelescopeTheta"); } @@ -72,6 +71,4 @@ fTotalNumSimulatedShowers = 0; fAllEvtsTriggered = kTRUE; - fAllEvtsAtSameTheta = kTRUE; - fTelescopeTheta = -1.; return kTRUE; @@ -81,15 +78,8 @@ { const Double_t energy = fMcEvt->GetEnergy(); - const Double_t impact = fMcEvt->GetImpact()/100.; - if (fTelescopeTheta < 0.) - fTelescopeTheta = fMcEvt->GetTelescopeTheta(); + Double_t TelescopeTheta = 180.*fMcEvt->GetTelescopeTheta()/TMath::Pi(); - if (fTelescopeTheta != fMcEvt->GetTelescopeTheta()) - { - fAllEvtsAtSameTheta = kFALSE; - } - - fCollArea->FillSel(energy, impact); + fCollArea->FillSel(energy, TelescopeTheta); return kTRUE; @@ -98,9 +88,4 @@ Bool_t MMcCT1CollectionAreaCalc::PostProcess() { - if ( ! fAllEvtsAtSameTheta ) - { - *fLog << err << dbginf << "ERROR: input data contain events at different zenith angles (not supported)... exiting." << endl; - return kFALSE; - } // // do the calculation of the effective area @@ -108,5 +93,5 @@ *fLog << inf << "Calculation Collection Area..." << endl; - fCollArea->CalcEfficiency(fTelescopeTheta); + fCollArea->CalcEfficiency(); return kTRUE; Index: trunk/MagicSoft/Mars/mmontecarlo/MMcCT1CollectionAreaCalc.h =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcCT1CollectionAreaCalc.h (revision 1822) +++ trunk/MagicSoft/Mars/mmontecarlo/MMcCT1CollectionAreaCalc.h (revision 1823) @@ -21,6 +21,4 @@ UInt_t fTotalNumSimulatedShowers; Bool_t fAllEvtsTriggered; - Bool_t fAllEvtsAtSameTheta; - Float_t fTelescopeTheta; public: