source: trunk/MagicSoft/Mars/mtemp/mmpi/macros/gridloop.C@ 5934

Last change on this file since 5934 was 4102, checked in by mazin, 21 years ago
*** empty log message ***
File size: 14.3 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cherenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17! Author(s): Daniel Mazin, 05/2004 <mailto:mazin@imppmu.mpg.de>
18!
19! Copyright: MAGIC Software Development, 2000-2004
20!
21!
22\* ======================================================================== */
23
24// Daniel Mazin 17.05.2004 mazin@mppmu.mpg.de
25// **********************************************************************************
26// this macro is used to produce false source plots
27// and it is a part from the macro falsesourcemethod.C.
28// See falsesourcemethod.C for more explanation
29// input: hillas parameters as vectors
30// output: several 2D plots in a root file +
31// ascii file with parameters like significance, Nexcess etc.
32// **********************************************************************************
33
34
35void gridloop()
36{
37
38 printf("inside of gridloop() ... \n");
39
40/* emty histograms */
41
42 histskyplot.Reset();
43 histNexOnOff.Reset();
44 histChi2Off.Reset();
45 histskyplotOnOff.Reset();
46 histNex.Reset();
47 histChi2.Reset();
48 histskyplot4.Reset();
49 histskyLiMa.Reset();
50 histsign.Reset();
51 histsign4.Reset();
52 histBerlin.Reset();
53
54/* end emty histograms */
55
56
57 TH1F histalphagrid;
58 histalphagrid.SetName("Alpha");
59 histalphagrid.SetTitle("Alpha");
60 histalphagrid.SetXTitle("Alpha (deg)");
61 histalphagrid.SetYTitle("Counts");
62 histalphagrid.SetDirectory(NULL);
63 histalphagrid.SetFillStyle(4000);
64 histalphagrid.UseCurrentStyle();
65
66 bins.SetEdges(18, 0.0, 90.);
67 bins.Apply(histalphagrid);
68
69 Double_t xcorr, ycorr, xtest, ytest, xtilde, ytilde, fconc1;
70
71 const Float_t binwidth1 = histalphagrid.GetBinWidth(1);
72 if (TMath::Abs( ALPHAMAX/binwidth1 - TMath::Nint(ALPHAMAX/binwidth1)) > TOLERANCE)
73 {
74 cout << "wrong binning for histogram: histalphagrid \n SORRY ... ABORT" << endl;
75 exit(1);
76 }
77 const Int_t numbbinOnMax = TMath::Nint(ALPHAMAX / binwidth1); // number bins in ON region
78
79 /* fit functions */
80 TF1 * fitsig = new TF1("fsig", "[0]", 0., ALPHAMAX);
81 fitsig->SetLineColor(4);
82 TF1 * fitbgpar = new TF1("fbgpar", "[0]*x*x + [1]", 30., 90.);
83 fitbgpar->SetLineColor(2);
84
85 TF1 * fitbg4 = new TF1("fbg4", "[0]*x*x*x*x + [1]*x*x + [2]", 30., 90.);
86 fitbg4->SetLineColor(3);
87 /* end fit functions */
88
89 Float_t Nex, Non, Noff, Sign;
90 Float_t xsource, ysource, xsournew, ysournew;
91
92 Double_t Nex4, Noff4, Sign4;
93 Double_t alpha, SignLiMa;
94
95 Double_t chisquarefit;
96
97 Double_t integon, integoff, normf, NoffOFF, SignOnOff, NexOnOff;
98 Float_t logsize, lgsize, lgsize2, tanbeta, beta;
99 const Float_t LOG3000 = log(3000.);
100
101
102
103 Double_t chi2par[GRIDBINS], apar[GRIDBINS], bpar[GRIDBINS]; // fit parameters from OFF sample
104 Float_t offbin1[GRIDBINS], offbin2[GRIDBINS], offbin3[GRIDBINS]; // Number of events in first 3 bins
105 Double_t xpos, ypos, binwidth2;
106 Int_t numgrid;
107
108
109 Int_t numbinsoff = 12; // FIXME: 12 should be changed according to the binning
110 Int_t nbinx, nbiny, gridpoint = 0;
111
112 Char_t psname[120], titelname[120], plot1name[120], plot2name[120], stringsig[120];
113
114 FILE * fp;
115
116 if (TYPEOPTION == kTRUE) // read fit parameters from the OFF sample
117 {
118 Char_t * datin = DATPARAMIN;
119
120 fp = fopen(datin, "r");
121
122 for (Int_t i=0; i<GRIDBINS; i++)
123 {
124 fscanf(fp,"%d %lf %lf %lf %lf %lf %lf %f %f %f",
125 &numgrid, &xpos, &ypos, &apar[i], &bpar[i],
126 &chi2par[i], &binwidth2,
127 &offbin1[i], &offbin2[i], &offbin3[i]);
128// cout << "numgrid: " << numgrid << ", apar " << apar[i] << endl;
129 }
130
131 fclose(fp);
132 }
133
134
135 Char_t * chout = DATPARAMOUT; // write parameters into file
136
137 fp = fopen(chout, "w");
138
139
140/* begin of the loop over all bins in the sky map */
141
142 for(Int_t gridx = 0; gridx < NUMSTEPS; gridx++)
143 {
144 xsource = MINXGRID + gridx * STEPGRID;
145//cout << xsource << endl;
146
147 for(Int_t gridy = 0; gridy < NUMSTEPS; gridy++)
148 {
149 ysource = MINYGRID + gridy * STEPGRID;
150//cout << ysource << endl;
151
152// empty histogram:
153// nbinx = histalphagrid.GetNbinsX();
154// for(Int_t ii = 0; ii< (nbinx+1); ii++) histalphagrid.SetBinContent(ii, 0.);
155// histalphagrid.SetEntries(0);
156 histalphagrid.Reset();
157
158
159// Bool_t anglephi = kFALSE;
160
161 /* loop over all events */
162 for(Int_t k=0; k< imnum; k++)
163 {
164
165 if (ROTOPTION == kTRUE) // derotate into sky coordinates
166 {
167 /* derotation : correct sky coordinates into camera coordinates */
168 xsournew = cosgam(k) * xsource - singam(k) * ysource;
169 ysournew = singam(k) * xsource + cosgam(k) * ysource;
170 /* end derotatiom */
171 }
172 else // do not derotate, plot into camera coordinates
173 {
174 xsournew = xsource;
175 ysournew = ysource;
176 }
177
178
179 /* calculate ALPHA und DIST acording to the source position */
180 tanbeta = (ypar(k) - ysournew) / (xpar(k) - xsournew);
181 beta = TMath::ATan(tanbeta);
182 alphapar(k) = (deltapar(k) - beta) * 180./ TMath::Pi();
183 distpar(k) = sqrt((ypar(k) - ysournew) * (ypar(k) - ysournew) +
184 (xpar(k) - xsournew) * (xpar(k) - xsournew));
185
186
187
188 if(alphapar(k) > 90.) alphapar(k) -= 180.;
189 if(alphapar(k) < -90.) alphapar(k) += 180.;
190
191// histalpha.Fill(alphapar(k),1.);
192// histdist.Fill( distpar(k),1.);
193
194 if(CUTOPTION==kFALSE)
195 {
196/* ****************************************************** */
197 /* static cuts */
198 if (lengthpar(k) > LENGTHMIN && lengthpar(k) < LENGTHMAX)
199 if (widthpar(k) > WIDTHMIN && widthpar(k) < WIDTHMAX)
200 if (sizepar(k) > SIZEMIN)
201 if (distpar(k) > DISTMIN && distpar(k) < DISTMAX)
202 if (leak1(k) < LEAKMAX)
203/* ****************************************************** */
204 {
205 alphapar(k) = TMath::Abs(alphapar(k));
206//printf("alpha passed the cuts: %5.2f \n", alphapar(k));
207//printf("histalphagrid: %5.2f \n", histalphagrid.GetNbinsX());
208 histalphagrid.Fill(alphapar(k),1.);
209
210 }
211
212 }
213/* ****************************************************** */
214/* dynamical cuts */
215 else // if(CUTOPTION==kTRUE)
216 {
217 logsize = log(sizepar(k));
218 lgsize = logsize-LOG3000;
219 lgsize2 = lgsize*lgsize;
220 if ( sizepar(k) > SIZEMIN )
221 if ( leak1(k) < LEAKMAX )
222 if ( lengthpar(k) > (LENGTHMINParA + LENGTHMINParB*lgsize + LENGTHMINParC*lgsize2) &&
223 lengthpar(k) < (LENGTHMAXParA + LENGTHMAXParB*lgsize + LENGTHMAXParC*lgsize2))
224 if ( widthpar(k) > (WIDTHMINParA + WIDTHMINParB*lgsize + WIDTHMINParC*lgsize2) &&
225 widthpar(k) < (WIDTHMAXParA + WIDTHMAXParB*lgsize + WIDTHMAXParC*lgsize2) )
226 if ( distpar(k) > (DISTMINParA + DISTMINParB*lgsize + DISTMINParC*lgsize2) &&
227 distpar(k) < (DISTMAXParA + DISTMAXParB*lgsize + DISTMAXParC*lgsize2) )
228 {
229 alphapar(k) = TMath::Abs(alphapar(k));
230 histalphagrid.Fill(alphapar(k),1.);
231 if(alphapar(k) < 8.) histBerlin.Fill(xsournew, ysournew);
232 }
233 }
234
235 }
236// histalphagrid.Fit("fsig","NR");
237//printf("content bin 1: %4.2f\n", histalphagrid.GetBinContent(1));
238 histalphagrid.Fit(fitsig,"NWR");
239 histalphagrid.Fit("fbg4","+NWR");
240 histalphagrid.Fit("fbgpar","+NWR");
241
242// cout << "OK" << endl;
243// calculate significance:
244 Non = 0.;
245 for(Int_t i=1; i<=numbbinOnMax;i++) Non += histalphagrid.GetBinContent(i);
246// Non = (fitsig->GetParameter(0)) * numbbinOnMax;
247//cout << "Non : " << Non << endl;
248// histalphagrid.Draw();
249
250 chisquarefit = fitbgpar->GetChisquare();
251 Noff = (1./3. * (fitbgpar->GetParameter(0)) * pow(ALPHAMAX,3.) +
252 (fitbgpar->GetParameter(1)) * ALPHAMAX) / binwidth1;
253
254//cout << "Noff : " << Noff << endl;
255
256 Nex = Non - Noff;
257
258 Noff4 = (1./5. * fitbg4->GetParameter(0) * pow(ALPHAMAX,5.) +
259 1./3. * (fitbg4->GetParameter(1)) * pow(ALPHAMAX,3.) +
260 fitbg4->GetParameter(2) * ALPHAMAX ) / binwidth1;
261
262//cout << "Noff4 : " << Noff4 << endl;
263
264 Nex4 = Non - Noff4;
265
266// calculate significance:
267
268// Sign4 = Nex4 / sqrt(Nex4 + 2.* Noff4);
269 if (Noff4<0.) Sign4 = 0.;
270 else Sign4 = LiMa17(Non,Noff4,1.);
271
272//cout << "Sign4 : " << Sign4 << endl;
273// Sign = Nex / sqrt(Nex + 2.* Noff);
274 if (Noff<0.) Sign = 0.;
275 else Sign = LiMa17(Non,Noff,1.);
276//cout << "Sign LiMa : " << Sign << endl;
277
278 /* calculate Noff from the OFF data */
279 if (TYPEOPTION == kTRUE)
280 {
281 // ON:
282
283 integon = 0.; // number of events between 30 and 90 degrees
284
285 for (Int_t ik = 0; ik < numbinsoff; ik++)
286 {
287 integon += histalphagrid.GetBinContent(7+ik); // FIXME: 7 should be changed according to the binning
288 }
289 // OFF:
290
291 integoff = ((1./3. * apar[gridpoint] * pow(90.,3.) + bpar[gridpoint] * 90.) -
292 (1./3. * apar[gridpoint] * pow(30.,3.) + bpar[gridpoint] * 30.)) / binwidth2;
293
294 normf = integoff / integon;
295
296 NoffOFF = (1./3. * apar[gridpoint] * pow(ALPHAMAX,3.) +
297 (bpar[gridpoint] * ALPHAMAX)) / binwidth2 / normf;
298
299 NexOnOff = Non - NoffOFF;
300
301 SignOnOff = NexOnOff / sqrt(NexOnOff + 2.* NoffOFF);
302
303 // calculate according to Li Ma:
304 if (NoffOFF<0.) SignLiMa=0.;
305 else SignLiMa = LiMa17(Non,NoffOFF*normf,1./normf);
306
307 }
308
309 cout << endl << " xsource: " << xsource << " ysource: " << ysource << endl;
310
311 cout << " Non : " << Non << " Noff : " << Noff << " Nex : " << Nex << endl;
312
313 if (TYPEOPTION == kTRUE)
314 cout << " ON and OFF : Non : " << Non << " Noff (estimated) : " << NoffOFF << " Nex : " << NexOnOff << endl;
315
316 cout << endl;
317 cout << " significance (parab): " << Sign << " sigma" << endl;
318 cout << " significance (4th order): " << Sign4 << " sigma" << endl;
319 if (TYPEOPTION == kTRUE)
320 {
321 cout << " significance (using ON and OFF): " << SignOnOff << " sigma" << endl;
322 cout << " significance (LiMa 17): " << SignLiMa << " sigma" << endl;
323 }
324
325
326 if (TYPEOPTION == kTRUE)
327 fprintf(fp,"%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf \n",
328 xsource, ysource, Sign, Non, Noff, Nex, SignLiMa,
329 SignOnOff, NoffOFF, NexOnOff, chi2par[gridpoint]);
330 else
331 fprintf(fp,"%d %lf %lf %lf %lf %lf %lf %lf %lf %lf\n",
332 gridpoint+1, xsource, ysource,
333 fitbgpar->GetParameter(0), fitbgpar->GetParameter(1),
334 fitbgpar->GetChisquare(), histalphagrid.GetBinWidth(1),
335 histalphagrid.GetBinContent(1), histalphagrid.GetBinContent(2),
336 histalphagrid.GetBinContent(3));
337
338
339 histskyplot.SetBinContent(gridx+1, gridy+1, Sign);
340 histNex.SetBinContent(gridx+1, gridy+1, Nex);
341 histChi2.SetBinContent(gridx+1, gridy+1, fitbgpar->GetChisquare() / 10.);
342
343 if (TYPEOPTION == kTRUE) histskyplotOnOff.SetBinContent(gridx+1, gridy+1, SignOnOff);
344 if (TYPEOPTION == kTRUE) histNexOnOff.SetBinContent(gridx+1, gridy+1, NexOnOff);
345 if (TYPEOPTION == kTRUE) histChi2Off.SetBinContent(gridx+1, gridy+1, chi2par[gridpoint] / 16.);
346
347 if (TYPEOPTION == kTRUE) histskyLiMa.SetBinContent(gridx+1, gridy+1, SignLiMa);
348 histskyplot4.SetBinContent(gridx+1, gridy+1, Sign4);
349 histsign.Fill(Sign,1.);
350 histsign4.Fill(Sign4,1.);
351// break;
352
353 if (TYPEOPTION == kTRUE)
354 {
355 TF1 * bgoff = new TF1("bgoff", parabfunc, 0., 90., 3);
356 bgoff->SetParameters(apar[gridpoint], bpar[gridpoint], normf);
357 bgoff->FixParameter(0, apar[gridpoint]);
358 bgoff->FixParameter(1, bpar[gridpoint]);
359 bgoff->FixParameter(2, normf);
360 bgoff->SetLineColor(8);
361 }
362
363// sprintf(stringsig,"S = %.2f sigma", Sign);
364 if (TYPEOPTION == kTRUE) sprintf(stringsig,"S = %.2f sigma", SignOnOff);
365/*
366 sprintf(psname,"/.magic/magicserv01/scratch/Daniel/plots/RootPlots/040215Mrk421_misp/all/alphaX%.2fY%.2fON.root", xsource, ysource);
367 sprintf(titelname,"alpha plot. assumed source position: x = %.2f y = %.2f", xsource, ysource);
368 TCanvas * cplot = new TCanvas();
369 histalphagrid.SetTitle(titelname);
370 histalphagrid.Draw();
371 bgoff->Draw("same");
372 leg = new TLegend(0.1,0.15,0.52,0.35);
373 leg->Draw();
374 leg->AddEntry(fitsig,"fit for ON region (0-20 deg)","l");
375 leg->AddEntry(fitbgpar,"fit for OFF region (30-90 deg)","l");
376 leg->AddEntry(bgoff,"fit from OFF data (0-90 deg)","l");
377 leg->SetHeader("Legend");
378 leg->SetFillColor(19);
379 leg->Draw();
380
381 text = new TPaveText(50,20,90,80);
382 text->AddText(0.5, 0.8, "calculated significance");
383 text->AddText(0.5, 0.5, "using ON and OFF:");
384 text->AddText(1, 0.2, stringsig);
385 text->SetTextSize(0.04);
386 text->Draw();
387
388 cplot->Modified();
389 cplot->Update();
390
391 cplot->SaveAs(psname);
392 delete cplot;
393*/
394 if (TYPEOPTION == kTRUE) delete bgoff;
395
396 gridpoint++;
397
398 }
399 }
400
401 /* close file */
402 fclose(fp);
403
404
405 TFile rootfile(ROOTPLOTNAME, "RECREATE",
406 "sky plots in all variations");
407
408
409// cout << " end of one of the grids 1 " << endl;
410
411 histskyplot.Write();
412
413 histNex.Write();
414
415 histChi2.Write();
416
417 histsign.Write();
418
419 histskyplot4.Write();
420
421 histsign4.Write();
422
423 histBerlin.Write();
424
425 if (TYPEOPTION == kTRUE)
426 {
427 histNexOnOff.Write();
428 histChi2Off.Write();
429 histskyplotOnOff.Write();
430 histskyLiMa.Write();
431 }
432
433 rootfile.Close();
434
435 delete fitsig;
436 delete fitbgpar;
437 delete fitbg4;
438
439 cout << " end of one of the grids " << endl;
440}
441
Note: See TracBrowser for help on using the repository browser.