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

Last change on this file since 4308 was 4234, checked in by rico, 20 years ago
Changelog
File size: 10.1 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("dist","MHillasSrc.fDist*0.6/189");
116 ton->SetAlias("conc","MNewImagePar.fConc");
117 ton->SetAlias("size","MHillas.fSize");
118 ton->SetAlias("event","MRawEvtHeader.fDAQEvtNumber");
119 ton->SetAlias("alpha","abs(MHillasSrc.fAlpha)");
120
121 toff->SetAlias("length","MHillas.fLength*0.6/189");
122 toff->SetAlias("width","MHillas.fWidth*0.6/189");
123 toff->SetAlias("dist","MHillasSrc.fDist*0.6/189");
124 toff->SetAlias("conc","MNewImagePar.fConc");
125 toff->SetAlias("size","MHillas.fSize");
126 toff->SetAlias("event","MRawEvtHeader.fDAQEvtNumber");
127 toff->SetAlias("alpha","abs(MHillasSrc.fAlpha)");
128
129
130 // general cut
131 char gencut[256];
132 char evecut[256];
133 sprintf(evecut,"event%%%d==0",onRate);
134 sprintf(gencut,"dist>%f && dist<%f",lowdistcut,uppdistcut);
135 cout << "General cut " << gencut << " (" << evecut << " for on-data)" << endl;
136
137 Bool_t isDiff = kFALSE;
138 if(nsizebins<0)
139 {
140 nsizebins*=-1;
141 isDiff=kTRUE;
142 }
143 else
144 nsizebins+=1;
145
146 // loop on size, width and length cuts
147 for(Int_t isb=0;isb<nsizebins;isb++) // loop on size
148 {
149 cout << isb << endl;
150 Float_t minsize = isDiff ?
151 TMath::Power(10,TMath::Log10(lowsizecut)+isb*(TMath::Log10(uppsizecut)-TMath::Log10(lowsizecut))/nsizebins) :
152 TMath::Power(10,TMath::Log10(lowsizecut)+isb*(TMath::Log10(uppsizecut)-TMath::Log10(lowsizecut))/(nsizebins-1));
153 Float_t maxsize = isDiff ? TMath::Power(10,TMath::Log10(lowsizecut)+(isb+1)*(TMath::Log10(uppsizecut)-TMath::Log10(lowsizecut))/nsizebins) : 9999999999.;
154
155 Float_t length=lowlgcut;
156 while(length<=upplgcut) // loop on length
157 {
158 Float_t width=lowwdcut;
159 while(width<=uppwdcut) // loop on width
160 {
161 TH1F* onhisto = new TH1F("OnHist" ,"Alpha Plot",nbins,0,90);
162 TH1F* offhisto = new TH1F("OffHist","Alpha Plot",nbins,0,90);
163
164 // define the width/length cut
165 char varcut[256];
166 sprintf(varcut,"size>%f && size<%f && length<%f && width<%f",minsize,maxsize,length,width);
167 cout << "Cutting " << varcut << endl;
168
169 // define the on/off data cut
170 char offcut[1024];
171 sprintf(offcut,"%s && %s",gencut,varcut);
172 char oncut[1024];
173 sprintf(oncut,"%s && %s",evecut,offcut);
174
175 // Fill the histos
176 ton->Draw("alpha>>OnHist",oncut);
177 toff->Draw("alpha>>OffHist",offcut);
178
179 // Normalize histos
180 const Int_t inibin = (Int_t)(frontiere/90.*nbins+1);
181 Float_t level=0;
182 Float_t leveloff=0;
183 for(Int_t ibin = inibin; ibin<=nbins;ibin++)
184 {
185 level+=onhisto->GetBinContent(ibin);
186 leveloff+=offhisto->GetBinContent(ibin);
187 }
188 offhisto->Sumw2(); // needed to compute correct errors after normalization
189 const Float_t norm = leveloff ? level/leveloff: 1;
190 cout << "Normalizing by factor " << norm <<endl;
191 offhisto->Scale(norm);
192
193 // compute significance/excess
194 Float_t sig=0,bg=0,esig=0,ebg=0;
195 Float_t significance=0;
196 const Int_t signbins = inibin-1;
197 for(Int_t ibin = 1; ibin<=signbins;ibin++)
198 {
199 sig += onhisto->GetBinContent(ibin);
200 esig += onhisto->GetBinError(ibin)*onhisto->GetBinError(ibin);
201 bg += offhisto->GetBinContent(ibin);
202 ebg += offhisto->GetBinError(ibin)*offhisto->GetBinError(ibin);
203 }
204 Float_t error= TMath::Sqrt(esig+ebg);
205 significance = error ? (sig-bg)/error : 0;
206
207 cout << "Excess: " << sig-bg << endl;
208 cout << "N bkg: " << bg << endl;
209 cout << "Significance: " << significance << endl;
210
211 // save the result in file
212 fout << minsize << '\t'
213 << maxsize << '\t'
214 << width << '\t'
215 << length << '\t'
216 << sig-bg << '\t'
217 << significance << '\t' << endl;
218
219 delete onhisto;
220 delete offhisto;
221
222 width+=wdstep;
223 }
224 length+=lgstep;
225 }
226 }
227 fout.close();
228}
229//-----------------------------------------------------------------------
230
231Bool_t readDatacards(TString& filename)
232{
233 ifstream ifun(filename.Data());
234 if(!ifun)
235 {
236 cout << "File " << filename << " not found" << endl;
237 return kFALSE;
238 }
239
240 TString word;
241
242 while(ifun >> word)
243 {
244 // skip comments
245 if(word[0]=='/' && word[1]=='/')
246 {
247 while(ifun.get()!='\n'); // skip line
248 continue;
249 }
250
251 // number of events
252 if(strcmp(word.Data(),"NEVENTS")==0)
253 ifun >> nmaxevents;
254
255 // number of events
256 if(strcmp(word.Data(),"ONRATE")==0)
257 ifun >> onRate;
258
259 // input file name (on data)
260 if(strcmp(word.Data(),"ONFILES")==0)
261 {
262 if(onFile.Length())
263 cout << "readDataCards Warning: overriding on-data file name" << endl;
264 ifun >> onFile;
265 }
266
267 // input file name (off data)
268 if(strcmp(word.Data(),"OFFFILES")==0)
269 {
270 if(offFile.Length())
271 cout << "readDataCards Warning: overriding off-data file name" << endl;
272 ifun >> offFile;
273 }
274
275 // output file name
276 if(strcmp(word.Data(),"OUTFILE")==0)
277 {
278 if(outputFile.Length())
279 cout << "readDataCards Warning: overriding output file name" << endl;
280 ifun >> outputFile;
281 }
282
283 // size cuts
284 if(strcmp(word.Data(),"SIZECUTS")==0)
285 {
286 ifun >> lowsizecut;
287 ifun >> uppsizecut;
288 ifun >> nsizebins;
289 }
290
291 // dist cuts
292 if(strcmp(word.Data(),"DISTCUTS")==0)
293 {
294 ifun >> lowdistcut ;
295 ifun >> uppdistcut ;
296 }
297
298 // length cuts
299 if(strcmp(word.Data(),"LENGTHCUTS")==0)
300 {
301 ifun >> lowlgcut;
302 ifun >> upplgcut;
303 ifun >> lgstep;
304 }
305
306 // width cuts
307 if(strcmp(word.Data(),"WIDTHCUTS")==0)
308 {
309 ifun >> lowwdcut;
310 ifun >> uppwdcut;
311 ifun >> wdstep;
312 }
313
314 }
315
316 // check compulsory values
317 if(!onFile.Length())
318 {
319 cout << "No on-data file name specified" << endl;
320 return kFALSE;
321 }
322
323 if(!offFile.Length())
324 {
325 cout << "No off-data file name specified" << endl;
326 return kFALSE;
327 }
328
329 if(!outputFile.Length())
330 {
331 cout << "No output file name specified" << endl;
332 return kFALSE;
333 }
334
335 // Dump read values
336 cout << "************************************************" << endl;
337 cout << "* Datacards read from file " << filename << endl;
338 cout << "************************************************" << endl;
339 cout << "Maximum number of input events: " << nmaxevents << endl;
340 cout << "On-data acceptance rate: 1/" << onRate << endl;
341 cout << "On-data input file name(s): " << onFile << endl;
342 cout << "Off-data input file name(s): " << offFile << endl;
343 cout << "Output file name: " << outputFile << endl;
344 cout << "Dist cuts (deg): [" << lowdistcut<< ","<<uppdistcut<< "]" << endl;
345 cout << "Scanning size range: [" << lowsizecut << "," << uppsizecut << "]"<< " in " << nsizebins << " steps" << endl;
346 cout << "Scanning length range (deg): [" << lowlgcut<< ","<<upplgcut<< "], with step (deg): "<< lgstep << endl;
347 cout << "Scanning width range (deg): [" << lowwdcut<< ","<<uppwdcut<< "], with step (deg): "<< wdstep << endl;
348 cout << "***********" << endl << endl;
349
350 return kTRUE;
351}
Note: See TracBrowser for help on using the repository browser.