Changeset 1823 for trunk/MagicSoft/Mars/mhist
- Timestamp:
- 03/13/03 22:25:00 (22 years ago)
- Location:
- trunk/MagicSoft/Mars/mhist
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/mhist/MHMcCT1CollectionArea.cc
r1821 r1823 52 52 // we set the energy range from 10 Gev to 30000 GeV (in log 4.5 orders 53 53 // 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. 56 55 // 57 56 fName = name ? name : "MHMcCT1CollectionArea"; … … 60 59 MBinning binsx; 61 60 MBinning binsy; 61 62 62 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); 64 66 65 67 fHistAll = new TH2D; 66 68 fHistSel = new TH2D; 67 fHistCol = new TH 1D;68 69 fHistCol = new TH2D; 70 69 71 MH::SetBinning(fHistAll, &binsx, &binsy); 70 72 MH::SetBinning(fHistSel, &binsx, &binsy); 73 MH::SetBinning(fHistCol, &binsx, &binsy); 74 71 75 72 76 fHistCol->SetName(fName); … … 75 79 76 80 fHistCol->SetTitle(fTitle); 77 fHistAll->SetTitle("All showers - Radiusvs Energy distribution");78 fHistSel->SetTitle("Selected showers - Radiusvs Energy distribution");81 fHistAll->SetTitle("All showers - Theta vs Energy distribution"); 82 fHistSel->SetTitle("Selected showers - Theta vs Energy distribution"); 79 83 80 84 fHistAll->SetDirectory(NULL); … … 83 87 84 88 fHistAll->SetXTitle("E [GeV]"); 85 fHistAll->SetYTitle(" r [m]");89 fHistAll->SetYTitle("theta [deg]"); 86 90 fHistAll->SetZTitle("N"); 87 91 88 92 fHistSel->SetXTitle("E [GeV]"); 89 fHistSel->SetYTitle(" r [m]");93 fHistSel->SetYTitle("theta [deg]"); 90 94 fHistSel->SetYTitle("N"); 91 95 92 96 fHistCol->SetXTitle("E [GeV]"); 93 fHistCol->SetYTitle("A [m^{2}]"); 97 fHistCol->SetYTitle("theta [deg]"); 98 fHistCol->SetZTitle("A [m^{2}]"); 94 99 } 95 100 … … 109 114 // Fill data into the histogram which contains the selected showers 110 115 // 111 void MHMcCT1CollectionArea::FillSel(Double_t energy, Double_t radius)112 { 113 fHistSel->Fill(energy, radius);116 void MHMcCT1CollectionArea::FillSel(Double_t energy, Double_t theta) 117 { 118 fHistSel->Fill(energy, theta); 114 119 } 115 120 … … 159 164 160 165 // 161 // This is necessary to get the expected b ahviour of DrawClone166 // This is necessary to get the expected behaviour of DrawClone 162 167 // 163 168 gROOT->SetSelectedPad(NULL); … … 190 195 // and set the 'ReadyToSave' flag 191 196 // 192 void MHMcCT1CollectionArea::CalcEfficiency( Float_t theta)197 void MHMcCT1CollectionArea::CalcEfficiency() 193 198 { 194 199 // … … 205 210 // 206 211 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++) 291 213 { 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) 296 299 continue; 297 300 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 } 311 344 } 312 345 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 generated327 // in the energy bin, it may happen that Ns (triggered showers) is328 // 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 345 346 SetReadyToSave(); 346 347 } 348 -
trunk/MagicSoft/Mars/mhist/MHMcCT1CollectionArea.h
r1821 r1823 15 15 TH2D *fHistAll; //! all simulated showers 16 16 TH2D *fHistSel; //! the selected showers 17 18 TH1D *fHistCol; // the collection area 17 TH2D *fHistCol; // the collection area 19 18 20 19 public: … … 27 26 void DrawSel(Option_t *option=""); 28 27 29 const TH 1D *GetHist() const { return fHistCol; }28 const TH2D *GetHist() const { return fHistCol; } 30 29 31 30 void Draw(Option_t *option=""); 32 31 TObject *DrawClone(Option_t *option="") const; 33 32 34 void CalcEfficiency( Float_t theta);33 void CalcEfficiency(); 35 34 36 35 ClassDef(MHMcCT1CollectionArea, 1) // Data Container to calculate Collection Area … … 38 37 39 38 #endif 39 40
Note:
See TracChangeset
for help on using the changeset viewer.