/* ======================================================================== *\ ! ! * ! * 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): A. Moralejo 3/2003 ! ! Copyright: MAGIC Software Development, 2000-2003 ! ! \* ======================================================================== */ ////////////////////////////////////////////////////////////////////////////// // // // MHMcCT1CollectionArea // // // ////////////////////////////////////////////////////////////////////////////// #include "MHMcCT1CollectionArea.h" #include #include #include "MMcEvt.hxx" #include "MH.h" #include "MBinning.h" #include "MParList.h" #include "MLog.h" #include "MLogManip.h" ClassImp(MHMcCT1CollectionArea); // -------------------------------------------------------------------------- // // Creates the three necessary histograms: // - selected showers (input) // - all showers (input) // - collection area (result) // MHMcCT1CollectionArea::MHMcCT1CollectionArea(const char *name, const char *title) { // // nbins, minEnergy, maxEnergy defaults: // we set the energy range from 100 Gev to 30000 GeV (in log, 3.5 orders // of magnitude) and for each order we take 10 subdivisions --> 35 xbins // we set the theta range from 12.5 to 48 deg, with 6 bins (the latter // choice has been done to make the bin centers as close as possible to // the actual zenith angles in the CT1 MC sample). // fName = name ? name : "MHMcCT1CollectionArea"; fTitle = title ? title : "Collection Area vs. log10 Energy"; fHistAll = new TH2D; fHistSel = new TH2D; fHistCol = new TH2D; fHistCol->SetName(fName); fHistAll->SetName("AllEvents"); fHistSel->SetName("SelectedEvents"); fHistCol->SetTitle(fTitle); fHistAll->SetTitle("All showers - Theta vs log10 Energy distribution"); fHistSel->SetTitle("Selected showers - Theta vs log10 Energy distribution"); fHistAll->SetDirectory(NULL); fHistSel->SetDirectory(NULL); fHistCol->SetDirectory(NULL); fHistAll->SetXTitle("log10 E [GeV]"); fHistAll->SetYTitle("theta [deg]"); fHistAll->SetZTitle("N"); fHistSel->SetXTitle("log10 E [GeV]"); fHistSel->SetYTitle("theta [deg]"); fHistSel->SetZTitle("N"); fHistCol->SetXTitle("log10 E [GeV]"); fHistCol->SetYTitle("theta [deg]"); fHistCol->SetZTitle("A [m^{2}]"); } // -------------------------------------------------------------------------- // // Delete the three histograms // MHMcCT1CollectionArea::~MHMcCT1CollectionArea() { delete fHistAll; delete fHistSel; delete fHistCol; } // -------------------------------------------------------------------------- // // Set the binnings and prepare the filling of the histograms // Bool_t MHMcCT1CollectionArea::SetupFill(const MParList *plist) { const MBinning* binsenergy = (MBinning*)plist->FindObject("BinningE"); const MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta"); if (!binsenergy || !binstheta) { *fLog << err << dbginf << "At least one MBinning not found... aborting."; *fLog << endl; return kFALSE; } SetBinning(fHistAll, binsenergy, binstheta); SetBinning(fHistSel, binsenergy, binstheta); SetBinning(fHistCol, binsenergy, binstheta); fHistAll->Sumw2(); fHistSel->Sumw2(); fHistCol->Sumw2(); return kTRUE; } // -------------------------------------------------------------------------- // // Fill data into the histogram which contains the selected showers // Bool_t MHMcCT1CollectionArea::Fill(const MParContainer *par, Double_t w) { MMcEvt &mcevt = *(MMcEvt*)par; fHistSel->Fill(log10(mcevt.GetEnergy()), kRad2Deg*mcevt.GetTelescopeTheta()); return kTRUE; } // -------------------------------------------------------------------------- // // Draw the histogram with all showers // void MHMcCT1CollectionArea::DrawAll(Option_t* option) { if (!gPad) MH::MakeDefCanvas(fHistAll); fHistAll->Draw(option); gPad->Modified(); gPad->Update(); } // -------------------------------------------------------------------------- // // Draw the histogram with the selected showers only. // void MHMcCT1CollectionArea::DrawSel(Option_t* option) { if (!gPad) MH::MakeDefCanvas(fHistSel); fHistSel->Draw(option); gPad->Modified(); gPad->Update(); } // -------------------------------------------------------------------------- // // Creates a new canvas and draws the histogram into it. // Be careful: The histogram belongs to this object and won't get deleted // together with the canvas. // TObject *MHMcCT1CollectionArea::DrawClone(Option_t* option) const { TCanvas &c = *MakeDefCanvas("CollArea", "Collection area plots", 600, 600); c.Divide(2,2); // // This is necessary to get the expected behaviour of DrawClone // gROOT->SetSelectedPad(NULL); c.cd(1); fHistCol->SetDirectory(NULL); fHistCol->DrawCopy(option); c.cd(2); fHistSel->SetDirectory(NULL); fHistSel->DrawCopy(option); c.cd(3); fHistAll->SetDirectory(NULL); fHistAll->DrawCopy(option); c.Modified(); c.Update(); return &c; } void MHMcCT1CollectionArea::Draw(Option_t* option) { if (!gPad) MH::MakeDefCanvas(fHistCol); fHistCol->Draw(option); gPad->Modified(); gPad->Update(); } // // Calculate the Efficiency (collection area) for the CT1 MC sample // and set the 'ReadyToSave' flag // void MHMcCT1CollectionArea::CalcEfficiency() { // // Here we estimate the total number of showers in each energy bin // from the known the energy range and spectral index of the generated // showers. This procedure is intended for the CT1 MC files. The total // number of generated events, collection area, spectral index etc will be // set here by hand, so make sure that the MC sample you are using is the // right one (check all these quantities in your files and compare with // what is written below. In some theta bins, there are two different // productions, with different energy limits but with the same spectral // slope. We account for this when calculating the original number of // events in each energy bin. // for (Int_t thetabin = 1; thetabin <= fHistAll->GetNbinsY(); thetabin++) { // 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 = pow(10.,fHistAll->GetXaxis()->GetBinLowEdge(i)); const Float_t e2 = pow(10.,fHistAll->GetXaxis()->GetBinLowEdge(i+1)); if (e1 < emin1 || e2 > emax2) continue; 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) { // // If energy is large, this case means that no or very few events // were generated at this energy bin. In this case we assign it // the effective area of the bin below it in energy. If energy is // below 1E4, it means that no events triggered -> eff area = 0 // if (fHistSel->GetXaxis()->GetBinLowEdge(ix) > 4.) { fHistCol->SetBinContent(ix, thetabin, fHistCol->GetBinContent(ix-1, thetabin)); fHistCol->SetBinError(ix, thetabin, fHistCol->GetBinError(ix-1, thetabin)); } 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 efferr = 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, efferr*area); } } SetReadyToSave(); }