source: trunk/MagicSoft/Mars/mtemp/mucm/classes/MHFlux.cc@ 6977

Last change on this file since 6977 was 5450, checked in by marcos, 22 years ago
*** empty log message ***
File size: 14.9 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Marcos Lopex 11/2004 <mailto:marcos@gae.ucm.es>
19!
20! Copyright: MAGIC Software Development, 2000-2002
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26// //
27// MHFlux //
28// //
29// 3D-histogram in alpha vs. E-est and Theta //
30// //
31//////////////////////////////////////////////////////////////////////////////
32
33#include "MHFlux.h"
34
35#include <TCanvas.h>
36#include <THStack.h>
37#include <TLegend.h>
38#include <TStyle.h>
39#include <TAxis.h>
40#include <TF1.h>
41#include <TGraphErrors.h>
42#include <TPaveText.h>
43
44#include "MHillasSrc.h"
45#include "MEnergyEst.h"
46#include "MPointingPos.h"
47#include "MRawRunHeader.h"
48
49#include "MHExcessEnergyTheta.h"
50#include "MHMcCollectionArea.h"
51#include "MHEffectiveOnTime.h"
52
53#include "MBinning.h"
54#include "MParList.h"
55
56#include "MLog.h"
57#include "MLogManip.h"
58
59ClassImp(MHFlux);
60
61using namespace std;
62
63// --------------------------------------------------------------------------
64//
65// Default Constructor. It sets name and title of the histogram.
66//
67MHFlux::MHFlux(const char *name, const char *title)
68 : fHist("","",10,0,1, 10,0,1), fAverageFlux("","",1,0,1)
69{
70 //
71 // set the name and title of this object
72 //
73 fName = name ? name : "MHFlux";
74 fTitle = title ? title : "Flux vs. E and Theta";
75
76
77 fHist.SetDirectory(NULL);
78 fHist.SetName("Flux vs. E and Theta");
79 fHist.SetTitle("Flux vs. E and Theta");
80 fHist.SetXTitle("E [GeV]");
81 fHist.SetYTitle("\\Theta [\\circ]");
82 fHist.SetZTitle("Flux [TeV^{-1} s^{-1} cm^{-2}]");
83
84 fAverageFlux.SetDirectory(NULL);
85 fAverageFlux.SetName("Average Flux");
86 fAverageFlux.SetTitle("Average Flux");
87 fAverageFlux.SetXTitle("E [GeV]");
88 fAverageFlux.SetYTitle("Flux [TeV^{-1} s^{-1} cm^{-2}]");
89}
90
91// --------------------------------------------------------------------------
92//
93// Set binnings and prepare filling of the histogram
94//
95Bool_t MHFlux::SetupFill(const MParList *plist)
96{
97
98
99 return kTRUE;
100}
101
102// --------------------------------------------------------------------------
103//
104// Look in the parlist for MMcEvt or MPointingPos depending on the run type.
105//
106Bool_t MHFlux::ReInit(MParList *pList)
107{
108
109 return kTRUE;
110}
111
112
113
114// --------------------------------------------------------------------------
115//
116// Fill the histogram
117//
118Bool_t MHFlux::Fill(const MParContainer *par, const Stat_t w)
119{
120
121 return kTRUE;
122}
123
124
125// --------------------------------------------------------------------------
126//
127// Calc
128//
129void MHFlux::Calc(MHExcessEnergyTheta* hExcess, MHMcCollectionArea* hColArea, MHEffectiveOnTime* hEffTime)
130{
131 const TH2D* hex = hExcess->GetHist();
132 const TH1D* hca = hColArea->GetHist();
133 const TH1D* heot = &hEffTime->GetHEffOnTheta();
134
135 TAxis* axisEnergy = hex->GetXaxis();
136 const TAxis* axisTheta = hex->GetYaxis();
137 const Int_t energyBins = hex->GetXaxis()->GetNbins();
138 const Int_t thetaBins = hex->GetYaxis()->GetNbins();;
139
140 MH::SetBinning(&fHist,axisEnergy, axisTheta);
141
142
143 //
144 // Calculate flux for each energy and theta
145 //
146 for (Int_t iy=1; iy<=thetaBins; iy++) // loop on theta
147 {
148 // Get Effective Time [sec] and its error
149 Double_t t = heot->GetBinContent(iy);
150 const Double_t dt = heot->GetBinError(iy);
151
152 if (t < 1e-3)
153 t = 0.0;
154
155 for (Int_t ix=1; ix<=energyBins; ix++)
156 {
157 const Double_t n = hex->GetBinContent(ix,iy);
158 const Double_t dn = hex->GetBinError(ix,iy);
159
160 // Get AreaEff and its error
161 Double_t energy = axisEnergy->GetBinCenter(ix);
162 Int_t bin = hca->GetXaxis()->FindBin(energy);
163
164 // Get NumberExcessEvents and its error
165 const Double_t a = hca->GetBinContent(bin)*1e4; //cm^2
166 const Double_t da = hca->GetBinError(bin) *1e4; //cm^2
167
168 // energy bin width in TeV
169 const Double_t en = axisEnergy->GetBinWidth(ix)*1e-3; //TeV
170
171 //
172 // Check that we can calculate the flux for the current bin
173 //
174 if (t==0)
175 cout << "No_Ton ";
176 if (a==0)
177 cout << "No_Aeff ";
178 if (n==0)
179 cout << "No_Events ";
180 if ((t == 0) || (a == 0) || (n == 0)) {
181 cout << endl;
182 continue;
183 }
184
185 //
186 // Flux calculation and its error
187 //
188 const Double_t flux = n/(en*t*a);
189
190 // error propagation formula
191 const Double_t errN = dn/(en*a*t);
192 const Double_t errA = da * n/(en*t*a*a);
193 const Double_t errT = dt * n/(en*a*t*t);
194 const Double_t error = sqrt(errN*errN + errA*errA + errT*errT);
195
196 cout << dn << " " << en << " " << a << " " << t << endl;
197
198 fHist.SetBinContent(ix,iy,flux);
199 fHist.SetBinError(ix,iy,error);
200
201 } //energy
202 } //theta
203 fHist.Print("all");
204}
205
206
207
208// --------------------------------------------------------------------------
209//
210// Draw the histogram
211//
212void MHFlux::Draw(Option_t *opt)
213{
214
215 // --------------------------------------------------------------------
216 //
217 // Draw lego plot of Flux vs. Energy and Theta
218 //
219 TCanvas *c1 = new TCanvas("Flux vs. E and Theta","Flux vs. E and Theta");
220 c1->SetLogx();
221 c1->SetLogz();
222 fHist.SetStats(0);
223 fHist.Draw("lego");
224
225
226 // --------------------------------------------------------------------
227 //
228 // Draw the Flux for each Theta bin
229 //
230 TCanvas *c2 = new TCanvas("Fluxes for each Theta bin","Fluxes for each Theta bin");
231 c2->SetLogx();
232 c2->SetLogy();
233 c2->SetGridx();
234 c2->SetGridy();
235
236
237 THStack* hs = new THStack("Fluxes for each Theta bin","Fluxes for each Theta bin");
238
239 TLegend * leg = new TLegend(0.73,0.65,0.89,0.89);
240
241 TAxis* yaxis = fHist.GetYaxis();
242 const Int_t nbiny = fHist.GetYaxis()->GetNbins();
243
244 for(Int_t iy=1; iy<=nbiny; iy++)
245 {
246
247 TH1D* h1 = fHist.ProjectionX(Form("%d",iy),iy,iy,"e"); //<----- Option e is very important, otherwise the errors are not copied
248
249 if(h1->GetEntries()==0)
250 continue;
251
252 h1->SetLineColor(iy);
253 hs->Add(h1,"e1");
254 leg->AddEntry(h1,Form("\\theta = %.0f",yaxis->GetBinCenter(iy)),"l");
255
256// TCanvas *c = new TCanvas();
257// c->SetLogx();
258// c->SetLogy();
259// h1->DrawCopy("e");
260// h1->Print("all");
261 }
262
263
264
265 // --------------------------------------------------------------------
266 //
267 // Calculate and Draw the Flux average on Theta bins
268 //
269 fAverageFlux.SetStats(0);
270
271 MH::SetBinning(&fAverageFlux,fHist.GetXaxis());
272
273 for(int ix=1; ix<=fHist.GetXaxis()->GetNbins(); ix++) // energy
274 {
275 Double_t sumw = 0;
276 Double_t sumcontents = 0;
277 Double_t sumerrors = 0;
278
279 for(int iy=1; iy<=fHist.GetYaxis()->GetNbins(); iy++) // theta
280 {
281 Double_t weight = fHist.GetYaxis()->GetBinWidth(iy);
282 sumw += weight;
283 sumcontents += fHist.GetBinContent(ix,iy)*weight;
284 sumerrors += fHist.GetBinError(ix,iy)*weight;
285 }
286
287 fAverageFlux.SetBinContent(ix,sumcontents/sumw);
288 fAverageFlux.SetBinError(ix,sumerrors/sumw);
289 }
290
291 // for(int ix=1; ix<=fHist.GetXaxis()->GetNbins(); ix++) // energy
292// {
293// Double_t sumw = 0;
294// Double_t sumcontents = 0;
295
296// cout << " energy bin "<<endl;
297// for(int iy=1; iy<=fHist.GetYaxis()->GetNbins(); iy++) // theta
298// {
299
300// Double_t bincontent = fHist.GetBinContent(ix,iy);
301// Double_t binerror = fHist.GetBinError(ix,iy);
302
303// if( bincontent == 0 || binerror == 0 )
304// continue;
305// cout << binerror << endl;
306
307// Double_t weight = 1/(binerror*binerror);
308// sumw += weight;
309// sumcontents += bincontent*weight;
310
311// cout << " theta bin " << fHist.GetBinContent(ix,iy)<< " " <<weight
312// << endl;
313// }
314// cout << "*****************" << sumcontents << " "<< sumw << endl;
315
316// if(sumcontents == 0 || sumw == 0 )
317// continue;
318
319// fAverageFlux.SetBinContent(ix,sumcontents/sumw);
320// fAverageFlux.SetBinError(ix,TMath::Sqrt(1/sumw));
321// }
322
323
324
325
326 fAverageFlux.SetMarkerStyle(8);
327 fAverageFlux.SetLineColor(6);
328 hs->Add(&fAverageFlux,"pe1");
329 leg->AddEntry(&fAverageFlux,"Average on Theta","l");
330
331
332 c2->cd();
333 hs->Draw("nostack");
334 leg->Draw();
335
336
337 TCanvas *c3 = new TCanvas("Average Flux","Average Flux");
338 c3->SetLogx();
339 c3->SetLogy();
340 c3->SetGridx();
341 c3->SetGridy();
342 fAverageFlux.Draw();
343 fAverageFlux.Print("all");
344
345 //
346 // Fix the Average Flux to a power law
347 //
348 TF1* fluxfit = new TF1("f1","[0]*pow(x,-[1])",90,1500);
349 fluxfit->SetParNames("f0","a");
350 fluxfit->SetParameter(0,5.10986e-05);
351 fluxfit->SetParameter(1,2.4);
352 fluxfit->SetTitle("Flux fit");
353 fluxfit->SetLineColor(27);
354 fluxfit->SetLineWidth(3);
355
356 fAverageFlux.Fit("f1","R");
357
358
359 //
360 // Draw the Crab spectrum measured by HEGRA between 500 GeV and 80 TeV
361 //
362 TF1* CrabFlux = new TF1("CrabFlux","[0]*pow(x/1000.,-[1])",350,2000);
363 CrabFlux->SetParameter(0,2.83e-11);
364 CrabFlux->SetParameter(1,2.62);
365 CrabFlux->SetLineStyle(2);
366 CrabFlux->SetLineColor(4);
367 CrabFlux->Draw("same");
368
369 //
370 // Draw formula
371 //
372 TPaveText* func = new TPaveText(0.16, 0.22, 0.67, 0.28,"NDC");
373 func->AddText(Form("#frac{dF}{dE} = %.2e * E^{-%.2f} [#frac{ph}{cm^{2} s TeV}]",fluxfit->GetParameter(0),fluxfit->GetParameter(1)));
374 func->SetFillStyle(0);
375 func->SetBorderSize(0);
376 func->Draw();
377
378
379// //
380// // Draw "Preliminary"
381// //
382// TPaveText* lab = new TPaveText(0.33, 0.83, 0.68, 0.89,"NDC");
383// lab->AddText("preliminary");
384// lab->SetTextColor(2);
385// lab->SetFillStyle(0);
386// lab->SetBorderSize(0);
387// lab->Draw();
388
389
390
391 // ---------------------------------------------------------------------
392 //
393 // Integral flux
394 //
395 TH1D *hIntegral = (TH1D*)fAverageFlux.Clone();
396 hIntegral->GetListOfFunctions()->Clear();
397
398 Int_t nbinsx = fAverageFlux.GetNbinsX();
399
400 for(int i=1; i<=nbinsx; i++)
401 {
402 cout <<"Integral Flux: Binwidth:" << hIntegral->GetBinWidth(i) << endl;
403
404 hIntegral->SetBinContent(i,hIntegral->GetBinContent(i)*hIntegral->GetBinWidth(i)*1e-3);
405 hIntegral->SetBinError(i,hIntegral->GetBinError(i)*hIntegral->GetBinWidth(i)*1e-3);
406 }
407
408
409 for(int i=nbinsx-1; i>=1; i--)
410 {
411 Double_t integralsofar = hIntegral->GetBinContent(i+1);
412 Double_t current = hIntegral->GetBinContent(i);
413
414 Double_t currentE = hIntegral->GetBinError(i);
415 Double_t Esofar = hIntegral->GetBinError(i+1);
416
417 hIntegral->SetBinContent(i,(current+integralsofar));
418 hIntegral->SetBinError(i,TMath::Sqrt(currentE*currentE+Esofar*Esofar));
419 }
420
421
422 hIntegral->SetTitle("Integral Flux");
423 hIntegral->SetXTitle("E [GeV]");
424 hIntegral->SetYTitle("Integral Flux [s^{-1} cm^{-2}]");
425
426
427 TCanvas *c20 = new TCanvas();
428 c20->SetLogx();
429 c20->SetLogy();
430 c20->SetGridx();
431 c20->SetGridy();
432
433 hIntegral->Draw();
434
435
436
437
438 // --------------------------------------------------------------------
439 //
440 // E^2 * Flux
441 //
442 TH1D *hEscaledFlux = (TH1D*)fAverageFlux.Clone();
443 hEscaledFlux->GetListOfFunctions()->Clear();
444
445 nbinsx = hEscaledFlux->GetNbinsX();
446
447 for(int i=1; i<=nbinsx; i++)
448 {
449
450 Double_t energy = hEscaledFlux->GetBinLowEdge(i)*1e-3; // TeV
451 Double_t Flux = hEscaledFlux->GetBinContent(i);
452 Double_t dFlux = hEscaledFlux->GetBinError(i);
453
454 hEscaledFlux->SetBinContent(i,energy*energy*Flux);
455 hEscaledFlux->SetBinError(i,energy*energy*dFlux);
456 }
457
458 TCanvas *c40 = new TCanvas();
459 c40->SetLogx();
460 c40->SetLogy();
461 c40->SetGridx();
462 c40->SetGridy();
463
464 hEscaledFlux->SetTitle("Escaled Flux");
465 hEscaledFlux ->SetXTitle("E [GeV]");
466 hEscaledFlux ->SetYTitle("E^{2}*Flux [TeV s^{-1} cm^{-2}]");
467
468 hEscaledFlux->Draw();
469
470
471// // -----------------------------------------------------------------
472// //
473// // Graph move a 30 %
474// //
475// TCanvas *c4 = new TCanvas();
476// c4->SetLogx();
477// c4->SetLogy();
478// c4->SetGridx();
479// c4->SetGridy();
480
481// Int_t nbins = fAverageFlux.GetNbinsX();
482// TArrayD x(nbins),y(nbins),dx(nbins),dy(nbins);
483
484// for(int i=1; i<=nbins; i++)
485// {
486// x[i-1] = fAverageFlux.GetXaxis()->GetBinCenter(i)*.7;
487// y[i-1] = fAverageFlux.GetBinContent(i);
488// dx[i-1] = fAverageFlux.GetXaxis()->GetBinWidth(i)*0.62;
489// dy[i-1] = fAverageFlux.GetBinError(i);
490// }
491
492// TGraphErrors* gr = new TGraphErrors(fAverageFlux.GetNbinsX(), x.GetArray(), y.GetArray(), dx.GetArray(), dy.GetArray());
493
494// gr->SetMarkerStyle(8);
495// gr->Draw("Ap");
496// gr->Print("all");
497
498// //
499// TF1* fluxfit2 = new TF1("f2","[0]*pow(x,-[1])",70,2500);
500// fluxfit2->SetParNames("f0","a");
501// fluxfit2->SetParameter(0,5.10986e-05);
502// fluxfit2->SetParameter(1,2.4);
503// fluxfit2->SetTitle("Flux fit");
504// fluxfit2->SetLineColor(27);
505// fluxfit2->SetLineWidth(3);
506
507
508// gr->Fit("f2","R");
509
510// gr->SetTitle("");
511// gr->SetMaximum(1e-3);
512// gr->SetMinimum(1e-12);
513
514// TLegend* leg2 = new TLegend(0.67,0.72,0.89,0.89);
515
516// leg2->AddEntry(gr,"MAGIC Sept. 2004","p");
517
518// gr->SetMarkerStyle(8);
519// // gr->SetLineColor(6);
520
521// // //
522// // // Draw the Crab spectrum measured by HEGRA between 500 GeV and 80 TeV
523// // //
524// // TF1* CrabFlux = new TF1("CrabFlux","[0]*pow(x/1000.,-[1])",350,2000);
525// // CrabFlux->SetParameter(0,2.83e-11);
526// // CrabFlux->SetParameter(1,2.62);
527// // CrabFlux->SetLineStyle(2);
528// // CrabFlux->SetLineColor(4);
529// // CrabFlux->Draw("same");
530
531// leg2->AddEntry(CrabFlux,"HEGRA ApJ 614","l");
532// leg2->Draw();
533// lab->Draw();
534
535// TPaveText* func2 = new TPaveText(0.16, 0.22, 0.67, 0.28,"NDC");
536// func2->AddText(Form("#frac{dF}{dE} = %.2e * E^{-%.2f} [#frac{ph}{cm^{2} s TeV}]",fluxfit2->GetParameter(0),fluxfit2->GetParameter(1)));
537// func2->SetFillStyle(0);
538// func2->SetBorderSize(0);
539// func2->Draw();
540
541
542// gr->GetHistogram()->SetXTitle("E [GeV]");
543// gr->GetHistogram()->SetYTitle("Flux [TeV^{-1} s^{-1} cm^{-2}]");
544
545// TH1F* h = gr->GetHistogram();
546// TCanvas *c12 = new TCanvas();
547// h->Draw();
548// cout << "Integral flux = "<< h->Integral("width") << endl;
549}
Note: See TracBrowser for help on using the repository browser.