source: trunk/MagicSoft/Mars/mtemp/mifae/programs/falseSource.cc@ 3973

Last change on this file since 3973 was 3973, checked in by rico, 22 years ago
*** empty log message ***
File size: 13.0 KB
Line 
1////////////////////////////////////////////////////////////////////////////////////
2//
3// _____False Source Method macro_____
4//
5// Take as input root files with hillas parameters, perform derotation (optional)
6// and generates the alpha histograms for On and Off for different positions of the
7// source. These histos can be ploted with signal.C and signalPoint.C
8//
9// Ester Aliu <aliu@ifae.es>
10// Oscar Blanch <blanch@ifae.es>
11// Javier Rico <jrico@ifae.es>
12////////////////////////////////////////////////////////////////////////////////////
13
14#include <fstream>
15#include <iostream>
16
17#include "TString.h"
18#include "TH1F.h"
19#include "TFile.h"
20
21#include "MHillasSrcCalc.h"
22#include "MSrcPosCam.h"
23#include "MHillasSrc.h"
24#include "MSrcRotate.h"
25#include "MParList.h"
26#include "MTaskList.h"
27#include "MHillas.h"
28#include "MReadTree.h"
29#include "MEvtLoop.h"
30#include "MLog.h"
31#include "MArgs.h"
32
33using namespace std;
34
35Bool_t readDatacards(TString& filename);
36void falseSource();
37
38//-----------------------------------------------------------------------------
39// declaration of variables read from datacards
40//-----------------------------------------------------------------------------
41
42TString onFile;
43TString offFile;
44TString outputFile;
45Bool_t kRotate=1;
46Double_t fRA= -1.;
47Double_t fDEC= -1.;
48ULong_t nmaxevents=999999999;
49
50// cuts (default values to be overriden by datacards)
51Float_t minsize = 0; // minimum size (# of photons)
52Float_t maxsize = 9999999; // maximum size (# of photons)
53Float_t mindist = 0; // minimum distance cut (deg)
54Float_t maxdist = 10; // maximum distance cut (deg)
55Float_t minwidth = 0.; // minimum width cut (deg)
56Float_t maxwidth = 10; // maximum width cut (deg)
57Float_t minlength = 0; // minimum length cut (deg)
58Float_t maxlength = 10; // maximum length cut (deg)
59Float_t maxXY = maxdist; // maximum X and Y distance from the camera center
60
61// binning
62Int_t binning = 21; // number of bins in false source search (one coordinate)
63Int_t nbin = 18; // number of bins for alpha plots
64Float_t field = 2.; // width of the examined field (degrees)
65
66//-----------------------------------------------------------------------------
67// constants
68//-----------------------------------------------------------------------------
69
70const TString defaultcard="falsesource.datacard";
71
72const Float_t alphamin = 0.; // minimum value in alpha (degrees)
73const Float_t alphamax = 90.; // maximum value in alpha (degrees)
74const Float_t conver = 189./0.6; // conversion factor degrees to mm
75
76//-----------------------------------------------------------------------------
77
78static void Usage()
79{
80 gLog <<endl;
81 gLog << "Usage is:" << endl;
82 gLog << " falseSource [-h] [-?] <datacards>" << endl << endl;
83 gLog << " <datacards>: datacards file name (dafault falsesource.datacards)" << endl;
84 gLog << " -?/-h: This help" << endl << endl;
85}
86
87//-----------------------------------------------------------------------------
88int main(int argc, char **argv)
89{
90 // evaluate arguments
91 MArgs arg(argc, argv);
92 if (arg.HasOnly("-?") || arg.HasOnly("-h") || arg.HasOnly("--help"))
93 {
94 Usage();
95 return -1;
96 }
97
98 TString datacard = arg.GetArgumentStr(0);
99 if(!datacard.Length())
100 datacard = defaultcard;
101
102 if(!readDatacards(datacard))
103 {
104 cout << "Error reading datacards. Stoping" << endl;
105 return -1;
106 }
107 falseSource();
108}
109
110void falseSource()
111{
112 // CUTS
113
114 // variable declaration
115 Float_t xpos;
116 Float_t ypos;
117 Char_t name[20];
118 Char_t title[50];
119 Float_t stepDegrees = field/(binning-1); // degrees
120 Float_t step = stepDegrees*conver; // mm
121 Int_t ival = (Int_t) ceil((Float_t)(binning-1)/2); // index bound
122
123 // Hillas variables
124 Float_t alpha;
125 Float_t size;
126 Float_t dist;
127 Float_t length;
128 Float_t width;
129 Float_t meanX;
130 Float_t meanY;
131
132 // create the histograms (empty)
133 TH1F *hOnAlpha[binning][binning];
134 TH1F *hOffAlpha[binning][binning];
135 for(Int_t i = -ival; i <= ival ; i++){
136 for(Int_t j = -ival; j <= ival ; j++){
137
138 sprintf(name,"hOnAlpha[%d][%d]", i, j);
139 sprintf(title,"Alpha-Plot(On data) (%d ,%d)", i, j);
140 hOnAlpha[i+ival][j+ival] = new TH1F(name, title, nbin, alphamin, alphamax);
141
142 sprintf(name,"hOffAlpha[%d][%d]", i, j);
143 sprintf(title,"Alpha-Plot(Off data) (%d ,%d)", i, j);
144 hOffAlpha[i+ival][j+ival] = new TH1F(name, title, nbin, alphamin, alphamax);
145
146 }
147 }
148
149
150
151 /******************************************************************
152 FILL ON-DATA HISTOS
153 *******************************************************************/
154
155 // source dependent hillas containers/tasks
156 MHillasSrcCalc* csrc1[binning][binning];
157 MSrcPosCam* source[binning][binning];
158 MHillasSrc* hillasSrc[binning][binning];
159 MSrcRotate* srcrotate[binning][binning];
160
161 // normal containers/classes
162 MParList plist;
163 MTaskList tlist;
164 MHillas hillas;
165
166 plist.AddToList(&tlist);
167 plist.AddToList(&hillas);
168
169 MReadTree read("Parameters", onFile);
170 read.DisableAutoScheme();
171 tlist.AddToList(&read);
172
173 Int_t k,l;
174 char sourceName[100];
175 char hillasSrcName[100];
176
177 // create the tasks/containers for the different source positions
178 for(Int_t i=-ival;i<=ival;i++)
179 {
180 xpos = step*i;
181 for(Int_t j=-ival;j<=ival;j++)
182 {
183 ypos = step*j;
184
185 // name the different containers
186 if (i<0 && j<0)
187 {
188 k = -i;
189 l = -j;
190 sprintf(sourceName, "MSrcPosCam_Min%dMin%d", k, l);
191 sprintf(hillasSrcName, "MHillasSrc_Min%dMin%d", k, l);
192 }
193
194 if (i<0 && j>=0)
195 {
196 k = -i;
197 sprintf(sourceName, "MSrcPosCam_Min%d%d", k, j);
198 sprintf(hillasSrcName, "MHillasSrc_Min%d%d", k, j);
199 }
200
201 if (i>=0 && j<0)
202 {
203 l = -j;
204 sprintf(sourceName, "MSrcPosCam_%dMin%d", i, l);
205 sprintf(hillasSrcName, "MHillasSrc_%dMin%d", i, l);
206 }
207
208 if (i>=0 && j>= 0)
209 {
210 sprintf(sourceName, "MSrcPosCam_%d%d", i, j);
211 sprintf(hillasSrcName, "MHillasSrc_%d%d", i, j);
212 }
213
214 // include containers in parameter list
215 source[i+ival][j+ival] = new MSrcPosCam(sourceName);
216 source[i+ival][j+ival]->SetXY(xpos, ypos);
217 plist.AddToList(source[i+ival][j+ival]);
218 source[i+ival][j+ival]->Print();
219
220 hillasSrc[i+ival][j+ival] = new MHillasSrc(hillasSrcName);
221 plist.AddToList(hillasSrc[i+ival][j+ival]);
222
223 // define the different tasks and include them in the task list
224 if(kRotate)
225 {
226 srcrotate[i+ival][j+ival] = new MSrcRotate(sourceName);
227 srcrotate[i+ival][j+ival]->SetRAandDEC(fRA,fDEC);
228 tlist.AddToList(srcrotate[i+ival][j+ival]);
229 }
230
231 TString HilName = "MHillas";
232 csrc1[i+ival][j+ival] = new MHillasSrcCalc(sourceName, hillasSrcName);
233 csrc1[i+ival][j+ival]->SetInput(HilName);
234
235 tlist.AddToList(csrc1[i+ival][j+ival]);
236
237 } // loop on j (y coordinate)
238 } // loop on i (x coordinate)
239
240
241 // Eventloop
242 MEvtLoop evtloop;
243 evtloop.SetParList(&plist);
244 // MProgressBar bar;
245 // evtloop.SetProgressBar(&bar);
246 if (!evtloop.PreProcess())
247 return;
248
249 // select the events passing the cuts
250 UInt_t nread=0;
251 while (tlist.Process())
252 {
253 for(Int_t i = -ival; i <= ival ; i++)
254 for(Int_t j = -ival; j <= ival ; j++)
255 {
256 width = hillas.GetWidth()/conver;
257 length = hillas.GetLength()/conver;
258 size = hillas.GetSize();
259 dist = hillasSrc[i+ival][j+ival]->GetDist()/conver;
260 meanX = hillas.GetMeanX()/conver;
261 meanY = hillas.GetMeanY()/conver;
262
263 if (width<maxwidth && length<maxlength && dist>mindist && dist<maxdist && size>minsize && meanX<maxXY && meanY<maxXY)
264 {
265 alpha = hillasSrc[i+ival][j+ival]->GetAlpha();
266 hOnAlpha[i+ival][j+ival]->Fill(abs(alpha));
267 }
268 }
269 if(++nread>nmaxevents) break;
270 }
271
272 // end
273 evtloop.PostProcess();
274 tlist.PrintStatistics();
275
276
277 /******************************************************************
278 FILL OFF-DATA HISTOS
279 *******************************************************************/
280
281 MReadTree read2("Parameters", offFile);
282 read2.DisableAutoScheme();
283
284 tlist.AddToListBefore(&read2, &read, "All");
285 tlist.RemoveFromList(&read);
286
287 if (!evtloop.PreProcess())
288 return;
289
290 nread=0;
291 while (tlist.Process())
292 {
293 for(Int_t i = -ival; i <= ival ; i++)
294 for(Int_t j = -ival; j <= ival ; j++)
295 {
296 width = hillas.GetWidth()/conver;
297 length = hillas.GetLength()/conver;
298 size = hillas.GetSize();
299 dist = hillasSrc[i+ival][j+ival]->GetDist()/conver;
300 meanX = hillas.GetMeanX()/conver;
301 meanY = hillas.GetMeanY()/conver;
302
303 if (width<maxwidth && length<maxlength && dist>mindist && dist<maxdist && size>minsize && meanX<maxXY && meanY<maxXY)
304 {
305 alpha = hillasSrc[i+ival][j+ival]->GetAlpha();
306 hOffAlpha[i+ival][j+ival]->Fill(abs(alpha));
307 }
308 }
309 if(++nread>nmaxevents) break;
310 }
311
312 evtloop.PostProcess();
313
314 // Save results
315 TFile *f = new TFile(outputFile, "RECREATE");
316
317 for(Int_t i = -ival; i <= ival ; i++){
318 for(Int_t j = -ival; j <= ival ; j++){
319 hOnAlpha[i+ival][j+ival]->Write();
320 hOffAlpha[i+ival][j+ival]->Write();
321 }
322 }
323
324 f->Close();
325
326}
327//-----------------------------------------------------------------------
328
329Bool_t readDatacards(TString& filename)
330{
331 ifstream ifun(filename.Data());
332 if(!ifun)
333 {
334 cout << "File " << filename << " not found" << endl;
335 return kFALSE;
336 }
337
338 TString word;
339
340 while(ifun >> word)
341 {
342 // skip comments
343 if(word[0]=='/' && word[1]=='/')
344 {
345 while(ifun.get()!='\n'); // skip line
346 continue;
347 }
348
349 // number of events
350 if(strcmp(word.Data(),"NEVENTS")==0)
351 ifun >> nmaxevents;
352
353 // on-data file name
354 if(strcmp(word.Data(),"ONFILES")==0)
355 {
356 if(onFile.Length())
357 cout << "readDataCards Warning: overriding on-data file name" << endl;
358 ifun >> onFile;
359 }
360
361 // off-data file name
362 if(strcmp(word.Data(),"OFFFILES")==0)
363 {
364 if(offFile.Length())
365 cout << "readDataCards Warning: overriding off-data file name" << endl;
366 ifun >> offFile;
367 }
368
369 // output file name
370 if(strcmp(word.Data(),"OUTFILE")==0)
371 {
372 if(outputFile.Length())
373 cout << "readDataCards Warning: overriding output file name" << endl;
374 ifun >> outputFile;
375 }
376
377 // rotation flag
378 if(strcmp(word.Data(),"ROTFLAG")==0)
379 ifun >> kRotate;
380
381 // source celestial coordinates
382 if(strcmp(word.Data(),"SRCCOORDS")==0)
383 {
384 ifun >> fRA;
385 ifun >> fDEC;
386 }
387
388 // field width
389 if(strcmp(word.Data(),"FIELD")==0)
390 ifun >> field;
391
392 // binning
393 if(strcmp(word.Data(),"BINNING")==0)
394 ifun >> binning;
395
396
397 // Number of bins in alpha plots
398 if(strcmp(word.Data(),"NBIN")==0)
399 ifun >> nbin;
400
401
402
403
404 // size cut
405 if(strcmp(word.Data(),"SIZECUT")==0)
406 {
407 ifun >> minsize;
408 ifun >> maxsize;
409 }
410
411 // dist cut
412 if(strcmp(word.Data(),"DISTCUT")==0)
413 {
414 ifun >> mindist;
415 ifun >> maxdist;
416 }
417
418 // width cut
419 if(strcmp(word.Data(),"WIDTHCUT")==0)
420 {
421 ifun >> minwidth;
422 ifun >> maxwidth;
423 }
424
425 // length cut
426 if(strcmp(word.Data(),"LENGTHCUT")==0)
427 {
428 ifun >> minlength;
429 ifun >> maxlength;
430 }
431
432 // maxX and maxY upper cut
433 if(strcmp(word.Data(),"CENTERCUT")==0)
434 ifun >> maxXY;
435 }
436
437 // check compulsory values
438 if(!onFile.Length())
439 {
440 cout << "No on-data file name specified" << endl;
441 return kFALSE;
442 }
443 if(!offFile.Length())
444 {
445 cout << "No off-data file name specified" << endl;
446 return kFALSE;
447 }
448 if(!outputFile.Length())
449 {
450 cout << "No output file name specified" << endl;
451 return kFALSE;
452 }
453
454
455 // Dump read values
456 cout << "************************************************" << endl;
457 cout << "* Datacards read from file " << filename << endl;
458 cout << "************************************************" << endl;
459 cout << "Maximum number of input events: " << nmaxevents << endl;
460 cout << "On-data file name(s): " << onFile << endl;
461 cout << "Off-data file name(s): " << offFile << endl;
462 cout << "Output file name: " << outputFile << endl;
463 cout << "De-rotation flag " << kRotate << endl;
464 cout << "Source celestial coordiantes (rad): RA = " << fRA << ", DEC = " << fDEC << endl;
465 cout << "Field width (degrees): " << field << endl;
466 cout << "Field binning: " << binning << endl;
467 cout << "Number of alpha plot bins: " << nbin << endl;
468 cout << "Size cuts (# of photons): ("<<minsize<<","<<maxsize<<")" << endl;
469 cout << "Dist cuts (degrees): ("<<mindist<<","<<maxdist<<")" << endl;
470 cout << "Length cuts (degrees): ("<<minlength<<","<<maxlength<<")" << endl;
471 cout << "Width cuts (degrees): ("<<minwidth<<","<<maxwidth<<")" << endl;
472 cout << "maxX and maxY upper cut (degrees): " << maxXY << endl;
473 cout << "***********" << endl << endl;
474
475 return kTRUE;
476}
Note: See TracBrowser for help on using the repository browser.