All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups
AddClusterProperties.cc
Go to the documentation of this file.
1 #include "AddClusterProperties.h"
2 #include "WeightedPoints3D.h"
3 #include <iostream>
4 
5 #include <EVENT/LCCollection.h>
6 #include <EVENT/ReconstructedParticle.h>
7 #include <EVENT/Cluster.h>
8 #include <IMPL/ClusterImpl.h>
9 #include <IMPL/ReconstructedParticleImpl.h>
10 
11 // ----- include for verbosity dependend logging ---------
12 #include "marlin/VerbosityLevels.h"
13 #include "marlin/Exceptions.h"
14 #include <vector>
15 
16 #ifdef MARLIN_USE_AIDA
17 #include <marlin/AIDAProcessor.h>
18 #include <AIDA/IHistogramFactory.h>
19 #include <AIDA/ICloud1D.h>
20 //#include <AIDA/IHistogram1D.h>
21 #endif // MARLIN_USE_AIDA
22 
23 
24 using namespace lcio ;
25 using namespace marlin ;
26 
27 
28 AddClusterProperties aAddClusterProperties ;
29 
30 
31 AddClusterProperties::AddClusterProperties() : Processor("AddClusterProperties") {
32 
33  // modify processor description
34  _description = "AddClusterProperties does whatever it does ..." ;
35 
36 
37  registerInputCollection( LCIO::RECONSTRUCTEDPARTICLE,
38  "PFOCollectionName" ,
39  "Name of the input PFO collection" ,
40  _PFOName ,
41  std::string("PandoraPFOs") ) ;
42 
43  registerInputCollection( LCIO::CLUSTER,
44  "ClusterCollection" ,
45  "Name of the Clusters input collection" ,
46  _clusterCollectionName ,
47  std::string("PandoraClusters") ) ;
48 
49 
50 }
51 
52 
53 
54 
55 
56 
57 
58 void AddClusterProperties::init() {
59 
60  streamlog_out(DEBUG) << " init called " << std::endl ;
61 
62  // usually a good idea to
63  printParameters() ;
64 
65  _nRun = 0 ;
66  _nEvt = 0 ;
67 
68 }
69 
70 
71 void AddClusterProperties::processRunHeader( LCRunHeader* /*run*/) {
72 
73  _nRun++ ;
74 }
75 
76 
77 
78 void AddClusterProperties::processEvent( LCEvent * evt ) {
79 
80 
81  int npoints_index=0;
82  int sum_wgt_index=0 ;
83  int sum_wgt2_index=0 ;
84  int sum_wgt4_index=0 ;
85  int ecal_index=0 ;
86  int hcal_index=0 ;
87  int yoke_index=0 ;
88  int lcal_index=0 ;
89  int lhcal_index=0 ;
90  int bcal_index=0 ;
91  StringVec shape_keys ;
92  LCCollection* clucol = NULL;
93  try{
94  clucol = evt->getCollection( _clusterCollectionName );
95  StringVec shapeParams ;
96  shapeParams = clucol->getParameters().getStringVals ("ClusterShapeParameters", shapeParams);
97  shapeParams.push_back("npoints") ;
98  shapeParams.push_back("sum_wgt") ;
99  shapeParams.push_back("sum_wgt^2") ;
100  shapeParams.push_back("sum_wgt^4") ;
101  clucol->parameters().setValues( "ClusterShapeParameters" , shapeParams ) ;
102  shape_keys = clucol->getParameters().getStringVals( std::string("ClusterShapeParameters"),shape_keys);
103  for ( unsigned kkk=0 ; kkk < shape_keys.size() ; kkk++ ) {
104  if ( shape_keys[kkk] == "npoints" ) { npoints_index = kkk ; }
105  if ( shape_keys[kkk] == "sum_wgt" ) { sum_wgt_index = kkk ; }
106  if ( shape_keys[kkk] == "sum_wgt^2" ) { sum_wgt2_index = kkk ; }
107  if ( shape_keys[kkk] == "sum_wgt^4" ) { sum_wgt4_index = kkk ; }
108  }
109  StringVec subDetectorNames;
110  subDetectorNames = clucol->getParameters().getStringVals ("ClusterSubdetectorNames", subDetectorNames);
111  for ( unsigned kkk=0 ; kkk < subDetectorNames.size() ; kkk++ ) {
112  if ( subDetectorNames[kkk] == "ecal" ) { ecal_index = kkk ; }
113  if ( subDetectorNames[kkk] == "hcal" ) { hcal_index = kkk ; }
114  if ( subDetectorNames[kkk] == "yoke" ) { yoke_index = kkk ; }
115  if ( subDetectorNames[kkk] == "lcal" ) { lcal_index = kkk ; }
116  if ( subDetectorNames[kkk] == "lhcal" ) { lhcal_index = kkk ; }
117  if ( subDetectorNames[kkk] == "bcal" ) { bcal_index = kkk ; }
118  }
119  }
120  catch( lcio::DataNotAvailableException& e )
121  {
122  streamlog_out(WARNING) << _clusterCollectionName << " collection not available" << std::endl;
123  clucol = NULL;
124  }
125 
126  if( clucol != NULL ){
127  int nclu = clucol->getNumberOfElements() ;
128  for(int i=0; i< nclu ; i++){
129  streamlog_out(DEBUG3) << "i cluster " << i << std::endl;
130  ClusterImpl* clu = dynamic_cast<ClusterImpl*>( clucol->getElementAt( i ) ) ;
131  CalorimeterHitVec hits=clu->getCalorimeterHits();
132  int np=clu->getCalorimeterHits().size();
133  if ( np == 0 ) {
134  streamlog_out(ERROR) << " Cluster with no hits. Input probably isnt a RECO-file, so just stop here ! " << std::endl;
135  throw StopProcessingException( dynamic_cast<marlin::Processor*>(this) );
136  }
137 
138  // get the hits:
139 
140  std::vector<double> ehit,xhit,yhit,zhit;
141  for (int ihit=0; ihit< np ; ihit++){
142  CalorimeterHit* calo_hit = dynamic_cast<CalorimeterHit*>(hits[ihit]);
143  ehit.push_back(calo_hit->getEnergy());
144  xhit.push_back(calo_hit->getPosition()[0]);
145  yhit.push_back(calo_hit->getPosition()[1]);
146  zhit.push_back(calo_hit->getPosition()[2]);
147  }
148 
149  // Analayse the cluster:
150 
151  WeightedPoints3D wgtp = WeightedPoints3D( int(hits.size()) , &ehit[0] , &xhit[0], &yhit[0], &zhit[0] );
152 
153  double* cog=wgtp.getCentreOfGravity();
154 
155  double *covv=wgtp.getCentreOfGravityErrors();
156  double cov[3][3]; for( int iii=0 ; iii<3 ; iii++ ) { for ( int jjj=0 ; jjj<3 ; jjj++ ) { cov[iii][jjj]=covv[jjj+iii*3] ; } } ;
157 
158  double *eval=wgtp.getEigenVal();
159  double *eval_err=wgtp.getEigenValErrors();
160 
161  double *evpv=wgtp.getEigenVecPolar();
162  double evp[2][3]; for( int iii=0 ; iii<2 ; iii++ ) { for ( int jjj=0 ; jjj<3 ; jjj++ ) { evp[iii][jjj]=evpv[jjj+iii*3] ; } } ;
163 
164  double *evpc=wgtp.getEigenVecCartesian();
165  double evc[2][3]; for( int iii=0 ; iii<2 ; iii++ ) { for ( int jjj=0 ; jjj<3 ; jjj++ ) { evc[iii][jjj]=evpc[jjj+iii*3] ; } } ;
166 
167  double* evpev=wgtp.getEigenVecPolarErrors();
168  double evpe[2][2][3]; for( int iii=0 ; iii<2 ; iii++ ) { for ( int jjj=0 ; jjj<2 ; jjj++ ) {for ( int kkk=0 ; kkk<3 ; kkk++ ) { evpe[iii][jjj][kkk]=evpev[kkk+jjj*3+iii*3*2] ; } } };
169 
170  double sum_wgtsqr = wgtp.getTotalSquaredWeight();
171  double sum_wgt4 = wgtp.getTotalQuarticWeight();
172  double sum_wgt = wgtp.getTotalWeight();
173  streamlog_out(DEBUG3) << "i Totals : " << sum_wgt << " " << clu->getEnergy() << std::endl;
174  // if ( wgtp ) delete wgtp ;
175 
176  // Pack/translate/cast/... :
177 
178  FloatVec PositionError ;
179  PositionError.push_back(cov[0][0]*sum_wgtsqr/(sum_wgt*sum_wgt));
180  PositionError.push_back(cov[0][1]*sum_wgtsqr/(sum_wgt*sum_wgt));
181  PositionError.push_back(cov[1][1]*sum_wgtsqr/(sum_wgt*sum_wgt));
182  PositionError.push_back(cov[0][2]*sum_wgtsqr/(sum_wgt*sum_wgt));
183  PositionError.push_back(cov[1][2]*sum_wgtsqr/(sum_wgt*sum_wgt));
184  PositionError.push_back(cov[2][2]*sum_wgtsqr/(sum_wgt*sum_wgt));
185  streamlog_out(DEBUG3) << "i PositionError : " << std::endl; streamlog_out(DEBUG3) << "i " ;
186  for ( int iii=0 ; iii < 6 ; iii++ ) { streamlog_out(DEBUG3) << PositionError[iii] << " " ;} ;streamlog_out(DEBUG3) << std::endl;
187 
188  float Position[3];
189  Position[0] = cog[0];
190  Position[1] = cog[1];
191  Position[2] = cog[2];
192  streamlog_out(DEBUG3) << "i Position : " << std::endl; streamlog_out(DEBUG3) << "i " ;
193  for ( int iii=0 ; iii < 3 ; iii++ ) { streamlog_out(DEBUG3) << Position[iii] << " " ;} ;streamlog_out(DEBUG3) << std::endl;
194 
195  float theta = evp[0][2];
196  float phi = evp[1][2];
197  streamlog_out(DEBUG3) << "i theta/phi : " << std::endl; streamlog_out(DEBUG3) << "i " ;
198  streamlog_out(DEBUG3) << theta << " " << phi << std::endl;
199 
200  FloatVec DirectionError ;
201  DirectionError.push_back(evpe[0][0][2]);
202  DirectionError.push_back(evpe[1][0][2]);
203  DirectionError.push_back(evpe[1][1][2]);
204  streamlog_out(DEBUG3) << "i DirectionError : " << std::endl; streamlog_out(DEBUG3) << "i " ;
205  for ( int iii=0 ; iii < 3 ; iii++ ) { streamlog_out(DEBUG3) << DirectionError[iii] << " " ;} ;streamlog_out(DEBUG3) << std::endl;
206 
207  FloatVec shape = clu->getShape() ;
208  shape.resize(shape_keys.size());
209  shape[npoints_index]=1.0*hits.size() ;
210  shape[sum_wgt_index]=sum_wgt;
211  shape[sum_wgt2_index]=sum_wgtsqr;
212  shape[sum_wgt4_index]=sum_wgt4;
213 
214 
215  float Eerror=clu->getEnergyError();
216  if ( Eerror == 0.0 ) {
217  // not set, so as per HLRWS:
218  float E=clu->getEnergy();
219  const FloatVec pec = clu->getSubdetectorEnergies();
220  float Eem = pec[ecal_index]+pec[lcal_index]+pec[bcal_index];
221  float Ehad = pec[hcal_index]+pec[yoke_index]+pec[lhcal_index];
222  Eerror=0.0;
223  if ( Eem/E < 0.95 ) {
224  Eerror=sqrt((0.6*sqrt(E))*(0.6*sqrt(E)) + (0.03*E)*(0.03*E));
225  } else {
226  Eerror=sqrt((0.17*sqrt(E))*(0.17*sqrt(E)) + (0.01*E)*(0.01*E));
227  }
228  clu->setEnergyError(Eerror);
229  streamlog_out(DEBUG3) << " energies : " << E << " " << Eem << " " << Ehad << " " << Eerror << " " << Eerror/sqrt(E) << std::endl;
230  }
231 
232  for ( int kkk = npoints_index ; kkk <=sum_wgt4_index ; kkk++ ) {if (shape_keys[kkk] != shape_keys[kkk] ) { streamlog_out(WARNING) << " shape_keys " << kkk << " is NaN " << std::endl ; }}
233  for ( int kkk = 0 ; kkk < 2 ; kkk++ ) {if ( Position[kkk] != Position[kkk] ) { streamlog_out(WARNING)<< " Position " << kkk << " is NaN " << std::endl ; }}
234  for ( int kkk = 0 ; kkk < 6 ; kkk++ ) {if ( PositionError[kkk] != PositionError[kkk] ) { streamlog_out(WARNING) << " PositionError " << kkk << " is NaN " << std::endl ; }}
235  for ( int kkk = 0 ; kkk < 3 ; kkk++ ) {if ( DirectionError[kkk] != DirectionError[kkk] ) { streamlog_out(WARNING) << " DirectionError " << kkk << " is NaN " << std::endl ; }}
236  if ( theta != theta ) { streamlog_out(WARNING) << " theta is Nan " << std::endl ; }
237  if ( phi != phi ) { streamlog_out(WARNING) << " phi is Nan " << std::endl ; }
238 
239  // add to cluster:
240 
241  clu->setShape( shape ) ;
242  clu->setPosition(Position);
243  clu->setPositionError(&PositionError[0]);
244  clu->setITheta(theta);
245  clu->setIPhi(phi);
246  clu->setDirectionError(&DirectionError[0]);
247 
248  if( streamlog_level(DEBUG2) ) {
249  // if ( streamlog::out.write<streamlog::DEBUG2>() ){
250  debuging(clucol,clu,cog,&cov[0][0],eval,eval_err,&evp[0][0],&evpe[0][0][0],&evc[0][0],np,sum_wgt,sum_wgtsqr,sum_wgt4);
251  }
252 
253  }
254  }
255 
256  // Now loop PFOs to add the error on P for neutral ones, from the cluster
257  // NB: not enough to just use the cluster COG + errors + error on E: Also must
258  // check if the neutral actually is a V0 or gamma-conversion, in which case the
259  // p should be from the fitted values (or at least the sum on the tracks)
260 
261  LCCollection* pfocol = NULL;
262  try{
263  pfocol = evt->getCollection( _PFOName );
264  }
265  catch( lcio::DataNotAvailableException& e )
266  {
267  streamlog_out(WARNING) << _PFOName << " collection not available" << std::endl;
268  pfocol = NULL;
269  }
270  if( pfocol != NULL ){
271  int npfo = pfocol->getNumberOfElements() ;
272  for(int i=0; i< npfo ; i++){
273  // get the PFO...
274  ReconstructedParticleImpl* part = dynamic_cast<ReconstructedParticleImpl*>( pfocol->getElementAt( i ) ) ;
275  // and its clusters ...
276  ClusterVec clusters=part->getClusters();
277  int nclu = part->getClusters().size();
278  TrackVec tracks=part->getTracks();
279  int ntrk = part->getTracks().size();
280  if ( part->getCharge() == 0 && nclu == 1 && ntrk == 0 ) {
281  // for now, only look at the clear-cut cases - see above
282  streamlog_out(DEBUG3) << "i normal, neutral pfo " << i << std::endl;
283  Cluster* clu= dynamic_cast<Cluster*>(clusters[0]) ;
284  const float* cogv = clu->getPosition();
285  float dist = 0.0;
286  for ( int iii=0 ; iii < 3 ; iii++ ) {
287  dist+=cogv[iii]*cogv[iii];
288  }
289  dist = sqrt(dist);
290  float Eclu = clu->getEnergy();
291  float E = part->getEnergy();
292  streamlog_out(DEBUG3) << " pfo/cluster E " << E << " / " << Eclu << std::endl;
293  float Eerror = clu->getEnergyError();
294  const double* mompfo=part->getMomentum();
295  double mom[4] ;
296  for ( int kkk=0 ; kkk<3 ; kkk++) {mom[kkk] = E*cogv[kkk]/dist;}
297  mom[3] = E ;
298  streamlog_out(DEBUG3) << "i clu-mom: " ; for (int kkk=0 ; kkk<4 ; kkk++ ) { streamlog_out(DEBUG3) << mom[kkk] << " " ;} streamlog_out(DEBUG3) << std::endl;
299  streamlog_out(DEBUG3) << "i pfo-mom: " ; for (int kkk=0 ; kkk<4 ; kkk++ ) { streamlog_out(DEBUG3) << mompfo[kkk] << " " ;} streamlog_out(DEBUG3) << std::endl;
300  part->setMomentum(mom);
301  FloatVec PositionError = clu->getPositionError();
302  double covmat[4][4];
303  int nnn = 0;
304  for ( int kkk=0 ; kkk < 3 ; kkk++ ) {
305  for ( int lll=0 ; lll<=kkk ; lll++) {
306  covmat[lll][kkk] = PositionError[nnn];
307  covmat[kkk][lll] = PositionError[nnn];
308  nnn++ ;
309  }
310  }
311  for ( int kkk=0 ; kkk < 4 ; kkk++ ) {
312  covmat[kkk][3] = 0. ;
313  covmat[3][kkk] = 0. ;
314  }
315  covmat[3][3] = Eerror*Eerror;
316  streamlog_out(DEBUG3) << "i covmat : " << std::endl;
317  for ( int kkk=0 ; kkk <4 ; kkk++ ) {
318  for ( int iii=0 ; iii <4 ; iii++ ) { streamlog_out(DEBUG3) << covmat[kkk][iii] << " " ;}
319  streamlog_out(DEBUG3) << std::endl;
320  }
321 
322  double dp_drE[4][4] ;
323  double prefact= E/(dist*dist*dist);
324  for ( int kkk=0 ; kkk < 4 ; kkk++ ) {
325  if ( kkk < 3 ) {
326  for ( int lll=0 ; lll< 3 ; lll++) {
327  if ( lll == kkk ) {
328  dp_drE[lll][kkk] = prefact*(dist*dist-cogv[kkk]*cogv[kkk]);
329  } else {
330  dp_drE[lll][kkk] = -prefact*cogv[kkk]*cogv[lll];
331  }
332  }
333  } else {
334  for ( int lll=0 ; lll< 3 ; lll++) {
335  dp_drE[lll][3] = 0 ;
336  }
337  for ( int lll=0 ; lll< 3 ; lll++) {
338  dp_drE[3][lll] = prefact*cogv[lll]*dist*dist/E ;
339  }
340  dp_drE[3][3] = prefact*dist*dist*dist/E;
341  }
342  }
343  streamlog_out(DEBUG3) << "i dp_drE : " << std::endl;
344  for ( int kkk=0 ; kkk <4 ; kkk++ ) {
345  for ( int iii=0 ; iii <4 ; iii++ ) { streamlog_out(DEBUG3) << dp_drE[kkk][iii] << " " ;}
346  streamlog_out(DEBUG3) << std::endl;
347  }
348  // propagate
349  double p_cov[4][4];
350  double tmp_mat [4][4];
351  for ( int iii =0 ; iii < 4 ; iii++ ) {
352  for ( int jjj=0 ; jjj < 4 ; jjj++ ) {
353  tmp_mat [iii][jjj] = 0. ;
354  for (int kkk=0 ; kkk < 4 ; kkk++ ) {
355  tmp_mat [iii][jjj] += covmat[iii][kkk]* dp_drE[kkk][jjj];
356  }
357  }
358  }
359  for ( int iii =0 ; iii < 4 ; iii++ ) {
360  for ( int jjj=0 ; jjj < 4 ; jjj++ ) {
361  p_cov[iii][jjj] = 0. ;
362  for (int kkk=0 ; kkk < 4 ; kkk++ ) {
363  p_cov[iii][jjj] += dp_drE[kkk][iii]*tmp_mat[kkk][jjj] ;
364  }
365  }
366  }
367  streamlog_out(DEBUG3) << "i pcov : " << std::endl;
368  for ( int kkk=0 ; kkk <4 ; kkk++ ) {
369  for ( int iii=0 ; iii <4 ; iii++ ) { streamlog_out(DEBUG3) << "i " << p_cov[kkk][iii] << " " ;}
370  streamlog_out(DEBUG3) << std::endl;
371  }
372  nnn=0;
373  float p_cov_v[10];
374  for ( int iii=0 ; iii < 4 ; iii++ ) {
375  for ( int jjj=0 ; jjj <= iii ; jjj++ ) {
376  p_cov_v[nnn] = p_cov[iii][jjj] ;
377  if ( p_cov_v[nnn] != p_cov_v[nnn] ) { streamlog_out(WARNING) << " p_cov_v " << nnn << " is NaN " << p_cov_v[nnn] << std::endl ; }
378  nnn++ ;
379  }
380  }
381 
382  part->setCovMatrix(p_cov_v);
383 
384  } else {
385  if ( part->getCharge() == 0 ) {
386  streamlog_out(DEBUG3) << " atypical neutral : nclu = " << nclu << " , ntrk = " << ntrk << std::endl;
387  }
388  }
389  }
390  }
391  //-- note: this will not be printed if compiled w/o MARLINDEBUG=1 !
392 
393  streamlog_out(DEBUG4) << " processing event: " << _nEvt << " which is event " << evt->getEventNumber()
394  << " in run: " << evt->getRunNumber() << std::endl ;
395 
396 
397 
398  _nEvt ++ ;
399 }
400 void AddClusterProperties::debuging(LCCollection* clucol ,ClusterImpl* clu,double* cog,double* cov,double* eval,double* eval_err,double* evp,double* evpe,double* evc,int np,double sum_wgt,double sum_wgtsqr,double sum_wgt4) {
401 
402  streamlog_out(DEBUG2) << " input: cog " << cog[0] << " " << cog[1] << " " << cog[2] << std::endl;
403  streamlog_out(DEBUG2) << " input: covariance mat " << cov[0*3+0] << " " << cov[0*3+1] << " " << cov[0*3+2] << std::endl;
404  streamlog_out(DEBUG2) << " input: covariance mat " << cov[1*3+0] << " " << cov[1*3+1] << " " << cov[1*3+2] << std::endl;
405  streamlog_out(DEBUG2) << " input: covariance mat " << cov[2*3+0] << " " << cov[2*3+1] << " " << cov[2*3+2] << std::endl;
406  streamlog_out(DEBUG2) << " input: eigenvals " << eval[0] << " " << eval[1] << " " << eval[2] << std::endl;
407  streamlog_out(DEBUG2) << " input: V(eigenvals) " << eval_err[0] << " " << eval_err[1] << " " << eval_err[2] << std::endl;
408  streamlog_out(DEBUG2) << " input: eigen vec 1 (p) " << evp[0*3+0] << " " << evp[1*3+0] << std::endl;
409  streamlog_out(DEBUG2) << " input: eigen vec 2 (p) " << evp[0*3+1] << " " << evp[1*3+1] << std::endl;
410  streamlog_out(DEBUG2) << " input: eigen vec 3 (p) " << evp[0*3+2] << " " << evp[1*3+2] << std::endl;
411  streamlog_out(DEBUG2) << " input: cov(eigen vec) 1 (p) " << evpe[0*2*3+0*3+0] << " " << evpe[1*2*3+1*3+0] << " " << evpe[1*2*3+0*3+0] << std::endl;
412  streamlog_out(DEBUG2) << " input: cov(eigen vec) 2 (p) " << evpe[0*2*3+0*3+1] << " " << evpe[1*2*3+1*3+1] << " " << evpe[1*2*3+0*3+1] << std::endl;
413  streamlog_out(DEBUG2) << " input: cov(eigen vec) 3 (p) " << evpe[0*2*3+0*3+2] << " " << evpe[1*2*3+1*3+2] << " " << evpe[1*2*3+0*3+2] << std::endl;
414  streamlog_out(DEBUG2) << " input: eigen vec 1 (c) " << evc[0*3+0] << " " << evc[1*3+0] << " " << evc[1*3+0] << std::endl;
415  streamlog_out(DEBUG2) << " input: eigen vec 2 (c) " << evc[0*3+1] << " " << evc[1*3+1] << " " << evc[1*3+1] << std::endl;
416  streamlog_out(DEBUG2) << " input: eigen vec 3 (c) " << evc[0*3+2] << " " << evc[1*3+2] << " " << evc[1*3+2] << std::endl;
417  streamlog_out(DEBUG2) << " input: E, np, sum, sum_sq, sum_4 " << clu->getEnergy() << " " << np << " " << sum_wgt << " " << sum_wgtsqr << " " << sum_wgt4 << std::endl;
418 
419 
420  // check reading back: get values now stored in clu
421  const float* cogv = clu->getPosition();
422  vector<double> cog_fr_clu ;
423  for (int iii=0; iii<3 ; iii++) { cog_fr_clu.push_back(cogv[iii]);}
424  FloatVec PositionError = clu->getPositionError();
425  StringVec shape_keys ;
426  shape_keys = clucol->getParameters().getStringVals( std::string("ClusterShapeParameters"),shape_keys);
427  FloatVec wgts_sumv = clu->getShape() ;
428  int npnt = 0 ;
429  double wgt_sum= 0.0;
430  double wgt_sq_sum= 0.0;
431  double wgt_4_sum= 0.0;
432  for ( unsigned kkk=0 ; kkk < shape_keys.size() ; kkk++ ) {
433  if ( shape_keys[kkk] == "npoints" ) { npnt = int(wgts_sumv[kkk]) ; }
434  if ( shape_keys[kkk] == "sum_wgt" ) { wgt_sum = wgts_sumv[kkk] ; }
435  if ( shape_keys[kkk] == "sum_wgt^2" ) { wgt_sq_sum = wgts_sumv[kkk] ; }
436  if ( shape_keys[kkk] == "sum_wgt^4" ) { wgt_4_sum = wgts_sumv[kkk] ; }
437  }
438 
439  double seen_covmat[3][3];
440  int nnn = 0;
441  for ( int kkk=0 ; kkk < 3 ; kkk++ ) {
442  for ( int lll=0 ; lll<=kkk ; lll++) {
443  seen_covmat[lll][kkk] = PositionError[nnn]*wgt_sum*wgt_sum/wgt_sq_sum;
444  seen_covmat[kkk][lll] = PositionError[nnn]*wgt_sum*wgt_sum/wgt_sq_sum;
445  nnn++ ;
446  }
447  }
448  vector<double> cov_fr_clu ;
449  for( int jjj=0 ; jjj<3 ; jjj++) {for (int iii=0; iii<3 ; iii++) { cov_fr_clu.push_back(seen_covmat[jjj][iii]);}}
450  // for (int iii=0; iii<6 ; iii++) { cov.push_back(PositionError[iii]);}
451  FloatVec DirectionError = clu->getDirectionError();
452  vector<double> thphcov;
453  for (int iii=0; iii<3 ; iii++) { thphcov.push_back(DirectionError[iii]);}
454  // a new WeightedPoints3D object, ctored with the cluster-shape params
455  WeightedPoints3D wgtp_readback = WeightedPoints3D( cog_fr_clu, cov_fr_clu , thphcov ,npnt, wgt_sum , wgt_sq_sum , wgt_4_sum);
456  // read back stuff just entered
457  double* cog_readback=wgtp_readback.getCentreOfGravity();
458  double* covv=wgtp_readback.getCentreOfGravityErrors();
459  int np_readback=wgtp_readback.getNumberOfPoints();
460  double sum_wgt_readback=wgtp_readback.getTotalWeight();
461  double sum_wgtsqr_readback=wgtp_readback.getTotalSquaredWeight();
462  double sum_wgt4_readback=wgtp_readback.getTotalQuarticWeight();
463  double cov_readback[3][3] ;
464  for( int iii=0 ; iii<3 ; iii++ ) { for ( int jjj=0 ; jjj<3 ; jjj++ ) { cov_readback[iii][jjj]=covv[jjj+iii*3] ; } } ;
465  //for( int iii=0 ; iii<3 ; iii++ ) { for ( int jjj=0 ; jjj<3 ; jjj++ ) { cov_readback[iii][jjj]=covv[jjj+iii*3]*sum_wgtsqr_readback/(sum_wgt_readback*sum_wgt_readback) ; } } ;
466  // now get things not actually entered, but calculated.
467  double* evpv=wgtp_readback.getEigenVecPolar();
468  double evp_readback[2][3] ;
469  for( int iii=0 ; iii<2 ; iii++ ) { for ( int jjj=0 ; jjj<3 ; jjj++ ) { evp_readback[iii][jjj]=evpv[jjj+iii*3] ; } } ;
470  double* evcv=wgtp_readback.getEigenVecCartesian();
471  double evc_readback[3][3] ;
472  for( int iii=0 ; iii<3 ; iii++ ) { for ( int jjj=0 ; jjj<3 ; jjj++ ) { evc_readback[iii][jjj]=evcv[jjj+iii*3] ; } } ;
473  double* eval_readback=wgtp_readback.getEigenVal();
474  double* eval_err_readback=wgtp_readback.getEigenValErrors();
475  // this should be a mix of input and calculated: [][][3] is entry, the others calcuated. The latter should be
476  // different from the exact values from the points, since the fourth moments are missing.
477  double* evpev=wgtp_readback.getEigenVecPolarErrors();
478  double evpe_readback[2][2][3] ;
479  for( int iii=0 ; iii<2 ; iii++ ) { for ( int jjj=0 ; jjj<2 ; jjj++ ) {for ( int kkk=0 ; kkk<3 ; kkk++ ) { evpe_readback[iii][jjj][kkk]=evpev[kkk+jjj*3+iii*3*2] ; } } };
480  streamlog_out(DEBUG2) << " read back: cog " << cog_readback[0] << " " << cog_readback[1] << " " << cog_readback[2] << std::endl;
481  streamlog_out(DEBUG2) << " read back: covariance mat " << cov_readback[0][0] << " " << cov_readback[0][1] << " " << cov_readback[0][2] << std::endl;
482  streamlog_out(DEBUG2) << " read back: covariance mat " << cov_readback[1][0] << " " << cov_readback[1][1] << " " << cov_readback[1][2] << std::endl;
483  streamlog_out(DEBUG2) << " read back: covariance mat " << cov_readback[2][0] << " " << cov_readback[2][1] << " " << cov_readback[2][2] << std::endl;
484  streamlog_out(DEBUG2) << " read back: eigenvals " << eval_readback[0] << " " << eval_readback[1] << " " << eval_readback[2] << std::endl;
485  streamlog_out(DEBUG2) << " read back: V(eigenvals) " << eval_err_readback[0] << " " << eval_err_readback[1] << " " << eval_err_readback[2] << std::endl;
486  streamlog_out(DEBUG2) << " read back: eigen vec 1 (p) " << evp_readback[0][0] << " " << evp_readback[1][0] << std::endl;
487  streamlog_out(DEBUG2) << " read back: eigen vec 2 (p) " << evp_readback[0][1] << " " << evp_readback[1][1] << std::endl;
488  streamlog_out(DEBUG2) << " read back: eigen vec 3 (p) " << evp_readback[0][2] << " " << evp_readback[1][2] << std::endl;
489  streamlog_out(DEBUG2) << " read back: cov(eigen vec) 1 (p) " << evpe_readback[0][0][0] << " " << evpe_readback[1][1][0] << " " << evpe_readback[1][0][0] << std::endl;
490  streamlog_out(DEBUG2) << " read back: cov(eigen vec) 2 (p) " << evpe_readback[0][0][1] << " " << evpe_readback[1][1][1] << " " << evpe_readback[1][0][1] << std::endl;
491  streamlog_out(DEBUG2) << " read back: cov(eigen vec) 3 (p) " << evpe_readback[0][0][2] << " " << evpe_readback[1][1][2] << " " << evpe_readback[1][0][2] << std::endl;
492  streamlog_out(DEBUG2) << " read back: eigen vec 1 (c) " << evc_readback[0][0] << " " << evc_readback[1][0] << " " << evc_readback[1][0] << std::endl;
493  streamlog_out(DEBUG2) << " read back: eigen vec 2 (c) " << evc_readback[0][1] << " " << evc_readback[1][1] << " " << evc_readback[1][1] << std::endl;
494  streamlog_out(DEBUG2) << " read back: eigen vec 3 (c) " << evc_readback[0][2] << " " << evc_readback[1][2] << " " << evc_readback[1][2] << std::endl;
495  streamlog_out(DEBUG2) << " read back: E, np, sum, sum_sq, sum_4 " << clu->getEnergy() << " " << np_readback << " " << sum_wgt_readback << " " << sum_wgtsqr_readback << " " << sum_wgt4_readback << std::endl;
496 }
497 
498 void AddClusterProperties::check( LCEvent * /*evt*/ ) {
499  // nothing to check here - could be used to fill checkplots in reconstruction processor
500 }
501 
502 
503 void AddClusterProperties::end(){
504 
505  // std::cout << "AddClusterProperties::end() " << name()
506  // << " processed " << _nEvt << " events in " << _nRun << " runs "
507  // << std::endl ;
508 
509 }
AddClusterProperties::processRunHeader
virtual void processRunHeader(LCRunHeader *run)
Called for every run.
Definition: AddClusterProperties.cc:71
AddClusterProperties::_clusterCollectionName
std::string _clusterCollectionName
Input collection name.
Definition: AddClusterProperties.h:71
AddClusterProperties::debuging
virtual void debuging(LCCollection *clucol, ClusterImpl *clu, double *cog, double *cov, double *eval, double *eval_err, double *evp, double *evpe, double *evc, int np, double sum_wgt, double sum_wgtsqr, double sum_wgt4)
Definition: AddClusterProperties.cc:400
AddClusterProperties
Example processor for marlin.
Definition: AddClusterProperties.h:33
AddClusterProperties::processEvent
virtual void processEvent(LCEvent *evt)
Called for every event - the working horse.
Definition: AddClusterProperties.cc:78
AddClusterProperties::end
virtual void end()
Called after data processing for clean up.
Definition: AddClusterProperties.cc:503
AddClusterProperties::init
virtual void init()
Called at the begin of the job before anything is read.
Definition: AddClusterProperties.cc:58
AddClusterProperties::check
virtual void check(LCEvent *evt)
Definition: AddClusterProperties.cc:498
sistrip::e
static const float e
Definition: PhysicalConstants.h:26
sistrip::StringVec
std::vector< std::string > StringVec
Definition: SiStripClus.h:56
aAddClusterProperties
AddClusterProperties aAddClusterProperties
Definition: AddClusterProperties.cc:28