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

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