source: trunk/MagicSoft/Mars/mtemp/mpadova/macros/wobblemap.C@ 6668

Last change on this file since 6668 was 6396, checked in by moralejo, 20 years ago
*** empty log message ***
File size: 29.4 KB
Line 
1void wobblemap(Float_t maxhadronness = 0.15, Float_t minsize=2000., Float_t maxsize = 1.e10)
2{
3
4 //// CHANGE NAME OF FILE TO NOT OVERWRITE FORMER FILES!!! ////////
5 TString outname = "Hists_";
6 if (maxsize > 1.e6)
7 outname += "above";
8
9 outname += Form("%4d", (Int_t)minsize);
10
11 if (maxsize < 1.e6)
12 {
13 outname += "to";
14 outname += Form("%4d", (Int_t)maxsize);
15 }
16
17 outname += "_had";
18
19 TString dummy = Form("%3.2f", (Double_t)maxhadronness); // cast is necessary!!
20 outname += dummy;
21 outname += ".root";
22
23 cout << outname << endl;
24
25 TFile* hfile = new TFile(outname,"recreate");
26
27 MGeomCamMagic mgeom;
28 Double_t mm2deg = mgeom.GetConvMm2Deg();
29
30 Double_t CrabDec = 22.014444444; // deg
31 Double_t CrabRA = 5.5755555555; // hours
32
33// Double_t CrabDec = 19.8197221; // deg
34// Double_t CrabRA = 8.895194; // hours // FALSE source for the 3EG0853 data!
35
36 Float_t onregion_radius = 0.25; // degree
37 Float_t onr2 = onregion_radius * onregion_radius; // deg2
38 Float_t datatime = 0.;
39
40 // CHANGE MC PARAMETERS AND FILE BELOW ACCORDING TO STUDIED DATA SETS!!!!
41 // Period 025:
42// TString spotsize = "2.0cm";
43// Float_t minMCtheta = 5.73; // lower edge of zbin 1
44// Float_t maxMCtheta = 23.79; // upper edge of zbin 8
45// Float_t n_gamma_produced = 1.e5 * 8. * 6. / 2.;
46 // #Events per file * #zbins * #files per zbin / 2 (test sample only)
47
48 // YOU ALSO HAVE TO CHANGE THE SCALING OF W, L AND THE EXPRESSION FOR DISP!
49
50 // Period 021
51 TString spotsize = "1.4cm";
52 Float_t minMCtheta = 15.20; // lower edge of zbin 4
53 Float_t maxMCtheta = 28.96; // upper edge of zbin 12
54 Float_t n_gamma_produced = 1.e5 * 9. * 6. / 2.;
55 // #Events per file * #zbins * #files per zbin / 2 (test sample only)
56
57
58
59 TChain *MCchain = new TChain("Events");
60 // MCchain->Add("/data2/magic/data/rootdata/Crab/Period025/2005_01_10/RandomForest_1.0_spot/star_S1000_mcgammas.root");
61
62 if (spotsize == "1.4cm")
63 MCchain->Add("/data2/magic/data/rootdata/Crab/Period021/2004_09_15/RandomForest_1.4_spot/star_S1000_mcgammas.root");
64 else if (spotsize == "2.0cm")
65 MCchain->Add("/data2/magic/data/rootdata/Crab/Period025/2005_01_10/RandomForest_2.0_spot/star_S1000_mcgammas.root");
66
67 //////////////////////////////////////////////////////////////////////////////////////////////////
68
69 TChain *Wchain = new TChain("Events");
70
71 Wchain->Add("/data2/magic/data/rootdata/Crab/Period021/2004_09_15/RandomForest_1.4_spot/star_S1000_20040915_CrabNebulaW1.root");
72 datatime += 1062.3;
73
74 Wchain->Add("/data2/magic/data/rootdata/Crab/Period021/2004_09_24/RandomForest_1.4_spot/star_S1000_20040924_CrabNebulaW1.root");
75 datatime += 1316.2;
76
77 Wchain->Add("/data2/magic/data/rootdata/Crab/Period021/2004_09_15/RandomForest_1.4_spot/star_S1000_20040915_CrabNebulaW2.root");
78 datatime += 1241.1;
79
80 Wchain->Add("/data2/magic/data/rootdata/Crab/Period021/2004_09_18/RandomForest_1.4_spot/star_S1000_20040918_CrabNebulaW2.root");
81 datatime += 1388.1;
82
83 Wchain->Add("/data2/magic/data/rootdata/Crab/Period021/2004_09_24/RandomForest_1.4_spot/star_S1000_20040924_CrabNebulaW2.root");
84 datatime += 767.0;
85
86
87// Wchain->Add("/data2/magic/data/rootdata/Crab/Period025/2005_01_10/RandomForest_2.0_spot/star_S1000_20050110_CrabNebulaW1.root");
88// datatime += 9469.7;
89
90// Wchain->Add("/data2/magic/data/rootdata/Crab/Period025/2005_01_11/RandomForest_2.0_spot/star_S1000_20050111_CrabNebulaW1.root");
91// datatime += 4400.1;
92
93// Wchain->Add("/data2/magic/data/rootdata/Crab/Period025/2005_01_11/RandomForest_2.0_spot/star_S1000_20050111_CrabNebulaW2.root");
94// datatime += 4673.1;
95
96
97 // November data: WE SKIP THEM, they are very few!
98 // Wchain->Add("/data2/magic/data/rootdata/Crab/Period023/2004_11_23/RandomForest_1.0_spot/star_S1000_20041123_CrabNebulaW1.root");
99 // datatime += 394.2;
100 // Wchain->Add("/data2/magic/data/rootdata/Crab/Period023/2004_11_23/RandomForest_1.0_spot/star_S1000_20041123_CrabNebulaW2.root");
101 // datatime += 0.0; // No data with theta below 28.96 deg!
102
103
104// Wchain->Add("/data2/magic/data/rootdata/3EG0853/Period024/2004_12_10/RandomForest_2.0_spot/star_S1000_20041210_3EG0853.root");
105// datatime += 4487.9;
106
107
108
109 TH2F *CamMap = new TH2F("CamMap", "Events incoming direction map ",
110 50, -2.5, 2.5,50, -2.5, 2.5);
111
112
113 TH2F *SkyMap = new TH2F("SkyMap", "Events incoming direction map ",
114 70, CrabRA-0.175, CrabRA+0.175,
115 70, CrabDec-2.8315, CrabDec+2.8315);
116
117 TH2F *ONandOFF = new TH2F("ONandOFF", "ON and OFF regions",
118 140, CrabRA-0.175, CrabRA+0.175,
119 140, CrabDec-2.8315, CrabDec+2.8315);
120
121 ////////////// BACKGROUND MODEL: //////////////////////////////////////////////////
122
123 // Radial distribution, dN/dr with trespect to camera center:
124 // SIZE > 2000
125 TF1* myfunc = new TF1("myfunc","x*([0]+[1]*log(x)+[2]*log(x)*log(x))", 0., 2.);
126 myfunc->SetParameters(2.06765e+02, -2.36253e+02, -2.13230e+01);
127
128 // SIZE 1000-2000
129 // TF1* myfunc = new TF1("myfunc","x*([0]+[1]*log(x)+[2]*log(x)*log(x)+[3]*log(x)**4)", 0., 2.);
130 // myfunc->SetParameters(6.45011e+02, -1.01606e+03, 1.63394e+02, -3.60749e+01);
131
132 TH1F *rndmrad = new TH1F("rndmrad","",200, 0., 2.);
133 rndmrad->FillRandom("myfunc",1e6);
134 rndmrad->Draw("e");
135
136 TH2F *BgMap = (TH2F*) SkyMap->Clone("BgMap");
137 TRandom dice;
138
139 ///////////////////////////////////////////////////////////////////////////////////
140
141
142 TH1F *Alpha = new TH1F("Alpha", "", 60, -90., 90.);
143 TH1F *AbsAlpha = new TH1F("AbsAlpha", "", 45, 0., 90.);
144
145 TH1F *WidthOn = new TH1F("WidthOn", "WIDTH", 120, 0., 0.60);
146 TH1F *WidthOff = (TH1F*) WidthOn->Clone("WidthOff");
147 WidthOn->Sumw2();
148 WidthOff->Sumw2();
149
150 TH1F *LengthOn = new TH1F("LengthOn", "LENGTH", 50, 0., 1.);
151 TH1F *LengthOff = (TH1F*) LengthOn->Clone("LengthOff");
152 LengthOn->Sumw2();
153 LengthOff->Sumw2();
154
155 TH1F *SWidthOn = new TH1F("SWidthOn", "SWIDTH", 75, 0.3, 5.3);
156 TH1F *SWidthOff = (TH1F*) SWidthOn->Clone("SWidthOff");
157 SWidthOn->Sumw2();
158 SWidthOff->Sumw2();
159
160 TH1F *SLengthOn = new TH1F("SLengthOn", "SLENGTH", 45, 0.2, 3.2);
161 TH1F *SLengthOff = (TH1F*) SLengthOn->Clone("SLengthOff");
162 SLengthOn->Sumw2();
163 SLengthOff->Sumw2();
164
165 TH1F *ConcOn = new TH1F("ConcOn", "CONC", 30, 0., 0.6);
166 TH1F *ConcOff = (TH1F*) ConcOn->Clone("ConcOff");
167 ConcOn->Sumw2();
168 ConcOff->Sumw2();
169
170 TH1F *M3LongOn = new TH1F("M3LongOn", "M3LONG", 80, -0.8, 0.8);
171 TH1F *M3LongOff = (TH1F*) M3LongOn->Clone("M3LongOff");
172 M3LongOn->Sumw2();
173 M3LongOff->Sumw2();
174
175 TH1F *DistOn = new TH1F("DistOn", "DIST", 40, 0., 2.);
176 TH1F *DistOff = (TH1F*) DistOn->Clone("DistOff");
177 DistOn->Sumw2();
178 DistOff->Sumw2();
179
180 TH1F *SizeOn = new TH1F("SizeOn", "LOG10SIZE", 60, 2.5, 5.5);
181 TH1F *SizeOff = (TH1F*) SizeOn->Clone("SizeOff");
182 SizeOn->Sumw2();
183 SizeOff->Sumw2();
184
185
186 TH1F *HadronnessOn = new TH1F("HadronnessOn", "Hadronness", 101, 0., 1.01);
187 TH1F *HadronnessOff = (TH1F*) HadronnessOn->Clone("HadronnessOff");
188 HadronnessOn->Sumw2();
189 HadronnessOff->Sumw2();
190
191 TH1F* r_source = new TH1F("r_source", "r_source", 50, 0., 2.5);
192
193
194 // Now the MC histograms:
195
196 TH1F *AlphaMC = (TH1F*) Alpha->Clone("AlphaMC");
197 TH1F *AbsAlphaMC = (TH1F*) AbsAlpha->Clone("AbsAlphaMC");
198 TH1F *WidthMC = (TH1F*) WidthOn->Clone("WidthMC");
199 TH1F *LengthMC = (TH1F*) LengthOn->Clone("LengthMC");
200 TH1F *SWidthMC = (TH1F*) SWidthOn->Clone("SWidthMC");
201 TH1F *SLengthMC = (TH1F*) SLengthOn->Clone("SLengthMC");
202
203 TH1F *ConcMC = (TH1F*) ConcOn->Clone("ConcMC");
204 TH1F *M3LongMC = (TH1F*) M3LongOn->Clone("M3LongMC");
205 TH1F *M3LongMC_all = (TH1F*) M3LongOn->Clone("M3LongMC_all");
206 TH1F *DistMC = (TH1F*) DistOn->Clone("DistMC");
207 TH1F *SizeMC = (TH1F*) SizeOn->Clone("SizeMC");
208
209
210 TH1F *HadronnessMC = (TH1F*) HadronnessOn->Clone("HadronnessMC");
211
212 TH1F *EnergyMC = new TH1F("EnergyMC", "MC, log_{10}E_{\gamma} (GeV)", 1000, 0., 2000.);
213 EnergyMC->Sumw2();
214
215 ////////////////////////////////////////////////////////////////
216
217 MObservatory mobserv;
218 MStarCamTrans* starcamtrans = new MStarCamTrans(mgeom, mobserv);
219
220
221 MHillas *mhil = new MHillas();
222 MHillasSrc *msrc = new MHillasSrc();
223 MHillasExt *mhilext = new MHillasExt();
224 MNewImagePar *mnew = new MNewImagePar();
225 MPointingPos *mpoint = new MPointingPos();
226 MHadronness *mhad = new MHadronness();
227 MPointingPos *mpoint = new MPointingPos();
228 MTime *mtime = new MTime();
229
230 Wchain->SetBranchStatus("*", 0);
231 Wchain->SetBranchStatus("MHillas.fSize", 1);
232 Wchain->SetBranchStatus("MHillas.fWidth", 1);
233 Wchain->SetBranchStatus("MHillas.fLength", 1);
234 Wchain->SetBranchStatus("MHillas.fMeanX", 1);
235 Wchain->SetBranchStatus("MHillas.fMeanY", 1);
236 Wchain->SetBranchStatus("MHillas.fCosDelta", 1);
237 Wchain->SetBranchStatus("MHillas.fSinDelta", 1);
238 Wchain->SetBranchStatus("MHillasSrc.fCosDeltaAlpha", 1);
239 Wchain->SetBranchStatus("MHillasExt.fM3Long", 1);
240 Wchain->SetBranchStatus("MNewImagePar.fNumCorePixels", 1);
241 Wchain->SetBranchStatus("MNewImagePar.fConc", 1);
242 Wchain->SetBranchStatus("MPointingPos.fZd",1);
243 Wchain->SetBranchStatus("MHadronness.fHadronness",1);
244 Wchain->SetBranchStatus("MPointingPos.fDec",1);
245 Wchain->SetBranchStatus("MPointingPos.fRa",1);
246 Wchain->SetBranchStatus("MTime.*",1);
247 Wchain->SetBranchStatus("MHillasSrc.fDist",1);
248
249 Wchain->SetBranchAddress("MHillas.", &mhil);
250 Wchain->SetBranchAddress("MHillasSrc.", &msrc);
251 Wchain->SetBranchAddress("MHillasExt.", &mhilext);
252 Wchain->SetBranchAddress("MNewImagePar.", &mnew);
253 Wchain->SetBranchAddress("MPointingPos.", &mpoint);
254 Wchain->SetBranchAddress("MHadronness.", &mhad);
255 Wchain->SetBranchAddress("MTime.", &mtime);
256
257
258 for (Int_t i = 0; i < Wchain->GetEntries(); i++) // event count starts in 0
259 {
260 if ((i+1)%10000 == 0)
261 cout << i+1 << " of " << Wchain->GetEntries() << endl;
262
263 Wchain->GetEvent(i);
264
265 Float_t size = mhil->GetSize();
266
267 if (size < minsize)
268 continue;
269 else if (size > maxsize)
270 continue;
271
272 Float_t log10size = log10(size);
273
274 Float_t ZA = mpoint.GetZd();
275 // Exclude events with no drive info. For Crab ZAmin = 6.8
276 if (ZA < minMCtheta)
277 continue;
278 else if (ZA > maxMCtheta)
279 continue; // MC up to ~29 deg loaded, but data may have lower upper end!
280
281
282 Float_t width = mhil->GetWidth()* mm2deg;;
283 Float_t length = mhil->GetLength()* mm2deg;;
284
285
286 Float_t disp;
287
288 // spot 1 cm:
289 // Float_t disp = (1.-(width/length))*(2.68576-0.830937*log10size+0.1227*log10size*log10size);
290
291
292 if (spotsize == "1.4cm")
293 disp = (1.-(width/length))*(3.64292-1.28713*log10size+0.179722*log10size*log10size);
294
295 else if (spotsize == "2.0cm")
296 disp = (1.-(width/length))*(4.11949-1.39848*log10size+0.183514*log10size*log10size);
297
298
299 Float_t xmean = mhil->GetMeanX()* mm2deg;;
300 Float_t ymean = mhil->GetMeanY()* mm2deg;;
301
302 Float_t cosdelta = mhil->GetCosDelta();
303 Float_t sindelta = mhil->GetSinDelta();
304
305 Float_t asylon = mhilext->GetM3Long();
306
307
308 Double_t sourcex, sourcey;
309
310 // Choose right orientation:
311
312 if (asylon < 0.)
313 {
314 sourcex = xmean+cosdelta*disp;
315 sourcey = ymean+sindelta*disp;
316 }
317 else
318 {
319 sourcex = xmean-cosdelta*disp;
320 sourcey = ymean-sindelta*disp;
321 }
322
323 Double_t SourceDec = 0;
324 Double_t SourceRA = 0.;
325 Double_t SourceHA = 0.;
326
327 Double_t CenterDec = mpoint->GetDec();
328 Double_t CenterRA = mpoint->GetRa();
329
330 Double_t gmst = mtime->GetGmst(); // Greenwich mean sidereal time [rad]
331
332 gmst *= 24. / (2.*TMath::Pi()); // Conversion to hours
333
334 Double_t Wlongitude = -1. * mobserv.GetLongitudeDeg()/15.; // h
335 // Warning: GetLongitude return is positive if location is E of Greenwich!!
336 Double_t lst = gmst - Wlongitude; // Local sidereal time [h]
337
338 Double_t CenterHA = lst - CenterRA;
339
340 if (CenterHA < 0.)
341 CenterHA += 24.;
342 else if (CenterHA > 24.)
343 CenterHA -= 24.;
344 // Hour angle of center of camera [h]
345
346
347 // Find Sky coordinates of reconstructed direction:
348
349 starcamtrans->Cel0CamToCel(CenterDec, CenterHA, sourcex/mm2deg,
350 sourcey/mm2deg, SourceDec, SourceHA);
351 SourceRA = lst - SourceHA; // Estimated source RA [h]
352 if (SourceRA < 0.)
353 SourceRA += 24.;
354 else if (SourceRA > 24.)
355 SourceRA -= 24.;
356
357
358 // Generate background
359
360 Float_t r = rndmrad->GetRandom();
361 Float_t phi = 2.*TMath::Pi()*dice.Rndm();
362 Float_t xbg = r * cos(phi);
363 Float_t ybg = r * sin(phi);
364 Double_t BGDec = 0.;
365 Double_t BGHA = 0.;
366 Double_t BGRA = 0.;
367 // Sky coordinates:
368 starcamtrans->Cel0CamToCel(CenterDec, CenterHA, xbg/mm2deg, ybg/mm2deg, BGDec, BGHA);
369 BGRA = lst - BGHA; // Estimated source RA [h]
370 if (BGRA < 0.)
371 BGRA += 24.;
372 else if (BGRA > 24.)
373 BGRA -= 24.;
374
375
376 // Now calculate Alpha with respect to true Crab position
377
378 Double_t CrabHA = lst - CrabRA;
379 if (CrabHA > 24.)
380 CrabHA -= 24.;
381 else if (CrabHA < 0.)
382 CrabHA += 24.;
383
384 Double_t xcrab, ycrab; // mm
385
386 starcamtrans->Cel0CelToCam(CenterDec, CenterHA, CrabDec, CrabHA,
387 xcrab, ycrab);
388
389 xcrab *= mm2deg;
390 ycrab *= mm2deg; // Convert to deg
391
392 // Calculate Alpha relative to Crab:
393
394 Float_t scalar = ((xmean-xcrab)*cosdelta + (ymean-ycrab)*sindelta) /
395 sqrt( pow(xmean-xcrab, 2.) + pow(ymean-ycrab, 2.) );
396
397 Float_t alpha = 180./TMath::Pi()*acos(scalar);
398
399 if (alpha > 90.)
400 alpha -= 180.;
401
402
403 // Hadronness cut
404
405 Float_t hadronness = mhad->GetHadronness();
406 if (hadronness > maxhadronness)
407 continue;
408
409 // WARNING!: If scaling is changed, do it also in the MC part !!
410
411// Float_t swidth = width/(-0.028235+(0.03231*log10size));
412// Float_t slength = length/(-0.13527+(0.09876*log10size));
413
414// Spot 1 cm:
415// Float_t swidth = width/(-0.0332763+(0.0320227*log10size));
416// Float_t slength = length/(-0.174697+(0.107548*log10size));
417
418 Float_t swidth, slength;
419
420 if (spotsize == "1.4cm")
421 {
422 swidth = width/(-0.0308984+(0.0327119*log10size));
423 slength = length/(-0.181605+(0.109609*log10size));
424 }
425 else if (spotsize == "2.0cm")
426 {
427 Float_t swidth = width/(-0.0279008+(0.0344538*log10size));
428 Float_t slength = length/(-0.186394+(0.111055*log10size));
429 }
430
431 Float_t conc = mnew->GetConc();
432 Float_t dist = msrc->GetDist()*mm2deg;
433
434
435 Float_t crab_from_center = sqrt(xcrab*xcrab+ycrab*ycrab);
436 Float_t source_from_center = sqrt(sourcex*sourcex+sourcey*sourcey);
437 // Skip events reconstructed further than 2 deg from center:
438
439 if (source_from_center > 2.)
440 continue;
441
442 // Fill histograms
443
444 CamMap->Fill(sourcex, sourcey);
445 SkyMap->Fill(SourceRA, SourceDec);
446
447 BgMap->Fill(BGRA, BGDec);
448
449
450 // Fill histogram of distance from reconstructed source position to camera
451 // center, but exclude the quadrant where the candidate source is.... (dirty trick)
452
453 if ( (xcrab*sourcex+ycrab*sourcey)/
454 source_from_center / crab_from_center < 0.707) // cos 45 deg
455 r_source->Fill(source_from_center);
456
457
458 // ON region
459
460 Float_t ondist2 = (sourcex-xcrab)*(sourcex-xcrab) +
461 (sourcey-ycrab)*(sourcey-ycrab);
462
463 Float_t offdist2_a = (sourcex+xcrab)*(sourcex+xcrab) +
464 (sourcey+ycrab)*(sourcey+ycrab);
465 Float_t offdist2_b = (sourcex-ycrab)*(sourcex-ycrab) +
466 (sourcey+xcrab)*(sourcey+xcrab);
467 Float_t offdist2_c = (sourcex+ycrab)*(sourcex+ycrab) +
468 (sourcey-xcrab)*(sourcey-xcrab);
469
470
471 if (ondist2 < onr2)
472 {
473 WidthOn->Fill(width);
474 LengthOn->Fill(length);
475 SWidthOn->Fill(swidth);
476 SLengthOn->Fill(slength);
477 ConcOn->Fill(conc);
478 DistOn->Fill(dist);
479
480
481 // 3rd moment with sign referred to source position:
482
483 // Float_t m3long = asylon * TMath::Sign(1., msrc->GetCosDeltaAlpha()) * mm2deg;
484 Float_t m3long = asylon * TMath::Sign(1., xmean-xcrab) * mm2deg;
485
486 M3LongOn->Fill(m3long);
487 SizeOn->Fill(log10size);
488
489 HadronnessOn->Fill(hadronness);
490 ONandOFF->Fill(SourceRA, SourceDec);
491 }
492
493 if (offdist2_a < onr2 || offdist2_b < onr2 || offdist2_c < onr2)
494 {
495 WidthOff->Fill(width, 1./3.);
496 LengthOff->Fill(length, 1./3.);
497 SWidthOff->Fill(swidth, 1./3.);
498 SLengthOff->Fill(slength, 1./3.);
499 ConcOff->Fill(conc, 1./3.);
500 SizeOff->Fill(log10size, 1./3.);
501
502 HadronnessOff->Fill(hadronness, 1./3.);
503 ONandOFF->Fill(SourceRA, SourceDec);
504
505 if (offdist2_a < onr2)
506 {
507 DistOff->Fill(sqrt((xmean+xcrab)*(xmean+xcrab)+(ymean+ycrab)*(ymean+ycrab)), 1./3.);
508
509 m3long = asylon * TMath::Sign(1.,xmean+xcrab) * mm2deg;
510 M3LongOff->Fill(m3long, 1./3.);
511 }
512 else if (offdist2_b < onr2)
513 {
514 DistOff->Fill(sqrt((xmean-ycrab)*(xmean-ycrab)+(ymean+xcrab)*(ymean+xcrab)), 1./3.);
515
516 m3long = asylon * TMath::Sign(1.,xmean-ycrab) * mm2deg;
517 M3LongOff->Fill(m3long, 1./3.);
518 }
519
520 else if (offdist2_c < onr2)
521 {
522 DistOff->Fill(sqrt((xmean+ycrab)*(xmean+ycrab)+(ymean-xcrab)*(ymean-xcrab)), 1./3.);
523
524 m3long = asylon * TMath::Sign(1.,xmean+ycrab) * mm2deg;
525 M3LongOff->Fill(m3long, 1./3.);
526 }
527 }
528
529
530 // Very rough cut on 3rd moment related to Crab position:
531 // (only for alpha plot)
532
533 if ( (asylon < 0. && xcrab < xmean) ||
534 (asylon > 0. && xcrab > xmean) )
535 continue;
536
537 Alpha->Fill(alpha);
538 AbsAlpha->Fill(abs(alpha));
539
540 }
541
542 // Convert y-axis units to Hz
543
544 WidthOn->Scale(1./datatime);
545 LengthOn->Scale(1./datatime);
546 SWidthOn->Scale(1./datatime);
547 ConcOn->Scale(1./datatime);
548 DistOn->Scale(1./datatime);
549 M3LongOn->Scale(1./datatime);
550 SizeOn->Scale(1./datatime);
551 SLengthOn->Scale(1./datatime);
552 HadronnessOn->Scale(1./datatime);
553
554 WidthOff->Scale(1./datatime);
555 LengthOff->Scale(1./datatime);
556 SWidthOff->Scale(1./datatime);
557 ConcOff->Scale(1./datatime);
558 DistOff->Scale(1./datatime);
559 M3LongOff->Scale(1./datatime);
560 SizeOff->Scale(1./datatime);
561 SLengthOff->Scale(1./datatime);
562 HadronnessOff->Scale(1./datatime);
563
564 TH1F* WidthDiff = (TH1F*) WidthOn->Clone("WidthDiff");
565 WidthDiff->Add(WidthOff, -1.);
566
567 TH1F* LengthDiff = (TH1F*) LengthOn->Clone("LengthDiff");
568 LengthDiff->Add(LengthOff, -1.);
569
570 TH1F* SWidthDiff = (TH1F*) SWidthOn->Clone("SWidthDiff");
571 SWidthDiff->Add(SWidthOff, -1.);
572
573 TH1F* SLengthDiff = (TH1F*) SLengthOn->Clone("SLengthDiff");
574 SLengthDiff->Add(SLengthOff, -1.);
575
576 TH1F* ConcDiff = (TH1F*) ConcOn->Clone("ConcDiff");
577 ConcDiff->Add(ConcOff, -1.);
578
579 TH1F* DistDiff = (TH1F*) DistOn->Clone("DistDiff");
580 DistDiff->Add(DistOff, -1.);
581
582 TH1F* M3LongDiff = (TH1F*) M3LongOn->Clone("M3LongDiff");
583 M3LongDiff->Add(M3LongOff, -1.);
584
585 TH1F* SizeDiff = (TH1F*) SizeOn->Clone("SizeDiff");
586 SizeDiff->Add(SizeOff, -1.);
587
588 TH1F* HadronnessDiff = (TH1F*) HadronnessOn->Clone("HadronnessDiff");
589 HadronnessDiff->Add(HadronnessOff, -1.);
590
591
592 TH2F* SkyMapRaw = SkyMap->Clone("SkyMapRaw");
593 TH2F* SkyMapSubtracted= SkyMap->Clone("SkyMapSubtracted");
594
595 SkyMapSubtracted->Add(BgMap, -1.*(SkyMap->Integral() -
596 datatime*WidthDiff->Integral()) /
597 BgMap->Integral());
598
599 SkyMap->Scale(1./datatime);
600 BgMap->Scale(1./datatime);
601
602 SkyMap->Add(BgMap, -1.*(SkyMap->Integral()-WidthDiff->Integral())/BgMap->Integral());
603
604
605 // gStyle->SetPalette(8,0); // Greyscale
606 gStyle->SetPalette(1,0); // Pretty palette
607 // gStyle->SetPalette(51,0); // Deep sea palette
608
609 CamMap->SetStats(0);
610 SkyMapRaw->SetStats(0);
611 SkyMapSubtracted->SetStats(0);
612 SkyMap->SetStats(0);
613
614 // CamMap->Draw("zcol");
615 // SkyMap->Draw("zcol");
616
617 // Alpha->Draw();
618 // AbsAlpha->Draw();
619
620// WidthOn->Draw("e");
621// WidthOff->Draw("same");
622
623// LengthOn->Draw("e");
624// LengthOff->Draw("histo,same");
625
626
627 ///////////////////// MONTE CARLO /////////////////////////
628
629 Double_t xcrab = 120.; //mm
630 Double_t ycrab = 0.;
631
632 xcrab *= mm2deg;
633 ycrab *= mm2deg; // Convert to deg
634
635
636 mhil = new MHillas();
637 mhilsrc = new MHillasSrc();
638 mhilext = new MHillasExt();
639 mnew = new MNewImagePar();
640 mpoint = new MPointingPos();
641 mhad = new MHadronness();
642 MMcEvt* mcevt = new MMcEvt();
643
644 MCchain->SetBranchStatus("*", 0);
645 MCchain->SetBranchStatus("MHillas.fSize", 1);
646 MCchain->SetBranchStatus("MHillas.fWidth", 1);
647 MCchain->SetBranchStatus("MHillas.fLength", 1);
648 MCchain->SetBranchStatus("MHillas.fMeanX", 1);
649 MCchain->SetBranchStatus("MHillas.fMeanY", 1);
650 MCchain->SetBranchStatus("MHillas.fCosDelta", 1);
651 MCchain->SetBranchStatus("MHillas.fSinDelta", 1);
652 MCchain->SetBranchStatus("MHillasSrc.fCosDeltaAlpha", 1);
653 MCchain->SetBranchStatus("MHillasSrc.fDist", 1);
654 MCchain->SetBranchStatus("MHillasExt.fM3Long", 1);
655 MCchain->SetBranchStatus("MNewImagePar.fNumCorePixels", 1);
656 MCchain->SetBranchStatus("MNewImagePar.fConc", 1);
657 MCchain->SetBranchStatus("MHadronness.fHadronness",1);
658 MCchain->SetBranchStatus("MMcEvt.fTelescopeTheta",1);
659 MCchain->SetBranchStatus("MMcEvt.fEnergy",1);
660
661 MCchain->SetBranchAddress("MHillas.", &mhil);
662 MCchain->SetBranchAddress("MHillasExt.", &mhilext);
663 MCchain->SetBranchAddress("MNewImagePar.", &mnew);
664 MCchain->SetBranchAddress("MHadronness.", &mhad);
665 MCchain->SetBranchAddress("MMcEvt.", &mcevt);
666 MCchain->SetBranchAddress("MHillasSrc.", &msrc);
667
668
669 // Original # gammas in the MC sample between 1 and 30 TeV (estimate!):
670 Float_t n_gamma_produced_1to30TeV =
671 n_gamma_produced*(pow(1000.,-1.6) - pow(30000.,-1.6)) /
672 (pow(10.,-1.6) - pow(30000.,-1.6));
673
674 // Correction due to the modification of the spectrum: -2.62 instead of -2.6
675 n_gamma_produced_1to30TeV *= (pow(1000.,-1.62) - pow(30000.,-1.62)) /
676 (pow(1000.,-1.6) - pow(30000.,-1.6));
677
678 Float_t evts2rate = 1.74e-11 / n_gamma_produced_1to30TeV; // cm-2 s-1
679 evts2rate *= TMath::Pi()*pow(300.e2, 2.); // * MC production area in cm2 => s-1
680
681 cout << "Processing Monte Carlo gammas...." << endl;
682
683 Int_t MC_gammas_in_off = 0;
684 Int_t MC_gammas_in_on = 0;
685
686 for (Int_t i = 0; i < MCchain->GetEntries(); i++) // event count starts in 0
687 {
688 if ((i+1)%10000 == 0)
689 cout << i+1 << " of " << MCchain->GetEntries() << endl;
690
691 MCchain->GetEvent(i);
692
693 Float_t size = mhil->GetSize();
694
695 if (size < minsize)
696 continue;
697 else if (size > maxsize)
698 continue;
699
700 Float_t log10size = log10(size);
701
702 Float_t ZA = mcevt->GetTelescopeTheta()*180./TMath::Pi();
703
704 // Exclude MC from zbin0. For Crab ZAmin = 6.8 deg. We do not make exactly that cut
705 // to simplify the flux calculation.
706
707 if (ZA < minMCtheta) // Corresponds to the upper theta limit of zbin0
708 continue;
709 else if (ZA > maxMCtheta)
710 continue;
711
712 Float_t width = mhil->GetWidth()* mm2deg;;
713 Float_t length = mhil->GetLength()* mm2deg;;
714
715 // Float_t disp = (1.-(width/length))*(0.9776+0.101062*log10(size));
716
717 // spot 1 cm:
718 // Float_t disp = (1.-(width/length))*(2.68576-0.830937*log10size+0.1227*log10size*log10size);
719
720 Float_t disp;
721
722 if (spotsize == "1.4cm")
723 disp = (1.-(width/length))*(3.64292-1.28713*log10size+0.179722*log10size*log10size);
724 else if (spotsize == "2.0cm")
725 disp = (1.-(width/length))*(4.11949-1.39848*log10size+0.183514*log10size*log10size);
726
727 Float_t xmean = mhil->GetMeanX()* mm2deg;;
728 Float_t ymean = mhil->GetMeanY()* mm2deg;;
729
730 Float_t cosdelta = mhil->GetCosDelta();
731 Float_t sindelta = mhil->GetSinDelta();
732
733 Float_t asylon = mhilext->GetM3Long();
734
735 Double_t sourcex, sourcey;
736
737 // Choose right orientation:
738
739 if (asylon < 0.)
740 {
741 sourcex = xmean+cosdelta*disp;
742 sourcey = ymean+sindelta*disp;
743 }
744 else
745 {
746 sourcex = xmean-cosdelta*disp;
747 sourcey = ymean-sindelta*disp;
748 }
749
750
751 // Calculate Alpha relative to Crab:
752
753 Float_t scalar = ((xmean-xcrab)*cosdelta + (ymean-ycrab)*sindelta) /
754 sqrt( pow(xmean-xcrab, 2.) + pow(ymean-ycrab, 2.) );
755
756 Float_t alpha = 180./TMath::Pi()*acos(scalar);
757
758 if (alpha > 90.)
759 alpha -= 180.;
760
761
762 // Hadronness cut
763
764 Float_t hadronness = mhad->GetHadronness();
765 if (hadronness > maxhadronness)
766 continue;
767
768// Float_t swidth = width/(-0.028235+(0.03231*log10size));
769// Float_t slength = length/(-0.13527+(0.09876*log10size));
770
771// Spot 1 cm:
772// Float_t swidth = width/(-0.0332763+(0.0320227*log10size));
773// Float_t slength = length/(-0.174697+(0.107548*log10size));
774
775 Float_t swidth, slength;
776
777 if (spotsize == "1.4cm")
778 {
779 swidth = width/(-0.0308984+(0.0327119*log10size));
780 slength = length/(-0.181605+(0.109609*log10size));
781 }
782 else if (spotsize == "2.0cm")
783 {
784 swidth = width/(-0.0279008+(0.0344538*log10size));
785 slength = length/(-0.186394+(0.111055*log10size));
786 }
787
788 Float_t conc = mnew->GetConc();
789 Float_t dist = msrc->GetDist()*mm2deg;
790
791 // ON region
792
793 Float_t ondist2 = (sourcex-xcrab)*(sourcex-xcrab) +
794 (sourcey-ycrab)*(sourcey-ycrab);
795
796 Float_t offdist2_a = (sourcex+xcrab)*(sourcex+xcrab) +
797 (sourcey+ycrab)*(sourcey+ycrab);
798 Float_t offdist2_b = (sourcex-ycrab)*(sourcex-ycrab) +
799 (sourcey+xcrab)*(sourcey+xcrab);
800 Float_t offdist2_c = (sourcex+ycrab)*(sourcex+ycrab) +
801 (sourcey-xcrab)*(sourcey-xcrab);
802
803
804 Float_t energy = mcevt->GetEnergy(); // GeV
805
806 // Calculate weight to account for true Crab spectrum:
807 // We take the (curved) shape of the spectrum BELOW 1 TeV from
808 // Astropart.Phys. 19 (2003) 339-350 (Fabrizio, Konopelko):
809 // dN/dE = C*(E/1TeV)^[-2.47-0.11*log(E/1TeV)]
810 // normalizing it above 1 TeV with the value from ApJ614
811
812 // Integrated Crab flux 1 - 30 TeV from Aharonian et al. 2004, ApJ 614:
813 // 1.74e-11 cm-2 s-1 (shape: power law with index -2.62)
814
815 Float_t mcweight;
816 if (energy > 1000.)
817 mcweight = pow(energy/1000., 2.6-2.62); // correct spectrum
818 else
819 {
820 mcweight = pow(energy/1000., (2.6-2.47-0.11*log(energy/1000.)));
821 }
822
823 mcweight *= evts2rate; // weight to get results in Hz
824
825 // 3rd moment with sign referred to source position, convert to degree:
826
827 // Float_t m3long = asylon * TMath::Sign(1., msrc->GetCosDeltaAlpha()) * mm2deg;
828
829 Float_t m3long = asylon * TMath::Sign(1., xmean-xcrab) * mm2deg;
830
831 M3LongMC_all->Fill(m3long, mcweight);
832
833 if (ondist2 < onr2)
834 {
835 WidthMC->Fill(width, mcweight);
836 LengthMC->Fill(length, mcweight);
837 SWidthMC->Fill(swidth, mcweight);
838 SLengthMC->Fill(slength, mcweight);
839
840 ConcMC->Fill(conc, mcweight);
841 DistMC->Fill(dist, mcweight);
842 M3LongMC->Fill(m3long, mcweight);
843 SizeMC->Fill(log10size, mcweight);
844
845 HadronnessMC->Fill(hadronness, mcweight);
846 EnergyMC->Fill(energy, mcweight);
847
848 MC_gammas_in_on++;
849 }
850
851 if (offdist2_a < onr2 || offdist2_b < onr2 || offdist2_c < onr2)
852 MC_gammas_in_off++;
853
854 // Very rough cut on 3rd moment related to Crab position:
855 // (only for alpha plot)
856
857 if ( (asylon < 0. && xcrab < xmean) ||
858 (asylon > 0. && xcrab > xmean) )
859 continue;
860
861 AlphaMC->Fill(alpha);
862 AbsAlphaMC->Fill(abs(alpha));
863
864 }
865
866 cout << "MC gammas in ON region: " << MC_gammas_in_on << endl;
867 cout << "MC gammas in OFF region: " << MC_gammas_in_off << " / 3" << endl;
868
869 WidthMC->SetLineColor(4);
870 LengthMC->SetLineColor(4);
871 SWidthMC->SetLineColor(4);
872 SLengthMC->SetLineColor(4);
873 HadronnessMC->SetLineColor(4);
874 ConcMC->SetLineColor(4);
875 DistMC->SetLineColor(4);
876 M3LongMC->SetLineColor(4);
877 M3LongMC_all->SetLineColor(2);
878 SizeMC->SetLineColor(4);
879
880 ///////////////////////// DRAW /////////////////////////////////
881
882
883 WidthDiff->Draw();
884 WidthMC->Draw("chisto,same");
885
886 cout << "Data rate: " << WidthDiff->Integral() << " Hz" << endl;
887 cout << "MC rate: " << WidthMC->Integral() << " Hz" << endl;
888 cout << "MC events: " << WidthMC->GetEntries() << endl;
889
890 hfile->Write();
891
892 return;
893}
894
895///////////////////////////////////////
896
897void Smooth(Char_t* filename="Hists_above2000_had0.15.root")
898{
899 TH2F* SkyMap = new TH2F();
900
901 TFile *f = new TFile(filename);
902 SkyMap->Read("SkyMapSubtracted");
903
904
905 Int_t binfactor = 5;
906 Int_t range = 11; // odd number both!
907
908 Int_t nx = binfactor*SkyMap->GetNbinsX();
909 Int_t ny = binfactor*SkyMap->GetNbinsY();
910
911 TH2F* SkyMap2 = new TH2F("SkyMap2", "Events incoming direction map ",
912 nx,
913 SkyMap->GetXaxis()->GetXmin(),
914 SkyMap->GetXaxis()->GetXmax(),
915 ny,
916 SkyMap->GetYaxis()->GetXmin(),
917 SkyMap->GetYaxis()->GetXmax());
918
919
920 // Parameters from fit to real Crab data
921
922 TF2* angres = new TF2("g",
923 "exp(-(x/0.01065)**2-(y/0.135684)**2)",
924 -SkyMap->GetXaxis()->GetBinWidth(1)*range/2.,
925 SkyMap->GetXaxis()->GetBinWidth(1)*range/2.,
926 -SkyMap->GetYaxis()->GetBinWidth(1)*range/2.,
927 SkyMap->GetYaxis()->GetBinWidth(1)*range/2. );
928
929 Int_t nbinsx = (angres->GetXmax() - angres->GetXmin()) /
930 (SkyMap2->GetXaxis()->GetBinWidth(1));
931
932 Int_t nbinsy = (angres->GetYmax() - angres->GetYmin()) /
933 (SkyMap2->GetYaxis()->GetBinWidth(1));
934
935 angres->SetNpx(nbinsx);
936 angres->SetNpy(nbinsy);
937
938 angres->Draw("lego1");
939
940
941 TH2F* hist = (TH2F*)angres->GetHistogram()->Clone("hist");
942
943 Float_t norm = 1./hist->Integral();
944
945 hist->Scale(norm);
946
947 cout << hist->Integral() << endl;
948
949 Float_t maxx = SkyMap2->GetXaxis()->GetXmax();
950 Float_t minx = SkyMap2->GetXaxis()->GetXmin();
951 Float_t maxy = SkyMap2->GetYaxis()->GetXmax();
952 Float_t miny = SkyMap2->GetYaxis()->GetXmin();
953
954 for (Int_t ix = 1; ix <= SkyMap->GetNbinsX(); ix++)
955 {
956 for (Int_t iy = 1; iy <= SkyMap->GetNbinsY(); iy++)
957 {
958 Float_t events = SkyMap->GetBinContent(ix, iy);
959
960 for (Int_t ix2 = 1; ix2 <= nbinsx; ix2++)
961 for (Int_t iy2 = 1; iy2 <= nbinsy; iy2++)
962 {
963 Float_t events2 = events * hist->GetBinContent(ix2, iy2);
964
965 Float_t newx = SkyMap->GetXaxis()->GetBinCenter(ix)
966 + hist->GetXaxis()->GetBinCenter(ix2);
967 Float_t newy = SkyMap->GetYaxis()->GetBinCenter(iy)
968 + hist->GetYaxis()->GetBinCenter(iy2);
969
970 SkyMap2->Fill(newx, newy, events2);
971 }
972 }
973 if (ix%10 == 0)
974 cout << ix << endl;
975 }
976
977 gStyle->SetPalette(1,0);
978 SkyMap2->SetStats(0);
979 SkyMap2->Draw("zcol");
980
981 return;
982}
Note: See TracBrowser for help on using the repository browser.