static void AEMarchingCubes(uint8_t* voxels, uint32_t stepX, uint32_t stepY, uint32_t stepZ, double isovalue, uint32_t w, uint32_t h, uint32_t d, AEMCTriangleArray* triangleArray) {
    // Run the processCube function on every cube in the grid
	for(uint32_t x = stepX; x < w-2*stepX; x += stepX) {
		for(uint32_t y = stepY; y < h-2*stepY; y += stepY) {
			for(uint32_t z = stepZ; z < d-2*stepZ; z += stepZ) {
				/*GRIDCELL c = {{
					{x,y,z, 
                        (double)(voxels[Index(x+stepX,y,z,w,h,d)]-voxels[Index(x-stepX,y,z,w,h,d)]) / -stepX,
                        (double)(voxels[Index(x,y+stepY,z,w,h,d)]-voxels[Index(x,y-stepY,z,w,h,d)]) / -stepY,
                        (double)(voxels[Index(x,y,z+stepZ,w,h,d)]-voxels[Index(x,y,z-stepZ,w,h,d)]) / -stepZ
                    },
					{x+stepX,y,z, 
                        (double)(voxels[Index(x+2*stepX,y,z,w,h,d)]-voxels[Index(x,y,z,w,h,d)]) / -stepX,
                        (double)(voxels[Index(x+stepX,y+stepY,z,w,h,d)]-voxels[Index(x+stepX,y-stepY,z,w,h,d)]) / -stepY,
                        (double)(voxels[Index(x+stepX,y,z+stepZ,w,h,d)]-voxels[Index(x+stepX,y,z-stepZ,w,h,d)]) / -stepZ
                    },
					{x+stepX,y,z+stepZ, 
                        (double)(voxels[Index(x+2*stepX,y,z+stepZ,w,h,d)]-voxels[Index(x,y,z+stepZ,w,h,d)]) / -stepX,
                        (double)(voxels[Index(x+stepX,y+stepY,z+stepZ,w,h,d)]-voxels[Index(x+stepX,y-stepY,z+stepZ,w,h,d)]) / -stepY,
                        (double)(voxels[Index(x+stepX,y,z+2*stepZ,w,h,d)]-voxels[Index(x+stepX,y,z,w,h,d)]) / -stepZ
                    },
					{x,y,z+stepZ, 
                        (double)(voxels[Index(x+stepX,y,z+stepZ,w,h,d)]-voxels[Index(x-stepX,y,z+stepZ,w,h,d)]) / -stepX,
                        (double)(voxels[Index(x,y+stepY,z+stepZ,w,h,d)]-voxels[Index(x,y-stepY,z+stepZ,w,h,d)]) / -stepY,
                        (double)(voxels[Index(x,y,z+2*stepZ,w,h,d)]-voxels[Index(x,y,z,w,h,d)]) / -stepZ
                    },
					{x,y+stepY,z, 
                        (double)(voxels[Index(x+stepX,y+stepY,z,w,h,d)]-voxels[Index(x-stepX,y+stepY,z,w,h,d)]) / -stepX,
                        (double)(voxels[Index(x,y+2*stepY,z,w,h,d)]-voxels[Index(x,y,z,w,h,d)]) / -stepY,
                        (double)(voxels[Index(x,y+stepY,z+stepZ,w,h,d)]-voxels[Index(x,y+stepY,z-stepZ,w,h,d)]) / -stepZ
                    },
					{x+stepX,y+stepY,z, 
                        (double)(voxels[Index(x+2*stepX,y+stepY,z,w,h,d)]-voxels[Index(x+stepX,y+stepY,z,w,h,d)]) / -stepX,
                        (double)(voxels[Index(x+stepX,y+2*stepY,z,w,h,d)]-voxels[Index(x+stepX,y,z,w,h,d)]) / -stepY,
                        (double)(voxels[Index(x+stepX,y+stepY,z+stepZ,w,h,d)]-voxels[Index(x+stepX,y+stepY,z-stepZ,w,h,d)]) / -stepZ
                    },
					{x+stepX,y+stepY,z+stepZ, 
                        (double)(voxels[Index(x+2*stepX,y+stepY,z+stepZ,w,h,d)]-voxels[Index(x,y+stepY,z+stepZ,w,h,d)]) / -stepX,
                        (double)(voxels[Index(x+stepX,y+2*stepY,z+stepZ,w,h,d)]-voxels[Index(x+stepX,y,z+stepZ,w,h,d)]) / -stepY,
                        (double)(voxels[Index(x+stepX,y+stepY,z+2*stepZ,w,h,d)]-voxels[Index(x+stepX,y+stepY,z,w,h,d)]) / -stepZ
                    },
					{x,y+stepY,z+stepZ, 
                        (double)(voxels[Index(x+stepX,y+stepY,z+stepZ,w,h,d)]-voxels[Index(x-stepX,y+stepY,z+stepZ,w,h,d)]) / -stepX,
                        (double)(voxels[Index(x,y+2*stepY,z+stepZ,w,h,d)]-voxels[Index(x,y,z+stepZ,w,h,d)]) / -stepY,
                        (double)(voxels[Index(x,y+stepY,z+2*stepZ,w,h,d)]-voxels[Index(x,y+stepY,z,w,h,d)]) / -stepZ
                    }
				},{
					voxels[Index(x,y,z,w,h,d)],
					voxels[Index(x+stepX,y,z,w,h,d)],
					voxels[Index(x+stepX,y,z+stepZ,w,h,d)],
					voxels[Index(x,y,z+stepZ,w,h,d)],
					voxels[Index(x,y+stepY,z,w,h,d)],
					voxels[Index(x+stepX,y+stepY,z,w,h,d)],
					voxels[Index(x+stepX,y+stepY,z+stepZ,w,h,d)],
					voxels[Index(x,y+stepY,z+stepZ,w,h,d)]
				}};*/
				GRIDCELL c={{
					{x,y,z,0,0,0},
					{x+stepX,y,z,0,0,0},
					{x+stepX,y,z+stepZ,0,0,0},
					{x,y,z+stepZ,0,0,0},
					{x,y+stepY,z,0,0,0},
					{x+stepX,y+stepY,z,0,0,0},
					{x+stepX,y+stepY,z+stepZ,0,0,0},
					{x,y+stepY,z+stepZ,0,0,0}
				},{
					voxels[Index(x,y,z,w,h,d)],
					voxels[Index(x+stepX,y,z,w,h,d)],
					voxels[Index(x+stepX,y,z+stepZ,w,h,d)],
					voxels[Index(x,y,z+stepZ,w,h,d)],
					voxels[Index(x,y+stepY,z,w,h,d)],
					voxels[Index(x+stepX,y+stepY,z,w,h,d)],
					voxels[Index(x+stepX,y+stepY,z+stepZ,w,h,d)],
					voxels[Index(x,y+stepY,z+stepZ,w,h,d)]
				}};
				///Polygonise(c, isovalue);
				TRIANGLE triangles[5];
				size_t count=Polygonise(c, isovalue, triangles);
				for (size_t i=0; i<count; i++) {
					AEArrayAddBytes(triangleArray, triangles+i);
				}
			}
		}
	}
}

