| 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 | // //
|
|---|
| 40 | //--------------------------------------------------------------------------//
|
|---|
| 41 |
|
|---|
| 42 | const Int_t numPixels = 577;
|
|---|
| 43 | Int_t nPixelsFitted; // 396[inners] 36[3rings] 60[4rings] 90[5rings]
|
|---|
| 44 | Bool_t isPixelsFitted[numPixels];
|
|---|
| 45 | Float_t z[numPixels],x[numPixels],y[numPixels],errorz[numPixels];
|
|---|
| 46 | Float_t chisquare;
|
|---|
| 47 |
|
|---|
| 48 | //______________________________________________________________________________
|
|---|
| 49 | void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
|
|---|
| 50 | {
|
|---|
| 51 | Int_t i;
|
|---|
| 52 |
|
|---|
| 53 | //calculate chisquare
|
|---|
| 54 | Double_t chisq = 0;
|
|---|
| 55 | Double_t delta;
|
|---|
| 56 | nPixelsFitted=0;
|
|---|
| 57 | for (i=1;i<numPixels; i++) {
|
|---|
| 58 | if (isPixelsFitted[i])
|
|---|
| 59 | {
|
|---|
| 60 | if (errorz[i] != 0.0)
|
|---|
| 61 | {
|
|---|
| 62 | delta = (z[i]-func(x[i],y[i],par))/errorz[i];
|
|---|
| 63 | chisq += delta*delta;
|
|---|
| 64 | nPixelsFitted++;
|
|---|
| 65 | }
|
|---|
| 66 | else
|
|---|
| 67 | cout << "This should never happen errorz[" << i << "] " << errorz[i] << endl;
|
|---|
| 68 | }
|
|---|
| 69 | }
|
|---|
| 70 | f = chisq;
|
|---|
| 71 | chisquare = chisq;
|
|---|
| 72 | }
|
|---|
| 73 |
|
|---|
| 74 | //______________________________________________________________________________
|
|---|
| 75 | Double_t func(float x,float y,Double_t *par)
|
|---|
| 76 | {
|
|---|
| 77 | 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]));
|
|---|
| 78 | return value;
|
|---|
| 79 | }
|
|---|
| 80 |
|
|---|
| 81 | Bool_t HandleInput()
|
|---|
| 82 | {
|
|---|
| 83 | TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
|
|---|
| 84 |
|
|---|
| 85 | while (1)
|
|---|
| 86 | {
|
|---|
| 87 | //
|
|---|
| 88 | // While reading the input process gui events asynchronously
|
|---|
| 89 | //
|
|---|
| 90 | timer.TurnOn();
|
|---|
| 91 | TString input = Getline("Type 'q' to exit, <return> to go on: ");
|
|---|
| 92 | timer.TurnOff();
|
|---|
| 93 |
|
|---|
| 94 | if (input=="q\n")
|
|---|
| 95 | return kFALSE;
|
|---|
| 96 |
|
|---|
| 97 | if (input=="\n")
|
|---|
| 98 | return kTRUE;
|
|---|
| 99 | };
|
|---|
| 100 |
|
|---|
| 101 | return kFALSE;
|
|---|
| 102 | }
|
|---|
| 103 |
|
|---|
| 104 | void pointspreadfunction(TString fname, TString clname = "NULL", Int_t userNumLines = 1000)
|
|---|
| 105 | {
|
|---|
| 106 |
|
|---|
| 107 | // First run over continuos light files to have a DC calibration
|
|---|
| 108 |
|
|---|
| 109 |
|
|---|
| 110 | Float_t currmean[numPixels];
|
|---|
| 111 | Float_t currsquaremean[numPixels];
|
|---|
| 112 | Float_t currrms[numPixels];
|
|---|
| 113 | Float_t meanofcurrmean = 0;
|
|---|
| 114 |
|
|---|
| 115 | for (Int_t swpix=0; swpix<numPixels; swpix++)
|
|---|
| 116 | {
|
|---|
| 117 | currmean[swpix] = 0.;
|
|---|
| 118 | currsquaremean[swpix] = 0.;
|
|---|
| 119 | currrms[swpix] = 0.;
|
|---|
| 120 | }
|
|---|
| 121 |
|
|---|
| 122 | Int_t numLines=0;
|
|---|
| 123 |
|
|---|
| 124 | //containers
|
|---|
| 125 | MGeomCamMagic geomcam;
|
|---|
| 126 | MCameraDC curr;
|
|---|
| 127 | MCameraDC currbis;
|
|---|
| 128 |
|
|---|
| 129 | const Float_t conv4mm2deg = geomcam.GetConvMm2Deg();
|
|---|
| 130 |
|
|---|
| 131 | if (clname != "NULL")
|
|---|
| 132 | {
|
|---|
| 133 | MParList plist0;
|
|---|
| 134 |
|
|---|
| 135 | MTaskList tlist0;
|
|---|
| 136 | plist0.AddToList(&tlist0);
|
|---|
| 137 |
|
|---|
| 138 | plist0.AddToList(&geomcam);
|
|---|
| 139 | plist0.AddToList(&curr);
|
|---|
| 140 |
|
|---|
| 141 | MReportFileRead read0(clname);
|
|---|
| 142 | read0.SetHasNoHeader();
|
|---|
| 143 | read0.AddToList("MReportCurrents");
|
|---|
| 144 |
|
|---|
| 145 | tlist0.AddToList(&read0);
|
|---|
| 146 |
|
|---|
| 147 | MEvtLoop evtloop0;
|
|---|
| 148 | evtloop0.SetParList(&plist0);
|
|---|
| 149 |
|
|---|
| 150 |
|
|---|
| 151 | if (!evtloop0.PreProcess())
|
|---|
| 152 | return;
|
|---|
| 153 |
|
|---|
| 154 | while (tlist0.Process())
|
|---|
| 155 | {
|
|---|
| 156 | for (Int_t swpix=0; swpix<numPixels; swpix++)
|
|---|
| 157 | {
|
|---|
| 158 | meanofcurrmean += curr[swpix];
|
|---|
| 159 | currmean[swpix] += curr[swpix];
|
|---|
| 160 | currsquaremean[swpix] += curr[swpix]*curr[swpix];
|
|---|
| 161 | }
|
|---|
| 162 | numLines++;
|
|---|
| 163 | }
|
|---|
| 164 |
|
|---|
| 165 | evtloop0.PostProcess();
|
|---|
| 166 |
|
|---|
| 167 | meanofcurrmean /= (numLines*numPixels);
|
|---|
| 168 | for (Int_t swpix=0; swpix<numPixels; swpix++)
|
|---|
| 169 | {
|
|---|
| 170 |
|
|---|
| 171 | currmean[swpix] /= numLines;
|
|---|
| 172 | currsquaremean[swpix] /= numLines;
|
|---|
| 173 | currrms[swpix] = sqrt(fabs(currsquaremean[swpix] - currmean[swpix]*currmean[swpix]));
|
|---|
| 174 |
|
|---|
| 175 | curr[swpix] = currmean[swpix];
|
|---|
| 176 | currbis[swpix] = currrms[swpix];
|
|---|
| 177 |
|
|---|
| 178 | currmean[swpix] /= meanofcurrmean;
|
|---|
| 179 | currrms[swpix] /= meanofcurrmean;
|
|---|
| 180 |
|
|---|
| 181 | // cout << "swpixel " << swpix << " dc fCalib " << currmean[swpix] << " +- " << currrms[swpix] << endl;
|
|---|
| 182 | }
|
|---|
| 183 |
|
|---|
| 184 |
|
|---|
| 185 | // curr.Print();
|
|---|
| 186 | /* MHCamera display0(geomcam);
|
|---|
| 187 | display0.SetPrettyPalette();
|
|---|
| 188 | display0.Draw();
|
|---|
| 189 |
|
|---|
| 190 | display0.SetCamContent(currbis);
|
|---|
| 191 | cout << "PSF>> DC mean values drawn" << endl;
|
|---|
| 192 | // Remove the comments if you want to go through the file
|
|---|
| 193 | // event-by-event:
|
|---|
| 194 | if (!HandleInput())
|
|---|
| 195 | break;
|
|---|
| 196 | */
|
|---|
| 197 | }
|
|---|
| 198 | else
|
|---|
| 199 | {
|
|---|
| 200 | for (Int_t swpix=0; swpix<numPixels; swpix++)
|
|---|
| 201 | {
|
|---|
| 202 | currmean[swpix] = 1.;
|
|---|
| 203 | currrms[swpix] = 0.2;
|
|---|
| 204 | }
|
|---|
| 205 |
|
|---|
| 206 | }
|
|---|
| 207 |
|
|---|
| 208 | // Now we can run over the dc data to extract the psf rms.
|
|---|
| 209 |
|
|---|
| 210 | const Int_t maxNumLines = 10000;
|
|---|
| 211 |
|
|---|
| 212 | Double_t ux[maxNumLines];
|
|---|
| 213 | Double_t uy[maxNumLines];
|
|---|
| 214 | Double_t sx[maxNumLines];
|
|---|
| 215 | Double_t sy[maxNumLines];
|
|---|
| 216 | Double_t chisqu[maxNumLines];
|
|---|
| 217 | Double_t time[maxNumLines];
|
|---|
| 218 |
|
|---|
| 219 | MParList plist;
|
|---|
| 220 |
|
|---|
| 221 | MGeomCamMagic geomcam;
|
|---|
| 222 | MCameraDC curr;
|
|---|
| 223 | MTaskList tlist;
|
|---|
| 224 |
|
|---|
| 225 | plist.AddToList(&geomcam);
|
|---|
| 226 | plist.AddToList(&curr);
|
|---|
| 227 | plist.AddToList(&tlist);
|
|---|
| 228 |
|
|---|
| 229 | MReportFileRead read(fname);
|
|---|
| 230 | read.SetHasNoHeader();
|
|---|
| 231 | read.AddToList("MReportCurrents");
|
|---|
| 232 |
|
|---|
| 233 | tlist.AddToList(&read);
|
|---|
| 234 |
|
|---|
| 235 | MEvtLoop evtloop;
|
|---|
| 236 | evtloop.SetParList(&plist);
|
|---|
| 237 |
|
|---|
| 238 | if (!evtloop.PreProcess())
|
|---|
| 239 | return;
|
|---|
| 240 |
|
|---|
| 241 | MHCamera display(geomcam);
|
|---|
| 242 | display.SetPrettyPalette();
|
|---|
| 243 | // display.SetInvDeepBlueSeaPalette()
|
|---|
| 244 | display.Draw();
|
|---|
| 245 | gPad->SetLogy();
|
|---|
| 246 | gPad->cd(1);
|
|---|
| 247 |
|
|---|
| 248 | Double_t amin,edm,errdef;
|
|---|
| 249 | Int_t nvpar,nparx,icstat;
|
|---|
| 250 |
|
|---|
| 251 | Double_t max,maxerror;
|
|---|
| 252 | Double_t xmean,xsigma,ymean,ysigma;
|
|---|
| 253 | Double_t xmeanerror,xsigmaerror,ymeanerror,ysigmaerror;
|
|---|
| 254 |
|
|---|
| 255 | // TVector xellipsecenter;
|
|---|
| 256 | // TVector yellipsecenter;
|
|---|
| 257 |
|
|---|
| 258 | TEllipse ellipse;
|
|---|
| 259 | ellipse.SetFillStyle(4000);
|
|---|
| 260 | ellipse.SetLineWidth(2);
|
|---|
| 261 | ellipse.SetLineColor(2);
|
|---|
| 262 |
|
|---|
| 263 | ellipse.Draw();
|
|---|
| 264 |
|
|---|
| 265 | Int_t nbinsxy = 80;
|
|---|
| 266 | Float_t minxy = -600*conv4mm2deg;
|
|---|
| 267 | Float_t maxxy = 600*conv4mm2deg;
|
|---|
| 268 | Float_t fromdegtobin = (maxxy-minxy)/nbinsxy;
|
|---|
| 269 |
|
|---|
| 270 | TH2D psfhist("psfhist","",nbinsxy,minxy,maxxy,nbinsxy,minxy,maxxy);
|
|---|
| 271 | psfhist->GetXaxis()->SetTitle("[deg]");
|
|---|
| 272 | psfhist->GetYaxis()->SetTitle("[deg]");
|
|---|
| 273 | psfhist->GetZaxis()->SetTitle("DC [uA]");
|
|---|
| 274 |
|
|---|
| 275 | TCanvas *psfcanvas = new TCanvas("psfcanvas","Point Spread Funtion 2D",200,20,900,700);
|
|---|
| 276 |
|
|---|
| 277 | tlist.Process();
|
|---|
| 278 |
|
|---|
| 279 | Int_t numLines=0;
|
|---|
| 280 | Float_t minDCStar = 6.0;
|
|---|
| 281 |
|
|---|
| 282 | Int_t numPixelsInStar = 0;
|
|---|
| 283 | Float_t maxDC = 0;
|
|---|
| 284 | Int_t swpixelmaxDC;
|
|---|
| 285 |
|
|---|
| 286 | Bool_t isPixelsFittedTmp[numPixels];
|
|---|
| 287 |
|
|---|
| 288 | for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
|
|---|
| 289 | isPixelsFittedTmp[swpixel] = kFALSE;
|
|---|
| 290 |
|
|---|
| 291 | for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
|
|---|
| 292 | {
|
|---|
| 293 | if(curr[swpixel] > maxDC)
|
|---|
| 294 | {
|
|---|
| 295 | swpixelmaxDC = swpixel;
|
|---|
| 296 | maxDC = curr[swpixelmaxDC];
|
|---|
| 297 | }
|
|---|
| 298 |
|
|---|
| 299 | if(curr[swpixel]>minDCStar)
|
|---|
| 300 | {
|
|---|
| 301 | numPixelsInStar++;
|
|---|
| 302 | isPixelsFitted[swpixel] = kTRUE;
|
|---|
| 303 | }
|
|---|
| 304 | else
|
|---|
| 305 | isPixelsFitted[swpixel] = kFALSE;
|
|---|
| 306 | }
|
|---|
| 307 |
|
|---|
| 308 | if (numPixelsInStar == 0)
|
|---|
| 309 | {
|
|---|
| 310 | cout << "PSF>> Warning: none pixel over minDCStar(" << minDCStar << ')' << endl;
|
|---|
| 311 | return;
|
|---|
| 312 | }
|
|---|
| 313 |
|
|---|
| 314 | //1st neighboor ring
|
|---|
| 315 | for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
|
|---|
| 316 | if (isPixelsFitted[swpixel])
|
|---|
| 317 | for(Int_t next=0; next<geomcam[swpixel].GetNumNeighbors(); next++)
|
|---|
| 318 | isPixelsFittedTmp[geomcam[swpixel].GetNeighbor(next)] = kTRUE;
|
|---|
| 319 |
|
|---|
| 320 | for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
|
|---|
| 321 | if (isPixelsFittedTmp[swpixel])
|
|---|
| 322 | isPixelsFitted[swpixel] = kTRUE;
|
|---|
| 323 |
|
|---|
| 324 | //2on neighboor ring
|
|---|
| 325 | for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
|
|---|
| 326 | if (isPixelsFitted[swpixel])
|
|---|
| 327 | for(Int_t next=0; next<geomcam[swpixel].GetNumNeighbors(); next++)
|
|---|
| 328 | isPixelsFittedTmp[geomcam[swpixel].GetNeighbor(next)] = kTRUE;
|
|---|
| 329 |
|
|---|
| 330 | for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
|
|---|
| 331 | if (isPixelsFittedTmp[swpixel])
|
|---|
| 332 | isPixelsFitted[swpixel] = kTRUE;
|
|---|
| 333 |
|
|---|
| 334 |
|
|---|
| 335 | //3rt neighboor ring
|
|---|
| 336 | for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
|
|---|
| 337 | if (isPixelsFitted[swpixel])
|
|---|
| 338 | for(Int_t next=0; next<geomcam[swpixel].GetNumNeighbors(); next++)
|
|---|
| 339 | isPixelsFittedTmp[geomcam[swpixel].GetNeighbor(next)] = kTRUE;
|
|---|
| 340 |
|
|---|
| 341 | for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
|
|---|
| 342 | if (isPixelsFittedTmp[swpixel])
|
|---|
| 343 | isPixelsFitted[swpixel] = kTRUE;
|
|---|
| 344 |
|
|---|
| 345 |
|
|---|
| 346 | for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
|
|---|
| 347 | {
|
|---|
| 348 | // cout << "isPixelsFitted[" << swpixel << "] ";
|
|---|
| 349 | // cout << (Int_t)isPixelsFitted[swpixel] << endl;
|
|---|
| 350 | curr[swpixel] = (Float_t)isPixelsFitted[swpixel];
|
|---|
| 351 | }
|
|---|
| 352 |
|
|---|
| 353 | /* MHCamera display0(geomcam);
|
|---|
| 354 | display0.SetPrettyPalette();
|
|---|
| 355 | display0.Draw();
|
|---|
| 356 |
|
|---|
| 357 | display0.SetCamContent(curr);
|
|---|
| 358 | cout << "PSF>> Fitted pixels drawn" << endl;
|
|---|
| 359 | // Remove the comments if you want to go through the file
|
|---|
| 360 | // event-by-event:
|
|---|
| 361 | if (!HandleInput())
|
|---|
| 362 | break;
|
|---|
| 363 | */
|
|---|
| 364 | TMinuit *gMinuit = new TMinuit(7); //initialize TMinuit with a maximum of 5 params
|
|---|
| 365 | gMinuit->SetFCN(fcn);
|
|---|
| 366 |
|
|---|
| 367 | Double_t arglist[10];
|
|---|
| 368 | Int_t ierflg = 0;
|
|---|
| 369 |
|
|---|
| 370 | arglist[0] = 1;
|
|---|
| 371 | gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
|
|---|
| 372 | // arglist[0] = -1;
|
|---|
| 373 | arglist[0] = 0;
|
|---|
| 374 | gMinuit->mnexcm("SET PRI", arglist ,1,ierflg);
|
|---|
| 375 |
|
|---|
| 376 | // Set starting values and step sizes for parameters
|
|---|
| 377 | Double_t vstart[5];
|
|---|
| 378 | Double_t step[5];
|
|---|
| 379 | Double_t lowlimit[5] = {minDCStar, -2., 0.05, -2, 0.05};
|
|---|
| 380 | Double_t uplimit[5] = {30., 2., 1., 2., 1.};
|
|---|
| 381 |
|
|---|
| 382 | vstart[0] = maxDC;
|
|---|
| 383 | vstart[1] = geomcam[swpixelmaxDC].GetX()*conv4mm2deg;
|
|---|
| 384 | vstart[2] = 30*TMath::Sqrt(numPixelsInStar/2)*conv4mm2deg;
|
|---|
| 385 | vstart[3] = geomcam[swpixelmaxDC].GetY()*conv4mm2deg;
|
|---|
| 386 | vstart[4] = 30*TMath::Sqrt(numPixelsInStar/2)*conv4mm2deg;
|
|---|
| 387 |
|
|---|
| 388 | for(Int_t i=0; i<5; i++)
|
|---|
| 389 | {
|
|---|
| 390 | if (vstart[i] != 0)
|
|---|
| 391 | step[i] = TMath::Abs(vstart[i]/sqrt(2));
|
|---|
| 392 | else
|
|---|
| 393 | step[i] = 1;
|
|---|
| 394 |
|
|---|
| 395 | cout << "DGB vstat[" << i << "] " << vstart[i] << " step[" << i << "] " << step[i] <<endl;
|
|---|
| 396 | }
|
|---|
| 397 | // cout << endl;
|
|---|
| 398 |
|
|---|
| 399 | gMinuit->mnparm(0, "max", vstart[0], step[0], lowlimit[0], uplimit[0],ierflg);
|
|---|
| 400 | gMinuit->mnparm(1, "xmean", vstart[1], step[1], lowlimit[1], uplimit[1],ierflg);
|
|---|
| 401 | gMinuit->mnparm(2, "xsigma", vstart[2], step[2], lowlimit[2], uplimit[2],ierflg);
|
|---|
| 402 | gMinuit->mnparm(3, "ymean", vstart[3], step[3], lowlimit[3], uplimit[3],ierflg);
|
|---|
| 403 | gMinuit->mnparm(4, "ysigma", vstart[4], step[4], lowlimit[4], uplimit[4],ierflg);
|
|---|
| 404 |
|
|---|
| 405 | numLines++;
|
|---|
| 406 |
|
|---|
| 407 | while (tlist.Process() && numLines < maxNumLines)
|
|---|
| 408 | {
|
|---|
| 409 | // curr.Print();
|
|---|
| 410 |
|
|---|
| 411 | for (Int_t swpixel=1; swpixel<577; swpixel++)
|
|---|
| 412 | {
|
|---|
| 413 | // cout << swpixel << '\t' << curr[swpixel] << '\t' << geomcam[swpixel].GetX() << '\t' << geomcam[swpixel].GetY() << endl;
|
|---|
| 414 |
|
|---|
| 415 | x[swpixel] = geomcam[swpixel].GetX()*conv4mm2deg;
|
|---|
| 416 | y[swpixel] = geomcam[swpixel].GetY()*conv4mm2deg;
|
|---|
| 417 | z[swpixel] = curr[swpixel]/currmean[swpixel];
|
|---|
| 418 | 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]);
|
|---|
| 419 |
|
|---|
| 420 |
|
|---|
| 421 | psfhist->SetBinContent((Int_t)((x[swpixel]+600*conv4mm2deg)/fromdegtobin),(Int_t)((y[swpixel]+600*conv4mm2deg)/fromdegtobin),z[swpixel]);
|
|---|
| 422 | }
|
|---|
| 423 |
|
|---|
| 424 | psfcanvas->cd(1);
|
|---|
| 425 | psfhist->Draw("lego2");
|
|---|
| 426 |
|
|---|
| 427 | // Now ready for minimization step
|
|---|
| 428 | arglist[0] = 500;
|
|---|
| 429 | arglist[1] = 1.;
|
|---|
| 430 | gMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);
|
|---|
| 431 |
|
|---|
| 432 | // Print results
|
|---|
| 433 | /* gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
|
|---|
| 434 | gMinuit->mnprin(3,amin);
|
|---|
| 435 | */
|
|---|
| 436 | gMinuit->GetParameter(0,max,maxerror);
|
|---|
| 437 | gMinuit->GetParameter(1,xmean,xmeanerror);
|
|---|
| 438 | gMinuit->GetParameter(2,xsigma,xsigmaerror);
|
|---|
| 439 | gMinuit->GetParameter(3,ymean,ymeanerror);
|
|---|
| 440 | gMinuit->GetParameter(4,ysigma,ysigmaerror);
|
|---|
| 441 |
|
|---|
| 442 | /* cout << endl;
|
|---|
| 443 | cout << "numLine " << numLines << endl;
|
|---|
| 444 | cout << "max \t" << max << " +- " << maxerror << endl;
|
|---|
| 445 | cout << "xmean \t" << xmean << " +- " << xmeanerror << endl;
|
|---|
| 446 | cout << "xsigma \t" << TMath::Abs(xsigma) << " +- " << xsigmaerror << endl;
|
|---|
| 447 | cout << "ymean \t" << ymean << " +- " << ymeanerror << endl;
|
|---|
| 448 | cout << "ysigma \t" << TMath::Abs(ysigma) << " +- " << ysigmaerror << endl;
|
|---|
| 449 | cout << "chisquare/ndof \t" << chisquare/(nPixelsFitted-5) << endl;
|
|---|
| 450 | */
|
|---|
| 451 |
|
|---|
| 452 | chisqu[numLines]=chisquare/(nPixelsFitted-5);
|
|---|
| 453 |
|
|---|
| 454 | if(chisqu[numLines]<100.)
|
|---|
| 455 | {
|
|---|
| 456 | ux[numLines]=xmean;
|
|---|
| 457 | uy[numLines]=ymean;
|
|---|
| 458 | sx[numLines]=TMath::Abs(xsigma);
|
|---|
| 459 | sy[numLines]=TMath::Abs(ysigma);
|
|---|
| 460 | time[numLines]=numLines;
|
|---|
| 461 |
|
|---|
| 462 | display.SetCamContent(curr);
|
|---|
| 463 | gPad->cd(1);
|
|---|
| 464 | // ellipse.DrawEllipse(xmean,ymean,TMath::Abs(xsigma),TMath::Abs(ysigma),0,360,0);
|
|---|
| 465 | ellipse.SetX1(xmean/conv4mm2deg);
|
|---|
| 466 | ellipse.SetY1(ymean/conv4mm2deg);
|
|---|
| 467 | ellipse.SetR1(TMath::Abs(xsigma)/conv4mm2deg);
|
|---|
| 468 | ellipse.SetR2(TMath::Abs(ysigma)/conv4mm2deg);
|
|---|
| 469 |
|
|---|
| 470 | gPad->Modified();
|
|---|
| 471 | gPad->Update();
|
|---|
| 472 |
|
|---|
| 473 | /* if (numLines%10==0)
|
|---|
| 474 | {
|
|---|
| 475 | char imagepath[60];
|
|---|
| 476 | sprintf(&imagepath,"/home/Javi/mnt_users/Data/CaCo/img/ZetaTau_%04d.gif",numLines/10);
|
|---|
| 477 | cout << imagepath << endl;
|
|---|
| 478 | gPad->SaveAs(imagepath);
|
|---|
| 479 | }
|
|---|
| 480 | */
|
|---|
| 481 | numLines++;
|
|---|
| 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 | }
|
|---|
| 495 | }
|
|---|
| 496 |
|
|---|
| 497 |
|
|---|
| 498 | evtloop.PostProcess();
|
|---|
| 499 |
|
|---|
| 500 | //Draw the ditributions of the sigmas the point spread function
|
|---|
| 501 |
|
|---|
| 502 | cout<<"Number of lines "<<numLines<<endl;
|
|---|
| 503 |
|
|---|
| 504 | gROOT->Reset();
|
|---|
| 505 | gStyle->SetCanvasColor(0);
|
|---|
| 506 | gStyle->SetCanvasBorderMode(0);
|
|---|
| 507 | gStyle->SetPadBorderMode(0);
|
|---|
| 508 | gStyle->SetFrameBorderMode(0);
|
|---|
| 509 | gStyle->SetOptStat(00000000);
|
|---|
| 510 |
|
|---|
| 511 | c1 = new TCanvas("c1","Time evolution & distributions",0,0,1200,850);
|
|---|
| 512 | c1->Divide(3,2);
|
|---|
| 513 |
|
|---|
| 514 | c1->cd(1);
|
|---|
| 515 |
|
|---|
| 516 | TMath math;
|
|---|
| 517 |
|
|---|
| 518 | Double_t minmeanx, maxmeanx;
|
|---|
| 519 | minmeanx = ux[math.LocMin(numLines,ux)];
|
|---|
| 520 | maxmeanx = ux[math.LocMax(numLines,ux)];
|
|---|
| 521 |
|
|---|
| 522 | Double_t minmeany, maxmeany;
|
|---|
| 523 | minmeany = uy[math.LocMin(numLines,uy)];
|
|---|
| 524 | maxmeany = uy[math.LocMax(numLines,uy)];
|
|---|
| 525 |
|
|---|
| 526 | Double_t minmean, maxmean;
|
|---|
| 527 | minmean = math.Min(minmeanx,minmeany);
|
|---|
| 528 | maxmean = math.Max(maxmeanx,maxmeany);
|
|---|
| 529 |
|
|---|
| 530 | Double_t diff;
|
|---|
| 531 | diff = maxmean - minmean;
|
|---|
| 532 | diff = 0.1*diff;
|
|---|
| 533 | minmean = minmean - diff;
|
|---|
| 534 | maxmean = maxmean + diff;
|
|---|
| 535 |
|
|---|
| 536 | Double_t mintime, maxtime;
|
|---|
| 537 | mintime = time[math.LocMin(numLines,time)];
|
|---|
| 538 | maxtime = time[math.LocMax(numLines,time)];
|
|---|
| 539 |
|
|---|
| 540 | TH2D *h1 = new TH2D("h1",fname,1,mintime-1,maxtime+1,1,minmean,maxmean);
|
|---|
| 541 | h1->GetXaxis()->SetTitle("Event number");
|
|---|
| 542 | h1->GetYaxis()->SetTitle("mean position (deg)");
|
|---|
| 543 | h1->Draw();
|
|---|
| 544 |
|
|---|
| 545 | TGraph *grtimeevolmeanx = new TGraph(numLines,time,ux);
|
|---|
| 546 | grtimeevolmeanx->SetMarkerColor(3);
|
|---|
| 547 | grtimeevolmeanx->SetMarkerStyle(20);
|
|---|
| 548 | grtimeevolmeanx->SetMarkerSize (0.4);
|
|---|
| 549 | grtimeevolmeanx->Draw("P");
|
|---|
| 550 |
|
|---|
| 551 | TGraph *grtimeevolmeany = new TGraph(numLines,time,uy);
|
|---|
| 552 | grtimeevolmeany->SetMarkerColor(6);
|
|---|
| 553 | grtimeevolmeany->SetMarkerStyle(24);
|
|---|
| 554 | grtimeevolmeany->SetMarkerSize (0.4);
|
|---|
| 555 | grtimeevolmeany->Draw("P");
|
|---|
| 556 |
|
|---|
| 557 | legmeanxy = new TLegend(0.8,0.85,0.95,0.95);
|
|---|
| 558 | legmeanxy.SetTextSize(0.03);
|
|---|
| 559 | legmeanxy.AddEntry(grtimeevolmeanx,"mean x","P");
|
|---|
| 560 | legmeanxy.AddEntry(grtimeevolmeany,"mean y","P");
|
|---|
| 561 | legmeanxy.Draw();
|
|---|
| 562 |
|
|---|
| 563 | c1->cd(2);
|
|---|
| 564 |
|
|---|
| 565 | TMath math;
|
|---|
| 566 |
|
|---|
| 567 | Double_t minsigmax, maxsigmax;
|
|---|
| 568 | minsigmax = sx[math.LocMin(numLines,sx)];
|
|---|
| 569 | maxsigmax = sx[math.LocMax(numLines,sx)];
|
|---|
| 570 |
|
|---|
| 571 | Double_t minsigmay, maxsigmay;
|
|---|
| 572 | minsigmay = sy[math.LocMin(numLines,sy)];
|
|---|
| 573 | maxsigmay = sy[math.LocMax(numLines,sy)];
|
|---|
| 574 |
|
|---|
| 575 | Double_t minsigma, maxsigma;
|
|---|
| 576 | minsigma = math.Min(minsigmax,minsigmay);
|
|---|
| 577 | maxsigma = math.Max(maxsigmax,maxsigmay);
|
|---|
| 578 |
|
|---|
| 579 | diff = maxsigma - minsigma;
|
|---|
| 580 | diff = 0.1*diff;
|
|---|
| 581 | minsigma = minsigma - diff;
|
|---|
| 582 | maxsigma = maxsigma + diff;
|
|---|
| 583 |
|
|---|
| 584 | TH2D *h2 = new TH2D("h2","",1,mintime-1,maxtime+1,1,minsigma,maxsigma);
|
|---|
| 585 | h2->GetXaxis()->SetTitle("Event number");
|
|---|
| 586 | h2->GetYaxis()->SetTitle("PSF Rms (deg)");
|
|---|
| 587 | h2->Draw();
|
|---|
| 588 |
|
|---|
| 589 | TGraph* grtimeevolsigmax= new TGraph(numLines,time,sx);
|
|---|
| 590 | grtimeevolsigmax->SetMarkerColor(3);
|
|---|
| 591 | grtimeevolsigmax->SetMarkerStyle(20);
|
|---|
| 592 | grtimeevolsigmax->SetMarkerSize (0.4);
|
|---|
| 593 | grtimeevolsigmax->Draw("P");
|
|---|
| 594 |
|
|---|
| 595 | TGraph* grtimeevolsigmay= new TGraph(numLines,time,sy);
|
|---|
| 596 | grtimeevolsigmay->SetMarkerColor(6);
|
|---|
| 597 | grtimeevolsigmay->SetMarkerStyle(24);
|
|---|
| 598 | grtimeevolsigmay->SetMarkerSize (0.4);
|
|---|
| 599 | grtimeevolsigmay->Draw("P");
|
|---|
| 600 |
|
|---|
| 601 | legsigmaxy = new TLegend(0.8,0.85,0.95,0.95);
|
|---|
| 602 | legsigmaxy.SetTextSize(0.03);
|
|---|
| 603 | legsigmaxy.AddEntry(grtimeevolsigmax,"sigma x","P");
|
|---|
| 604 | legsigmaxy.AddEntry(grtimeevolsigmay,"sigma y","P");
|
|---|
| 605 | legsigmaxy.Draw();
|
|---|
| 606 |
|
|---|
| 607 | c1->cd(3);
|
|---|
| 608 |
|
|---|
| 609 | Double_t minchisqu, maxchisqu;
|
|---|
| 610 |
|
|---|
| 611 | minchisqu = chisqu[math.LocMin(numLines,chisqu)];
|
|---|
| 612 | maxchisqu = chisqu[math.LocMax(numLines,chisqu)];
|
|---|
| 613 |
|
|---|
| 614 | diff = maxchisqu - minchisqu;
|
|---|
| 615 | diff = 0.1*diff;
|
|---|
| 616 | minchisqu = minchisqu - diff;
|
|---|
| 617 | maxchisqu = maxchisqu + diff;
|
|---|
| 618 |
|
|---|
| 619 | TH2D *h3 = new TH2D("h3","",1,mintime-1,maxtime+1,1,minchisqu,maxchisqu);
|
|---|
| 620 | h3->GetXaxis()->SetTitle("Event number");
|
|---|
| 621 | h3->Draw();
|
|---|
| 622 |
|
|---|
| 623 | TGraph * grtimeevolchisqu = new TGraph(numLines,time,chisqu);
|
|---|
| 624 | grtimeevolchisqu->SetMarkerColor(215);
|
|---|
| 625 | grtimeevolchisqu->SetMarkerStyle(20);
|
|---|
| 626 | grtimeevolchisqu->SetMarkerSize(0.4);
|
|---|
| 627 | grtimeevolchisqu->Draw("P");
|
|---|
| 628 |
|
|---|
| 629 | legchisqu = new TLegend(0.55,0.90,0.95,0.95);
|
|---|
| 630 | legchisqu.SetTextSize(0.03);
|
|---|
| 631 | legchisqu.AddEntry(grtimeevolchisqu,"chi square / ndof","P");
|
|---|
| 632 | legchisqu.Draw();
|
|---|
| 633 |
|
|---|
| 634 | //***************************************
|
|---|
| 635 |
|
|---|
| 636 | const Int_t nbins = 100;
|
|---|
| 637 |
|
|---|
| 638 | TH1D *xsigmahist = new TH1D("xsigmahist","",nbins,minsigma,maxsigma);
|
|---|
| 639 | TH1D *ysigmahist = new TH1D("ysigmahist","",nbins,minsigma,maxsigma);
|
|---|
| 640 | TH1D *xmeanhist = new TH1D("xmeanhist","",nbins,minmean,maxmean);
|
|---|
| 641 | TH1D *ymeanhist = new TH1D("ymeanhist","",nbins,minmean,maxmean);
|
|---|
| 642 | TH1D *chisquhist = new TH1D("chisquhist","",nbins,minchisqu,maxchisqu);
|
|---|
| 643 |
|
|---|
| 644 | for (Int_t i=0; i<numLines; i++)
|
|---|
| 645 | {
|
|---|
| 646 | xsigmahist->Fill(TMath::Abs(sx[i]));
|
|---|
| 647 | ysigmahist->Fill(TMath::Abs(sy[i]));
|
|---|
| 648 | xmeanhist->Fill(ux[i]);
|
|---|
| 649 | ymeanhist->Fill(uy[i]);
|
|---|
| 650 | chisquhist->Fill(chisqu[i]);
|
|---|
| 651 |
|
|---|
| 652 | }
|
|---|
| 653 |
|
|---|
| 654 | c1->cd(5);
|
|---|
| 655 |
|
|---|
| 656 | // TCanvas *sigmacanvas = new TCanvas("sigmacanvas","Point Spread Funtion RMS",100,10,800,600);
|
|---|
| 657 |
|
|---|
| 658 | TMath math;
|
|---|
| 659 | Double_t maxsigma;
|
|---|
| 660 | Int_t binmaxx, binmaxy;
|
|---|
| 661 | xsigmahist->SetLineColor(3);
|
|---|
| 662 | xsigmahist->SetLineWidth(2);
|
|---|
| 663 | xsigmahist->SetXTitle("RMS [deg]");
|
|---|
| 664 | binmaxx = xsigmahist->GetMaximumBin();
|
|---|
| 665 | binmaxx = xsigmahist->GetBinContent(binmaxx);
|
|---|
| 666 |
|
|---|
| 667 | ysigmahist->SetLineColor(6);
|
|---|
| 668 | ysigmahist->SetLineWidth(2);
|
|---|
| 669 | binmaxy = ysigmahist->GetMaximumBin();
|
|---|
| 670 | binmaxy = ysigmahist->GetBinContent(binmaxy);
|
|---|
| 671 |
|
|---|
| 672 | maxsigma = math.Max(binmaxx,binmaxy);
|
|---|
| 673 | maxsigma += 1;
|
|---|
| 674 |
|
|---|
| 675 | xsigmahist->SetMaximum(maxsigma);
|
|---|
| 676 | ysigmahist->SetMaximum(maxsigma);
|
|---|
| 677 | xsigmahist->DrawCopy();
|
|---|
| 678 | ysigmahist->DrawCopy("Same");
|
|---|
| 679 |
|
|---|
| 680 | TLegend *legendsigma = new TLegend(.3,.8,.95,.95);
|
|---|
| 681 | legendsigma->SetTextSize(0.03);
|
|---|
| 682 | char xsigmatitle[100];
|
|---|
| 683 | char ysigmatitle[100];
|
|---|
| 684 | sprintf(xsigmatitle,"PSF Rms on X axis -- %5.2f +/- %5.2f deg",xsigmahist->GetMean(),xsigmahist->GetRMS());
|
|---|
| 685 | sprintf(ysigmatitle,"PSF Rms on Y axis -- %5.2f +/- %5.2f deg",ysigmahist->GetMean(),ysigmahist->GetRMS());
|
|---|
| 686 | legendsigma->AddEntry(xsigmahist,xsigmatitle,"F");
|
|---|
| 687 | legendsigma->AddEntry(ysigmahist,ysigmatitle,"F");
|
|---|
| 688 | legendsigma->Draw();
|
|---|
| 689 |
|
|---|
| 690 | c1->cd(4);
|
|---|
| 691 |
|
|---|
| 692 | Double_t maxmean;
|
|---|
| 693 |
|
|---|
| 694 | xmeanhist->SetLineColor(3);
|
|---|
| 695 | xmeanhist->SetLineWidth(2);
|
|---|
| 696 | xmeanhist->SetXTitle("mean [deg]");
|
|---|
| 697 | binmaxx = xmeanhist->GetMaximumBin();
|
|---|
| 698 | binmaxx = xmeanhist->GetBinContent(binmaxx);
|
|---|
| 699 |
|
|---|
| 700 | ymeanhist->SetLineColor(6);
|
|---|
| 701 | ymeanhist->SetLineWidth(2);
|
|---|
| 702 | binmaxy = ymeanhist->GetMaximumBin();
|
|---|
| 703 | binmaxy = ymeanhist->GetBinContent(binmaxy);
|
|---|
| 704 |
|
|---|
| 705 | maxmean = math.Max(binmaxx,binmaxy);
|
|---|
| 706 | maxmean += 1;
|
|---|
| 707 |
|
|---|
| 708 | xmeanhist->SetMaximum(maxmean);
|
|---|
| 709 | ymeanhist->SetMaximum(maxmean);
|
|---|
| 710 | xmeanhist->DrawCopy();
|
|---|
| 711 | ymeanhist->DrawCopy("Same");
|
|---|
| 712 |
|
|---|
| 713 | TLegend *legendmean = new TLegend(.35,.8,.95,.95);
|
|---|
| 714 | legendmean->SetTextSize(0.03);
|
|---|
| 715 | char xmeantitle[100];
|
|---|
| 716 | char ymeantitle[100];
|
|---|
| 717 | sprintf(xmeantitle,"mean on X axis -- %5.2f +/- %5.2f deg",xmeanhist->GetMean(),xmeanhist->GetRMS());
|
|---|
| 718 | sprintf(ymeantitle,"mean on Y axis -- %5.2f +/- %5.2f deg",ymeanhist->GetMean(),ymeanhist->GetRMS());
|
|---|
| 719 | legendmean->AddEntry(xmeanhist,xmeantitle,"F");
|
|---|
| 720 | legendmean->AddEntry(ymeanhist,ymeantitle,"F");
|
|---|
| 721 | legendmean->Draw();
|
|---|
| 722 |
|
|---|
| 723 | //meancanvas->Modified();
|
|---|
| 724 | //meancanvas->Update();
|
|---|
| 725 |
|
|---|
| 726 | c1->cd(6);
|
|---|
| 727 |
|
|---|
| 728 | chisquhist->SetLineColor(215);
|
|---|
| 729 | chisquhist->SetLineWidth(2);
|
|---|
| 730 | chisquhist->SetXTitle("chi square / ndof");
|
|---|
| 731 | TAxis * axis = chisquhist->GetXaxis();
|
|---|
| 732 | axis->SetLabelSize(0.025);
|
|---|
| 733 | chisquhist->DrawCopy();
|
|---|
| 734 |
|
|---|
| 735 | TLegend *legendchisqu = new TLegend(.4,.85,.95,.95);
|
|---|
| 736 | legendchisqu->SetTextSize(0.03);
|
|---|
| 737 | char chisqutitle[100];
|
|---|
| 738 | sprintf(chisqutitle,"chi square / ndof -- %5.2f +/- %5.2f ",chisquhist->GetMean(),chisquhist->GetRMS());
|
|---|
| 739 | legendchisqu->AddEntry(chisquhist,chisqutitle,"F");
|
|---|
| 740 | legendchisqu->Draw();
|
|---|
| 741 |
|
|---|
| 742 |
|
|---|
| 743 | return;
|
|---|
| 744 |
|
|---|
| 745 | }
|
|---|
| 746 |
|
|---|
| 747 |
|
|---|
| 748 |
|
|---|
| 749 |
|
|---|
| 750 |
|
|---|
| 751 |
|
|---|
| 752 |
|
|---|
| 753 |
|
|---|
| 754 |
|
|---|
| 755 |
|
|---|
| 756 |
|
|---|
| 757 |
|
|---|
| 758 |
|
|---|
| 759 |
|
|---|
| 760 |
|
|---|
| 761 |
|
|---|
| 762 |
|
|---|
| 763 |
|
|---|
| 764 |
|
|---|
| 765 |
|
|---|
| 766 |
|
|---|
| 767 |
|
|---|
| 768 |
|
|---|
| 769 |
|
|---|
| 770 |
|
|---|
| 771 |
|
|---|
| 772 |
|
|---|
| 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 |
|
|---|
