Changeset 3643 for trunk/MagicSoft/Mars/manalysis
- Timestamp:
- 04/04/04 17:56:14 (21 years ago)
- Location:
- trunk/MagicSoft/Mars/manalysis
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/manalysis/MHPedestalCam.cc
r3446 r3643 28 28 // MHPedestalCam // 29 29 // // 30 // Hold the Pedestal information for all pixels in the camera // 31 // // 30 // Fills the extracted pedestals of MExtractedSignalCam into 31 // the MHGausEvents-classes MHPedestalPix for every: 32 // 33 // - Pixel, stored in the TObjArray's MHCalibrationCam::fHiGainArray 34 // or MHCalibrationCam::fHiGainArray, respectively, depending if 35 // MArrivalTimePix::IsLoGainUsed() is set. 36 // 37 // - Average pixel per AREA index (e.g. inner and outer for the MAGIC camera), 38 // stored in the TObjArray's MHCalibrationCam::fAverageHiGainAreas and 39 // MHCalibrationCam::fAverageHiGainAreas 40 // 41 // - Average pixel per camera SECTOR (e.g. sectors 1-6 for the MAGIC camera), 42 // stored in the TObjArray's MHCalibrationCam::fAverageHiGainSectors 43 // and MHCalibrationCam::fAverageHiGainSectors 44 // 45 // Every pedestal is directly taken from MExtractedSignalCam, filled by the 46 // appropriate extractor. 47 // 48 // The pedestals are filled into a histogram and an array, in order to perform 49 // a Fourier analysis (see MHGausEvents). The signals are moreover averaged on an 50 // event-by-event basis and written into the corresponding average pixels. 51 // 52 // The histograms are fitted to a Gaussian, mean and sigma with its errors 53 // and the fit probability are extracted. If none of these values are NaN's and 54 // if the probability is bigger than MHGausEvents::fProbLimit (default: 0.5%), 55 // the fit is declared valid. 56 // Otherwise, the fit is repeated within ranges of the previous mean 57 // +- MHGausEvents::fPickupLimit (default: 5) sigma (see MHGausEvents::RepeatFit()) 58 // In case this does not make the fit valid, the histogram means and RMS's are 59 // taken directly (see MHGausEvents::BypassFit()). 60 // 61 // Outliers of more than MHGausEvents::fPickupLimit (default: 5) sigmas 62 // from the mean are counted as Pickup events (stored in MHGausEvents::fPickup) 63 // 64 // The number of summed FADC slices is taken to re-normalize 65 // the result to a pedestal per pixel with the formulas (see MHPedestalPix::Renorm()): 66 // - Mean Pedestal / slice = Mean Pedestal / Number slices 67 // - Mean Pedestal Error / slice = Mean Pedestal Error / Number slices 68 // - Sigma Pedestal / slice = Sigma Pedestal / Sqrt (Number slices) 69 // - Sigma Pedestal Error / slice = Sigma Pedestal Error / Sqrt (Number slices) 70 // 71 // The class also fills arrays with the signal vs. event number, creates a fourier 72 // spectrum (see MHGausEvents::CreateFourierSpectrum()) and investigates if the 73 // projected fourier components follow an exponential distribution. 74 // 75 // This same procedure is performed for the average pixels. 76 // 77 // The following results are written into MPedestalCam: 78 // 79 // - MCalibrationPix::SetMean() 80 // - MCalibrationPix::SetMeanErr() 81 // - MCalibrationPix::SetSigma() 82 // - MCalibrationPix::SetSigmaErr() 83 // - MCalibrationPix::SetProb() 84 // - MCalibrationPix::SetNumPickup() 85 // 86 // For all averaged areas, the fitted sigma is multiplied with the square root of 87 // the number involved pixels in order to be able to compare it to the average of 88 // sigmas in the camera. 89 // 32 90 ///////////////////////////////////////////////////////////////////////////// 33 91 #include "MHPedestalCam.h" … … 46 104 #include "MPedestalPix.h" 47 105 106 #include "MGeomCam.h" 107 #include "MGeomPix.h" 108 48 109 ClassImp(MHPedestalCam); 49 110 … … 51 112 // -------------------------------------------------------------------------- 52 113 // 53 // Default constructor. Creates a MHPedestalPix object for each pixel 114 // Default constructor. 115 // 116 // Initializes: 117 // - fExtractHiGainSlices to 1. 118 // - fExtractLoGainSlices to 1. 54 119 // 55 120 MHPedestalCam::MHPedestalCam(const char *name, const char *title) 56 : fExtract Slices(0.)121 : fExtractHiGainSlices(1.), fExtractLoGainSlices(1.) 57 122 { 58 123 fName = name ? name : "MHPedestalCam"; 59 124 fTitle = title ? title : ""; 60 125 61 // 62 // loop over all Pixels and create two histograms 63 // one for the Low and one for the High gain 64 // connect all the histogram with the container fHist 65 // 66 fArray = new TObjArray; 67 fArray->SetOwner(); 68 69 } 126 } 127 128 70 129 71 130 // -------------------------------------------------------------------------- 72 131 // 73 // Delete the array conatining the pixel pedest information 74 // 75 MHPedestalCam::~MHPedestalCam() 76 { 77 delete fArray; 78 } 79 80 // -------------------------------------------------------------------------- 81 // 82 // Get i-th pixel (pixel number) 83 // 84 MHPedestalPix &MHPedestalCam::operator[](UInt_t i) 85 { 86 return *(MHPedestalPix*)(fArray->At(i)); 87 } 88 89 // -------------------------------------------------------------------------- 90 // 91 // Get i-th pixel (pixel number) 92 // 93 const MHPedestalPix &MHPedestalCam::operator[](UInt_t i) const 94 { 95 return *(MHPedestalPix*)(fArray->At(i)); 96 } 97 98 99 // -------------------------------------------------------------------------- 100 // 101 // Our own clone function is necessary since root 3.01/06 or Mars 0.4 102 // I don't know the reason 103 // 104 TObject *MHPedestalCam::Clone(const char *) const 105 { 106 107 const Int_t n = fArray->GetSize(); 108 109 // 110 // FIXME, this might be done faster and more elegant, by direct copy. 111 // 112 MHPedestalCam *cam = new MHPedestalCam; 113 114 cam->fArray->Expand(n); 115 116 for (int i=0; i<n; i++) 117 { 118 delete (*cam->fArray)[i]; 119 (*cam->fArray)[i] = (*fArray)[i]->Clone(); 120 } 121 122 return cam; 123 } 124 125 126 127 // -------------------------------------------------------------------------- 128 // 129 // To setup the object we get the number of pixels from a MGeomCam object 130 // in the Parameter list. 131 // MHPedestalPix sets its parameters to 0. (other than default which is -1.) 132 // 133 Bool_t MHPedestalCam::SetupFill(const MParList *pList) 134 { 132 // Searches pointer to: 133 // - MPedestalCam 134 // - MExtractedSignalCam 135 // 136 // Retrieves from MExtractedSignalCam: 137 // - number of used High Gain FADC slices (MExtractedSignalCam::GetNumUsedHiGainFADCSlices()) 138 // - number of used Low Gain FADC slices (MExtractedSignalCam::GetNumUsedLoGainFADCSlices()) 139 // 140 // Initializes, if empty to MGeomCam::GetNumPixels(): 141 // - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray 142 // 143 // Initializes, if empty to MGeomCam::GetNumAreas() for: 144 // - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas 145 // 146 // Initializes, if empty to MGeomCam::GetNumSectors() for: 147 // - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors 148 // 149 // Calls MHCalibrationCam::InitPedHists() for every entry in: 150 // - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray 151 // - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas 152 // - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors 153 // 154 // Sets Titles and Names for the Histograms 155 // - MHCalibrationCam::fAverageHiGainAreas 156 // - MHCalibrationCam::fAverageHiGainSectors 157 // 158 // Sets number of bins to MHCalibrationCam::fAverageNbins for: 159 // - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas 160 // - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors 161 // 162 Bool_t MHPedestalCam::ReInitHists(MParList *pList) 163 { 164 165 MExtractedSignalCam *signal = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam"); 166 if (!signal) 167 { 168 gLog << err << "Cannot find MExtractedSignalCam... abort." << endl; 169 return kFALSE; 170 } 171 135 172 136 173 fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam"); … … 142 179 } 143 180 144 fArray->Delete(); 181 Float_t sliceshi = signal->GetNumUsedHiGainFADCSlices(); 182 Float_t sliceslo = signal->GetNumUsedLoGainFADCSlices(); 183 184 if (sliceshi == 0.) 185 { 186 gLog << err << "Number of used signal slices in MExtractedSignalCam is zero ... abort." 187 << endl; 188 return kFALSE; 189 } 190 191 if (fExtractHiGainSlices != 0. && sliceshi != fExtractHiGainSlices ) 192 { 193 gLog << err << "Number of used High Gain signal slices changed in MExtractedSignalCam ... abort." 194 << endl; 195 return kFALSE; 196 } 197 198 if (fExtractLoGainSlices != 0. && sliceslo != fExtractLoGainSlices ) 199 { 200 gLog << err << "Number of used Low Gain signal slices changed in MExtractedSignalCam ... abort." 201 << endl; 202 return kFALSE; 203 } 204 205 fExtractHiGainSlices = sliceshi; 206 fExtractLoGainSlices = sliceslo; 207 208 const Int_t npixels = fGeom->GetNumPixels(); 209 const Int_t nsectors = fGeom->GetNumSectors(); 210 const Int_t nareas = fGeom->GetNumAreas(); 211 212 if (fHiGainArray->GetEntries()==0) 213 { 214 fHiGainArray->Expand(npixels); 215 for (Int_t i=0; i<npixels; i++) 216 { 217 (*fHiGainArray)[i] = new MHPedestalPix; 218 InitPedHists((MHPedestalPix&)(*this)[i],i,fExtractHiGainSlices); 219 } 220 } 221 222 if (fLoGainArray->GetEntries()==0) 223 { 224 fLoGainArray->Expand(npixels); 225 for (Int_t i=0; i<npixels; i++) 226 { 227 (*fLoGainArray)[i] = new MHPedestalPix; 228 InitPedHists((MHPedestalPix&)(*this)(i),i,fExtractLoGainSlices); 229 } 230 } 231 232 if (fAverageHiGainAreas->GetEntries()==0) 233 { 234 fAverageHiGainAreas->Expand(nareas); 235 236 for (Int_t j=0; j<nareas; j++) 237 { 238 (*fAverageHiGainAreas)[j] = 239 new MHPedestalPix("AverageHiGainArea", 240 "Average Pedestals area idx "); 241 242 InitPedHists((MHPedestalPix&)GetAverageHiGainArea(j),j,fExtractHiGainSlices); 243 244 GetAverageHiGainArea(j).GetHGausHist()->SetTitle("Pedestals average Area Idx "); 245 GetAverageHiGainArea(j).SetNbins(fAverageNbins); 246 } 247 } 248 249 if (fAverageLoGainAreas->GetEntries()==0) 250 { 251 fAverageLoGainAreas->Expand(nareas); 252 253 for (Int_t j=0; j<nareas; j++) 254 { 255 (*fAverageLoGainAreas)[j] = 256 new MHPedestalPix("AverageLoGainArea", 257 "Pedestals average Area idx "); 258 259 InitPedHists((MHPedestalPix&)GetAverageLoGainArea(j),j,fExtractLoGainSlices); 260 261 GetAverageLoGainArea(j).GetHGausHist()->SetTitle("Pedestals average Area Idx "); 262 GetAverageLoGainArea(j).SetNbins(fAverageNbins); 263 } 264 } 265 266 if (fAverageHiGainSectors->GetEntries()==0) 267 { 268 fAverageHiGainSectors->Expand(nsectors); 269 270 for (Int_t j=0; j<nsectors; j++) 271 { 272 (*fAverageHiGainSectors)[j] = 273 new MHPedestalPix("AverageHiGainSector", 274 "Pedestals average sector "); 275 276 InitPedHists((MHPedestalPix&)GetAverageHiGainSector(j),j,fExtractHiGainSlices); 277 278 GetAverageHiGainSector(j).GetHGausHist()->SetTitle("Pedestals average Sector "); 279 GetAverageHiGainSector(j).SetNbins(fAverageNbins); 280 } 281 } 282 283 if (fAverageLoGainSectors->GetEntries()==0) 284 { 285 fAverageLoGainSectors->Expand(nsectors); 286 287 for (Int_t j=0; j<nsectors; j++) 288 { 289 (*fAverageLoGainSectors)[j] = 290 new MHPedestalPix("AverageLoGainSector", 291 "Pedestals average sector "); 292 293 InitPedHists((MHPedestalPix&)GetAverageLoGainSector(j),j,fExtractLoGainSlices); 294 295 GetAverageLoGainSector(j).GetHGausHist()->SetTitle("Pedestals average Sector "); 296 GetAverageLoGainSector(j).SetNbins(fAverageNbins); 297 } 298 } 299 145 300 return kTRUE; 146 301 147 302 } 148 303 149 // -------------------------------------------------------------------------- 150 // 151 Bool_t MHPedestalCam::Fill(const MParContainer *par, const Stat_t w) 304 305 // ------------------------------------------------------------- 306 // 307 // If MBadPixelsPix::IsBad(): 308 // - calls MHGausEvents::SetExcluded() 309 // 310 // Calls: 311 // - MHGausEvents::InitBins() 312 // - MHGausEvents::ChangeHistId(i) 313 // - MHGausEvents::SetEventFrequency(fPulserFrequency) 314 // - MHPedestalPix::SetNSlices(nslices) 315 // 316 void MHPedestalCam::InitPedHists(MHPedestalPix &hist, const Int_t i, const Float_t nslices) 317 { 318 319 hist.InitBins(); 320 hist.ChangeHistId(i); 321 hist.SetEventFrequency(fPulserFrequency); 322 hist.SetNSlices(nslices); 323 } 324 325 326 327 // ------------------------------------------------------------------------------- 328 // 329 // Retrieves pointer to MExtractedSignalCam: 330 // 331 // Retrieves from MGeomCam: 332 // - number of pixels 333 // - number of pixel areas 334 // - number of sectors 335 // 336 // Fills HiGain or LoGain histograms (MHGausEvents::FillHistAndArray()), respectively 337 // with the signals MExtractedSignalCam::GetExtractedSignalHiGain() and 338 // MExtractedSignalCam::GetExtractedSignalLoGain(), respectively. 339 // 340 Bool_t MHPedestalCam::FillHists(const MParContainer *par, const Stat_t w) 152 341 { 153 342 … … 159 348 } 160 349 161 162 Float_t slices = signal->GetNumUsedHiGainFADCSlices(); 163 164 if (slices == 0.) 165 { 166 gLog << err << "Number of used signal slices in MExtractedSignalCam is zero ... abort." 167 << endl; 168 return kFALSE; 169 } 170 171 if (fExtractSlices != 0. && slices != fExtractSlices ) 172 { 173 gLog << err << "Number of used signal slices changed in MExtractedSignalCam ... abort." 174 << endl; 175 return kFALSE; 176 } 177 178 fExtractSlices = slices; 179 180 const Int_t n = signal->GetSize(); 181 182 if (fArray->GetEntries()==0) 183 { 184 fArray->Expand(n); 185 186 for (Int_t i=0; i<n; i++) 187 { 188 (*fArray)[i] = new MHPedestalPix; 189 MPedestalPix &pix = (*fPedestals)[i]; 190 pix.InitUseHists(); 191 MHPedestalPix &hist = (*this)[i]; 192 hist.ChangeHistId(i); 193 hist.InitBins(); 194 } 195 } 196 197 if (fArray->GetEntries() != n) 198 { 199 gLog << err << "ERROR - Size mismatch... abort." << endl; 200 return kFALSE; 201 } 202 203 for (Int_t i=0; i<n; i++) 204 { 350 const Int_t npixels = fGeom->GetNumPixels(); 351 const Int_t nareas = fGeom->GetNumAreas(); 352 const Int_t nsectors = fGeom->GetNumSectors(); 353 354 Float_t sumareahi [nareas], sumarealo [nareas]; 355 Float_t sumsectorhi[nsectors], sumsectorlo[nsectors]; 356 Int_t numareahi [nareas], numarealo [nareas]; 357 Int_t numsectorhi[nsectors], numsectorlo[nsectors]; 358 359 for (UInt_t j=0; j<nareas; j++) 360 { 361 sumareahi[j] = sumarealo[j] = 0.; 362 numareahi[j] = numarealo[j] = 0; 363 } 364 365 for (UInt_t j=0; j<nsectors; j++) 366 { 367 sumsectorhi[j] = sumsectorlo[j] = 0.; 368 numsectorhi[j] = numsectorlo[j] = 0; 369 } 370 371 372 for (Int_t i=0; i<npixels; i++) 373 { 374 375 MHGausEvents &histhi = (*this)[i]; 376 MHGausEvents &histlo = (*this)(i); 377 378 if (histhi.IsExcluded()) 379 continue; 205 380 206 381 const MExtractedSignalPix &pix = (*signal)[i]; 207 382 208 const Float_t signal = pix.GetExtractedSignalHiGain(); 209 210 MHPedestalPix &hist = (*this)[i]; 211 // 212 // Don't fill signal per slice, we get completely screwed up 213 // with the sigma. Better fill like it is and renorm later 214 // 215 hist.FillHistAndArray(signal); 383 const Float_t pedhi = pix.GetExtractedSignalHiGain(); 384 const Float_t pedlo = pix.GetExtractedSignalLoGain(); 385 386 const Int_t aidx = (*fGeom)[i].GetAidx(); 387 const Int_t sector = (*fGeom)[i].GetSector(); 388 389 histhi.FillHistAndArray(pedhi) ; 390 sumareahi [aidx] += pedhi; 391 numareahi [aidx] ++; 392 sumsectorhi[sector] += pedhi; 393 numsectorhi[sector] ++; 394 395 histlo.FillHistAndArray(pedlo); 396 sumarealo [aidx] += pedlo; 397 numarealo [aidx] ++; 398 sumsectorlo[sector] += pedlo; 399 numsectorlo[sector] ++; 400 401 } 402 403 for (UInt_t j=0; j<nareas; j++) 404 { 405 MHGausEvents &histhi = GetAverageHiGainArea(j); 406 histhi.FillHistAndArray(numareahi[j] == 0 ? 0. : sumareahi[j]/numareahi[j]); 407 408 MHGausEvents &histlo = GetAverageLoGainArea(j); 409 histlo.FillHistAndArray(numarealo[j] == 0 ? 0. : sumarealo[j]/numarealo[j]); 410 } 411 412 for (UInt_t j=0; j<nsectors; j++) 413 { 414 MHGausEvents &histhi = GetAverageHiGainSector(j); 415 histhi.FillHistAndArray(numsectorhi[j] == 0 ? 0. : sumsectorhi[j]/numsectorhi[j]); 416 417 MHGausEvents &histlo = GetAverageLoGainSector(j); 418 histlo.FillHistAndArray(numsectorlo[j] == 0 ? 0. : sumsectorlo[j]/numsectorlo[j]); 216 419 } 217 420 … … 219 422 } 220 423 221 Bool_t MHPedestalCam::Finalize()222 {223 for (Int_t i=0; i<fArray->GetSize(); i++)224 {225 226 MHPedestalPix &hist = (*this)[i];227 228 //229 // 1) Return if the charge distribution is already succesfully fitted230 // or if the histogram is empty231 //232 if (hist.IsGausFitOK() || hist.IsEmpty())233 continue;234 235 //236 // 2) Fit the Hi Gain histograms with a Gaussian237 //238 hist.FitGaus();239 hist.Renorm(fExtractSlices);240 hist.CreateFourierSpectrum();241 242 }243 return kTRUE;244 }245 246 424 // -------------------------------------------------------------------------- 425 // 426 // Calls: 427 // - MHCalibrationCam::FitHiGainArrays() with Bad Pixels flags 0 428 // - MHCalibrationCam::FitLoGainArrays() with Bad Pixels flags 0 429 // 430 Bool_t MHPedestalCam::FinalizeHists() 431 { 432 433 FitHiGainArrays((MCalibrationCam&)(*fCam),*fBadPixels, 434 MBadPixelsPix::kHiGainNotFitted, 435 MBadPixelsPix::kHiGainOscillating); 436 FitLoGainArrays((MCalibrationCam&)(*fCam),*fBadPixels, 437 MBadPixelsPix::kLoGainNotFitted, 438 MBadPixelsPix::kLoGainOscillating); 439 440 return kTRUE; 441 442 } 443 444 // ------------------------------------------------------------------ 247 445 // 248 446 // The types are as follows: … … 277 475 { 278 476 279 if (f Array->GetSize() <= idx)477 if (fHiGainArray->GetSize() <= idx) 280 478 return kFALSE; 281 479 … … 344 542 } 345 543 544 // -------------------------------------------------------------------------- 545 // 546 // Calls MHGausEvents::DrawClone() for pixel idx 547 // 346 548 void MHPedestalCam::DrawPixelContent(Int_t idx) const 347 549 { -
trunk/MagicSoft/Mars/manalysis/MHPedestalCam.h
r3642 r3643 32 32 void DrawPixelContent(Int_t idx) const; 33 33 34 ClassDef(MHPedestalCam, 1) // Storage Container for all pedestal information of the camera34 ClassDef(MHPedestalCam, 1) // Histogram class for Charge Camera Pedestals 35 35 }; 36 36 -
trunk/MagicSoft/Mars/manalysis/MHPedestalPix.h
r3642 r3643 11 11 private: 12 12 13 static const Int_t fgChargeNbins; // Default for fN Bins (now set to: 450 )13 static const Int_t fgChargeNbins; // Default for fNbins (now set to: 450 ) 14 14 static const Axis_t fgChargeFirst; // Default for fFirst (now set to: -0.5 ) 15 15 static const Axis_t fgChargeLast; // Default for fLast (now set to: 449.5) … … 31 31 void Renorm(); 32 32 33 ClassDef(MHPedestalPix, 1) // Histogram s for each calibrated pixel33 ClassDef(MHPedestalPix, 1) // Histogram class for Charge Pedestal Pixel 34 34 }; 35 35
Note:
See TracChangeset
for help on using the changeset viewer.