Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

source: trunk/MagicSoft/Mars/mtemp/mberlin/macros/SrcCorrect.C@ 6724

Visit:

Last change on this file since 6724 was 4137, checked in by hengsteb, 21 years ago
* empty log message *
File size: 19.1 KB

Line
1	#include "MParList.h"
2	#include "MTaskList.h"
3	#include "MEvtLoop.h"
4
5	#include "MReadTree.h"
6
7	#include "MHillasSrc.h"
8	#include "MSrcPosCam.h"
9
10	#include "MHillasCalc.h"
11	#include "MHillasSrcCalc.h"
12	#include "MTaskInteractive.h"
13	#include "MContinue.h"
14	#include "MRawRunHeader.h"
15	#include "MWriteRootFile.h"
16
17	#include "TH2F.h"
18	#include "TF2.h"
19	#include "TArrayF.h"
20	#include "TObjArray.h"
21	#include "TStyle.h"
22	#include "TFile.h"
23	#include "TMath.h"
24	#include "TString.h"
25	#include "TCanvas.h"
26	#include "TGraphErrors.h"
27	#include "TMarker.h"
28
29	#include "iostream.h"
30
31	//----------------------------------------------------------------
32	// MAIN INPUT/OUTPUT
33
34	#include "IOMkn421.h"
35	//----------------------------------------------------------------
36
37	Double_t rmean=0.3;
38	Double_t dtlimit=10; // time slice width [min]
39	//Double_t hlimit =0.85;
40	Double_t hlimit =0.9;
41
42	//----------------------------------------------------------------
43	// false source hist ---------------------------------------------
44	//----------------------------------------------------------------
45	const Int_t nx = 100; // no. of points in x-grid
46	const Int_t ny = 100; // no. of points in y-grid
47
48	Double_t alimit = 8.; // signal region in alpha-plot [degree]
49
50	Double_t xmin = -1.; // units in degree !!
51	Double_t xmax = 1.;
52	Double_t ymin = -1.;
53	Double_t ymax = 1.;
54
55	Double_t dx=(xmax-xmin)/Double_t(nx);
56	Double_t dy=(ymax-ymin)/Double_t(ny);
57
58	TH2F FSrcHist("FSrcHist","",nx,xmin-dx,xmax+dx,ny,ymin-dy,ymax+dy);
59	TObjArray FSrcArray(1);
60
61	TH1F APlot("APlot","",36,0,90);
62	TObjArray APlotArray(1);
63
64	TArrayF xpeak(100);
65	TArrayF ypeak(100);
66
67	TArrayF xpeakerr(100);
68	TArrayF ypeakerr(100);
69
70	Int_t runno = -1;
71	Int_t runno_old = -1;
72	Int_t runno_lo = -1;
73	Int_t runno_up = -1;
74
75	Double_t mjdend_old = -1.;
76	Double_t mjdstart_lo = -1.;
77	Double_t mjdstart_up = -1.;
78
79	Double_t mjd_lo = -1.;
80	Double_t mjd_up = -1.;
81
82	Double_t dt = 0.;
83	Int_t islice = 0;
84
85	TF2 *fpeak = NULL;
86	FILE *fsrcpos = NULL;
87
88	//****************************************************************************************
89	// first part: track the source in time-intervalls of duration dtlimit
90	//****************************************************************************************
91
92	MHillas *hil=NULL;
93	MRawRunHeader *run=NULL;
94
95	Double_t gaus2d(Double_t x,Double_t par)
96	{
97	Double_t N=par[0]; //integral
98	Double_t mx=par[1]; //mean_x
99	Double_t my=par[2]; //mean_y
100	Double_t sx=par[3]; //sigma_x
101	Double_t sy=par[4]; //sigma_y
102	Double_t rho=par[5];//correlation factor
103	Double_t b=par[6]; //constant background
104
105	Double_t z=(x[0]-mx)*(x[0]-mx)/sx/sx;
106	z+=(x[1]-my)*(x[1]-my)/sy/sy;
107	z-=2.rho(x[0]-mx)*(x[1]-my)/sx/sy;
108
109	Double_t val=N/(2.TMath::Pi()sxsysqrt(1-rho*rho));
110	val=exp( -z/(2.(1.-rho*rho)) );
111	val+=b;
112	return val;
113	}
114
115	Bool_t CalcPeakXY(Double_t xmean, Double_t ymean,
116	Double_t xmeanerr, Double_t ymeanerr)
117	{
118	Int_t mbin=FSrcHist.GetMaximumBin();
119	Int_t biny = mbin/(nx+2);
120	Int_t binx = mbin%(nx+2);
121
122	if(mbin!=FSrcHist.GetBin(binx,biny))
123	{
124	cout<<"Error extracting binx, biny from global bin!"<<endl;
125	cout<<" mbin="<<mbin<<" binfound="<<FSrcHist.GetBin(binx,biny)<<endl;
126	cout<<" binx="<<binx<<" biny="<<biny<<endl;
127
128	return kFALSE;
129	}
130
131	Double_t ndelta = rmean/FSrcHist.GetBinWidth(mbin);
132
133	Double_t nmax=FSrcHist.GetBinContent(binx,biny);
134
135	Int_t ilo=Int_t(binx-ndelta);
136	Int_t iup=Int_t(binx+ndelta);
137
138	Int_t jlo=Int_t(biny-ndelta);
139	Int_t jup=Int_t(biny+ndelta);
140
141	Double_t xsum = 0.; Double_t xsum2 = 0.;
142	Double_t ysum = 0.; Double_t ysum2 = 0.;
143	Double_t cnt = 0.;
144	//Double_t cov = 0.;
145
146	for(Int_t i=ilo;i<=iup;i++)
147	for(Int_t j=jlo;j<=jup;j++)
148	{
149	Double_t n=(Double_t)FSrcHist.GetBinContent(i,j);
150	if(n<hlimit*nmax)continue;
151	Double_t x=(Double_t)FSrcHist.GetXaxis()->GetBinCenter(i);
152	Double_t y=(Double_t)FSrcHist.GetYaxis()->GetBinCenter(j);
153
154	cnt+=n;
155
156	xsum +=n*x;
157	xsum2+=nxx;
158	ysum +=n*y;
159	ysum2+=nyy;
160	}
161
162	Double_t xbar=0.; Double_t xbarerr=0.;
163	Double_t ybar=0.; Double_t ybarerr=0.;
164
165	if(cnt>1)
166	{
167	xbar=xsum/cnt;
168	ybar=ysum/cnt;
169	}
170	if(cnt>2)
171	{
172	xbarerr=cnt(xsum2/cnt-xbarxbar)/(cnt-1.);
173	xbarerr=TMath::Sqrt(xbarerr);
174
175	ybarerr=cnt(ysum2/cnt-ybarybar)/(cnt-1.);
176	ybarerr=TMath::Sqrt(ybarerr);
177	}
178
179	/*
180	//fit not yet tested!!
181	Double_t b=FSrcHist.GetEntries()/nx/ny;
182	Double_t nn=bndeltandelta;
183
184	fpeak->SetParameter(0,nn); // integral
185	fpeak->SetParameter(1,*xmean); // mean_x
186	fpeak->SetParameter(2,*ymean); // mean_y
187	fpeak->SetParameter(3,0.15); // sigma_x
188	fpeak->SetParameter(4,0.15); // sigma_y
189	fpeak->SetParameter(5,0.); // correlation factor
190	fpeak->SetParameter(6,b); // constant background
191
192	FSrcHist.Fit("fpeak");
193	xbar=fpeak->GetParameter(1);
194	ybar=fpeak->GetParameter(2);
195	xbarerr=fpeak->GetParameter(3);
196	ybarerr=fpeak->GetParameter(4);
197	*/
198
199	*xmean=xbar;
200	*ymean=ybar;
201	*xmeanerr=xbarerr;
202	*ymeanerr=ybarerr;
203
204	return kTRUE;
205	}
206
207	void AddFSrcHist(Int_t i,Double_t delta)
208	{
209	Int_t n=FSrcArray.GetSize();
210	if(i>=n) FSrcArray.Expand(n+1);
211
212	TH2F htemp=(TH2F)FSrcHist.Clone();
213	TString name=Form("%d",i);
214	TString title=Form("TimeIntervall = %f min",delta);
215
216	htemp->SetName(name.Data());
217	htemp->SetTitle(title.Data());
218	htemp->SetDrawOption("colz");
219
220	FSrcArray[i] = htemp;
221
222	return;
223	}
224
225	void AddAPlot(Int_t i,Double_t delta)
226	{
227	Int_t n=APlotArray.GetSize();
228	if(i>=n) APlotArray.Expand(n+1);
229
230	TH1F htemp=(TH1F)APlot.Clone();
231	TString name=Form("%d",i);
232	TString title=Form("TimeIntervall = %f min",delta);
233
234	htemp->SetName(name.Data());
235	htemp->SetTitle(title.Data());
236	htemp->SetFillColor(16);
237
238	APlotArray[i] = htemp;
239
240	return;
241	}
242
243	//----------------------------------------------------------------
244	// Source Position Grid ------------------------------------------
245	//----------------------------------------------------------------
246
247	Int_t STrackPreProc(MParList *plist)
248	{
249	hil = (MHillas*) plist->FindObject("MHillas");
250	run = (MRawRunHeader*) plist->FindObject("MRawRunHeader");
251
252	return (hil && run) ? kTRUE : kFALSE;
253	}
254
255	Int_t STrackProc()
256	{
257	runno = run->GetRunNumber();
258	if(runno==0) return kTRUE;
259
260	Double_t mjdstart = (run->GetRunStart()).GetMjd();
261	Double_t mjdend = (run->GetRunEnd()).GetMjd();
262	Double_t t = (mjdend-mjdstart)24.60.;
263
264	Double_t dtrun;
265	if(runno!=runno_old)
266	{
267	if(dt<0.001) {runno_lo=runno;mjd_lo=mjdstart;mjdend_old=mjdstart;} // very first run
268
269	dtrun=(mjdend_old-mjd_lo)24.60.;
270	//if(dt>dtlimit)
271	if(dtrun>dtlimit)
272	{
273	runno_up=runno_old;
274	mjd_up=mjdend_old;
275
276	Double_t x; Double_t xerr;
277	Double_t y; Double_t yerr;
278
279	if(!CalcPeakXY(&x,&y,&xerr,&yerr)) return kFALSE;
280	AddFSrcHist(islice,dt);
281
282	xpeak[islice]=x;
283	ypeak[islice]=y;
284	xpeakerr[islice]=xerr;
285	ypeakerr[islice]=yerr;
286
287	printf("\n TimeSlice %d RunNo %d-%d Mjd %f-%f Dt=%f \n",islice,runno_lo,runno_up,mjd_lo,mjd_up,dt);
288	printf(" X=%f+-%f Y=%f+-%f \n\n",x,xerr,y,yerr);
289	fprintf(fsrcpos,"%d %d %d %f %f %f %f %f %f %f\n",islice,runno_lo,runno_up,x,y,xerr,yerr,mjd_lo,mjd_up,dt);
290
291	FSrcHist.Reset();
292	dt=0.;
293
294	runno_lo=runno;
295	mjd_lo=mjdstart;
296
297	islice++;
298	}
299
300	dt+=t;
301	printf("Runno=%d MjdStart=%f MjdEnde=%f dt=%f \n",runno,mjdstart,mjdend,dt);
302	runno_old=runno;
303	mjdend_old=mjdend;
304	}
305
306
307	MHillasSrc hillas_src;
308	MSrcPosCam srcpos;
309
310	Double_t length=hil->GetLength();
311	Double_t width=hil->GetWidth();
312
313	// size cut
314	if(hil->GetSize()<SizeLo) return kTRUE;
315
316	// cuts in scaled length, width
317	if((length>LengthUp) \|\| (length<LengthLo) \|\| (width>WidthUp) \|\| (width<WidthLo))
318	return kTRUE;
319
320	// cut in dist from camera center
321	srcpos.SetXY(0,0);
322	hillas_src.SetSrcPos(&srcpos);
323	hillas_src.Calc(hil);
324
325	Double_t dist=hillas_src.GetDist()*mm2deg;
326	if(dist>DistUp \|\| dist<DistLo) return kTRUE;
327
328	// false source grid
329	for (Int_t i=0;i<nx;i++)
330	for (Int_t j=0;j<ny;j++)
331	{
332	Double_t x=(xmax-xmin)*Double_t(i)/Double_t(nx-1)+xmin;
333	Double_t y=(ymax-ymin)*Double_t(j)/Double_t(ny-1)+ymin;
334
335	srcpos.SetXY(xdeg2mm,ydeg2mm);
336	hillas_src.SetSrcPos(&srcpos);
337	hillas_src.Calc(hil);
338
339	if(TMath::Abs(hillas_src.GetAlpha())<alimit)
340	FSrcHist.Fill(x,y);
341
342	}
343
344	return kTRUE;
345	}
346
347	Int_t STrackPostProc()
348	{
349	if(dt>0.001)
350	{
351	runno_up=runno_old;
352	mjd_up=mjdend_old;
353
354	Double_t x; Double_t xerr;
355	Double_t y; Double_t yerr;
356
357	if(!CalcPeakXY(&x,&y,&xerr,&yerr)) return kFALSE;
358	AddFSrcHist(islice,dt);
359
360	xpeak[islice]=x;
361	ypeak[islice]=y;
362	xpeakerr[islice]=xerr;
363	ypeakerr[islice]=yerr;
364
365	printf("\n TimeSlice %d RunNo %d-%d Mjd %f-%f Dt=%f \n",islice,runno_lo,runno_up,mjd_lo,mjd_up,dt);
366	printf(" X=%f+-%f Y=%f+-%f \n\n",x,xerr,y,yerr);
367	fprintf(fsrcpos,"%d %d %d %f %f %f %f %f %f %f\n",islice,runno_lo,runno_up,x,y,xerr,yerr,mjd_lo,mjd_up,dt);
368
369	islice++;
370	}
371
372	return kTRUE;
373	}
374
375	//----------------------------------------------------------------
376	//----------------------------------------------------------------
377
378	void SourceTrack(Int_t nent)
379	{
380
381	MParList plist;
382
383	MTaskList tlist;
384	plist.AddToList(&tlist);
385
386	MReadTree read("Events",fileHScaled);
387	read.DisableAutoScheme();
388
389	//TString runcuts("(MRawRunHeader.fRunNumber<17427)");
390	//runcuts+="\|\| (MRawRunHeader.fRunNumber>17428)";
391	//MContinue runsel(runcuts.Data());
392
393	MTaskInteractive strack;
394	strack.SetPreProcess(STrackPreProc);
395	strack.SetProcess(STrackProc);
396	strack.SetPostProcess(STrackPostProc);
397
398	tlist.AddToList(&read);
399	//tlist.AddToList(&runsel);
400	tlist.AddToList(&strack);
401
402	MEvtLoop evtloop;
403	evtloop.SetParList(&plist);
404
405	//------------------------------------------
406	// initializations
407
408	FSrcHist.SetDrawOption("colz");
409	FSrcHist.SetStats(kFALSE);
410
411	fpeak=new TF2("fpeak",gaus2d,-1,1,-1,1,6);
412	fsrcpos=fopen(fileSrcPos.Data(),"w");
413
414	//------------------------------------------
415	// standard eventloop
416
417	if (!evtloop.Eventloop(nent))
418	return;
419
420	tlist.PrintStatistics();
421
422	fclose(fsrcpos);
423
424	TFile file(fileFSrcPlot.Data(),"recreate");
425	FSrcArray.Write();
426	FSrcHist.Write();
427	file.Close();
428
429	const Int_t dim=islice;
430	TGraphErrors g(dim,xpeak.GetArray(),ypeak.GetArray(),
431	xpeakerr.GetArray(),ypeakerr.GetArray());
432
433	//------------------------------------------
434	// canvas with complete hist
435	TCanvas *c1=new TCanvas("c1","",400,800);
436	gStyle->SetPalette(1);
437
438	c1->Divide(2,1);
439	c1->cd(1);
440	gPad->Update();
441
442	FSrcHist.Reset();
443	for(Int_t i=0;i<dim;i++)
444	{
445	TH2F hptr = (TH2F) (FSrcArray[i]);
446	FSrcHist.Add(hptr);
447	}
448
449	TString label=Form("Size>%f ",SizeLo);
450	label+=Form("%f<Dist<%f ", DistLo,DistUp);
451	label+=Form("%f<Width<%f " ,WidthLo,WidthUp);
452	label+=Form("%f<Length<%f ",LengthLo,LengthUp);
453	label+=Form("\|Alpha\|<%f" ,alimit);
454
455	FSrcHist.SetTitle(label.Data());
456	FSrcHist.GetXaxis()->SetTitle("X");
457	FSrcHist.GetYaxis()->SetTitle("Y");
458
459	FSrcHist.DrawClone("colz");
460	g.DrawClone("pl");
461
462	c1->cd(2);
463	gPad->Update();
464	g.GetXaxis()->SetRangeUser(-1.,1.);
465	g.GetYaxis()->SetRangeUser(-1.,1.);
466	g.DrawClone("apl");
467
468	//------------------------------------------
469	// canvas with slides
470	TCanvas *c2=new TCanvas("c2","",800,800);
471	gStyle->SetPalette(1);
472
473	Int_t n=Int_t(TMath::Ceil(sqrt(double(dim))));
474
475	c2->Divide(n,n);
476	for(Int_t i=0;i<dim;i++)
477	{
478	Int_t col=i/n;
479	Int_t row=i%n;
480	c2->cd(row*n+col+1);
481
482	gPad->Update();
483	TH2F hptr = (TH2F) (FSrcArray[i]);
484
485	TMarker m;
486	m.SetMarkerStyle(2);
487	m.SetMarkerSize(4);
488	m.SetMarkerColor(1);
489	m.SetX(xpeak[i]);
490	m.SetY(ypeak[i]);
491
492	hptr->DrawClone("colz");
493	m.DrawClone();
494	}
495
496	return;
497	}
498
499	//****************************************************************************************
500	// second part: correct for moving source position
501	//****************************************************************************************
502
503	// initilizations
504
505	Double_t dtref=0.;
506	Int_t islice_old;
507	Double_t dtref_old=0;
508
509	Double_t xref = 0.; Double_t xreferr = 0.;
510	Double_t yref = 0.; Double_t yreferr = 0.;
511
512	MSrcPosCam *src=NULL;
513
514	//----------------------------------------------------------------
515	// Source Position Grid ------------------------------------------
516	//----------------------------------------------------------------
517	Int_t SCorrectPreProc(MParList *plist)
518	{
519	hil = (MHillas*) plist->FindObject("MHillas");
520	src = (MSrcPosCam*) plist->FindCreateObj("MSrcPosCam");
521	run = (MRawRunHeader*) plist->FindObject("MRawRunHeader");
522
523	return (hil && run && src) ? kTRUE : kFALSE;
524	}
525
526	Int_t SCorrectProc()
527	{
528	if(hil->GetSize()<=0.) return kTRUE;
529
530	Int_t runno = run->GetRunNumber();
531	if(runno==0) return kTRUE;
532
533	Double_t mjdstart = (run->GetRunStart()).GetMjd();
534	Double_t mjdend = (run->GetRunEnd()).GetMjd();
535	Double_t t=(mjdend-mjdstart)24.60.;
536
537	if(runno!=runno_old)
538	{
539	if((runno<runno_lo)\|\|(runno>runno_up))
540	{
541	dt=t;
542
543	fscanf(fsrcpos,"%d %d %d %lf %lf %lf %lf %lf %lf %lf\n",&islice,&runno_lo,&runno_up,&xref,&yref,&xreferr,&yreferr,&mjd_lo,&mjd_up,&dtref);
544	printf("\nApplying correction for runs %d - %d (MJD %f - %f), DtRef=%f :\n",runno_lo,runno_up,mjd_lo,mjd_up,dtref);
545	printf("Source pos set to X=%f Y=%f \n\n",xref,yref);
546
547	if(islice>0)
548	{
549	AddAPlot(islice_old,dtref_old);
550	APlot.Reset();
551	}
552	islice_old=islice;
553	dtref_old=dtref;
554	}
555	else
556	{
557	dt+=t;
558	}
559	printf("Runno=%d MjdStart=%f MjdEnde=%f Dt=%f \n",runno,mjdstart,mjdend,dt);
560	runno_old=runno;
561
562	}
563
564	src->SetXY(xrefdeg2mm,yrefdeg2mm);
565	src->SetReadyToSave();
566	// false source grid
567
568	MHillasSrc hillas_src;
569	MSrcPosCam srcpos;
570
571	Double_t length=hil->GetLength();
572	Double_t width =hil->GetWidth();
573
574	// size cut
575	if(hil->GetSize()<SizeLo) return kTRUE;
576
577	// cuts in scaled length, width
578	if((length>LengthUp) \|\| (length<LengthLo) \|\| (width>WidthUp) \|\| (width<WidthLo))
579	return kTRUE;
580
581	// cut in dist from camera center
582	srcpos.SetXY(0,0);
583	hillas_src.SetSrcPos(&srcpos);
584	hillas_src.Calc(hil);
585
586	Double_t dist=hillas_src.GetDist()*mm2deg;
587	if(dist>DistUp \|\| dist<DistLo) return kTRUE;
588
589	// fill 1-dim APlot for extracted source position
590	srcpos.SetXY(xrefdeg2mm,yrefdeg2mm);
591	hillas_src.SetSrcPos(&srcpos);
592	hillas_src.Calc(hil);
593	APlot.Fill(TMath::Abs(hillas_src.GetAlpha()));
594
595	// false source plot for corrected source movement
596	for (Int_t i=0;i<nx;i++)
597	for (Int_t j=0;j<ny;j++)
598	{
599	Double_t x=(xmax-xmin)*Double_t(i)/Double_t(nx-1)+xmin;
600	Double_t y=(ymax-ymin)*Double_t(j)/Double_t(ny-1)+ymin;
601
602	//printf("xref=%f yref=%f \n",xref,yref);
603	srcpos.SetXY((x+xref)deg2mm,(y+yref)deg2mm);
604	hillas_src.SetSrcPos(&srcpos);
605	hillas_src.Calc(hil);
606
607	if(TMath::Abs(hillas_src.GetAlpha())<alimit)
608	FSrcHist.Fill(x,y);
609	}
610
611	return kTRUE;
612	}
613
614	Int_t SCorrectPostProc()
615	{
616	AddAPlot(islice_old,dtref_old);
617
618	return kTRUE;
619	}
620
621	//----------------------------------------------------------------
622	//----------------------------------------------------------------
623
624	void SourceCorrect(Int_t nent)
625	{
626	MParList plist;
627	MTaskList tlist;
628	plist.AddToList(&tlist);
629
630	MReadTree read("Events", fileHScaled);
631	read.DisableAutoScheme();
632
633	MTaskInteractive scorrect;
634	scorrect.SetPreProcess(SCorrectPreProc);
635	scorrect.SetProcess(SCorrectProc);
636	scorrect.SetPostProcess(SCorrectPostProc);
637
638	MHillasSrcCalc scalc;
639
640	MWriteRootFile write(fileHScaledC.Data());
641	write.AddContainer("MSrcPosCam", "Events");
642	write.AddContainer("MHillas", "Events");
643	write.AddContainer("MHillasExt", "Events");
644	write.AddContainer("MHillasSrc", "Events");
645	write.AddContainer("MNewImagePar", "Events");
646	write.AddContainer("MRawRunHeader","Events");
647
648	tlist.AddToList(&read);
649	tlist.AddToList(&scorrect);
650	tlist.AddToList(&scalc);
651	tlist.AddToList(&write);
652
653	MEvtLoop evtloop;
654	evtloop.SetParList(&plist);
655
656	//------------------------------------------
657	// standard eventloop
658
659	fsrcpos=fopen(fileSrcPos.Data(),"r");
660	if(!fsrcpos)
661	{
662	cout<<"fsrcpos does not exist!"<<endl;
663	return;
664	}
665
666	if (!evtloop.Eventloop(nent))
667	return;
668
669	tlist.PrintStatistics();
670
671	fclose(fsrcpos);
672
673
674	gStyle->SetPalette(1);
675	TFile file(fileFSrcPlotC.Data(),"recreate");
676	FSrcHist.Write();
677	APlotArray.Write();
678	file.Close();
679
680	//------------------------------------------
681	// canvas with complete hist
682	TCanvas *c=new TCanvas("c","",800,400);
683	c->cd();
684
685	TString label=Form("Size>%f ",SizeLo);
686	label+=Form("%f<Dist<%f ", DistLo,DistUp);
687	label+=Form("%f<Width<%f " ,WidthLo,WidthUp);
688	label+=Form("%f<Length<%f ",LengthLo,LengthUp);
689	label+=Form("\|Alpha\|<%f" ,alimit);
690
691	FSrcHist.SetTitle(label.Data());
692	FSrcHist.GetXaxis()->SetTitle("X");
693	FSrcHist.GetYaxis()->SetTitle("Y");
694
695	FSrcHist.DrawCopy("colz");
696
697	//------------------------------------------
698	// canvas with alpha hists
699	TCanvas *c2=new TCanvas("c2","",800,800);
700	gStyle->SetPalette(1);
701
702	const Int_t dim=islice+1;
703	Int_t n=Int_t(TMath::Ceil(sqrt(double(dim))));
704
705	c2->Divide(n,n);
706	for(Int_t i=0;i<dim;i++)
707	{
708	Int_t col=i/n;
709	Int_t row=i%n;
710	c2->cd(row*n+col+1);
711
712	gPad->Update();
713	TH1F hptr = (TH1F) (APlotArray[i]);
714
715	hptr->DrawClone();
716	}
717
718	return;
719	}
720
721	//****************************************************************************************
722	// main
723	//****************************************************************************************
724
725	void SrcCorrect(Int_t icorrect,Int_t nent=-1)
726	{
727
728	if(!icorrect)
729	{
730	// (re-) initilizations
731	islice=0;
732	runno_old=-1;
733	runno_lo=-1;
734	runno_up=-1;
735
736	dt=0;
737
738	mjd_lo=-1;
739	mjd_up=-1;
740
741	hil=NULL; run=NULL;
742
743	cout<<"Starting source track reconstruction."<<endl<<flush;
744	FSrcHist.Reset();
745	SourceTrack(nent);
746
747	}else{
748
749	// (re-) initilizations
750	islice=0;
751	runno_old=-1;
752	runno_lo=-1;
753	runno_up=-1;
754
755	dt=0;
756
757	mjd_lo=-1;
758	mjd_up=-1;
759
760	hil=NULL; run=NULL;
761
762	cout<<"Starting source correct."<<endl<<flush;
763	FSrcHist.Reset();
764	SourceCorrect(nent);
765	}
766
767	return;
768	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: