269 fAu=1.49597870700e11;
270 fPc=3.08567758149e16;
280 fHbar=6.58211928e-22;
310 fMapTS.LoadUTCparameterFiles();
311 for (Int_t i=0; i<3; i++)
329 for (Int_t i=0; i<6; i++)
348 for (Int_t i=0; i<5; i++)
355 for (Int_t i=0; i<10; i++)
391 for (Int_t i=0; i<2; i++)
406 if (gROOT->GetListOfCanvases()->FindObject(
"NcAstrolab"))
delete fCanvas;
448 if (fTasks) fTasks->Clear();
467 size=t.fRefs->GetSize();
468 fRefs=new TObjArray(size);
469 for (Int_t i=0; i<size; i++)
471 NcSignal* sx=(NcSignal*)t.fRefs->At(i);
472 if (sx) fRefs->AddAt(sx->Clone(),i);
478 size=t.fSigs->GetSize();
479 fSigs=new TObjArray(size);
480 for (Int_t i=0; i<size; i++)
482 NcSignal* sx=(NcSignal*)t.fSigs->At(i);
483 if (sx) fSigs->AddAt(sx->Clone(),i);
494 for (Int_t ih=0; ih<2; ih++)
499 for (Int_t i=0; i<4; i++)
501 fMarkerSize[i]=t.fMarkerSize[i];
502 fMarkerStyle[i]=t.fMarkerStyle[i];
503 fMarkerColor[i]=t.fMarkerColor[i];
510 SetTimeScramble(t.fTscmode,t.fTscmin,t.fTscmax,t.fTscfunc);
522 SetPositionScramble(t.fRscmode,t.fDscmin,t.fDscmax,t.fDscfunc,t.fThetascmin,t.fThetascmax,t.fThetascfunc,t.fPhiscmin,t.fPhiscmax,t.fPhiscfunc);
531 if (fx) fNuAngle=(TF1*)fx->Clone();
561 TString name=GetName();
562 TString title=GetTitle();
563 printf(
" *%-s::Data*",ClassName());
564 if (name!=
"") printf(
" Name : %-s",name.Data());
565 if (title!=
"") printf(
" Title : %-s",title.Data());
570 printf(
" Position longitude : ");
PrintAngle(l,
"deg",u,3);
571 printf(
" latitude : ");
PrintAngle(b,
"deg",u,3);
572 printf(
" for detector ID : %-i \n",
fLabId);
573 printf(
" Local user reference frame orientation w.r.t X0=South, Y0=East and Z0=Zenith : \n");
574 printf(
" (Note : At the Poles (e.g. IceCube) South means the Greenwich meridian) \n");
575 TString saxes[3]={
"Local X-axis",
"Local Y-axis",
"Local Z-axis"};
576 for (Int_t i=0; i<5; i+=2)
589 if (utc && mode!=4)
Date(4);
593 printf(
" ------------------ Time scrambling ----------------- \n");
596 printf(
" *** Scrambling is applied only for off-source patch (background) data generation in MatchBurstData() *** \n");
597 printf(
" *** The actual stored source and measurement data are not modified *** \n");
602 printf(
" *** Each obtained time difference will be scrambled (mode %-i) *** \n",
fTscmode);
603 printf(
" *** The angular differences are not affected *** \n");
604 printf(
" --> Tailored for scrambling entries in a specific time window without affecting \n");
605 printf(
" the event selection based on angular separation w.r.t. a source. \n");
611 printf(
" *** Each measurement is stored with a scrambled timestamp (mode %-i) *** \n",
fTscmode);
612 printf(
" *** The data of the sources are not modified *** \n");
613 printf(
" *** Off-source (background) data : At each source matching, the measurement is given a new scrambled fake timestamp *** \n");
614 printf(
" --> Tailored for blind analyses without access to the unblinded measurement data \n");
618 printf(
" *** At each source matching, the measurement is given a new scrambled fake timestamp (mode %-i) *** \n",
fTscmode);
620 printf(
" --> Possible patterns in the distribution of sources on the sky will stay intact \n");
624 printf(
" *** At each measurement matching, the source is retrieved from the storage with a scrambled fake timestamp (mode %-i) *** \n",
fTscmode);
625 printf(
" *** The measurements are not affected --> Detection efficiency is unaltered *** \n");
626 printf(
" --> Both time and angular differences will be scrambled \n");
627 printf(
" --> Possible patterns in the distribution of sources on the sky will be washed out \n");
629 TString sx=
"offsets";
630 if (abs(
fTscmode)==1) sx=
"differences";
631 printf(
" Time %-s are randomly drawn from the interval [%-g,%-g] sec. \n",sx.Data(),
fTscmin,
fTscmax);
634 printf(
" Randomising TF1 function %-s is used. \n",
fTscfunc->GetName());
638 printf(
" Uniform randomisation is used. \n");
644 printf(
" ------------------ Position scrambling ------------------ \n");
647 printf(
" *** Scrambling is applied only for off-source patch (background) data generation in MatchBurstData() *** \n");
648 printf(
" *** The actual stored source and measurement data are not modified *** \n");
652 printf(
" *** The data of the sources and the measurement timestamps are not modified *** \n");
658 printf(
" *** Each obtained angular difference will be scrambled (mode %-i) *** \n",
fRscmode);
659 printf(
" *** The time differences are not affected *** \n");
660 printf(
" ---> Tailored for scrambling entries within a specific angular range w.r.t. a source,");
661 printf(
" without affecting the event selection based on the time difference w.r.t. a transient source. \n");
662 printf(
" Angular differences are randomly drawn from the interval [%-g,%-g] degrees. \n",
fDscmin,
fDscmax);
665 printf(
" Randomising TF1 function %-s is used. \n",
fDscfunc->GetName());
669 printf(
" Homogeneous solid angle randomisation is used. \n");
675 printf(
" *** Each measurement is stored with a scrambled fake local position (mode %-i) *** \n",
fRscmode);
676 printf(
" *** Off-source (background) data : At each source matching, the measurement is given a new scrambled fake local position *** \n");
677 printf(
" --> Tailored for blind analyses without access to the unblinded measurement data \n");
682 printf(
" *** At each source matching, the measurement is given a new scrambled fake local position (mode %-i) *** \n",
fRscmode);
687 printf(
" The local coordinates (r,theta,phi) of the measurement are modified to (r+dr,theta+dtheta,phi+dphi) with : \n");
689 printf(
" dr is randomly drawn from the interval [%-g,%-g] \n",
fDscmin,
fDscmax);
692 printf(
" Randomising TF1 function %-s is used. \n",
fDscfunc->GetName());
696 printf(
" Uniform randomisation is used. \n");
702 printf(
" Randomising TF1 function %-s is used. \n",
fThetascfunc->GetName());
706 printf(
" Uniform cos(theta) randomisation is used. \n");
709 printf(
" dphi is randomly drawn from the interval [%-g,%-g] degrees. \n",
fPhiscmin,
fPhiscmax);
712 printf(
" Randomising TF1 function %-s is used. \n",
fPhiscfunc->GetName());
716 printf(
" Uniform phi randomisation is used. \n");
721 printf(
" ------------------ \n");
724 Int_t iseed,cnt1,cnt2;
726 cout <<
" *** Current settings of the internal NcRandom randomiser : iseed=" << iseed <<
" cnt1=" << cnt1 <<
" cnt2=" << cnt2 << endl;
730 cout <<
" *** The internal NcRandom randomiser is currently not intialised ***" << endl;
731 cout <<
" Automatic initialisation will be performed with the actual timestamp at the first random number request." << endl;
732 cout <<
" This will ensure different random sequences for different NcAstrolab instances." << endl;
733 cout <<
" To obtain reproducible scrambled results, please invoke SetRandomiser() before the first SetSignal() invokation." << endl;
735 printf(
" ------------------ \n");
805 Double_t r=1,theta=0,phi=0;
815 Double_t p[3]={r,theta,phi};
816 fLabPos.SetPosition(p,
"sph",
"deg");
871 SetNameTitle(
"User",
"Virtual Lab for general use");
881 if (name==
"Greenwich")
884 SetNameTitle(
"Greenwich",
"The Royal Observatory in the UK");
897 SetNameTitle(
"Amanda",
"Antarctic Muon And Neutrino Detector Array");
908 SetNameTitle(
"IceCube",
"The South Pole Neutrino Observatory");
920 SetNameTitle(
"WSRT",
"The Westerbork Synthesis Radio Telescope");
931 SetNameTitle(
"Astron",
"The Netherlands Institute for Radio Astronomy");
942 SetNameTitle(
"ARA",
"The Askaryan Radio Array at the South Pole");
952 SetNameTitle(
"RNO-G",
"The Greenland Radio Neutrino Observatory at Summit Station");
958 Int_t ids[35]={11,12,13,14,15,16,17,21,22,23,24,25,26,27,33,34,35,36,37,44,45,46,47,54,55,56,57,64,65,66,67,74,75,76,77};
959 Float_t ls[35]={-38.5023,-38.4962,-38.4901,-38.4841,-38.4780,-38.4719,-38.4657,-38.4660,-38.4599,-38.4538,-38.4477,-38.4416,-38.4355,-38.4293,
960 -38.4175,-38.4114,-38.4053,-38.3991,-38.3930,-38.3751,-38.3689,-38.3627,-38.3566,-38.3388,-38.3326,-38.3264,-38.3202,-38.3025,
961 -38.2963,-38.2900,-38.2838,-38.2662,-38.2599,-38.2537,-38.2474};
962 Float_t bs[35]={72.58923,72.60009,72.61095,72.62181,72.63267,72.64353,72.65439,72.58741,72.59827,72.60912,72.61998,72.63084,72.64170,72.65256,
963 72.60729,72.61815,72.62901,72.63987,72.65073,72.61631,72.62717,72.63803,72.64889,72.61447,72.62533,72.63618,72.64704,72.61262,
964 72.62347,72.63433,72.64518,72.61076,72.62161,72.63247,72.64332};
969 for (Int_t i=0; i<35; i++)
991 SetNameTitle(
"ARCA",
"The KM3NeT/ARCA Neutrino detector near Sicily");
999 cout <<
" *" << ClassName() <<
"::SetExperiment* Unsupported experiment name : " << name.Data() << endl;
1000 printf(
" Experiment is set to %-s with detector identifier %-i \n",
fExperiment.Data(),
fLabId);
1050 Double_t pi=acos(-1.);
1052 Double_t offset=90.;
1053 if (u==
"rad") offset=pi/2.;
1056 fLabPos.GetPosition(p,
"sph",u);
1165 if (!
fRan)
return 0;
1167 iseed=
fRan->GetSeed();
1168 cnt1=
fRan->GetCnt1();
1169 cnt2=
fRan->GetCnt2();
1297 if (out==
"deg" || out==
"rad")
1301 printf(
"%*.*f %-s",5+ndig,ndig,b,out.Data());
1305 printf(
"%-.*f %-s",ndig,b,out.Data());
1310 Double_t epsilon=1.e-12;
1311 Int_t word=0,ddd=0,hh=0,mm=0,ss=0;
1322 s=fabs(b)-Double_t(ddd*10000+mm*100+ss);
1348 printf(
" -0d %02i' %0*.*f\"",mm,3+ndig,ndig,s);
1352 printf(
"%4id %02i' %0*.*f\"",ddd,mm,3+ndig,ndig,s);
1359 printf(
"-0d %-i' %-.*f\"",mm,ndig,s);
1363 printf(
"%-id %-i' %-.*f\"",ddd,mm,ndig,s);
1377 s=fabs(b)-Double_t(hh*10000+mm*100+ss);
1403 printf(
" -0h %02im %0*.*fs",mm,3+ndig,ndig,s);
1407 printf(
"%3ih %02im %0*.*fs",hh,mm,3+ndig,ndig,s);
1414 printf(
"-0h %-im %-.*fs",mm,ndig,s);
1418 printf(
"%-ih %-im %-.*fs",hh,mm,ndig,s);
1512 if (frame!=
"equ" && frame!=
"gal" && frame!=
"ecl" && frame!=
"hor" && frame!=
"icr" && frame!=
"loc")
return 0;
1514 if (frame==
"equ" && mode!=
"M" && mode!=
"m" && mode!=
"T" && mode!=
"t" && mode!=
"B" && mode!=
"b" && mode!=
"J" && mode!=
"j")
return 0;
1527 ts->
GetMJD(mjd,isec,ins);
1529 ts2.
SetMJD(mjd,isec,ins,ips);
1532 vec[0]=r->
GetX(1,
"sph",
"rad");
1533 vec[1]=r->
GetX(2,
"sph",
"rad");
1534 vec[2]=r->
GetX(3,
"sph",
"rad");
1545 sx=
SetSignal(&q,
"hor",mode,&ts2,jref,name,type);
1586 if (mode==
"T" || mode==
"t")
1593 if (mode==
"T" || mode==
"t" || mode==
"M" || mode==
"m")
1600 if (mode==
"B" || mode==
"b") te.
SetEpoch(1950,
"B");
1601 if (mode==
"J" || mode==
"j") te.
SetEpoch(2000,
"J");
1665 fSigs=
new TObjArray();
1669 size=
fSigs->GetSize();
1671 jlast=
fSigs->GetLast();
1674 if (jref>size)
fSigs->Expand(jref);
1678 if (jlast==jmax)
fSigs->Expand(size+1);
1696 sxsig->SetName(name);
1697 sxsig->SetTitle(
"Observed event in ICRS coordinates");
1706 fSigs->AddAt(sxsig,jref-1);
1715 fRefs=
new TObjArray();
1719 size=
fRefs->GetSize();
1721 jlast=
fRefs->GetLast();
1724 if (jref>size)
fRefs->Expand(jref);
1728 if (jlast==jmax)
fRefs->Expand(size+1);
1746 sxref->SetName(name);
1747 sxref->SetTitle(
"Reference event in ICRS coordinates");
1756 fRefs->AddAt(sxref,jref-1);
1761 if (
fRscmode !=2 || !type)
return sx;
1768 Int_t index=
fSigs->IndexOf(sx);
1805 Double_t pi=acos(-1.);
1810 Float_t temp=cosmin;
1814 Double_t cosang=
fRan->Uniform(cosmin,cosmax);
1815 dtheta=acos(cosang)*180./pi;
1831 if (vec[0]<=0) vec[0]=1e-20;
1894 if (name.Contains(
"*"))
return 0;
1904 ts->
Almanac(0,0,0,0,name,&lx,&bx,&rx);
1912 if (set && jref)
SetSignal(rx,lx,
"deg",bx,
"deg",
"ecl",ts,jref,
"M",name,type);
1917NcSignal*
NcAstrolab::SetSignal(Double_t d,Double_t a,TString au,Double_t b,TString bu,TString frame,
NcTimestamp* ts,Int_t jref,TString mode,TString name,Int_t type)
2006 Double_t vec[3]={0,0,0};
2012 if (mode!=
"M" && mode!=
"m" && mode!=
"T" && mode!=
"t" && mode!=
"B" && mode!=
"b" && mode!=
"J" && mode!=
"j")
return 0;
2046 if (frame==
"loc" || frame==
"pdir")
2067NcSignal*
NcAstrolab::SetSignal(Double_t d,Double_t a,TString au,Double_t b,TString bu,TString frame,TString s,Double_t e,Int_t jref,TString mode,TString name,Int_t type)
2154 NcSignal* sx=
SetSignal(d,a,au,b,bu,frame,&tx,jref,mode,name,type);
2241 if (!s)
return -999;
2308 TObjArray* arr=
fRefs;
2309 if (type) arr=
fSigs;
2316 n=arr->GetEntries();
2388 if (frame!=
"equ" && frame!=
"gal" && frame!=
"ecl" && frame!=
"hor" && frame!=
"icr" && frame!=
"loc")
return 0;
2390 if (frame==
"equ" && mode!=
"M" && mode!=
"m" && mode!=
"T" && mode!=
"t" && mode!=
"B" && mode!=
"b" && mode!=
"J" && mode!=
"j")
return 0;
2411 if (fabs(dt)>
fMaxDt)
return 0;
2415 TString name=sx->GetName();
2433 if (frame==
"equ" && mode!=
"J" && mode!=
"j")
2438 if (mode!=
"B" && mode!=
"b")
2450 if (mode==
"T" || mode==
"t")
Nutate(q,ts);
2543 Double_t xg=R*cos(B)*cos(L);
2544 Double_t yg=R*cos(B)*sin(L);
2545 Double_t zg=R*sin(B);
2549 ts->
Almanac(0,0,0,0,
"Earth*",&l0,&b0,&r0);
2554 Double_t x0=r0*cos(b0)*cos(l0);
2555 Double_t y0=r0*cos(b0)*sin(l0);
2556 Double_t z0=r0*sin(b0);
2564 Double_t Rh=sqrt(xh*xh+yh*yh+zh*zh);
2565 Double_t Bh=atan2(zh,sqrt(xh*xh+yh*yh));
2566 Double_t Lh=atan2(yh,xh);
2568 Double_t twopi=2.*acos(-1.);
2569 while (Lh<0) { Lh+=twopi; }
2570 while (Lh>twopi) { Lh-=twopi; }
2621 Double_t xh=R*cos(B)*cos(L);
2622 Double_t yh=R*cos(B)*sin(L);
2623 Double_t zh=R*sin(B);
2627 ts->
Almanac(0,0,0,0,
"Earth*",&l0,&b0,&r0);
2632 Double_t x0=r0*cos(b0)*cos(l0);
2633 Double_t y0=r0*cos(b0)*sin(l0);
2634 Double_t z0=r0*sin(b0);
2642 Double_t Rg=sqrt(xg*xg+yg*yg+zg*zg);
2643 Double_t Bg=atan2(zg,sqrt(xg*xg+yg*yg));
2644 Double_t Lg=atan2(yg,xg);
2646 Double_t twopi=2.*acos(-1.);
2647 while (Lg<0) { Lg+=twopi; }
2648 while (Lg>twopi) { Lg-=twopi; }
2747 if (frame==
"equ" || frame==
"gal" || frame==
"ecl" || frame==
"icr")
2868 if (j>=0) sx=
GetSignal(d,a,au,b,bu,frame,ts,j,mode,type);
2872NcSignal*
NcAstrolab::GetSignal(Double_t& d,Double_t& a,TString au,Double_t& b,TString bu,TString frame,TString s,Double_t e,Int_t jref,TString mode,Int_t type)
2953 if (s!=
"B" && s!=
"b" && s!=
"J" && s!=
"j")
return 0;
2962NcSignal*
NcAstrolab::GetSignal(Double_t& d,Double_t& a,TString au,Double_t& b,TString bu,TString frame,TString s,Double_t e,TString name,TString mode,Int_t type)
3047 if (j>=0) sx=
GetSignal(d,a,au,b,bu,frame,s,e,j,mode,type);
3074 if (jref<0)
return 0;
3082 if (!type && !
fRefs)
return 0;
3083 if (type && !
fSigs)
return 0;
3088 if (jref<=fSigs->GetSize()) sx=(
NcSignal*)
fSigs->At(jref-1);
3092 if (jref<=fRefs->GetSize()) sx=(
NcSignal*)
fRefs->At(jref-1);
3157 if (!
fRefs)
return 0;
3162 nrem=
fRefs->GetEntries();
3169 if (j>0 && j<=fRefs->GetSize())
3171 TObject* obj=
fRefs->RemoveAt(j-1);
3180 if (compress)
fRefs->Compress();
3239 TObjArray* arr=
fRefs;
3240 if (type) arr=
fSigs;
3242 if (!arr)
return nrem;
3247 nrem=arr->GetEntries();
3263 if (j>0 && j<=arr->GetSize())
3265 TObject* obj=arr->RemoveAt(j-1);
3277 Int_t n=arr->GetEntries();
3311 if (name==
"")
return nrem;
3353 if (name==
"")
return 0;
3361 TObjArray* arr=
fRefs;
3362 if (type) arr=
fSigs;
3369 for (Int_t j=1; j<=arr->GetSize(); j++)
3374 namex=sx->GetName();
3375 if (!namex.Contains(name))
continue;
3385 Int_t n=arr->GetEntries();
3408 if (name==
"")
return -1;
3412 TObjArray* arr=
fRefs;
3413 if (type) arr=
fSigs;
3415 if (!arr)
return -1;
3417 for (Int_t i=0; i<arr->GetSize(); i++)
3422 if (name==sx->GetName())
3449 TObjArray* arr=
fRefs;
3450 if (type) arr=
fSigs;
3452 if (!arr)
return -1;
3454 for (Int_t i=0; i<arr->GetSize(); i++)
3561 TString slha=
"LAHA";
3562 if (mode==
"M" || mode==
"m")
3567 else if ((mode==
"B" || mode==
"b" || mode==
"J" || mode==
"j") && emode==
"M")
3576 d=90.-r.
GetX(2,
"sph",
"deg");
3577 a=r.
GetX(3,
"sph",
"rad");
3578 if (mode==
"B" || mode==
"b") mode=
"B1950";
3579 if (mode==
"J" || mode==
"j") mode=
"J2000";
3580 cout <<
"Equatorial (" << mode.Data() <<
") | a :";
3581 if (!align) {cout <<
" ";}
PrintAngle(a,
"rad",
"hms",ndig,align);
3582 if (!align) {cout <<
" ";}
PrintAngle(a,
"rad",
"deg",ndig,align);
3584 if (!align) {cout <<
" ";}
PrintAngle(d,
"deg",
"dms",ndig,align);
3585 if (!align) {cout <<
" ";}
PrintAngle(d,
"deg",
"deg",ndig,align);
3586 cout <<
" | " << slha.Data() <<
" : ";
3588 cout <<
" ";
PrintAngle(lha,
"deg",
"deg",ndig,align);
3594 b=90.-r.
GetX(2,
"sph",
"deg");
3595 l=r.
GetX(3,
"sph",
"deg");
3596 cout <<
"Galactic | l :";
3597 if (!align) {cout <<
" ";}
PrintAngle(l,
"deg",
"deg",ndig,align);
3598 if (!align) {cout <<
" ";}
PrintAngle(l,
"deg",
"dms",ndig,align);
3600 if (!align) {cout <<
" ";}
PrintAngle(b,
"deg",
"deg",ndig,align);
3601 if (!align) {cout <<
" ";}
PrintAngle(b,
"deg",
"dms",ndig,align);
3602 cout <<
" | " << slha.Data() <<
" : ";
3604 cout <<
" ";
PrintAngle(lha,
"deg",
"deg",ndig,align);
3610 b=90.-r.
GetX(2,
"sph",
"deg");
3611 l=r.
GetX(3,
"sph",
"deg");
3612 cout <<
"ICRS | l :";
3613 if (!align) {cout <<
" ";}
PrintAngle(l,
"deg",
"deg",ndig,align);
3614 if (!align) {cout <<
" ";}
PrintAngle(l,
"deg",
"dms",ndig,align);
3616 if (!align) {cout <<
" ";}
PrintAngle(b,
"deg",
"deg",ndig,align);
3617 if (!align) {cout <<
" ";}
PrintAngle(b,
"deg",
"dms",ndig,align);
3618 cout <<
" | " << slha.Data() <<
" : ";
3620 cout <<
" ";
PrintAngle(lha,
"deg",
"deg",ndig,align);
3626 b=90.-r.
GetX(2,
"sph",
"deg");
3627 l=r.
GetX(3,
"sph",
"deg");
3628 cout <<
"Geocentric ecliptic | l :";
3629 if (!align) {cout <<
" ";}
PrintAngle(l,
"deg",
"deg",ndig,align);
3630 if (!align) {cout <<
" ";}
PrintAngle(l,
"deg",
"dms",ndig,align);
3632 if (!align) {cout <<
" ";}
PrintAngle(b,
"deg",
"deg",ndig,align);
3633 if (!align) {cout <<
" ";}
PrintAngle(b,
"deg",
"dms",ndig,align);
3634 cout <<
" | " << slha.Data() <<
" : ";
3636 cout <<
" ";
PrintAngle(lha,
"deg",
"deg",ndig,align);
3641 Double_t alt=90.-r.
GetX(2,
"sph",
"deg");
3642 Double_t azi=180.-r.
GetX(3,
"sph",
"deg");
3651 cout <<
"Horizontal | azi :";
3652 if (!align) {cout <<
" ";}
PrintAngle(azi,
"deg",
"deg",ndig,align);
3653 if (!align) {cout <<
" ";}
PrintAngle(azi,
"deg",
"dms",ndig,align);
3655 if (!align) {cout <<
" ";}
PrintAngle(alt,
"deg",
"deg",ndig,align);
3656 if (!align) {cout <<
" ";}
PrintAngle(alt,
"deg",
"dms",ndig,align);
3657 cout <<
" | " << slha.Data() <<
" : ";
3659 cout <<
" ";
PrintAngle(lha,
"deg",
"deg",ndig,align);
3664 Double_t theta=r.
GetX(2,
"sph",
"deg");
3665 Double_t phi=r.
GetX(3,
"sph",
"deg");
3666 cout <<
"Local-frame | phi :";
3667 if (!align) {cout <<
" ";}
PrintAngle(phi,
"deg",
"deg",ndig,align);
3668 if (!align) {cout <<
" ";}
PrintAngle(phi,
"deg",
"dms",ndig,align);
3669 cout <<
" | theta :";
3670 if (!align) {cout <<
" ";}
PrintAngle(theta,
"deg",
"deg",ndig,align);
3671 if (!align) {cout <<
" ";}
PrintAngle(theta,
"deg",
"dms",ndig,align);
3672 cout <<
" | " << slha.Data() <<
" : ";
3674 cout <<
" ";
PrintAngle(lha,
"deg",
"deg",ndig,align);
3679 TString name=sx->GetName();
3680 if (name !=
"") cout <<
" " << name.
Data();
3765 if (j>=0)
PrintSignal(frame,mode,ts,ndig,j,emode,type,align);
3844 if (nmax>0) width=log10(nmax);
3852 if (maxsig>maxj) maxj=maxsig;
3853 if (maxj>0) width=log10(maxj);
3861 if (mode==
"T" || mode==
"t") dform=-1;
3862 if ((mode==
"B" || mode==
"b" || mode==
"J" || mode==
"j") && emode==
"T") dform=-1;
3870 printf(
" *%-s::ListSignals* Name : %-s Title : %-s \n",ClassName(),GetName(),GetTitle());
3873 cout <<
" Timestamp of the measurement stored at index=" << j;
3877 cout <<
" *Scrambled* timestamp of the measurement stored at index=" << j;
3879 cout <<
" (Lab time offset w.r.t. UT : ";
PrintTime(
fToffset,12); cout <<
")" << endl;
3882 cout <<
" Corresponding location of this measurement" << endl;
3883 cout <<
" ";
PrintSignal(frame,mode,tx,ndig,j,emode,1); cout << endl;
3914 nstored=arr->GetEntries();
3917 for (Int_t i=1; i<=arr->GetSize(); i++)
3923 if (nmax>=0 && jlist>nmax)
break;
3926 namex=sx->GetName();
3927 if (name!=
"*" && !namex.Contains(name))
continue;
3931 printf(
" *%-s::ListSignals* Name : %-s Title : %-s \n",ClassName(),GetName(),GetTitle());
3935 cout <<
" User provided timestamp (Lab time offset w.r.t. UT : ";
PrintTime(
fToffset,12); cout <<
")";
3945 cout <<
" Current timestamp of the laboratory (Lab time offset w.r.t. UT : ";
PrintTime(
fToffset,12); cout <<
")";
3955 if (nmax<0 || nmax>=nstored)
3961 cout <<
" === All stored reference signals according to the above timestamp ===" << endl;
3965 cout <<
" === All stored reference signals according to their actual recorded timestamp ===" << endl;
3972 cout <<
" === All stored measurements according to their actual observation timestamp ===" << endl;
3976 cout <<
" === All stored measurements according to their *scrambled* observation timestamp ===" << endl;
3977 cout <<
" === Time scrambling was performed by adding dt from the interval [dtmin,dtmax] to their actual timestamp" << endl;
3978 cout <<
" === dtmin : " <<
fTscmin <<
" dtmax : " <<
fTscmax <<
" sec.";
3981 cout <<
" Randomising TF1 function " <<
fTscfunc->GetName() <<
" was used." << endl;
3985 cout <<
" Uniform randomisation was used." << endl;
3996 cout <<
" === The first " << nmax <<
" stored reference signals according to the above timestamp ===" << endl;
4000 cout <<
" === The first " << nmax <<
" stored reference signals according to their actual recorded timestamp ===" << endl;
4007 cout <<
" === The first " << nmax <<
" stored measurements according to their actual observation timestamp ===" << endl;
4011 cout <<
" === The first " << nmax <<
" stored measurements according to their *scrambled* observation timestamp ===" << endl;
4012 cout <<
" === Time scrambling was performed by adding dt from the interval [dtmin,dtmax] to their actual timestamp" << endl;
4013 cout <<
" === dtmin : " <<
fTscmin <<
" dtmax : " <<
fTscmax <<
" sec.";
4016 cout <<
" Randomising TF1 function " <<
fTscfunc->GetName() <<
" was used." << endl;
4020 cout <<
" Uniform randomisation was used." << endl;
4027 if (type==1 || (!type && j<0)) tx=0;
4028 printf(
" Index : %*d ",width,i);
PrintSignal(frame,mode,tx,ndig,i,emode,type,kTRUE); printf(
"\n");
4111 Double_t pi=acos(-1.);
4115 Double_t x=-16.6170;
4119 a*=pi/(180.*3600.*1000.);
4120 x*=pi/(180.*3600.*1000.);
4121 e*=pi/(180.*3600.*1000.);
4124 mat[0]=1.-0.5*(a*a+x*x);
4128 mat[4]=1.-0.5*(a*a+e*e);
4132 mat[8]=1.-0.5*(e*e+x*x);
4150 Double_t mat[9]={0,0,0,0,0,0,0,0,0};
4157 Double_t pi=acos(-1.);
4159 Double_t t=(ts->
GetJD()-2451545.0)/36525.;
4162 Double_t eps0=84381.406;
4163 Double_t psi=5038.481507*t-1.0790069*pow(t,2)-0.00114045*pow(t,3)+0.000132851*pow(t,4)
4164 -0.0000000951*pow(t,4);
4165 Double_t om=eps0-0.025754*t+0.0512623*pow(t,2)-0.00772503*pow(t,3)-0.000000467*pow(t,4)
4166 +0.0000003337*pow(t,5);
4167 Double_t chi=10.556403*t-2.3814292*pow(t,2)-0.00121197*pow(t,3)+0.000170663*pow(t,4)
4168 -0.0000000560*pow(t,5);
4171 eps0*=pi/(180.*3600.);
4172 psi*=pi/(180.*3600.);
4173 om*=pi/(180.*3600.);
4174 chi*=pi/(180.*3600.);
4176 Double_t s1=sin(eps0);
4177 Double_t s2=sin(-psi);
4178 Double_t s3=sin(-om);
4179 Double_t s4=sin(chi);
4180 Double_t c1=cos(eps0);
4181 Double_t c2=cos(-psi);
4182 Double_t c3=cos(-om);
4183 Double_t c4=cos(chi);
4185 mat[0]=c4*c2-s2*s4*c3;
4186 mat[1]=c4*s2*c1+s4*c3*c2*c1-s1*s4*s3;
4187 mat[2]=c4*s2*s1+s4*c3*c2*s1+c1*s4*s3;
4188 mat[3]=-s4*c2-s2*c4*c3;
4189 mat[4]=-s4*s2*c1+c4*c3*c2*c1-s1*c4*s3;
4190 mat[5]=-s4*s2*s1+c4*c3*c2*s1+c1*c4*s3;
4192 mat[7]=-s3*c2*c1-s1*c3;
4193 mat[8]=-s3*c2*s1+c3*c1;
4209 Double_t mat[9]={0,0,0,0,0,0,0,0,0};
4216 Double_t pi=acos(-1.);
4218 Double_t dpsi,deps,eps;
4219 ts->
Almanac(&dpsi,&deps,&eps);
4222 dpsi*=pi/(180.*3600.);
4223 deps*=pi/(180.*3600.);
4224 eps*=pi/(180.*3600.);
4226 Double_t s1=sin(eps);
4227 Double_t s2=sin(-dpsi);
4228 Double_t s3=sin(-(eps+deps));
4229 Double_t c1=cos(eps);
4230 Double_t c2=cos(-dpsi);
4231 Double_t c3=cos(-(eps+deps));
4237 mat[4]=c3*c2*c1-s1*s3;
4238 mat[5]=c3*c2*s1+c1*s3;
4240 mat[7]=-s3*c2*c1-s1*c3;
4241 mat[8]=-s3*c2*s1+c3*c1;
4266 Double_t vec[3]={1,0,0};
4303 fG.SetAngles(x.
GetX(2,
"sph",
"deg"),x.
GetX(3,
"sph",
"deg"),
4304 y.
GetX(2,
"sph",
"deg"),y.
GetX(3,
"sph",
"deg"),
4305 z.
GetX(2,
"sph",
"deg"),z.
GetX(3,
"sph",
"deg"));
4319 Double_t dpsi,deps,eps;
4320 ts->
Almanac(&dpsi,&deps,&eps);
4329 Double_t theta2=90.-eps;
4331 Double_t theta3=eps;
4334 fE.SetAngles(theta1,phi1,theta2,phi2,theta3,phi3);
4359 Double_t vec[3]={1,0,0};
4378 fH.SetAngles(x.
GetX(2,
"sph",
"deg"),x.
GetX(3,
"sph",
"deg"),
4379 y.
GetX(2,
"sph",
"deg"),y.
GetX(3,
"sph",
"deg"),
4380 z.
GetX(2,
"sph",
"deg"),z.
GetX(3,
"sph",
"deg"));
4417 fL.SetAngles(t1,p1,t2,p2,t3,p3);
4471 if (in==out)
return a;
4474 Double_t pi=acos(-1.);
4475 Double_t epsilon=1.e-12;
4476 Int_t word=0,ddd=0,hh=0,mm=0,ss=0;
4481 if (in==
"rad") b*=180./pi;
4483 if (in==
"hrs") b*=15.;
4492 s=b-Double_t(ddd*10000+mm*100+ss);
4493 b=Double_t(ddd)+Double_t(mm)/60.+(Double_t(ss)+s)/3600.;
4503 s=b-Double_t(hh*10000+mm*100+ss);
4504 b=15.*(Double_t(hh)+Double_t(mm)/60.+(Double_t(ss)+s)/3600.);
4512 if (out==
"rad") b*=pi/180.;
4514 if (out==
"hrs") b/=15.;
4545 b=Double_t(10000*ddd+100*mm+ss)+s;
4578 b=Double_t(10000*hh+100*mm+ss)+s;
4610 omega=(ph2-ph1)*(cos(th1)-cos(th2));
4611 if (omega<0) omega=-omega;
4656 if (mode==
"M" || mode==
"m")
GetSignal(d,a,
"deg",b,
"deg",
"equ",ts,jref,
"M",type);
4657 if (mode==
"A" || mode==
"a")
GetSignal(d,a,
"deg",b,
"deg",
"equ",ts,jref,
"T",type);
4872 MatchSignals(matches,da,au,dt,tu,mode,j,j,0,1,1,1);
4884 MatchSignals(matches,da,au,dt,tu,mode,1,0,0,1,1,1);
4907Double_t
NcAstrolab::GetSeparation(TString name1,TString name2,TString au,Double_t* dt,TString tu,Int_t mode,Double_t* diftheta,Double_t* difphi)
4986 Double_t da=
GetSeparation(i,j,au,dtx,tu,mode,0,diftheta,difphi);
4988Double_t dthetax,dphix;
4990printf(
" i=%-i j=%-i da=%-g dt=%-g %-s dtheta=%-g dphi=%-g %-s \n",i,j,da,dtx,tu.Data(),dthetax,dphix,au.Data());
4996Double_t
NcAstrolab::GetSeparation(Int_t i,Int_t j,TString au,Double_t& dt,TString tu,Int_t mode,Int_t bkgpatch,Double_t* diftheta,Double_t* difphi)
5048 if (diftheta) *diftheta=-999;
5049 if (difphi) *difphi=-999;
5051 if (!i || !j)
return dang;
5052 if ((i>0 || j>0) && !
fRefs)
return dang;
5053 if ((i<0 || j<0) && !
fSigs)
return dang;
5085 if (!si || !sj)
return dang;
5089 if (!ti || !tj)
return dang;
5136 TString name=sxref->GetName();
5153 if (Tscmode==1 || Tscmode==3 || (Tscmode<0 && bkgpatch))
5164 if (Tscmode==3 || (abs(Tscmode)>1 && bkgpatch)) txtmp.
AddSec(trndm);
5168 if (Tscmode==1 || (Tscmode==-1 && bkgpatch))
5172 trndm/=double(24*3600);
5192 name=sxref->GetName();
5196 GetSignal(rxevt,
"loc",
"T",txevt,idxevt,evttype);
5197 GetSignal(rxref,
"loc",
"T",&txtmp,idxref,reftype);
5200 Double_t pi=acos(-1.);
5209 if (Rscmode==1 || (Rscmode==-1 && bkgpatch))
5217 Float_t temp=cosmin;
5221 cosang=
fRan->Uniform(cosmin,cosmax);
5223 if (au==
"deg") dang*=180./pi;
5228 if (au==
"rad") dang*=pi/180.;
5233 if (Rscmode==3 || (abs(Rscmode)>1 && bkgpatch))
5268 Float_t temp=cosmin;
5272 cosang=
fRan->Uniform(cosmin,cosmax);
5273 dtheta=acos(cosang)*180./pi;
5289 if (vec[0]<=0) vec[0]=1e-20;
5302 dtheta=vec[1]-refvec[1];
5303 dphi=vec[2]-refvec[2];
5304 if (diftheta) *diftheta=dtheta;
5305 if (difphi) *difphi=dphi;
5396 MatchSignals(matches,da,au,dt,tu,mode,1,0,0,1,1,1);
5399 if (!nhits)
return 0;
5407 for (Int_t i=1; i<=nhits; i++)
5420 if (!jfill)
return 0;
5425void NcAstrolab::MatchSignals(
NcDevice& matches,Double_t da,TString au,Double_t dt,TString tu,Int_t mode,Int_t i1,Int_t i2,Int_t itype,Int_t j1,Int_t j2,Int_t jtype)
5543 TString name=
"Matches";
5544 TString title=
"Space and time matchings of NcAstrolab stored signals";
5545 matches.SetNameTitle(name.Data(),title.Data());
5547 if (tu==
"d") tux=
"days";
5548 if (tu==
"s") tux=
"sec";
5549 TString namedamin=
"psimin in ";
5551 TString namedtmin=
"dtmin in ";
5559 TString nameda=
"psi in ";
5561 TString namedt=
"t2-t1 in ";
5570 if ((!itype || !jtype) && !
fRefs)
5572 printf(
" *%-s::MatchSignals* Error: itype=%-i jtype=%-i but no reference signals are present. \n",ClassName(),itype,jtype);
5576 if ((itype || jtype) && !
fSigs)
5578 printf(
" *%-s::MatchSignals* Error: itype=%-i jtype=%-i but no measurements are present. \n",ClassName(),itype,jtype);
5593 if (i2<1 || i2>nrefs) i2=nrefs;
5597 if (i2<1 || i2>nsigs) i2=nsigs;
5601 if (j2<1 || j2>nrefs) j2=nrefs;
5605 if (j2<1 || j2>nsigs) j2=nsigs;
5608 if (i1<1 || j1<1 || i1>i2 || j1>j2)
5610 printf(
" *%-s::MatchSignals* Inconsistent parameters: i1=%-i i2=%-i itype=%-i j1=%-i j2=%-i jtype=%-i. \n",ClassName(),i1,i2,itype,j1,j2,jtype);
5614 Double_t dang,dtime;
5624 for (Int_t i=i1; i<=i2; i++)
5629 for (Int_t j=j1; j<=j2; j++)
5632 if (itype==jtype && i==j)
continue;
5639 if ((fabs(dang)<=da || da<0) && (fabs(dtime)<=dt || dt<0))
5645 name+=sx->GetName();
5649 title+=sx->GetName();
5651 data.SetNameTitle(name.Data(),title.Data());
5652 data.SetUniqueID(
id);
5662 if (first || fabs(dang)<dangmin)
5669 if (first || fabs(dtime)<fabs(dtmin))
5687void NcAstrolab::MatchSignals(
NcDevice& matches,TString name,Double_t da,TString au,Double_t dt,TString tu,Int_t mode,Int_t itype,Int_t j1,Int_t j2,Int_t jtype)
5810 printf(
" *%-s::MatchSignals* Object %-s not found for itype=%-i. \n",ClassName(),name.Data(),itype);
5814 MatchSignals(matches,da,au,dt,tu,mode,i,i,itype,j1,j2,jtype);
5903 if (tmax>tmin)
fTscfunc->SetRange(tmin,tmax);
6027 if (abs(mode)==1 && dmin<0) dmin=0;
6028 if (abs(mode)==1 && dmax>180) dmax=180;
6031 if (dmax<dmin) dmax=dmin;
6032 if (thmax<thmin) thmax=thmin;
6033 if (phimax<phimin) phimax=phimin;
6046 if (dmax>dmin)
fDscfunc->SetRange(dmin,dmax);
6070 if (phimax>phimin)
fPhiscfunc->SetRange(phimin,phimax);
6242 sx=
GetSignal(r,frame,mode,ts,jref,type);
6248 TString namesave=sx->GetName();
6255 Double_t pi=acos(-1.);
6257 if (frame==
"equ" || frame==
"gal" || frame==
"icr" || frame==
"ecl" || frame==
"loc")
6259 theta=(pi/2.)-r.
GetX(2,
"sph",
"rad");
6260 phi=r.
GetX(3,
"sph",
"rad");
6265 theta=(pi/2.)-r.
GetX(2,
"sph",
"rad");
6266 phi=pi-r.
GetX(3,
"sph",
"rad");
6277 if (frame==
"gal" || frame==
"icr" || frame==
"ecl")
6282 if (frame==
"hor" || frame==
"loc")
6298 if (proj==
"hamh" || proj==
"aith" || proj==
"merh" || proj==
"cylh" || proj==
"angh") hist=1;
6299 if (proj==
"UTh" || proj==
"LTh" || proj==
"GSTh" || proj==
"LSTh") hist=1;
6300 if (proj==
"UYh" || proj==
"LYh" || proj==
"GSYh" || proj==
"LSYh") hist=1;
6305 if (!
fCanvas || !(gROOT->GetListOfCanvases()->FindObject(
"NcAstrolab")))
fCanvas=
new TCanvas(
"NcAstrolab",
"Skymap");
6323 titleup+=
"Geocentric Equatorial (";
6325 if (mode==
"J") titleup+=
"2000";
6326 if (mode==
"B") titleup+=
"1950";
6329 if (frame==
"gal") titleup+=
"Heliocentric Galactic";
6330 if (frame==
"ecl") titleup+=
" Geocentric Ecliptic";
6331 if (frame==
"hor") titleup+=
" Standard Horizon";
6332 if (frame==
"icr") titleup+=
"Static Barycentric ICRS";
6335 titleup+=
" User defined Local";
6339 titleup+=
" Coordinates";
6340 titlelow=
"Projection : ";
6341 if (proj==
"ham" || proj==
"hamh") titlelow+=
"Hammer";
6342 if (proj==
"cyl" || proj==
"cylh") titlelow+=
"Cylindrical";
6343 if (proj==
"ait" || proj==
"aith") titlelow+=
"Aitoff";
6344 if (proj==
"mer" || proj==
"merh") titlelow+=
"Mercator";
6345 if (proj==
"ang" || proj==
"angh")
6347 titlelow+=
"sin(b) vs. l";
6351 titlelow+=
" Central Meridian : ";
6353 Int_t angmax,hmin,hmax,dmin,dmax;
6459 if (clr==1 || proj!=
fProj)
6481 Float_t xmargin=0.5;
6482 Float_t ymargin=0.3;
6483 fCanvas->Range(xlow-xmargin,ylow-ymargin,xup+xmargin,yup+ymargin);
6486 if (proj==
"ham" || proj==
"ait")
6489 TEllipse* outline=
new TEllipse(0,0,xup,yup);
6499 const Int_t nphi=13;
6500 Double_t gphiarr[nphi]={0,30,60,90,120,150,180,210,240,270,300,330,360};
6504 Float_t gstep=180./float(ndots);
6507 for (Int_t iph=0; iph<nphi; iph++)
6509 gphi=gphiarr[iph]*pi/180.;
6510 if (frame==
"hor") gphi=pi-gphi;
6512 for (Int_t ith=1; ith<ndots; ith++)
6514 gtheta=gtheta-(gstep*pi/180.);
6515 Project(gphi,gtheta,proj,xgrid,ygrid);
6526 Double_t gtharr[nth]={15,30,45,60,75,105,120,135,150,165};
6529 gstep=360./float(ndots);
6535 for (Int_t ith=0; ith<nth; ith++)
6537 gtheta=pi/2.-(gtharr[ith]*pi/180.);
6538 igs=int(90.-gtharr[ith]);
6539 if (frame==
"loc") igs=int(gtharr[ith]);
6544 for (Int_t iphi=1; iphi<ndots; iphi++)
6547 Project(gphi,gtheta,proj,xgrid,ygrid);
6563 if (proj==
"ham" || proj==
"ait")
6565 lgs->DrawLatex(xtext-0.25,ytext,gs.Data());
6569 lgs->DrawLatex(xtext-0.4,ytext-0.02,gs.Data());
6574 if (proj==
"ham" || proj==
"ait")
6576 lgs->DrawLatex(xtext-0.3,ytext-0.1,gs.Data());
6580 lgs->DrawLatex(xtext-0.4,ytext-0.02,gs.Data());
6586 TLine* line=
new TLine(xlow,0,xup,0);
6589 line=
new TLine(0,yup,0,ylow);
6594 TLatex* header=
new TLatex;
6596 header->SetTextAlign(21);
6597 header->DrawLatex(0,yup+0.2,titleup.Data());
6598 TLatex* footer=
new TLatex;
6600 footer->SetTextAlign(21);
6601 footer->DrawLatex(0,ylow-0.25,titlelow.Data());
6604 TLatex* left=
new TLatex;
6606 if (proj==
"ham" || proj==
"ait")
6608 left->DrawLatex(xlow-0.4,0,sleft.Data());
6612 left->DrawLatex(xlow-0.15,yup+0.05,sleft.Data());
6615 TLatex* right=
new TLatex;
6617 if (proj==
"ham" || proj==
"ait")
6619 right->DrawLatex(xup+0.1,0,sright.Data());
6623 right->DrawLatex(xup-0.1,yup+0.05,sright.Data());
6626 TLatex* up=
new TLatex;
6628 if (proj==
"ham" || proj==
"ait")
6630 up->DrawLatex(-0.1,yup+0.05,sup.Data());
6634 up->DrawLatex(-0.1,yup+0.05,scenter.Data());
6639 up->DrawLatex(xlow-0.4,yup-0.04,sup.Data());
6643 TLatex* low=
new TLatex;
6645 if (proj==
"ham" || proj==
"ait")
6647 low->DrawLatex(-0.15,ylow-0.15,slow.Data());
6651 if (proj!=
"ang") low->DrawLatex(xlow-0.4,ylow,slow.Data());
6659 sx=
SetSignal(1,0,
"deg",0,
"deg",
"gal",0,-1,
"J",
"GC",0);
6660 Int_t idx=
fRefs->IndexOf(sx);
6668 if (frame==
"equ" || frame==
"gal" || frame==
"icr" || frame==
"ecl" || frame==
"loc")
6670 thetagc=(pi/2.)-rgc.
GetX(2,
"sph",
"rad");
6671 phigc=rgc.
GetX(3,
"sph",
"rad");
6675 thetagc=(pi/2.)-rgc.
GetX(2,
"sph",
"rad");
6676 phigc=pi-rgc.
GetX(3,
"sph",
"rad");
6681 Project(phigc,thetagc,proj,xgc,ygc);
6704 if (frame==
"equ") xfac=6;
6705 if (proj==
"angh") yfac=1;
6709 if (clr==1 || proj!=
fProj)
6713 for (Int_t i=0; i<2; i++)
6723 fHist[type]->Reset();
6727 TString title=titleup;
6730 fHist[type]->SetNameTitle(
"SkyMap",title.Data());
6731 fHist[type]->GetXaxis()->SetTitle(
"Degrees from central Meridian");
6734 fHist[type]->SetBins(1000,-181,181,100,-1.1,1.1);
6735 fHist[type]->GetYaxis()->SetTitle(
"sin(b)");
6739 fHist[type]->SetBins(1000,-181,181,500,-91,91);
6740 fHist[type]->GetYaxis()->SetTitle(
"Projected Latitude in degrees");
6744 fHist[type]->GetXaxis()->SetTitle(
"Hours from central Meridian");
6747 fHist[type]->SetBins(200,-12.1,12.1,100,-1.1,1.1);
6748 if (frame==
"equ")
fHist[type]->GetYaxis()->SetTitle(
"sin(#delta)");
6752 fHist[type]->SetBins(200,-12.1,12.1,500,-91,91);
6753 if (frame==
"equ")
fHist[type]->GetYaxis()->SetTitle(
"Projected Declination in degrees");
6760 fHist[type]->GetYaxis()->SetTitle(
"sin(alt)=cos(zenith)");
6764 fHist[type]->GetYaxis()->SetTitle(
"Projected Altitude in degrees");
6771 fHist[type]->GetYaxis()->SetTitle(
"cos(#theta)=sin(b)");
6775 fHist[type]->GetYaxis()->SetTitle(
"Projected degrees from the equator");
6783 fHist[type]->Fill(x*xfac,theta*180./pi);
6785 else if (proj==
"hamh" || proj==
"aith" || proj==
"cylh" || proj==
"angh")
6787 fHist[type]->Fill(x*xfac,y*yfac);
6789 else if (proj==
"UTh" || proj==
"LTh" || proj==
"GSTh" || proj==
"LSTh")
6792 fHist[type]->SetBins(100,0,24,181,-90.5,90.5);
6796 TString title=
"Day view";
6809 if (proj==
"LTh") tmode=
"LMT";
6813 if (mode==
"T") tmode=
"GAST";
6818 if (mode==
"T") tmode=
"LAST";
6820 if (proj==
"UTh") title+=
";Universal Time";
6821 if (proj==
"LTh") title+=
";Local Time";
6822 if (proj==
"GSTh") title+=
";Greenwich Sidereal Time";
6823 if (proj==
"LSTh") title+=
";Local Sidereal Time";
6830 title+=
" in hours;";
6831 TString ytitle=
"Declination ";
6832 if (frame==
"equ" && mode==
"J") ytitle+=
"(J2000)";
6833 if (frame==
"equ" && mode==
"B") ytitle+=
"(B1950)";
6834 if (frame==
"equ" && mode==
"M") ytitle+=
"(Mean)";
6835 if (frame==
"equ" && mode==
"T") ytitle+=
"(True)";
6836 if (frame==
"gal") ytitle=
"Galactic latitude";
6837 if (frame==
"ecl") ytitle=
"Geocentric Ecliptic latitude";
6838 if (frame==
"hor") ytitle=
"Horizon altitude";
6839 if (frame==
"icr") ytitle=
"ICRS latitude";
6840 if (frame==
"loc") ytitle=
"Angle w.r.t. local frame equator";
6841 ytitle+=
" in degrees";
6843 fHist[type]->SetTitle(title);
6844 fHist[type]->SetStats(kFALSE);
6849 for (Int_t i=0; i<24; i++)
6851 if (tmode==
"UT") hour=tx.
GetUT();
6852 if (tmode==
"LMT") hour=tx.
GetLT(toffset);
6853 if (tmode==
"GMST") hour=tx.
GetGMST();
6854 if (tmode==
"GAST") hour=tx.
GetGAST();
6855 if (tmode==
"LMST") hour=tx.
GetLMST(toffset);
6856 if (tmode==
"LAST") hour=tx.
GetLAST(toffset);
6859 GetSignal(d,a,
"deg",b,
"deg",frame,&tx,jref,mode,type);
6861 if (frame==
"hor") b=90.-b;
6862 if (frame==
"loc") b=90.-a;
6864 fHist[type]->Fill(hour,b);
6872 else if (proj==
"UYh" || proj==
"LYh" || proj==
"GSYh" || proj==
"LSYh")
6875 fHist[type]->SetBins(1500,0,370,181,-90.5,90.5);
6882 if (proj==
"LYh") tmode=
"LMT";
6886 if (mode==
"T") tmode=
"GAST";
6891 if (mode==
"T") tmode=
"LAST";
6895 if (proj==
"LYh" || proj==
"LSYh") ts->
GetDayTimeString(tmode,0,toffset,0,&time);
6897 TString title=
"Year view";
6909 title+=
";Day of the year;";
6910 TString ytitle=
"Declination ";
6911 if (frame==
"equ" && mode==
"J") ytitle+=
"(J2000)";
6912 if (frame==
"equ" && mode==
"B") ytitle+=
"(B1950)";
6913 if (frame==
"equ" && mode==
"M") ytitle+=
"(Mean)";
6914 if (frame==
"equ" && mode==
"T") ytitle+=
"(True)";
6915 if (frame==
"gal") ytitle=
"Galactic latitude";
6916 if (frame==
"ecl") ytitle=
"Geocentric Ecliptic latitude";
6917 if (frame==
"hor") ytitle=
"Horizon altitude";
6918 if (frame==
"icr") ytitle=
"ICRS latitude";
6919 if (frame==
"loc") ytitle=
"Angle w.r.t. local frame equator";
6920 ytitle+=
" in degrees";
6922 fHist[type]->SetTitle(title);
6923 fHist[type]->SetStats(kFALSE);
6929 if (proj==
"UYh" || proj==
"GSYh")
6931 tx.
SetUT(year.Atoi(),1,1,
"00:00:00");
6934 if (proj==
"LYh" || proj==
"LSYh")
6936 tx.
SetLT(toffset,year.Atoi(),1,1,
"00:00:00");
6937 hour=ts->
GetLT(toffset);
6943 for (Int_t i=0; i<367; i++)
6947 GetSignal(d,a,
"deg",b,
"deg",frame,&tx,jref,mode,type);
6949 if (frame==
"hor") b=90.-b;
6950 if (frame==
"loc") b=90.-a;
6952 if (proj==
"UYh" || proj==
"GSYh") day=tx.GetDayOfYear();
6953 if (proj==
"LYh" || proj==
"LSYh") day=tx.GetDayOfYear(kFALSE,toffset);
6955 fHist[type]->Fill(day,b);
6965 if ((!type &&
fHist[1]) || (type &&
fHist[0]))
6967 fHist[type]->Draw(
"same");
6971 fHist[type]->Draw();
6974 if (proj==
"UTh" || proj==
"LTh" || proj==
"GSTh" || proj==
"LSTh" ||
6975 proj==
"UYh" || proj==
"LYh" || proj==
"GSYh" || proj==
"LSYh")
6983 if (proj==
"UTh" || proj==
"LTh" || proj==
"GSTh" || proj==
"LSTh") line=
new TLine(0,0,24,0);
6984 if (proj==
"UYh" || proj==
"LYh" || proj==
"GSYh" || proj==
"LSYh") line=
new TLine(0,0,370,0);
6987 line->SetLineWidth(3);
7090 GetSignal(d,a,
"deg",b,
"deg",frame,ts,name,mode,type);
7206 if (
fRefs && type<=0)
7216 if (!j || !tx) tx=ts;
7222 for (Int_t i=1; i<=
fRefs->GetSize(); i++)
7228 if (nmax>=0 && jdisp>nmax)
break;
7231 namex=sx->GetName();
7232 if (name!=
"*" && !namex.Contains(name))
continue;
7258 for (Int_t j=1; j<=
fSigs->GetSize(); j++)
7264 if (nmax>=0 && jdisp>nmax)
break;
7267 namex=sx->GetName();
7268 if (name!=
"*" && !namex.Contains(name))
continue;
7303 if (type<0 || type >3)
return;
7322 if (type<0 || type >3)
return;
7341 if (type<0 || type >3)
return;
7376 Double_t pi=acos(-1.);
7377 Double_t twopi=2.*pi;
7424 Double_t pi=acos(-1.);
7447 if (proj==
"ang" || proj==
"angh")
7465 Double_t pi=acos(-1.);
7482 Double_t k=1./sqrt(1.+cos(b)*cos(l/2.));
7483 x=2.*k*cos(b)*sin(l/2.);
7498 Double_t pi=acos(-1.);
7501 Double_t k=acos(cos(b)*cos(l/2.));
7504 x=4.*k*cos(b)*sin(l/2.)/(pi*sin(k));
7505 y=2.*k*sin(b)/(pi*sin(k));
7526 Double_t pi=acos(-1.);
7527 Double_t bcut=85.051*pi/180.;
7531 if (b > bcut) b=bcut;
7532 if (b < -bcut) b=-bcut;
7533 y=0.5*log((1.+sin(b))/(1.-sin(b)))/pi;
7610 if (name==
"Me")
fMe=value;
7611 if (name==
"Mmu")
fMmu=value;
7612 if (name==
"Mtau")
fMtau=value;
7620 if (name==
"Mp")
fMp=value;
7621 if (name==
"Mn")
fMn=value;
7622 if (name==
"MW")
fMW=value;
7623 if (name==
"GammaW")
fGammaW=value;
7624 if (name==
"MZ")
fMZ=value;
7625 if (name==
"GammaZ")
fGammaZ=value;
7626 if (name==
"AlphaEM")
fAlphaEM=value;
7627 if (name==
"Fermi")
fFermi=value;
7636 if (name==
"Boltz")
fBoltz=value;
7643 if (name==
"Gn")
fGn=value;
7644 if (name==
"Au")
fAu=value;
7645 if (name==
"Pc")
fPc=value;
7646 if (name==
"Hubble")
fHubble=value;
7647 if (name==
"OmegaM")
fOmegaM=value;
7648 if (name==
"OmegaR")
fOmegaR=value;
7649 if (name==
"OmegaL")
fOmegaL=value;
7650 if (name==
"OmegaB")
fOmegaB=value;
7651 if (name==
"OmegaC")
fOmegaC=value;
7681 if (name==
"SpeedC")
return fSpeedC;
7682 if (name==
"Qe")
return fQe;
7683 if (name==
"Me")
return fMe;
7684 if (name==
"Mmu")
return fMmu;
7685 if (name==
"Mtau")
return fMtau;
7686 if (name==
"Amu")
return fAmu;
7687 if (name==
"Mp")
return fMp;
7688 if (name==
"Mn")
return fMn;
7689 if (name==
"MW")
return fMW;
7690 if (name==
"GammaW")
return fGammaW;
7691 if (name==
"GammaZ")
return fGammaZ;
7692 if (name==
"AlphaEM")
return fAlphaEM;
7693 if (name==
"Fermi")
return fFermi;
7694 if (name==
"Planck")
return fPlanck;
7695 if (name==
"Boltz")
return fBoltz;
7696 if (name==
"Newton")
return fNewton;
7697 if (name==
"Gn")
return fGn;
7698 if (name==
"Au")
return fAu;
7699 if (name==
"Pc")
return fPc;
7700 if (name==
"Hubble")
return fHubble;
7701 if (name==
"OmegaM")
return fOmegaM;
7702 if (name==
"OmegaR")
return fOmegaR;
7703 if (name==
"OmegaL")
return fOmegaL;
7704 if (name==
"OmegaB")
return fOmegaB;
7705 if (name==
"OmegaC")
return fOmegaC;
7708 if (name==
"Hbar")
return fHbar;
7709 if (name==
"Hbarc")
return fHbarc;
7710 if (name==
"Hbarc2")
return fHbarc2;
7721 if (name==
"Jy")
return 1e-23;
7768 if (z<=0 ||
fHubble<=0)
return 0;
7772 TF1 f(
"f",
"1./sqrt([0]*pow((1.+x),4)+[1]*pow((1.+x),3)+[2])");
7779 Double_t dist=f.Integral(0,z);
7782 Double_t distm=dist*1e6*
fPc;
7786 if (u==
"Gpc") val=dist*1e-3;
7787 if (u==
"Mpc") val=dist;
7788 if (u==
"pc") val=dist*1e6;
7789 if (u==
"ly") val=dist*3.26156e6;
7790 if (u==
"m") val=distm;
7791 if (u==
"km") val=distm*1e-3;
7792 if (u==
"cm") val=distm*1e2;
7794 if (!t) val=val/(z+1.);
7956 if (z<=0 ||
fHubble<=0)
return 0;
7960 TF1 f(
"f",
"1./((1.+x)*sqrt([0]*pow((1.+x),4)+[1]*pow((1.+x),3)+[2]))");
7967 Double_t dist=f.Integral(0,z);
7970 Double_t distm=dist*1e6*
fPc;
7974 if (u==
"Gpc") val=dist*1e-3;
7975 if (u==
"Mpc") val=dist;
7976 if (u==
"pc") val=dist*1e6;
7977 if (u==
"ly") val=dist*3.26156e6;
7978 if (u==
"m") val=distm;
7979 if (u==
"km") val=distm*1e-3;
7980 if (u==
"cm") val=distm*1e2;
8033 if (z<0 ||
fHubble<=0)
return 0;
8035 TF1 f(
"f",
"sqrt([0]*pow((1.+x),4)+[1]*pow((1.+x),3)+[2])");
8042 Double_t H=f.Eval(z);
8045 Double_t Hm=H/(1e6*
fPc);
8049 if (u==
"Gpc") val=H/1e-3;
8050 if (u==
"Mpc") val=H;
8051 if (u==
"pc") val=H/1e6;
8052 if (u==
"ly") val=H/3.26156e6;
8054 if (u==
"km") val=Hm/1e-3;
8055 if (u==
"cm") val=Hm/1e2;
8090 if (Z<0 || N<0)
return 0;
8104 Double_t term1=a*ra;
8105 Double_t term2=b*pow(ra,2./3.);
8106 Double_t term3=s*pow((rn-rz),2.)/ra;
8107 Double_t term4=d*pow(rz,2.)/pow(ra,1./3.);
8113 if (oz && on) term5=delta/sqrt(ra);
8114 if (!oz && !on) term5=-delta/sqrt(ra);
8117 Double_t bnz=term1-term2-term3-term4-term5;
8127 Double_t mass=rz*
fMp+rn*
fMn-bnz;
8132 mass=2.013553212712*
fAmu;
8138 mass=3.0155007134*
fAmu;
8144 mass=3.0149322468*
fAmu;
8150 mass=4.001506179125*
fAmu;
8179 value=bnz/(
fAmu*ra);
8219 if (Z<=0 || A<1 || A<Z)
return -1;
8221 X0=716.4*A/(Z*(Z+1.)*(log(287./sqrt(Z))));
8224 if (rho==0) X0=X0/mN;
8225 if (rho>0) X0=X0/rho;
8269 if (sigma<=0 || rho<=0 || mode<0 || mode>4)
return -1;
8272 Double_t n=fabs(rho)/mN;
8276 lambda=1./(sigma*n);
8278 if (mode==1 || mode==4) lambda=lambda*n;
8279 if (mode==3) lambda=lambda*rho;
8319 if (x<0 || lambda<=0)
return -1;
8321 prob=1.-exp(-x/lambda);
8404 if (x<0 || lambda<=0)
return -1;
8406 prob=exp(-x/lambda);
8486 if (prob<=0 || prob>1 || lambda<=0)
return -1;
8488 x=-lambda*log(prob);
8526 if (prob<=0 || prob>1)
return -1;
8565 if (prob<=0 || prob>1 || lambda<=0)
return -1;
8606 if (prob<=0 || prob>1)
return -1;
8665 if (eprimgev) *eprimgev=0;
8666 if (alpha) *alpha=0;
8667 if (!mode || mode>3 || mode<-4 || !type || abs(type)>3)
return 0;
8669 const Double_t fnumuccn=6.77e-15;
8670 const Double_t fanumuccn=3.34e-15;
8674 const Double_t fnuetote=0.25+sinw2+4.*pow(sinw2,2.)/3.;
8675 const Double_t fanuetote=(1./12.)+(sinw2/3.)+4.*pow(sinw2,2.)/3.;
8676 Double_t fnumucce=1.;
8677 const Double_t fnumunce=0.25-sinw2+4.*pow(sinw2,2.)/3.;
8678 const Double_t fanumunce=(1./12.)-(sinw2/3.)+4.*pow(sinw2,2.)/3.;
8679 const Double_t f4=1./3.;
8682 const Double_t c0nu=-1.826;
8683 const Double_t c1nu=-17.31;
8684 const Double_t c2nunc=-6.448;
8685 const Double_t c2nucc=-6.406;
8686 const Double_t c3nu=1.431;
8687 const Double_t c4nunc=-18.61;
8688 const Double_t c4nucc=-17.91;
8689 const Double_t c0anu=-1.033;
8690 const Double_t c1anu=-15.95;
8691 const Double_t c2anunc=-7.296;
8692 const Double_t c2anucc=-7.247;
8693 const Double_t c3anu=1.569;
8694 const Double_t c4anunc=-18.30;
8695 const Double_t c4anucc=-17.72;
8697 Double_t rncnu=0.2261/0.7221;
8698 Double_t rncanu=0.1307/0.3747;
8701 Double_t ynucc[12]={0.483,0.477,0.472,0.426,0.332,0.237,0.250,0.237,0.225,0.216,0.208,0.205};
8702 Double_t ynunc[12]={0.474,0.470,0.467,0.428,0.341,0.279,0.254,0.239,0.227,0.217,0.210,0.207};
8703 Double_t yanucc[12]={0.333,0.340,0.354,0.345,0.301,0.266,0.249,0.237,0.225,0.216,0.208,0.205};
8704 Double_t yanunc[12]={0.350,0.354,0.368,0.358,0.313,0.273,0.253,0.239,0.227,0.217,0.210,0.207};
8706 Double_t loge=log10(egev);
8708 Int_t index=int(loge+0.5);
8709 if (index<1) index=1;
8710 if (index>12) index=12;
8715 if (abs(mode)==1) y=ynucc[index-1];
8716 if (abs(mode)==2) y=ynunc[index-1];
8717 if (mode==3) y=(ynucc[index-1]+rncnu*ynunc[index-1])/(1.+rncnu);
8718 if (mode==-3) y=(ynucc[index-1]+fnumunce*ynunc[index-1])/(1.+fnumunce);
8722 if (abs(mode)==1) y=yanucc[index-1];
8723 if (abs(mode)==2) y=yanunc[index-1];
8724 if (mode==3) y=(yanucc[index-1]+rncanu*yanunc[index-1])/(1.+rncanu);
8725 if (mode==-3) y=(yanucc[index-1]+fnumunce*yanunc[index-1])/(1.+fnumunce);
8727 *eprimgev=egev*(1.-y);
8734 *alpha=sqrt(2.e-3*mtarg/egev)*y*180./((1.-y)*acos(-1.));
8739 Double_t c0=0,c1=0,c2=0,c3=0,c4=0;
8755 if (type<0) fact=fanumuccn;
8759 fact=fnumuccn*rncnu;
8760 if (type<0) fact=fanumuccn*rncanu;
8804 Double_t lne=log(loge-c0);
8805 Double_t logsigma=c1+c2*lne+c3*pow(lne,2.)+c4/lne;
8807 xsec=pow(10.,logsigma);
8815 if (mode==-3 && type==-1 && egev>elow && egev<eup)
8817 xsec=5.02e-7/xscale;
8818 if (eprimgev) *eprimgev=0;
8819 if (alpha) *alpha=0;
8824 Double_t mlepton=
fMe;
8825 if (abs(type)==2) mlepton=
fMmu;
8826 if (abs(type)==3) mlepton=
fMtau;
8827 Double_t eth=1.e-3*(pow(mlepton,2.)-pow(
fMe,2.))/(2.*
fMe);
8834 if (eprimgev) *eprimgev=egev;
8835 if (alpha) *alpha=0;
8845 if (type>1) fact=fnumucce;
8849 if (type>1) fact=fnumunce;
8850 if (type<-1) fact=fanumunce;
8854 if (type==1) fact=fnuetote;
8855 if (type==-1) fact=fanuetote;
8856 if (type>1) fact=fnumucce+fnumunce;
8857 if (type<-1) fact=fanumunce;
8859 if (mode==-4 && type==-1) fact=f4;
8908 if (E<=0 || mode<0 || mode>2)
return value;
8914 Double_t mean=1.38711583;
8915 Double_t median=0.86842105;
8916 Double_t mpv=0.560150;
8917 Double_t sigma=0.226679;
8920 Double_t p=log10(E);
8921 Double_t scale=1./(pow(1.5,p)*sqrt(E));
8931 fNuAngle=
new TF1(
"NuAngle",
"TMath::Landau(x,[0],[1],1)",0,90);
8932 fNuAngle->SetTitle(
"Landau pdf;Neutrino-lepton opening angle in degrees;PDF");
8940 Double_t ang=
fNuAngle->GetRandom();
8944 Double_t fact=acos(-1.)/180.;
8954 if (mode==1) value=median;
8955 if (mode==2) value=ang;
8986 Double_t pi=acos(-1.);
8987 Double_t cosmax=cos(thetamin*pi/180.);
8988 Double_t cosmin=cos(thetamax*pi/180.);
8989 Double_t cost=
fRan->Uniform(cosmin,cosmax);
8990 Double_t theta=acos(cost)*180./pi;
8991 Double_t phi=
fRan->Uniform(phimin,phimax);
8996 Double_t err[3]={0,0,0};
9004 v.
SetVector(norm,theta,phi,
"sph",
"deg");
9038 Double_t norm=v.
GetX(1,
"sph",
"deg");
9039 Double_t theta=v.
GetX(2,
"sph",
"deg");
9040 Double_t phi=v.
GetX(3,
"sph",
"deg");
9041 Double_t err[3]={0,0,0};
9053 v.
SetVector(norm,theta,phi,
"sph",
"deg");
9060 m.SetAngles(90.+theta,phi,90,phi+90.,theta,phi);
9063 Double_t pi=acos(-1.);
9065 Double_t cosmin=cos(fabs(sigma)*pi/180.);
9067 Double_t phimax=360;
9071 cost=
fRan->Uniform(cosmin,cosmax);
9072 theta=acos(cost)*180./pi;
9076 theta=
fRan->Gauss(0.,sigma);
9078 phi=
fRan->Uniform(phimin,phimax);
9081 v.
SetVector(norm,theta,phi,
"sph",
"deg");
9116 Double_t norm=v.
GetX(1,
"sph",
"deg");
9117 Double_t theta=v.
GetX(2,
"sph",
"deg");
9118 Double_t phi=v.
GetX(3,
"sph",
"deg");
9119 Double_t err[3]={0,0,0};
9131 v.
SetVector(norm,theta,phi,
"sph",
"deg");
9138 m.SetAngles(90.+theta,phi,90,phi+90.,theta,phi);
9142 Double_t phimax=360;
9144 phi=
fRan->Uniform(phimin,phimax);
9147 v.
SetVector(norm,theta,phi,
"sph",
"deg");
9256 if (mode<0 || mode>3)
return hdx;
9258 if (!hx)
return hdx;
9260 if (nc<1)
return hdx;
9262 Int_t nenhx=hx->GetEntries();
9263 if (nenhx<=nc)
return hdx;
9265 Int_t idxbin=TMath::Nint(dxbin);
9266 if (idxbin<-2)
return hdx;
9271 if (dxmin>=0 && dxmax>=dxmin && dxbin>0)
9274 Double_t range=dxmax-dxmin;
9275 if (range>dxbin) nbins=TMath::Nint(range/dxbin);
9276 hdx.SetBins(nbins,dxmin,dxmax);
9280 Double_t binwidth=hdx.GetXaxis()->GetBinWidth(1);
9281 s.Form(
"Dx interval distribution between %-i consecutive entries (nc=%-i, mode=%-i);Dx interval;Counts per bin of size %-.3g",nc+1,nc,mode,binwidth);
9282 hdx.SetNameTitle(
"DxHistogram",s);
9293 Double_t x1,x2,deltax;
9295 Double_t deltaxmin=0;
9296 Double_t deltaxmax=0;
9297 Bool_t found=kFALSE;
9301 Int_t nbhx=hx->GetNbinsX();
9306 for (Int_t i=1; i<=nbhx; i++)
9310 xlow=hx->GetBinLowEdge(i);
9311 bsize=hx->GetBinWidth(i);
9313 x1=hx->GetBinCenter(i);
9314 if (mode==1 || mode==3) x1=
fRan->Uniform(xlow,xup);
9315 value=hx->GetBinContent(i);
9318 if (mode<2) nx1=TMath::Nint(value);
9324 if (nx1>1) jstart=i;
9326 for (Int_t j=jstart; j<=nbhx; j++)
9328 xlow=hx->GetBinLowEdge(j);
9329 bsize=hx->GetBinWidth(j);
9331 x2=hx->GetBinCenter(j);
9332 if (mode==1 || mode==3) x2=
fRan->Uniform(xlow,xup);
9333 value=hx->GetBinContent(j);
9336 if (mode<2) nx2=TMath::Nint(value);
9338 if (j==i) nx2=nx1-1;
9341 if (nx2<1)
continue;
9348 if (dxmin>=0 && dxmax>=dxmin && dxbin>0)
9354 if (!found || deltax<deltaxmin) deltaxmin=deltax;
9355 if (!found || deltax>deltaxmax) deltaxmax=deltax;
9367 if (!found)
return hdx;
9370 Int_t nen=hdx.GetEntries();
9374 if (!idxbin && dxbin<=0) dxbin=hx->GetBinWidth(1);
9377 dxbin=(hx->GetBinWidth(1))*fact;
9378 if (deltaxmin>0 && deltaxmin>dxbin) dxbin=deltaxmin;
9379 if (dxbin<=0) dxbin=hx->GetBinWidth(1);
9383 dxbin=hx->GetBinWidth(1);
9384 dxbin=dxbin*float(nc);
9392 if (mode==1 || mode==3)
9394 bsize=hx->GetBinWidth(1);
9395 dxmin=dxmin-2.*bsize;
9396 if (dxmin<0) dxmin=0;
9401 dxmax=deltaxmax+dxbin;
9403 if (mode==1 || mode==3)
9405 bsize=hx->GetBinWidth(1);
9406 dxmax=dxmax+2.*bsize;
9456 if (!hin)
return hout;
9458 Int_t nbins=hin->GetNbinsX();
9459 if (!nbins)
return hout;
9462 const TArrayD* xarr=hin->GetXaxis()->GetXbins();
9463 Int_t xsize=xarr->GetSize();
9466 Double_t xmin=hin->GetXaxis()->GetXmin();
9467 Double_t xmax=hin->GetXaxis()->GetXmax();
9468 hout.SetBins(nbins,xmin,xmax);
9472 const Double_t* xbins=xarr->GetArray();
9473 hout.SetBins(nbins,xbins);
9477 hout.SetNameTitle(
"DifHistogram",hin->GetTitle());
9480 TString sxin=hin->GetXaxis()->GetTitle();
9481 TString syin=hin->GetYaxis()->GetTitle();
9490 syout=f->GetExpFormula(
"p");
9501 sxin.ReplaceAll(
"^{10}log",
"");
9502 sxin.ReplaceAll(
"^{10}Log",
"");
9503 sxin.ReplaceAll(
"log10",
"");
9504 sxin.ReplaceAll(
"Log10",
"");
9505 sxin.ReplaceAll(
"log",
"");
9506 sxin.ReplaceAll(
"Log",
"");
9507 sxin.ReplaceAll(
"ln",
"");
9508 sxin.ReplaceAll(
"Ln",
"");
9513 syout.ReplaceAll(
"((",
"(");
9514 syout.ReplaceAll(
"))",
")");
9517 hout.GetXaxis()->SetTitle(sxout.Data());
9518 hout.GetYaxis()->SetTitle(syout.Data());
9528 for (Int_t i=1; i<=nbins; i++)
9530 x=hin->GetBinCenter(i);
9531 y=hin->GetBinContent(i);
9532 err=fabs(hin->GetBinError(i));
9533 width=hin->GetBinWidth(i);
9534 binlow=hin->GetBinLowEdge(i);
9538 if (width<=0)
continue;
9544 width=pow(10.,binup)-pow(10.,binlow);
9549 width=exp(binup)-exp(binlow);
9563 hout.SetBinContent(i,y);
9564 hout.SetBinError(i,err);
9600 hout.SetName(
"CountsHistogram");
9605 s=
"CountsHistogram;";
9614 s+=spec.GetXaxis()->GetTitle();
9618 hout.SetTitle(s.Data());
9621 Double_t step=(xmax-xmin)/
double(nbins);
9622 Double_t* xbins=
new Double_t[nbins+1];
9624 for (Int_t ibin=0; ibin<=nbins; ibin++)
9630 hout.SetBins(nbins,xbins);
9637 for (Int_t ibin=1; ibin<=nbins; ibin++)
9642 dx=xbins[ibin]-xbins[ibin-1];
9648 dx=pow(10.,xbins[ibin])-pow(10.,xbins[ibin-1]);
9649 N=spec.Eval(xval)*dx;
9654 dx=exp(xbins[ibin])-exp(xbins[ibin-1]);
9655 N=spec.Eval(xval)*dx;
9699 hout.SetName(
"CountsHistogram");
9700 Int_t nbins=hin.GetNbinsX();
9702 if (nbins<1)
return hout;
9704 TAxis* ax=hin.GetXaxis();
9705 Double_t xmin=ax->GetXmin();
9706 Double_t xmax=ax->GetXmax();
9707 hout.SetBins(nbins,xmin,xmax);
9714 TString tx=ax->GetTitle();
9721 if (mode==0) s+=
"x";
9722 if (mode==1) s+=
"^{10}Log(x)";
9723 if (mode==2) s+=
"Ln(x)";
9737 for (Int_t ibin=1; ibin<=nbins; ibin++)
9739 x=hin.GetBinCenter(ibin);
9740 xlow=hin.GetBinLowEdge(ibin);
9741 xup=xlow+hin.GetBinWidth(ibin);
9742 N=hin.GetBinContent(ibin);
9752 dx=pow(10.,xup)-pow(10.,xlow);
9758 dx=exp(xup)-exp(xlow);
9765 fval=fscale->Eval(xval);
9799 if (!hin || mode<1 || mode>2)
return hout;
9801 Int_t nbins=hin->GetNbinsX();
9802 if (!nbins)
return hout;
9805 const TArrayD* xarr=hin->GetXaxis()->GetXbins();
9806 Int_t xsize=xarr->GetSize();
9809 Double_t xmin=hin->GetXaxis()->GetXmin();
9810 Double_t xmax=hin->GetXaxis()->GetXmax();
9811 hout.SetBins(nbins,xmin,xmax);
9815 const Double_t* xbins=xarr->GetArray();
9816 hout.SetBins(nbins,xbins);
9820 hout.SetNameTitle(
"LogHistogram",hin->GetTitle());
9825 if (mode==2) s=
"Ln(";
9826 s+=hin->GetYaxis()->GetTitle();
9830 hout.GetXaxis()->SetTitle(hin->GetXaxis()->GetTitle());
9831 hout.GetYaxis()->SetTitle(s.Data());
9837 for (Int_t i=1; i<=nbins; i++)
9839 y=hin->GetBinContent(i);
9840 err=fabs(hin->GetBinError(i));
9857 hout.SetBinContent(i,y);
9859 hout.SetBinError(i,err);
9900 Int_t Noff=int(pars[0]);
9901 Double_t Toff=pars[1];
9902 Double_t bmax=pars[2];
9903 Double_t prec=pars[3];
9905 if (b<=0 || Noff<0 || Toff<=0)
return 0;
9907 Double_t rNoff=Noff;
9908 if (bmax<0) bmax=100.*rNoff/Toff;
9918 lnU=log(Toff)+rNoff*log(b*Toff)-b*Toff;
9921 lnD=math.
LnGamma(Noff+1,bmax*Toff,1);
9925 if (lnprob < -fabs(prec))
return 0;
9927 if (lnprob > fabs(prec)) lnprob=fabs(prec);
9979 Int_t Non=TMath::Nint(pars[0]);
9980 Double_t Ton=pars[1];
9981 Int_t Noff=TMath::Nint(pars[2]);
9982 Double_t Toff=pars[3];
9983 Double_t smax=pars[4];
9984 Double_t bmax=pars[5];
9985 Double_t prec=pars[6];
9987 if (s<0 || Non<0 || Ton<=0 || Noff<0 || Toff<=0)
return 0;
9990 if (smax<0) smax=100.*rNon/Ton;
9992 Double_t rNoff=Noff;
9993 if (bmax<0) bmax=100.*rNoff/Toff;
9998 Int_t ndim=Non+Noff+1;
9999 TArrayD lnfacN(ndim);
10003 for (Int_t i=1; i<ndim; i++)
10014 Double_t gammaP1=0;
10015 Double_t gammaP2=0;
10018 for(Int_t i=0; i<=Non; i++)
10023 gammaP1=math.
Gamma(Non+Noff+1-i,bmax*(Ton+Toff),0);
10024 gammaP2=math.
Gamma(i+1,smax*Ton,0);
10027 lnU=-s*Ton+ri*log(s)+ri*log(Ton+Toff)-lnfacN[i]-lnfacN[Non-i]+lnfacN[Non+Noff-i]-lnfacN[Non+Noff]+lnfacN[Non];
10029 if ((lnU > -fabs(prec)) && (lnU < fabs(prec))) sumU+=exp(lnU)*gammaP1;
10032 lnD=ri*log(Ton+Toff)-(ri+1.)*log(Ton)-lnfacN[i]-lnfacN[Non-i]+lnfacN[Non+Noff-i]+lnfacN[i]-lnfacN[Non+Noff]+lnfacN[Non];
10034 if ((lnD > -fabs(prec)) && (lnD < fabs(prec))) sumD+=exp(lnD)*gammaP1*gammaP2;
10037 if (sumD) prob=sumU/sumD;
10074 if (bmax<0) bmax=100.*Noff/Toff;
10079 pdf.SetParName(0,
"Noff");
10080 pdf.SetParName(1,
"Toff");
10081 pdf.SetParName(2,
"bmax");
10082 pdf.SetParName(3,
"prec");
10084 pdf.SetParameter(
"Noff",Noff);
10085 pdf.SetParameter(
"Toff",Toff);
10086 pdf.SetParameter(
"bmax",bmax);
10087 pdf.SetParameter(
"prec",prec);
10089 pdf.SetTitle(
"Bayesian posterior background rate PDF;Background rate B in Hz;p(B|Noff,Toff,I)");
10090 pdf.SetRange(0,bmax);
10135 if (smax<0) smax=100.*Non/Ton;
10141 if (bmax<0) bmax=100.*Noff/Toff;
10146 pdf.SetParName(0,
"Non");
10147 pdf.SetParName(1,
"Ton");
10148 pdf.SetParName(2,
"Noff");
10149 pdf.SetParName(3,
"Toff");
10150 pdf.SetParName(4,
"smax");
10151 pdf.SetParName(5,
"bmax");
10152 pdf.SetParName(6,
"prec");
10154 pdf.SetParameter(
"Non",Non);
10155 pdf.SetParameter(
"Ton",Ton);
10156 pdf.SetParameter(
"Noff",Noff);
10157 pdf.SetParameter(
"Toff",Toff);
10158 pdf.SetParameter(
"smax",smax);
10159 pdf.SetParameter(
"bmax",bmax);
10160 pdf.SetParameter(
"prec",prec);
10162 pdf.SetTitle(
"Bayesian posterior signal rate PDF;Signal rate S in Hz;p(S|Non,Ton,Noff,Toff,I)");
10163 pdf.SetRange(0,smax);
10187 if (p<=0 || p>100)
return 0;
10195 nu=pdf.GetQuantiles(1,ua,xa);
10214 if (p<=0 || p>100 || !his)
return 0;
10222 his->ComputeIntegral();
10225 nu=his->GetQuantiles(1,ua,xa);
10269 if (p<=0 || p>100 || n<2)
return 0;
10272 Double_t prec=1./double(n);
10275 Double_t* q=
new Double_t[n];
10276 Double_t* sumq=
new Double_t[n];
10278 for (Int_t i=0; i<n; i++)
10283 Int_t ncalc=pdf.GetQuantiles(n,q,sumq);
10295 Double_t xmode=pdf.GetMaximumX();
10297 Double_t diff,diffmin;
10298 diffmin=fabs(q[ncalc-1]-q[0]);
10299 for (Int_t i=0; i<ncalc; i++)
10301 diff=fabs(xmode-q[i]);
10310 Double_t xmin=q[0];
10311 Double_t xmax=q[ncalc-1];
10312 Double_t totint=pdf.Integral(xmin,xmax);
10327 Double_t ylow=pdf.Eval(q[ilow]);
10328 Double_t yup=pdf.Eval(q[iup]);
10329 Double_t frac=p/100.;
10330 if (frac>1) frac=1;
10331 Double_t credint=-1;
10332 while (credint<frac*totint)
10334 if (yup>ylow && iup<(ncalc-1))
10340 else if (ylow>yup && ilow>0)
10344 ylow=pdf.Eval(xlow);
10346 else if (iup<(ncalc-1))
10356 ylow=pdf.Eval(xlow);
10362 credint=pdf.Integral(xlow,xup);
10366 Double_t intfrac=credint/totint;
10446 if (p<=0 || p>100 || !his)
return 0;
10448 Int_t nbins=his->GetNbinsX();
10451 if (nbins<2)
return 0;
10454 his->ComputeIntegral();
10458 Double_t* q=
new Double_t[n];
10459 Int_t ncalc=his->GetQuantiles(n,q);
10470 Int_t imode=his->GetMaximumBin();
10473 Double_t totint=his->Integral(1,ncalc,
"width");
10487 Double_t ylow=his->GetBinContent(ilow);
10488 Double_t yup=his->GetBinContent(iup);
10489 Double_t frac=p/100.;
10490 if (frac>1) frac=1;
10491 Double_t credint=-1;
10492 while (credint<frac*totint)
10494 if (yup>ylow && iup<(ncalc-1))
10498 yup=his->GetBinContent(iup);
10500 else if (ylow>yup && ilow>0)
10504 ylow=his->GetBinContent(ilow);
10506 else if (iup<(ncalc-1))
10510 yup=his->GetBinContent(iup);
10516 ylow=his->GetBinContent(ilow);
10522 credint=his->Integral(ilow,iup,
"width");
10526 Double_t intfrac=credint/totint;
10665 if (!mode.Contains(
"M") && !mode.Contains(
"K") && !mode.Contains(
"P"))
return -1;
10666 if (mode.Contains(
"M") && (mode.Contains(
"K") || mode.Contains(
"P")))
return -1;
10667 if (mode.Contains(
"K") && (mode.Contains(
"M") || mode.Contains(
"P")))
return -1;
10668 if (mode.Contains(
"P") && (mode.Contains(
"M") || mode.Contains(
"K")))
return -1;
10670 if (!h1)
return -1;
10671 if (!h2 && !pdf)
return -1;
10672 if (h2 && pdf)
return -1;
10674 ULong64_t nrep=ULong64_t(nr);
10680 nrep=ULong64_t(1.e19);
10688 TAxis* xaxis=h1->GetXaxis();
10689 Double_t xmin1=xaxis->GetXmin();
10690 Double_t xmax1=xaxis->GetXmax();
10691 Double_t range1=xmax1-xmin1;
10692 Int_t nbins1=h1->GetNbinsX();
10693 Double_t nen1=h1->GetSumOfWeights();
10694 Double_t underflow1=h1->GetBinContent(0);
10695 Double_t overflow1=h1->GetBinContent(nbins1+1);
10696 if (mode.Contains(
"U")) nen1=nen1+underflow1;
10697 if (mode.Contains(
"O")) nen1=nen1+overflow1;
10699 if (nbins1<=0 || nen1<=0 || range1<=0)
10701 cout <<
" *" << ClassName() <<
"::KolmogorovTest* Histogram h1 is empty or has inconsistent data." << endl;
10702 cout <<
" h1 : nentries=" << nen1 <<
" nbins=" << nbins1 <<
" xmin=" << xmin1 <<
" xmax=" << xmax1 << endl;
10708 xaxis=h2->GetXaxis();
10709 Double_t xmin2=xaxis->GetXmin();
10710 Double_t xmax2=xaxis->GetXmax();
10711 Double_t range2=xmax2-xmin2;
10712 Int_t nbins2=h2->GetNbinsX();
10713 Double_t nen2=h2->GetSumOfWeights();
10715 if (nen2<=0 || range2<=0)
10717 cout <<
" *" << ClassName() <<
"::KolmogorovTest* Histogram h2 is empty or has inconsistent data." << endl;
10718 cout <<
" h2 : nentries=" << nen2 <<
" nbins=" << nbins2 <<
" xmin=" << xmin2 <<
" xmax=" << xmax2 << endl;
10722 Double_t prec=1e-6;
10723 if (nbins2!=nbins1 || fabs(xmin2-xmin1)>prec || fabs(xmax2-xmax1)>prec)
10725 cout <<
" *" << ClassName() <<
"::KolmogorovTest* Histograms h1 and h2 do not have the same binning." << endl;
10726 cout <<
" h1 : nbins=" << nbins1 <<
" xmin=" << xmin1 <<
" xmax=" << xmax1 << endl;
10727 cout <<
" h2 : nbins=" << nbins2 <<
" xmin=" << xmin2 <<
" xmax=" << xmax2 << endl;
10735 pdf->SetRange(xmin1,xmax1);
10736 pdf->SetNpx(nbins1);
10737 h2=(TH1*)pdf->GetHistogram()->Clone();
10738 h2->SetName(
"hpdf");
10740 for (Int_t i=0; i<=nbins1+1; i++)
10742 h2->SetBinError(i,0);
10748 if (mode.Contains(
"U")) s+=
"U";
10749 if (mode.Contains(
"O")) s+=
"O";
10750 if (mode.Contains(
"N") && !pdf) s+=
"N";
10755 Double_t d0=h2->KolmogorovTest(h1,s2.Data());
10758 if (mode.Contains(
"M")) s+=
"M";
10759 if (mode.Contains(
"I")) s+=
"D";
10762 if (mode.Contains(
"I"))
10766 cout <<
" *" << ClassName() <<
"::KolmogorovTest* Single sample KS-test results for execution mode "<< mode.Data() << endl;
10767 if (mode.Contains(
"N")) cout <<
" === For a single sample KS-test the mode=N is suppressed ===" << endl;
10771 cout <<
" *" << ClassName() <<
"::KolmogorovTest* Two sample KS-test results for execution mode "<< mode.Data() << endl;
10774 value=h1->KolmogorovTest(h2,s.Data());
10777 if (ksh) ksh->SetBins(101,0,1.01);
10783 if (mode.Contains(
"P"))
10785 htemp=(TH1*)h1->Clone();
10786 for (jrep=0; jrep<nrep; jrep++)
10789 for (Int_t ien=0; ien<nen1; ien++)
10791 xval=h2->GetRandom();
10794 dist=htemp->KolmogorovTest(h2,s2.Data());
10795 if (ksh) ksh->Fill(dist);
10797 if (dist>=d0) sumd++;
10800 if (ncut && sumd>=ncut)
break;
10803 value=double(sumd)/sumrep;
10804 if (nrx) *nrx=sumrep;
10805 if (mode.Contains(
"I"))
10807 cout <<
" P-value = " << value <<
" after " << sumrep <<
" pseudo experiments." << endl;
10811 if (mode.Contains(
"I")) cout <<
" Returned value = " << value << endl;
10816 TString xlabel=
"Dmax";
10817 TString ylabel=
"Counts after ";
10819 ylabel+=
" pseudo experiments";
10821 ksh->SetTitle(
"KS-test Dmax distribution from pseudo experiments");
10822 ksh->SetXTitle(xlabel.Data());
10823 ksh->SetYTitle(ylabel.Data());
10831 Float_t ymax=ksh->GetMaximum();
10833 TLine* vline=
new TLine(x,ymin,x,ymax);
10834 vline->SetLineStyle(2);
10835 vline->SetLineWidth(2);
10836 vline->SetLineColor(4);
10838 TString title=
"P-value : %-10.3g";
10839 TString sh=title.Format(title.Data(),value);
10841 TLegend* leg=
new TLegend(0.6,0.8,0.8,0.9);
10842 leg->SetFillColor(0);
10843 leg->SetHeader(sh.Data());
10844 leg->AddEntry(vline,
"Observed Dmax",
"L");
10846 TList* hlist=ksh->GetListOfFunctions();
10859 if (htemp)
delete htemp;
10894 TString title=
"Cumulative Distribution of histogram ";
10895 hcd.SetNameTitle(name.Data(),title.Data());
10897 if (!h)
return hcd;
10899 TAxis* xaxis=h->GetXaxis();
10900 TAxis* yaxis=h->GetYaxis();
10901 Double_t xmin=xaxis->GetXmin();
10902 Double_t xmax=xaxis->GetXmax();
10903 Double_t range=xmax-xmin;
10904 Int_t nbins=h->GetNbinsX();
10905 Double_t nen=h->GetSumOfWeights();
10906 TString nameh=h->GetName();
10907 TString xtitle=xaxis->GetTitle();
10908 TString ytitle=yaxis->GetTitle();
10910 hcd.SetNameTitle(name.Data(),title.Data());
10911 hcd.SetXTitle(xtitle.Data());
10912 hcd.SetYTitle(ytitle.Data());
10914 if (nbins<=0 || nen<=0 || range<=0)
return hcd;
10916 if(!(mode.Contains(
"F") || mode.Contains(
"B")) || (mode.Contains(
"F") && mode.Contains(
"B")))
return hcd;
10918 hcd.SetBins(nbins,xmin,xmax);
10920 if (mode.Contains(
"N")) title=
"Normalized ";
10921 if (mode.Contains(
"F")) title+=
"Forward ";
10922 if (mode.Contains(
"B")) title+=
"Backward ";
10923 title+=
"Cumulative Distribution of histogram ";
10925 hcd.SetNameTitle(name.Data(),title.Data());
10926 hcd.SetXTitle(xtitle.Data());
10927 hcd.SetYTitle(ytitle.Data());
10930 if (mode.Contains(
"N")) norm=nen;
10934 if (mode.Contains(
"F"))
10936 for (Int_t ibin=1; ibin<=nbins; ibin++)
10938 y=h->GetBinContent(ibin);
10940 hcd.SetBinContent(ibin,sum);
10945 for (Int_t ibin=nbins; ibin>=1; ibin--)
10947 y=h->GetBinContent(ibin);
10949 hcd.SetBinContent(ibin,sum);
10986 TString title=
"Cumulative Distribution Histogram of function ";
10987 hcd.SetNameTitle(name.Data(),title.Data());
10989 if (!f)
return hcd;
10992 Double_t xminold=f->GetXmin();
10993 Double_t xmaxold=f->GetXmax();
10995 f->SetRange(xmin,xmax);
10997 TH1* hf=(TH1*)f->GetHistogram();
11001 if (hcd.GetEntries()>0)
11004 if (mode.Contains(
"N")) title=
"Normalized ";
11005 if (mode.Contains(
"F")) title+=
"Forward ";
11006 if (mode.Contains(
"B")) title+=
"Backward ";
11007 title+=
"Cumulative Distribution Histogram of function ";
11010 TString namef=f->GetName();
11013 hcd.SetTitle(title.Data());
11016 f->SetRange(xminold,xmaxold);
11050 if (!dir || (frame!=
"equ" && frame!=
"gal" && frame!=
"ecl" && frame!=
"hor" && frame!=
"icr" && frame!=
"loc") ||
11051 (mode!=
"M" && mode!=
"T" && mode!=
"B" && mode!=
"J"))
11054 cout <<
" *" << ClassName() <<
"::InitDataNames* Invalid input encountered." << endl;
11055 cout <<
" dir=" << dir <<
" frame=" << frame <<
" mode=" << mode << endl;
11066 TString sdir=
"arrival";
11067 if (dir<0) sdir=
"moving";
11071 if (mode==
"M") frame=
"mean";
11072 if (mode==
"T") frame=
"true";
11073 if (mode==
"B") frame=
"B1950";
11074 if (mode==
"J") frame=
"J2000";
11075 frame+=
" equatorial";
11077 if (frame==
"gal") frame=
"galactic";
11078 if (frame==
"ecl") frame=
"ecliptic";
11079 if (frame==
"hor") frame=
"horizontal";
11080 if (frame==
"icr") frame=
"ICRS";
11081 if (frame==
"loc") frame=
"local";
11084 cout <<
" *" << ClassName() <<
"::InitDataNames* Prepared for input of " << sdir <<
" directions in " << frame <<
" coordinates." << endl;
11177 Bool_t error=kFALSE;
11179 if (obsname!=
"Name" && obsname!=
"Run" && obsname!=
"Event" && obsname!=
"Eventb" && obsname!=
"VetoLevel"
11180 && obsname!=
"DetId" && obsname!=
"Date" && obsname!=
"Tobs" && obsname!=
"Tstart" && obsname!=
"Tend"
11181 && obsname!=
"d" && obsname!=
"a" && obsname!=
"b" && obsname!=
"z"
11182 && obsname!=
"E" && obsname!=
"L" && obsname!=
"S" && obsname!=
"F" && obsname!=
"I" && obsname!=
"J"
11183 && obsname!=
"T90" && obsname!=
"T100"
11184 && obsname!=
"dsigma" && obsname!=
"csigma" && obsname!=
"zsigma"
11185 && obsname!=
"Esigma" && obsname!=
"Lsigma" && obsname!=
"Ssigma" && obsname!=
"Fsigma" && obsname!=
"Isigma"
11186 && obsname!=
"T90sigma" && obsname!=
"T100sigma") error=kTRUE;
11187 if (obsname==
"Date" && (units!=
"ddmmyyyy" && units!=
"yyyymmdd" && units!=
"mmddyyyy" && units!=
"yyyyddmm")) error=kTRUE;
11188 if ((obsname==
"Tobs" || obsname==
"Tstart" || obsname==
"Tend")
11189 && (units!=
"JD" && units!=
"MJD" && units!=
"TJD" && units!=
"hms" && units!=
"hrs")) error=kTRUE;
11190 if ((obsname==
"a" || obsname==
"b" || obsname==
"csigma")
11191 && (units!=
"rad" && units!=
"deg" && units!=
"dms" && units!=
"hms" && units!=
"hrs")) error=kTRUE;
11192 if (func!=
"none" && func!=
"Log" && func!=
"Ln") error=kTRUE;
11196 cout <<
" *" << ClassName() <<
"::SetDataNames* Invalid input encountered." << endl;
11197 cout <<
" obsname=" << obsname <<
" units=" << units <<
" func=" << func << endl;
11202 TObjString* obs=
new TObjString(obsname);
11203 TObjString* var=
new TObjString(varname);
11204 TObjString* u=
new TObjString(units);
11205 TObjString* f=
new TObjString(func);
11206 TObjString* val=
new TObjString(
"");
11224 TString sdir=
"undefined";
11228 TString frame=
"undefined";
11235 frame+=
" equatorial";
11244 cout <<
" *" << ClassName() <<
"::ListDataNames* Settings for input of " << sdir <<
" directions in " << frame <<
" coordinates." << endl;
11250 cout <<
" *** No settings were specified ***" << endl;
11254 cout <<
" *** The following " << n <<
" settings were specified ***" << endl;
11261 for (Int_t i=1; i<=n; i++)
11267 cout <<
" obsname=" << obs->GetString() <<
" varname=" << var->GetString() <<
" units=" << u->GetString() <<
" func=" << f->GetString() << endl;
11393 fBurstParameters->SetNameTitle(
"BurstParameters",
"Parameter settings for transient burst investigations");
11403 Double_t pi=acos(-1.);
11549 if (fTunits==0) Tfact=86400;
11550 if (fTunits==1) Tfact=3600;
11551 if (fTunits==3) Tfact=1e-9;
11552 if (fTunits==4) Tfact=1e-12;
11563 if (!fRecoangle) fAngresmin=fAngresfix;
11564 Float_t Minsigmatot=-1;
11565 if (fSumsigmas==-1) Minsigmatot=fAngresmin;
11566 if (fSumsigmas==0) Minsigmatot=fSigmamin;
11567 if (fSumsigmas==1) Minsigmatot=fSigmamin+fAngresmin;
11568 if (fSumsigmas==2) Minsigmatot= sqrt(fSigmamin*fSigmamin+fAngresmin*fAngresmin);
11569 name=
"Minsigmatot";
11574 if (!fRecoangle) fAngresmax=fAngresfix;
11575 Float_t Maxsigmatot=-1;
11576 if (fSumsigmas==-1) Maxsigmatot=fAngresmax;
11577 if (fSumsigmas==0) Maxsigmatot=fSigmamax;
11578 if (fSumsigmas==1) Maxsigmatot=fSigmamax+fAngresmax;
11579 if (fSumsigmas==2) Maxsigmatot= sqrt(fSigmamax*fSigmamax+fAngresmax*fAngresmax);
11580 name=
"Maxsigmatot";
11587 Float_t Dtwin=fTmax-fTmin;
11610 Float_t phimin=fRAmin;
11611 Float_t phimax=fRAmax;
11612 Float_t thmin=90.-fDeclmax;
11613 Float_t thmax=90.-fDeclmin;
11614 Float_t OmegaDecl=
GetSolidAngle(thmin,thmax,
"deg",phimin,phimax,
"deg");
11621 Float_t RbkgDecl=fBkgrate;
11624 RbkgDecl=fabs(fBkgrate)*OmegaDecl;
11631 Float_t NbkgHour=RbkgDecl*3600.;
11638 Float_t NbkgWin=RbkgDecl*fDtwin*Tfact;
11731 if (fTunits==1) tu=
"hours";
11732 if (fTunits==2) tu=
"sec";
11733 if (fTunits==3) tu=
"ns";
11734 if (fTunits==4) tu=
"ps";
11736 cout <<
" ========================= User provided source c.q. burst settings ===============================" << endl;
11739 printf(
" No limitation has been put on the number of sources to be accepted for analysis. \n");
11743 printf(
" Maximal number of sources to be accepted for analysis : %-i \n",fNmaxsrc);
11747 printf(
" No limitation has been put on the number of observed events to be accepted for analysis. \n");
11751 printf(
" Maximal number of observed events to be accepted for analysis : %-i \n",fNmaxevt);
11753 printf(
" Right ascension interval (J2000 in degrees) for source c.q. event position acceptance : [%-g,%-g] \n",fRAmin,fRAmax);
11754 printf(
" Declination interval (J2000 in degrees) for source c.q. event position acceptance : [%-g,%-g] \n",fDeclmin,fDeclmax);
11755 printf(
" Redshift interval for source acceptance : [%-g,%-g] \n",fabs(fZmin),fZmax);
11756 if (fZmin<0) printf(
" Random redshift values taken from z-distribution in case of unknown redshift \n");
11757 if (fAvgrbz>=0) printf(
" User defined average source redshift : %-g \n",fAvgrbz);
11758 printf(
" Event energy interval (in GeV) for event acceptance : [%-g,%-g] \n",fEmin,fEmax);
11759 printf(
" Position uncertainty interval (sigma in degrees) for source acceptance : [%-g,%-g] \n",fabs(fSigmamin),fSigmamax);
11760 if (fSigmamin<0) printf(
" Random sigma values taken from sigma-distribution when missing in loaded data \n");
11761 printf(
" Event angular resolution interval (sigma in degrees) for event acceptance : [%-g,%-g] \n",fAngresmin,fAngresmax);
11762 if (fDatype>0) printf(
" Sigma (combination) selection (-1=event sigma only 0=source sigma only 1=linear summation 2=quadratic summation) : %-i \n",fSumsigmas);
11765 if (fDatype<0) printf(
" Unrestricted angular search area. \n");
11766 if (!fDatype) printf(
" Fixed angular search circle around the source position : %-g degrees \n",fDawin);
11767 if (fDatype==1) printf(
" Fixed angular search circle (in combined max. source/event sigma) around the source position : %-g \n",fDawin);
11768 if (fDatype==2) printf(
" Variable angular search circle (in combined actual source/event sigma) around the source position : %-g \n",fDawin);
11772 if (fDatype<0) printf(
" Unrestricted angular search area. \n");
11773 if (!fDatype) printf(
" Fixed angular local zenith band above/below the source position : %-g degrees \n",fabs(fDawin));
11774 if (fDatype==1) printf(
" Fixed angular local zenith band (in combined max. source/event sigma) above/below the source position : %-g \n",fabs(fDawin));
11775 if (fDatype==2) printf(
" Variable angular local zenith band (in combined actual source/event sigma) above/below the source position : %-g \n",fabs(fDawin));
11777 printf(
" Number of requested background patches per source : %-i \n",fNbkg);
11780 printf(
" Duration interval (t90 in sec) for burst acceptance : [%-g,%-g] \n",fabs(fT90min),fT90max);
11781 if (fT90min<0) printf(
" Random values taken from T90-distribution in case T90 and T100 were missing \n");
11782 if (fAvgrbt90>0) printf(
" User defined average burst T90 duration : %-g sec. \n",fAvgrbt90);
11783 if (fTmax>fTmin) printf(
" Total search time window (in %-s) with the burst trigger at t=0 : [%-g,%-g] \n",tu.Data(),fTmin,fTmax);
11785 if (fTmax<=fTmin) printf(
" Search will be performed with no restrictions on the time window \n");
11786 if (fTbin<0) printf(
" Automatic time histogram binning with as mean number of bkg counts/bin : %-g \n",fabs(fTbin));
11787 if (!fTbin) printf(
" Variable time histogram binning with as size (in %-s) for the first time : %-g \n",tu.Data(),fVarTbin);
11792 printf(
" Time histogram bin size in average T90 units : %-g",fTbin);
11793 if (fAvgrbt90>0 || fNgrbs>0) printf(
" (=%-g sec)",fTbin*fabs(fAvgrbt90));
11798 printf(
" Time histogram bin size : %-g %-s \n",fTbin,tu.Data());
11803 printf(
" Automatic angular histogram binning with as mean number of bkg counts per bin : %-g \n",fabs(fAbin));
11807 printf(
" Angular histogram bin size : %-g degrees \n",fAbin);
11815 printf(
" Number of generated neutrinos per burst : %-g without statistical fluctuations \n",fabs(fGrbnu));
11819 printf(
" Maximum number of generated neutrinos per burst : %-g \n",fGrbnu);
11820 printf(
" The actual number of neutrinos may be less due to statistical fluctuations \n");
11822 printf(
" Energy interval (in GeV) for signal event generation at the source : [%-g,%-g] \n",fESigmin,fESigmax);
11825 printf(
" The signal neutrino energy at Earth will be corrected for redshift \n");
11829 printf(
" No redshift correction will be applied on the generated signal neutrino energy \n");
11833 printf(
" Neutrino production is assumed to be NOT coupled to the observed burst duration \n");
11834 printf(
" Mean decoupled time difference between burst gammas/GW and neutrinos : %-g sec. \n",fDtnu);
11838 printf(
" Neutrino production is assumed to be coupled to the observed burst duration \n");
11839 printf(
" Mean coupled time difference (in units of T90 w.r.t. trigger) between burst gammas/GW and neutrinos : %-g",fDtnu);
11843 printf(
" Sigma of mean time difference between burst gammas/GW and neutrinos : %-g sec. \n",fDtnus);
11847 printf(
" Sigma of mean time difference (in units of T90) between burst gammas/GW and neutrinos : %-g \n",fabs(fDtnus));
11849 printf(
" Default spectral index for a dN/dE=E^-alpha source induced signal energy spectrum within [%-g,%-g] GeV : %-g \n",fEmin,fEmax,fAlphasig);
11850 printf(
" --- The user may have changed the spectrum by invokation of the memberfunction MakeBurstEnergydist(). \n");
11851 printf(
" Default spectral index for a dN/dE=E^-alpha background energy spectrum within [%-g,%-g] GeV : %-g \n",fEmin,fEmax,fAlphabkg);
11852 printf(
" --- The user may have changed the spectrum by invokation of the memberfunction MakeBurstEnergydist(). \n");
11853 printf(
" Fixed event reco angular resolution (sigma in degrees), also used when no distribution (value) is available : %-g \n",fAngresfix);
11854 printf(
" Event reconstruction angular uncertainty selection (0=use fixed value 1=mean 2=median 3=draw from distribution) : %-i \n",fRecoangle);
11855 printf(
" Neutrino-lepton kinematic opening angle selection for CC interactions (0=none 1=mean 2=median 3=draw from pdf) : %-i \n",fKinangle);
11859 printf(
" Time resolution of the neutrino detector : %-g sec. \n",fTimres);
11860 printf(
" Area covered c.q. overlooked by the neutrino detector : %-g m^2 \n",fSensarea);
11861 if (fBkgrate) printf(
" User defined mean rate of background events for the specified declination interval (<0 : rate per steradian) : %-g Hz \n",fBkgrate);
11864 printf(
" ============================== Derived parameters ==================================== \n");
11865 printf(
" Combined source position and event reco angular uncertainty interval (sigma in degrees) : [%-g,%-g] \n",fMinsigmatot,fMaxsigmatot);
11866 printf(
" Solid angle coverage corresponding to the selected declination band : %-g steradian \n",fOmegaDecl);
11869 printf(
" Background event rate (from user setting) for the selected declination band : %-g Hz \n",fRbkgDecl);
11870 printf(
" Mean number of background events (from user setting) per hour from the selected declination band : %-g \n",fNbkgHour);
11871 printf(
" Mean number of background events (from user setting) in the time window from the selected declination band : %-g \n",fNbkgWin);
11875 printf(
" Number of bursts accepted for analysis : %-i \n",fNgrbs);
11876 printf(
" Median source redshift from the data sample : %-g \n",fabs(fAvgrbz));
11877 if (fAvgrbt90) printf(
" Median burst T90 duration from the data sample : %-g sec. \n",fabs(fAvgrbt90));
11878 printf(
" Median souce position uncertainty (sigma in degrees) from the data sample : %-g \n",fAvgrbsigma);
11880 if (fNevts>0) printf(
" Number of observed events accepted for analysis : %-i \n",fNevts);
11881 printf(
" ====================================================================================== \n");
11945 if (type==
"-" && src) type=
"GRB";
11946 if (type==
"-" && !src) type=
"EVT";
11952 cout <<
" *" << ClassName() <<
"::LoadInputData* No variables were specified for " << type <<
" data."<< endl;
11984 if (src && fZmin<0) zdist=
GetBurstZdist(
"LoadInputData()",type);
11988 if (src && fT90min<0) t90dist=
GetBurstT90dist(
"LoadInputData()",type);
11991 TH1* sigmaposdist=0;
11995 TChain data(tree.Data());
11996 data.Add(file.Data());
12005 TObjString* pval=0;
12016 TString Name,
Date,Tobs,Tstart,Tend;
12018 Float_t z,csigma,T90,T100,E;
12037 Int_t nent=data.GetEntries();
12043 fRefs=
new TObjArray(nent);
12048 size=
fRefs->GetSize();
12049 fRefs->Expand(size+nent);
12056 fSigs=
new TObjArray(nent);
12061 size=
fSigs->GetSize();
12062 fSigs->Expand(size+nent);
12067 for (Int_t ient=0; ient<nent; ient++)
12069 if (nmax>=0 && nnew>=nmax)
break;
12070 if (src && fNmaxsrc>=0 && (fNgrbs+nnew)>=fNmaxsrc)
break;
12071 if (!src && fNmaxevt>=0 && (fNevts+nnew)>=fNmaxevt)
break;
12073 data.GetEntry(ient);
12092 for (Int_t ivar=1; ivar<=nvars; ivar++)
12094 obsname=((TObjString*)
fDataNames.GetObject(ivar,1))->GetString();
12095 varname=((TObjString*)
fDataNames.GetObject(ivar,2))->GetString();
12096 units=((TObjString*)
fDataNames.GetObject(ivar,3))->GetString();
12097 func=((TObjString*)
fDataNames.GetObject(ivar,4))->GetString();
12099 pval=(TObjString*)
fDataNames.GetObject(ivar,5);
12101 if (!pval)
continue;
12103 lx=data.GetLeaf(varname);
12108 pval->SetString(
"-999");
12112 if (obsname==
"Name")
12116 Name=lxc->GetValueString();
12120 value=lx->GetValue();
12123 if (func==
"Log") value=pow(10,value);
12124 if (func==
"Ln") value=exp(value);
12127 if (obsname==
"a" || obsname==
"b" || obsname==
"csigma") value=
ConvertAngle(value,units,
"deg");
12130 if (units.IsFloat())
12139 pval->SetString(val);
12142 if (obsname==
"Date")
12146 if (units==
"ddmmyyyy") dmode=0;
12147 if (units==
"yyyymmdd") dmode=1;
12148 if (units==
"mmddyyyy") dmode=2;
12149 if (units==
"yyyyddmm") dmode=3;
12151 if (obsname==
"Tobs")
12154 if (units==
"JD") tobs.
SetJD(value);
12155 if (units==
"MJD") tobs.
SetMJD(value);
12156 if (units==
"TJD") tobs.
SetTJD(value);
12166 value=s+double(100*m+10000*h);
12170 if (obsname==
"Tstart")
12173 if (units==
"JD") tstart.
SetJD(value);
12174 if (units==
"MJD") tstart.
SetMJD(value);
12175 if (units==
"TJD") tstart.
SetTJD(value);
12185 value=s+double(100*m+10000*h);
12189 if (obsname==
"Tend")
12192 if (units==
"JD") tend.
SetJD(value);
12193 if (units==
"MJD") tend.
SetMJD(value);
12194 if (units==
"TJD") tend.
SetTJD(value);
12204 value=s+double(100*m+10000*h);
12209 if (obsname==
"d") d=value;
12210 if (obsname==
"a") a=value;
12211 if (obsname==
"b") b=value;
12212 if (obsname==
"z") z=value;
12213 if (obsname==
"csigma") csigma=value;
12214 if (obsname==
"T90") T90=value;
12215 if (obsname==
"T100") T100=value;
12216 if (obsname==
"E") E=value;
12223 if (a<-900 || b<-900)
continue;
12226 if (
fDataFrame==
"equ" &&
fDataMode==
"J" && (a<fRAmin || a>fRAmax || b<fDeclmin || b>fDeclmax))
continue;
12229 if (Tobs!=
"none" && Tobs!=
"set" &&
Date!=
"none") tobs.
SetUT(
Date,Tobs,dmode);
12230 if (Tstart!=
"none" && Tstart!=
"set" &&
Date!=
"none") tstart.
SetUT(
Date,Tstart,dmode);
12231 if (Tend!=
"none" && Tend!=
"set" &&
Date!=
"none") tend.
SetUT(
Date,Tend,dmode);
12234 if (Tobs==
"none" || (Tobs!=
"set" &&
Date==
"none"))
continue;
12237 tobs.GetDate(kTRUE,0,&yyyy,&mm,&dd);
12238 idate=dd+100*mm+10000*yyyy;
12246 jdate=idate%1000000;
12248 if (jdate<100000) grbname+=
"0";
12252 if (date1 && idate<date1)
continue;
12253 if (date2 && idate>date2)
continue;
12257 if (T90<=0) T90=T100;
12258 if (fT90min<0 && T90<0 && t90dist) T90=t90dist->GetRandom();
12260 if (T90<fabs(fT90min) || T90>fT90max)
continue;
12262 if (fZmin<0 && z<0 && zdist) z=zdist->GetRandom();
12264 if (z<fabs(fZmin) || z>fZmax)
continue;
12266 if (fSigmamin<0 && csigma<0 && sigmaposdist) csigma=sigmaposdist->GetRandom();
12268 if (csigma<fabs(fSigmamin) || csigma>fSigmamax)
continue;
12272 if (E<fEmin || E>fEmax)
continue;
12273 if (csigma<fAngresmin || csigma>fAngresmax)
continue;
12277 sx=
SetSignal(d,a,
"deg",b,
"deg",
fDataFrame,&tobs,-1,
fDataMode,grbname,iobs);
12283 GetSignal(d,a,
"deg",b,
"deg",
"equ",&tobs,jlast,
"J",iobs);
12286 if (a<fRAmin || a>fRAmax || b<fDeclmin || b>fDeclmax)
12299 for (Int_t ivar=1; ivar<=nvars; ivar++)
12301 obsname=((TObjString*)
fDataNames.GetObject(ivar,1))->GetString();
12302 val=((TObjString*)
fDataNames.GetObject(ivar,5))->GetString();
12305 if (obsname==
"Name")
continue;
12311 if (obsname==
"a" || obsname==
"b")
continue;
12314 if (obsname==
"Date") value=idate;
12315 if (obsname==
"Tobs") value=tobs.
GetMJD();
12316 if (obsname==
"Tstart") value=tstart.
GetMJD();
12317 if (obsname==
"TEnd") value=tend.
GetMJD();
12320 if (obsname==
"z") value=z;
12321 if (obsname==
"csigma") value=csigma;
12322 if (obsname==
"T90") value=T90;
12340 cout <<
" *" << ClassName() <<
"::LoadInputData* " << nnew <<
" new source(s) of type " << type
12341 <<
" stored from Tree:" << tree <<
" of file(s):" << file << endl;
12342 cout <<
" Total number of stored sources c.q. bursts : " << nstored << endl;
12349 cout <<
" *" << ClassName() <<
"::LoadInputData* " << nnew
12350 <<
" new observed event(s) were stored from Tree:" << tree <<
" of file(s):" << file << endl;
12351 cout <<
" Total number of stored events : " << nstored << endl;
12394 Double_t pi=acos(-1.);
12395 n=TMath::Nint(
double(n)*fOmegaDecl/(4.*pi));
12410 if (!zdist || !t90dist || !sigmaposdist)
12413 cout <<
" *" << ClassName() <<
"::GenBurstGCNdata* A distribution for random values is missing." << endl;
12418 Float_t thlow=90.-fDeclmax;
12419 Float_t thup=90.-fDeclmin;
12420 Float_t phimin=fRAmin;
12421 Float_t phimax=fRAmax;
12422 if (thlow<0) thlow=0;
12423 if (thup>180) thup=180;
12429 Float_t sigmagrb=0;
12431 Double_t thetagrb,phigrb,ragrb,decgrb;
12434 for (Int_t igrb=1; igrb<=n; igrb++)
12436 if (fNmaxsrc>=0 && (fNgrbs+ngen)>=fNmaxsrc)
break;
12439 if (fabs(fZmin)==fZmax) zgrb=fZmax;
12440 while (zgrb<fabs(fZmin) || zgrb>fZmax)
12442 zgrb=zdist->GetRandom();
12446 if (fabs(fT90min)==fT90max) t90grb=fT90max;
12447 while (t90grb<fabs(fT90min) || t90grb>fT90max)
12449 t90grb=t90dist->GetRandom();
12450 t90grb=pow(
float(10),t90grb);
12454 if (fabs(fSigmamin)==fSigmamax) sigmagrb=fSigmamax;
12455 while (sigmagrb<fabs(fSigmamin) || sigmagrb>fSigmamax)
12457 sigmagrb=sigmaposdist->GetRandom();
12462 thetagrb=rgrb.
GetX(2,
"sph",
"deg");
12463 phigrb=rgrb.
GetX(3,
"sph",
"deg");
12470 decgrb=90.-thetagrb;
12471 sx=
SetSignal(1,ragrb,
"deg",decgrb,
"deg",
"equ",0,-1,
"J",grbname);
12492 cout <<
" *" << ClassName() <<
"::GenBurstGCNdata* " << ngen <<
" new generated bursts with name " << name <<
" were stored." << endl;
12493 cout <<
" Total number of stored bursts : " << fNgrbs << endl;
12531 TChain data(tree.Data());
12532 data.Add(file.Data());
12534 Int_t nen=data.GetEntries();
12535 TLeaf* lx=data.FindLeaf(name.Data());
12539 cout <<
" *" << ClassName() <<
"::MakeBurstZdist* Missing information for tree variable:" << name << endl;
12540 cout <<
" of Tree:" << tree <<
" with " << nen <<
" entries in file:" << file << endl;
12549 TH1F* hz=
new TH1F(
"hz",
"Archival data of observed burst redshifts",nb,zmin,zmax);
12551 hz->GetXaxis()->SetTitle(
"Burst redshift");
12552 hz->GetYaxis()->SetTitle(
"Counts");
12557 TH1F* hd=
new TH1F(
"hd",
"Burst distances derived from the archival redshift data",nb,dmin,dmax);
12559 hd->GetXaxis()->SetTitle(
"Burst physical distance in Mpc");
12560 hd->GetYaxis()->SetTitle(
"Counts");
12570 for (Int_t ien=0; ien<nen; ien++)
12572 data.GetEntry(ien);
12574 lx=data.GetLeaf(name.Data());
12578 if (z<zmin || z>zmax)
continue;
12587 cout <<
" *" << ClassName() <<
"::MakeBurstZdist* " << nz <<
" archival z-values have been obtained from tree variable:" << name
12588 <<
" of Tree:" << tree <<
" in file(s):" << file << endl;
12605 if (type.Contains(
"GRB"))
12608 cout <<
" *" << ClassName() <<
"::GetBurstZdist* Called from " << name << endl;
12609 cout <<
" *** Archival observed redshift distribution not found. ***" << endl;
12610 cout <<
" A Landau fit from Swift GRB redshift data will be used to provide missing c.q. random z values." << endl;
12615 TF1 fz(
"fz",
"59.54*TMath::Landau(x,1.092,0.5203)");
12618 TH1* hfz=fz.GetHistogram();
12619 zdist=(TH1*)hfz->Clone();
12620 zdist->SetNameTitle(
"hZpdf",
"Landau fit for Swift GRB z data");
12621 zdist->GetXaxis()->SetTitle(
"GRB redshift");
12622 zdist->GetYaxis()->SetTitle(
"Counts");
12629 cout <<
" *" << ClassName() <<
"::GetBurstZdist* Called from " << name << endl;
12630 cout <<
" *** No redshift fit is available for source class " << type <<
" ***" << endl;
12670 TChain data(tree.Data());
12671 data.Add(file.Data());
12673 Int_t nen=data.GetEntries();
12674 TLeaf* lx=data.FindLeaf(name.Data());
12678 cout <<
" *" << ClassName() <<
"::MakeBurstT90dist* Missing information for tree variable:" << name << endl;
12679 cout <<
" of Tree:" << tree <<
" with " << nen <<
" entries in file:" << file << endl;
12688 TH1F* ht90=
new TH1F(
"ht90",
"Archival data of observed burst durations",nb,xmin,xmax);
12690 ht90->GetXaxis()->SetTitle(
"Burst duration ^{10}log(T90) in sec.");
12691 ht90->GetYaxis()->SetTitle(
"Counts");
12699 for (Int_t ien=0; ien<nen; ien++)
12701 data.GetEntry(ien);
12703 lx=data.GetLeaf(name.Data());
12706 t90=lx->GetValue();
12709 ht90->Fill(log10(t90));
12714 cout <<
" *" << ClassName() <<
"::MakeBurstT90dist* " << nt90 <<
" archival T90 values have been obtained from variable:" << name
12715 <<
" of Tree:" << tree <<
" in file(s):" << file << endl;
12732 if (type.Contains(
"GRB"))
12735 cout <<
" *" << ClassName() <<
"::GetBurstT90dist* Called from " << name << endl;
12736 cout <<
" *** Observational T90 distribution not found. ***" << endl;
12737 cout <<
" A double Gaussian fit from Fermi GRB T90 data will be used to provide missing c.q. random T90 values." << endl;
12742 TF1 ft(
"ft",
"44.39*TMath::Gaus(x,-0.131,0.481)+193.8*TMath::Gaus(x,1.447,0.4752)");
12745 TH1* hft=ft.GetHistogram();
12746 t90dist=(TH1*)hft->Clone();
12747 t90dist->SetNameTitle(
"hT90pdf",
"Double Gauss fit for Fermi t90 data");
12748 t90dist->GetXaxis()->SetTitle(
"GRB duration ^{10}log(T90) in sec.");
12749 t90dist->GetYaxis()->SetTitle(
"Counts");
12756 cout <<
" *" << ClassName() <<
"::GetBurstT90dist* Called from " << name << endl;
12757 cout <<
" *** No T90 fit is available for source class " << type <<
" ***" << endl;
12806 TChain data(tree.Data());
12807 data.Add(file.Data());
12809 Int_t nen=data.GetEntries();
12810 TLeaf* lx=data.FindLeaf(name.Data());
12814 cout <<
" *" << ClassName() <<
"::MakeBurstSigmaPosdist* Missing information for tree variable:" << name << endl;
12815 cout <<
" of Tree:" << tree <<
" with " << nen <<
" entries in file:" << file << endl;
12820 TH1* sigmaposdist=(TH1*)
fBurstHistos.FindObject(
"hsigmapos");
12824 TH1F* hsigmapos=
new TH1F(
"hsigmapos",
"Archival data of observed 1-sigma burst position uncertainties",nb,xmin,xmax);
12826 hsigmapos->GetXaxis()->SetTitle(
"Burst position uncertainty (sigma in degrees)");
12827 hsigmapos->GetYaxis()->SetTitle(
"Counts");
12831 TH1* hsigmapos=(TH1*)
fBurstHistos.FindObject(
"hsigmapos");
12834 Double_t sigmapos=0;
12835 for (Int_t ien=0; ien<nen; ien++)
12837 data.GetEntry(ien);
12839 lx=data.GetLeaf(name.Data());
12842 sigmapos=lx->GetValue();
12847 if (sigmapos<xmin || sigmapos>xmax)
continue;
12849 hsigmapos->Fill(sigmapos);
12853 cout <<
" *" << ClassName() <<
"::MakeBurstSigmaPosdist* " << nsigmapos <<
" archival sigmapos values have been obtained from variable:" << name
12854 <<
" of Tree:" << tree <<
" in file(s):" << file << endl;
12868 TH1* sigmaposdist=(TH1*)
fBurstHistos.FindObject(
"hsigmapos");
12871 if (type.Contains(
"GRB"))
12874 cout <<
" *" << ClassName() <<
"::GetBurstSigmaPosdist* Called from " << name << endl;
12875 cout <<
" *** Archival observed GRB position uncertainty distribution not found. ***" << endl;
12876 cout <<
" A Landau fit from observed GRB data will be used to provide missing c.q. random 1-sigma uncertainty values." << endl;
12878 sigmaposdist=(TH1*)
fBurstHistos.FindObject(
"hSigmaSourcePDF");
12881 TF1 fsigmapos(
"fsigmapos",
"245.2*TMath::Landau(x,-2.209,0.6721,1)");
12882 fsigmapos.SetRange(0,90);
12883 fsigmapos.SetNpx(10000);
12884 TH1* hfsigmapos=fsigmapos.GetHistogram();
12885 sigmaposdist=(TH1*)hfsigmapos->Clone();
12886 sigmaposdist->SetNameTitle(
"hSigmaSourcePDF",
"Landau fit for burst 1-sigma position uncertainty data");
12887 sigmaposdist->GetXaxis()->SetTitle(
"Burst position uncertainty (sigma in degrees)");
12888 sigmaposdist->GetYaxis()->SetTitle(
"Counts");
12895 cout <<
" *" << ClassName() <<
"::GetBurstSigmaPosdist* Called from " << name << endl;
12896 cout <<
" *** No position uncertainty fit is available for source class " << type <<
" ***" << endl;
12900 return sigmaposdist;
12931 if (Emin<=0) Emin=1e-10;
12933 if ((mode!=1 && mode!=2) || Emax<=Emin)
12935 cout <<
" *" << ClassName() <<
"::MakeBurstEnergydist* Inconsistent data: mode=" << mode <<
" Emin=" << Emin <<
" Emax=" << Emax << endl;
12940 Double_t xmin=log10(Emin);
12941 Double_t xmax=log10(Emax);
12943 TString sf=
"dN/dE=";
12944 sf+=spec.GetExpFormula(
"p");
12945 if (mode==1) sf+=
" background";
12946 if (mode==2) sf+=
" signal";
12947 sf+=
" distribution";
12948 sf.ReplaceAll(
"x",
"E");
12950 if (mode==1) sh=
"Background energy distribution;^{10}Log(E) in GeV;";
12951 if (mode==2) sh=
"Source induced signal energy distribution;^{10}Log(E) in GeV;";
12957 if (mode==1) dist=(TH1*)
fBurstHistos.FindObject(
"hBkgEpdf");
12958 if (mode==2) dist=(TH1*)
fBurstHistos.FindObject(
"hSigEpdf");
12971 TH1F* hBkgEpdf=(TH1F*)his.Clone();
12972 hBkgEpdf->SetName(
"hBkgEpdf");
12979 TH1F* hSigEpdf=(TH1F*)his.Clone();
12980 hSigEpdf->SetName(
"hSigEpdf");
12986 cout <<
" *" << ClassName() <<
"::MakeBurstEnergydist* Created a " << sf
12987 <<
" on [Emin,Emax]=[" << Emin <<
"," << Emax <<
"] GeV" << endl;
13019 TF1 spec(
"spec",
"pow(x,[0])");
13020 spec.SetParameter(0,-alpha);
13084 if (Emin<=0) Emin=1e-10;
13086 if ((abs(mode)!=1 && abs(mode)!=2) || Emax<=Emin)
13088 cout <<
" *" << ClassName() <<
"::MakeBurstEnergydist* Inconsistent data: mode=" << mode <<
" Emin=" << Emin <<
" Emax=" << Emax << endl;
13093 Double_t xmin=log10(Emin);
13094 Double_t xmax=log10(Emax);
13097 TChain data(tree.Data());
13098 data.Add(file.Data());
13100 Int_t nen=data.GetEntries();
13102 if (!nen || !data.FindLeaf(name1.Data()) || !data.FindLeaf(name2.Data()))
13104 cout <<
" *" << ClassName() <<
"::MakeBurstEnergydist* Missing information for tree variable:" << name1
13105 <<
" and/or tree variable:" << name2 << endl;
13106 cout <<
" of Tree:" << tree <<
" with " << nen <<
" entries in file:" << file << endl;
13112 if (abs(mode)==1) flag=TMath::Nint(
fBurstParameters->GetSignal(
"PDFbkgE"));
13113 if (abs(mode)==2) flag=TMath::Nint(
fBurstParameters->GetSignal(
"PDFsigE"));
13115 if (flag) mode=abs(mode);
13119 if (abs(mode)==1) Edist=(TH1*)
fBurstHistos.FindObject(
"hBkgEpdf");
13120 if (abs(mode)==2) Edist=(TH1*)
fBurstHistos.FindObject(
"hSigEpdf");
13122 if (mode>0 && Edist)
13137 TH1F* hBkgEpdf=
new TH1F(
"hBkgEpdf",
"Archival data of observed background energies",nb,xmin,xmax);
13139 hBkgEpdf->GetXaxis()->SetTitle(
"^{10}log(Energy) in GeV");
13140 hBkgEpdf->GetYaxis()->SetTitle(
"Counts");
13147 TH1F* hSigEpdf=
new TH1F(
"hSigEpdf",
"Archival data of observed signal energies",nb,xmin,xmax);
13149 hSigEpdf->GetXaxis()->SetTitle(
"^{10}log(Energy) in GeV");
13150 hSigEpdf->GetYaxis()->SetTitle(
"Counts");
13158 if (abs(mode)==1) Edist=(TH1*)
fBurstHistos.FindObject(
"hBkgEpdf");
13159 if (abs(mode)==2) Edist=(TH1*)
fBurstHistos.FindObject(
"hSigEpdf");
13165 for (Int_t ien=0; ien<nen; ien++)
13167 data.GetEntry(ien);
13169 lx=data.GetLeaf(name1.Data());
13172 logE=lx->GetValue();
13174 lx=data.GetLeaf(name2.Data());
13177 dec=lx->GetValue();
13182 if (dec>=fDeclmin && dec<=fDeclmax)
13190 if (abs(mode)==1) smode=
"archival background";
13191 if (abs(mode)==2) smode=
"archival signal";
13195 cout <<
" *" << ClassName() <<
"::MakeBurstEnergydist* A new " << smode
13196 <<
" energy distribution has been created." << endl;
13200 cout <<
" *" << ClassName() <<
"::MakeBurstEnergydist* Statistics of the existing " << smode
13201 <<
" energy distribution have been increased." << endl;
13204 cout <<
" " << nE <<
" energy values have been obtained from variable:" << name1
13205 <<
" of Tree:" << tree <<
" in file(s):" << file << endl;
13208void NcAstrolab::MakeBurstRecoAngresdist(TString file,TString tree,TString name1,TString name2,TString ua,TString name3,TString ud,Double_t Emin,Double_t Emax,Int_t nbe,Int_t nba)
13261 if (Emin<=0) Emin=1e-10;
13265 cout <<
" *" << ClassName() <<
"::MakeBurstRecoAngresdist* Inconsistent data: Emin=" << Emin <<
" Emax=" << Emax << endl;
13270 Double_t xmin=log10(Emin);
13271 Double_t xmax=log10(Emax);
13274 TChain data(tree.Data());
13275 data.Add(file.Data());
13277 Int_t nen=data.GetEntries();
13279 if (!nen || !data.FindLeaf(name1.Data()) || !data.FindLeaf(name2.Data()) || !data.FindLeaf(name3.Data()))
13281 cout <<
" *" << ClassName() <<
"::MakeBurstRecoAngresdist* Missing information for tree variable:" << name1
13282 <<
" and/or tree variable:" << name2 <<
" and/or tree variable:" << name3 << endl;
13283 cout <<
" of Tree:" << tree <<
" with " << nen <<
" entries in file:" << file << endl;
13288 TH2* Angresdist=(TH2*)
fBurstHistos.FindObject(
"hAngresE");
13292 TH2F* hAngresE=
new TH2F(
"hAngresE",
"Archival data of observed reconstruction angle resolution vs. energy",
13293 nbe,xmin,xmax,nba,0,180.1);
13295 hAngresE->GetXaxis()->SetTitle(
"^{10}log(Energy) in GeV");
13296 hAngresE->GetYaxis()->SetTitle(
"Angular resolution in degrees");
13300 TH2* hAngresE=(TH2*)
fBurstHistos.FindObject(
"hAngresE");
13307 for (Int_t ien=0; ien<nen; ien++)
13309 data.GetEntry(ien);
13311 lx=data.GetLeaf(name1.Data());
13314 logE=lx->GetValue();
13316 lx=data.GetLeaf(name2.Data());
13319 dang=lx->GetValue();
13324 lx=data.GetLeaf(name3.Data());
13327 dec=lx->GetValue();
13332 if (dec>=fDeclmin && dec<=fDeclmax)
13334 hAngresE->Fill(logE,dang);
13339 cout <<
" *" << ClassName() <<
"::MakeBurstRecoAngresdist* " << nE <<
" archival entries have been obtained for variables:" << name2
13340 <<
" vs. " << name1 <<
" of Tree:" << tree <<
" in file(s):" << file << endl;
13366 TH1* hSigEpdf=(TH1*)
fBurstHistos.FindObject(
"hSigEpdf");
13368 if (!hSigEpdf)
return E;
13370 Int_t nbins=hSigEpdf->GetNbinsX();
13371 Int_t nentries=hSigEpdf->GetEntries();
13373 if (nbins<=0 || nentries<=0)
return E;
13375 TAxis* xaxis=hSigEpdf->GetXaxis();
13377 if (!xaxis)
return E;
13379 Double_t xlow=xaxis->GetBinLowEdge(1);
13380 Double_t xup=xaxis->GetBinUpEdge(nbins);
13382 Double_t logEmin=0;
13389 logEmin=log10(Emin);
13392 Double_t logEmax=0;
13399 logEmax=log10(Emax);
13402 if (logEmax<=logEmin || logEmin>=xup || logEmax<=xlow)
return E;
13404 while (E<logEmin || E>logEmax)
13406 E=hSigEpdf->GetRandom();
13409 E=pow(
float(10),E);
13437 TH1* hBkgEpdf=(TH1*)
fBurstHistos.FindObject(
"hBkgEpdf");
13439 if (!hBkgEpdf)
return E;
13441 Int_t nbins=hBkgEpdf->GetNbinsX();
13442 Int_t nentries=hBkgEpdf->GetEntries();
13444 if (nbins<=0 || nentries<=0)
return E;
13446 TAxis* xaxis=hBkgEpdf->GetXaxis();
13448 if (!xaxis)
return E;
13450 Double_t xlow=xaxis->GetBinLowEdge(1);
13451 Double_t xup=xaxis->GetBinUpEdge(nbins);
13453 Double_t logEmin=0;
13460 logEmin=log10(Emin);
13463 Double_t logEmax=0;
13470 logEmax=log10(Emax);
13473 if (logEmax<=logEmin || logEmin>=xup || logEmax<=xlow)
return E;
13475 while (E<logEmin || E>logEmax)
13477 E=hBkgEpdf->GetRandom();
13480 E=pow(
float(10),E);
13513 if (Amin<0) Amin=0;
13528 TH2* hAngresE=(TH2*)
fBurstHistos.FindObject(
"hAngresE");
13539 Int_t nbins=hAngresE->GetNbinsX();
13540 Int_t nentries=hAngresE->GetEntries();
13542 if (nbins<=0 || nentries<=0)
return dang;
13544 TAxis* xaxis=hAngresE->GetXaxis();
13546 if (!xaxis)
return dang;
13548 Double_t xlow=xaxis->GetBinLowEdge(1);
13549 Double_t xup=xaxis->GetBinUpEdge(nbins);
13551 Double_t logEmin=0;
13558 logEmin=log10(Emin);
13561 Double_t logEmax=0;
13568 logEmax=log10(Emax);
13571 if (logEmax<logEmin || logEmin>=xup || logEmax<=xlow)
return dang;
13573 Int_t ilow=xaxis->FindBin(logEmin);
13574 Int_t iup=xaxis->FindBin(logEmax);
13576 TH1D* hproj=hAngresE->ProjectionY(
"hproj",ilow,iup);
13578 if (!hproj)
return dang;
13580 nbins=hproj->GetNbinsX();
13581 nentries=hproj->GetEntries();
13583 if (nbins<=0 || nentries<=0)
return dang;
13585 if (fRecoangle==1) dang=hproj->GetMean();
13595 xaxis=hproj->GetXaxis();
13597 if (!xaxis)
return dang;
13599 xlow=xaxis->GetBinLowEdge(1);
13600 xup=xaxis->GetBinUpEdge(nbins);
13602 if (Amax<=Amin || Amin>=xup || Amax<=xlow)
return dang;
13605 while (dang<Amin || dang>Amax)
13607 dang=hproj->GetRandom();
13611 if (hproj)
delete hproj;
13642 printf(
"\n *%-s::GenBurstSignals* Error : [Tmin,Tmax]=[%-g,%-g] whereas Tmin<Tmax is required. \n",ClassName(),fTmin,fTmax);
13701 TH1* hSigmaReco=(TH1*)
fBurstHistos.FindObject(
"hSigmaReco");
13703 TH1* hSigEzcor=(TH1*)
fBurstHistos.FindObject(
"hSigEzcor");
13714 Float_t sigmagrb=0;
13717 Double_t thetagrb,phigrb;
13718 Double_t dmu,thetamu,phimu;
13724 Float_t dangmaxon=0;
13725 Float_t dangmaxoff=0;
13726 Float_t thlow,thup;
13727 Float_t ranlow,ranup;
13731 Float_t solidangle=0;
13732 Float_t ramu,decmu;
13735 Double_t sigmareco=0;
13736 Float_t sigmatot=0;
13738 Int_t fixedwinset=0;
13746 for (Int_t igrb=0; igrb<fNgrbs; igrb++)
13756 GetSignal(dgrb,thetagrb,
"deg",phigrb,
"deg",
"loc",tx,igrb+1);
13760 if (!fDatype) dangmax=fabs(fDawin);
13764 if (fMaxsigmatot>0) dangmax=fabs(fDawin*fMaxsigmatot);
13769 if (!fSumsigmas) dangmax=fabs(fDawin*sigmagrb);
13773 if (fSumsigmas==-1) sigmatot=fAngresfix;
13774 if (fSumsigmas==1) sigmatot=sigmagrb+fAngresfix;
13775 if (fSumsigmas==2) sigmatot=sqrt(sigmagrb*sigmagrb+fAngresfix*fAngresfix);
13776 if (sigmatot>=0) dangmax=fabs(fDawin*sigmatot);
13781 if (dangmax<0) fixedwinset=0;
13784 if (fDatype==1 && fMaxsigmatot<0) dangmax=-1;
13787 name=
"fixedwinset";
13804 if (fixedwinset && dangmax>=0)
13808 thlow=thetagrb-0.5*dangmax;
13809 thup=thetagrb+0.5*dangmax;
13821 sx->
SetSignal(solidangle*
float(fNbkg),
"OmegaOff");
13826 if (fSumsigmas==0 || fSumsigmas==1 || fSumsigmas==2) dangmaxon=fabs(fDawin*sigmagrb);
13827 dangmaxoff=dangmaxon;
13832 for (Int_t bkgpatch=0; bkgpatch<=fNbkg; bkgpatch++)
13834 nmu=int(
fRan->Poisson(fNbkgWin));
13835 for (Int_t imu=0; imu<nmu; imu++)
13840 dt=
fRan->Uniform(ranlow,ranup);
13845 thlow=90.-fDeclmax;
13848 decmu=90.-rmu.
GetX(2,
"sph",
"deg");
13849 ramu=rmu.
GetX(3,
"sph",
"deg");
13852 SetSignal(1,ramu,
"deg",decmu,
"deg",
"equ",&tmu,fNgrbs+1,
"J",
"bkgtemp",0);
13853 GetSignal(dmu,thetamu,
"deg",phimu,
"deg",
"loc",&tmu,fNgrbs+1);
13858 dang=fabs(thetagrb-thetamu);
13866 if (fDatype>=0 && fixedwinset && dang>dangmax)
continue;
13871 if (E<0 || E<fEmin || E>fEmax)
continue;
13875 if (sigmareco<0) sigmareco=fAngresfix;
13877 if (sigmareco<fAngresmin || sigmareco>fAngresmax)
continue;
13879 if (hSigmaReco) hSigmaReco->Fill(sigmareco);
13882 if (fSumsigmas==-1) sigmatot=sigmareco;
13883 if (fSumsigmas==0) sigmatot=sigmagrb;
13884 if (fSumsigmas==1) sigmatot=sigmagrb+sigmareco;
13885 if (fSumsigmas==2) sigmatot=sqrt(sigmagrb*sigmagrb+sigmareco*sigmareco);
13891 if (sigmatot>=0) dangmax=fabs(fDawin*sigmatot);
13894 if (fDatype>=0 && dang>dangmax)
continue;
13898 if (dangmax>dangmaxon) dangmaxon=dangmax;
13905 if (dangmax>dangmaxoff) dangmaxoff=dangmax;
13921 if (fGrbnu>=0 || nmugrb<
int(fabs(fGrbnu)*
float(fNgrbs)))
13927 nmu=int(fabs(fGrbnu));
13928 if (!nmu &&
fRan->Uniform()<fabs(fGrbnu)) nmu=1;
13929 for (Int_t imu=0; imu<nmu; imu++)
13936 dt=
fRan->Gauss(fDtnu,fabs(fDtnus)*t90grb);
13940 dt=
fRan->Gauss(fDtnu,fDtnus);
13948 dt=
fRan->Gauss(fDtnu*t90grb,fabs(fDtnus)*t90grb);
13952 dt=
fRan->Gauss(fDtnu*t90grb,fDtnus);
13955 if (fTimres>0) dt=
fRan->Gauss(dt,fTimres);
13966 if (hSigE) hSigE->Fill(E);
13969 if (fEzcor) E=E/(zgrb+1.);
13971 if (hSigEzcor) hSigEzcor->Fill(E);
13973 if (E<0 || E<fEmin || E>fEmax)
continue;
13978 Int_t mode=fKinangle-1;
13985 if (sigmareco<0) sigmareco=fAngresfix;
13987 if (sigmareco<fAngresmin || sigmareco>fAngresmax)
continue;
13989 if (hSigmaReco) hSigmaReco->Fill(sigmareco);
13997 if (fSumsigmas==-1) sigmatot=sigmareco;
13998 if (fSumsigmas==0) sigmatot=sigmagrb;
13999 if (fSumsigmas==1) sigmatot=sigmagrb+sigmareco;
14000 if (fSumsigmas==2) sigmatot=sqrt(sigmagrb*sigmagrb+sigmareco*sigmareco);
14006 if (sigmatot>=0) dangmax=fabs(fDawin*sigmatot);
14009 if (fDatype>=0 && dang>dangmax)
continue;
14011 if (dangmax>dangmaxon) dangmaxon=dangmax;
14021 if (fixedwinset)
continue;
14026 thlow=thetagrb-0.5*dangmaxon;
14027 thup=thetagrb+0.5*dangmaxon;
14031 thlow=thetagrb-0.5*dangmaxoff;
14032 thup=thetagrb+0.5*dangmaxoff;
14034 sx->
SetSignal(dangmaxoff,
"dangmaxOff");
14035 sx->
SetSignal(solidangle*
float(fNbkg),
"OmegaOff");
14046 sx->
SetSignal(dangmaxoff,
"dangmaxOff");
14047 sx->
SetSignal(solidangle*
float(fNbkg),
"OmegaOff");
14058 name=
"SolidangleOn";
14061 name=
"SolidangleOff";
14064 for (Int_t igrb=1; igrb<=fNgrbs; igrb++)
14097 Double_t pi=acos(-1.);
14118 Float_t fEnergyBkg=fEnergyOn-fEnergySig;
14122 if (fTunits==1) tu=
"hours";
14123 if (fTunits==2) tu=
"sec";
14124 if (fTunits==3) tu=
"ns";
14125 if (fTunits==4) tu=
"ps";
14130 Double_t TminOn=fTmin;
14131 Double_t TmaxOn=fTmax;
14132 Double_t TminOff=fTmin;
14133 Double_t TmaxOff=fTmax;
14149 Float_t DtwinOn=TmaxOn-TminOn;
14150 Float_t DtwinOff=TmaxOff-TminOff;
14151 Float_t NbkgWinOn=fRbkgDecl*DtwinOn*fTfact;
14152 Float_t NbkgWinOff=fRbkgDecl*DtwinOff*fTfact;
14161 Double_t binsize=0;
14162 Double_t binsizecos=0;
14164 Double_t abinsizeOn=0;
14166 Double_t abinsizeOff=0;
14169 TH1F* hOnSourceZ=0;
14170 TH1F* hOffSourceZ=0;
14171 TH1F* hOnSigSourceZ=0;
14172 TH1F* hOnSigmaSource=0;
14173 TH1F* hOffSigmaSource=0;
14174 TH1F* hOnSigSigmaSource=0;
14175 TH1F* hOnSigmaReco=0;
14176 TH1F* hOffSigmaReco=0;
14177 TH1F* hOnSigSigmaReco=0;
14178 TH1F* hOnSigmaComb=0;
14179 TH1F* hOffSigmaComb=0;
14180 TH1F* hOnSigSigmaComb=0;
14184 title.Form(
"On-Source object redshifts with matching event(s);Redshift;Counts");
14185 hOnSourceZ=
new TH1F(
"hOnSourceZ",title,nbins,1,0);
14186 hOnSourceZ->SetBuffer(non);
14188 title.Form(
"On-Source object position uncertainties with matching event(s);Object position angular uncertainty (sigma in degrees);Counts");
14189 hOnSigmaSource=
new TH1F(
"hOnSigmaSource",title,nbins,1,0);
14190 hOnSigmaSource->SetBuffer(non);
14192 title.Form(
"On-source event reconstruction uncertainties in the final sample;Event angular reconstruction uncertainty (sigma in degrees);Counts");
14193 hOnSigmaReco=
new TH1F(
"hOnSigmaReco",title,nbins,1,0);
14194 hOnSigmaReco->SetBuffer(non);
14196 title.Form(
"On-source combined object position and event reconstruction uncertainty;Combined object position and event reco angular uncertainty (sigma in degrees);Counts");
14197 hOnSigmaComb=
new TH1F(
"hOnSigmaComb",title,nbins,1,0);
14198 hOnSigmaComb->SetBuffer(non);
14203 title.Form(
"Off-Source object redshifts with matching event(s);Redshift;Counts");
14204 hOffSourceZ=
new TH1F(
"hOffSourceZ",title,nbins,1,0);
14205 hOffSourceZ->SetBuffer(noff);
14207 title.Form(
"Off-Source object position uncertainties with matching event(s);Object position angular uncertainty (sigma in degrees);Counts");
14208 hOffSigmaSource=
new TH1F(
"hOffSigmaSource",title,nbins,1,0);
14209 hOffSigmaSource->SetBuffer(noff);
14211 title.Form(
"Off-source event reconstruction uncertainties in the final sample;Event angular reconstruction uncertainty (sigma in degrees);Counts");
14212 hOffSigmaReco=
new TH1F(
"hOffSigmaReco",title,nbins,1,0);
14213 hOffSigmaReco->SetBuffer(noff);
14215 title.Form(
"Off-source combined object position and event reconstruction uncertainty;Combined object position and event reco angular uncertainty (sigma in degrees);Counts");
14216 hOffSigmaComb=
new TH1F(
"hOffSigmaComb",title,nbins,1,0);
14217 hOffSigmaComb->SetBuffer(noff);
14222 title.Form(
"On-Source object redshifts with matching simulated signal event(s);Redshift;Counts");
14223 hOnSigSourceZ=
new TH1F(
"hOnSigSourceZ",title,nbins,1,0);
14224 hOnSigSourceZ->SetBuffer(nsig);
14226 title.Form(
"On-Source object position uncertainties with matching simulated signal event(s);Object position angular uncertainty (sigma in degrees);Counts");
14227 hOnSigSigmaSource=
new TH1F(
"hOnSigSigmaSource",title,nbins,1,0);
14228 hOnSigSigmaSource->SetBuffer(nsig);
14230 title.Form(
"On-source simulated signal event reconstruction uncertainties in the final sample;Event angular reconstruction uncertainty (sigma in degrees);Counts");
14231 hOnSigSigmaReco=
new TH1F(
"hOnSigSigmaReco",title,nbins,1,0);
14232 hOnSigSigmaReco->SetBuffer(nsig);
14234 title.Form(
"On-source combined object position and simulated signal event reconstruction uncertainty;Combined object position and event reco angular uncertainty (sigma in degrees);Counts");
14235 hOnSigSigmaComb=
new TH1F(
"hOnSigSigmaComb",title,nbins,1,0);
14236 hOnSigSigmaComb->SetBuffer(nsig);
14246 title.Form(
"On-source reconstructed event energy in the final sample;Event energy in GeV;Counts");
14247 if (!mode) title.ReplaceAll(
"reconstructed",
"simulated");
14248 hOnE=
new TH1F(
"hOnE",title,nbins,1,0);
14252 title.Form(
"Off-source reconstructed event energy in the final sample;Event energy in GeV;Counts");
14253 if (!mode) title.ReplaceAll(
"reconstructed",
"simulated");
14254 hOffE=
new TH1F(
"hOffE",title,nbins,1,0);
14258 title.Form(
"On-source simulated signal event energy in the final sample;Event energy in GeV;Counts");
14259 hOnSigE=
new TH1F(
"hOnSigE",title,nbins,1,0);
14265 if (
fTscmode>0) scrt=
"(scrambled)";
14274 TH1F* hOnSigcosa=0;
14280 nbins=int((AmaxOn-AminOn)/binsize);
14284 nbins=int(((AmaxOn-AminOn)/180.)*NbkgWinOn*
float(fNgrbs)/fabs(fAbin));
14285 if (nbins) binsize=(AmaxOn-AminOn)/
float(nbins);
14289 title.Form(
"Reconstructed %-s opening angle of on-source events in time window;Opening angle (degrees);Counts per %-.3g degrees",scrp.Data(),binsize);
14290 hOna=
new TH1F(
"hOna",title,nbins+1,AminOn,AmaxOn+binsize);
14295 binsize=AmaxOn/float(nbins);
14296 title.Form(
"Reconstructed %-s opening angle of on-source events in time window;Opening angle (degrees);Counts per %-.3g degrees",scrp.Data(),binsize);
14297 hOna=
new TH1F(
"hOna",title,nbins+2,AminOn-binsize,AmaxOn+binsize);
14300 abinsizeOn=binsize;
14301 if (nbins) binsizecos=(cos(AminOn*pi/180.)-cos(AmaxOn*pi/180.))/
float(nbins);
14304 title.Form(
"Reconstructed %-s cos(opening angle) of on-source events in time window;cos(opening angle);Counts per %-.3g",scrp.Data(),binsizecos);
14305 hOnCosa=
new TH1F(
"hOnCosa",title,nbins+1,cos(AmaxOn*pi/180.),cos(AminOn*pi/180.)+binsizecos);
14310 binsizecos=fabs(cos(AmaxOn*pi/180.))/float(nbins);
14311 title.Form(
"Reconstructed %-s cos(opening angle) of on-source events in time window;cos(opening angle);Counts per %-.3g",scrp.Data(),binsizecos);
14312 hOnCosa=
new TH1F(
"hOnCosa",title,nbins+2,cos(AmaxOn*pi/180.)-binsizecos,cos(AminOn*pi/180.)+binsizecos);
14317 title.Form(
"Reconstructed opening angle of on-source simulated signal events in time window;Opening angle (degrees);Counts per %-.3g degrees",binsize);
14318 hOnSiga=
new TH1F(
"hOnSiga",title,nbins+2,AminOn-binsize,AmaxOn+binsize);
14319 title.Form(
"Reconstructed cos(opening angle) of on-source simulated signal events in time window;cos(opening angle);Counts per %-.3g",binsizecos);
14320 hOnSigcosa=
new TH1F(
"hOnSigcosa",title,nbins+2,cos(AmaxOn*pi/180.)-binsizecos,cos(AminOn*pi/180.)+binsizecos);
14330 nbins=int((AmaxOff-AminOff)/binsize);
14334 nbins=int(((AmaxOff-AminOff)/180.)*NbkgWinOff*
float(fNgrbs)*
float(fNbkg)/fabs(fAbin));
14335 if (nbins) binsize=(AmaxOff-AminOff)/
float(nbins);
14339 title.Form(
"Reconstructed opening angle of off-source events in time window;Opening angle (degrees);Counts per %-.3g degrees",binsize);
14340 hOffa=
new TH1F(
"hOffa",title,nbins+1,AminOff,AmaxOff+binsize);
14345 binsize=AmaxOff/float(nbins);
14346 title.Form(
"Reconstructed opening angle of off-source events in time window;Opening angle (degrees);Counts per %-.3g degrees",binsize);
14347 hOffa=
new TH1F(
"hOffa",title,nbins+2,AminOff-binsize,AmaxOff+binsize);
14350 abinsizeOff=binsize;
14351 if (nbins) binsizecos=(cos(AminOff*pi/180.)-cos(AmaxOff*pi/180.))/
float(nbins);
14354 title.Form(
"Reconstructed cos(opening angle) of off-source events in time window;cos(opening angle);Counts per %-.3g",binsizecos);
14355 hOffCosa=
new TH1F(
"hOffCosa",title,nbins+1,cos(AmaxOff*pi/180.),cos(AminOff*pi/180.)+binsizecos);
14360 binsizecos=fabs(cos(AmaxOff*pi/180.))/float(nbins);
14361 title.Form(
"Reconstructed cos(opening angle) of off-source events in time window;cos(opening angle);Counts per %-.3g",binsizecos);
14362 hOffCosa=
new TH1F(
"hOffCosa",title,nbins+2,cos(AmaxOff*pi/180.)-binsizecos,cos(AminOff*pi/180.)+binsizecos);
14385 addtitle.Form(
" (=%-g*<T90>)",fTbin);
14386 if (fTbin<0) addtitle.Form(
" (~%-g counts/bin)",fabs(fTbin));
14390 nbins=int(NbkgWinOn*
float(fNgrbs)/fabs(fTbin));
14391 Int_t temp=int(NbkgWinOff*
float(fNgrbs)*
float(fNbkg)/fabs(fTbin));
14392 if (temp>nbins) nbins=temp;
14393 if (nbins) binsize=DtwinOn/float(nbins);
14398 if (fTbint90) binsize=fTbin*fabs(fAvgrbt90)/fTfact;
14399 nbins=DtwinOn/binsize;
14403 title.Form(
"Arrival times %-s of on-source events in time window;Event arrival time (in %-s) w.r.t. burst trigger;Counts per %-.3g %-s",scrt.Data(),tu.Data(),binsize,tu.Data());
14404 if (fTbint90 || fTbin<0) title+=addtitle;
14405 hOnt=
new TH1F(
"hOnt",title,nbins+2,TminOn-binsize,TmaxOn+binsize);
14406 title.Form(
"Arrival time %-s vs. reconstructed %-s opening angle of on-source events in time window;Opening angle (degrees);Event arrival time (in %-s) w.r.t. burst trigger",scrt.Data(),scrp.Data(),tu.Data());
14407 hOnta=
new TH2F(
"hOnta",title,nabinsOn+2,AminOn-abinsizeOn,AmaxOn+abinsizeOn,nbins+2,TminOn-binsize,TmaxOn+binsize);
14408 title.Form(
"Redshift corrected arrival times %-s of on-source events in time window;Event arrival time (in %-s) w.r.t. burst trigger;Counts per %-.3g %-s",scrt.Data(),tu.Data(),binsize,tu.Data());
14409 if (fTbint90 || fTbin<0) title+=addtitle;
14410 hOnZt=
new TH1F(
"hOnZt",title,nbins+2,TminOn-binsize,TmaxOn+binsize);
14411 title.Form(
"Redshift corrected %-s arrival time vs. reconstructed %-s opening angle of on-source events in time window;Opening angle (degrees);Event arrival time (in %-s) w.r.t. burst trigger",scrt.Data(),scrp.Data(),tu.Data());
14412 hOnZta=
new TH2F(
"hOnZta",title,nabinsOn+2,AminOn-abinsizeOn,AmaxOn+abinsizeOn,nbins+2,TminOn-binsize,TmaxOn+binsize);
14416 title.Form(
"Arrival times of on-source simulated signal events in time window;Event arrival time (in %-s) w.r.t. burst trigger;Counts per %-.3g %-s",tu.Data(),binsize,tu.Data());
14417 if (fTbint90 || fTbin<0) title+=addtitle;
14418 hOnSigt=
new TH1F(
"hOnSigt",title,nbins+2,TminOn-binsize,TmaxOn+binsize);
14419 title.Form(
"Arrival time vs. reconstructed opening angle of on-source simulated signal events in time window;Opening angle (degrees);Event arrival time (in %-s) w.r.t. burst trigger",tu.Data());
14420 hOnSigta=
new TH2F(
"hOnSigta",title,nabinsOn+2,AminOn-abinsizeOn,AmaxOn+abinsizeOn,nbins+2,TminOn-binsize,TmaxOn+binsize);
14421 title.Form(
"Redshift corrected arrival times of on-source simulated signal events in time window;Event arrival time (in %-s) w.r.t. burst trigger;Counts per %-.3g %-s",tu.Data(),binsize,tu.Data());
14422 if (fTbint90 || fTbin<0) title+=addtitle;
14423 hOnSigZt=
new TH1F(
"hOnSigZt",title,nbins+2,TminOn-binsize,TmaxOn+binsize);
14424 title.Form(
"Redshift corrected arrival time vs. reconstructed opening angle of on-source simulated signal events in time window;Opening angle (degrees);Event arrival time (in %-s) w.r.t. burst trigger",tu.Data());
14425 hOnSigZta=
new TH2F(
"hOnSigZta",title,nabinsOn+2,AminOn-abinsizeOn,AmaxOn+abinsizeOn,nbins+2,TminOn-binsize,TmaxOn+binsize);
14431 if (fTbint90) binsize=fTbin*fabs(fAvgrbt90)/fTfact;
14432 nbins=DtwinOff/binsize;
14436 title.Form(
"Arrival times of off-source events in time window;Event arrival time (in %-s) w.r.t. burst trigger;Counts per %-.3g %-s",tu.Data(),binsize,tu.Data());
14437 if (fTbint90 || fTbin<0) title+=addtitle;
14438 hOfft=
new TH1F(
"hOfft",title,nbins+2,TminOff-binsize,TmaxOff+binsize);
14439 title.Form(
"Arrival time vs. reconstructed opening angle of off-source events in time window;Opening angle (degrees);Event arrival time (in %-s) w.r.t. burst trigger",tu.Data());
14440 hOffta=
new TH2F(
"hOffta",title,nabinsOff+2,AminOff-abinsizeOff,AmaxOff+abinsizeOff,nbins+2,TminOff-binsize,TmaxOff+binsize);
14441 title.Form(
"Redshift corrected arrival times of off-source events in time window;Event arrival time (in %-s) w.r.t. burst trigger;Counts per %-.3g %-s",tu.Data(),binsize,tu.Data());
14442 if (fTbint90 || fTbin<0) title+=addtitle;
14443 hOffZt=
new TH1F(
"hOffZt",title,nbins+2,TminOff-binsize,TmaxOff+binsize);
14444 title.Form(
"Redshift corrected arrival time vs. reconstructed opening angle of off-source events in time window;Opening angle (degrees);Event arrival time (in %-s) w.r.t. burst trigger",tu.Data());
14445 hOffZta=
new TH2F(
"hOffZta",title,nabinsOff+2,AminOff-abinsizeOff,AmaxOff+abinsizeOff,nbins+2,TminOff-binsize,TmaxOff+binsize);
14450 Double_t* binarr=0;
14451 Int_t nbx=int(DtwinOn/fVarTbin);
14452 Float_t gamma=fabs(fAvgrbz)+1.;
14453 Float_t* bins=
new Float_t[nbx];
14455 Float_t xlow=0,xup=0,size=fVarTbin;
14456 for (Int_t i=0; i<nbx-1; i++)
14459 if (xup>DtwinOn/2.)
14470 binarr=
new Double_t[2*nbins-1];
14471 for (Int_t j=nbins; j>0; j--)
14473 binarr[nbins-j]=-bins[j-1];
14474 binarr[nbins+j-2]=bins[j-1];
14479 title.Form(
"Arrival times %-s of on-source events in time window;Event arrival time (in %-s) w.r.t. burst trigger;Counts per time bin",scrt.Data(),tu.Data());
14480 hOnt=
new TH1F(
"hOnt",title,nbins,binarr);
14483 title.Form(
"Arrival time %-s vs. reconstructed %-s opening angle of on-source events in time window;Opening angle (degrees);Event arrival time (in %-s) w.r.t. burst trigger",scrt.Data(),scrp.Data(),tu.Data());
14484 hOnta=
new TH2F(
"hOnta",title,nabinsOn,AminOn,AmaxOn,nbins,binarr);
14489 title.Form(
"Arrival times of off-source events in time window;Event arrival time (in %-s) w.r.t. burst trigger;Counts per time bin",tu.Data());
14490 hOfft=
new TH1F(
"hOfft",title,nbins,binarr);
14493 title.Form(
"Arrival time vs. reconstructed opening angle of off-source events in time window;Opening angle (degrees);Event arrival time (in %-s) w.r.t. burst trigger",tu.Data());
14494 hOffta=
new TH2F(
"hOffta",title,nabinsOff,AminOff,AmaxOff,nbins,binarr);
14499 title.Form(
"Arrival times of on-source simulated signal events in time window;Event arrival time (in %-s) w.r.t. burst trigger;Counts per time bin",tu.Data());
14500 hOnSigt=
new TH1F(
"hOnSigt",title,nbins,binarr);
14503 title.Form(
"Arrival time vs. reconstructed opening angle of on-source simulated signal events in time window;Opening angle (degrees);Event arrival time (in %-s) w.r.t. burst trigger",tu.Data());
14504 hOnSigta=
new TH2F(
"hOnSigta",title,nabinsOn,AminOn,AmaxOn,nbins,binarr);
14524 if (hOnSourceZ) hOnSourceZ->Fill(value1);
14525 if (hOnSigmaSource) hOnSigmaSource->Fill(value2);
14526 if (hOnSigmaReco) hOnSigmaReco->Fill(value3);
14527 if (hOnSigmaComb) hOnSigmaComb->Fill(value4);
14536 if (hOnE) hOnE->Fill(value1);
14537 if (hOna) hOna->Fill(value2);
14538 if (hOnCosa) hOnCosa->Fill(cos(value2*pi/180.));
14539 if (hOnt) hOnt->Fill(value3);
14540 if (hOnta) hOnta->Fill(value2,value3);
14541 if (hOnZt) hOnZt->Fill(value4);
14542 if (hOnZta) hOnZta->Fill(value2,value4);
14552 if (hOffSourceZ) hOffSourceZ->Fill(value1);
14553 if (hOffSigmaSource) hOffSigmaSource->Fill(value2);
14554 if (hOffSigmaReco) hOffSigmaReco->Fill(value3);
14555 if (hOffSigmaComb) hOffSigmaComb->Fill(value4);
14564 if (hOffE) hOffE->Fill(value1);
14565 if (hOffa) hOffa->Fill(value2);
14566 if (hOffCosa) hOffCosa->Fill(cos(value2*pi/180.));
14567 if (hOfft) hOfft->Fill(value3);
14568 if (hOffta) hOffta->Fill(value2,value3);
14569 if (hOffZt) hOffZt->Fill(value4);
14570 if (hOffZta) hOffZta->Fill(value2,value4);
14580 if (hOnSigSourceZ) hOnSigSourceZ->Fill(value1);
14581 if (hOnSigSigmaSource) hOnSigSigmaSource->Fill(value2);
14582 if (hOnSigSigmaReco) hOnSigSigmaReco->Fill(value3);
14583 if (hOnSigSigmaComb) hOnSigSigmaComb->Fill(value4);
14592 if (hOnSigE) hOnSigE->Fill(value1);
14593 if (hOnSiga) hOnSiga->Fill(value2);
14594 if (hOnSigcosa) hOnSigcosa->Fill(cos(value2*pi/180.));
14595 if (hOnSigt) hOnSigt->Fill(value3);
14596 if (hOnSigta) hOnSigta->Fill(value2,value3);
14597 if (hOnSigZt) hOnSigZt->Fill(value4);
14598 if (hOnSigZta) hOnSigZta->Fill(value2,value4);
14608 title.Form(
"Bayesian blocks for on-source arrival times %-s;Event arrival time (in %-s) w.r.t. burst trigger;Event rate (%-s^{-1}) for %-i stacked time windows",scrt.Data(),tu.Data(),tu.Data(),fNgrbs);
14609 title.ReplaceAll(
"days^",
"day^");
14610 title.ReplaceAll(
"hours^",
"hour^");
14611 hOnBBt.SetNameTitle(
"hOnBBt",title);
14613 title.Form(
"Bayesian blocks for redshift corrected on-source arrival times %-s;Event arrival time (in %-s) w.r.t. burst trigger;Event rate (%-s^{-1}) for %-i stacked time windows",scrt.Data(),tu.Data(),tu.Data(),fNgrbs);
14614 title.ReplaceAll(
"days^",
"day^");
14615 title.ReplaceAll(
"hours^",
"hour^");
14616 hOnBBzt.SetNameTitle(
"hOnBBzt",title);
14623 title.Form(
"Bayesian blocks for off-source arrival times;Event arrival time (in %-s) w.r.t. burst trigger;Event rate (%-s^{-1}) for %-i stacked time windows",tu.Data(),tu.Data(),fNgrbs*fNbkg);
14624 title.ReplaceAll(
"days^",
"day^");
14625 title.ReplaceAll(
"hours^",
"hour^");
14626 hOffBBt.SetNameTitle(
"hOffBBt",title);
14628 title.Form(
"Bayesian blocks for redshift corrected off-source arrival times;Event arrival time (in %-s) w.r.t. burst trigger;Event rate (%-s^{-1}) for %-i stacked time windows",tu.Data(),tu.Data(),fNgrbs*fNbkg);
14629 title.ReplaceAll(
"days^",
"day^");
14630 title.ReplaceAll(
"hours^",
"hour^");
14631 hOffBBzt.SetNameTitle(
"hOffBBzt",title);
14639 axis=hOnBBt.GetXaxis();
14640 xmin=axis->GetXmin();
14641 xmax=axis->GetXmax();
14642 axis=hOffBBt.GetXaxis();
14643 temp=axis->GetXmin();
14644 if (temp<xmin) xmin=temp;
14645 temp=axis->GetXmax();
14646 if (temp>xmax) xmax=temp;
14651 title.Form(
"Event rate (%-s^{-1}) scaled per time window",tu.Data());
14652 title.ReplaceAll(
"days^",
"day^");
14653 title.ReplaceAll(
"hours^",
"hour^");
14654 nb1=BB.
Rebin(&hOnBBt,&hOnUBBt,kFALSE,0,xmin,xmax);
14655 nb2=BB.
Rebin(&hOffBBt,&hOffUBBt,kFALSE,0,xmin,xmax);
14656 if (nb2>nb1) BB.
Rebin(&hOnBBt,&hOnUBBt,kFALSE,nb2,xmin,xmax);
14657 hOnUBBt.SetName(
"hOnUBBt");
14658 hOffUBBt.SetName(
"hOffUBBt");
14659 if (fNgrbs>1) hOnUBBt.Scale(1./
float(fNgrbs));
14660 if (fNgrbs*fNbkg>1) hOffUBBt.Scale(1./
float(fNgrbs*fNbkg));
14661 hOnUBBt.SetYTitle(title);
14662 hOffUBBt.SetYTitle(title);
14663 BB.
Divide(&hOnUBBt,&hOffUBBt,&hRatUBBt,kFALSE,1);
14664 hRatUBBt.SetName(
"hRatUBBt");
14665 hRatUBBt.SetYTitle(
"Ratio");
14667 nbins=hOnBBzt.GetNbinsX();
14668 if (hOffBBzt.GetNbinsX()>nbins) nbins=hOffBBzt.GetNbinsX();
14669 axis=hOnBBzt.GetXaxis();
14670 xmin=axis->GetXmin();
14671 xmax=axis->GetXmax();
14672 axis=hOffBBzt.GetXaxis();
14673 temp=axis->GetXmin();
14674 if (temp<xmin) xmin=temp;
14675 temp=axis->GetXmax();
14676 if (temp>xmax) xmax=temp;
14681 nb1=BB.
Rebin(&hOnBBzt,&hOnUBBzt,kFALSE,0,xmin,xmax);
14682 nb2=BB.
Rebin(&hOffBBzt,&hOffUBBzt,kFALSE,0,xmin,xmax);
14683 if (nb2>nb1) BB.
Rebin(&hOnBBzt,&hOnUBBzt,kFALSE,nb2,xmin,xmax);
14684 hOnUBBzt.SetName(
"hOnUBBzt");
14685 hOffUBBzt.SetName(
"hOffUBBzt");
14686 if (fNgrbs>1) hOnUBBzt.Scale(1./
float(fNgrbs));
14687 if (fNgrbs*fNbkg>1) hOffUBBzt.Scale(1./
float(fNgrbs*fNbkg));
14688 hOnUBBzt.SetYTitle(title);
14689 hOffUBBzt.SetYTitle(title);
14690 BB.
Divide(&hOnUBBzt,&hOffUBBzt,&hRatUBBzt,kFALSE,1);
14691 hRatUBBzt.SetName(
"hRatUBBzt");
14692 hRatUBBzt.SetYTitle(
"Ratio");
14699 if (hOffSigmaSource)
fBurstHistos.Add(hOffSigmaSource);
14700 if (hOnSigSigmaSource)
fBurstHistos.Add(hOnSigSigmaSource);
14703 if (hOnSigSigmaReco)
fBurstHistos.Add(hOnSigSigmaReco);
14706 if (hOnSigSigmaComb)
fBurstHistos.Add(hOnSigSigmaComb);
14722 if (hOnt && hOnBBt.GetEntries())
fBurstHistos.Add(hOnBBt.Clone());
14723 if (hOfft && hOffBBt.GetEntries())
fBurstHistos.Add(hOffBBt.Clone());
14724 if (hOnt && hOnUBBt.GetEntries())
fBurstHistos.Add(hOnUBBt.Clone());
14725 if (hOfft && hOffUBBt.GetEntries())
fBurstHistos.Add(hOffUBBt.Clone());
14726 if (hOnt && hOfft && hRatUBBt.GetEntries())
fBurstHistos.Add(hRatUBBt.Clone());
14733 if (hOnZt && hOnBBzt.GetEntries())
fBurstHistos.Add(hOnBBzt.Clone());
14734 if (hOffZt && hOffBBzt.GetEntries())
fBurstHistos.Add(hOffBBzt.Clone());
14735 if (hOnZt && hOnUBBzt.GetEntries())
fBurstHistos.Add(hOnUBBzt.Clone());
14736 if (hOffZt && hOffUBBzt.GetEntries())
fBurstHistos.Add(hOffUBBzt.Clone());
14737 if (hOnZt && hOffZt && hRatUBBzt.GetEntries())
fBurstHistos.Add(hRatUBBzt.Clone());
14741 for (Int_t ih=0; ih<nh; ih++)
14745 hx->BufferEmpty(1);
14747 name=hx->GetName();
14749 if (name==
"hTdiff")
14751 binsize=hx->GetBinWidth(1);
14752 s.Form(
"Counts per %-.3g %-s",binsize,tu.Data());
14753 axis=hx->GetYaxis();
14754 if (axis) axis->SetTitle(s);
14757 if (name==
"hAdiff" || name.Contains(
"Sigma"))
14759 binsize=hx->GetBinWidth(1);
14760 s.Form(
"Counts per %-.3g degrees",binsize);
14761 axis=hx->GetYaxis();
14762 if (axis) axis->SetTitle(s);
14765 if (name==
"hOnE" || name==
"hOffE" || name==
"hOnSigE")
14767 binsize=hx->GetBinWidth(1);
14768 s.Form(
"Counts per %-.3g GeV",binsize);
14769 axis=hx->GetYaxis();
14770 if (axis) axis->SetTitle(s);
14777 Float_t nsigOn=nsig;
14778 Float_t nbkgOn=nOn-nsigOn;
14779 Float_t Ton=DtwinOn*fTfact*float(fNgrbs);
14780 Float_t Toff=DtwinOff*fTfact*float(fNgrbs)*float(fNbkg);
14781 Float_t rateOn=nOn/(DtwinOn*fTfact);
14782 Float_t rateOff=nOff/(DtwinOff*fTfact);
14785 Float_t TwinOn=DtwinOn*fTfact;
14786 Float_t TwinOff=DtwinOff*fTfact;
14787 Float_t AvSolidangleOn=0;
14788 Float_t AvSolidangleOff=0;
14795 Float_t AvNsigOn=0;
14796 Float_t AvNbkgOn=0;
14797 Float_t AvRsigOn=0;
14798 Float_t AvRbkgOn=0;
14799 Float_t AvEsigOn=0;
14800 Float_t AvEbkgOn=0;
14801 Float_t scale=fNgrbs;
14804 AvSolidangleOn=fSolidangleOn/scale;
14807 AvEon=fEnergyOn/scale;
14808 AvNsigOn=nsigOn/scale;
14809 AvNbkgOn=nbkgOn/scale;
14810 AvRsigOn=nsigOn/Ton;
14811 AvRbkgOn=nbkgOn/Ton;
14812 AvEsigOn=fEnergySig/scale;
14813 AvEbkgOn=fEnergyBkg/scale;
14815 scale=fNgrbs*fNbkg;
14818 AvSolidangleOff=fSolidangleOff/scale;
14821 AvEoff=fEnergyOff/scale;
14825 printf(
"\n *%-s::MakeBurstDataStats* Statistics of the stacked observed event samples. \n",ClassName());
14826 printf(
" Integrated on-source exposure time of the %-i stacked time windows : %-g %-s \n",fNgrbs,Ton/fTfact,tu.Data());
14827 printf(
" Integrated off-source exposure time of the %-i*%-i stacked time windows : %-g %-s \n",fNgrbs,fNbkg,Toff/fTfact,tu.Data());
14828 printf(
" Total accumulated on-source solid angle : %-g sr in %-i stacked patches --> Average per patch : %-g sr \n",fSolidangleOn,fNgrbs,AvSolidangleOn);
14829 printf(
" Total accumulated off-source solid angle : %-g sr in %-i*%-i stacked patches --> Average per patch : %-g sr \n",fSolidangleOff,fNgrbs,fNbkg,AvSolidangleOff);
14830 if (fSensarea>0) printf(
" Area covered c.q. overlooked by the detector sensors : %-g m^2. \n",fSensarea);
14833 printf(
" *On source* Total number of recorded on-source events : %-g --> Average per patch : %-g events. \n",nOn,AvNon);
14834 printf(
" --- Average values per on-source patch --- \n");
14835 printf(
" Steady event rate during the time window : %-g Hz",AvRon);
14836 if (AvSolidangleOn) printf(
" --> %-g Hz sr^-1",AvRon/AvSolidangleOn);
14840 printf(
" Particle fluence : %-g cm^-2",AvNon/fSensarea);
14841 if (AvSolidangleOn) printf(
" --> %-g cm^-2 sr^-1",AvNon/(fSensarea*AvSolidangleOn));
14843 printf(
" Particle flux : %-g cm^-2 s^-1",AvRon/fSensarea);
14844 if (AvSolidangleOn) printf(
" --> Intensity : %-g cm^-2 s^-1 sr^-1",AvRon/(fSensarea*AvSolidangleOn));
14847 printf(
" Cumulated energy : %-g GeV",AvEon);
14848 if (AvSolidangleOn) printf(
" --> %-g GeV sr^-1",AvEon/AvSolidangleOn);
14850 printf(
" Power : %-g GeV/s",AvEon/TwinOn);
14851 if (AvSolidangleOn) printf(
" --> %-g GeV s^-1 sr^-1",AvEon/(TwinOn*AvSolidangleOn));
14855 printf(
" Energy fluence : %-g GeV cm^-2",AvEon/fSensarea);
14856 if (AvSolidangleOn) printf(
" --> %-g GeV cm^-2 sr^-1",AvEon/(fSensarea*AvSolidangleOn));
14858 printf(
" Energy flux : %-g GeV cm^-2 s^-1",AvEon/(fSensarea*TwinOn));
14859 if (AvSolidangleOn) printf(
" --> Intensity : %-g GeV cm^-2 s^-1 sr^-1",AvEon/(fSensarea*TwinOn*AvSolidangleOn));
14865 printf(
" *Off source* Total number of recorded off-source (background) events : %-g --> Average per patch : %-g events. \n",nOff,AvNoff);
14866 printf(
" --- Average values per off-source patch --- \n");
14867 printf(
" Steady event rate during the time window : %-g Hz",AvRoff);
14868 if (AvSolidangleOff) printf(
" --> %-g Hz sr^-1",AvRoff/AvSolidangleOff);
14872 printf(
" Particle fluence : %-g cm^-2",AvNoff/fSensarea);
14873 if (AvSolidangleOff) printf(
" --> %-g cm^-2 sr^-1",AvNoff/(fSensarea*AvSolidangleOff));
14875 printf(
" Particle flux : %-g cm^-2 s^-1",AvRoff/fSensarea);
14876 if (AvSolidangleOff) printf(
" --> Intensity : %-g cm^-2 s^-1 sr^-1",AvRoff/(fSensarea*AvSolidangleOff));
14879 printf(
" Cumulated energy : %-g GeV",AvEoff);
14880 if (AvSolidangleOff) printf(
" --> %-g GeV sr^-1",AvEoff/AvSolidangleOff);
14882 printf(
" Power : %-g GeV/s",AvEoff/TwinOff);
14883 if (AvSolidangleOff) printf(
" --> %-g GeV s^-1 sr^-1",AvEoff/(TwinOn*AvSolidangleOff));
14887 printf(
" Energy fluence : %-g GeV cm^-2",AvEoff/fSensarea);
14888 if (AvSolidangleOff) printf(
" --> %-g GeV cm^-2 sr^-1",AvEoff/(fSensarea*AvSolidangleOff));
14890 printf(
" Energy flux : %-g GeV cm^-2 s^-1",AvEoff/(fSensarea*TwinOff));
14891 if (AvSolidangleOff) printf(
" --> Intensity : %-g GeV cm^-2 s^-1 sr^-1",AvEoff/(fSensarea*TwinOff*AvSolidangleOff));
14899 printf(
" -(Unknown)- Total number of injected on-source signal events : %-i \n",nmugrb);
14900 printf(
" Total number of recorded on-source signal events : %-g",nsigOn);
14901 if (fNgrbs) printf(
" --> Average per patch : %-g events.",nsigOn/
float(fNgrbs));
14903 printf(
" Total number of recorded on-source bkg events : %-g",nbkgOn);
14904 if (fNgrbs) printf(
" --> Average per patch : %-g events.",nbkgOn/
float(fNgrbs));
14907 printf(
" --- Average signal values per on-source patch --- \n");
14908 printf(
" Steady event rate during the time window : %-g Hz",AvRsigOn);
14909 if (AvSolidangleOn) printf(
" --> %-g Hz sr^-1",AvRsigOn/AvSolidangleOn);
14913 printf(
" Particle fluence : %-g cm^-2",AvNsigOn/fSensarea);
14914 if (AvSolidangleOn) printf(
" --> %-g cm^-2 sr^-1",AvNsigOn/(fSensarea*AvSolidangleOn));
14916 printf(
" Particle flux : %-g cm^-2 s^-1",AvRsigOn/fSensarea);
14917 if (AvSolidangleOn) printf(
" --> Intensity : %-g cm^-2 s^-1 sr^-1",AvRsigOn/(fSensarea*AvSolidangleOn));
14920 printf(
" Cumulated energy : %-g GeV",AvEsigOn);
14921 if (AvSolidangleOn) printf(
" --> %-g GeV sr^-1",AvEsigOn/AvSolidangleOn);
14923 printf(
" Power : %-g GeV/s",AvEsigOn/TwinOn);
14924 if (AvSolidangleOn) printf(
" --> %-g GeV s^-1 sr^-1",AvEsigOn/(TwinOn*AvSolidangleOn));
14928 printf(
" Energy fluence : %-g GeV cm^-2",AvEsigOn/fSensarea);
14929 if (AvSolidangleOn) printf(
" --> %-g GeV cm^-2 sr^-1",AvEsigOn/(fSensarea*AvSolidangleOn));
14931 printf(
" Energy flux : %-g GeV cm^-2 s^-1",AvEsigOn/(fSensarea*TwinOn));
14932 if (AvSolidangleOn) printf(
" --> Intensity : %-g GeV cm^-2 s^-1 sr^-1",AvEsigOn/(fSensarea*TwinOn*AvSolidangleOn));
14935 printf(
" --- Average background values per on-source patch --- \n");
14936 printf(
" Steady event rate during the time window : %-g Hz",AvRbkgOn);
14937 if (AvSolidangleOn) printf(
" --> %-g Hz sr^-1",AvRbkgOn/AvSolidangleOn);
14941 printf(
" Particle fluence : %-g cm^-2",AvNbkgOn/fSensarea);
14942 if (AvSolidangleOn) printf(
" --> %-g cm^-2 sr^-1",AvNbkgOn/(fSensarea*AvSolidangleOn));
14944 printf(
" Particle flux : %-g cm^-2 s^-1",AvRbkgOn/fSensarea);
14945 if (AvSolidangleOn) printf(
" --> Intensity : %-g cm^-2 s^-1 sr^-1",AvRbkgOn/(fSensarea*AvSolidangleOn));
14948 printf(
" Cumulated energy : %-g GeV",AvEbkgOn);
14949 if (AvSolidangleOn) printf(
" --> %-g GeV sr^-1",AvEbkgOn/AvSolidangleOn);
14951 printf(
" Power : %-g GeV/s",AvEbkgOn/TwinOn);
14952 if (AvSolidangleOn) printf(
" --> %-g GeV s^-1 sr^-1",AvEbkgOn/(TwinOn*AvSolidangleOn));
14956 printf(
" Energy fluence : %-g GeV cm^-2",AvEbkgOn/fSensarea);
14957 if (AvSolidangleOn) printf(
" --> %-g GeV cm^-2 sr^-1",AvEbkgOn/(fSensarea*AvSolidangleOn));
14959 printf(
" Energy flux : %-g GeV cm^-2 s^-1",AvEbkgOn/(fSensarea*TwinOn));
14960 if (AvSolidangleOn) printf(
" --> Intensity : %-g GeV cm^-2 s^-1 sr^-1",AvEbkgOn/(fSensarea*TwinOn*AvSolidangleOn));
15124 if (fTunits==1) tu=
"hrs";
15125 if (fTunits==2) tu=
"s";
15126 if (fTunits==3) tu=
"ns";
15127 if (fTunits==4) tu=
"ps";
15133 TString name=
"Matches";
15134 TString title=
"Space and time matchings of NcAstrolab stored signals";
15135 matches.SetNameTitle(name,title);
15137 if (tu==
"d") tux=
"days";
15138 if (tu==
"hrs") tux=
"hours";
15139 if (tu==
"s") tux=
"sec";
15140 TString namedamin=
"psimin in deg";
15141 TString namedtmin=
"dtmin in ";
15149 TString nameda=
"psi in deg";
15150 TString namedt=
"t2-t1 in ";
15159 if ((!itype || !jtype) && !
fRefs)
15161 printf(
" *%-s::MatchBurstData* Error: itype=%-i jtype=%-i but no reference signals are present. \n",ClassName(),itype,jtype);
15165 if ((itype || jtype) && !
fSigs)
15167 printf(
" *%-s::MatchBurstData* Error: itype=%-i jtype=%-i but no measurements are present. \n",ClassName(),itype,jtype);
15177 if (itype) itype=1;
15178 if (jtype) jtype=1;
15181 if (i2<1 || i2>nrefs) i2=nrefs;
15185 if (i2<1 || i2>nsigs) i2=nsigs;
15189 if (j2<1 || j2>nrefs) j2=nrefs;
15193 if (j2<1 || j2>nsigs) j2=nsigs;
15196 if (i1<1 || j1<1 || i1>i2 || j1>j2)
15198 printf(
" *%-s::MatchBurstData* Inconsistent parameters: i1=%-i i2=%-i itype=%-i j1=%-i j2=%-i jtype=%-i. \n",ClassName(),i1,i2,itype,j1,j2,jtype);
15202 if (fDatype==1 && fMaxsigmatot<0)
15204 printf(
" *%-s::MatchBurstData* Incompatible parameter settings Datype=%-i Maxsigmatot=%-g \n",ClassName(),fDatype,fMaxsigmatot);
15205 printf(
" === No matching analysis will be performed === \n");
15215 if (
fTscmode>0) scrt=
"(scrambled)";
15220 title.Form(
"Time difference %-s between events and bursts for the full selected dataset;Tevent-Tburst in %-s;Counts",scrt.Data(),tux.Data());
15221 TH1F* hTdiff=
new TH1F(
"hTdiff",title,nbins,1,0);
15222 if (ni && nj) hTdiff->SetBuffer(ni*nj);
15226 title.Form(
"Angular separation %-s between events and bursts for the full selected dataset;Opening angle in degrees;Counts",scrp.Data());
15227 TH1F* hAdiff=
new TH1F(
"hAdiff",title,nbins,1,0);
15228 if (ni && nj) hAdiff->SetBuffer(ni*nj);
15234 Double_t dang,dtime,diftheta;
15243 Float_t sigmagrb=0;
15244 Float_t sigmareco=0;
15245 Float_t sigmatot=0;
15247 Double_t dangmin=0;
15251 Bool_t first=kTRUE;
15252 Double_t dangmax=-1;
15256 Float_t thetagrb=0;
15257 Float_t solidangle=0;
15262 for (Int_t bkgpatch=0; bkgpatch<=fNbkg; bkgpatch++)
15266 for (Int_t i=i1; i<=i2; i++)
15269 if (!sxi)
continue;
15273 if (itype && !fRecoangle) sigmai=fAngresfix;
15275 for (Int_t j=j1; j<=j2; j++)
15278 if (itype==jtype && i==j)
continue;
15281 if (!sxj)
continue;
15285 if (jtype && !fRecoangle) sigmaj=fAngresfix;
15313 if (fSumsigmas==-1) sigmatot=sigmareco;
15314 if (fSumsigmas==0) sigmatot=sigmagrb;
15315 if (fSumsigmas==1) sigmatot=sigmagrb+sigmareco;
15316 if (fSumsigmas==2) sigmatot=sqrt(sigmagrb*sigmagrb+sigmareco*sigmareco);
15319 if (!fDatype) dangmax=fabs(fDawin);
15323 if (fMaxsigmatot>0) dangmax=fabs(fDawin*fMaxsigmatot);
15327 if (sigmatot>=0) dangmax=fabs(fDawin*sigmatot);
15331 if (fDatype==1 && fMaxsigmatot<0) dangmax=-1;
15340 dang=
GetSeparation(ix,jx,
"deg",dtime,tu,1,bkgpatch,&diftheta);
15344 dang=
GetSeparation(ix,jx,
"deg",dtime,
"s",1,bkgpatch,&diftheta);
15351 if (hTdiff) hTdiff->Fill(dtime);
15354 if (hAdiff) hAdiff->Fill(dang);
15362 thlow=thetagrb-0.5*dangmax;
15363 thup=thetagrb+0.5*dangmax;
15372 sxgrb->
SetSignal(solidangle,
"OmegaOn");
15383 thlow=thetagrb-0.5*dangmax;
15384 thup=thetagrb+0.5*dangmax;
15393 sxgrb->
SetSignal(solidangle,
"OmegaOff");
15401 if (diftheta<thlow || diftheta>thup)
continue;
15405 if (fabs(dang)>dangmax)
continue;
15409 if (fTmax>fTmin && (dtime<fTmin || dtime>fTmax))
continue;
15422 name+=sxi->GetName();
15424 title+=sxj->GetName();
15426 data.SetNameTitle(name,title);
15427 data.SetUniqueID(
id);
15441 thlow=thetagrb-0.5*dangmax;
15442 thup=thetagrb+0.5*dangmax;
15468 thlow=thetagrb-0.5*dangmax;
15469 thup=thetagrb+0.5*dangmax;
15478 if (solidangle>sxgrb->
GetSignal(
"OmegaOff")) sxgrb->
SetSignal(solidangle,
"OmegaOff");
15492 if (first || fabs(dang)<dangmin)
15494 dangmin=fabs(dang);
15499 if (first || fabs(dtime)<fabs(dtmin))
15524 name=
"SolidangleOn";
15527 name=
"SolidangleOff";
15552 for (Int_t k=k1; k<=k2; k++)
15555 if (!sxgrb)
continue;
15557 solidangle=sxgrb->
GetSignal(
"OmegaOn");
15560 solidangle=sxgrb->
GetSignal(
"OmegaOff");
15686 printf(
" *%-s::MatchBurstData* Object %-s not found for itype=%-i. \n",ClassName(),name.Data(),itype);
15725 if (fTunits==1) tu=
"hours";
15726 if (fTunits==2) tu=
"sec";
15727 if (fTunits==3) tu=
"ns";
15728 if (fTunits==4) tu=
"ps";
15740 fBurstOnReco.SetNameTitle(
"BurstOnReco",
"On-source reco data");
15742 fBurstOnReco.SetNames(
"zburst",
"sigmaburst",
"sigmareco",
"sigmacomb");
15743 fBurstOnMatch.SetNameTitle(
"BurstOnMatch",
"On-source matching data");
15746 fBurstSigReco.SetNameTitle(
"BurstSigReco",
"Simulated signal reco data");
15748 fBurstSigReco.SetNames(
"zburst",
"sigmaburst",
"sigmareco",
"sigmacomb");
15749 fBurstSignal.SetNameTitle(
"BurstSignal",
"Simulated signal events");
15752 fBurstOffReco.SetNameTitle(
"BurstOffReco",
"Off-source reco data");
15754 fBurstOffReco.SetNames(
"zburst",
"sigmaburst",
"sigmareco",
"sigmacomb");
15755 fBurstOffMatch.SetNameTitle(
"BurstOffMatch",
"Off-source matching data");
15781 Int_t srcbufsize=10000;
15782 if (fNgrbs && fNgrbs<srcbufsize) srcbufsize=fNgrbs;
15783 Int_t evtbufsize=10000;
15784 if (fNevts && fNevts<evtbufsize) evtbufsize=fNevts;
15787 title.Form(
"Redshifts for the selected source sample;Redshift;Counts");
15788 TH1F* hSourceZ=
new TH1F(
"hSourceZ",title,nbins,1,0);
15789 hSourceZ->SetBuffer(srcbufsize);
15793 title.Form(
"Distances for the selected source sample derived from the redshifts;Physical distance in Mpc;Counts");
15794 TH1F* hSourceD=
new TH1F(
"hSourceD",title,nbins,xmin,xmax);
15795 hSourceD->SetBuffer(srcbufsize);
15799 if (fabs(fT90min)>0) xmin=log10(fabs(fT90min));
15801 if (fT90max>0) xmax=log10(fT90max);
15804 nbins=TMath::Nint(range/binsize);
15811 title.Form(
"Burst durations for the selected source sample;Burst duration ^{10}log(T90) in sec.;Counts");
15812 TH1F* hBurstT90=
new TH1F(
"hBurstT90",title,nbins,xmin,xmax);
15816 title.Form(
"Position uncertainties for the selected source sample;Position angular uncertainty (sigma in degrees);Counts");
15817 TH1F* hSigmaSource=
new TH1F(
"hSigmaSource",title,nbins,1,0);
15818 hSigmaSource->SetBuffer(srcbufsize);
15825 title.Form(
"Reconstructed energy for the full selected real event sample;Reconstructed event energy in GeV;Counts");
15826 hEreco=
new TH1F(
"hEreco",title,nbins,1,0);
15827 hEreco->SetBuffer(evtbufsize);
15836 title.Form(
"Injected signal energy at the source;Event energy in GeV;Counts");
15837 hSigE=
new TH1F(
"hSigE",title,nbins,1,0);
15839 title.Form(
"(Redshift corrected) injected signal energy arriving at Earth;Event energy in GeV;Counts");
15840 hSigEzcor=
new TH1F(
"hSigEzcor",title,nbins,1,0);
15845 TH1F* hSigmaReco=0;
15847 title.Form(
"Event reconstruction uncertainties for the full selected sample;Event angular reconstruction uncertainty (sigma in degrees);Counts");
15848 hSigmaReco=
new TH1F(
"hSigmaReco",title,nbins,1,0);
15849 hSigmaReco->SetBuffer(evtbufsize);
15856 Float_t sigmagrb=0;
15864 for (Int_t i=1; i<=nsig; i++)
15875 hSourceZ->Fill(zgrb);
15876 hSourceD->Fill(dgrb);
15877 if (t90grb>0) hBurstT90->Fill(log10(t90grb));
15878 hSigmaSource->Fill(sigmagrb);
15880 if (fAvgrbz<0) zsample.
Enter(zgrb);
15881 if (fAvgrbt90<0) t90sample.
Enter(t90grb);
15882 sigmasample.
Enter(sigmagrb);
15886 if (hSourceZ->GetEntries())
fBurstHistos.Add(hSourceZ);
15887 if (hSourceD->GetEntries())
fBurstHistos.Add(hSourceD);
15888 if (hBurstT90->GetEntries())
fBurstHistos.Add(hBurstT90);
15889 if (hSigmaSource->GetEntries())
fBurstHistos.Add(hSigmaSource);
15905 Float_t fAvgrbsigma=sigmasample.
GetMedian(1);
15908 Float_t sigmareco=0;
15911 for (Int_t i=1; i<=nsig; i++)
15920 if (hSigmaReco) hSigmaReco->Fill(sigmareco);
15921 if (hEreco) hEreco->Fill(Ereco);
15925 if (hSigmaReco->GetEntries())
fBurstHistos.Add(hSigmaReco);
15973 Int_t nmu=int(fabs(fGrbnu)*
float(fNgrbs));
15979 Float_t sigmagrb=0;
15984 Double_t thetagrb=0;
15988 Float_t dangmaxOn=0;
15991 Double_t sigmareco=0;
15992 Float_t sigmatot=0;
15996 Float_t solidangle=0;
15997 Int_t fixedwinset=0;
16000 TH1* hSigEzcor=(TH1*)
fBurstHistos.FindObject(
"hSigEzcor");
16005 jgrb=int(
fRan->Uniform(0.,
float(fNgrbs)));
16006 if (jgrb==0) jgrb=1;
16012 GetSignal(dgrb,thetagrb,
"deg",phigrb,
"deg",
"loc",tx,jgrb);
16017 fixedwinset=TMath::Nint(sx->
GetSignal(
"fixedwinset"));
16033 dt=
fRan->Gauss(fDtnu,fabs(fDtnus)*t90grb);
16037 dt=
fRan->Gauss(fDtnu,fDtnus);
16045 dt=
fRan->Gauss(fDtnu*t90grb,fabs(fDtnus)*t90grb);
16049 dt=
fRan->Gauss(fDtnu*t90grb,fDtnus);
16052 if (fTimres>0) dt=
fRan->Gauss(dt,fTimres);
16063 if (hSigE) hSigE->Fill(E);
16066 if (fEzcor) E=E/(zgrb+1.);
16068 if (hSigEzcor) hSigEzcor->Fill(E);
16070 if (E<0 || E<fEmin || E>fEmax)
continue;
16075 Int_t mode=fKinangle-1;
16082 if (sigmareco<0) sigmareco=fAngresfix;
16084 if (sigmareco<fAngresmin || sigmareco>fAngresmax)
continue;
16092 if (fSumsigmas==-1) sigmatot=sigmareco;
16093 if (fSumsigmas==0) sigmatot=sigmagrb;
16094 if (fSumsigmas==1) sigmatot=sigmagrb+sigmareco;
16095 if (fSumsigmas==2) sigmatot=sqrt(sigmagrb*sigmagrb+sigmareco*sigmareco);
16101 if (sigmatot>=0) dangmax=fabs(fDawin*sigmatot);
16104 if (fDatype>=0 && dang>dangmax)
continue;
16113 if (fixedwinset)
continue;
16118 thlow=thetagrb-0.5*dangmax;
16119 thup=thetagrb+0.5*dangmax;
16128 if (dangmax>dangmaxOn) sx->
SetSignal(dangmax,
"dangmaxOn");
16129 if (solidangle>OmegaOn) sx->
SetSignal(solidangle,
"OmegaOn");
16172 if (fNgrbs<=0 || fNbkg<=0 || Non<=0 || Noff<=0)
16174 printf(
" \n *%-s::GetBurstBayesianSignalRate* \n",ClassName());
16175 if (fNgrbs<=0 || Non<=0) printf(
" === No on-source data available === \n");
16176 if (fNbkg<=0 || Noff<=0) printf(
" === No off-source data available === \n");
16184 if (fTunits==1) tu=
"hours";
16185 if (fTunits==2) tu=
"sec";
16186 if (fTunits==3) tu=
"ns";
16187 if (fTunits==4) tu=
"ps";
16195 if (fAvSolidangleOff) Ra=fAvSolidangleOn/fAvSolidangleOff;
16205 Double_t rmode=fsigrpdf.GetMaximumX();
16206 Double_t bkgrmode=fbkgrpdf.GetMaximumX();
16213 fbkgrpdf.SetRange(0,3.*Noff/Toff);
16214 fbkgrpdf.SetNpx(n);
16215 TH1* hBkgRatePDF=(TH1*)fbkgrpdf.GetHistogram()->Clone();
16216 hBkgRatePDF->SetName(
"hBkgRatePDF");
16218 fsigrpdf.SetRange(0,3.*Non/Ton);
16219 fsigrpdf.SetNpx(n);
16220 TH1* hSigRatePDF=(TH1*)fsigrpdf.GetHistogram()->Clone();
16221 hSigRatePDF->SetName(
"hSigRatePDF");
16224 printf(
"\n *%-s::GetBurstBayesianSignalRate* Credible interval [rlow,rup] for p=%-g%% with a precision of 1/%-i \n",ClassName(),p,n);
16232 if (fRecoangle && fSumsigmas && fDatype==2)
16234 printf(
" === Warning: Variable angular cones based on actual track reconstruction sigmas were used. \n");
16235 printf(
" Large variations between the on-source and off-source background counts may be present. \n");
16238 printf(
" The %-g%% credible interval from the Bayesian posterior signal pdf : [%-g,%-g] Hz \n",100.*frac,rlow,rup);
16239 printf(
" Modes of the on-source posterior PDFs : Signal=%-g Hz Background=%-g Hz",rmode,bkgrmode);
16240 if (bkgrmode) printf(
" Signal/Background=%-g",rmode/bkgrmode);
16242 cout <<
" The following signal and background rate PDF histograms have been generated :" << endl;
16243 cout <<
" ... " << hSigRatePDF->GetName() <<
" : " << hSigRatePDF->GetTitle() << endl;
16244 cout <<
" ... " << hBkgRatePDF->GetName() <<
" : " << hBkgRatePDF->GetTitle() << endl;
16248 Float_t TwinOn=Ton/float(fNgrbs);
16250 printf(
" Integrated on-source exposure time of the %-i stacked time windows : %-g %-s \n",fNgrbs,Ton/fTfact,tu.Data());
16251 printf(
" Integrated off-source exposure time of the %-i*%-i stacked time windows : %-g %-s \n",fNgrbs,fNbkg,Toff/fTfact,tu.Data());
16252 printf(
" Total accumulated on-source solid angle : %-g sr in %-i stacked patches --> Average per patch : %-g sr \n",fSolidangleOn,fNgrbs,fAvSolidangleOn);
16253 printf(
" Total accumulated off-source solid angle : %-g sr in %-i*%-i stacked patches --> Average per patch : %-g sr \n",fSolidangleOff,fNgrbs,fNbkg,fAvSolidangleOff);
16254 if (fSensarea>0) printf(
" Area covered c.q. overlooked by the detector sensors : %-g m^2. \n",fSensarea);
16257 printf(
" Total number of recorded on-source events : %-g --> Average per patch : %-g events. \n",Non,Non/
float(fNgrbs));
16258 printf(
" Total number of recorded off-source events : %-g --> Average per patch : %-g events. \n",Noff,Noff/
float(fNgrbs*fNbkg));
16260 printf(
" --- Average values per on-source patch based on the posterior PDFs --- \n");
16262 printf(
" *Lower bound* Steady signal event rate during the time window : %-g Hz",rlow);
16263 if (fAvSolidangleOn) printf(
" --> %-g Hz sr^-1",rlow/fAvSolidangleOn);
16265 printf(
" Expected number of signal events : %-g",rlow*TwinOn);
16266 if (fAvSolidangleOn) printf(
" --> %-g events/sr",rlow*TwinOn/fAvSolidangleOn);
16270 printf(
" Expected signal particle fluence : %-g cm^-2",rlow*TwinOn/fSensarea);
16271 if (fAvSolidangleOn) printf(
" --> : %-g cm^-2 sr^-1",rlow*TwinOn/(fSensarea*fAvSolidangleOn));
16273 printf(
" Expected signal particle flux : %-g cm^-2 s^-1",rlow/fSensarea);
16274 if (fAvSolidangleOn) printf(
" --> Intensity : %-g cm^-2 s^-1 sr^-1",rlow/(fSensarea*fAvSolidangleOn));
16278 printf(
" *Mode* Steady signal event rate during the time window : %-g Hz",rmode);
16279 if (fAvSolidangleOn) printf(
" --> %-g Hz sr^-1",rmode/fAvSolidangleOn);
16281 printf(
" Expected number of signal events : %-g",rmode*TwinOn);
16282 if (fAvSolidangleOn) printf(
" --> %-g events/sr",rmode*TwinOn/fAvSolidangleOn);
16286 printf(
" Expected signal particle fluence : %-g cm^-2",rmode*TwinOn/fSensarea);
16287 if (fAvSolidangleOn) printf(
" --> : %-g cm^-2 sr^-1",rmode*TwinOn/(fSensarea*fAvSolidangleOn));
16289 printf(
" Expected signal particle flux : %-g cm^-2 s^-1",rmode/fSensarea);
16290 if (fAvSolidangleOn) printf(
" --> Intensity : %-g cm^-2 s^-1 sr^-1",rmode/(fSensarea*fAvSolidangleOn));
16294 printf(
" *Upper bound* Steady signal event rate during the time window : %-g Hz",rup);
16295 if (fAvSolidangleOn) printf(
" --> %-g Hz sr^-1",rup/fAvSolidangleOn);
16297 printf(
" Expected number of signal events : %-g",rup*TwinOn);
16298 if (fAvSolidangleOn) printf(
" --> %-g events/sr",rup*TwinOn/fAvSolidangleOn);
16302 printf(
" Expected signal particle fluence : %-g cm^-2",rup*TwinOn/fSensarea);
16303 if (fAvSolidangleOn) printf(
" --> : %-g cm^-2 sr^-1",rup*TwinOn/(fSensarea*fAvSolidangleOn));
16305 printf(
" Expected signal particle flux : %-g cm^-2 s^-1",rup/fSensarea);
16306 if (fAvSolidangleOn) printf(
" --> Intensity : %-g cm^-2 s^-1 sr^-1",rup/(fSensarea*fAvSolidangleOn));
16310 printf(
" *Background* Steady background event rate during the time window : %-g Hz",bkgrmode);
16311 if (fAvSolidangleOn) printf(
" --> %-g Hz sr^-1",bkgrmode/fAvSolidangleOn);
16313 printf(
" Expected number of background events : %-g",bkgrmode*TwinOn);
16314 if (fAvSolidangleOn) printf(
" --> %-g events/sr",bkgrmode*TwinOn/fAvSolidangleOn);
16318 printf(
" Expected background particle fluence : %-g cm^-2",bkgrmode*TwinOn/fSensarea);
16319 if (fAvSolidangleOn) printf(
" --> : %-g cm^-2 sr^-1",bkgrmode*TwinOn/(fSensarea*fAvSolidangleOn));
16321 printf(
" Expected background particle flux : %-g cm^-2 s^-1",bkgrmode/fSensarea);
16322 if (fAvSolidangleOn) printf(
" --> Intensity : %-g cm^-2 s^-1 sr^-1",bkgrmode/(fSensarea*fAvSolidangleOn));
16326 return hSigRatePDF;
16350 if (Non<=0 || Noff<=0)
16352 printf(
" \n *%-s::GetBurstLiMaSignificance* \n",ClassName());
16353 if (Non<=0) printf(
" === No on source data available === \n");
16354 if (Noff<=0) printf(
" === No off source data available === \n");
16362 if (fAvSolidangleOff) Ra=fAvSolidangleOn/fAvSolidangleOff;
16371 printf(
"\n *%-s::GetBurstLiMaSignificance* The Li-Ma signal significance is : %-g \n",ClassName(),sigma);
16379 if (fRecoangle && fSumsigmas && fDatype==2)
16381 printf(
" === Warning: Variable angular cones based on actual track reconstruction sigmas were used. \n");
16382 printf(
" Large variations between the on-source and off-source background counts may be present. \n");
16462 TString text=
"none";
16463 if (type==
"time") text=
"arrival time";
16464 if (type==
"BBtime") text=
"Bayesian Block event rate";
16465 if (type==
"BBrat") text=
"normalized On-source/Off-source Bayesian Block event rate ratio";
16466 if (type==
"angle") text=
"opening angle";
16467 if (type==
"cosa") text=
"cos(opening angle)";
16468 if (type==
"dt") text=
"arrival time interval";
16472 text.ReplaceAll(
"arrival time",
"redshift corrected arrival time");
16473 text.ReplaceAll(
"event rate",
"redshift corrected event rate");
16479 cout <<
" *" << ClassName() <<
"::GetBurstBayesianPsiStatistics* Unknown statistics type : " << type << endl;
16484 cout <<
" *" << ClassName() <<
"::GetBurstBayesianPsiStatistics* Analysis of " << text <<
" statistics" << endl;
16489 if (fTunits==1) tu=
"hours";
16490 if (fTunits==2) tu=
"sec";
16491 if (fTunits==3) tu=
"ns";
16492 if (fTunits==4) tu=
"ps";
16498 Double_t psitot=-1, psibkg=-1, psirat=-1;
16500 Float_t psimintot=-1, psimaxtot=-1, psifractot=0;
16501 Float_t psiminbkg=-1, psimaxbkg=-1, psifracbkg=0;
16502 Float_t psiminrat=-1, psimaxrat=-1, psifracrat=0;
16503 Double_t nrxtot=-1, nrxbkg=-1, nrxrat=-1;
16504 Double_t pvaluetot=-1, pvaluebkg=-1, pvaluerat=-1;
16529 if (!tot) printf(
" === No on source data available === \n");
16530 if (!bkg) printf(
" === No off source data available === \n");
16532 if (!tot && !bkg)
return;
16534 if (tot) psitot=math.
PsiValue(tot,0,0,freq);
16535 if (bkg) psibkg=math.
PsiValue(bkg,0,0,freq);
16536 psidif=psitot-psibkg;
16542 if (psitot<psimintot) psimintot=psitot;
16544 if (psimaxtot>psimintot) psifractot=(psimaxtot-psitot)/(psimaxtot-psimintot);
16549 if (psibkg<psiminbkg) psiminbkg=psibkg;
16551 if (psimaxbkg>psiminbkg) psifracbkg=(psimaxbkg-psibkg)/(psimaxbkg-psiminbkg);
16570 rtot=(TH1F*)hPsiOn->Clone();
16577 if (!zcor) rtot=
new TH1F(
"hPsiOnt",
"Psi distr. for bkg hypothesis of on-source arrival time data",100,psimintot-1.,psimaxtot+1.);
16578 if (zcor) rtot=
new TH1F(
"hPsiOnZt",
"Psi distr. for bkg hypothesis of redshift corrected on-source arrival time data",100,psimintot-1.,psimaxtot+1.);
16584 rbkg=(TH1F*)hPsiOff->Clone();
16591 if (!zcor) rbkg=
new TH1F(
"hPsiOfft",
"Psi distr. for bkg hypothesis of off-source arrival time data",100,psiminbkg-1.,psimaxbkg+1.);
16592 if (zcor) rbkg=
new TH1F(
"hPsiOffZt",
"Psi distr. for bkg hypothesis of redshift corrected off-source arrival time data",100,psiminbkg-1.,psimaxbkg+1.);
16598 pvaluetot=math.
PsiPvalue(-1,nr,tot,0,0,freq,0,rtot,ncut,&nrxtot);
16603 pvaluebkg=math.
PsiPvalue(-1,nr,bkg,0,0,freq,0,rbkg,ncut,&nrxbkg);
16607 cout <<
" The following randomised Psi histograms have been generated :" << endl;
16608 if (rtot) cout <<
" ... " << rtot->GetName() <<
" : " << rtot->GetTitle() << endl;
16609 if (rbkg) cout <<
" ... " << rbkg->GetName() <<
" : " << rbkg->GetTitle() << endl;
16616 if (type==
"BBtime")
16629 if (!tot) printf(
" === No on source data available === \n");
16630 if (!bkg) printf(
" === No off source data available === \n");
16632 if (!tot && !bkg)
return;
16637 for (Int_t i=1; i<=tot->GetNbinsX(); i++)
16639 tot->AddBinContent(i,1);
16645 for (Int_t i=1; i<=bkg->GetNbinsX(); i++)
16647 bkg->AddBinContent(i,1);
16651 if (tot) psitot=math.
PsiValue(tot,0,0,freq);
16652 if (bkg) psibkg=math.
PsiValue(bkg,0,0,freq);
16653 psidif=psitot-psibkg;
16659 if (psitot<psimintot) psimintot=psitot;
16661 if (psimaxtot>psimintot) psifractot=(psimaxtot-psitot)/(psimaxtot-psimintot);
16666 if (psibkg<psiminbkg) psiminbkg=psibkg;
16668 if (psimaxbkg>psiminbkg) psifracbkg=(psimaxbkg-psibkg)/(psimaxbkg-psiminbkg);
16682 hPsiOff=(TH1F*)
fBurstHistos.FindObject(
"hPsiOffBBzt");
16687 rtot=(TH1F*)hPsiOn->Clone();
16694 if (!zcor) rtot=
new TH1F(
"hPsiOnBBt",
"Psi distr. for bkg hypothesis of on-source event rate Bayesian Block data",100,psimintot-1.,psimaxtot+1.);
16695 if (zcor) rtot=
new TH1F(
"hPsiOnBBzt",
"Psi distr. for bkg hypothesis of redshift corrected on-source event rate Bayesian Block data",100,psimintot-1.,psimaxtot+1.);
16701 rbkg=(TH1F*)hPsiOff->Clone();
16708 if (!zcor) rbkg=
new TH1F(
"hPsiOffBBt",
"Psi distr. for bkg hypothesis of off-source event rate Bayesian Block data",100,psiminbkg-1.,psimaxbkg+1.);
16709 if (zcor) rbkg=
new TH1F(
"hPsiOffBBzt",
"Psi distr. for bkg hypothesis of redshift corrected off-source event rate Bayesian Block data",100,psiminbkg-1.,psimaxbkg+1.);
16715 pvaluetot=math.
PsiPvalue(-1,nr,tot,0,0,freq,0,rtot,ncut,&nrxtot);
16720 pvaluebkg=math.
PsiPvalue(-1,nr,bkg,0,0,freq,0,rbkg,ncut,&nrxbkg);
16724 cout <<
" The following randomised Psi histograms have been generated :" << endl;
16725 if (rtot) cout <<
" ... " << rtot->GetName() <<
" : " << rtot->GetTitle() << endl;
16726 if (rbkg) cout <<
" ... " << rbkg->GetName() <<
" : " << rbkg->GetTitle() << endl;
16732 for (Int_t i=1; i<=tot->GetNbinsX(); i++)
16734 tot->AddBinContent(i,-1);
16739 for (Int_t i=1; i<=bkg->GetNbinsX(); i++)
16741 bkg->AddBinContent(i,-1);
16762 printf(
" === No data available === \n");
16766 psirat=math.
PsiValue(rat,0,0,freq);
16767 psidif=psitot-psibkg;
16771 if (psirat<psiminrat) psiminrat=psirat;
16773 if (psimaxrat>psiminrat) psifracrat=(psimaxrat-psirat)/(psimaxrat-psiminrat);
16780 hPsiRat=(TH1F*)
fBurstHistos.FindObject(
"hPsiRatUBBt");
16784 hPsiRat=(TH1F*)
fBurstHistos.FindObject(
"hPsiRatUBBzt");
16789 rrat=(TH1F*)hPsiRat->Clone();
16794 if (!zcor) rrat=
new TH1F(
"hPsiRatUBBt",
"Psi distr. for bkg hypothesis of on/off source event rate Bayesian Block data",100,psiminrat-1.,psimaxrat+1.);
16795 if (zcor) rrat=
new TH1F(
"hPsiRatUBBzt",
"Psi distr. for bkg hypothesis of redshift corrected on/off source event rate Bayesian Block data",100,psiminrat-1.,psimaxrat+1.);
16798 pvaluerat=math.
PsiPvalue(-1,nr,rat,0,0,freq,0,rrat,ncut,&nrxrat);
16801 cout <<
" The following randomised Psi histograms have been generated :" << endl;
16802 if (rrat) cout <<
" ... " << rrat->GetName() <<
" : " << rrat->GetTitle() << endl;
16814 if (!tot) printf(
" === No on source data available === \n");
16815 if (!bkg) printf(
" === No off source data available === \n");
16817 if (!tot && !bkg)
return;
16819 TF1 pdfa(
"pdfa",
"sin(x*acos(-1.)/180.)");
16820 if (tot) psitot=math.
PsiValue(tot,0,&pdfa,freq);
16821 if (bkg) psibkg=math.
PsiValue(bkg,0,&pdfa,freq);
16822 psidif=psitot-psibkg;
16828 if (psitot<psimintot) psimintot=psitot;
16830 if (psimaxtot>psimintot) psifractot=(psimaxtot-psitot)/(psimaxtot-psimintot);
16835 if (psibkg<psiminbkg) psiminbkg=psibkg;
16837 if (psimaxbkg>psiminbkg) psifracbkg=(psimaxbkg-psibkg)/(psimaxbkg-psiminbkg);
16848 rtot=(TH1F*)hPsiOn->Clone();
16853 if (tot) rtot=
new TH1F(
"hPsiOna",
"Psi distr. for bkg hypothesis of on-source opening angle data",100,psimintot-1.,psimaxtot+1.);
16858 rbkg=(TH1F*)hPsiOff->Clone();
16863 if (bkg) rbkg=
new TH1F(
"hPsiOffa",
"Psi distr. for bkg hypothesis of off-source opening angle data",100,psiminbkg-1.,psimaxbkg+1.);
16868 pvaluetot=math.
PsiPvalue(-1,nr,tot,0,&pdfa,freq,0,rtot,ncut,&nrxtot);
16873 pvaluebkg=math.
PsiPvalue(-1,nr,bkg,0,&pdfa,freq,0,rbkg,ncut,&nrxbkg);
16877 cout <<
" The following randomised Psi histograms have been generated :" << endl;
16878 if (rtot) cout <<
" ... " << rtot->GetName() <<
" : " << rtot->GetTitle() << endl;
16879 if (rbkg) cout <<
" ... " << rbkg->GetName() <<
" : " << rbkg->GetTitle() << endl;
16891 if (!tot) printf(
" === No on source data available === \n");
16892 if (!bkg) printf(
" === No off source data available === \n");
16894 if (!tot && !bkg)
return;
16896 if (tot) psitot=math.
PsiValue(tot,0,0,freq);
16897 if (bkg) psibkg=math.
PsiValue(bkg,0,0,freq);
16898 psidif=psitot-psibkg;
16904 if (psitot<psimintot) psimintot=psitot;
16906 if (psimaxtot>psimintot) psifractot=(psimaxtot-psitot)/(psimaxtot-psimintot);
16911 if (psibkg<psiminbkg) psiminbkg=psibkg;
16913 if (psimaxbkg>psiminbkg) psifracbkg=(psimaxbkg-psibkg)/(psimaxbkg-psiminbkg);
16920 hPsiOff=(TH1F*)
fBurstHistos.FindObject(
"hPsiOffCosa");
16924 rtot=(TH1F*)hPsiOn->Clone();
16929 if (tot) rtot=
new TH1F(
"hPsiOnCosa",
"Psi distr. for bkg hypothesis of on-source cos(opening angle) data",100,psimintot-1.,psimaxtot+1.);
16934 rbkg=(TH1F*)hPsiOff->Clone();
16939 if (bkg) rbkg=
new TH1F(
"hPsiOffCosa",
"Psi distr. for bkg hypothesis of off-source cos(opening angle) data",100,psiminbkg-1.,psimaxbkg+1.);
16944 pvaluetot=math.
PsiPvalue(-1,nr,tot,0,0,freq,0,rtot,ncut,&nrxtot);
16949 pvaluebkg=math.
PsiPvalue(-1,nr,bkg,0,0,freq,0,rbkg,ncut,&nrxbkg);
16953 cout <<
" The following randomised Psi histograms have been generated :" << endl;
16954 if (rtot) cout <<
" ... " << rtot->GetName() <<
" : " << rtot->GetTitle() << endl;
16955 if (rbkg) cout <<
" ... " << rbkg->GetName() <<
" : " << rbkg->GetTitle() << endl;
16971 Int_t ntot=hOndt.GetEntries();
16972 Int_t nbkg=hOffdt.GetEntries();
16974 if (!ntot) printf(
" === No on source data available === \n");
16975 if (!nbkg) printf(
" === No off source data available === \n");
16977 if (!ntot && !nbkg)
return;
16979 if (ntot) psitot=math.
PsiValue(&hOndt,0,&pdfdttot,freq);
16980 if (nbkg) psibkg=math.
PsiValue(&hOffdt,0,&pdfdtbkg,freq);
16981 psidif=psitot-psibkg;
16985 psimintot=math.
PsiExtreme(&hOndt,0,&pdfdttot,-2);
16986 if (psitot<psimintot) psimintot=psitot;
16987 psimaxtot=math.
PsiExtreme(&hOndt,0,&pdfdttot,-1);
16988 if (psimaxtot>psimintot) psifractot=(psimaxtot-psitot)/(psimaxtot-psimintot);
16992 psiminbkg=math.
PsiExtreme(&hOffdt,0,&pdfdtbkg,-2);
16993 if (psibkg<psiminbkg) psiminbkg=psibkg;
16994 psimaxbkg=math.
PsiExtreme(&hOffdt,0,&pdfdtbkg,-1);
16995 if (psimaxbkg>psiminbkg) psifracbkg=(psimaxbkg-psibkg)/(psimaxbkg-psiminbkg);
17005 nametot.Form(
"hPsiOndt%-i",ndt);
17006 namebkg.Form(
"hPsiOffdt%-i",ndt);
17010 nametot.Form(
"hPsiOnZdt%-i",ndt);
17011 namebkg.Form(
"hPsiOffZdt%-i",ndt);
17015 TH1F* hPsiOff=(TH1F*)
fBurstHistos.FindObject(namebkg);
17020 rtot=(TH1F*)hPsiOn->Clone();
17025 if (!zcor) title.Form(
"Psi distr. for bkg hypothesis of on-source dt data for n=%-i",ndt);
17026 if (zcor) title.Form(
"Psi distr. for bkg hypothesis of redshift corrected on-source dt data for n=%-i",ndt);
17027 if (ntot) rtot=
new TH1F(nametot,title,100,psimintot-1.,psimaxtot+1.);
17032 rbkg=(TH1F*)hPsiOff->Clone();
17037 if (!zcor) title.Form(
"Psi distr. for bkg hypothesis of off-source dt data for n=%-i",ndt);
17038 if (zcor) title.Form(
"Psi distr. for bkg hypothesis of redshift corrected off-source dt data for n=%-i",ndt);
17039 if (nbkg) rbkg=
new TH1F(namebkg,title,100,psiminbkg-1.,psimaxbkg+1.);
17044 pvaluetot=math.
PsiPvalue(-1,nr,&hOndt,0,&pdfdttot,freq,0,rtot,ncut,&nrxtot);
17049 pvaluebkg=math.
PsiPvalue(-1,nr,&hOffdt,0,&pdfdtbkg,freq,0,rbkg,ncut,&nrxbkg);
17053 cout <<
" The following randomised Psi histograms have been (re)generated :" << endl;
17054 if (rtot) cout <<
" ... " << rtot->GetName() <<
" : " << rtot->GetTitle() << endl;
17055 if (rbkg) cout <<
" ... " << rbkg->GetName() <<
" : " << rbkg->GetTitle() << endl;
17060 cout <<
" *** Observed Psi values (in dB) for the hypothesis of no burst signal ***" << endl;
17061 if (psitot>=0) cout <<
" For the on-source stacked patches : psi = " << psitot << endl;
17062 if (psibkg>=0) cout <<
" For the off-source stacked patches : psi = " << psibkg << endl;
17063 if (psitot>=0 && psibkg>=0) cout <<
" --> Difference between observed on-source and off-source psi values : " << psidif << endl;
17064 if (psirat>=0) cout <<
" For the on/off source event rate Bayesian Blocks : psi = " << psirat << endl;
17065 cout <<
" *** Extreme Psi values for the case of pure background ***" << endl;
17066 if (psimintot>=0) cout <<
" For on-source psimin : " << psimintot <<
" psimax : " << psimaxtot;
17067 if (psifractot>0) printf(
" (psimax-psi)/range : %-g",psifractot);
17068 if (psimintot>=0 || psifractot>0) printf(
"\n");
17069 if (psiminbkg>=0) cout <<
" For off-source psimin : " << psiminbkg <<
" psimax : " << psimaxbkg;
17070 if (psifracbkg>0) printf(
" (psimax-psi)/range : %-g",psifracbkg);
17071 if (psiminbkg>=0 || psifracbkg>0) printf(
"\n");
17072 if (psiminrat>=0) cout <<
" For on/off source event rate Bayesian Blocks psimin : " << psiminrat <<
" psimax : " << psimaxrat;
17073 if (psifracrat>0) printf(
" (psimax-psi)/range : %-g",psifracrat);
17074 if (psiminrat>=0 || psifracrat>0) printf(
"\n");
17078 cout <<
" *** P-values of the observed on-source and off-source psi values ***" << endl;
17079 if (nrxtot>0) cout <<
" For the on-source stacked patches : P-value = " << pvaluetot <<
" Used number of randomisations : " << nrxtot << endl;
17080 if (nrxbkg>0) cout <<
" For the off-source stacked patches : P-value = " << pvaluebkg <<
" Used number of randomisations : " << nrxbkg << endl;
17081 if (nrxrat>0) cout <<
" For the on/off source event rate Bayesian Blocks : P-value = " << pvaluerat <<
" Used number of randomisations : " << nrxrat << endl;
17124 TString text=
"none";
17125 if (type==
"time") text=
"arrival time";
17126 if (type==
"BBtime") text=
"Bayesian Block event rate";
17127 if (type==
"BBrat") text=
"normalized On-source/Off-source Bayesian Block event rate ratio";
17128 if (type==
"angle") text=
"opening angle";
17129 if (type==
"cosa") text=
"cos(opening angle)";
17130 if (type==
"dt") text=
"arrival time interval";
17134 text.ReplaceAll(
"arrival time",
"redshift corrected arrival time");
17135 text.ReplaceAll(
"event rate",
"redshift corrected event rate");
17141 cout <<
" *" << ClassName() <<
"::GetBurstChi2Statistics* Unknown statistics type : " << type << endl;
17146 cout <<
" *" << ClassName() <<
"::GetBurstChi2Statistics* Analysis of " << text <<
" statistics" << endl;
17176 if (!tot) printf(
" === No on source data available === \n");
17177 if (!bkg) printf(
" === No off source data available === \n");
17179 if (!tot && !bkg)
return;
17181 if (tot) chitot=math.
Chi2Value(tot,0,0,&ndftot);
17182 if (bkg) chibkg=math.
Chi2Value(bkg,0,0,&ndfbkg);
17188 if (type==
"BBtime")
17201 if (!tot) printf(
" === No on source data available === \n");
17202 if (!bkg) printf(
" === No off source data available === \n");
17204 if (!tot && !bkg)
return;
17206 if (tot) chitot=math.
Chi2Value(tot,0,0,&ndftot);
17207 if (bkg) chibkg=math.
Chi2Value(bkg,0,0,&ndfbkg);
17226 printf(
" === No data available === \n");
17230 chirat=math.
Chi2Value(rat,0,0,&ndfrat);
17241 if (!tot) printf(
" === No on source data available === \n");
17242 if (!bkg) printf(
" === No off source data available === \n");
17244 if (!tot && !bkg)
return;
17246 TF1 pdf(
"pdf",
"sin(x*acos(-1.)/180.)");
17247 if (tot) chitot=math.
Chi2Value(tot,0,&pdf,&ndftot);
17248 if (bkg) chibkg=math.
Chi2Value(bkg,0,&pdf,&ndfbkg);
17259 if (!tot) printf(
" === No on source data available === \n");
17260 if (!bkg) printf(
" === No off source data available === \n");
17262 if (!tot && !bkg)
return;
17264 if (tot) chitot=math.
Chi2Value(tot,0,0,&ndftot);
17265 if (bkg) chibkg=math.
Chi2Value(bkg,0,0,&ndfbkg);
17280 Int_t ntot=hOndt.GetEntries();
17281 Int_t nbkg=hOffdt.GetEntries();
17283 if (!ntot) printf(
" === No on source data available === \n");
17284 if (!nbkg) printf(
" === No off source data available === \n");
17286 if (!ntot && !nbkg)
return;
17288 if (ntot) chitot=math.
Chi2Value(&hOndt,0,&pdfdttot,&ndftot);
17289 if (nbkg) chibkg=math.
Chi2Value(&hOffdt,0,&pdfdtbkg,&ndfbkg);
17293 Float_t chidif=chitot-chibkg;
17294 cout <<
" *** Observed Chi-squared values for the hypothesis of no burst signal ***" << endl;
17295 if (chitot>0) cout <<
" For the on-source stacked patches : chi2 = " << chitot <<
" ndf = " << ndftot << endl;
17296 if (chibkg>0) cout <<
" For the off-source stacked patches : chi2 = " << chibkg <<
" ndf = " << ndfbkg << endl;
17297 if (chitot>0 && chibkg>0) cout <<
" --> Difference between observed on-source and off-source chi2 values : " << chidif << endl;
17298 if (chirat>0) cout <<
" For the on/off source event rate Bayesian Blocks : chi2 = " << chirat <<
" ndf = " << ndfrat << endl;
17301 Float_t sigmatot=-1;
17303 Float_t sigmabkg=-1;
17305 Float_t sigmarat=-1;
17310 sigmatot=math.
Chi2Pvalue(chitot,ndftot,0,1);
17315 sigmabkg=math.
Chi2Pvalue(chibkg,ndfbkg,0,1);
17320 sigmarat=math.
Chi2Pvalue(chirat,ndfrat,0,1);
17323 cout <<
" *** P-values of the observed on-source and off-source chi2 values ***" << endl;
17324 if (ptot>=0) cout <<
" For the on-source stacked patches : P-value = " << ptot <<
" (" << sigmatot <<
" sigma)" << endl;
17325 if (pbkg>=0) cout <<
" For the off-source stacked patches : P-value = " << pbkg <<
" (" << sigmabkg <<
" sigma)" << endl;
17326 if (prat>=0) cout <<
" For the on/off source event rate Bayesian Blocks : P-value = " << prat <<
" (" << sigmarat <<
" sigma)" << endl;
17356 if (fTunits==1) tu=
"hours";
17357 if (fTunits==2) tu=
"sec";
17358 if (fTunits==3) tu=
"ns";
17359 if (fTunits==4) tu=
"ps";
17370 nametot.Form(
"hOndt%-i",ndt);
17371 namebkg.Form(
"hOffdt%-i",ndt);
17376 nametot.Form(
"hOnZdt%-i",ndt);
17377 namebkg.Form(
"hOffZdt%-i",ndt);
17384 Double_t deltatbin=0;
17391 nOndt=sOndt.
GetN();
17395 hOndt=
new TH1F(nametot,
"histo",10000,1,0);
17397 hOndt->SetBuffer(nOndt);
17399 for (Int_t i=1; i<=nOndt; i++)
17404 hOndt->BufferEmpty(1);
17407 deltatbin=hOndt->GetXaxis()->GetBinWidth(1);
17408 if (!zcor) title.Form(
"Time intervals between %-i consecutive events in the on-source time window;dt in %-s;Counts per %-.3g %-s",(ndt+1),tu.Data(),deltatbin,tu.Data());
17409 if (zcor) title.Form(
"Time intervals between %-i consecutive events in the redshift corrected on-source time window;dt in %-s;Counts per %-.3g %-s",(ndt+1),tu.Data(),deltatbin,tu.Data());
17410 hOndt->SetTitle(title);
17419 nOffdt=sOffdt.
GetN();
17423 hOffdt=
new TH1F(namebkg,
"histo",10000,1,0);
17425 hOffdt->SetBuffer(nOffdt);
17427 for (Int_t i=1; i<=nOffdt; i++)
17429 hOffdt->Fill(sOffdt.
GetEntry(i,1));
17432 hOffdt->BufferEmpty(1);
17435 deltatbin=hOffdt->GetXaxis()->GetBinWidth(1);
17436 if (!zcor) title.Form(
"Time intervals between %-i consecutive events in the off-source time window;dt in %-s;Counts per %-.3g %-s",(ndt+1),tu.Data(),deltatbin,tu.Data());
17437 if (zcor) title.Form(
"Time intervals between %-i consecutive events in the redshift corrected off-source time window;dt in %-s;Counts per %-.3g %-s",(ndt+1),tu.Data(),deltatbin,tu.Data());
17438 hOffdt->SetTitle(title);
17444 if (hOndt || hOffdt) cout <<
" The following arrival time interval (dt) histograms have been generated :" << endl;
17445 if (hOndt) cout <<
" ... " << hOndt->GetName() <<
" : " << hOndt->GetTitle() << endl;
17446 if (hOffdt) cout <<
" ... " << hOffdt->GetName() <<
" : " << hOffdt->GetTitle() << endl;
17448 if (hOndt) hisdtOn=*hOndt;
17449 if (hOffdt) hisdtOff=*hOffdt;
17465 Double_t twin=tmax-tmin;
17467 if (twin) rateOn=nevt/twin;
17473 if (twin) rateOff=nevt/twin;
17486 nametot.Form(
"hpdfOndt%-i",ndt);
17487 namebkg.Form(
"hpdfOffdt%-i",ndt);
17491 nametot.Form(
"hpdfOnZdt%-i",ndt);
17492 namebkg.Form(
"hpdfOffZdt%-i",ndt);
17495 TH1* hpdfOffdt=(TH1*)
fBurstHistos.FindObject(namebkg);
17497 Double_t xmaxfdt=1;
17498 Double_t deltatmax=0;
17501 deltatmax=hOndt->GetXaxis()->GetXmax();
17506 deltatmax=hOffdt->GetXaxis()->GetXmax();
17507 if (deltatmax>xmaxfdt) xmaxfdt=deltatmax;
17510 pdfdtOn.SetRange(0,xmaxfdt);
17511 pdfdtOn.SetNpx(npx);
17512 pdfdtOff.SetRange(0,xmaxfdt);
17513 pdfdtOff.SetNpx(npx);
17516 if (!hpdfOndt && nOndt)
17518 hpdfOndt=(TH1*)pdfdtOn.GetHistogram()->Clone();
17519 hpdfOndt->SetName(nametot.Data());
17520 title.Form(
"dt in %-s",tu.Data());
17521 hpdfOndt->GetXaxis()->SetTitle(title);
17525 if (!hpdfOffdt && nOffdt)
17527 hpdfOffdt=(TH1*)pdfdtOff.GetHistogram()->Clone();
17528 hpdfOffdt->SetName(namebkg.Data());
17529 hpdfOffdt->GetXaxis()->SetTitle(title);
17533 if (hpdfOndt || hpdfOffdt) cout <<
" The following arrival time interval (dt) PDFs have been generated :" << endl;
17534 if (hpdfOndt) cout <<
" ... " << hpdfOndt->GetName() <<
" : " << hpdfOndt->GetTitle() << endl;
17535 if (hpdfOffdt) cout <<
" ... " << hpdfOffdt->GetName() <<
" : " << hpdfOffdt->GetTitle() << endl;
17552 cout <<
" =============== The following " << nh <<
" histograms have been generated ===============" << endl;
17553 for (Int_t ih=0; ih<nh; ih++)
17557 cout <<
" " << hx->GetName() <<
" : " << hx->GetTitle() << endl;
17559 cout <<
" ===============================================================================" << endl;
17575 TFile fout(filename.Data(),
"RECREATE",
"NcAstrolab analysis results");
17579 for (Int_t ih=0; ih<nh; ih++)
17589 cout <<
" *" << ClassName() <<
"::WriteBurstHistograms* All generated histograms have been written to file " << filename << endl;
17604 if (gROOT->IsBatch())
17606 printf(
"\n *%-s::SkyMapPanel* GUI is not available in batch mode. \n",ClassName());
17625 GetDate(kTRUE,0,&year,&month,&day);
17627 fMapDate.Form(
"%02i-%02i-%-i",day,month,year);
17638 TString names[9]={
"User",
"IceCube",
"RNO-G",
"ARA",
"Amanda",
"WSRT",
"Astron",
"Greenwich",
"ARCA"};
17639 for (Int_t i=1; i<=9; i++)
17667 fSkyMapPanel->Connect(
"CloseWindow()",
"NcAstrolab",
this,
"MapClose()");
17670 TGCompositeFrame* frames[4]={0,0,0,0};
17671 TGLayoutHints* layouts[4]={0,0,0,0};
17674 frames[0]=
new TGCompositeFrame(
fSkyMapPanel,1,1,kHorizontalFrame|kSunkenFrame);
17675 layouts[0]=
new TGLayoutHints(kLHintsExpandX,border,border,0,0);
17680 frames[1]=
new TGCompositeFrame(
fSkyMapPanel,1,1,kHorizontalFrame|kSunkenFrame);
17681 layouts[1]=
new TGLayoutHints(kLHintsExpandX,border,border,0,0);
17686 frames[2]=
new TGCompositeFrame(
fSkyMapPanel,1,1,kHorizontalFrame|kSunkenFrame);
17687 layouts[2]=
new TGLayoutHints(kLHintsExpandX,border,border,0,0);
17691 frames[3]=
new TGCompositeFrame(
fSkyMapPanel,1,1,kHorizontalFrame|kSunkenFrame);
17692 layouts[3]=
new TGLayoutHints(kLHintsExpandX,border,border,0,0);
17697 for (Int_t i=0; i<4; i++)
17699 if (frames[i])
fSkyMapPanel->AddFrame(frames[i],layouts[i]);
17720 if (!frame)
return;
17722 TGGroupFrame* panel=
new TGGroupFrame(frame,
"Lab longitude, latitude, experiment site and detector ID",kHorizontalFrame);
17723 panel->SetTitlePos(TGGroupFrame::kCenter);
17724 frame->AddFrame(panel);
17729 fMapLabLBI[0]->Connect(
"TextChanged(const char*)",
"NcAstrolab",
this,
"MapLocl(const char*)");
17736 fMapLabLBI[1]->Connect(
"TextChanged(const char*)",
"NcAstrolab",
this,
"MapLocb(const char*)");
17742 fMapLabU->Connect(
"Selected(Int_t)",
"NcAstrolab",
this,
"MapUloc(Int_t)");
17753 fMapLabE->Connect(
"Selected(Int_t)",
"NcAstrolab",
this,
"MapExperiment(Int_t)");
17755 TString names[9]={
"User",
"IceCube",
"RNO-G",
"ARA",
"Amanda",
"WSRT",
"Astron",
"Greenwich",
"ARCA"};
17756 for (Int_t i=1; i<=9; i++)
17766 fMapLabLBI[2]=
new TGNumberEntryField(panel,-1,
fLabId,TGNumberFormat::kNESInteger);
17767 fMapLabLBI[2]->SetToolTipText(
"The (optional) detector element ID (0=global)");
17768 fMapLabLBI[2]->Connect(
"TextChanged(const char*)",
"NcAstrolab",
this,
"MapLocId(const char*)");
17773 TGTextButton* enter=
new TGTextButton(panel,
"Enter");
17774 enter->SetToolTipText(
"Enter the provided data");
17775 enter->Connect(
"Clicked()",
"NcAstrolab",
this,
"MapLocEnter()");
17776 TGLayoutHints* Lenter=
new TGLayoutHints(kLHintsCenterX,10,0,0,0);
17777 panel->AddFrame(enter,Lenter);
17814 TString s[4]={
"deg",
"dms",
"hms",
"rad"};
17826 TString s[9]={
"User",
"IceCube",
"RNO-G",
"ARA",
"Amanda",
"WSRT",
"Astron",
"Greenwich",
"ARCA"};
17856 SetNameTitle(
"User",
"Virtual Lab for general use");
17858 printf(
"\n *** Settings adopted for a virtual lab for general use. *** \n");
17875 printf(
"\n *** Lab settings adopted for the %-s location. *** \n",
fMapLabExpName.Data());
17880 for (Int_t i=0; i<6; i++)
17882 val.Form(
"%-.2f",
fAxes[i]);
17895 if (!frame)
return;
17897 TGGroupFrame* panel=
new TGGroupFrame(frame,
"Timestamp to be used for Entries, List and Map",kHorizontalFrame);
17898 panel->SetTitlePos(TGGroupFrame::kCenter);
17899 frame->AddFrame(panel);
17903 fMapTSdatetime->SetToolTipText(
"Date/Time as dd-mm-yyyy/hh:mm:ss.sss or MJD");
17904 fMapTSdatetime->Connect(
"TextChanged(const char*)",
"NcAstrolab",
this,
"MapDateTime(const char*)");
17911 fMapTStimetype->Connect(
"Selected(Int_t)",
"NcAstrolab",
this,
"MapTimeType(Int_t)");
17931 TGTextButton* labTS=
new TGTextButton(panel,
"Store as Lab TS");
17932 labTS->SetToolTipText(
"Store the current selection as Lab timestamp");
17933 labTS->Connect(
"Clicked()",
"NcAstrolab",
this,
"MapLabTS()");
17934 TGLayoutHints* LlabTS=
new TGLayoutHints(kLHintsCenterX,10,0,0,0);
17935 panel->AddFrame(labTS,LlabTS);
17957 TString s[13]={
"UTC",
"LMT",
"UT1",
"MJD",
"JD",
"TJD",
"Unix",
"GPS",
"TAI",
"TT",
"SysClock",
"Lab",
"EntryName"};
17978 fMapTS.GetMJD(mjd,sec,ns);
17979 Int_t ps=
fMapTS.GetPs();
17980 this->
SetMJD(mjd,sec,ns,ps,
"U");
17982 printf(
"\n *** Lab timestamp modified *** \n");
17993 if (!frame)
return;
17995 TGGroupFrame* panel=
new TGGroupFrame(frame,
"Local frame axes orientations w.r.t. X0=South, Y0=East, Z0=Zenith)",kHorizontalFrame);
17996 panel->SetTitlePos(TGGroupFrame::kCenter);
17997 frame->AddFrame(panel);
18001 fMapLabLframe[0]->SetToolTipText(
"Local X-axis zenith angle in deg");
18007 fMapLabLframe[1]->SetToolTipText(
"Local X-axis phi angle in deg");
18013 fMapLabLframe[2]->SetToolTipText(
"Local Y-axis zenith angle in deg");
18019 fMapLabLframe[3]->SetToolTipText(
"Local Y-axis phi angle in deg");
18025 fMapLabLframe[4]->SetToolTipText(
"Local Z-axis zenith angle in deg");
18031 fMapLabLframe[5]->SetToolTipText(
"Local Z-axis phi angle in deg");
18038 TGTextButton* enter=
new TGTextButton(panel,
"Enter");
18039 enter->SetToolTipText(
"Enter the provided data");
18040 enter->Connect(
"Clicked()",
"NcAstrolab",
this,
"MapLabLframeEnter()");
18041 TGLayoutHints* Lenter=
new TGLayoutHints(kLHintsLeft,10,0,0,0);
18042 panel->AddFrame(enter,Lenter);
18059 printf(
"\n *** Local frame will NOT be changed for experiment site %-s *** \n",
fMapLabExpName.Data());
18063 for (Int_t i=0; i<6; i++)
18065 val.Form(
"%-.2f",
fAxes[i]);
18071 for (Int_t i=0; i<6; i++)
18087 if (!frame)
return;
18089 TGGroupFrame* panel=
new TGGroupFrame(frame,
"Informative output",kHorizontalFrame);
18090 panel->SetTitlePos(TGGroupFrame::kCenter);
18091 frame->AddFrame(panel);
18094 TGComboBox* cinfo=
new TGComboBox(panel);
18095 cinfo->Connect(
"Selected(Int_t)",
"NcAstrolab",
this,
"MapCinfo(Int_t)");
18096 cinfo->AddEntry(
"Lab",1);
18097 cinfo->AddEntry(
"TSbox",2);
18098 cinfo->AddEntry(
"Entry",3);
18099 cinfo->AddEntry(
"Nstore",4);
18100 cinfo->Resize(60,20);
18101 panel->AddFrame(cinfo);
18102 cinfo->Select(1,kTRUE);
18106 TGComboBox* tinfo=
new TGComboBox(panel);
18107 tinfo->Connect(
"Selected(Int_t)",
"NcAstrolab",
this,
"MapTinfo(Int_t)");
18108 tinfo->AddEntry(
"LAT/LAST",1);
18109 tinfo->AddEntry(
"LMT/LMST",2);
18110 tinfo->AddEntry(
"Julian",3);
18111 tinfo->AddEntry(
"UT1",4);
18112 tinfo->AddEntry(
"UTC",5);
18113 tinfo->AddEntry(
"TAI",6);
18114 tinfo->AddEntry(
"GPS",7);
18115 tinfo->AddEntry(
"TT",8);
18116 tinfo->AddEntry(
"Unix",9);
18117 tinfo->Resize(85,20);
18118 panel->AddFrame(tinfo);
18119 tinfo->Select(1,kTRUE);
18123 TGComboBox* uinfo=
new TGComboBox(panel);
18124 uinfo->Connect(
"Selected(Int_t)",
"NcAstrolab",
this,
"MapUinfo(Int_t)");
18125 uinfo->AddEntry(
"deg",1);
18126 uinfo->AddEntry(
"dms",2);
18127 uinfo->AddEntry(
"hms",3);
18128 uinfo->AddEntry(
"rad",4);
18129 uinfo->Resize(50,20);
18130 panel->AddFrame(uinfo);
18131 uinfo->Select(1,kTRUE);
18135 TGTextEntry* mapiname=
new TGTextEntry(panel,
"");
18136 mapiname->SetToolTipText(
"Stored entry name for info");
18137 mapiname->Connect(
"TextChanged(const char*)",
"NcAstrolab",
this,
"MapIname(const char*)");
18138 mapiname->SetAlignment(kTextRight);
18139 mapiname->Resize(100,20);
18140 panel->AddFrame(mapiname);
18143 TGTextButton* info=
new TGTextButton(panel,
"Info");
18144 info->Connect(
"Clicked()",
"NcAstrolab",
this,
"MapInfo()");
18145 info->SetToolTipText(
"Provide info on the specified item");
18146 TGLayoutHints* Linfo=
new TGLayoutHints(kLHintsLeft,10,0,0,0);
18147 panel->AddFrame(info,Linfo);
18158 TString s[4]={
"Lab",
"TS box",
"EntryName",
"Nstore"};
18188 TString s[4]={
"deg",
"dms",
"hms",
"rad"};
18213 TString modes[9]={
"LAT/LAST",
"LMT/LMST",
"Julian",
"UT1",
"UTC",
"TAI",
"GPS",
"TT",
"Unix"};
18220 Int_t ntot=nref+nmeas;
18221 printf(
"\n *** Info about the stored entries *** \n");
18222 printf(
" Ntotal=%-i Nrefs=%-i Nmeas=%-i \n",ntot,nref,nmeas);
18228 printf(
"\n *** Info about the current Lab settings *** \n");
18233 printf(
"\n *** Info for the Lab timestamp *** \n");
18241 printf(
" %-s \n",datetime.Data());
18247 printf(
"\n *** Info for the timestamp in the user selection box *** \n");
18258 printf(
" Unix time : %-.12f \n",
fMapTS.GetUnixTime());
18263 printf(
" %-s \n",datetime.Data());
18275 printf(
"\n *** Info for entry : %-s *** \n",
fMapIname.Data());
18292 printf(
" Unix time : %-.12f \n",tx2.
GetUnixTime());
18297 printf(
" %-s \n",datetime.Data());
18304 printf(
"\n *** No entry found with name : %-s *** \n",
fMapIname.Data());
18317 if (!frame)
return;
18319 TGGroupFrame* panel=
new TGGroupFrame(frame,
"Entries in (a,b) coordinates",kHorizontalFrame);
18320 panel->SetTitlePos(TGGroupFrame::kCenter);
18321 frame->AddFrame(panel);
18324 TGNumberEntryField* ea=
new TGNumberEntryField(panel,-1,0);
18325 ea->SetToolTipText(
"Angle a");
18326 ea->Connect(
"TextChanged(const char*)",
"NcAstrolab",
this,
"MapEa(const char*)");
18327 ea->Resize(100,20);
18328 panel->AddFrame(ea);
18332 TGComboBox* ua=
new TGComboBox(panel);
18333 ua->Connect(
"Selected(Int_t)",
"NcAstrolab",
this,
"MapUa(Int_t)");
18334 ua->AddEntry(
"deg",1);
18335 ua->AddEntry(
"dms",2);
18336 ua->AddEntry(
"hms",3);
18337 ua->AddEntry(
"rad",4);
18338 ua->AddEntry(
"hrs",5);
18340 panel->AddFrame(ua);
18341 ua->Select(1,kTRUE);
18345 TGNumberEntryField* eb=
new TGNumberEntryField(panel,-1,0);
18346 eb->SetToolTipText(
"Angle b");
18347 eb->Connect(
"TextChanged(const char*)",
"NcAstrolab",
this,
"MapEb(const char*)");
18348 eb->Resize(100,20);
18349 panel->AddFrame(eb);
18353 TGComboBox* ub=
new TGComboBox(panel);
18354 ub->Connect(
"Selected(Int_t)",
"NcAstrolab",
this,
"MapUb(Int_t)");
18355 ub->AddEntry(
"deg",1);
18356 ub->AddEntry(
"dms",2);
18357 ub->AddEntry(
"hms",3);
18358 ub->AddEntry(
"rad",4);
18359 ub->AddEntry(
"hrs",5);
18361 panel->AddFrame(ub);
18362 ub->Select(1,kTRUE);
18365 TGComboBox* ecoords=
new TGComboBox(panel);
18366 ecoords->Connect(
"Selected(Int_t)",
"NcAstrolab",
this,
"MapEcoord(Int_t)");
18367 ecoords->AddEntry(
"(ra,dec) (J2000)",1);
18368 ecoords->AddEntry(
"(ra,dec) (Mean)",2);
18369 ecoords->AddEntry(
"(ra,dec) (True)",3);
18370 ecoords->AddEntry(
"(ra,dec) (B1950)",4);
18371 ecoords->AddEntry(
"Galactic (l,b)",5);
18372 ecoords->AddEntry(
"Ecliptic (l,b)",6);
18373 ecoords->AddEntry(
"Horizon (azi,zen)",7);
18374 ecoords->AddEntry(
"ICR (l,b)",8);
18375 ecoords->AddEntry(
"Local (theta,phi)",9);
18376 ecoords->AddEntry(
"pdir (theta,phi)",10);
18377 ecoords->Resize(125,20);
18378 panel->AddFrame(ecoords);
18379 ecoords->Select(1,kTRUE);
18383 TGComboBox* etypes=
new TGComboBox(panel);
18384 etypes->Connect(
"Selected(Int_t)",
"NcAstrolab",
this,
"MapEtype(Int_t)");
18385 etypes->AddEntry(
"Meas",1);
18386 etypes->AddEntry(
"Ref",2);
18387 etypes->Resize(55,20);
18388 panel->AddFrame(etypes);
18389 etypes->Select(1,kTRUE);
18393 TGTextEntry* ename=
new TGTextEntry(panel,
"");
18394 ename->SetToolTipText(
"The (optional) entry name");
18395 ename->Connect(
"TextChanged(const char*)",
"NcAstrolab",
this,
"MapEname(const char*)");
18396 ename->SetAlignment(kTextRight);
18397 ename->Resize(100,20);
18398 panel->AddFrame(ename);
18401 TGTextButton* enter=
new TGTextButton(panel,
"Enter");
18402 enter->SetToolTipText(
"Enter the provided entry");
18403 enter->Connect(
"Clicked()",
"NcAstrolab",
this,
"MapEnter()");
18404 TGLayoutHints* Lenter=
new TGLayoutHints(kLHintsCenterX,10,0,0,0);
18405 panel->AddFrame(enter,Lenter);
18408 TGTextButton* remove=
new TGTextButton(panel,
"Remove");
18409 remove->SetToolTipText(
"Remove the entry specified by type and name pattern (name=* means all)");
18410 remove->Connect(
"Clicked()",
"NcAstrolab",
this,
"MapRemove()");
18411 TGLayoutHints* Lremove=
new TGLayoutHints(kLHintsCenterX,10,0,0,0);
18412 panel->AddFrame(remove,Lremove);
18435 TString s[5]={
"deg",
"dms",
"hms",
"rad",
"hrs"};
18459 TString s[5]={
"deg",
"dms",
"hms",
"rad",
"hrs"};
18471 TString system[10]={
"equ",
"equ",
"equ",
"equ",
"gal",
"ecl",
"hor",
"icr",
"loc",
"pdir"};
18472 TString mode[4]={
"J",
"M",
"T",
"B"};
18514 SetSignal(1,
fMapEa,
fMapEua,
fMapEb,
fMapEub,
fMapEcoord,&
fMapTS,-1,
fMapEmode,
fMapEname,
fMapEtype);
18516 printf(
"\n *** Specified entry stored *** \n");
18534 printf(
"\n *** Number of specified entries have been removed : %-i *** \n",nrem);
18545 if (!frame)
return;
18547 TGGroupFrame* panel=
new TGGroupFrame(frame,
"Map/List options for the (a,b) coordinates",kHorizontalFrame);
18548 panel->SetTitlePos(TGGroupFrame::kCenter);
18549 frame->AddFrame(panel);
18552 TGVerticalFrame* render=
new TGVerticalFrame(panel,1,1);
18553 panel->AddFrame(render);
18556 TGGroupFrame* coordsys=
new TGGroupFrame(render,
"Coordinate system",kHorizontalFrame);
18557 coordsys->SetTitlePos(TGGroupFrame::kCenter);
18558 render->AddFrame(coordsys);
18559 TGComboBox* dcoords=
new TGComboBox(coordsys);
18560 dcoords->Connect(
"Selected(Int_t)",
"NcAstrolab",
this,
"MapDcoord(Int_t)");
18561 dcoords->AddEntry(
"Equatorial (J2000)",1);
18562 dcoords->AddEntry(
"Equatorial (Mean)",2);
18563 dcoords->AddEntry(
"Equatorial (True)",3);
18564 dcoords->AddEntry(
"Equatorial (B1950)",4);
18565 dcoords->AddEntry(
"Galactic",5);
18566 dcoords->AddEntry(
"Ecliptic",6);
18567 dcoords->AddEntry(
"Horizon",7);
18568 dcoords->AddEntry(
"ICR",8);
18569 dcoords->AddEntry(
"Local",9);
18570 dcoords->Resize(140,20);
18571 coordsys->AddFrame(dcoords);
18572 dcoords->Select(1,kTRUE);
18576 TGGroupFrame* mapview=
new TGGroupFrame(render,
"Map representation",kHorizontalFrame);
18577 mapview->SetTitlePos(TGGroupFrame::kCenter);
18578 render->AddFrame(mapview);
18579 TGComboBox* projs=
new TGComboBox(mapview);
18580 projs->Connect(
"Selected(Int_t)",
"NcAstrolab",
this,
"MapProj(Int_t)");
18581 projs->AddEntry(
"Hammer projection",1);
18582 projs->AddEntry(
"Aitoff projection",2);
18583 projs->AddEntry(
"Mercator projection",3);
18584 projs->AddEntry(
"sin(b) vs. a",4);
18585 projs->AddEntry(
"b vs. a",5);
18586 projs->AddEntry(
"b vs. UT (0-24 hrs)",6);
18587 projs->AddEntry(
"b vs. LT (0-24 hrs)",7);
18588 projs->AddEntry(
"b vs. GST (0-24 hrs)",8);
18589 projs->AddEntry(
"b vs. LST (0-24 hrs)",9);
18590 projs->AddEntry(
"b vs. Day at UT",10);
18591 projs->AddEntry(
"b vs. Day at LT",11);
18592 projs->AddEntry(
"b vs. Day at GST",12);
18593 projs->AddEntry(
"b vs. Day at LST",13);
18594 projs->Resize(150,20);
18595 mapview->AddFrame(projs);
18596 projs->Select(1,kTRUE);
18600 TGGroupFrame* meridian=
new TGGroupFrame(render,
"Meridian ordering",kHorizontalFrame);
18601 meridian->SetTitlePos(TGGroupFrame::kCenter);
18602 render->AddFrame(meridian);
18605 TGComboBox* mermode=
new TGComboBox(meridian);
18606 mermode->Connect(
"Selected(Int_t)",
"NcAstrolab",
this,
"MapMerMode(Int_t)");
18607 mermode->AddEntry(
"Auto",1);
18608 mermode->AddEntry(
"---->",2);
18609 mermode->AddEntry(
"<----",3);
18610 mermode->Resize(55,20);
18611 meridian->AddFrame(mermode);
18612 mermode->Select(1,kTRUE);
18616 TGNumberEntryField* merc=
new TGNumberEntryField(meridian,-1,0);
18617 merc->SetToolTipText(
"Central meridian position");
18618 merc->Connect(
"TextChanged(const char*)",
"NcAstrolab",
this,
"MapMerC(const char*)");
18619 merc->Resize(65,20);
18620 meridian->AddFrame(merc);
18624 TGComboBox* meruc=
new TGComboBox(meridian);
18625 meruc->Connect(
"Selected(Int_t)",
"NcAstrolab",
this,
"MapMerUc(Int_t)");
18626 meruc->AddEntry(
"deg",1);
18627 meruc->AddEntry(
"dms",2);
18628 meruc->AddEntry(
"hms",3);
18629 meruc->AddEntry(
"rad",4);
18630 meruc->Resize(50,20);
18631 meridian->AddFrame(meruc);
18632 meruc->Select(1,kTRUE);
18636 TGVButtonGroup* options=
new TGVButtonGroup(panel,
"Options");
18637 options->SetTitlePos(TGGroupFrame::kCenter);
18638 options->Connect(
"Clicked(Int_t)",
"NcAstrolab",
this,
"MapDoptions(Int_t)");
18639 TGCheckButton* bhist=
new TGCheckButton(options,
"Hist");
18640 bhist->SetToolTipText(
"Project data in a histogram");
18641 TGCheckButton* bclr=
new TGCheckButton(options,
"Clr");
18642 bclr->SetToolTipText(
"Clear display before drawing");
18643 TGCheckButton* bref=
new TGCheckButton(options,
"Ref");
18644 bref->SetToolTipText(
"Display reference signals");
18645 TGCheckButton* bmeas=
new TGCheckButton(options,
"Meas");
18646 bmeas->SetToolTipText(
"Display measured signals");
18647 TGCheckButton* brefts=
new TGCheckButton(options,
"RefTS");
18648 brefts->SetToolTipText(
"Display/list each reference signal using its actual recorded timestamp");
18649 panel->AddFrame(options);
18651 options->SetButton(1,kFALSE);
18652 options->SetButton(2,kTRUE);
18653 options->SetButton(3,kTRUE);
18654 options->SetButton(4,kTRUE);
18655 options->SetButton(5,kFALSE);
18663 TGNumberEntryField* nmax=
new TGNumberEntryField(options,-1,-1,TGNumberFormat::kNESInteger);
18664 nmax->SetToolTipText(
"Max. number of Drawn/Listed entries per type (-1=all)");
18665 nmax->Connect(
"TextChanged(const char*)",
"NcAstrolab",
this,
"MapNmax(const char*)");
18666 nmax->Resize(40,20);
18667 options->AddFrame(nmax);
18671 TGTextEntry* sname=
new TGTextEntry(options,
"*");
18672 sname->SetToolTipText(
"The requested entry name pattern (*=all)");
18673 sname->Connect(
"TextChanged(const char*)",
"NcAstrolab",
this,
"MapDname(const char*)");
18674 sname->SetAlignment(kTextRight);
18675 sname->Resize(100,20);
18676 TGLayoutHints* Lsname=
new TGLayoutHints(kLHintsLeft,0,0,5,0);
18677 options->AddFrame(sname,Lsname);
18681 TGNumberEntryField* ndigs=
new TGNumberEntryField(options,-1,1,TGNumberFormat::kNESInteger);
18682 ndigs->SetToolTipText(
"Number of digits for the listed entries");
18683 ndigs->Connect(
"TextChanged(const char*)",
"NcAstrolab",
this,
"MapNdigs(const char*)");
18684 ndigs->Resize(40,20);
18685 options->AddFrame(ndigs);
18689 TGVerticalFrame* others=
new TGVerticalFrame(panel,1,1);
18690 panel->AddFrame(others);
18693 TGGroupFrame* markers=
new TGGroupFrame(others,
"Marker settings",kHorizontalFrame);
18694 markers->SetTitlePos(TGGroupFrame::kCenter);
18695 others->AddFrame(markers);
18698 TGNumberEntryField* marksize=
new TGNumberEntryField(markers,-1,1);
18699 marksize->SetToolTipText(
"Marker size");
18700 marksize->Connect(
"TextChanged(const char*)",
"NcAstrolab",
this,
"MapMarkSize(const char*)");
18701 marksize->Resize(40,20);
18702 markers->AddFrame(marksize);
18706 TGComboBox* markstyle=
new TGComboBox(markers);
18707 markstyle->Connect(
"Selected(Int_t)",
"NcAstrolab",
this,
"MapMarkStyle(Int_t)");
18708 markstyle->AddEntry(
"Dot",1);
18709 markstyle->AddEntry(
"Star",2);
18710 markstyle->AddEntry(
"Square",3);
18711 markstyle->AddEntry(
"Utriangle",4);
18712 markstyle->AddEntry(
"Dtriangle",5);
18713 markstyle->AddEntry(
"Diamond",6);
18714 markstyle->AddEntry(
"Cross",7);
18715 markstyle->AddEntry(
"Ast",8);
18716 markstyle->AddEntry(
"Plus",9);
18717 markstyle->AddEntry(
"Times",10);
18718 markstyle->AddEntry(
"Circle",11);
18719 markstyle->AddEntry(
"oStar",12);
18720 markstyle->AddEntry(
"oSquare",13);
18721 markstyle->AddEntry(
"oUtriangle",14);
18722 markstyle->AddEntry(
"oDtriangle",15);
18723 markstyle->AddEntry(
"oDiamond",16);
18724 markstyle->AddEntry(
"oCross",17);
18725 markstyle->Resize(90,20);
18726 markers->AddFrame(markstyle);
18727 markstyle->Select(1,kTRUE);
18731 TGComboBox* markcolor=
new TGComboBox(markers);
18732 markcolor->Connect(
"Selected(Int_t)",
"NcAstrolab",
this,
"MapMarkColor(Int_t)");
18733 markcolor->AddEntry(
"Black",1);
18734 markcolor->AddEntry(
"Red",2);
18735 markcolor->AddEntry(
"Blue",3);
18736 markcolor->AddEntry(
"Green",4);
18737 markcolor->AddEntry(
"Yellow",5);
18738 markcolor->AddEntry(
"Magenta",6);
18739 markcolor->AddEntry(
"Cyan",7);
18740 markcolor->AddEntry(
"Orange",8);
18741 markcolor->AddEntry(
"Violet",9);
18742 markcolor->AddEntry(
"Pink",10);
18743 markcolor->AddEntry(
"Azure",11);
18744 markcolor->AddEntry(
"Spring",12);
18745 markcolor->AddEntry(
"Teal",13);
18746 markcolor->AddEntry(
"Gray",14);
18747 markcolor->AddEntry(
"White",15);
18748 markcolor->Resize(80,20);
18749 markers->AddFrame(markcolor);
18750 markcolor->Select(3,kTRUE);
18754 TGComboBox* marktype=
new TGComboBox(markers);
18755 marktype->Connect(
"Selected(Int_t)",
"NcAstrolab",
this,
"MapMarkType(Int_t)");
18756 marktype->AddEntry(
"Ref",1);
18757 marktype->AddEntry(
"Meas",2);
18758 marktype->AddEntry(
"GC",3);
18759 marktype->AddEntry(
"Grid",4);
18760 marktype->Resize(55,20);
18761 markers->AddFrame(marktype);
18762 marktype->Select(2,kTRUE);
18766 TGGroupFrame* solar=
new TGGroupFrame(others,
"Selection of Solar system reference objects",kHorizontalFrame);
18767 solar->SetTitlePos(TGGroupFrame::kCenter);
18768 others->AddFrame(solar);
18770 TGButtonGroup* Ssolar=
new TGButtonGroup(solar,0,3,5);
18771 Ssolar->SetTitlePos(TGGroupFrame::kCenter);
18772 Ssolar->Connect(
"Clicked(Int_t)",
"NcAstrolab",
this,
"MapSolar(Int_t)");
18773 TGCheckButton* bsun=
new TGCheckButton(Ssolar,
"Sun");
18774 bsun->SetToolTipText(
"Enter/Remove the Sun as a reference entry");
18775 TGCheckButton* bmoon=
new TGCheckButton(Ssolar,
"Moon");
18776 bmoon->SetToolTipText(
"Enter/Remove the Moon as a reference entry");
18777 TGCheckButton* bmercury=
new TGCheckButton(Ssolar,
"Mercury");
18778 bmercury->SetToolTipText(
"Enter/Remove Mercury as a reference entry");
18779 TGCheckButton* bvenus=
new TGCheckButton(Ssolar,
"Venus");
18780 bvenus->SetToolTipText(
"Enter/Remove Venus as a reference entry");
18781 TGCheckButton* bearth=
new TGCheckButton(Ssolar,
"Earth");
18782 bearth->SetToolTipText(
"Enter/Remove the Earth as a reference entry");
18783 TGCheckButton* bmars=
new TGCheckButton(Ssolar,
"Mars");
18784 bmars->SetToolTipText(
"Enter/Remove Mars as a reference entry");
18785 TGCheckButton* bjupiter=
new TGCheckButton(Ssolar,
"Jupiter");
18786 bjupiter->SetToolTipText(
"Enter/Remove Jupiter as a reference entry");
18787 TGCheckButton* bsaturn=
new TGCheckButton(Ssolar,
"Saturn");
18788 bsaturn->SetToolTipText(
"Enter/Remove Saturn as a reference entry");
18789 TGCheckButton* buranus=
new TGCheckButton(Ssolar,
"Uranus");
18790 buranus->SetToolTipText(
"Enter/Remove Uranus as a reference entry");
18791 TGCheckButton* bneptune=
new TGCheckButton(Ssolar,
"Neptune");
18792 bneptune->SetToolTipText(
"Enter/Remove Neptune as a reference entry");
18793 solar->AddFrame(Ssolar);
18795 Ssolar->SetButton(1,kFALSE);
18796 Ssolar->SetButton(2,kFALSE);
18797 Ssolar->SetButton(3,kFALSE);
18798 Ssolar->SetButton(4,kFALSE);
18799 Ssolar->SetButton(5,kFALSE);
18800 Ssolar->SetButton(6,kFALSE);
18801 Ssolar->SetButton(7,kFALSE);
18802 Ssolar->SetButton(8,kFALSE);
18803 Ssolar->SetButton(9,kFALSE);
18804 Ssolar->SetButton(10,kFALSE);
18805 for (Int_t i=0; i<10; i++)
18811 TGVerticalFrame* comms=
new TGVerticalFrame(solar,1,1);
18812 TGLayoutHints* Lcomms=
new TGLayoutHints(kLHintsLeft,10,0,15,0);
18813 solar->AddFrame(comms,Lcomms);
18816 TGTextButton* enter=
new TGTextButton(comms,
"Enter");
18817 enter->SetToolTipText(
"Enter the selected Solar system objects as reference signals");
18818 enter->Connect(
"Clicked()",
"NcAstrolab",
this,
"MapEnterSolar()");
18819 TGLayoutHints* Lenter=
new TGLayoutHints(kLHintsCenterX,0,0,10,15);
18820 comms->AddFrame(enter,Lenter);
18823 TGTextButton* remove=
new TGTextButton(comms,
"Remove");
18824 remove->SetToolTipText(
"Remove the selected Solar system objects from the reference signals");
18825 remove->Connect(
"Clicked()",
"NcAstrolab",
this,
"MapRemoveSolar()");
18826 TGLayoutHints* Lremove=
new TGLayoutHints(kLHintsCenterX,0,0,0,0);
18827 comms->AddFrame(remove,Lremove);
18838 TString system[9]={
"equ",
"equ",
"equ",
"equ",
"gal",
"ecl",
"hor",
"icr",
"loc"};
18839 TString mode[4]={
"J",
"M",
"T",
"B"};
18854 TString s[13]={
"ham",
"ait",
"mer",
"ang",
"cyl",
"UTh",
"LTh",
"GSTh",
"LSTh",
"UYh",
"LYh",
"GSYh",
"LSYh"};
18855 TString sh[13]={
"hamh",
"aith",
"merh",
"angh",
"cylh",
"UTh",
"LTh",
"GSTh",
"LSTh",
"UYh",
"LYh",
"GSYh",
"LSYh"};
18858 for (Int_t k=0; k<13; k++)
18873 Int_t modes[3]={0,1,-1};
18897 TString s[4]={
"deg",
"dms",
"hms",
"rad"};
18977 Int_t styles[17]={8,29,21,22,23,33,34,31,2,5,24,30,25,26,32,27,28};
18992 Int_t colors[15]={kBlack,kRed,kBlue,kGreen,kYellow,kMagenta,kCyan,kOrange,kViolet,kPink,kAzure,kSpring,kTeal,kGray,kWhite};
19038 TString names[10]={
"Sun",
"Moon",
"Mercury",
"Venus",
"Earth",
"Mars",
"Jupiter",
"Saturn",
"Uranus",
"Neptune"};
19044 for (Int_t i=0; i<10; i++)
19048 printf(
"\n *** Selected Solar system object(s) entered *** \n");
19059 TString names[10]={
"Sun",
"Moon",
"Mercury",
"Venus",
"Earth",
"Mars",
"Jupiter",
"Saturn",
"Uranus",
"Neptune"};
19067 for (Int_t i=0; i<10; i++)
19076 printf(
"\n *** Number of Solar system object that have been removed : %-i *** \n",nrem);
19087 if (!frame)
return;
19089 TGGroupFrame* panel=
new TGGroupFrame(frame,
"Commands",kVerticalFrame);
19090 panel->SetTitlePos(TGGroupFrame::kCenter);
19091 frame->AddFrame(panel);
19093 TGTextButton* list=
new TGTextButton(panel,
"List");
19094 list->Connect(
"Clicked()",
"NcAstrolab",
this,
"MapList()");
19095 list->SetToolTipText(
"List the selected entries");
19096 TGLayoutHints* Llist=
new TGLayoutHints(kLHintsCenterX,0,0,10,10);
19097 panel->AddFrame(list,Llist);
19099 TGTextButton* map=
new TGTextButton(panel,
"Map");
19100 map->Connect(
"Clicked()",
"NcAstrolab",
this,
"MapDraw()");
19101 map->SetToolTipText(
"Display the selected entries");
19102 TGLayoutHints* Lmap=
new TGLayoutHints(kLHintsCenterX,0,0,10,10);
19103 panel->AddFrame(map,Lmap);
19105 TGTextButton* close=
new TGTextButton(panel,
"Close");
19106 close->Connect(
"Clicked()",
"NcAstrolab",
this,
"MapClose()");
19107 close->SetToolTipText(
"Close this panel window");
19108 TGLayoutHints* Lclose=
new TGLayoutHints(kLHintsCenterX,0,0,10,10);
19109 panel->AddFrame(close,Lclose);
19111 TGTextButton* exit=
new TGTextButton(panel,
"Exit");
19112 exit->Connect(
"Clicked()",
"NcAstrolab",
this,
"MapExit()");
19113 exit->SetToolTipText(
"Exit this ROOT session");
19114 TGLayoutHints* Lexit=
new TGLayoutHints(kLHintsCenterX,0,0,10,10);
19115 panel->AddFrame(exit,Lexit);
19203 gApplication->Terminate(0);
19247 fMapTS.SetUnixTime(val);
19261 tx->
GetMJD(imjd,isec,ins);
19263 fMapTS.SetMJD(imjd,isec,ins,ips);
19267 printf(
"\n *** No stored entry with name %-s found --> Lab TS will be used *** \n",name.Data());
19270 fMapTS.SetMJD(imjd,isec,ins,ips);
19319 if (strlen(name)) lab->SetName(name);
Handling of 3-vectors in various reference frames.
Double_t GetX(Int_t i, TString f, TString u="rad")
void GetErrors(Double_t *e, TString f, TString u="rad") const
void SetVector(Double_t *v, TString f, TString u="rad")
virtual void Load(Nc3Vector &q)
Nc3Vector Cross(Nc3Vector &q) const
void SetErrors(Double_t *e, TString f, TString u="rad")
virtual Double_t GetOpeningAngle(Nc3Vector &q, TString u="rad")
void GetVector(Double_t *v, TString f, TString u="rad") const
Nc3Vector GetPrimed(TRotMatrix *m) const
Nc3Vector GetUnprimed(TRotMatrix *m) const
Virtual lab to provide (astro)physical parameters, treat data and relate observations with astrophysi...
void Precess(Nc3Vector &r, NcTimestamp *ts1, NcTimestamp *ts2)
void MapEname(const char *text)
void Nutate(Nc3Vector &r, NcTimestamp *ts)
void MakeBurstRecoAngresdist(TString file, TString tree, TString name1, TString name2, TString ua, TString name3, TString ud, Double_t Emin, Double_t Emax, Int_t nbe=100, Int_t nba=1000)
Double_t fMapMerC
! The GUI entered central meridian location for the Map
void ProjectMercator(Double_t l, Double_t b, Double_t &x, Double_t &y)
Double_t fMeridian
! Central meridian (in rad) for the sky display
Double_t GetTargetThickness(Double_t prob, Double_t lambda) const
virtual void LabLocationPanel(TGCompositeFrame *frame)
void GetBurstBayesianPsiStatistics(TString type, Double_t nr=-1, Int_t ncut=10, Int_t ndt=2, Bool_t zcor=kFALSE, Int_t freq=0)
Double_t GetCredibleInterval(TF1 pdf, Double_t p, Double_t &xlow, Double_t &xup, Int_t n=1000)
void SetRandomiser(Int_t iseed, Int_t cnt1=0, Int_t cnt2=0, NcTimestamp *ts=0)
Double_t GetUpperLimit(TF1 pdf, Double_t p)
TH1F GetDxHistogram(TH1 *hx, Int_t nc, Double_t dxbin=-1, Double_t dxmin=-1, Double_t dxmax=-1, Int_t mode=1, Double_t fact=1)
Double_t fMapMarkSize
! The GUI entered marker size for the Map
void HelioToGeocentric(Double_t &R, Double_t &B, Double_t &L, NcTimestamp *ts, TString Bu="deg", TString Lu="deg")
void MapMarkType(Int_t i)
Int_t GetTimeScramble(Double_t *tmin=0, Double_t *tmax=0, TF1 *frndm=0)
TRotMatrix fP
! Matrix for precession correction
void MapEb(const char *text)
void MatchSignals(NcDevice &matches, Double_t da, TString au, Double_t dt, TString tu, Int_t mode=1, Int_t i1=1, Int_t i2=0, Int_t itype=0, Int_t j1=1, Int_t j2=0, Int_t jtype=1)
TString fMapEcoord
! The GUI entered coordinate system of the entry
Int_t SetSolarSystem(TString name, NcTimestamp *ts, Int_t type=0)
Double_t GetLightTravelDistance(Double_t z, TString u="Mpc") const
Double_t fMapLabLocB
! The GUI entered Lab latitude
TString fMapEua
! The GUI entered angular units of a
void SetBurstParameter(TString name, Double_t value)
void SmearPosition(Nc3Vector &v, Double_t sigma)
Double_t GetLightTravelTime(Double_t z) const
Int_t GetPositionScramble(Double_t *dmin=0, Double_t *dmax=0, TF1 *df=0, Double_t *thmin=0, Double_t *thmax=0, TF1 *thf=0, Double_t *phimin=0, Double_t *phimax=0, TF1 *phif=0)
void ProjectCylindrical(Double_t l, Double_t b, Double_t &x, Double_t &y)
void MapExperiment(Int_t i)
void MapMarkSize(const char *text)
TRotMatrix fG
! Matrix for conversion of equatorial to galactic coordinates
TString fMapTime
! The GUI entered time
Double_t GetBurstRecoAngres(Double_t Emin=-1, Double_t Emax=-1, Double_t Amin=0, Double_t Amax=999) const
TRotMatrix fB
! The frame bias matrix for conversion of ICRS to J2000 coordinates
virtual void EntriesPanel(TGCompositeFrame *frame)
TH1 * GetBurstBayesianSignalRate(Double_t p, Double_t &rlow, Double_t &rup, Int_t n=1000)
Double_t ConvertAngle(Double_t a, TString in, TString out) const
void SetNmatrix(NcTimestamp *ts)
void ProjectHammer(Double_t l, Double_t b, Double_t &x, Double_t &y)
Double_t GetDifference(Int_t jref, TString au, Double_t &dt, TString tu, Int_t mode=1, Int_t *ia=0, Int_t *it=0)
void ListSignals(TString frame, TString mode, Int_t ndig=1, TString emode="T", Int_t nmax=10, Int_t j=-1, Int_t type=-1, NcTimestamp *ts=0, TString name="*")
Double_t GetSourceAttributes(NcSignal *s, Float_t *z=0, Float_t *T90=0)
TString fProj
! Projection which is currently in use
Double_t GetLuminosityDistance(Double_t z, TString u="Mpc") const
NcDevice * fBurstParameters
TRotMatrix fH
! Matrix for conversion of equatorial to horizontal coordinates
Int_t RemoveSignal(Int_t j, Int_t type, Int_t compress)
Double_t GetBurstBackgroundEnergy(Double_t Emin=-1, Double_t Emax=-1) const
void SetLabPosition(Nc3Vector &r)
TString fMapTimeType
! The GUI entered time type
TString GetExperiment() const
void MapDoptions(Int_t i)
Double_t GetLabTimeOffset() const
void GenBurstGCNdata(Int_t n, TString name="GRB", Bool_t scale=kFALSE)
Double_t GetComovingDistance(Double_t z, TString u="Mpc") const
TH2 * fHist[2]
! Temp. histograms for the sky display
TGComboBox * fMapTStimetype
! The GUI TS time type selection box
TString fMapCinfo
! The GUI selected info category
TArrayI * fIndices
! Storage indices of the matching reference signals
TH1 * GetBurstSigmaPosdist(TString name, TString type)
TString fMapProj
! The GUI selected projection for the Map
Bool_t fMapLabTS
! The GUI selection to use the Lab timestamp for the List/Map
void GetLocalFrame(Float_t arr[6])
void LoadInputData(Bool_t src, TString file, TString tree, Int_t date1=0, Int_t date2=0, Int_t nmax=-1, TString type="-")
void InitBurstHistograms(Int_t mode)
TString fMapLabLocU
! The GUI entered Lab location angular units
TCanvas * fCanvas
! The canvas for the skymap
void MapDname(const char *text)
void MapIname(const char *text)
NcAstrolab(const char *name="User", const char *title="Virtual Lab for general use")
Double_t GetBurstLiMaSignificance() const
virtual void SkyMapPanel()
Int_t GetSignalIndex(TString name, Int_t type=0)
TGTextEntry * fMapTSdatetime
! The GUI TS date/time specification
void MapLocb(const char *text)
TRotMatrix fN
! Matrix for nutation correction
void ListBurstParameters() const
TGNumberEntryField * fMapLabLBI[3]
! The GUI number entries for the Lab location specification
void MapEa(const char *text)
Double_t GetMeanFreePath(Double_t sigma, Double_t rho, Int_t mode) const
NcPosition GetLabPosition() const
void SetPositionScramble(Int_t mode, Double_t dmin, Double_t dmax, TF1 *df=0, Double_t thmin=0, Double_t thmax=0, TF1 *thf=0, Double_t phimin=0, Double_t phimax=0, TF1 *phif=0)
void InitDataNames(Int_t dir, TString frame, TString mode="J")
Double_t GetHubbleParameter(Double_t z, TString u="Mpc") const
Int_t fMapMarkStyle
! The GUI selected marker style for the Map
Int_t SetSourceAttributes(NcSignal *s, Double_t sigmapos, TString u, Double_t z=-999, Double_t T90=-999)
Double_t GetMaxDt() const
void SetTimeScramble(Int_t mode, Double_t tmin, Double_t tmax, TF1 *frndm=0)
TString fMapDcoord
! The GUI selected coordinate system for the Map/List
void SetPmatrix(NcTimestamp *ts)
Double_t fMapEb
! The GUI entered b coordinate of an entry
TString fMapDate
! The GUI entered date
TF1 GetSignalRatePDF(Double_t Non, Double_t Ton, Double_t Noff, Double_t Toff, Double_t Ra=1, Double_t Re=1, Double_t smax=-1, Double_t bmax=-1, Double_t prec=709)
Double_t fMapLabLocL
! The GUI entered Lab longitude
Int_t fMapMarkType
! The GUI selected entry type to apply the marker attributes on
void MapMarkStyle(Int_t i)
void Project(Double_t l, Double_t b, TString proj, Double_t &x, Double_t &y)
void MapDateTime(const char *text)
TString fMapDmode
! The GUI selected coordinate system mode for the Map/List
Double_t GetShieldingThickness(Double_t prob, Double_t lambda) const
TString fMapEub
! The GUI entered angular units of b
void MapMarkColor(Int_t i)
Double_t GetBackgroundRateProb(Double_t *vars, Double_t *pars)
Double_t GetInteractionProbability(Double_t x, Double_t lambda) const
void MakeBurstSigmaPosdist(TString file, TString tree, TString name, TString u, Int_t nb=900, Float_t xmin=0, Float_t xmax=90)
void MapNdigs(const char *text)
Int_t GetNRefSignals(Int_t mode=0) const
TH1F GetDifHistogram(TH1 *hin, Int_t mode, TString s="", TF1 *f=0) const
Int_t fMapEtype
! The GUI entered entry type
TGComboBox * fMapLabU
! The GUI Lab location angular unit selection box
void SetDataNames(TString obsname, TString varname, TString units="1", TString func="none")
TString fMapDname
! The GUI entered name pattern for entries to be shown in the Map/List
TGNumberEntryField * fMapLabLframe[6]
! The GUI number entries for the local frame specification
void SetGmatrix(TString mode)
Double_t GetNeutrinoXsection(Int_t mode, Int_t type, Double_t egev, Double_t xscale=1, Double_t *eprimgev=0, Double_t *alpha=0) const
TH1 * GetBurstT90dist(TString name, TString type)
NcDevice * GetBurstParameters()
void ListBurstHistograms() const
void SetLT(Int_t y, Int_t m, Int_t d, Int_t hh, Int_t mm, Int_t ss, Int_t ns=0, Int_t ps=0)
void MapLocl(const char *text)
void MakeBurstDataStats(Int_t mode, Int_t nmugrb=0)
void SetCentralMeridian(Int_t mode=0, Double_t phi=0, TString u="deg")
Int_t fGal
! Type indicator for fG values (0=uninitialised 1=B1950 2=J2000)
Int_t fMapLabId
! The GUI entered Lab detector Id
void SetEmatrix(NcTimestamp *ts)
Int_t RemoveRefSignal(Int_t j, Int_t compress)
Int_t fMapMarkColor
! The GUI selected marker color for the Map
void MakeBurstT90dist(TString file, TString tree, TString name, Int_t nb=50, Float_t xmin=-5, Float_t xmax=5)
void DisplaySignal(TString frame, TString mode, NcTimestamp *ts, Int_t j=-1, TString proj="ham", Int_t clr=0, TString name="")
virtual void TimestampPanel(TGCompositeFrame *frame)
void SetHmatrix(NcTimestamp *ts)
void PrintSignal(TString frame, TString mode, NcTimestamp *ts, Int_t ndig, Int_t jref=0, TString emode="T", Int_t type=0, Bool_t align=kFALSE)
void SetLabTimeOffset(Double_t dt)
TH1F GetCumulHistogram(TH1 *h, TString name, TString mode="F") const
virtual void MapListOptionsPanel(TGCompositeFrame *frame)
TObjArray * fMarkers
! Temp. array to hold the markers for the signal display
TString fMapDateTime
! The GUI entered datetime
Double_t GetSolidAngle(Double_t thetamin, Double_t thetamax, TString tu, Double_t phimin, Double_t phimax, TString pu) const
void SetMarkerColor(Int_t color, Int_t type)
Int_t fMapNdigs
! The GUI selected number of digits for the List output
TRotMatrix fE
! Matrix for conversion of equatorial to ecliptic coordinates
Double_t GetProperDistance(Double_t z, TString u="Mpc", Int_t t=1) const
TString fMapIname
! The GUI selected entry name for the info
Double_t fMapEa
! The GUI entered a coordinate of an entry
Bool_t fMapSolar[10]
! The GUI selection of solar system objects
void GetBurstChi2Statistics(TString type, Int_t ndt=2, Bool_t zcor=kFALSE)
Double_t GetBurstSignalEnergy(Double_t Emin=-1, Double_t Emax=-1) const
void MapLocId(const char *text)
Bool_t fMapDoptions[5]
! The GUI Map/List options (histo, clr, ref, meas, refTS)
Int_t RemoveSignals(TString name, Int_t type, Int_t compress)
void SetPhysicalParameter(TString name, Double_t value)
void MatchBurstData(NcDevice &matches, Int_t i1=1, Int_t i2=0, Int_t itype=0, Int_t j1=1, Int_t j2=0, Int_t jtype=1)
void ShiftPosition(Nc3Vector &v, Double_t angle)
void MakeBurstZdist(TString file, TString tree, TString name, Int_t nb=200, Float_t zmin=0, Float_t zmax=20)
TArrayI * MatchRefSignal(Double_t da, TString au, Double_t dt, TString tu, Int_t mode=1)
void WriteBurstHistograms(TString filename)
void GeoToHeliocentric(Double_t &R, Double_t &B, Double_t &L, NcTimestamp *ts, TString Bu="deg", TString Lu="deg")
Double_t GetSignalRateProb(Double_t *vars, Double_t *pars)
Double_t GetNuclearMass(Int_t Z, Int_t N, Int_t mode=1) const
virtual void InfoPanel(TGCompositeFrame *frame)
void BurstCompensate(Int_t &nmugrb)
void SetLocalFrame(Double_t t1, Double_t p1, Double_t t2, Double_t p2, Double_t t3, Double_t p3)
void Data(Int_t mode=1, TString u="deg", Bool_t utc=kTRUE)
void ProjectAitoff(Double_t l, Double_t b, Double_t &x, Double_t &y)
void SetMarkerSize(Float_t size, Int_t type)
void DisplaySignals(TString frame, TString mode, NcTimestamp *ts, TString proj="ham", Int_t clr=0, Int_t nmax=-1, Int_t j=-1, Int_t type=-1, TString name="*")
NcSignal * GetSignal(Double_t &d, Double_t &a, TString au, Double_t &b, TString bu, TString frame, NcTimestamp *ts, Int_t jref, TString mode="T", Int_t type=0)
virtual void LabLocalFramePanel(TGCompositeFrame *frame)
TString fMapMerUc
! The GUI selected angular units for the central meridian location
void SetExperiment(TString name, Int_t id=0)
TH1 * GetBurstZdist(TString name, TString type)
TString fMapUinfo
! The GUI selected angular units for the Lab info
TString fMapEname
! The GUI entered name of the entry
Int_t GetLabDetectorId() const
Double_t GetPhysicalParameter(TString name) const
virtual TObject * Clone(const char *name="") const
Double_t GetNeutrinoAngle(Double_t E, TString u, Int_t mode, TF1 *f=0)
TGMainFrame * fSkyMapPanel
! The main frame for the SkyMapPanel GUI
Int_t GetNsignals(Int_t type, Int_t mode=0) const
void SetMarkerStyle(Int_t style, Int_t type)
NcSignal * SetSignal(Double_t d, Double_t a, TString au, Double_t b, TString bu, TString frame, NcTimestamp *ts, Int_t jref, TString mode="T", TString name="", Int_t type=0)
Int_t fMapMerMode
! The GUI selected meridian orientation for the Map
void PrintAngle(Double_t a, TString in, TString out, Int_t ndig=1, Bool_t align=kFALSE) const
void MapMerC(const char *text)
Int_t fMapTinfo
! The GUI selected mode for the timestamp info
void SetMaxDt(Double_t s)
TH1F GetLogHistogram(TH1 *hin, Int_t mode, TString s="") const
NcTimestamp fMapTS
! The GUI entered timestamp to be used for the List/Map
Double_t GetPhysicalDistance(Double_t z, TString u="Mpc", Int_t t=1) const
void MapTimeType(Int_t i)
Double_t KolmogorovTest(TString mode, TH1 *h1, TH1 *h2=0, TF1 *pdf=0, Double_t nr=1000, TH1F *ksh=0, Int_t ncut=0, Double_t *nrx=0, Int_t mark=1)
TGComboBox * fMapLabE
! The GUI Lab experiment site selection box
virtual void CommandPanel(TGCompositeFrame *frame)
Double_t GetSurvivalProbability(Double_t x, Double_t lambda) const
Double_t GetHourAngle(TString mode, NcTimestamp *ts, Int_t jref=0, Int_t type=0)
Double_t GetRadiationLength(Double_t Z, Double_t A, Double_t rho=-1) const
void MakeBurstEnergydist(Int_t mode, TString file, TString tree, TString name1, TString name2, TString u, Double_t Emin, Double_t Emax, Int_t nb=1000)
Int_t fMapNmax
! The GUI selected max. number of signals of each type to Map/List
Int_t fBias
! Initialisation flag for fB values (0=uninitialised 1=initialised)
void RandomPosition(Nc3Vector &v, Double_t thetamin, Double_t thetamax, Double_t phimin, Double_t phimax)
TString fMapLabExpName
! The GUI entered Lab experimental site
TF1 GetBackgroundRatePDF(Double_t Noff, Double_t Toff, Double_t bmax=-1, Double_t prec=709)
TString fMapEmode
! The GUI entered coordinate system mode of the entry
NcRandom * GetRandomiser(Int_t &iseed, Int_t &cnt1, Int_t &cnt2) const
void GetBurstDtDistributions(Int_t ndt, TH1F &hisdtOn, TF1 &pdfdtOn, TH1F &hisdtOff, TF1 &pdfdtOff, Bool_t zcor)
void MapNmax(const char *text)
TH1F GetCountsHistogram(TF1 &spec, Int_t nbins, Double_t xmin, Double_t xmax, Int_t mode, TString s="") const
Double_t GetSeparation(TString name1, TString name2, TString au, Double_t *dt=0, TString tu="s", Int_t mode=1, Double_t *diftheta=0, Double_t *difphi=0)
void AddNamedSlot(TString s)
Int_t GetSlotIndex(TString name, Int_t opt=0) const
(Bayesian) Block treatment of sequential data.
Int_t Rebin(TH1 *hin, TH1 *hout, Bool_t scale, Int_t nbins=0, Double_t xmin=0, Double_t xmax=-1)
Int_t Divide(TH1 *h1, TH1 *h2, TH1 *hout, Bool_t scale, Double_t c, Double_t d=0)
Double_t GetBlocks(TH1 *hin, Double_t fpr, TH1 *hout, Int_t ntrig=0)
Signal (Hit) handling of a generic device.
virtual TObject * Clone(const char *name="") const
NcSignal * GetHit(Int_t j) const
virtual void Reset(Int_t mode=0)
Various mathematical tools for scientific analysis.
Double_t PsiValue(Int_t m, Int_t *n, Double_t *p=0, Int_t f=0) const
Double_t LnGamma(Double_t z) const
Double_t Chi2Pvalue(Double_t chi2, Int_t ndf, Int_t sides=0, Int_t sigma=0, Int_t mode=1) const
Double_t LiMaSignificance(Double_t Non, Double_t Ton, Double_t Noff, Double_t Toff, Double_t Ra=1, Double_t Re=1) const
Double_t Gamma(Double_t z) const
Double_t PsiExtreme(Double_t n, Int_t m, Double_t *p=0, Int_t k=0) const
Double_t PsiPvalue(Double_t psi0, Double_t nr, Double_t n, Int_t m, Double_t *p=0, Int_t f=0, Double_t *na=0, TH1F *psih=0, Int_t ncut=0, Double_t *nrx=0, Int_t mark=1)
TF1 PoissonDtDist(Double_t r, Int_t n) const
Double_t Chi2Value(Int_t m, Int_t *n, Double_t *p=0, Int_t *ndf=0) const
Handling of positions (with timestamps) in various reference frames.
void GetPosition(Double_t *r, TString f, TString u="rad", Float_t s=-1) const
void SetPosition(Double_t *r, TString f, TString u="rad")
NcTimestamp * GetTimestamp()
void SetTimestamp(NcTimestamp &t)
Generate universal random numbers and sequences on all common machines.
Sampling and statistics tools for various multi-dimensional data samples.
Double_t GetMedian(Int_t i)
void SetStoreMode(Int_t mode=1, Int_t nmax=0, Int_t i=0)
Double_t GetEntry(Int_t i, Int_t j, Int_t mode=0, Int_t k=0)
Generic handling of (extrapolated) detector signals.
virtual void Reset(Int_t mode=0)
virtual void SetSignal(Double_t sig, Int_t j=1)
virtual Float_t GetSignal(Int_t j=1, Int_t mode=0) const
virtual void Data(TString f="car", TString u="rad") const
Double_t GetEpoch(TString mode)
Double_t GetMJD(Int_t y, Int_t m, Int_t d, Int_t hh, Int_t mm, Int_t ss, Int_t ns) const
void SetLT(Double_t dt, Int_t y, Int_t m, Int_t d, Int_t hh, Int_t mm, Int_t ss, Int_t ns=0, Int_t ps=0, TString utc="A", Int_t leap=0, Double_t dut=0)
Double_t GetJD(Int_t y, Int_t m, Int_t d, Int_t hh, Int_t mm, Int_t ss, Int_t ns) const
void GetUT(Int_t &hh, Int_t &mm, Int_t &ss, Int_t &ns, Int_t &ps)
TString GetDayTimeString(TString mode, Int_t ndig=0, Double_t offset=0, TString *date=0, TString *time=0, Bool_t full=kTRUE)
Double_t Almanac(Double_t *dpsi=0, Double_t *deps=0, Double_t *eps=0, Double_t *dl=0, TString name="", Double_t *el=0, Double_t *eb=0, Double_t *dr=0, Double_t *value=0, Int_t j=0)
Double_t GetLMST(Double_t offset)
Double_t GetLT(Double_t offset)
void GetGMST(Int_t &hh, Int_t &mm, Int_t &ss, Int_t &ns, Int_t &ps)
void SetEpoch(Double_t e, TString mode, TString utc="U", Int_t leap=0, Double_t dut=0)
void SetUT(Int_t y, Int_t m, Int_t d, Int_t hh, Int_t mm, Int_t ss, Int_t ns=0, Int_t ps=0, TString utc="A", Int_t leap=0, Double_t dut=0)
void PrintTime(Double_t h, Int_t ndig=1) const
Double_t GetLAST(Double_t offset)
TTree * LoadUTCparameterFiles(TString leapfile="$(NCFS)/IERS/leap.txt", TString dutfile="$(NCFS)/IERS/dut.txt")
void SetTJD(Int_t tjd, Int_t sec, Int_t ns, Int_t ps=0, TString utc="U", Int_t leap=0, Double_t dut=0)
void SetMJD(Int_t mjd, Int_t sec, Int_t ns, Int_t ps=0, TString utc="U", Int_t leap=0, Double_t dut=0)
void SetJD(Int_t jd, Int_t sec, Int_t ns, Int_t ps=0, TString utc="U", Int_t leap=0, Double_t dut=0)
void Add(Int_t d, Int_t s, Int_t ns, Int_t ps=0)
void Convert(Double_t date, Int_t &days, Int_t &secs, Int_t &ns) const
Int_t GetDifference(NcTimestamp *t, Int_t &days, Int_t &sec, Int_t &ns, Int_t &ps, TString type="UT")
void AddSec(Double_t seconds)
void Date(Int_t mode=3, Double_t offset=0)