Attachment 'KITcalibrator.cc'

Download

   1 #include "KITcalibrator.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 KITcalibrator aKITcalibrator ;
  11 KITcalibrator::KITcalibrator() : Processor("KITcalibrator") {
  12   
  13   // modify processor description
  14   _description = "KITcalibrator does whatever it does ..." ;
  15   
  16   registerProcessorParameter( "ECAL_Collection",
  17 			      "Name of the CalorimeterHit Collection for ECAL ",
  18 			      _Ecal_col,
  19 			      std::string("ECAL"));
  20   registerProcessorParameter( "Core_Collection",
  21 			      "Name of the Cluster Collection for Cores ",
  22 			      _CoreCollection,
  23 			      std::string("CORE"));
  24 
  25   registerProcessorParameter( "Cleaning",
  26 			      "To do the cleaning on hits or not ",
  27 			      _ToClean,
  28 			      std::string("YES"));
  29   registerProcessorParameter( "TopologicalCut",
  30 			      "At which number of neighbors to put the threshold, condition is < so you need to put N+1 ",
  31 			       _CleanCut,
  32 			      (int)5);
  33 
  34   registerProcessorParameter( "NumberOfLevels",
  35 			      "Number of levels for central loop ",
  36 			       _Nlevels,
  37 			      (int)10);
  38   vector<float> mipstep;
  39   mipstep.push_back(0.1);
  40   mipstep.push_back(1.5);
  41   mipstep.push_back(2.5);
  42   mipstep.push_back(4.0);
  43   mipstep.push_back(6.0);
  44   mipstep.push_back(9.0);
  45   mipstep.push_back(16.0);
  46   mipstep.push_back(26.0);
  47   mipstep.push_back(41.0);
  48   mipstep.push_back(65.0);
  49 
  50   registerProcessorParameter( "Levels",
  51 			      "Levels for central loop in MIP ",
  52 			      _mipstep,
  53 			      mipstep);
  54 
  55   registerProcessorParameter( "MinHit0",
  56 			      "Minimal Number of hits for ground level cluster ",
  57 			      _MinHit0,
  58 			      (int)4);
  59 
  60   registerProcessorParameter( "MinHitSplit",
  61 			      "Minimal Number of hits for i-th level cluster ",
  62 			      _MinHitSplit,
  63 			      (int)2);
  64 
  65   registerProcessorParameter( "Rcut",
  66 			      "Fluctuation suprresion cut",
  67 			      _rCut,
  68 			      (double)0.4);
  69   registerProcessorParameter( "Distcut",
  70 			      "Square of distance cut for merging ",
  71 			      _distCut,
  72 			      (double)35.0);
  73   registerProcessorParameter( "Coscut",
  74 			      "Cosine of the angle for merging ",
  75 			      _cosCut,
  76 			      (double)0.95);
  77   
  78   registerProcessorParameter( "ROOT_NAME",
  79 			      "Name of the root output file ",
  80 			      _rootTreeName,
  81 			      std::string("KITcalibrator"));
  82 
  83   registerProcessorParameter( "DebugLevel",
  84 			      "limits the amount of information written to std out (0 - none, 9 - maximal information)",
  85 			      _debugLevel,
  86 			      int(0) );
  87   
  88   registerProcessorParameter( "Ncluster",
  89 			      "expected N cluster from cluster cheater",
  90 			      _desiredNCluster,
  91 			      int(1) );
  92 }
  93 
  94 
  95 void KITcalibrator::init() { 
  96      
  97   printParameters() ;
  98   
  99   _rootTreeName = _rootTreeName+".root";
 100   _rootFile = new TFile( _rootTreeName.c_str(), "RECREATE");
 101   _myTree   = new TTree("tree","statistics");
 102   
 103   _myTree->Branch("_N_trueCluster",&_NtrueCluster,"_N_trueCluster/I");
 104   _myTree->Branch("_true_clusterEnergy",  &_trueClusterEnergy,"_true_clusterEnergy/D");
 105   _myTree->Branch("_energy_level0_cluster",&_level0_ClusterEnergy,"_energy_level0_cluster/D");
 106   _myTree->Branch("_energy_level1_cluster",&_level1_ClusterEnergy,"_energy_level1_cluster/D");
 107   _myTree->Branch("_energy_level2_cluster",&_level2_ClusterEnergy,"_energy_level2_cluster/D");
 108   _myTree->Branch("_energy_level3_cluster",&_level3_ClusterEnergy,"_energy_level3_cluster/D");
 109   _myTree->Branch("_energy_level4_cluster",&_level4_ClusterEnergy,"_energy_level4_cluster/D");
 110   _myTree->Branch("_energy_level5_cluster",&_level5_ClusterEnergy,"_energy_level5_cluster/D");
 111   _myTree->Branch("_energy_level6_cluster",&_level6_ClusterEnergy,"_energy_level6_cluster/D");
 112   _myTree->Branch("_energy_level7_cluster",&_level7_ClusterEnergy,"_energy_level7_cluster/D");
 113   _myTree->Branch("_energy_level8_cluster",&_level8_ClusterEnergy,"_energy_level8_cluster/D");
 114   _myTree->Branch("_energy_level9_cluster",&_level9_ClusterEnergy,"_energy_level9_cluster/D");
 115 	 
 116   
 117   CreateCalibration_LDCPrime02Sc(&coreCalibVector);
 118 
 119   MarlinCED::init(this) ;
 120   _nRun = 0 ;
 121   _nEvt = 0 ;   
 122 }
 123 
 124 void KITcalibrator::processRunHeader( LCRunHeader* run) { 
 125       
 126   _nRun++ ;
 127 } 
 128 
 129 void KITcalibrator::processEvent( LCEvent * evt ) {  
 130   _NtrueCluster = 0;
 131   _trueClusterEnergy      = 0.0;
 132   _level0_ClusterEnergy   = 0.0;
 133   _level1_ClusterEnergy   = 0.0;
 134   _level2_ClusterEnergy   = 0.0;
 135   _level3_ClusterEnergy   = 0.0;
 136   _level4_ClusterEnergy   = 0.0;
 137   _level5_ClusterEnergy   = 0.0;
 138   _level6_ClusterEnergy   = 0.0;
 139   _level7_ClusterEnergy   = 0.0;
 140   _level8_ClusterEnergy   = 0.0;
 141   _level9_ClusterEnergy   = 0.0;
 142 
 143 
 144   LCCollection* colm = evt->getCollection("MCParticle");
 145   MCParticle* mc = dynamic_cast<MCParticle*>(colm->getElementAt(0));
 146   
 147   myPGdb::ZONE zon ;
 148   double zone[3];
 149   zone[0] = mc->getEndpoint()[0];  
 150   zone[1] = mc->getEndpoint()[1];  
 151   zone[2] = mc->getEndpoint()[2];
 152   
 153   Point3D dsf2(zone[0], zone[1], zone[2]);
 154   zon = myPGDB.get_zone(dsf2);
 155 
 156   bool nijeivica = true;
 157   if ( zon==myPGdb::ECAL1_BAR || zon==myPGdb::ECAL2_BAR ) { //if in ecal barrel
 158     if(myPGDB[myPGdb::ECAL1_BAR].z_outer - fabs(zone[2]) < 200.0)
 159       nijeivica = false;
 160   }
 161   
 162   //check availability of geometrical information
 163   if ( (zon==myPGdb::ECAL1_BAR || zon==myPGdb::ECAL2_BAR || zon==myPGdb::ECAL1_CAP || zon==myPGdb::ECAL2_CAP) && nijeivica ) { //
 164     try {  
 165       
 166       LCCollection* ecalCol = evt->getCollection(_Ecal_col.c_str()) ;
 167       LCCollection* true_col = evt->getCollection("TrueClusters");
 168 
 169       CellIDDecoder<CalorimeterHit>CDECAL = CellIDDecoder<CalorimeterHit>(ecalCol);
 170       
 171       _NtrueCluster = 0;
 172       if ( true_col!=0 ) {
 173 	//_NtrueCluster = true_col->getNumberOfElements();
 174 	
 175 	for (int colCounter=0; colCounter<true_col->getNumberOfElements(); colCounter++) {
 176 	  
 177 	  Cluster* cl = dynamic_cast<Cluster*>(true_col->getElementAt(colCounter));
 178 	  _trueClusterEnergy = cl->getEnergy();
 179 	  if ( cl->getEnergy()>0.2 )
 180 	    _NtrueCluster++;
 181 	}
 182       }
 183 
 184       if ( _debugLevel >0 )
 185 	std::cout<<"KIT calibrator found "<<_NtrueCluster<<" TrueClusters"<<std::endl;
 186     
 187       if ( _NtrueCluster == _desiredNCluster && ecalCol != 0 ) {
 188 	
 189 	unsigned int nelem = ecalCol->getNumberOfElements(); 
 190 	// array with vectors of superhit pointers
 191 	vector<mySuperhit*> calo[_Nlevels]; 
 192 	// creating all superhits, the vector in calo[0] is filled with the superhits; applies mip calibration
 193 	CreateAllShits(ecalCol,CDECAL,calo); 
 194 	
 195 	if ( _debugLevel >0 )
 196 	  std::cout<<"KIT calibrator created mySuperhits "<<std::endl;
 197 	
 198 	// shift ECAL hits to equal distances
 199 	TotalPrecalc(calo,CDECAL,nelem,_CleanCut);  
 200 	if ( _debugLevel >0 )
 201 	  std::cout<<"KIT calibrator performed precalculations "<<std::endl;
 202 	
 203 	// setting the parameters of the alghorithm
 204 	vector <PROTSEED> photonSeed;
 205 	CoreCut2 Ccut;
 206 	Ccut.rCut    = _rCut;
 207 	Ccut.distCut = _distCut;  
 208 	Ccut.cosCut  = _cosCut;
 209 	Ccut.MinHit0 = (unsigned int) _MinHit0;
 210 	Ccut.MinHitSplit = (unsigned int) _MinHitSplit;
 211 
 212 	//temporary storage of cluster candidates
 213 	const unsigned int myNLevels = _Nlevels;
 214 	vector<myTmpcl*> clusterVector[myNLevels];
 215 	double Esum[myNLevels];	  
 216 
 217 	// finding cores.
 218 	if ( _ToClean=="YES" || _ToClean=="yes" || _ToClean=="Yes") {
 219 	  if(_debugLevel >0)
 220 	    cout << " topological cleaning " << endl;
 221 	  FindMyCores( &(calo[4]), clusterVector, &photonSeed, _Nlevels, _mipstep, Ccut );
 222 	} else {
 223 	  FindMyCores( &(calo[0]), clusterVector, &photonSeed, _Nlevels, _mipstep, Ccut );
 224 	}
 225 	
 226 	if ( photonSeed.size() != 0 ) {
 227 	  if ( photonSeed[0].active == true ) {
 228 	    if ( _debugLevel >0 )
 229 	      std::cout<<"KIT calibrator found active photon seed"<<std::endl;
 230 	    for ( int levelCounter=0; levelCounter<_Nlevels; levelCounter++ ) { //has to be adjusted if more than 10 cluster levels are required!
 231 	      if ( _debugLevel >2 ) {
 232 		std::cout<<"from "<<nelem<<" ECAL hits: "<<clusterVector[levelCounter].size()
 233 			 <<" are level "<<levelCounter
 234 			 <<" with Ehit >  "<<_mipstep[levelCounter]
 235 			 <<" mip "<<std::endl;
 236 	      }
 237 
 238 	      Esum[levelCounter] = 0; 
 239 	      for(unsigned int vectorCounter=0; vectorCounter<clusterVector[levelCounter].size(); vectorCounter++) {
 240 		//FIXME this energy is way too high >> photonSeed[seedCounter].cl->getEnergy()
 241 		Esum[levelCounter] += clusterVector[levelCounter][vectorCounter]->getEnergy();  
 242 	      }
 243 	      
 244 	      switch(levelCounter) {
 245 	      case 0:
 246 		_level0_ClusterEnergy = Esum[levelCounter];
 247 		_myTree->Fill();
 248 		break;
 249 	      case 1:
 250 		_level1_ClusterEnergy = Esum[levelCounter];
 251 		_myTree->Fill();
 252 		break;
 253 	      case 2:
 254 		_level2_ClusterEnergy = Esum[levelCounter];
 255 		_myTree->Fill();
 256 		break;
 257 	      case 3:
 258 		_level3_ClusterEnergy = Esum[levelCounter];
 259 		_myTree->Fill();
 260 		break;
 261 	      case 4:
 262 		_level4_ClusterEnergy = Esum[levelCounter];
 263 		_myTree->Fill();
 264 		break;
 265 	      case 5:
 266 		_level5_ClusterEnergy = Esum[levelCounter];
 267 		_myTree->Fill();
 268 		break;
 269 	      case 6:
 270 		_level6_ClusterEnergy = Esum[levelCounter];
 271 		_myTree->Fill();
 272 		break;
 273 	      case 7:
 274 		_level7_ClusterEnergy = Esum[levelCounter];
 275 		_myTree->Fill();
 276 		break;
 277 	      case 8:
 278 		_level8_ClusterEnergy = Esum[levelCounter];
 279 		_myTree->Fill();
 280 		break;
 281 	      case 9:
 282 		_level9_ClusterEnergy = Esum[levelCounter];
 283 		_myTree->Fill();
 284 		break;
 285 	      }//switch
 286 	    }//all level  //seedCounter
 287 	  }//photonSeed.active
 288 	}//photonSeed.size=1
 289 
 290 	// for strong memory and nice dreams ..
 291 	for ( unsigned int levelCounter=0; levelCounter<myNLevels; levelCounter++ ) {
 292 	  if ( clusterVector[levelCounter].size()!=0 ) 
 293 	    for ( unsigned int vectorCounter=0; vectorCounter<clusterVector[levelCounter].size(); vectorCounter++ )
 294 	      delete clusterVector[levelCounter][vectorCounter]; 
 295 	}
 296 	
 297 	for ( unsigned int zoneCounter=0; zoneCounter<2; zoneCounter++ ) {        
 298 	  if(calo[zoneCounter].size()!=0)
 299 	    for( unsigned int caloCounter=0; caloCounter<calo[zoneCounter].size(); caloCounter++)
 300 	      delete (calo[zoneCounter])[caloCounter];
 301 	}
 302       }//ecal collection not empty && desired nCluster
 303     } catch ( DataNotAvailableException &e ) {}
 304   }//zone
 305   _nEvt ++ ;
 306 }
 307 
 308 void KITcalibrator::end(){
 309   
 310   _myTree->Write();
 311   delete _myTree;
 312   _rootFile->Write();
 313   _rootFile->Close();
 314   delete _rootFile;
 315 }

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.