void compare_patches(double eps, int ndim, double *patch1, int lo1[], int hi1[],
                     int dims1[],double *patch2, int lo2[], int hi2[], 
                     int dims2[])
					           
{
	int i,j, elems=1;	
	int subscr1[MAXDIMS], subscr2[MAXDIMS];
        double diff,max;
    int offset1 = 0, offset2 = 0;

	for(i=0;i<ndim;i++){   /* count # of elements & verify consistency of both patches */
		int diff = hi1[i]-lo1[i];
		assert(diff == (hi2[i]-lo2[i]));
		assert(diff < dims1[i]);
		assert(diff < dims2[i]);
		elems *= diff+1;
		subscr1[i]= lo1[i];
		subscr2[i]=lo2[i];
	}

	
	/* compare element values in both patches */ 
	for(j=0; j< elems; j++){ 
		int idx1, idx2;
		
		idx1 = Index(ndim, subscr1, dims1);	 /* calculate element Index from a subscript */
		idx2 = Index(ndim, subscr2, dims2);

		if(j==0){
			offset1 =idx1;
			offset2 =idx2;
		}
		idx1 -= offset1;
		idx2 -= offset2;
		

                diff = patch1[idx1] - patch2[idx2];
                max  = MAX(ABS(patch1[idx1]),ABS(patch2[idx2]));
                if(max == 0. || max <eps) max = 1.; 

		if(eps < ABS(diff)/max){
			char msg[48];
			sprintf(msg,"(proc=%d):%f",me,patch1[idx1]);
			print_subscript("ERROR: a",ndim,subscr1,msg);
			sprintf(msg,"%f\n",patch2[idx2]);
			print_subscript(" b",ndim,subscr2,msg);
                        printf("\nA = %f B = %f\n", patch1[idx1], patch2[idx2]);
                        fflush(stdout);
                        sleep(1);
                        fprintf(stderr, "Bailing out\n");
                        MPI_Abort(MPI_COMM_WORLD, 1);
		}

		{ /* update subscript for the patches */
		   update_subscript(ndim, subscr1, lo1,hi1, dims1);
		   update_subscript(ndim, subscr2, lo2,hi2, dims2);
		}
	}
				 
	

	/* make sure we reached upper limit */
	/*for(i=0;i<ndim;i++){ 
		assert(subscr1[i]==hi1[i]);
		assert(subscr2[i]==hi2[i]);
	}*/ 
}

