source: branches/MarsISDCBranchBasedOn17887/macros/tutorials/weights.C@ 18477

Last change on this file since 18477 was 7159, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 3.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 Lopez, 10/2003 <mailto:marcos@gae.ucm.es>
19!
20! Copyright: MAGIC Software Development, 2000-2004
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27// weights.C
28// =========
29//
30// This macro shows how to use the class MMcWeightEnergySpecCalc
31// to convert the energy spectrum of the MC showers generated with Corsika,
32// to a different one.
33//
34//////////////////////////////////////////////////////////////////////////////
35
36Double_t myfunction(Double_t *x, Double_t *par)
37{
38 Double_t xx = x[0];
39
40 return pow(xx,-2)*exp(-xx/20);
41}
42
43void weights(TString filename="/up1/data/Magic-MC/CameraAll/Gammas/zbin0/Gamma_zbin0_0_7_1000to1009_w0-4:4:2.root")
44{
45
46 //
47 // PartList
48 //
49 MParList parlist;
50 MTaskList tasklist;
51
52 MHMcEnergyImpact h1("h1");
53 MHMcEnergyImpact h2("h2");
54 parlist.AddToList(&h1);
55 parlist.AddToList(&h2);
56
57 MBinning binsenergy("BinningEnergy");
58 binsenergy.SetEdgesLog(100, 1, 1e5);
59 parlist.AddToList(&binsenergy);
60
61 MBinning binsimpact("BinningImpact");
62 binsimpact.SetEdges(100, 0, 450);
63 parlist.AddToList(&binsimpact);
64
65 parlist.AddToList(&tasklist);
66
67
68 //
69 // TaskList
70 //
71 MReadMarsFile reader("Events", filename);
72 reader.EnableBranch("fEnergy");
73 reader.EnableBranch("fImpact");
74
75
76 // -------------------------------------------------------------
77 //
78 // Option 1. Just change the slope of the MC power law spectrum
79 //
80 //MMcWeightEnergySpecCalc wcalc(-2.0); //<-- Uncomment me
81
82 //
83 // Option 2. A completely differente specturm pass as a TF1 function
84 // e.g. spectrum with exponential cutoff
85 //
86 //TF1 spec("spectrum","pow(x,[0])*exp(-x/[1])"); //<-- Uncomment me
87 //spec->SetParameter(0,-2.0); //<-- Uncomment me
88 //spec->SetParameter(1,50); //<-- Uncomment me
89 //MMcWeightEnergySpecCalc wcalc(spec); //<-- Uncomment me
90
91 //
92 // Option 3. A completely differente specturm pass as a cahr*
93 //
94 //char* func = "pow(x,-2)"; //<-- Uncomment me
95 //MMcWeightEnergySpecCalc wcalc(func); //<-- Uncomment me
96
97 //
98 // Option 4. A completely differente specturm pass as a c++ function
99 //
100 MMcWeightEnergySpecCalc wcalc((void*)myfunction); //<-- Uncomment me
101 //
102 //-------------------------------------------------------------
103
104 MFillH hfill(&h1,"MMcEvt");
105 MFillH hfill2(&h2,"MMcEvt");
106 hfill2.SetWeight("MWeight");
107
108 tasklist.AddToList(&reader);
109 tasklist.AddToList(&wcalc);
110 tasklist.AddToList(&hfill);
111 tasklist.AddToList(&hfill2);
112
113 //
114 // EventLoop
115 //
116 MEvtLoop magic;
117 magic.SetParList(&parlist);
118
119 if (!magic.Eventloop())
120 return;
121
122 tasklist.PrintStatistics();
123 parlist.Print();
124
125 //
126 // Draw the Results
127 //
128 TCanvas *c = new TCanvas();
129 c->SetLogy();
130 c->SetLogx();
131
132 TH1D* hist1 = (h1->GetHist())->ProjectionX();
133 TH1D* hist2 = (h2->GetHist())->ProjectionX();
134 hist2->SetLineColor(2);
135
136 hist1->DrawClone();
137 hist2->DrawClone("same");
138}
Note: See TracBrowser for help on using the repository browser.