Changeset 2920
- Timestamp:
- 01/26/04 19:32:58 (21 years ago)
- Location:
- trunk/MagicSoft/Mars
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/Changelog
r2919 r2920 11 11 in arrays: gsNCells, gsNTrigPixels, 12 12 gsNPixInCell, gsNLutInCell, gsNPixInLut, fNumPixCell. 13 - Added method MMcTriggerLvl2::GetCellCompactPixel(int cell, MGeomCam *fCam) 13 - Added method 14 MMcTriggerLvl2::GetCellCompactPixel(int cell, MGeomCam *fCam) 14 15 which computes compact pixels into a given L2T macrocell. 15 - Added method MMcTriggerLvl2::CalcBiggerCellPseudoSize() which 16 computes fCellPseudoSize, the maximum Pseudo Size into L2T macrocells 17 - Added method MMcTriggerLvl2::GetCellPseudoSize() const which 18 returns fCellPseudoSize 19 - Added method MMcTriggerLvl2::IsPixelInCell(Int_t pixel, Int_t cell), 16 - Added method 17 MMcTriggerLvl2::CalcBiggerCellPseudoSize() 18 which computes fCellPseudoSize, the maximum Pseudo Size into L2T 19 macrocells 20 - Added method 21 MMcTriggerLvl2::GetCellPseudoSize() const 22 which returns fCellPseudoSize 23 - Added method 24 MMcTriggerLvl2::IsPixelInCell(Int_t pixel, Int_t cell), 20 25 which controls whether a pixel belongs to a given L2T cell. 21 - Added method MMcTriggerLvl2::GetMaxCell() const which returns fMaxCell, the cell with 22 the maximum fCellPseudoSize. 23 24 25 26 2004/01/26: Markus Gaug 26 - Added method 27 MMcTriggerLvl2::GetMaxCell() const 28 which returns fMaxCell, the cell with the maximum 29 fCellPseudoSize. 30 31 32 33 2004/01/26: Markus Gaug, Michele Doro 27 34 28 35 * manalysis/MArrivalTime.[h,cc], manalysis/MArrivalTimeCalc.[h,cc]: … … 33 40 34 41 * mcalib/MHCalibrationBlindPixel.[h,cc]: 35 - force mu_0 in Blind Pixel Fit to be around 0 in fKPoisson4 42 - force mu_{0} in Blind Pixel Fit to be around 0 in fKPoisson4 43 - implement combined Polya fit and Michele's back-scattered electron 44 fit 36 45 37 46 -
trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc
r2913 r2920 260 260 { 261 261 262 gStyle->SetOptFit( 0);262 gStyle->SetOptFit(1); 263 263 gStyle->SetOptStat(111111); 264 264 … … 366 366 break; 367 367 case kEPolya: 368 break;369 368 fSinglePheFit = new TF1("SinglePheFit",&fPolya,min,max,8); 369 break; 370 370 case kEMichele: 371 371 break; … … 381 381 const Double_t mu_1_guess = mu_0_guess + 50.; 382 382 const Double_t si_1_guess = si_0_guess + si_0_guess; 383 // Michele 384 const Double_t lambda_1cat_guess = 0.5; 385 const Double_t lambda_1dyn_guess = 0.5; 386 const Double_t mu_1cat_guess = mu_0_guess + 50.; 387 const Double_t mu_1dyn_guess = mu_0_guess + 20.; 388 const Double_t si_1cat_guess = si_0_guess + si_0_guess; 389 const Double_t si_1dyn_guess = si_0_guess; 390 // Polya 391 const Double_t excessPoisson_guess = 0.5; 392 const Double_t delta1_guess = 8.; 393 const Double_t delta2_guess = 5.; 394 const Double_t electronicAmp_guess = 0.0025; 395 383 396 // 384 397 // Initialize boundaries and start parameters … … 389 402 case kEPoisson4: 390 403 fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm); 391 fSinglePheFit->SetParNames("#lambda","#mu_ 0","#mu_1","#sigma_0","#sigma_1","Area");404 fSinglePheFit->SetParNames("#lambda","#mu_{0}","#mu_{1}","#sigma_{0}","#sigma_{1}","Area"); 392 405 fSinglePheFit->SetParLimits(0,0.,1.); 393 406 fSinglePheFit->SetParLimits(1,-2.,2.); … … 400 413 case kEPoisson6: 401 414 fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm); 402 fSinglePheFit->SetParNames("#lambda","#mu_ 0","#mu_1","#sigma_0","#sigma_1","Area");415 fSinglePheFit->SetParNames("#lambda","#mu_{0}","#mu_{1}","#sigma_{0}","#sigma_{1}","Area"); 403 416 fSinglePheFit->SetParLimits(0,0.,1.); 404 417 fSinglePheFit->SetParLimits(1,min,(max-min)/1.5); … … 410 423 411 424 case kEPolya: 425 fSinglePheFit->SetParameters(lambda_guess, excessPoisson_guess, 426 delta1_guess,delta2_guess, 427 electronicAmp_guess, 428 si_0_guess, 429 norm, mu_0_guess); 430 fSinglePheFit->SetParNames("#lambda","b_{tot}", 431 "#delta_{1}","#delta_{2}", 432 "amp_{e}","#sigma_{0}", 433 "Area", "#mu_{0}"); 434 fSinglePheFit->SetParLimits(0,0.,1.); 435 fSinglePheFit->SetParLimits(1,0.,1.); 436 fSinglePheFit->SetParLimits(2,6.,12.); 437 fSinglePheFit->SetParLimits(3,3.,8.); 438 fSinglePheFit->SetParLimits(4,0.,0.005); 439 fSinglePheFit->SetParLimits(5,min,(max-min)/1.5); 440 fSinglePheFit->SetParLimits(6,norm-0.1,norm+0.1); 441 fSinglePheFit->SetParLimits(7,-35.,15.); 412 442 break; 413 443 414 444 case kEMichele: 445 446 415 447 break; 416 448 … … 450 482 fSigma1Err = fSinglePheFit->GetParError(4); 451 483 break; 484 case kEPolya: 485 fLambda = fSinglePheFit->GetParameter(0); 486 fMu0 = fSinglePheFit->GetParameter(7); 487 fMu1 = 0.; 488 fSigma0 = fSinglePheFit->GetParameter(5); 489 fSigma1 = 0.; 490 491 fLambdaErr = fSinglePheFit->GetParError(0); 492 fMu0Err = fSinglePheFit->GetParError(7); 493 fMu1Err = 0.; 494 fSigma0Err = fSinglePheFit->GetParError(5); 495 fSigma1Err = 0.; 452 496 default: 453 497 break; … … 479 523 480 524 // Perform the cross-check fitting only the pedestal: 481 fSinglePhePedFit = new TF1("GausPed","gaus",rmin, fMu0);525 fSinglePhePedFit = new TF1("GausPed","gaus",rmin,0.); 482 526 fHBlindPixelCharge->Fit(fSinglePhePedFit,opt); 483 527 -
trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.h
r2913 r2920 71 71 72 72 //Getters 73 const Double_t GetLambda() const { return fLambda; }74 const Double_t GetLambdaCheck() const { return fLambdaCheck;}75 const Double_t GetMu0() const { return fMu0; }76 const Double_t GetMu1() const { return fMu1; }77 const Double_t GetSigma0() const { return fSigma0; }78 const Double_t GetSigma1() const { return fSigma1; }79 80 const Double_t GetLambdaErr() const { return fLambdaErr; }81 const Double_t GetLambdaCheckErr() 82 const Double_t GetMu0Err() const { return fMu0Err; }83 const Double_t GetMu1Err() const { return fMu1Err; }84 const Double_t GetSigma0Err() const { return fSigma0Err; }85 const Double_t GetSigma1Err() const { return fSigma1Err; }86 87 const Double_t GetChiSquare() const { return fChisquare; }88 const Double_t GetProb() const { return fProb;}89 const Int_t GetNdf() const { return fNdf;}90 91 const Double_t GetMeanTime() const { return fMeanTime;}92 const Double_t GetMeanTimeErr() const { return fMeanTimeErr; }93 const Double_t GetSigmaTime() const { return fSigmaTime; }94 const Double_t GetSigmaTimeErr() const { return fSigmaTimeErr;}95 96 const Bool_t IsFitOK() { return fFitOK;}97 98 const TH1I *GetHBlindPixelChargevsN() const { return fHBlindPixelChargevsN;}99 const TH1F *GetHBlindPixelPSD() const { return fHBlindPixelPSD;}73 const Double_t GetLambda() const { return fLambda; } 74 const Double_t GetLambdaCheck() const { return fLambdaCheck; } 75 const Double_t GetMu0() const { return fMu0; } 76 const Double_t GetMu1() const { return fMu1; } 77 const Double_t GetSigma0() const { return fSigma0; } 78 const Double_t GetSigma1() const { return fSigma1; } 79 80 const Double_t GetLambdaErr() const { return fLambdaErr; } 81 const Double_t GetLambdaCheckErr() const { return fLambdaCheckErr; } 82 const Double_t GetMu0Err() const { return fMu0Err; } 83 const Double_t GetMu1Err() const { return fMu1Err; } 84 const Double_t GetSigma0Err() const { return fSigma0Err; } 85 const Double_t GetSigma1Err() const { return fSigma1Err; } 86 87 const Double_t GetChiSquare() const { return fChisquare; } 88 const Double_t GetProb() const { return fProb; } 89 const Int_t GetNdf() const { return fNdf; } 90 91 const Double_t GetMeanTime() const { return fMeanTime; } 92 const Double_t GetMeanTimeErr() const { return fMeanTimeErr; } 93 const Double_t GetSigmaTime() const { return fSigmaTime; } 94 const Double_t GetSigmaTimeErr() const { return fSigmaTimeErr; } 95 96 const Bool_t IsFitOK() const { return fFitOK; } 97 98 const TH1I *GetHBlindPixelChargevsN() const { return fHBlindPixelChargevsN; } 99 const TH1F *GetHBlindPixelPSD() const { return fHBlindPixelPSD; } 100 100 101 101 // Draws … … 131 131 inline static Double_t fFitFuncMichele(Double_t *x, Double_t *par) 132 132 { 133 134 const Double_t qe = 0.247; 135 136 Double_t lambda = par[0]; 137 138 Double_t b1 = par[1]; 139 Double_t b2 = par[2]; 140 141 Double_t sum = 0.; 133 134 Double_t lambda1cat = par[0]; 135 Double_t lambda1dyn = par[1]; 136 Double_t mu0 = par[2]; 137 Double_t mu1cat = par[3]; 138 Double_t mu1dyn = par[4]; 139 Double_t sigma0 = par[5]; 140 Double_t sigma1cat = par[6]; 141 Double_t sigma1dyn = par[7]; 142 143 Double_t sumcat = 0.; 144 Double_t sumdyn = 0.; 142 145 Double_t arg = 0.; 143 146 144 return TMath::Exp(-1.*lambda)*par[5]*sum; 145 147 if (mu1cat < mu0) 148 return fNoWay; 149 150 if (sigma1cat < sigma0) 151 return fNoWay; 152 153 // if (sigma1cat < sigma1dyn) 154 // return NoWay; 155 156 //if (mu1cat < mu1dyn) 157 // return NoWay; 158 159 // if (lambda1cat < lambda1dyn) 160 // return NoWay; 161 162 Double_t mu2cat = (2.*mu1cat)-mu0; 163 Double_t mu2dyn = (2.*mu1dyn)-mu0; 164 Double_t mu3cat = (3.*mu1cat)-(2.*mu0); 165 Double_t mu3dyn = (3.*mu1dyn)-(2.*mu0); 166 167 Double_t sigma2cat = TMath::Sqrt((2.*sigma1cat*sigma1cat) - (sigma0*sigma0)); 168 Double_t sigma2dyn = TMath::Sqrt((2.*sigma1dyn*sigma1dyn) - (sigma0*sigma0)); 169 Double_t sigma3cat = TMath::Sqrt((3.*sigma1cat*sigma1cat) - (2.*sigma0*sigma0)); 170 Double_t sigma3dyn = TMath::Sqrt((3.*sigma1dyn*sigma1dyn) - (2.*sigma0*sigma0)); 171 172 Double_t lambda2cat = lambda1cat*lambda1cat; 173 Double_t lambda2dyn = lambda1dyn*lambda1dyn; 174 Double_t lambda3cat = lambda2cat*lambda1cat; 175 Double_t lambda3dyn = lambda2dyn*lambda1dyn; 176 177 // k=0: 178 arg = (x[0] - mu0)/sigma0; 179 sumcat = TMath::Exp(-0.5*arg*arg)/sigma0; 180 sumdyn =sumcat; 181 182 // k=1cat: 183 arg = (x[0] - mu1cat)/sigma1cat; 184 sumcat += lambda1cat*TMath::Exp(-0.5*arg*arg)/sigma1cat; 185 // k=1dyn: 186 arg = (x[0] - mu1dyn)/sigma1dyn; 187 sumdyn += lambda1dyn*TMath::Exp(-0.5*arg*arg)/sigma1dyn; 188 189 // k=2cat: 190 arg = (x[0] - mu2cat)/sigma2cat; 191 sumcat += 0.5*lambda2cat*TMath::Exp(-0.5*arg*arg)/sigma2cat; 192 // k=2dyn: 193 arg = (x[0] - mu2dyn)/sigma2dyn; 194 sumdyn += 0.5*lambda2dyn*TMath::Exp(-0.5*arg*arg)/sigma2dyn; 195 196 197 // k=3cat: 198 arg = (x[0] - mu3cat)/sigma3cat; 199 sumcat += 0.1666666667*lambda3cat*TMath::Exp(-0.5*arg*arg)/sigma3cat; 200 // k=3dyn: 201 arg = (x[0] - mu3dyn)/sigma3dyn; 202 sumdyn += 0.1666666667*lambda3dyn*TMath::Exp(-0.5*arg*arg)/sigma3dyn; 203 204 sumcat = TMath::Exp(-1.*lambda1cat)*sumcat; 205 sumdyn = TMath::Exp(-1.*lambda1dyn)*sumdyn; 206 207 return par[8]*(sumcat+sumdyn)/2.; 208 146 209 } 147 210 … … 339 402 340 403 } 404 405 inline static Double_t fPolya(Double_t *x, Double_t *par) 406 { 407 408 const Double_t QEcat = 0.247; // mean quantum efficiency 409 const Double_t sqrt2 = 1.4142135623731; 410 const Double_t sqrt3 = 1.7320508075689; 411 const Double_t sqrt4 = 2.; 412 413 const Double_t lambda = par[0]; // mean number of photons 414 415 const Double_t excessPoisson = par[1]; // non-Poissonic noise contribution 416 const Double_t delta1 = par[2]; // amplification first dynode 417 const Double_t delta2 = par[3]; // amplification subsequent dynodes 418 419 const Double_t electronicAmpl = par[4]; // electronic amplification and conversion to FADC charges 420 421 const Double_t pmtAmpl = delta1*delta2*delta2*delta2*delta2*delta2; // total PMT gain 422 const Double_t A = 1. + excessPoisson - QEcat 423 + 1./delta1 424 + 1./delta1/delta2 425 + 1./delta1/delta2/delta2; // variance contributions from PMT and QE 426 427 const Double_t totAmpl = QEcat*pmtAmpl*electronicAmpl; // Total gain and conversion 428 429 const Double_t mu0 = par[7]; // pedestal 430 const Double_t mu1 = totAmpl; // single phe position 431 const Double_t mu2 = 2*totAmpl; // double phe position 432 const Double_t mu3 = 3*totAmpl; // triple phe position 433 const Double_t mu4 = 4*totAmpl; // quadruple phe position 434 435 const Double_t sigma0 = par[5]; 436 const Double_t sigma1 = electronicAmpl*pmtAmpl*TMath::Sqrt(QEcat*A); 437 const Double_t sigma2 = sqrt2*sigma1; 438 const Double_t sigma3 = sqrt3*sigma1; 439 const Double_t sigma4 = sqrt4*sigma1; 440 441 const Double_t lambda2 = lambda*lambda; 442 const Double_t lambda3 = lambda2*lambda; 443 const Double_t lambda4 = lambda3*lambda; 444 445 //-- calculate the area---- 446 Double_t arg = (x[0] - mu0)/sigma0; 447 Double_t sum = TMath::Exp(-0.5*arg*arg)/sigma0; 448 449 // k=1: 450 arg = (x[0] - mu1)/sigma1; 451 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1; 452 453 // k=2: 454 arg = (x[0] - mu2)/sigma2; 455 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2; 456 457 // k=3: 458 arg = (x[0] - mu3)/sigma3; 459 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3; 460 461 // k=4: 462 arg = (x[0] - mu4)/sigma4; 463 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4; 464 465 return TMath::Exp(-1.*lambda)*par[6]*sum; 466 } 467 468 341 469 342 470 ClassDef(MHCalibrationBlindPixel, 1) // Histograms from the Calibration Blind Pixel
Note:
See TracChangeset
for help on using the changeset viewer.