source: branches/MarsMoreSimulationTruth/mtemp/mifae/library/MIslandsCalc.cc@ 20115

Last change on this file since 20115 was 5466, checked in by aliu, 20 years ago
*** empty log message ***
File size: 19.2 KB
Line 
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): Ester Aliu, 2/2004 <aliu@ifae.es>
19|
20! Last Update: 7/2004
21!
22!
23! Copyright: MAGIC Software Development, 2000-2004
24!
25!
26\* ======================================================================== */
27
28/////////////////////////////////////////////////////////////////////////////
29//
30// MIslandsCalc
31//
32// The Island Calc task calculates some islands parameters for each
33// of the events such as:
34//
35// - fPixNum // number of pixels in the island
36// - fSigToNoise // signal to noise of the island
37// - fTimeSpread // "time" of the island
38// - fMeanX // mean X position of the island
39// - fMeanY // mean Y position of the island
40// - fDist // dist between an island and the continent
41// - fLength // major axis of the island ellipse
42// - fWidth // minor axis of the island ellipse
43// - fDistL // dist divided by lenght of the larger island
44// - fDistW // dist divided by width of the larger island
45// - fDistS // dist divided by size of the larger island
46//
47// - fPixList // list of pixels in the island (TArrayI)
48// - fPeakPulse // mean arrival time of the pixels in the island (TArrayF)
49//
50// Input Containers:
51// MGeomCam
52// MCerPhotEvt
53// MPedestalCam
54// MArrivalTimeCam
55//
56// Output Containers:
57// MIslands
58// MImgIsland
59//
60/////////////////////////////////////////////////////////////////////////////
61#include "MIslandsCalc.h"
62
63#include <stdlib.h> // atof
64#include <fstream> // ofstream, SavePrimitive
65
66#include "MLog.h"
67#include "MLogManip.h"
68
69#include "MIslands.h"
70#include "MImgIsland.h"
71
72#include "MParList.h"
73
74#include "MGeomPix.h"
75#include "MGeomCam.h"
76
77#include "MCerPhotPix.h"
78#include "MCerPhotEvt.h"
79
80#include "MPedestalCam.h"
81#include "MPedestalPix.h"
82
83#include "MArrivalTimeCam.h"
84#include "MArrivalTimePix.h"
85
86ClassImp(MIslandsCalc);
87
88
89using namespace std;
90
91// --------------------------------------------------------------------------
92//
93// Default constructor.
94//
95MIslandsCalc::MIslandsCalc(const char* name, const char* title)
96 : fIsl(NULL)
97{
98 fName = name ? name : "MIslandsCalc";
99 fTitle = title ? title : "Calculate island parameters";
100}
101
102
103// --------------------------------------------------------------------------
104Int_t MIslandsCalc::PreProcess (MParList *pList)
105{
106 fCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
107 if (!fCam)
108 {
109 *fLog << dbginf << "MGeomCam not found (no geometry information available)... aborting." << endl;
110 return kFALSE;
111 }
112
113 fEvt = (MCerPhotEvt*)pList->FindObject(AddSerialNumber("MCerPhotEvt"));
114 if (!fEvt)
115 {
116 *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl;
117 return kFALSE;
118 }
119
120 fPed = (MPedestalCam*)pList->FindObject(AddSerialNumber("MPedestalCam"));
121 if (!fPed)
122 {
123 *fLog << dbginf << "MPedestalCam not found... aborting." << endl;
124 return kFALSE;
125 }
126
127 fTime = (MArrivalTimeCam*)pList->FindObject(AddSerialNumber("MArrivalTimeCam"));
128 if (!fTime)
129 {
130 *fLog << dbginf << "MArrivalTimeCam not found... aborting." << endl;
131 return kFALSE;
132 }
133
134 if (strlen(fIslName) > 0)
135 {
136 fIsl = (MIslands*)pList->FindCreateObj("MIslands", AddSerialNumber(fIslName));
137 }
138 else
139 {
140 fIsl = (MIslands*)pList->FindCreateObj(AddSerialNumber("MIslands"));
141 }
142 if (!fIsl)
143 return kFALSE;
144
145 return kTRUE;
146}
147
148
149Int_t MIslandsCalc::Process(){
150
151 fIsl->GetList()->Delete();
152 IslandPar();
153
154 return kTRUE;
155}
156
157
158Int_t MIslandsCalc::IslandPar(){
159
160 //calculates all the island parameters
161
162 const Int_t nPix=fCam->GetNumPixels();
163 const Int_t nVect=50;
164 Int_t numisl;
165
166 Int_t** vect;
167 vect = new Int_t*[nVect];
168 for(Int_t i=0;i<nVect;i++)
169 vect[i]= new Int_t[nPix];
170
171 Int_t* num;
172 num = new Int_t[nVect];
173
174 if (fIslandAlgorithm == 1)
175 Calc1(numisl,nVect,nPix,vect,num);
176 if (fIslandAlgorithm == 2)
177 Calc2(numisl,nVect,nPix,vect,num);
178
179 //set the number of islands in one event
180 fIsl->SetIslNum(numisl);
181
182 //cout << "code numisl " << numisl << endl;
183 //examine each island...
184
185 Float_t noise;
186 Float_t signal;
187
188 Int_t PixelNumIsl[numisl];
189 Float_t SigToNoise[numisl];
190 Float_t time[nPix];
191 Float_t timeVariance[numisl];
192 Float_t meanX[numisl];
193 Float_t meanY[numisl];
194 Float_t length[numisl];
195 Float_t width[numisl];
196 Float_t dist[numisl];
197 Float_t distL[numisl];
198 Float_t distW[numisl];
199 Float_t distS[numisl];
200
201
202 Float_t size[numisl], sizeLargeIsl, distance, alpha, alphaW, sizetot;
203
204 sizeLargeIsl = 0;
205 alphaW = 0;
206 sizetot = 0;
207
208 for(Int_t i = 1; i<=numisl ; i++)
209 {
210
211 MImgIsland *imgIsl = new MImgIsland;
212
213 imgIsl->InitSize(num[i]);
214
215 Int_t n = 0;
216
217 Float_t minTime = 10000.;
218 Float_t maxTime = 0.;
219
220 Float_t minX = 10000.;
221 Float_t maxX = 0.;
222
223 Float_t minY = 10000.;
224 Float_t maxY = 0.;
225
226 signal = 0;
227 noise = 0;
228
229 size[i-1] = 0;
230 meanX[i-1] = 0;
231 meanY[i-1] = 0;
232 dist[i-1] = 0;
233
234 PixelNumIsl[i-1] = 0;
235 timeVariance[i-1] = 0;
236
237 for(Int_t idx=0 ; idx<nPix ; idx++)
238 {
239 MCerPhotPix *pix = fEvt->GetPixById(idx);
240 if(!pix) continue;
241 const MGeomPix &gpix2 = (*fCam)[pix->GetPixId()];
242 const MPedestalPix &ped = (*fPed)[idx];
243 const MArrivalTimePix &timepix = (*fTime)[idx];
244 const Float_t nphot = pix->GetNumPhotons();
245
246 if (vect[i][idx]==1){
247
248 PixelNumIsl[i-1]++;
249 signal += nphot * (fCam->GetPixRatio(idx));
250 noise += pow(ped.GetPedestalRms(),2);
251
252 size[i-1] += nphot;
253 if (i == 1)
254 sizeLargeIsl += nphot;
255
256 meanX[i-1] += nphot * gpix2.GetX();
257 meanY[i-1] += nphot * gpix2.GetY();
258
259 time[i-1] = timepix.IsLoGainUsed() ? timepix.GetArrivalTimeLoGain() : timepix.GetArrivalTimeHiGain();
260
261 imgIsl->SetPixList(PixelNumIsl[i-1]-1,pix->GetPixId());
262 imgIsl->SetPeakPulse(PixelNumIsl[i-1]-1,time[i-1]);
263
264 //calculates the time spread only for core pixels
265 if (fEvt->IsPixelCore(idx)){
266
267 if (time[i-1] > maxTime){
268 maxTime = time[i-1];
269 maxX = gpix2.GetX();
270 maxY = gpix2.GetY();
271 }
272 if (time[i-1] < minTime){
273 minTime = time[i-1];
274 minX = gpix2.GetX();
275 minY = gpix2.GetY();
276 }
277
278 }
279 n++;
280 }
281 }
282
283 meanX[i-1] /= size[i-1];
284 meanY[i-1] /= size[i-1];
285
286 dist[i-1] = TMath::Power(meanX[i-1]-meanX[0],2) + TMath::Power(meanY[i-1]-meanY[i-1],2);
287 dist[i-1] = TMath::Sqrt(dist[i-1]);
288
289 //timeVariance[i-1] = (maxTime-minTime);
290
291 if (maxX!=minX && maxY!=minY)
292 timeVariance[i-1] = (maxTime-minTime)/sqrt(TMath::Power(maxX-minX,2) + TMath::Power(maxY-minY,2));
293 else
294 timeVariance[i-1] = -1;
295
296 //noise = 0, in the case of MC w/o noise
297 if (noise == 0) noise = 1;
298
299 SigToNoise[i-1]= (Float_t)signal/(Float_t)sqrt(noise);
300
301 imgIsl->SetPixNum(PixelNumIsl[i-1]);
302 imgIsl->SetSigToNoise(SigToNoise[i-1]);
303 imgIsl->SetTimeSpread(timeVariance[i-1]);
304 imgIsl->SetMeanX(meanX[i-1]);
305 imgIsl->SetMeanY(meanY[i-1]);
306 imgIsl->SetDist(dist[i-1]);
307 imgIsl->SetSizeIsl(size[i-1]);
308
309
310 // sanity check: if one island has 2 or less pixels
311 if (num[i]>2){
312
313
314 // calculate width and lenght of each island
315
316 Double_t corrxx=0; // [m^2]
317 Double_t corrxy=0; // [m^2]
318 Double_t corryy=0; // [m^2]
319
320 for(Int_t idx=0 ; idx<nPix ; idx++){
321
322 MCerPhotPix *pix = fEvt->GetPixById(idx);
323 if(!pix) continue;
324 const MGeomPix &gpix3 = (*fCam)[pix->GetPixId()];
325 const Float_t nphot = pix->GetNumPhotons();
326
327 if (vect[i][idx]==1){
328
329 const Float_t dx = gpix3.GetX() - meanX[i-1]; // [mm]
330 const Float_t dy = gpix3.GetY() - meanY[i-1]; // [mm]
331
332 corrxx += nphot * dx*dx; // [mm^2]
333 corrxy += nphot * dx*dy; // [mm^2]
334 corryy += nphot * dy*dy; // [mm^2]
335
336 }
337 }
338
339 const Double_t d0 = corryy - corrxx;
340 const Double_t d1 = corrxy*2;
341 const Double_t d2 = d0 + TMath::Sqrt(d0*d0 + d1*d1);
342 const Double_t tand = d2 / d1;
343 const Double_t tand2 = tand*tand;
344
345 const Double_t s2 = tand2+1;
346 const Double_t s = TMath::Sqrt(s2);
347
348 const Double_t CosDelta = 1.0/s; // need these in derived classes
349 const Double_t SinDelta = tand/s; // like MHillasExt
350
351 const Double_t axis1 = (tand2*corryy + d2 + corrxx)/s2/size[i-1];
352 const Double_t axis2 = (tand2*corrxx - d2 + corryy)/s2/size[i-1];
353
354 //
355 // fLength^2 is the second moment along the major axis of the ellipse
356 // fWidth^2 is the second moment along the minor axis of the ellipse
357 //
358 // From the algorithm we get: fWidth <= fLength is always true
359 //
360 // very small numbers can get negative by rounding
361 //
362 length[i-1] = axis1<0 ? 0 : TMath::Sqrt(axis1); // [mm]
363 width[i-1] = axis2<0 ? 0 : TMath::Sqrt(axis2); // [mm]
364
365
366 // alpha calculation
367
368 const Double_t mx = meanX[i-1]; // [mm]
369 const Double_t my = meanY[i-1]; // [mm]
370
371 //FIXME: xpos, ypos from MSrcPos
372 const Double_t xpos = 0.;
373 const Double_t ypos = 0.;
374
375 const Double_t sx = mx - xpos; // [mm]
376 const Double_t sy = my - ypos; // [mm]
377
378 const Double_t sd = SinDelta; // [1]
379 const Double_t cd = CosDelta; // [1]
380
381 //
382 // Distance from source position to center of ellipse.
383 // If the distance is 0 distance, Alpha is not specified.
384 // The calculation has failed and returnes kFALSE.
385 //
386 distance = TMath::Sqrt(sx*sx + sy*sy); // [mm]
387 if (distance==0){
388
389 for (Int_t l = 0; l< nVect; l++)
390 delete [] vect[l];
391
392 delete [] vect;
393 delete [] num;
394
395 return 1;
396 }
397
398 //
399 // Calculate Alpha and Cosda = cos(d,a)
400 // The sign of Cosda will be used for quantities containing
401 // a head-tail information
402 //
403 // *OLD* const Double_t arg = (sy-tand*sx) / (dist*sqrt(tand*tand+1));
404 // *OLD* fAlpha = asin(arg)*kRad2Deg;
405 //
406 const Double_t arg2 = cd*sx + sd*sy; // [mm]
407 if (arg2==0){
408
409 for (Int_t l = 0; l< nVect; l++)
410 delete [] vect[l];
411
412 delete [] vect;
413 delete [] num;
414
415 return 2;
416 }
417
418 const Double_t arg1 = cd*sy - sd*sx; // [mm]
419
420 //
421 // Due to numerical uncertanties in the calculation of the
422 // square root (dist) and arg1 it can happen (in less than 1e-5 cases)
423 // that the absolute value of arg exceeds 1. Because this uncertainty
424 // results in an Delta Alpha which is still less than 1e-3 we don't care
425 // about this uncertainty in general and simply set values which exceed
426 // to 1 saving its sign.
427 //
428 const Double_t arg = arg1/distance;
429 alpha = TMath::Abs(arg)>1 ? TMath::Sign(90., arg) : TMath::ASin(arg)*TMath::RadToDeg(); // [deg]
430
431 alphaW += alpha*size[i-1];
432 sizetot += size[i-1];
433
434 ////////////////////////////////////////
435
436 distL[i-1]=dist[i-1]/length[0];
437 distW[i-1]=dist[i-1]/width[0];
438 distS[i-1]= dist[i-1]/size[0];
439
440 imgIsl->SetLength(length[i-1]);
441 imgIsl->SetWidth(width[i-1]);
442
443 imgIsl->SetDistL(distL[i-1]);
444 imgIsl->SetDistW(distW[i-1]);
445 imgIsl->SetDistS(distS[i-1]);
446
447 imgIsl->SetAlpha(alpha);
448
449 }
450
451 fIsl->GetList()->Add(imgIsl);
452
453 }
454
455 fIsl->SetAlphaW(alphaW/sizetot);
456 //fIsl->SetReadyToSave();
457
458 for (Int_t l = 0; l< nVect; l++)
459 delete [] vect[l];
460
461 delete [] vect;
462 delete [] num;
463
464 return kTRUE;
465}
466
467//------------------------------------------------------------------------------------------
468void MIslandsCalc::Calc1(Int_t& numisl, const Int_t nv, const Int_t npix, Int_t** vect, Int_t* num){
469
470
471 /////////////////////////////
472 //
473 // ALGORITHM # 1
474 // counts the number of islands as you can see in
475 // the event display after doing the std image cleaning
476 //
477 /////////////////////////////
478
479 Int_t sflag;
480 Int_t control = 0;
481
482 Int_t nvect = 0;
483
484 numisl = 0;
485
486 Int_t zeros[nv];
487
488 for(Int_t m = 0; m < nv ; m++)
489 for(Int_t n = 0; n < npix ; n++)
490 vect[m][n] = 0;
491
492 for(Int_t n = 0; n < nv ; n++)
493 zeros[n] = 0;
494
495 //cout << "new event" <<endl;
496 MCerPhotPix *pix;
497
498 // Loop over used pixels only
499 TIter Next(*fEvt);
500
501 //after the interpolation of the pixels, these can be disordered by index. This is important for this algorithm of calculating islands
502 fEvt->Sort();
503
504 while ((pix=static_cast<MCerPhotPix*>(Next())))
505 {
506 const Int_t idx = pix->GetPixId();
507
508 const MGeomPix &gpix = (*fCam)[idx];
509 const Int_t nnmax = gpix.GetNumNeighbors();
510
511 if( fEvt->IsPixelUsed(idx))
512 {
513 //cout <<idx <<endl;
514 sflag = 0;
515
516 for(Int_t j=0; j < nnmax ; j++)
517 {
518 const Int_t idx2 = gpix.GetNeighbor(j);
519
520 if (idx2 < idx)
521 {
522 for(Int_t k = 1; k <= nvect; k++)
523 {
524 if (vect[k][idx2] == 1)
525 {
526 sflag = 1;
527 vect[k][idx] = 1;
528 }
529 }
530 }
531 }
532
533 if (sflag == 0)
534 {
535 nvect++;
536 vect[nvect][idx] = 1;
537 }
538
539 }
540 }
541
542 numisl = nvect;
543
544
545 // Repeated Chain Corrections
546
547 Int_t jmin = 0;
548
549 for(Int_t i = 1; i <= nvect; i++){
550 control=0;
551 for(Int_t j = i+1; j <= nvect; j++){
552 control = 0;
553 for(Int_t k = 0; k < npix; k++){
554 if (vect[i][k] == 1 && vect[j][k] == 1){
555 control = 1;
556 k=npix;
557 }
558 }
559 if (control == 1){
560 for(Int_t k = 0; k < npix; k++){
561 if(vect[j][k] == 1) vect[i][k] = 1;
562 vect[j][k] = 0;
563 zeros[j] = 1;
564 }
565 numisl = numisl-1;
566 }
567 }
568
569 for(Int_t j = 1; j <= i-1; j++){
570 for(Int_t k = 0; k < npix; k++){
571 if (vect[i][k] == 1 && vect[j][k] == 1){
572 control = 2;
573 jmin=j;
574 k=npix;
575 j=i;
576 }
577 }
578
579 if (control == 2){
580 for (Int_t k = 0; k < npix; k++){
581 if(vect[i][k]==1) vect[jmin][k]=1;
582 vect[i][k] = 0;
583 zeros[i] = 1;
584 }
585 numisl = numisl-1;
586 }
587 }
588 }
589
590 Int_t pixMAX = 0;
591 Int_t idMAX = 1;
592 Int_t l = 1;
593 Int_t numpixels;
594
595 for(Int_t i = 1; i<= nvect ; i++)
596 {
597 numpixels = 0;
598
599 if (zeros[i] == 0)
600 {
601 for(Int_t k=0; k<npix; k++)
602 {
603 vect[l][k] = vect[i][k];
604 if (vect[l][k] == 1)
605 numpixels++;
606 }
607
608 num[l] = numpixels;
609
610 if (numpixels>pixMAX)
611 {
612 pixMAX = numpixels;
613 idMAX = l;
614 }
615 l++;
616 }
617 }
618
619
620 //the larger island will correspond to the 1st component of the vector
621
622 num[nvect+1] = num[1];
623 num[1] = num[idMAX];
624 num[idMAX] = num[nvect+1];
625
626
627 for(Int_t k = 0; k<npix; k++)
628 {
629 vect[nvect+1][k] = vect[1][k];
630 vect[1][k] = vect[idMAX][k];
631 vect[idMAX][k] = vect[nvect+1][k];
632 }
633}
634
635//------------------------------------------------------------------------------------------
636
637void MIslandsCalc::Calc2(Int_t& numisl, const Int_t nv, const Int_t npix, Int_t** vect, Int_t* num){
638
639
640 /////////////////////////////
641 //
642 // ALGORITHM # 2
643 // counts the number of islands considering as the same
644 // islands the ones separated for 2 or less pixels
645 //
646 /////////////////////////////
647
648 Int_t sflag;
649 Int_t control;
650
651 Int_t nvect = 0;
652 numisl = 0;
653
654 Int_t zeros[nv];
655
656 Int_t kk[npix];
657
658 for(Int_t m = 0; m < nv ; m++)
659 for(Int_t n = 0; n < npix ; n++)
660 vect[m][n] = 0;
661
662 for(Int_t n = 0; n < nv ; n++)
663 zeros[n] = 0;
664
665 for(Int_t n = 0; n < npix ; n++)
666 kk[n] = 0;
667
668 MCerPhotPix *pix;
669
670 //after the interpolation of the pixels, these can be disordered by index. This is important for this algorithm of calculating islands
671 fEvt->Sort();
672
673 // 1st loop over used pixels only
674 TIter Next0(*fEvt);
675
676 while ((pix=static_cast<MCerPhotPix*>(Next0())))
677 {
678 const Int_t idx = pix->GetPixId();
679
680 const MGeomPix &gpix = (*fCam)[idx];
681 const Int_t nnmax = gpix.GetNumNeighbors();
682
683 if( fEvt->IsPixelUsed(idx))
684 {
685 kk[idx] = 1 ;
686 for(Int_t j=0; j< nnmax; j++)
687 {
688 kk[gpix.GetNeighbor(j)] = 1;
689 }
690 }
691
692 }
693
694
695 //2nd loop over all pixels
696 TIter Next(*fEvt);
697
698 while ((pix=static_cast<MCerPhotPix*>(Next())))
699 {
700 const Int_t idx = pix->GetPixId();
701
702 const MGeomPix &gpix = (*fCam)[idx];
703 const Int_t nnmax = gpix.GetNumNeighbors();
704
705 if ( kk[idx] > 0)
706 {
707 sflag = 0;
708
709 for(Int_t j=0; j < nnmax ; j++)
710 {
711 const Int_t idx2 = gpix.GetNeighbor(j);
712
713 if (idx2 < idx)
714 {
715 for(Int_t k = 1; k <= nvect; k++)
716 {
717 if (vect[k][idx2] == 1)
718 {
719 sflag = 1;
720 vect[k][idx] = 1;
721 }
722 }
723 }
724 }
725
726 if (sflag == 0)
727 {
728 nvect++;
729 vect[nvect][idx] = 1;
730 }
731
732 }
733 }
734
735 numisl = nvect;
736
737 // Repeated Chain Corrections
738
739 Int_t jmin = 0;
740
741 for(Int_t i = 1; i <= nvect; i++){
742 control=0;
743 for(Int_t j = i+1; j <= nvect; j++){
744 control = 0;
745 for(Int_t k = 0; k < npix; k++){
746 if (vect[i][k] == 1 && vect[j][k] == 1){
747 control = 1;
748 k=npix;
749 }
750 }
751 if (control == 1){
752 for(Int_t k = 0; k < npix; k++){
753 if(vect[j][k] == 1) vect[i][k] = 1;
754 vect[j][k] = 0;
755 zeros[j] = 1;
756 }
757 numisl = numisl-1;
758 }
759 }
760
761 for(Int_t j = 1; j <= i-1; j++){
762 for(Int_t k = 0; k < npix; k++){
763 if (vect[i][k] == 1 && vect[j][k] == 1){
764 control = 2;
765 jmin=j;
766 k=npix;
767 j=i;
768 }
769 }
770
771 if (control == 2){
772 for (Int_t k = 0; k < npix; k++){
773 if(vect[i][k]==1) vect[jmin][k]=1;
774 vect[i][k] = 0;
775 zeros[i] = 1;
776 }
777 numisl = numisl-1;
778 }
779 }
780 }
781
782 Int_t l = 1;
783 Int_t numpixels;
784 Int_t pixMAX = 0;
785 Int_t idMAX = 1;
786
787 for(Int_t i = 1; i<= nvect ; i++)
788 {
789 numpixels = 0;
790
791 if (zeros[i] == 0)
792 {
793 for(Int_t k = 0; k<npix; k++)
794 {
795 vect[l][k] = vect[i][k];
796 if (vect[l][k] == 1)
797 numpixels++;
798 }
799
800 num[l] = numpixels;
801
802 if (numpixels>pixMAX)
803 {
804 pixMAX = numpixels;
805 idMAX = l;
806 }
807 l++;
808 }
809 }
810
811
812 //the larger island will correspond to the 1st component of the vector
813
814 num[nvect +1] = num[1];
815 num[1] = num[idMAX];
816 num[idMAX]=num[nvect+1];
817
818 for(Int_t k = 0; k<npix; k++)
819 {
820 vect[nvect+1][k] = vect[1][k];
821 vect[1][k] = vect[idMAX][k];
822 vect[idMAX][k] = vect[nvect+1][k];
823 }
824
825}
826
Note: See TracBrowser for help on using the repository browser.