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

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