/* ======================================================================== *\ ! ! * ! * This file is part of MARS, the MAGIC Analysis and Reconstruction ! * Software. It is distributed to you in the hope that it can be a useful ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. ! * It is distributed WITHOUT ANY WARRANTY. ! * ! * Permission to use, copy, modify and distribute this software and its ! * documentation for any purpose is hereby granted without fee, ! * provided that the above copyright notice appear in all copies and ! * that both that copyright notice and this permission notice appear ! * in supporting documentation. It is provided "as is" without express ! * or implied warranty. ! * ! ! ! Author(s): Thomas Bretz, 12/2000 ! Author(s): Harald Kornmayer, 1/2001 ! Author(s): Nadia Tonello, 4/2003 ! ! Copyright: MAGIC Software Development, 2000-2003 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // MImgCleanStd // // The Image Cleaning task selects the pixels you use for the Hillas // parameters calculation. // // There are two methods to make the selection: the standard one, as done // in the analysis of CT1 data, and the democratic one, as suggested by // W.Wittek. The number of photo-electrons of a pixel is compared with the // pedestal RMS of the pixel itself (standard method) or with the average // RMS of the inner pixels (democratic method). // In both cases, the possibility to have a camera with pixels of // different area is taken into account. // The too noisy pixels can be recognized and eventally switched off // (Unmap: set blind pixels to UNUSED) separately, using the // MBlindPixelCalc Class. In the MBlindPixelCalc class there is also the // function to replace the value of the noisy pixels with the interpolation // of the content of the neighbors (SetUseInterpolation). // // Example: // ... // MBlindPixelCalc blind; // blind.SetUseInterpolation(); // blind.SetUseBlindPixels(); // // MImgCleanStd clean; // ... // tlist.AddToList(&blind); // tlist.AddToList(&clean); // // Look at the MBlindPixelCalc Class for more details. // // Starting point: default values ---------------------------------------- // // When an event is read, before the image cleaning, all the pixels that // are in MCerPhotEvt are set as USED and NOT CORE. All the pixels belong // to RING number 1 (like USED pixels). // Look at MCerPhotPix.h to see how these informations of the pixel are // stored. // The default cleaning METHOD is the STANDARD one and the number of the // rings around the CORE pixel it analyzes is 1. Look at the Constructor // of the class in MImgCleanStd.cc to see (or change) the default values. // // Example: To modify this setting, use the member functions // SetMethod(MImgCleanStd::kDemocratic) and SetCleanRings(UShort_t n). // // MImgCleanStd:CleanStep1 ----------------------------------------------- // // The first step of cleaning defines the CORE pixels. The CORE pixels are // the ones which contain the informations about the core of the electro- // magnetic shower. // The ratio (A_0/A_i) is calculated from fCam->GetPixRatio(i). A_0 is // the area of the central pixel of the camera, A_i is the area of the // examined pixel. In this way, if we have a MAGIC-like camera, with the // outer pixels bigger than the inner ones, the level of cleaning in the // two different regions is weighted. // This avoids problems of deformations of the shower images. // The signal S_i and the pedestal RMS Prms_i of the pixel are called from // the object MCerPhotPix. // If (default method = kStandard) //Begin_Html //       //End_Html // the pixel is set as CORE pixel. L_1 (n=1) is called "first level of // cleaning" (default: fCleanLvl1 = 3). // All the other pixels are set as UNUSED and belong to RING 0. // After this point, only the CORE pixels are set as USED, with RING // number 1. // // MImgCleanStd:CleanStep2 ---------------------------------------------- // // The second step of cleaning looks at the isolated CORE pixels and sets // them to UNUSED. An isolated pixel is a pixel without CORE neighbors. // At the end of this point, we have set as USED only CORE pixels with at // least one CORE neighbor. // // MImgCleanStd:CleanStep3 ---------------------------------------------- // // The third step of cleaning looks at all the pixels (USED or UNUSED) that // surround the USED pixels. // If the content of the analyzed pixel survives at the second level of // cleaning, i.e. if //Begin_Html //       //End_Html // the pixel is set as USED. L_2 (n=2) is called "second level of cleaning" // (default:fCleanLvl2 = 2.5). // // When the number of RINGS to analyze is 1 (default value), only the // pixels that have a neighbor CORE pixel are analyzed. // // There is the option to decide the number of times you want to repeat // this procedure (number of RINGS analyzed around the core pixels = n). // Every time the level of cleaning is the same (fCleanLvl2) and the pixel // will belong to ring r+1, 1 < r < n+1. This is described in // MImgCleanStd:CleanStep4 . // // Dictionary and member functions --------------------------------------- // // Here there is the detailed description of the member functions and of // the terms commonly used in the class. // // // STANDARD CLEANING: // ================= // This is the method used for the CT1 data analysis. It is the default // method of the class. // The number of photo-electrons of a pixel (S_i) is compared to the // pedestal RMS of the pixel itself (Prms_i). To have the comparison to // the same photon density for all the pixels, taking into account they // can have different areas, we have to keep in mind that the number of // photons that hit each pixel, goes linearly with the area of the pixel. // The fluctuations of the LONS are proportional to sqrt(A_i), so when we // compare S_i with Prms_i, only a factor sqrt(A_0/A_i) is missing to // have the same (N.photons/Area) threshold for all the pixels. // // !!WARNING: if noise independent from the // pixel size (example: electronic noise) is introduced, // then the noise fluctuations are no longer proportional // to sqrt(A_i), and then the cut value (for a camera with // pixels of different sizes) resulting from the above // procedure would not be proportional to pixel size as we // intend. In that case, democratic cleaning is preferred. // // If //Begin_Html //       //End_Html // the pixel survives the cleaning and it is set as CORE (when L_n is the // first level of cleaning, fCleanLvl1) or USED (when L_n is the second // level of cleaning, fCleanLvl2). // // Example: // // MImgCleanStd clean; // //creates a default Cleaning object, with default setting // ... // tlist.AddToList(&clean); // // add the image cleaning to the main task list // // // DEMOCRATIC CLEANING: // =================== // You use this cleaning method when you want to compare the number of // photo-electons of each pixel with the average pedestal RMS of the // inner pixels (for the MAGIC camera they are the smaller ones): //Begin_Html //       //End_Html // In this case, the simple ratio (A_0/A_i) is used to weight the level of // cleaning, because both the inner and the outer pixels (that in MAGIC // have a different area) are compared to the same pedestal RMS, coming // from the inner pixels. // // Make sure that you used a class calculating the MPedPhotCam which also // updated the contents of the mean values (Recalc) correctly. // // // PROBABILITY CLEANING // ==================== // This method takes signal height (over signal noise) and arrival time // into account. Instead of comparing signal/Noise with cleaning level // one and two, we calculate // - P_ped: The probability that a signal is a pedestal (using // the signal height and the pedestal) For this probability the // same algorithm like in kScaled is used (which is a standard // cleaning which scales the noise with the mean noise of pixels // with the same size) // - P_sig: The probability that the signal corresponds to the pixel // with the highest signal. For this probability we use the // arrival time only. // // The cleaning now is done in levels of Probability (eg. 0.2, 0.05) // The meaning of the cleaning levels is essentially the same (the same cleaning // algorithm is used) but the cleaning is not done in levels of signal/noise // but in level of this probability. // // This probability is calculated as (1-P_ped)*P_sig // // Example: // // MImgCleanStd clean(0.2, 0.05); // clean.SetMethod(MImgCleanStd::kProbability); // // // ABSOLUTE CLEANING // ================= // This method takes signal height (photons) times area ratio // ad the cleaning levels. // // The cleaning now is done in these levels (eg. 16, 20) // The meaning of the cleaning levels is essentially the same (the same cleaning // algorithm is used) but the cleaning is not done in different 'units' // // Example: // // MImgCleanStd clean(20, 16); // clean.SetMethod(MImgCleanStd::kAbsolulte); // // // Member Function: SetMethod() // ============================ // When you call the MImgCleanStd task, the default method is kStandard. // // If you want to switch to the kDemocratic method you have to // call this member function. // // Example: // // MImgCleanStd clean; // //creates a default Cleaning object, with default setting // // clean.SetMethod(MImgCleanStd::kDemocratic); // //now the method of cleaning is changed to Democratic // // // FIRST AND SECOND CLEANING LEVEL // =============================== // When you call the MImgCleanStd task, the default cleaning levels are // fCleanLvl1 = 3, fCleanLvl2 = 2.5. You can change them easily when you // create the MImgCleanStd object. // // Example: // // MImgCleanStd clean(Float_t lvl1,Float_t lvl2); // //creates a default cleaning object, but the cleaning levels are now // //lvl1 and lvl2. // // RING NUMBER // =========== // The standard cleaning procedure is such that it looks for the // informations of the boundary part of the shower only on the first // neighbors of the CORE pixels. // There is the possibility now to look not only at the firs neighbors // (first ring),but also further away, around the CORE pixels. All the new // pixels you can find with this method, are tested with the second level // of cleaning and have to have at least an USED neighbor. // // They will be also set as USED and will be taken into account during the // calculation of the image parameters. // The only way to distinguish them from the other USED pixels, is the // Ring number, that is bigger than 1. // // Example: You can decide how many rings you want to analyze using: // // MImgCleanStd clean; // //creates a default cleaning object (default number of rings =1) // clean.SetCleanRings(UShort_t r); // //now it looks r times around the CORE pixels to find new pixels with // //signal. // // // Input Containers: // MGeomCam // MPedPhotCam // MCerPhotEvt // // Output Containers: // MCerPhotEvt // ///////////////////////////////////////////////////////////////////////////// #include "MImgCleanStd.h" #include // atof #include // ofstream, SavePrimitive #include #include // TGFrame #include // TGLabel #include // TGTextEntry #include "MLog.h" #include "MLogManip.h" #include "MParList.h" #include "MCameraData.h" #include "MGeomPix.h" #include "MGeomCam.h" #include "MCerPhotPix.h" #include "MCerPhotEvt.h" #include "MGGroupFrame.h" // MGGroupFrame ClassImp(MImgCleanStd); using namespace std; enum { kImgCleanLvl1, kImgCleanLvl2 }; static const TString gsDefName = "MImgCleanStd"; static const TString gsDefTitle = "Task to perform image cleaning"; const TString MImgCleanStd::gsNamePedPhotCam="MPedPhotCam"; // default name of the 'MPedPhotCam' container // -------------------------------------------------------------------------- // // Default constructor. Here you can specify the cleaning method and levels. // If you don't specify them the 'common standard' values 3.0 and 2.5 (sigma // above mean) are used. // Here you can also specify how many rings around the core pixels you want // to analyze (with the fixed lvl2). The default value for "rings" is 1. // MImgCleanStd::MImgCleanStd(const Float_t lvl1, const Float_t lvl2, const char *name, const char *title) : fCleaningMethod(kStandard), fCleanLvl1(lvl1), fCleanLvl2(lvl2), fCleanRings(1), fNamePedPhotCam(gsNamePedPhotCam) { fName = name ? name : gsDefName.Data(); fTitle = title ? title : gsDefTitle.Data(); } // -------------------------------------------------------------------------- // // The first step of cleaning defines the CORE pixels. All the other pixels // are set as UNUSED and belong to RING 0. // After this point, only the CORE pixels are set as USED, with RING // number 1. // // NT 28/04/2003: now the option to use the standard method or the // democratic method is implemented: // // kStandard: This method looks for all pixels with an entry (photons) // that is three times bigger than the noise of the pixel // (default: 3 sigma, clean level 1) // // kDemocratic: this method looks for all pixels with an entry (photons) // that is n times bigger than the noise of the mean of the // inner pixels (default: 3 sigmabar, clean level 1) // // void MImgCleanStd::CleanStep1() { const TArrayD &data = fData->GetData(); // // check the number of all pixels against the noise level and // set them to 'unused' state if necessary // MCerPhotPix *pix; // Loop over all pixels MCerPhotEvtIter Next(fEvt, kFALSE); while ((pix=static_cast(Next()))) if (!pix->IsPixelUnmapped() && data[pix->GetPixId()] <= fCleanLvl1) pix->SetPixelUnused(); } // -------------------------------------------------------------------------- // // Check if the survived pixel have a neighbor, that also // survived, otherwise set pixel to unused (removes pixels without // neighbors). // void MImgCleanStd::CleanStep2() { MCerPhotPix *pix; // Loop over used pixels only TIter Next(*fEvt); while ((pix=static_cast(Next()))) { // get pixel id of this entry const Int_t idx = pix->GetPixId(); // check for 'used' neighbors of this pixel const MGeomPix &gpix = (*fCam)[idx]; const Int_t nnmax = gpix.GetNumNeighbors(); Bool_t hasNeighbor = kFALSE; //loop on the neighbors to check if they are used for (Int_t j=0; jIsPixelUsed(idx2)) { hasNeighbor = kTRUE; break; } } if (hasNeighbor == kFALSE) pix->SetPixelUnused(); } // // now we declare all pixels that survive as CorePixels // Next.Reset(); while ((pix=static_cast(Next()))) { if (pix->IsPixelUsed()) pix->SetPixelCore(); } } void MImgCleanStd::CleanStep3b(MCerPhotPix &pix) { const Int_t idx = pix.GetPixId(); // // check if the pixel's next neighbor is a core pixel. // if it is a core pixel set pixel state to: used. // MGeomPix &gpix = (*fCam)[idx]; const Int_t nnmax = gpix.GetNumNeighbors(); for (Int_t j=0; jIsPixelCore(idx2)) continue; pix.SetPixelUsed(); break; } } // -------------------------------------------------------------------------- // // NT: Add option "rings": default value = 1. // Look n (n>1) times for the boundary pixels around the used pixels. // If a pixel has more than 2.5 (clean level 2.5) sigma, // it is declared as used. // // If a value<2 for fCleanRings is used, no CleanStep4 is done. // void MImgCleanStd::CleanStep4(UShort_t r, MCerPhotPix &pix) { // // Skip events that have already a defined status; // if (pix.GetRing() != 0) return; // // check if the pixel's next neighbor is a used pixel. // if it is a used pixel set pixel state to: used, // and tell to which ring it belongs to. // const Int_t idx = pix.GetPixId(); MGeomPix &gpix = (*fCam)[idx]; const Int_t nnmax = gpix.GetNumNeighbors(); for (Int_t j=0; jGetPixById(idx2); if (!npix || !npix->IsPixelUsed() || npix->GetRing()>r-1 ) continue; pix.SetRing(r); break; } } // -------------------------------------------------------------------------- // // Look for the boundary pixels around the core pixels // if a pixel has more than 2.5 (clean level 2.5) sigma, and // a core neigbor, it is declared as used. // void MImgCleanStd::CleanStep3() { const TArrayD &data = fData->GetData(); for (UShort_t r=1; r(NextAll()))) { // // if pixel is a core pixel or unmapped, go to the next pixel // if (pix->IsPixelCore() || pix->IsPixelUnmapped()) continue; if (data[pix->GetPixId()] <= fCleanLvl2) continue; if (r==1) CleanStep3b(*pix); else CleanStep4(r, *pix); } } } // -------------------------------------------------------------------------- // // Check if MEvtHeader exists in the Parameter list already. // if not create one and add them to the list // Int_t MImgCleanStd::PreProcess (MParList *pList) { fCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam")); if (!fCam) { *fLog << err << "MGeomCam not found (no geometry information available)... aborting." << endl; return kFALSE; } fEvt = (MCerPhotEvt*)pList->FindObject(AddSerialNumber("MCerPhotEvt")); if (!fEvt) { *fLog << err << "MCerPhotEvt not found... aborting." << endl; return kFALSE; } fPed = (MPedPhotCam*)pList->FindObject(AddSerialNumber(fNamePedPhotCam), "MPedPhotCam"); if (!fPed) { *fLog << err << "MPedPhotCam not found... aborting." << endl; return kFALSE; } *fLog << all << "Name of MPedPhotCam container = " << fNamePedPhotCam << endl; fTime = (MArrivalTime*)pList->FindObject(AddSerialNumber("MArrivalTime")); if (!fTime && fCleaningMethod==kProbability) *fLog << warn << "MArrivalTime not found... probability cleaning done with signal only!" << endl; fData = (MCameraData*)pList->FindCreateObj(AddSerialNumber("MCameraData")); if (!fData) return kFALSE; Print(); return kTRUE; } // -------------------------------------------------------------------------- // // Cleans the image. // Int_t MImgCleanStd::Process() { switch (fCleaningMethod) { case kStandard: fData->CalcCleaningLevel(*fEvt, *fPed, *fCam); break; case kScaled: fData->CalcCleaningLevel2(*fEvt, *fPed, *fCam); break; case kDemocratic: fData->CalcCleaningLevelDemocratic(*fEvt, *fPed, *fCam); break; case kProbability: fData->CalcCleaningProbability(*fEvt, *fPed, *fCam, fTime); break; case kAbsolute: fData->CalcCleaningAbsolute(*fEvt, *fCam); break; default: break; } #ifdef DEBUG *fLog << all << "CleanStep 1" << endl; #endif CleanStep1(); #ifdef DEBUG *fLog << all << "CleanStep 2" << endl; #endif CleanStep2(); // For speed reasons skip the rest of the cleaning if no // action will be taken! if (fCleanLvl1>fCleanLvl2) { #ifdef DEBUG *fLog << all << "CleanStep 3" << endl; #endif CleanStep3(); } #ifdef DEBUG *fLog << all << "Calc Islands" << endl; #endif // Takes roughly 10% of the time fEvt->CalcIslands(*fCam); #ifdef DEBUG *fLog << all << "Done." << endl; #endif return kTRUE; } // -------------------------------------------------------------------------- // // Print descriptor and cleaning levels. // void MImgCleanStd::Print(Option_t *o) const { *fLog << all << GetDescriptor() << " using "; switch (fCleaningMethod) { case kDemocratic: *fLog << "democratic"; break; case kStandard: *fLog << "standard"; break; case kScaled: *fLog << "scaled"; break; case kProbability: *fLog << "probability"; break; case kAbsolute: *fLog << "absolute"; break; } *fLog << " cleaning initialized with level " << fCleanLvl1 << " and " << fCleanLvl2; *fLog << " (CleanRings=" << fCleanRings << ")" << endl; } // -------------------------------------------------------------------------- // // Create two text entry fields, one for each cleaning level and a // describing text line. // void MImgCleanStd::CreateGuiElements(MGGroupFrame *f) { // // Create a frame for line 3 and 4 to be able // to align entry field and label in one line // TGHorizontalFrame *f1 = new TGHorizontalFrame(f, 0, 0); TGHorizontalFrame *f2 = new TGHorizontalFrame(f, 0, 0); /* * --> use with root >=3.02 <-- * TGNumberEntry *fNumEntry1 = new TGNumberEntry(frame, 3.0, 2, M_NENT_LVL1, kNESRealOne, kNEANonNegative); TGNumberEntry *fNumEntry2 = new TGNumberEntry(frame, 2.5, 2, M_NENT_LVL1, kNESRealOne, kNEANonNegative); */ TGTextEntry *entry1 = new TGTextEntry(f1, "****", kImgCleanLvl1); TGTextEntry *entry2 = new TGTextEntry(f2, "****", kImgCleanLvl2); // --- doesn't work like expected (until root 3.02?) --- fNumEntry1->SetAlignment(kTextRight); // --- doesn't work like expected (until root 3.02?) --- fNumEntry2->SetAlignment(kTextRight); entry1->SetText("3.0"); entry2->SetText("2.5"); entry1->Associate(f); entry2->Associate(f); TGLabel *l1 = new TGLabel(f1, "Cleaning Level 1"); TGLabel *l2 = new TGLabel(f2, "Cleaning Level 2"); l1->SetTextJustify(kTextLeft); l2->SetTextJustify(kTextLeft); // // Align the text of the label centered, left in the row // with a left padding of 10 // TGLayoutHints *laylabel = new TGLayoutHints(kLHintsCenterY|kLHintsLeft, 10); TGLayoutHints *layframe = new TGLayoutHints(kLHintsCenterY|kLHintsLeft, 5, 0, 10); // // Add one entry field and the corresponding label to each line // f1->AddFrame(entry1); f2->AddFrame(entry2); f1->AddFrame(l1, laylabel); f2->AddFrame(l2, laylabel); f->AddFrame(f1, layframe); f->AddFrame(f2, layframe); f->AddToList(entry1); f->AddToList(entry2); f->AddToList(l1); f->AddToList(l2); f->AddToList(laylabel); f->AddToList(layframe); } // -------------------------------------------------------------------------- // // Process the GUI Events comming from the two text entry fields. // Bool_t MImgCleanStd::ProcessMessage(Int_t msg, Int_t submsg, Long_t param1, Long_t param2) { if (msg!=kC_TEXTENTRY || submsg!=kTE_ENTER) return kTRUE; TGTextEntry *txt = (TGTextEntry*)FindWidget(param1); if (!txt) return kTRUE; Float_t lvl = atof(txt->GetText()); switch (param1) { case kImgCleanLvl1: fCleanLvl1 = lvl; *fLog << "Cleaning level 1 set to " << lvl << " sigma." << endl; return kTRUE; case kImgCleanLvl2: fCleanLvl2 = lvl; *fLog << "Cleaning level 2 set to " << lvl << " sigma." << endl; return kTRUE; } return kTRUE; } // -------------------------------------------------------------------------- // // Implementation of SavePrimitive. Used to write the call to a constructor // to a macro. In the original root implementation it is used to write // gui elements to a macro-file. // void MImgCleanStd::StreamPrimitive(ofstream &out) const { out << " MImgCleanStd " << GetUniqueName() << "("; out << fCleanLvl1 << ", " << fCleanLvl2; if (fName!=gsDefName || fTitle!=gsDefTitle) { out << ", \"" << fName << "\""; if (fTitle!=gsDefTitle) out << ", \"" << fTitle << "\""; } out << ");" << endl; if (fCleaningMethod!=kStandard) { out << " " << GetUniqueName() << ".SetMethod(MImgCleanStd::k"; switch (fCleaningMethod) { case kScaled: out << "Scaled"; break; case kDemocratic: out << "Democratic"; break; case kProbability: out << "Probability"; break; case kAbsolute: out << "Absolute"; break; default: break; } out << ");" << endl; } if (fCleanRings!=1) out << " " << GetUniqueName() << ".SetCleanRings(" << fCleanRings << ");" << endl; if (gsNamePedPhotCam!=fNamePedPhotCam) out << " " << GetUniqueName() << ".SetNamePedPhotCam(\"" << fNamePedPhotCam << "\");" << endl; } // -------------------------------------------------------------------------- // // Read the setup from a TEnv, eg: // MImgCleanStd.CleanLevel1: 3.0 // MImgCleanStd.CleanLevel2: 2.5 // MImgCleanStd.CleanMethod: Standard, Scaled, Democratic, Probability, Absolute // MImgCleanStd.CleanRings: 1 // Int_t MImgCleanStd::ReadEnv(const TEnv &env, TString prefix, Bool_t print) { Bool_t rc = kFALSE; if (IsEnvDefined(env, prefix, "CleanRings", print)) { rc = kTRUE; SetCleanRings(GetEnvValue(env, prefix, "CleanRings", fCleanRings)); } if (IsEnvDefined(env, prefix, "CleanLevel1", print)) { rc = kTRUE; fCleanLvl1 = GetEnvValue(env, prefix, "CleanLevel1", fCleanLvl1); } if (IsEnvDefined(env, prefix, "CleanLevel2", print)) { rc = kTRUE; fCleanLvl2 = GetEnvValue(env, prefix, "CleanLevel2", fCleanLvl2); } if (IsEnvDefined(env, prefix, "CleanMethod", print)) { rc = kTRUE; TString s = GetEnvValue(env, prefix, "CleanMethod", ""); s.ToLower(); if (s.BeginsWith("standard")) SetMethod(kStandard); if (s.BeginsWith("scaled")) SetMethod(kScaled); if (s.BeginsWith("democratic")) SetMethod(kDemocratic); if (s.BeginsWith("probability")) SetMethod(kProbability); if (s.BeginsWith("absolute")) SetMethod(kAbsolute); } return rc; }