| 1 | /* ======================================================================== *\
|
|---|
| 2 | !
|
|---|
| 3 | ! *
|
|---|
| 4 | ! * This file is part of MARS, the MAGIC Analysis and Reconstruction
|
|---|
| 5 | ! * Software. It is distributed to you in the hope that it can be a useful
|
|---|
| 6 | ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
|
|---|
| 7 | ! * It is distributed WITHOUT ANY WARRANTY.
|
|---|
| 8 | ! *
|
|---|
| 9 | ! * Permission to use, copy, modify and distribute this software and its
|
|---|
| 10 | ! * documentation for any purpose is hereby granted without fee,
|
|---|
| 11 | ! * provided that the above copyright notice appear in all copies and
|
|---|
| 12 | ! * that both that copyright notice and this permission notice appear
|
|---|
| 13 | ! * in supporting documentation. It is provided "as is" without express
|
|---|
| 14 | ! * or implied warranty.
|
|---|
| 15 | ! *
|
|---|
| 16 | !
|
|---|
| 17 | ! Author(s): Javier Lopez, 12/2003 <mailto:jlopez@ifae.es>
|
|---|
| 18 | ! Author(s): Alex Armada, 1/2004 <mailto:armada@ifae.es>
|
|---|
| 19 | !
|
|---|
| 20 | ! Copyright: MAGIC Software Development, 2000-2003
|
|---|
| 21 | !
|
|---|
| 22 | !
|
|---|
| 23 | \* ======================================================================== */
|
|---|
| 24 |
|
|---|
| 25 |
|
|---|
| 26 | //------------------------------------------------------------------------- //
|
|---|
| 27 | // //
|
|---|
| 28 | // This macro fits the dc signal of a star using a two dimension gaussian //
|
|---|
| 29 | // for each dc measurement. Then the values of parameters of the fit are //
|
|---|
| 30 | // stored in histograms and shown at the end of the macro. //
|
|---|
| 31 | // //
|
|---|
| 32 | // USAGE: //
|
|---|
| 33 | // It has two arguments, //
|
|---|
| 34 | // 1- The first one is the dc file with the tracked star //
|
|---|
| 35 | // 2- The second one is a continuos light file used to intercalibrate //
|
|---|
| 36 | // the gain of the photomultipliers. //
|
|---|
| 37 | // (It's possible not to use the calibration file and then the gain //
|
|---|
| 38 | // of the pmts are supouse to be the same for all of them. //
|
|---|
| 39 | // 3- The third argument is just the number of dc measurements you want //
|
|---|
| 40 | // analize. If you put 0 it just stops after each fit and show you //
|
|---|
| 41 | // results. //
|
|---|
| 42 | // //
|
|---|
| 43 | //--------------------------------------------------------------------------//
|
|---|
| 44 |
|
|---|
| 45 | const Int_t numPixels = 577;
|
|---|
| 46 | Int_t nPixelsFitted;
|
|---|
| 47 | Bool_t isPixelsFitted[numPixels];
|
|---|
| 48 | Float_t z[numPixels],x[numPixels],y[numPixels],errorz[numPixels];
|
|---|
| 49 | Float_t chisquare;
|
|---|
| 50 |
|
|---|
| 51 | //______________________________________________________________________________
|
|---|
| 52 | //
|
|---|
| 53 | // Function used by Minuit to do the fit
|
|---|
| 54 | //
|
|---|
| 55 | void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
|
|---|
| 56 | {
|
|---|
| 57 | Int_t i;
|
|---|
| 58 |
|
|---|
| 59 | //calculate chisquare
|
|---|
| 60 | Double_t chisq = 0;
|
|---|
| 61 | Double_t delta;
|
|---|
| 62 | nPixelsFitted=0;
|
|---|
| 63 | for (i=1;i<numPixels; i++) {
|
|---|
| 64 | if (isPixelsFitted[i])
|
|---|
| 65 | {
|
|---|
| 66 | if (errorz[i] != 0.0)
|
|---|
| 67 | {
|
|---|
| 68 | delta = (z[i]-func(x[i],y[i],par))/errorz[i];
|
|---|
| 69 | chisq += delta*delta;
|
|---|
| 70 | nPixelsFitted++;
|
|---|
| 71 | }
|
|---|
| 72 | else
|
|---|
| 73 | cout << "This should never happen errorz[" << i << "] " << errorz[i] << endl;
|
|---|
| 74 | }
|
|---|
| 75 | }
|
|---|
| 76 | f = chisq;
|
|---|
| 77 | chisquare = chisq;
|
|---|
| 78 | }
|
|---|
| 79 |
|
|---|
| 80 | //______________________________________________________________________________
|
|---|
| 81 | //
|
|---|
| 82 | // The 2D gaussian fucntion used to fit the spot of the star
|
|---|
| 83 | //
|
|---|
| 84 | Double_t func(float x,float y,Double_t *par)
|
|---|
| 85 | {
|
|---|
| 86 | Double_t value=par[0]*TMath::exp(-(x-par[1])*(x-par[1])/(2*par[2]*par[2]))*TMath::exp(-(y-par[3])*(y-par[3])/(2*par[4]*par[4]));
|
|---|
| 87 | return value;
|
|---|
| 88 | }
|
|---|
| 89 |
|
|---|
| 90 | Bool_t HandleInput()
|
|---|
| 91 | {
|
|---|
| 92 | TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
|
|---|
| 93 |
|
|---|
| 94 | while (1)
|
|---|
| 95 | {
|
|---|
| 96 | //
|
|---|
| 97 | // While reading the input process gui events asynchronously
|
|---|
| 98 | //
|
|---|
| 99 | timer.TurnOn();
|
|---|
| 100 | TString input = Getline("Type 'q' to exit, <return> to go on: ");
|
|---|
| 101 | timer.TurnOff();
|
|---|
| 102 |
|
|---|
| 103 | if (input=="q\n")
|
|---|
| 104 | return kFALSE;
|
|---|
| 105 |
|
|---|
| 106 | if (input=="\n")
|
|---|
| 107 | return kTRUE;
|
|---|
| 108 | };
|
|---|
| 109 |
|
|---|
| 110 | return kFALSE;
|
|---|
| 111 | }
|
|---|
| 112 |
|
|---|
| 113 |
|
|---|
| 114 | void pointspreadfunction(TString fname, TString clname = "NULL", Int_t userNumLines = 1000)
|
|---|
| 115 | {
|
|---|
| 116 |
|
|---|
| 117 |
|
|---|
| 118 | Float_t currmean[numPixels];
|
|---|
| 119 | Float_t currsquaremean[numPixels];
|
|---|
| 120 | Float_t currrms[numPixels];
|
|---|
| 121 | Float_t meanofcurrmean = 0;
|
|---|
| 122 |
|
|---|
| 123 | for (Int_t swpix=0; swpix<numPixels; swpix++)
|
|---|
| 124 | {
|
|---|
| 125 | currmean[swpix] = 0.;
|
|---|
| 126 | currsquaremean[swpix] = 0.;
|
|---|
| 127 | currrms[swpix] = 0.;
|
|---|
| 128 | }
|
|---|
| 129 |
|
|---|
| 130 | Int_t numLines=0;
|
|---|
| 131 |
|
|---|
| 132 | //containers
|
|---|
| 133 | MGeomCamMagic geomcam;
|
|---|
| 134 | MCameraDC curr;
|
|---|
| 135 | MCameraDC currbis;
|
|---|
| 136 |
|
|---|
| 137 | const Float_t conv4mm2deg = geomcam.GetConvMm2Deg();
|
|---|
| 138 |
|
|---|
| 139 |
|
|---|
| 140 | if (clname != "NULL")
|
|---|
| 141 | {
|
|---|
| 142 |
|
|---|
| 143 | //
|
|---|
| 144 | // First run over continuos light files to have a DC calibration
|
|---|
| 145 | //
|
|---|
| 146 |
|
|---|
| 147 | MParList plist0;
|
|---|
| 148 |
|
|---|
| 149 | MTaskList tlist0;
|
|---|
| 150 | plist0.AddToList(&tlist0);
|
|---|
| 151 |
|
|---|
| 152 | plist0.AddToList(&geomcam);
|
|---|
| 153 | plist0.AddToList(&curr);
|
|---|
| 154 |
|
|---|
| 155 | MReportFileRead read0(clname);
|
|---|
| 156 | read0.SetHasNoHeader();
|
|---|
| 157 | read0.AddToList("MReportCurrents");
|
|---|
| 158 |
|
|---|
| 159 | tlist0.AddToList(&read0);
|
|---|
| 160 |
|
|---|
| 161 | MEvtLoop evtloop0;
|
|---|
| 162 | evtloop0.SetParList(&plist0);
|
|---|
| 163 |
|
|---|
| 164 |
|
|---|
| 165 | if (!evtloop0.PreProcess())
|
|---|
| 166 | return;
|
|---|
| 167 |
|
|---|
| 168 | while (tlist0.Process())
|
|---|
| 169 | {
|
|---|
| 170 | for (Int_t swpix=0; swpix<numPixels; swpix++)
|
|---|
| 171 | {
|
|---|
| 172 | meanofcurrmean += curr[swpix];
|
|---|
| 173 | currmean[swpix] += curr[swpix];
|
|---|
| 174 | currsquaremean[swpix] += curr[swpix]*curr[swpix];
|
|---|
| 175 | }
|
|---|
| 176 | numLines++;
|
|---|
| 177 | }
|
|---|
| 178 |
|
|---|
| 179 | evtloop0.PostProcess();
|
|---|
| 180 |
|
|---|
| 181 | meanofcurrmean /= (numLines*numPixels);
|
|---|
| 182 | for (Int_t swpix=0; swpix<numPixels; swpix++)
|
|---|
| 183 | {
|
|---|
| 184 |
|
|---|
| 185 | currmean[swpix] /= numLines;
|
|---|
| 186 | currsquaremean[swpix] /= numLines;
|
|---|
| 187 | currrms[swpix] = sqrt(fabs(currsquaremean[swpix] - currmean[swpix]*currmean[swpix]));
|
|---|
| 188 |
|
|---|
| 189 | curr[swpix] = currmean[swpix];
|
|---|
| 190 | currbis[swpix] = currrms[swpix];
|
|---|
| 191 |
|
|---|
| 192 | currmean[swpix] /= meanofcurrmean;
|
|---|
| 193 | currrms[swpix] /= meanofcurrmean;
|
|---|
| 194 |
|
|---|
| 195 | }
|
|---|
| 196 |
|
|---|
| 197 |
|
|---|
| 198 | /* MHCamera display0(geomcam);
|
|---|
| 199 | display0.SetPrettyPalette();
|
|---|
| 200 | display0.Draw();
|
|---|
| 201 |
|
|---|
| 202 | // curr.Print();
|
|---|
| 203 | display0.SetCamContent(currbis);
|
|---|
| 204 | cout << "PSF>> DC mean values drawn" << endl;
|
|---|
| 205 | // Remove the comments if you want to go through the file
|
|---|
| 206 | // event-by-event:
|
|---|
| 207 | if (!HandleInput())
|
|---|
| 208 | break;
|
|---|
| 209 | */
|
|---|
| 210 | }
|
|---|
| 211 | else
|
|---|
| 212 | {
|
|---|
| 213 | // If you don't use the continuous light this is the currrms[] array
|
|---|
| 214 | // is the error associated to the currents TMinuit will use.
|
|---|
| 215 | for (Int_t swpix=0; swpix<numPixels; swpix++)
|
|---|
| 216 | {
|
|---|
| 217 | currmean[swpix] = 1.;
|
|---|
| 218 | currrms[swpix] = 0.2;
|
|---|
| 219 | }
|
|---|
| 220 |
|
|---|
| 221 | }
|
|---|
| 222 |
|
|---|
| 223 | //
|
|---|
| 224 | // Now we can run over the dc data to extract the psf.
|
|---|
| 225 | //
|
|---|
| 226 | const Int_t maxNumLines = 10000;
|
|---|
| 227 |
|
|---|
| 228 | Double_t ux[maxNumLines];
|
|---|
| 229 | Double_t uy[maxNumLines];
|
|---|
| 230 | Double_t sx[maxNumLines];
|
|---|
| 231 | Double_t sy[maxNumLines];
|
|---|
| 232 | Double_t chisqu[maxNumLines];
|
|---|
| 233 | Double_t time[maxNumLines];
|
|---|
| 234 |
|
|---|
| 235 | MParList plist;
|
|---|
| 236 |
|
|---|
| 237 | MGeomCamMagic geomcam;
|
|---|
| 238 | MCameraDC curr;
|
|---|
| 239 | MTaskList tlist;
|
|---|
| 240 |
|
|---|
| 241 | plist.AddToList(&geomcam);
|
|---|
| 242 | plist.AddToList(&curr);
|
|---|
| 243 | plist.AddToList(&tlist);
|
|---|
| 244 |
|
|---|
| 245 | MReportFileRead read(fname);
|
|---|
| 246 | read.SetHasNoHeader();
|
|---|
| 247 | read.AddToList("MReportCurrents");
|
|---|
| 248 |
|
|---|
| 249 | tlist.AddToList(&read);
|
|---|
| 250 |
|
|---|
| 251 | MEvtLoop evtloop;
|
|---|
| 252 | evtloop.SetParList(&plist);
|
|---|
| 253 |
|
|---|
| 254 | if (!evtloop.PreProcess())
|
|---|
| 255 | return;
|
|---|
| 256 |
|
|---|
| 257 | MHCamera display(geomcam);
|
|---|
| 258 | display.SetPrettyPalette();
|
|---|
| 259 | display.Draw();
|
|---|
| 260 | gPad->SetLogy();
|
|---|
| 261 | gPad->cd(1);
|
|---|
| 262 |
|
|---|
| 263 | Double_t amin,edm,errdef;
|
|---|
| 264 | Int_t nvpar,nparx,icstat;
|
|---|
| 265 |
|
|---|
| 266 | Double_t max,maxerror;
|
|---|
| 267 | Double_t xmean,xsigma,ymean,ysigma;
|
|---|
| 268 | Double_t xmeanerror,xsigmaerror,ymeanerror,ysigmaerror;
|
|---|
| 269 |
|
|---|
| 270 | TEllipse ellipse;
|
|---|
| 271 | ellipse.SetFillStyle(4000);
|
|---|
| 272 | ellipse.SetLineWidth(2);
|
|---|
| 273 | ellipse.SetLineColor(2);
|
|---|
| 274 |
|
|---|
| 275 | ellipse.Draw();
|
|---|
| 276 |
|
|---|
| 277 | Int_t nbinsxy = 80;
|
|---|
| 278 | Float_t minxy = -600*conv4mm2deg;
|
|---|
| 279 | Float_t maxxy = 600*conv4mm2deg;
|
|---|
| 280 | Float_t fromdegtobin = (maxxy-minxy)/nbinsxy;
|
|---|
| 281 |
|
|---|
| 282 | TH2D psfhist("psfhist","",nbinsxy,minxy,maxxy,nbinsxy,minxy,maxxy);
|
|---|
| 283 | psfhist->GetXaxis()->SetTitle("[deg]");
|
|---|
| 284 | psfhist->GetYaxis()->SetTitle("[deg]");
|
|---|
| 285 | psfhist->GetZaxis()->SetTitle("DC [uA]");
|
|---|
| 286 |
|
|---|
| 287 | TCanvas *psfcanvas = new TCanvas("psfcanvas","Point Spread Funtion 2D",200,20,900,700);
|
|---|
| 288 |
|
|---|
| 289 |
|
|---|
| 290 | //
|
|---|
| 291 | // Using the first dc measurement we search the pixels which contains the star and define
|
|---|
| 292 | // an area to be fitted by Minuit which is 3 rings of neightbords around the peak of the star.
|
|---|
| 293 | //
|
|---|
| 294 |
|
|---|
| 295 | tlist.Process();
|
|---|
| 296 |
|
|---|
| 297 | Int_t numLines=0;
|
|---|
| 298 | Float_t minDCStar = 6.0;
|
|---|
| 299 |
|
|---|
| 300 | Int_t numPixelsInStar = 0;
|
|---|
| 301 | Float_t maxDC = 0;
|
|---|
| 302 | Int_t swpixelmaxDC;
|
|---|
| 303 |
|
|---|
| 304 | Bool_t isPixelsFittedTmp[numPixels];
|
|---|
| 305 |
|
|---|
| 306 | for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
|
|---|
| 307 | isPixelsFittedTmp[swpixel] = kFALSE;
|
|---|
| 308 |
|
|---|
| 309 | for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
|
|---|
| 310 | {
|
|---|
| 311 | if(curr[swpixel] > maxDC)
|
|---|
| 312 | {
|
|---|
| 313 | swpixelmaxDC = swpixel;
|
|---|
| 314 | maxDC = curr[swpixelmaxDC];
|
|---|
| 315 | }
|
|---|
| 316 |
|
|---|
| 317 | if(curr[swpixel]>minDCStar)
|
|---|
| 318 | {
|
|---|
| 319 | numPixelsInStar++;
|
|---|
| 320 | isPixelsFitted[swpixel] = kTRUE;
|
|---|
| 321 | }
|
|---|
| 322 | else
|
|---|
| 323 | isPixelsFitted[swpixel] = kFALSE;
|
|---|
| 324 | }
|
|---|
| 325 |
|
|---|
| 326 | if (numPixelsInStar == 0)
|
|---|
| 327 | {
|
|---|
| 328 | cout << "PSF>> Warning: none pixel over minDCStar(" << minDCStar << ')' << endl;
|
|---|
| 329 | return;
|
|---|
| 330 | }
|
|---|
| 331 |
|
|---|
| 332 | //1st neighboor ring
|
|---|
| 333 | for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
|
|---|
| 334 | if (isPixelsFitted[swpixel])
|
|---|
| 335 | for(Int_t next=0; next<geomcam[swpixel].GetNumNeighbors(); next++)
|
|---|
| 336 | isPixelsFittedTmp[geomcam[swpixel].GetNeighbor(next)] = kTRUE;
|
|---|
| 337 |
|
|---|
| 338 | for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
|
|---|
| 339 | if (isPixelsFittedTmp[swpixel])
|
|---|
| 340 | isPixelsFitted[swpixel] = kTRUE;
|
|---|
| 341 |
|
|---|
| 342 | //2on neighboor ring
|
|---|
| 343 | for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
|
|---|
| 344 | if (isPixelsFitted[swpixel])
|
|---|
| 345 | for(Int_t next=0; next<geomcam[swpixel].GetNumNeighbors(); next++)
|
|---|
| 346 | isPixelsFittedTmp[geomcam[swpixel].GetNeighbor(next)] = kTRUE;
|
|---|
| 347 |
|
|---|
| 348 | for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
|
|---|
| 349 | if (isPixelsFittedTmp[swpixel])
|
|---|
| 350 | isPixelsFitted[swpixel] = kTRUE;
|
|---|
| 351 |
|
|---|
| 352 |
|
|---|
| 353 | //3rt neighboor ring
|
|---|
| 354 | for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
|
|---|
| 355 | if (isPixelsFitted[swpixel])
|
|---|
| 356 | for(Int_t next=0; next<geomcam[swpixel].GetNumNeighbors(); next++)
|
|---|
| 357 | isPixelsFittedTmp[geomcam[swpixel].GetNeighbor(next)] = kTRUE;
|
|---|
| 358 |
|
|---|
| 359 | for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
|
|---|
| 360 | if (isPixelsFittedTmp[swpixel])
|
|---|
| 361 | isPixelsFitted[swpixel] = kTRUE;
|
|---|
| 362 |
|
|---|
| 363 |
|
|---|
| 364 | for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
|
|---|
| 365 | curr[swpixel] = (Float_t)isPixelsFitted[swpixel];
|
|---|
| 366 |
|
|---|
| 367 | /* MHCamera display0(geomcam);
|
|---|
| 368 | display0.SetPrettyPalette();
|
|---|
| 369 | display0.Draw();
|
|---|
| 370 |
|
|---|
| 371 | display0.SetCamContent(curr);
|
|---|
| 372 | cout << "PSF>> Fitted pixels drawn" << endl;
|
|---|
| 373 | // Remove the comments if you want to go through the file
|
|---|
| 374 | // event-by-event:
|
|---|
| 375 | if (!HandleInput())
|
|---|
| 376 | break;
|
|---|
| 377 | */
|
|---|
| 378 |
|
|---|
| 379 | // Minuit initialization
|
|---|
| 380 |
|
|---|
| 381 | TMinuit *gMinuit = new TMinuit(7); //initialize TMinuit with a maximum of 5 params
|
|---|
| 382 | gMinuit->SetFCN(fcn);
|
|---|
| 383 |
|
|---|
| 384 | Double_t arglist[10];
|
|---|
| 385 | Int_t ierflg = 0;
|
|---|
| 386 |
|
|---|
| 387 | arglist[0] = 1;
|
|---|
| 388 | gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
|
|---|
| 389 | // arglist[0] = -1;
|
|---|
| 390 | arglist[0] = 0;
|
|---|
| 391 | gMinuit->mnexcm("SET PRI", arglist ,1,ierflg);
|
|---|
| 392 |
|
|---|
| 393 | // Set starting values and step sizes for parameters
|
|---|
| 394 | Double_t vstart[5];
|
|---|
| 395 | Double_t step[5];
|
|---|
| 396 | Double_t lowlimit[5] = {minDCStar, -2., 0.05, -2, 0.05};
|
|---|
| 397 | Double_t uplimit[5] = {30., 2., 1., 2., 1.};
|
|---|
| 398 |
|
|---|
| 399 | vstart[0] = maxDC;
|
|---|
| 400 | vstart[1] = geomcam[swpixelmaxDC].GetX()*conv4mm2deg;
|
|---|
| 401 | vstart[2] = 30*TMath::Sqrt(numPixelsInStar/2)*conv4mm2deg;
|
|---|
| 402 | vstart[3] = geomcam[swpixelmaxDC].GetY()*conv4mm2deg;
|
|---|
| 403 | vstart[4] = 30*TMath::Sqrt(numPixelsInStar/2)*conv4mm2deg;
|
|---|
| 404 |
|
|---|
| 405 | for(Int_t i=0; i<5; i++)
|
|---|
| 406 | {
|
|---|
| 407 | if (vstart[i] != 0)
|
|---|
| 408 | step[i] = TMath::Abs(vstart[i]/sqrt(2));
|
|---|
| 409 | else
|
|---|
| 410 | step[i] = uplimit[i]/2;
|
|---|
| 411 | }
|
|---|
| 412 |
|
|---|
| 413 | gMinuit->mnparm(0, "max", vstart[0], step[0], lowlimit[0], uplimit[0],ierflg);
|
|---|
| 414 | gMinuit->mnparm(1, "xmean", vstart[1], step[1], lowlimit[1], uplimit[1],ierflg);
|
|---|
| 415 | gMinuit->mnparm(2, "xsigma", vstart[2], step[2], lowlimit[2], uplimit[2],ierflg);
|
|---|
| 416 | gMinuit->mnparm(3, "ymean", vstart[3], step[3], lowlimit[3], uplimit[3],ierflg);
|
|---|
| 417 | gMinuit->mnparm(4, "ysigma", vstart[4], step[4], lowlimit[4], uplimit[4],ierflg);
|
|---|
| 418 |
|
|---|
| 419 | while (tlist.Process() && numLines < maxNumLines)
|
|---|
| 420 | {
|
|---|
| 421 |
|
|---|
| 422 | for (Int_t swpixel=1; swpixel<577; swpixel++)
|
|---|
| 423 | {
|
|---|
| 424 |
|
|---|
| 425 | x[swpixel] = geomcam[swpixel].GetX()*conv4mm2deg;
|
|---|
| 426 | y[swpixel] = geomcam[swpixel].GetY()*conv4mm2deg;
|
|---|
| 427 | z[swpixel] = curr[swpixel]/currmean[swpixel];
|
|---|
| 428 | errorz[swpixel] = TMath::Sqrt((curr[swpixel]*currrms[swpixel]/(currmean[swpixel]*currmean[swpixel]))*(curr[swpixel]*currrms[swpixel]/(currmean[swpixel]*currmean[swpixel]))+(0.1)/currmean[swpixel]*(0.1)/currmean[swpixel]);
|
|---|
| 429 |
|
|---|
| 430 |
|
|---|
| 431 | psfhist->SetBinContent((Int_t)((x[swpixel]+600*conv4mm2deg)/fromdegtobin),(Int_t)((y[swpixel]+600*conv4mm2deg)/fromdegtobin),z[swpixel]);
|
|---|
| 432 | }
|
|---|
| 433 |
|
|---|
| 434 | psfcanvas->cd(1);
|
|---|
| 435 | psfhist->Draw("lego2");
|
|---|
| 436 |
|
|---|
| 437 | // Now ready for minimization step
|
|---|
| 438 | arglist[0] = 500;
|
|---|
| 439 | arglist[1] = 1.;
|
|---|
| 440 | gMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);
|
|---|
| 441 |
|
|---|
| 442 | // Print results
|
|---|
| 443 | /* gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
|
|---|
| 444 | gMinuit->mnprin(3,amin);
|
|---|
| 445 | */
|
|---|
| 446 | gMinuit->GetParameter(0,max,maxerror);
|
|---|
| 447 | gMinuit->GetParameter(1,xmean,xmeanerror);
|
|---|
| 448 | gMinuit->GetParameter(2,xsigma,xsigmaerror);
|
|---|
| 449 | gMinuit->GetParameter(3,ymean,ymeanerror);
|
|---|
| 450 | gMinuit->GetParameter(4,ysigma,ysigmaerror);
|
|---|
| 451 |
|
|---|
| 452 | /* cout << endl;
|
|---|
| 453 | cout << "numLine " << numLines << endl;
|
|---|
| 454 | cout << "max \t" << max << " +- " << maxerror << endl;
|
|---|
| 455 | cout << "xmean \t" << xmean << " +- " << xmeanerror << endl;
|
|---|
| 456 | cout << "xsigma \t" << TMath::Abs(xsigma) << " +- " << xsigmaerror << endl;
|
|---|
| 457 | cout << "ymean \t" << ymean << " +- " << ymeanerror << endl;
|
|---|
| 458 | cout << "ysigma \t" << TMath::Abs(ysigma) << " +- " << ysigmaerror << endl;
|
|---|
| 459 | cout << "chisquare/ndof \t" << chisquare/(nPixelsFitted-5) << endl;
|
|---|
| 460 | */
|
|---|
| 461 |
|
|---|
| 462 | chisqu[numLines]=chisquare/(nPixelsFitted-5);
|
|---|
| 463 |
|
|---|
| 464 | if(chisqu[numLines]<100.)
|
|---|
| 465 | {
|
|---|
| 466 | ux[numLines]=xmean;
|
|---|
| 467 | uy[numLines]=ymean;
|
|---|
| 468 | sx[numLines]=TMath::Abs(xsigma);
|
|---|
| 469 | sy[numLines]=TMath::Abs(ysigma);
|
|---|
| 470 | time[numLines]=numLines;
|
|---|
| 471 |
|
|---|
| 472 | display.SetCamContent(curr);
|
|---|
| 473 | gPad->cd(1);
|
|---|
| 474 | ellipse.SetX1(xmean/conv4mm2deg);
|
|---|
| 475 | ellipse.SetY1(ymean/conv4mm2deg);
|
|---|
| 476 | ellipse.SetR1(TMath::Abs(xsigma)/conv4mm2deg);
|
|---|
| 477 | ellipse.SetR2(TMath::Abs(ysigma)/conv4mm2deg);
|
|---|
| 478 |
|
|---|
| 479 | gPad->Modified();
|
|---|
| 480 | gPad->Update();
|
|---|
| 481 |
|
|---|
| 482 | // Remove the comments if you want to go through the file
|
|---|
| 483 | //event-by-event:
|
|---|
| 484 | if (userNumLines>0)
|
|---|
| 485 | {
|
|---|
| 486 | if (numLines>userNumLines)
|
|---|
| 487 | break;
|
|---|
| 488 | }
|
|---|
| 489 | else
|
|---|
| 490 | {
|
|---|
| 491 | if (!HandleInput())
|
|---|
| 492 | break;
|
|---|
| 493 | }
|
|---|
| 494 | numLines++;
|
|---|
| 495 | }
|
|---|
| 496 | }
|
|---|
| 497 |
|
|---|
| 498 |
|
|---|
| 499 | evtloop.PostProcess();
|
|---|
| 500 |
|
|---|
| 501 | //
|
|---|
| 502 | // Draw the ditributions of the sigmas the point spread function
|
|---|
| 503 | //
|
|---|
| 504 |
|
|---|
| 505 | cout<<"PSF>> Number of lines "<<numLines<<endl;
|
|---|
| 506 |
|
|---|
| 507 | gROOT->Reset();
|
|---|
| 508 | gStyle->SetCanvasColor(0);
|
|---|
| 509 | gStyle->SetCanvasBorderMode(0);
|
|---|
| 510 | gStyle->SetPadBorderMode(0);
|
|---|
| 511 | gStyle->SetFrameBorderMode(0);
|
|---|
| 512 | gStyle->SetOptStat(00000000);
|
|---|
| 513 |
|
|---|
| 514 | //
|
|---|
| 515 | // Find in the file name the date, run and project name to put it in the title
|
|---|
| 516 | //
|
|---|
| 517 |
|
|---|
| 518 | Size_t pos = fname.Last('/');
|
|---|
| 519 | TString iRun = TString(fname(pos+24,5));
|
|---|
| 520 | TString iYear = TString(fname(pos+4,4));
|
|---|
| 521 | TString iMonth = TString(fname(pos+9,2));
|
|---|
| 522 | TString iDay = TString(fname(pos+12,2));
|
|---|
| 523 |
|
|---|
| 524 | TString iHour = TString(fname(pos+15,2));
|
|---|
| 525 | TString iMin = TString(fname(pos+18,2));
|
|---|
| 526 | TString iSec = TString(fname(pos+21,2));
|
|---|
| 527 |
|
|---|
| 528 | Size_t poslast = fname.Last('.');
|
|---|
| 529 | Size_t posfirst = poslast-1;
|
|---|
| 530 | while (fname[posfirst] != '_')
|
|---|
| 531 | posfirst--;
|
|---|
| 532 |
|
|---|
| 533 | TString iSource = TString(fname(posfirst+1,poslast-posfirst-1));
|
|---|
| 534 |
|
|---|
| 535 |
|
|---|
| 536 | char str[100];
|
|---|
| 537 |
|
|---|
| 538 | // sprintf(str,"Date %s/%s/%s Run %s Source %s",iYear.Data(),iMonth.Data(),iDay.Data(),iRun.Data(),iSource.Data());
|
|---|
| 539 | sprintf(str,"Date %s/%s/%s StartTime %s:%s:%s Run %s Source %s",iYear.Data(),iMonth.Data(),iDay.Data(),iHour.Data(),iMin.Data(),iSec.Data(),iRun.Data(),iSource.Data());
|
|---|
| 540 |
|
|---|
| 541 | c1 = new TCanvas("c1",str,0,0,1200,850);
|
|---|
| 542 | // c1 = new TCanvas("c1","Time evolution & distributions",0,0,1200,850);
|
|---|
| 543 | c1->Divide(3,2);
|
|---|
| 544 |
|
|---|
| 545 | c1->cd(1);
|
|---|
| 546 |
|
|---|
| 547 | TMath math;
|
|---|
| 548 |
|
|---|
| 549 | Double_t minmeanx, maxmeanx;
|
|---|
| 550 | minmeanx = ux[math.LocMin(numLines,ux)];
|
|---|
| 551 | maxmeanx = ux[math.LocMax(numLines,ux)];
|
|---|
| 552 |
|
|---|
| 553 | Double_t minmeany, maxmeany;
|
|---|
| 554 | minmeany = uy[math.LocMin(numLines,uy)];
|
|---|
| 555 | maxmeany = uy[math.LocMax(numLines,uy)];
|
|---|
| 556 |
|
|---|
| 557 | Double_t minmean, maxmean;
|
|---|
| 558 | minmean = math.Min(minmeanx,minmeany);
|
|---|
| 559 | maxmean = math.Max(maxmeanx,maxmeany);
|
|---|
| 560 |
|
|---|
| 561 | Double_t diff;
|
|---|
| 562 | diff = maxmean - minmean;
|
|---|
| 563 | diff = 0.1*diff;
|
|---|
| 564 | minmean = minmean - diff;
|
|---|
| 565 | maxmean = maxmean + diff;
|
|---|
| 566 |
|
|---|
| 567 | Double_t mintime, maxtime;
|
|---|
| 568 | mintime = time[math.LocMin(numLines,time)];
|
|---|
| 569 | maxtime = time[math.LocMax(numLines,time)];
|
|---|
| 570 |
|
|---|
| 571 | TH2D *h1 = new TH2D("h1","",1,mintime-1,maxtime+1,1,minmean,maxmean);
|
|---|
| 572 | h1->GetXaxis()->SetTitle("Event number");
|
|---|
| 573 | h1->GetYaxis()->SetTitle("mean position (deg)");
|
|---|
| 574 | h1->Draw();
|
|---|
| 575 |
|
|---|
| 576 | TGraph *grtimeevolmeanx = new TGraph(numLines,time,ux);
|
|---|
| 577 | grtimeevolmeanx->SetMarkerColor(3);
|
|---|
| 578 | grtimeevolmeanx->SetMarkerStyle(20);
|
|---|
| 579 | grtimeevolmeanx->SetMarkerSize (0.4);
|
|---|
| 580 | grtimeevolmeanx->Draw("P");
|
|---|
| 581 |
|
|---|
| 582 | TGraph *grtimeevolmeany = new TGraph(numLines,time,uy);
|
|---|
| 583 | grtimeevolmeany->SetMarkerColor(6);
|
|---|
| 584 | grtimeevolmeany->SetMarkerStyle(24);
|
|---|
| 585 | grtimeevolmeany->SetMarkerSize (0.4);
|
|---|
| 586 | grtimeevolmeany->Draw("P");
|
|---|
| 587 |
|
|---|
| 588 | legmeanxy = new TLegend(0.8,0.85,0.95,0.95);
|
|---|
| 589 | legmeanxy.SetTextSize(0.03);
|
|---|
| 590 | legmeanxy.AddEntry(grtimeevolmeanx,"mean x","P");
|
|---|
| 591 | legmeanxy.AddEntry(grtimeevolmeany,"mean y","P");
|
|---|
| 592 | legmeanxy.Draw();
|
|---|
| 593 |
|
|---|
| 594 | c1->cd(2);
|
|---|
| 595 |
|
|---|
| 596 | TMath math;
|
|---|
| 597 |
|
|---|
| 598 | Double_t minsigmax, maxsigmax;
|
|---|
| 599 | minsigmax = sx[math.LocMin(numLines,sx)];
|
|---|
| 600 | maxsigmax = sx[math.LocMax(numLines,sx)];
|
|---|
| 601 |
|
|---|
| 602 | Double_t minsigmay, maxsigmay;
|
|---|
| 603 | minsigmay = sy[math.LocMin(numLines,sy)];
|
|---|
| 604 | maxsigmay = sy[math.LocMax(numLines,sy)];
|
|---|
| 605 |
|
|---|
| 606 | Double_t minsigma, maxsigma;
|
|---|
| 607 | minsigma = math.Min(minsigmax,minsigmay);
|
|---|
| 608 | maxsigma = math.Max(maxsigmax,maxsigmay);
|
|---|
| 609 |
|
|---|
| 610 | diff = maxsigma - minsigma;
|
|---|
| 611 | diff = 0.1*diff;
|
|---|
| 612 | minsigma = minsigma - diff;
|
|---|
| 613 | maxsigma = maxsigma + diff;
|
|---|
| 614 |
|
|---|
| 615 | TH2D *h2 = new TH2D("h2","",1,mintime-1,maxtime+1,1,minsigma,maxsigma);
|
|---|
| 616 | h2->GetXaxis()->SetTitle("Event number");
|
|---|
| 617 | h2->GetYaxis()->SetTitle("PSF Rms (deg)");
|
|---|
| 618 | h2->Draw();
|
|---|
| 619 |
|
|---|
| 620 | TGraph* grtimeevolsigmax= new TGraph(numLines,time,sx);
|
|---|
| 621 | grtimeevolsigmax->SetMarkerColor(3);
|
|---|
| 622 | grtimeevolsigmax->SetMarkerStyle(20);
|
|---|
| 623 | grtimeevolsigmax->SetMarkerSize (0.4);
|
|---|
| 624 | grtimeevolsigmax->Draw("P");
|
|---|
| 625 |
|
|---|
| 626 | TGraph* grtimeevolsigmay= new TGraph(numLines,time,sy);
|
|---|
| 627 | grtimeevolsigmay->SetMarkerColor(6);
|
|---|
| 628 | grtimeevolsigmay->SetMarkerStyle(24);
|
|---|
| 629 | grtimeevolsigmay->SetMarkerSize (0.4);
|
|---|
| 630 | grtimeevolsigmay->Draw("P");
|
|---|
| 631 |
|
|---|
| 632 | legsigmaxy = new TLegend(0.8,0.85,0.95,0.95);
|
|---|
| 633 | legsigmaxy.SetTextSize(0.03);
|
|---|
| 634 | legsigmaxy.AddEntry(grtimeevolsigmax,"sigma x","P");
|
|---|
| 635 | legsigmaxy.AddEntry(grtimeevolsigmay,"sigma y","P");
|
|---|
| 636 | legsigmaxy.Draw();
|
|---|
| 637 |
|
|---|
| 638 | c1->cd(3);
|
|---|
| 639 |
|
|---|
| 640 | Double_t minchisqu, maxchisqu;
|
|---|
| 641 |
|
|---|
| 642 | minchisqu = chisqu[math.LocMin(numLines,chisqu)];
|
|---|
| 643 | maxchisqu = chisqu[math.LocMax(numLines,chisqu)];
|
|---|
| 644 |
|
|---|
| 645 | diff = maxchisqu - minchisqu;
|
|---|
| 646 | diff = 0.1*diff;
|
|---|
| 647 | minchisqu = minchisqu - diff;
|
|---|
| 648 | maxchisqu = maxchisqu + diff;
|
|---|
| 649 |
|
|---|
| 650 | TH2D *h3 = new TH2D("h3","",1,mintime-1,maxtime+1,1,minchisqu,maxchisqu);
|
|---|
| 651 | h3->GetXaxis()->SetTitle("Event number");
|
|---|
| 652 | h3->Draw();
|
|---|
| 653 |
|
|---|
| 654 | TGraph * grtimeevolchisqu = new TGraph(numLines,time,chisqu);
|
|---|
| 655 | grtimeevolchisqu->SetMarkerColor(215);
|
|---|
| 656 | grtimeevolchisqu->SetMarkerStyle(20);
|
|---|
| 657 | grtimeevolchisqu->SetMarkerSize(0.4);
|
|---|
| 658 | grtimeevolchisqu->Draw("P");
|
|---|
| 659 |
|
|---|
| 660 | legchisqu = new TLegend(0.55,0.90,0.95,0.95);
|
|---|
| 661 | legchisqu.SetTextSize(0.03);
|
|---|
| 662 | legchisqu.AddEntry(grtimeevolchisqu,"chi square / ndof","P");
|
|---|
| 663 | legchisqu.Draw();
|
|---|
| 664 |
|
|---|
| 665 | //***************************************
|
|---|
| 666 |
|
|---|
| 667 | const Int_t nbins = 100;
|
|---|
| 668 |
|
|---|
| 669 | TH1D *xsigmahist = new TH1D("xsigmahist","",nbins,minsigma,maxsigma);
|
|---|
| 670 | TH1D *ysigmahist = new TH1D("ysigmahist","",nbins,minsigma,maxsigma);
|
|---|
| 671 | TH1D *xmeanhist = new TH1D("xmeanhist","",nbins,minmean,maxmean);
|
|---|
| 672 | TH1D *ymeanhist = new TH1D("ymeanhist","",nbins,minmean,maxmean);
|
|---|
| 673 | TH1D *chisquhist = new TH1D("chisquhist","",nbins,minchisqu,maxchisqu);
|
|---|
| 674 |
|
|---|
| 675 | for (Int_t i=0; i<numLines; i++)
|
|---|
| 676 | {
|
|---|
| 677 | xsigmahist->Fill(TMath::Abs(sx[i]));
|
|---|
| 678 | ysigmahist->Fill(TMath::Abs(sy[i]));
|
|---|
| 679 | xmeanhist->Fill(ux[i]);
|
|---|
| 680 | ymeanhist->Fill(uy[i]);
|
|---|
| 681 | chisquhist->Fill(chisqu[i]);
|
|---|
| 682 |
|
|---|
| 683 | }
|
|---|
| 684 |
|
|---|
| 685 | c1->cd(5);
|
|---|
| 686 |
|
|---|
| 687 | TMath math;
|
|---|
| 688 | Double_t maxsigma;
|
|---|
| 689 | Int_t binmaxx, binmaxy;
|
|---|
| 690 | xsigmahist->SetLineColor(3);
|
|---|
| 691 | xsigmahist->SetLineWidth(2);
|
|---|
| 692 | xsigmahist->SetXTitle("RMS [deg]");
|
|---|
| 693 | binmaxx = xsigmahist->GetMaximumBin();
|
|---|
| 694 | binmaxx = xsigmahist->GetBinContent(binmaxx);
|
|---|
| 695 |
|
|---|
| 696 | ysigmahist->SetLineColor(6);
|
|---|
| 697 | ysigmahist->SetLineWidth(2);
|
|---|
| 698 | binmaxy = ysigmahist->GetMaximumBin();
|
|---|
| 699 | binmaxy = ysigmahist->GetBinContent(binmaxy);
|
|---|
| 700 |
|
|---|
| 701 | maxsigma = math.Max(binmaxx,binmaxy);
|
|---|
| 702 | maxsigma += 1;
|
|---|
| 703 |
|
|---|
| 704 | xsigmahist->SetMaximum(maxsigma);
|
|---|
| 705 | ysigmahist->SetMaximum(maxsigma);
|
|---|
| 706 | xsigmahist->DrawCopy();
|
|---|
| 707 | ysigmahist->DrawCopy("Same");
|
|---|
| 708 |
|
|---|
| 709 | TLegend *legendsigma = new TLegend(.3,.8,.95,.95);
|
|---|
| 710 | legendsigma->SetTextSize(0.03);
|
|---|
| 711 | char xsigmatitle[100];
|
|---|
| 712 | char ysigmatitle[100];
|
|---|
| 713 | sprintf(xsigmatitle,"PSF Rms on X axis -- %5.2f +/- %5.2f deg",xsigmahist->GetMean(),xsigmahist->GetRMS());
|
|---|
| 714 | sprintf(ysigmatitle,"PSF Rms on Y axis -- %5.2f +/- %5.2f deg",ysigmahist->GetMean(),ysigmahist->GetRMS());
|
|---|
| 715 | legendsigma->AddEntry(xsigmahist,xsigmatitle,"F");
|
|---|
| 716 | legendsigma->AddEntry(ysigmahist,ysigmatitle,"F");
|
|---|
| 717 | legendsigma->Draw();
|
|---|
| 718 |
|
|---|
| 719 | c1->cd(4);
|
|---|
| 720 |
|
|---|
| 721 | Double_t maxmean;
|
|---|
| 722 |
|
|---|
| 723 | xmeanhist->SetLineColor(3);
|
|---|
| 724 | xmeanhist->SetLineWidth(2);
|
|---|
| 725 | xmeanhist->SetXTitle("mean [deg]");
|
|---|
| 726 | binmaxx = xmeanhist->GetMaximumBin();
|
|---|
| 727 | binmaxx = xmeanhist->GetBinContent(binmaxx);
|
|---|
| 728 |
|
|---|
| 729 | ymeanhist->SetLineColor(6);
|
|---|
| 730 | ymeanhist->SetLineWidth(2);
|
|---|
| 731 | binmaxy = ymeanhist->GetMaximumBin();
|
|---|
| 732 | binmaxy = ymeanhist->GetBinContent(binmaxy);
|
|---|
| 733 |
|
|---|
| 734 | maxmean = math.Max(binmaxx,binmaxy);
|
|---|
| 735 | maxmean += 1;
|
|---|
| 736 |
|
|---|
| 737 | xmeanhist->SetMaximum(maxmean);
|
|---|
| 738 | ymeanhist->SetMaximum(maxmean);
|
|---|
| 739 | xmeanhist->DrawCopy();
|
|---|
| 740 | ymeanhist->DrawCopy("Same");
|
|---|
| 741 |
|
|---|
| 742 | TLegend *legendmean = new TLegend(.35,.8,.95,.95);
|
|---|
| 743 | legendmean->SetTextSize(0.03);
|
|---|
| 744 | char xmeantitle[100];
|
|---|
| 745 | char ymeantitle[100];
|
|---|
| 746 | sprintf(xmeantitle,"mean on X axis -- %5.2f +/- %5.2f deg",xmeanhist->GetMean(),xmeanhist->GetRMS());
|
|---|
| 747 | sprintf(ymeantitle,"mean on Y axis -- %5.2f +/- %5.2f deg",ymeanhist->GetMean(),ymeanhist->GetRMS());
|
|---|
| 748 | legendmean->AddEntry(xmeanhist,xmeantitle,"F");
|
|---|
| 749 | legendmean->AddEntry(ymeanhist,ymeantitle,"F");
|
|---|
| 750 | legendmean->Draw();
|
|---|
| 751 |
|
|---|
| 752 | //meancanvas->Modified();
|
|---|
| 753 | //meancanvas->Update();
|
|---|
| 754 |
|
|---|
| 755 | c1->cd(6);
|
|---|
| 756 |
|
|---|
| 757 | chisquhist->SetLineColor(215);
|
|---|
| 758 | chisquhist->SetLineWidth(2);
|
|---|
| 759 | chisquhist->SetXTitle("chi square / ndof");
|
|---|
| 760 | TAxis * axis = chisquhist->GetXaxis();
|
|---|
| 761 | axis->SetLabelSize(0.025);
|
|---|
| 762 | chisquhist->DrawCopy();
|
|---|
| 763 |
|
|---|
| 764 | TLegend *legendchisqu = new TLegend(.4,.85,.95,.95);
|
|---|
| 765 | legendchisqu->SetTextSize(0.03);
|
|---|
| 766 | char chisqutitle[100];
|
|---|
| 767 | sprintf(chisqutitle,"chi square / ndof -- %5.2f +/- %5.2f ",chisquhist->GetMean(),chisquhist->GetRMS());
|
|---|
| 768 | legendchisqu->AddEntry(chisquhist,chisqutitle,"F");
|
|---|
| 769 | legendchisqu->Draw();
|
|---|
| 770 |
|
|---|
| 771 |
|
|---|
| 772 | return;
|
|---|
| 773 |
|
|---|
| 774 | }
|
|---|
| 775 |
|
|---|
| 776 |
|
|---|
| 777 |
|
|---|
| 778 |
|
|---|
| 779 |
|
|---|
| 780 |
|
|---|
| 781 |
|
|---|
| 782 |
|
|---|
| 783 |
|
|---|
| 784 |
|
|---|
| 785 |
|
|---|
| 786 |
|
|---|
| 787 |
|
|---|
| 788 |
|
|---|
| 789 |
|
|---|
| 790 |
|
|---|
| 791 |
|
|---|
| 792 |
|
|---|
| 793 |
|
|---|
| 794 |
|
|---|
| 795 |
|
|---|
| 796 |
|
|---|
| 797 |
|
|---|
| 798 |
|
|---|
| 799 |
|
|---|
| 800 |
|
|---|
| 801 |
|
|---|
| 802 |
|
|---|
| 803 |
|
|---|
| 804 |
|
|---|
| 805 |
|
|---|
| 806 |
|
|---|
| 807 |
|
|---|
| 808 |
|
|---|
| 809 |
|
|---|
| 810 |
|
|---|
| 811 |
|
|---|
| 812 |
|
|---|
| 813 |
|
|---|
| 814 |
|
|---|
| 815 |
|
|---|
| 816 |
|
|---|
| 817 |
|
|---|
| 818 |
|
|---|
| 819 |
|
|---|
| 820 |
|
|---|
| 821 |
|
|---|
| 822 |
|
|---|
| 823 |
|
|---|
| 824 |
|
|---|
| 825 |
|
|---|
| 826 |
|
|---|
| 827 |
|
|---|
| 828 |
|
|---|
| 829 |
|
|---|
| 830 |
|
|---|
| 831 |
|
|---|
| 832 |
|
|---|
| 833 |
|
|---|
| 834 |
|
|---|
| 835 |
|
|---|
| 836 |
|
|---|
| 837 |
|
|---|
| 838 |
|
|---|
| 839 |
|
|---|
| 840 |
|
|---|
| 841 |
|
|---|
| 842 |
|
|---|
| 843 |
|
|---|
| 844 |
|
|---|
| 845 |
|
|---|
| 846 |
|
|---|
| 847 |
|
|---|
| 848 |
|
|---|
| 849 |
|
|---|
| 850 |
|
|---|
| 851 |
|
|---|
| 852 |
|
|---|
| 853 |
|
|---|
