Attachment 'KIT.cc'

Download

   1 #include "KITdev.h"
   2 
   3 using namespace lcio ;
   4 using namespace marlin ;
   5 using namespace std;
   6 using namespace IMPL;
   7 using namespace EVENT;
   8 using namespace UTIL;
   9                          	
  10 
  11 /*
  12 Photon finding Algorythm by Predrag Krstonosic.
  13 depends on Superhits, KITutil (for calibration), 
  14 Physics Geometry Database myPGDB (for pre-calculation needed for NN Clustering),
  15 MarlinMath, B_Util
  16 
  17 reads in CalorimeterHits, and writes out reconstructed Photons
  18 developed and optimized for low energetic photons coming from the IP in LDC-00_02
  19 
  20 Adopted to LDCPrime_02Sc by Nanda Wattimena
  21 */
  22 
  23 KITdev aKITdev ;
  24 KITdev::KITdev() : Processor("KITdev") {
  25   
  26   _description = "KITdev: an emShowerFinder for LDC" ;
  27   
  28   registerProcessorParameter( "ECAL_Collection",
  29 			      "Name of the CalorimeterHit Collection for ECAL ",
  30 			      _Ecal_col,
  31 			      std::string("ECAL"));
  32   registerProcessorParameter( "HCAL_Collection",
  33 			      "Name of the CalorimeterHit Collection for HCAL ",
  34 			      _Hcal_col,
  35 			      std::string("HCAL"));
  36 
  37   registerProcessorParameter( "Core_Collection",
  38 			      "Name of the Cluster Collection for Cores ",
  39 			      _CoreCollection,
  40 			      std::string("CORE"));
  41 
  42   registerProcessorParameter( "detector_Model",
  43 			      "chose detector geometry: LDC00, LDC0105, LDCPrime", //TODO: use this information!
  44 			      _detectorModel,
  45 			      std::string("LDCPrime"));
  46 
  47   registerProcessorParameter( "Cleaning",
  48 			      "Apply topological cleaning?",
  49 			      _ToClean,
  50 			      std::string("YES"));
  51   registerProcessorParameter( "TopologicalCut",
  52 			      "if so, require >= N+1 neighbors",
  53 			       _CleanCut,
  54 			      (int)5);
  55 
  56   registerProcessorParameter( "NumberOfLevels",
  57 			      "Number of levels for central loop ",
  58 			       _Nlevels,
  59 			      (int)10);
  60   //defult values of energy per hit [mip], at which to start a new level
  61   vector<float> mipstep; 
  62   mipstep.push_back(0.1);
  63   mipstep.push_back(1.5);
  64   mipstep.push_back(2.5);
  65   mipstep.push_back(4.0);
  66   mipstep.push_back(6.0);
  67   mipstep.push_back(9.0);
  68   mipstep.push_back(16.0);
  69   mipstep.push_back(26.0);
  70   mipstep.push_back(41.0);
  71   mipstep.push_back(65.0);
  72 
  73   registerProcessorParameter( "Levels",
  74 			     "level thresholds (energy per hit [MIP]) ",
  75 			      _mipstep,
  76 			       mipstep);
  77 
  78   registerProcessorParameter( "MinHit0",
  79 			     "Minimal Number of hits for ground level cluster ",
  80 			      _MinHit0,
  81 			      (int)4);
  82 
  83   registerProcessorParameter( "MinHitSplit",
  84 			      "Minimal Number of hits for i-th level cluster ",
  85 			      _MinHitSplit,
  86 			      (int)2);
  87 
  88   registerProcessorParameter( "rCut",
  89 			      "Fluctuation suprresion cut",
  90 			      _rCut,
  91 			      (double)0.4);
  92   registerProcessorParameter( "distCut",
  93 			      "Square of distance cut [mm] for merging ",
  94 			      _distCut,
  95 			      (double)35.0);
  96   registerProcessorParameter( "cosCut",
  97 			      "Cosine of the angle for merging ",
  98 			      _cosCut,
  99 			      (double)0.95);
 100 
 101   registerProcessorParameter( "probabilityDensityCut",
 102 			      "cut on the probability density to assign hits to shower cores",
 103 			      _probabilityDensityCut,
 104 			      (double)1E-3);
 105 
 106 
 107   registerProcessorParameter( "energyDeviationCut",
 108 			      "cut on energy deviation of em shower candidates from estimated energy ( abs( (Ecluster - Eestimated)/Eestimated ) < cut )",
 109 			      _energyDeviationCut,
 110 			      (double)0.25);
 111 
 112 
 113   /******  debugging tools: ******/
 114   registerProcessorParameter( "DebugLevel",
 115 			      "limits the amount of information written to std out (0 - none, 9 - maximal information)",
 116 			      _debugLevel,
 117 			      int(0) );
 118 }
 119 
 120 
 121 
 122 
 123 void KITdev::init() { 
 124      
 125    printParameters() ;
 126   _nRun = 0 ;
 127   _nEvt = 0 ;   
 128   _nSkippedEvt = 0 ;   
 129   
 130   // FIXME: hard coded cell id's for old Mokka (e.g. Mokka v5.4) versions)
 131   CellIDDecoder<CalorimeterHit>::setDefaultEncoding("M:3,S-1:3,I:9,J:9,K-1:6");
 132 
 133 
 134   _Ehit_ecal    = 0.0;
 135   _Ehit_hcal    = 0.0;
 136 
 137   _Et_miss      = 0.0;
 138   _Et_x         = 0.0;
 139   _Et_y         = 0.0;
 140   
 141   _NhitsEcal    = 0;
 142   _NhitsHcal    = 0;
 143   _NhitsCalos   = 0;
 144 
 145   _trueDist    = 0;
 146   _trueEgamma  = 0;
 147   _truePhi     = 0;
 148   _trueTheta   = 0;
 149 
 150   _trueEneutralino  = 0;
 151 
 152 }
 153 
 154 
 155 void KITdev::processRunHeader( LCRunHeader* run) { 
 156       
 157   ++_nRun;
 158 } 
 159 
 160 void KITdev::processEvent( LCEvent * evt ) {  
 161   
 162 
 163   try {   
 164     //open all kind of collections
 165     LCCollection* ecalColl = evt->getCollection(_Ecal_col.c_str()) ;
 166     LCCollection* hcalColl = evt->getCollection(_Hcal_col.c_str()); 
 167 
 168     LCCollection* trueColl = evt->getCollection("TrueClusters");
 169     LCCollection* mcColl   = evt->getCollection("MCParticle");
 170     LCCollection* LCRColl  = evt->getCollection("TrueClusterToMCP");
 171 
 172 
 173     std::vector<float> trueDistance;
 174     std::vector<float> trueGammaEnergy;
 175     std::vector<float> trueGammaStartX;
 176     std::vector<float> trueGammaStartY;
 177     std::vector<float> trueGammaStartZ;
 178     std::vector<float> trueCogX;
 179     std::vector<float> trueCogY;
 180     std::vector<float> trueCogZ;
 181     std::vector<float> trueGammaDirX;
 182     std::vector<float> trueGammaDirY;
 183     std::vector<float> trueGammaDirZ;
 184     std::vector<float> trueGammaPosX;
 185     std::vector<float> trueGammaPosY;
 186     std::vector<float> trueGammaPosZ;
 187     std::vector<float> trueNeutralinoEnergy;
 188     std::vector<float> trueNeutralinoDirX;
 189     std::vector<float> trueNeutralinoDirY;
 190     std::vector<float> trueNeutralinoDirZ;
 191     std::vector<float> trueThetaDir;
 192     std::vector<float> truePhiDir;
 193     std::vector<float> closestDistance;
 194     std::vector<int> nClusterHits;
 195     std::vector<float> phiPos;
 196     std::vector<float> thetaPos;
 197 
 198     /* get energy and vertex from photon that initiated the cluster (MCtruth) */
 199     _Ntrcl = 0;
 200     if ( trueColl!=0 ) {
 201       _Ntrcl = trueColl->getNumberOfElements();
 202       
 203 
 204       for (int i=0; i<_Ntrcl; i++) {
 205 
 206 	Cluster* cluster = dynamic_cast<Cluster*>(trueColl->getElementAt(i));
 207 	LCRelationNavigator* nav = new LCRelationNavigator(LCRColl);
 208 	const LCObjectVec& relMCParticlesToCluster = nav->getRelatedToObjects(cluster);	
 209 	
 210 	if ( relMCParticlesToCluster.size() > 1 ) 
 211 	  std::cout << "Error: More than one MCParticle related to cluster." << std::endl;
 212 
 213 	// FIXME: here is only the 0th contribution taken into account
 214 	MCParticle* mcpOfCluster = dynamic_cast<MCParticle*>(relMCParticlesToCluster.at(0));
 215 	MCParticle* mcp = dynamic_cast<MCParticle*>(mcColl->getElementAt(0));
 216 
 217 	int mcpOfClusterID = mcpOfCluster->getPDG();
 218 	int parentID = mcp->getPDG();
 219 	int nDaughters = mcp->getDaughters().size();
 220 	int daughterID[nDaughters];
 221 
 222 	//get true gamma information
 223 	if(mcpOfClusterID == 22) {
 224 	_trueEgamma      = mcpOfCluster->getEnergy();
 225 	_trueGammaPos[0] = mcpOfCluster->getVertex()[0];
 226 	_trueGammaPos[1] = mcpOfCluster->getVertex()[1];
 227 	_trueGammaPos[2] = mcpOfCluster->getVertex()[2];
 228 	
 229 	_trueDist = sqrt(  pow(_trueGammaPos[0],2) + pow(_trueGammaPos[1],2) + pow(_trueGammaPos[2],2) );
 230 	
 231 	_trueGammaDir[0] = mcpOfCluster->getMomentum()[0];
 232 	_trueGammaDir[1] = mcpOfCluster->getMomentum()[1];
 233 	_trueGammaDir[2] = mcpOfCluster->getMomentum()[2];
 234 	
 235 	double norm = sqrt( pow(_trueGammaDir[0],2) + pow(_trueGammaDir[1],2) + pow(_trueGammaDir[2],2) );
 236 	
 237 	_trueGammaDir[0] /= norm;
 238 	_trueGammaDir[1] /= norm;
 239 	_trueGammaDir[2] /= norm;
 240 	
 241 	_truePhi   = atan2(_trueGammaDir[1],_trueGammaDir[0]) *180/M_PI;
 242 	_trueTheta = acos(_trueGammaDir[2]) *180/M_PI;
 243 	}
 244 	//get true neutralino information
 245 	if( parentID == 1000022) {
 246 	  _trueEneutralino = mcp->getEnergy();
 247 	  
 248 	  _trueNeutralinoDir[0] = mcp->getMomentum()[0];
 249 	  _trueNeutralinoDir[1] = mcp->getMomentum()[1];
 250 	  _trueNeutralinoDir[2] = mcp->getMomentum()[2];
 251 	  
 252 	  double norm = sqrt( pow(_trueNeutralinoDir[0],2) + pow(_trueNeutralinoDir[1],2) + pow(_trueNeutralinoDir[2],2) );
 253 	  
 254 	  _trueNeutralinoDir[0] /= norm;
 255 	  _trueNeutralinoDir[1] /= norm;
 256 	  _trueNeutralinoDir[2] /= norm;
 257 
 258 	  if (_debugLevel >1 ) {
 259 	    std::cout<<"/*** MC truth ***/"<<std::endl;
 260 	    std::cout<<"found neutralino "<<std::endl;
 261 	    std::cout<<"       - with energy:   "<<_trueEneutralino<<" GeV"<<std::endl;
 262 	    std::cout<<"       - with momentum: "<<_trueNeutralinoDir[0]<<", "<<_trueNeutralinoDir[1]<<", "<<_trueNeutralinoDir[2]<<std::endl;
 263 	    std::cout<<"/*** *** *** ***/"<<std::endl;
 264 	  }
 265 	}//is neutralino
 266 
 267 	else {
 268 	  if (_debugLevel >1 ) 
 269 	    std::cout<<"!!! Did not find neutralino !!!"<<std::endl;
 270 	  
 271 	  _trueEneutralino      = 0;	  
 272 	  _trueNeutralinoDir[0] = 0;
 273 	  _trueNeutralinoDir[1] = 0;
 274 	  _trueNeutralinoDir[2] = 0;
 275 	}
 276 
 277 	
 278 	if ( _debugLevel > 0 ) {
 279 	  std::cout<<"/*** MC truth ***/"<<std::endl;
 280 	  std::cout<<"found cluster "<<std::endl;
 281 	  std::cout<<"       - with energy: "<<_trueEgamma<<" GeV"<<std::endl;
 282 	  std::cout<<"       - at position: "<<_trueGammaPos[0]<<", "<<_trueGammaPos[1]<<", "<<_trueGammaPos[2]<<std::endl;
 283 	  std::cout<<"       - distance to IP: "<<_trueDist<<" mm"<<std::endl;
 284 	  std::cout<<"       - momentum: "<<_trueGammaDir[0]<<", "<<_trueGammaDir[1]<<", "<<_trueGammaDir[2]<<std::endl;
 285 	  std::cout<<"       - phi: "<<_truePhi<<" degree"<<std::endl;
 286 	  std::cout<<"       - theta: "<<_trueTheta<<" degree"<<std::endl;
 287 	  std::cout<<"/*** *** *** ***/"<<std::endl;
 288 	}
 289    
 290 	
 291 	trueDistance.push_back(_trueDist);
 292 	trueGammaEnergy.push_back(_trueEgamma);
 293 	trueGammaStartX.push_back(_trueGammaPos[0]);
 294 	trueGammaStartY.push_back(_trueGammaPos[1]);
 295 	trueGammaStartZ.push_back(_trueGammaPos[2]);
 296 	trueCogX.push_back(_trueCog[0]);
 297 	trueCogY.push_back(_trueCog[1]);
 298 	trueCogZ.push_back(_trueCog[2]);
 299 	trueGammaDirX.push_back(_trueGammaDir[0]);
 300 	trueGammaDirY.push_back(_trueGammaDir[1]);
 301 	trueGammaDirZ.push_back(_trueGammaDir[2]);
 302 	trueThetaDir.push_back(_trueTheta);	    
 303 	truePhiDir.push_back(_truePhi);	  
 304 
 305 	trueNeutralinoEnergy.push_back(_trueEneutralino);  	
 306 	trueNeutralinoDirX.push_back(_trueNeutralinoDir[0]);
 307 	trueNeutralinoDirY.push_back(_trueNeutralinoDir[1]);
 308 	trueNeutralinoDirZ.push_back(_trueNeutralinoDir[2]);
 309       }//all true cluster 
 310       
 311       if ( _debugLevel > 0 ) 
 312 	std::cout<<"KITdev found "<<_Ntrcl<<" true cluster in event "<<_nEvt<<std::endl;
 313     }//true cluster collection not empty	
 314   
 315     
 316     CellIDDecoder<CalorimeterHit>CDECAL = CellIDDecoder<CalorimeterHit>(ecalColl);
 317     CellIDDecoder<CalorimeterHit>CDHCAL = CellIDDecoder<CalorimeterHit>(hcalColl);
 318     
 319     LCCollectionVec * clusterColl = new LCCollectionVec(LCIO::CLUSTER);
 320     
 321     //set everything to 0
 322     _E_estimated = 0.0;
 323     
 324     _Esum_ecal   = 0.0;
 325     _Esum_hcal   = 0.0;
 326     _Esum_calos  = 0.0;
 327         
 328     _Et_miss = 0.0;
 329     _Et_x    = 0.0;
 330     _Et_y    = 0.0;
 331     
 332     _phiDegree   = 0;
 333     _thetaDegree = 0;
 334     _dca         = 0;
 335     _nPhotons    = 0;
 336       
 337   
 338     unsigned int nHcalElements;
 339     unsigned int nEcalElements;
 340     if ( ecalColl!=0 ) {
 341       if ( _debugLevel > 0 )
 342 	std::cout<<"KITdev found collection "<<_Ecal_col.c_str()<<std::endl;
 343       
 344       nEcalElements = ecalColl->getNumberOfElements();
 345       
 346       //get ECAL hits and energy in this event
 347       _NhitsEcal = nEcalElements;
 348 	
 349      
 350       for (unsigned int i=0; i<nEcalElements; i++) {
 351 	CalorimeterHit* ecalHit = dynamic_cast<CalorimeterHit*>(ecalColl->getElementAt(i));
 352 
 353 	unsigned int layer = CDECAL(ecalHit)["K-1"];
 354 	
 355 	_Ehit_ecal   = ecalHit->getEnergy();
 356 	_Esum_ecal  += ecalHit->getEnergy();
 357 	_Esum_calos += ecalHit->getEnergy();
 358 	
 359 	//calculate transverse energy in ECAL
 360 	float normalisation = sqrt ( pow(ecalHit->getPosition()[0],2) + pow(ecalHit->getPosition()[1],2) + pow(ecalHit->getPosition()[2],2));
 361 	_Et_x += ecalHit->getEnergy() * ecalHit->getPosition()[0] / normalisation; 
 362 	_Et_y += ecalHit->getEnergy() * ecalHit->getPosition()[1] / normalisation;
 363 	
 364 	if ( layer <= myPGDB[myPGdb::ECAL1_BAR].max_lay ) {
 365 	  
 366 	  _Ehit_ecal_mip  = _Ehit_ecal / myPGDB[myPGdb::ECAL1_BAR].mip_whole; //mip_whole=mip_vis*e_coeff;
 367 	  
 368 	  if ( _debugLevel > 8 )
 369 	    std::cout<<myPGDB[myPGdb::ECAL1_BAR].mip_vis<<" * "<<myPGDB[myPGdb::ECAL1_BAR].e_coeff<<"= mip_whole SF1: "<<myPGDB[myPGdb::ECAL1_BAR].mip_whole<<std::endl;
 370 	  
 371 	}
 372 	else {
 373 	  _Ehit_ecal_mip  = _Ehit_ecal / myPGDB[myPGdb::ECAL2_BAR].mip_whole;
 374 	  
 375 	  if ( _debugLevel > 8 )
 376 	    std::cout<<myPGDB[myPGdb::ECAL2_BAR].mip_vis<<" * "<<myPGDB[myPGdb::ECAL2_BAR].e_coeff<<"= mip_whole SF1: "<<myPGDB[myPGdb::ECAL2_BAR].mip_whole<<std::endl;
 377 	  
 378 	}
 379 	
 380       }//all ECAL elem.
 381       if ( _debugLevel > 0 )
 382 	std::cout<<" with "<<_NhitsEcal<<" hits and "<<_Esum_ecal<<"[GeV] energy"<<std::endl;
 383       
 384       //get HCAL hits and energy in this event
 385       if ( hcalColl!=0 ) {
 386 	if ( _debugLevel > 0 )
 387 	  std::cout<<"KITdev found collection "<<_Hcal_col.c_str()<<std::endl;
 388 	
 389 	nHcalElements = hcalColl->getNumberOfElements(); 
 390 	_NhitsHcal = nHcalElements;
 391 	  
 392 	for (unsigned int i=0; i<nHcalElements; i++) {
 393 	  CalorimeterHit* hcalHit = dynamic_cast<CalorimeterHit*>(hcalColl->getElementAt(i));
 394 	    
 395 	  _Esum_hcal += hcalHit->getEnergy();
 396 	  _Ehit_hcal = hcalHit->getEnergy();
 397 	  _Esum_calos += hcalHit->getEnergy();
 398 	   
 399 	  //calculate transverse energy in HCAL
 400 	  float normalisation = sqrt ( pow(hcalHit->getPosition()[0],2) + pow(hcalHit->getPosition()[1],2) + pow(hcalHit->getPosition()[2],2));
 401 	  _Et_x += hcalHit->getEnergy() * hcalHit->getPosition()[0] / normalisation; 
 402 	  _Et_y += hcalHit->getEnergy() * hcalHit->getPosition()[1] / normalisation;
 403 	}
 404 	_NhitsCalos = nHcalElements + nEcalElements;	
 405 
 406 	if ( _debugLevel > 0 )
 407 	  std::cout<<" with "<<_NhitsHcal<<" hits and "<<_Esum_hcal<<"[GeV] energy"<<std::endl;
 408       }//hcal col != 0
 409       
 410       //since energy conservation should be valid, the total transverse energy should be 0 and thus Et = - Et_miss
 411       _Et_miss = sqrt ( pow(_Et_x,2) + pow(_Et_y,2) );      
 412     
 413       //save parameters for later processors
 414       evt->parameters().setValue("_Et_miss",(float)_Et_miss);
 415       evt->parameters().setValue("_Esum_ECAL",(float)_Esum_ecal);
 416       evt->parameters().setValue("_Nhits_ECAL",(int)nEcalElements);
 417       evt->parameters().setValue("_Esum_HCAL",(float)_Esum_hcal);
 418       evt->parameters().setValue("_Nhits_HCAL",(int)nHcalElements);
 419       evt->parameters().setValue("_Esum_calos",(float)_Esum_calos);
 420       evt->parameters().setValue("_Nhits_calos",(int)_NhitsCalos);
 421       evt->parameters().setValue("_trueNcluster",(int)_Ntrcl);
 422 
 423       vector<Superhit2*> calo[10];  // array with vectors of superhit pointers
 424       // creating all superhits, the vector in calo[0] is filled with the superhits; applies mip calibration
 425       CreateAllShits2(ecalColl,CDECAL,calo,_debugLevel); 
 426       
 427       // precalculation, shifting cells to equal distance for NN finder -> KITutil
 428       TotalPrecalc2(calo,CDECAL,nEcalElements,_CleanCut,_debugLevel);  
 429 
 430 
 431       // setting the parameters of the alghorithm
 432       vector <PROTSEED2> photonSeed;
 433       CoreCut Ccut;
 434       Ccut.rCut        = _rCut;
 435       Ccut.distCut     = _distCut;
 436       Ccut.cosCut      = _cosCut;
 437       Ccut.MinHit0     = (unsigned int) _MinHit0;
 438       Ccut.MinHitSplit = (unsigned int) _MinHitSplit;
 439 	
 440       //temporary storage of cluster candidates
 441       const int myNLevels  = _Nlevels;
 442       vector<Tmpcl2*> clusterVector[myNLevels];
 443 	
 444       // apply topological cleaning 
 445       // -> remove hits with < _CleanCut neighbouring hits
 446       if ( _ToClean=="YES" || _ToClean=="yes" || _ToClean=="Yes") {
 447 	if ( _debugLevel > 0 )
 448 	  std::cout << " topological cleaning " << std::endl;
 449 	
 450 	//main loop: level wise clustering -> KITutil 
 451 	FindCores2( &(calo[4]), clusterVector, &photonSeed, _Nlevels, _mipstep, Ccut, _debugLevel );
 452       } else {
 453 	FindCores2( &(calo[0]), clusterVector, &photonSeed, _Nlevels, _mipstep, Ccut, _debugLevel );
 454       }
 455       
 456       // container to store information to assign hits to photon candidates
 457       std::vector<ECALHitWithAttributes> ECALHitsWithAttributes;
 458 	
 459       //calibration to get from measured cluster energy to estimated photon energy
 460       //is hardcoded for different levels and energies in KITutil.cc
 461       std::vector<CoreCalib2> coreCalibration;
 462       CreateCalibration_LDCPrime02Sc(&coreCalibration);
 463       //CreateCalibration_LDC00(&coreCalibration);
 464       int nClusters = 0;
 465       for ( unsigned int seedCounter=0; seedCounter<photonSeed.size(); seedCounter++) {
 466 	if ( photonSeed[seedCounter].active==true ) {
 467 	  if ( _debugLevel > 1 )
 468 	    std::cout<<"KIT found photon seed"<<std::endl;
 469 
 470 	  //find cluster center
 471 	  double centerPosition[3];
 472 	  double directionOfCluster[3];
 473 	  centerPosition[0] = photonSeed[seedCounter].cl->getCenter()[0];
 474 	  centerPosition[1] = photonSeed[seedCounter].cl->getCenter()[1];
 475 	  centerPosition[2] = photonSeed[seedCounter].cl->getCenter()[2];
 476 	  
 477 	  // already normalised cluster direction
 478 	  photonSeed[seedCounter].cl->findInertia();
 479 	  directionOfCluster[0] = photonSeed[seedCounter].cl->direction[0]; 
 480 	  directionOfCluster[1] = photonSeed[seedCounter].cl->direction[1];
 481 	  directionOfCluster[2] = photonSeed[seedCounter].cl->direction[2];
 482 	    
 483 	  double coreEnergy = 0.0;
 484 	  double startPosition[3];
 485 	  double startPositionMinProjection[3];
 486 	  double startPositionMaxProjection[3];
 487 	  double minProjectionAlongMainPrincipalAxis = DBL_MAX;
 488 	  double maxProjectionAlongMainPrincipalAxis = 0.0;
 489 
 490 	  //sum over hit in this cluster
 491 	  for ( unsigned int hitCounter=0; hitCounter<photonSeed[seedCounter].cl->hits.size(); hitCounter++ ) {
 492 	    
 493 	    //sum up core energy
 494 	    coreEnergy += photonSeed[seedCounter].cl->hits[hitCounter]->chit->getEnergy();
 495 	    
 496 	    //get cluster geometry
 497 	    double vCenterToHit[3];		
 498 	    double projectionAlongMainPrincipalAxis = 0.0;
 499 	    
 500 	    for ( unsigned int positionCounter=0; positionCounter<3; positionCounter++ ) {
 501 	      
 502 	      vCenterToHit[positionCounter] = centerPosition[positionCounter] //vector from center to hit
 503 		- photonSeed[seedCounter].cl->hits[hitCounter]->chit->getPosition()[positionCounter];
 504 	      
 505 	      projectionAlongMainPrincipalAxis += directionOfCluster[positionCounter]*vCenterToHit[positionCounter];
 506 	    }//all positions
 507 	    
 508 	    if ( projectionAlongMainPrincipalAxis < minProjectionAlongMainPrincipalAxis ) {
 509 	      
 510 	      minProjectionAlongMainPrincipalAxis = projectionAlongMainPrincipalAxis;
 511 	      
 512 	      for ( unsigned int positionCounter=0; positionCounter<3; positionCounter++ ) {
 513 		startPositionMinProjection[positionCounter] = centerPosition[positionCounter] 
 514 		  - 2.0*minProjectionAlongMainPrincipalAxis*directionOfCluster[positionCounter];
 515 	      }
 516 	    }
 517 		
 518 	    if ( projectionAlongMainPrincipalAxis > maxProjectionAlongMainPrincipalAxis ) {
 519 	      
 520 	      maxProjectionAlongMainPrincipalAxis = projectionAlongMainPrincipalAxis;
 521 	      for ( unsigned int positionCounter=0; positionCounter<3; positionCounter++ ) { 
 522 		startPositionMaxProjection[positionCounter] = centerPosition[positionCounter] 
 523 		  - 2.0*maxProjectionAlongMainPrincipalAxis*directionOfCluster[positionCounter];
 524 	      }
 525 	    }
 526 
 527 	    //get cluster start
 528 	    double rStartPositionMinProjection = sqrt( pow(startPositionMinProjection[0],2) 
 529 						       + pow(startPositionMinProjection[1],2) 
 530 						       + pow(startPositionMinProjection[2],2));
 531 	    double rStartPositionMaxProjection = sqrt( pow(startPositionMaxProjection[0],2) 
 532 						       + pow(startPositionMaxProjection[1],2) 
 533 							 + pow(startPositionMaxProjection[2],2));
 534 	      
 535 	    if ( rStartPositionMinProjection < rStartPositionMaxProjection ) {
 536 	      
 537 	      startPosition[0] = startPositionMinProjection[0];
 538 	      startPosition[1] = startPositionMinProjection[1];
 539 	      startPosition[2] = startPositionMinProjection[2];
 540 	    }
 541 	    else {
 542 
 543 	      startPosition[0] = startPositionMaxProjection[0];
 544 	      startPosition[1] = startPositionMaxProjection[1];
 545 	      startPosition[2] = startPositionMaxProjection[2];
 546 	    }
 547 	  }//all hits in cluster
 548 	    
 549 	  _xStart_cluster = startPosition[0];
 550 	  _yStart_cluster = startPosition[1];
 551 	  _zStart_cluster = startPosition[2];
 552 	  
 553 	  //estimate photon energy depending on the cluster energy, level and detector calibration
 554 	  
 555 	  double photonEnergy = giveMeEEstimate2(photonSeed[seedCounter].level,photonSeed[seedCounter].cl->getEnergy(),coreCalibration);
 556 
 557 	  if ( photonEnergy > 0.0 ) {	
 558 	    
 559 	    Photon2* photonFinder = new Photon2(photonEnergy,centerPosition,startPosition);
 560 	    
 561 	    // loop over all ECAL entries, to collect remaining hits
 562 	    for ( unsigned int ecalHitCounter = 0; ecalHitCounter < nEcalElements; ++ecalHitCounter) {
 563 	      
 564 	      double probabilityAndDistance[2];
 565 	      CalorimeterHit* ECALHit = dynamic_cast<CalorimeterHit*>(ecalColl->getElementAt(ecalHitCounter));
 566 	      
 567 	      photonFinder->Prob(ECALHit,_probabilityDensityCut,probabilityAndDistance);
 568 		
 569 	      bool isECALHitAlreadyAssigned = false;
 570 	      int indexOfECALHitAlreadyAssigned = 0;
 571 
 572 	      for ( unsigned int k = 0; k < ECALHitsWithAttributes.size(); ++k) {
 573 		if ( ECALHitsWithAttributes.at(k).ECALHit == ECALHit) {
 574 		  
 575 		  isECALHitAlreadyAssigned = true;
 576 		  indexOfECALHitAlreadyAssigned = k;
 577 		  break;
 578 		}  
 579 	      }
 580 	      //avoid double counting
 581 	      if ( !isECALHitAlreadyAssigned ) {		    
 582 		if ( probabilityAndDistance[0] > 0.00001 ) {
 583 		  
 584 		  ECALHitWithAttributes hitWithAttributes;
 585 		  
 586 		  hitWithAttributes.ECALHit = ECALHit;
 587 		  hitWithAttributes.relatedCores.push_back(&(photonSeed[seedCounter]));
 588 		  hitWithAttributes.probabilitiesForThisECALHit.push_back(probabilityAndDistance[0]);
 589 		  hitWithAttributes.distancesToCoresForThisECALHit.push_back(probabilityAndDistance[1]);
 590 		  hitWithAttributes.estimatedEnergyPerCore.push_back(photonEnergy);
 591 		  
 592 		  ECALHitsWithAttributes.push_back(hitWithAttributes);
 593 		}
 594 	      }
 595 	      else if ( probabilityAndDistance[0] > 0.00001 ) { 
 596 		  
 597 		ECALHitsWithAttributes.at(indexOfECALHitAlreadyAssigned).relatedCores.push_back(&(photonSeed[seedCounter]));
 598 		ECALHitsWithAttributes.at(indexOfECALHitAlreadyAssigned).probabilitiesForThisECALHit.push_back(probabilityAndDistance[0]);
 599 		ECALHitsWithAttributes.at(indexOfECALHitAlreadyAssigned).distancesToCoresForThisECALHit.push_back(probabilityAndDistance[1]);
 600 		ECALHitsWithAttributes.at(indexOfECALHitAlreadyAssigned).estimatedEnergyPerCore.push_back(photonEnergy);
 601 		
 602 		if ( _debugLevel > 4 ) 
 603 		  std::cout<<" All hits assigned "<<std::endl;	  
 604 	      
 605 	      }
 606 	    }//ecalHitCounter	
 607 	    
 608 	    estimatedE.push_back(photonEnergy);
 609 	    _nPhotons++;
 610 			
 611 	    if ( _debugLevel > 0 ) 
 612 	    std::cout<<" found photon with energy "<<photonEnergy<<std::endl;
 613 	    delete photonFinder;
 614 	    photonFinder = 0;
 615 	  }//photon energy > 0
 616 	  else
 617 	    std::cout<<" energy == 0 -> no photon"<<std::endl;	  
 618 	}//photonSeed is active
 619       }//all cluster
 620 
 621       //******* Collect the remaining hits ********//
 622 	
 623       //map hit to clusters
 624       std::map<PROTSEED2*,ClusterImpl*> mapCoreToCluster;
 625       std::map<PROTSEED2*,double> mapCoreToEstimatedEnergy;
 626 
 627       for ( unsigned int selectedHitsCounter = 0; selectedHitsCounter < ECALHitsWithAttributes.size(); ++selectedHitsCounter ) {
 628 	  
 629 	PROTSEED2* coreToAddHitTo;
 630 	double estimatedEnergyOfCoreToAddHitTo = 0.0;
 631 	    
 632 
 633 	  if ( ECALHitsWithAttributes.at(selectedHitsCounter).relatedCores.size() == 1 ) {
 634 	    
 635 	    coreToAddHitTo = ECALHitsWithAttributes.at(selectedHitsCounter).relatedCores.at(0);	
 636 	    estimatedEnergyOfCoreToAddHitTo = ECALHitsWithAttributes.at(selectedHitsCounter).estimatedEnergyPerCore.at(0);
 637 	  } else {
 638 	      
 639 	    unsigned int indexOfMaxProbability = 0;
 640 	    double maxProbability = 0.0;
 641 	    
 642 	    for ( unsigned int selectedCoresCounter = 0; selectedCoresCounter  <  ECALHitsWithAttributes.at(selectedHitsCounter).relatedCores.size(); ++selectedCoresCounter) {
 643 	      if ( ECALHitsWithAttributes.at(selectedHitsCounter).probabilitiesForThisECALHit.at(selectedCoresCounter) > maxProbability ) {
 644 	      
 645 		indexOfMaxProbability = selectedCoresCounter;
 646 		maxProbability =
 647 		  ECALHitsWithAttributes.at(selectedHitsCounter).probabilitiesForThisECALHit.at(selectedCoresCounter); 
 648 	      }
 649 	    }//selected cores
 650 	  
 651 	    coreToAddHitTo = ECALHitsWithAttributes.at(selectedHitsCounter).relatedCores.at(indexOfMaxProbability);
 652 	    estimatedEnergyOfCoreToAddHitTo = 
 653 	      ECALHitsWithAttributes.at(selectedHitsCounter).estimatedEnergyPerCore.at(indexOfMaxProbability); 
 654 	  }
 655 	    std::map<PROTSEED2*,ClusterImpl*>::iterator position = mapCoreToCluster.find(coreToAddHitTo);
 656 	    
 657 	    if ( position ==  mapCoreToCluster.end() ) { 
 658 	      
 659 	      mapCoreToCluster[coreToAddHitTo] = new ClusterImpl();
 660 	      mapCoreToEstimatedEnergy[coreToAddHitTo] = estimatedEnergyOfCoreToAddHitTo;
 661 	      position = mapCoreToCluster.find(coreToAddHitTo);
 662 	      
 663 	      if ( _debugLevel > 4 ) 
 664 		std::cout<<" Mapped Cores to Clusters "<<std::endl;	  
 665 	    }	
 666 	    
 667 	    (*position).second->addHit(ECALHitsWithAttributes.at(selectedHitsCounter).ECALHit,(float)1.0); 
 668 	}//selected ecal hits
 669 	
 670 
 671 	if ( _debugLevel > 0 ) 
 672 	  std::cout<<" Core to Cluster Map size "<<mapCoreToCluster.size()<<std::endl;	  
 673 
 674 	int iteratorCounter = 0;
 675 
 676 	for ( std::map<PROTSEED2*,ClusterImpl*>::iterator clusterIterator = mapCoreToCluster.begin(); clusterIterator != mapCoreToCluster.end(); clusterIterator++ ) {
 677        	    
 678 	  ClusterImpl* cluster = (*clusterIterator).second;	  
 679 	  double estimatedEnergyOfCluster = mapCoreToEstimatedEnergy[(*clusterIterator).first];
 680 	    
 681 	  // set energy and position of cluster
 682 	  unsigned int caloHits = cluster->getCalorimeterHits().size();
 683 	    
 684 	  float* xPos      = new float[caloHits];
 685 	  float* yPos      = new float[caloHits];
 686 	  float* zPos      = new float[caloHits];
 687 	  float* hitEnergy = new float[caloHits];
 688 	  float energy     = 0.0;
 689 	  
 690 	  for ( unsigned int caloHitsCounter = 0; caloHitsCounter < caloHits; ++caloHitsCounter ) {
 691 	      
 692 	    xPos[caloHitsCounter] = cluster->getCalorimeterHits().at(caloHitsCounter)->getPosition()[0];
 693 	    yPos[caloHitsCounter] = cluster->getCalorimeterHits().at(caloHitsCounter)->getPosition()[1];
 694 	    zPos[caloHitsCounter] = cluster->getCalorimeterHits().at(caloHitsCounter)->getPosition()[2];
 695 	    hitEnergy[caloHitsCounter] = cluster->getCalorimeterHits().at(caloHitsCounter)->getEnergy();
 696 	    
 697 	    energy += cluster->getCalorimeterHits().at(caloHitsCounter)->getEnergy(); 
 698 	  }
 699 	
 700 	  // cut of EM candidtate clusters which differ more than _energyDeviationCut from enery estimate
 701 	  double energyDeviation = fabs( (energy-estimatedEnergyOfCluster)/estimatedEnergyOfCluster ); 
 702 	  //FIXME: This accounts only for the leading cluster!
 703 	  
 704 	  // simple cut on the starting layer of the shower
 705 	  // FIXME: this is not working properly for the edges of the ECAL and HCAL => need pseudo layer !!!
 706 	  const unsigned int cutOnLayer = 10; // FIXME: remove hard-coded cut	
 707 	  const unsigned int contOnNumberOfHitsInLayerCut = 4; // FIXME: remove hard-coded cut
 708 	  unsigned int nHitsInLayerCut = 0;
 709 	  bool hitWithinCutRangeFound  = false;
 710 	  
 711 	  //	  double energyDeviation = 0;
 712 	  
 713 	  
 714 	  _E_estimated = estimatedE.at(iteratorCounter);
 715 	  //	  energyDeviation = fabs( (energy - _E_estimated) / _E_estimated );
 716 	  
 717 	  if (_debugLevel > 0) {
 718 	    std::cout<<"energyDeviation:  (" <<energy<<" - "<<_E_estimated<<") / "<<_E_estimated
 719 		     <<" = "<<energyDeviation<<std::endl;
 720 	  }
 721 	  
 722 	  for ( unsigned int caloHitsCounter = 0; caloHitsCounter < cluster->getCalorimeterHits().size(); ++caloHitsCounter ) {
 723 	    
 724 	    unsigned int layer = CDECAL(cluster->getCalorimeterHits().at(caloHitsCounter))["K-1"];
 725 	    if ( layer <= cutOnLayer ) {
 726 	      ++nHitsInLayerCut;
 727 	    }
 728 
 729 	    if ( nHitsInLayerCut >= contOnNumberOfHitsInLayerCut ) {//if >= 4 layers hit within the first 10 layers
 730 	      hitWithinCutRangeFound = true;
 731 	      break;
 732 	    }
 733 
 734 	    //COMMENT IN && HITWITHINRANGE
 735 	    if ( energyDeviation < _energyDeviationCut ) {
 736 	      //	      if ( hitWithinCutRangeFound ) {
 737 	      // flag hits 
 738 	      for ( unsigned int j = 0; j < cluster->getCalorimeterHits().size(); ++j ) {
 739 		cluster->getCalorimeterHits().at(j)->ext<isPartOfEMShowerCandidate>() = 1;
 740 		//clusterColl->addHit(cluster->getCalorimeterHits().at(j));
 741 	      }
 742 	      ClusterShapes* clusterShape = new ClusterShapes(caloHits,hitEnergy,xPos,yPos,zPos);
 743 	      
 744 	      float position[3];
 745 	      position[0] = clusterShape->getCentreOfGravity()[0];
 746 	      position[1] = clusterShape->getCentreOfGravity()[1];
 747 	      position[2] = clusterShape->getCentreOfGravity()[2];
 748 	      
 749 	      float normalisation = sqrt( pow(position[0],2) + pow(position[1],2) + pow(position[2],2) );
 750 	      float normPos[3];
 751 	      normPos[0] = position[0]/normalisation;
 752 	      normPos[1] = position[1]/normalisation;
 753 	      normPos[2] = position[2]/normalisation;
 754 
 755 	      float _phiPos = atan2(normPos[1],normPos[0]); 
 756 	      float _thetaPos = acos(normPos[2]);
 757 
 758 	      //get main cluster axis
 759 	      float mainClusterAxis[3];
 760 	      mainClusterAxis[0] = clusterShape->getEigenVecInertia()[0];
 761 	      mainClusterAxis[1] = clusterShape->getEigenVecInertia()[1];
 762 	      mainClusterAxis[2] = clusterShape->getEigenVecInertia()[2];
 763 	      //get cluster angle
 764 	      float phiCluster = atan2(clusterShape->getEigenVecInertia()[1],clusterShape->getEigenVecInertia()[0]);
 765 	      float thetaCluster = acos(clusterShape->getEigenVecInertia()[2]); 
 766 
 767 	      std::vector<float> shapeParameter(3) ;
 768 
 769 	      shapeParameter[0] = mainClusterAxis[0]; 
 770 	      shapeParameter[1] = mainClusterAxis[1]; 
 771 	      shapeParameter[2] = mainClusterAxis[2]; 
 772 
 773 	      cluster->setPosition(position); 
 774 	      cluster->setEnergy(energy);
 775 	      cluster->setITheta(thetaCluster); 
 776 	      cluster->setIPhi(phiCluster);
 777 	      cluster->setShape(shapeParameter);
 778 	      clusterColl->addElement(cluster);
 779 
 780 	      nClusters++;
 781 
 782 	      nClusterHits.push_back(cluster->getCalorimeterHits().size());
 783 	      phiPos.push_back(_phiPos);
 784 	      thetaPos.push_back(_thetaPos);
 785 
 786 	      if(energy == 0.)
 787 		std::cout<<"WARNING: found cluster with 0 GeV energy!"<<std::endl;
 788 
 789 
 790 	      _phiDegree = phiCluster*180/M_PI;
 791 	      _thetaDegree = thetaCluster*180/M_PI;
 792 	      
 793 	      //glue cluster direction to cog
 794 	      float photonLine[3];
 795 	      for (int axisCounter=0; axisCounter<3; axisCounter++) {
 796 		photonLine[axisCounter] = position[axisCounter] + mainClusterAxis[axisCounter]; //g=cog+k*EV
 797 	      }
 798 	      
 799 	      //calculate distance of closest approach between gamma and IP:
 800 	      // d = | cog x EV_1 | 
 801 	      float dcaVec[3];
 802 	      dcaVec[0] = position[1]*mainClusterAxis[2] - position[2]*mainClusterAxis[1];
 803 	      dcaVec[1] = position[2]*mainClusterAxis[0] - position[0]*mainClusterAxis[2];
 804 	      dcaVec[2] = position[0]*mainClusterAxis[1] - position[1]*mainClusterAxis[0];
 805 	      
 806 	      _dca = sqrt(  pow(dcaVec[0],2) + pow(dcaVec[1],2)+ pow(dcaVec[2],2) );
 807 	      
 808 	      closestDistance.push_back(_dca);
 809 
 810 	      if ( _debugLevel > 0 ) {
 811 		std::cout<<" ************************* "<<std::endl;
 812 		std::cout<<"Final cluster with energy: "<<energy
 813 			 <<" at position "<<position[0]
 814 			 <<" , "<<position[1]
 815 			 <<" , "<<position[2]<<std::endl;
 816 		std::cout<<"Cluster direction phi:   "<<_phiDegree<<std::endl;
 817 		std::cout<<"                  theta: "<<_thetaDegree<<std::endl;
 818 		std::cout<<"        main axis: "<<mainClusterAxis[0]
 819 			 <<" , "<<mainClusterAxis[1]
 820 			 <<" , "<<mainClusterAxis[2]<<std::endl;
 821 		std::cout<<"minimal distance from IP:   "<<_dca<<std::endl;
 822 		std::cout<<" ************************* "<<std::endl;
 823 	      }//DEBUG
 824 		
 825 		
 826 	      delete clusterShape;
 827 	      clusterShape = 0; 
 828 	      break;
 829 	      //	      }//hit within range cut
 830 	    //	      else std::cout << "hit not within range cut " << std::endl;
 831 	      
 832 	    }//energy devitation cut
 833 	    else if (_debugLevel >0) std::cout << "energy deviation "<<energyDeviation<<" > " << _energyDeviationCut << std::endl;
 834 	    
 835 
 836 	  }//all calorimeter hits
 837 
 838 	  delete[] xPos;
 839 	  xPos = 0;
 840 	  delete[] yPos;
 841 	  yPos = 0;
 842 	  delete[] zPos;
 843 	  zPos = 0;
 844 	  delete[] hitEnergy;
 845 	  hitEnergy = 0;
 846 	  
 847 	  iteratorCounter++;
 848 	}//in cluster map 
 849 	
 850 	//save parameters for rootTreeWriter
 851 	clusterColl->parameters().setValues("_dca",closestDistance);    
 852 	clusterColl->parameters().setValues("_nClusterHits",nClusterHits);
 853 	clusterColl->parameters().setValues("_phiPos",phiPos);
 854 	clusterColl->parameters().setValues("_thetaPos",thetaPos);
 855 	clusterColl->parameters().setValues("_trueDist",trueDistance);
 856 	clusterColl->parameters().setValues("_trueGammaEnergy",trueGammaEnergy);
 857 	clusterColl->parameters().setValues("_trueGammaPosX",trueGammaPosX);
 858 	clusterColl->parameters().setValues("_trueGammaPosY",trueGammaPosY);
 859 	clusterColl->parameters().setValues("_trueGammaPosZ",trueGammaPosZ);
 860 	clusterColl->parameters().setValues("_trueGammaDirX",trueGammaDirX);
 861 	clusterColl->parameters().setValues("_trueGammaDirY",trueGammaDirY);
 862 	clusterColl->parameters().setValues("_trueGammaDirZ",trueGammaDirZ);
 863 	clusterColl->parameters().setValues("_trueNeutralinoEnergy",trueNeutralinoEnergy);
 864 	clusterColl->parameters().setValues("_trueNeutralinoDirX",trueNeutralinoDirX);
 865 	clusterColl->parameters().setValues("_trueNeutralinoDirY",trueNeutralinoDirY);
 866 	clusterColl->parameters().setValues("_trueNeutralinoDirZ",trueNeutralinoDirZ);
 867 	clusterColl->parameters().setValues("_trueThetaDir",trueThetaDir);
 868 	clusterColl->parameters().setValues("_truePhiDir",truePhiDir);
 869 
 870 	mapCoreToCluster.clear();
 871 	ECALHitsWithAttributes.clear();
 872 	estimatedE.clear();
 873 
 874 	  // for strong memory and nice dreams clean up again
 875 	for ( int levelCounter=0; levelCounter<_Nlevels; levelCounter++ ) {
 876 	  if ( clusterVector[levelCounter].size()!=0 )
 877 	    for ( unsigned int calibrationCounter=0; calibrationCounter<clusterVector[levelCounter].size(); calibrationCounter++ ) {
 878 	      delete clusterVector[levelCounter][calibrationCounter];
 879 	    }
 880 	}
 881 	  
 882 	for ( unsigned int caloCounter=0; caloCounter<2; caloCounter++ ) {        
 883 	  if ( calo[caloCounter].size()!=0 )
 884 	    for ( unsigned int caloSize=0; caloSize<calo[caloCounter].size(); caloSize++ )
 885 	      delete (calo[caloCounter][caloSize]);
 886 	}
 887       }//ecalCollection not empty      
 888 	
 889       evt->addCollection(clusterColl,_CoreCollection.c_str());
 890       evt->parameters().setValue("_nGamma",(int)_nPhotons);
 891 
 892   } catch ( DataNotAvailableException &e ) {
 893     std::cout << "KITdev:: missing collection in event " << _nEvt << std::endl; 
 894   }
 895   
 896   _nEvt++ ;
 897 
 898 }//process event
 899 
 900 void KITdev::end(){
 901 
 902 }

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.