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.You are not allowed to attach a file to this page.