| 1 | #include "MCerPhotEvt.h"
|
|---|
| 2 |
|
|---|
| 3 | #include <math.h>
|
|---|
| 4 | #include <fstream.h>
|
|---|
| 5 |
|
|---|
| 6 | #include <TCanvas.h>
|
|---|
| 7 | #include <TClonesArray.h>
|
|---|
| 8 |
|
|---|
| 9 | #include "MLog.h"
|
|---|
| 10 | #include "MGeomCam.h"
|
|---|
| 11 | #include "MCamNeighbor.h"
|
|---|
| 12 | #include "MCamDisplay.h"
|
|---|
| 13 | #include "MHexagon.h"
|
|---|
| 14 |
|
|---|
| 15 | ClassImp(MCerPhotEvt)
|
|---|
| 16 |
|
|---|
| 17 | MCerPhotEvt::MCerPhotEvt(const char *name, const char *title ) : fType(0), fNbPixels(0)
|
|---|
| 18 | {
|
|---|
| 19 | // the default constructor
|
|---|
| 20 |
|
|---|
| 21 |
|
|---|
| 22 | *fName = name ? name : "MCerPhotEvt";
|
|---|
| 23 | *fTitle = name ? name : "(Number of Photon)-Event Information";
|
|---|
| 24 |
|
|---|
| 25 | fPixels = new TClonesArray ("MCerPhotPix", 577) ;
|
|---|
| 26 |
|
|---|
| 27 | //
|
|---|
| 28 | // FIXME: is this really necessary?
|
|---|
| 29 | //
|
|---|
| 30 | fPixels->Clear();
|
|---|
| 31 | }
|
|---|
| 32 |
|
|---|
| 33 | void MCerPhotEvt::Draw(Option_t* option)
|
|---|
| 34 | {
|
|---|
| 35 | // FIXME!!!
|
|---|
| 36 | //
|
|---|
| 37 |
|
|---|
| 38 | MCamDisplay *disp = new MCamDisplay(fType) ;
|
|---|
| 39 |
|
|---|
| 40 | disp->DrawPhotNum(this) ;
|
|---|
| 41 | }
|
|---|
| 42 |
|
|---|
| 43 | Int_t MCerPhotEvt::GetNbPixels()
|
|---|
| 44 | {
|
|---|
| 45 | return fNbPixels;
|
|---|
| 46 | }
|
|---|
| 47 |
|
|---|
| 48 | void MCerPhotEvt::AddPixel(Int_t id, Float_t nph, Float_t err)
|
|---|
| 49 | {
|
|---|
| 50 | //
|
|---|
| 51 | // add a new pixel to the list and increase the number
|
|---|
| 52 | // of valid pixels in the list by one
|
|---|
| 53 | //
|
|---|
| 54 | (*fPixels)[fNbPixels++] = new MCerPhotPix( id, nph, err);
|
|---|
| 55 | }
|
|---|
| 56 |
|
|---|
| 57 | void MCerPhotEvt::Clear(Option_t *)
|
|---|
| 58 | {
|
|---|
| 59 | //
|
|---|
| 60 | // reset counter and delete netries in list.
|
|---|
| 61 | //
|
|---|
| 62 | fNbPixels = 0 ;
|
|---|
| 63 | fPixels->Clear() ;
|
|---|
| 64 | }
|
|---|
| 65 |
|
|---|
| 66 | void MCerPhotEvt::Print(Option_t *)
|
|---|
| 67 | {
|
|---|
| 68 | const Int_t entries = fPixels->GetEntries();
|
|---|
| 69 |
|
|---|
| 70 | *fLog << "MCerPhotEvt::Print()" << endl
|
|---|
| 71 | << "Number of Pixels: " << fNbPixels
|
|---|
| 72 | << "(" << entries << ")"
|
|---|
| 73 | << endl ;
|
|---|
| 74 |
|
|---|
| 75 | for (Int_t il=0; il<entries; il++ )
|
|---|
| 76 | (*this)[il].Print();
|
|---|
| 77 | }
|
|---|
| 78 |
|
|---|
| 79 | void MCerPhotEvt::CleanLevel1()
|
|---|
| 80 | {
|
|---|
| 81 | //
|
|---|
| 82 | // This method looks for all pixels with an entry (photons)
|
|---|
| 83 | // that is three times bigger than the noise of the pixel
|
|---|
| 84 | //
|
|---|
| 85 |
|
|---|
| 86 | const Int_t entries = fPixels->GetEntries();
|
|---|
| 87 |
|
|---|
| 88 | //
|
|---|
| 89 | // check the number of all pixels against the noise level and
|
|---|
| 90 | // set them to 'unused' state if necessary
|
|---|
| 91 | //
|
|---|
| 92 | for (Int_t il=0; il<entries; il++ )
|
|---|
| 93 | {
|
|---|
| 94 | MCerPhotPix &pix = (*this)[il];
|
|---|
| 95 |
|
|---|
| 96 | const Float_t entry = pix.GetNumPhotons();
|
|---|
| 97 | const Float_t noise = pix.GetErrorPhot();
|
|---|
| 98 |
|
|---|
| 99 | if (entry < 3 * noise )
|
|---|
| 100 | pix.SetPixelUnused();
|
|---|
| 101 | }
|
|---|
| 102 | }
|
|---|
| 103 |
|
|---|
| 104 | void MCerPhotEvt::CleanLevel2()
|
|---|
| 105 | {
|
|---|
| 106 | //
|
|---|
| 107 | // check if the survived pixel have a neighbor, that also
|
|---|
| 108 | // survived
|
|---|
| 109 | //
|
|---|
| 110 |
|
|---|
| 111 | const Int_t entries = fPixels->GetEntries();
|
|---|
| 112 |
|
|---|
| 113 | for (Int_t il=0; il<entries; il++)
|
|---|
| 114 | {
|
|---|
| 115 | //
|
|---|
| 116 | // get entry il from list
|
|---|
| 117 | //
|
|---|
| 118 | MCerPhotPix &pix = (*this)[il];
|
|---|
| 119 |
|
|---|
| 120 | //
|
|---|
| 121 | // check if pixel is in use, if not goto next pixel in list
|
|---|
| 122 | //
|
|---|
| 123 | if (!pix.IsPixelUsed())
|
|---|
| 124 | continue;
|
|---|
| 125 |
|
|---|
| 126 | //
|
|---|
| 127 | // get pixel id of this entry
|
|---|
| 128 | //
|
|---|
| 129 | const Int_t id = pix.GetPixId() ;
|
|---|
| 130 |
|
|---|
| 131 | //
|
|---|
| 132 | // count number of next neighbors of this pixel which
|
|---|
| 133 | // state is 'used'
|
|---|
| 134 | //
|
|---|
| 135 | Int_t itest = 0 ;
|
|---|
| 136 | for (Int_t in=0 ; in < 6; in++ )
|
|---|
| 137 | {
|
|---|
| 138 | const Int_t id2 = fNN.GetNN(id, in) ;
|
|---|
| 139 |
|
|---|
| 140 | if (id2 < 0)
|
|---|
| 141 | continue;
|
|---|
| 142 |
|
|---|
| 143 | if (IsPixelUsed(id2))
|
|---|
| 144 | itest++ ;
|
|---|
| 145 | }
|
|---|
| 146 |
|
|---|
| 147 | //
|
|---|
| 148 | // check if no next neighbor has the state 'used'
|
|---|
| 149 | // set this pixel to 'unused', too.
|
|---|
| 150 | //
|
|---|
| 151 | if (itest==0)
|
|---|
| 152 | pix.SetPixelUnused();
|
|---|
| 153 | }
|
|---|
| 154 |
|
|---|
| 155 | //
|
|---|
| 156 | // now we declare all pixels that survive as CorePixels
|
|---|
| 157 | //
|
|---|
| 158 | for (Int_t il=0; il<entries; il++)
|
|---|
| 159 | {
|
|---|
| 160 | MCerPhotPix &pix = (*this)[il];
|
|---|
| 161 |
|
|---|
| 162 | if (pix.IsPixelUsed())
|
|---|
| 163 | pix.SetCorePixel();
|
|---|
| 164 | }
|
|---|
| 165 |
|
|---|
| 166 | }
|
|---|
| 167 |
|
|---|
| 168 | void MCerPhotEvt::CleanLevel3()
|
|---|
| 169 | {
|
|---|
| 170 | //
|
|---|
| 171 | // Look for the boundary pixels around the core pixels
|
|---|
| 172 | // if a pixel has more than 2.5 sigma, and a core neigbor
|
|---|
| 173 | // it is declared as used.
|
|---|
| 174 | //
|
|---|
| 175 | const Int_t entries = fPixels->GetEntries();
|
|---|
| 176 |
|
|---|
| 177 | for (Int_t il=0; il<entries; il++)
|
|---|
| 178 | {
|
|---|
| 179 | //
|
|---|
| 180 | // get pixel as entry il from list
|
|---|
| 181 | //
|
|---|
| 182 | MCerPhotPix &pix = (*this)[il];
|
|---|
| 183 |
|
|---|
| 184 | //
|
|---|
| 185 | // if pixel is a core pixel go to the next pixel
|
|---|
| 186 | //
|
|---|
| 187 | if (pix.IsCorePixel())
|
|---|
| 188 | continue;
|
|---|
| 189 |
|
|---|
| 190 | //
|
|---|
| 191 | // check the num of photons against the noise level
|
|---|
| 192 | //
|
|---|
| 193 | const Float_t entry = pix.GetNumPhotons();
|
|---|
| 194 | const Float_t noise = pix.GetErrorPhot();
|
|---|
| 195 |
|
|---|
| 196 | if (entry <= 2.5 * noise )
|
|---|
| 197 | continue;
|
|---|
| 198 |
|
|---|
| 199 | //
|
|---|
| 200 | // get pixel id of this entry
|
|---|
| 201 | //
|
|---|
| 202 | const Int_t id = pix.GetPixId();
|
|---|
| 203 |
|
|---|
| 204 | //
|
|---|
| 205 | // check if the pixel's next neighbor is a core pixel.
|
|---|
| 206 | // if it is a core pixel set pixel state to: used.
|
|---|
| 207 | //
|
|---|
| 208 | for (Int_t in=0; in<6 ; in++)
|
|---|
| 209 | {
|
|---|
| 210 | const Int_t id2 = fNN.GetNN(id, in);
|
|---|
| 211 |
|
|---|
| 212 | if (id2 <0)
|
|---|
| 213 | continue;
|
|---|
| 214 |
|
|---|
| 215 | if (!IsPixelCore(id2))
|
|---|
| 216 | continue;
|
|---|
| 217 |
|
|---|
| 218 | pix.SetPixelUsed();
|
|---|
| 219 |
|
|---|
| 220 | break ;
|
|---|
| 221 | }
|
|---|
| 222 | }
|
|---|
| 223 | }
|
|---|
| 224 |
|
|---|
| 225 |
|
|---|
| 226 | Bool_t MCerPhotEvt::IsPixelExisting(Int_t id)
|
|---|
| 227 | {
|
|---|
| 228 | //
|
|---|
| 229 | // Checks if in the pixel list is an entry with pixel id
|
|---|
| 230 | //
|
|---|
| 231 | const Int_t entries = fPixels->GetEntries();
|
|---|
| 232 |
|
|---|
| 233 | for (Int_t il=0; il<entries; il++)
|
|---|
| 234 | {
|
|---|
| 235 | if (id == (*this)[il].GetPixId())
|
|---|
| 236 | return kTRUE ;
|
|---|
| 237 | }
|
|---|
| 238 |
|
|---|
| 239 | return kFALSE ;
|
|---|
| 240 | }
|
|---|
| 241 |
|
|---|
| 242 | Bool_t MCerPhotEvt::IsPixelUsed(Int_t id)
|
|---|
| 243 | {
|
|---|
| 244 | //
|
|---|
| 245 | // Checks if in the pixel list is an entry with pixel id
|
|---|
| 246 | //
|
|---|
| 247 | const Int_t entries = fPixels->GetEntries();
|
|---|
| 248 |
|
|---|
| 249 | for (Int_t il=0; il<entries; il++ )
|
|---|
| 250 | {
|
|---|
| 251 | MCerPhotPix &pix = (*this)[il];
|
|---|
| 252 |
|
|---|
| 253 | if (id == pix.GetPixId() && pix.IsPixelUsed())
|
|---|
| 254 | return kTRUE ;
|
|---|
| 255 | }
|
|---|
| 256 |
|
|---|
| 257 | return kFALSE ;
|
|---|
| 258 | }
|
|---|
| 259 |
|
|---|
| 260 | Bool_t MCerPhotEvt::IsPixelCore(Int_t id)
|
|---|
| 261 | {
|
|---|
| 262 | //
|
|---|
| 263 | // Checks if in the pixel list is an entry with pixel id
|
|---|
| 264 | //
|
|---|
| 265 | const Int_t entries = fPixels->GetEntries();
|
|---|
| 266 |
|
|---|
| 267 | for (Int_t il=0; il<entries; il++ )
|
|---|
| 268 | {
|
|---|
| 269 | MCerPhotPix &pix = (*this)[il];
|
|---|
| 270 |
|
|---|
| 271 | if ( id == pix.GetPixId() && pix.IsCorePixel())
|
|---|
| 272 | return kTRUE ;
|
|---|
| 273 | }
|
|---|
| 274 |
|
|---|
| 275 | return kFALSE ;
|
|---|
| 276 | }
|
|---|
| 277 |
|
|---|
| 278 | Float_t MCerPhotEvt::GetMinNumPhotons()
|
|---|
| 279 | {
|
|---|
| 280 | //
|
|---|
| 281 | // get the minimum number of photons of all valid pixels in the list
|
|---|
| 282 | //
|
|---|
| 283 | if (fNbPixels <= 0)
|
|---|
| 284 | return -5. ;
|
|---|
| 285 |
|
|---|
| 286 | Float_t minval = (*this)[0].GetNumPhotons();
|
|---|
| 287 |
|
|---|
| 288 | Float_t testval;
|
|---|
| 289 | for (Int_t i=1 ; i<fNbPixels; i++ )
|
|---|
| 290 | {
|
|---|
| 291 | testval = (*this)[i].GetNumPhotons();
|
|---|
| 292 |
|
|---|
| 293 | if (testval < minval)
|
|---|
| 294 | minval = testval;
|
|---|
| 295 | }
|
|---|
| 296 |
|
|---|
| 297 | return minval;
|
|---|
| 298 | }
|
|---|
| 299 |
|
|---|
| 300 | Float_t MCerPhotEvt::GetMaxNumPhotons()
|
|---|
| 301 | {
|
|---|
| 302 | //
|
|---|
| 303 | // get the maximum number of photons of all valid pixels in the list
|
|---|
| 304 | //
|
|---|
| 305 | if (fNbPixels <= 0)
|
|---|
| 306 | return 50.;
|
|---|
| 307 |
|
|---|
| 308 | Float_t maxval = (*this)[0].GetNumPhotons();
|
|---|
| 309 |
|
|---|
| 310 | Float_t testval;
|
|---|
| 311 | for (Int_t i=1; i<fNbPixels; i++)
|
|---|
| 312 | {
|
|---|
| 313 | testval = (*this)[i].GetNumPhotons();
|
|---|
| 314 |
|
|---|
| 315 | if (testval > maxval)
|
|---|
| 316 | maxval = testval;
|
|---|
| 317 | }
|
|---|
| 318 | return maxval;
|
|---|
| 319 | }
|
|---|
| 320 |
|
|---|