source: trunk/MagicSoft/Mars/mtemp/mifae/programs/optimizeCuts.cc@ 4211

Last change on this file since 4211 was 4098, checked in by rico, 22 years ago
*** empty log message ***
File size: 9.0 KB
Line 
1////////////////////////////////////////////////////////////////////////////////////
2//
3// _____ Optimize cuts _____
4//
5// Take as input root files with hillas parameters and compute the optimal cuts by
6// mapping the width/length plane (only upper cuts so far)
7//
8// Jose Flix <jflix@ifae.es>
9// Javier Rico <jrico@ifae.es>
10////////////////////////////////////////////////////////////////////////////////////
11
12#include <fstream>
13#include <iostream>
14
15#include "TString.h"
16#include "TChain.h"
17#include "TH1F.h"
18
19#include "MArgs.h"
20#include "MLog.h"
21
22using namespace std;
23
24Bool_t readDatacards(TString& filename);
25void optimizeCuts();
26
27
28//-----------------------------------------------------------------------------
29// declaration of variables read from datacards
30//-----------------------------------------------------------------------------
31
32TString onFile;
33TString offFile;
34TString outputFile;
35ULong_t nmaxevents=999999999;
36Float_t lowdistcut=0;
37Float_t uppdistcut=10;
38Float_t sizecut=1800;
39UInt_t onRate=1;
40Float_t lowlgcut=0.01;
41Float_t upplgcut=0.3;
42Float_t lgstep=0.005;
43Float_t lowwdcut=0.01;
44Float_t uppwdcut=0.3;
45Float_t wdstep=0.005;
46
47//-----------------------------------------------------------------------------
48// constants
49//-----------------------------------------------------------------------------
50
51const TString defaultcard="optimizecuts.datacard";
52const Float_t conver = 189./0.6; // conversion factor degrees to mm
53const Int_t nbins=18; //number of bins in alpha plots
54const Float_t frontiere=20.; // [deg] value below (above) wich signal (background is computed (normalized)
55
56//-----------------------------------------------------------------------------
57
58static void Usage()
59{
60 gLog <<endl;
61 gLog << "Usage is:" << endl;
62 gLog << " optimizeCuts [-h] [-?] <datacards>" << endl << endl;
63 gLog << " <datacards>: datacards file name (default " << defaultcard <<")" << endl;
64 gLog << " -?/-h: This help" << endl << endl;
65}
66
67//-----------------------------------------------------------------------------
68int main(int argc, char **argv)
69{
70 // evaluate arguments
71 MArgs arg(argc, argv);
72 if (arg.HasOnly("-?") || arg.HasOnly("-h") || arg.HasOnly("--help"))
73 {
74 Usage();
75 return -1;
76 }
77
78 TString datacard = arg.GetArgumentStr(0);
79 if(!datacard.Length())
80 datacard = defaultcard;
81
82 if(!readDatacards(datacard))
83 {
84 cout << "Error reading datacards. Stoping" << endl;
85 return -1;
86 }
87 optimizeCuts();
88}
89
90//-----------------------------------------------------------------------------
91void optimizeCuts()
92{
93 // create the TChains for ON and OFF data
94 TChain* ton= new TChain("Parameters");
95 TChain* toff= new TChain("Parameters");
96 if(!ton->Add(onFile))
97 {
98 cout << "On-data file not found or not valid format" << endl;
99 return;
100 }
101 if(!toff->Add(offFile))
102 {
103 cout << "Off-data file not found or not valid format" << endl;
104 return;
105 }
106
107 ofstream fout;
108 fout.open(outputFile.Data());
109
110 // define aliases
111 ton->SetAlias("length","MHillas.fLength*0.6/189");
112 ton->SetAlias("width","MHillas.fWidth*0.6/189");
113 ton->SetAlias("dist","MHillasSrc.fDist*0.6/189");
114 ton->SetAlias("conc","MNewImagePar.fConc");
115 ton->SetAlias("size","MHillas.fSize");
116 ton->SetAlias("event","MRawEvtHeader.fDAQEvtNumber");
117 ton->SetAlias("alpha","abs(MHillasSrc.fAlpha)");
118
119 toff->SetAlias("length","MHillas.fLength*0.6/189");
120 toff->SetAlias("width","MHillas.fWidth*0.6/189");
121 toff->SetAlias("dist","MHillasSrc.fDist*0.6/189");
122 toff->SetAlias("conc","MNewImagePar.fConc");
123 toff->SetAlias("size","MHillas.fSize");
124 toff->SetAlias("event","MRawEvtHeader.fDAQEvtNumber");
125 toff->SetAlias("alpha","abs(MHillasSrc.fAlpha)");
126
127
128 // general cut
129 char gencut[256];
130 char evecut[256];
131 sprintf(evecut,"event%%%d==0",onRate);
132 sprintf(gencut,"size>%f && dist>%f && dist<%f",sizecut,lowdistcut,uppdistcut);
133 cout << "General cut " << gencut << " (" << evecut << " for on-data)" << endl;
134
135 // loop on widht and length cuts
136 Float_t length=lowlgcut;
137 while(length<=upplgcut)
138 {
139 Float_t width=lowwdcut;
140 while(width<=uppwdcut)
141 {
142 TH1F* onhisto = new TH1F("OnHist" ,"Alpha Plot",nbins,0,90);
143 TH1F* offhisto = new TH1F("OffHist","Alpha Plot",nbins,0,90);
144
145 // define the width/length cut
146 char varcut[256];
147 sprintf(varcut,"length<%f && width<%f",length,width);
148 cout << "Cutting " << varcut << endl;
149
150 // define the on/off data cut
151 char offcut[1024];
152 sprintf(offcut,"%s && %s",gencut,varcut);
153 char oncut[1024];
154 sprintf(oncut,"%s && %s",evecut,offcut);
155
156 // Fill the histos
157 ton->Draw("alpha>>OnHist",oncut);
158 toff->Draw("alpha>>OffHist",offcut);
159
160 // Normalize histos
161 const Int_t inibin = (Int_t)(frontiere/90.*nbins+1);
162 Float_t level=0;
163 Float_t leveloff=0;
164 for(Int_t ibin = inibin; ibin<=nbins;ibin++)
165 {
166 level+=onhisto->GetBinContent(ibin);
167 leveloff+=offhisto->GetBinContent(ibin);
168 }
169 offhisto->Sumw2(); // needed to compute correct errors after normalization
170 const Float_t norm = leveloff ? level/leveloff: 1;
171 cout << "Normalizing by factor " << norm <<endl;
172 offhisto->Scale(norm);
173
174 // compute significance/excess
175 Float_t sig=0,bg=0,esig=0,ebg=0;
176 Float_t significance=0;
177 const Int_t signbins = inibin-1;
178 for(Int_t ibin = 1; ibin<=signbins;ibin++)
179 {
180 sig += onhisto->GetBinContent(ibin);
181 esig += onhisto->GetBinError(ibin)*onhisto->GetBinError(ibin);
182 bg += offhisto->GetBinContent(ibin);
183 ebg += offhisto->GetBinError(ibin)*offhisto->GetBinError(ibin);
184 }
185 Float_t error= TMath::Sqrt(esig+ebg);
186 significance = error ? (sig-bg)/error : 0;
187
188 cout << "Excess: " << sig-bg << endl;
189 cout << "N bkg: " << bg << endl;
190 cout << "Significance: " << significance << endl;
191
192 // save the result in file
193 fout << width << '\t'
194 << length << '\t'
195 << sig-bg << '\t'
196 << significance << '\t' << endl;
197
198 delete onhisto;
199 delete offhisto;
200
201 width+=wdstep;
202 }
203 length+=lgstep;
204 }
205
206 fout.close();
207}
208//-----------------------------------------------------------------------
209
210Bool_t readDatacards(TString& filename)
211{
212 ifstream ifun(filename.Data());
213 if(!ifun)
214 {
215 cout << "File " << filename << " not found" << endl;
216 return kFALSE;
217 }
218
219 TString word;
220
221 while(ifun >> word)
222 {
223 // skip comments
224 if(word[0]=='/' && word[1]=='/')
225 {
226 while(ifun.get()!='\n'); // skip line
227 continue;
228 }
229
230 // number of events
231 if(strcmp(word.Data(),"NEVENTS")==0)
232 ifun >> nmaxevents;
233
234 // number of events
235 if(strcmp(word.Data(),"ONRATE")==0)
236 ifun >> onRate;
237
238 // input file name (on data)
239 if(strcmp(word.Data(),"ONFILES")==0)
240 {
241 if(onFile.Length())
242 cout << "readDataCards Warning: overriding on-data file name" << endl;
243 ifun >> onFile;
244 }
245
246 // input file name (off data)
247 if(strcmp(word.Data(),"OFFFILES")==0)
248 {
249 if(offFile.Length())
250 cout << "readDataCards Warning: overriding off-data file name" << endl;
251 ifun >> offFile;
252 }
253
254 // output file name
255 if(strcmp(word.Data(),"OUTFILE")==0)
256 {
257 if(outputFile.Length())
258 cout << "readDataCards Warning: overriding output file name" << endl;
259 ifun >> outputFile;
260 }
261
262 // size cut
263 if(strcmp(word.Data(),"SIZECUT")==0)
264 ifun >> sizecut;
265
266 // dist cuts
267 if(strcmp(word.Data(),"DISTCUTS")==0)
268 {
269 ifun >> lowdistcut ;
270 ifun >> uppdistcut ;
271 }
272
273 // length cuts
274 if(strcmp(word.Data(),"LENGTHCUTS")==0)
275 {
276 ifun >> lowlgcut;
277 ifun >> upplgcut;
278 ifun >> lgstep;
279 }
280
281 // width cuts
282 if(strcmp(word.Data(),"WIDTHCUTS")==0)
283 {
284 ifun >> lowwdcut;
285 ifun >> uppwdcut;
286 ifun >> wdstep;
287 }
288
289 }
290
291 // check compulsory values
292 if(!onFile.Length())
293 {
294 cout << "No on-data file name specified" << endl;
295 return kFALSE;
296 }
297
298 if(!offFile.Length())
299 {
300 cout << "No off-data file name specified" << endl;
301 return kFALSE;
302 }
303
304 if(!outputFile.Length())
305 {
306 cout << "No output file name specified" << endl;
307 return kFALSE;
308 }
309
310 // Dump read values
311 cout << "************************************************" << endl;
312 cout << "* Datacards read from file " << filename << endl;
313 cout << "************************************************" << endl;
314 cout << "Maximum number of input events: " << nmaxevents << endl;
315 cout << "On-data acceptance rate: 1/" << onRate << endl;
316 cout << "On-data input file name(s): " << onFile << endl;
317 cout << "Off-data input file name(s): " << offFile << endl;
318 cout << "Output file name: " << outputFile << endl;
319 cout << "Size cut: " << sizecut << endl;
320 cout << "Dist cuts (deg): [" << lowdistcut<< ","<<uppdistcut<< "]" << endl;
321 cout << "Scanning length range (deg): [" << lowlgcut<< ","<<upplgcut<< "], with step (deg): "<< lgstep << endl;
322 cout << "Scanning width range (deg): [" << lowwdcut<< ","<<uppwdcut<< "], with step (deg): "<< wdstep << endl;
323 cout << "***********" << endl << endl;
324
325 return kTRUE;
326}
Note: See TracBrowser for help on using the repository browser.