source: branches/Mars_McMismatchStudy/macros/trigrate.C@ 19116

Last change on this file since 19116 was 2377, checked in by tbretz, 21 years ago
*** empty log message ***
File size: 6.2 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): Thomas Bretz, 12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2001
21!
22\* ======================================================================== */
23
24/*
25 ! Modified 4/7/2002, Abelardo Moralejo:
26 ! Added one optional input parameter: a camera .root file containing
27 ! pure NSB events. One such file is generated running the camera over an
28 ! "empty" reflector file, with the NSB option on, and asking the camera
29 ! program (see camera program manual) to do the simulation even if no
30 ! photoelectron from the shower arrives at the camera. One also needs to
31 ! write to the output file all the events, both triggered and untriggered
32 ! (see again camera manual). These nsb camera files must contain the same
33 ! trigger conditions as the proton file.
34 !
35 ! If no nsb file is supplied, the macro will assume no triggers from
36 ! pure NSB fluctuations.
37 */
38
39Float_t GetNSBEvents(TString name, Float_t *BgR, int dim)
40{
41 Int_t numnsbevents;
42
43 TFile bgfile(name);
44 TTree *events = (TTree*) bgfile.Get("Events");
45
46 TH1F h("h","",5,.5,5.5);
47
48 UInt_t from = dim>0 ? 1 : -dim;
49 UInt_t to = dim>0 ? dim : -dim;
50
51 cout << endl << "Calculating NSB triggers from file " << name << "..." << endl << endl;
52
53 for (UInt_t i = from; i <= to; i++)
54 {
55 TString cond("MMcTrig;");
56 cond += i;
57 cond += ".";
58 TBranch *b = events->GetBranch(cond);
59 MMcTrig *mctrig;
60 b->SetAddress(&mctrig);
61
62 Int_t tottrig = 0.;
63
64 UInt_t imax = b->GetEntries();
65
66 for (UInt_t iev = 0; iev < imax; iev++)
67 {
68 b->GetEvent(iev);
69 tottrig += mctrig->GetFirstLevel();
70 }
71 // Set total number of L1 triggers:
72 BgR[dim>0? i-1: 0] = (Float_t) tottrig;
73
74 numnsbevents = (Float_t) imax;
75 }
76
77 return numnsbevents;
78}
79
80
81void trigrate(int dim=0, char *filename = "data/camera.root",
82 char *nsbfile = NULL)
83{
84 // The dim parameter has three possible values:
85 // dim : = 0 -> root file with 1 trigger condition.
86 // > 0 -> number of trigger condition to be analised
87 // in multi conditon file.
88 // < 0 -> selects the -dim trigger condition.
89 //
90 // first we have to create our empty lists
91 //
92 MParList parlist;
93 MTaskList tasklist;
94
95 //
96 // Setup the parameter list.
97 // - we do not need to create any other container. All of them
98 // are created automatically without loss - we don't have to
99 // access them-
100 //
101 // - we need to create MHMcRate only. The other containers
102 // are created automatically without loss - we don't have to
103 // access them-
104 // - MHMcRate must be created by us because we need the pointer
105 // to it and if it would get created automatically it would also be
106 // deleted automatically
107 // - Actually, depending on using a single trigger option MonteCarlo
108 // file or a multyple trigger option, a MHMcRate or an array of
109 // MHMcRate are needed.
110 //
111 parlist.AddToList(&tasklist);
112
113 //
114 // You don't have to add the MHMcRate container here by hand.
115 // But if you want to print or display these containers later on
116 // it is necessary (Rem: No printing or displaying is done in this
117 // macro yet)
118 //
119 UInt_t from = dim>0 ? 1 : -dim;
120 UInt_t to = dim>0 ? dim : -dim;
121
122 Int_t num = to-from+1;
123
124 TObjArray hists(MParList::CreateObjList("MHMcRate", from, to));
125 hists.SetOwner();
126
127 //
128 // Check if the list really contains the right number of histograms
129 //
130 if (hists.GetEntriesFast() != num)
131 return;
132
133 // Set for each MHMcRate object the trigger condition number in the
134 // camera file (for the case of camera files with several conditions,
135 // produced with the trigger_loop option.
136 //
137
138 if (dim < 0)
139 ((MHMcRate*)(hists[0]))->SetTriggerCondNum((Short_t)(-dim));
140 if (dim > 0)
141 for (Short_t i = from ; i <= to; i++)
142 ((MHMcRate*)(hists[i-1]))->SetTriggerCondNum(i);
143
144 //
145 // Add the histograms to the paramater list.
146 //
147 parlist.AddToList(&hists);
148
149 //
150 // Setup out tasks:
151 // - First we have to read the events
152 // - Then we can calculate rates, for what the number of
153 // triggered showers from a empty reflector file for the
154 // analised trigger conditions should be set (BgR[])
155 //
156
157 MReadMarsFile reader("Events", filename);
158 tasklist.AddToList(&reader);
159
160 // Now we have to build the BgR array, containing the number
161 // of triggers (may be more than 1 trigger/event!) from the
162 // numnsbevents simulated in the nsbfile (3rd input parameter).
163 // If no nsbfile is supplied, we assume no triggers from NSB
164
165 Float_t* BgR = new Float_t[num];
166 memset(BgR, 0, num*sizeof(Float_t));
167
168 Float_t numnsbevents = 5.e4; // some default value.
169 if (nsbfile)
170 numnsbevents = GetNSBEvents(nsbfile, BgR, dim);
171
172 cout << "Number of Trigger conditions: " << num << endl;
173
174 MMcTriggerRateCalc rate(dim, BgR, numnsbevents);
175
176 tasklist.AddToList(&rate);
177
178 //
179 // set up the loop for the processing
180 //
181 MEvtLoop magic;
182 magic.SetParList(&parlist);
183
184 //
185 // Start to loop over all events
186 //
187 MProgressBar bar;
188 magic.SetProgressBar(&bar);
189 if (!magic.Eventloop())
190 return;
191
192 delete BgR;
193
194 hists.Print();
195
196 if (num > 1)
197 {
198 gStyle->SetOptStat(0);
199 rate->Draw();
200
201 TFile f("ratehists.root", "recreate");
202
203 rate->GetHist(2)->Write();
204 rate->GetHist(3)->Write();
205 rate->GetHist(4)->Write();
206 }
207
208
209}
Note: See TracBrowser for help on using the repository browser.