Attachment 'ClusterWriteEngine.cc'

Download

   1 #include "KIT_ClusterWriteEngine.hh"
   2 
   3 #include <cfloat>
   4 #include <iostream>
   5 
   6 #include "EVENT/Cluster.h"
   7 
   8 #include <string>
   9 
  10 using namespace lcio;
  11 using namespace std;
  12 using namespace marlin;
  13 
  14 #define DDEBUG(name) std::cout << __FILE__ <<","<<__LINE__ << "; " << #name<<": " << name << std::endl;
  15 #define IDEBUG(name) std::cout << __FILE__ <<","<<__LINE__ << "; " << #name <<" at " << &name << std::endl;
  16 
  17 #define INVALIDF (-FLT_MAX)
  18 #define INVALIDI INT_MIN
  19 
  20 /* This engine can be run within the RootTreeWriter
  21  * it will write all available information for clusters
  22  * it was tested for root 5.16.00
  23  */
  24 namespace marlin
  25 {
  26   void KIT_ClusterWriteEngine::registerParameters()
  27   {
  28     _hostProcessor.relayRegisterProcessorParameter("KIT_ClusterWriteEngine_prefix",
  29 						   "KIT_ClusterWriteEngine prefix to tree variables,e.g. cluster",
  30 						   _prefix,
  31 						   std::string("cluster"));
  32 
  33     _hostProcessor.relayRegisterInputCollection(LCIO::RECONSTRUCTEDPARTICLE,_engineName+"_InCol",
  34 						"Name of input collection",
  35 						_clusterColName, std::string("clusterColl")  );
  36 
  37   }
  38 
  39   void KIT_ClusterWriteEngine::registerBranches( TTree* hostTree )
  40   {
  41     if ( _prefix.size() > 0 )
  42       if ( _prefix[ _prefix.length()-1 ] != '_' )
  43 	_prefix += "_";
  44     //declare all engine branches:
  45     hostTree->Branch(string(_prefix+"iEvt").c_str(), &_clusterFill.iEvt, string(_prefix+"iEvt/I").c_str());
  46     hostTree->Branch(string(_prefix+"nCluster").c_str(), &_clusterFill.nCluster, string(_prefix+"nCluster/I").c_str());
  47     hostTree->Branch(string(_prefix+"nGamma").c_str(), &_clusterFill.nGamma, string(_prefix+"nGamma/I").c_str());
  48     hostTree->Branch(string(_prefix+"energy").c_str(),&_clusterFill.energy, string(_prefix+"energy["+_prefix+"nCluster]/F").c_str());
  49     hostTree->Branch(string(_prefix+"etMiss").c_str(),&_clusterFill.etMiss, string(_prefix+"etMiss/F").c_str());
  50     hostTree->Branch(string("calos_energySum").c_str(),&_clusterFill.esumCalos, string("calos_energySum/F").c_str());
  51     hostTree->Branch(string("calos_nHits").c_str(),&_clusterFill.nHitsCalos, string("calos_nHits/I").c_str());
  52     hostTree->Branch(string("ECAL_energySum").c_str(),&_clusterFill.esumECAL, string("ECAL_energySum/F").c_str());
  53     hostTree->Branch(string("ECAL_nHits").c_str(),&_clusterFill.nHitsECAL, string("ECAL_nHits/I").c_str());
  54     hostTree->Branch(string("HCAL_energySum").c_str(),&_clusterFill.esumHCAL, string("HCAL_energySum/F").c_str());
  55     hostTree->Branch(string("HCAL_nHits").c_str(),&_clusterFill.nHitsHCAL, string("HCAL_nHits/I").c_str());
  56     hostTree->Branch(string(_prefix+"position").c_str(),&_clusterFill.clusterPos, string(_prefix+"position["+_prefix+"nCluster][3]/F").c_str());
  57     hostTree->Branch(string(_prefix+"thetaDir").c_str(),&_clusterFill.thetaDir, string(_prefix+"thetaDir["+_prefix+"nCluster]/F").c_str());
  58     hostTree->Branch(string(_prefix+"phiDir").c_str(),&_clusterFill.phiDir, string(_prefix+"phiDir["+_prefix+"nCluster]/F").c_str());
  59     hostTree->Branch(string(_prefix+"cosThetaPos").c_str(),&_clusterFill.cosThetaPos, string(_prefix+"cosThetaPos["+_prefix+"nCluster]/F").c_str());
  60     hostTree->Branch(string(_prefix+"phiPos").c_str(),&_clusterFill.phiPos, string(_prefix+"phiPos["+_prefix+"nCluster]/F").c_str());
  61     hostTree->Branch(string(_prefix+"dca").c_str(),&_clusterFill.dca, string(_prefix+"dca["+_prefix+"nCluster]/F").c_str());
  62     hostTree->Branch(string(_prefix+"clusterAxis").c_str(),&_clusterFill.clusterAxis, string(_prefix+"clusterAxis["+_prefix+"nCluster][3]/F").c_str());
  63 
  64 
  65     hostTree->Branch(string(_prefix+"nHits").c_str(), &_clusterFill.nHitsPerCluster, string(_prefix+"nHits["+_prefix+"nCluster]/I").c_str());
  66 
  67 
  68     hostTree->Branch(string(_prefix+"trueDist").c_str(), &_clusterFill.trueDist, string(_prefix+"trueDist["+_prefix+"nCluster]/F").c_str());
  69     hostTree->Branch(string(_prefix+"trueGammaEnergy").c_str(), &_clusterFill.trueGammaEnergy, string(_prefix+"trueGammaEnergy["+_prefix+"nCluster]/F").c_str());
  70     hostTree->Branch(string(_prefix+"trueGammaDirX").c_str(), &_clusterFill.trueGammaDirX, string(_prefix+"trueGammaDirX["+_prefix+"nCluster]/F").c_str());
  71     hostTree->Branch(string(_prefix+"trueGammaDirY").c_str(), &_clusterFill.trueGammaDirY, string(_prefix+"trueGammaDirY["+_prefix+"nCluster]/F").c_str());
  72     hostTree->Branch(string(_prefix+"trueGammaDirZ").c_str(), &_clusterFill.trueGammaDirZ, string(_prefix+"trueGammaDirZ["+_prefix+"nCluster]/F").c_str());
  73     hostTree->Branch(string(_prefix+"trueGammaPosX").c_str(), &_clusterFill.trueGammaPosX, string(_prefix+"trueGammaPosX["+_prefix+"nCluster]/F").c_str());
  74     hostTree->Branch(string(_prefix+"trueGammaPosY").c_str(), &_clusterFill.trueGammaPosY, string(_prefix+"trueGammaPosY["+_prefix+"nCluster]/F").c_str());
  75     hostTree->Branch(string(_prefix+"trueGammaPosZ").c_str(), &_clusterFill.trueGammaPosZ, string(_prefix+"trueGammaPosZ["+_prefix+"nCluster]/F").c_str());
  76     hostTree->Branch(string(_prefix+"trueNeutralinoEnergy").c_str(), &_clusterFill.trueNeutralinoEnergy, string(_prefix+"trueNeutralinoEnergy["+_prefix+"nCluster]/F").c_str());
  77     hostTree->Branch(string(_prefix+"trueNeutralinoDirX").c_str(), &_clusterFill.trueNeutralinoDirX, string(_prefix+"trueNeutralinoDirX["+_prefix+"nCluster]/F").c_str());
  78     hostTree->Branch(string(_prefix+"trueNeutralinoDirY").c_str(), &_clusterFill.trueNeutralinoDirY, string(_prefix+"trueNeutralinoDirY["+_prefix+"nCluster]/F").c_str());
  79     hostTree->Branch(string(_prefix+"trueNeutralinoDirZ").c_str(), &_clusterFill.trueNeutralinoDirZ, string(_prefix+"trueNeutralinoDirZ["+_prefix+"nCluster]/F").c_str());
  80     hostTree->Branch(string(_prefix+"trueThetaDir").c_str(), &_clusterFill.trueThetaDir, string(_prefix+"trueThetaDir["+_prefix+"nCluster]/F").c_str());
  81     hostTree->Branch(string(_prefix+"truePhiDir").c_str(), &_clusterFill.truePhiDir, string(_prefix+"truePhiDir["+_prefix+"nCluster]/F").c_str());
  82     hostTree->Branch(string(_prefix+"trueNcluster").c_str(), &_clusterFill.trueNcluster, string(_prefix+"trueNcluster/I").c_str());
  83     // init not yet initialised variables
  84     _ievt = 0;
  85 
  86   }//register branches
  87 
  88   //for each event do:
  89   void KIT_ClusterWriteEngine::fillVariables( EVENT::LCEvent* evt ) {
  90     
  91     LCCollection* clusterCol;
  92 
  93 
  94     try {
  95       clusterCol = evt->getCollection( _clusterColName );     // open cluster collection;
  96       
  97       _nCluster = clusterCol->getNumberOfElements();          //# cluster per event
  98                       
  99       vector<float> dcaVec; 
 100       clusterCol->getParameters().getFloatVals("_dca",dcaVec);
 101       vector<float> clusterAxisVec; 
 102       clusterCol->getParameters().getFloatVals("_clusterAxis",clusterAxisVec);
 103       vector<int> nHitsVec; 
 104       clusterCol->getParameters().getIntVals("_nClusterHits",nHitsVec);
 105       vector<float> phiPosVec; 
 106       clusterCol->getParameters().getFloatVals("_phiPos",phiPosVec);
 107       vector<float> thetaPosVec; 
 108       clusterCol->getParameters().getFloatVals("_thetaPos",thetaPosVec);
 109       vector<float> trueDistVec; 
 110       clusterCol->getParameters().getFloatVals("_trueDist",trueDistVec);
 111       vector<float> trueGammaEnergyVec; 
 112       clusterCol->getParameters().getFloatVals("_trueGammaEnergy",trueGammaEnergyVec);
 113       vector<float> trueGammaDirXVec; 
 114       clusterCol->getParameters().getFloatVals("_trueGammaDirX",trueGammaDirXVec);
 115       vector<float> trueGammaDirYVec; 
 116       clusterCol->getParameters().getFloatVals("_trueGammaDirY",trueGammaDirYVec);
 117       vector<float> trueGammaDirZVec; 
 118       clusterCol->getParameters().getFloatVals("_trueGammaDirZ",trueGammaDirZVec);
 119       vector<float> trueGammaPosXVec; 
 120       clusterCol->getParameters().getFloatVals("_trueGammaPosX",trueGammaPosXVec);
 121       vector<float> trueGammaPosYVec; 
 122       clusterCol->getParameters().getFloatVals("_trueGammaPosY",trueGammaPosYVec);
 123       vector<float> trueGammaPosZVec; 
 124       clusterCol->getParameters().getFloatVals("_trueGammaPosZ",trueGammaPosZVec);
 125       vector<float> trueNeutralinoEnergyVec; 
 126       clusterCol->getParameters().getFloatVals("_trueNeutralinoEnergy",trueNeutralinoEnergyVec);
 127       vector<float> trueNeutralinoDirXVec; 
 128       clusterCol->getParameters().getFloatVals("_trueNeutralinoDirX",trueNeutralinoDirXVec);
 129       vector<float> trueNeutralinoDirYVec; 
 130       clusterCol->getParameters().getFloatVals("_trueNeutralinoDirY",trueNeutralinoDirYVec);
 131       vector<float> trueNeutralinoDirZVec; 
 132       clusterCol->getParameters().getFloatVals("_trueNeutralinoDirZ",trueNeutralinoDirZVec);
 133       vector<float> trueThetaDirVec; 
 134       clusterCol->getParameters().getFloatVals("_trueThetaDir",trueThetaDirVec);
 135       vector<float> truePhiDirVec; 
 136       clusterCol->getParameters().getFloatVals("_truePhiDir",truePhiDirVec);
 137 
 138       //fill cluster specific parameters
 139       for ( int clusterCounter=0; clusterCounter < _nCluster; clusterCounter++ ) {
 140 	Cluster* cluster = dynamic_cast<Cluster*>(clusterCol->getElementAt(clusterCounter));
 141 	
 142 	//check wether collection is empty
 143 	if (cluster == 0) {
 144 	  std::cout<<"WARNING: did not find cluster in event "<<_ievt<<std::endl;
 145 	  std::cout<<"\t KIT_ClusterWriteEngine will skip event "<<_ievt<<std::endl;
 146 	  return;
 147 	}
 148 
 149 	_nHitsPerCluster = nHitsVec[clusterCounter];
 150 	_dca             = dcaVec[clusterCounter];	  
 151 	_phiPos          = phiPosVec[clusterCounter];
 152 	_thetaPos        = thetaPosVec[clusterCounter];
 153 	_cosThetaPos     = cos(_thetaPos);
 154 	_energy   = cluster->getEnergy();
 155  	_thetaDir = cluster->getITheta();
 156 	_phiDir   = cluster->getIPhi();
 157 
 158 	/* DEBUG
 159 	std::cout<<" cluster "<<clusterCounter
 160 		 <<" with "<<_nHitsPerCluster<<" hits "
 161 		 <<" and "<<_energy<<" GeV "
 162 		 <<" at distance "<<_dca<<std::endl; 
 163 	*/
 164 
 165 	_trueDist             = trueDistVec[clusterCounter];
 166 	_trueGammaEnergy      = trueGammaEnergyVec[clusterCounter];
 167 	_trueGammaDirX        = trueGammaDirXVec[clusterCounter];
 168 	_trueGammaDirY        = trueGammaDirYVec[clusterCounter];
 169 	_trueGammaDirZ        = trueGammaDirZVec[clusterCounter];
 170 	_trueGammaPosX        = trueGammaPosXVec[clusterCounter];
 171 	_trueGammaPosY        = trueGammaPosYVec[clusterCounter];
 172 	_trueGammaPosZ        = trueGammaPosZVec[clusterCounter];
 173 	_trueNeutralinoEnergy = trueNeutralinoEnergyVec[clusterCounter];
 174 	_trueNeutralinoDirX   = trueNeutralinoDirXVec[clusterCounter];
 175 	_trueNeutralinoDirY   = trueNeutralinoDirYVec[clusterCounter];
 176 	_trueNeutralinoDirZ   = trueNeutralinoDirZVec[clusterCounter];
 177 	_trueThetaDir         = trueThetaDirVec[clusterCounter];
 178 	_truePhiDir           = truePhiDirVec[clusterCounter];
 179 
 180 	/* DEBUG
 181 	std::cout<<"true cluster distance "<<_trueDist
 182 		 <<" true energy "<<_trueGammaEnergy<<std::endl;
 183 	*/
 184  	const float* _clusterPos =  cluster->getPosition();
 185 
 186 
 187 	//tranform angle from rad to degree
 188 	_phiDir   = _phiDir*180/M_PI;
 189 	_thetaDir = _thetaDir*180/M_PI;
 190 	_phiPos   = _phiPos*180/M_PI;
 191 	_thetaPos = _thetaPos*180/M_PI;
 192 
 193 	vector<float> clusterShape = cluster->getShape(); 
 194 
 195 	float _clusterAxis[3];
 196 	_clusterAxis[0]   = clusterShape[0];   //cluster main principle axis 
 197 	_clusterAxis[1]   = clusterShape[1];   //2nd cluster principle axis
 198 	_clusterAxis[2]   = clusterShape[2];   //3rd cluster principle axis
 199 
 200 	//fill branches
 201 	_clusterFill.energy[clusterCounter]          = _energy;
 202 	_clusterFill.nHitsPerCluster[clusterCounter] = _nHitsPerCluster;
 203 	_clusterFill.clusterPos[clusterCounter][0]   = _clusterPos[0];
 204 	_clusterFill.clusterPos[clusterCounter][1]   = _clusterPos[1];
 205 	_clusterFill.clusterPos[clusterCounter][2]   = _clusterPos[2];
 206 	_clusterFill.thetaDir[clusterCounter]        = _thetaDir; 
 207 	_clusterFill.phiDir[clusterCounter]          = _phiDir; 
 208 	_clusterFill.cosThetaPos[clusterCounter]     = _cosThetaPos; 
 209 	_clusterFill.phiPos[clusterCounter]          = _phiPos; 
 210 	_clusterFill.dca[clusterCounter]             = _dca;
 211 	_clusterFill.clusterAxis[clusterCounter][0]  = _clusterAxis[0];
 212 	_clusterFill.clusterAxis[clusterCounter][1]  = _clusterAxis[1];
 213 	_clusterFill.clusterAxis[clusterCounter][2]  = _clusterAxis[2];
 214 
 215 	_clusterFill.trueDist[clusterCounter]        = _trueDist;
 216 	_clusterFill.trueGammaEnergy[clusterCounter] = _trueGammaEnergy;
 217 	_clusterFill.trueGammaDirX[clusterCounter]   = _trueGammaDirX;
 218 	_clusterFill.trueGammaDirY[clusterCounter]   = _trueGammaDirY;
 219 	_clusterFill.trueGammaDirZ[clusterCounter]   = _trueGammaDirZ;
 220 	_clusterFill.trueGammaPosX[clusterCounter]   = _trueGammaPosX;
 221 	_clusterFill.trueGammaPosY[clusterCounter]   = _trueGammaPosY;
 222 	_clusterFill.trueGammaPosZ[clusterCounter]   = _trueGammaPosZ;
 223 	_clusterFill.trueNeutralinoEnergy[clusterCounter] = _trueNeutralinoEnergy;
 224 	_clusterFill.trueNeutralinoDirX[clusterCounter]   = _trueNeutralinoDirX;
 225 	_clusterFill.trueNeutralinoDirY[clusterCounter]   = _trueNeutralinoDirY;
 226 	_clusterFill.trueNeutralinoDirZ[clusterCounter]   = _trueNeutralinoDirZ;
 227 	_clusterFill.trueThetaDir[clusterCounter]    = _trueThetaDir;
 228 	_clusterFill.truePhiDir[clusterCounter]      = _truePhiDir;
 229 
 230       }// end loop over all cluster
 231 	
 232       //add some event specific parameters
 233       _etMiss       = evt->getParameters().getFloatVal("_Et_miss");
 234       _esum_calos   = evt->getParameters().getFloatVal("_Esum_calos");
 235       _nhits_calos  = evt->getParameters().getIntVal("_Nhits_calos");
 236       _esum_ECAL    = evt->getParameters().getFloatVal("_Esum_ECAL");
 237       _nhits_ECAL   = evt->getParameters().getIntVal("_Nhits_ECAL");
 238       _esum_HCAL    = evt->getParameters().getFloatVal("_Esum_HCAL");
 239       _nhits_HCAL   = evt->getParameters().getIntVal("_Nhits_HCAL");
 240       _trueNcluster = evt->getParameters().getIntVal("_trueNcluster");
 241       _nGamma       = evt->getParameters().getIntVal("_nGamma");
 242 
 243       _clusterFill.nCluster     = _nCluster;             
 244       _clusterFill.etMiss       = _etMiss; 	
 245       _clusterFill.esumCalos    = _esum_calos; 		
 246       _clusterFill.nHitsCalos   = _nhits_calos; 	
 247       _clusterFill.esumECAL     = _esum_ECAL; 		
 248       _clusterFill.nHitsECAL    = _nhits_ECAL; 	
 249       _clusterFill.esumHCAL     = _esum_HCAL; 		
 250       _clusterFill.nHitsHCAL    = _nhits_HCAL; 	
 251       _clusterFill.trueNcluster = _trueNcluster; 	
 252       _clusterFill.nGamma       = _nGamma; 	
 253       _clusterFill.iEvt         = _ievt;
 254 
 255       _ievt ++;
 256 
 257     }//try
 258     
 259     catch ( DataNotAvailableException err ) {
 260       cout <<  "KIT_ClusterWriteEngine WARNING: Collection "<< _clusterColName
 261 	   << " not available in event "<< evt->getEventNumber() << endl;
 262     }//catch
 263 
 264   }//fillVariables
 265 
 266 }//namespace marlin
 267 

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 14:46:21, 14.1 KB) [[attachment:ClusterWriteEngine.cc]]
  • [get | view] (2010-03-23 14:47:51, 4.6 KB) [[attachment:ClusterWriteEngine.hh]]
  • [get | view] (2010-03-23 14:57:17, 17.4 KB) [[attachment:eventSelection.c]]
  • [get | view] (2010-03-23 15:40:25, 7.3 KB) [[attachment:lifetime.c]]
 All files | Selected Files: delete move to page copy to page

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