| 1 | /* ======================================================================== *\ | 
|---|
| 2 | ! | 
|---|
| 3 | ! * | 
|---|
| 4 | ! * This file is part of MARS, the MAGIC Analysis and Reconstruction | 
|---|
| 5 | ! * Software. It is distributed to you in the hope that it can be a useful | 
|---|
| 6 | ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. | 
|---|
| 7 | ! * It is distributed WITHOUT ANY WARRANTY. | 
|---|
| 8 | ! * | 
|---|
| 9 | ! * Permission to use, copy, modify and distribute this software and its | 
|---|
| 10 | ! * documentation for any purpose is hereby granted without fee, | 
|---|
| 11 | ! * provided that the above copyright notice appear in all copies and | 
|---|
| 12 | ! * that both that copyright notice and this permission notice appear | 
|---|
| 13 | ! * in supporting documentation. It is provided "as is" without express | 
|---|
| 14 | ! * or implied warranty. | 
|---|
| 15 | ! * | 
|---|
| 16 | ! | 
|---|
| 17 | ! | 
|---|
| 18 | !   Author(s): Markus Gaug   02/2004 <mailto:markus@ifae.es> | 
|---|
| 19 | ! | 
|---|
| 20 | !   Copyright: MAGIC Software Development, 2000-2004 | 
|---|
| 21 | ! | 
|---|
| 22 | ! | 
|---|
| 23 | \* ======================================================================== */ | 
|---|
| 24 | ///////////////////////////////////////////////////////////////////////////// | 
|---|
| 25 | // | 
|---|
| 26 | // MHCalibrationPulseTimeCam | 
|---|
| 27 | // | 
|---|
| 28 | // Fills the extracted signals of MExtractedSignalCam into the MHCalibrationPix-classes | 
|---|
| 29 | // MHCalibrationPulseTimeHiGainPix and MHCalibrationPulseTimeLoGainPix for every: | 
|---|
| 30 | // | 
|---|
| 31 | // - Pixel, stored in the TOrdCollection's MHCalibrationCam::fHiGainArray and | 
|---|
| 32 | //   MHCalibrationCam::fLoGainArray | 
|---|
| 33 | // | 
|---|
| 34 | // - Average pixel per AREA index (e.g. inner and outer for the MAGIC camera), | 
|---|
| 35 | //   stored in the TOrdCollection's MHCalibrationCam::fAverageHiGainAreas and | 
|---|
| 36 | //   MHCalibrationCam::fAverageLoGainAreas | 
|---|
| 37 | // | 
|---|
| 38 | // - Average pixel per camera SECTOR (e.g. sectors 1-6 for the MAGIC camera), | 
|---|
| 39 | //   stored in the TOrdCollection's MHCalibrationCam::fAverageHiGainSectors and | 
|---|
| 40 | //   MHCalibrationCam::fAverageLoGainSectors | 
|---|
| 41 | // | 
|---|
| 42 | // Every signal is taken from MExtractedSignalCam and filled into a histogram and | 
|---|
| 43 | // an array, in order to perform a Fourier analysis (see MHGausEvents). | 
|---|
| 44 | // The signals are moreover averaged on an event-by-event basis and written into | 
|---|
| 45 | // the corresponding average pixels. | 
|---|
| 46 | // | 
|---|
| 47 | // Additionally, the (FADC slice) position of the maximum is stored in an Absolute | 
|---|
| 48 | // Arrival Time histogram. This histogram serves for a rough cross-check if the | 
|---|
| 49 | // signal does not lie at or outside the edges of the extraction window. | 
|---|
| 50 | // | 
|---|
| 51 | // The PulseTime histograms are fitted to a Gaussian, mean and sigma with its errors | 
|---|
| 52 | // and the fit probability are extracted. If none of these values are NaN's and | 
|---|
| 53 | // if the probability is bigger than MHGausEvents::fProbLimit (default: 0.5%), | 
|---|
| 54 | // the fit is declared valid. | 
|---|
| 55 | // Otherwise, the fit is repeated within ranges of the previous mean | 
|---|
| 56 | // +- MHCalibrationPix::fPickupLimit (default: 5) sigma (see MHCalibrationPix::RepeatFit()) | 
|---|
| 57 | // In case this does not make the fit valid, the histogram means and RMS's are | 
|---|
| 58 | // taken directly (see MHCalibrationPix::BypassFit()) and the following flags are set: | 
|---|
| 59 | // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainNotFitted ) or | 
|---|
| 60 | // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainNotFitted ) and | 
|---|
| 61 | // - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnreliableRun   ) | 
|---|
| 62 | // | 
|---|
| 63 | // Outliers of more than MHCalibrationPix::fPickupLimit (default: 5) sigmas | 
|---|
| 64 | // from the mean are counted as Pickup events (stored in MHCalibrationPix::fPickup) | 
|---|
| 65 | // | 
|---|
| 66 | // Unless more than fNumHiGainSaturationLimit (default: 1%) of the overall FADC | 
|---|
| 67 | // slices show saturation, the following flag is set: | 
|---|
| 68 | // - MCalibrationPulseTimePix::SetHiGainSaturation(); | 
|---|
| 69 | // In that case, the calibration constants are derived from the low-gain results. | 
|---|
| 70 | // | 
|---|
| 71 | // If more than fNumLoGainSaturationLimit (default: 1%) of the overall | 
|---|
| 72 | // low-gain FADC slices saturate, the following flags are set: | 
|---|
| 73 | // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainSaturation ) and | 
|---|
| 74 | // - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnsuitableRun    ) | 
|---|
| 75 | // | 
|---|
| 76 | // The class also fills arrays with the signal vs. event number, creates a fourier | 
|---|
| 77 | // spectrum and investigates if the projected fourier components follow an exponential | 
|---|
| 78 | // distribution. In case that the probability of the exponential fit is less than | 
|---|
| 79 | // MHGausEvents::fProbLimit (default: 0.5%), the following flags are set: | 
|---|
| 80 | // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainOscillating ) or | 
|---|
| 81 | // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainOscillating ) and | 
|---|
| 82 | // - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnreliableRun     ) | 
|---|
| 83 | // | 
|---|
| 84 | // This same procedure is performed for the average pixels. | 
|---|
| 85 | // | 
|---|
| 86 | // The following results are written into MCalibrationPulseTimeCam: | 
|---|
| 87 | // | 
|---|
| 88 | // - MCalibrationPix::SetHiGainSaturation() | 
|---|
| 89 | // - MCalibrationPix::SetHiGainMean() | 
|---|
| 90 | // - MCalibrationPix::SetHiGainMeanErr() | 
|---|
| 91 | // - MCalibrationPix::SetHiGainSigma() | 
|---|
| 92 | // - MCalibrationPix::SetHiGainSigmaErr() | 
|---|
| 93 | // - MCalibrationPix::SetHiGainProb() | 
|---|
| 94 | // - MCalibrationPix::SetHiGainNumPickup() | 
|---|
| 95 | // | 
|---|
| 96 | // For all averaged areas, the fitted sigma is multiplied with the square root of | 
|---|
| 97 | // the number involved pixels in order to be able to compare it to the average of | 
|---|
| 98 | // sigmas in the camera. | 
|---|
| 99 | // | 
|---|
| 100 | ///////////////////////////////////////////////////////////////////////////// | 
|---|
| 101 | #include "MHCalibrationPulseTimeCam.h" | 
|---|
| 102 | #include "MHCalibrationCam.h" | 
|---|
| 103 |  | 
|---|
| 104 | #include "MLog.h" | 
|---|
| 105 | #include "MLogManip.h" | 
|---|
| 106 |  | 
|---|
| 107 | #include "MParList.h" | 
|---|
| 108 |  | 
|---|
| 109 | #include "MHCalibrationPix.h" | 
|---|
| 110 |  | 
|---|
| 111 | #include "MCalibrationIntensityCam.h" | 
|---|
| 112 | #include "MCalibrationChargeCam.h" | 
|---|
| 113 | #include "MCalibrationChargePix.h" | 
|---|
| 114 |  | 
|---|
| 115 | #include "MGeomCam.h" | 
|---|
| 116 | #include "MGeomPix.h" | 
|---|
| 117 |  | 
|---|
| 118 | #include "MBadPixelsIntensityCam.h" | 
|---|
| 119 | #include "MBadPixelsCam.h" | 
|---|
| 120 | #include "MBadPixelsPix.h" | 
|---|
| 121 |  | 
|---|
| 122 | #include "MRawEvtData.h" | 
|---|
| 123 | #include "MRawRunHeader.h" | 
|---|
| 124 | #include "MRawEvtPixelIter.h" | 
|---|
| 125 |  | 
|---|
| 126 | #include "MExtractedSignalCam.h" | 
|---|
| 127 | #include "MExtractedSignalPix.h" | 
|---|
| 128 |  | 
|---|
| 129 | #include "MArrayI.h" | 
|---|
| 130 | #include "MArrayD.h" | 
|---|
| 131 |  | 
|---|
| 132 | #include <TOrdCollection.h> | 
|---|
| 133 | #include <TPad.h> | 
|---|
| 134 | #include <TVirtualPad.h> | 
|---|
| 135 | #include <TCanvas.h> | 
|---|
| 136 | #include <TStyle.h> | 
|---|
| 137 | #include <TF1.h> | 
|---|
| 138 | #include <TLatex.h> | 
|---|
| 139 | #include <TLegend.h> | 
|---|
| 140 | #include <TGraph.h> | 
|---|
| 141 | #include <TEnv.h> | 
|---|
| 142 |  | 
|---|
| 143 | ClassImp(MHCalibrationPulseTimeCam); | 
|---|
| 144 |  | 
|---|
| 145 | using namespace std; | 
|---|
| 146 |  | 
|---|
| 147 | const Byte_t  MHCalibrationPulseTimeCam::fgSaturationLimit   = 254; | 
|---|
| 148 | const Byte_t  MHCalibrationPulseTimeCam::fgLowerSignalLimit  = 100; | 
|---|
| 149 | const Int_t   MHCalibrationPulseTimeCam::fgNumPixelsRequired = 2; | 
|---|
| 150 | const Int_t   MHCalibrationPulseTimeCam::fgHiGainNbins =  40; | 
|---|
| 151 | const Axis_t  MHCalibrationPulseTimeCam::fgHiGainFirst =  -0.5; | 
|---|
| 152 | const Axis_t  MHCalibrationPulseTimeCam::fgHiGainLast  =  19.5; | 
|---|
| 153 | const Float_t MHCalibrationPulseTimeCam::fgProbLimit   = 0.001; | 
|---|
| 154 | const TString MHCalibrationPulseTimeCam::gsHistName    = "PulseTime"; | 
|---|
| 155 | const TString MHCalibrationPulseTimeCam::gsHistTitle   = "Extracted Times"; | 
|---|
| 156 | const TString MHCalibrationPulseTimeCam::gsHistXTitle  = "Time [FADC slices]"; | 
|---|
| 157 | const TString MHCalibrationPulseTimeCam::gsHistYTitle  = "Nr. events"; | 
|---|
| 158 | const TString MHCalibrationPulseTimeCam::fgReferenceFile = "mjobs/signalref.rc"; | 
|---|
| 159 | // -------------------------------------------------------------------------- | 
|---|
| 160 | // | 
|---|
| 161 | // Default Constructor. | 
|---|
| 162 | // | 
|---|
| 163 | // Sets: | 
|---|
| 164 | // - all pointers to NULL | 
|---|
| 165 | // | 
|---|
| 166 | // - fNbins to fgHiGainNbins | 
|---|
| 167 | // - fFirst to fgHiGainFirst | 
|---|
| 168 | // - fLast  to fgHiGainLast | 
|---|
| 169 | // | 
|---|
| 170 | // - fHistName   to gsHistName | 
|---|
| 171 | // - fHistTitle  to gsHistTitle | 
|---|
| 172 | // - fHistXTitle to gsHistXTitle | 
|---|
| 173 | // - fHistYTitle to gsHistYTitle | 
|---|
| 174 | // | 
|---|
| 175 | // - fSaturationLimit to fgSaturationLimit | 
|---|
| 176 | // - fLowerSignalLimit to fgLowerSignalLimit | 
|---|
| 177 | // - fNumPixelsRequired to fgNumPixelsRequired | 
|---|
| 178 | // | 
|---|
| 179 | MHCalibrationPulseTimeCam::MHCalibrationPulseTimeCam(const char *name, const char *title) | 
|---|
| 180 | : fBadPixels(NULL) | 
|---|
| 181 | { | 
|---|
| 182 |  | 
|---|
| 183 | fName  = name  ? name  : "MHCalibrationPulseTimeCam"; | 
|---|
| 184 | fTitle = title ? title : "Class to fill the extracted pulse times for cosmics "; | 
|---|
| 185 |  | 
|---|
| 186 | SetNbins(fgHiGainNbins); | 
|---|
| 187 | SetFirst(fgHiGainFirst); | 
|---|
| 188 | SetLast (fgHiGainLast ); | 
|---|
| 189 |  | 
|---|
| 190 | SetProbLimit(fgProbLimit); | 
|---|
| 191 |  | 
|---|
| 192 | SetHistName  (gsHistName  .Data()); | 
|---|
| 193 | SetHistTitle (gsHistTitle .Data()); | 
|---|
| 194 | SetHistXTitle(gsHistXTitle.Data()); | 
|---|
| 195 | SetHistYTitle(gsHistYTitle.Data()); | 
|---|
| 196 |  | 
|---|
| 197 | SetReferenceFile(); | 
|---|
| 198 | SetLoGain(kFALSE); | 
|---|
| 199 | SetOscillations(kFALSE); | 
|---|
| 200 |  | 
|---|
| 201 | SetSaturationLimit(); | 
|---|
| 202 | SetLowerSignalLimit(); | 
|---|
| 203 | SetNumPixelsRequired(); | 
|---|
| 204 |  | 
|---|
| 205 | fInnerRefTime = 5.; | 
|---|
| 206 | fOuterRefTime = 5.; | 
|---|
| 207 | } | 
|---|
| 208 | // -------------------------------------------------------------------------- | 
|---|
| 209 | // | 
|---|
| 210 | // Creates new MHCalibrationPulseTimeCam only with the averaged areas: | 
|---|
| 211 | // the rest has to be retrieved directly, e.g. via: | 
|---|
| 212 | //     MHCalibrationPulseTimeCam *cam = MParList::FindObject("MHCalibrationPulseTimeCam"); | 
|---|
| 213 | //  -  cam->GetAverageSector(5).DrawClone(); | 
|---|
| 214 | //  -  (*cam)[100].DrawClone() | 
|---|
| 215 | // | 
|---|
| 216 | TObject *MHCalibrationPulseTimeCam::Clone(const char *) const | 
|---|
| 217 | { | 
|---|
| 218 |  | 
|---|
| 219 | MHCalibrationPulseTimeCam *cam = new MHCalibrationPulseTimeCam(); | 
|---|
| 220 |  | 
|---|
| 221 | // | 
|---|
| 222 | // Copy the data members | 
|---|
| 223 | // | 
|---|
| 224 | cam->fColor                  = fColor; | 
|---|
| 225 | cam->fRunNumbers             = fRunNumbers; | 
|---|
| 226 | cam->fPulserFrequency        = fPulserFrequency; | 
|---|
| 227 | cam->fFlags                  = fFlags; | 
|---|
| 228 | cam->fNbins                  = fNbins; | 
|---|
| 229 | cam->fFirst                  = fFirst; | 
|---|
| 230 | cam->fLast                   = fLast; | 
|---|
| 231 | cam->fReferenceFile          = fReferenceFile; | 
|---|
| 232 |  | 
|---|
| 233 | if (!IsAverageing()) | 
|---|
| 234 | return cam; | 
|---|
| 235 |  | 
|---|
| 236 | const Int_t navhi   =  fAverageHiGainAreas->GetSize(); | 
|---|
| 237 |  | 
|---|
| 238 | for (int i=0; i<navhi; i++) | 
|---|
| 239 | cam->fAverageHiGainAreas->AddAt(GetAverageHiGainArea(i).Clone(),i); | 
|---|
| 240 |  | 
|---|
| 241 | return cam; | 
|---|
| 242 | } | 
|---|
| 243 |  | 
|---|
| 244 | // -------------------------------------------------------------------------- | 
|---|
| 245 | // | 
|---|
| 246 | // Gets the pointers to: | 
|---|
| 247 | // - MRawEvtData | 
|---|
| 248 | // | 
|---|
| 249 | Bool_t MHCalibrationPulseTimeCam::SetupHists(const MParList *pList) | 
|---|
| 250 | { | 
|---|
| 251 |  | 
|---|
| 252 | fBadPixels = (MBadPixelsCam*)pList->FindObject("MBadPixelsCam"); | 
|---|
| 253 | if (!fBadPixels) | 
|---|
| 254 | { | 
|---|
| 255 | *fLog << warn << GetDescriptor() << "MBadPixelsCam not found... " << endl; | 
|---|
| 256 | } | 
|---|
| 257 |  | 
|---|
| 258 | return kTRUE; | 
|---|
| 259 | } | 
|---|
| 260 |  | 
|---|
| 261 | // -------------------------------------------------------------------------- | 
|---|
| 262 | // | 
|---|
| 263 | // Gets or creates the pointers to: | 
|---|
| 264 | // - MExtractedSignalCam | 
|---|
| 265 | // - MCalibrationPulseTimeCam or MCalibrationIntensityPulseTimeCam | 
|---|
| 266 | // - MBadPixelsCam | 
|---|
| 267 | // | 
|---|
| 268 | // Initializes the number of used FADC slices from MExtractedSignalCam | 
|---|
| 269 | // into MCalibrationPulseTimeCam and test for changes in that variable | 
|---|
| 270 | // | 
|---|
| 271 | // Calls: | 
|---|
| 272 | // - InitHiGainArrays() | 
|---|
| 273 | // | 
|---|
| 274 | // Sets: | 
|---|
| 275 | // - fSumhiarea   to nareas | 
|---|
| 276 | // - fSumloarea   to nareas | 
|---|
| 277 | // - fSumhisector to nsectors | 
|---|
| 278 | // - fSumlosector to nsectors | 
|---|
| 279 | // | 
|---|
| 280 | Bool_t MHCalibrationPulseTimeCam::ReInitHists(MParList *pList) | 
|---|
| 281 | { | 
|---|
| 282 |  | 
|---|
| 283 | MExtractedSignalCam *signal = | 
|---|
| 284 | (MExtractedSignalCam*)pList->FindObject(AddSerialNumber("MExtractedSignalCam")); | 
|---|
| 285 | if (!signal) | 
|---|
| 286 | { | 
|---|
| 287 | *fLog << err << "MExtractedSignalCam not found... abort." << endl; | 
|---|
| 288 | return kFALSE; | 
|---|
| 289 | } | 
|---|
| 290 |  | 
|---|
| 291 | if (!InitCams(pList,"PulseTime")) | 
|---|
| 292 | return kFALSE; | 
|---|
| 293 |  | 
|---|
| 294 | const Int_t npixels  = fGeom->GetNumPixels(); | 
|---|
| 295 | const Int_t nsectors = fGeom->GetNumSectors(); | 
|---|
| 296 | const Int_t nareas   = fGeom->GetNumAreas(); | 
|---|
| 297 |  | 
|---|
| 298 | InitHiGainArrays(npixels,nareas,nsectors); | 
|---|
| 299 |  | 
|---|
| 300 | fSumhiarea  .Set(nareas); | 
|---|
| 301 | fSumhisector.Set(nsectors); | 
|---|
| 302 |  | 
|---|
| 303 | return kTRUE; | 
|---|
| 304 | } | 
|---|
| 305 |  | 
|---|
| 306 | void MHCalibrationPulseTimeCam::InitHiGainArrays(const Int_t npixels, const Int_t nareas, const Int_t nsectors) | 
|---|
| 307 | { | 
|---|
| 308 |  | 
|---|
| 309 | if (fHiGainArray->GetSize()==0) | 
|---|
| 310 | { | 
|---|
| 311 | for (Int_t i=0; i<npixels; i++) | 
|---|
| 312 | { | 
|---|
| 313 | fHiGainArray->AddAt(new MHCalibrationPix(Form("%sHiGainPix%04d",fHistName.Data(),i), | 
|---|
| 314 | Form("%s High Gain Pixel %4d",fHistTitle.Data(),i)),i); | 
|---|
| 315 |  | 
|---|
| 316 | MHCalibrationPix &pix = (*this)[i]; | 
|---|
| 317 | pix.SetNbins(fNbins*2); | 
|---|
| 318 | pix.SetFirst(fFirst); | 
|---|
| 319 | pix.SetLast (fLast); | 
|---|
| 320 |  | 
|---|
| 321 | MBadPixelsPix &bad = fIntensBad ? (*fIntensBad)[i] : (*fBadPixels)[i]; | 
|---|
| 322 | InitHists(pix,bad,i); | 
|---|
| 323 |  | 
|---|
| 324 | if (fCam) | 
|---|
| 325 | (*fCam)[i].SetPixId(i); | 
|---|
| 326 | } | 
|---|
| 327 | } | 
|---|
| 328 |  | 
|---|
| 329 | if (!IsAverageing()) | 
|---|
| 330 | return; | 
|---|
| 331 |  | 
|---|
| 332 | if (fAverageHiGainAreas->GetSize()==0) | 
|---|
| 333 | { | 
|---|
| 334 | for (Int_t j=0; j<nareas; j++) | 
|---|
| 335 | { | 
|---|
| 336 | fAverageHiGainAreas->AddAt(new MHCalibrationPix(Form("%sHiGainArea%d",fHistName.Data(),j), | 
|---|
| 337 | Form("%s High Gain Area Idx %d",fHistTitle.Data(),j)),j); | 
|---|
| 338 |  | 
|---|
| 339 | MHCalibrationPix &pix = GetAverageHiGainArea(j); | 
|---|
| 340 |  | 
|---|
| 341 | pix.SetNbins(fNbins*2); | 
|---|
| 342 | pix.SetFirst(fFirst); | 
|---|
| 343 | pix.SetLast (fLast); | 
|---|
| 344 |  | 
|---|
| 345 | InitHists(pix,fIntensCam ? fIntensCam->GetAverageBadArea(j) : fCam->GetAverageBadArea(j),j); | 
|---|
| 346 | } | 
|---|
| 347 | } | 
|---|
| 348 |  | 
|---|
| 349 | if (fAverageHiGainSectors->GetSize()==0) | 
|---|
| 350 | { | 
|---|
| 351 | for (Int_t j=0; j<nsectors; j++) | 
|---|
| 352 | { | 
|---|
| 353 | fAverageHiGainSectors->AddAt(new MHCalibrationPix(Form("%sHiGainSector%02d",fHistName.Data(),j), | 
|---|
| 354 | Form("%s High Gain Sector %02d",fHistTitle.Data(),j)),j); | 
|---|
| 355 | MHCalibrationPix &pix = GetAverageHiGainSector(j); | 
|---|
| 356 |  | 
|---|
| 357 | pix.SetNbins(fNbins); | 
|---|
| 358 | pix.SetFirst(fFirst); | 
|---|
| 359 | pix.SetLast (fLast); | 
|---|
| 360 |  | 
|---|
| 361 | InitHists(pix,fIntensCam ? fIntensCam->GetAverageBadSector(j) : fCam->GetAverageBadSector(j),j); | 
|---|
| 362 | } | 
|---|
| 363 | } | 
|---|
| 364 | } | 
|---|
| 365 |  | 
|---|
| 366 | // -------------------------------------------------------------------------- | 
|---|
| 367 | // | 
|---|
| 368 | // Retrieves from MExtractedSignalCam: | 
|---|
| 369 | // - first used LoGain FADC slice | 
|---|
| 370 | // | 
|---|
| 371 | // Retrieves from MGeomCam: | 
|---|
| 372 | // - number of pixels | 
|---|
| 373 | // - number of pixel areas | 
|---|
| 374 | // - number of sectors | 
|---|
| 375 | // | 
|---|
| 376 | // For all TOrdCollection's (including the averaged ones), the following steps are performed: | 
|---|
| 377 | // | 
|---|
| 378 | // 1) Fill PulseTimes histograms (MHGausEvents::FillHistAndArray()) with: | 
|---|
| 379 | // - MExtractedSignalPix::GetExtractedSignalHiGain(); | 
|---|
| 380 | // - MExtractedSignalPix::GetExtractedSignalLoGain(); | 
|---|
| 381 | // | 
|---|
| 382 | Bool_t MHCalibrationPulseTimeCam::FillHists(const MParContainer *par, const Stat_t w) | 
|---|
| 383 | { | 
|---|
| 384 |  | 
|---|
| 385 | MRawEvtData *data = (MRawEvtData*)par; | 
|---|
| 386 | if (!data) | 
|---|
| 387 | { | 
|---|
| 388 | *fLog << err << "No argument in MHCalibrationPulseTimeCam::Fill... abort." << endl; | 
|---|
| 389 | return kFALSE; | 
|---|
| 390 | } | 
|---|
| 391 |  | 
|---|
| 392 | const UInt_t nareas   = fGeom->GetNumAreas(); | 
|---|
| 393 | const UInt_t nsectors = fGeom->GetNumSectors(); | 
|---|
| 394 |  | 
|---|
| 395 | fSumhiarea  .Reset(); | 
|---|
| 396 | fSumhisector.Reset(); | 
|---|
| 397 | fAverageAreaNum.Reset(); | 
|---|
| 398 | fAverageSectorNum.Reset(); | 
|---|
| 399 |  | 
|---|
| 400 | MRawEvtPixelIter pixel(data); | 
|---|
| 401 | while (pixel.Next()) | 
|---|
| 402 | { | 
|---|
| 403 |  | 
|---|
| 404 | const Int_t i = pixel.GetPixelId(); | 
|---|
| 405 |  | 
|---|
| 406 | if (fBadPixels) | 
|---|
| 407 | { | 
|---|
| 408 | MBadPixelsPix &bad = (*fBadPixels)[i]; | 
|---|
| 409 | if (bad.IsUnsuitable()) | 
|---|
| 410 | continue; | 
|---|
| 411 | } | 
|---|
| 412 |  | 
|---|
| 413 | Byte_t *start = pixel.GetHiGainSamples(); | 
|---|
| 414 | Byte_t *end   = start + pixel.GetNumHiGainSamples(); | 
|---|
| 415 | Byte_t *p     = start; | 
|---|
| 416 | Byte_t max    = 0; | 
|---|
| 417 | Int_t  maxpos = 0; | 
|---|
| 418 |  | 
|---|
| 419 | while (p < end) | 
|---|
| 420 | { | 
|---|
| 421 | if ((*p > max) && (*p < fSaturationLimit)) | 
|---|
| 422 | { | 
|---|
| 423 | max = *p; | 
|---|
| 424 | maxpos = p-start-1; | 
|---|
| 425 | } | 
|---|
| 426 | p++; | 
|---|
| 427 | } | 
|---|
| 428 |  | 
|---|
| 429 | start   = pixel.GetLoGainSamples(); | 
|---|
| 430 | end     = start + pixel.GetNumLoGainSamples(); | 
|---|
| 431 | p       = start; | 
|---|
| 432 |  | 
|---|
| 433 | while (p < end) | 
|---|
| 434 | { | 
|---|
| 435 | if ((*p > max) && (*p < fSaturationLimit)) | 
|---|
| 436 | { | 
|---|
| 437 | max = *p; | 
|---|
| 438 | maxpos = p-start+pixel.GetNumHiGainSamples() - 1; | 
|---|
| 439 | } | 
|---|
| 440 | p++; | 
|---|
| 441 | } | 
|---|
| 442 |  | 
|---|
| 443 | if (max < fLowerSignalLimit) | 
|---|
| 444 | continue; | 
|---|
| 445 |  | 
|---|
| 446 | const Float_t time = (Float_t)maxpos; | 
|---|
| 447 |  | 
|---|
| 448 | (*this)[i].FillHist(time); | 
|---|
| 449 |  | 
|---|
| 450 | const Int_t aidx   = (*fGeom)[i].GetAidx(); | 
|---|
| 451 | const Int_t sector = (*fGeom)[i].GetSector(); | 
|---|
| 452 |  | 
|---|
| 453 | fSumhiarea[aidx]      += time; | 
|---|
| 454 | fSumhisector[sector]  += time; | 
|---|
| 455 |  | 
|---|
| 456 | fAverageAreaNum[aidx]++; | 
|---|
| 457 | fAverageSectorNum[sector]++; | 
|---|
| 458 | } | 
|---|
| 459 |  | 
|---|
| 460 | for (UInt_t j=0; j<nareas; j++) | 
|---|
| 461 | { | 
|---|
| 462 |  | 
|---|
| 463 | const Int_t npix = fAverageAreaNum[j]; | 
|---|
| 464 |  | 
|---|
| 465 | if (npix > fNumPixelsRequired) | 
|---|
| 466 | { | 
|---|
| 467 | if (IsOscillations()) | 
|---|
| 468 | GetAverageHiGainArea(j).FillHistAndArray(fSumhiarea[j]/npix); | 
|---|
| 469 | else | 
|---|
| 470 | GetAverageHiGainArea(j).FillHist(fSumhiarea[j]/npix); | 
|---|
| 471 |  | 
|---|
| 472 | } | 
|---|
| 473 | } | 
|---|
| 474 |  | 
|---|
| 475 | for (UInt_t j=0; j<nsectors; j++) | 
|---|
| 476 | { | 
|---|
| 477 |  | 
|---|
| 478 | const Int_t npix = fAverageSectorNum[j]; | 
|---|
| 479 |  | 
|---|
| 480 | if (npix > 0) | 
|---|
| 481 | { | 
|---|
| 482 | if (IsOscillations()) | 
|---|
| 483 | GetAverageHiGainSector(j).FillHistAndArray(fSumhisector [j]/npix); | 
|---|
| 484 | else | 
|---|
| 485 | GetAverageHiGainSector(j).FillHist(fSumhisector [j]/npix); | 
|---|
| 486 | } | 
|---|
| 487 | } | 
|---|
| 488 |  | 
|---|
| 489 | return kTRUE; | 
|---|
| 490 | } | 
|---|
| 491 |  | 
|---|
| 492 | // -------------------------------------------------------------------------- | 
|---|
| 493 | // | 
|---|
| 494 | // For all TOrdCollection's (including the averaged ones), the following steps are performed: | 
|---|
| 495 | // | 
|---|
| 496 | // 1) Returns if the pixel is excluded. | 
|---|
| 497 | // 2) Tests saturation. In case yes, set the flag: MCalibrationPix::SetHiGainSaturation() | 
|---|
| 498 | //    or the flag: MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainSaturated ) | 
|---|
| 499 | // 3) Store the absolute arrival times in the MCalibrationPulseTimePix's. If flag | 
|---|
| 500 | //    MCalibrationPix::IsHiGainSaturation() is set, the Low-Gain arrival times are stored, | 
|---|
| 501 | //    otherwise the Hi-Gain ones. | 
|---|
| 502 | // 4) Calls to MHCalibrationCam::FitHiGainArrays() and MCalibrationCam::FitLoGainArrays() | 
|---|
| 503 | //    with the flags: | 
|---|
| 504 | //    - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainNotFitted ) | 
|---|
| 505 | //    - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainNotFitted ) | 
|---|
| 506 | //    - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainOscillating ) | 
|---|
| 507 | //    - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainOscillating ) | 
|---|
| 508 | // | 
|---|
| 509 | Bool_t MHCalibrationPulseTimeCam::FinalizeHists() | 
|---|
| 510 | { | 
|---|
| 511 |  | 
|---|
| 512 | *fLog << endl; | 
|---|
| 513 |  | 
|---|
| 514 | MCalibrationCam *calcam  = fIntensCam ? fIntensCam->GetCam() : fCam; | 
|---|
| 515 | // | 
|---|
| 516 | // Perform the fitting for the High Gain (done in MHCalibrationCam) | 
|---|
| 517 | // | 
|---|
| 518 | for (Int_t i=0; i<fHiGainArray->GetSize(); i++) | 
|---|
| 519 | { | 
|---|
| 520 |  | 
|---|
| 521 | MHCalibrationPix &hist = (*this)[i]; | 
|---|
| 522 |  | 
|---|
| 523 | if (hist.IsExcluded()) | 
|---|
| 524 | continue; | 
|---|
| 525 |  | 
|---|
| 526 | MCalibrationPix &pix    = (*calcam)[i]; | 
|---|
| 527 | CalcHists(hist,pix); | 
|---|
| 528 | } | 
|---|
| 529 |  | 
|---|
| 530 | if (!IsAverageing()) | 
|---|
| 531 | return kTRUE; | 
|---|
| 532 |  | 
|---|
| 533 | for (Int_t j=0; j<fAverageHiGainAreas->GetSize(); j++) | 
|---|
| 534 | { | 
|---|
| 535 |  | 
|---|
| 536 | MHCalibrationPix &hist = GetAverageHiGainArea(j); | 
|---|
| 537 | MCalibrationPix  &pix  = calcam->GetAverageArea(j); | 
|---|
| 538 | CalcHists(hist,pix); | 
|---|
| 539 | } | 
|---|
| 540 |  | 
|---|
| 541 | for (Int_t j=0; j<fAverageHiGainSectors->GetSize(); j++) | 
|---|
| 542 | { | 
|---|
| 543 | MHCalibrationPix &hist = GetAverageHiGainSector(j); | 
|---|
| 544 | MCalibrationPix  &pix  = calcam->GetAverageSector(j); | 
|---|
| 545 | CalcHists(hist,pix); | 
|---|
| 546 | } | 
|---|
| 547 |  | 
|---|
| 548 | return kTRUE; | 
|---|
| 549 | } | 
|---|
| 550 |  | 
|---|
| 551 | void MHCalibrationPulseTimeCam::CalcHists(MHCalibrationPix &hist, MCalibrationPix &pix) | 
|---|
| 552 | { | 
|---|
| 553 |  | 
|---|
| 554 |  | 
|---|
| 555 | if (hist.IsEmpty() || hist.IsOnlyOverflow() || hist.IsOnlyUnderflow()) | 
|---|
| 556 | { | 
|---|
| 557 | *fLog << warn << GetDescriptor() << ": Only over- or underflow in " << hist.GetName() << endl; | 
|---|
| 558 | return; | 
|---|
| 559 | } | 
|---|
| 560 |  | 
|---|
| 561 | hist.BypassFit(); | 
|---|
| 562 |  | 
|---|
| 563 | pix.SetHiGainMean       ( hist.GetMean()      ); | 
|---|
| 564 | pix.SetHiGainMeanVar    ( hist.GetMeanErr() * hist.GetMeanErr()   ); | 
|---|
| 565 | pix.SetHiGainRms        ( hist.GetHistRms()   ); | 
|---|
| 566 | pix.SetHiGainSigma      ( hist.GetSigma()     ); | 
|---|
| 567 | pix.SetHiGainSigmaVar   ( hist.GetSigmaErr()* hist.GetSigmaErr()  ); | 
|---|
| 568 |  | 
|---|
| 569 | if (IsDebug()) | 
|---|
| 570 | { | 
|---|
| 571 | *fLog << dbginf << GetDescriptor() << ": ID " << GetName() | 
|---|
| 572 | << " "<<pix.GetPixId() | 
|---|
| 573 | << " Mean: "         << hist.GetMean    () | 
|---|
| 574 | << " MeanErr: "      << hist.GetMeanErr () | 
|---|
| 575 | << " MeanSigma: "    << hist.GetSigma   () | 
|---|
| 576 | << " MeanSigmaErr: " << hist.GetSigmaErr() | 
|---|
| 577 | << " Prob: "         << hist.GetProb    () | 
|---|
| 578 | << endl; | 
|---|
| 579 | } | 
|---|
| 580 | return; | 
|---|
| 581 | } | 
|---|
| 582 |  | 
|---|
| 583 |  | 
|---|
| 584 | // -------------------------------------------------------------------------- | 
|---|
| 585 | // | 
|---|
| 586 | // Calls MHCalibrationPix::DrawClone() for pixel idx | 
|---|
| 587 | // | 
|---|
| 588 | void MHCalibrationPulseTimeCam::DrawPixelContent(Int_t idx) const | 
|---|
| 589 | { | 
|---|
| 590 | (*this)[idx].DrawClone(); | 
|---|
| 591 | } | 
|---|
| 592 |  | 
|---|
| 593 |  | 
|---|
| 594 | // ----------------------------------------------------------------------------- | 
|---|
| 595 | // | 
|---|
| 596 | // Default draw: | 
|---|
| 597 | // | 
|---|
| 598 | // Displays the averaged areas, both High Gain and Low Gain | 
|---|
| 599 | // | 
|---|
| 600 | // Calls the Draw of the fAverageHiGainAreas and fAverageLoGainAreas objects with options | 
|---|
| 601 | // | 
|---|
| 602 | void MHCalibrationPulseTimeCam::Draw(const Option_t *opt) | 
|---|
| 603 | { | 
|---|
| 604 |  | 
|---|
| 605 | const Int_t nareas = fAverageHiGainAreas->GetSize(); | 
|---|
| 606 | if (nareas == 0) | 
|---|
| 607 | return; | 
|---|
| 608 |  | 
|---|
| 609 | TString option(opt); | 
|---|
| 610 | option.ToLower(); | 
|---|
| 611 |  | 
|---|
| 612 | TVirtualPad *pad = gPad ? gPad : MH::MakeDefCanvas(this); | 
|---|
| 613 | pad->SetBorderMode(0); | 
|---|
| 614 | pad->Divide(1,nareas); | 
|---|
| 615 |  | 
|---|
| 616 | // | 
|---|
| 617 | // Loop over inner and outer pixels | 
|---|
| 618 | // | 
|---|
| 619 | for (Int_t i=0; i<nareas;i++) | 
|---|
| 620 | { | 
|---|
| 621 |  | 
|---|
| 622 | pad->cd(i+1); | 
|---|
| 623 |  | 
|---|
| 624 | MHCalibrationPix &hipix = GetAverageHiGainArea(i); | 
|---|
| 625 | DrawDataCheckPixel(hipix,i ? fOuterRefTime : fInnerRefTime); | 
|---|
| 626 | } | 
|---|
| 627 | } | 
|---|
| 628 |  | 
|---|
| 629 | // ----------------------------------------------------------------------------- | 
|---|
| 630 | // | 
|---|
| 631 | // Draw the average pixel for the datacheck: | 
|---|
| 632 | // | 
|---|
| 633 | // Displays the averaged areas, both High Gain and Low Gain | 
|---|
| 634 | // | 
|---|
| 635 | // Calls the Draw of the fAverageHiGainAreas and fAverageLoGainAreas objects with options | 
|---|
| 636 | // | 
|---|
| 637 | void MHCalibrationPulseTimeCam::DrawDataCheckPixel(MHCalibrationPix &pix, const Float_t refline) | 
|---|
| 638 | { | 
|---|
| 639 |  | 
|---|
| 640 | gPad->SetTicks(); | 
|---|
| 641 |  | 
|---|
| 642 | TH1F *hist = pix.GetHGausHist(); | 
|---|
| 643 |  | 
|---|
| 644 | TH1F *null = new TH1F("Null",hist->GetTitle(),100, | 
|---|
| 645 | pix.GetFirst() > 0. ? pix.GetFirst() : 0., | 
|---|
| 646 | pix.GetLast() > pix.GetFirst() | 
|---|
| 647 | ? pix.GetLast() : pix.GetFirst()*2.); | 
|---|
| 648 |  | 
|---|
| 649 | null->SetMaximum(1.1*hist->GetMaximum()); | 
|---|
| 650 | null->SetDirectory(NULL); | 
|---|
| 651 | null->SetBit(kCanDelete); | 
|---|
| 652 | null->SetStats(kFALSE); | 
|---|
| 653 | // | 
|---|
| 654 | // set the labels bigger | 
|---|
| 655 | // | 
|---|
| 656 | TAxis *xaxe = null->GetXaxis(); | 
|---|
| 657 | TAxis *yaxe = null->GetYaxis(); | 
|---|
| 658 | xaxe->CenterTitle(); | 
|---|
| 659 | yaxe->CenterTitle(); | 
|---|
| 660 | xaxe->SetTitleSize(0.06); | 
|---|
| 661 | yaxe->SetTitleSize(0.076); | 
|---|
| 662 | xaxe->SetTitleOffset(0.6); | 
|---|
| 663 | yaxe->SetTitleOffset(0.65); | 
|---|
| 664 | xaxe->SetLabelSize(0.06); | 
|---|
| 665 | yaxe->SetLabelSize(0.06); | 
|---|
| 666 | xaxe->SetTitle(hist->GetXaxis()->GetTitle()); | 
|---|
| 667 | yaxe->SetTitle(hist->GetYaxis()->GetTitle()); | 
|---|
| 668 |  | 
|---|
| 669 | null->Draw(); | 
|---|
| 670 | hist->Draw("same"); | 
|---|
| 671 |  | 
|---|
| 672 | DisplayRefLines(null,refline); | 
|---|
| 673 |  | 
|---|
| 674 | return; | 
|---|
| 675 |  | 
|---|
| 676 | } | 
|---|
| 677 |  | 
|---|
| 678 | void  MHCalibrationPulseTimeCam::DisplayRefLines(const TH1F *hist, const Float_t refline) const | 
|---|
| 679 | { | 
|---|
| 680 |  | 
|---|
| 681 | TGraph *gr    = new TGraph(2); | 
|---|
| 682 | gr->SetPoint(0,refline,0.); | 
|---|
| 683 | gr->SetPoint(1,refline,hist->GetMaximum()); | 
|---|
| 684 | gr->SetBit(kCanDelete); | 
|---|
| 685 | gr->SetLineColor(106); | 
|---|
| 686 | gr->SetLineStyle(2); | 
|---|
| 687 | gr->SetLineWidth(3); | 
|---|
| 688 | gr->Draw("L"); | 
|---|
| 689 |  | 
|---|
| 690 | TLegend *leg = new TLegend(0.8,0.35,0.99,0.65); | 
|---|
| 691 | leg->SetBit(kCanDelete); | 
|---|
| 692 | leg->AddEntry(gr,"Reference","l"); | 
|---|
| 693 |  | 
|---|
| 694 | leg->Draw(); | 
|---|
| 695 | } | 
|---|
| 696 |  | 
|---|
| 697 | Int_t MHCalibrationPulseTimeCam::ReadEnv(const TEnv &env, TString prefix, Bool_t print) | 
|---|
| 698 | { | 
|---|
| 699 |  | 
|---|
| 700 | Bool_t rc = kFALSE; | 
|---|
| 701 |  | 
|---|
| 702 | if (MHCalibrationCam::ReadEnv(env,prefix,print)) | 
|---|
| 703 | rc = kTRUE; | 
|---|
| 704 |  | 
|---|
| 705 | if (IsEnvDefined(env, prefix, "SaturationLimit", print)) | 
|---|
| 706 | { | 
|---|
| 707 | SetSaturationLimit(GetEnvValue(env, prefix, "SaturationLimit", fSaturationLimit)); | 
|---|
| 708 | rc = kTRUE; | 
|---|
| 709 | } | 
|---|
| 710 |  | 
|---|
| 711 | if (IsEnvDefined(env, prefix, "LowerSignalLimit", print)) | 
|---|
| 712 | { | 
|---|
| 713 | SetLowerSignalLimit(GetEnvValue(env,prefix,"LowerSignalLimit",fLowerSignalLimit)); | 
|---|
| 714 | rc = kTRUE; | 
|---|
| 715 | } | 
|---|
| 716 |  | 
|---|
| 717 | if (IsEnvDefined(env, prefix, "NumPixelsRequired", print)) | 
|---|
| 718 | { | 
|---|
| 719 | SetNumPixelsRequired(GetEnvValue(env,prefix,"NumPixelsRequired",fNumPixelsRequired)); | 
|---|
| 720 | rc = kTRUE; | 
|---|
| 721 | } | 
|---|
| 722 |  | 
|---|
| 723 | if (IsEnvDefined(env, prefix, "ReferenceFile", print)) | 
|---|
| 724 | { | 
|---|
| 725 | SetReferenceFile(GetEnvValue(env,prefix,"ReferenceFile",fReferenceFile.Data())); | 
|---|
| 726 | rc = kTRUE; | 
|---|
| 727 | } | 
|---|
| 728 |  | 
|---|
| 729 | TEnv refenv(fReferenceFile); | 
|---|
| 730 |  | 
|---|
| 731 | fInnerRefTime = refenv.GetValue("InnerRefTime",fInnerRefTime); | 
|---|
| 732 | fOuterRefTime = refenv.GetValue("OuterRefTime",fOuterRefTime); | 
|---|
| 733 |  | 
|---|
| 734 | return rc; | 
|---|
| 735 | } | 
|---|