3 #include "TPaletteAxis.h"
89 if (resolution<=0) {cout <<
"gskymap::gskymap: Error sky resolution must be >0" << endl;
exit(1);}
91 if(gSystem->Getenv(
"CWB_GWAT")!=NULL) {
92 worldMapPath=
TString(gSystem->Getenv(
"CWB_GWAT"))+
"/data";
99 char sname[64];
sprintf(sname,
"gskymap_%d",
int(gRandom->Rndm(13)*1.e9));
101 this->projection = projection;
102 this->coordinate = coordinate;
103 this->resolution = resolution;
106 wtopx=35; wtopy=46;
ww=817; wh=472;
112 if (coordinate.CompareTo(
"CWB")==0) {phMin=0;phMax=360;thMin=0;thMax=180;}
113 if (coordinate.CompareTo(
"CELESTIAL")==0) {phMin=0;phMax=360;thMin=-90;thMax=90;}
115 if(h2!=NULL)
delete h2;
117 h2 =
new TH2D(nameh,
"gskymap",
int(360*resolution), phMin, phMax,
int(180*resolution), thMin, thMax);
118 h2->SetStats(kFALSE);
119 h2->SetTitleFont(12);
121 h2->SetFillColor(kWhite);
123 h2->GetXaxis()->SetNdivisions(-604);
124 h2->GetXaxis()->SetLabelFont(42);
125 h2->GetXaxis()->SetLabelOffset(0.014);
126 h2->GetXaxis()->SetTitleOffset(1.2);
127 h2->GetYaxis()->SetTitleOffset(0.8);
128 h2->GetYaxis()->SetNdivisions(-602);
129 h2->GetYaxis()->SetLabelFont(42);
130 h2->GetYaxis()->SetLabelOffset(0.01);
131 h2->GetZaxis()->SetLabelFont(42);
132 h2->GetZaxis()->SetNoExponent(
false);
134 h2->GetXaxis()->SetTitleFont(42);
135 if ((coordinate.CompareTo(
"CWB")==0)||(coordinate.CompareTo(
"GEOGRAPHIC")==0)) {
136 h2->GetXaxis()->SetTitle(
"#phi (deg)");
137 }
else if (coordinate.CompareTo(
"CELESTIAL")==0) {
138 h2->GetXaxis()->SetTitle(
"RA (deg)");
140 {cout <<
"gskymap::gskymap: Error : coordinate system not declared" << endl;
exit(1);}
142 h2->GetXaxis()->CenterTitle(
true);
144 h2->GetYaxis()->SetTitleFont(42);
145 if ((coordinate.CompareTo(
"CWB")==0)||(coordinate.CompareTo(
"GEOGRAPHIC")==0)) {
146 h2->GetYaxis()->SetTitle(
"#theta (deg)");
147 }
else if (coordinate.CompareTo(
"CELESTIAL")==0) {
148 h2->GetYaxis()->SetTitle(
"DEC (deg)");
150 {cout <<
"gskymap::gskymap: Error : coordinate system not declared" << endl;
exit(1);}
152 h2->GetYaxis()->CenterTitle(
true);
154 h2->GetZaxis()->SetTitleOffset(0.6);
155 h2->GetZaxis()->SetTitleFont(42);
156 h2->GetZaxis()->SetTitle(zAxisTitle);
157 h2->GetZaxis()->CenterTitle(
true);
159 if(this->
size()) FillData();
176 canvas->ToggleEventStatus();
177 canvas->SetGridx(
false);
178 canvas->SetGridy(
false);
180 canvas->SetFillColor(kWhite);
181 canvas->SetLeftMargin(0.08);
182 canvas->SetBottomMargin(0.13);
183 canvas->SetBorderMode(0);
191 gStyle->SetFrameBorderMode(0);
194 gStyle->SetTitleH(0.050);
195 gStyle->SetTitleW(0.95);
196 gStyle->SetTitleY(0.98);
197 gStyle->SetTitleFont(12,
"D");
198 gStyle->SetTitleColor(kBlue,
"D");
199 gStyle->SetTextFont(12);
200 gStyle->SetTitleFillColor(kWhite);
201 gStyle->SetLineColor(kWhite);
202 gStyle->SetNumberContours(256);
203 gStyle->SetMarkerStyle(7);
204 gStyle->SetMarkerSize(2);
205 gStyle->SetCanvasColor(kWhite);
206 gStyle->SetStatBorderSize(1);
215 if(
h2!=NULL)
delete h2;
253 if (!in.good()) {cout <<
"gskymap::FillData: Error Opening File : " << fname << endl;
exit(1);}
255 in.getline(iline,1024);
256 in.getline(iline,1024);
257 if(
TString(iline).Contains(
"skymap")==0)
258 {cout <<
"gskymap::FillData: " << fname <<
" is not a skymap file" << endl;
exit(1);}
259 in.getline(iline,1024);
261 TObjString* tsms = (TObjString*)tok->At(1);
262 if(tsms->GetString().IsAlpha()==1)
263 {cout <<
"gskymap::FillData: resolution is not well defined " << endl;
exit(1);}
264 double sms = tsms->GetString().Atof();
265 in.getline(iline,1024);
267 TObjString* ttheta_1 = (TObjString*)tok->At(1);
268 if(ttheta_1->GetString().IsAlpha()==1)
269 {cout <<
"gskymap::FillData: theta_1 is not well defined " << endl;
exit(1);}
270 double theta_1 = ttheta_1->GetString().Atof();
271 in.getline(iline,1024);
273 TObjString* ttheta_2 = (TObjString*)tok->At(1);
274 if(ttheta_2->GetString().IsAlpha()==1)
275 {cout <<
"gskymap::FillData: theta_2 is not well defined " << endl;
exit(1);}
276 double theta_2 = ttheta_2->GetString().Atof();
277 in.getline(iline,1024);
279 TObjString* tphi_1 = (TObjString*)tok->At(1);
280 if(tphi_1->GetString().IsAlpha()==1)
281 {cout <<
"gskymap::FillData: phi_1 is not well defined " << endl;
exit(1);}
282 double phi_1 = tphi_1->GetString().Atof();
283 in.getline(iline,1024);
285 TObjString* tphi_2 = (TObjString*)tok->At(1);
286 if(tphi_2->GetString().IsAlpha()==1)
287 {cout <<
"gskymap::FillData: phi_2 is not well defined " << endl;
exit(1);}
288 double phi_2 = tphi_2->GetString().Atof();
289 in.getline(iline,1024);
291 TObjString* tlenght = (TObjString*)tok->At(1);
292 if(tlenght->GetString().IsAlpha()==1)
293 {cout <<
"gskymap::FillData: lenght is not well defined " << endl;
exit(1);}
294 int lenght = tlenght->GetString().Atoi();
302 {cout <<
"gskymap::FillData: lenght is not consistent" << endl;
exit(1);}
306 in.getline(iline,1024);
307 if (!in.good())
break;
311 if (tok->GetEntries()==1)
continue;
313 TObjString* tindex = (TObjString*)tok->At(0);
314 TObjString* tvalue = (TObjString*)tok->At(1);
316 if (tindex->GetString().IsAlpha()==1)
317 {cout <<
"gskymap::FillData: index " << tindex->GetString().Data() <<
" is not well defined " << endl;
exit(1);}
321 int index = tindex->GetString().Atoi();
322 double value = tvalue->GetString().Atof();
324 if(index>=0&&index<L)
327 {cout <<
"gskymap::FillData: index " << index <<
" not allowed - max " << L-1 << endl;}
346 double*
phi =
new double[
size];
348 double* binc =
new double[
size];
353 phi[
k]=(double)
j/(2*resolution);
354 theta[
k]=(double)
i/(2*resolution);
389 if(phi[
i]<0) phi[
i]+=360.;
390 if(phi[
i]>360) phi[
i]-=360.;
391 if(theta[
i]<0) theta[
i]+=180.;
392 if(theta[
i]>180) theta[
i]-=180.;
395 if(phi[
i]<-180) phi[
i]+=360.;
396 if(phi[
i]>180) phi[
i]-=360.;
397 if(theta[
i]<-90) theta[
i]+=180.;
398 if(theta[
i]>90) theta[
i]-=180.;
403 if(phi[
i]<0) phi[
i]+=360.;
404 if(phi[
i]>360) phi[
i]-=360.;
405 if(theta[
i]<-90) theta[
i]+=180.;
406 if(theta[
i]>90) theta[
i]-=180.;
427 h2->SetBinContent(
int(x*
resolution)+1,
int(y*resolution)+1,binc[
i]);
454 TGaxis::SetMaxDigits(3);
462 gStyle->SetPalette(1,0);
468 gStyle->SetPalette(1,0);
477 TPaletteAxis *
palette = (TPaletteAxis*)
h2->GetListOfFunctions()->FindObject(
"palette");
479 palette->SetX1NDC(0.91);
480 palette->SetX2NDC(0.933);
481 palette->SetTitleOffset(0.92);
482 palette->GetAxis()->SetTickSize(0.01);
486 Double_t rlonait[360*4],rlatait[360*4];
494 for (kk=0;kk<=360;kk+=15){
497 for (k=0;k<180*4;k++){
513 if (
coordinate.CompareTo(
"CWB")==0) {x+=180;y+=90;}
515 if (
coordinate.CompareTo(
"CELESTIAL")==0) {x+=180;x=360-
x;}
518 double distance = klen>0 ? sqrt(pow(x-rlonait[klen-1],2)+pow(y-rlatait[klen-1],2)) : mdistance;
519 if(distance>=mdistance || k==(180*4-1)) klen++;
521 TPolyLine *pL =
new TPolyLine(klen,rlonait,rlatait);
523 if(kk==0 || kk==360) pL->SetLineStyle(1);
else pL->SetLineStyle(3);
533 for (kk=15;kk<180;kk+=15){
536 for (k=0;k<360*4;k++){
552 if (
coordinate.CompareTo(
"CWB")==0) {x+=180;y+=90;}
554 if (
coordinate.CompareTo(
"CELESTIAL")==0) {x+=180;x=360-
x;}
557 double distance = klen>0 ? sqrt(pow(x-rlonait[klen-1],2)+pow(y-rlatait[klen-1],2)) : mdistance;
558 if(distance>=mdistance || k==(360*4-1)) klen++;
560 TPolyLine *pL =
new TPolyLine(klen,rlonait,rlatait);
562 if(kk==0 || kk==180) pL->SetLineStyle(1);
else pL->SetLineStyle(3);
618 for (
int l=0;
l<
L;
l++) {
629 for (
int l=0;
l<
L;
l++) {
636 }
else if (
projection.CompareTo(
"sinusoidal")==0) {
639 }
else if (
projection.CompareTo(
"parabolic")==0) {
652 TPolyLine *pL =
new TPolyLine(P,x.
data,y.
data);
671 if (
coordinate.CompareTo(
"CWB")==0) {phi+=180;
if(phi>180) phi=phi-360;}
672 if (
coordinate.CompareTo(
"CELESTIAL")==0) {phi+=180;
if(phi>180) phi=phi-360;phi=360-
phi;}
687 TMarker *mP =
new TMarker(x,y, 1);
688 mP->SetMarkerColor(kGray+1);
698 TGaxis::SetMaxDigits();
719 DrawMarker(phi, theta, marker, msize, mcolor);
753 TMarker *mP =
new TMarker(x,y, marker);
754 mP->SetMarkerSize(msize);
755 mP->SetMarkerColor(mcolor);
776 DrawText(phi, theta, text, tsize, tcolor);
810 TText *tP =
new TText(x,y,text);
811 tP->SetTextSize(tsize);
812 tP->SetTextColor(tcolor);
838 Double_t alpha2 = (l/2)*TMath::DegToRad();
839 Double_t
delta = b*TMath::DegToRad();
840 Double_t
r2 = TMath::Sqrt(2.);
843 Double_t denom = TMath::Sqrt(1. + cdec*
TMath::Cos(alpha2));
844 x = cdec*TMath::Sin(alpha2)*2.*r2/denom;
845 y = TMath::Sin(delta)*r2/denom;
846 x *= TMath::RadToDeg()/
f;
847 y *= TMath::RadToDeg()/
f;
868 Al = l*cos(b*TMath::DegToRad());
886 Al = l*(2.*
TMath::Cos(2*b*TMath::DegToRad()/3) - 1);
887 Ab = 180*TMath::Sin(b*TMath::DegToRad()/3);
904 const Int_t
NRGBs = 5;
905 const Int_t
NCont = 255;
907 if (
fabs(paletteId)==1) {
908 Double_t
stops[
NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
909 Double_t
red[
NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
910 Double_t
green[
NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
911 Double_t
blue[
NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 };
913 TColor::CreateGradientColorTable(NRGBs, stops, blue, green, red, NCont);
915 TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
918 if (
fabs(paletteId)==2) {
919 Double_t
stops[
NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
920 Double_t
red[
NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
922 Double_t
green[
NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
924 Double_t
blue[
NRGBs] = { 1.00, 1.00, 0.00, 0.00, 0.00 };
926 TColor::CreateGradientColorTable(NRGBs, stops, blue, green, red, NCont);
928 TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
931 if (
fabs(paletteId)==3) {
932 Double_t
stops[
NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
933 Double_t
red[
NRGBs] = { 0.00, 0.09, 0.18, 0.09, 0.00 };
934 Double_t
green[
NRGBs] = { 0.01, 0.02, 0.39, 0.68, 0.97 };
935 Double_t
blue[
NRGBs] = { 0.17, 0.39, 0.62, 0.79, 0.97 };
937 TColor::CreateGradientColorTable(NRGBs, stops, blue, green, red, NCont);
939 TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
942 if (
fabs(paletteId)==4) {
943 Double_t
stops[
NRGBs] = { 0.00, 0.50, 0.75, 0.875, 1.00 };
944 Double_t
red[
NRGBs] = { 1.00, 1.00, 1.00, 1.00, 1.00 };
945 Double_t
green[
NRGBs] = { 1.00, 0.75, 0.50, 0.25, 0.00 };
946 Double_t
blue[
NRGBs] = { 0.00, 0.00, 0.00, 0.00, 0.00 };
948 TColor::CreateGradientColorTable(NRGBs, stops, blue, green, red, NCont);
950 TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
953 if (
fabs(paletteId)==5) {
954 Double_t
stops[
NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
955 Double_t
red[
NRGBs] = { 1.00, 0.84, 0.61, 0.34, 0.00 };
956 Double_t
green[
NRGBs] = { 1.00, 0.84, 0.61, 0.34, 0.00 };
957 Double_t
blue[
NRGBs] = { 1.00, 0.84, 0.61, 0.34, 0.00 };
959 TColor::CreateGradientColorTable(NRGBs, stops, blue, green, red, NCont);
961 TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
964 if (
fabs(paletteId)==57) {
965 Double_t
stops[9] = { 0., 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0};
966 Double_t
red[9] = { 0.2082, 0.0592, 0.0780, 0.0232, 0.1802, 0.5301, 0.8186, 0.9956, 0.9764};
967 Double_t
green[9] = { 0.1664, 0.3599, 0.5041, 0.6419, 0.7178, 0.7492, 0.7328, 0.7862, 0.9832};
968 Double_t
blue[9] = { 0.5293, 0.8684, 0.8385, 0.7914, 0.6425, 0.4662, 0.3499, 0.1968, 0.0539};
970 TColor::CreateGradientColorTable(9, stops, blue, green, red, NCont);
972 TColor::CreateGradientColorTable(9, stops, red, green, blue, NCont);
976 gStyle->SetNumberContours(NCont);
1002 if (!in.good()) {cout <<
"gskymap::ReadWorlMapCoastLine: Error Opening File : " << ifileName << endl;
exit(1);}
1007 wm_lon =
new double[
size];
1008 wm_lat =
new double[
size];
1011 in.getline(iline,1024);
1012 if (!in.good())
break;
1016 if (tok->GetEntries()==1)
continue;
1018 TObjString*
tra = (TObjString*)tok->At(0);
1019 TObjString*
tdec = (TObjString*)tok->At(1);
1021 if (tdec->GetString().IsAlpha()==1)
continue;
1024 double ra = tra->GetString().Atof();
1025 double dec = tdec->GetString().Atof();
1027 if (
fabs(dec)<0.001)
continue;
1034 cout << cnt << endl;
1111 if(
canvas==NULL) {cout <<
"gskymap::Print: Error - canvas not initialized " << endl;
exit(1);}
1113 TGaxis::SetMaxDigits(3);
1133 pname.ReplaceAll(
".PNG",
".png");
1136 TGaxis::SetMaxDigits();
1164 if((fp = fopen(fname,
"w")) == NULL) {
1165 cout <<
"netevent::DumpSkyMap() error: cannot open file " << fname <<
". \n";
1172 x = double(
int(1000*x))/1000.;
1194 for (
int i=1;
i<=h2->GetNbinsX();
i++) {
1196 double xvalue = h2->GetXaxis()->GetBinCenter(
i)+h2->GetXaxis()->GetBinWidth(
i)/2.;
1197 int nbin=
int(h2->GetNbinsX()/4.);
1199 sprintf(xlabel,
"%d",
int(360-xvalue));
1201 h2->GetXaxis()->SetBinLabel(
i,xlabel);
1202 h2->GetXaxis()->SetLabelSize(0.06);
1203 h2->GetXaxis()->LabelsOption(
"h");
1204 h2->GetXaxis()->SetNdivisions(604);
1205 h2->GetXaxis()->SetLabelFont(42);
1206 h2->GetXaxis()->SetLabelOffset(0.010);
1207 h2->GetXaxis()->SetTitleOffset(1.5);
1234 cout <<
"gskymap::LoadObject : Error - input file don't contains object gskymap" << endl;
std::vector< TPolyLine * > gdL
fprintf(stdout,"start=%f duration=%f rate=%f\n", x.start(), x.size()/x.rate(), x.rate())
std::vector< TPolyLine * > gridxL
void CwbToCelestial(double ilongitude, double ilatitude, double &olongitude, double &olatitude, double gps=0)
void SetPlotStyle(int paletteId=1)
cout<< "skymap size : "<< L<< endl;for(int l=0;l< L;l++) sm.set(l, l);sm > const_cast< char * >("skymap.dat")
char * watversion(char c='s')
cout<< endl;cout<< "ts size = "<< ts.size()<< " ts rate = "<< ts.rate()<< endl;tf.Forward(ts, wdm);int levels=tf.getLevel();cout<< "tf size = "<< tf.size()<< endl;double dF=tf.resolution();double dT=1./(2 *dF);cout<< "rate(hz) : "<< RATE<< "\t layers : "<< nLAYERS<< "\t dF(hz) : "<< dF<< "\t dT(ms) : "<< dT *1000.<< endl;int itime=TIME_PIXEL_INDEX;int ifreq=FREQ_PIXEL_INDEX;int index=(levels+1)*itime+ifreq;double time=itime *dT;double freq=(ifreq >0)?ifreq *dF:dF/4;cout<< endl;cout<< "PIXEL TIME = "<< time<< " sec "<< endl;cout<< "PIXEL FREQ = "<< freq<< " Hz "<< endl;cout<< endl;wavearray< double > x
void DrawMarker(double phi, double theta, int marker, Size_t msize=1, Color_t tcolor=1)
void DrawText(double phi, double theta, TString text, double tsize=0.04, Color_t tcolor=1)
void Draw(int dpaletteId=0, Option_t *option="colfz")
std::vector< vectorD > value
gSM SetOptions(cwb_antpat_projection, COORDINATES, RESOLUTION/2)
std::vector< TText * > axisT
double RA2phi(double ph, double gps)
void ReverseXAxis(TH2D *h2)
void CelestialToCwb(double ilongitude, double ilatitude, double &olongitude, double &olatitude, double gps=0)
void ProjectParabolic(Double_t l, Double_t b, Double_t &Al, Double_t &Ab)
double fabs(const Complex &x)
void ProjectHammer(Double_t l, Double_t b, Double_t &Al, Double_t &Ab)
void LoadObject(const char *file, const char *name="gskymap")
void CwbToGeographic(double ilongitude, double ilatitude, double &olongitude, double &olatitude)
void DumpObject(const char *file, const char *name="gskymap")
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
void GalacticToEquatorial(double ilongitude, double ilatitude, double &olongitude, double &olatitude)
void set(size_t i, double a)
param: sky index param: value to set
void ProjectSinusoidal(Double_t l, Double_t b, Double_t &Al, Double_t &Ab)
void DumpSkyMap(char *fname)
gskymap() gskymap(double sms, double t1=0., double t2=180., double p1=0., double p2=360.) gskymap(char *ifile changed)
double get(size_t i)
param: sky index
void GeographicToCwb(double ilongitude, double ilatitude, double &olongitude, double &olatitude)
Color_t colorGalacticDisk
int ReadWorlMapCoastLine(double *&wm_lon, double *&wm_lat)
void Print(TString pname)
std::vector< TPolyLine * > gridyL
std::vector< TMarker * > wmM