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

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