source: trunk/MagicSoft/Mars/mtemp/meth/macros/CorrectSrcPos.C@ 4993

Last change on this file since 4993 was 4416, checked in by stark, 21 years ago
*** empty log message ***
File size: 10.5 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 Cerenkov 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!
18! Author(s): Luisa Sabrina Stark Schneebeli, 07/2004 <stark@particle.phys.ethz.ch>
19!
20! Copyright: MAGIC Software Development, 2000-2004
21!
22! This macro is partially derived from findstars.C
23\* ======================================================================== */
24
25
26Bool_t HandleInput()
27{
28 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
29 while (1)
30 {
31 //
32 // While reading the input process gui events asynchronously
33 //
34 timer.TurnOn();
35 TString input = Getline("Type 'q' to exit, <return> to go on: ");
36 timer.TurnOff();
37
38 if (input=="q\n")
39 return kFALSE;
40
41 if (input=="\n")
42 return kTRUE;
43 };
44
45 return kFALSE;
46}
47
48
49void CorrectSrcPos(const TString filename="dc_2004_04_21_22_*_Mrk421.root", const TString directory="/scratch/Period016/cacodata/2004_04_22/", const UInt_t numEvents = 0)
50{
51
52 //
53 // Create a empty Parameter List and an empty Task List
54 // The tasklist is identified in the eventloop by its name
55 //
56 MParList plist;
57
58 MTaskList tlist;
59 plist.AddToList(&tlist);
60
61
62 MGeomCamMagic geomcam;
63 MCameraDC dccam;
64 MStarLocalCam starcam;
65
66 plist.AddToList(&geomcam);
67 plist.AddToList(&dccam);
68 plist.AddToList(&starcam);
69
70 //
71 // Now setup the tasks and tasklist:
72 // ---------------------------------
73 //
74
75 // Reads the trees of the root file and the analysed branches
76 MReadReports read;
77 read.AddTree("Currents");
78 read.AddFile(directory+filename); // after the reading of the trees!!!
79 read.AddToBranchList("MReportCurrents.*");
80
81 MGeomApply geomapl;
82 TString continuoslightfile =
83 "/scratch/Period016/cacodata/2004_04_16/dc_2004_04_16_04_46_18_22368_Off3c279-2CL100.root";
84 MCalibrateDC dccal;
85 dccal.SetFileName(continuoslightfile);
86
87
88 const Int_t numblind = 12;
89 const Short_t x[numblind] = { 8, 27,
90 507, 508, 509, 510, 511,
91 543, 559, 560, 561, 567};
92 const TArrayS blindpixels(numblind,(Short_t*)x);
93 Float_t ringinterest = 100; //[mm]
94 Float_t tailcut = 3.5;
95 UInt_t integratedevents = 5;
96
97 MFindStars findstars;
98 // findstars.SetBlindPixels(blindpixels);
99 findstars.SetRingInterest(ringinterest);
100 findstars.SetDCTailCut(tailcut);
101 findstars.SetNumIntegratedEvents(integratedevents);
102 findstars.SetMinuitPrintOutLevel(0);
103
104 findstars.SetSourceRaH(11); // Mrk421
105 findstars.SetSourceRaM(4);
106 findstars.SetSourceRaS(26.0);
107 findstars.SetSourceDecD(38);
108 findstars.SetSourceDecM(12);
109 findstars.SetSourceDecS(36.0);
110
111
112 // prints
113 MPrint pdc("MCameraDC");
114 MPrint pstar("MStarLocalCam");
115
116 tlist.AddToList(&geomapl);
117 tlist.AddToList(&read);
118 tlist.AddToList(&dccal);
119 // tlist.AddToList(&pdc, "Currents");
120 tlist.AddToList(&findstars, "Currents");
121 // tlist.AddToList(&pstar, "Currents");
122
123 //
124 // Create and setup the eventloop
125 //
126 MEvtLoop evtloop;
127 evtloop.SetParList(&plist);
128
129 const Int_t numpar = 3;
130 Float_t starpos[10*numpar];
131 Float_t starposcat[10*numpar];
132
133 if (numEvents > 0)
134 {
135 if (!evtloop.Eventloop(numEvents))
136 return;
137 }
138 else
139 {
140 if (!evtloop.PreProcess())
141 return;
142
143 MHCamera display(geomcam);
144 display.SetPrettyPalette();
145
146 TCanvas *c = new TCanvas("c", "FindStars", 100, 100,600,600);
147
148 TH2D *h0 = new TH2D("h0","",1,-2,2,1,-2,2);
149 h0->GetXaxis()->SetTitle("(deg)");
150 h0->GetYaxis()->SetTitle("(deg)");
151 h0->Draw();
152
153 display.Draw();
154 gPad->cd(1);
155 //display.Draw("same");
156 starcam.Draw();
157
158 UShort_t year;
159 Byte_t month, day;
160
161 MTime *evttime = (MTime*)plist->FindObject("MTimeCurrents");
162
163 TList *starcatalog = new TList;
164 TList *stardc = new TList;
165 TVector3 *res = new TVector3;
166
167 Int_t maxdist = 100;
168 UInt_t numevents=0;
169 Int_t numEvt=0;
170 const Int_t maxnumEvt = 10000;
171 Int_t col = 2;
172 const Double_t kDeg2mm = 189/0.6;
173
174 Double_t resx[maxnumEvt];
175 Double_t resy[maxnumEvt];
176 Double_t time[maxnumEvt];
177
178
179//
180// Find in the file name the date, run and project name
181
182 Size_t pos = filename.Last('/');
183 TString iRun = TString(filename(pos+24,5));
184 TString iYear = TString(filename(pos+4,4));
185 TString iMonth = TString(filename(pos+9,2));
186 TString iDay = TString(filename(pos+12,2));
187
188 TString iHour = TString(filename(pos+15,2));
189 TString iMin = TString(filename(pos+18,2));
190 TString iSec = TString(filename(pos+21,2));
191
192 Size_t poslast = filename.Last('.');
193 Size_t posfirst = poslast-1;
194 while (filename[posfirst] != '_')
195 posfirst--;
196
197 TString iSource = TString(filename(posfirst+1,poslast-posfirst-1));
198 char str[100];
199 char str2[100];
200 sprintf(str,"Date %s %s %s Run %s Source %s",iYear.Data(),iMonth.Data(),iDay.Data(),iRun.Data(),iSource.Data());
201 sprintf(str2,"CorrectSrcPos_%s%s%s_%s_%s.txt",iYear.Data(),iMonth.Data(),iDay.Data(),iRun.Data(),iSource.Data());
202
203 UInt_t evtyear, evtmonth, evtday, evthour, evtmin, evtsec;
204 UShort_t evtmsec;
205
206 ofstream fposout1;
207 // fposout1.open (str2, ofstream::app);
208 fposout1.open (str2, ofstream::trunc);
209
210 while (tlist.Process())
211 {
212 if (numevents >= integratedevents)
213 {
214 evttime->GetTime(evthour, evtmin, evtsec,evtmsec);
215 evttime->Print();
216 evthour = evttime->Hour();
217 evtmin = evttime->Min();
218 evtsec = evttime->Sec();
219
220 // display.SetCamContent(dccam);
221 display.SetCamContent(findstars.GetDisplay());
222 gPad->Modified();
223 gPad->Update();
224 Int_t numStars = starcam.GetNumStars();
225 starcam.Print("bkmaxposrotsizechi");
226
227 res = findstars.GetCorrSrcPos();
228 cout<<"Corrected source position = "<<res->X()<<" , "<<res->Y()<<endl;
229 resx[numEvt] = res->X();
230 resy[numEvt] = res->Y();
231 time[numEvt] = numEvt;
232 numEvt ++;
233
234 // write the corrected source position and the star positions to file
235 if (!fposout1)
236 cout << "Could not open output file to write the center of the fit." << endl;
237 else{
238 for ( Int_t i = 0; i < numStars; i++)
239 {
240 if (i == 0){
241 fposout1 << str<<" "<<evthour<< " "<<evtmin<<" "<<evtsec<<" "<<evtmsec;
242 fposout1 <<" "<< res->X()<< " " << res->Y();
243 }
244 fposout1 <<" "<< starcam[i].GetMeanX()<< " " << starcam[i].GetMeanY();
245 }
246 }
247 fposout1 << endl;
248 numevents = 0;
249
250 // Remove the comments if you want to go through the file
251 // event-by-event:
252 //if (!HandleInput())
253 // break;
254 }
255 numevents++;
256 }
257 fposout1.close();
258
259 // Print source position
260 gROOT->Reset();
261 gStyle->SetCanvasColor(0);
262 gStyle->SetCanvasBorderMode(0);
263 gStyle->SetPadBorderMode(0);
264 gStyle->SetFrameBorderMode(0);
265 gStyle->SetOptStat(00000000);
266
267 TCanvas *c1 = new TCanvas("c1", "SourcePosition", 200, 200,600,1200);
268 c1->Divide(1,3);
269
270 TMath math;
271
272 Double_t minresx, maxresx;
273 minresx = resx[math.LocMin(numEvt,resx)];
274 maxresx = resx[math.LocMax(numEvt,resx)];
275
276 Double_t minresy, maxresy;
277 minresy = resy[math.LocMin(numEvt,resy)];
278 maxresy = resy[math.LocMax(numEvt,resy)];
279
280 Double_t minres, maxres;
281 minres = math.Min(minresx,minresy);
282 maxres = math.Max(maxresx,maxresy);
283
284 Double_t diff;
285 diff = maxres - minres;
286 diff = 0.1*diff;
287 minres = minres - diff;
288 maxres = maxres + diff;
289
290 Double_t mintime, maxtime;
291 mintime = time[math.LocMin(numEvt,time)];
292 maxtime = time[math.LocMax(numEvt,time)];
293
294 c1->cd(1);
295 TH2D *h1 = new TH2D("h1","",1,mintime-1,maxtime+1,1,minres,maxres);
296 h1->GetXaxis()->SetTitle("(evt.nr)");
297 h1->GetYaxis()->SetTitle("(mm)");
298 h1->Draw();
299
300 TGraph *grtimeevolx = new TGraph(numEvt,time,resx);
301 grtimeevolx->SetMarkerColor(215);
302 grtimeevolx->SetMarkerStyle(20);
303 grtimeevolx->SetMarkerSize (0.4);
304 grtimeevolx->Draw("P");
305
306 TGraph *grtimeevoly = new TGraph(numEvt,time,resy);
307 grtimeevoly->SetMarkerColor(6);
308 grtimeevoly->SetMarkerStyle(24);
309 grtimeevoly->SetMarkerSize (0.4);
310 grtimeevoly->Draw("P");
311
312 legxy = new TLegend(0.8,0.85,0.95,0.95);
313 legxy.SetTextSize(0.03);
314 legxy.AddEntry(grtimeevolx,"src pos x","P");
315 legxy.AddEntry(grtimeevoly,"src pos y","P");
316 legxy.Draw();
317
318 c1->cd(2);
319 TH2D *h2 = new TH2D("h2","",1,mintime-1,maxtime+1,1,minresx,maxresx);
320 h2->GetXaxis()->SetTitle("(evt.nr)");
321 h2->GetYaxis()->SetTitle("(mm)");
322 h2->Draw();
323
324 TGraph *grtimeevolx = new TGraph(numEvt,time,resx);
325 grtimeevolx->SetMarkerColor(215);
326 grtimeevolx->SetMarkerStyle(20);
327 grtimeevolx->SetMarkerSize (0.4);
328 grtimeevolx->Draw("P");
329
330 legxy = new TLegend(0.8,0.85,0.95,0.95);
331 legxy.SetTextSize(0.03);
332 legxy.AddEntry(grtimeevolx,"src pos x","P");
333 legxy.Draw();
334
335 c1->cd(3);
336 TH2D *h3 = new TH2D("h3","",1,mintime-1,maxtime+1,1,minresy,maxresy);
337 h3->GetXaxis()->SetTitle("(evt.nr)");
338 h3->GetYaxis()->SetTitle("(mm)");
339 h3->Draw();
340
341 TGraph *grtimeevoly = new TGraph(numEvt,time,resy);
342 grtimeevoly->SetMarkerColor(6);
343 grtimeevoly->SetMarkerStyle(24);
344 grtimeevoly->SetMarkerSize (0.4);
345 grtimeevoly->Draw("P");
346
347 legxy = new TLegend(0.8,0.85,0.95,0.95);
348 legxy.SetTextSize(0.03);
349 legxy.AddEntry(grtimeevoly,"src pos y","P");
350 legxy.Draw();
351
352 evtloop.PostProcess();
353 }
354
355 tlist.PrintStatistics();
356
357}
358
359
Note: See TracBrowser for help on using the repository browser.