Changeset 14166
- Timestamp:
- 06/13/12 12:46:17 (12 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
fact/tools/marsmacros/singlepe.C
r14165 r14166 3 3 #include <TSystem.h> 4 4 #include <TF1.h> 5 #include <TProfile.h> 6 #include <TProfile2D.h> 5 7 #include <TMath.h> 6 8 #include <TFitResultPtr.h> … … 28 30 #include "MHCamera.h" 29 31 #include "MGeomCamFACT.h" 32 #include "MRawRunHeader.h" 30 33 31 34 using namespace std; … … 36 39 struct Single 37 40 { 38 float fSignal;39 float fTime;41 float fSignal; 42 float fTime; 40 43 }; 41 44 42 45 class MSingles : public MParContainer, public MCamEvent 43 46 { 47 Int_t fIntegrationWindow; 48 44 49 vector<vector<Single>> fData; 45 50 46 51 public: 47 MSingles(const char *name=NULL, const char *title=NULL) 52 MSingles(const char *name=NULL, const char *title=NULL) : fIntegrationWindow(0) 48 53 { 49 54 fName = name ? name : "MSingles"; … … 60 65 UInt_t GetNumPixels() const { return fData.size(); } 61 66 67 void SetIntegrationWindow(Int_t w) { fIntegrationWindow = w; } 68 Int_t GetIntegrationWindow() const { return fIntegrationWindow; } 69 62 70 Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const 63 71 { 72 return kTRUE; 64 73 } 65 74 void DrawPixelContent(Int_t num) const { } … … 68 77 }; 69 78 79 class MH2F : public TH2F 80 { 81 public: 82 Int_t Fill(Double_t x,Double_t y) 83 { 84 Int_t binx, biny, bin; 85 fEntries++; 86 binx = TMath::Nint(x+1); 87 biny = fYaxis.FindBin(y); 88 bin = biny*(fXaxis.GetNbins()+2) + binx; 89 AddBinContent(bin); 90 if (fSumw2.fN) ++fSumw2.fArray[bin]; 91 if (biny == 0 || biny > fYaxis.GetNbins()) { 92 if (!fgStatOverflows) return -1; 93 } 94 ++fTsumw; 95 ++fTsumw2; 96 fTsumwy += y; 97 fTsumwy2 += y*y; 98 return bin; 99 } 100 ClassDef(MH2F, 1); 101 }; 102 103 class MProfile2D : public TProfile2D 104 { 105 public: 106 Int_t Fill(Double_t x, Double_t y, Double_t z) 107 { 108 Int_t bin,binx,biny; 109 fEntries++; 110 binx =TMath::Nint(x+1); 111 biny =fYaxis.FindBin(y); 112 bin = GetBin(binx, biny); 113 fArray[bin] += z; 114 fSumw2.fArray[bin] += z*z; 115 fBinEntries.fArray[bin] += 1; 116 if (fBinSumw2.fN) fBinSumw2.fArray[bin] += 1; 117 if (biny == 0 || biny > fYaxis.GetNbins()) { 118 if (!fgStatOverflows) return -1; 119 } 120 ++fTsumw; 121 ++fTsumw2; 122 fTsumwy += y; 123 fTsumwy2 += y*y; 124 fTsumwz += z; 125 fTsumwz2 += z*z; 126 return bin; 127 } 128 ClassDef(MProfile2D, 1); 129 }; 130 131 70 132 71 133 class MHSingles : public MH 72 134 { 73 MSingles *fSingles; 74 75 TH2F fHist; 135 MH2F fSignal; 136 MH2F fTime; 137 MProfile2D fPulse; 138 139 UInt_t fNumEntries; 140 141 MSingles *fSingles; //! 142 MPedestalSubtractedEvt *fEvent; //! 76 143 77 144 public: 78 MHSingles() : f Singles(0)145 MHSingles() : fNumEntries(0), fSingles(0), fEvent(0) 79 146 { 80 147 fName = "MHSingles"; 81 148 82 fHist.SetName("Name"); 83 fHist.SetTitle("Title"); 84 fHist.SetXTitle("X title"); 85 fHist.SetYTitle("Y title"); 86 fHist.SetDirectory(NULL); 149 fSignal.SetName("Signal"); 150 fSignal.SetTitle("Title"); 151 fSignal.SetXTitle("X title"); 152 fSignal.SetYTitle("Y title"); 153 fSignal.SetDirectory(NULL); 154 155 fTime.SetName("Time"); 156 fTime.SetTitle("Title"); 157 fTime.SetXTitle("X title"); 158 fTime.SetYTitle("Y title"); 159 fTime.SetDirectory(NULL); 160 161 fPulse.SetName("Pulse"); 162 fPulse.SetTitle("Title"); 163 fPulse.SetXTitle("X title"); 164 fPulse.SetYTitle("Y title"); 165 fPulse.SetDirectory(NULL); 87 166 } 88 167 … … 96 175 } 97 176 177 fEvent = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt"); 178 if (!fEvent) 179 { 180 *fLog << /* err << */ "MPedestalSubtractedEvt not found... abort." << endl; 181 return kFALSE; 182 } 183 184 fNumEntries = 0; 185 98 186 return kTRUE; 99 187 } 100 Bool_t ReInit(MParList *) 101 { 188 Bool_t ReInit(MParList *plist) 189 { 190 const MRawRunHeader *header = (MRawRunHeader*)plist->FindObject("MRawRunHeader"); 191 if (!header) 192 { 193 *fLog << /* err << */ "MRawRunHeader not found... abort." << endl; 194 return kFALSE; 195 } 196 197 const Int_t w = fSingles->GetIntegrationWindow(); 198 102 199 MBinning binsx, binsy; 103 200 binsx.SetEdges(fSingles->GetNumPixels(), -0.5, fSingles->GetNumPixels()-0.5); 104 binsy.SetEdges(300, -25, 600-25); 105 106 MH::SetBinning(fHist, binsx, binsy); 201 binsy.SetEdges(300, -7.5*w, (60-7.5)*w); 202 203 MH::SetBinning(fSignal, binsx, binsy); 204 205 const UShort_t roi = header->GetNumSamples(); 206 207 // Correct for DRS timing!! 208 MBinning binst(roi, -0.5, roi-0.5); 209 MH::SetBinning(fTime, binsx, binst); 210 211 MBinning binspy(2*roi, -roi-0.5, roi-0.5); 212 MH::SetBinning(fPulse, binsx, binspy); 107 213 108 214 return kTRUE; … … 111 217 Int_t Fill(const MParContainer *par, const Stat_t weight=1) 112 218 { 219 const Float_t *ptr = fEvent->GetSamples(0); 220 const Short_t roi = fEvent->GetNumSamples(); 221 113 222 for (unsigned int i=0; i<fSingles->GetNumPixels(); i++) 114 223 { … … 116 225 117 226 for (unsigned int j=0; j<result.size(); j++) 118 fHist.Fill(i, result[j].fSignal); 119 } 227 { 228 fSignal.Fill(i, result[j].fSignal); 229 fTime.Fill (i, result[j].fTime); 230 231 if (!ptr) 232 continue; 233 234 const Short_t tm = floor(result[j].fTime); 235 236 for (int s=0; s<roi; s++) 237 fPulse.Fill(i, s-tm, ptr[s]); 238 } 239 240 ptr += roi; 241 } 242 243 fNumEntries++; 120 244 121 245 return kTRUE; 122 246 } 123 247 124 TH2 *GetHist() { return &fHist; } 248 TH2 *GetSignal() { return &fSignal; } 249 TH2 *GetTime() { return &fTime; } 250 TH2 *GetPulse() { return &fPulse; } 251 252 UInt_t GetNumEntries() const { return fNumEntries; } 125 253 126 254 void Draw(Option_t *opt) … … 130 258 AppendPad(""); 131 259 132 pad->Divide( 1,1);260 pad->Divide(2,2); 133 261 134 262 pad->cd(1); 135 136 263 gPad->SetBorderMode(0); 137 264 gPad->SetFrameBorderMode(0); 138 139 fHist.Draw("colz"); 265 fSignal.Draw("colz"); 266 267 pad->cd(2); 268 gPad->SetBorderMode(0); 269 gPad->SetFrameBorderMode(0); 270 fTime.Draw("colz"); 271 272 pad->cd(3); 273 gPad->SetBorderMode(0); 274 gPad->SetFrameBorderMode(0); 275 fPulse.Draw("colz"); 140 276 } 141 277 … … 143 279 }; 144 280 281 MStatusDisplay *d = new MStatusDisplay; 282 145 283 MSingles *fSingles; 284 285 TH1F *fluct1 = new TH1F("", "", 100, -50, 50); 286 TH1F *fluct2 = new TH1F("", "", 100, -50, 50); 287 288 TGraph weights("weights.txt"); 146 289 147 290 Int_t PreProcess(MParList *plist) … … 158 301 } 159 302 303 d->AddTab("Fluct"); 304 fluct2->Draw(); 305 fluct1->Draw("same"); 306 307 fSingles->SetIntegrationWindow(20); 308 309 //if (weights.GetN()==0) 310 // return kFALSE; 311 160 312 return kTRUE; 161 313 } … … 163 315 Int_t Process() 164 316 { 165 166 for (int i=0; i<fEvent->GetNumPixels(); i++) 167 { 168 vector<Single> &result = fSingles->GetVector(i); 317 const UInt_t roi = fEvent->GetNumSamples(); 318 319 const size_t navg = 10; 320 const float threshold = 5; 321 const UInt_t start_skip = 20; 322 const UInt_t end_skip = 10; 323 const UInt_t integration_size = 10*3; 324 325 vector<float> val(roi-navg);//result of first sliding average 326 for (int pix=0; pix<fEvent->GetNumPixels(); pix++) 327 { 328 vector<Single> &result = fSingles->GetVector(pix); 169 329 result.clear(); 170 330 171 const UInt_t roi = fEvent->GetNumSamples(); 172 173 const Float_t *ptr = fEvent->GetSamples(i); 174 175 //--------------------------- 176 // ...Extract pulses here.... 177 //--------------------------- 178 179 //sliding average parameter (center weight = 3) 180 int avg1 = 8; 181 vector<float> Vslide(roi-2*avg1);//result of first sliding average 331 const Float_t *ptr = fEvent->GetSamples(pix); 182 332 183 333 //Sliding window 184 for (UInt_t k=avg1; k<roi-avg1; k++) 185 { 186 Vslide[k-avg1]=ptr[k]; 187 for(unsigned int offset=0; offset<avg1; offset++) 188 { 189 Vslide[k-avg1] += ptr[k+offset]; 190 Vslide[k-avg1] += ptr[k-offset]; 191 } 192 Vslide[k-avg1] /= float(2*avg1+1); 193 } 194 195 //peak finding via threshold 196 float threshold = 5; 197 Single single; 198 UInt_t start_sample = 5; 199 UInt_t end_sample = roi-2*avg1-21; 200 UInt_t min_dist = 20; 201 UInt_t integration_size = 5; 202 UInt_t integration_delay = 8; 203 204 for (UInt_t k=start_sample; k<end_sample; k++) 205 { 206 //search for threshold crossings 207 if( (Vslide[k-start_sample]<threshold) 208 &&(Vslide[k-1]<threshold) 209 &&(Vslide[k]>threshold) 210 &&(Vslide[k+5]>threshold) 211 ) 212 { 213 //average baseline before crossing 214 //double baseline = (Vslide[k-start_sample]+Vslide[k-start_sample+1]+Vslide[k-start_sample-1])/3.; 215 single.fSignal = 0; 216 single.fTime = k+integration_delay; 217 for(UInt_t l=integration_delay; l<integration_delay+integration_size; l++) 218 { 219 single.fSignal+=Vslide[k+l]; 220 } 221 //single.fSignal-=baseline;//baseline subtraction 222 result.push_back(single); 223 k+=min_dist; 224 } 334 for (size_t i=0; i<roi-navg; i++) 335 { 336 val[i] = 0; 337 for (size_t j=i; j<i+navg; j++) 338 val[i] += ptr[j]; 339 val[i] /= navg; 340 } 341 342 //peak finding via threshold 343 for (UInt_t i=start_skip; i<roi-navg-end_skip-10; i++) 344 { 345 //search for threshold crossings 346 if (val[i+0]>threshold || 347 val[i+4]>threshold || 348 349 val[i+5]<threshold || 350 val[i+9]<threshold) 351 continue; 352 353 // Evaluate arrival time more precisely!!! 354 // Get a better integration window 355 356 Single single; 357 single.fSignal = 0; 358 single.fTime = i; 359 360 // Crossing of the threshold is at 5 361 for (UInt_t j=i+5; j<i+5+integration_size; j++) 362 single.fSignal += ptr[j]; 363 364 result.push_back(single); 365 366 i += 5+integration_size; 225 367 } 226 368 } … … 233 375 } 234 376 235 Double_t 236 spek_fit( 237 Double_t* x, 238 Double_t* par 239 ); 240 241 int singlepe_final() 377 Double_t fcn(Double_t *x, Double_t *par); 378 379 int singlepe_final2() 242 380 { 243 381 // ====================================================== 244 382 245 const char *drsfile = "/fact/raw/2012/0 2/25/20120225_004.drs.fits.gz";383 const char *drsfile = "/fact/raw/2012/06/01/20120601_003.drs.fits.gz"; 246 384 247 385 MDirIter iter; … … 258 396 // map file to use (get that from La Palma!) 259 397 const char *map = usemap ? "/scratch_nfs/FACTmap111030.txt" : NULL; 260 gStyle->SetPalette(1,0); 398 261 399 // ------------------------------------------------------ 262 400 … … 295 433 // ====================================================== 296 434 297 MStatusDisplay *d = new MStatusDisplay;435 //MStatusDisplay *d = new MStatusDisplay; 298 436 299 437 MBadPixelsCam badpixels; … … 316 454 //MDrsCalibrationTime timecam; 317 455 318 gStyle->SetOptFit(kTRUE);456 MH::SetPalette(); 319 457 320 458 // ====================================================== … … 326 464 327 465 MSingles singles; 466 MRawRunHeader header; 328 467 329 468 MParList plist1; … … 332 471 plist1.AddToList(&badpixels); 333 472 plist1.AddToList(&singles); 473 plist1.AddToList(&header); 334 474 335 475 MEvtLoop loop1("DetermineCalConst"); … … 371 511 return 11; 372 512 513 gStyle->SetOptFit(kTRUE); 514 373 515 MHSingles* h = (MHSingles*) plist1.FindObject("MHSingles"); 374 516 if (!h) 375 517 return 0; 376 518 377 TH2 *hist = h->GetHist(); 519 TH2 *hsignal = h->GetSignal(); 520 TH2 *htime = h->GetTime(); 521 TH2 *hpulse = h->GetPulse(); 522 523 d->AddTab("Time"); 524 TH1 *ht = htime->ProjectionY(); 525 ht->Scale(1./1440); 526 ht->Draw(); 527 528 d->AddTab("Pulse"); 529 TH1 *hp = hpulse->ProjectionY(); 530 hp->Scale(1./1440); 531 hp->Draw(); 532 378 533 379 534 // ------------------ Spectrum Fit --------------- … … 381 536 // access the spektrumhistogram by the pointer *hist 382 537 383 hist->SetTitle("Amplitude Spectrum vs. Pixel ID");384 hist->SetXTitle("Pixel SoftID");385 hist->SetYTitle("Amplitude [mV]");386 hist->SetZTitle("counts");387 388 538 MGeomCamFACT fact; 389 539 MHCamera hcam(fact); 390 540 391 struct fit392 {393 float fAmplitude;394 float fGain;395 float f2GainRMS;396 float fCrosstlk;397 float fOffset;398 };399 fit fitParameters[singles.GetNumPixels()];400 401 TCanvas &canv1 = d->AddTab("Fits");402 TCanvas &canv2 = d->AddTab("Distributions");403 canv2.Divide(3,2);404 TCanvas &gain_canv = d->AddTab("Gain Distribution");405 406 541 //built an array of Amplitude spectra 407 TH1D* hAmplSpek[singles.GetNumPixels()]; 408 TString histotitle; 409 410 TF1 gaus4("g4","gaus"); 411 gaus4.SetLineColor(2); 412 413 TH1F hAmplitude( 414 "hAmplitude", 415 "Amplitude of Single Preak", 416 200, 417 50, 418 250 419 ); 420 hAmplitude.GetXaxis()->SetTitle("Amplitude (from Integral 5 Slices) [mV]"); 421 hAmplitude.GetYaxis()->SetTitle("counts"); 422 423 TH1F hGain( 424 "hGain", 425 "Gain", 426 40*5, 427 30, 428 70 429 ); 430 hGain.GetXaxis()->SetTitle("Gain (from Integral 5 Slices) [mV]"); 431 hGain.GetYaxis()->SetTitle("counts"); 432 433 TH1F hGainRMS( 434 "hSigma", 435 "Sigma of gauss peak in gainfit", 436 20*5, 437 0, 438 40 439 ); 440 hGainRMS.GetXaxis()->SetTitle("Sigma (from Integral 5 Slices) [mV]"); 441 hGainRMS.GetYaxis()->SetTitle("counts"); 442 443 TH1F hCrosstlk( 444 "hCrosstlk", 445 "Parameter proportional to Crosstalk", 446 30*5, 447 0, 448 4 449 ); 450 hCrosstlk.GetXaxis()->SetTitle("~Crosstalk (from Integral 5 Slices) [mV]"); 451 hCrosstlk.GetYaxis()->SetTitle("counts"); 452 453 TH1F hOffset( 454 "hOffset", 455 "X-Offset of Spectrum", 456 100*2, 457 0, 458 4 459 ); 460 hOffset.GetXaxis()->SetTitle(" X-Offset of Spectrum (from Integral 5 Slices) [mV]"); 461 hOffset.GetYaxis()->SetTitle("counts"); 462 463 TH1F hFitProb( 464 "hFitProb", 465 "Fit Probability", 466 100*2, 467 0, 468 100 469 ); 470 hFitProb.GetXaxis()->SetTitle("fit-probability [%]"); 471 hFitProb.GetYaxis()->SetTitle("counts"); 542 TH1F hAmplitude("Rate", "Average number of dark counts per event", 543 100, 0, 10); 544 TH1F hGain ("Gain", "Gain distribution", 545 50, 100, 250); 546 TH1F hGainRMS ("RelSigma", "Distribution of rel. Sigma", 547 100, 0, 30); 548 TH1F hCrosstlk ("Crosstalk", "Distribution of crosstalk probability", 549 100, 0, 25); 550 TH1F hOffset ("Offset", "Baseline offset distribution", 551 50, -7.5, 2.5); 552 TH1F hFitProb ("FitProb", "FitProb distribution", 553 50, 0, 100); 554 TH1F hSum1 ("Sum1", "Sum", 555 100, 0, 2000); 556 TH1F hSum2 ("Sum2", "Sum", 557 100, 0, 10); 472 558 473 559 // define fit function for Amplitudespectrum 474 TF1 fspek_fit("fspek_fit", spek_fit, 10, 600, 5); 475 fspek_fit.SetParNames("Ampl 1.Max", "Gain", "Sigma", "Ampl 2.Max", "Offset"); 476 477 // fit parameters 478 Double_t par[5]; 560 TF1 func("spektrum", fcn, 0, 1000, 5); 561 func.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset"); 562 func.SetLineColor(kBlue); 479 563 480 564 // Map for which pixel shall be plotted and which not 481 565 TArrayC usePixel(1440); 566 memset(usePixel.GetArray(), 1, 1440); 482 567 483 568 // List of Pixel that should be ignored in camera view 484 Int_t doNotUse[12]= {485 424,486 923,487 1208,488 583,489 830,490 1399,491 113,492 115,493 354,494 423,495 1195,496 1393 497 };569 usePixel[424] = 0; 570 usePixel[923] = 0; 571 usePixel[1208] = 0; 572 usePixel[583] = 0; 573 usePixel[830] = 0; 574 usePixel[1399] = 0; 575 usePixel[113] = 0; 576 usePixel[115] = 0; 577 usePixel[354] = 0; 578 usePixel[423] = 0; 579 usePixel[1195] = 0; 580 usePixel[1393] = 0; 581 582 int count_ok = 0; 498 583 499 584 // Begin of Loop over Pixels 500 585 for (UInt_t pixel = 0; pixel < singles.GetNumPixels(); pixel++) 501 586 { 502 usePixel[pixel] = 0; 503 504 //define projection histogram title 505 histotitle = "hPixelPro_"; 506 histotitle += pixel; 587 if (usePixel[pixel]==0) 588 continue; 507 589 508 590 //Projectipon of a certain Pixel out of Ampl.Spectrum 509 hAmplSpek[pixel] = hist->ProjectionY( histotitle , pixel+1, pixel+1); 510 511 //set histogram and Axis title 512 histotitle = "Amplitude Spectrum of Pixel #"; 513 histotitle += pixel; 514 hAmplSpek[pixel]->SetTitle(histotitle); 515 516 // define appearance of fit 517 fspek_fit.SetLineColor(kBlue); 518 fspek_fit.SetLineStyle(2); 519 520 par[0] = hAmplSpek[pixel]->GetBinContent(hAmplSpek[pixel]->GetMaximumBin()); 591 TH1D *hist = hsignal->ProjectionY("proj", pixel+1, pixel+1); 592 hist->SetDirectory(0); 593 594 const Int_t maxbin = hist->GetMaximumBin(); 595 const Double_t binwidth = hist->GetBinWidth(maxbin); 596 const Double_t ampl = hist->GetBinContent(maxbin); 597 598 double fwhm = 0; 599 for (int i=1; i<maxbin; i++) 600 if (hist->GetBinContent(maxbin-i)+hist->GetBinContent(maxbin+i)<ampl*2) 601 { 602 fwhm = (2*i+1)*binwidth; 603 break; 604 } 605 606 if (fwhm==0) 607 { 608 cout << "Could not determine start value for sigma." << endl; 609 continue; 610 } 611 612 // Amplitude, Gain, Sigma, Crosstalk Prob, Shift 521 613 // par[1] + par[4] is the position of the first maximum 522 par[1] = 0.8*hAmplSpek[pixel]->GetBinCenter(hAmplSpek[pixel]->GetMaximumBin()); 523 par[2] = 0.2*par[1]; 524 par[3] = 0.1*par[0]; 525 par[4] = 0.2*hAmplSpek[pixel]->GetBinCenter(hAmplSpek[pixel]->GetMaximumBin()); 526 fspek_fit.SetParameters(par); 614 Double_t par[5] = 615 { 616 ampl, hist->GetBinCenter(maxbin), fwhm*1.4, 0.1, -10 617 }; 618 619 func.SetParameters(par); 620 621 const double min = par[1]-fwhm*1.4; 622 const double max = par[1]*5; 527 623 528 624 //Fit Pixels spectrum 529 canv1.cd(); 530 TFitResultPtr fitStatus(hAmplSpek[ pixel ]->Fit(&fspek_fit,"+S","",0.8*par[1],5*par[1])); 531 double fit_prob = fitStatus->Prob(); 532 533 //cout << "pixel: " << pixel << "\t status: " << fitStatus << endl; 534 535 fitParameters[pixel].fAmplitude = fspek_fit.GetParameter(0); 536 fitParameters[pixel].fGain = fspek_fit.GetParameter(1); 537 fitParameters[pixel].f2GainRMS = fspek_fit.GetParameter(2); 538 fitParameters[pixel].fCrosstlk = fspek_fit.GetParameter(3)/fspek_fit.GetParameter(0); 539 fitParameters[pixel].fOffset = fspek_fit.GetParameter(4); 540 541 gain_canv.cd(); 542 gPad->SetLogy(1); 543 hGain.Fill(fitParameters[pixel].fGain); 544 hGain.Fit(&gaus4); 545 hGain.DrawCopy(); 546 547 canv2.cd(2); 548 gPad->SetLogy(1); 549 hGain.DrawCopy(); 550 551 canv2.cd(3); 552 gPad->SetLogy(1); 553 hGainRMS.Fill(fitParameters[pixel].f2GainRMS); 554 hGainRMS.DrawCopy(); 555 556 canv2.cd(1); 557 gPad->SetLogy(1); 558 hAmplitude.Fill(fitParameters[pixel].fAmplitude); 559 hAmplitude.DrawCopy(); 560 561 canv2.cd(4); 562 gPad->SetLogy(1); 563 hCrosstlk.Fill(fitParameters[pixel].fCrosstlk); 564 hCrosstlk.DrawCopy(); 565 566 canv2.cd(5); 567 gPad->SetLogy(1); 568 hOffset.Fill(fitParameters[pixel].fOffset); 569 hOffset.DrawCopy(); 570 571 canv2.cd(6); 572 gPad->SetLogy(1); 625 const TFitResultPtr rc = hist->Fit(&func, "N0QS", "", min, max); 626 627 const bool ok = int(rc)==0; 628 629 const double fit_prob = rc->Prob(); 630 const float fGain = func.GetParameter(1); 631 const float fAmplitude = func.Integral(0, 7*fGain)/h->GetNumEntries(); //func.GetParameter(0); 632 const float f2GainRMS = func.GetParameter(2); 633 const float fCrosstlk = func.GetParameter(3); 634 const float fOffset = func.GetParameter(4)/singles.GetIntegrationWindow(); 635 636 hAmplitude.Fill(fAmplitude); 637 hGain.Fill(fGain); 638 //hGainRMS.Fill(f2GainRMS); 639 hCrosstlk.Fill(fCrosstlk*100); 640 hOffset.Fill(fOffset); 573 641 hFitProb.Fill(100*fit_prob); 574 hFitProb.DrawCopy(); 575 576 577 hcam.SetBinContent(pixel+1, fitParameters[pixel].fGain); 578 if (!fitStatus) 579 { 580 usePixel[pixel] = 1; 581 cout << "Pixel: " << pixel << " success!" << endl << endl; 582 } 583 else 584 { 585 cout << "Pixel: " << pixel << " failed!" << endl << endl; 586 } 587 588 } 589 // TCanvas &canvBroken[12]; 590 // TString name; 591 for (int i = 0; i < 12; i++) 592 { 593 // name = "Pixel"; 594 // name += doNotUse[i]; 595 // canvBroken[i] = d->AddTab(name); 596 // hAmplSpek[doNotUse[ i ]]->Draw(); 597 598 usePixel[doNotUse[ i ]] = 0; 642 hGainRMS.Fill(100*f2GainRMS/fGain); 643 644 hcam.SetBinContent(pixel+1, fGain); 645 646 usePixel[pixel] = ok; 647 648 if (!ok) 649 { 650 cout << "Pixel " << setw(4) << pixel << ": failed!" << endl; 651 continue; 652 } 653 654 for (int b=1; b<=hist->GetNbinsX(); b++) 655 { 656 hSum1.Fill(hist->GetBinCenter(b), hist->GetBinContent(b)); 657 hSum2.Fill((hist->GetBinCenter(b)-fOffset)/fGain, hist->GetBinContent(b)); 658 } 659 660 if (count_ok==0) 661 { 662 TCanvas &c = d->AddTab(Form("Pix%d", pixel)); 663 c.cd(); 664 hist->SetName(Form("Pix%d", pixel)); 665 hist->DrawCopy()->SetDirectory(0); 666 //func.DrawCopy("same"); 667 } 668 669 count_ok++; 670 671 delete hist; 599 672 } 600 673 … … 606 679 hcam.DrawCopy(); 607 680 608 609 canv1.cd(); 681 TCanvas &c11 = d->AddTab("SumHist"); 682 gPad->SetLogy(); 683 hSum1.DrawCopy(); 684 685 TCanvas &c12 = d->AddTab("GainHist"); 686 gPad->SetLogy(); 687 hSum2.DrawCopy(); 688 689 TCanvas &c2 = d->AddTab("Result"); 690 c2.Divide(3,2); 691 c2.cd(1); 692 gPad->SetLogy(); 693 hAmplitude.DrawCopy(); 694 695 c2.cd(2); 696 gPad->SetLogy(); 697 hGain.DrawCopy(); 698 699 c2.cd(3); 700 gPad->SetLogy(); 701 hGainRMS.DrawCopy(); 702 703 c2.cd(4); 704 gPad->SetLogy(); 705 hCrosstlk.DrawCopy(); 706 707 c2.cd(5); 708 gPad->SetLogy(); 709 hOffset.DrawCopy(); 710 711 c2.cd(6); 712 gPad->SetLogy(); 713 hFitProb.DrawCopy(); 610 714 611 715 /* 612 TString title = "-- Calibrated signal #"; 613 title += seq.GetSequence(); 614 title += " ("; 615 title += drsfile; 616 title += ") --"; 617 d->SetTitle(title, kFALSE); 618 619 TString path; 620 path += Form("%s/20%6d_%03d-calibration.root", outpath, 621 seq.GetSequence()/1000, seq.GetSequence()%1000); 622 623 d->SaveAs(path); 624 */ 625 716 TCanvas &c3 = d->AddTab("Separation"); 717 gPad->SetLogy(); 718 hSep.DrawCopy(); 719 */ 626 720 return 0; 627 721 } 628 722 629 Double_t 630 spek_fit( 631 Double_t* x, 632 Double_t* par 633 ) 634 { 635 Double_t arg = 0; 636 Double_t cross = 1; 637 Double_t peak = 0; 638 // Double_t arg2 = 0; 639 // Double_t arg3 = 0; 640 // Double_t arg4 = 0; 641 642 // if (par[0] != 0) cross = par[3]/par[0]; 643 // if (par[2] != 0) arg = (x[0] - par[1] + par[4])/par[2]; 644 // if (par[2] != 0) arg2 = (x[0] - 2*par[1] + par[4])/par[2]; 645 // if (par[2] != 0) arg3 = (x[0] - 3*par[1] + par[4])/par[2]; 646 // if (par[2] != 0) arg4 = (x[0] - 4*par[1] + par[4])/par[2]; 647 // if (par[0] = 0) return 0; 648 // if (par[2] = 0) return 0; 649 650 if (par[0] != 0) cross = par[3]/par[0]; 723 Double_t fcn(Double_t *xx, Double_t *par) 724 { 725 Double_t x = xx[0]; 726 727 Double_t ampl = par[0]; 728 Double_t gain = par[1]; 729 Double_t sigma = par[2]; 730 Double_t cross = par[3]; 731 Double_t shift = par[4]; 732 733 Double_t xTalk = 1; 734 735 Double_t y = 0; 651 736 for (int order = 1; order < 5; order++) 652 737 { 653 arg = (x[0] - order*par[1] - par[4] )/par[2]; 654 Double_t xTalk = TMath::Power(cross, order-1); 655 // Double_t gauss = TMath::Exp(-0.5 * arg * arg/order); 656 Double_t gauss = TMath::Exp(-0.5 * arg * arg); 657 peak += xTalk*par[0]*gauss; 658 } 659 660 // peak += par[0]*TMath::Exp(-0.5 * arg * arg); 661 // peak += cross*par[0]*TMath::Exp(-0.5 * arg2 * arg2 /2); 662 // peak += cross*cross*par[0]*TMath::Exp(-0.5 * arg3 * arg3 /3); 663 // peak += cross*cross*cross*par[0]*TMath::Exp(-0.5 * arg4 * arg4 /4); 664 665 return peak; 738 Double_t arg = (x - order*gain - shift)/sigma; 739 740 Double_t gauss = TMath::Exp(-0.5 * arg * arg/order); 741 742 y += xTalk*gauss; 743 744 xTalk *= cross; 745 } 746 747 return ampl*y; 666 748 }
Note:
See TracChangeset
for help on using the changeset viewer.