Changeset 1748
- Timestamp:
- 02/07/03 13:26:50 (22 years ago)
- Location:
- trunk/MagicSoft/Mars
- Files:
-
- 4 added
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/manalysis/MSigmabar.cc
r1682 r1748 16 16 ! 17 17 ! 18 ! Author(s): Robert Wagner 10/2002 <mailto:magicsoft@rwagner.de> 18 ! Author(s): Robert Wagner 10/2002 <mailto:magicsoft@rwagner.de> 19 ! Wolfgang Wittek 01/2003 <mailto:wittek@mppmu.mpg.de> 19 20 ! 20 21 ! Copyright: MAGIC Software Development, 2002 … … 27 28 // MSigmabar // 28 29 // // 29 // This is the storage container to hold information about the mean sigma // 30 // (aka Sigmabar) of all pedestals // 30 // This is the storage container to hold information about // 31 // "average" pedestal sigmas // 32 // // 33 // In calculating averages all sigmas are normalized to the area of pixel 0// 31 34 // // 32 35 ///////////////////////////////////////////////////////////////////////////// … … 37 40 #include "MLog.h" 38 41 #include "MLogManip.h" 39 42 #include "MParList.h" 40 43 #include "MGeomCam.h" 41 44 #include "MPedestalCam.h" 45 #include "MPedestalPix.h" 42 46 #include "MGeomPix.h" 43 #include "MPedestalPix.h" 47 #include "MCerPhotEvt.h" 48 #include "MCerPhotPix.h" 44 49 45 50 ClassImp(MSigmabar); 46 51 52 // -------------------------------------------------------------------------- 53 // 47 54 MSigmabar::MSigmabar(const char *name, const char *title) 48 55 { … … 50 57 fTitle = title ? title : "Storage container for Sigmabar"; 51 58 52 fSigmabar = 0.0;53 fSigmabarInner = 0.0;54 fSigmabarOuter = 0.0;55 fRatioA = 0.0;56 59 57 60 fCalcPixNum=kTRUE; 58 61 } 59 62 63 // -------------------------------------------------------------------------- 64 // 60 65 MSigmabar::~MSigmabar() 61 66 { 62 67 // do nothing special. 63 68 } 64 65 69 66 70 // -------------------------------------------------------------------------- … … 68 72 // Actual calculation of sigmabar. This is done for each of the six sectors 69 73 // separately due to their possibly different HV behavior. Also inner and 70 // outer pixels are treated separately 74 // outer pixels are treated separately. 71 75 // 72 76 // Preliminary! Works for CT1 test, for real MAGIC crosschecks still have … … 74 78 // determination of sector to which a respective pixel belongs 75 79 // 76 Bool_t MSigmabar::Calc(const MGeomCam &geom, const MPedestalCam &ped) 77 { 78 const UInt_t npix = ped.GetSize(); 80 Float_t MSigmabar::Calc(const MGeomCam &geom, const MPedestalCam &ped, 81 const MCerPhotEvt &evt) 82 { 83 fSigmabar = 0.0; 84 fSigmabarInner = 0.0; 85 fSigmabarOuter = 0.0; 86 87 79 88 Float_t innerSquaredSum[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; 80 89 Float_t outerSquaredSum[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; 81 Int_t innerPixels[6] = {0,0,0,0,0,0}; 82 Int_t outerPixels[6] = {0,0,0,0,0,0}; 90 Int_t innerPixels[6] = {0,0,0,0,0,0}; 91 Int_t outerPixels[6] = {0,0,0,0,0,0}; 92 83 93 Int_t currentSector; 84 94 Float_t angle; 85 95 86 for (UInt_t i=0; i<npix;i++) 87 { 88 const MGeomPix gpix = geom[i]; 89 angle=((6*atan2(gpix.GetX(),gpix.GetY())/(2*TMath::Pi()))-0.5); 90 if (angle<0) angle+=6; 91 currentSector=(Int_t)angle; 92 geom.GetPixRatio(i) == 1 ? innerPixels[currentSector]++ : outerPixels[currentSector]++; 93 94 const MPedestalPix pix = ped[i]; 95 geom.GetPixRatio(i) == 1 ? innerSquaredSum[currentSector]+=(pix.GetSigma()*pix.GetSigma()) : outerSquaredSum[currentSector]+=((pix.GetSigma()*pix.GetSigma()) / geom.GetPixRatio(i)); 96 97 // Get once and forever the ratio of areas outer/inner pixels 98 if (fRatioA && (geom.GetPixRatio(i) != 1)) fRatioA=geom.GetPixRatio(i); 99 } 100 101 // Overall Sigma 96 // sum up sigma**2 for each sector, separately for inner and outer region; 97 // all pixels are renormalized to the area of pixel 0 98 // 99 // consider all pixels with Cherenkov photon information 100 // and require "Used" 101 // 102 103 const UInt_t npix = evt.GetNumPixels(); 104 105 for (UInt_t i=0; i<npix; i++) 106 { 107 MCerPhotPix &cerpix = evt.operator[](i); 108 if (!cerpix.IsPixelUsed()) 109 continue; 110 111 if ( cerpix.GetNumPhotons() == 0 ) 112 { 113 *fLog << "MSigmabar::Calc(); no.of photons is 0 for used pixel" 114 << endl; 115 continue; 116 } 117 118 Int_t j = cerpix.GetPixId(); 119 Double_t Area = geom.GetPixRatio(j); 120 121 const MGeomPix &gpix = geom[j]; 122 const MPedestalPix &pix = ped[j]; 123 124 //angle = 6.0*atan2(gpix.GetX(),gpix.GetY()) / (2.0*TMath::Pi()) - 0.5; 125 //if (angle<0.0) angle+=6.0; 126 127 angle = 6.0*atan2(gpix.GetY(),gpix.GetX()) / (2.0*TMath::Pi()); 128 if (angle<0.0) angle+=6.0; 129 currentSector=(Int_t)angle; 130 131 // count only those pixels which have a sigma != 0.0 132 Float_t sigma = pix.GetMeanRms(); 133 134 if ( sigma != 0.0 ) 135 { 136 if (Area < 1.5) 137 { 138 innerPixels[currentSector]++; 139 innerSquaredSum[currentSector]+= sigma*sigma / Area; 140 } 141 else 142 { 143 outerPixels[currentSector]++; 144 outerSquaredSum[currentSector]+= sigma*sigma / Area; 145 } 146 } 147 } 148 102 149 fSigmabarInner=0; fInnerPixels=0; 103 150 fSigmabarOuter=0; fOuterPixels=0; 104 151 for (UInt_t i=0; i<6; i++) { 105 fSigmabarInner +=innerSquaredSum[i];106 fInnerPixels +=innerPixels[i];107 fSigmabarOuter +=outerSquaredSum[i];108 fOuterPixels +=outerPixels[i];152 fSigmabarInner += innerSquaredSum[i]; 153 fInnerPixels += innerPixels[i]; 154 fSigmabarOuter += outerSquaredSum[i]; 155 fOuterPixels += outerPixels[i]; 109 156 } 110 157 111 fSigmabarInner/=fInnerPixels; 112 if (fSigmabarOuter != 0) fSigmabarOuter/=fOuterPixels; 113 fSigmabar=sqrt(fSigmabarInner + fSigmabarOuter/( fOuterPixels==0 ? 1 : fRatioA)); 114 fSigmabarInner=sqrt(fSigmabarInner); 115 fSigmabarOuter=sqrt(fSigmabarOuter); 158 159 // this is the sqrt of the average sigma**2; 160 // 161 if ( (fInnerPixels+fOuterPixels) > 0) 162 fSigmabar=sqrt( (fSigmabarInner + fSigmabarOuter) 163 / (fInnerPixels + fOuterPixels) ); 164 165 if (fInnerPixels > 0) fSigmabarInner /= fInnerPixels; 166 if (fOuterPixels > 0) fSigmabarOuter /= fOuterPixels; 167 168 // this is the sqrt of the average sigma**2 169 // for the inner and outer pixels respectively 170 // 171 fSigmabarInner = sqrt( fSigmabarInner ); 172 fSigmabarOuter = sqrt( fSigmabarOuter ); 116 173 174 117 175 for (UInt_t i=0; i<6; i++) { 118 fSigmabarInnerSector[i]=innerSquaredSum[i]/innerPixels[i]; 119 fSigmabarOuterSector[i]=outerSquaredSum[i]/( outerSquaredSum[i]==0 ? 1: outerPixels[i] ); 120 121 fSigmabarSector[i]=sqrt(fSigmabarInnerSector[i] + fSigmabarOuterSector[i]*fRatioA); 122 fSigmabarInnerSector[i]=sqrt(fSigmabarInnerSector[i]); 123 fSigmabarOuterSector[i]=sqrt(fSigmabarOuterSector[i]); 176 fSigmabarSector[i] = 0.0; 177 fSigmabarInnerSector[i] = 0.0; 178 fSigmabarOuterSector[i] = 0.0; 179 180 if ( (innerPixels[i]+outerPixels[i]) > 0) 181 fSigmabarSector[i] = sqrt( (innerSquaredSum[i] + outerSquaredSum[i]) 182 / (innerPixels[i] + outerPixels[i] ) ); 183 if ( innerPixels[i] > 0) 184 fSigmabarInnerSector[i] = innerSquaredSum[i] / innerPixels[i]; 185 if ( outerPixels[i] > 0) 186 fSigmabarOuterSector[i] = outerSquaredSum[i] / outerPixels[i]; 187 188 189 fSigmabarInnerSector[i] = sqrt( fSigmabarInnerSector[i] ); 190 fSigmabarOuterSector[i] = sqrt( fSigmabarOuterSector[i] ); 124 191 } 125 192 126 // Did all our calculations work? fOuterPixels==0 could happen, however.127 return (fInnerPixels!=0); 128 } 129 130 193 return (fSigmabar); 194 } 195 196 // -------------------------------------------------------------------------- 197 // 131 198 void MSigmabar::Print(Option_t *) const 132 199 { 133 *fLog << all; 134 *fLog << "Sigmabar Overall " << fSigmabar; 200 *fLog << endl; 201 *fLog << "Total number of inner pixels is " << fInnerPixels << endl; 202 *fLog << "Total number of outer pixels is " << fOuterPixels << endl; 203 *fLog << endl; 204 205 *fLog << "Sigmabar Overall : " << fSigmabar << " "; 135 206 *fLog << " Sectors: "; 136 for (Int_t i=0;i<6;i++) *fLog << fSigmabarSector[i] << " "; 137 *fLog << endl; 138 *fLog << "Sigmabar Inner " << fSigmabarInner; 207 for (Int_t i=0;i<6;i++) *fLog << fSigmabarSector[i] << ", "; 208 *fLog << endl; 209 210 211 *fLog << "Sigmabar Inner : " << fSigmabarInner << " "; 139 212 *fLog << " Sectors: "; 140 for (Int_t i=0;i<6;i++) *fLog << fSigmabarInnerSector[i] << " "; 141 *fLog << endl; 142 *fLog << "Sigmabar Outer " << fSigmabarOuter; 213 for (Int_t i=0;i<6;i++) *fLog << fSigmabarInnerSector[i] << ", "; 214 *fLog << endl; 215 216 217 *fLog << "Sigmabar Outer : " << fSigmabarOuter << " "; 143 218 *fLog << " Sectors: "; 144 for (Int_t i=0;i<6;i++) *fLog << fSigmabarOuterSector[i] << " "; 145 *fLog << endl; 146 *fLog << "Total number of inner pixels found is " << fInnerPixels << endl; 147 *fLog << "Total number of outer pixels found is " << fOuterPixels << endl; 148 *fLog << "Ratio of areas outer/inner pixels found is " << fRatioA << endl; 149 *fLog << endl; 150 } 219 for (Int_t i=0;i<6;i++) *fLog << fSigmabarOuterSector[i] << ", "; 220 *fLog << endl; 221 222 } 223 // -------------------------------------------------------------------------- 224 225 226 227 -
trunk/MagicSoft/Mars/manalysis/MSigmabar.h
r1693 r1748 4 4 #ifndef MARS_MParContainer 5 5 #include "MParContainer.h" 6 #endif 7 8 #ifndef MARS_MParList 9 #include "MParList.h" 6 10 #endif 7 11 … … 14 18 #endif 15 19 20 class MCerPhotEvt; 21 16 22 class MSigmabar : public MParContainer 17 23 { 18 24 private: 19 Float_t fSigmabar; // Sigmabar (mean standard deviation) of pedestal 25 Float_t fSigmabar; // Sigmabar ("average" RMS) of pedestal 26 Float_t fSigmabarInner; // --only for inner pixels 27 Float_t fSigmabarOuter; // --only for outer pixels 28 20 29 Float_t fSigmabarSector[6]; // --for the 6 sectors of the camera 21 30 Float_t fSigmabarInnerSector[6]; 22 31 Float_t fSigmabarOuterSector[6]; 23 Float_t fSigmabarInner; // --only for inner pixels24 Float_t fSigmabarOuter; // --only for outer pixels25 UInt_t fInnerPixels; // Overall number of inner pixels26 UInt_t fOuterPixels; // Overall number of outer pixels27 Float_t fRatioA; // Ratio of areas (outer/inner pixels)28 Bool_t fCalcPixNum;29 32 33 UInt_t fInnerPixels; // Overall number of inner pixels 34 UInt_t fOuterPixels; // Overall number of outer pixels 35 36 Bool_t fCalcPixNum; 37 30 38 public: 31 39 … … 46 54 // void SetSigmabarOuter(Float_t f) { fSigmabarOuter = f; } 47 55 48 Bool_t Calc(const MGeomCam &geom, const MPedestalCam &ped);56 Float_t Calc(const MGeomCam &geom, const MPedestalCam &ped, const MCerPhotEvt &evt); 49 57 50 58 ClassDef(MSigmabar, 1) // Storage Container for Sigmabar … … 53 61 #endif 54 62 63 64 65 66 -
trunk/MagicSoft/Mars/manalysis/MSigmabarCalc.cc
r1737 r1748 110 110 if (!fMcEvt) 111 111 { 112 *fLog << err << dbginf << "MHSigmabarTheta : MMcEvt not found... aborting." << endl; 112 *fLog << err << dbginf << "MHSigmabaCalc : MMcEvt not found... aborting." << endl; 113 return kFALSE; 114 } 115 116 fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt"); 117 if (!fEvt) 118 { 119 *fLog << err << dbginf << "MHSigmabarCalc : MCerPhotEvt not found... aborting." << endl; 113 120 return kFALSE; 114 121 } … … 126 133 Bool_t MSigmabarCalc::Process() 127 134 { 128 Bool_t rc = fSig->Calc(*fCam, *fPed);129 fSigmabarMax = TMath::Max((Double_t) fSig->GetSigmabar(), fSigmabarMax);130 fSigmabarMin = TMath::Min((Double_t) fSig->GetSigmabar(), fSigmabarMin);135 Float_t rc = fSig->Calc(*fCam, *fPed, *fEvt); 136 fSigmabarMax = TMath::Max((Double_t)rc, fSigmabarMax); 137 fSigmabarMin = TMath::Min((Double_t)rc, fSigmabarMin); 131 138 132 if (fMcEvt->GetTelescopeTheta()*kRad2Deg < 90)139 if (fMcEvt->GetTelescopeTheta()*kRad2Deg < 120) 133 140 fThetaMax = TMath::Max(fMcEvt->GetTelescopeTheta()*kRad2Deg, fThetaMax); 134 141 fThetaMin = TMath::Min(fMcEvt->GetTelescopeTheta()*kRad2Deg, fThetaMin); 135 142 136 return rc;143 return kTRUE; 137 144 } 138 145 -
trunk/MagicSoft/Mars/manalysis/MSigmabarCalc.h
r1737 r1748 43 43 MSigmabarParam *fSigParam; 44 44 MMcEvt *fMcEvt; 45 MCerPhotEvt *fEvt; 45 46 void Reset(); 46 47 … … 58 59 #endif 59 60 61 62
Note:
See TracChangeset
for help on using the changeset viewer.