main()
{     
      FILE *f;
      char filename[50];
      printf("input the name of the file:");
      scanf("%s",filename);//the file's name can be absolute path,include the .txt example: D:\\file1\\file2\\test.txt
      f=fopen(filename,"r");
      if (f==NULL)//if can not open the file,stop the program
      {
         printf("OMG,the file is not existed!\n");
         system("pause");
         return;
         }
      get_the_map(f);//get the information of map
      getchar();
      printf("input the start street name:");
      gets(start_street);
      printf("input the start house number:");
      scanf("%d",&start_house);
      getchar();
      printf("input the destination street name:");
      gets(destination_street);
      printf("input the destination house number:");
      scanf("%d",&destination_house);
      
      Insert_start_destination(street_num,start_house,destination_house,start_street,destination_street);
      if (start_unique=='\0' || destination_unique=='\0')//if we don't get the unique number for start and destination nodes,we stop.
      {
           system("pause");
           return 0;
           }
      if (start_unique == destination_unique)//if the start node has the same coordinate as the destination node.stop.
      {
           printf("the start point is the same as destination!\n");
           system("pause");
           return 0;
           }
      int k=1;
      int i,j;
      Dijkstra();//Dijkstra is to store each node's father unique numbers
      struct search_path BFS_array[50];
      /*the destination node has his father nodes,and these father nodes also have their own father nodes.
        So it's a tree,and I use the BFS method to find every shortest path.
      */
      int label;
      label=Index(destination_unique);
      BFS_array[1].unique=destination_unique;
      i=1;
      int from,present;//"from" and "present" are labels in the array 
      from=2;
      j=from;
      present=from+adj_list[label].father_num;
      while (i<=adj_list[label].father_num)//begin at the destination,so that store the destination's father firstly
      {
            BFS_array[j].unique=adj_list[label].father_unique_num[i];
            BFS_array[j].pre_label=1;
            j++;
            i++;
            }
      get_BFS_list(BFS_array,from,present);//"from" to store the present first father's position,and "present" to store the array's present position
      i=1;
      int counter=0;
      while(BFS_array[i].unique!='\0')//"counter" is to get the total number of the shortest pathes
      {
            if (BFS_array[i].unique==start_unique)
                  counter++;
            i++;
      }
      struct each_path shortest_ways[counter+1];
      get_the_pathes(shortest_ways,BFS_array);
      int ways=counter;
      for (i=1;i<=ways;i++)//for each shortest path,store each node's distance between he and his father
      {
          path_link present;
          present=&shortest_ways[i];
          float pre_x,pre_y;
          if (present->unique_num==1)//if the unique number is 1,do some special operations
          {
                pre_x=start_x;
                pre_y=start_y;
          }
          else
          {
                present->x=present->unique_num/100;
                present->y=present->unique_num%100;
                pre_x=present->x;
                pre_y=present->y;
          }
          present=present->next_node;
          while(present!=NULL)
          {
                float X,Y;
                if (present->unique_num==2)//if the next node is destination node,also do some special operations
                {
                      X=pre_x-destination_x;
                      Y=pre_y-destination_y;
                      present->distance=sqrt(pow(X,2)+pow(Y,2))*mile;
                }
                else
                {
//.........这里部分代码省略.........

//.........这里部分代码省略.........

		u32 address, vaddress, size, offset, unknown;
		char name[512], container[512];
		if (four_columns)
		{
			// sometimes there is no unknown number, and sometimes it is because it is an entry of something else
			if (length > 37 && line[37]==' ')
			{
				unknown = 0;
				sscanf(line, "%08x %08x %08x %08x %511s", &address, &size, &vaddress, &offset, name);
				char *s = strstr(line, "(entry of ");
				if (s)
				{
					sscanf(s + 10, "%511s", container);
					char *s2 = (strchr(container, ')'));
					if (s2 && container[0]!='.')
					{
						s2[0] = '\0';
						strcat(container, "::");
						strcat(container, name);
						strcpy(name, container);
					}
				}
			}
			else
			{
				sscanf(line, "%08x %08x %08x %08x %i %511s", &address, &size, &vaddress, &offset, &unknown, name);
			}
		}
		// some entries in the table have a function name followed by " (entry of " followed by a container name, followed by ")"
		// instead of a space followed by a number followed by a space followed by a name
		else if (length > 27 && line[27] != ' ' && strstr(line, "(entry of "))
		{
			unknown = 0;
			sscanf(line, "%08x %08x %08x %511s", &address, &size, &vaddress, name);
			char *s = strstr(line, "(entry of ");
			if (s)
			{
				sscanf(s + 10, "%511s", container);
				char *s2 = (strchr(container, ')'));
				if (s2 && container[0] != '.')
				{
					s2[0] = '\0';
					strcat(container, "::");
					strcat(container, name);
					strcpy(name, container);
				}
			}
		}
		else
		{
			sscanf(line, "%08x %08x %08x %i %511s", &address, &size, &vaddress, &unknown, name);
		}

		const char *namepos = strstr(line, name);
		if (namepos != nullptr) //would be odd if not :P
			strcpy(name, namepos);
		name[strlen(name) - 1] = 0;

		// we want the function names only .... TODO: or do we really? aren't we wasting information here?
		for (size_t i = 0; i < strlen(name); i++)
		{
			if (name[i] == ' ') name[i] = 0x00;
			if (name[i] == '(') name[i] = 0x00;
		}

		// Check if this is a valid entry.
		if (strcmp(name, ".text") != 0 && strcmp(name, ".init") != 0 && strlen(name) > 0)
		{
			vaddress |= 0x80000000;
			bool good = !bad;
			if (!good)
			{
				// check for BLR before function
				u32 opcode = PowerPC::HostRead_Instruction(vaddress - 4);
				if (opcode == 0x4e800020)
				{
					// check for BLR at end of function
					opcode = PowerPC::HostRead_Instruction(vaddress + size - 4);
					if (opcode == 0x4e800020)
						good = true;
				}
			}
			if (good)
			{
				++good_count;
				AddKnownSymbol(vaddress | 0x80000000, size, name); // ST_FUNCTION
			}
			else
			{
				++bad_count;
			}
		}
	}

	Index();
	if (bad)
		SuccessAlertT("Loaded %d good functions, ignored %d bad functions.", good_count, bad_count);
	return true;
}

void cSkeleton<_DataType>::ConnectingFlaggedVoxels()
{
	int				i, j, k, loc[3], NextPos_i[3], NumContinuousVoxels;
	unsigned char	VoxelFlag_uc;


	printf ("Connecting Flagged Voxels ... \n");
	fflush(stdout);


#ifdef	DEBUG_CONNECTING
//	int		MaxNumContinuousVoxels = 0;
#endif

	for (k=1; k<Depth_mi-1; k++) {
		for (j=1; j<Height_mi-1; j++) {
			for (i=1; i<Width_mi-1; i++) {
	
				loc[0] = Index(i, j, k);
				if (VoxelFlags_muc[loc[0]]==FLAG_NONUNIFORM) {


#ifdef	DEBUG_CONNECTING
					printf ("Start (%d %d %d)\n", i, j, k);
#endif
					NextPos_i[0] = i;
					NextPos_i[1] = j;
					NextPos_i[2] = k;

					NumContinuousVoxels = 0;
					do {

						NextPos_i[0] = (int)((float)NextPos_i[0] + GVFDistance_mf[loc[0]*3 + 0] + 0.5);
						NextPos_i[1] = (int)((float)NextPos_i[1] + GVFDistance_mf[loc[0]*3 + 1] + 0.5);
						NextPos_i[2] = (int)((float)NextPos_i[2] + GVFDistance_mf[loc[0]*3 + 2] + 0.5);

						loc[0] = Index (NextPos_i[0], NextPos_i[1], NextPos_i[2]);
						NumContinuousVoxels++;
						
#ifdef	DEBUG_CONNECTING
						printf ("(%d %d %d)", NextPos_i[0], NextPos_i[1], NextPos_i[2]);
						DisplayFlag(VoxelFlags_muc[loc[0]]);
						printf (",  ");
#endif
						
						VoxelFlag_uc = VoxelFlags_muc[loc[0]];
						if (VoxelFlag_uc == FLAG_CONNECTED || 
							VoxelFlag_uc == FLAG_NONUNIFORM || 
							VoxelFlag_uc == FLAG_EMPTY) break;
						
						VoxelFlags_muc[loc[0]] = FLAG_CONNECTED;
						
					} while (VoxelFlag_uc != FLAG_NONUNIFORM);

				}
			}
		}
	}

	for (i=0; i<WHD_mi; i++) {
		if (VoxelFlags_muc[i]==FLAG_NONUNIFORM) VoxelFlags_muc[i] = FLAG_CONNECTED;
	}

	int NumConnectedVoxels = 0;
	for (i=0; i<WHD_mi; i++) {
		if (VoxelFlags_muc[i]==FLAG_CONNECTED) NumConnectedVoxels++;
	}
	printf ("Num. Connected Voxels = %d\n", NumConnectedVoxels);
	printf ("Num. Connected / NumSegmented = %f %%\n", (float)NumConnectedVoxels/NumSegmentedVoxels_mi*100.0);
	printf ("Num. Connected / Total Num. Voxels = %f %%\n", (double)NumConnectedVoxels/WHD_mi*100.0);
	fflush (stdout);
	
}

void FitterUtilsSimultaneousExpOfPolyTimesX::fit(bool wantplot, bool constPartReco,
      double fracPartReco_const,
      ofstream& out, TTree* t, bool update, string plotsfile)
{

   //***************Get the PDFs from the workspace

   TFile fw(workspacename.c_str());   
   RooWorkspace* workspace = (RooWorkspace*)fw.Get("workspace");
   RooRealVar *B_plus_M = workspace->var("B_plus_M");
   RooRealVar *misPT = workspace->var("misPT");
   RooRealVar *l1Kee = workspace->var("l1Kee");
   RooRealVar *l2Kee = workspace->var("l2Kee");
   RooRealVar *l3Kee = workspace->var("l3Kee");
   RooRealVar *l4Kee = workspace->var("l4Kee");
   RooRealVar *l5Kee = workspace->var("l5Kee");
   RooRealVar *l1KeeGen = workspace->var("l1KeeGen");
   RooRealVar *l2KeeGen = workspace->var("l2KeeGen");
   RooRealVar *l3KeeGen = workspace->var("l3KeeGen");
   RooRealVar *l4KeeGen = workspace->var("l4KeeGen");
   RooRealVar *l5KeeGen = workspace->var("l5KeeGen");
   RooRealVar *fractionalErrorJpsiLeak = workspace->var("fractionalErrorJpsiLeak");


   RooRealVar l1Kemu(*l1Kee);
   l1Kemu.SetName("l1Kemu"); l1Kemu.SetTitle("l1Kemu");    
   RooRealVar l2Kemu(*l2Kee);
   l2Kemu.SetName("l2Kemu"); l2Kemu.SetTitle("l2Kemu");    
   RooRealVar l3Kemu(*l3Kee);
   l3Kemu.SetName("l3Kemu"); l3Kemu.SetTitle("l3Kemu");    
   RooRealVar l4Kemu(*l4Kee);
   l4Kemu.SetName("l4Kemu"); l4Kemu.SetTitle("l4Kemu");    
   RooRealVar l5Kemu(*l5Kee);
   l5Kemu.SetName("l5Kemu"); l5Kemu.SetTitle("l5Kemu");    


   RooHistPdf *histPdfSignalZeroGamma = (RooHistPdf *) workspace->pdf("histPdfSignalZeroGamma");
   RooHistPdf *histPdfSignalOneGamma = (RooHistPdf *) workspace->pdf("histPdfSignalOneGamma");
   RooHistPdf *histPdfSignalTwoGamma = (RooHistPdf *) workspace->pdf("histPdfSignalTwoGamma");
   RooHistPdf *histPdfPartReco = (RooHistPdf *) workspace->pdf("histPdfPartReco");
   RooHistPdf *histPdfJpsiLeak(0);
   if(nGenJpsiLeak>0) histPdfJpsiLeak = (RooHistPdf *) workspace->pdf("histPdfJpsiLeak");

   //Here set in the Kemu PDF the parameters that have to be shared

   RooExpOfPolyTimesX* combPDF = new RooExpOfPolyTimesX("combPDF", "combPDF",  *B_plus_M, *misPT,  *l1Kee, *l2Kee, *l3Kee, *l4Kee, *l5Kee);
   RooExpOfPolyTimesX* KemuPDF = new RooExpOfPolyTimesX("kemuPDF", "kemuPDF",  *B_plus_M, *misPT,  l1Kemu, *l2Kee, *l3Kee, *l4Kee, *l5Kee);



   RooWorkspace* workspaceGen = (RooWorkspace*)fw.Get("workspaceGen");
   RooDataSet* dataGenSignalZeroGamma = (RooDataSet*)workspaceGen->data("dataGenSignalZeroGamma");
   RooDataSet* dataGenSignalOneGamma = (RooDataSet*)workspaceGen->data("dataGenSignalOneGamma");
   RooDataSet* dataGenSignalTwoGamma = (RooDataSet*)workspaceGen->data("dataGenSignalTwoGamma");
   RooDataSet* dataGenPartReco = (RooDataSet*)workspaceGen->data("dataGenPartReco");
   RooDataSet* dataGenComb = (RooDataSet*)workspaceGen->data("dataGenComb");
   RooDataSet* dataGenKemu = (RooDataSet*)workspaceGen->data("dataGenKemu");
   RooDataSet* dataGenJpsiLeak(0);
   if(nGenJpsiLeak>0) dataGenJpsiLeak = (RooDataSet*)workspaceGen->data("dataGenJpsiLeak");


   if(wantplot)
   {
      //**************Must get the datasets

      RooDataSet* dataSetSignalZeroGamma = (RooDataSet*)workspace->data("dataSetSignalZeroGamma");
      RooDataSet* dataSetSignalOneGamma = (RooDataSet*)workspace->data("dataSetSignalOneGamma");
      RooDataSet* dataSetSignalTwoGamma = (RooDataSet*)workspace->data("dataSetSignalTwoGamma");
      RooDataSet* dataSetPartReco = (RooDataSet*)workspace->data("dataSetPartReco");
      RooDataSet* dataSetComb = (RooDataSet*)workspace->data("dataSetComb");
      RooDataSet* dataSetJpsiLeak = (RooDataSet*)workspace->data("dataSetJpsiLeak");

      //**************Plot all the different components

      cout<<"dataGenSignalZeroGamma: "<<dataGenSignalZeroGamma<<endl;
      PlotShape(*dataSetSignalZeroGamma, *dataGenSignalZeroGamma, *histPdfSignalZeroGamma, plotsfile, "cSignalZeroGamma", *B_plus_M, *misPT);
      PlotShape(*dataSetSignalOneGamma, *dataGenSignalOneGamma, *histPdfSignalOneGamma, plotsfile, "cSignalOneGamma", *B_plus_M, *misPT);
      PlotShape(*dataSetSignalTwoGamma, *dataGenSignalTwoGamma, *histPdfSignalTwoGamma, plotsfile, "cSignalTwoGamma", *B_plus_M, *misPT);
      PlotShape(*dataSetPartReco, *dataGenPartReco, *histPdfPartReco, plotsfile, "cPartReco", *B_plus_M, *misPT);
      PlotShape(*dataSetComb, *dataGenComb, *combPDF, plotsfile, "cComb", *B_plus_M, *misPT);
      if(nGenJpsiLeak>1) PlotShape(*dataSetJpsiLeak, *dataGenJpsiLeak, *histPdfJpsiLeak, plotsfile, "cJpsiLeak", *B_plus_M, *misPT);
   }

   //***************Merge datasets

   RooDataSet* dataGenTot(dataGenPartReco);
   dataGenTot->append(*dataGenSignalZeroGamma);
   dataGenTot->append(*dataGenSignalOneGamma);
   dataGenTot->append(*dataGenSignalTwoGamma);
   dataGenTot->append(*dataGenComb);
   if(nGenJpsiLeak>0) dataGenTot->append(*dataGenJpsiLeak);

   //**************Create index category and join samples

   RooCategory category("category", "category");
   category.defineType("Kee");
   category.defineType("Kemu");

   RooDataSet dataGenSimultaneous("dataGenSimultaneous", "dataGenSimultaneous", RooArgSet(*B_plus_M, *misPT), Index(category), Import("Kee", *dataGenTot), Import("Kemu", *dataGenKemu));

//.........这里部分代码省略.........

LUA_API void lua_rawgeti (lua_State *L, int index, int n) {
  StkId o = Index(L, index);
  LUA_ASSERT(ttype(o) == LUA_TTABLE, "table expected");
  *L->top = *luaH_getnum(hvalue(o), n);
  api_incr_top(L);
}

LUA_API void lua_settable (lua_State *L, int index) {
  StkId t = Index(L, index);
  StkId top = L->top;
  luaV_settable(L, t, top-2);
  L->top = top-2;  /* pop index and value */
}

bool ElevationMap::fuseAll()
{
  ROS_DEBUG("Requested to fuse entire elevation map.");
  boost::recursive_mutex::scoped_lock scopedLock(fusedMapMutex_);
  return fuse(Index(0, 0), fusedMap_.getSize());
}

LUA_API void lua_gettable (lua_State *L, int index) {
  StkId t = Index(L, index);
  StkId top = L->top;
  *(top-1) = *luaV_gettable(L, t);
  L->top = top;  /* tag method may change top */
}

// Column major grid dissection order for the sparse 
// QR decomposition algorithm.
void grid_dissection_order( const grid2 *g0, int *perm )
{
  int width = Width(*g0);
  int height = Height(*g0);
  int size = width*height;

  int mid, x, y, i, *p;
  grid2 g1, g2, gs;
  g1 = g2 = gs = *g0;

  if( width >= 4 || height >= 4 ) {
#ifdef DEBUG_DISSECTION_ORDER
    printf( "dissecting %dx%d grid based at %d:\n", width, height, g0->base ); 
#endif
    if( width < height ) {
      mid = height / 2;

      // 1st grid
      g1.ran[1][1] = g0->ran[1][0] + mid-1;

      // 2nd grid
      g2.base = g0->base + Size(g1);
      g2.ran[1][0] = g0->ran[1][0] + mid+1;
      
      // Separating grid
      gs.base = g2.base + Size(g2);
      gs.ran[1][0] = g0->ran[1][0] + mid-1;
      gs.ran[1][1] = g0->ran[1][0] + mid+1;

#ifdef DEBUG_DISSECTION_ORDER
      printf( "\tg1: base = %d,   \tran_x = [%d,%d),   \tran_y = [%d,%d)\n", 
	      g1.base, g1.ran[0][0], g1.ran[0][1], 
	      g1.ran[1][0], g1.ran[1][1] ); 
      printf( "\tgs: base = %d,   \tran_x = [%d,%d),   \tran_y = [%d,%d)\n", 
	      gs.base, gs.ran[0][0], gs.ran[0][1], 
	      gs.ran[1][0], gs.ran[1][1] ); 
      printf( "\tg2: base = %d,   \tran_x = [%d,%d),   \tran_y = [%d,%d)\n", 
	      g2.base, g2.ran[0][0], g2.ran[0][1], 
	      g2.ran[1][0], g2.ran[1][1] ); 
      printf( "\n" );
      printf( "Mapping Nodes:\n" );
#endif

      for( i=g0->base, p=perm+i; i < size; ++i, ++p ) {
	Node(*g0,i,x,y);	
#ifdef DEBUG_DISSECTION_ORDER
	printf( "\tNode %d @ (%d,%d) -> ", i, x, y );
	if( y >= g2.ran[1][0] ) printf("%d (g2)\n", Index(g2,x,y));
	else if( y >= gs.ran[1][0] ) printf("%d (gs)\n", Index(gs,x,y));
	else printf("%d (g1)\n", Index(g1,x,y));
#endif
	if( y >= g2.ran[1][0] ) {*p = Index(g2,x,y); continue;}
	if( y >= gs.ran[1][0] ) {*p = Index(gs,x,y); continue;}

	*p = Index(g1,x,y);
      } // for( i, p )
      
      grid_dissection_order( &g1, perm );
      grid_dissection_order( &g2, perm );
    } else { // width >= height
      mid = width / 2;

      // 1st grid
      g1.ran[0][1] = g0->ran[0][0] + mid-1;

      // 2nd grid
      g2.base = g0->base + Size(g1);
      g2.ran[0][0] = g0->ran[0][0] + mid+1;
      
      // Separating grid
      gs.base = g2.base + Size(g2);
      gs.ran[0][0] = g0->ran[0][0] + mid-1;
      gs.ran[0][1] = g0->ran[0][0] + mid+1;

#ifdef DEBUG_DISSECTION_ORDER
      printf( "\tg1: base = %d,   \tran_x = [%d,%d),   \tran_y = [%d,%d)\n", 
	      g1.base, g1.ran[0][0], g1.ran[0][1], 
	      g1.ran[1][0], g1.ran[1][1] ); 
      printf( "\tgs: base = %d,   \tran_x = [%d,%d),   \tran_y = [%d,%d)\n", 
	      gs.base, gs.ran[0][0], gs.ran[0][1], 
	      gs.ran[1][0], gs.ran[1][1] ); 
      printf( "\tg2: base = %d,   \tran_x = [%d,%d),   \tran_y = [%d,%d)\n", 
	      g2.base, g2.ran[0][0], g2.ran[0][1], 
	      g2.ran[1][0], g2.ran[1][1] ); 
      printf( "\n" );
      printf( "Mapping Nodes:\n" );
#endif

      i = Index(*g0,gs.ran[0][0],gs.ran[1][0]);
      for( p=perm+i; i < size; ++i, ++p ) {
	Node(*g0,i,x,y);
#ifdef DEBUG_DISSECTION_ORDER
	printf( "\tNode %d @ (%d,%d) -> ", i, x, y );
	if( x >= g2.ran[0][0] ) printf("%d (g2)\n", Index(g2,x,y));
	else if( x >= gs.ran[0][0] ) printf("%d (gs)\n", Index(gs,x,y));
	else printf("%d (g1)\n", *p);
#endif
	
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  std::cout<<" Data File 11 loaded"<<std::endl;

  RooDataSet* MCDataAll= new RooDataSet("MCDataAll","MCDataAll",Args);

  MCDataAll->append(*MCData12);
  MCDataAll->append(*MCData11);
  
  RooPlot* massFrame = MCBMass.frame(Title("Data Import Check"),Bins(50));
  MCDataAll->plotOn(massFrame);
  
  RooPlot *BDTFrame = BDT_response.frame(Title("BDT Cut Check"),Bins(50));
  MCDataAll->plotOn(BDTFrame);
  TCanvas C;
  C.Divide(2,1);
  C.cd(1);
  massFrame->Draw();
  C.cd(2);
  BDTFrame->Draw();
  C.SaveAs("ImportChecks.eps");

  //________________________________MAKE MCROODATACATEGORIES__________________________________

  RooDataSet* MCBData=(RooDataSet*)MCDataAll->reduce(RooArgSet(MCBMass));
  MCBData->Print("v");
  
  RooDataSet* MCEtaData=(RooDataSet*)MCDataAll->reduce(RooArgSet(MCEtaMass));
  MCEtaData->Print("v");

  RooCategory MCMassType("MCMassType","MCMassType") ;
  MCMassType.defineType("B") ;
  MCMassType.defineType("Eta") ;
  
  // Construct combined dataset in (x,sample)
  RooDataSet MCcombData("MCcombData","MC combined data",Args,Index(MCMassType),Import("B",*MCBData),Import("Eta",*MCEtaData));

  
  //=============================================== MC FIT MODEL===================================
  
  RooRealVar Mean("Mean","Mean",5279.29,5276.0,5284.00);
  RooRealVar Sigma("Sigma","Sigma",20.54,17.0,24.8);
  RooRealVar LAlpha("LAlpha","LAlpha",-1.064,-2.5,0.0);
  RooRealVar RAlpha("RAlpha","RAlpha",1.88,0.0,5.0);
  RooRealVar LN("LN","LN",13.0,0.0,40.0);
  RooRealVar RN("RN","RN",2.56,0.0,6.0);

  RooCBShape CBLeft("CBLeft","CBLeft",MCBMass,Mean,Sigma,LAlpha,LN);
  
  RooCBShape CBRight("CBRight","CBRight",MCBMass,Mean,Sigma,RAlpha,RN);

  RooRealVar FitFraction("FitFraction","FitFraction",0.5,0.0,1.0);
  RooAddPdf DCB("DCB","DCB",RooArgList(CBRight,CBLeft),FitFraction);

  RooRealVar SignalYield("SignalYield","SignalYield",4338.0,500.0,10000.0);
  //  RooExtendPdf ExtDCB("ExtDCB","ExtDCB",DCB,SignalYield);
  
  //==============================ETA DCB ++++++++++++++++++++++++++++++
  
  RooRealVar MCEtamean("MCEtamean","MCEtamean",958.0,955.0,960.0);
  RooRealVar MCEtasigma("MCEtasigma","MCEtasigma",9.16,8.0,14.0);
  RooRealVar EtaLAlpha("EtaLAlpha","EtaLAlpha",-1.45,-5.0,1.0);
  RooRealVar EtaRAlpha("EtaRAlpha","EtaRAlpha",1.76,0.0,4.0);
  RooRealVar EtaLN("EtaLN","EtaLN",0.1,0.0,20.0);
  RooRealVar EtaRN("EtaRN","EtaRN",0.1,0.0,20.0);

  RooCBShape EtaCBLeft("EtaCBLeft","EtaCBLeft",MCEtaMass,MCEtamean,MCEtasigma,EtaLAlpha,EtaLN);
  

//.........这里部分代码省略.........
    u32 address, vaddress, size, offset, alignment;
    char name[512], container[512];
    if (four_columns)
    {
      // sometimes there is no alignment value, and sometimes it is because it is an entry of
      // something else
      if (length > 37 && line[37] == ' ')
      {
        alignment = 0;
        sscanf(line, "%08x %08x %08x %08x %511s", &address, &size, &vaddress, &offset, name);
        char* s = strstr(line, "(entry of ");
        if (s)
        {
          sscanf(s + 10, "%511s", container);
          char* s2 = (strchr(container, ')'));
          if (s2 && container[0] != '.')
          {
            s2[0] = '\0';
            strcat(container, "::");
            strcat(container, name);
            strcpy(name, container);
          }
        }
      }
      else
      {
        sscanf(line, "%08x %08x %08x %08x %i %511s", &address, &size, &vaddress, &offset,
               &alignment, name);
      }
    }
    // some entries in the table have a function name followed by " (entry of " followed by a
    // container name, followed by ")"
    // instead of a space followed by a number followed by a space followed by a name
    else if (length > 27 && line[27] != ' ' && strstr(line, "(entry of "))
    {
      alignment = 0;
      sscanf(line, "%08x %08x %08x %511s", &address, &size, &vaddress, name);
      char* s = strstr(line, "(entry of ");
      if (s)
      {
        sscanf(s + 10, "%511s", container);
        char* s2 = (strchr(container, ')'));
        if (s2 && container[0] != '.')
        {
          s2[0] = '\0';
          strcat(container, "::");
          strcat(container, name);
          strcpy(name, container);
        }
      }
    }
    else
    {
      sscanf(line, "%08x %08x %08x %i %511s", &address, &size, &vaddress, &alignment, name);
    }

    const char* namepos = strstr(line, name);
    if (namepos != nullptr)  // would be odd if not :P
      strcpy(name, namepos);
    name[strlen(name) - 1] = 0;
    if (name[strlen(name) - 1] == '\r')
      name[strlen(name) - 1] = 0;

    // Check if this is a valid entry.
    if (strlen(name) > 0)
    {
      // Can't compute the checksum if not in RAM
      bool good = !bad && PowerPC::HostIsInstructionRAMAddress(vaddress) &&
                  PowerPC::HostIsInstructionRAMAddress(vaddress + size - 4);
      if (!good)
      {
        // check for BLR before function
        PowerPC::TryReadInstResult read_result = PowerPC::TryReadInstruction(vaddress - 4);
        if (read_result.valid && read_result.hex == 0x4e800020)
        {
          // check for BLR at end of function
          read_result = PowerPC::TryReadInstruction(vaddress + size - 4);
          good = read_result.valid && read_result.hex == 0x4e800020;
        }
      }
      if (good)
      {
        ++good_count;
        if (section_name == ".text" || section_name == ".init")
          AddKnownSymbol(vaddress, size, name, Symbol::Type::Function);
        else
          AddKnownSymbol(vaddress, size, name, Symbol::Type::Data);
      }
      else
      {
        ++bad_count;
      }
    }
  }

  Index();
  if (bad)
    SuccessAlertT("Loaded %d good functions, ignored %d bad functions.", good_count, bad_count);
  return true;
}

LUA_API void lua_rawset (lua_State *L, int index) {
  StkId t = Index(L, index);
  LUA_ASSERT(ttype(t) == LUA_TTABLE, "table expected");
  *luaH_set(L, hvalue(t), L->top-2) = *(L->top-1);
  L->top -= 2;
}

LUA_API void lua_rawseti (lua_State *L, int index, int n) {
  StkId o = Index(L, index);
  LUA_ASSERT(ttype(o) == LUA_TTABLE, "table expected");
  *luaH_setint(L, hvalue(o), n) = *(L->top-1);
  L->top--;
}

TObject *luaA_index (lua_State *L, int index) {
  return Index(L, index);
}

void turn_by_turn(struct each_path *shortest_ways,int ways)//"ways" save the number of the shorest ways,shortest_ways[] save all the shortest pathes
{
      int i;
      printf("/---------------------Next Path-------------------/\n");  
      for (i=1;i<=ways;i++)
      {
          path_link present,next,future;
          present=&shortest_ways[i];
          next=present->next_node;
          future=next->next_node;
          next->section_distance=next->distance;
          
          float x1,yl,x2,y2,x3,y3;
          if (present->unique_num==1)//if the presetn node is start node,take the x and y value in special ways
          {
                x1=start_x;
                yl=start_y;
          }
          else if (present->unique_num!=1)
          {
                x1=present->x;
                yl=present->y;
          }
          if (next->unique_num==2)//if the next node is destination node,take the x and y value in special ways
          {
                x2=destination_x;
                y2=destination_y;
          }
          else if (next->unique_num!=2)
          {
                x2=next->x;
                y2=next->y;
          }
          if (x1>x2 && yl>y2)//diceide the first direction
                printf("Head SouthWest ");
          else if (x1<x2 && yl>y2)
                printf("Head SouthEast ");
          else if (x1>x2 && yl<y2)
                printf("Head NorthWest ");
          else if (x1<x2 && yl<y2)
                printf("Head NorthEest ");
          else if (x1>x2 && yl==y2)
                printf("Head West ");
          else if (x1<x2 && yl==y2)
                printf("Head Eest ");
          else if (x1==x2 && yl<y2)
                printf("Head North ");
          else if (x1==x2 && yl>y2)
                printf("Head South ");
                
          int present_label,future_label,next_label;
          int j,k;
          char start_street[20];
          present_label=Index(present->unique_num);
          next_label=Index(next->unique_num);
          for (j=1;j<=adj_list[present_label].present_num;j++)//find the start_street name
          {
                int get_name=0;
                for (k=1;k<=adj_list[next_label].present_num;k++)
                if (strcmp(adj_list[present_label].street_name[j],adj_list[next_label].street_name[k])==0)
                {
                      get_name=1;
                      strcpy(start_street,adj_list[present_label].street_name[j]);
                      break;
                }     
                if (get_name==1)
                break;
          }
          printf("on %s.\n",start_street);//start_street must be the same street for both first node and the second node
          if (future==NULL)//if  future==NULL means they are very closed.
          {
                float distance,miles,decimal;//to find out the miles in 1 decimal place and the foot in integer
                int integer,feet;
                
                distance=next->section_distance;
                distance=distance*10;
                integer=distance/1;
                decimal=distance-integer;
                miles=(float)integer/10;   
                feet=5280*decimal/10;
                if (miles==0)
                     printf("stop on street %s,house number is %d,in %d foot.\n",destination_street,destination_house,feet);
                else if (feet==0)
                     printf("stop on street %s,house number is %d,in %.1f miles.\n",destination_street,destination_house,miles);
                else if (miles>0 && feet>0)
                     printf("stop on street %s,house number is %d,in %.1f miles and %d foot.\n",destination_street,destination_house,miles,feet);
          }
          while (future!=NULL)
          {
                present_label=Index(present->unique_num);
                future_label=Index(future->unique_num);
                next_label=Index(next->unique_num);
                
                int same_street=0;//to check whether the future node has the same as the previous two nodes
                char same_street_name[20];
                for (j=1;j<=adj_list[present_label].present_num;j++)
                {
                      int get_name=0;
                      for (k=1;k<=adj_list[next_label].present_num;k++)
                          if (strcmp(adj_list[present_label].street_name[j],adj_list[next_label].street_name[k])==0)
//.........这里部分代码省略.........

// This one can load both leftover map files on game discs (like Zelda), and mapfiles
// produced by SaveSymbolMap below.
bool PPCSymbolDB::LoadMap(const std::string& filename)
{
	File::IOFile f(filename, "r");
	if (!f)
		return false;

	bool started = false;

	char line[512];
	while (fgets(line, 512, f.GetHandle()))
	{
		if (strlen(line) < 4)
			continue;

		char temp[256];
		sscanf(line, "%255s", temp);

		if (strcmp(temp, "UNUSED")==0) continue;
		if (strcmp(temp, ".text")==0)  {started = true; continue;};
		if (strcmp(temp, ".init")==0)  {started = true; continue;};
		if (strcmp(temp, "Starting")==0) continue;
		if (strcmp(temp, "extab")==0) continue;
		if (strcmp(temp, ".ctors")==0) break; //uh?
		if (strcmp(temp, ".dtors")==0) break;
		if (strcmp(temp, ".rodata")==0) continue;
		if (strcmp(temp, ".data")==0) co