// TEST SCRIPT FOR GetContourList Function // Author: Josh de Bever // CSI Medical Physics Group // The University of Western Ontario // London, Ontario, Canada // Date: Aug. 26, 2004 CERNtest(){ gROOT->Reset(); gSystem->Load("libphysics"); const Double_t PI = TMath::Pi(); TCanvas* c = new TCanvas("c","Sample",0,0,700,600); Int_t i, j, pti, TotalConts; Int_t nZsamples = 80; Int_t nPhiSamples = 80; Double_t HofZwavelength = 4.0; // 4meters Double_t dZ = HofZwavelength/(Double_t)(nZsamples - 1); Double_t dPhi = 2*PI/(Double_t)(nPhiSamples - 1); TArrayD z(nZsamples); TArrayD HofZ(nZsamples); TArrayD phi(nPhiSamples); TArrayD FofPhi(nPhiSamples); TGraph2D *streamFn = new TGraph2D(); // Discritize Z and Phi Values for ( i = 0; i < nZsamples; i++) { z[i] = (i)*dZ - HofZwavelength/2.0; HofZ[i] = SawTooth(z[i], HofZwavelength); printf("%d: Z = % f H(Z) = % f\n", i, z[i], HofZ[i]); } for(Int_t i=0; i < nPhiSamples; i++){ phi[i] = (i)*dPhi; FofPhi[i] = sin(phi[i]); } printf("Input Functions Calculated\n"); // Load Graph2D for (Int_t i = 0; i < nZsamples; i++) { for(Int_t j = 0; j < nPhiSamples; j++){ streamFn->SetPoint(pti , z[i], phi[j], HofZ[i]*FofPhi[j]); pti++; } } gStyle->SetPalette(1); streamFn->Draw("surf1"); // ************************************************************************* // THIS IS WHERE THE ACTION HAPPENS, WILL TAKE A LONG TIME // TO GET FIRST CONTOUR, THEN DO THE REST VERY FAST // *********************************************************************** TotalConts = DoGetContour(streamFn, 10); printf("\n\n\tExtracted %d Graphs as Contours\n", TotalConts); } // ********************************************************************************* // ********************************************************************************* Int_t DoGetContour(TGraph2D *streamFn, Int_t nContours) { TList* contour = NULL; // List of Graphs for one contour Int_t nTotalConts = 0; // total # of graphs in all contours Int_t i = 0; // contour counter Double_t level = 0; // current level Double_t max = 0, min = 0; Double_t spacing = 0; if (streamFn->GetN() < 1) return -1 ; // Points haven't been properly calculated max = streamFn->GetZmax(); min = streamFn->GetZmin(); spacing = (max - min)/(nContours); // Get Contours from stream function for(i = 0; i < nContours; i++){ level = (max - 0.5*spacing) - i*spacing; printf("\n\n******************************\n Level = %f \n", level); printf("...GETTING CONTOUR LIST\n" ); // Get list of Graphs making up this contour level contour = streamFn->GetContourList(level); printf("\tThis Level Has %d Contours\n", contour->GetSize();); printf("******************************\n"); nTotalConts += contour->GetSize(); } // END FOR return nTotalConts; // success } // ********************************************************************************* // ********************************************************************************* Double_t SawTooth(Double_t x, Double_t WaveLen){ /* This function is specific to a sawtooth function with period WaveLen symetric about x = 0 and with amplitude = 1. Each segment is 1/4 of the wavelength. | /\ | / \ | / \ | / \ /--------\--------/------------ |\ / | \ / | \ / | \/ */ Double_t y; if ( (x < -WaveLen/2) || (x > WaveLen/2)) y = -99999999; // ERROR X OUT OF BOUNDS, Not elegant if (x <= -WaveLen/4){ y = x + 2.0; } else if ( (x > -WaveLen/4) && (x <= WaveLen/4)){ y = -x ; } else if( ( x > WaveLen/4) && (x <= WaveLen/2)){ y = x - 2.0; } return y; }