| 1 | #include <iomanip>
|
|---|
| 2 |
|
|---|
| 3 | Double_t GetResidual(Double_t fRawEl, Double_t fRawAz,
|
|---|
| 4 | Double_t fStarEl, Double_t fStarAz)
|
|---|
| 5 | {
|
|---|
| 6 | fRawEl *= TMath::DegToRad();
|
|---|
| 7 | fRawAz *= TMath::DegToRad();
|
|---|
| 8 | fStarEl *= TMath::DegToRad();
|
|---|
| 9 | fStarAz *= TMath::DegToRad();
|
|---|
| 10 | Double_t del = fRawEl-fStarEl;
|
|---|
| 11 | Double_t daz = fRawAz-fStarAz;
|
|---|
| 12 | Double_t dphi2 = daz/2.;
|
|---|
| 13 | Double_t cos2 = cos(dphi2)*cos(dphi2);
|
|---|
| 14 | Double_t sin2 = sin(dphi2)*sin(dphi2);
|
|---|
| 15 | Double_t d = cos(del)*cos2 - cos(fRawEl+fStarEl)*sin2;
|
|---|
| 16 |
|
|---|
| 17 | Double_t dist = acos(d);
|
|---|
| 18 |
|
|---|
| 19 | return dist * TMath::RadToDeg();
|
|---|
| 20 | }
|
|---|
| 21 |
|
|---|
| 22 | void DrawMarker(TVirtualPad *pad, Double_t r0, Double_t phi0, Double_t r1, Double_t phi1)
|
|---|
| 23 | {
|
|---|
| 24 | TView *view = pad->GetView();
|
|---|
| 25 |
|
|---|
| 26 | if (!view)
|
|---|
| 27 | {
|
|---|
| 28 | cout << "No View!" << endl;
|
|---|
| 29 | return;
|
|---|
| 30 | }
|
|---|
| 31 |
|
|---|
| 32 | TMarker mark0;
|
|---|
| 33 | TMarker mark1;
|
|---|
| 34 | mark0.SetMarkerStyle(kStar);
|
|---|
| 35 | mark1.SetMarkerStyle(kStar);
|
|---|
| 36 | mark1.SetMarkerColor(kRed);
|
|---|
| 37 |
|
|---|
| 38 | r0 /= 90;
|
|---|
| 39 | r1 /= 90;
|
|---|
| 40 | phi0 *= TMath::DegToRad();
|
|---|
| 41 | phi1 *= TMath::DegToRad();
|
|---|
| 42 |
|
|---|
| 43 | Double_t x0[3] = { r0*cos(phi0), r0*sin(phi0), 0};
|
|---|
| 44 | Double_t x1[3] = { r1*cos(phi1), r1*sin(phi1), 0};
|
|---|
| 45 |
|
|---|
| 46 | mark0.DrawMarker(x0[0], x0[1]);
|
|---|
| 47 | mark1.DrawMarker(x1[0], x1[1]);
|
|---|
| 48 |
|
|---|
| 49 | return;
|
|---|
| 50 | Double_t y0[3], y1[3];
|
|---|
| 51 |
|
|---|
| 52 | view->WCtoNDC(x0, y0);
|
|---|
| 53 | view->WCtoNDC(x1, y1);
|
|---|
| 54 |
|
|---|
| 55 | mark0.DrawMarker(y0[0], y0[1]);
|
|---|
| 56 | mark1.DrawMarker(y1[0], y1[1]);
|
|---|
| 57 | }
|
|---|
| 58 |
|
|---|
| 59 | int fill(const char *fname, TGraph *g, TH1 *h)
|
|---|
| 60 | {
|
|---|
| 61 | /*
|
|---|
| 62 | TH2F h2res1("Res2D1", " Dataset positions on the sky ", 32, 0, 360, 10, 0, 90);
|
|---|
| 63 | h2res1.SetBit(TH1::kNoStats);
|
|---|
| 64 | h2res1.DrawCopy("surf1pol");
|
|---|
| 65 | gPad->Modified();
|
|---|
| 66 | gPad->Update();
|
|---|
| 67 | gPad->SetTheta(90);
|
|---|
| 68 | gPad->SetPhi(-90);
|
|---|
| 69 |
|
|---|
| 70 | DrawMarker(gPad, 45, 0, 0, 0);
|
|---|
| 71 | gPad->Modified();
|
|---|
| 72 | gPad->Update();
|
|---|
| 73 |
|
|---|
| 74 | return;
|
|---|
| 75 |
|
|---|
| 76 | */
|
|---|
| 77 | ifstream fin(fname);
|
|---|
| 78 |
|
|---|
| 79 | cout << "Reading " << setw(23) << fname << "..." << flush;
|
|---|
| 80 |
|
|---|
| 81 | while (1)
|
|---|
| 82 | {
|
|---|
| 83 | TString str;
|
|---|
| 84 | str.ReadLine(fin);
|
|---|
| 85 | if (!fin)
|
|---|
| 86 | break;
|
|---|
| 87 |
|
|---|
| 88 | if (str.Contains("#"))
|
|---|
| 89 | continue;
|
|---|
| 90 |
|
|---|
| 91 | Float_t alt, az, dalt, daz, mjd;
|
|---|
| 92 | sscanf(str.Data(), "%f %f %*f %*f %*f %*f %f %f %f",
|
|---|
| 93 | &alt, &az, &dalt, &daz, &mjd);
|
|---|
| 94 |
|
|---|
| 95 | if (dalt==0/* || GetResidual(alt, az, alt+dalt, az+daz)>0.1*/)
|
|---|
| 96 | continue;
|
|---|
| 97 |
|
|---|
| 98 | //cout << dalt << " " << daz << " " << GetResidual(alt, az, alt+dalt, az+daz) << endl;
|
|---|
| 99 |
|
|---|
| 100 | mjd -= 53140.097505;
|
|---|
| 101 |
|
|---|
| 102 | Double_t res = GetResidual(alt, az, alt-dalt, az-daz);
|
|---|
| 103 |
|
|---|
| 104 | g[0].SetPoint(g[0].GetN(), g[0].GetN(), fabs(dalt));
|
|---|
| 105 | g[1].SetPoint(g[1].GetN(), g[1].GetN(), fabs(daz));
|
|---|
| 106 | g[2].SetPoint(g[2].GetN(), g[2].GetN(), res);
|
|---|
| 107 |
|
|---|
| 108 | h->Fill(res);
|
|---|
| 109 | }
|
|---|
| 110 |
|
|---|
| 111 | cout << "done (" << setw(3) << (Int_t)h->GetEntries() << "/";
|
|---|
| 112 | cout << setw(4) << g[0].GetN() << ") " << flush;
|
|---|
| 113 |
|
|---|
| 114 | return g[0].GetN();
|
|---|
| 115 | }
|
|---|
| 116 |
|
|---|
| 117 | struct Description_t
|
|---|
| 118 | {
|
|---|
| 119 | const char *fName;
|
|---|
| 120 | const char *fTitle;
|
|---|
| 121 | const char *fFile;
|
|---|
| 122 | };
|
|---|
| 123 |
|
|---|
| 124 | const Int_t counts = 22;
|
|---|
| 125 | Description_t desc[counts] =
|
|---|
| 126 | {
|
|---|
| 127 | // No good pointing model has been applied yet
|
|---|
| 128 | // {"0404", "TPoints Residuals 4/2004" , "tpoint/tpoint0404.txt"},
|
|---|
| 129 | // {"0405", "TPoints Residuals 5/2004" , "tpoint/tpoint0405.txt"},
|
|---|
| 130 | // {"04081", "TPoints Residuals 8/2004-1" , "tpoint/tpoint0408-1.txt"},
|
|---|
| 131 | {"04082", "TPoints Residuals 8/2004-2" , "tpoint/tpoint0408-2.txt"},
|
|---|
| 132 | {"0409", "TPoints Residuals 9/2004" , "tpoint/tpoint0409.txt"},
|
|---|
| 133 | // Culmination tests
|
|---|
| 134 | // {"0410", "TPoints Residuals 10/2004" , "tpoint/tpoint0410.txt"},
|
|---|
| 135 | {"0411", "TPoints Residuals 11/2004" , "tpoint/tpoint0411.txt"},
|
|---|
| 136 | {"0412", "TPoints Residuals 12/2004" , "tpoint/tpoint0412.txt"},
|
|---|
| 137 | // Worse pointing due to realignment of the mirror
|
|---|
| 138 | {"0503", "TPoints Residuals 3/2005" , "tpoint/tpoint0503.txt"},
|
|---|
| 139 | {"0504", "TPoints Residuals 4/2005" , "tpoint/tpoint0504.txt"},
|
|---|
| 140 | {"05051", "TPoints Residuals 5/2005-1" , "tpoint/tpoint0505-1.txt"},
|
|---|
| 141 | // Mirror alignment has been fixed
|
|---|
| 142 | {"05052", "TPoints Residuals 5/2005-2" , "tpoint/tpoint0505-2.txt"},
|
|---|
| 143 | // Fixes to pointing model due to fixing a screw
|
|---|
| 144 | {"0506", "TPoints Residuals 6/2005" , "tpoint/tpoint0506.txt"},
|
|---|
| 145 |
|
|---|
| 146 | {"0508", "TPoints Residuals 8/2005" , "tpoint/tpoint0508.txt"},
|
|---|
| 147 | {"0509", "TPoints Residuals 9/2005" , "tpoint/tpoint0509.txt"},
|
|---|
| 148 | // Quick-and-dirty mirror alignment
|
|---|
| 149 | {"0510", "TPoints Residuals 10/2005" , "tpoint/tpoint0510.txt"},
|
|---|
| 150 | // New mirror alignment after Tenerife meeting
|
|---|
| 151 | {"0511", "TPoints Residuals 11/2005" , "tpoint/tpoint0511.txt"},
|
|---|
| 152 |
|
|---|
| 153 | {"0512", "TPoints Residuals 12/2005" , "tpoint/tpoint0512.txt"},
|
|---|
| 154 | {"0601", "TPoints Residuals 1/2006" , "tpoint/tpoint0601.txt"},
|
|---|
| 155 | {"0603", "TPoints Residuals 3/2006" , "tpoint/tpoint0603.txt"},
|
|---|
| 156 | {"0604", "TPoints Residuals 4/2006" , "tpoint/tpoint0604.txt"},
|
|---|
| 157 | //{"0605", "TPoints Residuals 5/2006" , "tpoint/tpoint0605.txt"},
|
|---|
| 158 | //{"0606", "TPoints Residuals 6/2006" , "tpoint/tpoint0606.txt"},
|
|---|
| 159 | {"0701", "TPoints Residuals 1/2007" , "tpoint/tpoint0701.txt"},
|
|---|
| 160 | {"0702", "TPoints Residuals 2/2007" , "tpoint/tpoint0702.txt"},
|
|---|
| 161 | {"0703", "TPoints Residuals 3/2007" , "tpoint/tpoint0703.txt"},
|
|---|
| 162 | {"0704", "TPoints Residuals 4/2007" , "tpoint/tpoint0704.txt"},
|
|---|
| 163 | {"0705", "TPoints Residuals 5/2007" , "tpoint/tpoint0705.txt"},
|
|---|
| 164 | };
|
|---|
| 165 |
|
|---|
| 166 | void plot()
|
|---|
| 167 | {
|
|---|
| 168 | TGraph g[3];
|
|---|
| 169 |
|
|---|
| 170 | MBinning bins(10, 0, 0.2);
|
|---|
| 171 |
|
|---|
| 172 | TArrayI n(counts);
|
|---|
| 173 |
|
|---|
| 174 | TH1F hx[counts];
|
|---|
| 175 | for (int i=0; i<counts; i++)
|
|---|
| 176 | {
|
|---|
| 177 | hx[i].SetNameTitle(desc[i].fName, desc[i].fTitle);
|
|---|
| 178 | hx[i].SetDirectory(0);
|
|---|
| 179 | bins.Apply(hx[i]);
|
|---|
| 180 |
|
|---|
| 181 | cout << setw(2) << i << ": " << flush;
|
|---|
| 182 | n[i] = fill(desc[i].fFile, g, &hx[i]);
|
|---|
| 183 | cout << " Mean: " << setw(5) << setprecision(2) << hx[i].GetMean() << " deg +/- " << hx[i].GetRMS()<<endl;
|
|---|
| 184 | }
|
|---|
| 185 |
|
|---|
| 186 | g[0].SetMarkerColor(kGreen);
|
|---|
| 187 | g[1].SetMarkerColor(kMagenta);
|
|---|
| 188 | g[2].SetMarkerColor(kBlack);
|
|---|
| 189 | //g[2].SetLineColor(kBlack);
|
|---|
| 190 | g[0].SetMarkerStyle(kFullDotMedium);
|
|---|
| 191 | g[1].SetMarkerStyle(kFullDotMedium);
|
|---|
| 192 | g[2].SetMarkerStyle(kFullDotLarge);
|
|---|
| 193 |
|
|---|
| 194 | // --------- First Canvas ----------
|
|---|
| 195 |
|
|---|
| 196 | new TCanvas("Vs Time", "");
|
|---|
| 197 |
|
|---|
| 198 | TObject *obj[4];
|
|---|
| 199 |
|
|---|
| 200 | for (int i=0; i<3; i++)
|
|---|
| 201 | {
|
|---|
| 202 | g[i].SetFillColor(kWhite);
|
|---|
| 203 | g[i].SetLineColor(kWhite);
|
|---|
| 204 | obj[i] = g[i].Clone();
|
|---|
| 205 | obj[i]->SetBit(kCanDelete);
|
|---|
| 206 | obj[i]->Draw(i==0?"AP":"P");
|
|---|
| 207 | }
|
|---|
| 208 |
|
|---|
| 209 | TLegend leg(0.905, 0.86, 0.99, 0.99);
|
|---|
| 210 | leg.AddEntry(obj[0], " \\Delta\\theta");
|
|---|
| 211 | leg.AddEntry(obj[1], " \\Delta\\phi");
|
|---|
| 212 | leg.AddEntry(obj[2], " \\Delta");
|
|---|
| 213 | leg.DrawClone()->SetBit(kCanDelete);
|
|---|
| 214 |
|
|---|
| 215 | TLine l;
|
|---|
| 216 | l.SetLineColor(kBlue);
|
|---|
| 217 | for (int i=0; i<n.GetSize(); i++)
|
|---|
| 218 | l.DrawLine(n[i], 0, n[i], 0.2);
|
|---|
| 219 |
|
|---|
| 220 | // --------- Second Canvas ----------
|
|---|
| 221 |
|
|---|
| 222 | new TCanvas("Distrib", "");
|
|---|
| 223 |
|
|---|
| 224 | Double_t max=0;
|
|---|
| 225 | for (int i=0; i<n.GetSize(); i++)
|
|---|
| 226 | {
|
|---|
| 227 | if (hx[i].GetEntries()==0)
|
|---|
| 228 | {
|
|---|
| 229 | cout << "Skip #" << i << endl;
|
|---|
| 230 | continue;
|
|---|
| 231 | }
|
|---|
| 232 |
|
|---|
| 233 | hx[i].Scale(1./hx[i].GetEntries());
|
|---|
| 234 | max = TMath::Max(max, hx[i].GetMaximum());
|
|---|
| 235 | }
|
|---|
| 236 | for (int i=0; i<n.GetSize(); i++)
|
|---|
| 237 | {
|
|---|
| 238 | hx[i].SetMaximum(max*1.05);
|
|---|
| 239 | if (i<6)
|
|---|
| 240 | hx[i].SetLineColor(kRed+i);
|
|---|
| 241 | else
|
|---|
| 242 | hx[i].SetMarkerStyle(kPlus+i-6);
|
|---|
| 243 |
|
|---|
| 244 |
|
|---|
| 245 | }
|
|---|
| 246 |
|
|---|
| 247 | for (int i=0; i<counts; i++)
|
|---|
| 248 | hx[i].DrawCopy(i==0?"LP":"LPsame");
|
|---|
| 249 |
|
|---|
| 250 | return;
|
|---|
| 251 |
|
|---|
| 252 | for (int i=0; i<n.GetSize(); i++)
|
|---|
| 253 | {
|
|---|
| 254 | cout << "Mean: " << Form("%.3f +- %.3f", hx[i].GetMean(), hx[i].GetRMS());
|
|---|
| 255 | cout << " (Overflows=" << hx[i].GetBinContent(hx[i].GetNbinsX()+1)*hx[i].GetEntries() << ")" << endl;
|
|---|
| 256 | }
|
|---|
| 257 | }
|
|---|
