source: trunk/MagicSoft/Mars/mtemp/findTelAxisFromStars.C@ 5151

Last change on this file since 5151 was 4868, checked in by wittek, 20 years ago
*** empty log message ***
File size: 11.5 KB
Line 
1void ReadSetup(TString fname, MAstroCamera &cam)
2{
3 MMcConfigRunHeader *config=0;
4 MGeomCam *geom=0;
5
6 TFile file(fname);
7 TTree *tree = (TTree*)file.Get("RunHeaders");
8 tree->SetBranchAddress("MMcConfigRunHeader", &config);
9 if (tree->GetBranch("MGeomCam"))
10 tree->SetBranchAddress("MGeomCam", &geom);
11 tree->GetEntry(0);
12
13 cam.SetMirrors(*config->GetMirrors());
14 cam.SetGeom(*geom);
15}
16
17//--------------------------------------------------------------------------
18//
19// findTelAxisFromStars
20//
21// macro to
22// - determine the positions of stars in the camera from the DC currents
23// - and to find the expected position of the source in the camera
24//
25
26void findTelAxisFromStars(const TString filename="20040422_23213_D_Mrk421_E.root", const TString directory="/.magic/magicserv01/MAGIC/rootdata/2004_04_22/", const UInt_t numEvents = 0)
27{
28 gLog.SetNoColors();
29
30 //-------------------------------------------------------------------
31 // book some histograms which are to be filled in the event loop
32 //
33
34 TH2D *aberr = new TH2D("OpticalAberr", "Opt_Aberr",
35 100, 0.0, 650.0, 100, 0.0, 650.0);
36 aberr->SetXTitle("Ideal distance from center [mm]");
37 aberr->SetYTitle("Expected distance from center [mm]");
38
39 TProfile *reldiff = new TProfile("RelativeDiff", "Rel_Diff",
40 100, 0.0, 650.0, -0.2, 0.2, "S");
41 reldiff->SetXTitle("Ideal distance from center [mm]");
42 reldiff->SetYTitle("(R - R0)/R0");
43
44 //-------------------------------------------------------------------
45
46
47 MParList plist;
48 MTaskList tlist;
49 plist.AddToList(&tlist);
50
51 MGeomCamMagic geomcam;
52 MCameraDC dccam;
53 MStarCam starcam;
54
55 //$$$$$$$$$$$$$$$$ ww
56 MObservatory obs;
57 MStarCam sourcecam;
58 sourcecam.SetName("MSourceCam");
59
60 MHTelAxisFromStars htelaxis;
61
62 plist.AddToList(&obs);
63 plist.AddToList(&htelaxis);
64 //$$$$$$$$$$$$$$$$ ww
65
66
67 plist.AddToList(&geomcam);
68 plist.AddToList(&dccam);
69 plist.AddToList(&starcam);
70 plist.AddToList(&sourcecam);
71
72 // Reads the trees of the root file and the analysed branches
73 MReadReports read;
74 read.AddTree("Currents");
75 read.AddTree("Drive"); // If you do not include Drive info, the MFindStars class
76 // will not use the star catalog method
77 read.AddFile(directory+filename); // after the reading of the trees!!!
78 read.AddToBranchList("MReportCurrents.*");
79 read.AddToBranchList("MReportDrive.*");
80
81 MGeomApply geomapl;
82 // TString continuoslightfile =
83 // "/data/MAGIC/Period016/rootdata/2004_04_16/20040416_22368_P_Off3c279-2CL100_E.root";
84
85 TString continuoslightfile =
86"/.magic/magicserv01/MAGIC/rootdata/2004_04_16/20040416_22368_P_Off3c279-2CL100_E.root";
87
88
89 Float_t mindc = 0.9; //[uA]
90
91 MCalibrateDC dccal;
92 dccal.SetFileName(continuoslightfile);
93 dccal.SetMinDCAllowed(mindc);
94
95 const Int_t numblind = 5;
96 const Short_t x[numblind] = { 47, 124, 470, 475, 571};
97 const TArrayS blindpixels(numblind,(Short_t*)x);
98 Float_t ringinterest = 100; //[mm]
99 Float_t tailcut = 2.5;
100 UInt_t integratedevents = 1;
101
102 // We need the MAGIC mirror geometry from a MC header:
103 //TString geometryfile = "/mcdata/standard/camera/NSB_013/Gamma/Gamma_zbin9_90_7_1480to1489_w0.root";
104 TString geometryfile =
105"/.magic/magicserv01/MAGIC/mcdata/standard/camera/NSB_013/Gamma/Gamma_zbin9_90_7_1480to1489_w0.root";
106
107
108
109
110 // We need the Bright Star Catalog:
111 //TString catalogfile = "/home/rwagner/bsc5.dat";
112 TString catalogfile = "mtemp/mmpi/macros/bsc5.dat";
113
114 //$$$$$$$$$$$$$$$$ ww
115 MSourceDirections sdirs;
116 sdirs.SetGeometryFile(geometryfile);
117
118 const Double_t ra = MAstro::Hms2Rad(11, 4, 26);
119 const Double_t dec = MAstro::Dms2Rad(38, 12, 36);
120 sdirs.SetRaDec(ra,dec);
121 sdirs.AddDirection(ra,dec,1,"Mkn 421");
122
123 const Double_t ra = MAstro::Hms2Rad(11, 0, 50);
124 const Double_t dec = MAstro::Dms2Rad(39, 12, 44);
125 sdirs.AddDirection(ra,dec,1,"My_UMa 49");
126
127 const Double_t ra = MAstro::Hms2Rad(11, 4, 31);
128 const Double_t dec = MAstro::Dms2Rad(38, 14, 29);
129 sdirs.AddDirection(ra,dec,1,"My_UMa 51");
130
131 // dummy source dirtections
132 const Double_t ra = MAstro::Hms2Rad(11, 4, 31);
133 const Double_t dec = MAstro::Dms2Rad(38, 20, 0);
134 sdirs.AddDirection(ra,dec,1,"Dummy1");
135
136 const Double_t ra = MAstro::Hms2Rad(11, 4, 31);
137 const Double_t dec = MAstro::Dms2Rad(38, 20, 0);
138 sdirs.AddDirection(ra,dec,1,"Dummy2");
139
140 const Double_t ra = MAstro::Hms2Rad(11, 4, 31);
141 const Double_t dec = MAstro::Dms2Rad(38, 30, 0);
142 sdirs.AddDirection(ra,dec,1,"Dummy3");
143
144 const Double_t ra = MAstro::Hms2Rad(11, 4, 31);
145 const Double_t dec = MAstro::Dms2Rad(38, 40, 0);
146 sdirs.AddDirection(ra,dec,1,"Dummy4");
147
148 const Double_t ra = MAstro::Hms2Rad(11, 4, 31);
149 const Double_t dec = MAstro::Dms2Rad(38, 50, 0);
150 sdirs.AddDirection(ra,dec,1,"Dummy5");
151
152 const Double_t ra = MAstro::Hms2Rad(11, 4, 31);
153 const Double_t dec = MAstro::Dms2Rad(39, 0, 0);
154 sdirs.AddDirection(ra,dec,1,"Dummy6");
155
156 const Double_t ra = MAstro::Hms2Rad(11, 4, 31);
157 const Double_t dec = MAstro::Dms2Rad(39, 10, 0);
158 sdirs.AddDirection(ra,dec,1,"Dummy7");
159
160 const Double_t ra = MAstro::Hms2Rad(11, 4, 31);
161 const Double_t dec = MAstro::Dms2Rad(39, 20, 0);
162 sdirs.AddDirection(ra,dec,1,"Dummy8");
163
164 const Double_t ra = MAstro::Hms2Rad(11, 4, 31);
165 const Double_t dec = MAstro::Dms2Rad(39, 30, 0);
166 sdirs.AddDirection(ra,dec,1,"Dummy9");
167
168 const Double_t ra = MAstro::Hms2Rad(11, 4, 31);
169 const Double_t dec = MAstro::Dms2Rad(39, 40, 0);
170 sdirs.AddDirection(ra,dec,1,"Dummy10");
171
172 const Double_t ra = MAstro::Hms2Rad(11, 4, 31);
173 const Double_t dec = MAstro::Dms2Rad(39, 50, 0);
174 sdirs.AddDirection(ra,dec,1,"Dummy11");
175
176 const Double_t ra = MAstro::Hms2Rad(11, 4, 31);
177 const Double_t dec = MAstro::Dms2Rad(40, 0, 0);
178 sdirs.AddDirection(ra,dec,1,"Dummy12");
179
180 const Double_t ra = MAstro::Hms2Rad(11, 4, 31);
181 const Double_t dec = MAstro::Dms2Rad(40, 10, 0);
182 sdirs.AddDirection(ra,dec,1,"Dummy13");
183
184 const Double_t ra = MAstro::Hms2Rad(11, 4, 31);
185 const Double_t dec = MAstro::Dms2Rad(40, 20, 0);
186 sdirs.AddDirection(ra,dec,1,"Dummy14");
187
188
189 sdirs.SetRadiusFOV(3);
190 //$$$$$$$$$$$$$$$$ ww
191
192 MFindStars findstars;
193 findstars.SetBlindPixels(blindpixels);
194 findstars.SetRingInterest(ringinterest);
195 findstars.SetDCTailCut(tailcut);
196 findstars.SetNumIntegratedEvents(integratedevents);
197 findstars.SetMinuitPrintOutLevel(-1);
198 findstars.SetGeometryFile(geometryfile);
199 findstars.SetBSCFile(catalogfile);
200 const Double_t ra = MAstro::Hms2Rad(11, 4, 26);
201 const Double_t dec = MAstro::Dms2Rad(38, 12, 36);
202 findstars.SetRaDec(ra,dec);
203 findstars.SetLimMag(8);
204 findstars.SetRadiusFOV(1.5);
205
206
207
208 //$$$$$$$$$$$$$$$$ ww
209 Int_t InputType = 0;
210 //findstars.SetUseCorrelatedGauss(kFALSE);
211
212 MTelAxisFromStars telaxis;
213 //telaxis.FixRotationAngleAt(-1.0);
214 //telaxis.FixScaleFactorAt(-1.0);
215 telaxis.SetInputType(InputType);
216
217 MFillH fillhisto("MHTelAxisFromStars[MHTelAxisFromStars]","");
218 htelaxis.SetInputType(InputType);
219 //$$$$$$$$$$$$$$$$ ww
220
221
222 tlist.AddToList(&geomapl);
223 tlist.AddToList(&read);
224 tlist.AddToList(&dccal);
225 tlist.AddToList(&findstars, "Currents");
226
227 //$$$$$$$$$$$$$$$$ ww
228 tlist.AddToList(&sdirs);
229 tlist.AddToList(&telaxis);
230 tlist.AddToList(&fillhisto);
231 //$$$$$$$$$$$$$$$$ ww
232
233
234
235 // The following lines you only need if in addition you want to display
236 // independent MAstroCamera output
237 //
238 // TString fname = "/mcdata/standard/camera/NSB_013/Gamma/Gamma_zbin9_90_7_1480to1489_w0.root";
239 // MObservatory magic1;
240 // const Double_t ra = MAstro::Hms2Rad(11, 4, 26); //Mkn421
241 // const Double_t dec = MAstro::Dms2Rad(38, 12, 36);
242 //
243 // MAstroCamera stars;
244 // ReadSetup(fname, stars);
245 // stars.SetLimMag(9);
246 // stars.SetRadiusFOV(3);
247 // stars.SetRaDec(ra, dec);
248 // stars.ReadBSC("/home/rwagner/bsc5.dat");
249 // stars.SetObservatory(magic1);
250
251 MEvtLoop evtloop;
252 evtloop.SetParList(&plist);
253
254 if (!evtloop.PreProcess())
255 return;
256
257 MHCamera display(geomcam);
258 display.SetPrettyPalette();
259 display.Draw();
260 gPad->cd(1);
261 starcam.Draw();
262
263 UInt_t numevents=0;
264
265 while (tlist.Process())
266 {
267 //gLog << "---------------------------------------------" << endl;
268 //gLog << "Macro : MStarPos content = " << endl;
269 //starcam.Print("namepossizchierr");
270 //gLog << "---------------------------------------------" << endl;
271 //gLog << "Macro : MSourcePos content = " << endl;
272 //sourcecam.Print("namepossizchierr");
273 //gLog << "---------------------------------------------" << endl;
274
275 numevents++;
276 if (numevents%integratedevents==0)
277 {
278 display.SetCamContent(findstars.GetDisplay());
279 gPad->Modified();
280 gPad->Update();
281 // This line prints the results:
282 // starcam.Print();
283
284 // This is how to access the TList of stars:
285 // TList* starlist = starcam.GetList();
286
287 // This is how to iterate over stars found:
288 // TIter Next(starlist);
289 // MStarPos* star;
290 // UInt_t starnum = 0;
291 // cout << filename << " ";
292 // cout << "Iterating over list" << endl;
293 // while ( (star=(MStarPos*)Next()) )
294 // {
295 // cout << "star[" << starnum << "] ";
296 // cout << star->GetMeanX() << " "
297 // << star->GetMeanY() << " ";
298 // starnum++;
299 // }
300 // cout << endl;
301 }
302
303 // fill the histograms
304 // loop over stars in "MSourceCam"
305 TList* sourcelist = sourcecam.GetList();
306 TIter Nextsource(sourcelist);
307 MStarPos* source;
308 while ( (source=(MStarPos*)Nextsource()) )
309 {
310 Double_t R = sqrt( source->GetXExp()*source->GetXExp()
311 + source->GetYExp()*source->GetYExp() );
312
313 Double_t R0= sqrt( source->GetXIdeal()*source->GetXIdeal()
314 + source->GetYIdeal()*source->GetYIdeal() );
315
316 //cout << "R0, R = " << R0 << ", " << R << endl;
317
318 aberr->Fill(R0, R, 1.0);
319 reldiff->Fill(R0, (R-R0)/R0), 1.0);
320 }
321
322 //--------------------------------------
323 MTime time;
324 time.Set(2004, 4, 22, 21, 51, 15);
325
326 //superimpose star picture
327 // stars.SetTime(time);
328 // TObject *o = stars.Clone();
329 // o->SetBit(kCanDelete);
330 // o->Draw();
331
332 // wait after each event
333 //if (!HandleInput())
334 // break;
335
336 }
337
338
339
340 evtloop.PostProcess();
341 tlist.PrintStatistics();
342
343 plist.Print();
344
345 //$$$$$$$$$$$$$$$$ ww
346
347 //gLog << "Event loop finished; call DrawClone of MHTelAxisFromStars" << endl;
348
349 //TObject *srccam = plist.FindObject("MSourceCam");
350 //gLog << "srccam = " << srccam << endl;
351
352 TObject *obj = plist.FindObject("MHTelAxisFromStars");
353 if (obj)
354 {
355 //obj->Print();
356 //obj->Dump();
357 obj->DrawClone();
358 }
359 else
360 gLog << "address of MHTelAxisFromStars container is zero" << endl;
361
362 //$$$$$$$$$$$$$$$$ ww
363
364 // plot the histograms
365 TCanvas *canv = new TCanvas("Optical_Aberration", "Optical aberration",
366 0, 0, 600, 300);
367 canv->Divide(2,1);
368 canv->SetBorderMode(0);
369
370 canv->cd(1);
371 aberr->Draw();
372
373 canv->cd(2);
374 reldiff->Draw();
375
376}
377
378//----------------------------------------------------------------------
379//
380//
381
382Bool_t HandleInput()
383{
384 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
385
386 while (1)
387 {
388 //
389 // While reading the input process gui events asynchronously
390 //
391 timer.TurnOn();
392 TString input = Getline("Type 'q' to exit, <return> to go on: ");
393 timer.TurnOff();
394
395 if (input=="q\n")
396 return kFALSE;
397
398 if (input=="\n")
399 return kTRUE;
400 };
401
402 return kFALSE;
403}
404//=========================================================================
405
406
407
408
409
410
411
412
413
Note: See TracBrowser for help on using the repository browser.