Index: /trunk/MagicSoft/Mars/mtemp/mifae/Changelog =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/Changelog (revision 4285) +++ /trunk/MagicSoft/Mars/mtemp/mifae/Changelog (revision 4286) @@ -18,4 +18,16 @@ -*-*- END OF LINE -*-*- + + 2004/05/27 Ester Aliu + * programs/makeHillas.cc + - add the possibility of using more than one algorithm to + calculate the islands and use different algorithms for counting + islands + + * library/MIslands.[h,cc],MIslandCalc.cc + - add a new island cleaning which consists in removing all the + islands except the larger one + - add a new algorithm of counting islands consisting in consider + as the same islands those islands separated by 2 or less pixels 2004/06/02 Javier Rico Index: /trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandCalc.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandCalc.cc (revision 4285) +++ /trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandCalc.cc (revision 4286) @@ -4,5 +4,5 @@ ! Author(s): Ester Aliu, 2/2004 ! -! Last Update: 5/2004 +! Last Update: 6/2004 ! \* ======================================================================== */ @@ -124,5 +124,26 @@ Int_t MIslandCalc::Process() { - + + if (fIslandAlgorithm == 1) + Calc1(); + + if (fIslandAlgorithm == 2) + Calc2(); + return kTRUE; + +} + + +Int_t MIslandCalc::Calc1(){ + + + ///////////////////////////// + // + // ALGORITHM # 1 + // counts the number of algorithms as you can see in + // the event display after doing the std image cleaning + // + ///////////////////////////// + Float_t noise; Float_t signal; @@ -333,7 +354,252 @@ fIsl->SetReadyToSave(); - return 1; - + return kTRUE; } +Int_t MIslandCalc::Calc2(){ + + + ///////////////////////////// + // + // ALGORITHM # 2 + // counts the number of islands considering as the same + // islands the ones separated for 2 or less pixels + // + ///////////////////////////// + + Float_t noise; + Float_t signal; + + Int_t npix = 577; + + Int_t sflag; + Int_t control; + + Int_t nvect = 0; + Int_t fIslNum = 0; + + Int_t vect[50][577]; + + Int_t zeros[50]; + + Int_t kk[577]; + + for(Int_t m = 0; m < 50 ; m++) + for(Int_t n = 0; n < npix ; n++) + vect[m][n] = 0; + + for(Int_t n = 0; n < 50 ; n++) + zeros[n] = 0; + + for(Int_t n = 0; n < 577 ; n++) + kk[n] = 0; + + MCerPhotPix *pix; + + //1st loop over all pixels + MCerPhotEvtIter Next0(fEvt, kFALSE); + + while ((pix=static_cast(Next0()))) + { + const Int_t idx = pix->GetPixId(); + + const MGeomPix &gpix = (*fCam)[idx]; + const Int_t nnmax = gpix.GetNumNeighbors(); + + if( fEvt->IsPixelUsed(idx)) + { + kk[idx] = 1 ; + for(Int_t j=0; j< nnmax; j++) + { + kk[gpix.GetNeighbor(j)] = 1; + } + } + + } + + //2nd loop over all pixels + MCerPhotEvtIter Next(fEvt, kFALSE); + + while ((pix=static_cast(Next()))) + { + const Int_t idx = pix->GetPixId(); + + const MGeomPix &gpix = (*fCam)[idx]; + const Int_t nnmax = gpix.GetNumNeighbors(); + + if ( kk[idx] > 0) + { + sflag = 0; + + for(Int_t j=0; j < nnmax ; j++) + { + const Int_t idx2 = gpix.GetNeighbor(j); + + if (idx2 < idx) + { + for(Int_t k = 1; k <= nvect; k++) + { + if (vect[k][idx2] == 1) + { + sflag = 1; + vect[k][idx] = 1; + } + } + } + } + + if (sflag == 0) + { + nvect++; + vect[nvect][idx] = 1; + } + + } + } + + fIslNum = nvect; + + + // Repeated Chain Corrections + + for(Int_t i = 1; i <= nvect; i++) + { + for(Int_t j = i+1; j <= nvect; j++) + { + control = 0; + for(Int_t k = 0; k < npix; k++) + { + + if (vect[i][k] == 1 && vect[j][k] == 1) + { + control = 1; + break; + } + } + if (control == 1) + { + for(Int_t k = 0; k < npix; k++) + { + if(vect[j][k] == 1) + vect[i][k] = 1; + vect[j][k] = 0; + zeros[j] = 1; + } + fIslNum = fIslNum-1; + } + + } + } + + + Int_t l = 1; + + for(Int_t i = 1; i<= nvect ; i++) + { + if (zeros[i] == 0) + { + for(Int_t k = 0; kSetIslNum(fIslNum); + + //examine each island... + Int_t fPixNum[fIslNum]; + Float_t fSigToNoise[fIslNum]; + Float_t time[577]; + Float_t timeVariance[fIslNum]; + + //reset the "sets" functions + if (fIslNum <1) + fIsl->SetIslNum(0); + + for(Int_t i = 0; i<20 ;i++) + { + fIsl->SetPixNum(i,-1); + fIsl->SetSigToNoise(i,-1); + fIsl->SetTimeSpread(i,-1); + + for(Int_t idx = 0; idxSetIslId(idx, -1); + fIsl->SetArrivalTime(i, idx, -1 ); + } + } + + + for(Int_t i = 1; i<=fIslNum ; i++) + { + Int_t n = 0; + Int_t ncore = 0; + + Float_t MIN = 10000; + Float_t MAX = 0; + + signal = 0; + noise = 0; + fPixNum[i-1] = 0; + timeVariance[i-1] = 0; + + for(Int_t idx=0 ; idxGetPixById(idx); + const MPedestalPix &ped = (*fPed)[idx]; + const MArrivalTimePix &timepix = (*fTime)[idx]; + + if (pix == NULL) break; + + if (vect[i][idx] ==1 && fEvt->IsPixelUsed(idx)){ + + fPixNum[i-1]++; + signal += pix->GetNumPhotons() * (fCam->GetPixRatio(idx)); + noise += pow(ped.GetPedestalRms(),2); + + time[n] = timepix.IsLoGainUsed() ? timepix.GetArrivalTimeLoGain() : timepix.GetArrivalTimeHiGain(); + + + if (fEvt->IsPixelCore(idx)){ + + if (time[n] > MAX) + MAX = time[n]; + if (time[n] < MIN) + MIN = time[n]; + + ncore++; + } + + fIsl->SetIslId(idx, i-1); + fIsl->SetArrivalTime(i-1, idx, time[n]); + + n++; + } + + } + + // Float_t mean = timeVariance[i-1]/ncore; + + timeVariance[i-1] = 0; + + timeVariance[i-1] = (MAX - MIN)/ncore; + timeVariance[i-1] = (MAX - MIN)/ncore; + + fSigToNoise[i-1]= (Float_t)signal/(Float_t)sqrt(noise); + + fIsl->SetPixNum(i-1,fPixNum[i-1]); + fIsl->SetSigToNoise(i-1,fSigToNoise[i-1]); + fIsl->SetTimeSpread(i-1, timeVariance[i-1]); + + } + + fIsl->SetReadyToSave(); + + return 1; +} Index: /trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandCalc.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandCalc.h (revision 4285) +++ /trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandCalc.h (revision 4286) @@ -31,7 +31,11 @@ TString fIslName; // name of the 'MIslands' container - + + Int_t fIslandAlgorithm; + Int_t PreProcess(MParList *plist); Int_t Process(); + Int_t Calc1(); // algorithm of counting islands #1 + Int_t Calc2(); // algorithm of counting islands #2 public: @@ -39,4 +43,6 @@ void SetOutputName(TString outname) { fIslName = outname; } + void SetAlgorithm(Int_t m) {fIslandAlgorithm = m;} + ClassDef(MIslandCalc, 0) // task doing the image cleaning }; Index: /trunk/MagicSoft/Mars/mtemp/mifae/programs/makeHillas.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/programs/makeHillas.cc (revision 4285) +++ /trunk/MagicSoft/Mars/mtemp/mifae/programs/makeHillas.cc (revision 4286) @@ -65,8 +65,8 @@ // initial and final time slices to be used in signal extraction -const Byte_t hifirst = 0; -const Byte_t hilast = 13; +const Byte_t hifirst = 1; +const Byte_t hilast = 14; const Byte_t lofirst = 3; -const Byte_t lolast = 12; +const Byte_t lolast = 14; // declaration of variables read from datacards @@ -86,4 +86,5 @@ Float_t lnew = 40; Int_t kmethod = 1; +Int_t kalgorithm = 1; Int_t nfiles = 0; @@ -238,5 +239,8 @@ isl2.SetName("MIslands2"); - if (islflag == 1 || islflag == 2) + MIslands isl3; + isl3.SetName("MIslands3"); + + if (islflag == 1 || islflag == 2 || islflag == 3) { plist4.AddToList(&timecam); @@ -247,4 +251,9 @@ { plist4.AddToList(&isl2); + } + + if (islflag == 3) + { + plist4.AddToList(&isl3); } @@ -276,4 +285,5 @@ MIslandCalc island; island.SetOutputName("MIslands1"); + island.SetAlgorithm(kalgorithm); MBadPixelsTreat interpolatebadpixels; @@ -286,4 +296,8 @@ MIslandCalc island2; island2.SetOutputName("MIslands2"); + island2.SetAlgorithm(kalgorithm); + + MIslandCalc island3; + island3.SetOutputName("MIslands3"); @@ -307,8 +321,11 @@ write->AddContainer("MSrcPosCam" , "Parameters"); - if (islflag == 1 || islflag == 2) + if (islflag == 1 || islflag == 2 || islflag == 3) write->AddContainer("MIslands1" , "Parameters"); if (islflag == 2) write->AddContainer("MIslands2" , "Parameters"); + if (islflag == 3) + write->AddContainer("MIslands3" , "Parameters"); + if(display) @@ -319,4 +336,6 @@ if (islflag == 2) disphillas->SetIslandsName("MIslands2"); + if (islflag == 3) + disphillas->SetIslandsName("MIslands3"); } @@ -329,5 +348,5 @@ tlist4.AddToList(&clean); - if (islflag == 1 || islflag == 2) + if (islflag == 1 || islflag == 2 || islflag == 3) { tlist4.AddToList(&timecalc); @@ -340,5 +359,12 @@ tlist4.AddToList(&island2); } - + + if (islflag == 3) + { + tlist4.AddToList(&islclean); + tlist4.AddToList(&island3); + } + + //tlist4.AddToList(&blind2); tlist4.AddToList(&hcalc); @@ -475,4 +501,7 @@ { ifun >> islflag; + + // if (islflag == 1 || islflag == 2) + ifun >> kalgorithm; } @@ -522,5 +551,5 @@ cout << "Cleaning level: ("<