Changeset 2012 for trunk/MagicSoft/Mars
- Timestamp:
- 04/25/03 10:58:26 (22 years ago)
- Location:
- trunk/MagicSoft/Mars
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/Changelog
r2010 r2012 15 15 * mhistmc/MHMcEnergyMigration.h: 16 16 - fixed the includes 17 18 * mgui/MCamDisplay.cc: 19 - changed autoscaling (max<1:max=1 --> max==min:max=min+1) 20 21 * manalysis/MBlindPixelCalc.cc: 22 - interpolate: take pixel area into account 23 24 * mhist/MHSigmaTheta.h: 25 - removed nonsense GetSigmaThetaByName(const TString name) 26 - removed nonsense GetSigmaPixThetaByName(const TString name) 27 - removed nonsense GetDiffPixThetaByName(const TString name) 28 29 * manalysis/MPadSchweizer.cc: 30 - fixed naming 31 - fixed usage of operators 32 - added some const qualifiers 33 - replaced 'int OK' by 'Bool_t ok' 34 - fixed wrong usage floating point value 0 17 35 18 36 -
trunk/MagicSoft/Mars/manalysis/MPadSchweizer.cc
r2001 r2012 280 280 //*fLog << "Entry MPadSchweizer::Process();" << endl; 281 281 282 Int_t rc ;282 Int_t rc=0; 283 283 284 284 const UInt_t npix = fEvt->GetNumPixels(); 285 285 286 Double_t SigbarOld; 287 288 289 SigbarOld = fSigmabar->Calc(*fCam, *fPed, *fEvt); 286 fSigmabar->Calc(*fCam, *fPed, *fEvt); 290 287 //*fLog << "before padding : " << endl; 291 288 //fSigmabar->Print(""); … … 308 305 // Calculate average sigma of the event 309 306 // 310 SigbarOld = fSigmabar->Calc(*fCam, *fPed, *fEvt);307 Double_t sigbarold = fSigmabar->Calc(*fCam, *fPed, *fEvt); 311 308 //fSigmabar->Print(""); 312 309 313 if ( SigbarOld > 0.0)310 if (sigbarold > 0) 314 311 { 315 312 //*fLog << "MPadSchweizer::Process(); Sigmabar of event to be padded is > 0 : " … … 322 319 } 323 320 324 Double_t Theta= kRad2Deg*fMcEvt->GetTelescopeTheta();321 const Double_t theta = kRad2Deg*fMcEvt->GetTelescopeTheta(); 325 322 // *fLog << "Theta = " << Theta << endl; 326 323 … … 330 327 // generate a Sigmabar according to the histogram fHSigmaTheta 331 328 // 332 Double_t Sigmabar;333 Int_t binNumber = fHSigmaTheta->GetXaxis()->FindBin( Theta);329 Double_t sigmabar=0; 330 Int_t binNumber = fHSigmaTheta->GetXaxis()->FindBin(theta); 334 331 335 332 if ( binNumber < 1 || binNumber > fHSigmaTheta->GetNbinsX() ) 336 333 { 337 334 *fLog << "MPadSchweizer::Process(); binNumber out of range : Theta, binNumber = " 338 << Theta << ", " << binNumber << "; Skip event " << endl;335 << theta << ", " << binNumber << "; Skip event " << endl; 339 336 // event cannot be padded; skip event 340 337 … … 345 342 else 346 343 { 347 TH1D* fHSigma =344 TH1D* hsigma = 348 345 fHSigmaTheta->ProjectionY("", binNumber, binNumber, ""); 349 if ( fHSigma->GetEntries() == 0.0 )346 if ( hsigma->GetEntries() == 0 ) 350 347 { 351 348 *fLog << "MPadSchweizer::Process(); projection for Theta bin " 352 349 << binNumber << " has no entries; Skip event " << endl; 353 350 // event cannot be padded; skip event 354 delete fHSigma;351 delete hsigma; 355 352 356 353 rc = 3; … … 360 357 else 361 358 { 362 Sigmabar = fHSigma->GetRandom();359 sigmabar = hsigma->GetRandom(); 363 360 364 361 //*fLog << "Theta, bin number = " << Theta << ", " << binNumber 365 362 // << ", Sigmabar = " << Sigmabar << endl; 366 363 } 367 delete fHSigma;364 delete hsigma; 368 365 } 369 366 … … 375 372 376 373 // Skip event if target Sigmabar is <= SigbarOld 377 if ( Sigmabar <= SigbarOld)374 if (sigmabar <= sigbarold) 378 375 { 379 376 *fLog << "MPadSchweizer::Process(); target Sigmabar is less than SigbarOld : " 380 << Sigmabar << ", " << SigbarOld << ", Skip event" << endl;377 << sigmabar << ", " << sigbarold << ", Skip event" << endl; 381 378 382 379 rc = 4; … … 393 390 // - using a fixed value (F2excess) for the excess noise factor 394 391 395 Double_t elNoise2; // [photons]396 392 Double_t F2excess = 1.3; 397 Double_t lambdabar; // [photons] 398 399 400 Int_t binTheta = fHDiffPixTheta->GetXaxis()->FindBin( Theta);393 394 const Double_t sigmabar2 = sigmabar*sigmabar; 395 396 Int_t binTheta = fHDiffPixTheta->GetXaxis()->FindBin(theta); 401 397 if (binTheta != binNumber) 402 398 { … … 410 406 // otherwise the electronic noise of an individual pixel (elNoise2Pix) 411 407 // may become negative 412 TH1D* fHNoise = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, 413 0, 9999, ""); 414 Double_t RMS = fHNoise->GetRMS(1); 415 delete fHNoise; 416 elNoise2 = TMath::Min(RMS, Sigmabar*Sigmabar - SigbarOld*SigbarOld); 408 TH1D* hnoise = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, 409 0, 9999, ""); 410 411 const Double_t elNoise2 = TMath::Min(hnoise->GetRMS(1), 412 sigmabar2 - sigbarold*sigbarold); // [photons] 413 414 delete hnoise; 415 417 416 //*fLog << "elNoise2 = " << elNoise2 << endl; 418 417 419 lambdabar = (Sigmabar*Sigmabar - SigbarOld*SigbarOld - elNoise2) 420 / F2excess; 418 const Double_t lambdabar = (sigmabar2 - sigbarold*sigbarold - elNoise2) 419 / F2excess; // [photons] 420 421 421 // This value of lambdabar is the same for all pixels; 422 422 // note that lambdabar is normalized to the area of pixel 0 423 424 423 425 424 //---------- start loop over pixels --------------------------------- … … 428 427 // pad only pixels - which are used (before image cleaning) 429 428 // 430 Double_t Sig = 0.0;431 Double_t Sigma2 = 0.0;432 Double_t Diff = 0.0;433 Double_t addSig2 = 0 .0;434 Double_t elNoise2Pix = 0 .0;429 Double_t sig = 0; 430 Double_t sigma2 = 0; 431 Double_t diff = 0; 432 Double_t addSig2 = 0; 433 Double_t elNoise2Pix = 0; 435 434 436 435 437 436 for (UInt_t i=0; i<npix; i++) 438 437 { 439 MCerPhotPix &pix = fEvt->operator[](i);438 MCerPhotPix &pix = (*fEvt)[i]; 440 439 if ( !pix.IsPixelUsed() ) 441 440 continue; … … 449 448 450 449 Int_t j = pix.GetPixId(); 451 Double_t Area = fCam->GetPixRatio(j);452 453 MPedestalPix &ppix = fPed->operator[](j);450 Double_t area = fCam->GetPixRatio(j); 451 452 MPedestalPix &ppix = (*fPed)[j]; 454 453 Double_t oldsigma = ppix.GetMeanRms(); 455 454 … … 460 459 461 460 Int_t count; 462 Int_t OK;463 TH1D* fHDiff;464 TH1D* fHSig;461 Bool_t ok=kFALSE; 462 TH1D* hDiff; 463 TH1D* hSig; 465 464 466 465 switch (fPadFlag) … … 468 467 case 1 : 469 468 // throw the Sigma for this pixel from the distribution fHDiffPixTheta 470 fHDiff = 471 fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, 469 hDiff = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, 472 470 binPixel, binPixel, ""); 473 if ( fHDiff->GetEntries() == 0.0 )471 if ( hDiff->GetEntries() == 0 ) 474 472 { 475 473 *fLog << "MPadSchweizer::Process(); projection for Theta bin " … … 483 481 484 482 count = 0; 485 OK = 0;486 483 for (Int_t m=0; m<20; m++) 487 484 { 488 485 count += 1; 489 Diff = fHDiff->GetRandom();486 diff = hDiff->GetRandom(); 490 487 // the following condition ensures that elNoise2Pix > 0.0 491 if ( ( Diff+Sigmabar*Sigmabar-oldsigma*oldsigma/Area-492 lambdabar*F2excess) > 0 .0)488 if ( (diff+sigmabar2-oldsigma*oldsigma/area- 489 lambdabar*F2excess) > 0 ) 493 490 { 494 OK = 1;491 ok = kTRUE; 495 492 break; 496 493 } 497 494 } 498 if ( OK == 0)495 if (!ok) 499 496 { 500 *fLog << "Theta, j, count, Sigmabar, Diff = " << Theta << ", "501 << j << ", " << count << ", " << Sigmabar << ", "502 << Diff << endl;503 Diff = lambdabar*F2excess + oldsigma*oldsigma/Area504 - Sigmabar*Sigmabar;497 *fLog << "Theta, j, count, Sigmabar, Diff = " << theta << ", " 498 << j << ", " << count << ", " << sigmabar << ", " 499 << diff << endl; 500 diff = lambdabar*F2excess + oldsigma*oldsigma/area 501 - sigmabar2; 505 502 } 506 delete fHDiff;507 Sigma2 = Diff + Sigmabar*Sigmabar;503 delete hDiff; 504 sigma2 = diff + sigmabar2; 508 505 break; 509 506 510 507 case 2 : 511 508 // throw the Sigma for this pixel from the distribution fHSigmaPixTheta 512 fHSig = 513 fHSigmaPixTheta->ProjectionZ("", binTheta, binTheta, 509 hSig = fHSigmaPixTheta->ProjectionZ("", binTheta, binTheta, 514 510 binPixel, binPixel, ""); 515 if ( fHSig->GetEntries() == 0.0 )511 if ( hSig->GetEntries() == 0 ) 516 512 { 517 513 *fLog << "MPadSchweizer::Process(); projection for Theta bin " … … 525 521 526 522 count = 0; 527 OK = 0;528 523 for (Int_t m=0; m<20; m++) 529 524 { 530 525 count += 1; 531 Sig = fHSig->GetRandom();532 Sigma2 = Sig*Sig/Area;526 sig = hSig->GetRandom(); 527 sigma2 = sig*sig/area; 533 528 // the following condition ensures that elNoise2Pix > 0.0 534 if ( ( Sigma2-oldsigma*oldsigma/Area-lambdabar*F2excess) > 0.0 )529 if ( (sigma2-oldsigma*oldsigma/area-lambdabar*F2excess) > 0.0 ) 535 530 { 536 OK = 1;531 ok = kTRUE; 537 532 break; 538 533 } 539 534 } 540 if ( OK == 0)535 if (!ok) 541 536 { 542 *fLog << "Theta, j, count, Sigmabar, Sig = " << Theta << ", "543 << j << ", " << count << ", " << Sigmabar << ", "544 << Sig << endl;545 Sigma2 = lambdabar*F2excess + oldsigma*oldsigma/Area;537 *fLog << "Theta, j, count, Sigmabar, Sig = " << theta << ", " 538 << j << ", " << count << ", " << sigmabar << ", " 539 << sig << endl; 540 sigma2 = lambdabar*F2excess + oldsigma*oldsigma/area; 546 541 } 547 delete fHSig;542 delete hSig; 548 543 break; 549 544 } … … 551 546 //--------------------------------- 552 547 // get the additional sigma^2 for this pixel (due to the padding) 553 addSig2 = Sigma2*Area - oldsigma*oldsigma;548 addSig2 = sigma2*area - oldsigma*oldsigma; 554 549 555 550 556 551 //--------------------------------- 557 552 // get the additional electronic noise for this pixel 558 elNoise2Pix = addSig2 - lambdabar*F2excess* Area;553 elNoise2Pix = addSig2 - lambdabar*F2excess*area; 559 554 560 555 … … 562 557 // throw actual number of additional NSB photons (NSB) 563 558 // and its RMS (sigmaNSB) 564 Double_t NSB0 = gRandom->Poisson(lambdabar* Area);565 Double_t arg = NSB0*(F2excess-1 .0) + elNoise2Pix;559 Double_t NSB0 = gRandom->Poisson(lambdabar*area); 560 Double_t arg = NSB0*(F2excess-1) + elNoise2Pix; 566 561 Double_t sigmaNSB0; 567 562 568 if (arg >= 0 .0)563 if (arg >= 0) 569 564 { 570 565 sigmaNSB0 = sqrt( arg ); … … 572 567 else 573 568 { 574 *fLog << "MPadSchweizer::Process(); argument of sqrt < 0 .0 : "569 *fLog << "MPadSchweizer::Process(); argument of sqrt < 0 : " 575 570 << arg << endl; 576 571 sigmaNSB0 = 0.0000001; … … 581 576 // smear NSB0 according to sigmaNSB0 582 577 // and subtract lambdabar because of AC coupling 583 Double_t NSB = gRandom->Gaus(NSB0, sigmaNSB0) - lambdabar* Area;578 Double_t NSB = gRandom->Gaus(NSB0, sigmaNSB0) - lambdabar*area; 584 579 585 580 //--------------------------------- … … 590 585 pix.SetNumPhotons( newphotons ); 591 586 592 fHNSB->Fill( NSB/sqrt(Area) ); 593 fHPhotons->Fill( oldphotons/sqrt(Area), newphotons/sqrt(Area) ); 587 const Double_t areasqrt = sqrt(area); 588 589 fHNSB->Fill( NSB/areasqrt ); 590 fHPhotons->Fill( oldphotons/areasqrt, newphotons/areasqrt ); 594 591 595 592 … … 604 601 605 602 fHSigmaPedestal->Fill( oldsigma, newsigma ); 606 fHSigPixTh-> Fill( Theta, (Double_t) j, newsigma );603 fHSigPixTh-> Fill( theta, (Double_t) j, newsigma ); 607 604 } 608 605 //---------- end of loop over pixels --------------------------------- … … 610 607 // Calculate Sigmabar again and crosscheck 611 608 612 Double_t SigbarNew = fSigmabar->Calc(*fCam, *fPed, *fEvt);609 Double_t sigbarnew = fSigmabar->Calc(*fCam, *fPed, *fEvt); 613 610 //*fLog << "after padding : " << endl; 614 611 //fSigmabar->Print(""); 615 612 616 fHSigbarTheta->Fill( Theta, SigbarNew );613 fHSigbarTheta->Fill( theta, sigbarnew ); 617 614 618 615 // this loop is only for filling the histogram fHDiffPixTh 619 616 for (UInt_t i=0; i<npix; i++) 620 617 { 621 MCerPhotPix &pix = fEvt->operator[](i);618 MCerPhotPix &pix = (*fEvt)[i]; 622 619 if ( !pix.IsPixelUsed() ) 623 620 continue; … … 631 628 632 629 Int_t j = pix.GetPixId(); 633 Double_t Area = fCam->GetPixRatio(j);634 635 MPedestalPix &ppix = fPed->operator[](j);630 Double_t area = fCam->GetPixRatio(j); 631 632 MPedestalPix &ppix = (*fPed)[j]; 636 633 Double_t newsigma = ppix.GetMeanRms(); 637 634 638 fHDiffPixTh->Fill( Theta, (Double_t) j,639 newsigma*newsigma/ Area-SigbarNew*SigbarNew);635 fHDiffPixTh->Fill( theta, (Double_t) j, 636 newsigma*newsigma/area-sigbarnew*sigbarnew); 640 637 } 641 638
Note:
See TracChangeset
for help on using the changeset viewer.