Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 2011) +++ trunk/MagicSoft/Mars/Changelog (revision 2012) @@ -15,4 +15,22 @@ * mhistmc/MHMcEnergyMigration.h: - fixed the includes + + * mgui/MCamDisplay.cc: + - changed autoscaling (max<1:max=1 --> max==min:max=min+1) + + * manalysis/MBlindPixelCalc.cc: + - interpolate: take pixel area into account + + * mhist/MHSigmaTheta.h: + - removed nonsense GetSigmaThetaByName(const TString name) + - removed nonsense GetSigmaPixThetaByName(const TString name) + - removed nonsense GetDiffPixThetaByName(const TString name) + + * manalysis/MPadSchweizer.cc: + - fixed naming + - fixed usage of operators + - added some const qualifiers + - replaced 'int OK' by 'Bool_t ok' + - fixed wrong usage floating point value 0 Index: trunk/MagicSoft/Mars/manalysis/MPadSchweizer.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPadSchweizer.cc (revision 2011) +++ trunk/MagicSoft/Mars/manalysis/MPadSchweizer.cc (revision 2012) @@ -280,12 +280,9 @@ //*fLog << "Entry MPadSchweizer::Process();" << endl; - Int_t rc; + Int_t rc=0; const UInt_t npix = fEvt->GetNumPixels(); - Double_t SigbarOld; - - - SigbarOld = fSigmabar->Calc(*fCam, *fPed, *fEvt); + fSigmabar->Calc(*fCam, *fPed, *fEvt); //*fLog << "before padding : " << endl; //fSigmabar->Print(""); @@ -308,8 +305,8 @@ // Calculate average sigma of the event // - SigbarOld = fSigmabar->Calc(*fCam, *fPed, *fEvt); + Double_t sigbarold = fSigmabar->Calc(*fCam, *fPed, *fEvt); //fSigmabar->Print(""); - if (SigbarOld > 0.0) + if (sigbarold > 0) { //*fLog << "MPadSchweizer::Process(); Sigmabar of event to be padded is > 0 : " @@ -322,5 +319,5 @@ } - Double_t Theta = kRad2Deg*fMcEvt->GetTelescopeTheta(); + const Double_t theta = kRad2Deg*fMcEvt->GetTelescopeTheta(); // *fLog << "Theta = " << Theta << endl; @@ -330,11 +327,11 @@ // generate a Sigmabar according to the histogram fHSigmaTheta // - Double_t Sigmabar; - Int_t binNumber = fHSigmaTheta->GetXaxis()->FindBin(Theta); + Double_t sigmabar=0; + Int_t binNumber = fHSigmaTheta->GetXaxis()->FindBin(theta); if ( binNumber < 1 || binNumber > fHSigmaTheta->GetNbinsX() ) { *fLog << "MPadSchweizer::Process(); binNumber out of range : Theta, binNumber = " - << Theta << ", " << binNumber << "; Skip event " << endl; + << theta << ", " << binNumber << "; Skip event " << endl; // event cannot be padded; skip event @@ -345,12 +342,12 @@ else { - TH1D* fHSigma = + TH1D* hsigma = fHSigmaTheta->ProjectionY("", binNumber, binNumber, ""); - if ( fHSigma->GetEntries() == 0.0 ) + if ( hsigma->GetEntries() == 0 ) { *fLog << "MPadSchweizer::Process(); projection for Theta bin " << binNumber << " has no entries; Skip event " << endl; // event cannot be padded; skip event - delete fHSigma; + delete hsigma; rc = 3; @@ -360,10 +357,10 @@ else { - Sigmabar = fHSigma->GetRandom(); + sigmabar = hsigma->GetRandom(); //*fLog << "Theta, bin number = " << Theta << ", " << binNumber // << ", Sigmabar = " << Sigmabar << endl; } - delete fHSigma; + delete hsigma; } @@ -375,8 +372,8 @@ // Skip event if target Sigmabar is <= SigbarOld - if (Sigmabar <= SigbarOld) + if (sigmabar <= sigbarold) { *fLog << "MPadSchweizer::Process(); target Sigmabar is less than SigbarOld : " - << Sigmabar << ", " << SigbarOld << ", Skip event" << endl; + << sigmabar << ", " << sigbarold << ", Skip event" << endl; rc = 4; @@ -393,10 +390,9 @@ // - using a fixed value (F2excess) for the excess noise factor - Double_t elNoise2; // [photons] Double_t F2excess = 1.3; - Double_t lambdabar; // [photons] - - - Int_t binTheta = fHDiffPixTheta->GetXaxis()->FindBin(Theta); + + const Double_t sigmabar2 = sigmabar*sigmabar; + + Int_t binTheta = fHDiffPixTheta->GetXaxis()->FindBin(theta); if (binTheta != binNumber) { @@ -410,16 +406,19 @@ // otherwise the electronic noise of an individual pixel (elNoise2Pix) // may become negative - TH1D* fHNoise = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, - 0, 9999, ""); - Double_t RMS = fHNoise->GetRMS(1); - delete fHNoise; - elNoise2 = TMath::Min(RMS, Sigmabar*Sigmabar - SigbarOld*SigbarOld); + TH1D* hnoise = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, + 0, 9999, ""); + + const Double_t elNoise2 = TMath::Min(hnoise->GetRMS(1), + sigmabar2 - sigbarold*sigbarold); // [photons] + + delete hnoise; + //*fLog << "elNoise2 = " << elNoise2 << endl; - lambdabar = (Sigmabar*Sigmabar - SigbarOld*SigbarOld - elNoise2) - / F2excess; + const Double_t lambdabar = (sigmabar2 - sigbarold*sigbarold - elNoise2) + / F2excess; // [photons] + // This value of lambdabar is the same for all pixels; // note that lambdabar is normalized to the area of pixel 0 - //---------- start loop over pixels --------------------------------- @@ -428,14 +427,14 @@ // pad only pixels - which are used (before image cleaning) // - Double_t Sig = 0.0; - Double_t Sigma2 = 0.0; - Double_t Diff = 0.0; - Double_t addSig2 = 0.0; - Double_t elNoise2Pix = 0.0; + Double_t sig = 0; + Double_t sigma2 = 0; + Double_t diff = 0; + Double_t addSig2 = 0; + Double_t elNoise2Pix = 0; for (UInt_t i=0; ioperator[](i); + MCerPhotPix &pix = (*fEvt)[i]; if ( !pix.IsPixelUsed() ) continue; @@ -449,7 +448,7 @@ Int_t j = pix.GetPixId(); - Double_t Area = fCam->GetPixRatio(j); - - MPedestalPix &ppix = fPed->operator[](j); + Double_t area = fCam->GetPixRatio(j); + + MPedestalPix &ppix = (*fPed)[j]; Double_t oldsigma = ppix.GetMeanRms(); @@ -460,7 +459,7 @@ Int_t count; - Int_t OK; - TH1D* fHDiff; - TH1D* fHSig; + Bool_t ok=kFALSE; + TH1D* hDiff; + TH1D* hSig; switch (fPadFlag) @@ -468,8 +467,7 @@ case 1 : // throw the Sigma for this pixel from the distribution fHDiffPixTheta - fHDiff = - fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, + hDiff = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, binPixel, binPixel, ""); - if ( fHDiff->GetEntries() == 0.0 ) + if ( hDiff->GetEntries() == 0 ) { *fLog << "MPadSchweizer::Process(); projection for Theta bin " @@ -483,35 +481,33 @@ count = 0; - OK = 0; for (Int_t m=0; m<20; m++) { count += 1; - Diff = fHDiff->GetRandom(); + diff = hDiff->GetRandom(); // the following condition ensures that elNoise2Pix > 0.0 - if ( (Diff+Sigmabar*Sigmabar-oldsigma*oldsigma/Area- - lambdabar*F2excess) > 0.0 ) + if ( (diff+sigmabar2-oldsigma*oldsigma/area- + lambdabar*F2excess) > 0 ) { - OK = 1; + ok = kTRUE; break; } } - if (OK == 0) + if (!ok) { - *fLog << "Theta, j, count, Sigmabar, Diff = " << Theta << ", " - << j << ", " << count << ", " << Sigmabar << ", " - << Diff << endl; - Diff = lambdabar*F2excess + oldsigma*oldsigma/Area - - Sigmabar*Sigmabar; + *fLog << "Theta, j, count, Sigmabar, Diff = " << theta << ", " + << j << ", " << count << ", " << sigmabar << ", " + << diff << endl; + diff = lambdabar*F2excess + oldsigma*oldsigma/area + - sigmabar2; } - delete fHDiff; - Sigma2 = Diff + Sigmabar*Sigmabar; + delete hDiff; + sigma2 = diff + sigmabar2; break; case 2 : // throw the Sigma for this pixel from the distribution fHSigmaPixTheta - fHSig = - fHSigmaPixTheta->ProjectionZ("", binTheta, binTheta, + hSig = fHSigmaPixTheta->ProjectionZ("", binTheta, binTheta, binPixel, binPixel, ""); - if ( fHSig->GetEntries() == 0.0 ) + if ( hSig->GetEntries() == 0 ) { *fLog << "MPadSchweizer::Process(); projection for Theta bin " @@ -525,25 +521,24 @@ count = 0; - OK = 0; for (Int_t m=0; m<20; m++) { count += 1; - Sig = fHSig->GetRandom(); - Sigma2 = Sig*Sig/Area; + sig = hSig->GetRandom(); + sigma2 = sig*sig/area; // the following condition ensures that elNoise2Pix > 0.0 - if ( (Sigma2-oldsigma*oldsigma/Area-lambdabar*F2excess) > 0.0 ) + if ( (sigma2-oldsigma*oldsigma/area-lambdabar*F2excess) > 0.0 ) { - OK = 1; + ok = kTRUE; break; } } - if (OK == 0) + if (!ok) { - *fLog << "Theta, j, count, Sigmabar, Sig = " << Theta << ", " - << j << ", " << count << ", " << Sigmabar << ", " - << Sig << endl; - Sigma2 = lambdabar*F2excess + oldsigma*oldsigma/Area; + *fLog << "Theta, j, count, Sigmabar, Sig = " << theta << ", " + << j << ", " << count << ", " << sigmabar << ", " + << sig << endl; + sigma2 = lambdabar*F2excess + oldsigma*oldsigma/area; } - delete fHSig; + delete hSig; break; } @@ -551,10 +546,10 @@ //--------------------------------- // get the additional sigma^2 for this pixel (due to the padding) - addSig2 = Sigma2*Area - oldsigma*oldsigma; + addSig2 = sigma2*area - oldsigma*oldsigma; //--------------------------------- // get the additional electronic noise for this pixel - elNoise2Pix = addSig2 - lambdabar*F2excess*Area; + elNoise2Pix = addSig2 - lambdabar*F2excess*area; @@ -562,9 +557,9 @@ // throw actual number of additional NSB photons (NSB) // and its RMS (sigmaNSB) - Double_t NSB0 = gRandom->Poisson(lambdabar*Area); - Double_t arg = NSB0*(F2excess-1.0) + elNoise2Pix; + Double_t NSB0 = gRandom->Poisson(lambdabar*area); + Double_t arg = NSB0*(F2excess-1) + elNoise2Pix; Double_t sigmaNSB0; - if (arg >= 0.0) + if (arg >= 0) { sigmaNSB0 = sqrt( arg ); @@ -572,5 +567,5 @@ else { - *fLog << "MPadSchweizer::Process(); argument of sqrt < 0.0 : " + *fLog << "MPadSchweizer::Process(); argument of sqrt < 0 : " << arg << endl; sigmaNSB0 = 0.0000001; @@ -581,5 +576,5 @@ // smear NSB0 according to sigmaNSB0 // and subtract lambdabar because of AC coupling - Double_t NSB = gRandom->Gaus(NSB0, sigmaNSB0) - lambdabar*Area; + Double_t NSB = gRandom->Gaus(NSB0, sigmaNSB0) - lambdabar*area; //--------------------------------- @@ -590,6 +585,8 @@ pix.SetNumPhotons( newphotons ); - fHNSB->Fill( NSB/sqrt(Area) ); - fHPhotons->Fill( oldphotons/sqrt(Area), newphotons/sqrt(Area) ); + const Double_t areasqrt = sqrt(area); + + fHNSB->Fill( NSB/areasqrt ); + fHPhotons->Fill( oldphotons/areasqrt, newphotons/areasqrt ); @@ -604,5 +601,5 @@ fHSigmaPedestal->Fill( oldsigma, newsigma ); - fHSigPixTh-> Fill( Theta, (Double_t) j, newsigma ); + fHSigPixTh-> Fill( theta, (Double_t) j, newsigma ); } //---------- end of loop over pixels --------------------------------- @@ -610,14 +607,14 @@ // Calculate Sigmabar again and crosscheck - Double_t SigbarNew = fSigmabar->Calc(*fCam, *fPed, *fEvt); + Double_t sigbarnew = fSigmabar->Calc(*fCam, *fPed, *fEvt); //*fLog << "after padding : " << endl; //fSigmabar->Print(""); - fHSigbarTheta->Fill( Theta, SigbarNew ); + fHSigbarTheta->Fill( theta, sigbarnew ); // this loop is only for filling the histogram fHDiffPixTh for (UInt_t i=0; ioperator[](i); + MCerPhotPix &pix = (*fEvt)[i]; if ( !pix.IsPixelUsed() ) continue; @@ -631,11 +628,11 @@ Int_t j = pix.GetPixId(); - Double_t Area = fCam->GetPixRatio(j); - - MPedestalPix &ppix = fPed->operator[](j); + Double_t area = fCam->GetPixRatio(j); + + MPedestalPix &ppix = (*fPed)[j]; Double_t newsigma = ppix.GetMeanRms(); - fHDiffPixTh->Fill( Theta, (Double_t) j, - newsigma*newsigma/Area-SigbarNew*SigbarNew); + fHDiffPixTh->Fill( theta, (Double_t) j, + newsigma*newsigma/area-sigbarnew*sigbarnew); }