Attachment 'KITutil.cc'

Download

   1 #include "KITutil_LDC00.h"
   2 
   3 
   4 // Will be called by KIT.cc
   5 // requires inout from myPGDB & MarlinMath
   6 
   7 //Shits = superHits
   8 void CreateAllShits2(LCCollection* colt,CellIDDecoder<CalorimeterHit>& id,vector<Superhit2*>* calo, int debug)
   9 {
  10   unsigned int n_elements = colt->getNumberOfElements();
  11   int superHitCounterSampling1 = 0;
  12   int superHitCounterSampling2 = 0;
  13 
  14   for ( unsigned int elementCounter=0; elementCounter<n_elements; elementCounter++) {     
  15     CalorimeterHit* caloHit = dynamic_cast<CalorimeterHit*>(colt->getElementAt(elementCounter)); 
  16     unsigned int layer = (unsigned int)id(caloHit)["K-1"];
  17     Superhit2 *superHit;
  18 
  19     //requires information about detector geometry from Physics Geometry database myPGDB
  20     //for 2 different ECAL samplings in barrel ECAL1 & ECAL2
  21     if ( layer<=myPGDB[myPGdb::ECAL1_BAR].max_lay ) {//maxlayer from gear-file
  22       superHit = new Superhit2(myPGDB[myPGdb::ECAL1_BAR].mip_whole,caloHit); //mip_whole is hardcoded via e_coeff in myPGDB == energy of a mip?
  23       superHitCounterSampling1++;
  24     } else {
  25       superHit = new Superhit2(myPGDB[myPGdb::ECAL2_BAR].mip_whole,caloHit);
  26       superHitCounterSampling2++;
  27     }
  28     calo[0].push_back(superHit);
  29   }
  30   if (debug > 3 ) {
  31     std::cout<<"DEBUG KITutil:  *  from "<< n_elements <<" calorimeter hits created "
  32 	     << superHitCounterSampling1 <<" superhits in sampling 1 and "
  33 	     << superHitCounterSampling2 <<" superhits in sampling 2"<<std::endl;
  34   }
  35 }
  36 
  37 //precalculations does nearest neighbour clustering 
  38 //by shifting cells to equal distances in x, y & z
  39 //called by KIT
  40 void TotalPrecalc2(vector<Superhit2*>* calo,CellIDDecoder<CalorimeterHit>& id,unsigned int nelem,int Ncut,int debug)
  41 {
  42   if ( debug > 3 )
  43     std::cout<<"DEBUG KITutil: - precalculation for "<<nelem<<" ECAL hits "<<std::endl;
  44 
  45   vector<Superhit2*> tmpsh[9]; //fields 0 to 7 are superhits in stave 0 to 7, field 8 holds the endcaps
  46   
  47   int barrelCounter = 0;
  48   int endcapCounter = 0;
  49 
  50   for ( unsigned int hitCounter=0; hitCounter<nelem; hitCounter++ ) {
  51     
  52     CalorimeterHit* chh=(calo[0])[hitCounter]->chit;
  53     int stave  = id(chh)["S-1"]; //which octant
  54     int module = id(chh)["M"];   //which module
  55     int layer  = id(chh)["K-1"]; //which layer
  56     
  57     // some remaining superhit members initialised
  58     (calo[0])[hitCounter]->S=stave;
  59     (calo[0])[hitCounter]->M=module;
  60     (calo[0])[hitCounter]->K=layer;
  61     (calo[0])[hitCounter]->isCalorimeter=1; //is ecal
  62     
  63     if ( module !=0 && module!=6 ) { //if hit not in endcaps
  64       if ( layer!=0 ) {	    //and not in first layer
  65 	tmpsh[stave].push_back((calo[0])[hitCounter]);	      
  66       } else {
  67 	
  68 	//for hits in first layer, phi = angle between hit and y axis
  69 	//the following section determines if the first layer hit is a neighbour of the next stave. 
  70 	//If so it is assigned to both stave vectors
  71 
  72 	double phi = atan2(-(calo[0])[hitCounter]->chit->getPosition()[0],(calo[0])[hitCounter]->chit->getPosition()[1]);
  73 	if ( phi< 0.0 ) {
  74 	  if ( stave!=0 ) phi = 2.0*M_PI + phi; //for staves 1 - 7; M_PI = pi = 3.14
  75 	  else phi = fabs(phi);           
  76 	}
  77 	
  78 	double diff = phi - M_PI_4*stave; //M_PI_4 = pi/4 = 0.79
  79 	// 45 degrees time is the angle between to corners of the octagon
  80 	// now the hit lies in stave 0 and if fabs(phi)<45 DEG, it is not overlapping, see below
  81 
  82 	//The hit is written to its original stave anyway
  83 	tmpsh[stave].push_back((calo[0])[hitCounter]);
  84 
  85 	if ( fabs(diff)<M_PI_4/2.0 ) { // hit in layer 0, but not in the interesting "overlap" area -> stored a 2nd time ?
  86 	  tmpsh[stave].push_back((calo[0])[hitCounter]);
  87 	} else { // If its in the overlap area, put it in both stove collections
  88 	  
  89 	  (calo[0])[hitCounter]->connect = true;
  90 	  
  91 	  tmpsh[stave].push_back((calo[0])[hitCounter]);
  92 	  
  93 	  int stavem1 = stave - 1;
  94 	  // store hit in additional stave
  95 	  if ( stavem1>=0 )
  96 	    tmpsh[stavem1].push_back((calo[0])[hitCounter]);
  97 	  else
  98 	    tmpsh[7].push_back((calo[0])[hitCounter]);
  99 	}
 100       } // special treatment for first layer (= layer 0)
 101       barrelCounter++;
 102 
 103     }   // hit in barrel
 104     else { //hits in endcap
 105       tmpsh[8].push_back((calo[0])[hitCounter]);
 106       endcapCounter++;
 107     }
 108   }//all elements
 109   if ( debug > 3 )
 110     std::cout<<"DEBUG KITutil:  * found "<<endcapCounter<<" hits in ECAL endcap "<<" and "<<barrelCounter<<" hits in ECAL barrel"<<std::endl;
 111 
 112   // Next section: find nearest neighbours
 113   // CAREFULL: detector geometry dependend !!!      
 114  
 115   double celx = myPGDB[myPGdb::ECAL1_BAR].cell_size; // cell_size has to be checked
 116   double distcut = 0.98*(2.0*celx)*(2.0*celx); //where does this come from?
 117   // hits in barrel
 118   for ( unsigned int staveCounter=0; staveCounter<8; staveCounter++ ) {
 119     if(tmpsh[staveCounter].size()!=0)
 120       Precalc2(tmpsh[staveCounter],myPGDB[myPGdb::ECAL1_BAR].r_inner,
 121 	       myPGDB[myPGdb::ECAL1_CAP].z_inner,
 122 	       myPGDB[myPGdb::ECAL1_BAR].cell_size,distcut,true,staveCounter,1,id,debug); 
 123   } // isCalorimeter = 1 for ecal, 2 for hcal, 0 else
 124   
 125   if ( tmpsh[8].size()!=0 ) // hits in end caps: tmpsh[8]
 126     Precalc2(tmpsh[8],myPGDB[myPGdb::ECAL1_BAR].r_inner,
 127 	     myPGDB[myPGdb::ECAL1_CAP].z_inner,
 128 	     myPGDB[myPGdb::ECAL1_BAR].cell_size,distcut,false,1,1,id,debug); 
 129   
 130   // rearrange calo[0]
 131   for ( unsigned int ii=0; ii<calo[0].size(); ii++) 
 132     calo[0][ii]->top = calo[0][ii]->neighbours.size();
 133 
 134   for ( unsigned int ii=0; ii<(calo[0]).size(); ii++) { 
 135     if(debug > 3) 
 136       std::cout<<"DEBUG KITutil: - found hit with "<<(calo[0])[ii]->top<<" neighbours"<<std::endl;
 137       
 138     if ( (calo[0])[ii]->top !=0 ) {//if hit has neighbours
 139       if ( (calo[0])[ii]->top >= 1 && (calo[0])[ii]->top < Ncut ) {
 140 	(calo[2]).push_back((calo[0])[ii]);// if hit with less than required neighbours -> no cluster 
 141       } else { //for topological cleaning 
 142 	(calo[4]).push_back((calo[0])[ii]);  // hit with required neighbours -> topological cleaning possible
 143       }
 144     } else {  // if calo[0] hit has no neighbours at all
 145       (calo[6]).push_back((calo[0])[ii]); //-> no cluster
 146     }
 147   }
 148 }
 149 
 150 
 151 
 152 void Precalc2 ( vector< Superhit2* >& shvec,double r, double z, double cell, double dist,bool isCalorimeter,int stave,int module,CellIDDecoder<CalorimeterHit>& id,int debug ) {
 153   
 154   unsigned int nPoints = shvec.size();
 155   ANNpointArray dataPoints;
 156   ANNpoint      queryPoint;      //Annpoint: double[3]
 157   ANNidxArray   nnIndex;
 158   ANNdistArray  distance;
 159   ANNkd_tree*   kdTree;
 160   
 161   queryPoint = annAllocPt(3); // 3d point 
 162   dataPoints = annAllocPts(nPoints,3);
 163   nnIndex    = new ANNidx[36];
 164   distance   = new ANNdist[36];
 165   float xxx[3];          // transformed position in GridTransform
 166   float Ecalradius = r;
 167   float Ecalhalfz = z;
 168   float Ecalcellsize = cell;
 169   
 170   if ( debug>3 )
 171     std::cout<<"DEBUG KITutil: shifting hits to new positions"<<std::endl;
 172 
 173   for ( unsigned int hitCounter=0; hitCounter<nPoints; hitCounter++ ) {  
 174 
 175     GridTransform2(shvec[hitCounter]->chit,Ecalradius,Ecalhalfz,Ecalcellsize,xxx,isCalorimeter,stave,module,id);
 176     //dataPoints are filled with transformed coordinates for search alg. 
 177     dataPoints[hitCounter][0]   = xxx[0];  
 178     dataPoints[hitCounter][1]   = xxx[1];  
 179     dataPoints[hitCounter][2]   = xxx[2];  
 180     //point[] is the transformed position in class SuperHit
 181     shvec[hitCounter]->point[0] = xxx[0];
 182     shvec[hitCounter]->point[1] = xxx[1];
 183     shvec[hitCounter]->point[2] = xxx[2];
 184   }
 185   
 186   //search for nearest neighbours with ANN code, assumed to work correctly see ANN.cpp
 187   kdTree = new ANNkd_tree(dataPoints,nPoints,3);                                                       
 188   ANNdist trd = dist;
 189   
 190   for ( unsigned int hitCounter=0; hitCounter<nPoints; hitCounter++ ) {
 191     //number of neighbours of each hit  	      
 192     unsigned int N_NearestNeighbours = (unsigned int)kdTree->annkFRSearch(dataPoints[hitCounter],trd,26,nnIndex,distance,0.0);
 193     if ( N_NearestNeighbours>26 ) 
 194       N_NearestNeighbours = 26;
 195     if( isCalorimeter ) { // inside calorimeter: isCalorimeter(ecal)=1, isCalorimeter(hcal)=2 
 196       if( shvec[hitCounter]->connect ) { // hit in barrel and connected to two staves
 197 	for ( unsigned int neighbourCounter=0; neighbourCounter<N_NearestNeighbours; neighbourCounter++ ) {
 198 	  if ( shvec[hitCounter]->neighbours.end() 
 199 	       == find(shvec[hitCounter]->neighbours.begin(),shvec[hitCounter]->neighbours.end(), shvec[nnIndex[neighbourCounter]]) ) {
 200 	    //add neighbour to SuperHit
 201 	    shvec[hitCounter]->neighbours.push_back(shvec[nnIndex[neighbourCounter]]);
 202 	    shvec[hitCounter]->top++;
 203 	    
 204 	  }
 205 	  if ( shvec[nnIndex[neighbourCounter]]->neighbours.end() 
 206 	       == find(shvec[nnIndex[neighbourCounter]]->neighbours.begin(),shvec[nnIndex[neighbourCounter]]->neighbours.end(), shvec[hitCounter]) ) {		
 207 	    shvec[nnIndex[neighbourCounter]]->neighbours.push_back(shvec[hitCounter]);
 208 	    shvec[nnIndex[neighbourCounter]]->top++;
 209 	  }
 210 	}//all neighbours
 211 	
 212       } else if ( N_NearestNeighbours !=0 ) {
 213 	
 214 	shvec[hitCounter]->top = N_NearestNeighbours;		   
 215 	for ( unsigned int neighbourCounter=0; neighbourCounter<N_NearestNeighbours ; neighbourCounter++ )
 216 	  shvec[hitCounter]->neighbours.push_back( shvec[nnIndex[neighbourCounter]]);
 217       }
 218     } else if ( N_NearestNeighbours !=0 ) {
 219       
 220       shvec[hitCounter]->top=N_NearestNeighbours;
 221       for ( unsigned int neighbourCounter=0; neighbourCounter<N_NearestNeighbours; neighbourCounter++ )
 222 	shvec[hitCounter]->neighbours.push_back( shvec[nnIndex[neighbourCounter]]);
 223     }	        
 224   }//all hits
 225   
 226   delete [] nnIndex;			     
 227   delete [] distance;
 228   annDeallocPts(dataPoints);
 229   delete kdTree;
 230   annClose();   
 231 }
 232 
 233 
 234 //called by precalc2
 235 void GridTransform2 ( CalorimeterHit* caloHit,float& radius, float& halfz, float& cellsize,float*X,bool isCalorimeter,int stave,int module,CellIDDecoder<CalorimeterHit>& id ) {
 236   unsigned int layer;
 237   double coss;
 238   double sinn;
 239   
 240   if ( !isCalorimeter ) {	// !calorimeter -> tracker??
 241     stave  = id(caloHit)["S-1"];
 242     module = id(caloHit)["M"]; //module is back to 6 or 1 for endcap? 
 243     layer  = id(caloHit)["K-1"];
 244     
 245   } else {
 246     
 247     int stave2 = id(caloHit)["S-1"];               
 248     // The following calcs are done for the secondary stave the hit is assigned to, where it floats "freely"in space until now.
 249     if ( stave2!=stave ) {
 250       //define rotation matrix, hits are rotated to stave 0
 251       coss = cos(M_PI_4*stave);
 252       sinn = sin(M_PI_4*stave);
 253       //rotate stave to position 0 (only y-axis)
 254       float tmp_y = coss*caloHit->getPosition()[1] - sinn*caloHit->getPosition()[0];
 255       //divide distance of hit from ecal face by sampling layer thickness->layer
 256       layer = (unsigned int) ((tmp_y-radius)/myPGDB[myPGdb::ECAL1_BAR].sampling); // largest int number: 2.3->2, therefore n_sampl-1 in next line
 257 
 258       if ( layer> myPGDB[myPGdb::ECAL1_BAR].n_sampl - 1 ) { // last layer of first ECAl structure
 259 	layer = myPGDB[myPGdb::ECAL1_BAR].n_sampl - 1
 260 	  + (unsigned int)(( tmp_y-radius- myPGDB[myPGdb::ECAL1_BAR].n_sampl* 
 261 			     myPGDB[myPGdb::ECAL1_BAR].sampling)/myPGDB[myPGdb::ECAL2_BAR].sampling); 
 262       }
 263       
 264     } else { // (stave2 == stave) original layer for hits in their original stave 
 265       layer = id(caloHit)["K-1"]; 
 266     }
 267     //Now we know of every hit's (new) layer
 268   }
 269   // out of if 
 270   coss = cos(M_PI_4*stave);
 271   sinn = sin(M_PI_4*stave);
 272   
 273   if ( module !=0 && module!=6 ) {
 274     float tmp_x = coss*caloHit->getPosition()[0] + sinn*caloHit->getPosition()[1]; // rotate, find new y position
 275     float tmp_y = radius + (layer+1)*cellsize; // center of layer position. 
 276     // y using the variable cellsize, the hits are equally spaced in 3d
 277     
 278     //rotate back
 279     X[0] = coss*tmp_x - sinn*tmp_y;
 280     X[1] = coss*tmp_y + sinn*tmp_x;
 281     X[2] = caloHit->getPosition()[2];
 282 		 
 283   } else { // endcap
 284     X[0] = caloHit->getPosition()[0];
 285     X[1] = caloHit->getPosition()[1];
 286     X[2] = (-1+module/3)*(halfz+layer*cellsize);
 287   }  
 288 }
 289 
 290 
 291 // searching for cluster cores
 292 // wil be called by KIT 
 293 // secal1 = superhits in calo[4] with topological cleaning
 294 // secal1 = superhits in calo[0] without topological cleaning
 295 void FindCores2 (Shitvec2* secal1, Tmpclvec2* clusterVector , vector <PROTSEED2> * photonSeed2, unsigned int Nlevels, vector<float> mipstep, CoreCut Ccut,int debug )
 296 { 
 297   if (debug >2 )
 298     std::cout<<" ***** started main loop: searching for cluster cores  in "<<secal1->size()<<" hits ***"<<std::endl;
 299 
 300   typedef pair<int,int> test;
 301   int psl_plun_global=0;
 302   
 303   double Diagby2       = myPGDB[myPGdb::ECAL1_BAR].cell_size*sqrt(2.0)/2.0;
 304   double barrel_rInner = myPGDB[myPGdb::ECAL1_BAR].r_inner; //inner ECAL radius
 305   double barrel_zMax   = myPGDB[myPGdb::ECAL1_CAP].z_inner; //maximum z barrel 
 306   double Xatf[3];
 307   vector<Superhit2*> mySuperhits[Nlevels]; //lcio does not allow for more than 16 levels
 308   
 309   for ( unsigned int hitCounter = 0; hitCounter<secal1->size(); hitCounter++ ) {	      
 310     int myLevel = 0;
 311 
 312     //search for cluster level thresholds: energy per hit > mipstep[level]
 313     for ( unsigned int levelCounter=0; levelCounter<Nlevels; levelCounter++ ) {
 314       if ( (*secal1)[hitCounter]->mip > mipstep[levelCounter] ) 
 315 	myLevel = levelCounter;
 316       else 
 317 	break;
 318     }  
 319     //write back hits into correct levels
 320     for ( int levelCounter=myLevel; levelCounter>=0; levelCounter-- ) {
 321       mySuperhits[levelCounter].push_back((*secal1)[hitCounter]);
 322       if ( debug >3 )
 323 	std::cout<<"DEBUG KITutil:  * writing hits of level "<<levelCounter<<" to sorted superhits vector"<<std::endl;
 324     }
 325   }//all super hits
 326   //temporary cluster vector
 327   vector<Tmpcl2*>  testvector;
 328   //assigning neighbouring hits from superHitVector to temporary cluster for each level
 329   for ( unsigned int levelCounter=0; levelCounter<Nlevels; levelCounter++ ) { 
 330     cluster5(&mySuperhits[levelCounter],&clusterVector[levelCounter],debug);  
 331   }
 332 
 333   for ( unsigned int levelCounter=0; levelCounter<Nlevels; levelCounter++ ) {
 334     
 335     if ( clusterVector[levelCounter].size()!=0 ) {
 336       for ( unsigned int clusterCounter=0; clusterCounter<clusterVector[levelCounter].size(); clusterCounter++ ) {
 337 	clusterVector[levelCounter][clusterCounter]->calcCenter(); //energy weighted cluster center
 338 	if (debug >4 )
 339 	  std::cout<<"calculated cluster center for level "<<levelCounter<<std::endl;
 340 	clusterVector[levelCounter][clusterCounter]->findInertia();//cluster start
 341 	if (debug >4 )
 342 	  std::cout<<"calculated cluster start for level "<<levelCounter<<std::endl;
 343 	
 344       }
 345     }
 346   }//all levels
 347 	
 348   int myLevel = Nlevels - 1;  
 349   for ( unsigned int levelCounter=0; levelCounter<Nlevels; levelCounter++ ) {
 350     if ( clusterVector[levelCounter].size()==0 ) {
 351       myLevel = levelCounter - 1;
 352       break; 
 353     }//no hits in this cluster level go one level down
 354     myLevel = levelCounter;
 355   }
 356   psl_plun_global = myLevel;
 357   
 358   vector < PROTSEED2 > photonSeed;
 359 
 360   for (unsigned int clusterCounter=0; clusterCounter<clusterVector[0].size(); clusterCounter++ ) {      
 361     if ( clusterVector[0][clusterCounter]->hits.size()>Ccut.MinHit0 ) {
 362       //yet another protseed
 363       PROTSEED2 km;
 364       km.cl     = clusterVector[0][clusterCounter];
 365       km.level  = 0;
 366       km.X[0]   = 0.0;
 367       km.X[1]   = 0.0;
 368       km.X[2]   = 0.0;
 369       km.active = true;
 370       photonSeed.push_back(km);
 371     }
 372   }
 373 	       
 374 	    
 375   for ( unsigned int seedCounter=0; seedCounter<photonSeed.size(); seedCounter++ ) { 		 
 376 
 377     Xatf[0] = 0.0;
 378     Xatf[1] = 0.0;
 379     Xatf[2] = 0.0;
 380     
 381     //search for cluster closest to IP
 382     LineCaloIntersectD2(photonSeed[seedCounter].cl->getCenter(), photonSeed[seedCounter].cl->direction,barrel_rInner,barrel_zMax,Xatf,debug);//cluster center, cluster direction, inner barrel radius, inner z barrel, 
 383     
 384     if (debug >4 ) {
 385       std::cout<<"calculated cluster - IP intersect for seed "<<seedCounter
 386 	       <<" of "<<photonSeed.size()<<std::endl;
 387       std::cout<<"- cluster center: "<<photonSeed[seedCounter].cl->getCenter()<<std::endl;
 388       std::cout<<"- cluster direction: "<<photonSeed[seedCounter].cl->direction<<std::endl;
 389       std::cout<<"- ECAL1 inner barrel: "<<barrel_rInner<<std::endl;
 390       std::cout<<"- ECAL1 max z: "<<barrel_zMax<<std::endl;
 391     }
 392 
 393     //assign this cluster shifted by its depth in ECAL to photonSeed
 394     photonSeed[seedCounter].X[0] = Xatf[0];
 395     photonSeed[seedCounter].X[1] = Xatf[1];
 396     photonSeed[seedCounter].X[2] = Xatf[2];
 397     photonSeed2->push_back(photonSeed[seedCounter]); 
 398 		 
 399     if (debug >5 )
 400       std::cout<<"assign shifted hits to new photonSeed of size"<<photonSeed2->size()<<std::endl;
 401   }//all seed entries
 402   
 403 
 404   //search for splitted clusters
 405   for ( int levelCounter=1; levelCounter<=myLevel; levelCounter++) {
 406     
 407     for ( unsigned int seed2Counter=0; seed2Counter<photonSeed2->size(); seed2Counter++) {   
 408       if ( (*photonSeed2)[seed2Counter].active ) {
 409 
 410 	Tmpclvec2 tempCluster; //clusters inside clusters
 411 	ClusterInCluster2((*photonSeed2)[seed2Counter].cl,clusterVector[levelCounter],tempCluster, debug);
 412 	
 413 	//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
 414 	//FIXME: CAREFULL photonSeed.Size increases 1 for each call of ClusterInCluster2 ! -> Infinite loop !
 415 	//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
 416 	if (debug >4 )
 417 	  std::cout<<"after ClusterInCluster2:  photonSeed size = "<<photonSeed2->size()<<std::endl;
 418 	
 419 	if ( tempCluster.size()<=1 )  {
 420 	  if (debug >4 )
 421 	    std::cout<<"did not find cluster in cluster "<<std::endl;
 422 	  // 1-in-1 nothing to do 
 423 	} else { //found cluster tempCluster in cluster photonSeed2
 424 	         // Split the split in two categories 
 425 		 // clear one (larger distance) and not so clear
 426 	  
 427 	  if (debug >4 )
 428 	    std::cout<<"level "<<levelCounter<<": found cluster in cluster "<<std::endl;
 429 
 430 	  vector < int > largerDistanceCluster;
 431 	  vector < int > smallerDistanceCluster;
 432 	  //check included clusters
 433 	  for ( unsigned int clusterCounter1=0; clusterCounter1<tempCluster.size(); clusterCounter1++ ) { 
 434 	    for ( unsigned int clusterCounter2=0; clusterCounter2<tempCluster.size(); clusterCounter2++ ) {
 435 	      if ( clusterCounter1!=clusterCounter2 ) {
 436 		if ( tempCluster[clusterCounter1]->hits.size()>=Ccut.MinHitSplit //more than minimum numgber of hits in included cluster
 437 		     && tempCluster[clusterCounter2]->hits.size()>=Ccut.MinHitSplit ) {
 438 		  
 439 		  double XIP[3];
 440 		  XIP[0] = 0.0;  
 441 		  XIP[1] = 0.0;  
 442 		  XIP[2] = 0.0;
 443 		  
 444 		  
 445 		  if ( (LinePointDistance2(XIP,(*photonSeed2)[seed2Counter].cl->center,tempCluster[clusterCounter2]->center) > Diagby2 
 446 			|| LinePointDistance2(XIP,(*photonSeed2)[seed2Counter].cl->center,tempCluster[clusterCounter1]->center) > Diagby2 ) ) {//distance of tempCluster to photonSeed2 > cellSize *sqrt(2) / 2 ->clear cluster split 
 447 		    if ( find(largerDistanceCluster.begin(),largerDistanceCluster.end(),(int)clusterCounter2)==largerDistanceCluster.end() )
 448 		      
 449 		      largerDistanceCluster.push_back(clusterCounter2);
 450 
 451 		    if ( find(largerDistanceCluster.begin(),largerDistanceCluster.end(),(int)clusterCounter1)==largerDistanceCluster.end() )
 452 
 453 		      largerDistanceCluster.push_back(clusterCounter1);
 454 		  } else { //small distance between clusters -> split unclear
 455 		    if (debug >5 )
 456 		      std::cout<<"not so clearly seperated small distant cluster -> Split unclear"<<std::endl;
 457 		    
 458 		    if (find(smallerDistanceCluster.begin(),smallerDistanceCluster.end(),(int)clusterCounter2)==smallerDistanceCluster.end() )
 459 
 460 		      smallerDistanceCluster.push_back(clusterCounter2);
 461 
 462 		    if (find(smallerDistanceCluster.begin(),smallerDistanceCluster.end(),(int)clusterCounter1)==smallerDistanceCluster.end() )
 463 
 464 		      smallerDistanceCluster.push_back(clusterCounter1);
 465 		  }
 466 		}//found cluster in cluster
 467 	      }//clusterCounters don't agree
 468 	    }//all tempCluster entries, first loop
 469 	  }//all tempCluster entries, second loop
 470 	  
 471 	  if ( largerDistanceCluster.size()!=0 ) { //if clearly seperated cluster exists
 472 
 473 	    if (debug >5 )
 474 	      std::cout<<"clearly seperated large distant cluster "<<std::endl;
 475 
 476 	    if (debug >6 )
 477 	      std::cout<<"for clearly seperated large distant clusters:  photonSeed size = "<<photonSeed2->size()<<std::endl;	    
 478 	    int clusterSize1 = 0;
 479 	    int clusterSize3 = (*photonSeed2)[seed2Counter].cl->hits.size();
 480 	    for ( unsigned int vectorEntry=0; vectorEntry<largerDistanceCluster.size(); vectorEntry++ )
 481 	      clusterSize1 += tempCluster[largerDistanceCluster[vectorEntry]]->hits.size();
 482 	    
 483 	    double ratio = ((double)clusterSize1) / ((double)clusterSize3);
 484 	    (*photonSeed2)[seed2Counter].active = false;
 485 
 486 	    if ( ratio > Ccut.rCut ) { //require minimum ratio of cluster entries to supress fluctuations
 487 	      for ( unsigned int clusterEntry=0; clusterEntry<smallerDistanceCluster.size(); clusterEntry++ ) { //have a look at the not so clear splitted clusters
 488 
 489 		if (debug >5 )
 490 		  std::cout<<"cluster ratio larger than"<<Ccut.rCut<<" hits, loop over small distance clusters"<<std::endl;
 491 
 492 		Xatf[0] = 0.0;
 493 		Xatf[1] = 0.0;
 494 		Xatf[2] = 0.0;
 495 
 496 		PROTSEED2 clusterSeed;
 497 		
 498 		clusterSeed.cl       = tempCluster[smallerDistanceCluster[clusterEntry]];
 499 		clusterSeed.level    = levelCounter;
 500 		LineCaloIntersectD2(tempCluster[smallerDistanceCluster[clusterEntry]]->getCenter(),tempCluster[smallerDistanceCluster[clusterEntry]]->direction,barrel_rInner,barrel_zMax,Xatf,debug);
 501 		clusterSeed.X[0]     = Xatf[0];
 502 		clusterSeed.X[1]     = Xatf[1];
 503 		clusterSeed.X[2]     = Xatf[2];
 504 		clusterSeed.active   = true;
 505 		bool split = true;
 506 	
 507 		for ( unsigned int seed2Counter=0; seed2Counter<photonSeed2->size(); seed2Counter++ ) {
 508 		  if ( (*photonSeed2)[seed2Counter].active ) {
 509 		    if ( (*photonSeed2)[seed2Counter].cl==clusterSeed.cl )
 510 		      split=false;
 511 		    if (debug >5 )
 512 		      std::cout<<"do not split cluster"<<std::endl;
 513 		  }
 514 		}
 515 		if (split) {
 516 		  photonSeed2->push_back(clusterSeed);
 517 		  if (debug >5 )
 518 		    std::cout<<"splitting cluster"<<std::endl;
 519 		}
 520 	      }
 521 	    }//ratio > ratio cut
 522 		      
 523 	    
 524 	    for ( unsigned int clusterEntry=0; clusterEntry<largerDistanceCluster.size(); clusterEntry++ ) { 
 525 			  
 526 	      Xatf[0] = 0.0;
 527 	      Xatf[1] = 0.0;
 528 	      Xatf[2] = 0.0;
 529 	      
 530 	      PROTSEED2 clusterSeed;
 531 			 
 532 	      if (debug >4 )
 533 		std::cout<<"cluster larger than "<<Ccut.rCut<<" hits, loop over large distance clusters"<<std::endl;
 534 
 535 	      clusterSeed.cl    = tempCluster[largerDistanceCluster[clusterEntry]];
 536 	      clusterSeed.level = myLevel;
 537 	      LineCaloIntersectD2(tempCluster[largerDistanceCluster[clusterEntry]]->getCenter(),tempCluster[largerDistanceCluster[clusterEntry]]->direction,barrel_rInner,barrel_zMax,Xatf,debug);
 538 	      clusterSeed.X[0]   = Xatf[0];
 539 	      clusterSeed.X[1]   = Xatf[1];
 540 	      clusterSeed.X[2]   = Xatf[2];
 541 	      clusterSeed.active = true;
 542 	      bool split = true;
 543 	      for ( unsigned int seed2Counter=0; seed2Counter<photonSeed2->size(); seed2Counter++ ) {
 544 		if ( (*photonSeed2)[seed2Counter].active ) {
 545 		  if ( (*photonSeed2)[seed2Counter].cl==clusterSeed.cl )
 546 		    split = false;
 547 		}
 548 	      }
 549 	      if ( split )
 550 		photonSeed2->push_back(clusterSeed);
 551 	    }
 552 	  }//found split clusterSeed at large distance
 553 	  if ( debug >6 ) std::cout<<"exit large distance cluster"<<std::endl;
 554 	}//found cluster in cluster
 555 	if ( debug >6 ) std::cout<<"exit cluster in cluster"<<std::endl;
 556       }// if active 
 557       if ( debug >6 ) std::cout<<"exit photonseed.active"<<std::endl;
 558       if ( debug >6 ) std::cout<<"for seed "<<seed2Counter<<" of "<< photonSeed2->size()<<std::endl;
 559     }//in photonSeed2
 560     if ( debug >6 ) std::cout<<"exit loop over photonseed"<<std::endl;
 561   }//all hit energy threshold levels
 562   if ( debug >6 ) std::cout<<"exit loop over threhold limits"<<std::endl;
 563 
 564   vector <test> mergeVector;
 565   
 566   for ( unsigned int hitCounter=0; hitCounter<photonSeed2->size(); hitCounter++ ) {
 567     if( (*photonSeed2)[hitCounter].active ) {		      
 568 	
 569       double Xa[3];
 570       ModuleNormal2((*photonSeed2)[hitCounter].cl->center,barrel_zMax,Xa); //shift hits in barrel
 571       if (debug >3 )
 572 	std::cout<<" start merging: shifting barrel hits"<<std::endl;
 573       
 574       //double loop over all hits?
 575       for ( unsigned int hitCounter2=0; hitCounter2<photonSeed2->size(); hitCounter2++) {
 576 
 577 	//if hitcounters do not agree, get difference between these hits
 578 	if ( hitCounter2!=hitCounter && (*photonSeed2)[hitCounter2].active ) {
 579 	  
 580 	  double dir_diff[3];
 581 	  dir_diff[0] = (*photonSeed2)[hitCounter].cl->center[0]
 582 	    - (*photonSeed2)[hitCounter2].cl->center[0];
 583 	  dir_diff[1] = (*photonSeed2)[hitCounter].cl->center[1]
 584 	    - (*photonSeed2)[hitCounter2].cl->center[1];
 585 	  dir_diff[2] = (*photonSeed2)[hitCounter].cl->center[2]
 586 	    - (*photonSeed2)[hitCounter2].cl->center[2];    
 587 			    
 588 	  //MERGE CONDITION for clusters:
 589 	  //D_cl_cl2 returns smallest distance between cluster1 & cluster2
 590 	  //Dot2 returns angle between clusters center & split candidate ??
 591 	  if ( D_cl_cl2((*photonSeed2)[hitCounter].cl,(*photonSeed2)[hitCounter2].cl) < Ccut.distCut //minimum distance between clusters < distance cut
 592 	       && fabs( Dot2((*photonSeed2)[hitCounter2].cl->center,dir_diff)) > Ccut.cosCut) { //minimum angle between clusters
 593 	    
 594 	    test cluster1,cluster2; //pair of cluster1 & cluster2
 595 	    cluster1.first  = hitCounter;
 596 	    cluster1.second = hitCounter2;
 597 	    cluster2.first  = hitCounter2;
 598 	    cluster2.second = hitCounter;
 599 	    //write result in test vector
 600 	    if ( find(mergeVector.begin(),mergeVector.end(),cluster1)==mergeVector.end() 
 601 		 && find(mergeVector.begin(),mergeVector.end(),cluster2)==mergeVector.end() ) {
 602 	      mergeVector.push_back(cluster1);
 603 	      
 604 	      if (debug >3 )
 605 		std::cout<<" merging clusters"<<std::endl;
 606 	    }
 607 	  }//if merge condition
 608 	}//hitCounters do not agree
 609       }//all photonSeed entries
 610     }//photonSeed active
 611   }//all photonSeed entries
 612 
 613 
 614   //	  if(testVector.size()!=0) //if testvector not empty -> differences in hits
 615   if ( mergeVector.size()!=0) {
 616            
 617     int benefit[mergeVector.size()];
 618     for ( unsigned int vectorCounter=0; vectorCounter<mergeVector.size(); vectorCounter++) {
 619       benefit[vectorCounter] = 0;
 620     }
 621     for ( unsigned int vectorCounter=0; vectorCounter<mergeVector.size(); vectorCounter++ ) {
 622       if ( benefit[vectorCounter]==0 ) { 
 623 
 624 	vector <int> trinityVector; //holds entries from merged clusters
 625 	trinityVector.push_back(mergeVector[vectorCounter].first);
 626 	trinityVector.push_back(mergeVector[vectorCounter].second);
 627 	benefit[vectorCounter] = 1;
 628 	for(unsigned int vectorCounter2=0; vectorCounter2<mergeVector.size(); vectorCounter2++ ) {
 629 		     
 630 	  if ( (find(trinityVector.begin(),trinityVector.end(),mergeVector[vectorCounter2].first) != trinityVector.end() 
 631 	      && find(trinityVector.begin(),trinityVector.end(),mergeVector[vectorCounter2].second) == trinityVector.end() ) 
 632 	       || (find(trinityVector.begin(),trinityVector.end(),mergeVector[vectorCounter2].second) != trinityVector.end() 
 633 		   && find(trinityVector.begin(),trinityVector.end(),mergeVector[vectorCounter2].first) == trinityVector.end() ) ) {
 634 	      
 635 	    if ( find(trinityVector.begin(),trinityVector.end(),mergeVector[vectorCounter2].first) == trinityVector.end() ) {
 636 	      if ( benefit[vectorCounter2] !=1 ) { //->vectorCounter1 != vectorCounter
 637 				 
 638 		benefit[vectorCounter2] = 1;
 639 		if ( (*photonSeed2)[mergeVector[vectorCounter2].first].active )
 640 		  trinityVector.push_back(mergeVector[vectorCounter2].first);
 641 	      }
 642 	      
 643 	    } else if ( benefit[vectorCounter2] != 1 ) { 
 644 	      if( (*photonSeed2)[mergeVector[vectorCounter2].second].active)
 645 		trinityVector.push_back(mergeVector[vectorCounter2].second);
 646 	      benefit[vectorCounter2] = 1;
 647 	    }  
 648 	  } 
 649 	}
 650 	
 651 		 
 652 	int levelCheck = 1;
 653 	for ( unsigned int vectorEntry=1; vectorEntry<trinityVector.size(); vectorEntry++) {
 654 	  if( (*photonSeed2)[trinityVector[0]].level != (*photonSeed2)[trinityVector[vectorEntry]].level )
 655 	    levelCheck++;
 656 	}
 657 	
 658 	if ( levelCheck==1 ) { //all cluster entries in same level
 659 	  
 660 	  Tmpcl2* cl = new Tmpcl2();
 661 	  for ( unsigned int vectorEntry=0; vectorEntry<trinityVector.size(); vectorEntry++ ) {
 662 	    if ( (*photonSeed2)[trinityVector[vectorEntry]].active ) {
 663 	      for ( unsigned int hitCounter=0; hitCounter<(*photonSeed2)[trinityVector[vectorEntry]].cl->hits.size(); hitCounter++ ) {	
 664 		cl->hits.push_back((*photonSeed2)[trinityVector[vectorEntry]].cl->hits[hitCounter]);
 665 	      }
 666 	      (*photonSeed2)[trinityVector[vectorEntry]].active = false;
 667 	    }
 668 	  }
 669 	  if ( cl->hits.size()!=0 ) {//found merged cluster
 670 	    
 671 	    cl->calcCenter();
 672 	    if (debug >4 )
 673 	      std::cout<<" calculate center of merged cluster"<<std::endl;
 674 	    cl->findInertia();
 675 	    if (debug >4 )
 676 	      std::cout<<" calculate start of merged cluster"<<std::endl;
 677 	    clusterVector[(*photonSeed2)[trinityVector[0]].level].push_back(cl);
 678 	    
 679 	    Xatf[0] = 0.0;
 680 	    Xatf[1] = 0.0;
 681 	    Xatf[2] = 0.0;
 682 		
 683 	    PROTSEED2 clusterSeed; 
 684 	    clusterSeed.cl     = cl;
 685 	    clusterSeed.level  = (*photonSeed2)[trinityVector[0]].level;
 686 	    LineCaloIntersectD2(cl->getCenter(),cl->direction,barrel_rInner,barrel_zMax,Xatf,debug);
 687 	    clusterSeed.X[0]   = Xatf[0];
 688 	    clusterSeed.X[1]   = Xatf[1];
 689 	    clusterSeed.X[2]   = Xatf[2];
 690 	    clusterSeed.active = true;
 691 		
 692 	    bool split = true;
 693 	    for ( unsigned int seedEntry=0; seedEntry<photonSeed2->size(); seedEntry++ ) {
 694 	      if ( (*photonSeed2)[seedEntry].active ) {
 695 		if( (*photonSeed2)[seedEntry].cl==clusterSeed.cl)
 696 		  split = false;
 697 	      }
 698 	    }
 699 	    if(split)
 700 	      photonSeed2->push_back(clusterSeed);
 701 	    
 702 	  } else {//no merged clusters
 703 	    delete cl;
 704 	  }
 705 	} else {//levels do not agree
 706 	  
 707 	  if (debug >5 )
 708 	    std::cout<<" clusters to merge have different levels -> create new vector"<<std::endl;
 709 	  vector <int > newVector;
 710 	  for ( unsigned int j=0; j<trinityVector.size(); j++ ) {//write trinityVector into newVector
 711 	    if ( find(newVector.begin(),newVector.end(),(*photonSeed2)[trinityVector[j]].level)==newVector.end() )
 712 	      newVector.push_back((*photonSeed2)[trinityVector[j]].level);
 713 	  }
 714 	  
 715 	  if ( trinityVector.size()==2 ) {//should always be, else trinityVector empty??
 716 	    if ( (*photonSeed2)[trinityVector[0]].cl->getEnergy() >=(*photonSeed2)[trinityVector[1]].cl->getEnergy() ) {
 717 	      (*photonSeed2)[trinityVector[1]].active = false;
 718 	    } else {
 719 	      (*photonSeed2)[trinityVector[0]].active = false;
 720 	    }
 721 	  } else {//trinityVector.size !=2
 722                           
 723 	    vector <int> selectedHits;
 724 	    int sumedUp[trinityVector.size()];
 725 	    for ( unsigned int trinityEntry=0; trinityEntry<trinityVector.size(); trinityEntry++ ) {
 726 	      sumedUp[trinityEntry] = 0;
 727 	    }
 728 	    for ( unsigned int newVectorEntry=0; newVectorEntry<newVector.size(); newVectorEntry++ ) {
 729 	      vector<int> sumUpVector;
 730 	      for ( unsigned int trinityEntry=0; trinityEntry<trinityVector.size(); trinityEntry++) {
 731 		if ( sumedUp[trinityEntry]==0 ) {
 732 		  if ( newVector[newVectorEntry]==(*photonSeed2)[trinityVector[trinityEntry]].level) {
 733 		    sumUpVector.push_back(trinityVector[trinityEntry]);
 734 		    sumedUp[trinityEntry] = 1;
 735 		  }
 736 		} 
 737 	      }
 738 
 739 	      if ( sumUpVector.size()>1 ) {
 740 		
 741 		Tmpcl2* cl = new Tmpcl2();
 742 		for ( unsigned int sumUpVectorEntry=0; sumUpVectorEntry<sumUpVector.size(); sumUpVectorEntry++ ) {
 743 		  (*photonSeed2)[sumUpVector[sumUpVectorEntry]].active = false;
 744 
 745 		  for ( unsigned int hitCounter=0; hitCounter<(*photonSeed2)[sumUpVector[sumUpVectorEntry]].cl->hits.size(); hitCounter++ )
 746 		    cl->hits.push_back((*photonSeed2)[sumUpVector[sumUpVectorEntry]].cl->hits[hitCounter]);
 747 		  if (debug >5 )
 748 		    std::cout<<" writing new, merged cluster"<<std::endl;
 749 		} 
 750 		
 751 		cl->calcCenter();
 752 		cl->findInertia();
 753 		clusterVector[newVector[newVectorEntry]].push_back(cl);
 754 		
 755 
 756 		Xatf[0] = 0.0;
 757 		Xatf[1] = 0.0;
 758 		Xatf[2] = 0.0;
 759 		
 760 		PROTSEED2 clusterSeed;
 761 		    
 762 		clusterSeed.cl     = cl;
 763 		clusterSeed.level  = (*photonSeed2)[newVector[newVectorEntry]].level;
 764 		LineCaloIntersectD2(cl->getCenter(),cl->direction,barrel_rInner,barrel_zMax,Xatf,debug);
 765 		clusterSeed.X[0]   = Xatf[0];
 766 		clusterSeed.X[1]   = Xatf[1];
 767 		clusterSeed.X[2]   = Xatf[2];
 768 		clusterSeed.active = true;
 769 			
 770 		bool split = true;
 771 		for ( unsigned int im=0; im<photonSeed2->size(); im++ ) {
 772 		  if ( (*photonSeed2)[im].active ) {
 773 		    if ( (*photonSeed2)[im].cl==clusterSeed.cl )
 774 		      split=false;
 775 		  }
 776 		}
 777 		if ( split ) {
 778 		  photonSeed2->push_back(clusterSeed);
 779 		  selectedHits.push_back(photonSeed2->size()-1);
 780 		} else {
 781 		  delete cl;
 782 		}
 783 	      } else {
 784 		selectedHits.push_back(sumUpVector[0]);
 785 	      }			   
 786 	    }//all newVector entries
 787 
 788 	    if ( selectedHits.size()>1 ) {
 789 	      if (debug >5 )
 790 		std::cout<<" pairing merged cluster"<<std::endl;
 791 	      typedef pair <double,int> DpairI;
 792 	      vector <DpairI> selectedPairsVector;
 793 	      for ( unsigned int j=0; j<selectedHits.size(); j++ ) {
 794 
 795 		DpairI tmp;
 796 		tmp.first  = (*photonSeed2)[selectedHits[j]].cl->getEnergy();			    
 797 		tmp.second = selectedHits[j];
 798 		selectedPairsVector.push_back(tmp);
 799 	      }
 800 	      
 801 	      sort(selectedPairsVector.begin(),selectedPairsVector.end());
 802 	      
 803 	      for ( unsigned int j=0; j<selectedHits.size()-1; j++ ) {
 804 		(*photonSeed2)[selectedPairsVector[j].second].active=false;
 805 	      }  
 806 	    }//found selected Hits       
 807 	  }//trinityVector.size !=2
 808 	}//levels do not agree
 809       }//benefit==0  
 810     } // mergeVector entries
 811   }//merge clusters
 812   if (debug >2 )
 813     std::cout<<" ***** FindCores2 is done ***"<<std::endl;
 814 }//FindCores2
 815 
 816 
 817 
 818 //called by FindCores2, asigning neighbouring hits to cluster
 819 void cluster5 ( vector<Superhit2*>* superHitVector, vector<Tmpcl2*>* clusterVector, int debug ) 
 820 {
 821   if ( (*superHitVector).size()!=0 ) { 
 822     if ( debug > 5 )
 823       std::cout<<"DEBUG KITutil:  * clustering neighbouring hits of same level "<<std::endl;
 824 
 825     unsigned int superHitVectorSize = (*superHitVector).size();
 826     
 827     for ( unsigned int hitCounter=0; hitCounter<superHitVectorSize; hitCounter++ ) {
 828       (*superHitVector)[hitCounter]->is_assigned = false;
 829       //fill superHitVector into temporary clusterVector
 830       if ( (*superHitVector)[hitCounter]->top==0 ) {
 831 	Tmpcl2* tempCluster = new Tmpcl2();  
 832 	tempCluster->hits.push_back((*superHitVector)[hitCounter]);
 833 	(*superHitVector)[hitCounter]->is_assigned = true;
 834 	clusterVector->push_back(tempCluster);
 835 	continue;
 836 	
 837       } else if ( (*superHitVector)[hitCounter]->is_assigned==false ) {	  
 838 	vector<Superhit2*> sshv;
 839 	stack <Superhit2*> superhitStack;
 840 	sshv.push_back((*superHitVector)[hitCounter]);
 841 	
 842 	//assigning neighbouring hits
 843 	for ( unsigned int neighbourCounter=0; neighbourCounter<(*superHitVector)[hitCounter]->neighbours.size(); neighbourCounter++ ) {
 844 	  if ( !((*superHitVector)[hitCounter]->neighbours[neighbourCounter]->is_assigned) 
 845 	       && superHitVector->end() != find(superHitVector->begin(),superHitVector->end(),(*superHitVector)[hitCounter]->neighbours[neighbourCounter]) ) {
 846 	    
 847 	    sshv.push_back((*superHitVector)[hitCounter]->neighbours[neighbourCounter]);
 848 	    (*superHitVector)[hitCounter]->neighbours[neighbourCounter]->is_assigned = true;
 849 	    superhitStack.push((*superHitVector)[hitCounter]->neighbours[neighbourCounter]);
 850 	  }
 851 	}//all neighbours
 852 	  
 853 	while ( !superhitStack.empty() ) {
 854 	  
 855 	  Superhit2* superHit = superhitStack.top();
 856 	  superhitStack.pop();
 857 	    
 858 	  for ( unsigned int neighbourCounter=0; neighbourCounter<superHit->neighbours.size(); neighbourCounter++ ) {  
 859 	    if ( superHit->neighbours[neighbourCounter]!=0 ) {
 860 	      if ( ! superHit->neighbours[neighbourCounter]->is_assigned &&
 861 		   superHitVector->end() != find(superHitVector->begin(),superHitVector->end(),superHit->neighbours[neighbourCounter]) ) 
 862 		if ( sshv.end() == find(sshv.begin(),sshv.end(),superHit->neighbours[neighbourCounter]) ) {
 863 		  sshv.push_back(superHit->neighbours[neighbourCounter]);
 864 		  superHit->neighbours[neighbourCounter]->is_assigned=true;
 865 		  
 866 		  superhitStack.push(superHit->neighbours[neighbourCounter]);
 867 		}
 868 	    } // neighbours not 0
 869 	  } // all neighbours  
 870 	  if ( debug >5 )
 871 	    std::cout<<"DEBUG KITutil:  * found neighbours "<<superHit->neighbours.size()<<std::endl;
 872 	} // stack not empty
 873 	
 874 	if ( sshv.size()!=0 ) {
 875 	  Tmpcl2* tempCluster = new Tmpcl2();
 876 	  for ( unsigned int sshvCounter=0; sshvCounter<sshv.size(); sshvCounter++ ) {
 877 	    tempCluster->hits.push_back(sshv[sshvCounter]);
 878 	  }
 879 	  clusterVector->push_back(tempCluster);
 880 	}
 881       } // no superhit assigned
 882     } // over all hits
 883 
 884     for ( unsigned int superHitVectorCounter=0; superHitVectorCounter<superHitVector->size(); superHitVectorCounter++ )
 885       if ( (*superHitVector)[superHitVectorCounter]->is_assigned == false ) {
 886 	Tmpcl2* tempCluster = new Tmpcl2();
 887 	tempCluster->hits.push_back((*superHitVector)[superHitVectorCounter]);
 888 	clusterVector->push_back(tempCluster);
 889 	if ( debug >5 )
 890 	  std::cout<<"DEBUG KITutil:  * found cluster "<<tempCluster<<std::endl;
 891       }
 892   }//superhit vector size not 0
 893 }//cluster5
 894 
 895 
 896 Superhit2::~Superhit2 () {}
 897 
 898 Superhit2::Superhit2 ( float E, CalorimeterHit* chit_in ) {
 899   chit = chit_in;
 900   for ( unsigned int positionCounter=0; positionCounter<3; positionCounter++ ) {
 901     point[positionCounter]  = 0.0;   
 902     pointt[positionCounter] = chit_in->getPosition()[positionCounter];
 903   }
 904 
 905   mip = chit->getEnergy() / E; //mip calibration
 906   connect     = false;
 907   is_assigned = false;
 908   mipE = E;
 909   top  = 0;
 910   cl   = 0; 
 911   
 912   S   = 0; //stave
 913   M   = 0; //module
 914   K   = 0; //layer
 915   isCalorimeter = 0; // ecal=1, hcal=2, else=0
 916   
 917 }
 918 
 919 mySuperhit::~mySuperhit () {}
 920 
 921 mySuperhit::mySuperhit ( float E, CalorimeterHit* chit_in ) {
 922   chit = chit_in;
 923   for ( unsigned int positionCounter=0; positionCounter<3; positionCounter++ ) {
 924     point[positionCounter]  = 0.0;   
 925     pointt[positionCounter] = chit_in->getPosition()[positionCounter];
 926   }
 927 
 928   mip = chit->getEnergy() / E; //mip calibration
 929   connect     = false;
 930   is_assigned = false;
 931   mipE = E;
 932   top  = 0;
 933   cl   = 0; 
 934   
 935   S   = 0; //stave
 936   M   = 0; //module
 937   K   = 0; //layer  
 938 }
 939 
 940 
 941 double* myTmpcl:: getCenter()
 942 {
 943   return center;
 944 }
 945 
 946 
 947 myTmpcl::myTmpcl(){
 948   energy=0.0;
 949 }
 950 myTmpcl::~myTmpcl(){}
 951 
 952 double myTmpcl::getEnergy(){
 953 
 954   energy =0.0;
 955   for ( unsigned int i=0; i<hits.size(); i++ ) {
 956     energy+=hits[i]->chit->getEnergy();
 957   }
 958   return energy;
 959 }
 960 
 961 //calculating energy weighted center of cluster
 962 void myTmpcl::calcCenter()
 963 {
 964   center[0] = 0.0;
 965   center[1] = 0.0;
 966   center[2] = 0.0;
 967   double energySum = 0.0;
 968   double energy    = 0.0;
 969   for ( unsigned int hitCounter=0; hitCounter<hits.size(); hitCounter++ ) {
 970     energy = hits[hitCounter]->chit->getEnergy();
 971     
 972     center[0] += hits[hitCounter]->pointt[0] * energy;
 973     center[1] += hits[hitCounter]->pointt[1] * energy;
 974     center[2] += hits[hitCounter]->pointt[2] * energy;
 975     energySum += energy;
 976   } 
 977   center[0] = center[0] / energySum;
 978   center[1] = center[1] / energySum;
 979   center[2] = center[2] / energySum;
 980 }
 981 
 982 
 983 //finding cluster start
 984 void myTmpcl::findInertia()
 985 {
 986 
 987   double aIne[3][3];  
 988   for ( int xPosCounter(0); xPosCounter<3; ++xPosCounter )
 989     for ( int yPosCounter(0); yPosCounter < 3; ++yPosCounter ) {
 990       aIne[xPosCounter][yPosCounter] = 0.0;
 991     }
 992 
 993   
 994   for ( unsigned int hitCounter=0;  hitCounter<hits.size();  hitCounter++) {
 995 
 996     //distance hit - cluster center
 997     double distX   = hits[hitCounter]->pointt[0] - center[0];
 998     double distY   = hits[hitCounter]->pointt[1] - center[1];
 999     double distZ   = hits[hitCounter]->pointt[2] - center[2];
1000 
1001     double eHit = hits[hitCounter]->chit->getEnergy();
1002     eHit = 1.0;
1003 
1004     aIne[0][0] += eHit * (distY*distY + distZ*distZ); 
1005     aIne[1][1] += eHit * (distX*distX + distZ*distZ);
1006     aIne[2][2] += eHit * (distX*distX + distY*distY);
1007 
1008     aIne[0][1] -= eHit * distX*distY;
1009     aIne[0][2] -= eHit * distX*distZ;
1010     aIne[1][2] -= eHit * distY*distZ;
1011   }
1012 
1013   for (int i(0); i<2; ++i ) {
1014     for (int j=i+1; j<3; ++j ) {
1015       aIne[j][i] = aIne[i][j];
1016     }
1017   }
1018  
1019 
1020   gsl_matrix_view aMatrix = gsl_matrix_view_array((double*)aIne,3,3);
1021   gsl_vector* aVector     = gsl_vector_alloc(3);
1022   gsl_matrix* aEigenVec   = gsl_matrix_alloc(3,3);
1023   gsl_eigen_symmv_workspace* wa = gsl_eigen_symmv_alloc(3);
1024   gsl_eigen_symmv(&aMatrix.matrix,aVector,aEigenVec,wa);
1025   gsl_eigen_symmv_free(wa);
1026   gsl_eigen_symmv_sort(aVector,aEigenVec,GSL_EIGEN_SORT_ABS_ASC);
1027 
1028   double norm = 0.0;
1029   for ( int i(0); i < 3; i++ ) {
1030     inteigen[i] = gsl_vector_get(aVector,i);
1031     norm += inteigen[i] * inteigen[i];
1032     for (int j(0); j<3; j++ ) {
1033       inteigenvec[i+3*j] = gsl_matrix_get(aEigenVec,i,j);
1034     }
1035   }
1036   norm = sqrt(norm);
1037  
1038 
1039   for ( int i(0); i<3; i++ ) 
1040     inteigen[i] = inteigen[i]/norm;
1041 
1042   direction[0] = inteigenvec[0];
1043   direction[1] = inteigenvec[1];
1044   direction[2] = inteigenvec[2];
1045 
1046   gsl_vector_free(aVector);
1047   gsl_matrix_free(aEigenVec);
1048 }//find inertia
1049 
1050 
1051 double* Tmpcl2:: getCenter()
1052 {
1053   return center;
1054 }
1055 
1056 
1057 Tmpcl2::Tmpcl2(){
1058   energy=0.0;
1059 }
1060 Tmpcl2::~Tmpcl2(){}
1061 
1062 double Tmpcl2::getEnergy(){
1063 
1064   energy =0.0;
1065   for ( unsigned int i=0; i<hits.size(); i++ ) {
1066     energy+=hits[i]->chit->getEnergy();
1067   }
1068   return energy;
1069 }
1070 
1071 //calculating energy weighted center of cluster
1072 void Tmpcl2::calcCenter()
1073 {
1074   center[0] = 0.0;
1075   center[1] = 0.0;
1076   center[2] = 0.0;
1077   double energySum = 0.0;
1078   double energy    = 0.0;
1079   for ( unsigned int hitCounter=0; hitCounter<hits.size(); hitCounter++ ) {
1080     energy = hits[hitCounter]->chit->getEnergy();
1081     
1082     center[0] += hits[hitCounter]->pointt[0] * energy;
1083     center[1] += hits[hitCounter]->pointt[1] * energy;
1084     center[2] += hits[hitCounter]->pointt[2] * energy;
1085     energySum += energy;
1086   } 
1087   center[0] = center[0] / energySum;
1088   center[1] = center[1] / energySum;
1089   center[2] = center[2] / energySum;
1090 }
1091 
1092 
1093 //finding cluster start
1094 void Tmpcl2::findInertia()
1095 {
1096 
1097   double aIne[3][3];  
1098   for ( int xPosCounter(0); xPosCounter<3; ++xPosCounter )
1099     for ( int yPosCounter(0); yPosCounter < 3; ++yPosCounter ) {
1100       aIne[xPosCounter][yPosCounter] = 0.0;
1101     }
1102 
1103   
1104   for ( unsigned int hitCounter=0;  hitCounter<hits.size();  hitCounter++) {
1105 
1106     //distance hit - cluster center
1107     double distX   = hits[hitCounter]->pointt[0] - center[0];
1108     double distY   = hits[hitCounter]->pointt[1] - center[1];
1109     double distZ   = hits[hitCounter]->pointt[2] - center[2];
1110 
1111     double eHit = hits[hitCounter]->chit->getEnergy();
1112     eHit = 1.0;
1113 
1114     aIne[0][0] += eHit * (distY*distY + distZ*distZ); 
1115     aIne[1][1] += eHit * (distX*distX + distZ*distZ);
1116     aIne[2][2] += eHit * (distX*distX + distY*distY);
1117 
1118     aIne[0][1] -= eHit * distX*distY;
1119     aIne[0][2] -= eHit * distX*distZ;
1120     aIne[1][2] -= eHit * distY*distZ;
1121   }
1122 
1123   for (int i(0); i<2; ++i ) {
1124     for (int j=i+1; j<3; ++j ) {
1125       aIne[j][i] = aIne[i][j];
1126     }
1127   }
1128  
1129 
1130   gsl_matrix_view aMatrix = gsl_matrix_view_array((double*)aIne,3,3);
1131   gsl_vector* aVector     = gsl_vector_alloc(3);
1132   gsl_matrix* aEigenVec   = gsl_matrix_alloc(3,3);
1133   gsl_eigen_symmv_workspace* wa = gsl_eigen_symmv_alloc(3);
1134   gsl_eigen_symmv(&aMatrix.matrix,aVector,aEigenVec,wa);
1135   gsl_eigen_symmv_free(wa);
1136   gsl_eigen_symmv_sort(aVector,aEigenVec,GSL_EIGEN_SORT_ABS_ASC);
1137 
1138   double norm = 0.0;
1139   for ( int i(0); i < 3; i++ ) {
1140     inteigen[i] = gsl_vector_get(aVector,i);
1141     norm += inteigen[i] * inteigen[i];
1142     for (int j(0); j<3; j++ ) {
1143       inteigenvec[i+3*j] = gsl_matrix_get(aEigenVec,i,j);
1144     }
1145   }
1146   norm = sqrt(norm);
1147  
1148 
1149   for ( int i(0); i<3; i++ ) 
1150     inteigen[i] = inteigen[i]/norm;
1151 
1152   direction[0] = inteigenvec[0];
1153   direction[1] = inteigenvec[1];
1154   direction[2] = inteigenvec[2];
1155 
1156   gsl_vector_free(aVector);
1157   gsl_matrix_free(aEigenVec);
1158 }//find inertia
1159 
1160 
1161 //find shower start in first ECAL layer?
1162 //where does 400 come from?
1163 void LineCaloIntersectD2( double* X1,double* dir,double& r_inner,double& zmax,double* X,int debug)
1164 {  //X1: cluster center; dir: cluster direction; r_inner: inner barrel radius, zmax: maximum barrel z-value; X: 0 
1165 
1166   if(debug >4)
1167     std::cout<<"---- FindCores2 called LineCaloIntersect -----"<<std::endl;
1168 
1169   double X2[3]; //cluster center + direction 
1170   //distinguish between upper and lower barrel half
1171   X2[0] = X1[0] + dir[0] * 400.0;
1172   X2[1] = X1[1] + dir[1] * 400.0;
1173   X2[2] = X1[2] + dir[2] * 400.0;
1174 
1175   LineCaloIntersect2(X2,X1,r_inner,zmax,X,debug);
1176 
1177   if ( X[0]==0.0 ) {
1178     
1179     X2[0] = X1[0] - dir[0] * 400.0;
1180     X2[1] = X1[1] - dir[1] * 400.0;
1181     X2[2] = X1[2] - dir[2] * 400.0;
1182     
1183     LineCaloIntersect2(X2,X1,r_inner,zmax,X,debug);
1184       
1185   }
1186   if(debug >4)
1187     std::cout<<"---- Finished LineCaloIntersect ----"<<std::endl;
1188   
1189 }
1190 
1191 
1192 // CAREFULL hardcoded numbers!!
1193 // written for LDC00_02 -> assumes 8fold detector geometry + 2 endcaps!
1194 //searches for cluster closest to IP
1195 void LineCaloIntersect2(double* X1, double* X2,double& r_inner,double& zmax,double* X, int debug)
1196 {
1197   double  n[10][3];
1198   
1199   double sqrt2_half = M_SQRT2/2.0; //= sqrt(2)/2 = 0.71
1200    
1201   //n[stave][xyz] ->stave geometry
1202   n[0][0] = 0.0; 
1203   n[0][1] = 1.0; 
1204   n[0][2] = 0.0;
1205   
1206   n[1][0] = sqrt2_half; 
1207   n[1][1] = sqrt2_half;
1208   n[1][2] = 0.0;
1209 
1210   n[2][0] = 1.0; 
1211   n[2][1] = 0.0;
1212   n[2][2] = 0.0;
1213   
1214   n[3][0] = -1.0; 
1215   n[3][1] = 0.0;
1216   n[3][2] = 0.0;
1217   
1218   n[4][0] = 0.0; 
1219   n[4][1] = -1.0;
1220   n[4][2] = 0.0;
1221 
1222   n[5][0] = -sqrt2_half; 
1223   n[5][1] = -sqrt2_half;
1224   n[5][2] = 0.0;
1225  
1226   n[6][0] = sqrt2_half; 
1227   n[6][1] = -sqrt2_half;
1228   n[6][2] = 0.0;
1229 
1230   n[7][0] = -sqrt2_half; 
1231   n[7][1] = sqrt2_half;
1232   n[7][2] = 0.0;
1233 
1234   //endcaps	
1235   n[8][0] = 0.0; 
1236   n[8][1] = 0.0;
1237   n[8][2] = 1.0;
1238 
1239   n[9][0] = 0.0; 
1240   n[9][1] = 0.0;
1241   n[9][2] = -1.0;
1242  
1243  
1244 
1245   if ( abs(X2[2])<zmax ) { //hit in barrel:
1246 
1247    
1248   if(debug >5)
1249     std::cout<<"---- LineCaloIntersect2 found hit in barrel ----"<<std::endl;
1250 
1251     vector < vector <double> > test;
1252     for ( unsigned int staveCounter=0; staveCounter<8; staveCounter++ ) {
1253       
1254       double tmp = - (n[staveCounter][0]*n[staveCounter][0] + n[staveCounter][1]*n[staveCounter][1] + n[staveCounter][2]*n[staveCounter][2]) * r_inner
1255 	+ n[staveCounter][0]*X1[0] + n[staveCounter][1]*X1[1] + n[staveCounter][2]*X1[2];
1256       
1257       double tmp2 = n[staveCounter][0]*(X1[0]-X2[0]) 
1258 	+ n[staveCounter][1]*(X1[1]-X2[1]) 
1259 	+ n[staveCounter][2]*(X1[2]-X2[2]);
1260 	    
1261       double t = -2.0;     
1262       if ( tmp2!=0.0 ) {
1263 	t = tmp/tmp2; 
1264 	//for stave[0]: t = (-r_inner+clusterCenterY)/(400*clusterDirectionY)
1265       }
1266 
1267       if (t>0.0 && t<=1.0 ) {
1268 	//for stave[0]: if (clusterCenterY>r_inner && -r_inner+clusterCenterY<400*clusterDirectionY) ->hit inside ECAL
1269 	vector <double> tmpi;
1270 	tmpi.push_back(X1[0]+(X2[0]-X1[0])*t); //2*clusterCenterX-r_inner ->shifts cluster center 2*depper into the ECAL stave 
1271 	tmpi.push_back(X1[1]+(X2[1]-X1[1])*t);
1272 	tmpi.push_back(X1[2]+(X2[2]-X1[2])*t);
1273 	test.push_back(tmpi);
1274       }
1275     }//all staves
1276              
1277     if (test.size()==1 ) {
1278       
1279       X[0] = test[0][0];
1280       X[1] = test[0][1];
1281       X[2] = test[0][2];
1282       double rts = sqrt(X[0]*X[0] + X[1]*X[1]);
1283       if ( rts< ((r_inner+2.0)*1.0823922002924) ) { //why +2 ?? why *1.08...??
1284 	return; //check wether hit still in ECAL??
1285       } else {
1286 	X[0] = 0.0;
1287 	X[1] = 0.0;
1288 	X[2] = 0.0;
1289 	return;//if so, seteverything back
1290       } 
1291     } //if test positive for exactly one stave
1292     else {
1293 
1294       double rmin = 1e+22;
1295       unsigned int imin = 999;
1296 
1297       //loop over all positive staves
1298       for ( unsigned int i=0; i<test.size(); i++ ) {
1299 	double Xa[3];
1300 	Xa[0] = test[i][0]; 
1301 	Xa[1] = test[i][1]; 
1302 	Xa[2] = test[i][2];
1303 	
1304 	double tmpr = sqrt(Xa[0]*Xa[0] + Xa[1]*Xa[1]);
1305 	if( tmpr<rmin ) { 
1306 	  rmin = tmpr; //search smallest cluster distance to IP
1307 	  imin = i;    
1308 	}
1309       }
1310       
1311       if ( imin < test.size() ) {
1312 	X[0] = test[imin][0]; //return shifted cluster position of stave with cluster closest to ECAL surface
1313 	X[1] = test[imin][1];
1314 	X[2] = test[imin][2];
1315 	return;
1316       } else { //did not find stave smallest cluster distance
1317 	X[0] = 0.0;
1318 	X[1] = 0.0;
1319 	X[2] = 0.0;
1320 	return;
1321       }
1322     }// testsize != 1
1323    
1324   } else { // endcap start
1325     if(debug >5)
1326       std::cout<<"---- LineCaloIntersect2 found hit in endcap ----"<<std::endl;
1327       
1328       double n[8][3];
1329       if ( X2[2]>0.0 ) {
1330 	n[0][0] = 0.0; 
1331 	n[0][1] = 0.0;
1332 	n[0][2] = 1.0;
1333       } else {
1334 	n[0][0] = 0.0; 
1335 	n[0][1] = 0.0;
1336 	n[0][2] = -1.0;
1337       }
1338 
1339       
1340       double tmp  = -(n[0][0]*n[0][0] + n[0][1]*n[0][1] + n[0][2]*n[0][2]) * zmax
1341 	+ n[0][0]*X1[0] + n[0][1]*X1[1] + n[0][2]*X1[2];
1342       double tmp2 = n[0][0]*(X1[0]-X2[0]) 
1343 	+ n[0][1]*(X1[1]-X2[1]) 
1344 	+ n[0][2]*(X1[2]-X2[2]);
1345       double t    = 0.0;
1346       
1347       if ( tmp2!=0.0 ) t=tmp/tmp2;
1348       
1349       if ( t>0.0 && t<=1.0 ) {
1350 	X[0] = X1[0]+(X2[0]-X1[0])*t ;
1351 	X[1] = X1[1]+(X2[1]-X1[1])*t ;
1352 	X[2] = X1[2]+(X2[2]-X1[2])*t ; 
1353 	return;
1354       } else {
1355 	X[0] = 0.0;
1356 	X[1] = 0.0;
1357 	X[2] = 0.0;
1358 	return;
1359       }
1360     }//endcap end
1361 }
1362 void ClusterInCluster(myTmpcl* cl, vector<myTmpcl*>& clv,vector<myTmpcl*>& clout)
1363 {
1364   sort(cl->hits.begin(),cl->hits.end());//sort hits in first cluster
1365   for (unsigned int i=0; i<clv.size(); i++) {//loop over 2nd cluster
1366     sort(clv[i]->hits.begin(),clv[i]->hits.end());//sort hits in 2nd cluster
1367     if ( includes(cl->hits.begin(),cl->hits.end(),clv[i]->hits.begin(),clv[i]->hits.end()) )
1368       //every element of clusterVector[i] included in cl 
1369       clout.push_back(clv[i]);//return included cluster
1370       }
1371   //  std::cout<<"ClusterInCluster loop over "<<clv.size() <<" hits of clusterVector included in photonSeed2 -> filled tempCluster"<<std::endl;
1372 }
1373 
1374 void ClusterInCluster2(Tmpcl2* cl, vector<Tmpcl2*>& clusterVector, int debug)
1375 {
1376  
1377    sort(cl->hits.begin(),cl->hits.end());
1378 
1379    for ( unsigned int vectorCounter=0; vectorCounter<clusterVector.size(); vectorCounter++ ) {
1380 
1381      sort(clusterVector[vectorCounter]->hits.begin(),clusterVector[vectorCounter]->hits.end());
1382 
1383      if ( includes(cl->hits.begin(),cl->hits.end(),clusterVector[vectorCounter]->hits.begin(),clusterVector[vectorCounter]->hits.end()) ) { 
1384        cl->daughters.push_back(clusterVector[vectorCounter]);
1385        clusterVector[vectorCounter]->parents.push_back(cl);
1386      }
1387    }   
1388 }
1389 
1390 
1391 //called by FindCores2
1392 //returns distance of included cluster hit to surrounding cluster hit ??
1393 double LinePointDistance2( double* X1, double* X2, double* X0){
1394 
1395   double tmp1 = X1[1]*X2[0] - X1[0]*X2[1] - X1[1]*X0[0] + X2[1]*X0[0] + X1[0]*X0[1] - X2[0]*X0[1]; 
1396   double tmp2 = X1[2]*X2[0] - X1[0]*X2[2] - X1[2]*X0[0] + X2[2]*X0[0] + X1[0]*X0[2] - X2[0]*X0[2]; 
1397   double tmp3 = X1[2]*X2[1] - X1[1]*X2[2] - X1[2]*X0[1] + X2[2]*X0[1] + X1[1]*X0[2] - X2[1]*X0[2]; 
1398 
1399   double tmp4 = (X1[0]-X2[0])*(X1[0]-X2[0]) + (X1[1]-X2[1])*(X1[1]-X2[1]) + (X1[2]-X2[2])*(X1[2]-X2[2]);
1400   
1401   double tmp5 = tmp1*tmp1 + tmp2*tmp2 + tmp3*tmp3;
1402   tmp5 = tmp5/tmp4;
1403   return sqrt(tmp5);
1404 }
1405 
1406 //called by findCores2
1407 void ModuleNormal2(double* X1,double& zmax, double* X0)
1408 {
1409   if ( abs(X1[2])<zmax ) { //in barrel
1410     
1411     //for hits in first layer, phi = angle between hit and y axis
1412     //the following section determines if the first layer hit is a neighbour of the next stave. 
1413     //If so it is assigned to both stave vectors -> very similar to precalculation
1414 
1415     double phi  = atan2(X1[0],X1[1]);
1416     double phi0 = M_PI_4/2.0; //pi/4/2 = 0.39
1417    
1418     int i;
1419     if ( phi>=0. ) 
1420       i = (int)((phi+phi0) / M_PI_4);
1421     else 
1422       i = (int)((phi-phi0) / M_PI_4);
1423 
1424     double phi_p = i * M_PI_4;
1425     if ( i!=0 ) {
1426 
1427       X0[0] = sin(phi_p);
1428       X0[1] = cos(phi_p);
1429       X0[2] = 0.0;
1430     } else {
1431 
1432       X0[0] = 0.0;
1433       X0[1] = 1.0;
1434       X0[2] = 0.0;
1435     }
1436   } else {
1437     if ( X1[2]>0. ) {
1438 
1439       X0[0] = 0.0;
1440       X0[1] = 0.0;
1441       X0[2] = 1.0;
1442     } else {
1443  
1444      X0[0] = 0.0;
1445      X0[1] = 0.0;
1446      X0[2] = -1.0;
1447     }
1448   }
1449 }
1450 
1451 double D_cl_cl(myTmpcl* cl1,myTmpcl* cl2) 
1452 {
1453 
1454  vector <double> dist;
1455  
1456  for ( unsigned int cluster1entry=0; cluster1entry<cl1->hits.size(); cluster1entry++ ) {
1457  
1458    double x1 = cl1->hits[cluster1entry]->pointt[0];
1459    double y1 = cl1->hits[cluster1entry]->pointt[1];
1460    double z1 = cl1->hits[cluster1entry]->pointt[2]; 
1461  
1462    for ( unsigned int cluster2entry=0; cluster2entry<cl2->hits.size(); cluster2entry++ ) {
1463      
1464      double tmp = ( (x1 - cl2->hits[cluster2entry]->pointt[0]) * (x1 - cl2->hits[cluster2entry]->pointt[0])
1465 		    + (y1 - cl2->hits[cluster2entry]->pointt[1]) * (y1 - cl2->hits[cluster2entry]->pointt[1])
1466 		    + (z1 - cl2->hits[cluster2entry]->pointt[2]) * (z1 - cl2->hits[cluster2entry]->pointt[2]) ) ;     
1467      dist.push_back(tmp);
1468    }
1469  }
1470  sort(dist.begin(),dist.end());
1471  return sqrt(dist[0]); //return smallest distance between cluster1 & cluster2
1472 }
1473 
1474 //calculate distance between 2 clusters called by FindCores
1475 double D_cl_cl2(Tmpcl2* cl1,Tmpcl2* cl2) 
1476 {
1477 
1478  vector <double> dist;
1479  
1480  for ( unsigned int cluster1entry=0; cluster1entry<cl1->hits.size(); cluster1entry++ ) {
1481  
1482    double x1 = cl1->hits[cluster1entry]->pointt[0];
1483    double y1 = cl1->hits[cluster1entry]->pointt[1];
1484    double z1 = cl1->hits[cluster1entry]->pointt[2]; 
1485  
1486    for ( unsigned int cluster2entry=0; cluster2entry<cl2->hits.size(); cluster2entry++ ) {
1487      
1488      double tmp = ( (x1 - cl2->hits[cluster2entry]->pointt[0]) * (x1 - cl2->hits[cluster2entry]->pointt[0])
1489 		    + (y1 - cl2->hits[cluster2entry]->pointt[1]) * (y1 - cl2->hits[cluster2entry]->pointt[1])
1490 		    + (z1 - cl2->hits[cluster2entry]->pointt[2]) * (z1 - cl2->hits[cluster2entry]->pointt[2]) ) ;     
1491      dist.push_back(tmp);
1492    }
1493  }
1494  sort(dist.begin(),dist.end());
1495  return sqrt(dist[0]); //return smallest distance between cluster1 & cluster2
1496 }
1497 
1498 // called by FindCores
1499 inline double Dot2(double* X1,double* X2)
1500 {
1501   double n1 = sqrt(X1[0]*X1[0] + X1[1]*X1[1] + X1[2]*X1[2]);
1502   double n2 = sqrt(X2[0]*X2[0] + X2[1]*X2[1] + X2[2]*X2[2]);
1503 
1504   return (X1[0]*X2[0] + X1[1]*X2[1] + X1[2]*X2[2]) / (n1*n2);
1505 }
1506 
1507 
1508 void ClusterInCluster2(Tmpcl2* cl, vector<Tmpcl2*>& clusterVector, vector<Tmpcl2*>& clout, int debug)
1509 {
1510   //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
1511   //FIXME: The main problem seems to be this sort routine, but why?
1512   //This should not touch the content, only the order of the cluster!
1513   //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
1514 
1515   //sort(cl->hits.begin(),cl->hits.end());//sort hits in first cluster
1516   for (unsigned int i=0; i<clusterVector.size(); i++) {//loop over 2nd cluster
1517     sort(clusterVector[i]->hits.begin(),clusterVector[i]->hits.end());//sort hits in 2nd cluster
1518     if ( includes(cl->hits.begin(),cl->hits.end(),clusterVector[i]->hits.begin(),clusterVector[i]->hits.end()) )
1519       //every element of clusterVector[i] included in cl 
1520       clout.push_back(clusterVector[i]);//return included cluster
1521       }
1522   if ( debug >2)
1523     std::cout<<"ClusterInCluster loop over "<<clusterVector.size() <<" hits of clusterVector included in photonSeed2 -> filled tempCluster"<<std::endl;
1524 
1525 }
1526 
1527 
1528 Photon2::~Photon2(){}
1529 
1530 Photon2:: Photon2(double E_in, double* direction, double* beginning)
1531 {
1532  //photon energy
1533  double E_e = E_in;
1534 
1535  dir[0] = direction[0];
1536  dir[1] = direction[1];
1537  dir[2] = direction[2];
1538 
1539  start[0] = beginning[0];
1540  start[1] = beginning[1];
1541  start[2] = beginning[2];
1542 
1543  //my Particle Geometrical database finds hits in ECAL barrel
1544  myPGdb::ZONE zones ;
1545  zones = myPGdb::ECAL1_BAR;
1546  //get ECAL parameters 
1547  Z        = myPGDB[zones].Zeff;
1548  x0eff    = myPGDB[zones].x0eff;
1549  sampling = myPGDB[zones].sampling;
1550  Eceff    = myPGDB[zones].Eceff;
1551  eprime   = myPGDB[zones].eprime;
1552  Rm       = myPGDB[zones].Rmeff;
1553 
1554  z1 = 0.0251 + 0.00319  * log(E_e*1000.0);
1555  z2 = 0.1162 - 0.000381 * Z;
1556  k1 = 0.659  - 0.00309  * Z;
1557  k2 = 0.645;
1558  k3 = -2.59;
1559  k4 = 0.3585 + 0.0421  * log(E_e*1000.0);
1560  p1 = 2.632  - 0.00094 * Z; 
1561  p2 = 0.401  + 0.00187 * Z;
1562  p3 = 1.313  - 0.0686  * log(E_e*1000.0);
1563  y  = E_e * 1000.0 / Eceff;
1564 
1565  Fs = x0eff/sampling;
1566 
1567  Thom = log(y) - 0.858;
1568  Tsam = Thom - 0.59/Fs - 0.53*(1.0-eprime);
1569  alfahom = 0.21 + (0.492+2.38/Z) * log(y);
1570  alfasam = alfahom - 0.444/Fs;
1571  betasam = 0.5;
1572 
1573 }
1574 
1575 
1576 void Photon2::Prob(CalorimeterHit* ch,double cut,double* out)
1577 {
1578  
1579   double X[3];
1580   PointOnLine22(start, dir,ch->getPosition(),X);
1581   
1582   double TP[3];  
1583   TP[0] = ch->getPosition()[0] - start[0];
1584   TP[1] = ch->getPosition()[1] - start[1];
1585   TP[2] = ch->getPosition()[2] - start[2];
1586 
1587   // distance must have direction  
1588   if ( Dot2(TP,dir)<0.0 ) {
1589     out[0] = 0.0;
1590     out[1] = 0.0;
1591     return ;
1592   }
1593 
1594   double pos_X0 = sqrt( (start[0]-X[0])*(start[0]-X[0])+
1595 			(start[1]-X[1])*(start[1]-X[1])+
1596 			(start[2]-X[2])*(start[2]-X[2]))/x0eff;
1597   double tau =  pos_X0/Tsam; 
1598 
1599   double pos[3]; 
1600   pos[0] = ch->getPosition()[0];
1601   pos[1] = ch->getPosition()[1];
1602   pos[2] = ch->getPosition()[2];
1603   
1604   double radius = LinePointDistance2(start,dir,pos);
1605 
1606   // average radial profiles 
1607   double RChom = z1 + z2*tau;
1608   double RThom = k1 * (exp(k3*(tau-k2)) + exp(k4*(tau-k2)));
1609   double phom  = p1 * exp((p2-tau)/p3 - exp((p2-tau)/p3));
1610   
1611   // average radial profiles sampling
1612   double RCsam = RChom - 0.0203*(1.0-eprime) + 0.0397*exp(-tau)/Fs;
1613   double RTsam = RThom - 0.14*(1.0-eprime) - 0.495*exp(-tau)/Fs;
1614   double psam  = phom  + (1.0-eprime)*(0.348-0.642*exp(-pow(tau-1.0,2.0))/Fs);
1615   
1616   double rb = radius/Rm;  
1617   
1618 if ( rb<20.0 && pos_X0<35.0 ) {
1619   double tmp = pow(betasam*pos_X0,alfasam-1.0)*betasam*exp(-betasam*pos_X0)/gsl_sf_gamma(alfasam);
1620   tmp = Ee*tmp*(psam*2.0*rb*RCsam*RCsam/pow(rb*rb+RCsam*RCsam,2.0) +
1621 		  (1-psam)*(2.0*rb*RTsam*RTsam/pow(rb*rb+RTsam*RTsam,2.0)))/rb;
1622   if ( tmp>cut ) {
1623     out[0] = tmp;
1624     out[1] = radius;
1625     
1626     //    std::cout<<"set probability to "<<tmp<<std::endl;
1627     
1628   } else {
1629     out[0] = 0.0;
1630     out[1] = 0.0;
1631     
1632     //    std::cout<<"probability "<<tmp<<" < "<<cut<<" -> set to 0."<<std::endl;
1633   }
1634 } else {
1635     out[0] = 0.0;
1636     out[1] = 0.0;
1637     
1638     //    std::cout<<"radius "<<rb<<" < 20 || pos_X0 "<<pos_X0<<" < 35 -> set to 0."<<std::endl;
1639   }
1640   
1641 }
1642 
1643 
1644 
1645 void PointOnLine3(const double* X1,const double* X2,const float* X0,double* Xline)
1646 {
1647 
1648   double tmp1 = (X1[0]-X0[0]) * (X2[0]-X1[0]) + (X1[1]-X0[1]) * (X2[1]-X1[1]) + (X1[2]-X0[2]) * (X2[2]-X1[2]);
1649   double tmp4 = (X1[0]-X2[0]) * (X1[0]-X2[0]) + (X1[1]-X2[1]) * (X1[1]-X2[1]) + (X1[2]-X2[2]) * (X1[2]-X2[2]);
1650   double t = -tmp1/tmp4;
1651 
1652 
1653  Xline[0] = X1[0] + (X2[0]-X1[0])*t;
1654  Xline[1] = X1[1] + (X2[1]-X1[1])*t;
1655  Xline[2] = X1[2] + (X2[2]-X1[2])*t;
1656  
1657 }
1658 void PointOnLine22(const double* Xstart,const double* dir,const float* X0,double* Xline)
1659 {
1660   double X1[3];
1661   double X2[3];
1662   X1[0] = Xstart[0] - dir[0]*100.0;
1663   X1[1] = Xstart[1] - dir[1]*100.0;
1664   X1[2] = Xstart[2] - dir[2]*100.0;
1665 
1666   X2[0] = Xstart[0] + dir[0]*100.0;
1667   X2[1] = Xstart[1] + dir[1]*100.0;
1668   X2[2] = Xstart[2] + dir[2]*100.0;
1669 
1670 
1671   PointOnLine3(X1,X2,X0,Xline);
1672 
1673 }
1674 
1675 double giveMeEEstimate2(int level,double Ecore, vector<CoreCalib2> cc)
1676 {
1677   for ( unsigned int i=0; i<cc.size(); i++ ) { 
1678     //    std::cout<<cc[i].level<<" == "<<level<<" ? && "<<cc[i].Emax<<" > "<<Ecore<<" ?"<<std::endl;
1679     if ( cc[i].level==level && cc[i].Emax>Ecore ) {  //FIXME: does find more than one appropriate calibration! 
1680       /*
1681       std::cout << "For level " << level 
1682 		<< " core with energy: " << Ecore 
1683 		<< " : estimated photon energy = cc[i].b: " << cc[i].b 
1684 		<< " + cc[i].a: " << cc[i].a 
1685 		<< " * Ecore = " << cc[i].b + cc[i].a*Ecore << std::endl;
1686       */
1687       if ( cc[i].b==0.0 || cc[i].a==0.0 )
1688 	return Ecore;                      //if no valid calibration, return uncalibrated value
1689       else
1690 	return cc[i].b + cc[i].a*Ecore;
1691     }
1692     else 
1693       return Ecore;                        //if no valid calibration, return uncalibrated value
1694   }
1695   return 0.0;
1696 }
1697 
1698 
1699 // !!! has to be adopted for different detector models !!!
1700 // relates measured energy in cluster to estimated true energy of the cluster, 
1701 // depending on the cluster level
1702 // numbers = energy cut in the enery per hit spectrum per level, for various energies
1703 // E_true = aa*E_cluster + b; E_min < E_cluster < E_max
1704 // can be obtained with KALIBRATOR.cc
1705 
1706 void CreateCalibration_LDC00(vector<CoreCalib2>* cc)
1707 {
1708 
1709   
1710   double aEnom[9]  = {0.5,1.0,2.0,3.0,5.0,8.0,12.0,20.0,40.0};
1711   double aa[9][10] = {{0.283,0.357,0.337,0.398,0.000,0.000,0.000,0.000,0.000,0.000},
1712 		      {0.703,0.346,0.814,0.785,0.939,0.000,0.000,0.000,0.000,0.000},
1713 		      {0.505,0.821,1.000,1.112,1.527,1.860,0.000,0.000,0.000,0.000},
1714 		      {0.809,0.892,1.999,1.554,2.170,2.870,0.000,0.000,0.000,0.000},
1715 		      {1.075,0.827,1.290,2.350,2.123,3.958,4.814,0.000,0.000,0.000},
1716 		      {1.140,1.310,1.765,3.652,4.030,5.600,6.670,0.000,0.000,0.000},
1717 		      {0.760,1.845,3.444,2.850,5.150,6.760,10.72,0.000,0.000,0.000},
1718 		      {3.770,4.460,2.270,5.930,8.730,11.00,17.40,17.77,0.000,0.000},
1719 		      {1.890,3.560,5.120,7.320,9.650,11.71,21.74,0.000,0.000,0.000}};
1720 
1721   double bb[9][10] = {{0.730,0.595,0.772,0.683,0.000,0.000,0.000,0.000,0.000,0.000},
1722 		      {0.493,1.079,0.458,0.592,0.447,0.000,0.000,0.000,0.000,0.000},
1723 		      {0.941,0.813,0.766,0.782,0.559,0.405,0.000,0.000,0.000,0.000},
1724 		      {0.892,0.910,0.550,0.809,0.593,0.347,0.000,0.000,0.000,0.000},
1725 		      {0.918,1.015,0.967,0.800,0.954,0.510,0.356,0.000,0.000,0.000},
1726 		      {0.974,0.986,0.970,0.778,0.784,0.593,0.509,0.000,0.000,0.000},
1727 		      {1.032,0.971,0.869,0.986,0.835,0.742,0.339,0.000,0.000,0.000},
1728 		      {0.891,0.888,1.025,0.884,0.782,0.620,0.314,0.361,0.000,0.000},
1729 		      {1.020,0.992,0.977,0.962,0.945,0.943,0.724,0.000,0.000,0.000}};
1730   
1731 double aEmin[9][10] = {{0.250,0.200,0.150,0.060,0.000,0.000,0.000,0.000,0.000,0.000},
1732 		       {0.600,0.500,0.350,0.200,0.100,0.000,0.000,0.000,0.000,0.000},
1733 		       {1.500,1.200,1.100,0.900,0.650,0.300,0.000,0.000,0.000,0.000},
1734 		       {2.400,2.100,1.900,1.600,1.300,0.600,0.000,0.000,0.000,0.000},
1735 		       {4.000,4.100,3.900,3.200,2.900,2.400,1.100,0.000,0.000,0.000},
1736 		       {7.000,6.800,6.200,5.500,5.000,4.200,2.600,0.000,0.000,0.000},
1737 		       {10.50,10.00,9.500,8.600,8.400,8.300,5.200,0.000,0.000,0.000},
1738 		       {18.50,18.20,17.40,16.00,15.50,13.00,11.50,9.000,0.000,0.000},
1739 		       {36.00,35.00,34.00,33.50,32.50,31.00,28.00,0.000,0.000,0.000}};
1740  
1741  double aEmax[9][10]={{0.600,0.500,0.350,0.300,0.000,0.000,0.000,0.000,0.000,0.000},
1742 		      {1.200,1.100,1.000,0.900,0.700,0.000,0.000,0.000,0.000,0.000},
1743 		      {2.500,2.200,1.900,1.800,1.400,1.200,0.000,0.000,0.000,0.000},
1744 		      {3.300,3.100,3.000,2.700,2.400,2.000,0.000,0.000,0.000,0.000},
1745 		      {5.600,5.500,5.200,4.800,4.300,3.900,2.900,0.000,0.000,0.000},
1746 		      {9.000,8.800,8.300,7.600,7.200,6.400,5.200,0.000,0.000,0.000},
1747 		      {13.50,13.00,12.50,11.70,11.10,11.00,8.300,0.000,0.000,0.000},
1748 		      {22.20,22.80,21.50,20.50,19.00,16.30,15.50,13.50,0.000,0.000},
1749 		      {44.00,43.50,42.50,41.50,39.20,37.50,33.50,0.000,0.000,0.000}};
1750  
1751  for ( unsigned int levelCounter=0; levelCounter<10; levelCounter++ ) {
1752    for ( unsigned int energyCounter=0; energyCounter<12; energyCounter++ ) {
1753      if ( aEmin[energyCounter][levelCounter]!=0.0 && aEmax[energyCounter][levelCounter]!=0.0 ) {
1754 
1755        CoreCalib2 cck;
1756        cck.level = levelCounter;
1757        cck.Enom  = aEnom[energyCounter];
1758        cck.a     = aa[energyCounter][levelCounter];
1759        cck.b     = bb[energyCounter][levelCounter];
1760        cck.Emin  = aEmin[energyCounter][levelCounter];
1761        cck.Emax  = aEmax[energyCounter][levelCounter];
1762        
1763        cc->push_back(cck);
1764      }
1765    }	      
1766  }
1767  
1768 }
1769 
1770 
1771 void CreateCalibration_LDCPrime02Sc(vector<CoreCalib2>* cc)
1772 {
1773   //  fit straight line: E_true = aa*E_reco + bb for each energy and level  
1774  
1775   double aEnom[13]  = {0.5,1.0,2.0,3.0,5.0,8.0,12.0,20.0,40.0,70.0,100.0,150.0,200.0};   //nominal photon energy  
1776 
1777   double aa[13][10] = {{0.390,0.386,0.875,0.327,0.000,0.000,0.000,0.000,0.000,0.000},
1778 		       {0.574,0.437,0.365,0.380,0.378,0.000,0.000,0.000,0.000,0.000},
1779 		       {0.693,0.644,0.735,0.701,0.663,0.386,0.000,0.000,0.000,0.000},
1780 		       {0.992,1.003,0.703,0.920,0.705,0.666,0.000,0.000,0.000,0.000},
1781 		       {0.942,0.894,0.857,0.685,0.851,0.719,0.000,0.000,0.000,0.000},
1782 		       {0.850,0.876,1.036,0.806,0.746,0.532,0.000,0.000,0.000,0.000},
1783 		       {0.911,0.901,0.851,0.847,0.796,0.693,0.468,0.464,0.000,0.000},
1784 		       {0.832,1.005,1.006,0.920,0.975,0.865,0.655,0.441,0.472,0.000},
1785 		       {0.969,0.899,0.854,0.939,0.858,0.762,0.716,0.378,0.417,0.000},
1786 		       {1.056,1.073,1.046,0.894,0.879,1.127,1.140,1.199,0.308,0.650},
1787 		       {0.738,1.080,2.044,1.089,1.032,0.862,0.885,0.758,1.294,0.422},
1788 		       {0.879,0.879,0.879,0.879,0.879,0.879,0.879,0.879,0.000,0.000},
1789 		       {0.836,0.836,0.836,0.836,0.836,0.836,0.386,0.836,0.000,0.000}};
1790 
1791   double bb[13][10] = {{0.412,0.436,0.383,0.471,0.000,0.000,0.000,0.000,0.000,0.000},
1792 		       {0.690,0.784,0.870,0.886,0.926,0.000,0.000,0.000,0.000,0.000},
1793 		       {1.095,1.209,1.155,1.281,1.422,1.748,0.000,0.000,0.000,0.000},
1794 		       {0.894,0.965,1.675,1.454,1.986,2.251,0.000,0.000,0.000,0.000},
1795 		       {1.466,1.804,2.109,2.900,2.760,3.442,0.000,0.000,0.000,0.000},
1796 		       {2.658,2.722,2.135,3.841,4.589,5.984,0.000,0.000,0.000,0.000},
1797 		       {3.158,3.588,4.480,5.106,6.226,7.719,9.918,10.47,0.000,0.000},
1798 		       {6.274,3.616,4.306,6.645,7.095,9.971,14.18,17.14,17.78,0.000},
1799 		       {5.938,9.247,11.98,10.65,14.84,20.11,25.15,33.79,34.32,0.000},
1800 		       {3.345,3.711,7.038,18.97,22.79,13.86,21.53,-6.09,62.73,54.05},
1801 		       {33.77,3.409,-81.3,8.968,17.70,36.57,43.36,60.02,35.92,85.20},
1802 		       {29.30,29.30,29.30,29.30,29.30,29.30,29.30,29.30,0.000,0.000},
1803 		       {46.13,46.13,46.13,46.13,46.13,46.13,46.13,46.13,0.000,0.000}};
1804   
1805   double aEmin[13][10] = {{0.250,0.200,0.150,0.060,0.000,0.000,0.000,0.000,0.000,0.000},
1806 			  {0.600,0.500,0.350,0.200,0.100,0.000,0.000,0.000,0.000,0.000},
1807 			  {1.500,1.200,1.100,0.900,0.650,0.300,0.000,0.000,0.000,0.000},
1808 			  {2.400,2.100,1.900,1.600,1.300,0.600,0.000,0.000,0.000,0.000},
1809 			  {4.000,4.100,3.900,3.200,2.900,2.400,0.000,0.000,0.000,0.000},
1810 			  {4.279,3.848,3.501,4.618,2.439,1.695,0.000,0.000,0.000,0.000},
1811 			  {9.189,6.657,6.144,5.337,4.811,3.780,2.749,0.954,0.000,0.000},
1812 			  {14.29,13.64,12.69,11.56,10.76,8.683,6.213,3.714,3.040,0.000},
1813 			  {31.50,30.82,28.61,28.23,25.75,22.29,18.14,13.61,9.857,0.000},
1814 			  {57.10,56.41,51.97,51.54,48.47,32.96,39.48,28.08,25.26,23.19},
1815 			  {62.30,84.58,71.61,68.98,59.43,62.63,53.12,33.60,35.83,34.79},
1816 			  {134.9,134.9,134.9,194.9,134.9,134.9,134.9,134.9,0.000,0.000},
1817 			  {182.7,182.7,182.7,182.7,182.7,165.0,135.0,112.0,0.000,0.000}};
1818  
1819   double aEmax[13][10]={{0.600,0.500,0.350,0.300,0.000,0.000,0.000,0.000,0.000,0.000},
1820 			{1.200,1.100,1.000,0.900,0.700,0.000,0.000,0.000,0.000,0.000},
1821 			{2.500,2.200,1.900,1.800,1.400,1.200,0.000,0.000,0.000,0.000},
1822 			{3.300,3.100,3.000,2.700,2.400,2.000,0.000,0.000,0.000,0.000},
1823 			{5.600,5.500,5.200,4.800,4.300,3.900,0.000,0.000,0.000,0.000},
1824 			{8.464,8.442,8.051,6.774,6.945,6.120,0.000,0.000,0.000,0.000},
1825 			{12.79,12.45,12.01,11.35,10.21,9.137,7.234,6.474,0.000,0.000},
1826 			{20.40,19.61,19.38,18.38,16.55,15.58,12.83,11.91,9.463,0.000},
1827 			{40.94,39.63,39.58,36.69,34.23,32.34,27.23,24.58,22.26,0.000},
1828 			{71.53,70.61,73.92,66.78,63.47,71.53,49.70,50.04,42.17,34.28},
1829 			{136.1,100.0,106.0,139.5,140.0,92.21,85.62,91.04,72.17,132.5},
1830 			{152.6,152.6,152.6,152.6,152.6,152.6,152.6,152.6,0.000,0.000},
1831 			{204.1,204.1,204.1,204.1,204.1,204.1,204.1,204.1,0.000,0.000}};
1832 
1833   for ( unsigned int levelCounter=0; levelCounter<10; levelCounter++ ) {
1834     for ( unsigned int energyCounter=0; energyCounter<12; energyCounter++ ) {
1835       if ( aEmin[energyCounter][levelCounter]!=0.0 && aEmax[energyCounter][levelCounter]!=0.0 ) {
1836 
1837 	CoreCalib2 cck;
1838 	cck.level = levelCounter;
1839 	cck.Enom  = aEnom[energyCounter];
1840 	cck.a     = aa[energyCounter][levelCounter];
1841 	cck.b     = bb[energyCounter][levelCounter];
1842 	cck.Emin  = aEmin[energyCounter][levelCounter];
1843 	cck.Emax  = aEmax[energyCounter][levelCounter];
1844        
1845 	cc->push_back(cck);
1846       }
1847     }	      
1848   }
1849   
1850 }//creatCalibration_LDCPrime02Sc
1851 

Attached Files

To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.
  • [get | view] (2010-03-23 09:57:46, 33.4 KB) [[attachment:KIT.cc]]
  • [get | view] (2010-03-23 09:58:09, 4.3 KB) [[attachment:KIT.h]]
  • [get | view] (2010-03-23 10:20:16, 10.0 KB) [[attachment:KITcalibrator.cc]]
  • [get | view] (2010-03-23 12:53:20, 2.9 KB) [[attachment:KITcalibrator.h]]
  • [get | view] (2010-03-23 09:55:11, 66.1 KB) [[attachment:KITutil.cc]]
  • [get | view] (2010-03-23 09:55:36, 8.2 KB) [[attachment:KITutil.h]]
  • [get | view] (2010-03-23 09:53:52, 22.9 KB) [[attachment:PDGB.cc]]
  • [get | view] (2010-03-23 09:54:13, 8.2 KB) [[attachment:PDGB.h]]
  • [get | view] (2010-03-23 10:26:19, 3.2 KB) [[attachment:calibrateKIT.xml]]
  • [get | view] (2010-03-23 09:33:22, 3.9 KB) [[attachment:marlinKIT.xml]]
 All files | Selected Files: delete move to page copy to page

You are not allowed to attach a file to this page.